{'fraction_non_alphanumeric': 0.03656649271189491, 'fraction_numerical': 0.015224041749145222, 'mean_word_length': 4.6429731925264015, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of spontaneous black hole scalarization and mergers in Einstein-scalar-Gauss-Bonnet gravity

We study the dynamics of black holes in Einstein-scalar-Gauss-Bonnet theories that exhibit spontaneous black hole scalarization using recently introduced methods for solving the full, non-perturbative equations of motion. For one sign of the coupling parameter, non-spinning vacuum black holes are unstable to developing scalar hair, while for the other, instability only sets in for black holes with sufficiently large spin. We study scalarization in both cases, demonstrating that there is a range of parameter space where the theory maintains hyperbolic evolution and for which the instability saturates in a scalarized black hole that is stable without symmetry assumptions. However, this parameter space range is significantly smaller than the range for which stationary scalarized black hole solutions exist. We show how different choices for the subleading behavior of the Gauss-Bonnet coupling affect the dynamics of the instability and the final state, or lack thereof. Finally, we present mergers of binary black holes and demonstrate the imprint of the scalar hair in the gravitational radiation.

We study the dynamics of black holes in Einstein-scalar-Gauss-Bonnet theories that exhibit spontaneous black hole scalarization using recently introduced methods for solving the full, nonperturbative equations of motion. For one sign of the coupling parameter, non-spinning vacuum black holes are unstable to developing scalar hair, while for the other, instability only sets in for black holes with sufficiently large spin. We study scalarization in both cases, demonstrating that there is a range of parameter space where the theory maintains hyperbolic evolution and for which the instability saturates in a scalarized black hole that is stable without symmetry assumptions. However, this parameter space range is significantly smaller than the range for which stationary scalarized black hole solutions exist. We show how different choices for the subleading behavior of the Gauss-Bonnet coupling affect the dynamics of the instability and the final state, or lack thereof. Finally, we present mergers of binary black holes and demonstrate the imprint of the scalar hair in the gravitational radiation.

I. INTRODUCTION
In recent years, the observations of black holes (BHs) through electromagnetic and gravitational waves have furnished new opportunities to test our understanding of gravity (see, e.g. [1][2][3][4][5][6][7][8][9][10][11]). However, in order to perform model selection tests of General Relativity (GR) with these observations, one needs accurate predictions for modified gravity theories in the strong field and dynamical regime, which in these cases, is an outstanding theoretical problem.
A interesting class of theories to test against GR is Einstein-scalar-Gauss-Bonnet (ESGB) gravity. Variants of ESBG give rise to BH solutions with scalar hair, hence can differ qualitatively from GR in the strong field regime (e.g. in BH mergers) while still passing weak field tests. Here we focus on versions of ESGB where GR solutions with vanishing scalar field remain solutions of the modified theory, but which in some circumstances are unstable to perturbations in the scalar field (this contrasts with linearly-coupled ESGB, studied extensively [12][13][14][15][16][17][18][19][20][21], where stationary BHs always form a scalar cloud.) Recently, stationary scalarized BH solutions to the full equations of motion were constructed in this class of theories [22][23][24][25][26][27][28]. These studies show that for particular ranges of mass and spin, set by the magnitude and sign of the Gauss-Bonnet (GB) coupling, the scalarized BHs can differ significantly from their GR counterparts (e.g. the BH spacetime can have 20% of its mass in the scalar cloud), likely impacting gravitational wave observations. These scalarized BH solutions are plausibly the end state of the linear scalarization instability of vacuum BHs of sufficiently small mass (in comparison to the GB coupling length scale) [22,23] or sufficiently high spin (re- * weast@perimeterinstitute.ca † lloydripley@gmail.com ferred to as spin-induced scalarization) [26], which would provide a formation channel. However, the scalarization process is poorly understood, particularly in the spinning BH case, where not even the stability of the scalarized solutions is known. The dynamics in these theories have only been studied by treating the scalar field as a test on a GR background, both for isolated [26,29,30] and binary BHs [31] (with the expection of a nonlinear study of Ref. [32] which considered a different ESGB variant than considered here in spherical symmetry). A major challenge in studying the nonlinear dynamics for these theories has been in finding a well-posed scheme for the ESGB equations of motion (EOM). Here, we build on the methods of Ref. [33], where we demonstrated the feasibility of finding full binary BH solution solutions in linearly-coupled ESGB gravity using the modified generalized harmonic (MGH) formulation [34,35]. Using these methods, we study the nonlinear development and saturation of the scalarization instability for two variants of ESGB gravity, working with initially nearly vacuum (binary) BH solutions with a small scalar field perturbation. We find a range of parameters where this initial data leads to the formation of a stationary scalarized BH. However, we also find that in a significant portion of the parameter space, the predictability of the theory (i.e. the hyperbolicity of the equations of motion) breaks down during scalarization, even when stationary scalarized BH solutions exist. By studying binary BH mergers, we show that even when restricted to this parameter space, BH scalarization can have a significant impact on the resulting gravitational waves.

II. SPONTANEOUS BH SCALARIZATION IN ESGB GRAVITY
We first briefly review the theories we consider, and the heuristic arguments behind why BHs may be dynamically unstable to scalar field perturbations in these theories.

arXiv:2105.08571v2 [gr-qc] 7 Sep 2021
The action for ESGB gravity is where G = R 2 −4R µν R µν +R αµβν R αµβν is the GB scalar, and we use geometric units with G = c = 1 here and throughout. ESGB gravity appears in the low-energy effective actions for certain string theories [36,37] and more generally captures the leading order scalar-tensor interactions expected in an effective gradient expansion of the Einstein equations [34,38]. We will consider two different classes of scalar GB coupling that allow for spontaneous BH scalarization: where λ, σ, and λ e are constants. Note that these agree to leading order in φ 2 when λ = λ e . We choose these because scalarized BH solutions with these couplings have been constructed and studied [22][23][24][25][26][27][28]. The first captures the leading and first subleading term for a coupling invariant under φ → −φ, while the second is one particular higher-order completion with this property.
To review the idea behind spontaneous scalarization, we only need to consider the scalar field equation of motion: where, expanding around φ = 0, β (φ) = β (0)φ+O(φ 2 ). Provided β G > 0, this term will act like a tachyonic mass, and for small enough BH masses, there will be a linear instability which could potentially give rise to a stable scalarized BH solution. There are two possibilities for G and β to have the same sign. For non-spinning and slowly spinning BHs, G > 0 everywhere exterior to the BH, so in order to see spontaneous scalarization one needs β > 0 [22][23][24][25]. For rapidly enough rotating BHs, G is no longer positive definite, which allows for spin induced spontaneous scalarization if β < 0 [26][27][28].

III. METHODOLOGY
We numerically evolve the full ESBG EOM using the MGH formulation [34,35] as described in Ref. [33]. We use similar choices for the gauge, numerical parameters, etc. as in Ref. [33], except that we find the scalarized BHs evolved here also benefit from the addition of long wavelength constraint damping obtained by setting ρ = −0.5 in Eq. (2) of Ref. [33].
For initial data, we start from single or binary vacuum BH solutions (the latter constructed as in Ref. [39]) with a small Gaussian scalar perturbation centered on the BH(s). For most cases presented here, we use an initial amplitude of φ 0 = 0.01, though we have verified smaller amplitudes give the same results, and that the error induced by not solving the constraint equations including the perturbation is negligible. See the Supplementary Material (which cites the references [22-25, 32, 33, 40]) for details on resolution, convergence, and the exact form of the initial perturbation.
We use many of the same diagnostics as in Ref. [33], which we briefly review here. We measure the gravitational radiation by extracting the Newman-Penrose scalar Ψ 4 , and use this to calculate an associated gravitational wave luminosity P GW . We also measure the flux of energy in the scalar field P SF .
During the evolution, we track any apparent horizons present at a given time, and measure their areas and associated angular momentum J BH . From this we compute a BH mass M BH via the Christodoulou formula. We will refer to the mass that lies outside the BH horizon(s)which, to a good approximation, can be attributed to the scalar cloud-as M φ ≡ M − M BH , where M is the global Arnowitt-Deser-Misner (ADM) mass of the spacetime, and similarly define an angular momentum J φ ≡ J − J BH . Finally, we compute the scalar charge from the asymptotic behavior of the scalar field at large r: φ = Q SF /r + O(1/r 2 ), where we have fixed that φ → 0 at spatial infinity with our initial conditions.
In addition to evolutions with the MGH formulation, we also present evolutions of BHs in spherical symmetry, using the formalism and code described in Ref. [32]. We do this in order to determine exactly where the hyperbolicity of ESGB breaks down in spherical symmetry for various choices of coupling, which we can compare to the MGH evolutions without symmetry assumptions. Spherically symmetric spacetimes are not only computationally less expensive (and thus we are able to systematically scan the parameter space), but there is less gauge ambiguity in determining when the EOM are hyperbolic (although see Ref. [41] for a recently introduced formalism for generic backgrounds).

IV. BLACK HOLE SCALARIZATION AND SATURATION
We begin by considering the scalarization of isolated non-spinning and spinning BHs using our MGH code [33]. Our main result is that we find a range of parameters for both signs of β (0) where the scalarization instability stability saturates in the formation of a stable BH with scalar hair containing up to a few percent of the total mass.
Considering first β (0) > 0, we show the dynamics of BH scalarization for several cases in Fig. 1. Following an exponential growth phase, where the BH develops scalar charge at the expense of losing mass, we find that the instability eventually saturates and settles to a nearly stationary scalarized BH solution. (We note that the area of BHs can decrease in theories like ESGB that violate the Scalarization of BHs with positive GB coupling. Top: The scalar charge measured at a large distance as a function of look-back time. Bottom: The amount of mass (solid lines) and, for the cases with spin, angular momentum (dashed lines) going from the BH into the scalar cloud and radiation. The legend in the top panel applies to the bottom as well.
Null Convergence Condition, see e.g. [18].) Increasing λ e results in higher instability rates and more massive scalar clouds at saturation. At fixed coupling, considering nonzero BH spin decreases the instability rate and cloud mass. In addition to results with the exponential coupling [Eq. (A1b)], in Fig. 1 we also show one case with the coupling given by Eq. (A1a) and λ = −σ/10 = 1.4M 2 .
(We note that with σ = 0, spherical scalarized BHs are radially unstable [42].) In this case, a similar amount of mass goes into the scalar cloud compared to the exponential coupling with λ e = 0.78M 2 , but the instability happens much faster and initially overshoots, e.g., the final scalar charge (see top panel of Fig. 1).
We were unable to obtain hyperbolic evolutions through saturation for non-spinning BHs with positive λ e or λ = −σ/10 much higher than the above mentioned values. As we discuss in the next section, this is because we are approaching the regime where the asymptotic hyperbolicity of the theory breaks down.
We also consider spin-induced BH scalarization, which occurs when β (0) < 0, finding similar results. In Fig. 2, we show several different cases where the initial BH spin ranges from a BH = 0.8 to a BH = 0.95 with different values of the coupling given by Eq. (A1a) (with σ = 0 unless otherwise noted). We again find a range of parameters where the scalarization instability saturates and leads to the formation of a stationary BH solution with up to a few percent of the mass of the BH converted to scalar hair. For higher values of spin, the instability sets in at lower values of the coupling. We find the gravitational waves from the scalarization process to be negligible, but the scalar radiation increases with the instability rate (see top panel of Fig. 2; the time dependence of Q SF from these cases is similar to the β (0) > 0 cases). Though not shown, we also considered an exponential coupling up to λ e = 1.05M 2 for a BH = 0.9 and found similar results. Again, we were not able to obtain numerical evolutions for significantly higher values of the coupling than shown in Fig. 2 for the given BH spin values. That is, the breakdown in the evolution occurs for lower coupling values for higher BH spins.

V. LIMITS ON SCALARIZED BHS FROM HYPERBOLICITY
For the case of spherically symmetric BHs (i.e. using the code in Ref. [32]), we will now explicitly demonstrate that the reason we are unable to follow the dynamical formation of BHs with scalar hair containing more than a few percent of the total mass is because the hyperbolicity of the theory breaks down (and we conjecture that something similar happens beyond spherical symmetry).
For the polynomical coupling, Eq. (A1a), we fix λ/M 2 , and then find the value of σ for which the theory becomes elliptic. As the former is made larger, the latter must be made more and more negative in order to control the magnitude of the coupling β at saturation. We plot the dividing line in parameter space in Fig. 3. We find that for a given λ 0.84 the minimum absolute value of σ that is necessary for hyperbolic evolution through scalarization is The bottom panel of Fig. 3 shows that the maximum amount of energy liberated from the BH and put into the scalar cloud by the instability is always 5%, and the maximum falls of like M φ 8.5 × 10 −2 × λ −3/2 M 4 for λ/M 2 > 1.5. This may be connected to the fact that, as illustrated above, larger values of λ tend to initially overshoot the scalarized solution (as opposed to smaller values of λ which smoothly saturate) and hence may more easily violate hyperbolicity dynamically.
For the case of positive exponential coupling, Eq. (A1b), we find the dividing coupling between hyperbolic and elliptic evolution for spherical scalarization is λ e = 0.834M 2 . This is within ∼ 10% of the value one would obtain by neglecting O φ 6 terms in this coupling and using Eq. (4) with σ = −3λ. By contrast, the analysis of Ref. [42] concluded that scalarized BH solutions with up to λ e ≈ 8.55M 2 were radially stable (and similar results were found for the quadratic-quartic coupling [24,25].) Comparing these results to the MGH evolutions that do not explicitly enforce spherical symmetry, we see that the values of λ e used for the former are within 10% of the maximum value that retains hyperbolicity in spherical symmetry (and within ∼ 15% for λ = −0.1σ). This difference is likely just due to the more limited numerical resolution used for the MGH evolutions, as approaching extremality, the narrowing region between the elliptic region and the horizon becomes more and more difficult to resolve.
Going beyond spherical symmetry to the case of spininduced scalarization, we do not have any definitive results on the breakdown of hyperbolicity (only positive results establishing hyperbolicity for a range of parameters). However, given the above, we can conjecture that the reason we were not able to evolve significantly larger couplings is that, for this case as well, elliptic regions develop outside the BH horizon during scalarization, even in the regime where stationary scalarized BH solutions exist.

VI. HEAD ON COLLISIONS OF SCALARIZED BHS
Given that we find that scalarized BHs with up to a few percent of the total mass in the scalar cloud can form in the regime where the theory is hyperbolic, it is interesting to consider how this will affect a binary BH merger. Here we focus on the case of the (axisymmetric) headon collision of an equal-mass binary, considering both a non-spinning binary (with λ > 0) and a spinning binary (with λ < 0).
In Fig. 4, we plot the total gravitational wave and scalar field luminosity for head-on black hole collisions. As described above, we begin with two nearly-vacuum The gravitational (top) and scalar (bottom) luminosity for the head-on merger of equal mass BHs. For the λ ≥ 0 cases (red, black, and blue curves), the BHs are nonspinning, while for the λ < 0 case (green curve; multiplied by ×10 2 to be visible), each member of the binary has dimensionless spin aBH = 0.8 initially. We see that the gravitational wave luminosity is larger when two scalarized BHs collide, as compared to two vacuum BHs in the theory.
BH solutions, with a small scalar field perturbation outside their horizons. With the couplings we consider, the BHs subsequently spontaneously scalarize, and reach saturation well before they collide. For λ = −0.1σ m 2 (where m is the mass of one of the binary constituents) we find that the gravitational waves from the scalarized BH collisions have noticeably larger amplitude compared to the GR case, and that the scalar luminosity is comparable to the gravitational wave luminosity. Due to the increased radiation with larger λ, the merger also happens slightly faster for the same initial separation/velocity.
We also show a case with λ = −0.1σ = −2.55m 2 where (unlike the above mentioned case) the binary constituents are both spinning with a BH = 0.8. Here, despite the fact that ∼ 0.9% of the mass is a spin-induced scalar cloud prior to merger (compared to ∼ 1.7% and ∼ 2.8% for the above cases with λ/m 2 = 1 and 1.4, respectively), the scalar radiation from merger is several orders of magnitude smaller than the positive λ cases, and the gravitational wave luminosity (not shown) does not noticeably differ from the GR case. Hence, it appears that larger BH spins are necessary to have a strong impact on the binary BH collision for spin-induced scalarization (though for an inspiral, the impact is likely greater).

VII. DISCUSSION AND CONCLUSION
In this work we have studied full, non-perturbative solutions to several ESGB theories that exhibit spontaneous BH scalarization. We studied theories where the GB coupling was either a single-parameter exponential function, or a two parameter quadratic function of φ 2 , and for both choices of sign, leading to either mass or spin-induced scalarization. Though these theories have a multidimensional parameter space, which we have not fully explored, we can infer several general results. We have shown that the end state of the nonlinear evolution of the linear scalarization instability of vacuum nonspinning and spinning BHs [22][23][24][25] results in the formation of a stable, scalarized BH, for a range of coupling parameters and BH masses and spins. However, we find that for a given set of couplings, at large enough curvature scales, the theory can lose hyperbolicity and, in contrast to the linearly coupled theory [18,33,43], this breakdown occurs at much lower values compared to the maximum values where stationary solutions can constructed. In particular, while stationary solutions can be constructed with over 20% of the total mass is attributable to the scalar cloud, here we did not find any cases where this was greater than 5%.
Thus, in the quadratic (plus higher order) coupled ESGB theories considered here, as for the linearly coupled case, there is a minimum mass for a stable BH to form. In the former case, vacuum BHs below this limit are unstable, while in the latter case they are explicitly non-stationary. But in either case, their subsequent evolution will generically result in a breakdown of the predictability of the theory. Nevertheless, the valid range is still interesting from the point of view of potentially impacting the gravitational wave signal of a binary BH merger, as we have demonstrated for some example headon binary collisions.
We note that, while here we studied the scalarization instability starting from stationary BHs, instead of following the dynamical formation of an unstable BH from collapsing matter (see Ref. [30] for such a calculation in the spherically symmetric, test-field limit), we expect our results to apply to such cases as well. This is because BHs generically form with some matter distribution falling within its Schwarzschild radius at relativistic speeds, while the scalarization instability timescale is generally much longer than the BH light-crossing time.
Finally, for larger GB couplings, the scalarization instability growth rate is faster and there is a tendency to initially overshoot the final stationary solution, resulting in more scalar radiation. However, maintaining hyperbolicity during the evolution also requires larger higher order corrections that reduce the coupling at large field values, and the energy scale at which this breakdown occurs becomes smaller for larger couplings.
This work also shows another example where the methods for numerically evolving Horndeski modified gravity theories of Ref. [33], based on the MGH formula-tion [34,35], work at larger couplings, where the deviations from GR are significant. For future work, it would interesting to study the predicted gravitational wave signal from scalarized binary BH mergers in these theories.

VIII. ACKNOWLEDGMENTS
We thank Thomas Sotiriou for discussions on ESGB gravity and spontaneous scalarization, and Max Corman for comments on our article. W.E. acknowledges support from an NSERC Discovery grant. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Insti-tute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science. This research was enabled in part by support provided by SciNet (www.scinethpc.ca/) and Compute Canada (www.computecanada.ca). Some of the simulations presented in this article were performed on computational resources managed and supported by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's High Performance Computing Center and Visualization Laboratory at Princeton University.
Here we present convergence test results for the modified generalized harmonic code [33] and spherically symmetric code [32] that we use, as well as compare results with different symmetries imposed.
In Fig. 5, we compare the spontaneous scalarization of a spinning BH with a BH = 0.9 and GB coupling λ = −σ/10 = 1.1M 2 for different numerical resolutions and initial scalar perturbations. The medium resolution shown here uses 7 levels of mesh refinement with a refinement ratio of 2 : 1, and has a grid spacing of dx ≈ 0.01M on the finest level, which covers the BH. This is the default resolution for the results presented in this work. For the initial perturbation, we use a Gaussian profile for the scalar that is centered on the BH: For most cases, we use an initial amplitude of φ 0 = 0.01. However, in Fig. 5 we also show a case with φ 0 = 0.005 to demonstrate that this is sufficiently far in the linear regime of the instability, and that the constraint violation due to not including the backreaction of the scalar in the initial data is smaller than the truncation error of the highest resolution studies we considered. We find that the constraints converge to zero at the expected fourth order. Finally, in the top panel of Fig. 5, we also show a case where we do not enforce axisymmetry (in contrast to the  The scalar field energy (as measured from the canonical stress-energy tensor; top panel) and violation of the MGH constraints (bottom panel) for the scalarization of a BH with dimensionless spin aBH = 0.9 and coupling λ = −σ/10 = 1.1M 2 . The initial scalar perturbation has amplitude either φ0 = 0.01 or 0.005 for the axisymmetric cases, or φ0 = 0.015 (and is equatorially offset from the BH) for the fully 3D case. A constant offset has been applied to the time shown on the horizontal axis, so that the curves align when ESF/M = 10 −4 . After a brief transient phase, we find the linear growth and nonlinear saturation to be the same. We show results from three different numerical resolutions, where the medium and high resolutions correspond to, respectively, ×4/3 and ×8/3 as many points in each linear dimension compared to the low resolution. The decrease in the constraints with resolution is consistent with approximately fourth order convergence.
other cases), and where we choose the center of the initial Gaussian perturbation to be offset from the BH by M/4 in the equatorial plane, in order to explicitly break axisymmetry. We find that scalarization and saturation proceeds as in the axisymmetric case.
For the spherical symmetric evolutions, we do exactly solve the constraint equations for ESGB gravity when constructing perturbed initial data. As in Ref. [32] we use initial data that consisted of a small amplitude scalar field bump profile outside of a Schwarzschild BH: where r u > r l > 2M , where M is the initial Misner-Sharp mass of the BH (see [32] for more discussion). We consider r l = 2.4M , r u = 5.2M , and a 0 = 10 −2 , so that the initial scalar field seed only adds δM φ /M 10 3 to the total mass of the spacetime. In spherical symmetry, we compute the global mass by extracting the Misner-Sharp mass at spatial infinity [32] (we note that in spherical symmetry the Misner-Sharp mass is equal to the ADM mass at spatial infinity [40]). Note that we need a 0 = 0, as the Schwarzschild solution is a stationary solution to the theories we consider. In Fig. 6, we demonstrate convergence of the spherically symmetric code. FIG. 6. We plot the value of the scalar field at the BH horizon, and the norm of the rr component of the Einstein equations in spherical symmetry, for the case of exponential coupling Eq. (A1b), with λ/M 2 = 0.8. The slow decrease in the value of φ over time is due to the slow increase (due to truncation error) of the BH mass. The low, medium, and high resolutions are with radial grid resolution nx = 2 12 + 1, nx = 2 13 + 1, and nx = 2 14 + 1, respectively. We see that we can both resolve the scalar field dynamics, and that the residuals are converging between second and fourth order, which is consistent with the second order method we use to solve for the metric variable and fourth order stencils we use to evolve the scalar field equations of motion (for more details on the numerical implementation, see Ref. [32]).