Dynamical descalarization in binary black hole mergers

Scalar fields coupled to the Gauss-Bonnet invariant can undergo a tachyonic instability, leading to spontaneous scalarization of black holes. Studies of this effect have so far been restricted to single black hole spacetimes. We present the first results on dynamical scalarization in head-on collisions and quasicircular inspirals of black hole binaries with numerical relativity simulations. We show that black hole binaries can either form a scalarized remnant or dynamically descalarize by shedding off its initial scalar hair. The observational implications of these findings are discussed.

Introduction.− Despite the elegance of Einstein's theory, it presents several shortcomings: explaining the late-time acceleration of the Universe and providing a consistent theory of quantum gravity or the presence of spacetime singularities [e.g., in black holes (BHs) ]. Candidate theories (of quantum gravity) that remedy these shortcomings typically predict the coupling to additional fields or higher curvature corrections [1]. Binary BHs, their gravitational wave (GW) emission, and the first GW detections by the LIGO-Virgo Collaboration [2,3] offer unique insights into the nonlinear regime of gravity that unfolds during the BHs' inspiral and merger and enable new precision tests of gravity [4,5]. So far, these tests have been parametrized null tests against General Relativity (GR) [6,7] or used a mapping between these parameters and those of specific theories [8][9][10]. To do the latter, however, requires GW predictions in specific theories.
One of the most compelling beyond-GR theories, scalar Gauss-Bonnet (sGB) gravity introduces a dynamical scalar field coupled to the Gauss-Bonnet invariant. sGB gravity emerges in the low-energy limit of quantum gravity paradigms such as string theory [11], through a dimensional reduction of Lovelock gravity [12] and is the simplest model that contains higher curvature operators. The most studied class of sGB gravity with a dilatonic or linear coupling to the scalar field gives rise to hairy BHs [13][14][15][16][17][18][19]. This theory, however, has been strongly constrained with GW observations from binary BHs [9].
We turn our attention to another interesting class of sGB gravity that is both unconstrained by GW observations and gives rise to (spontaneously) scalarized BHs [20,21]. Spontaneous scalarization is a familiar concept in beyond-GR theories; e.g., it is well established for neutron stars in scalar-tensor theories [22,23]. In such theories, the neutron star matter itself can induce a tachyonic instability that spontaneously scalarizes the star [24]. When placed in a binary system, initially unscalarized neutron stars can scalarize dynamically near their merger or a scalarized neutron star can induce a scalar field in their unscalarized companion [25][26][27][28]. In sGB gravity, it is the spacetime curvature itself that in-duces scalarization of BHs [20,21], although this has only been shown for isolated BHs so far. In this Letter we investigate, for the first time, dynamical scalarization in binary BHs. We concentrate on head-on collisions of BHs, but also present the first binary BH inspiral study. Before doing so, it is convenient to first review the basics of sGB gravity and spontaneous BH scalarization.
The latter class still admits all vacuum (BH) solutions of GR together with Φ = Φ 0 = const. In fact, if f (Φ 0 )G < 0 these solutions are unique due to a no-hair theorem [21]. A linear stability study of these Φ 0 = const. solutions around a Schwarzschild BH reveals that this condition is a requirement for the absence of a tachyonic instability (m 2 eff > 0) for the scalar field perturbations [21]. If the effective mass m 2 eff < 0, a tachyonic instability is triggered and the the sGB scalar field is excited and spontaneously scalarizes the BH. This linear instability [34] is quenched at the nonlinear level, resulting in a scalarized BH as end-state [35]. The simplest theory that admits scalarized BHs is described by the quadratic coupling f (Φ) =β 2 Φ 2 , whereβ 2 = const. [36]. The relevant parameter in this theory is the dimensionless constant β 2 = (α GB /m 2 )β 2 , where m is the characteristic mass of the system. The onset of scalarization is fully determined by the scalar's linear dynamics on a given GR background. For a Schwarzschild BH of mass m, for which G 0 everywhere, scalarization first occurs for a spherically symmetric scalar field if β 2 = β c ∼ 1.45123, a result in agreement with nonlinear calculations [20,21]. For values below β c the scalar perturbation decays monotonically at late times (we call them "subcritical"), precisely at β c the scalar field forms a bound state around the BH ("critical"), and above it the scalar field growths exponentially with time ("supercritical"). This result was recently generalized to Kerr BHs, where spin-induced scalarization can take place for β 2 < 0, for dimensionless spin parameters χ 0.5 [37][38][39][40]. Nonlinear rotating scalarized BH solutions in sGB gravity were found for both positive [41,42] and negative values of β 2 [43,44]. So far studies of scalarization in sGB gravity focused on single BHs. We advance these studies to BH binaries, and expand upon [45], focusing on the quadratic theory f (Φ) =β 2 Φ 2 , as discussed next.
Numerical methods and simulations.− We investigate BH scalarization in the decoupling limit, i.e., we numerically evolve the scalar field on a time-dependent background in vacuum GR that represents binary BH spacetimes. Unless stated otherwise, we follow the approach of [45] and refer to it for details. We foliate the spacetime into spatial hypersurfaces with 3-metric γ ij and extrinsic cur- being the Lie derivative along the shift vector β i , and α is the lapse function. We write Einstein's equations as a Cauchy problem and adopt the Baumgarte-Shapiro-Shibata-Nakamura formulation [46,47] of the time evolution equations complemented with the moving-puncture gauge conditions [48,49]. We prepare Brill-Lindquist initial data [50,51] for head-on collisions or Bowen-York initial data [52,53] for a quasicircular BH binary.
To evolve the scalar field, we introduce its momentum K Φ = −α −1 d t Φ, and write its field equation (2) as where D i is the covariant derivative associated with γ ij , K = γ ij K ij , f = 2β 2 Φ, and G is the Gauss-Bonnet invariant of the background spacetime. We set the system's total mass to unit, i.e., M = m 1 + m 2 = 1, where m = m 1,2 is the component's mass. The scalar field is initialized either as a spherically symmetric Gaussian shell (G) located at r 0 = 12M and with width σ = 1M as in [45] or as a bound state (B) around each binary component, Here, = m + 2r, and c 1 = 3.68375, c 2 = 4.972416, c 3 = 4.972416 · 10 2 are fitting constants to reproduce the numerical results in [21].
We perform our numerical simulations with Canuda [45,[54][55][56], coupled to the open-source Einstein Toolkit [57,58]. We extended the implementation of [45,55,56] to general coupling functions f , including the quadratic coupling. We employ the method of lines with fourth-order finite difference stencils to realize spatial derivatives and a fourth-order Runge-Kutta time integrator. We use box-in-box mesh refinement provided by Carpet [59]. The numerical grid contains seven refinement levels, with the outer boundary located at 256M and a grid spacing of dx = 1.0M on the outer mesh. To assess the numerical accuracy of our simulations we evolved case (b) in Fig. 1 with additional resolutions dx = 0.9M and dx = 0.8M . We find second-order convergence and a relative discretization error of ∆Φ 00 /Φ 00 0.5%, where Φ 00 is the = m = 0 multipole of the scalar field. We present the corresponding convergence plot for the scalar monopole and for the gravitational wave = 2, m = 0 mode in Fig. 5 of the Supplemental Material.
Results.− We performed a large set of BH head-on collisions with varying mass ratio q = m 1 /m 2 1, total mass M = 1 and initial separation d = 25M , considering both initial data in Eq. (4). The BHs merge at t M ∼ 179.5M , as estimated from the peak of the = 2, m = 0 multipole of the gravitational waveform. To guide our choices of β 2 , we recall that the critical coupling for the fundamental mode is β 2,c = β c (m/M ) 2 with β c ∼ 1.45123, and m denotes either the individual BHs' mass m 1,2 or the total mass M . For example, for an equal-mass binary with m 1 = m 2 = M/2, the critical coupling for the individual holes is β 2,c = β c /4 = 0.36275 and that of the final hole is approximately β f 2,c = β c where we neglected the small mass loss in the form of GWs during the collision [60,61].
Here we present a selection of our results, illustrated in Fig. 1, to highlight our most important findings. An expanded discussion will be presented in a companion paper [62]. We vary the initial state by setting the coupling parameter β 2 such that (a) none of the BHs are initially scalarized, (b) the smaller-mass BH initially carries a bound-state scalar field, both BHs carry initially s sss ss wheres and s stand for initial or final states that are either nonscalarized or scalarized respectively. Each diagram is labeled by the initial data (Gaussian shell "G" or bound state "B"), the mass ratio q = m1/m2 (1 or 1/2) and the coupling parameter β2. In Fig. 2 we show the = m = 0 scalar field multipole extracted on a sphere of fixed radius r ex = 50M , as a function of time, and we present snapshots of the scalar's profile in the Supplemental Material. In case (a), the scalar perturbation is not supported at all (since m eff = 0) and, indeed, after a brief interaction at early times it decays already before the BHs collide. In cases (b) and (c) we find a constant scalar field before the BHs collide, that is consistent with a bound state around the individual (q = 1) or smaller-mass BH (q = 1/2). After the merger the scalar field decays since the curvature (and thus m eff ) decreases and the system no longer supports a bound state -the final BH dynamically descalarizes. In case (d), the scalar field grows exponentially before the merger because it is supercritical for the individual BHs and settles to a constant in time that is consistent with a bound state around the final BH.
In Fig. 3 we show two-dimensional snapshots of the scalar field and spacetime curvature for case (b) which illustrates the dynamical descalarization phenomenon [63]. The color map is shared among all panels and shows the amplitude of log 10 |Φ|, while the curves are isocurvature levels of G M 4 = {1, 10 −1 , 10 −2 , 10 −3 }. Initially, at t = 1M , both BHs (whose locations are revealed by the isocurvature levels) are surrounded by nontrivial scalar field profiles given by Eq. (4). At t = 50M , the smaller BH hosts a bound state scalar that is dragged along the hole's motion, inducing scalar dipole radiation that would impact the GWs emitted. In contrast, the scalar field around the larger BH disperses because its curvature is too small to sustain a bound state for a coupling of β 2 = 0.36281. The system thus evolves as a s +s process in the notation of Fig. 1. At t = 160M , the BHs are about to merge, as indicated by the two lobes in the isocurvature contours, the curvature of the combined system decreases and the scalar field starts dissipating. At t = 182M , which is shortly after the collision, the system has descalarized since for the final BH β f 2,c > β 2 . We also simulated the inspiral of an equal-mass, nonspinning BH binary with initial separation of d = 10M , β 2 = 0.36281, and bound state scalar field initial data. This corresponds to an initial configuration in which both BHs are scalarized, and then, after merger, the remnant is not scalarized, which is analogous to case (c) of Fig. 1 in the head-on case. In Fig. 4, we show the gravitational quadrupole waveform (bottom panel), as characterized by the = m = 2 mode of the Newman-Penrose scalar Ψ 4 , together with the scalar field's monopole (top) and quadrupole (middle). The scalar's monopole Φ 00 exhibits the distinctive signature of descalarization: the increase in the field's amplitude during the inspiral of scalarized BHs is followed by a complete dissipation of the scalar field after the merger (t M ∼ 917M ) as the curvature of the remnant BH no longer supports a bound state. In addition, the dynamics of the BH binary sources scalar quadrupole radiation (of the initially spherically symmetric scalar). The field's amplitude grows exponentially during the inspiral and decays after the BHs have merged. The origin of this excitation is not direct scalarization of  the = 2 scalar bound state, but due to the inspiral of two scalarized (or "hairy") BHs. This interpretation is further supported by the observation that the phase of the = m = 2 scalar mode is driven by the binary's orbital frequency. We also observed this for the = m = 4 mode and expect it to happen for all even = m modes. For q = 1, the odd = m modes are suppressed due to symmetry, whereas they would be excited in the general case q = 1. The descalarization during the merger is reminiscent of the decrease in scalar charge observed in the shift-symmetric theory [45], however, with the striking difference that here the remnant BH is a rotating GR solution.
Discussions.− We presented the first numerical relativity simulations of the scalar field dynamics in binary BH spacetimes in quadratic sGB gravity [21]. We found that the interplay between mass ratio q and β 2 can result in different scenarios for the scalar field dynamics. Most notably, it can lead to a dynamical descalarization of the binary, which we observed in both head-on and quasicircular inspiral simulations. Here we focused on β 2 0, but the case β 2 < 0 would be particularly interesting to study in inspiral simulations. More specifically, the spinning remnant of a binary BH merger typically has a dimensionless spin χ ∼ 0.7 [64], sufficient to trigger a spin-induced tachyonic instability of the scalar field [37]. This is currently under study [62]. It would be interesting to frame this effect within the effective field theory (EFT) of [65] or in a post-Newtonian framework [66][67][68], although the latter may not be suitable for the modeling of a nonlinear dynamical scalarization process.
The scalar excitations we have discovered during binary BH coalescence in this class of sGB theories have important implications to GW observations and tests of GR. In particular, the scalar excitations will drain the binary of energy as they propagate away from the system, the monopole scalar piece inducing dipole losses, and the quadrupole piece correcting the quadrupole GW losses of GR, which, based on [69], are expected to only have the same "plus" and "cross" polarizations. This enhanced dissipation of energy and angular momentum, in turn, will force the binary to inspiral faster than in GR, and therefore, leave an imprint in the GWs emitted through corrections to the rate at which the GW frequency increases during the inspiral. This GW phase shift will enable us to project bounds on sGB gravity that are similar in spirit but complementary to the analysis of [45]. In fact, because the merger leaves behind a "bald" Kerr black hole due to dynamical descalarization, the (scalar) energy flux is, in general, larger as compared to shiftsymmetric sGB, where the remnant black hole always retains some of its hair. This suggests that strong observational bounds might be placed on this theory.
Having worked in the decoupling limit, a question nat-urally arises: what would we expect in the fully nonlinear regime of sGB gravity? It is known that nonlinear effects set an upper bound on the scalar field magnitude at the BH horizon [29], so that the domain of existence of scalarized BHs exhibits a very narrow bandlike structure in the phase space spanned by BH mass and coupling β 2 ; see Fig. 2 of [21]. This means that case (d) would only occur for sufficiently small mass ratios such that both the initial binary and its final state remain in band. In general, however, comparable mass BH binaries could undergo ans +s → s process, in which two unscalarized BHs would merge, forming BH within the scalarization band, i.e., a dynamical BH scalarization. The descalarization of the BH remnant would also impact the GW emission during the ringdown. Specifically, the waveforms in Fig. 4 show that the ringdown time scales of scalar and tensorial modes are comparable. This suggests that one should expect to see the imprint of the descalarization onto the quasinormal mode spectra of the Kerr black hole in the nonlinear case. Performing these studies in practice would require a general, wellposed formulation of the time evolution equations outside the EFT approach [45,70], small values of the coupling parameter [71,72], or spherical symmetry [35,[73][74][75].
Finding such a formulation has proven challenging [76][77][78][79], although first results in this direction were presented in [80]. Our work motivates and paves the way for future studies of nonperturbative, beyond-GR effects in BH binaries, with potential implications to tests of GR with GW astronomy. -SUPPLEMENTAL MATERIAL -

CONVERGENCE PLOTS
We assess the discretization error of our simulations by exemplarily running the head-on collision of equal-mass black holes that initially carry a bound-state scalar field with coupling parameter β 2 = 0.36281 at three different resolutions dx c = 1.0M , dx m = 0.9M and dx f = 0.8M . Here M is the system's total mass, which we set to unit. This setup corresponds to case (c) in Fig. 1 of the main text. Focusing on the scalar field monopole (Φ 00 ) and the gravitational quadrupole (Ψ 4,20 ) we compute the differences between the course and medium, and medium and high resolution runs. For Φ 00 , we rescaled the latter difference by the convergence factor Q 2 = 1.12, as shown in the left panel in Fig. 5, indicating second-order convergence. For Ψ 4,20 , we rescaled the latter difference by Q 4 = 1.39, as shown in the right panel, indicating fourth-order convergence. Computing the relative difference ∆Φ 00 /Φ 00 between the coarsest resolution simulation with dx c = 1.0M and the second-order Richardson extrapolation, we find a numerical error of ∆Φ 00 /Φ 00 0.5% as stated in the main text.  Figure 6 presents the scalar field profile along the collision axis x/M at different instances throughout the evolution before, near and after the merger of the BHs. In case (a), the scalar field is below the critical value to form any bound state configurations and, indeed, after a brief interaction at early times it decays already before the BHs collide. In cases (b) and (c), the scalar field forms a bound state that is anchored around the individual (q = 1) or smaller-mass BH (q = 1/2). As the BHs approach each other, the scalar field follows their dynamics and moves along the collision course with only small adjustments to its spatial configuration. After the BHs merge, the critical value β 2,c to form a bound state increases, i.e., the BH can no longer support a scalar bound state. Consequently, the configuration becomes subcritical and the scalar field is depleted, indicating dynamical descalarization of the BH binary. Finally, case (d) is set up such that the final configuration is near critical to form a bound state, always leading to a supercritical setup before merger. Indeed, we observe that the scalar field grows (exponentially), before settling to a constant-in-time radial profile after the merger. This rapid growth is due to the fact that β 2 ∼ 1.45123 is four times larger than the critical scalarization value for the initial BHs.