Black hole hair formation in shift-symmetric generalised scalar-tensor gravity

A linear coupling between a scalar field and the Gauss-Bonnet invariant is the only known interaction term between a scalar and the metric that: respects shift symmetry; does not lead to higher order equations; inevitably introduces black hole hair in asymptotically flat, 4-dimensional spacetimes. Here we focus on the simplest theory that includes such a term and we explore the dynamical formation of scalar hair. In particular, we work in the decoupling limit that neglects the backreaction of the scalar onto the metric and evolve the scalar configuration numerically in the background of a Schwarzschild black hole or a collapsing dust star described by the Oppenheimer-Snyder solution. For all types of initial data that we consider, the scalar relaxes at late times to the known, static, analytic configuration that is associated with a hairy, spherically symmetric black hole. This suggests that the corresponding black hole solutions are indeed endpoints of collapse.


Introduction
A century after black holes and gravitational waves were first predicted as solutions to Einstein's equations, the LIGO Scientific and VIRGO collaborations reported the first direct observations of gravitational waves originating from coalescing black-hole binaries [1][2][3]. This remarkable discovery can also be considered as the first direct observation of black holes and has opened up an entirely new chapter in understanding and probing gravity in its strong-field regime [4,5]. It is possible that electromagnetic and gravitational wave observations of astrophysical black holes will reveal deviations from the (perturbed) Kerr geometry and allow us to infer the existence of a new fundamental field [6][7][8][9][10].
In general relativity, no-hair theorems have established that black holes are surprisingly simple objects, parametrized fully by only three (global) charges: their mass M , angular momentum J, and electromagnetic charge Q [11][12][13][14][15]; see e.g. Refs. [10,16,17] for recent reviews on the topic. It is well-known that black holes can have hair in the presence of Yang-Mills fields [18][19][20]. However, our focus here will be extensions of general relativity that involve a scalar field. In such theories no-hair theorems still exist. They are essentially a consequence of the fact that the equation where 2 is the curved spacetime d'Alembertian, admits only the trivial solution Φ = constant in an asymptotically flat region of spacetime that has a Killing horizon as an inner boundary [21]. This leads to the conclusion that stationary black hole solutions in scalar-tensor theories are the same as in general relativity. This result has been extended to scalars with nonlinear self-interactions in Ref. [22]. By means of field redefinitions and conformal transformations the applicability of the proof extends to the widest class of scalar-tensor theories that are quadratic in derivatives. Recent pedagogical reviews on no-hair theorems involving scalar fields can be found in Refs. [23,24]. As discussed in detail in Ref. [23], no-hair theorems rely on a number of assumptions such as: asymptotic flatness and stationarity of the spacetime, absence of matter, and the requirement that additional fields exhibit the same symmetries as the metric. The validity of these assumptions can be disputed. It is known that black holes can develop scalar hair if they have matter in their vicinity [25,26], if the scalar is complex and has a time-dependent phase [27][28][29], or if the asymptotics are cosmological or anti-de Sitter [30][31][32][33][34].
Still, the most obvious way to evade no-hair theorems is to consider a broader class of scalar-tensor theories in which the action contains terms with more than two derivatives. Horndeski [35] has pinned-down the most general scalar-tensor theory that leads to second-order field equations. The action coincides with that of generalized galileons, which have recently received much attention in cosmology; see e.g. Ref. [36] for a mathematical introduction and references therein for phenomenological applications. There is no no-hair theorem that applies to this general class of theories. Instead, a counter example has been known for quite some time. As shown in Ref. [37], an exponential coupling between the scalar and the Gauss-Bonnet invariant G = R 2 − 4 R ab R ab + R abcd R abcd leads to hairy black hole solutions. Such a coupling is known to be present in the low energy effective action of heterotic string theory [37][38][39][40] and can also arise in a dimensional reduction of Lovelock gravity [41,42]. Several studies of black holes in theories where a scalar is coupled to the Gauss-Bonnet invariant have followed [43][44][45][46][47].
An interesting subclass of Horndeski theories consists of the subset that satisfies shift symmetry, Φ → Φ + constant, as this symmetry prevents the scalar from acquiring a mass. This symmetry also excludes the exponential coupling e Φ G that led to the hairy solution of Ref. [37]. ‡ Indeed, it has been shown in Ref. [48] that in shift-symmetric Horndeski theories static, spherically symmetric, asymptotically flat black holes cannot have hair. However, as pointed out in Ref. [49], there is a shift-symmetric coupling that manages to circumvent this no-hair theorem. Since the Gauss-Bonnet invariant G is a total divergence, the linear coupling ΦG is invariant under shifts up to a boundary term. Assuming it has a canonical kinetic term, the scalar satisfies the equation of motion where λ is a coupling constant. Hence, the scalar is sourced by G, which contains the Kretschmann scalar R abcd R abcd . Since the latter in general does not vanish in a black hole spacetime, the scalar will be forced to have a nontrivial configuration. Note that another interesting way to circumvent the no-hair theorem of Ref. [48] is to allow Φ to have a linear dependence on Killing time [49][50][51]. However, we will not consider this option here. A solution to (2) that describes the scalar profile of a hairy, static, spherically symmetric, black hole has been obtained in Ref. [49] working perturbatively in the coupling λ. This matches the solutions found earlier in Refs. [44,52] using the same technique, but working with a more general coupling and applying a weak-field approximation for the scalar field as well. In Ref. [53] instead, a nonperturbative, numerical solution has been presented and compared in detail with the perturbative one. This numerical solution resembles strongly the one found in Ref. [37] for the exponential coupling. All of these solutions are static and are expected to be endpoints of gravitational collapse. Our main focus here is to present a first, preliminary exploration of whether this is indeed the case.
Our motivation is threefold: (i) These solutions constitute a two-parameter family, parametrized by the mass and the scalar charge of the black hole. However, generically the scalar is singular on the event horizon, unless the scalar charge and the mass satisfy a bond. Imposing regularity selects a one-parameter family and it is an open question whether solutions within this family are dynamically selected during collapse. (ii) Stellar configurations in the theory in question have been shown to have vanishing scalar monopole, i.e. the asymptotic fall-off for the scalar is necessarily faster than r −1 [54]. In contrast, in the known black hole solutions the scalar does exhibit an r −1 falloff. This implies that this monopolar component should develop during collapse. (iii) A scalar-tensor theory that evades no-hair theorems is expected to lead to detectable deviations from general relativity in the strong field regime. The first step towards confronting its prediction with observations is to understand black hole formation and evolution.
Our exploration will be preliminary because we will resort to the decoupling limit, i.e. we will neglect the scalar field's backreaction onto the geometry. This approximation reduces the problem to solving (2) on a background spacetime that is a solution to Einstein's equations, potentially with matter. We will consider two different backgrounds: a Schwarzschild black hole, previously considered in Ref. [55], and an Oppenheimer-Snyder spacetime [56], which is the simplest model of stellar collapse. The evolution of the scalar will correspond to the formation of scalar hair on these spacetimes.
The rest of this paper is organized as follows. In Section 2 we define the theory with an action and derive equations of motion. We also discuss the decoupling limit and perturbative solutions. In Section 3 we formulate the problem in a way suitable for numerical methods and in Section 4 we present our numerical results. Section 5 contains our conclusions.

Action and equations of motion
The action that we will consider here reads where κ = 16πG, S Ψ denotes the matter action, λ and µ are coupling constants, and G is the Gauss-Bonnet invariant The coupling µ might appear redundant, as it could be absorbed in a redefinition of Φ, but we choose to keep it for reasons that will become apparent in the next section. In the following we will employ geometric units G = 1 and c = 1.
Varying the action with respect to the metric g ab and the scalar field Φ yields their field equations If the matter stress-energy tensor is non-zero, as is the case for the collapsing, homogeneous dust star, we complement these by the conservation of the energy-momentum tensor and the continuity equation where E and u a are the rest-mass density and velocity of the dust. The canonical scalar field energy-momentum tensor is and the correction due to the Gauss-Bonnet term is The action can be straightforwardly generalised by introducing a potential for Φ, by generalising the coupling between Φ and G, etc., but here we will focus on the simplest case that inevitably leads to hairy black holes. Any of these generalisations would break shift symmetry for Φ.

The decoupling limit
We are interested in the dynamical development of scalar hair for black holes, so ideally we would like to study the evolution of the scalar field and its imprint on the spacetime geometry during gravitational collapse of a star and the formation of a black hole. However, for the sake of simplicity, we will consider a simpler problem, namely the evolution of the scalar field and hair formation in a given spacetime background.
This approximation can be formally derived from the original theory as a decoupling limit. Consider the field equations (5) and (6). By taking the limit µ → 0 one can turn off the backreaction of the scalar field on the metric and is left with together with (8). That is, the field equations for the metric reduce to Einstein's equations in the presence of matter while the scalar's equation of motion remains unaffected.

The nature and role of λ
It is important to stress that there are two distinct ways to view the theory (3) depending on the status of the coupling λ. If λ is taken to be a usual coupling constant, the action can be taken as exact and studied as a classical theory of gravity. If λ is instead considered to double as a book keeping parameter of an expansion, the action can be taken to describe some effective theory. In this case the theory is known to order λ only, and hence one can only trust solutions to this order.
To make this more concrete and rigorous let us define the dimensionless parameter ε = λ/l 2 , where l is a characteristic length scale, and consider the small coupling limit -as opposed to decoupling -where ε 1. In the effective action scenario one has to work perturbatively in ε. Thus, the solutions will be of the form where the pair (g ab , Φ 0 ) constitutes an exact, potentially dynamical solution of the system (5) and (6) for ε = 0 (or λ = 0), Here, G ab , T ab and 2 (0) denote the Einstein tensor, the canonical scalar field energymomentum tensor and the d'Alembertian constructed from the background fields (ḡ ab , Φ 0 ). One can then use the expansion to generate a solution at O (ε) by solving the equations where quantities with superscript (0) are constructed from the background metricḡ ab , and G (1) ab and T (1) ab are the Einstein tensor and the scalar's stress tensor to first order. Higher order corrections should be discarded because the theory is only known to O (ε).
Note that for this discussion µ has been taken to be O(1), as generically the decoupling limit has nothing to do with the small coupling limit we are discussing here. In fact, solutions with nontrivial Φ 0 will have nonvanishing T (0) ab and hence the scalar will have nonvanishing backreaction on the spacetime already at zeroth order in ε. Notably, stationary, asymptotically flat, black-hole spacetimes do have trivial Φ 0 .
To see this one needs to first consider (15) and (16). These are effectively the equations of general relativity coupled to a scalar field. Hence, provided that the scalar shares the symmetries of the metric, no-hair theorems [21] apply and dictate that the only vacuum solution is Φ 0 = constant and the spacetime is described by the Kerr geometry. With Φ 0 = constant, (17) and (18) become exactly the same as (11) and (12). Hence, the full solution at decoupling will match the leading order solution at small coupling for stationary, asymptotically flat, black hole spacetimes.
Another point we wish to clarify in this section is the role of λ within the decoupling limit. Consider the transformation Φ → λΦ. At the level of the action (3), this transformation allows one to effectively set λ to 1 by simply redefining µ. This does not affect the process of taking the decoupling limit, and hence λ becomes a redundant coupling at decoupling. The same can be seen at the level of the field equations. At decoupling µ → 0, λ and Φ are entirely absent from (5). The transformation Φ → λΦ makes λ drop out from (6) as well. Clearly, instead of generating solutions for different values of the coupling constant one can select a specific λ and then obtain the remaining solutions simply by rescaling Φ. Hence, from now on we will just set the dimensionless coupling λ/M 2 = 1.

The late-time behaviour of the scalar field
The fact that the full solution at decoupling matches the leading order solution at small coupling for any stationary, asymptotically flat, black hole spacetime is particularly relevant to our work. Static, spherically symmetric, asymptotically flat solutions to the theory in action (3) have been studied in Refs. [49,53]. The small coupling solution is known analytically to quadratic order and it is unique. The leading order part of this solution, i.e. the scalar configuration on a Schwarzschild background, will be the exact, static, asymptotically flat solution at decoupling.
Below, we will use this scalar field profile to benchmark our numerical simulations in the decoupling limit at late times, when the field has settled down to a time-independent state. As explained in more detail in Section 3.2, we numerically evolve the background spacetimes using puncture coordinates [57][58][59][60] denoted as (t, r, θ, φ). At late times, these evolutions yield the well-known trumpet slices of the Schwarzschild spacetime [61][62][63]. However, because the metric functions in this slicing are not known in analytic form, here we instead employ isotropic coordinates (t S , ρ, θ, φ). The two coordinates systems agree within 0.1% at late times and for radii ρ ≥ 10M and r ≥ 10M , as we have explicitly verified in Appendix A.1. Hence, it is convenient to have the scalar profile of the analytically known static solution in isotropic coordinates. Instead of starting from the solution as given in Ref. [49] and perform a coordinate transformation, we prefer to rederive the solution in the desired coordinate system.
The Schwarzschild metric in isotropic coordinates is given by where the conformal factor and lapse function are and the horizon corresponds to ρ H = M/2. The scalar field equation (12) then reads Direct integration leads to a solution with two integration constants. Fixing them by demanding (i) regularity at the horizon, and (ii) Due to shift symmetry we can always set Φ ∞ = 0 without loss of generality. This agrees with the solution given in Ref. [49] after applying the coordinate transformation r = ψ 2 ρ, wherer denotes the areal radius coordinate.

Spacetime split revisited
Since we plan to numerically evolve the system of field equations (11) and (12), we will perform the ADM-York decomposition common in numerical relativity [64][65][66][67]. To this end we foliate the 4-dimensional manifold (M, g ab ) into a set of spatial hypersurfaces (Σ t , γ ij ) labelled by a time parameter t. We introduce the unit timelike vector n a orthogonal to the hypersurfaces with norm n a n a = −1. It can be expressed in terms of the lapse function α and shift vector β i as The 3-metric γ ab = g ab + n a n b acts as a projection operator with γ a b n b = 0 by construction. Then, the line element takes the form In the following we denote the covariant derivative and Riemann tensor with respect to the 3-metric as D i and R i jkl , while the extrinsic curvature is where L n is the Lie-derivative along n a .

Background spacetime
Dynamics in the decoupling limit boil down to solving the scalar's equation (12) on a given background. In vacuum the Schwarzschild solution is the unique solution of (11) under the assumption of spherical symmetry. Hence, without adding any matter fields one can simply study the evolution of the scalar on the Schwarzschild geometry. Simple as this setup might be, it can still capture the most important aspects of the problem that we are trying to understand here, as one can still use it to model the formation of a nontrivial scalar configuration on a black hole spacetime. Moreover, one can check if the evolution indeed has the static solution discussed in Sec. 2.4 as an endpoint. Hence, this will be the first case of background we will consider, revisiting the results of Ref. [55]. Once matter fields are included one needs to solve (11) together with (8). The solution of the system (8) and (11) can be taken to represent a star collapsing to form a black hole, while the solution of (12) will represent the time evolution of Φ during and after the formation of the black hole. Here we will use the simplest spacetime that can be thought of as representing idealised stellar collapse: the Oppenheimer-Snyder solution [56]. This is an analytic model of a homogeneous dust star collapsing into a black hole in spherical symmetry. Despite the fact that modelling matter as dust neglects important phenomena, such as the effect of radiation for instance, it seems to be an adequate approximation for our purposes. Since we are working in the decoupling limit, details regarding the structure of the matter configuration should not be particularly important for the behaviour of the scalar field, which is what is of interest here.
In its exterior the Oppenheimer-Snyder collapse is described by the Schwarzschild solution in Schwarzschild coordinates (t S ,r, θ, φ) with its surface located atr B . The interior of the star is given by a closed Friedmann metric Let us denote the surface of the star as χ B in the coordinate system (τ, χ, θ, φ). These can be related to the conformal time η ∈ (−π, 0) via Continuity of the metric is ensured by matching the circumference on the boundary between the star's interior and exterior regions, namely where ψ is the conformal factor and ρ the isotropic radial coordinate. The initial scale factor a B and value of χ B can be expressed in terms of the areal radiusr , Although both the Schwarzschild and the Oppenheimer-Snyder solution are known explicitly in the coordinate systems used above, these forms are not particularly suitable for our numerical simulations. In each case, one needs to introduce a foliation that penetrates the black hole horizon and at the same time allows us to continue the simulations after a black hole forms. To achieve this, instead of attempting to explicitly rewrite the solutions in a suitable foliation, we prefer to generate them numerically using 1 + log-slicing condition [68]. Details on the numerical evolution of the background spacetimes, including the prescription of initial data, can be found in Appendix B.

Scalar field evolution
In order to evolve the scalar field equation (12) in any of the two backgrounds we first re-write it as a time-evolution problem. Therefore, we introduce the scalar's conjugate momentum This definition immediately provides an evolution equation for the scalar field whereas the 3 + 1-decomposition of (12) yields the momentum's evolution. Hence, one has the set of equations where L β is the Lie-derivative along the shift vector, D i is the covariant derivative w.r.t. the 3-metric, and K is the trace of the extrinsic curvature. Since we are working in the decoupling limit, the Gauss-Bonnet invariant G depends only on the background geometry.
The system of evolution equations (33) determines the scalar field dynamics in 3+1 dimensions. They need to be supplemented with a set of initial conditions (Φ, Π)| t=0 , and we will specify two different types. Initial Data 1: The first set of data is for the trivial field configuration This simple setup already leads to interesting results. In particular, it demonstrates excellently that the scalar has to develop a nontrivial profile even if it is assumed to be trivial initially, as it is sourced by the Gauss-Bonnet invariant. Initial Data 2: The second type of initial data is a scalar field cloud anchored around the compact object given by where A 0 , r 0 and σ are the amplitude, location and width of the Gaussian. Σ(θ, φ) determines the angular distribution of Π 0 and is defined as a superposition of spherical harmonics. We focus on two specific choices, namely a spherically symmetric or "monopole" configuration with Σ(θ, φ) = Σ 00 ≡ Y 00 and a dipole configuration with Σ(θ, φ) = Σ 11 ≡ Y 1−1 − Y 11 . Unless denoted otherwise we will always set the dimensionless amplitude A 0 /M = 1 since it only leads to a re-scaling of the scalar in the decoupling limit.

Implementation
We have implemented the field equations in the decoupling limit (11) and (12) as part of the Lean code [69]. Originally based only on the Cactus Computational toolkit [70,71] and the Carpet mesh refinement package [72,73], Lean has now been adapted to the Einstein Toolkit [74][75][76]. We refer the interested reader to Ref. [77] for more details about the upgraded infrastructure. Lean has been extended to evolve additional bosonic fields coupled to gravity in Refs. [78][79][80].
To accomplish our present project we have not only incorporated new thorns into Lean to evolve the field equations in the decoupling limit but also new thorns capable of evolving the Oppenheimer-Snyder collapse. The details of this implementation are discussed in Appendix B. We analytically prescribe initial data for the background spacetime and the scalar fields. We carry out simulations of the general relativity background using the χ-version of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation [66,81,82] together with puncture coordinates [57][58][59][60]. We apply the method-of-lines to perform the evolutions, where spatial derivatives are typically approximated by fourth-or sixth-order finite difference stencils, and we use the fourthorder Runge-Kutta time integrator. In order to track the black-hole formation during the Oppenheimer-Snyder collapse and to obtain information about the black hole's properties we employ the apparent horizon finder AHFinderDirect [83,84]. At the outer boundary we employ Sommerfeld, i.e. radiative, boundary conditions as implemented in the Einstein Toolkit [74][75][76].
Our numerical domain typically contains 7 refinement levels, with the outer boundary located at 120M and resolution h/M = 1.0 on the outermost grid. This translates into the grid setup {(120, 24, 12, 6, 3, 1.5, 0.6), M/64} , in the notation of Section IIE of Ref. [69], with resolution h I /M = 1/64 on the innermost refinement level. While this setup is typically sufficent to obtain accurate numerical results, we found it necessary to push the outer boundary to 240M when we used Initial Data 2 and Π 0 was dipolar. In order to estimate the numerical error we have performed benchmark tests against the analytic solution (22) as well as convergence tests. The analysis, described in Appendix A and illustrated in figures 2 and A2, reveals a discretization error of ∆Φ/Φ 2% after an evolution time of about t/M ∼ 200. At late times, the numerical solution agrees with the analytic one within less than |Φ/Φ ana − 1| 1% for radii r/M ≥ 10.0.
To analyze the formation of nontrivial scalar hair we consider both the field's radial profile as well as its multipolar components extracted on spheres of fixed radii r ex as a function of time. In particular, we perform a multipole decomposition where Y lm (θ, φ) are the standard spherical harmonics. We also compute the canonical scalar field energy density E SF = T (Φ) ab n a n b , where the energy-momentum tensor T (Φ) ab is given in (9), and the quantity E GB = λG GB ab n a n b , which can be interpreted as the contribution of the Gauss-Bonnet coupling to the scalar's total energy density. Though we have verified that both quantities remain finite and smooth in all cases, we do not present them or discuss them in detail below. The scalar profile is always smooth and these quantities do not provide any additional information in the decoupling limit. §

Scalar field dynamics around Schwarzschild black holes
The specific choices of initial data for a scalar field in a Schwarzschild black hole background are summarized in table 1 and include both Initial Data 1 and Initial Data 2 with a spherically symmetric or dipole scalar configuration. We illustrate the time § Note that E GB /E SF ∼ λ. We have argued that λ is a redundant coupling within the decoupling limit approximation, but beyond decoupling the value of λ will determine the relative importance of the two contributions in the backreaction the scalar will have on the metric and control potential deviations between decoupling and small coupling solutions.  Figure 1. Radial scalar field profile, multiplied by the radius, at different instances of time in a Schwarzschild background for various initial data. Top-left: the field and its time derivative have been chosen to vanish initially but the scalar still develops a nontrivial profile as it is sourced by the Kretschmann scalar. Top-right: the field vanishes initially and the derivative Π 0 is given as a spherically symmetric Gaussian shell with parameters Σ(θ, φ) = Σ 00 , r 0 /M = 10 and σ/M = 1 in (35). Bottom-left: initially vanishing scalar with Π 0 given as a dipolar Gaussian shell with parameters Σ(θ, φ) = Σ 11 , r 0 /M = 10 and σ/M = 1 in (35). The type of data breaks spherical symmetry. We present the profiles along the θ = 0 axis. During the evolution the scalar field sheds off its dipole moment through quasi-normal ringing as shown in the right panel of figure 3 and settles down to a spherical profile. Bottom-right: Comparison of early (t/M = 10) and late time (t/M = 300) profiles in Schwarzschild geometry for various initial configurations. In all cases, at late times the field converges to the known analytic solution (22) with an asymptotic fall-off r|Φ| = constant, independently of the initial field content, as it is sourced by the Kretschmann scalar. evolution of the scalar field profile in detail for the various characteristic cases in figure  1. All four panels actually present the radial profile rescaled by the radius, r|Φ|, as this illustrates clearly the asymptotic behaviour. r|Φ| always remained smooth throughout the evolution. The presence of apparent kinks in the plots is due to the fact that we plot the absolute value of Φ and use a logarithmic scale. In all cases, the solutions approach r|Φ| = constant for large radii at late times. This agrees well with the leading order behaviour Φ ∼ 2λ M r + O 1 r 2 expected from the analytic solution (22). SBH ID 1 SBH ID 2: Σ 00 (r 0 = 6, σ = 1) SBH ID 2: Σ 00 (r 0 = 6, σ = 2) SBH ID 2: Σ 00 (r 0 = 10, σ = 1) SBH ID 2: Σ 11 (r 0 = 10, σ = 1) OS ID 1 OS ID 2: Σ 00 (r 0 = 6, σ = 2) OS ID 2: Σ 00 (r 0 = 10, σ = 1) Figure 2. The fractional deviation between the late-time numerical profile and the static, analytic solutions |Φ/Φ ana − 1| at late times t/M = 300 for different types of initial data. We see that the deviation remains below |Φ/Φ ana −1| 1.0% independent of the initial scalar field configuration. "SBH" refers to the Schwarzschild background whereas "OS" stands for the Oppenheimer-Snyder background.
In fact, irrespective of the choice of initial data the solution always converges to the known static, analytic scalar profile of (22) at late times. This can be seen in the the bottom-right panel of figure 1 where we show two time instances, one at early times and one at late times, for all the different initial data we considered. It can also be seen in more detail in figure 2, where we show the fractional deviation between the static, analytic solutions and the late-time numerical profile for all the cases we studied.
In the left panel of figure 3 we illustrate the time evolution of the l = m = 0 mode, constructed by projecting the scalar field onto the corresponding spherical harmonic as in (37), and measured at a fixed coordinate radius r ex /M = 40, for different initial configurations. The scalar field dynamics at early times are dominated by the specific Table 1. List of selected simulations performed in the background of a Schwarzschild black hole with mass parameter M = 1 and dimensionless coupling constant λ/M 2 = 1. The scalar field is initially either trivial, i.e. Initial Data 1, or Π 0 is given as a Gaussian shell with angular distribution Σ lm , located at r 0 /M and with width σ/M , i.e. Initial Data 2. initial setup, but after t/M ∼ 100 all types of data converge to the same solution. The case of an initially dipolar scalar field is of particular interest. The spherical component follows exactly that of an initially trivial field, while at the same time the scalar sheds off its dipolar component. This is depicted in the right panel of figure 3 where we present the l = m = 1 multipole of the scalar field extracted at r ex /M = 40. After the early time response, we clearly see the quasi-normal ringdown followed by the latetime tail. Moreover, we estimate the ringdown frequency of the numerical data to be M ω 11 = 0.292 − ı0.097. This is in excellent agreement, within 0.6%, with predictions in general relativity [85]. The oscillatory ringdown phase is followed by a power-law decay Φ ∼ t −5.2 . This tail, computed from our time-domain data, agrees within 4% with the theoretical prediction Φ ∼ t −(2l+3) = t −5 for l = 1 in general relativity [86][87][88].

Scalar field dynamics in Oppenheimer-Snyder background
Next we evolve the scalar's field equation (12) in the Oppenheimer-Snyder background.
In practice, we evolve the (8) and (11)  The different types of initial configurations, as indicated in the legend, determine the evolution at early times. As can be seen in all figures, after the stellar collapse the scalar eventually approaches the known analytic solution in Schwarzschild spacetime (22) and exhibits an r −1 asymptotic fall-off independent of the initial data.
in table 2.
In the top panels of figure 4 we present the radial profile multiplied by the radius, r|Φ|, at different instances in time for a scalar that is initially entirely trivial in Oppenheimer-Snyder backgrounds with r B /M = 5 and r B /M = 10 respectively. For the bottom-left panel of figure 4 we have used Initial Data 2, with Π 0 being a spherically symmetric Gaussian field in an Oppenheimer-Snyder background with r B /M = 5. In Table 2.
List of selected simulations performed in the background of an Oppenheimer-Snyder dust collapse with initial surface radius r B /M and mass parameter M = 1. The scalar field is initially either trivial, i.e. Initial Data 1, or Π 0 is given as a Gaussian shell with angular distribution Σ lm , located at r 0 /M and with width σ/M , i.e. Initial Data 2. Already during the stellar phase the Gauss-Bonnet invariant forces the scalar to develop a nontrivial profile even if it is trivial initially. After a horizon forms and the exterior spacetime settles to a Schwarzschild black hole, we recover the behaviour already discussed in the previous section 4.2: for all types of initial data the scalar approaches r|Φ| = constant asymptotically. More specifically, as shown in figure 2, the entire configuration approaches the known, static, analytic solution in which the spacetime is described by a Schwarzschild black hole and the scalar profile is that of (22).
In the bottom-right panel of figure 4 we present the l = m = 0 multipole of the scalar for various field configurations evolved in the background of a collapsing dust star with initial radius r B /M = 5. During the stellar phase, we clearly observe the excitation of a nontrivial scalar configuration, even in the case of trivial initial setup, induced by the Gauss-Bonnet invariant. After the collapse that occurs at around t AH /M = 21.25 in this spacetime, however, the time evolution of the scalar becomes insensitive to the original geometry and exhibits the same behaviour as in Schwarzschild. In particular, the field again converges to the same hairy black hole solution regardless of the inital setup.

Discussion
We have investigated the dynamical formation of scalar hair in the simplest theory that fashions a linear coupling between a scalar field and the Gauss-Bonnet invariant. This coupling is known to yield black hole hair in stationary configurations. In order to simplify our analysis we have worked in the decoupling limit, where the backreaction of the scalar onto the the spacetime geometry is neglected. This reduces the problem to solving the scalar's equations of motion in a background that is a solution to Einstein's equations.
We have considered two types of backgrounds, a Schwarzschild black hole (see also Ref. [55]) and the Oppenheimer-Snyder solution that describes the collapse of a dust star. We have explored several choices of initial data, including the case of a trivial scalar with vanishing time derivatives, and nontrivial cases where the initial scalar configuration is spherically symmetric or dipolar. In all cases the scalar configuration eventually relaxes to the known, analytic, static configuration. Beyond decoupling this configuration corresponds to a hairy black hole. Although not a rigorous mathematical proof, this is a strong indication that this solution is indeed the spherically symmetric endpoint of stellar collapse.
It should be stressed that the known static configuration of (22) is not the unique solution of (12). As discussed in Refs. [49,53], there exists a 2-parameter family of solutions that generically diverge on the horizon. Imposing regularity on the horizon implies a bond between the two parameters and selects a 1-parameter subclass. Although this appears to be a reasonable condition, one cannot know a priori if it constitutes tuning or if dynamical evolution naturally leads to this subclass. Our results imply the latter and, hence, clearly suggest that collapse will lead to the formation of hairy black holes in a theory where a scalar field couples linearly to the Gauss-Bonnet invariant.
Among the cases of initial data we studied, the one where the time derivative of the scalar field is initially given by a dipolar Gaussian shell is of particular interest because it does not respect spherical symmetry. Since our backgrounds are spherically symmetric, within the decoupling approximation, the Gauss-Bonnet invariant fails to source nonspherical contributions. Hence, they decay with rates predicted by general relativity. For example, we have seen that a dipole field loses its dipole mode via quasi-normal ringing with frequencies matching those predicted in general relativity. At the same time, the field does develop a spherical profile that converges to the known, analytic, static solution.
As is clear in our simulations, the dynamical behaviour of the scalar does not differ significantly when we switch from a Schwarzschild to an Oppenheimer-Snyder background. The early time behavior is affected, especially at small radii, and this can be attributed to the fact that the matter inside the star contributes to the curvature tensor, and effectively sources the scalar through the Gauss-Bonnet invariant. Our simulations clearly show that the scalar develops a nontrivial profile immediately during the stellar phase and well before an apparent horizon forms. This is inevitable because it is sourced strongly by the Gauss-Bonnet invariant that does not vanish at any stage of the evolution. It is worth pointing out that this is by no means in contradiction with the result of Ref. [54], where it has been shown that stationary solutions of (12) in an asymptotically flat spacetime without a horizon will have vanishing monopole. Firstly, our solutions are not stationary and when they approach stationarity at late times a horizon has already formed. Hence the result of Ref. [54] is not applicable here. Moreover, a vanishing monopole does not imply that the scalar configuration is trivial but only that its asymptotic fall-off is faster than r −1 . Remarkably, this is indeed the case for our solutions before they reach stationarity.
Our work has a number of exciting extensions. Within the decoupling limit, there are two natural next steps: to use a more realistic stellar collapse model than Oppenheimer-Snyder and to relax the symmetry assumptions of the background so as to allow for rotating black holes. Work in both directions is underway. It is also important to go beyond the decoupling approximation, as this would allow one to calculate the effect that the scalar field configuration has on the spacetime. As has been discussed in Ref. [53] for the static black hole case, the metric configuration changes significantly in the interior of the horizon once the scalar field's backreaction is taken into account. In that case the singularity has finite area and black holes have a mimimum mass [53]. Hence, it would be interesting to explore the dynamical formation of black holes beyond the decoupling approximation.
this is to compare directly the lapse function α, the shift vector β i , and the 3-metric γ ij as obtained by our simulations at late times with the same components as one can read them off the Schwarzschild metric in isotropic coordinates. In all cases we found agreement to within 0.1% for ρ/M ≥ 10 and r/M ≥ 10. An illustration for the case of the lapse is given in figure A1. This justifies using isotropic coordinates to perform the comparison between our numerical solutions and the known, static, analytic solution for the scalar profile.  Figure A1.
Top: radial profile of the lapse function at t/M = 200 for the numerically evolved Schwarzschild black hole (red dashed line) and Oppenheimer-Snyder collapse (blue dashed-dotted line) using the puncture gauge and its analytic value for a Schwarzschild black hole in isotropic coordinates (black solid line). Bottom: deviation between the numerically computed lapse function and its analytic value. As expected, the deviations are significant near the black hole and they drop below 0.1% in the far region. Note that the axis label "r/M " really stands for both the (dimensionless) isotropic radial coordinate ρ/M and the (dimentsionless) puncture radial coordinate r/M .

Appendix A.2. Convergence analysis
We estimate the numerical error by performing a convergence analysis exemplarily for evolutions of Initial Data 1 in (34)  We present the convergence plots for monopole mode Φ 00 extracted at r ex /M = 40 in figure A2. Specifically, we show the difference between the coarse and medium, and medium and high resolution runs, where we have rescaled the latter by Q 4 = 2.1 in the Schwarzschild case and Q 2 = 1.2 for the Oppenheimer-Snyder evolution with r B /M = 5.0 indicating, respectively, 4 th and 2 nd order convergence. We estimate the numerical error in the scalar field to be about ∆Φ 00 /Φ 00 2% after an evolution time of t/M ∼ 200.
Convergence analysis of the scalar field initialized by Initial Data 1 and evolved in the background of a Schwarzschild black hole (top) and Oppenheimer-Snyder collapse (bottom), measured at r ex /M = 40. The rescaling of the difference between the medium and high resolution (red dashed lines) indicate, respectively, 4 th and 2 nd order convergence.

Appendix A.3. Benchmark tests
As we have shown in the main body of the text, the scalar field numerically evolves towards the known, static, analytic solution given in (22) for all cases we have studied. Considering this as the expected behaviour, we can employ it to benchmark our numerical solution at late times. For this purpose, the fractional deviation between the numerical and analytic solutions depicted in figure 2 can be reinterpreted as a relative error. For all cases, this relative error is |Φ/Φ ana − 1| 1.0% at late times t/M = 300.

Appendix B. Evolution of the background spacetime
In order to generate the background spacetimes numerically we employ a coordinate gauge that allows for a smooth evolution across the horizon, namely puncture coordinates [57][58][59][60]. We use (standard) numerical relativity techniques that have been established over the last decade and details can be found, e.g., in Refs. [66,68,89]. We will briefly summarize them here assuming the presence of dust, i.e., a homogeneous, pressure-less perfect fluid. We recover the black hole evolution for vanishing matter energy-momentum tensor and matter quantities.
In order to numerically evolve the background spacetime it is convenient to perform a spacetime split as described in Sec. 3.1. Recall, that we foliate the 4-dimensional spacetime into a set of 3-dimensional spatial hypersurfaces whose geometry is encoded in the 3-metric γ ij , its embedding is described by the extrinsic curvature K ij defined in (26), and we introduce the unit normal vector n a orthogonal to the spatial slices. In the presence of dust we additionally decompose its velocity field according to where u a n a = −w and its spatial components are v a for which v a n a = 0 by construction. The normalization u a u a = −1 implies w 2 = 1+v i v i . It has proven convenient to redefine E * OS = wE OS . Performing the ADM-York decomposition of the conservation and continuity equations (B.2) yields the evolution equations for the energy density and velocity field of the collapsing dust shell The 3 + 1 split of Einstein's equations (B.1) yields the (gravity) evolution equations and constraint equations As indicated before, we recover those for a black hole by setting the energy density E * OS = 0.

Appendix B.3. Initial data
Let us first focus on the derivation of suitable initial data (γ ij , K ij , α, β i )| t=0 that are complemented by the appropriate matter quantities in the case of the Oppenheimer-Snyder spacetime. To construct initial configurations of the background geometry, we need to solve the constraints (B.9) and (B.10). Therefore, we start by performing the York-Lichnerowicz conformal decomposition [90,91] of the metric γ ij = ψ 4γ ij ,γ ij = Diag 1, ρ 2 , ρ 2 sin 2 θ , (B.11) where ψ andγ ij are the conformal factor and metric. After applying this decomposition, the spatial line element becomes dl 2 = ψ 4 dρ 2 + ρ 2 dΩ 2 . (B.12) Comparing with (19) we observe that this is nothing else but writing the initial spatial slices in isotropic coordinates (ρ, θ, φ). Bearing in mind the definition of the extrinsic curvature (26) we see immediately that K ij = 0 initially and, hence, the momentum constraint (B.10) is satisfied trivially. Instead, the conformal factor will be specific to the particular spacetime under consideration and is constructed by solving the Hamiltonian constraint (B.9). Before we derive it for each of the cases below, let us provide the last piece of information to complete our (more generic) initial conditions, namely those for the gauge functions.
Instead of taking the lapse function in isotropic coordinates, we initialize it either as α = 1 or as the pre-collapsed lapse, α = ψ −2 , that has proven necessary for numerically stable simulations of black-hole spacetimes [66]. The shift vector is β i = 0. The gauge functions will adjust themselves to puncture coordinates by virtue of their evolution equations (B.17) and (B.18).
Let us now derive the conformal factor for each of the background spacetimes.

. Evolution equations
To follow the time development of the background numerically we adopt a free-evolution scheme, i.e., we solve the constraints (B.9) and (B.10) only for the initial data, which is then evolved. Throughout the evolution we monitor the constraints and verify that they remain satisfied within the numerical accuracy. In practice, we evolve (B.5) -(B.8) using the BSSN formulation of Einstein's equations [81,82] which is known to yield numerically stable evolutions. In this approach, the dynamical variables are given by χ = γ −1/3 ,γ ij = χγ ij ,Γ i =γ jkΓi jk , The evolution equations are further modified by appropriate constraint addition, and their explicit form can be found, e.g., in Ref. [66].
The system of evolution PDEs (B.5) -(B.8) is closed by a suitable choice of coordinate conditions. In particular, we employ puncture coordinates [57][58][59][60] ∂ t α = β k ∂ k α − 2αK , (B.17) where we set the parameters to η β = 1/M and ζ Γ = 3/4. We illustrate the evolution of the Oppenheimer-Snyder collapse in figure B1 where we depict the lapse function and energy density at different instances in time. The results are in good agreement with those presented in Ref. [68].