---
abstract: 'We study scalar fields in a black hole background and show that, when the scalar is suitably coupled to curvature, rapid rotation can induce a tachyonic instability. This instability, which is the hallmark of spontaneous scalarization in the linearized regime, is expected to be quenched by nonlinearities and endow the black hole with scalar hair. Hence, our results demonstrate the existence of a broad class of theories that share the same stationary black hole solutions with general relativity at low spins, but which exhibit black hole hair at sufficiently high spins ($a/M\gtrsim 0.64$). This result has clear implications for tests of general relativity and the nature of black holes with gravitational and electromagnetic observations.'
author:
- Alexandru Dima
- Enrico Barausse
- Nicola Franchini
- 'Thomas P. Sotiriou'
bibliography:
- 'master.bib'
title: 'Spin-induced black hole spontaneous scalarization'
---

***Introduction:*** Direct and indirect detections leave little doubt that black holes (BH) exist in nature . In general relativity (GR) the mass and the spin of an astrophysical BH fully determine its properties. An electric charge is also technically allowed, but is expected to be paltry for astrophysical BHs, see e.g. [@Barausse:2014tra]. Any other quantity, *hair* in jargon, is not necessary according to no-hair theorems [@Israel:1967wq; @Carter:1971zc; @Robinson:1975bv]. Future gravitational wave detectors will finally allow us to confront theorems and observations with unprecedented precision [@Berti:2016lat; @Barausse:2016eii; @Toubiana:2020vtf], improving upon current observations, which are perfectly compatible with hairless BHs [@LIGOScientific:2019fpa; @Isi:2019aib].

It is tempting to interpret an absence of BH hair as a vindication of GR minimally coupled to the Standard Model. However, new fundamental fields can be more elusive. It is illustrative to consider scalar fields: no-hair theorems exist for stationary BHs in scalar-tensor theories [@Hawking:1972qk; @Sotiriou:2011dz], and static, spherically symmetric and slowly rotating BHs in shift-symmetric generalized (Horndeski) scalar-tensor theories [@Hui:2012qt; @Sotiriou:2013qea].[^1] In fact, it turns out that there is a single coupling term in the Horndeski class that gives rise to hair: a linear coupling between the scalar and the Gauss-Bonnet (GB) invariant [@Sotiriou:2013qea; @Sotiriou:2014pfa], given by $$\label{eq:GBinv}
\mathcal{G}=R^{\mu\nu\rho\sigma}R_{\mu\nu\rho\sigma}-4~R^{\mu\nu}R_{\mu\nu}+R^2.$$ Considering that the Horndeski class contains all actions for a massless scalar nonminimally coupled to gravity that yield second order equations upon variation, absence of hair actually seems to be the norm rather than the exception for scalar fields. Indeed, known hairy BH solutions circumvent theorems by evading one or more of their assumptions, see [*e.g.*]{} [@Sotiriou:2013qea; @Babichev:2013cya; @Cardoso:2013fwa; @Cardoso:2013opa; @Herdeiro:2014goa; @Antoniou:2017acq; @Sotiriou:2015pka].

A further complication in attempting to detect new fields through BH hair is the possibility that, even within the context of the same theory, only certain BHs might actually exhibit it. This was realized only recently, as the first models of BH [*scalarization*]{} appeared in the literature [@Silva:2017uqg; @Doneva:2017bvd]. For concreteness, consider the action $$\label{eq:gbaction}
S=\frac{1}{2}\int d^4x \sqrt{-g}\left(R-\frac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi+
f(\phi)\mathcal{G}\right),$$ where $f$ is some function of $\phi$, and where we have also set (as in the rest of this paper) $8\pi G = c =1$. Varying the action with respect to $\phi$ yields $$\label{eq:ScalarEq}
{\mathop{}\!\mathbin\Box}\phi = - f'(\phi)\mathcal{G},$$ where $f'(\phi)\equiv df/d\phi$. Assume that $f'(\phi_0)=0$, for some constant $\phi_0$. Then solutions with $\phi=\phi_0$ are admissible and they are also solutions of GR. A no-hair theorem [@Silva:2017uqg] ensures that they are unique if they are stationary, provided that $f''(\phi)\mathcal{G} < 0$.

The fact that GR BHs are stationary solutions to this theory is not sufficient to conclude that there are no observable deviations from GR, as the perturbations over these solutions do not generally obey the GR field equations [@Barausse:2008xv]. These perturbations may even grow unstable, thus rendering the GR solutions irrelevant. Indeed, one can think of $-f''\mathcal{G}$ as the (square of the) mass of the scalar perturbation on a fixed background. Hence, the condition above ensures that this effective (squared) mass is positive. If the condition is violated and the effective (squared) mass becomes sufficiently negative, the GR solutions suffer a tachyonic instability and the scalar develops a nontrivial profile.

A similar scalarization effect was shown to occur for neutron stars in a different class of scalar-tensor theories more that 25 years ago [@Damour:1993hw], and is triggered when the star compactness reaches a critical threshold. Related “dynamical” scalarization effects [@ST1; @ST2; @ST3; @ST4] are present in the same theories for neutron star binaries, whenever their separation is sufficiently small (or the binary’s “compactness” sufficiently large). However, in the class of theories considered in [@Damour:1993hw; @ST1; @ST2; @ST3; @ST4], scalarization is not present without matter, and BHs are vacuum solutions.[^2]

Black hole scalarization is fairly well understood. It starts as a linear tachyonic instability and, as such, its onset is controlled only by terms that contribute to linear perturbations around GR solutions. In this sense, action  with $f(\phi)=\eta\phi^2/2$ is sufficient to study the onset of scalarization [@Silva:2017uqg; @Andreou:2019ikc]. As the instability develops and the scalar grows, nonlinear terms become increasingly important and eventually quench the instability. Hence, the endpoint and properties of the scalarized solutions are actually controlled by the nonlinear interactions of the scalar [@Silva:2018qhn; @Macedo:2019sem]. A characteristic example is that in models with different nonlinear interactions, scalarized solutions have different stability properties [@Blazquez-Salcedo:2018jnn; @Silva:2018qhn].

Here we will focus exclusively on the onset of scalarization, so we will restrict our attention to quadratic scalar GB (qsGB) gravity, i.e. $f(\phi)=\eta\phi^2/2$ (without loss of generality [@Andreou:2019ikc]). The effective (squared) mass of the scalar on a fixed background is then $$\label{meff}
    \mu^2_{\text{\scshape\tiny{eff}}} = -\eta \mathcal{G}.$$ For the Schwarzschild solution, one has $\mathcal{G}=48 M^2/r^6$, which is always positive and decreasing with $r$, and which yields the horizon value $\mathcal{G}(r=2M)=3/(4M^4)$. Hence, a tachyonic instability only occurs for $\eta>0$, and the instability is expected to be more violent for smaller masses.[^3] This is why the focus in the literature so far has been on $\eta>0$ (or the equivalent condition in more complicated models). However, for a Kerr BH of mass $M$ and spin parameter $a$ in Boyer-Lindquist $(t,r,\theta,\varphi)$ one has $$\label{eq:gbterm}
\mathcal{G}_{\text{\scshape\tiny{Kerr}}}=\frac{48 M^2}{(r^2+\chi^2)^6}\left(r^6-15r^4\chi^2+15r^2\chi^4-\chi^6\right)$$ where, for brevity, $\chi\equiv a\cos\theta$. Clearly, $\mathcal{G}_{\text{\scshape\tiny{Kerr}}}$ is not monotonic, and can even become negative close to the horizon. This explains the results of [@Cunha:2019dwb; @Collodel:2019kkx], where it was shown that rotation suppresses scalarization for $\eta>0$.

In this Letter we focus on $\eta<0$, which yields a real effective mass $\mu_{\text{\scshape\tiny{eff}}}$ for low BH spins, but which can yield an imaginary $\mu_{\text{\scshape\tiny{eff}}}$ for high spins. We investigate the behavior of linear scalar perturbations to the GR solution by evolving Eq.  on a Kerr background, with the goal of assessing for what BH spins and couplings $\eta$ the perturbations become unstable. Indeed, at least two possible instability mechanisms may be at play in Eq. . The first is the tachyonic instability associated to spontaneous scalarization, mentioned above. The second could be a superradiant instability, which is known to exist at high spins for constant real masses [@Damour:1976; @Zouros:1979iw; @Detweiler:1980uk; @Dolan:2007mj], and potentially also for non-constant effective masses [@Dima:2020rzg] such as the one of Eq. . Superradiance occurs when bosonic waves with non-vanishing angular momentum are amplified when scattered by a spinning BH, at the expense of the rotational energy of the BH, which as a result spins down. For massive bosons, superradiant scattering can develop into an instability because the field is confined near the BH by its own mass.

It should be stressed that, in principle, both instabilites could be present. However, they have distinct features (timescales, the angular momenta involved, dependence on the BH spin). We show below that the tachyonic instability is by far the dominant effect for $\eta<0$. More broadly, our results strongly suggest that there exist theories in which scalarization occurs only for rapidly rotating BHs.

***Methodology:*** For $f(\phi)=0$ and over a Kerr background, Eq.  separates into ordinary differential equations when $\phi$ is decomposed onto a basis of spheroidal harmonics. However, the choice $f(\phi)=\eta \phi^2/2$ yields an intrinsically non-separable equation. We therefore resort to a time-domain numerical integration of this equation, by using techniques akin to those presented in [@Dima:2020rzg; @Dolan:2007mj], to which we refer for more details.

In brief, the idea is to project Eq.  onto a basis of *spherical*[^4] harmonics $\mathbf{Y}_{lm}$, which yields 1+1 evolutions equations (in $t$ and $r$) for the components of the scalar field, $$\label{eq:decomp}
\psi_{lm}(t,r)\equiv \int  \phi\, r \mathbf{Y}_{lm}\mathbf{Y}^*_{lm}d\Omega\equiv \langle l m | \phi\, r  | l m\rangle\,.$$ These equations are coupled and given explicitly by $$\begin{aligned}
\label{eq:perteq}
&\left[(r^2+a^2)^2-a^2\Delta(1- c^m_{ll})\right]\ddot{\psi}_{l}+a^2\Delta(c_{l,l+2}^m\ddot{\psi}_{l+2}+\nonumber\\
&+c_{l,l-2}^m\ddot{\psi}_{l-2})+4iamMr\dot{\psi}_l \nonumber\\
&-(r^2+a^2)^2\psi''_l-(2iam(r^2+a^2)-2a^2\frac{\Delta}{r})\psi'_l\nonumber\\
&+\Delta [ l(l+1)+\frac{2M}{r}-\frac{a^2}{r^2}+\frac{2iam}{r} ]\psi_l\nonumber\\
&+\sum_{j}\Delta \langle l m | \mu^2_{\text{\scshape\tiny{eff}}} (r^2+\chi^2) | j m \rangle\psi_j=0 \,,\\
&\Delta\equiv r^2-2 M r +a^2    \,,\\
\label{eq:coscoup}
&c^m_{jl}\equiv \langle l m | \cos^2\!\theta | j m\rangle \nonumber\\&= \frac{\delta_{lj}}{3}+\frac{2}{3}\sqrt{\frac{2j+1}{2l+1}}\langle j,2,m,0|l,m\rangle \cdot \langle j,2,0,0|l,0\rangle\,,\end{aligned}$$ where $\langle j_1,j_2,m_1,m_2|j_3,m_3\rangle$ are the Clebsch-Gordan coefficients [@AngMom]. Note that the evolution of modes of different $m$ decouples because of the axisymmetry of the problem. Morever, because of reflection symmetry with respect to the origin, even-$l$ and odd-$l$ modes also decouple: the evolution of a mode $(l,m)$ is coupled to that of all the modes $(l+2k,m)$, with $k=1,2,3,\ldots$.

To numerically evolve the system , we discretize the spatial grid and use a method of lines. By integrating in time using a fourth order explicit Runge-Kutta time-step inside the computational grid (as done e.g in [@Dima:2020rzg]), it becomes apparent that the equations are stiff for large $\eta$, and that the numerical integration becomes unstable. To overcome this problem, we have used an Implicit-Explicit (IMEX) Runge-Kutta solver with adaptive time step, namely the IMEX-SSP3(3,3,2) and IMEX-SSP(4,3,3) schemes of [@Palenzuela:2008sf]. Note that implicit methods [@PareschiRusso], while effective at dealing with stiff problems, are typically less accurate and more computationally expensive. However, implicit-explicit algorithms, by employing explicit steps for the non-stiff terms and implicit steps only for the stiff ones, can tackle stiff problems with limited computational overhead. We successfully compared our code to results from both frequency-domain techniques [@Berti:2009kk] and similar time-domain codes [@Dolan:2012yt]. Our implementation was also tested by analyzing the convergence of the results (and their overall robustness) vs time-step and spatial-grid resolution.

***Results:*** To investigate the possible presence of an instability, we evolve the scalar field by integrating the system given by Eq. , with $l$ ranging from 0 to $l_{\rm max}=23$ and $|m|\leq l$, and with Gaussian initial conditions for each mode $\psi_{lm}$. The results are robust against the choice of the cutoff $l_{\rm max}$ – as long as that is sufficiently large – and initial conditions, which only affect the early transient evolution of the scalar and not the unstable growth phase, if present. We consider BH spins $a/M \sim 0.6$ – $0.999$ and qsGB coupling $\abs{\eta}/ M^2 \sim 0.1$ – $10^{5}$.

From the simulations showing an exponential scalar growth, we extract the instability timescale $\tau$ of the reconstructed field $|\phi|=\left(\sum_{lm} |\psi_{lm}|^2 \right)^{1/2}\propto \exp(t/\tau)$ by fitting the time evolution of the scalar’s amplitude after the initial transient. The contours in Fig. \[fig:tacReconstructed\] show $\tau^{-1}$ as a function of $a/M$ and $\abs{\eta} /M^2$. The instability becomes stronger as either the spin or the coupling increases. Moreover, there appears to be a minimum spin $a_{\rm min}\simeq 0.64$, below which the instability does not occur. The thick green line denotes the combinations of parameters for which the instability disappears (i.e. $\tau\to\infty$). To understand how the different modes affect the instability threshold, we perform the following test. We select a specific mode, and forcefully set to zero all the couplings that it has with other modes in Eq. , therefore overriding the non-separability of the problem. This ensures that if the instability is due to this mode, it will not leak into the others through the couplings, and viceversa. We then let the system evolve, choosing Gaussian initial data for the chosen mode only. This test reveals that the instability around the critical line is dominated by the $\psi_{10}$ mode, as can be seen in Fig. \[fig:tacReconstructed\].

![Instability timescale $\tau$ (color code) for the reconstructed field as a function of spin and GB coupling. The instability threshold is shown by a solid green line. Note that the unstable region coincides with that of the $l=1,m=0$ mode, whose boundary is indistinguishable from the green solid line, but is larger than the unstable region when the $l=1,m=0$ mode is suppressed (blue dashed line).[]{data-label="fig:tacReconstructed"}](insta_main.png)

![ Same as in Fig. \[fig:tacReconstructed\], but for even modes with $m=0$ only. The instability threshold is shown by a solid green line, while that for the $l=m=0$ mode is shown by a red dash-dotted line. The blue dashed line represents the instability threshold for all the higher even modes ($l\geq2$, $m=0$). The small mismatch between the green solid and blue dashed lines is an artifact of the numerical interpolation.[]{data-label="fig:tacSpherical"}](insta_sph.png)

Even and odd parity modes automatically decouple in Eq. . In Fig. \[fig:tacSpherical\], we show the instability timescales for even modes with $m=0$. The instability occurs in a smaller parameter region than in Fig. \[fig:tacReconstructed\], which was dominated by the $l=1$ mode instability. Another striking difference with Fig. \[fig:tacReconstructed\] is the kink in the critical line separating stable and unstable regions. To explain this feature, we have performed the same, brute force decoupling test discussed above, but this time for the $l=m=0$ mode. The red dash-dotted line in Fig. \[fig:tacSpherical\] corresponds to the instability threshold for this mode, and it coincides at high spins/small couplings with the critical line for all even modes (green solid line). Moreover, suppressing the coupling of the spherical mode $l=m=0$ to the GB invariant and setting its initial conditions to zero, we evolved the non-spherical even modes ($l\geq2$, $m=0$). The result of this test is also shown in Fig. \[fig:tacSpherical\] as a blue dashed line, which matches the solid green line for lower spins/higher couplings. These two tests therefore explain completely the kink in Fig \[fig:tacSpherical\].

We now proceed to examine whether the instability is dominantly tachyonic or powered by superradiance. The growth times, as shown in Figs. \[fig:tacReconstructed\] and \[fig:tacSpherical\], can be as small as $\sim 0.02 M$. This seems to favor a tachyonic origin, as superradiance acts on longer timescales (see e.g. [@Dima:2020rzg; @Dolan:2007mj]). Moreover, the fact that the spherical mode $l=m=0$ can be unstable bodes ill for superradiance, as this mode can never satisfy the superradiance condition $\omega< m \Omega$ (with $\omega$ and $\Omega$ respectively the wave and horizon angular frequencies).

One may naively expect the spherical mode $l=m=0$ not to suffer from a tachyonic instability either, since $\mu^2_{\text{\scshape\tiny{eff}}}=-\eta \mathcal{G}$ is positive everywhere in a Schwarzschild spacetime when $\eta<0$ (as considered here). However, the (squared) effective mass for the $l=m=0$ mode is actually $-\eta \langle00| \mathcal{G}_{\rm Kerr}|00\rangle$, which only matches the naive estimate $-\eta \mathcal{G}_{\rm Schwarzschild}$ at leading order in spin, correcting it by terms ${\cal O}(a^2)$. This explains, in particular, why no instability occurs at low spins in Figs. \[fig:tacReconstructed\] and \[fig:tacSpherical\].

To further confirm the tachyonic nature of the instabilities, we have conducted the following test. We re-ran our simulations with the (squared) effective mass replaced by its absolute value, $\mu^2_{\text{\scshape\tiny{eff}}} \rightarrow |\mu^2_{\text{\scshape\tiny{eff}}}|$. This is enough to suppress the instabilities, and further shows that the latter were due to the change of sign of the GB invariant close to the horizon.

![Energy flux $\mathcal{F}_E$ through the BH horizon vs time, for $a=0.99 M$. The blue, orange and magenta lines correspond respectively to $\eta=-10 M^2$, to a tachyonic mass $\mu M=i$, and to a constant, real mass $\mu M=0.42$. The inset zooms on the constant, real mass flux (of which we show a moving average to decrease the oscillations caused by the dynamics). That flux is negative, signaling energy extraction from the BH, as expected for superradiant instabilities.[]{data-label="fig:fluxes"}](fluxes.png)

One can also look at the scalar fluxes through the event horizon after the initial transient. In Fig. \[fig:fluxes\], we compare the scalar field’s energy flux through the horizon for $\eta=-10 M^2$ (blue) vs the same fluxes for minimally coupled scalar fields with imaginary (orange) and real (magenta) constant masses. Clearly, the flux for a scalar coupled to the GB invariant resembles more closely the tachyonic (i.e. imaginary mass) scalar field evolution, both in timescale and sign. Note that the constant, real mass case, whose evolution is unstable due to superradiance, shows a slower growth and negative energy fluxes. The latter are indeed the hallmark of a superradiant instability, which [*removes*]{} rotational energy and angular momentum from the BH.

The most plausible explanation for why Kerr BHs in qsGB do not suffer from superradiant instabilities seems to be the rapid falloff of the GB invariant (thus of the effective mass) at large distances, $\mathcal{G}(r\rightarrow \infty) \sim 1/r^6$. Scalar perturbations with a position-dependent mass were studied in [@Dima:2020rzg], which showed that a steep decay of the mass with distance quenches the superradiant instability. This happens because the effective potential for scalar perturbations does not develop wells, and thus quasi-bound states, unless the mass remains relatively constant till at least $r\sim 2-3 M$ [@Dima:2020rzg].

***Conclusions:*** We have shown that a coupling, with a suitable sign, between a scalar and the GB invariant can lead to an instability triggered by rapid rotation. We have also demonstrated that this instability is not related to superradiance, but is instead tachyonic in nature. Nonlinear effects, which our approach does not capture, are expected to quench that instability and lead to a BH with scalar hair. The process is analogous to the more conventional spontaneous scalarization, but the threshold is controlled by the black hole rotation instead of its curvature.

The action that we use is sufficient for studying the onset of the instability that we have found for BHs. However, the endpoint of this instability, and hence the amount of hair a BH would carry, will strongly depend on nonlinear (self)interactions.[^5] There is no obvious reason to believe that this instability is restricted to BHs, and it could well affect rapidly rotating stars as well. Hence, our results demonstrate that there is a broad class of theories where rotation might control deviations from GR. Our findings also have clear implications for searches of new physics in the strong-field regime. Black hole scalar hair induces vacuum dipole gravitational emission, which is potentially observable in the low frequency inspiral of binary system by gravitational wave interferometers [@Barausse:2016eii; @Toubiana:2020vtf], deviations from GR in the spectrum of the gravitational wave ringdown [@Berti:2016lat] or in the electromagnetic spectrum of accretion disks [@Bambi:2011jq], and it may also impact the black hole shadow observed by the Event Horizon Telescope [@Akiyama:2019cqa].

We stress that we are not aware of any observational upper bounds on $\eta$, which we therefore allow here to reach very high values, for illustrative purposes and in order to excite higher modes. Note that slowly rotating black holes in qsGB would be identical to their GR counterpart. Compact stars can scalarize for $\eta<0$ [@Silva:2017uqg] and hence yield constraints. However, this effect could easily be quenched by adding a coupling between the scalar field and the Ricci scalar [@Andreou:2019ikc; @Ventagli:2020rnx]. The latter might be necessary to get a sensible cosmology [@Antoniou:2020nax], and would have no effect for black holes, thus leaving our analysis unaffected.

[[**[*Acknowledgments.*]{}**]{}]{} We thank Carlos Palenzuela and Miguel Bezares for insightful advice on technical aspects of the IMEX schemes and their validation, and Hector O. Silva for useful discussions on black hole scalarization. A.D., E.B. and N.F. acknowledge financial support provided under the European Union’s H2020 ERC Consolidator Grant “GRavity from Astrophysical to Microscopic Scales” grant agreement no. GRAMS-815673. T.P.S. acknowledges partial support from the STFC Consolidated Grant No. ST/P000703/1. We also acknowledge networking support by the COST Action GWverse Grant No. CA16104.

[^1]: No-hair theorems also exist for stars in shift-symmetric scalar tensor theories [@starsNohair1; @starsNohair2; @Lehebel:2017fag; @Yagi:2015oca].

[^2]: Black holes can scalarize if they have matter in their vicinity [@Cardoso:2013fwa; @Cardoso:2013opa], but the densities necessary to obtain a measurable effect are probably astrophysically unrealistic.

[^3]: Note that in curved spacetimes $\mu^2_{\text{\scshape\tiny{eff}}}$ can be somewhat negative without necessarily developing a tachyonic instability.

[^4]: There is no advantage in using spheroidal harmonics, for which analytic expressions are unavailable, as they do not lead to a separable equation anyway.

[^5]: Stationary scalarized black hole solutions that constitute the endpoint of the instabilty will be presented elsewhere [@prep].
