Spherical collapse in scalar-Gauss-Bonnet gravity: taming ill-posedness with a Ricci coupling

We study spherical collapse of a scalar cloud in scalar-Gauss-Bonnet gravity - a theory in which black holes can develop scalar hair if they are in a certain mass range. We show that an additional quadratic coupling of the scalar field to the Ricci scalar can mitigate loss of hyperbolicity problems that have plagued previous numerical collapse studies and instead lead to well-posed evolution. This suggests that including specific additional interactions can be a successful strategy for tackling well-posedness problems in effective field theories of gravity with nonminimally coupled scalars. Our simulations also show that spherical collapse leads to black holes with scalar hair when their mass is below a mass threshold and above a minimum mass bound and that above the mass threshold the collapse leads to black holes without hair, in line with results in the static case and perturbative analyses. For masses below the minimum mass bound we find that the scalar cloud smoothly dissipates, leaving behind flat space.

Introduction.-Astrophysicalphenomena in which gravity is at its strongest regime, such as the collapse of a compact object or the merger of coalescing binaries, cannot be adequately described with perturbative techniques.To model them in General Relativity (GR), one instead needs to formulate Einstein's equations as an initial value problem (IVP) and solve them numerically.Testing gravity in its nonlinear regime, to search for new physics beyond GR or the Standard Model of particle physics that could manifest, requires formulating IVPs and performing numerical simulations in such beyond-GR scenarios.To obtain predictions a well-posed IVP formulation should exist, i.e., the solution would be unique and depend continuously on the initial conditions.
A major obstacle is that in many interesting such scenarios, it is not obvious how to obtain a well-posed formulation of the IVP because the new physics drastically changes the structure of the field equations as partial differential equations.Scalar fields provide a characteristic example.No-hair theorems [1][2][3][4][5][6] dictate that scalars cannot leave an imprint on black holes (BHs) in most cases.All the known counterexamples (e.g.[7][8][9][10][11][12][13]) require coupling the scalar to the Gauss-Bonnet (GB) invariant, G = R µνρσ R µνρσ − 4R µν R µν + R 2 , where R µνρσ , R µν , and R being the Riemann tensor, Ricci tensor, and Ricci scalar respectively.Einstein's equations are quasi-linear, i.e. linear in the second derivatives of the metric, and this is a key property in what regards their well-posedness as an IVP.In contrast, G is clearly quadratic in the curvature and hence quadratic in the second derivatives of the metric.This places well-posedness in jeopardy.
One can circumvent this problem by working perturbatively in the coupling constants that control the deviations from GR. Assuming that the solutions are continu-ously connected to GR as the coupling goes to zero, one can generate solutions order by order [14][15][16][17][18][19].One disadvantage of this approach is that secular growth can drive it out of its range of validity [20].Another is that the perturbative treatment in the coupling is not suitable for capturing effects that involve nonlinearities in the new fields.
This last point is clearly demonstrated in the phenomenon of scalarization [21], first introduced in Ref. [22] for neutron stars.It was shown more recently that scalarization can affect BHs as well if a scalar field exhibits a suitable coupling with the GB invariant G [11,12].Scalarization can be understood as a linear tachyonic instability of the scalar that is eventually nonlinearly quenched.The instability occurs in GR spacetimes describing compact objects, it is controlled by interactions that are quadratic in the scalar [23] and, for BHs, it appears below a mass [11,12] or above a spin [24,25] threshold.The endpoint is a compact object with a nontrivial scalar configuration, whose properties are determined by nonlinear interactions of the scalar [26][27][28].Hence, the dynamics of theories exhibiting scalarization cannot be fully captured working perturbatively in the coupling constant [29][30][31][32][33].
Indeed, the study of the IVP for scalars nonminimally coupled to gravity, and in particular for the broad class of theories that lead to second-order equations described by the Horndeski action [34], has received a lot of attention recently [35][36][37][38][39][40][41][42][43].It was established in [44] that an appropriate formulation exists that renders the IVP well-posed in these theories in the weakly coupled regime, i.e., when the nonminimal couplings of the scalar remain small.Numerical studies, restricted so far to theories where the scalar couples only to the GB invariant (see [45] for an exception), have verified this result for small values of the coupling constant that controls this coupling [31,[46][47][48].For larger values of the coupling constant however well-posedness is eventually lost, as the equations tend to change character from hyperbolic to elliptic in a region of spacetime [31,[37][38][39][46][47][48][49].
Viewing these theories as nonlinear effective field theories (EFTs), and hence as products of some truncation of a more "fundamental" theory, a promising strategy could be to try to employ a method inspired by viscous relativistic hydrodynamics to render them predictive [50].How to import such a method to gravity theories is currently being explored [51][52][53][54][55][56][57][58].An interesting alternative would be to study whether additional couplings of the scalar, which one could expect to be there in an EFT, could actually improve the theories behaviour and lead to well-posed evolution.Indeed, it has been shown in Ref. [59] that an additional coupling of the type ϕ 2 R leads to an improvement of the hyperbolic nature of the equations for linear perturbations in scalarization theories that contain ϕ 2 G couplings.Interestingly, that same curvature interaction, ϕ 2 R, has been shown to resolve radial stability problems for scalarized BHs [28,59], to help evade binary pulsar constraints by suppressing neutron star scalarization [60], and to render GR a cosmological attractor without the need for the fine-tuning of the initial conditions [61], thereby making scalarization models compatible with late time cosmological dynamics.
Motivated by the above, we consider here a theory with quadratic scalar couplings to both the GB invariant and the Ricci scalar and study gravitational collapse of a scalar cloud.We show that the Ricci coupling does indeed improve significantly the dynamical behaviour of the theory and allows it to be predictive for values of the GB coupling that otherwise would have yielded an illposed IVP.Our numerical simulations allow us to track the formation of a scalarized BH for suitable initial data.For data that would have led to a formation of a BHs with mass smaller than the existence threshold of the theory we instead see the scalar cloud smoothly dispersing.The theory.-Weconsider the following action: where g = det(g µν ), X = −∇ µ ϕ∇ µ ϕ/2, and we use units of G = c = 1.The parameter β is a dimensionless coupling constant while α is of dimension length squared.
Varying action (1) with respect to the metric yields, where and with respect to the scalar fields, gives where □ := ∇ µ ∇ µ , G µν is the Einstein tensor, and ϵ σγλτ is the Levi-civita totally anti-symmetric tensor density.
Evolution in spherical symmetry.-Inthis paper, we focus on the study of the IVP in spherical symmetry.We follow closely the prescription given in [38].
To this end, we consider a spherically-symmetric background, with the following ansatz in polar coordinates (t, r, θ, φ) and same symmetries for the scalar field ϕ = ϕ(t, r).We introduce new variables With this choice of variables and the ansatz of Eq. ( 5), the system (2)-( 4) is reduced to three time-evolution equations for ϕ, P and Q, and two radial constraints for A and B. This system of equations might not always admit a well-posed IVP.We will say that the system is wellposed if there exists a unique solution that depends continuously on its initial data.This will be the case if the system is strongly hyperbolic, i.e. its principal part (the highest derivative terms) is diagonalizable with real eigenvalues [62,63].
In numerical considerations, the characteristic speeds constitute a useful diagnostic tool as they convey information regarding the character of the evolution equations.In spherical symmetry, they are given by [38] where C, and D depend on the functions A, B, ϕ, P , and Q, and their derivation can be found in the Appendix.
The system is hyperbolic if D > 0, otherwise, the system is elliptic (D < 0) or parabolic (D = 0).Initial data.-As for initial data (ID), we are free to specify the values of ϕ and P at t = 0, whereas by definition we must have Q(0, r) = ∂ r ϕ(0, r).We use two types of initial data, type-I ID is a static Gaussian pulse while type-II ID is an approximately in-going pulse where a 0 , r 0 and w 0 are constants.
Numerical implementation.-Weemploy a fully constrained evolution scheme.The equations are discretized over the domain [0, T ] × [0, R] for some choice of R, and T .For a given resolution N , the radial step ∆r is given by ∆r = R/N , and the time step is defined by setting the Courant parameter to λ = 0.25 i.e., ∆t = λ∆r.To impose regularity at the centre we take a staggered grid and perform the following expansion and solve for A 0 , A 2 , A 4 , B 2 , B 4 by using the initial data.Moreover, the boundary conditions at r = 0 require that ∂ r A = 0 and ∂ r B = 0.One can see from the equations that by choosing B 0 = 0 this is automatically satisfied.
Without such treatment at the origin, it was only possible to evolve the system for a small subset of coupling constants.At the outer boundary of the domain, we impose approximately outgoing boundary conditions for ϕ, Q, P , while the metric functions are completely specified by solving the constraint.Given the initial data, we solve the constraint equations for the metric functions on the discretized domain [0, R] using a fourth-order Runge-Kutta (RK4) scheme.This scheme requires knowledge of the value of ϕ, Q, and P at intermediate virtual points and we obtain such information by employing a fifthorder Lagrangian interpolator.Once we have the solution, we set A(t, r) → A(t, r) − A(t, R) by utilizing the remaining gauge freedom.This guarantees that at the outer boundary of the domain, the time function t measures the proper time of a static observer.After obtaining the metric functions, we integrate the variables ϕ, Q, and P in time using the method of lines with an RK4 scheme and a sixth-order Kreiss-Olliger (KO) dissipation term.We use a fourth-order finite differences operator satisfying summation by parts to discretize radial derivatives [64].At the initial step we compute the initial Misner-Sharp mass [65] as M = R 1 − e −2B(0,R) /2, which we use for the rescaling of the output.We keep track of the formation of an apparent horizon by checking the largest radius for which exp(A − B) falls below a chosen tolerance, indicating the formation of a BH.
Results.-We first study how the inclusion of the Ricci term improves the hyperbolic nature of the evolution system for a specific value of α/M 2 = 0.25.In Fig. 1 we plot how the maximum of c − and the minimum of c + vary with time for different values of β, type-II ID with r 0 = 25, w 0 = 6 and for two different amplitudes a 0 = 0.01, 0.016.A vanishing or rather small value of β leads to the formation of a naked elliptic region (NER), i.e. a region of spacetime where the characteristic speeds change sign signalling that the character of the equations has changed from hyperbolic to elliptic without being shielded by an apparent horizon [66].On the other hand, for a sufficiently large coupling to the Ricci scalar, the plot suggests that this coupling heals the problem We show the maximum of c − (solid lines) and the minimum of c + (dashed lines) in space for different values of β.In the top panel, for β = 0, 0.35 the characteristic speeds approach zero and will eventually cross it and change sign.Therefore, the equations change character and we cannot follow the evolution further, but for β = 0.525 the characteristic speeds do not cross zero and the end state of the evolution is flat geometry.In the bottom panel, we observe the formation of an apparent horizon for β = 0.5.The curves indeed asymptote to zero, unlike the β = 0, 0.2 cases.
and allows for the evolution to continue for later times.
Motivated by these results, we explore the parameter space more thoroughly.We start with fixed α/M 2 = 0.25 and type-II ID with r 0 = 25, w 0 = 6, and varied the amplitude a 0 in the range [0.2, 2] × 10 −2 and the Ricci coupling β in the range [0, 1].For each of the cases considered, we establish whether the outcome of the evolution is flat space, NER, or the development of an apparent horizon.Outcomes are summarized in Fig. 2. For β = 0 and low enough initial amplitudes, the theory is predictive, and the outcome is flat spacetime.For certain larger amplitudes an apparent horizon forms.However, we encounter a NER both in the transition between these two outcomes and when we increase the amplitude further.Remarkably, increasing the value of β eventually removes the NER in all cases.
On a few occasions, labelled by an orange star in Fig. 2, it was not possible to conclude whether the outcome is with r 0 = 25, w 0 = 6, and a varying amplitude a 0 .The end state of the evolution is indicated by green squares for flat spacetimes, red crosses for NERs, and black dots for BHs.Orange stars denote cases for which it is difficult to conclude if the end state is a BH or a NER (in polar coordinates).For large enough β the Ricci coupling "cures" the loss of hyperbolicity for both flat space and BH end states.a formation of an apparent horizon or a NER due to the steep gradients in the metric sector.This is due to the fact that polar coordinates are not horizon penetrating and it does not affect our previous conclusions regarding the effect of the Ricci coupling.We have performed additional simulations with different values of α/M 2 and for the two different types of ID that confirm the positive effect this coupling has on hyperbolicity in the case of spherical collapse.A summary of the other simulations performed can be found in table I in the Appendix .
Since for large enough β we can track the evolution, we delve a bit deeper into the properties of the end state.We present three representative cases in Fig. 3, where we show an early and a late snapshot of the evolution of the scalar field ϕ and the metric combination exp (A − B), which vanishes on the apparent horizon, for type-II ID with fixed r 0 = 25, and w 0 = 6 for three different cases.(The fact that exp (A − B) vanishes also at smaller radii, once an apparent horizon forms, is due to the use of coordinates that are not horizon-penetrating.)For each of them, the amplitude a 0 is chosen such that it would produce an apparent horizon in GR.From the analysis of static solutions in Ref. [28], we know that, for large enough values of β, BHs below some threshold mass and down to some minimum mass will exhibit hair.BHs over the threshold mass will have no hair.Below that minimum mass, the Schwarzschild metric is unstable and no hairy BHs appear to exist either.The parameters of the three panels of Fig. 3, top to bottom, have been chosen FIG.3: Scalar field and e A−B , which vanishes at the apparent horizon, as a function of radius for two different time instances for type-II ID with fixed r 0 = 25, and w 0 = 6.ϕ static is the static solution.The black dotted line in the first two panels is the location of the apparent horizon.Top panel: is for α/M 2 = 0.75, β = 1.5 with a 0 = 0.015.Snapshots with {t 1 , t 2 }/M = {22, 88} show the formation of an apparent horizon where the quantity e A−B → 0, with the late time behavior maintaining a scalar profile with a 1/r fall off.The blue dashed curve shows the static solution for the corresponding α/M 2 at each instance (as the mass might decrease during the evolution).Middle panel: is for α/M 2 = 0.25, β = 1 with a 0 = 0.02 and {t 1 , t 2 }/M = {11, 28}.In this case, an apparent horizon forms but the late-time behavior does not support a scalar profile.Bottom panel: is for α/M 2 = 1.25, β = 2 with a 0 = 0.01 and {t 1 , t 2 }/M = {39, 97}.Here, the end state is flat geometry for which the scalar field dissipates to infinity.such that the mass associated with the ID corresponds to each of the three cases respectively.As can be seen from the plots, the end state of evolution is in agreement with the analysis of the static solutions.When applicable, we provide a comparison with a static and spherically symmetric profile, obtained by integrating numerically the equations ( 2)-(4), as described in [28].Note that the lack of perfect overlap between the scalar profile and the static configuration is in part due to the limitations of using polar coordinates.Remarkably, in the case where the static limit predicts that no (stable) BH can exist (bottom panel), we find no impediment in the time evolution and the scalar field dissipates to infinity leaving a flat background.Conclusions.-Wehave investigated the effect of nonminimally coupling the scalar field quadratically to the Ricci scalar on the well-posedness of the IVP in sGB gravity, in the case of spherical collapse.Our results show that this additional coupling, for large enough values of the corresponding coupling constant, can mitigate against the formation of a NER, which signals a loss of hyperbolicity and plagued earlier numerical simulations.This demonstrates that including specific additional interactions -other than those that are essential for having interesting phenomenology, such as BH hair -can be a successful strategy for tackling well-posedness problems in EFTs of gravity with nonminimally coupled scalars.
There are three important limitations to our results, which also highlight interesting future directions.The first one is spherical symmetry.It is likely that the coupling to the Ricci scalar might not be sufficient to cure ill-posedness in a less symmetric setup and that a broader range of couplings would need to be explored.3 + 1 simulations would also allow us to explore numerically the non-radial stability of scalarized black holes [67,68].The second limitation is that we only considered the collapse of a scalar cloud.Work is already underway to generalize our results to the case of stellar collapse.The third limitation is that our numerical implementation involved coordinates that are not horizon penetrating, and hence are not ideal for probing the properties of the end state when the latter is a BH.Interestingly, our simulations did allow us to confirm the expectations arising from combining static results and linear perturbations (e.g.[11,12,21,28,59]: that spherical collapse will lead to BHs with scalar hair when their mass is below a mass threshold and above a minimum mass bound and that above the mass threshold collapse leads to BHs without hair.Remarkably, in simulations that would form a BH below the minimum mass bound, where stable BHs are not know to exist, the scalar cloud smoothly dissipated, leaving behind flat space.We are currently exploring all of these cases in more detail in a numerical implementation that uses horizonpenetrating coordinates, which should allow us to trace the evolution past the formation of an apparent horizon and probe the end state in more detail.

APPENDIX
Characteristic speeds.-Given a system of first-order partial differential equations V I (x, u, ∂u) = 0,1 with the index I counting the number of equations, x µ being the spacetime coordinates, and u J denoting the field content.The principal symbol is defined as follows [38,44,63,69] for a given covector ξ µ .The covector ξ µ is called characteristic if it satisfies the characteristic equation The system of equations is well-posed if all solutions of the characteristic equation are real.In spherical symmetry, the principal symbol matrix can be written as [38] where where E Q = 0, E P = 0 are the evolution equations, and C A = 0, C B = 0 are the constraint equations.The characteristic speed in spherical symmetry is defined as c := − ξt ξr , if ξ µ satisfies the characteristic equation then the corresponding values of the characteristic speeds are given by where, In both panels, we show the convergence of the Misner-Sharp mass where we plot the absolute difference between the low and medium resolutions, and between the medium and high resolutions for type I-ID with r 0 = 25, w 0 = 6.The latter is rescaled by a factor of c 4 = 16 for both plots showing fourth-order convergence.In the top panel, we have a 0 = 0.3, and α/M 2 = 1.25, β = 0.5 with the end state being flat spacetime.In the bottom panel, we have a 0 = 0.4, and α/M 2 = 0.025, β = 0.9, which forms an apparent horizon indicated by the black dotted line.
present two cases, one for which the end state is a flat background, and another where we observe the formation of an apparent horizon.We use three different resolutions {∆r Low , ∆r Mid , ∆r High } = {0.0625,0.03125, 0.015625}.In Fig. 4 we display the convergence behavior of the Misner-Sharp mass.In both panels, we show fourth-order convergence up until the formation of the apparent horizon (in the second panel) which is expected.Exploring the parameter space.-Intable I, we summarize some additional simulations we performed.For each case, we specify what type of ID we used, the amplitude of the pulse a 0 (we fix for all the cases r 0 = 25 and w 0 = 6), the value of the coupling constants α/M 2 and β, the Misner-Sharp mass M at t = 0 that we use to rescale lengths, and finally the outcome of the simulation.Among the outcomes, we distinguish the following cases.We dub flat all the cases in which the scalar field disperses to infinity; NER when the discriminant D becomes negative, signalling the appearance of a NER; BH and sBH when we can trace the formation of an apparent horizon from the condition exp(A − B) → 0. In the former, the scalar field approaches zero in the late time stages of the evolution, while in the latter, it is consistent with stationary scalarized BHs.
Finally, we note that in some cases, which encompass large values of α/M 2 and β, we encountered difficulty in imposing the boundary conditions at the centre.A solution to this problem was to increase the order of the Taylor expansion at the centre to select accurate boundary conditions for A and B as shown in equations ( 11)- (13).We also noticed that increasing the strength of the KO dissipation and decreasing the order of the stencils of finite differences from fourth to second order mitigated this effect.

5 FIG. 1 :
FIG.1:Both plots are for α/M 2 = 0.25 and type-II ID with r 0 = 25, w 0 = 6.The plot on the top is for a 0 = 0.01, and the one on the bottom is for a 0 = 0.016.We show the maximum of c − (solid lines) and the minimum of c + (dashed lines) in space for different values of β.In the top panel, for β = 0, 0.35 the characteristic speeds approach zero and will eventually cross it and change sign.Therefore, the equations change character and we cannot follow the evolution further, but for β = 0.525 the characteristic speeds do not cross zero and the end state of the evolution is flat geometry.In the bottom panel, we observe the formation of an apparent horizon for β = 0.5.The curves indeed asymptote to zero, unlike the β = 0, 0.2 cases.

FIG. 2 :
FIG.2: Scatter plot for α/M 2 = 0.25 and type-II ID with r 0 = 25, w 0 = 6, and a varying amplitude a 0 .The end state of the evolution is indicated by green squares for flat spacetimes, red crosses for NERs, and black dots for BHs.Orange stars denote cases for which it is difficult to conclude if the end state is a BH or a NER (in polar coordinates).For large enough β the Ricci coupling "cures" the loss of hyperbolicity for both flat space and BH end states.
FIG. 4:In both panels, we show the convergence of the Misner-Sharp mass where we plot the absolute difference between the low and medium resolutions, and between the medium and high resolutions for type I-ID with r 0 = 25, w 0 = 6.The latter is rescaled by a factor of c 4 = 16 for both plots showing fourth-order convergence.In the top panel, we have a 0 = 0.3, and α/M 2 = 1.25, β = 0.5 with the end state being flat spacetime.In the bottom panel, we have a 0 = 0.4, and α/M 2 = 0.025, β = 0.9, which forms an apparent horizon indicated by the black dotted line.

TABLE I :
This table is for r 0 = 25, and w 0 = 6.M is the initial Misner-Sharp mass, and the outcome indicates the end state of the evolution.It is evident that the addition of the Ricci coupling yields a well-posed IVP.