{'fraction_non_alphanumeric': 0.038109431552054505, 'fraction_numerical': 0.014716840536512668, 'mean_word_length': 4.265097379851478, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Shrouded black holes in Einstein-Gauss-Bonnet gravity

We study black holes in a modified gravity scenario involving a scalar field quadratically coupled to the Gauss-Bonnet invariant. The scalar is assumed to be in a spontaneously broken phase at spatial infinity due to a bare Higgs-like potential. For a proper choice of sign, the non-minimal coupling to gravity leads to symmetry restoration near the black hole horizon, prompting the development of the scalar wall in its vicinity. The wall thickness depends on the bare mass of the scalar and can be much smaller than the Schwarzschild radius. In a weakly coupled regime, the quadratic coupling to the Gauss-Bonnet invariant effectively becomes linear, and no walls are formed. We find approximate analytical solutions for the scalar field in the test field regime, and obtain numerically static black hole solutions within this setup. We discuss cosmological implications of the model and show that it is fully consistent with the existence of an inflationary stage, unlike the spontaneous scalarization scenario assuming the opposite sign of the non-minimal coupling to gravity. Our model predicts the speed of gravitational waves to be extremely close to unity, - in a comfortable agreement with the observation of the GW170817 event and its electromagnetic counterpart.


Introduction
The growing number of black holes discovered by LIGO, VIRGO, and KAGRA collaborations [1,2,3,4] greatly improves the possibilities for strong field tests of gravity, which have been mainly limited to binary pulsars until the recent past. Such strong field tests are crucial, if a given modification of gravity is indistinguishable from General Relativity (GR) in the weak field regime, i.e., in the Solar system.
In this work we focus on scalar-tensor theories of gravity involving a scalar field ϕ coupled to the Gauss-Bonnet invariant R 2 GB . Such couplings have received a lot of attention in the literature: in particular, black hole solutions have been identified for a dilatonic coupling of the form ∼ e ϕ R 2 GB [5,6,7] as well as for the linear coupling ∼ ϕR 2 GB [8,9,10,11]. Recently, it has been shown that the quadratic scalar-Gauss-Bonnet coupling ∼ ϕ 2 R 2 GB leads to spontaneous scalarization of black holes and neutron stars for a certain choice of the sign of the coupling [12,13,14,15]. In a nutshell, scalarization occurs due to the instability of the trivial solution ϕ = 0 reproducing GR.
As a result, the scalar acquires a non-trivial profile in the vicinity of compact objects, which in turn provokes potentially testable departures from GR. Note that historically the first example of spontaneous scalarization has been given in Ref. [16], assuming the scalar coupling to the Ricci curvature, in which case scalarization takes place only for neutron stars. This model of scalarization is subject to stringent observational constraints [17] clearly demonstrating how deviations from GR caused by scalarization can be efficiently tested. Implications of scalar-Gauss-Bonnet coupling for black hole mergers [18,19,20,21] and binary pulsars [22] are currently under active investigation.
We find new solutions for static black holes in a model with the quadratic coupling ∼ ϕ 2 R 2 GB assuming a non-trivial bare (i.e., non-vanishing upon switching off gravity) potential V (ϕ). Studies of scalarization performed so far do not include the potential V (ϕ), or assume that it has one trivial minimum, which can be set at the origin ϕ = 0 with no loss of generality [23,24,25,26,27]. Unlike the previous studies, in this paper we consider a bare potential V (ϕ), which spontaneously breaks Here v is the expectation value of the scalar field ϕ at spatial infinity 1 . For a proper choice of the sign, the role of the scalar-Gauss-Bonnet coupling is to restore the symmetry in the vicinity of the black hole. Consequently, for a sufficiently large coupling to the Gauss-Bonnet invariant, the scalar ϕ sits close to zero near the horizon, and the overall scalar profile resembles that of a domain wall, which becomes more pronounced as one increases the bare mass of the field ϕ. We confirm this picture of shrouded black holes by running the full system of equations of motion involving the scalar field and metric potentials in Section 4. Before that, in Section 3 we treat the scalar profile in the test field approximation, in which case not only do numerical simulations become less cumbersome, but one can also estimate the profile analytically. In a weak coupling regime, the departure of the scalar from its expectation value v is always small, and one effectively recovers the model with a linear coupling to the Gauss-Bonnet invariant. In our case, the sign of the scalar-Gauss-Bonnet coupling is opposite compared to the scenario without the symmetry breaking potential [12,15] 2 .. This is advantageous in view of the following problem: in the spontaneous scalarization scenario of Refs. [12,15] the scalar-Gauss-Bonnet coupling mimics a large tachyonic mass in the accelerated Universe leading to a runaway solution inconsistent with existence of inflation [33]. With our choice of the sign of the coupling no instability during inflation occurs. Nevertheless, we face possible issues with post-inflationary cosmological evolution, albeit significantly alleviated. Indeed as we show in Section 5 there is an instability during the radiation-dominated stage. This instability, however, develops very slowly and thus does not spoil standard cosmological picture, provided that the reheating temperature is bounded as T reh 1 GeV. In Section 5, we also discuss stringent limits on model parameters following from the possibility of domain wall formation and Dark Matter overproduction, and point out the ways of avoiding these constraints. Finally, we show there that the speed of gravitational waves is well within the existing bound imposed from the simultaneous observation of the gravitational wave signal GW170817 [34] and its electromagnetic counterpart [35].

The model
We consider the following scalar-tensor model: where R 2 GB = R 2 − 4R µν R µν + R µνλρ R µνλρ is the Gauss-Bonnet invariant and is a parameter with dimension of length. We work in units 8πG N = 1, where G N is the Newton's constant, and we assume the mostly positive metric signature (−, +, +.+). The bare scalar potential V (ϕ) is taken to be of the form where v is a vacuum expectation value of the field ϕ, h 2 is the quartic self-interaction coupling constant, and −µ 2 is the tachyonic bare mass squared of the field ϕ. We assume that the scalar ϕ resides in the minimum of its potential V (ϕ) spontaneously breaking Z 2 -symmetry, i.e., (We choose the minimum with a positive sign with no loss of generality.) Note that the potential (2) fulfils V (ϕ ∞ ) = 0, so that it does not contribute to the cosmological constant Λ (cf. Ref. [36]). Neglecting Λ, we therefore have asymptotically Minkowski spacetime. In this work, we assume a particular sign of the scalar-Gauss-Bonnet coupling: With this choice, which is opposite compared to that of Refs. [12,15], Z 2 -symmetry, spontaneously broken at spatial infinity, can be restored in the strong gravity regime, in the vicinity of a compact object. Once this regime is realized, the scalar field acquires a non-trivial configuration in the form of a shell that shrouds a black hole. The main goal in this work is to find static black hole solutions in the model (1). Phenomenology of the model in the weak field regime, i.e., in the Solar system, has been considered in Ref. [37]. The bound on the model parameters inferred there using the Shapiro time delay measurements by the Cassini spacecraft reads 2 v 1.6 · 10 20 m 2 . This bound is very weak allowing to be orders of magnitude above the Schwarzschild radius of Solar mass objects.
We are interested in the situation where the scalar-Gauss-Bonnet coupling, negligible at spatial infinity, eventually overcomes the bare mass of the field ϕ, as one approaches the horizon. Otherwise, the field ϕ remains in the spontaneously broken phase. It is convenient to define the crossover radius r cross , where the scalar-Gauss-Bonnet coupling is of the order of the bare mass term ∼ µ 2 ϕ 2 , namely Given that R 2 GB (r) ≈ 12r 2 S /r 6 for r r S , we get Hereafter, the subscript S denotes values on the Schwarzschild radius. Formation of a scalar wall is possible provided that r cross r S , which translates into the bound on µ: or Note that for Solar mass objects and not much larger than r S , the scalar ϕ is ultra-light. So small masses µ are not at all exotic in physics, e.g., they are characteristic for axions and axion-like particles.

Test field approximation
Let us start by constructing a solution for the scalar field ϕ in the test field approximation, i.e., neglecting deviations of gravitational potentials from their GR values. Such a simplification is justified in theoretically well-motivated ranges of parameters r S (perturbative regime) and/or v 1. Furthermore, we will see that certain non-perturbative effects, which may have non-trivial phenomenological consequences, arise already in the test field approximation.
As we focus on static black holes in the present work, the line element in Schwarzschild coordinates is given by ds 2 = −e ν dt 2 + e λ dr 2 + r 2 · dθ 2 + sin 2 θdϕ 2 .
GR values of the gravitational potentials ν and λ are inferred from In the test field approximation, the equation of motion for the scalar ϕ reads where V ,ϕ ≡ dV dϕ and the prime denotes the derivative with respect to the radial coordinate r. Note that Eq. (11) is invariant under the change where L is a constant. The regularity of the solution on the boundary requires that all the terms ∼ 1/ 1 − r S r in Eq. (11) cancel out. This gives At infinity the boundary condition for the field ϕ is given by Eq. (3).
We solve Eq. (11) with the boundary conditions (3) and (13) using both the shooting method and the relaxation technique. The results are shown in Figs. 1 and 2 in terms of the field ϕ and its derivative ϕ , respectively. At this point, we have switched to dimensionless variables: Recall that we assume 8πG N = 1. Qualitatively, the scalar field profile is different in three cases: In all three cases, we assume that the crossover radius is larger than the Schwarzschild radius, so that Eq. (7) is fulfilled. Case I corresponds to the perturbative regime (with 2 being a perturbative parameter), i.e., the field ϕ remains close to its bare expectation value v all the way from infinity till the horizon surface, and one has ϕ → v in the limit → 0. Thus, for r S , one can effectively replace the quadratic coupling ∼ ϕ 2 R 2 GB with a linear one ∼ vϕR 2 GB . In cases II and III, there is a non-perturbative dependence of the scalar field solutions on the coupling constant 2 , manifested in particular by the strong exponential suppression of the field ϕ near the horizon surface. This confirms the main idea underlying the model: for sufficiently large 2 , the scalar-Gauss-Bonnet coupling leads to the symmetry restoration in the vicinity of a black hole. As a result, we observe a wall -the scalar field configuration interpolating between the broken and unbroken phase at spatial infinity and near the horizon, respectively. Cases II and III differ by the ratio of the crossover radius r cross and r ∼ µ −1 . Namely, the sphere r ∼ r cross is inside the sphere r ∼ µ −1 in case II, and vice versa in case III. In practice, this difference means that the walls are much steeper in case III. As µ increases, the wall becomes sharper and moves closer to the horizon.
In the remainder of this section, we discuss the analytical solution for the field ϕ in the test field regime. Details of the derivation of this solution can be found in the Appendix, and below we summarize the main results.
Case I. For r S , the departure of the field ϕ (including its value at the horizon ϕ S ) from the expectation value at inifnity v is very small and hence the effect of the Gauss-Bonnet term is rather mild, -in an agreement with the numerical result shown in Figs. 1 and 2. In particular, the scalar field profile at distances r r S is described by where A is an order one constant. Note that the coupling constant 2 enters the expression (16) in a perturbative fashion. This confirms that the regime r S corresponds to the weak field limit. Case II. In the limit of large and relatively small µ, the scalar has an exponentially suppressed value at the horizon ϕ S . At relatively small values of r, the scalar exponentially grows as  (14). The scalar field is normalized to its expectation value v at spatial infinity. The sets of parameters are chosen to cover the three cases outlined in Eq. (15): the orange line corresponds to case I (perturbative regime), the green line stands for case II, where the walls first appear; the walls become steeper in the large µ regime corresponding to case III depicted with the dark blue, aquamarine, and red lines.
We observe a non-perturbative dependence of the scalar profile on the coupling 2 . At r ∼ √ r S the exponential growth transitions to an exponentially slow relaxation towards the vacuum value v at infinity, Overall, the scalar profile has a configuration of a wall, which for r S r µ −1 can be described by the following analytical expression: where K 1/4 is a modified Bessel function of the order 1/4 of the second kind. Note that Eq. (19) reduces to expressions (17) and (18) in given regimes. The wall configuration can be characterized  (14). In the large regime, the scalar field derivative has peaks, which match the locations of the walls shown in Fig. 1; a wall width is estimated as a peak width. For increasing µ and fixed , a wall moves towards a horizon and its width shrinks, in accordance with Eqs. (22) and (23). by its location r wall and the width δ wall . It is not difficult to see that both are estimated as r wall ∼ δ wall ∼ √ r S , since √ r S is the only dimensionful parameter entering the expression (19). Case III. A similar story occurs for relatively large µ: the field ϕ starts from a tiny value ϕ S at the horizon and grows exponentially for r S r r cross , For r r cross , the growth is shut off: As µ increases, with held constant, the radius r cross shrinks, and consequently the scalar profile becomes more and more shallow. Eventually, when µ hits the upper bound in Eq. (7), one gets ϕ v. This could be expected from the beginning: when the bare potential V (ϕ) dominates over the scalar-Gauss-Bonnet coupling term all the way down to the horizon surface, the field ϕ remains in the broken phase. The configuration described by Eqs. (20) and (21) corresponds to a wall, which is much steeper than in case II, -in an agreement with the numerical results shown in Figs. 1 and 2. The location of the wall coincides with the peak of the derivative ϕ of the scalar profile given by Eqs. (20) and (21): see Eq. (6). The naive estimate for the wall width δ wall obtained by extrapolating expressions (20) and (21) to the wall center at r ∼ r cross , reads δ wall ∼ µ −1 . A more accurate analysis of the field behaviour around the location of the wall presented in the end of the Appendix gives a slightly different estimate: Note that µ can reach values as large as µ ∼ /r 2 S (for larger values the wall is not formed). In this limiting case, the wall is located very close to the horizon surface, and its width is a small fraction of the Schwarzschild radius.

Numerical analysis of the full system of equations
In this section, we go beyond the test field approximation and numerically solve the full set of metric and the field ϕ equations. As we are interested in static black hole solutions, we assume an ansatz for the metric of the form (9); however, in contrast to the previous section, here the metric is considered to be dynamical. For the metric potentials ν and λ the equations of motion are given by and The equation of motion for the field ϕ reads  1) and (2) obtained by numerically solving the full system of Eqs. (24), (25), and (27) are shown for different sets of model parameters. The scalar field ϕ is normalized to its expectation value v at spatial infinity. Dimensionless variables of Eq. (14) have been used.
Note that only three of the above four equations are independent due to the Bianchi identities.
Boundary conditions for gravitational potentials ν and λ at the horizon surface are given by The boundary condition for the scalar ϕ at spatial infinity is given by Eq. (3). To find the boundary condition for the scalar ϕ at the horizon surface, one performs the expansion near the horizon: where a 1 and b 1 are constants. Substituting Eq. (29) into Eq. (24) and keeping only the leading terms in 1/(r − r S ), we obtain Now substituting expansions (29) into Eq. (26), we get  1) and (2) is shown relative to its GR values e ν GR for the same sets of model parameters as in Fig. 3.
Finally, substituting Eqs. (29) and (31) into Eq. (27), we obtain Note that in the limit V (ϕ) → 0, we consistently reproduce the result of Refs. [13,15]. We solve the system of Eqs. (24), (25), and (27) with the boundary conditions (3), (28), and (32) numerically using both the shooting and relaxation methods. For the shooting technique we integrate from the horizon and use ϕ S as the shooting parameter. The value of ϕ S is adjusted such that the scalar field approaches the value v at infinity. On the other hand, when using the relaxation method, one specifies the boundary condition (32) directly, while the correct value of ϕ S is found automatically within the method itself.
The results of numerical simulations are demonstrated in Figs. 3, 4, and 5, which show the scalar profile, and the metric functions e ν and e λ relative to their GR values around Schwarzschild black holes, respectively. In particular, Fig. 3 confirms the formation of walls in the large regime identified in the test field approximation, see Fig. 1. The presence of a wall has a sizeable impact on the metric leading to a considerable deviation from GR potentials even for v 1. The deviation becomes more profound when the wall gets steep, while the parameters and v are kept fixed - There is a limitation on the model parameters, for which regular static black hole solutions exist. For a large enough backreaction of the scalar field, the equations of motion become illdefined at some r. More precisely, when the system of equations is written in a canonical form, the coefficient(s) of higher-order derivative variables are equal to zero at some r. To see how the problem manifests, one eliminates ν in Eq. (27) by making use of Eq. (26). As a result the equation for the field ϕ takes the form where by '...' we denote the terms with lower derivatives in ϕ. Clearly, no regular solution exists, if the coefficient in front of the second derivative ϕ is zero at some r. Using the Newtonian values of the gravitational potentials, i.e, assuming that the singularity in the equations of motion takes place reasonably far away from the horizon, one can express the condition for the existence of black hole solutions as follows: In the most interesting case of large and µ, the quantity on the l.h.s. reaches the maximum at r slightly larger than the crossover radius, where ϕ v, while the derivative ϕ is already negligible, i.e., ϕ 0. Using these values and Eq. (6), we arrive at the following constraint on the model parameters: This constraint limits the backreaction of the field ϕ on the gravitational potentials and hence departures from the GR metric, which nevertheless can be significant in the vicinity of the wall.

Cosmology
In this section, we briefly discuss the cosmological implications of the model given by Eqs. (1), (2), and (4). We show that consistency of the model with ΛCDM cosmology leads to important limitations on the early Universe evolution, in particular the reheating temperature T reh . Furthermore, possible problems at preheating as well as issues with domain walls and Dark Matter overproduction may invoke an extension of the model. (In)stability in the early Universe. On a cosmological background, the scalar-Gauss-Bonnet coupling feeds into the time-dependent mass of the field ϕ. This time-dependent mass can be tachyonic at certain epochs and reach very large values in the early Universe. As a result, one faces the risk of instability grossly altering the standard cosmological history [33]. Indeed, in the FLRW Universe expanding with the scalar factor a(t) the Gauss-Bonnet invariant can be written as where H ≡ȧ/a is the Hubble rate. Hence, for 2 < 0 assumed in Refs. [12,15] the field ϕ has a positive effective mass squared in the decelerating Universe and tachyonic mass in the accelerating Universe. The latter fact has dramatic consequences for inflation, when expansion occurs with a nearly constant Hubble rate H infl and the effective mass squared is given by Thus, for 2 < 0 there is a huge instability, which is inconsistent with the existence of an inflationary stage 3 . Recall that formation of a scalar wall takes place for large values of | | greatly exceeding H −1 infl , unless one assumes an extremely low scale inflation. In our case 2 > 0, and the situation is reversed: the instability develops in the decelerating Universe, and vice versa. In particular, due to the large positive mass squared (37), during inflation the field ϕ, including its perturbations, quickly relax to nearly zero values. This is true for any Fourier mode ϕ(k, t) with wavenumbers k/a · H 2 infl , in which case the mode decays as 1/a 3/2 . So, for these modes by the end of inflation one estimates whereā is the scale factor defined from k/ā ∼ · H 2 infl . For cosmologically interesting wavenumbers, these values ϕ(k, t e ) are very small even given that H −1 infl . Note that we are not interested in modes with higher k, such that the inequality k/a > l · H 2 infl is satisfied throughout inflation, since they are stable during the entire evolution of the Universe. Despite this important advantage over the case with 2 < 0, in our scenario the catastrophic instability may still take place during the radiation-dominated stage and possibly preheating, when the Universe is decelerating. Below we demonstrate this for the radiation-dominated stage and propose a solution to the problem. The equation of motion for the field ϕ at the times of interest is given byφ where we have neglected the potential term and used a(t) ∝ t 1/2 . The general solution of this equation is given by where I 1/4 and K 1/4 are modified Bessel functions of the order 1/4 of the first and second kind, respectively. We are interested in the growing part of the general solution ∼ C 2 . The coefficient C 2 is defined from the initial conditions for the field ϕ and its derivativeφ at the onset of the radiation-dominated stage, or the end of reheating, which we denote by t reh . At the times t t reh the solution for the field ϕ reads where we assumed /t reh 1. As it follows, even for tiny initial ϕ reh andφ reh and generically large /t reh , the field ϕ almost instantly reaches huge values, which are fatal for the existence of the radiation-dominated stage.
One obvious solution to this problem is to ensure that the exponent in Eq. (41) is cancelled by the exponentially suppressed amplitude of the scalar at t reh , i.e., to assume that /t reh ln (M Pl /ϕ reh ) 100. (In this section, we do not set 8πG N = 1, and restore the Planck mass M Pl ≈ 2.44 · 10 18 GeV.) The latter inequality here is fulfilled provided that the field ϕ, strongly decaying during inflation, does not experience a significant growth at preheating (see the comment at the end of this paragraph) 4 . The condition /t reh 100 has strong consequences for the reheating temperature of the Universe T reh . To demonstrate this, we use the expression for the Hubble rate during the radiation-dominated stage: where g * (T ) is the number of relativistic degrees of freedom. On the other hand, the Hubble rate at radiation-domination is given by H = 1/(2t). Equating the latter to the expression (43) and using /t reh 100, we obtain from Eq. (43): For the interesting values of of the order of the Schwarzschild radius of the Sun, the above constraint reads T reh 3 GeV .
We have set g * (T reh ) 30. Note that the upper bound on T reh here is completely reasonable: it is three orders of magnitude above the lower bound, T reh > 4 MeV at 95% CL [39]. Still, such low reheating temperatures imply that there is a long preheating stage separating the inflationary and radiation dominated eras. As the equation of state during preheating corresponds to deceleration, one again expects a runaway solution for the field ϕ. Whether or not this runaway can be achieved in a controllable way, remains an open issue. Let us make one comment in passing. The constraint (45) is valid independently of a particular choice of the potential V (ϕ) and only assumes 2 > 0. Namely, it is also applicable in the model with V (ϕ) = 0. While the latter scenario requires 2 < 0 for scalarization of static black holes, it has been shown in Refs. [29,30,31] that rapidly rotating Kerr black holes can develop scalar hair also in the case 2 > 0.
Domain walls. Formation of domain walls is common in models exhibiting spontaneous breaking of discrete symmetries [40], which is the case we are considering. As discussed above, the field ϕ quickly relaxes to zero during inflation. Thus, during the subsequent evolution the field ϕ occupies vacua with positive and negative expectation values, i.e., ϕ = ±v in causally disconnected patches with equal probability. Neighbouring regions, where the field ϕ takes opposite expectation values, 4 With this assumption, one has ϕ reh (k) ∼ ϕ(k, te). To estimate ϕ(k, te) from Eq. (38) and hence ϕ reh , we note that only the modes with the wavenumber k/a reh ∼ · H 2 reh or lower may experience instability during the radiationdominated stage. Such modes exponentially decay during inflation starting from the moment, when k a ∼ · H 2 infl . Hence, for the unstable mode with the maximum possible wavenumber one hasā ∼ (H reh /H infl ) 2 · a reh . The modes with lower k, which are also unstable, have smaller initial amplitudes; therefore we concentrate on the particular mode satisfying k/a reh ∼ · H 2 reh . Consequently, using Eq. (38), we find: Taking the values H infl ∼ 10 13 GeV, T reh 1 GeV, assuming ln( H reh ) < 100 and not extremely large a reh /ae, we get ln(M Pl /ϕ reh ) ∼ 100 -the value quoted in the main text. are separated by domain walls. Typically, the latter pose problems for ΛCDM model, because they tend to overclose the Universe and affect cosmic microwave background measurements in unacceptable manner. These considerations set a strong limit on the domain wall tension [40,41]: which can be converted into For example, assuming the expectation values v 0.1 M Pl , which are relevant from the viewpoint of shell formation, we obtain µ 2 · 10 −33 eV. With such a strong limit, the potential V (ϕ) can be neglected for the study of black holes of astrophysical size. This is not an option of our major interest in the present work. The restrictive bound (47) can be avoided, if one allows for a small explicit breaking of Z 2 -symmetry, in which case domain walls are dissolved at some point of cosmological history [40,42]. We assume that the small breaking of Z 2 -symmetry does not affect the analysis of black hole solutions, but this remains to be shown explicitly.
Dark Matter. Note that the field ϕ generically experiences oscillations around a minimum of the broken phase. The oscillations contribute to Dark Matter, and this suggests another way of constraining the model. They start at the time t * , when µ H * . Using Eq. (43), we obtain for the temperature of the Universe T * at this moment: At the onset of oscillations, their amplitude is naturally estimated as µ/h and then redshifts ∝ 1/a 3/2 . Thus, we can estimate the energy density of the Dark Matter component due to the field ϕ as The requirement that Dark Matter is not overproduced can be expressed as ρ osc (t eq ) ρ rad (t eq ), where t eq is the time of the matter-radiation equality, and ρ rad (t eq ) is the radiation energy density at this time. Using this inequality, entropy conservation in the comoving volume, which gives and Eq. (48), expressing the quartic self-interaction constant h 2 through the expectation value v = µ/h, we obtain where we have set g * (T * ) ≈ 3.4. The above bound would place a strong constraint on the range of parameters of our model. One can avoid the problem by adding a decay channel for oscillations of the field ϕ into lighter degrees of freedom. We postpone further discussion of the problem for future studies. Speed of gravitational waves. Finally, let us ensure that our model is consistent with the observations of the gravitational wave signal GW170817 [34] and the complementary gamma ray burst [35], which constrain the speed of gravitational waves to be |c T − 1| 10 −15 [43] disfavouring a number of dark energy and modified gravity scenarios [44,45,46]. In our model, the speed of gravitational waves is given by [47] The field ϕ is assumed to be in the minimum of spontaneously broken phase slowly drifting due to the Gauss-Bonnet coupling: Using Eq. (36), we estimate and |φ| ∼ H|φ|. The speed of gravitational waves is measured in the late Universe, when the Hubble rate is extremely small in the range of parameters interesting for black hole physics, i.e., H −1 and H µ. Hence, we haveφ 2 Hϕ|φ| ∼ ϕ|φ| and thus Taking, e.g., ∼ µ −1 ∼ r S ∼ 3 km (Schwarzschild radius of the Solar mass object), v ∼ 0.1M Pl , and the present day Hubble rate H 0 ≈ 70 km/(s · Mpc), we obtain Note that the speed of gravity in the vicinity of a black hole deviates from the speed of light by a greater amount, however, the region of space where this happens is negligibly small in comparison to the distance from the source of the event GW170817. We conclude that for any reasonable values of , µ, and v, the speed of gravitational waves is well within the existing observational bound 5 .

Discussions
In the present work we have constructed static black hole solutions in a model described by Eqs. (1), (2), and (4). Introducing the potential V (ϕ) is natural in the quadratic scalar-Gauss-Bonnet model, as the latter does not respect shift symmetry. Compared to previous studies, we considered the case in which the potential V (ϕ) spontaneously breaks Z 2 -symmetry, so that the scalar ϕ acquires a non-zero vacuum expectation value. Such a choice of potential is very common in particle physics, with the Higgs field being the most well-known example. For r S , the solution is in the perturbative regime, meaning that the departure of the scalar field ϕ from its expectation value is small and goes to zero in the limit 2 → 0. While the effect of the Gauss-Bonnet coupling is rather mild in this regime, the case of small is nevertheless interesting at least for two reasons. First, it allows for the perturbative treatment of black hole mergers, cf. Ref. [18]. Second, for r S we effectively recover the linear Gauss-Bonnet coupling ∼ vϕR 2 GB . Remarkably, the recovery occurs in the regime, where the linear coupling can be trusted according to the recent causality constraints [49]. Note that the coupling ∼ vϕR 2 GB is commonly discussed in the context of shift-symmetric models, while in our case the shift-symmetry is explicitly broken by the potential V (ϕ).
For r S , which is the most interesting case, the scalar profile around the black hole depends non-perturbatively on the coupling constant 2 . The Z 2 -symmetry, spontaneously broken at spatial infinity, is restored in the vicinity of the black hole, where the scalar field takes on exponentially small values. This gives rise to the main effect identified in this work: a black hole is shrouded by a wall for sufficiently large masses µ. Namely, the scalar field experiences a steep growth in a thin shell separating the domains, where it is in the broken phase and in the unbroken phase. This effect is already manifested in the test field approximation (see Figs. 1 and 2) and holds upon taking dynamics of the metric into account (see Fig. 3). As it can be seen from Figs. 4 and 5, the thinner the wall is, the larger the deviation from the Schwarzschild metric is for a fixed coupling constant 2 and expectation value v. Checking stability of black hole solutions we obtained in this paper, deserves a separate investigation in the line of Refs. [50,51]. In the future, it would be also interesting to see how the walls are manifested in the case of Kerr black holes and neutron stars. Furthermore, it is worth investigating implications of the walls for strong field tests of GR: black hole mergers and binary pulsars.
Finally, let us mention that besides black hole physics, early Universe cosmology represents a promising playground for studying GR in the large curvature regime. In the present work we briefly considered cosmological implications of the model at hand in Section 5. Unlike the commonly considered scenario of spontaneous scalarization with 2 < 0, no problem with catastrophic tachyonic instability at inflation arises in our case. The instability appears during the radiation-dominated stage, but it is harmless provided that the reheating temperature is constrained to be sufficiently low, T reh 1 GeV. As a result, one expects a rather long stage of preheating, which can be highly non-trivial due to the presence of the scalar-Gauss-Bonnet coupling and thus worth investigating. In addition, consistency with the ΛCDM model must be carefully studied in view of possible issues with Dark Matter overproduction and formation of domain walls. These challenges may necessitate a modification of our scenario. We leave this task for future work. It is also worthwhile noting, that in our model the speed of gravitational waves is extremely close to the speed of light, Eq. (56), being well within the present bound from the event GW170817 and its complimentary gamma ray burst. Acknowledgments S.R. is indebted to Leonardo Trombetta for useful discussions. The work of W. E. and S. R. is supported by the Czech Science Foundation GAČR, project 20-16531Y.

Appendix. Analytical solutions in test field approximation
In this Appendix, we construct an approximate analytical solution of Eq. (11) with the boundary conditions (3) and (13). First, we notice that in the range of masses µ interesting for our scenario and for ϕ S v, we have Thanks to the above inequality the last term in Eq. (13) can be dropped, so that the latter simplifies to which we assume in what follows. For generic model parameters, there are three separate ranges of r, where the solution for the field ϕ has a distinct behaviour: i) close to the horizon, ii) intermediate range, i.e., far from the horizon but inside the crossover radius given by Eq. (6), iii) above the crossover radius. For each range of r we find the general analytical solution, and then construct a full solution by matching the analytical solutions on the boundaries between the regions.
In the region of small r, i.e., close to the black hole horizon, Eq. (11) can be written as where we again took into account the inequality (57). The solution of the above equation respecting the boundary condition (58) reads where I 0 is the modified Bessel function of the first kind of the zeroth order.
In the intermediate region r S r r cross , one can replace the gravitational potential ν in Eq. (11) by its Newtonian value ν −r S /r and neglect the potential term with respect to the Gauss-Bonnet term using (6). As a result, we obtain ϕ + 2 r ϕ − 24 2 r 2 S r 6 ϕ = 0 .
The general solution of this equation reads where Γ is a Gamma function; I 1/4 and K 1/4 are modified Bessel functions of the order 1/4 of the first and second kinds, respectively. Finally, at very large distances r r cross , the equation of motion for the field ϕ is given bỹ whereφ ≡ ϕ − v and the Gauss-Bonnet term was neglected as compared to the potential term.
With the boundary condition ϕ → v in the limit r → ∞, one gets from the above equation: The general expressions (61), (64), and (67) contain arbitrary constants C 1 , C 2 , and D as well as the value of the scalar field at the horizon ϕ S . They are fixed by matching the solutions for different ranges of r. In what follows, we treat separately three cases outlined in Eq. (15). Case I. In this case, the arguments of the modified Bessel functions are small in the intermediate region r S r r cross . Therefore, we can use the following limiting expressions: and rewrite the general solution (64) as Matching the latter to Eq. (67) at r ∼ r cross gives and Now using the near horizon solution (61) and the small argument expansion of the modified Bessel function I 0 : Note that the coefficient in front of the second term here turns into zero at r cross , so that one can perform the expansion: 144l 2 r 2 S r 7 cross · (r − r cross ) .
The general solution of Eq. (88), which can be written as a linear combination of Airy functions, is not particularly illuminating. Without invoking the solution, however, one can deduce that there is only one dimensionful parameter in this regime given by the coefficient in front of (r − r cross ) in Eq. (89). This coefficient is thus identified as the inverse wall width, i.e., Using Eq. (6), we get the estimate (23) from the main text.