The onset of spontaneous scalarization in generalised scalar-tensor theories

In gravity theories that exhibit spontaneous scalarization, astrophysical objects are identical to their general relativistic counterpart until they reach a certain threshold in compactness or curvature. Beyond this threshold, they acquire a non-trivial scalar configuration, which also affects their structure. The onset of scalarization is controlled only by terms that contribute to linear perturbation around solutions of general relativity. The complete set of these terms has been identified for generalized scalar-tensor theories. Stepping on this result, we study the onset on scalarization in generalized scalar-tensor theories and determine the relevant thresholds in terms of the contributing coupling constants and the properties of the compact object.


I. INTRODUCTION
Gravitational wave astronomy now offers a way to observe the inspiral and coalescence of pairs of compact astrophysical objects, either black holes or neutron stars [1,2].This sheds light on two unexplored corners of gravity.First, during the late stages of the inspiral, the objects move at speeds close to the speed of light, allowing to test the dynamical regime of gravity.Second, the coalescence of compact objects involves high curvatures and non-linearity plays a significant role.Strong field effects can emerge progressively as the curvature increases.They can also appear abruptly, in a phase transition process.The study of these phase transitions has a rich history, particularly in the context of scalar-tensor gravity.It was originally investigated by Damour and Esposito-Farèse (DEF) [3] in the context of neutron stars, and later dubbed spontaneous scalarization (in analogy with magnetization).
The DEF model relies on a specific coupling between the scalar field and the metric.The model admits several branches of solutions, among which all general relativity (GR) solutions with a trivial scalar field.At low curvatures, GR solutions are stable against scalar perturbations.However, in the strong field regime, the scalar field acquires a tachyonic effective mass that destabilizes the solution [4].This instability is ultimately quenched by non-linear terms, resulting in a neutron star with a nontrivial scalar profile [5].The DEF model was studied thoroughly, especially in the light of binary pulsar constraints (e.g., [6][7][8]).These constraints have severely constrained the model in its original formulation, but they can be straightforwardly evaded by giving the scalar field a bare mass [9].
It has recently been shown that the DEF coupling is not unique in generating scalarization.In particular, scalar-Gauss-Bonnet gravity exhibits similar effects, both for neutron stars and black holes [10,11].A generic property of these models is that linear terms in the field equations determine the onset of scalarization, while nonlinear terms keep the instability under control and deter-mine the endpoint of the process [12][13][14].It is also possible to extend the concepts of scalarization to different field contents (e.g., [15][16][17]).At this point, it is legitimate to wonder what is the most generic model that can trigger scalarization.Within the context of Horndeski theory, this question was answered in [18].Horndeski theory offers a robust framework of models which retain second-order field equations, and are therefore free of Ostrogradsky ghost [19].Reference [18] listed all the terms that can play a role in the onset of scalarization.In this paper, we determine quantitatively the bounds on the couplings that allow scalarization to take place.We combine the effect of all relevant couplings simultaneously and also examine how the structure of the compact object affects the threshold of the instability.
The paper is organized as follows.In sec.II, we describe the framework that allows us to identify the onset of an instability.Then, in sec.III, we investigate the effect of the parameters that can generate an effective mass.Section IV is devoted to the influence of the background on the instability (mass of the object, equation of state).In sec.V, we focus on the parameter that deforms the effective metric in which perturbations propagate.Finally, we conclude in sec.VI.

II. SETUP
As mentioned above, in the framework of Horndeski theory, all the terms that can trigger a tachyonic instability (or play a role in its onset) were listed in [18].The minimal action containing all these terms reads where X = −∇ µ φ∇ µ φ/2, κ = 8πG/c 4 and G is the Gauss-Bonnet invariant: arXiv:2006.01153v2 [gr-qc] 12 Jun 2020 m φ is the bare mass of the scalar field, while γ, β and α parametrize deviations with respect to GR; β is dimensionless, while α and γ have the dimension of a length squared.S M denotes the matter action, where matter is assumed to couple minimally to the metric only: we are working in the so-called Jordan frame.Note that β is defined in order to match the notation used in the case of the (linearized) DEF model.The precise relation to the original DEF model is discussed in detail in Ref. [18].
In general, the action in eq. ( 1) admits all GR solutions with a trivial scalar field configuration φ = 0.One can obtain from this action any theory within the Horndeski class that admits GR solution with a constant scalar (not necessarily zero) and exhibits spontaneous scalarization: first one can shift the scalar by a suitable constant and then supplement the action with the desired nonlinear interactions.Our goal will be to investigate whether GR solutions with φ = 0 are stable or not, by focusing on the perturbations of the scalar field.The scalar field equation associated with the action (1) reads where the effective metric and the mass term are respectively It is clear from these equations that β and α generate an effective mass for φ in a curved background.On the other hand, γ determines the effective metric that defines the d'Alembertian which acts on the scalar field perturbations.Equation ( 3) is linear by construction, since we kept in the action only the terms that contribute linearly.This approach is valid when focusing on the onset of the scalarization, when linear terms dominate.In order to determine the final state of the process, one needs to include non-linear terms as well.Additionally, we work in the decoupling limit where we only perturb the scalar field.These perturbations will eventually back-react onto the metric, and a consistent analysis (beyond the onset of the instability) should thus include metric perturbations.GR vacuum solutions are Ricci flat.Therefore, it is immediate from action (1) that only the Gauss-Bonnet coupling α controls the scalarization of vacuum solutions (in particular electrically neutral black holes).Since we are interested in the combined effect of the parameters involved in (1), we focus on neutron stars, where matter is present under the form of a perfect fluid.We further consider a static and spherically symmetric background spacetime: The metric functions h and f are determined as solutions of an equivalent to the Tolman-Oppenheimer-Volkoff system of equations [20,21], together with the pressure P and energy density of the perfect fluid that composes the star.These equations can be found in appendix A.
In this framework, one has to specify some equation of state P ( ).We use two equations of state, SLy and MPA1 [22], both favored by LIGO-Virgo tidal measurements [1], which seem to prefer soft equations of state.We work in units where c = 1, G = 1 and M = 1.Thanks to spherical symmetry, we can decompose the scalar perturbation on the basis of spherical harmonics: We will focus on the breather mode, = m = 0, which is the first one to exhibit instability when it is present.In order to make the scalar field equation more transparent, we rescale this mode according to φ00 (t, r) = K(r)σ(t, r) with and we trade off the radial coordinate r for a new one, r * , defined through Equation ( 3) then takes the following form: where V eff depends on the parameters of the action (1) as well as on the background geometry.Its full expression can be found in appendix B. We further focus on exponentially growing perturbations: σ(t, r * ) = σ(r * )e ωt , with ω > 0.1 Equation ( 10) then boils down to a Schrödinger equation: where V eff is clearly an effective potential, and −(ω/c) 2 plays the role of the energy of the perturbation.The existence of a bound state for V eff with 'energy' E 0 < 0 implies the existence of an instability, with characteristic growth rate ω = c √ −E 0 .Our strategy will thus be the following: we start from values of the parameters for which the theory reduces to GR with a minimally coupled scalar field and hence there cannot be any instability.
We gradually increase the parameters, thus progressively deforming the potential.Whenever a bound state appears for V eff , we identify it with a new unstable mode for φ.By continuity, when deforming the potential, a new bound state will appear with a vanishing energy, E 0 = 0. Therefore, we solve eq. ( 11) for ω = 0, while scanning the parameters β, α and γ.This is less intuitive than choosing a set of parameters and scanning ω, but the final result is equivalent, and the procedure is easier to implement.
Equation (11) will admit a solution for any set of values of the parameters.Among these solutions, we identify as a bound state those with dσ/dr * → 0 for r * → +∞. 2hysically, this is necessary for the scalar perturbation to be localized in space.In terms of quantum mechanics, this corresponds to the fact that K(r * )σ(r * ) has to be square integrable.Note that, contrary to the naive expectation from eq. ( 11), it is K(r * )σ(r * ) that should be interpreted as the wave function, rather than σ(r * ) itself; this is very similar to the 1/r rescaling of the wave functions that allows to solve for the bound states of 3D spherically symmetric quantum wells (see e.g., secs.14-16 of [23]).

III. CHANGING THE EFFECTIVE MASS
As mentioned above, the main terms that can contribute to the effective mass m 2 are the bare mass term of the scalar field, the coupling between φ and the Ricci scalar, and the coupling between φ and the Gauss-Bonnet invariant.These terms are parametrized by three constants: m φ , β and α.Although the terms proportional to γ can affect the instability threshold, they contribute only as a multiplicative constant.Therefore, we will set γ = 0 throughout this section, and explore the role of this parameter in full detail in Sec.V. Several works already investigated the influence of each of the parameters m φ , β and α separately (e.g., [3,4,[10][11][12][13][14][24][25][26]).We study how varying several parameters simultaneously affects the threshold; in this way, we explore much wider regions of the parameter space.Most of our results are presented as 2D plots, where we freeze all parameters but two, and show the stable/unstable regions.

A. Coupling to curvature invariants
We first consider a vanishing bare mass, m φ = 0.The model is then parametrized by β and α only.The background we consider is a neutron star described by the SLy equation of state.We choose its central energy density to FIG. 1. Stable and unstable regions in the (β, α) space for a light star (M = 1.12 M , SLy equation of state).In the 2D plot, each line is labeled according to the number of nodes n of the corresponding unstable mode.Inside the white region, where the point (β, α) = (0, 0) lies, the GR solution is stable.Every line crossed while moving away from the origin corresponds to the appearance of a new unstable mode; any point in parameter space that lies within a grey region corresponds to an unstable solution.The lower panel shows |dσ/dr| at rmax when varying β in the same range as the 2D plot, with α = 0; it can be understood as a cut in the (β, α) plane along the β axis (each cusp corresponding to a line-crossing in the 2D plot).Similarly, the left panel shows a cut along the α axis.
be ρ c = 8.1 × 10 17 kg/m 3 , so that its gravitational mass is M = 1.12 M , which corresponds to the bottom of the mass range for observed neutron stars [27,28].The radius of the star is then R s = 11.7 km.The results are summarized in Fig. 1.The white area corresponds to the region of the parameter space where the background solution is stable.A new unstable mode appears when crossing each line while moving away from the origin.The lines are labeled with the number of nodes n of the associated mode, ranging from 0 to infinity.The n = 0 line is the boundary of the stable region.Any choice of parameters beyond this line will make the GR solution unstable.The left and bottom panels shows cuts in the (β, α) plane, along the α and β axes respectively; these panels actually reproduce known results, e.g. of [4] and [10].Notably, the scalarization threshold when only β is presents takes the well-known order of magnitude, β = −5.42.
To understand better the shape of the plot, especially along the β and α axes, we plot in Fig. 2 the Ricci scalar and the Gauss-Bonnet invariant.The Ricci scalar is always positive on this background.This is due to the fact FIG. 2. Ricci scalar and Gauss-Bonnet invariant for a light star (M = 1.12 M , SLy equation of state).The radial coordinate is rescaled by the radius of the star Rs = 11.7 km.The left panel shows that the Ricci scalar is non-negative everywhere; correspondingly, only β < 0 can lead to an instability.On the other hand, the Gauss-Bonnet invariant, shown in the right panel, is negative in the core of the star and positive towards its surface, leading to instabilities both for α < 0 and α > 0.
that we consider a relatively light neutron star, and due to the following relation: When the medium is not too dense, P and R > 0. We will see how this changes for a very dense star in sec.IV.As a consequence, only negative β can generate a negative effective mass.On the other hand, as shown in Fig. 2, the Gauss-Bonnet invariant is positive in some regions and negative in others.This is enough to trigger an instability when α becomes very negative or very positive, which is indeed what is observed in Fig. 1.
We notice, as expected, that the point (0, 0) is always inside the stable region.The tachyonic instability does not appear right away when β < 0 or α = 0, due to the curvature of spacetime.

B. Effect of the bare scalar mass
We now consider how the presence of a bare mass affects the results of the previous section.Fig. 3 shows the region of stability in the (β, α) plane when m φ = 1 in the system of units that we used, i.e., scalar particles have a mass of 1.33 × 10 −10 eV.The range of parameters in Fig. 3 is the same as in Fig. 1 in order to allow comparison.When one zooms out, Fig. 3 looks very similar to Fig. 1 (i.e., the same instability pattern remains valid, but it appears for higher values of the parameters).As can be seen from comparing Figs. 1 and 3, the stable region is widened in all directions.As expected, the presence of a bare mass stabilizes the solution (similarly, if the bare mass is tachyonic, the stability region shrinks).Therefore, even a very light bare mass for the scalar field is able to shield neutron stars from scalarization.This was already noted, e.g., in [9], where the effect of a bare mass in a cosmological setup is also discussed.FIG. 3. Stable and unstable regions in the (β, α) space for a bare mass of 1.33 × 10 −10 eV.As expected, the stable region is enlarged with respect to Fig. 1.Note that the range of the plot for β and α is the same as in Fig. 1 to facilitate comparison.

IV. CHANGING THE PROPERTIES OF THE STAR
Let us now examine how changing the background affects the stability.We first consider a more massive star, and then a different equation of state.

A. Mass of the star
We first increase the mass of the neutron star (and thus its compactness at the same time).We choose a central density of ρ c = 3.4 × 10 18 kg/m 3 , which corresponds to a mass of M = 2.04 M .This is the heaviest spherically symmetric star we can produce with the SLy equation of state; beyond this mass, the solutions become unstable already within GR.The results are displayed in Fig. 4. Since the curvature of the background increased with respect to Fig. 1, it is not surprising that the stable region shrunk.A more specific feature is that very positive values of β now also lead to an instability.This is due to the fact that the Ricci scalar now becomes negative towards the center of the star, as can be seen in Fig. 5.In terms of eq. ( 12), the energy density is no longer dominant with respect to the pressure in the extremely pressurized core of the neutron star, allowing R < 0. This effect was already noted in [29], and further studied in [30,31].

B. Equation of state
We now consider a different stellar model, the MPA1 equation of state [22].In order to compare with previous results, we keep the same mass as in Sec.III, M = 1.12 M .This corresponds to a central density of ρ c = 6.3 × 10 17 kg/m 3 .The stability region is shown FIG. 4. Stable and unstable regions in the (β, α) space for a heavy neutron star (M = 2.04 M , SLy equation of state).The vertical range for α is reduced with respect to Fig. 1 for readability.The wedge of stability in Fig. 1 has narrowed down to an island around (β, α) = (0, 0). in Fig. 6.By comparing Figs. 1 and 6, it is clear that the different choice of equation of state does not affect significantly the position of the stable and unstable regions.
Although the equation of state influences only mildly stars of similar mass, it can have indirect effects.Indeed, a softer equation of state will lead to a smaller pressure for a given energy density; thus, for a soft equation of state, the Ricci scalar ( 12) could remain positive in all configurations, discarding configurations like Fig. 4. Similarly, the equation of state can affect the range allowed for the parameter γ, as we will see in more detail in Sec.V.

V. CHANGING THE EFFECTIVE METRIC
We now return to the parameter γ and examine its role.As mentioned above, the terms controlled by γ cannot generate a tachyonic instability in the absence of the other couplings controlled by β, α or m φ .Nonetheless, FIG. 6. Stable and unstable regions in the (β, α) space for the MPA1 equation of state (M = 1.12 M ).The range for β and α is the same as in Fig. 1.We note that the choice of equation of state does not affect significantly the results.
choosing γ beyond a certain range leads to loss of hyperbolicity in the scalar field equation.Also the potential V eff does depend on γ, as can be seen from eq. (B1).We analyze these aspects below, before studying numerically the combined effect of γ and the other parameters.

A. Hyperbolicity
The parameter γ brings an additional contribution to the effective metric experienced by scalar perturbations, eq. ( 4).If this contribution exceeds a certain threshold and becomes dominant, the effective metric becomes elliptic, rendering the background unstable.We emphasize that this instability is distinct from the usual tachyonic instability associated with scalarization.Depending on the circumstances (see also below), it is either a gradient or a ghost instability.It is possible that this instability can be quenched nonlinearly, as is the case for conventional scalarization.Indeed, scalarization through a ghost instability has already been proposed in [32].Since here we are not including terms that are nonlinear in the scalar in our analysis, we cannot follow the development and potential quenching.
One can therefore view the analysis that follows in two different ways.Taking the conservative viewpoint, one may restrict to the well-controlled tachyonic scalarization.In this framework, our results allow to set bounds on the parameter γ.In a more open-minded perspective, which certainly deserves further investigation in the future, we are setting bounds beyond which ghost or gradient scalarization can be triggered.
Given that we are working on a GR background, we can make use of Einstein equations.The inverse of gµν for a perfect fluid is then FIG. 7. Hyperbolicity of the effective metric in the (M, γ) plane (SLy equation of state).The effective metric is hyperbolic only within the white region.The mass ranges from the putative lowest neutron star mass [27,28] to the maximal mass achieved with the SLy equation of state.The range of hyperbolicity narrows down when increasing the mass of the star.
where u µ is the 4-velocity of the fluid.Over the spherically symmetric background that we study, the determinant of the effective metric reads On a given background, the pressure and energy density are maximal at the center of the star, where we label their value as P c and c .The determinant of the physical metric g is always negative.Thus, the effective metric loses hyperbolicity either when κγ becomes larger than 1/P c or when κγ becomes more negative than −1/ c .To summarize, hyperbolicity of the effective metric requires that Reference [33] already noted the existence of the upper bound in a similar context.In the numerical analysis below, we will take the conservative approach and restrict to this range.For a given equation of state, the bounds in eq. ( 15) can be reformulated in terms of the mass of the star, M .This is shown in Fig. 7.The range in which gµν is hyperbolic closes up around γ = 0 when increasing the mass.The limits presented in Fig. 7 depend on the equation of state, but only mildly.In the framework of tachyonic scalarization, we can use these limits to put an absolute bound on γ.For the SLy equation of state, − 18.7 km 2 < γ < 34.9 km 2 .
These bounds should only be taken as order of magnitude estimates; they have been established in the decoupling limit and rely on a specific equation of state.However, a more detailed study would allow to put very precise constraints on γ.We are not aware of previously established bounds on this parameter (except [33]).Note that the analysis presented in this paragraph is entirely independent on the other parameters, which do not play a role in the effective metric.

B. Integral of the effective potential
Before moving to the numerical results, we out an difference between black holes and stars.In the case of black holes, the scalar field equation can also be brought to the form (11).In both cases (black hole and neutron star), we impose that dσ/dr * vanishes at large r * so that φ 2 is finite and the perturbation initially contains a finite amount of energy.The difference between neutron stars and black holes appears on the other side of the r * range.In the case of neutron star, r * goes down to 0, where φ = K σ ∼ φ 0 /r * for some constant φ 0 , unless σ(0) = 0. We do not want φ to diverge at the center of the star, so we impose σ(0) = 0.In the black hole case on the other hand, r * goes down to −∞ and nothing particular happens at r * = 0.The condition for the perturbation to be physical is then that dσ/dr * vanishes for r * → −∞, similarly to what happens at large r * .
It is then possible to establish an exact sufficient condition for the presence of an instability in the black hole case.Indeed, the function σ then respects the hypotheses of the theorem established in ref. [34]: if +∞ −∞ V eff (r * )dr * < 0, then V eff admits at least one bound state.
However, in the case of a star, due to the different boundary conditions, +∞ 0 V eff (r * )dr * can become negative -even infinitely negative -while V eff does not possess any bound state.This is equivalent to the 3Dspherically symmetric quantum well; consider such a well with depth −V 0 and width a.It admits a bound state only for V 0 a 2 > N (N is some constant that depends on the mass of the quantum particles), while the integral of the potential is −V 0 a < 0. Choosing the scaling V 0 = N/2/a 2 , the integral of the potential is then becoming infinitely negative for a → 0, while no bound state exists.
Hereinafter, V eff (r * )dr * indeed becomes infinitely negative in the limit where γ approaches one of the bounds of eq. ( 15).However, this does not necessarily mean that infinitely many bound states develop close to these boundaries.This seems to be true close to the upper bound, but not close to the lower one, as we will see in Fig. 8.

C. Numerical analysis
We now vary γ systematically in the range allowed by eq. ( 15), and β and α in similar ranges as in the previous ).It is also natural that the lines β = 0 and α = 0 (vertical axes in figs.8 and 9) correspond to stable configurations, because γ cannot create a tachyonic effective mass on its own.The boundary of the stable region is rather insensitive to the parameter γ; it does evolve slightly with the value of γ, but this can be seen only when zooming around smaller values of β and α with respect to figs.8 and 9. Close to both bounds of the γ range, V eff (r * )dr * diverges, but as explained in the previous section, this does not necessarily mean that infinitely many bound states should appear (or even that one bound state exists).Indeed, nothing particular happens when approaching the lower bound, while unstable FIG. 10.Summary of the stability (n = 0) contours in the (β, α) space for comparison.The white region corresponds to the region of stability for a heavy star, M = 2.04 M , using the SLy equation of state.The light grey region (together with the white region) is the stability region for M = 1.12 M , still with the SLy equation of state.The intermediate grey region (together with the previous paler regions) corresponds to a mass M = 1.12 M and the MPA1 equation of state.The dark grey region (again, with paler regions) corresponds to a mass M = 1.12 M and the SLy equation of state, together with a bare scalar mass of 1.33 × 10 −10 eV.Finally, the black region is unstable for all the previous models.
modes pile up when approaching the upper bound; both scenarios are allowed.

VI. CONCLUSION
We have investigated exhaustively the effect of all the terms that can play a role in the onset of spontaneous scalarization, within the context of Horndeski theory.Our analysis has identified the role of each term but has also revealed their combined effects.When taking each terms separately, our results agree with previous results regarding scalarization thresholds.More generally and when all terms are present, we have found that a very small bare mass suffices to stabilize GR solutions, and that the scalarization thresholds are only mildly sensitive to the choice of equation of state.
Our analysis allowed us to explore, for the first time, the multi-dimensional parameter space and provide scalarization thresholds that depend on more than one coupling.In Fig. 10 we presented in a single plot a summary of most of the stability contours presented in this paper.We expect that our multi-parameter analysis will be particularly useful in building viable scalarization models.Observational constraints on these models will come from different compact objects, such as black holes and neutron stars, and also from low curvature observations.The Ricci scalar and the Gauss-Bonnet invariant, which appear in the terms that contribute to the onset of scalarization, scale differently with curvature and are not sign-definite.Hence, by choosing the sign and mag-nitude of different couplings constants in the minimal action of eq. ( 1) one can construct models in which only certain types of compact objects are scalarized and also control the behaviour of the model at lower curvatures.A characteristic example in this direction has recently been discussed in Ref. [35]: a suitable choice of parameters yielded a black hole scalarization model with GR as a cosmological attractor.
One of the striking features revealed by our analysis is the role of the effective metric in which scalar perturbations propagate.It is controlled by a single coupling constant, γ.There exists a threshold beyond which the effective metric loses hyperbolicity.In the framework of tachyonic scalarization, we interpret this threshold as an absolute bound on the parameter γ.It is restricted to a rather narrow range (roughly, it should remain small with respect to the characteristic length of curvature).Therefore, it has a very limited effect on the threshold of tachyonic scalarization.On the other hand, the loss of hyperbolicity can be seen as an alternative instability that could lead to scalarization, in line with what was proposed in Ref. [32] in a more restricted setup.It would be interesting to understand if such an instability can indeed be controlled and give rise to a sensible scalarization process.
More generally, the next natural step is to take into account non-linearities in the scalar field, which are excluded from the minimal action (1) by definition.They determine the end state of any scalarization instability and control the deviations of scalarized compact objects from their GR counterparts.A detailed study would encompass the existence of scalarized solutions, their stability, their transition from the GR branch to the scalarized branch, and the associated observational signatures and constraints.We leave these for future work.

FIG. 5 .
FIG.5.Ricci scalar and Gauss-Bonnet invariant for a heavy star (M = 2.04 M , SLy equation of state).Now, the left panel shows that R becomes negative inside the core of the star, allowing instabilities for both signs of β.The Gauss-Bonnet invariant (right panel) still behaves as in the case of a light star, Fig.2.

FIG. 8 .
FIG.8.Stable and unstable regions in the (β, γ) space (M = 1.12 M , SLy equation of state).The region of stability is rather unaffected by a change in γ.The bound states pile up close to the upper bound for γ.