Scalarization is a mechanism that endows strongly selfgravitating bodies such as neutron stars and black holes with a scalarfield configuration. It resembles a phase transition in that the scalar configuration appears only when a certain quantity that characterizes the compact object for example its compactness or spin is beyond a threshold. A critical and comprehensive review of scalarization including the mechanism itself theories that exhibit it its manifestation in neutron stars black holes and their binaries potential extension to other fields and a thorough discussion of future perspectives is provided.
['General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena', 'High Energy Physics - Theory']
Spontaneous scalarization
{'Spontaneous scalarization': 'Daniela. D. Doneva ∗  \nTheoretical Astrophysics, Eberhard Karls University of Tubingen, 72076 Tubingen, Germany INRNE - Bulgarian Academy of Sciences, 1784 Sofia, Bulgaria  \nFethi M. Ramazano˘glu †  \nDepartment of Physics, Ko¸c University, Rumelifeneri Yolu, 34450 Sariyer, Istanbul, Turkey  \nHector O. Silva ‡  \nMax Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Muhlenberg 1, D-14476 Potsdam, Germany  \nThomas P. Sotiriou §  \nNottingham Centre of Gravity and School of Mathematical Sciences and School of Physics and Astronomy, University of Nottingham, University Park, Nottingham, NG7 2RD, United Kingdom  \nStoytcho S. Yazadjiev ¶  \nTheoretical Astrophysics,Eberhard Karls University of Tubingen, 72076 Tubingen, Germany Department of Theoretical Physics, Faculty of Physics, Sofia University, 1164 Sofia, Bulgaria Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Academy Georgi Bonchev Street 8, 1113 Sofia, Bulgaria  \nScalarization is a mechanism that endows strongly self-gravitating bodies, such as neutron stars and black holes, with a scalar-field configuration. It resembles a phase transition in that the scalar configuration appears only when a certain quantity that characterizes the compact object, for example, its compactness or spin, is beyond a threshold. A critical and comprehensive review of scalarization, including the mechanism itself, theories that exhibit it, its manifestation in neutron stars, black holes and their binaries, potential extension to other fields, and a thorough discussion of future perspectives, is provided.', 'I. Introduction': '2  \nII.  \n- VI. Open problems and future perspectives  \n41  \n- A. Scalarization and cosmology  \n42  \n- B. Dynamical evolution  \n42  \n- C. Model building in and beyond scalarization  \n43  \n- D. Observational prospects  \n44  \nAcknowledgments  \n44  \nReferences  \n45', 'I. INTRODUCTION': "Exploring the nature of gravity in the strong curvature regime has seen a recent surge of interest. This is expected to intensify, as it is driven by current and future observations of compact objects: black holes (BHs) and neutron stars (NSs). In particular, gravitational waves (GWs) produced by coalescing compact binaries have by now been routinely detected by the LIGO-VirgoKagra (LVK) Collaboration (Abbott et al. , 2019, 2021a,b). These observations enabled us to probe the highly dynamical and strong-field regime of general relativity (GR) for the first time. They have enabled one to perform new tests of GR and to constrain modifications thereof in a hitherto unexplored regime. Future spaceborne and ground-based GW observatories have testing GR and the standard model (SM) among their key priorities (Arun et al. , 2022; Barausse et al. , 2020; Kalogera et al. , 2021; Sathyaprakash et al. , 2019). At the same time, there is a new suite of electromagnetic observations that probe NSs with unprecedented sensitivity and timing resolution (Arzoumanian et al. , 2014; Gendreau and Arzoumanian, 2017; Gendreau et al. , 2012). On other fronts, the precision timing of binary pulsars has continually improved (Kramer et al. (2021)), measurements of the motion of stars at the Galactic Center are becoming more precise (Abuter et al. , 2018, 2020; Do et al. , 2019), and we have witnessed breakthroughs in supermassive BH imaging (Akiyama et al. , 2019).  \nAn interesting prospect is that these observations may reveal the existence of some otherwise elusive new fundamental fields, which could be an ingredient of new physics beyond the SM or beyond GR; see Clifton et al. (2012), Yunes and Siemens (2013), Berti et al. (2015), and Barack et al. (2019). For such fields to have remained undetected, there has to exist a mechanism to suppress them when gravity is weak. For scalar fields, which are ubiquitous in extensions of the SM and of GR, a possible realization of such a mechanism was first proposed by Damour and Esposito-Far'ese (1993) and named spontaneous scalarization . They showed that a specific type of nonminimal coupling between scalar field and gravity (or matter, after a field redefinition) leads to a theory that is indistinguishable from GR in weak-field gravitational experiments and yet predicts order unity deviations  \nfrom general-relativistic expectations in the strong-gravity regime of NSs.  \nAs today, the first model of spontaneous scalarization came about at a time in which gravitational experiments were producing new data from a then unexplored regime of gravity: the slow-velocity but strong-field regime of binary pulsars discovered by Hulse and Taylor (1975); see Damour (2015). This discovery inaugurated a new arena to test GR and its contenders (Damour and Taylor, 1992; Taylor and Weisberg, 1982; Taylor et al. , 1993). Meanwhile, slow velocities and weak-field Solar System tests had reached an accuracy that made it questionable whether viable theories that predict deviations from GR that are measurable with binary pulsars can exist (Will, 2014). Spontaneous scalarization settled this question and provided further motivation for the use of binary-pulsar observations to test GR. By now these observations have ruled out the Damour-Esposito-Far'ese (DEF) scalarization model (Antoniadis et al. , 2013; Kramer et al. , 2021; Zhao et al. , 2022).  \nIn recent years, however, spontaneous scalarization has received renewed interest. This is due to the realization that vacuum BH solutions of GR can also scalarize when the scalar field or fields couple suitably to the spacetime curvature (Doneva and Yazadjiev, 2018b; Silva et al. , 2018). This development also showed that the earlier DEF model is part of a much broader class of theories (Andreou et al. , 2019) that exhibit what resembles a phase transition in the strong field: once a quantity that describes a compact object, such as its compactness (Damour and Esposito-Far'ese, 1993; Doneva and Yazadjiev, 2018b; Silva et al. , 2018) or spin (Dima et al. , 2020), exceeds a certain threshold, the scalar field switches from a trivial constant configuration to a nontrivial one, and large deviations from GR appear. Conversely, one can think of this as deviations from GR getting severely 'screened' as soon as one crosses the same threshold in the opposite direction. It is this phase transition behavior that distinguishes scalarization from other models in which deviation from GR is induced and controlled by coupling to curvature, such as more general scalar-tensor theories with linear (Sotiriou and Zhou, 2014a,b; Yunes and Stein, 2011) or exponential (Kanti et al. , 1996) couplings to the GaussBonnet invariant.  \nSince the advent of GW astronomy, the broader class of theories that exhibit spontaneous scalarization have played a role similar to the one that the DEF model has played for binary-pulsar observations. These theories provide a putative explanation of why we have not detected new fundamental fields with existing observations, but we might still uncover them with high precision observations of astrophysical systems with specific characteristics.  \nThe aim of this review is to summarize, in a unified manner, the current status of this field. In Sec. II, we start by providing the theoretical background of the scalarization mechanism and its various subcases, following a  \npedagogical, rather than historical, approach. Next we discuss in more detail the literature on scalarization, first of NSs in Sec. III and then of BHs in Sec. IV. In so doing, we present the state of the art of our understanding of the consequences of scalarization for various situations of observational interest. In Sec. V, we discuss attempts to generalize scalarization to other field types. In Sec. VI, we outline open issues and summarize future perspectives in the field. Unless stated otherwise, we use geometrical units G = 1 = c and employ the mostly plus metric signature convention.", '1. Tachyonic instability and nonlinear quenching': "Before we discuss spontaneous scalarization in the context of gravity, it is instructive to review the dynamics of a real scalar field φ with a quartic self-interaction in Minkowski spacetime. The Lagrangian for this field is  \nL = 1 2 η µν ∂ µ φ∂ µ φ + V ( φ ) , (1)  \nwhere η µν is the Minkowski metric  \nV ( φ ) = 1 2 µ 2 φ 2 + 1 4 λφ 4 , (2)  \nµ is the bare mass, and λ is a coupling constant. The scalar then satisfies the following field equation:  \n□ η φ -µ 2 φ -λφ 3 = 0 , (3)  \nwhere □ η = η µν ∂ µ ∂ ν is the flat-spacetime d'Alembertian. φ = 0 is a solution of this equation. Consider now small perturbations δφ around φ = 0. By linearizing Eq. (3), we find that δφ obeys  \n□ η δφ -µ 2 δφ = 0 . (4)  \nThe corresponding dispersion relation is ω 2 = k 2 + µ 2 , where ω is the frequency and k is the wavenumber. If µ 2 > 0, the relevant solutions to this equation are plane waves and the perturbations decay. If instead µ 2 < 0, one encounters a tachyonic instability and the perturbations with small wave number exhibit exponential growth.  \nThis exponential growth seems catastrophic at first sight, but it does not have to be. As φ grows, the previously used linear approximation will quickly become invalid, and the nonlinear self-interaction λφ 3 will become important. It will be this interaction that will determine the end point of the instability. Assume that λ > 0 (and µ 2 < 0), in which case the potential has the well-known 'Mexican hat' shape. Equation (3) will then admit a second solution with constant φ , which we denote as φ min , as  \nit will correspond to the minimum of the potential. Equation (3) implies that φ 2 min = -µ 2 /λ . Thus, the tachyonic instability simply drives the scalar field away from the unstable local maximum of the potential and toward a stable minimum. This is sometimes referred to as tachyon condensation and is associated to a phase transition of the system.  \nThe key message from this simple example is that linearized perturbations around the unstable maximum capture the onset of the tachyonic instability but they are oblivious to the shape of the rest of the potential and hence cannot determine the end point. Nonlinear interactions, represented in this specific case by the φ 4 term in the potential (2), eventually quench the instability and drive the field to a different, stable configuration.", '2. Tachyonic instability in curved spacetime': 'In Sec. II.A.1 we considered a scalar0field that exhibited a tachyonic instability in flat spacetime. The generalization to curved spacetime is simple. If we promote the Minkowski metric to some general curved background described by a metric g µν , Eq. (4) becomes  \n□ δφ -µ 2 δφ = 0 , (5)  \nwhere □ = g µν ∇ µ ∇ ν , with ∇ µ the covariant derivative. The key difference here is that in curved space µ 2 < 0 is no longer sufficient for having a tachyonic instability.  \nTo see this, we take g µν to be the Schwarzschild metric, which can describe either a nonrotating BH or the exterior spacetime of a nonrotating NS in GR. The line element is  \nd s 2 = -( 1 -2 M r ) d t 2 + ( 1 -2 M r ) -1 d r 2 + r 2 dΩ 2 , (6)  \nwhere dΩ 2 = d θ 2 + sin 2 θ d ϕ 2 and M is the mass of the compact object. Because the spacetime is static and spherically symmetric, we can decompose the scalar perturbation δφ into spherical harmonics Y ℓm ( θ, ϕ ) and assume a harmonic time dependence  \nδφ = ∑ ℓm ψ ℓm ( r ) r Y ℓm ( θ, ϕ ) e -iωt , (7)  \nand by substitution into Eq. (5) we obtain a Schrodingerlike equation  \nd 2 ψ ℓm d r 2 ∗ + [ ω 2 -V eff ( r ) ] ψ ℓm = 0 , (8)  \nwhere we introduced the tortoise coordinate r ∗ defined as d r/ d r ∗ = 1 -2 M/r and V eff is an effective potential given by  \nV eff = ( 1 -2 M r )[ ℓ ( ℓ +1) r 2 + 2 M r 3 + µ 2 ] , (9)  \nwhich encodes information about the background, curved spacetime. To have an instability in a Schwarzschild BH spacetime it is sufficient but not necessary that  \n∫ ∞ -∞ d r ∗ V eff ( r ) ⩽ 0 , (10)  \n̸  \nwhere r ∗ = -∞ corresponds to the horizon radius in tortoise coordinates. For M = 0 (flat spacetime) this condition is always satisfied when µ 2 < 0 and Eq. (8) yields the same dispersion relation that we discussed in Sec. II.A. The situation is different for M = 0. Although it is not immediately obvious from upon inspection of Eq. (9), it turns out that µ 2 would have to be sufficiently negative for the tachyonic instability to occur. The main lesson here is that, in curved spacetime, the threshold for the tachyonic instability to happen depends on the spacetime; we will return to how one can determine it later. Note that, although we previously used Schwarzschild spacetime as an example, one can rederive Eq. (8) for a general static, spherically symmetric background provided that r ∗ is chosen appropriately.', '3. Scalarization and gravity': "̸  \nWe have thus far considered a scalar field with a negative bare mass squared, which is not well motivated. However, fields can acquire an effective mass squared µ 2 eff in specific situations due to their coupling to other fields. As an example, consider a scalar field that is nonminimally coupled to gravity and a term φ 2 R is present in the action, where R is the Ricci scalar. We then expect a contribution proportional to φR to the scalar's field equation; hence, the Ricci scalar contributes to the field's effective mass, that is, µ 2 eff ∝ µ 2 + R . Assume that the scalar field has no bare mass (i.e., µ = 0) and that the coupling to R is the only contribution to its effective mass µ eff . In flat spacetime, scalar-field perturbations would then be massless, whereas in curved spacetimes (with R = 0) they would be massive and in general also position dependent. Moreover, the sign of µ 2 eff would be controlled by the sign of R in this case. Hence, in some situations it would be possible for µ 2 eff to become sufficiently negative in some spacetime region and cause the scalar field to become tachyonically unstable despite this being impossible in flat spacetime.  \nJust as in the flat-spacetime example of the scalar with negative µ 2 and quartic interactions in Sec. II.A.1, this tachyonic instability does not have to be catastrophic. It can simply signal that the scalar needs to transition to a different configuration once curvature exceeds some threshold. The instability implies that the scalar field will grow, nonlinearities will become important, and, if they can quench the instability, then one can end up with a stable, different configuration, for both the scalar field and the spacetime.  \nThis is precisely the idea behind spontaneous scalarization, 1 first proposed in Damour and Esposito-Far'ese (1993). In a given generalized scalar-tensor theory, a configuration with a constant scalar field and a metric that solves Einstein's equations describes all gravitating systems except some that exhibit strong gravity. In the latter case, curvature becomes significant enough to render the constant scalar configurations tachyonically unstable. The tachyonic instability is eventually quenched by nonlinear effects, and there is a stable configuration with a nontrivial scalar and a spacetime that is no longer a solution to Einstein's equations.", '4. Strong-field phase transitions and weak-field screening': 'We have not yet shown that the mechanism of spontaneous scalarization, as described heuristically earlier, can be at play within a consistent gravity theory: we do so in Sec. II.B. Nonetheless, assuming that the proposal can be successfully implemented in some model, the following key observations can already be made:  \n- · Scalarization is a sharp transition to a new configuration that can differ significantly from the GR configuration for the same object, even when one is very near the threshold of the tachyonic instability. This is intuitive when one thinks of scalarization as a linear instability quenched by nonlinearity: even for a mild instability (large timescale) the scalar has to grow significantly for nonlinear effects to become important and manage to stop further growth.\n- · The onset of instability can be controlled by curvature couplings. In Sec. II.A.3, we considered as an example a coupling between the scalar field and the Ricci scalar R , but one can envisage couplings with other curvature invariants, as we see in in Sec. II.B. Hence, there can be models in which scalarization will occur only in the strong-field regime (where curvature can become large), while objects that exhibit weak gravity will show no deviation from GR (because the curvature is small). Combined with the previous point, this suggests that spontaneous scalarization can be thought of as a strong-field phase transition, whereby a field that is dormant in the weak field transitions to a nontrivial configuration in the strong field. Alternatively, one can think of scalarization in the reverse way: as a screening mechanism that forces a scalar field to transition to a trivial configuration in the weak field and hence explain why this field has managed to remain undetected so far.  \n- · The previous argument is based on the rather naive expectation that curvature invariants are a good measure of how strong the gravitational interaction is. As an example of the failure of this expectation, recall that for a Schwarzschild BH the Ricci scalar is zero; however, other curvature scalars, such as the Kretschmann scalar R µνρσ R µνρσ , are nonzero. In general, curvature invariants can have a complicated dependence on the properties of compact objects. Hence, more work is needed to understand what controls the threshold of the tachyonic instability and the onset of scalarization. This is addressed in Sec. II.C.', '1. Tachyonic instability and the minimal action': "We now turn our attention to the gravity theories that can exhibit spontaneous scalarization. As discussed, at the perturbative level the hallmark of spontaneous scalarization is a tachyonic instability. This begs the following question: Can we construct a minimal gravity theory which can have scalar-field perturbations that are tachyonically unstable? To do so, consider a gravity theory with a metric g µν and a scalar field φ . Assume that the theory is such that the following hold:  \n- A.1 Spacetimes that are solutions of Einstein's equations, potentially with a cosmological constant, and a constant scalar field are admissible solutions of this theory as well.\n- A.2 All terms in the action are at least quadratic in φ .\n- A.3 The field equations are second-order partial differential equations.  \nUnder these requirements, the equation governing the dynamics of scalar perturbations δφ on GR spacetimes can be cast in the form  \ng µν eff ∇ (0) µ ∇ (0) ν δφ -µ 2 eff δφ +NLC = 0 , (11)  \nwhere g µν eff , ∇ (0) and µ 2 eff are all computed in the background spacetime g (0) µν , and NLC denotes nonlinear corrections. In Eq. 11 g eff µν is an effective metric that can differ from g (0) µν for certain types of nonminimal couplings between the metric and the scalar field and µ 2 eff may contain not only a bare mass term but also other contributions.  \nIf one neglects nonlinearities and assumes that g eff is nondegenerate and has a Lorentzian signature, Eq. (11) becomes a curved-spacetime version of Eq. (4). This means that one can identify all theories with a single scalar field that are expected to lead to spontaneous scalarization by considering which couplings between a scalar and the metric can contribute to g eff and µ 2 eff while  \nstill satisfying the aforementioned assumptions A.1 -A.3 . The benefit of taking into account all possible such terms is that it would allow one to fully explore the onset of scalarization and identify a class of gravity theories that result in a scalarized spacetime.  \nAssumption A.2 appears to be essential to avoid having a source term in Eq. (11). We return to this shortly and show that it is a redundant assumption. But we first consider assumption A.3 . It ensures that there are no unwanted degrees of freedom, as would generically be the case if the equations contained higher-order derivatives (exceptions can exist, most notably in cases where field redefinitions can reduce the differential order of the equations). This assumption does limit the possibilities of the terms that one can consider. For example, a coupling term of the type φ 2 R µνλσ R µνλσ in the action would contribute to µ 2 eff but leads to higher-order field equations (in the absence of suitable counterterms). To deal with this potential pitfall, one can follow the lines of Andreou et al. (2019) and start with the Horndeski action (Deffayet et al. , 2009; Horndeski, 1974), also known as generalized scalar-tensor theory, 2  \nS = 1 16 πG 5 ∑ i =2 ∫ d 4 x √ -g L i + S m [Ψ m ; g µν ] (12)  \nand  \nL 2 = G 2 ( φ, X ) , (13a)  \nL 4 = G 4 ( φ, X ) R + G 4 X [ ( □ φ ) 2 -( φ µν ) 2 ] , (13c)  \nL 3 = -G 3 ( φ, X ) □ φ, (13b)  \nL 5 = G 5 ( φ, X ) G µν φ µν -1 6 G 5 X [ ( □ φ ) 3 -3 □ φ ( φ µν ) 2 +2( φ µν ) 3 ] , (13d)  \nwhere G 2 , G 3 , G 4 , and G 5 are arbitrary functions of the scalar field φ and its kinetic term X = -∇ µ φ ∇ µ φ/ 2. In Eq. (13) G iX = ∂G i /∂X ( i = 4 and 5), G µν is the Einstein tensor, and the notation φ µν = ∇ µ ∇ ν φ was introduced so, for example, ( φ µν ) 2 = φ µν φ µν = ∇ µ ∇ ν φ ∇ µ ∇ ν φ . Finally, S m is the matter action, with matter fields collectively denoted by Ψ m . This is the most general action for a metric and a scalar field that leads to second-order field equations in four dimensions upon direct variation; see Kobayashi (2019) for a review. We assume for the moment that matter couples minimally to the metric only. This means that the choice of fields g µν and φ correspond to the so-called Jordan frame; we return to this issue later.  \nImposing assumptions A.1 and A.2 on the action (12) places restrictions on the G i functions, as we later see.  \nWe refer to Andreou et al. (2019) for a detailed discussion. For our purposes, it is sufficient to say that, by perturbing around an arbitrary spacetime that is assumed to be a solution of Einstein's equations with a constant scalar field, we can identify all of the terms that contribute, at the linear level, to g µν eff and µ 2 eff , as defined in Eq. (11). These terms amount to the following action:  \nS min = 1 16 πG ∫ d 4 x √ -g [ R -1 2 ( γ 1 + γ 2 R ) ∇ µ φ ∇ µ φ + γ 2 R µν ∇ µ φ ∇ ν φ -1 2 µ 2 φ φ 2 -1 4 βφ 2 R + 1 2 αφ 2 G -2Λ] + S m [Ψ m ; g µν ] , (14)  \nwhere G is the Gauss-Bonnet invariant, which is defined in terms of the Riemann tensor and its familiar contractions as  \nG = R µνρσ R µνρσ -4 R µν R µν + R 2 , (15)  \nand where α , β , γ i , and µ 2 φ can be expressed in terms of the G i functions and their derivatives evaluated in the background configuration (Andreou et al. , 2019). 3 We refer to this action, in a slight abuse of terminology, as the minimal action for scalarization, in the sense that it contains all of the terms that contribute to the onset of scalarization manifesting as a tachyonic instability. As such, it can be used to study and understand what triggers scalarization and to determine the relevant instability thresholds.  \nBefore we go further, we examine what happens if we decide to drop assumption A.2 altogether, but still impose assumptions A.1 and A.3 . Working along the same lines as before, one arrives at a different set of terms composing the action  \nS ' min = 1 16 πG ∫ d 4 x √ -g [ R -1 2 γ ' 1 + γ ' 2 R φ ∇ µ φ ∇ µ φ + γ ' 2 φ R µν ∇ µ φ ∇ ν φ + τφ + η φR + λφ G -2Λ] + S m [Ψ m ; g µν ] . (16)  \nAs before, τ , η , λ , and γ ' i can be expressed in terms of the G i functions and their derivatives evaluated in the background configuration. It might seem counterintuitive that an action containing terms linear in φ leads to perturbation equations that are in the form of Eq. (11), which has no source terms. This is due to the presence of terms that are nonanalytic in φ .  \nThe first impression may be that abandoning assumption A.2 has given rise to a second minimal action for scalarization. However, action (16) is just a field redefinition away from action (14). Indeed, one can start with Eq. (16), introduce the redefinition 4 φ → φ 2 , and obtain action (14), with the following mapping of parameters: γ 1 = 4 γ ' 1 , µ 2 φ = -4 τ , β = -η , and α = λ (Andreou et al. , 2019). This equivalence demonstrates (i) that assumption A.2 is redundant and (ii) that up to field redefinitons Eq. (14) is sufficient to capture all terms that contribute to the onset of scalarization and satisfy assumptions A.1 and A.3 .  \nOne can see by inspection that the γ i terms in Eq. (14) will contribute to g µν eff , while the rest of the terms will contribute to µ 2 eff . Hence, if the latter vanish, the former cannot trigger scalarization by themselves, as the effective mass would vanish. Nevertheless, the γ i terms will affect the threshold of the tachyonic instability we associate with scalarization; cf. the discussion about the tachyonic instability in curved spacetime in Sec. II.A.2. Additionally, µ 2 φ corresponds to the bare mass of the scalar field, so it is expected to be positive. We then conclude that the terms that are expected to trigger scalarization in the strong-field regime are only the couplings of φ to R and G . In fact, we see shortly that these are indeed the terms present in the known models of scalarization.  \nTo summarize, the minimal action (14) can be used to study the onset of spontaneous scalarization triggered by a nonminimal coupling to gravity. It contains all possible terms that contribute to the associated tachyonic instability at the linearized level, so it could be used to study the threshold and onset of this instability in full generality. As previously discussed, as the instability progresses it is expected to be quenched nonlinearly, and the end point will be a scalarized configuration. The terms in Eq. (14) can contribute nonlinearly as well, but one could add a plethora of other nonlinear interactions, ranging from scalar self-interactions, for example, φ 4 , to nonminimal coupling terms that do not contribute to linear perturbations around curved spacetime with constant scalar, for instane, φ 4 G . That is, one can start with Eq. (14), or even a subset of terms therein, and construct different scalarization models. Models that differ only by terms that are not in Eq. (14) will have the same behavior regard to the onset of scalarization, and hence the configurations that one expects to not scalarize, but they can differ in the properties of scalarized solutions (Andreou et al. , 2019; Macedo et al. , 2019; Minamitsuji and Ikeda, 2019b; Silva et al. , 2019). This is further discussed alogn with specific examples in Secs. IV.A.1 and IV.A.2.  \nBefore proceeding to discuss more specific known models, we return to the issue of the coupling to matter. We have thus far assumed that matter couples minimally to the metric only. This assumption is sufficient to ensure that the theory is compatible with the weak equivalence principle (WEP) (Will, 2018). To satisfy the WEP it is sufficient to have matter couple minimally to some metric, but this does not need to be the same metric (or choice of other fields) for which the theory has second-order field equations. However, it is known that a disformal transformation (Bekenstein, 1993) of the form  \ng µν → C ( φ ) [ g µν + D ( φ ) ∇ µ φ ∇ ν φ ] (17)  \nleaves the Horndeski action (12) formally invariant (Bettoni and Liberati, 2013; Zumalac'arregui and Garc'ıaBellido, 2014). It was shown by Andreou et al. (2019) that coupling matter minimally to a metric that is related to g µν by such a disformal transformation as done, Minamitsuji and Silva (2016), amounts to a redefinition of γ 2 in the linearized theory around spacetimes that are solutions of Einstein's equations. Hence, such a coupling would be redundant when the onset of scalarization is studied using the minimal action (14). One could also entertain the idea of coupling matter to some composite metric ¯ g µν that depends on both g µν and φ in a different manner than in Eq. (17). In such a case, it is likely that assumption A.3 would be violated and one would have to start the analysis presented here with a generalization of the action (14).", "2. Damour-Esposito-Far'ese Model": "As mentioned, the concept of scalarization as we described it was first discussed by Damour and EspositoFar'ese (1993). They considered the theory  \nS = 1 16 πG ∗ ∫ d 4 x √ -g [ R -2 ∇ µ φ ∇ µ φ ] + S m [Ψ m ; A 2 ( φ ) g µν ] . (18)  \nEquation (18) is said to be written in the Einstein frame, which means that, contrary to our previous assumptions and conventions, the scalar field is coupled minimally to gravity and has a canonical kinetic term. The coupling with matter field Ψ m is through the function A 2 ( φ ). In Eq. (18) G ∗ carries a subscript, as it is not generally equal to G used thus far. Variation of the action (18) with respect to φ yields the field equation  \n□ φ = -4 πG ∗ α ( φ ) T, (19)  \nwhere  \nα ( φ ) = d ln A ( φ ) / d φ (20)  \nand T = g µν T µν is the trace of the matter energymomentum tensor in the Einstein frame defined as T µν =  \n2( -g ) -1 / 2 δS m /δg µν . We see that α ( φ ) controls the coupling strength between the scalar field and matter.  \n̸  \nIf α ( φ 0 ) = 0 for some constant scalar-field value φ 0 , the constant scalar configuration with T = 0 will be an admissible solution of the theory. It then follows from the generalized Einstein's equations  \nR µν = 2 ∇ µ φ ∇ ν φ +8 πG ∗ ( T µν -1 2 g µν T ) (21)  \nthat these will be solutions of GR since the first term on the right-hand side vanishes.  \nAt the same time, if we perturb Eq. (19) linearly in φ in a fixed background metric that is a solution of GR and compare with Eq. (11) we find that β ( φ 0 ) = (d α/ d φ ) φ = φ 0 and T determine the value and sign of the effective mass square of the perturbations, namely,  \nµ 2 eff = -4 πG ∗ β ( φ 0 ) T. (22)  \nFor stars, one generally has T < 0. Hence, for a negative sign of β ( φ 0 ) and the right magnitude of both quantities, the scalar can develop a tachyonic instability around a spacetime that describes stars in GR, as previously discussed and as studied in detail by Harada (1997). 5 It was shown by Novak (1998a) that this instability is quenched by nonlinearities and that the outcome is a NS with a nontrivial scalar-field configuration. These scalarized NSs were shown by Damour and EspositoFar'ese (1993) to have properties, such as their masses M and radii R , that can be dramatically different from their GR counterparts.  \nIn much of the literature considering scalarization, the function α ( φ ) is taken to have the form 6  \nα = α 0 + β 0 φ ≡ α DEF , (23)  \nwhere α 0 and β 0 are dimensionless constants. Sometimes this choice, rather that the more general action of Eq. (18), is referred to as the DEF model. The constant α 0 is then assumed to vanish to allow for constant φ solutions [cf. Eq. (19)], or it is assumed to be small. In the latter case, all stars will carry some nontrivial scalar field, but by tuning down α 0 any deviation from GR would be undetectable until scalarization kicks in. In its original formulation, the DEF model did not include a bare mass or self-interactions for the scalar, but a potential V ( φ ) can be added to the action and this option has been considered in the literature by Chen et al. (2015); Popchev (2015); and Ramazano˘glu and Pretorius (2016), as we see in Sec. III.  \nThus far it appears that the DEF model is not covered by our minimal action (14), because of our earlier assumption that the scalar does not couple to the matter. However, by defining A 2 ( φ ) g µν as a new metric and rewriting the action (18) in terms of that new metric, the scalar field is no longer coupled to matter. This is referred to as the Jordan frame. It was shown by Andreou et al. (2019) that, at linearized level and after a suitable scalar-field redefinition, the DEF model is equivalent to the action  \nS = 1 16 πG ∫ d 4 x √ -g [( 1 -1 4 β 0 φ 2 ) R -1 2 ∇ µ φ ∇ µ φ ] + S m [Ψ m ; g µν ] . (24)  \nEquation (24) is indeed a particular case of the action (14) in which γ 1 = 1, γ 2 = α = Λ = 0, and β = β 0 . Hence, the DEF model, with what regards to the onset of the tachyonic instability that leads to scalarization, is captured by the minimal action (14) and corresponds to one of the two couplings to curvature that can trigger scalarization.  \nBefore moving on, we mention that the original formulation of the DEF model, which leads directly to Eq. (19), suggests that it is the coupling to matter that controls and triggers scalarization. Indeed, when T = 0, as is the case for BHs, Eq. (19) becomes □ φ = 0 and admits only constant φ solutions for stationary and asymptotically flat configurations by virtue of a no-hair theorem by Hawking (1972) [this remains true when one includes a potential; see Sotiriou and Faraoni (2012)]. However, our previous analysis and the correspondence between the DEF model and action (24) at the linearized level makes it clear that the DEF model is part of a broader class of theories in which scalarization is present and controlled by the couplings to curvature, rather than matter, and this observation has been crucial for the development of models that exhibit BH scalarization. It is the fact that R and T are related through the trace of the theory's generalized Einstein equation [cf. Eq. (21)] that allows for both interpretations in the DEF model.", '3. Scalar-Gauss-Bonnet gravity': "It was first shown by Doneva and Yazadjiev (2018b) and Silva et al. (2018) that theories described by the scalar-Gauss-Bonnet (sGB) action  \nS = 1 16 πG ∫ d 4 x √ -g [ R -1 2 ∇ µ φ ∇ µ φ + f ( φ ) G ] + S m [Ψ m ; g µν ] (25)  \ncan exhibit BH scalarization provided that (d f/ d φ ) φ = φ 0 = 0 for some constant φ 0 . This is an existence condition for constant φ configurations that are solutions of GR. As proven by Silva et al. (2018), BH  \nsolutions of GR are unique solutions to the theory (25) provided that (d 2 f/ d φ 2 ) φ = φ 0 G < 0. To understand this, we can proceed as follows. By varying the action with respect to φ , we find  \n□ φ + f ,φ ( φ ) G = 0 , f ,φ ( φ ) = d f/ d φ. (26)  \nOnce again, we can consider linear perturbation of φ on a fixed background and compare with Eq. (11). We find that (d 2 f/ d φ 2 ) φ = φ 0 G plays the role of an effective mass square for the scalar perturbations,  \nµ 2 eff = -(d 2 f/ d φ 2 ) φ = φ 0 G . (27)  \nHence, violating the condition (d 2 f/ d φ 2 ) φ = φ 0 G < 0 is necessary, but not sufficient, to develop a tachyonic instability that can lead to scalarization.  \nDoneva and Yazadjiev (2018b) chose f ( φ ) to be f ( φ ) = λ 2 (1 -e -3 / 2 φ 2 ) / 12, whereas Silva et al. (2018) focused on f ( φ ) = αφ 2 / 2. Note that, in the linearized theory around φ = φ 0 and provided that the condition (d f/ d φ ) φ = φ 0 = 0 is satisfied, any choice of f ( φ ) is equivalent to f ( φ ) = αφ 2 / 2. Hence, all scalarization models described by action (25) are captured by the minimal action (14) and, in particular, by the coupling between the scalar and the Gauss-Bonnet invariant in regard to the onset of scalarization. However, different choices of f ( φ ) will exhibit different behavior in the nonlinear regime and hence scalarized BHs will generally have different properties.  \nIndeed, it was shown by Bl'azquez-Salcedo et al. (2018) that the static, spherically symmetric scalarized BHs that were found by Silva et al. (2018) for f ( φ ) = αϕ 2 / 2 are unstable against radial perturbations, unlike their for counterparts for f ( φ ) = λ 2 (1 -e -3 / 2 φ 2 ) / 12 found by Doneva and Yazadjiev (2018b). Silva et al. (2019) later demonstrated by examining the case f ( φ ) = αφ 2 / 2 + ξφ 4 that it is indeed the nonlinearity in φ that controls the stability of scalarized BHs. It was further shown by Macedo et al. (2019) and Minamitsuji and Ikeda (2019b) that a quartic scalar-field self-interaction would be sufficient to make scalarized BHs stable for the f ( φ ) = αφ 2 / 2 case. These results, which are discussed in Sec. IV, are a clear demonstration that although the onset of scalarization can be described fully using action (14), the end point of the tachyonic instability and the properties of the scalarized configurations will depend to the nonlinear interaction between the scalar and curvature and are thus model dependent.  \nWe remark that models described by action (25) lead to scalarization of compact stars as well for certain regions of their parameter spaces; see Doneva and Yazadjiev (2018a) and Silva et al. (2018).  \nAs we have seen, the DEF model and the sGB models of scalarization correspond respectively, to the φ 2 R and φ 2 G terms (in addition to R and the canonical kinetic term for the scalar field) in the minimal action (14). We have also argued that these two terms are the only terms that can trigger scalarization as a tachyonic instability around a spacetime that is a solution of GR. These facts together suggest considering the following action (Antoniou et al. , 2021a,b):  \nS = 1 16 πG ∫ d 4 x √ -g [ R -1 2 ∇ µ φ ∇ µ φ -1 4 βφ 2 R + 1 2 αφ 2 G ] + S m [Ψ m ; g µν ] . (28)  \nThis theory is interesting from the perspective of an effective field theory (EFT). The previously considered terms can be seen as part of an EFT in which the scalar field enjoys reflection symmetry (i.e., invariance under φ → -φ ), while shift symmetry (i.e., invariance under φ → φ +const . ) is broken only by the coupling to the curvature scalars. This theory is not a complete EFT, as there are other operators compatible with these symmetries, such as φ 4 R and G µν ∇ µ φ ∇ ν φ . Nonetheless, the theory is phenomenologically interesting for various reasons.  \nTo begin, it has GR with a constant scalar field as a late-time cosmic attractor for β > 0 (Antoniou et al. , 2021a). To appreciate why this is important, recall that one can think of scalarization in terms of a tachyonic instability of compact object configurations that are solutions of GR with φ = φ 0 . Below the threshold of this instability, these configurations are expected to be stable and exhibit no deviation from GR. The attractive feature of scalarization is that weakly gravitating systems will belong in this category and hence scalarization can be a form of weak-field screening of the scalar field. However, this argument assumes that φ = φ 0 everywhere in the Universe and deviates from this value only due to scalarization. If there were another reason for φ to evolve away from φ 0 , this would make even weakly gravitating objects develop a nontrivial scalar configuration. Cosmic evolution can indeed cause such evolution, as shown by Damour and Nordtvedt (1993), to the extent that weakly gravitating systems would become sufficiently scalarized to make the DEF model (Anderson et al. (2016)) and sGB scalarization models (Anson et al. (2019a) and Franchini and Sotiriou (2020)) fail weak-field and cosmological tests of gravity without severely fine-tuning the initial conditions for cosmic evolution.  \nThe theory described by action (28) provides an elegant solution to this problem (Antoniou et al. , 2021a). As pointed out by Damour and Nordtvedt (1993), the DEF model has GR as a cosmic attractor for β 0 > 0, whereas scalarization requires β 0 to be sufficiently negative. It was  \nalready mentioned that, at the linearized level around a GR background, the DEF model is equivalent to action (24), so it is reasonable to expect that the φ 2 R term in action (28) would tend to drive the scalar field to a constant in late-time cosmology. Moreover, this term should be dominant over the φ 2 G term at low curvatures, while the latter should dominate at high curvatures and trigger scalarization. It was indeed shown by Antoniou et al. (2021a) that cosmic evolution in the model of action (28) tracks the cosmic evolution of GR from radiation domination onward, and that φ is driven to φ 0 rapidly during matter domination.  \nIt was also shown by Antoniou et al. (2021b) and Ventagli et al. (2021) that for β > 0 one can have a range of values for α in which BHs scalarize but NSs do not. This is interesting because the strongest constraints on scalarization thus far, and specifically the DEF model, are based on binary pulsars (Kramer et al. , 2021; Zhao et al. , 2022). These constraints can be evaded if scalarization is limited to BHs. Antoniou et al. (2021b) also provided strong indication that scalarized BHs should be radially stable for β > 0, which is not the case for β = 0. Both of these issues are discussed in Sec. IV.  \nIn summary, adding a βφ 2 R term with β > 0 to the simplest sGB scalarization model addresses a series of concerns. This term would be there anyway in an EFT, as it has lower mass dimensions than φ 2 G . These considerations should act as a reminder that scalarization theories are currently still toy models in need of a completion.", '5. Tensor-multiscalar theories': "Thus we have considered theories with a single field, but one could also study scalarization in models with multiple scalar fields. These tensor-multiscalar theories were studied by Damour and Esposito-Far'ese (1992) and in their simplest form they are described by the action  \nS = 1 16 πG ∗ ∫ d 4 x √ -g [ R -2 γ ab ( φ ) ∇ µ φ a ∇ µ φ b -4 V ( φ )] + S m [Ψ m ; A 2 ( φ ) g µν ] , (29)  \nwhere φ denotes a multiplet of N scalar fields { φ 1 , . . . , φ N } and γ ab ( φ ) is a target-space metric that mixes their kinetic terms. Equation 29 can be seen as a generalization of action (18) to multiple fields. We return to this theory in Sec. III, as it was studied in the context of scalarized NSs by Horbatsch et al. (2015) and Doneva and Yazadjiev (2020a). However, we emphasize that one could consider further generalizations that involve, for example, couplings to the Gauss-Bonnet invariant or further derivative interactions between the scalar fields.", 'C. Types of scalarization': 'One can classify types of scalarization based on which property of the compact object triggers scalarization and controls the threshold of the tachyonic instability. As we saw in our discussion of the model (14), it is the couplings between the scalar and curvature invariants that control the onset of scalarization. Therefore, if one thinks of the onset of scalarization as an instability around a spacetime of GR, what controls the onset reduces to how curvature invariants depend on the properties of the object that curves spacetime.', '1. Induced by compactness': 'For static, spherically symmetric BHs, there is a straightforward answer: in GR the exterior is described by the Schwarzschild spacetime, so R = 0 and G = 48 M 2 /r 6 , where M is the mass of the BH and r is the areal radial coordinate. Hence, for a given M , G scales monotonically with r and whether the scalar will develop a negative enough effective mass square outside the horizon is controlled by M , with lighter BHs being susceptible to the scalarization instability. The opposite trend is present in NS scalarization in the DEF model: the trace of the energy-momentum tensor of the star is what controls the onset, and for most equations of state (EOS) its magnitude shows a steady increase as the radius decreases. The tachyonic instability is stronger at the center of the star, where the density is higher, and this in turn correlates with the total mass of the star, with heavier stars scalarizing.  \nTo unify these two pictures, one can observe that BHs get more compact, in the sense of average mean density, as their mass decreases, whereas for NSs it is the other way around for most EOS. Hence, how compact an object is plays a key role in scalarization and can determine whether the object will be affected by the tachyonic instability that triggers scalarization.  \nAlthough we colloquially refer to compactness as one of the triggers of scalarization in the rest of the review, there are some caveats that we ought to mention and on which we further elaborate. First, compactness for BHs is often defined in the literature as the mass over the horizon radius. By that definition, unlike with the more colloquial meaning of the term, Schwarzschild BHs of any mass have the same compactness (but not the same curvature near the horizon). Second, for NSs the relation between compactness and scalarization is highly dependent on the EOS, and there are cases where the tachyonic instability is not strongest at the center of the star. Finally, as we next discuss, other properties of compact objects also affect the value of curvature invariants, and hence can also control scalarization.', '2. Induced by spin': 'Rapidly rotating NSs in the context of scalarization were first studied by Doneva et al. (2013). They showed that for DEF-like models rapid rotation can enhance scalarization by increasing the parameter space where scalarization can occur. Conversely, it is also known that larger spin leads to weaker scalarization for BHs in sGB gravity with d 2 f/ d ϕ 2 < 0 (Cunha et al. , 2019).  \nIn suitable theories, spin can in itself induce a tachyonic instability that triggers scalarization (Dima et al. , 2020), and spinning scalarized BHs in these theories were subsequently constructed explicitly by Herdeiro et al. (2021b) and Berti et al. (2021). The fact that spin can trigger BH scalarization (rather than just controlling the scalar charge) is important because it opens up the possibility that only rapidly rotating BHs carry scalar charge, irrespective of their mass.  \nThese results, as well as further work on scalarized rotating compact objects that is discussed in Secs. III and IV, show that spin, a ubiquitous property of astrophysical objects, can play an important role in the amount of scalar charge that a scalarized object carries.', '3. Induced by matter or coupling to other fields': 'Thus far we have linked scalarization to a tachyonic instability at the perturbative level (although it is fundamentally a nonperturbative effect), and we have focused on models in which that instability is controlled by the nonminimal coupling to gravity. However, as we saw, in the DEF model one can think of µ 2 eff as being controlled either by the trace of the stress-energy tensor of matter T [see Eq. (22)] or by the Ricci scalar R [see the action (24)] The latter interpretation has the advantage of providing a unified framework of scalarization of BHs and stars as linked to a nonminimal coupling to gravity, which is the perspective that we followed. However, the former interpretation highlights that µ 2 eff in Eq. (11) could instead be attributed to any type of coupling between φ and another matter field.  \nFor example, Stefanov et al. (2008) considered a scalar field coupled to nonlinear electrodynamics, while Herdeiro et al. (2018) focused on Einstein-Maxwell-scalar theory with the addition of the coupling e -λφ 2 F µν F µν , where F µν is the usual Faraday tensor. In both cases it was shown that electrically charged BHs can develop scalar hair through scalarization. Further work in this direction is summarized in Sec. IV.C.  \nWhen µ 2 eff in Eq. (11) is thought of as being introduced by a coupling to matter, one is also led to consider whether surrounding matter, such as an accretion disk, a companion, or the Galaxy, could scalarize a BH even in models where BHs cannot scalarize in vacuum. It was shown by Cardoso et al. (2013a,b) that this can indeed  \noccur in the DEF model.', '4. Dynamical scalarization': "Thus far we have discussed the scalarization of isolated compact objects, but what happens when they form a binary? As we alluded to, when we embed a NS in an ambient scalar-field environment [which is the case for the nonvanishing α 0 in Eq. (23)], stars will always carry some small scalar charge and they can still nonperturbatively develop large charge values when they scalarize. This is called induced scalarization (Salgado et al. , 1998), and binaries provide a natural scenario for it to occur. Imagine two NSs, each with its own compactness, such that one is scalarized and the other is not. As the system inspirals, at some point the nonscalarized NS will start experiencing the presence of the scalar field sourced by its companion, and induced scalarization will then take place (Barausse et al. , 2013).  \nAnother, perhaps more dramatic, scenario is that of dynamical scalarization . In this case, two nonscalarized NSs can become scalarized once their orbital separation becomes sufficiently small. Qualitatively this can be quantified by some measure of an 'effective compactness' of the binary that, only for an isolated NS, can trigger scalarization once it reaches a certain threshold (Palenzuela et al. , 2014; Shibata et al. , 2014; Taniguchi et al. , 2015). In a quasicircular binary this effective compactness only increases (it scales inversely with the orbital separation), but that does not have to be the case in an eccentric orbit. In such cases, the effective compactness oscillates in time, being largest when the NSs are closest. This leads to a transient dynamical scalarization of the system, whereby the two NSs continuously scalarize and descalarize as the system inspirals.  \nWhat about BHs? In sGB models the scalar field is sourced by the Gauss-Bonnet invariant, and therefore a binary composed of two scalarized BHs would in general result in an unscalarized BH remnant since the latter has a larger mass and therefore a smaller spacetime curvature. That is, the system descalarizes (Silva et al. , 2021b). However, depending on the initial nonscalarized BHs' spins and masses, one can have cases where the remnant scalarizes due to its large spin, i.e., there can be dynamical spin-induced scalarization (Elley et al. , 2022).  \nAs we have seen, compact binaries lead to new manifestations of scalarization. We further discuss these in Secs. III and IV.", '5. Beyond scalarization': 'Thus far we have discussed scalarizations as a linear, tachyonic instability for a scalar field that is then quenched nonlinearly and leads to a nontrivial scalar configuration.  \nThere are many ways to extend this paradigm and yet keep the key outcome: to have fields that undergo what resembles a phase transition - a (sharp change from a trivial to a nontrivial configuration) in the strong-field regime.  \nOne direction is to generalize the mechanism to different fields, such as vectors, tensors, and spinors, generating models of spontaneous vectorization , tensorization or spinorization (Ramazano˘glu, 2017, 2018a,b,c). Another approach would be to construct models in which the transition is triggered not by a tachyonic instability but by some other linear instability (Ramazano˘glu, 2018a). Both of these directions are discussed in Sec. V. A third direction comes from the possibility that the transition might not be triggered by a linear instability, but might instead be a fully nonlinear effect (Doneva and Yazadjiev, 2022). It was shown that if one chooses f ( φ ) = exp( βφ 4 ) in the action (25), one obtains a theory in which scalar perturbations are massless around Kerr BHs, and hence there cannot be any linear tachyonic instability. Yet, stable scalarized solutions still exist for certain masses and spins.', 'D. Quantum aspects and classical analogs': "While most of this review deals with scalarization from a classical field theory perspective, it is noteworthy to comment on the quantum aspects of this phenomenon. In particular, Lima et al. (2010) studied a quantum scalar field nonminimally coupled to gravity living in the classical background of a compact star spacetime (i.e., they worked in the semiclassical approximation). They showed that some field modes can go through an exponential growth causing the vacuum expectation value of the field operator ˆ φ 2 to grow, ultimately causing the vacuum expectation value of the energy-momentum tensor ˆ T µν of the field itself to grow in the same manner; see also Lima and Vanzella (2010). In this sense, one can say that classical curved spacetimes can 'awake the vacuum state' of quantum fields. But what is the connection with spontaneous scalarization? The model studied by Lima et al. (2010) is simply Eq. (24) for a massive scalar field: the Jordanframe equivalent of the DEF model for a massive scalar. In the parameter space spanned by the scalar-field-Ricciscalar coupling constant and stellar compactness relevant to our discussion, the regions where the instability occurs agree with the classical prediction to where scalarization should happen, as shown by Pani et al. (2011b). Thus, one can think of the classical 'scalar-field perturbations,' which we have so often spoken about to be seeded by quantum field fluctuations (Landulfo et al. , 2015). We refer the interested reader to Landulfo et al. (2012); Lima et al. (2013); Mendes et al. (2014a,b); and Santiago et al. (2016) for other works on the semiclassical approach.  \nFinally, while our review focuses on astrophysical implications of scalarization, we remark that the realization  \nof this effect in condensed matter systems has also been studied. More specifically Ribeiro and Vanzella (2020) devised a classical analog that exhibits this phenomenon based on the nonlinear optics of metamaterials, and this could in principle be observed experimentally.", 'III. NEUTRON STAR SCALARIZATION': "Scalarization of compact objects was first considered in the context of NSs by Damour and Esposito-Far'ese (1993). As discussed in Sec. II, in this case the nonzero trace of the energy-momentum tensor of the nuclear matter acts as a source of the scalar field and evades the no-hair theorems existing for vacuum BHs in certain scalar-tensor theories. Since then, the scalarized NSs in the DEF model have attracted significant attention, and they are perhaps the most studied compact objects beyond GR. Over the years, NS scalarization has also been examined in other scalar-tensor theories beyond the DEF model.  \nIn this section, we present the developments in the field over the past three decades. In Sec. III.A, we start with the DEF model and its generalizations, placing a special emphasis on the constraints coming from binarypulsar observations. In Sec. III.B, we proceed to the dynamics of such compact objects, both isolated ones and those in binaries, discussing their stability in addition to GW emission. In Sec. III.C, we discuss the various astrophysical implications of scalarized NSs. Finally, in Sec. III.D we review NSs scalarization beyond the DEF model.", "A. Equilibrium neutron stars in the Damour-Esposito-Far'ese model": 'In this section, we review the equilibrium properties of NSs in the DEF theory, both in its original form and in extensions of the theory. The strongest constraints on the former come from binary-pulsar constraints, which we also discuss here.', "1. The original Damour-Esposito-Far'ese model and binary-pulsar constraints": "a. Static neutron stars: As we saw in Sec. II.B.2, the effective mass squared of scalar-field perturbations in the DEF model is given by Eq. (22), namely,  \nµ 2 eff = -4 πG ∗ β ( φ 0 ) T,  \nwhere we recall that β ( φ 0 ) is the derivative of the scalarmatter coupling α ( φ ) evaluated at a constant background scalar-field value φ 0 . Thus, the condition for scalarization µ 2 eff < 0 can be satisfied when the trace of the energymomentum tensor T and β ( φ 0 ) have the same sign. For  \nrealistic NSs, which are modeled as a perfect fluid with pressure p and energy density ε , we normally have  \nT = 3 p -ε < 0 , (30)  \nand hence scalarization can happen when β ( φ 0 ) < 0. In particular, for the coupling function (23) studied by Damour and Esposito-Far'ese (1993), we have β ( φ 0 ) = β 0 being a constant. Thus, most studies focus on the case where both T and β 0 are negative, but we remark that at sufficiently high densities some EOSs predict a positive sign of T , resulting in scalarization when β 0 > 0; see e.g., Mendes (2015) and Podkowka et al. (2018). For most NS models, p and ε are related by a barotropic EOS ε = ε ( p ).  \nIn Fig. 1 we show some properties of nonrotating scalarized NSs for α 0 = 0 and some illustrative values of β 0 < 0. Here in Sec. III.A.1, we consider the case of a vanishing scalar-field potential as in the original DEF model. The influence of a nonzero potential is further discussed later. The EOS is taken to be that of Akmal, Pandharipande, and Ravenhall (APR) (Akmal et al. , 1998). In the left panel of Fig. 1, we show the Einstein-frame ArnowittDeser-Misner (ADM) mass of a sequence of NS solutions, parametrized by the energy density at the center of the star ε c .  \nWhen α 0 = 0, NS solutions of GR are also solutions of the DEF model; they are indicated by solid black lines. These solutions are characterized by having a zero scalar field. We see that when a specific critical central energydensity value is reached (say, ε c , 1 ), the GR sequence becomes unstable and a new branch of stable solutions with a nontrivial scalar field (i.e., scalarized stars) bifurcates from it. In our example, the value of ε c , 1 depends only on β 0 and on the EOS. The scalarized branch merges again with the GR branch at a second bifurcation point at a larger energy density ε c , 2 . Hence, scalarized NSs exist only in the range ε c ∈ [ ε c , 1 , ε c , 2 ]. We see that the larger | β 0 | becomes, the more dramatic the deviations in the NS mass relative to GR are. In addition, the range of ε c in which scalarization can happen increases. This is also shown in the middle panel of Fig. 1, where we plot the ADM mass as a function of the radius R , and in the right panel of Fig. 1, where we show the scalar charge D as a function of the ADM mass. The scalar charge is defined in terms of an expansion at spatial infinity of the scalar field φ , namely,  \nφ = φ 0 + D/r + O ( r -2 ) , (31)  \nwhere φ 0 is the cosmological background value of φ , often assumed to be zero for simplicity. The right panel of Fig. 1 shows how D has a small magnitude near the bifurcation point where scalarization kicks in, grows monotonically with M , and then approaches zero again onc e M approaches the mass corresponding to the second bifurcation point at ε c , 2 .  \nM  \n/  \nM  \n2.5  \nM  \n/  \nM  \nFigure 1 Some properties of NSs in the DEF model with α 0 = 0, β 0 < 0, and using the APR EOS. We show the mass M as a function of the central energy density ε c (left panel), the mass M as a function of the radius R (middle panel), and the scalar charge D as a function of M (right panel). The solid black lines correspond to the GR solutions, which are also solutions in the DEF model. We use differently colored lines to represent the scalarized branches for β 0 = -5 . 0 , -4 . 8, and -4 . 5. NSs past the maximum mass in the M -ε c plane are unstable to radial oscillations. These sequences of stars are shown as dotted lines.  \n<!-- image -->  \nFor which values of β 0 does scalarization occur? Recall that the effective potential V eff needs to be sufficiently negative in order to support at least one tachyonic mode; cf. Eqs. (9) and (10). For typical NS densities, Eq. (22) implies that scalarization exists for β 0 ≲ -4. Damour and Esposito-Far'ese (1993) made this estimate using a Newtonian approximation and confirmed it by integrating the fully relativistic equations of stellar equilibrium. Subsequent works refined the threshold for scalarization to β 0 ≲ -4 . 35 and also showed that this bound is not very sensitive to the EOS (Altaha Motahar et al. , 2017; Novak, 1998b; Silva et al. , 2015). We note that, similar to the GR branch of NSs, scalarized solutions in the original DEF model are stable up to the maximum mass of the corresponding branch, 7 and the stable solutions are generally energetically favorable over the GR NSs. This is further discussed in Sec. III.B.1.  \nLet us briefly comment on the exact definition of mass in scalar-tensor theories (shown in Fig. 1 and elsewhere). In contrast to GR, the definition of mass in scalar-tensor theories is subtle, due to the fact that these theories violate the strong equivalence principle. This results in the appearance of different possible masses as a measure of the total energy of the star (Lee, 1974; Scheel et al. , 1995a,b; Whinnett, 1999; Yazadjiev, 1999). These works showed that only the so-called tensor mass has natural energylike properties. For example, the tensor mass is positive definite, it decreases monotonically by the emission of GWs and it is well defined even in dynamical spacetimes (Lee, 1974; Scheel et al. , 1995a,b). In addition, only the tensor mass leads to a physically acceptable picture since it peaks at the same point as the particle number, a property crucial for the stability of the static stars (Whinnett, 1999; Yazadjiev, 1999). Therefore, the tensor mass,  \nwhich is defined as the ADM mass in the Einstein frame, should be taken as the physical mass. As a matter of fact, though, for most of the commonly used coupling functions the Jordan-frame and the Einstein-frame NS masses are identical.  \n̸  \nAfter the discovery of the phenomenon, scalarization was examined for a larger set of parameters and in more detail by Damour and Esposito-Far'ese (1996a). They showed that the presence of some externally imposed scalar-field background φ 0 , as well as considering α 0 = 0, smoothens the transition to a scalarized state. This is what we described as induced scalarization in Sec. II.C.4. (Strictly speaking, we do not have pure scalarization when α 0 = 0, since GR is no longer a solution of the field equations.) Spontaneous scalarization was further studied by Salgado et al. (1998), who considered the problem in the Jordan frame and performed an approximate Newtonian analysis of the system. They showed that scalarization can also be associated with the fact that the effective gravitational constant in scalar-tensor theories decreases for large scalar fields. It was further argued by Whinnett and Torres (2004) that scalarization leads to violation of the weak energy condition in the inner regions of NSs, which can cause instabilities. It was later demonstrated by Salgado et al. (2004) that this is not a general feature of scalar-tensor theories and that there are subclasses of the theory where the weak energy condition is easily satisfied.  \n̸  \nAnother consequence of scalarization is that for sufficiently large β 0 the maximum allowed mass for NSs increases compared to GR 8 , which can have various observa-  \n/s8857  \n2.5  \n2  \n1.5  \n1  \nEOS APR  \nGR (  \n= -4.5  \n/s106  \n= 0  \n= -4.8  \n0  \n/s32  \n0  \n= -5.0  \n0  \n/s98  \n/s98  \nstable branch  \n/s98  \nStability (line pattern):  \nunstable branch  \n11  \n12  \n13  \n)  \nR  \n[km]  \nM  \n/  \nD  \n/s8857  \n1  \n0.8  \n0.6  \n0.4  \n0.2  \n0  \n1  \n1.5  \n2  \n2.5  \nEOS APR  \n= -4.5  \n= -4.8  \n0  \n/s32  \n0  \n= -5.0  \n0  \n/s98  \n/s98  \n/s98  \nM  \n/M  \n/s8857  \ntional consequences. This problem was studied by Sotani and Kokkotas (2017), who took the effects of various microphysics parameters into account. Empirical relations were derived for the maximum mass of scalarized NSs that are parametrized with respect to the nuclear saturation parameters and the maximum sound velocity in the core.  \nUntil now we have discussed works that considered negative β 0 . However, as we mentioned, when the trace of the energy-momentum tensor is positive, scalarization can occur for β 0 > 0 as well, and this scenario introduces qualitative differences relative to our story thus far (Mendes, 2015). In particular, for the commonly used coupling function α ( φ ) = β 0 φ of the DEF theory, scalarized stars are not stable for high values of β 0 ( β 0 ≫ 1), which is further discussed in Sec. III.B.2.  \nThe calculation of NS parameters spanning the DEF theory parameter space provides a challenging technical task. In particular, tests of this theory against binarypulsar observations (described in Sec. III.A.1.b) require knowledge of the scalar charges (and their derivatives with respect to the star's mass) for a large catalog of EOSs. This calculation was performed most extensively in the works by Anderson and Yunes (2019) and Guo et al. (2021), who provide the results in tabulated form or through surrogate models. Yagi and Stepniczka (2021) [see also Horbatsch and Burgess (2011)] computed scalar charges in the DEF model analytically using a combination of perturbative weak-field expansion and Pad'e resummation and found excellent agreement with numerical calculations.  \nSpontaneous scalarization for the case of several nearby compact objects was considered by Cardoso et al. (2020). Their analytical analysis showed that, even though an isolated body might be below the threshold for scalarization, a collection of such bodies could develop a nonzero scalar field while maintaining average compactness much below the scalarization limit.  \nWe have already discussed spontaneous scalarization of NSs, but other compact objects can also scalarize. The case of BHs is considered in Sec. IV, but as an additional nonvacuum example we now discuss the case of boson stars [cf. Liebling and Palenzuela (2023) for a review], which were also shown to scalarize (Whinnett, 2000). In such systems, one has a complex scalar field as a matter source for the boson star and an additional real scalar field responsible for scalarization similar to the DEF model. The dynamics of this process was examined by Alcubierre et al. (2010) who showed that nonlinear development of the scalar field is observed in the absence of self-interactions in the complex scalar field. Ruiz et al. (2012) studied spontaneous and induced scalarization starting with initial data corresponding to stable boson stars in GR. They showed that a strong emission of scalar radiation occurs during the scalarization process.  \nb. Observational constraints from binary pulsars: To date binary pulsars have set the best constraints on scalarized NSs in the DEF models (Kramer et al. , 2021; Zhao et al. , 2022). This is because determining different characteristics of these systems through pulsar timing can be made extremely precise by accumulating yearslong observations. One of such important observables is the rate at which the binary's orbit decreases by energy loss through GW emission. In contrast to GR, the DEF model has an additional scalar degree of freedom that leads to a new channel of energy loss. Thus, the shrinking of the orbit should happen faster. The energy flux of the scalar-dipole radiation that gives the dominant contribution is given by (Damour and Esposito-Far'ese, 1992),  \nF dipole scalar = A (1 + q 1 q 2 ) 3 ( q 2 -q 1 ) 2 , (32)  \nwhere A is a function depending on the properties of the binary, such as its total and reduced masses, and its eccentricity; q 2 and q 1 are the scalar charges of each binary component normalized with respect to their masses q i = D i /M i ( i = 1 and 2).  \nBinary-pulsar observations were considered in the context of scalarized NSs for the first time by Damour and Esposito-Far'ese (1996a). They calculated the gravitational form factors (also known as sensitivities) of slowly rotating NSs, which form the set of coupling constants appearing in the post-Keplerian description of the binary in scalar-tensor theory; Horbatsch and Burgess (2012) for a summary. Only a few such binary systems were known at the time, and β 0 was constrained to be greater than -5 using polytropic EOSs. These results were later refined to include realistic EOSs and a limit of β 0 > -4 . 5 was derived by Damour and Esposito-Far'ese (1998). In addition, an estimate was made that it would be difficult for LIGO and VIRGO to improve β 0 bounds for most EOS possibilities [some exceptions were pointed out by Sampson et al. (2014) and Shao et al. (2017)]. The reason is that even though the merger events observed using such GW detectors can lead to much stronger scalar-dipole radiation, they are inferior in accuracy compared to the radio observations of binary pulsars, leading to weaker overall constraints. Next-generation detectors such as the Cosmic Explorer and the Einstein Telescope, though, will be able to improve the bounds on scalar-dipole radiation (Sampson et al. , 2014; Shao et al. , 2017).  \nAn example of how constraints can be imposed on DEF models is presented in Fig. 2, which illustrates how the model fails the double pulsar test for specific values of the theory parameters. In this case, the failure is due to the additional energy loss from scalar GWs predicted by the DEF model, predominantly the dipolar contribution.  \nWith the advances in observational astronomy, more pulsars in binary systems suitable for constraining the scalar-dipole radiation have been discovered, a complete and up-to-date list can be found in Freire (2022). Consequently, observational bounds on the scalar-tensor gravity  \nFigure 2 Mass-mass diagram for the double pulsar PSR J07373039A/B for the DEF model with α 0 = 5 × 10 -4 , β 0 = -4, and the assumption that the NSs are described by EOS MPA1. Various measured post-Keplerian parameters are shown as different curves, with the width indicating the measurement uncertainty in each parameter. For the point α 0 = 5 × 10 -4 and β 0 = -4 of the parameter space of the DEF model to be consistent with observations, all curves would have to intersect in a region of the mass-mass plane. This is not the case, as can be seen with the ˙ ω (periastron advance) and ˙ P b (change of orbital period) curves. For further details and the definitions of the other post-Keplerian parameters see Kramer et al. (2021). From Kramer et al. (2021).  \n<!-- image -->  \n/circledot  \nparameter β 0 for such systems have been widely discussed in the literature (Antoniadis et al. , 2013; Chiba, 2022; Esposito-Far'ese, 2004; Freire et al. , 2012; Kramer et al. , 2021; Shao et al. , 2017; Shibata et al. , 2014; Voisin et al. , 2020; Wex, 2014; Zhao et al. , 2022). The strongest current limit comes from Zhao et al. (2022), who practically closed the scalarization window for the original DEF model; i.e., the possibility for scalarization is ruled out in this theory. As we later discuss, though, there are a number of other well-motivated models where scalarization is still possible or cannot be constrained at all by binary-pulsar observations. These include theories with a massive scalar field, tensor-multiscalar theories, or even the standard DEF model when rapid rotation of NSs, which enhances the effect of scalarization, is considered. Furthermore, the theoretical and numerical approaches developed for the study of the DEF model are still applicable to these generalized theories in most situations. Thus, we spend considerable time here on the aspects of the DEF model despite its original form being essentially ruled out.  \nConstraints on scalarization with β 0 > 0 using pulsartiming observations were investigated by Mendes and Ottoni (2019). Owing to the fact that the scalar charge  \nM  \n/  \nM  \nFigure 3 Mass-radius relation for NSs in the DEF model for β 0 = 0 , -4 . 5 and -4 . 8 and employing the APR EOS. Nonrotating NSs are shown as solid curves, whereas stars rotating at the mass-shedding limit are depicted as dashed lines. As in Fig. 1, stars unstable to radial oscillations are shown as dotted lines. The radius of rotating stars refers to their equatorial radius.  \n<!-- image -->  \nis suppressed as β 0 increases while the range of masses allowing spontaneous scalarization decreases, it turns out that only weak constraints can be imposed by binarypulsar observations in this part of parameter space.  \nc. Rotating scalarized neutron stars: Thus far we have commented only on static NS models. All observed NSs are at least slowly rotating, and some dynamical processes such as NS mergers or stellar core collapse can produce relatively long-lived rapidly rotating supramassive protoneutron stars. Hence, the inclusion of rotation in NS physics is an inseparable part of the goal to explore their astrophysical implications.  \nDamour and Esposito-Far'ese (1996a) were the first to study slowly rotating scalarized NSs to leading order in rotation frequency O (Ω) using the formalism of Hartle (1967) and Hartle and Thorne (1968), which allowed them to calculate the NS moment of inertia; see also Sotani (2012). It is interesting that in this case there is an exact analytical solution for the NS exterior (Damour and Esposito-Far'ese, 1996a). Static and slowly rotating NSs for a wide range of realistic EOSs, including examples with hyperons or quark matter, were considered by Altaha Motahar et al. (2017). The extension to second order in the rotational frequency O (Ω 2 ) was made by Pani and Berti (2014), who used the extension to calculate rotational corrections to the stellar radius and mass, and also its quadrupole moment.  \nRapidly uniformly rotating scalarized stars, without approximation, were obtained by Doneva et al. (2013).  \nThey showed that for a fixed β 0 the maximum deviation from GR that is achieved at the mass-shedding limit is considerably larger than in the static case, and the range of central energy densities where scalarization is possible is significantly broadened. This can be seen in Fig. 3, where sequences of static scalarized NSs are compared to NSs rotating at the Kepler limit for the same values of β 0 .  \nThere are a number of factors leading to differences from the static case. The first, more intuitive one is that the rotational energy of the star also acts as a source for the scalar field, and thus can change the onset and degree of scalarization. Meanwhile, the rapidly rotating models tend to be less compact, which can reduce the degree of scalarization. The large deviations from GR compared to the static case, conversely, are due mainly to the fact that scalarized stars can sustain much larger angular momentum before reaching the Kepler limit. This is a nonlinear effect that could not be normally caught in the slow-rotation approximation.  \nA natural consequence of the aforementioned rotation effects is that the minimum | β 0 | where scalarization is possible changes compared to the static case. Thus, for the same EOS II (Diaz Alonso and Ibanez Cabanell, 1985) used by Damour and Esposito-Far'ese (1996a), scalarization happens for β 0 < -3 . 9 in rotating stars (Doneva et al. , 2013) compared to β 0 < -4 . 35 in the static case (Damour and Esposito-Far'ese, 1996a). Therefore, one can conclude that, although binary-pulsar observations seem to rule out DEF scalarization for static or slowly rotating NSs (Zhao et al. , 2022), there is still an observationally viable range of β 0 where only rapidly rotating NSs can scalarize. This is potentially relevant for binary mergers and stellar core collapse, where such rapidly rotating NSs can form.  \nA caveat in the previous argument is that one does not expect the star to rotate uniformly in these extremely dynamical events. Such differential rotation was first studied by Doneva et al. (2018b), who adopted a simple rotation law that can still capture some of the main properties of the merger remnants, especially a few tens of milliseconds after the merger of the binary (Bauswein and Stergioulas, 2017).  \nWhen scalarization is considered, larger values of the maximum mass as well as of the angular momentum can be achieved for supramassive NSs compared to GR. Moreover, the scalar field causes rapidly rotating models to be less quasitoroidal than their general-relativistic counterparts. This can have direct astrophysical implications, especially for binary NS mergers, where the maximum possible mass and angular momentum that a NS can sustain are crucial for determining the merger outcome and the lifetime of the merger remnant (in case a supramassive  \nM  \n/  \nM  \n6  \nFigure 4 Mass as a function of radius for the APR4 EOS. The results for different values of the coupling constant β 0 and mass of the scalar field µ , given in eV, are plotted. The potential is assumed to have the form V ( φ ) = ( µ 2 / 2) φ 2 , with the scalar field defined though the action (25). From Yazadjiev et al. (2016).  \n<!-- image -->  \nor hypermassive NS forms 9 ); see Bauswein et al. (2013); Bauswein and Stergioulas (2017); Takami et al. (2014); and Weih et al. (2018).", '2. Massive scalar field': 'A key property of scalar-tensor theories that was neglected in the previously discussed studies is the possibility of having a nonzero scalar-field potential. The simplest case is to take a potential that leads to a nonzero scalarfield mass µ , but more complicated potentials, such as those with self-interaction terms, can be considered as well; see Eq. (2). Although this seems like a simple extension, it has a dramatic effect on the observational properties of NSs, especially on the GW emission. The reason lies in the different asymptotic behavior of the scalar field. In the case of zero potential, the scalar field decreases as D/r at infinity according to Eq. (31). This leads to a nonzero scalar charge D and thus nonzero scalar-dipole radiation. In the presence of nonzero scalar-field mass µ , though, the scalar field tends exponentially to zero after some characteristic distance related to its Compton wavelength λ φ = 2 π/µ as discussed by Ramazano˘glu and Pretorius (2016). Hence, the scalar field is effectively  \nconfined to a characteristic radius and its scalar charge is zero. If the orbital separation between the two objects in a binary-pulsar system is much larger that λ φ , the dynamics will not be directly altered by the scalar field and there is no significant emission of scalar-dipole radiation (Alsing et al. , 2012). Since λ φ is controlled by µ in the simplest case of a scalar-field potential, one can reconcile the DEF model (for arbitrarily small β 0 ) with binary pulsar constraints by giving the scalar field a mass µ ≳ 10 -16 eV.  \nScalarized NSs in massive scalar-tensor theories were first studied in the static case (Chen et al. , 2015; Popchev, 2015; Ramazano˘glu and Pretorius, 2016) and later extended to slow (Yazadjiev et al. , 2016) and rapid rotation (Doneva and Yazadjiev, 2016). The inclusion of a quartic self-interaction term to the potential was considered by Staykov et al. (2018). These works showed that the mass of the scalar field and the self-interaction have similar effects on the scalar field around NSs and that they both suppress scalarization. The quartic interaction by itself cannot affect the range of central energy densities where scalarized solutions exist, because it is a nonlinear contribution to the linearized scalar-field equation of motion [recall Eq. (11)]. In contrast, the mass term shrinks the domain of existence of scalarized NSs and, for large enough masses, no scalarization is possible at all. This is evident in Fig. 4, where NS mass is plotted as a function of its radius for different combinations of β 0 and µ . The massive scalar-field solutions are confined between the zero scalar-field mass models (the original DEF models) and the GR ones, corresponding loosely speaking to µ →∞ . Note that for µ ≳ 10 -16 eV the solutions are almost indistinguishable from the massless DEF model. For this reason the latter represents an upper limit on the possible deviations from GR in massive scalar-tensor theories.  \nThe exponential asymptotic behavior of the massive scalar field brings computational challenges to the construction of scalarized NSs, which led to new numerical approaches (Rosca-Mead et al. , 2020a) and also facilitated the construction of NSs for highly negative β 0 and large scalar-field masses. Rosca-Mead et al. (2020a) showed that for sufficiently negative β 0 qualitative changes in the strongly scalarized branch of solutions are possible. For example, the maximum of the scalar field can be located away from the stellar center. In their most extreme form, these solutions are composed of a highly compact NS model surrounded by a scalar-field shell. Also, Tuna et al. (2022) showed that some scalarized solutions in this part of the ( β 0 , µ )-parameter space gave indications of metastability: they were stable to small perturbations but had lower binding energy than their GR counterparts.  \nAn extension of these results to other forms of coupling functions and scalar-field potentials is the asymmetron model (Chen et al. , 2015; Morisaki and Suyama, 2017). It is interesting because of the fact that the asymmetron  \nmodel realizes proper cosmic evolution, and it can also account for the cold dark matter. Chen et al. (2015) and Morisaki and Suyama (2017) focused especially on large scalar-field masses spanning several orders of magnitude and having a Compton wavelength shorter than 10km, which is the typical size of a neutron star.', '3. Incorporating further physics': "New aspects of the original DEF model were recently studied with the inclusion of different physical details. For instance, Silva et al. (2015) studied the presence of anisotropic pressure of nuclear matter for both static and slowly rotating NSs. The motivation for this comes from the fact that some theoretical considerations, for instance, with magnetic fields or within the Skyrme model (a lowenergy EFT of quantum chromodynamics) (Nelmes and Piette, 2012) suggest that at high densities the NS EOS might have a significant degree of anisotropy (Herrera and Santos, 1997). In such a case, the effects of scalarization increase (decrease) when the tangential pressure is larger (smaller) than the radial pressure. The threshold value of β 0 for the development of scalarization, which in the isotropic case is β 0 < -4 . 35, can be increased due to the presence of anisotropy, thus widening the range of parameters in which scalarization is possible.  \nAnother astrophysically interesting extension of scalarization is to include the magnetic fields. According to observations and modeling, NS magnetic field values can span from 10 8 to 10 12 G for standard 'old' pulsars, ranging from 10 16 G at the surface of some magnetars to, hypothetically, as high as 10 17 -10 18 G in the cores of newly formed protoneutron stars. Such strong magnetic fields impact the properties of scalarized NSs, including their magnetic deformability, maximum mass, and range of scalarization, as studied by Soldateschi et al. (2020). They found a magnetically induced spontaneous scalarization whose essence is the following: strong toroidal magnetic fields can support descalarized configurations and, if the star's magnetic field decreases during some nonideal magnetohydrodynamical process, the star can undergo a rapid growth of the scalar field; i.e., it scalarizes. The magnetic quadrupolar deformations of scalarized NSs and the related GWs produced by rotating magnetars were studied by Soldateschi et al. (2021)  \nAnother interesting extension to the standard DEF model is related to challenging the idea that the fundamental physics remains unchanged in the star's interior, which is a common assumption when a nuclear matter EOS is constructed. This was studied by Coates et al. (2017), who considered two models in which the mass of the photon had a different value in the interior and the vicinity of a compact star compared to the mass measured by experiments performed in a weak-gravity regime. The first model is based on a Proca-like mass with an effec-  \ntive mass term dependent on φ . The second model can be thought of as a gravitational Higgs mechanism where the Higgs potential is replaced by the scalar-gravity coupling. In both cases the scalar field undergoes spontaneous scalarization, thus acquiring a nontrivial profile if the compactness passes a certain threshold, providing a mass to the photon by coupling to it in an appropriate manner. Although the focus of Coates et al. (2017) was on the electromagnetic field as a proof of principle, these results can be extended to other fields of the standard model. The signatures of such a gravitational Higgs mechanism on the behavior of magnetic field of NSs in Einstein-Maxwell theory was studied by Krall et al. (2020).", 'B. Dynamics of scalarized neutron stars and binary mergers': 'The dynamics of isolated NSs can be studied by solving the full nonlinear field equations of scalar-tensor gravity, which is often a challenging task. Instead, one usually first approaches the problem by linearizing the field equations around a background solution and then analyzing the resulting linearized dynamics. The study of nonlinear dynamics is then done when necessary and feasible. We follow this sequence in the Secs. III.B.1-III.B.4. We then review what happens when NSs in scalar-tensor theories are placed in binary systems.', '1. Linearized dynamics': "The studies of linearized dynamics concern the stability of scalarized stars and the analysis of their quasinormal mode (QNM) spectrum. The latter involves the study of different classes of NS oscillation modes that are tied to the emission of GWs. See Kokkotas and Schmidt (1999) for a review.  \na. Stability: We have seen that scalarized NSs coexist with their nonscalarized counterparts as solutions within the DEF model. Which of these branches of solutions is the one realizable in nature? One way of answering this question consists of calculating the fractional binding energy M 0 /M -1, where M 0 is the star's baryonic mass. Through this calculation, Damour and Esposito-Far'ese (1993) showed that scalarized NSs are usually energetically favored over the GR ones. This is evident in Fig. 5, where we plot the binding energy as a function of the baryonic mass. At constant baryonic mass, the scalarized solutions (solid curve) have larger binding energies relative to the GR solutions (dashed curve) and are thus energetically favorable. In addition, we can see a cusp at the maximum of the mass for both the scalarized and nonscalarized NSs. This suggests a change of stability. Since the branches beyond the maximum of the mass have a lower binding  \n1  \n-  \nM  \n/  \nM  \nFigure 5 Fractional binding energy M 0 /M -1 as a function of the baryonic mass M 0 for scalarized NSs in the DEF model with β 0 = -5 using the APR EOS. Inset: enlargement of the scalarized branch showing the formation of a cusp at the maximum of the mass in the upper-right corner.  \n<!-- image -->  \nenergy (dotted curve in the inset), they are unstable.  \nThe aforementioned stability analysis relied on the bulk properties of the star. A rigorous complementary approach considers the linear perturbations of the star. The first step in this direction was taken by Harada (1997), who studied scalar-field perturbation in the background of a NS with a zero (or constant) scalar field within the DEF model. Harada (1997) studied the perturbation equations in the frequency domain and showed that the GR solution becomes unstable after a specific critical central energy density. This is the point where the scalarized solutions branch out from the GR ones; see Fig. 1. Harada (1998) reached similar conclusions but worked in the context of catastrophe theory. The radial stability of scalarized NSs was also studied by Mendes and Ortiz (2018). They considered metric and scalar-field perturbations for both signs of β 0 . They found that scalarized NSs are stable against linear perturbations and that instability takes place past the point of maximum mass.  \nb. Gravitational waves from perturbed neutron stars: Linearized perturbations are also helpful for studying the NS oscillation modes directly related to GW emission. We first note that GWs in scalar-tensor gravity can carry additional polarizations compared to GR. In particular, one can have breathing modes in addition to the standard 'plus' and 'cross' polarizations of GR (Will, 2018). Moreover, radial perturbations in scalar-tensor gravity can excite GWs, contrary to what happens in GR. These perturbations source monopole scalar waves that result in a nonvanishing contribution to the perturbed Jordanframe Riemann tensor (linearized around a Minkowski  \n<!-- image -->  \n<!-- image -->  \n(b)  \n(  ) ( F ) ( H 1) (a) (b) (c) Figure 6 (a) Frequency and (b) inverse damping time of the three lowest-frequency radial modes of stellar models in scalar-tensor theories for the coupling function A ( φ ) = [ cosh ( √ 3 β 0 φ )] 1 / (3 β 0 ) , denoted by M1, and A ( φ ) = e β 0 φ 2 / 2 , denoted by M2 with β 0 = -5, as functions of compactness. Lines (points) represent values computed with frequency-domain (time-domain) techniques. F stands for the fundamental radial oscillation mode (with no nodes) and H 1 for its first overtone (with one node), while φ denotes the fundamental scalar mode (with no nodes). (c) Depiction of the fundamental scalar mode showing its real and imaginary parts. The employed EOS is a two-phase polytrope consisting of a stiff core and a soft crust. From Mendes and Ortiz (2018).  \n<!-- image -->  \n(  )  \n(  \nF  \n)  \n( H 1) ( F , H 1) background) in the transverse-traceless gauge (Damour and Esposito-Far'ese, 1992; Novak and Ibanez, 2000). For DEF-like scalar-tensor theories, this requires α 0 = 0. This contribution is then linked to the existence of a breathing polarization mode of the GW. We can then conclude that radially oscillating scalarized NSs in scalar-tensor theories with α 0 = 0 will emit GWs. Their amplitude is connected to α 0 that controls the generation of tensorial waves from the dynamics of the scalar field. Thus, larger α 0 leads to stronger coupling and stronger excitation of the breathing modes.  \n̸  \nThe first study of radial oscillations of NSs in the DEF model was performed by Sotani (2014) in the Cowling approximation (Cowling, 1941; McDermott et al. , 1983). In the original GR version of this approximation, the spacetime is held fixed while only the fluid is perturbed. In scalar-tensor gravity, the scalar-field perturbations are also often neglected. Despite its limitations, the Cowling approximation can actually capture well the qualitative features of the neutron star oscillation spectrum. The full problem (i.e., when both the metric and the scalarfield perturbations are taken into account) was addressed by Mendes and Ortiz (2018). They found a new family of modes named scalar modes that have no counterpart in GR. The results show that they have distinct frequencies and damping times compared to the fluid radial oscillation modes. In Fig. 6 we show the frequencies and damping times of the fundamental fluid radial oscillation mode and its first overtone, as well as the fundamental scalar radial mode, as functions of the stellar compactness. We see in Fig. 6(c) that scalarized NSs become unstable at the maximum mass, and it is the scalar mode that is responsible for the instability. This contrasts with GR, where it is the fundamental radial fluid mode that becomes  \n̸  \nunstable at the maximum of the mass.  \n-mode trivial sol. Not only do the radial oscillations of scalarized NSs differ significantly from those in the GR case, but the nonradial modes that are related to the tensorial gravitational wave emission can be strongly influenced by the scalar field. The behavior of nonradial oscillations modes, which are sources of the usual tensor polarizations of GWs, can also be strongly affected by the scalar field. These perturbations can be classified as 'polar' or 'axial' depending on how they behave under parity transformations (Regge and Wheeler, 1957; Thorne and Campolattaro, 1967).  \nIn particular, scalar-field perturbations are of the polar type, meaning that they couple only to polar perturbations of the fluid and metric. In GR, these perturbations are the most efficient GW sources. Sotani and Kokkotas (2004) were the first to study the nonradial oscillations of scalarized NSs. They used the Cowling approximation, in which only polar-parity fluid perturbations are dynamical. This simplifies the problem considerably and Sotani and Kokkotas (2004) calculated the fundamental f -mode frequency and the pressure p -mode frequency. In a follow-up work, Sotani and Kokkotas (2005) went beyond the Cowling approximation and derived equations for both axial and polar perturbations that included metric perturbations. They analyzed only the simpler axial perturbations (discussed later). A study of polar perturbations became possible only after new techniques were developed in GR by Kruger and Kokkotas (2020a,b). These techniques were then used in scalar-tensor theory by Kruger and Doneva (2021), who also considered self-interacting massive scalar fields. 10 Kruger and Doneva (2021) found that  \n]  \nz  \nH  \nk  \n[  \ny  \nc  \nn  \ne  \nu  \nq  \ne  \nr  \nF  \nFigure 7 f -mode frequencies (for ℓ = 2) as a function of the average NS density √ M/R 3 for several values of β . The SLy EOS is used (Douchin and Haensel, 2001). From Kruger and Doneva (2021).  \n<!-- image -->  \nthe scalar field leaves clear imprints on the oscillation frequencies of NSs. An extension of these results for rapidly rotating NSs, but now back to the Cowling approximation, was made by Yazadjiev et al. (2017). They also studied the Chandrasekhar-Friedman-Schutz instability driven by rotation (Chandrasekhar, 1970; Friedman and Schutz, 1978).  \nR M/ An interesting property of NS oscillations in GR is the almost linear and EOS-independent relation between the quadrupole ( ℓ = 2) f -mode frequency and the average density of the NS (Andersson and Kokkotas, 1998). Sotani and Kokkotas (2004) and Kruger and Doneva (2021) showed that this linear scaling is dramatically broken by scalarization. We show this in Fig. 7, where we also see that the deviations from GR increase with an increase of | β 0 | .  \nAnother class of modes that can be attributed, loosely speaking, to the 'oscillations' of the spacetime itself are the axial spacetime w modes. These modes are somehow easier to calculate (without approximations) because, in this case, the perturbations of the fluid and the scalar field are zero. The axial modes of scalarized NSs were considered for the first time by Sotani and Kokkotas (2005), who calculated the frequencies and the damping times of different classes of w modes. Extensions of these results to a variety of realistic EOSs including nuclear, hyperonic, and hybrid matter, were carried out by Altaha Motahar et al. (2018), while the case of massive self-interacting scalar fields was studied by Altaha Motahar et al. (2019).  \nNumerical calculations showed that the effect of scalarization is stronger on the damping times than the effect on the frequencies, and in general the values of both are lower than their GR values. In addition, EOS-independent relations between w -mode properties known to exist in GR can also be obtained in scalar-tensor theory.  \nA class of NS modes related to the crustal torsional oscillation was studied by Silva et al. (2014) for the DEF model in the Cowling approximation. These oscillations probably follow the giant flares in soft gamma-ray repeaters (Israel et al. , 2005; Strohmayer and Watts, 2005, 2006) and are associated with motions in the NS crust. Silva et al. (2014) found that, for values of β 0 consistent with binary-pulsar constraints at the time, the effect of scalarization on the torsional oscillation frequencies is smaller than the uncertainties in the microphysics modeling of the crust.", '2. Nonlinear stability and collapse to a black hole': 'While the stability and oscillations of a scalarized NS can be studied perturbatively, the formation of scalar hair starting from a GR NS solution is a fully nonlinear process. More precisely the initial exponential growth of the scalar field starting from an unstable GR solution can be modeled by linearized dynamics, but the subsequent saturation of the scalar field to an equilibrium value is a nonlinear phenomenon.  \nThe transition from a nonscalarized to a scalarized NS was first considered by Novak (1998a) [see also Degollado et al. (2020)], who studied the nonlinear evolution in spherical symmetry. The numerical simulations of Novak (1998a) showed that the scalar field for an unstable GR NS first grows exponentially and then saturates, with the system saturating to an equilibrium scalarized end state. 11 A consequence of the work by Novak (1998a) was the proof of nonlinear stability of scalarized NSs. That is, the numerical simulations showed that the scalarized NSs have stable evolution if the full nonlinear system of field equations (in spherical symmetry) is considered. Although effort was not strictly dedicated to a stability analysis, a fully nonlinear evolution of scalarized NSs was performed in a series of papers that we discuss later. This provides strong support for the stability of these objects under the most general circumstances.  \nNovak (1998b) studied the collapse of a scalarized NS to a BH in spherical symmetry. Since BH no-hair theorems include DEF models [see Herdeiro and Radu (2015) for a review], the resulting BH will be bald and the scalar  \nfield has to be radiated away during collapse. 12 This takes place in the form of scalar waves, in analogy to what we discussed in the case of radial NS oscillations in Sec. III.B.1.  \nOur previous discussion covered the β 0 < 0 case. Palenzuela and Liebling (2016) and Mendes and Ortiz (2016) studied the end state of the tachyonic instability in scalartensor theories for representative coupling functions with β 0 > 0 and realistic EOSs. This was done both through an energy balance analysis of the existing equilibrium configurations and by nonlinear dynamical simulations. They found that (contrary to the β 0 < 0 case) the final state of the instability is highly sensitive to the details of the coupling function, varying from gravitational collapse to spontaneous scalarization. They also found that in the original DEF model [where α ( φ ) = β 0 φ with α 0 = 0] scalarized solutions can become unstable compared to the GR ones when β 0 ≫ 1. However, stability can be recovered for all values of β 0 by considering different coupling functions. This is the case for coupling functions with bounded values, such as α ( φ →∞ ) = α ∞ for some constant α ∞ . This distinction could give rise to novel astrophysical tests for determining the detailed form of the coupling.', '3. Stellar core collapse': "Once we know that equilibrium scalarized NS solutions exist in the DEF model and they are stable against linear perturbations, the next step is to study how they are formed. Isolated NSs can form after the collapse of the core of a massive star during a supernova explosion. A NS on the other hand can collapse to a BH if a threshold for the mass and the angular momentum is reached. All these are highly dynamical nonlinear processes that can have strong observational signatures in both the electromagnetic and GW signals.  \nThe first study of a degenerate stellar core collapse (more specifically with white dwarf initial data) to a scalarized NS through a bounce and the formation of a shock was performed by Novak and Ibanez (2000). The simulations were done in spherical symmetry, which allowed them to calculate only the resulting gravitational monopolar radiation. They found that the emitted breathing modes can be potentially detected by LIGO or VIRGO. Notably the emitted signal will be substantially different than the collapse of a NS to a BH, allowing one to distinguish between them.  \nGerosa et al. (2016) studied the problem of spherically symmetric core collapse in the DEF model in further detail. Two types of initial data were used: collapse of a  \n3  \nFigure 8 Scalar wave rφ during a stellar core collapse extracted at 5 × 10 4 km. The legend lists deviations from the fiducial parameters µ = 10 -14 eV , α 0 = 10 -2 , β 0 = -20 , Γ 1 = 1 . 3 , Γ 2 = 2 . 5, and Γ th = 1 . 35, where µ is the scalar-field mass and Γ 1 , Γ 2 , and Γ th are parameters of the polytropic EOS and the thermal contribution. See Sperhake et al. (2017) for more details. From Sperhake et al. (2017).  \n<!-- image -->  \nstellar iron core and collapse of 'realistic' NS progenitors that were obtained from computations of stellar evolution. Depending on the theory parameters three possible outcomes of the core collapse are possible - collapse to a GR NS, collapse to a scalarized NS, or collapse to a shortlived protoneutron star followed by a nonscalarized BH formation. It was in the last case that the most prominent GW signal with a clear signature from the presence of nontrivial scalar fields was observed. While the fluid dynamics during the collapse is only weakly affected by the scalar field, the converse is not true: the scalar radiation depends strongly on the specifics of the matter collapse as well as on the choice of the coupling parameters α 0 and β 0 .  \nThe inclusion of a scalar-field mass to the core-collapse simulations has some interesting consequences (Sperhake et al. , 2017). The waveform rφ during such events is shown in Fig. 8. The mass itself does not have a large influence on the dynamics, but it allows the theory to be reconciled with binary-pulsar observations for a much larger range of α 0 and β 0 , as discussed in Sec. III.A.2. Thus, a dramatic increase in the radiated GW signal observable even with the existing LVK detectors was predicted.  \nA prominent feature is that we expect to receive an inverse chirp signal from such events that will last for years, with a near monochromatic signature on timescales of ∼ 1 month. In fact, the inverse chirp signal is connected to the scalar field's mass and, consequently, to the dispersion relation. Thus, the scalar wave burst that is emitted during the collapse is effectively stretched in time with decreasing amplitude. Extensions of these results to the case of a self-interacting scalar-field potential were made by Cheong and Li (2019) and Rosca-Mead et al. (2019). Constraints on the theory based on the  \nscalar-field evolution in the Einstein frame were imposed by Geng et al. (2020). The problem was considered in greater detail by Rosca-Mead et al. (2020b), who found the three possible scenarios of the collapse outcome in scalar-tensor theory with sufficiently negative β 0 . These constitute the formation of a BH following multiple NS stages, the multistage formation of a strongly scalarized NS, and the single-stage formation of a strongly scalarized NS. Rosca-Mead et al. (2019) found that the resulting GW signal can reach a signal-to-noise ratio of over 20 for the existing GW detectors, which has the potential to put strong constraints on the theory.", '4. Dynamical scalarization and neutron star mergers': "Another highly dynamical and nonlinear process that has important astrophysical, and especially GW implications is binary NS mergers. This problem can be more challenging to solve than the stellar core collapse because it is not possible by construction to apply certain approximations, such as spherical symmetry. That is why the binary merger dynamics in the DEF model was addressed shortly thereafter, first by Barausse et al. (2013) and then by Shibata et al. (2014). The overall conclusions of both works are similar and they reside in the fact that, even if one or both of the NSs are not scalarized before the merger, they can develop nonzero scalar fields during the inspiral. This was called dynamical scalarization. The main advances made by Shibata et al. (2014) constitute using several realistic EOS with consistently derived bounds on the parameter β 0 , as well as developing an initial data solver for scalarized binary NSs. The overall results demonstrated a significant change of the inspiral GW signal at the moment of dynamical scalarization and afterward. The reason is that the inspiral is accelerated due to the scalar-dipole radiation and the total number of GW cycles is significantly decreased with respect to the general-relativistic case.  \nThe actual merger and the postmerger phases were addressed only by Shibata et al. (2014), who showed that dynamical scalarization can happen not only in the inspiral phase but also during the merger since a massive compact star is formed in this process. The evolution during the NS merger of the maximum values of the restmass density ρ max and scalar field φ max for several binary NSs with total mass M = 2 . 7 M ⊙ and several values of β 0 is shown in Fig. 9. Two EOSs are considered in the two panels and the cosmological value of the scalar field is taken to be φ 0 = 10 -5 in order to satisfy constraints from the binary-pulsar observations (Antoniadis et al. , 2013; Freire et al. , 2012). Some of the binary NSs are not scalarized at the beginning of the evolution (the models with φ max = 0 at the beginning), but the values of β 0 are chosen such that all of them develop scalar fields at a certain point of the evolution. The actual merger of  \nϕ  \nϕ  \n<!-- image -->  \nH4,  \n=5x10  \n-5  \nFigure 9 Evolution of the maximum values of the rest-mass density ρ max and scalar field φ max during NS merger for several models with a total mass M = 2 . 7 M ⊙ with the APR4 EOS (top panel) and the H4 EOS (bottom panel). The merger sets in at the time as the maximum density steeply increases. We note that for β 0 ≥ -4 . 2 with the APR4 EOS and for β 0 = -4 . 75 with the H4 EOS, the scalarization had already occurred at t = 0. From Shibata et al. (2014).  \n<!-- image -->  \nthe two NSs is marked by the rapid increase of φ max . Afterward, either a hypermassive or supramassive NS is formed or the merger remnant collapses to a bald BH and the scalar field is radiated away. In agreement with the studies of equilibrium differentially rotating NSs discussed in Sec. III.A.1, the scalarized merger remnant can sustain a larger mass without collapsing to a BH (Doneva et al. , 2018b). This is evident in the bottom panel of Fig. 9, where collapse to a BH is observed in pure GR, while a supramassive NS is formed after the merger in scalartensor gravity with small enough β 0 .  \nThe quasiperiodic oscillations of the merger remnant were examined by Shibata et al. (2014), showing a clear distinction compared to the GR case. Such oscillations were studied in a series of papers in GR (see e.g. Bauswein and Janka (2012); Bauswein et al. (2012, 2014); Clark et al. (2014); Hotokezaka et al. (2013); Maione et al. (2016); Rezzolla and Takami (2016); Takami et al. (2014, 2015)), mainly as a tool to determine the nuclear matter EOS  \nϕ  \nfrom the postmerger GW signal. The observed differences with the scalarized case can potentially be used to discriminate between GR and modified gravity theories. This is especially the case because scalarized merger remnants are produced only if specific initial conditions related to the mass of the merging compact objects and the specifics of the EOS are met (Shibata et al. , 2014).  \nTaniguchi et al. (2015) took a different approach in which they calculated quasiequilibrium sequences of binary NSs at different time instants instead of performing time evolution. This approach is well known in GR (Gourgoulhon et al. , 2001; Taniguchi and Gourgoulhon, 2003; Taniguchi and Shibata, 2010), and it is assumed that the characteristic time of the system to settle to equilibrium is much smaller than the inspiral timescale. This is supposed to give a relatively accurate picture of the binary evolution even close to the merger. The results are in agreement with Shibata et al. (2014), while the small deviations relative to Barausse et al. (2013) are probably due to the fact that GR initial data were used in the latter. Taniguchi et al. (2015) showed that the absolute value of the binary binding energy is smaller than in GR. In addition, the GW cycles prior to the merger were significantly reduced compared to pure GR once scalarization kicks in, and the effect is considerably stronger than the one due to tidal interactions.  \nThe post-Newtonian (PN) approximation has also been a valued tool to model NS inspiral in the DEF model. Initial work has focused on tensor-multiscalar theories and simple scalar-tensor theories that do not allow scalarization (Damour and Esposito-Far'ese, 1992, 1996b; Lang, 2014; Mirshekari and Will, 2013). This approach cannot be immediately applied to the case of dynamical scalarization, since this is a nonperturbative effect that is absent from the weak-field regime. Palenzuela et al. (2014) were the first to address this issue. They used the equations of motion at 2.5PN order, derived in scalar-tensor gravity by Mirshekari and Will (2013), modified in such a way that the changes in the stars' scalar charges are taken into account. More specifically, they solved a system of nonlinear algebraic equations at each step of the orbital evolution to compute the scalar charges. The resulting inspiral evolution was found to be in agreement with the numerical simulations of Barausse et al. (2013). The approach of Palenzuela et al. (2014) is computationally inexpensive, and this allowed them the to study large portions of the parameter space, including (un)equalmass and eccentric binaries. This approach can also be used to efficiently generate inspiral GW templates in the DEF model. A further step forward taken by Sennett and Buonanno (2016) introduced a methodology called post-Dickean expansion. Their main improvement with respect to Palenzuela et al. (2014) is a 1PN extension of the feedback mechanism. The post-Dickean expansion was compared against the quasiequilibrium calculations of Taniguchi et al. (2015), and it was shown that this can  \naccurately predict the onset and magnitude of dynamical scalarization.  \nSennett et al. (2017) took a different approach in the perturbative study of the inspiral that constitutes an analytical model of dynamical scalarization using an effective action. The motivation was to cure two deficiencies in the previous PN studies. In an effective action approach, the nonlinear scalarization process is reduced to a pair of cubic equations that have a closed-form solution depending on the binary separation when dynamical scalarization happens and on the magnitude of the developed scalar charge. This simplifies the problem relative to Palenzuela et al. (2014) and Sennett and Buonanno (2016). In addition, the effective action approach allows one to construct a simple two-body Hamiltonian that can be used to compute the binary's binding energy. Sennett et al. (2017) used this Hamiltonian, in combination with Landau's theory of phase transitions, to interpret dynamical scalarization as a second-order phase transition. Khalil et al. (2019) extended these results to general theories admitting scalarization for either BHs or NSs and are valid for adiabatic (quasistationary) and quasicircular orbits. Khalil et al. (2022) took it a step further, where these approximations were dropped and the dynamical evolution around the phase transition to the scalarized regime was studied. The results showed that in some cases assuming a quasistationary evolution might not be accurate enough even for quasicircular binaries.  \nLast, Ponce et al. (2015) studied the effect of the scalar field on the electromagnetic radiation emitted in the merger of magnetized NSs. They found that deviations in the emitted electromagnetic flux due to scalarization are not negligible yet are challenging to measure. However, if combined with GW observations, constraints on scalar-tensor theory can in principle be placed, showing the usefulness of multimessenger astronomy.", "C. Astrophysical implications of scalarized neutron stars in the Damour-Esposito-Far'ese model": 'In this section, we further discuss the astrophysical implications of scalarized NSs. Many aspects have already been covered, such as the binary-pulsar observations, NS oscillations, stellar core collapse, and binary mergers, as they naturally appeared in the presentation. Here, we shed further light on the possible astrophysical implications of scalarization, trying to be as complete as possible in two main areas. We first discuss the astrophysical implications directly related to electromagnetic observables. Afterward, we address the problem of universal relations for scalarized NS models.  \na. Electromagnetic observations: Scalarized NSs in the DEF model and its extensions were studied in a vari-  \nety of astrophysical scenarios in an attempt to probe the existence of the scalar field. Since NSs are often surrounded by accretion disks, it is natural to study the effect of a nontrivial scalar field on disk properties. The simplest model is called the thin disk model, in which particles are assumed to move on geodesics around the central compact object. Any small perturbation acting upon these particles will lead to oscillations around their equilibrium orbit with some characteristic epicyclic frequencies. Epicyclic and orbital frequencies are thought to be, in one way or another, related to the interpretation of different accretion disk properties. An example is the quasiperiodic oscillations (QPOs) observed in the spectrum of low-mass x-ray binaries. For this reason, some works calculated these frequencies in the spacetime of scalarized NSs. DeDeo and Psaltis (2004) did this for nonrotating NSs, Staykov et al. (2019) included slow rotation and also a mass to the scalar field, while Doneva et al. (2014b) studied rapidly rotating NSs. The conclusion of these works is that, if one takes into account current observational constraints, only the cases of rapid rotation and massive scalar fields leave room for a significant effect of the scalar field. Yet, a question that deserves further work is how one can disentangle uncertainties in the NS EOS and in the disk model from modifications to GR.  \nAnother prospective way to infer information about accreting compact objects is through the shape of the Fe Kα fluorescent line at 6.4 keV. Since it is emitted from the inner regions of the accretion disk, it carries traces of the underlying spacetime geometry. This line can be observed for both accreting NSs and BHs, and the accuracy of the observations is expected to be improved through future x-ray missions such as Athena. Bucciantini and Soldateschi (2020) calculated the shape of this line by taking the light propagation around a scalarized NS into account. They argued that the influence of GR modification both on the intensity of the low-energy tails and on the position of the high-energy edge of the line are potentially observable in the future.  \nAccretion onto NSs can also trigger a process called the gravitational phase transition, which was named in analogy with matter phase transitions from confined nuclear matter to deconfined quark matter (Kuan et al. , 2022). This process can happen when the maximum mass of the scalarized NSs is smaller than the maximum mass of the zero scalar-field (GR) solution. The idea then is that if a scalarized NS close to this maximum mass accretes some matter, it may pass beyond the stability point. In GR this would cause a collapse to a BH. In scalar-tensor theory, the star will radiate its scalar hair and evolve toward the zero scalar-field (i.e., the GR) branch. A significant amount of GWs can be produced in this process, and they will potentially be detectable with the next generation of GW detectors.  \nThe NS surface as an emitter of x-ray radiation is also an important probe in strong-field gravity because the  \nelectromagnetic radiation is emitted from a region with large spacetime curvature. In this regard, the simplest observable is the gravitational redshift of surface atomic lines. The redshift carries information about the NS mass, radius, spin, and, in scalar-tensor gravity, the scalar field. DeDeo and Psaltis (2003) showed that scalarization had a significant effect on the redshift, but only for negative enough values of β 0 . Such values are already ruled out by binary-pulsar observations. One could entertain the idea that a small scalar-field mass could allow for large deviations in redshift relative to GR, while reconciling the theory with binary-pulsar constraints (Doneva and Yazadjiev, 2016; Popchev, 2015; Ramazano˘glu and Pretorius, 2016). A study of this problem has not yet been performed.  \nThe x-ray pulse profiles emitted by hot spots at NS surfaces have also received considerable attention. Observations of these signals allow for a relatively clean inference of NS masses and radii (Watts et al. , 2016). This potential was met with observations (Bogdanov et al. , 2019a; Miller et al. , 2019; Riley et al. , 2019) by the Neutron Star Interior Composition Explorer (NICER) (Arzoumanian et al. , 2014). Sotani (2017) and Silva and Yunes (2019b) developed pulse-profile models for NSs in a massless DEF model. Xu et al. (2020) considered the case of massive scalar fields and Hu et al. (2021) studied the positive β 0 case. Silva and Yunes (2019a) made the first study of what constraints can be placed on the DEF model with pulse-profile observations and found they can be competitive with binary-pulsar observations. However, their work is rather simple and not at the level of realism found in the analysis of real data. This remains an important avenue for future work.  \nAt last, Tuna et al. (2022) used NS mass and radius measurements (Bogdanov et al. , 2016; Ozel et al. , 2016) to constrain the massive extensions of the DEF model. They obtained a weak lower bound β 0 ≳ -20 for scalar-field masses µ ≲ 2 × 10 -11 eV; see Fig. 10. This is significant since no other bound is known for µ ≫ 10 -16 eV. These results show that large scalar masses enable agreement with observations, even for extremely negative β 0 , and demonstrate the difficulty to constrain scalarization when the scalar field is massive.  \nb. Universal relations: One of the largest obstacles when using NS observations to test modified theories of gravity is the uncertainty in the NS EOS, which remains unknown at high densities (Baym et al. , 2018; Lattimer and Prakash, 2016). In general modifications to NS properties predicted by different EOSs are degenerate with changes to the underlying gravity theory used to model these stars. One way to break this degeneracy is to consider relations between different NS properties that depend weakly on the EOS (Doneva and Pappas, 2018; Yagi and Yunes, 2017). In fact, in Sec. III.B.1.b we already met one such  \n<!-- image -->  \nFigure 10 The posterior probability density for β 0 ∈ [ -100 , 0] for various assumed values of the scalar-field mass µ in a massive version of the DEF-model that was introduced by Ramazano˘glu and Pretorius (2016) (marginalized over various EOSs). One can obtain the bound β 0 ≳ -20 for µ ≲ 2 × 10 -11 eV (left panel), but there is no effective bound for higher masses, at least not in the interval β 0 ∈ [ -100 , 0] (right panel). No bound is possible for µ ≳ 10 -10 eV using these data since scalarization does not occur for such high µ values within the considered β 0 interval. From Tuna et al. (2022).  \n<!-- image -->  \nexample relating the QNM frequencies and damping times to the mean density of the NSs. Here we focus on other EOS-independent relations connecting various equilibrium properties of scalarized NSs.  \nA class of universal relations that has attracted considerable attention is the I-Love-Q relations, which connect the normalized moment of inertia, the tidal Love number, and the quadrupole moment of the NSs (Yagi and Yunes, 2013a,b). If any two quantities in the I-Love-Q trio are measured independently, one can constrain deviations from GR by determining whether these observations agree with the GR I-Love-Q relation within observational accuracy. Silva et al. (2021a) applied this strategy to constrain dynamical Chern-Simons gravity (Alexander and Yunes, 2009; Jackiw and Pi, 2003) using observational data from LVK and NICER. The application to scalarization models has not yet been done.  \nThe I-Love-Q relations for scalarized NSs were studied in slow (Pani and Berti, 2014) and rapid rotations (Doneva et al. , 2014a), including the case of massive scalar fields (Doneva and Yazadjiev, 2016; Hu et al. , 2021). In Fig. 11 we show the relation between the normalized moment of inertia ¯ I and the quadrupole moment ¯ Q for sequences of scalarized NSs with a fixed rotational parameter α = M Ω / 2 π . We see that for small absolute values of β , for example, β 0 = -4 . 5, the deviations from GR are practically negligible for slow rotation (such as α = 0 . 32) while they are enhanced by rapid rotation. As we have seen, massive scalar fields allow for much smaller values β 0 while remaining consistent with observations. Hence, these scalarization models allow for larger differences be-  \n/s32  \nFigure 11 Comparison between ¯ I -¯ Q relations for GR and scalar-tensor theory for several values of β 0 , where ¯ I = I/M 3 and ¯ Q = QM/J 2 . The presented data are restricted to APR EOS in order to provide better visibility, but the deviations for other EOSs are small, typically below roughly 2%. Sequences for different values of the normalized rotational parameter α = M Ω / 2 π are given, with Ω the angular velocity of the star. The dash-dotted lines with β = -4 . 5 appear for each value of α , but they are limited to the low ¯ Q region and deviate only slightly from GR. From Doneva et al. (2014a)  \n<!-- image -->  \n.  \ntween the scalarized NS universal relations and the GR ones (Doneva and Yazadjiev, 2016).  \nThe moment of inertia I and the quadrupole moment Q are among the leading-order multipole moments in the asymptotic expansion of the metric functions at infinity. There are infinitely many other higher-order multipole moments, however, and it is interesting to explore whether similar universal relations hold for them. This was first addressed in GR (Pappas and Apostolatos, 2014; Yagi et al. , 2014), where such universal relations were derived for the higher multipole moments when proper normalization was applied. Pappas and Sotiriou (2015) and Pappas et al. (2019) extended these results to scalarized NSs. This required the generalization of the multipole moment formalism of Geroch (1970) and Hansen (1974) to scalartensor gravity. Pappas et al. (2019) showed that future observations of QPOs of low-mass x-ray binaries can in principle be used to measure different NS properties and distinguish different gravity theories.  \nAnother class of universal relations studied in the context of the DEF model connects the NS moment of inertia and compactness. This was endorsed by Lattimer and Schutz (2005) [see also Breu and Rezzolla (2016)], who used this relation to argue that pulsar-timing observations could lead to a measurement of the moment of inertia to within 10%. These relations were generalized to scalarized NSs with massless scalar-field potential by Altaha Mo-  \n/s32  \ntahar et al. (2017), and later to massive self-interacting scalar fields by Popchev et al. (2019). As with the I-LoveQ relations, only the massive scalar-field case leads to large deviations from GR when observational constraints are taken into account.  \nOfengeim (2020) proposed connecting the physical parameters of static NSs, such as mass, radius, central energy density, pressure, and sound speed, at the maximum-mass point for a given EOS. This is based on the observations that the nuclear matter EOS (without phase transitions) can be well parametrized by only two parameters (Lindblom, 2010). As a result, Ofengeim (2020) derived multiple constraints on the nuclear matter EOS. Danchev and Doneva (2021) generalized these relations to scalarized NSs. They showed how sensitive these relations can be to the underlying theory of gravity. Thus, all of the EOS restrictions derived in GR should be taken with care since the possibility for GR modifications is rarely taken into account when interpreting observational data.', "D. Extended scalar-tensor theories beyond the Damour-Esposito-Far'ese model": "Thus far we have discussed scalarization in the DEF model or other DEF-inspired models. However, as we saw in Sec. II.B, spontaneous scalarization is also possible for other modified theories of gravity. Here we discuss the theories where scalarized NS solutions were obtained.  \na. Scalar-tensor theories with disformal coupling Minamitsuji and Silva (2016) studied NS scalarization in scalartensor theories with disformal coupling of the form of Eq. (17). More specifically, instead of the standard conformal transformation between the Einstein-frame metric ( g µν ) and the Jordan-frame metric (˜ g µν )  \ng µν = A 2 ( φ )˜ g µν , (33)  \nwhere A ( φ ) is a function of the scalar field related to the Jordan-frame coupling between the scalar field and the Ricci scalar, we have a more general and complicated transformation that also involves scalar-field derivatives,  \ng µν = A 2 ( φ ) [ ˜ g µν +Λ B 2 ( φ ) ∇ µ φ ∇ ν φ ] , (34)  \nwhere B ( φ ) is a function of the scalar field and Λ is a parameter of dimensions of [length] 2 . The motivation behind this modification of the conformal factor comes from Bekenstein (1993) who aimed to find the most general coupling constructed from the metric g µν and the scalar field φ that respected causality and the WEP. More recently Andreou et al. (2019) showed that such disformal coupling is actually equivalent to a kinetic coupling with the Ricci scalar and thus is contained in the minimal action (14); see Sec. II.B.1 for a discussion.  \nMinamitsuji and Silva (2016) showed that, for negative values of the disformal coupling parameter Λ, scalarization can be suppressed, while for large positive values of Λ the stellar structure equation becomes singular. Thus, regular NS solutions cannot be found in the latter case, which can be used to impose the upper limit Λ ≲ 100 km 2 . They also explored the universal relations between the moment of inertia and the compactness of NSs, and they determined the range of parameters where these relations can deviate significantly from GR.  \nb. Scalar-Gauss-Bonnet theory: Another well-studied modified gravity theory that allows for NS scalarization is the scalar-Gauss-Bonnet gravity discussed in Sec. II.B. Spacetime curvature itself, as well as matter, can be the source of the scalar field in this theory. Scalarized NS solutions in scalar-Gauss-Bonnet gravity were studied for the first time by Doneva and Yazadjiev (2018a) and Silva et al. (2018). The solutions exhibit much different features than to the DEF model while being qualitatively similar to NS solutions in other classes of Gauss-Bonnet theories where scalarization is not possible (Kleihaus et al. , 2016; Pani et al. , 2011a; Saffer et al. , 2019). For example, for a given central energy density the stellar mass is always smaller than in GR. Moreover, there is normally only one bifurcation point at small central energy densities, and afterward the branches of solutions are terminated because of violation of the regularity conditions.  \nIn Fig. 12, we show the branches of scalarized NS solutions for a coupling function f ( φ ) = -λ 2 (2 β ) -1 [1 -exp( -βφ 2 / 4)]. Here λ is a parameter with dimensions of [length] -1 that controls the coupling strength in the action (25). Note that the sign of f ( φ ) is the opposite of that discussed in Sec. II.B.3, and thus cannot lead to BH spontaneous scalarization, since in vacuum µ 2 eff > 0, as defined in Eq. (27). However, owing the presence of matter, NSs can scalarize. Scalarized NS solutions were also found for the more 'conventional' sign of f ( φ ), namely, f ( φ ) = λ 2 (2 β ) -1 [1 -exp( -βφ 2 / 4)] (Doneva and Yazadjiev, 2018a). As in the DEF model, it was demonstrated that binary pulsars strongly constrain the coupling parameters of the theory, thereby leaving a small window for scalarization (Danchev et al. , 2022).  \nLittle is currently known about the astrophysical implications of these stars. Kuan et al. (2021a) considered the spherically symmetric core collapse of a noncompact star either to a protoneutron star or to a BH in Gauss-Bonnet theory. They also proposed a realistic physical mechanism for the formation of isolated scalarized BHs and NSs The complexity of the problem is greatly increased with respect to GR, though, and there are fundamental difficulties such as the loss of hyperbolicity of the evolution equations for certain ranges of parameters, effectively limiting the maximum possible scalar field to relatively low values; see also East and Ripley (2021a,b) and Ripley and  \nM  \n/  \nM  \nFigure 12 Scalarized NSs in scalar-Gauss-Bonnet gravity. We show the mass as a function of the radius for GR (trivial branch with φ = 0) as well as several scalarized branches with different values of the scalar-field coupling parameters λ and β . From Doneva and Yazadjiev (2018a).  \n<!-- image -->  \nPretorius (2020) for the case of BHs. Core-collapse simulations show that the remnant in scalar-Gauss-Bonnet theory can be rich, with (de)scalarization happening at the intermediate or final stages of the collapse, depending on the properties of the progenitor and the theory parameters. Since breathing modes are absent from this theory, the effect on GW emission can be estimated only if one drops the assumption of spherical symmetry. This is something that has not yet been done in any modified gravity theories due to the complexity of the problem.  \nc. Ricci-Gauss-Bonnet model: One can further modify the action in scalar-Gauss-Bonnet gravity to include additional terms. This is the case for the Ricci-Gauss-Bonnet model discussed in Sec. II.B.4. This model has advantages such as the possibility of reconciling scalarization with cosmology. Here we focus on the NS solutions within this theory that were considered by Ventagli et al. (2021). They conducted a thorough exploration of the theory parameter space in order to find the sectors where NS solutions exist. Since their scalar charge is nonzero, in these sectors one can put severe constraints on the theory based on binary-pulsar observations, as in the DEF model; see Sec. III.A.1.b.  \nThe free parameters in the theory, as evident from the action (28), are β and α which control the coupling strength to the Ricci scalar and the Gauss-Bonnet invariant, respectively. In Fig. 13 the existence of NS solutions is shown in a two-dimensional plot spanning the parameters α and β in a broad range for a fixed central energy density and the MPA1 EOS (Muther et al. , 1987). The white area Fig. 13 corresponds to the region of the parameter space where the GR solution is stable against  \n2.5  \nFigure 13 Regions of existence of scalarized solutions in the ( α, β ) plane for the MPA1 EOS with central energy density ε c = 6 . 3 × 10 17 kg / m 3 . In GR, a star with this choice of ε c and EOS is light, with M GR = 1 . 12 M ⊙ . The wedgelike red (wedgelike and upward-curving blue) region is the region where scalarized solutions with 0 (1) scalar-field nodes exist. The gray contours obtained by Ventagli et al. (2020) are superimposed and represent the lines beyond which GR solutions with the same density are unstable to scalar perturbations with 0 , 1 , 2 , etc., nodes. We see that the region where there are scalarized solutions with n nodes is included in the region where the GR solutions are unstable to scalar perturbations with n nodes but much smaller. The dashed boundary of the blue region corresponds to a breakdown of the integration inside the star. From Ventagli et al. (2021).  \n<!-- image -->  \nscalar perturbations. Throughout the parameter space that is spanned by α and β , a new unstable mode appears every time that one crosses a black line. These unstable modes can be labeled by the number of scalar-field nodes (denoted by n = 0 , 1 , 2 and 3 in the figure). Any point in the parameter space that lies within a gray region corresponds to a configuration where the GR neutron star solution is unstable. The red (blue) area corresponds to the region where scalarized solutions with n = 0 ( n = 1) nodes exist.  \nIt is evident that scalarized solutions exist in only part of the parameter space. In the β > 0 and smallα range (the white region in Fig. 13) the GR neutron stars are stable. This is the region that is most noteworthy from the point of view of BH scalarization and where we can reconcile scalarization with cosmology (Antoniou et al. , 2021a,b); see also Sec. IV.A.2. Thus, including a Ricci coupling to the scalar-Gauss-Bonnet action seems to allow for BH scalarization while evading the binary-pulsar constraints.  \nFigure 13 is for only one EOS and a specific central energy density. Different EOSs and central densities can lead to significant deformations of the existence and stability regions. The general point made by Ventagli et al. (2021), however, is that the stability (white) region always survives.  \nd. Tensor-multiscalar theories: Another class of gravity theories where NS scalarization was considered is the tensor-multiscalar theories (TMST)s (Damour and Esposito-Far'ese, 1992), whose basics were discussed in Sec. II.B. In this theory, the gravitational interaction is mediated by the spacetime metric g µν and N scalar fields φ a , which take values in a coordinate patch of an N -dimensional Riemannian target manifold E N with a positive-definite metric γ ab ( φ ) defined on it (Damour and Esposito-Far'ese, 1992; Horbatsch et al. , 2015).  \nThe main features of TMST are the inclusion of more than one scalar field and a structure called the targetspace metric. In that sense, there are two directions to go in order to obtain scalarized NSs. The first one is to consider a mixture of several more or less equivalent scalar fields; this was the approach followed by Horbatsch et al. (2015). They, the authors examined the case of two scalar fields in the form of a complex scalar with a maximally symmetric target-space metric γ ab . Doneva and Yazadjiev (2020a) studied scalarization in tensor-multiscalar gravity in a more complicated setup, namely, when γ ab is a threedimensional maximally symmetric space together with a nontrivial map φ : spacetime → target space. While the solutions given by Horbatsch et al. (2015) can be viewed as a generalization of the DEF model to multiple scalar fields, the scalarized NSs given by Doneva and Yazadjiev (2020a) have some distinct properties and they can be considered more as a limiting case of the topological NSs discovered by Doneva and Yazadjiev (2020b). More specifically the scalar field is zero in the stellar center, while only its first derivative is zero in the DEF model. The scalar charge for the scalarized stars in TMST is zero as well. This automatically reconciles this theory with the binary-pulsar observations (due to the absence of scalar-dipole radiation) while still allowing for large deviations from GR. Another interesting property is the fact that there is nonuniqueness of the scalarized solutions with respect to the central energy density, something that had never been observed before, at least for scalarized NSs. Figure 14 illustrates this property. Interestingly, Kuan et al. (2021b) showed that all of these solutions are stable until the maximum NS, mass regardless of the nonuniqueness.  \ne. Other scalarization models: Apart from the previously discussed models of scalarized NSs, other works have attempted to consider scalarization in more exotic theories of gravity, or even different types of scalarization.  \nFor example, Azri and Nasri (2021) considered NS scalarization in what is called scalar-connection gravity, where gravity is mediated by a scalar field and the connection. Other types of scalarization can happen in standard scalar-tensor theories if we allow for different scalar-field potentials. Minamitsuji and Tsujikawa (2023) showed this when considering the poten-  \nM  \n/  \nM  \n3  \nFigure 14 Tensor-multiscalar theories: the mass as a function of the central energy density (left panel) and as a function of the stellar radius (right panel) for scalarized NSs in tensormultiscalar theories where the target-space metric represents a three-dimensional maximally symmetric space and we have a nontrivial map φ : spacetime → target space. Sequences of scalarized solutions for different combinations of parameters are given. Observe the appearance of nonuniqueness of the solutions with a nonzero scalar field for certain combinations of the parameters (for example, the small central energydensity region indicated by the dotted line). From Doneva and Yazadjiev (2020a).  \n<!-- image -->  \ntial V ( φ ) = m 2 φ f 2 B [1 + cos( φ/f B )] of a pseudo-NambuGoldstone boson. Here ( µ φ , f B ) are constants with dimensions of length and the scalar field has a negative mass µ 2 = -m 2 φ . Minamitsuji and Tsujikawa (2023) showed that in this case the scalar field sits at its vacuum expectation value far from the source, while inside a NS a symmetry restoration can take place, resulting in a new type of scalarization. This theory also has an advantage in that it the avoids cosmological instabilities present in certain other scalarization models.", 'IV. BLACK-HOLE SCALARIZATION': 'In this section, we discuss the spontaneous scalarization of BHs in three parts. To begin, in Sec. IV.A we consider models in which BHs can scalarize in vacuum due to couplings between the scalar field and curvature scalars. Next, in Sec. IV.B we discuss models in which scalarization is induced by the presence of extra fields (such as gauge or matter fields) in the BH spacetime. Finally, in Sec. IV.C we review a selection of other models of BH scalarization.', 'A. Black-hole scalarization: Vacuum spacetimes': 'As we saw in Sec. II, some gravity theories admit the same vacuum BH solutions as GR, yet can give rise to new branches of solutions with scalar hair once certain  \nconditions are met. The prototypical example is described by scalar-Gauss-Bonnet theories, whose action in the absence of matter is given by Eq. (25), namely,  \nS = 1 16 πG ∫ d 4 x √ -g [ R -1 2 g µν ∂ ν φ∂ µ φ + f ( φ ) G ] .  \nAs discussed in Sec. I, the first models exhibiting BH scalarization in vacuum were proposed by Doneva and Yazadjiev (2018b) and Silva et al. (2018). These models remain the most studied BH scalarization models to date in the literature. For this reason, they will be our main focus in this section. 13  \nAs explained in Sec. II.B.3, in theories described by the action (25), a no-hair theorem guarantees that the BH solutions of GR are unique to the theory, as long as (d 2 f/ d φ 2 ) φ = φ 0 G < 0, for some constant φ 0 (Silva et al. , 2018). If this condition is violated, scalar-field perturbations can become tachyonic unstable and the endpoint is expected to be a nonlinear, scalarized BH (Ripley and Pretorius, 2020). Scalarized BHs have been shown to form dynamically, as outcomes from the core collapse of an initially unscalarized star (Kuan et al. , 2021a).  \nHere we first review our understanding of isolated scalarized BHs in scalar-Gauss-Bonnet theories, focusing on their properties and their stability. We then review what is understood about when they are found in binaries. Finally, we give an overview of vacuum BH solutions in models that generalize the action (25).', '1. Scalarized black holes': "Spontaneously scalarized BH solutions were first found by Doneva and Yazadjiev (2018b) and Silva et al. (2018). What distinguishes these works is the choice of the coupling function f ( φ ) that couples the scalar field to the Gauss-Bonnet invariant [cf. Eq. (25)],  \nf = ( λ 2 / 12) [1 -ε exp( -3 φ 2 / 2)] , (35a)  \nf = ( η/ 8) ε φ 2 , (35b)  \nwhere λ 2 and η , positive by definition, are coupling constants with dimensions of [length] 2 and ε = ± 1, which should not be confused with the energy density of Sec. III. The two coupling functions agree in the smallφ limit (i.e., when φ ≪ 1) and therefore result in the same prediction for the onset of scalarization of GR BHs. This threshold can be found by searching for bound state solutions, i.e., time-independent solutions of the linearized field equation  \nFigure 15 Results of the numerical integration of the equation of a monopole ( ℓ = m = 0) scalar perturbation δφ on a Schwarzschild BH background. Top panel: asymptotic value of | d σ/ d r | (with σ ≡ r δφ ) evaluated at a point r max ≫ r h , far from the event horizon r h = 2 M , as a function of η/M 2 . Cusps mark when bound state solutions form. They signal the transition between stable Schwarzschild BHs and scalarized solutions. Bottom-left panel: effective potential V eff (see Sec. II.A.2) for ( η/M ) 2 = 0 and 5. In the latter case, V eff develops a negative region and can support bound states. Bottom-right panel: radial profiles of δφ for the first three bound states corresponding to η/M 2 = 2 . 902, 19 . 50, and 50 . 93. These profiles have 0, 1, and 2 nodes, respectively. The properties of these solutions are shared between the models of Doneva and Yazadjiev (2018b) and Silva et al. (2018). From Silva et al. (2018).  \n<!-- image -->  \nfor φ around a GR BH, which are regular at the event horizon r h and that vanish at spatial infinity. Hence, the determination of the scalarization threshold reduces to a boundary value problem. In the simplest case of a Schwarzschild BH of mass M , the dimensionless quantity η/M 2 plays the role of the eigenvalue and can be determined numerically with standard techniques such as the shooting method. The smallest eigenvalue gives the threshold for the formation of the 'fundamental' scalarized BH solution (in the sense that the radial profile of the scalar field has n = 0 nodes). The other eigenvalues give the thresholds for the formation of 'excited states,' that is, solutions with one or more nodes. In Fig. 15 we show the results of such a calculation in the spacetime of a Schwarzschild BH and ε = 1, as done by Silva et al. (2018). We remark that no bound states can form for ε = -1, because the effective mass of scalar-field perturbations is positive; therefore, the effective potential is positive definite. See the discussion in Sec. II.A.2.  \nWhile both models in Eqs. (35) agree in their prediction of the threshold for scalarization, they differ significantly in their prediction of the properties of the nonlinear, scalarized solution, as discussed in Sec. II.A. In Fig. 16 we show the branches of scalarized BHs, in a parameter space spanned by the dimensionless scalar charge and BH mass. The solutions were obtained by solving the full system of  \n/s108  \n/  \nD  \nFigure 17 Rotating scalarized BHs in the Gaussian model (35a). The shaded regions represent the domain of existence of scalarized rotating BHs in the plane spanned by the dimensionless scalar charge ( Q s /λ ) and spin ( j = a/M ). These regions are bound by the Kerr family of solutions Q s /λ = 0 (solid blue line) and, in the nonrotating limit, by the Schwarzschild solution (black dot-dashed line), and also when the regularity condition for the scalar field is violated (dotted curve). The inset shows the case for positive coupling constant studied by Cunha et al. (2019), while the main panel shows nonlinear spin-induced scalarized BHs. From Herdeiro et al. (2021b).  \n<!-- image -->  \nFigure 16 Scalar charge D normalized to the coupling parameter λ as a function of the normalized BH mass M/λ for sequences of scalarized BH solutions. The GR solutions have zero scalar field and thus zero D and lie along the x axis, which is depicted by the horizontal magenta line. The dashed line corresponds to solutions with nodeless scalar-field profiles, which was found to be the only stable branch (Bl'azquez-Salcedo et al. , 2018). The curves with different colors correspond to scalarized solutions with scalar fields having different numbers of nodes. Solutions with one and two nodes branch off the GR solution around M/λ ≃ 0 . 15 and M/λ ≃ 0 . 24, respectively. The scalarized branches are symmetric with respect to the y axis due to the reflection symmetry ( φ →-φ ) in the theory. Adapted from Doneva and Yazadjiev (2018b).  \n<!-- image -->  \nfield equations; see Doneva and Yazadjiev (2018b) and Juli'e et al. (2022) for detailed discussions. The results are for the Gaussian coupling (35a) with ε = 1. We show the fundamental, zero-node scalar-field solution (dashed red curves) and also the excited scalarized BHs (solid green and blue curves). We see that the branches of excited BHs can terminate at finite masses M . This is due to a violation of an inequality at the BH horizon that must be satisfied for the scalar field to be real valued and is typical of scalar-Gauss-Bonnet theories; see e.g., Kanti et al. (1996) and Antoniou et al. (2018). The parameters in the coupling function (35a) are chosen in such a way that this disappearance of BH solutions does not affect the fundamental branch. For other couplings, such as Eq. (35b), even this branch can violate the regularity condition shortly after bifurcation; see Doneva et al. (2018a) and Silva et al. (2019). Hence, the domain of existence of scalarized solutions in the Gaussian model (35a) is larger with respect to that of the 'quadratic' model (35b).  \nAsymptotically flat, spinning scalarized BHs in the Gaussian models (with ε = 1) and quadratic theories were explored, respectively, by Cunha et al. (2019) and Collodel et al. (2020). In discussing these solutions, we consider a plane spanned by the dimensionless scalar charge Q s /λ and dimensionless spin parameter j ≡ J/M 2 , where J is the angular momentum of the BH. Figure 17 shows the domain of existence of scalarized rotating BH solutions. In particular, we focus on the inset, which corresponds to  \nε = 1 in Eq. (35a). We see that as the spin increases the existence domain of the solutions (the shaded region) becomes smaller. This result can be understood in terms of the Gauss-Bonnet invariant not being positive definite for the Kerr metric (Cherubini et al. , 2002). More specifically, a Kerr BH develops increasingly large regions where the Gauss-Bonnet invariant is negative as the spin increases and hence suppresses scalarization (or triggers scalarization, when ε = -1, as discussed in Sec. IV.A.3.) As a consequence, the deviations relative to GR predictions on physical observables (for example, the location of the innermost stable circular orbit of massive particles or the BH shadow) are smaller for rapidly rotating BHs. See Fig. 18 for an example. Collodel et al. (2020) studied the quadratic coupling (35b) and reached similar conclusions. They obtained radially excited solutions (i.e., scalar-field profiles with nodes in the radial direction) and also notes the existence of 'angularly' excited rotating BHs (i.e., scalar-field profiles with nodes in the angular directions).", '2. Stability of scalarized black holes and implications for model building': "Having established the existence of scalarized BHs, the natural next task is to study their stability. A first indication comes from the study of the thermodynamical properties of such BHs. Doneva and Yazadjiev (2018b) found that the scalarized solutions belonging to the fundamental branch in the Gaussian model are thermody-  \nFigure 18 Gravitational lensing and shadow produced by rotating BHs under similar observation conditions. Left panel: scalarized BH with mass M/λ ≈ 0 . 237 and spin j = 0 . 24. Right panel: a comparable Kerr BH. From Cunha et al. (2019).  \n<!-- image -->  \nFigure 19 Scalar-charge-mass diagram for scalarized BHs with a quadratic scalar-Gauss-Bonnet model with a scalarfield potential V ( φ ) = ( µ 2 / 2) φ 2 + ( λ/ 2) φ 4 . Dimensionless quantities obtained by factors of η 1 / 2 are used: BH charge ˆ Q = Q/η 1 / 2 and mass ˆ M/η 1 / 2 , and scalar-field parameters ˆ µ = µη 1 / 2 and ˆ λ = λη 1 / 2 . The vertical line marks the threshold for scalarization ˆ M th . BHs with ˆ M ⩾ ˆ M th are radially unstable, while the Schwarzschild solution is stable. For large ˆ λ , we can obtain solutions with ˆ M < ˆ M th and form two branches. The ones that are radially unstable solutions are denoted with a dashed line (upper branch) and those stable to radial perturbations are denoted with a solid line (lower branch). The dot marks are the marginally stable solutions, which correspond to a minimum mass but maximally scalar charged BHs for fixed ˆ µ and ˆ λ . From Macedo et al. (2019).  \n<!-- image -->  \nnamically favored relative to the Kerr solution (35a), in the sense that they have a smaller Wald entropy (Wald, 1993). This suggests the stability of the former; i.e., in a dynamical process BHs in this model would evolve toward a state of smaller entropy. In contrast, the quadratic coupling (35b) model always results in scalarized BHs with higher entropy relative to their GR counterparts. These observations remain true in the case of rotating BHs (Collodel et al. , 2020; Herdeiro et al. , 2021b).  \nA complementary question is whether scalarized BHs  \nFigure 20 Like Fig. 19, but now in a 'mixed' model where a quadratic coupling between the scalar field and the Ricci tensor is added to the quadratic scalar-Gauss-Bonnet interaction. As in Fig. 19, we show the charge ˆ Q = Q/α 1 / 2 and mass ˆ M/α 1 / 2 made dimensionless by division by α 1 / 2 , for different values of scalar-Ricci coupling β . See Eq. (28) for the theory's action. The left-bending black lines correspond to values of β for which all scalarized BHs have masses below the GR instability mass threshold, while the right-bending red lines mark values of β that lead to scalarized BH masses that are larger than the GR mass threshold. Past the turning points (marked with blue dots), the BHs are expected to be unstable to radial perturbations. From Antoniou et al. (2021b).  \n<!-- image -->  \nare stable to linear perturbations. The first step in this direction was carried out by Bl'azquez-Salcedo et al. (2018), who considered radial scalar-field and metric perturbations. They showed that Schwarzschild BHs became unstable at the scalarization threshold and that the stability of the fundamental (i.e., the zero-node) scalarized solution depended on which model in Eq. (35) was considered: in the quadratic case all solutions were unstable, whereas stable solutions existed in the Gaussian model.  \nThis important difference between the two models can be traced back to the higher-than-quadratic terms in the scalar field in the theory's action. More concretely Silva et al. (2019) and Minamitsuji and Ikeda (2019a) showed that the inclusion of a quartic interaction, say, ( ζ/ 16) φ 4 G , is sufficient to yield stable BH solutions as long as ζ is negative and satisfies ζ/η < -0 . 8. These two conditions are satisfied by the Gaussian model. Indeed, expanding Eqs. (35a) to order O ( φ 4 ), we find that ζ/η = -3 / 2, which helps to explain the stability of BHs in this model.  \nIt was later shown by Macedo et al. (2019) [see also Macedo (2020)] that the original quadratic model can be made stable under radial perturbations by making the scalar field massive and self-interacting, V ( φ ) ∝ µ 2 φ 2 + λφ 4 , as discussed in Fig. 19. From an EFT perspective, these terms are of lower order than the quartic scalar-Gauss-Bonnet interaction and hence should be included. A similar analysis was also performed for the Ricci-Gauss-Bonnet model by Antoniou et al. (2021b) and is shown in Fig. 20. See also Doneva et al. (2019) for a study of scalarization with massive scalar fields.  \nThe nonlinear stability of BHs in this model was also studied in the time domain by Ripley and Pretorius (2020). They found evidence for regions in mass-coupling parameter space in which the end state of the radial instability of the Schwarzschild BH is a stable scalarized BH. For larger couplings however, regions where the time-evolution equations change character from hyperbolic to elliptic appear outside the BH horizon. This signals a regime in which the theory does not have a well-posed initial-value problem. See Hilditch (2013) and Ripley (2022) for discussions. These results are in agreement with earlier findings by Bl'azquez-Salcedo et al. (2018), but in the linear regime. The radial stability of scalarized BHs in the Ricci-Gauss-Bonnet model introduced in Sec. II.B.4 was analyzed by Antoniou et al. (2022). They found that although the Ricci scalar does not affect the scalarization threshold, at the nonlinear level it can stabilize the BHs in the quadratic model.  \nThe problem of nonradial perturbations of scalarized BH solutions, which is relevant in the context of testing such theories with ringdown GW signals, was addressed by Bl'azquez-Salcedo et al. (2020a,b). However, the stability of rotating scalarized BHs remains an open problem, even in the slow rotating limit.", '3. Spin-induced scalarization': 'One property of the Gauss-Bonnet invariant of the Schwarzschild metric is that it is positive definite everywhere and, as a consequence, spontaneous scalarization can occur only for positive values of the coupling constant. Does this always have to be the case? In the case of the Kerr metric, it is known that regions where G < 0 outside the outer horizon are possible if the BH spins with a/M ⩾ 0 . 5 (Cherubini et al. , 2002), suggesting the possibility of a spin-induced scalarization when the coupling constant is negative, i.e., ε = -1 in Eqs. (35).  \nAlthough conceptually simple, the problem is not straightforward to analyze because (i) the analytical form of G forbids the separation of variables as done in Sec. I, and thus the determination of the scalarization threshold will in principle require more sophisticated numerical methods than the Schwarzschild case, and (ii) the Kerr metric (with a/M ⩾ 0 . 5) has regions where G can be either negative or positive, where in the latter regions, for negative coupling constants, the effective mass may trigger a superradiant instability (Brito et al. , 2015). To overcome these difficulties, Dima et al. (2020) reduced the (2+1)-dimensional nonseparable scalar-field equation to a (1 + 1)-dimensional system of equations that are coupled through spherical harmonic multipole indices l and m . They then performed time-domain numerical integration of the scalar field borrowing methods that were previously developed to study superradiant instability (Dolan, 2013). Dima et al. (2020) showed that the tachyonic instability  \nFigure 21 Instability timescale τ for spin-induced scalarization in the dimensionless spin a/M and dimensionless coupling constant -η/M 2 plane. The instability threshold for the total reconstructed field is shown as a solid green line, while the threshold when the m = 0 modes are excluded is shown as a blue dotted line. The red dashed line corresponds to the instability threshold for the m = 0 odd modes, while the dot-dashed cyan line marks the instability threshold for the spherical mode ℓ = m = 0. From Dima et al. (2020).  \n<!-- image -->  \nis the dominant one and charted the parameter space in which spin-induced scalarization would occur. This is shown in Fig. 21. Hod (2020) worked analytically in the infinitely large coupling limit ( η/M 2 →∞ ) and confirmed the expectation that a/M ≳ 0 . 5 is the minimal necessary spin value for which spin-induced scalarization occurs. The corrections for large but finite coupling were obtained by Hod (2022).  \nNonlinear spin-induced scalarized BH solutions were obtained by Herdeiro et al. (2021b) and by Berti et al. (2021). They confirmed the existence of scalarized BH solutions in the parameter space region in which linear theory predicts the spin-induced tachyonic instability. In Fig. 17 we show the domain of the spin-induced scalarized BHs in the Gaussian model (35a). We see that scalarized solutions exist only for j ≳ 0 . 5, and that their scalar charge increases with spin. Moreover, there are solutions that can violate the Kerr bound j ≤ 1 on BH spins.  \nA approach complementary to that of Dima et al. (2020) was used by Doneva et al. (2020a) and Doneva and Yazadjiev (2021a). They evolved the scalar-field equation in 2+1 dimensions, i.e., without doing a multipolar decomposition. The influence of a nonvanishing scalar-field mass was explored by Doneva et al. (2020b). A similar timedomain study that instead used the hyperboloidal foliation method (Zenginoglu, 2008) was performed by Zhang et al. (2020) and a model with both the Pontryagin and Gauss-Bonnet invariants was studied by Myung and Zou  \n(2021). Annulli et al. (2022) studied the effect of strong magnetic fields on spin-induced scalarization. They found that the magnetic field suppressed this effect; i.e., it pushed the scalarization threshold to larger values of the dimensionless spin j of the BH. The case of spin-induced scalarization has been much less studied in terms on its stability and on how additional nonlinear interactions affect the scalar charge.', '4. Scalarized black holes in binary systems': "Most work in BH scalarization has focused on isolated BHs, but to be able to confront these models against GW observations one has to turn to the two-body problem. What phenomenology would we expect in BH binaries? To answer this question we first consider work that explored the nonlinear regime of the late inspiral and merger of BH binaries, a regime that relies on numerical relativity . Next we consider work that focused on the inspiral alone, a regime that can be modeled with PN theory.  \nThe first work in this area was done by Silva et al. (2021b), who studied the scalar-field dynamics (i.e., in the decoupling limit) in head-on collisions and nonspinning, quasicircular inspirals of binary BH spacetimes obtained from numerical relativity simulations. This work adapted the methods developed by Witek et al. (2019) for shiftsymmetric scalar-Gauss-Bonnet theory to theories which that allow for spontaneous scalarization. In particular, because they were interested at phenomenology near the scalarization threshold, Silva et al. (2021b) adopted the quadratic model (35b), with a coupling strength η such that both (or one of them, depending on the mass ratio) binary components or the remnant BH can carry a scalarfield bound state through the simulation. They showed that BH binaries can either form a scalarized remnant or dynamically descalarize by shedding off the initial scalar hair (i.e., the scalar bound state configuration) depending on the values of the coupling constants and the mass ratio between the two BHs. An example of dynamical descalarization is shown in Fig. 22. Dynamical descalarization was also shown to occur in nonlinear scalarization models (see Sec. II.C.5) by Doneva et al. (2022) that simulated head-on BH collisions and worked in the decoupling limit.  \nA natural question that follows is: What would happen when ε = -1, i.e., the case in which spin-induced scalarization occurs? This is relevant from an observational point of view because BH remnants of binary BH coalescences have typical dimensionless spins of j ∼ 0 . 7. This value comfortably meets the criteria for spin-induced scalarization to happen. This question was explored by Elley et al. (2022), who found through a suite of numerical relativity simulations that: (i) spin-induced dynamical descalarization can happen when the remnant BH spin is smaller than that of the initially spinning (and scalarized) binary components, and (ii) spin-induced dynamical  \nscalarization in an initially nonspinning BH binary system can result in a scalarized rotating BH.  \nIn the latter case, if the value of the coupling constant is sufficiently large the scalar field may affect the late inspiral and ringdown of the binary. For values of the coupling constant near the scalarization threshold for the remnant BH, the tachyonic instability can be delayed to some time ≈ 100 M after the binary merger, and can thus hide the scalar field throughout the binary's inspiral. This was termed stealth scalarization by Elley et al. (2022); Fig. 23 displays an example.  \nA restriction of the decoupling limit is that it does not capture the backreaction of the scalar field onto the spacetime metric, and thus does not allow one to derive the modifications to the gravitational waves in the previously mentioned binary setups. These difficulties can be partially overcome by employing a formulation of the theory's field equations in a generalized harmonic gauge that was first developed by Kov'acs and Reall (2020a,b); see also Ripley (2022). This formulation in principle allows one to evolve the complete scalar-field and metric dynamics in scalar-Gauss-Bonnet gravity. However, the first binary BH coalescence simulations using this formalism in scalar-Gauss-Bonnet theories that allow for BH scalarization (East and Ripley, 2021a,b) demonstrated that evolution ceases to be well posed for larger values of the coupling. How different gauge choices and different models perform in this respect has not been well explored.  \nFranchini et al. (2022), motivated to go beyond the small-coupling approximation, explored the fixing-theequations approach to perform numerical simulations, an approach inspired by work in dissipative relativistic hydrodynamics (Cayuso et al. , 2017). The main idea is to modify the evolution equations of the theory such that the short wavelength modes (associated with strong coupling and the breakdown of well-posedness) with respect to some chosen scale are somewhat controlled, while longer wavelength modes become tractable. Franchini et al. (2022) were able to use the formalism to study scalar-field collapse in spherical symmetry in the scalarization model of Silva et al. (2019). This allowed them to evolve the system past situations where the original equations fail due to the appearance of regions in spacetime where the evolution changed character from hyperbolic to elliptic, which prevented evolution of the system in time (Ripley and Pretorius, 2020). Moreover, the 'fixed-equation' system evolves toward the static BHs predicted by the original system of equations, thus indicating that it may provide a suitable 'completion' of the original theory. Whether this approach remains valid in a situation with less symmetry, for example the case of BH binaries, remains to be explored.  \nProgress has also been made in modeling compact binary dynamics in the regime of wider separation and lower velocities, i.e., in the domain of validity of PN theory. In particular, Shiralilou et al. (2021, 2022) extended to a gen-  \n|  \n|  \nFigure 22 Dynamical descalarization in binary BH head-on collisions. Scalar-field and Gauss-Bonnet dynamics on the x -y plane of the BHs with initial mass ratio q = m 1 /m 2 = 1 / 2. We show the scalar-field amplitude log 10 | Φ | (color map) together with the isocurvature levels of the Gauss-Bonnet invariant at the beginning of the evolution (top-left panel), during the BH approach (top-right panel), shortly before the collision (bottom-left panel) and shortly after the merger (bottom-right panel). These levels correspond to 1 M -4 (solid line), 10 -1 M -4 (dashed line), 10 -2 M -4 (dot-dashed line) and 10 -3 M -4 (dotted line). From Silva et al. (2021b).  \n<!-- image -->  \nFigure 23 Snapshots of the scalar field, here labeled Φ, and the Gauss-Bonnet invariant G in the x -z plane llustrating the stealth scalarization. The color map indicates the amplitude of the scalar field. The curves represent isocurvature contours for the following values of | G | : 1 M -4 (solid line), 10 -1 M -4 (dashed line), 10 -2 M -4 (dot-dashed line), and 10 -3 M -4 (dotted line). Black (red) lines correspond to positive (negative) values of G . We show the inspiral (top-left panel), half an orbit before merger (top-right panel), formation of the first common apparent horizon (bottom-left panel), and at 200 M after the merger. From Elley et al. (2022).  \n<!-- image -->  \neral f ( φ ) coupling the calculations of Yagi et al. (2012), which applied for the shift-symmetric scalar-Gauss-Bonnet theory, i.e., f ( φ ) ∼ φ . These works calculated the leadingorder corrections due to curvature nonlinearities in the GW and scalar waveforms, finding that corrections due to the Gauss-Bonnet term appear at 1PN order in GWs. They also obtained the GW polarization and phasing. In  \naddition, Juli'e and Berti (2019) argued that during the adiabatic inspiral of two BHs the Wald entropy of each BH is constant. This allows a precise definition of the sensitivities of BHs in scalar-Gauss-Bonnet gravity in terms of the variation of the ADM mass as a function of the ambient scalar field that the BH is embedded in at fixed Wald entropy. This parallels the notion that the baryonic mass of NSs is constant during the inspiral (Damour and Esposito-Far'ese, 1996a). They also derived the two-body Lagrangian at 1PN order. Juli'e et al. (2022) used this result in conjunction with the numerical calculation of sensitivities of spontaneously scalarized BHs to study binaries. They found that in principle BHs can evolve toward a situation in which the inner singularity approaches the event horizon of the BH before the merger. This suggest the possibility that, depending on the model of scalarGauss-Bonnet gravity and the binary parameters, the field equations might simply not be able to predict a full binary evolution ranging from inspiral to merger.", 'B. Black-hole scalarization in the presence of matter': "Stefanov et al. (2008) suggested the earliest model of BH scalarization, which consisted of the DEF model coupled to nonlinear electromagnetism (Stefanov et al. , 2007a,b). They introduced a coupling between the scalar field and the Born-Infeld Lagrangian and found that BH solutions with scalar hair branch off from the generalrelativistic sequence of solutions. These BHs violate the no-scalar-hair theorems of Mayo and Bekenstein (1996) and Sen and Banerjee (2001), which are valid for charged, nonrotating BHs in scalar-tensor theory. The bifurcation point occurs where the GR sequence of solutions becomes unstable to scalar perturbations (Doneva et al. , 2010).  \nA variety of scalarization models have been studied in Einstein-Maxwell-scalar theory following Herdeiro et al. (2018), who introduced a coupling between the scalar  \nfield and the Lorentz invariant F 2 = F µν F µν of the form exp( -αφ 2 ) F 2 , where α is a dimensionless constant. In this model BHs can scalarize, with F 2 playing the role of the trace of the energy-momentum tensor in the NS scalarization in DEF-like models, or the Gauss-Bonnet invariant in BH scalarization models.  \nThe simpler nature of the model allowed Herdeiro et al. (2018) to perform numerical relativity simulations to study the scalarization process in the time domain by expanding upon already available code bases (Sanchis-Gual et al. , 2016a,b). This allowed Herdeiro et al. (2018) to show that the end point of the instability of the Reissner-Nordstrom BH was indeed a BH with scalar hair. Additional studies of this model were done by Astefanesei et al. (2019); Fernandes et al. (2019b); Herdeiro et al. (2021a); Luo et al. (2022); Niu et al. (2022); and Xiong et al. (2022). Within this model, Myung and Zou (2019b) studied the instability of the Reissner-Nordstrom solution, while Bl'azquezSalcedo et al. (2021); Myung and Zou (2018, 2019c) analyzed the stability of the scalarized solutions. Fernandes et al. (2019a) and Herdeiro and Oliveira (2020) explored other variations on this theme, including scalar-axion couplings, i.e., ∝ h ( φ ) F µν ˜ F µν , where ˜ F µν = (1 / 2) ϵ αβ µν F αβ . In addition, Ikeda et al. (2019) considered the coupling between the double dual of the Riemann tensor and F µν .  \nZhang (2021) studied an Einstein-Maxwell-scalar model, with a coupling between scalar and Maxwell fields of the form exp( βφ 4 ) F 2 . In this model, the Reissner-Nordstrom solution is linear stable against scalar perturbations, but it can become unstable under large, nonlinear perturbations. Zhang (2021) observed a critical phenomena separating unscalarized and scalarized BHs. This establishes an interesting connection between nonlinear scalarization with the result of critical phenomena in gravitational collapse, which was first observed by Choptuik (1993). Similar investigations were carried out in Einstein-Maxwellscalar the Reissner-Nordstrom-anti-de Sitter BHs (Zhang et al. , 2022) and in the Ricci-Gauss-Bonnet model (Liu et al. , 2023a,b).  \nCardoso et al. (2013a,b) explored a relevant scenario from an astrophysical perspective: matter in the vicinity of a BH could trigger scalarization and in turn endow the BH with scalar hair. They showed that this 'matterinduced' scalarization (see Sec. II.C.3) is in principle possible for the DEF model. The viability of this proposal in realistic matter configurations, such as an accretion disk or a dark matter halo, has not yet been explored. Matterinduced scalarization seems unlikely to happen in light of the severe constraints on the DEF model. However, this might not be the case for other scalarization models, and the topic deserves further investigation.  \nAfter it was understood that in scalar-Gauss-Bonnet theories the Schwarzschild solution can scalarize, it was natural to ask if the same can happen to its charged counterpart, i.e., the Reissner-Nordstrom solution. Doneva et al. (2018a) analyzed this and found the existence of two  \nbifurcation points, one at larger masses where the scalarized solutions bifurcated from the Reissner-Nordstrom one, and one at smaller masses where the scalar charge of the solution decreases again to zero and the branch merges again with the GR one. Scalarized charged BHs in the fundamental mode are also thermodynamically favored over the Reissner-Nordstrom solution. Brihaye and Hartmann (2019) showed that for the near-extremal ReissnerNordstrom case, scalarization can happen for either sign of the scalar-to-Gauss-Bonnet coupling constant (35b). Other works varying the choice of the coupling function f between the scalar field (or the axion field) and the Maxwell invariant and studying the scalarization of dyonic BHs include Fernandes et al. (2019b), Fernandes et al. (2019a), Bl'azquez-Salcedo et al. (2020c) and Bl'azquezSalcedo et al. (2021).  \nHerdeiro and Radu (2019) put forward their motivation for studying these models. They noted that quantum effects can break the scale invariance and vanishing energy-momentum trace properties of electro-vacuum GR BHs. As a consequence, even the simplest nonminimally scalar-field coupling ξφ 2 R can result in BH scalarization when these quantum effects are taken into account. Within these models, they discussed the scalarization of the Reissner-Nordstrom solution and a noncommutative generalization of the Schwarzschild solution. They found that the resulting scalarized BHs are in general entropically favored over the GR solutions.", 'C. Variations of the curvature-induced scalarization model': "In light of the nontrivial effect of rotation on scalarization, it is natural to consider what happens when the Gauss-Bonnet invariant is replaced by the Pontryagin invariant ∗ RR as the 'curvature source' to which the scalar field is couple. The Pontryagin density is known to vanish in spherically symmetric spacetimes (such as that of a Schwarzschild BH) and becomes nonzero in nonspherically symmetric spacetimes (such as that of a Kerr BH). In this sense, all scalarized BHs in theories that replace the Gauss-Bonnet invariant with the Pontryagin density are necessarily 'spin induced'; a characteristic shared with BH solutions in dynamical Chern-Simons gravity [see e.g., Konno et al. (2009) and Yunes and Pretorius (2009), and Alexander and Yunes (2009) for a review] which features a linear coupling between the scalar field and the Pontryagin invariant. The coupling between a scalar field and the Pontryagin density leads to equations that are of third order in derivatives [see Delsate et al. (2015) and Motohashi and Suyama (2012)], which is not the case for a coupling with the Gauss-Bonnet invariant. Thus, although both of these couplings can be seen as part of an EFT, they lead to distinct challenges when it comes to nonlinear evolution. Myung and Zou (2019a) studied the combined effect of Gauss-Bonnet and Pontryagin densities. In the  \ntest field limit, the scalar-field dynamics with an effective mass proportional to ϕ 2 ∗ RR was also studied in a Kerr background (Doneva and Yazadjiev, 2021b; Gao et al. , 2019) and in the Schwarzschild-Newman-Unti-Tamburino background (Brihaye et al. , 2019).  \nAnother extension of the original Gauss-Bonnet spontaneous scalarization model involves the consideration of n > 1 scalar fields φ a , whose interaction is determined by their 'target space' γ ab ( φ c ), an n -dimensional Riemannian manifold (Damour and Esposito-Far'ese, 1992; Horbatsch et al. , 2015). These are similar to the tensormultiscalar DEF models for NS scalarization presented in Sec. III.D. In these models, the scalar-field dynamics is determined by the quantity γ ab ( φ ) g ρσ ∂ ρ φ a ∂ σ φ b . Doneva et al. (2020c) numerically obtained a scalarized BH for the case n = 3 and maximally symmetric target-space geometries. An important feature of these solutions is that the scalar fields decay asymptotically as 1 /r 2 (i.e., the BHs do not have a monopole scalar charge). This shows that BHs in these theories will not emit dipole scalar radiation when they are placed in binaries.", 'V. GENERALIZATIONS OF SCALARIZATION TO OTHER FIELDS AND INSTABILITIES': 'We identified the underlying reason for spontaneous scalarization in its various forms to be a tachyonic instability. However, the scalar nature of the field did not play a special role in the mechanism. This suggests that other fields, such as vectors might also spontaneously develop nontrivial configurations around compact objects when they exhibit suitable couplings to curvature. In this section, we investigate this idea, which leads to the so-called spontaneous tensorization theories.  \nA second type of generalization of spontaneous scalarization considers new types of instabilities, as opposed to new types of fields. For instance, instead of replacing the scalar field with, say, a vector, we replace the tachyonic instability with, say, a ghost instability. We later see that a surprising key result of spontaneous tensorization is that these two paths are intimately connected. Namely, even if we intend only to have a theory of tachyons living on nonscalar fields, ghost and gradient instabilities necessarily arise in almost all models.', 'A. Spontaneous vectorization': 'What happens if we replace the scalar field of spontaneous scalarization with a vector? In scalarization, we started with the most general context in Sec. II.B.1, considering all the allowed couplings to the metric. Here we follow the reverse path, starting with individual models of vector-tensor theories and considering the more encompassing theory later. This exposition is preferable  \nbecause the study of specific models that are straightforward generalizations of known scalarization models reveals some pathologies and provides guidance for further model building.  \nThe idea and the first specific model of vectorization were introduced in analogy with a massive version of the DEF model in Eq. (18). Consider the vector-tensor theory action (Beltran Jimenez et al. , 2013; Ramazano˘glu, 2017)  \nS = 1 16 πG ∗ ∫ d 4 x √ -g [ R -F µν F µν -2 µ 2 X X 2 ] + S m [ Ψ m ; A 2 ( X 2 ) g µν ] , (36)  \nwhere F µν = ∇ µ X ν - ∇ ν X µ and X 2 = g µν X µ X ν = X µ X µ . The vector field X µ has the canonical kinetic term in this frame, and the matter degrees of freedom couple to the metric ˜ g µν = A 2 ( X 2 ) g µν , which is conformally scaled with respect to g µν . The vector field equation is  \n∇ ρ F ρµ = [ -8 πG ∗ A 4 Λ ˜ T + µ 2 X ] X µ , (37)  \nwhere Λ = d ln A / d x , T µν is the energy-momentum tensor in the Einstein frame, and ˜ T is the trace of the stressenergy tensor in the Jordan frame, i.e., with respect to the metric ˜ g µν .  \nEquation (37) is that of a massive vector (Proca) field where the expression inside the square brackets acts as an effective mass µ 2 eff . In parallel with the DEF model, when A has an appropriate form, for instance, A = exp( β 0 X 2 / 2) with sufficiently negative β 0 , µ 2 eff becomes negative. Furthermore, for dense enough objects such as NSs this occurs for order-of-unity values of β 0 . The expectation is that the vector field will grow from arbitrarily small perturbations around X µ = 0 due to this tachyonic behavior, which can be called spontaneous vectorization, in exact analogy with spontaneous scalarization. However, we later see that there are many subtle points in this narrative.  \nEven though we presented spontaneous vectorization in the Einstein frame with a nonminimal coupling to matter, it can be converted to a theory of minimal matter coupling and vector-curvature couplings as in spontaneous scalarization (Ramazano˘glu, 2019a). Even more directly, there are theories of spontaneous vectorization that are purely conceived through curvature coupling, with the first example being the Hellings-Nordtvedt theory (Hellings and Nordtvedt, 1973) studied by Annulli et al. (2019)  \nS = 1 16 πG ∫ d 4 x √ -g [ R -F µν F µν -Ω X 2 R -ηX µ X ν R µν ] + S m [Ψ m ; g µν ] , (38)  \nwhere Ω and η are coupling constants. Scalarization through coupling to the Gauss-Bonnet term also has a vector analog as in the action (Ramazano˘glu, 2019b)  \nS = 1 16 πG ∫ d 4 x √ -g [ R -F µν F µν + f ( X 2 ) G ] + S m [Ψ m ; g µν ] , (39)  \nwhich has a similar structure as Eq. (25). The action (39) leads to BH vectorization as well, unlike the other models we have seen thus far.  \nThe models of vectorization that we have examined are merely specific examples, and any theory of scalarization can be generalized to vector fields in principle. In this regard, Brihaye and Verbin (2020) studied the coupling of a vector field to a scalar, and Oliveira and Pombo (2021) and Brihaye et al. (2022) studied the coupling between two vector fields to obtain vectorized compact objects. Kase et al. (2020) considered generalized Proca theories (Heisenberg, 2014; Tasinato, 2014) with various couplings of a massive vector field and showed that spontaneous vectorization occurred in those models as well. Ramazano˘glu (2018b) investigated models where the effective mass of the vector field was generated by the Higgs mechanism, thus preserving the gauge symmetry. On a separate front, conformal scaling of the metric in the matter action (36) can be generalized to disformal transformations, and spontaneous vectorization still occurs (Minamitsuji, 2021; Ramazano˘glu and Unluturk, 2019). In terms of approximate solutions, Garcia-Saenz et al. (2022) recently calculated the Schwarzschild QNMs of nonminimally coupled vector fields. There have also been efforts to study all these phenomena in a unified manner in compact binaries, using more generic EFT tools (Khalil et al. , 2019).  \nTo summarize, coupling of the vector field to any nonvanishing term in the action can be considered for spontaneous vectorization (Ramazano˘glu, 2019b), and most options have been considered in at least a preliminary sense. The case of all possible couplings, an analog of the minimal action of scalarization in Sec. II.B.1, was studied by Garcia-Saenz et al. (2021). They also revealed some with the fundamental problems of vectorization, as we later discuss more broadly.  \nVectorization models strongly resemble scalarization models in terms of their actions, but this can be misleading. One major difference is in the number of additional degrees of freedom. Scalar fields contribute a single new degree of freedom irrespective of whether they are part of spontaneous scalarization provided that the field equations are of second order in derivatives. For vectorization, this is not the case. The action (37) and all of the vectortensor theory actions that we have discussed break the gauge freedom found in massless vector fields. 14 As a result, the vectors of vectorization models carry 3 degrees of freedom instead of the 2 found in electromagnetism. This is not an immediate reason of concern, since minimally coupled massive vectors, known as Proca fields (Proca,  \n1936), also break the gauge symmetry, and they still provide a classical field theory with well-behaved dynamics. However, as we later see, the extra degree of freedom is far more problematic in vectorization.  \nAlthough the vectorization models are designed to have a linear tachyonic instability in analogy with scalarization, the presence of the third degree of freedom, whose dynamics is not obvious, casts doubt on the success of this design. Namely, whether vectorization can indeed be described as a tachyonic instability of the vector modes that is then quenched nonlinearly is not apparent. This concern is amplified by the fact that vectorized compact objects in the models studied thus far seem to have some striking differences with respect to their scalarized counterparts.  \nAs an example, consider vectorized BHs in the theory described by the action (39). Entropy can be a measure of stability for BHs, where higher entropy solutions are favored in terms of stability. Entropies of spherically symmetric vectorized BHs were numerically computed by Barton et al. (2021), and GR BHs were shown to be entropically favored over vectorized ones in all cases. The verdict of stability is not definite without time evolution, but this is a strong indication for the instability of vectorized BHs. In contrast, scalarized BHs in scalar-Gauss-Bonnet theories, such as Eq. (25) can be entropically favored over those in GR for appropriate coupling functions (Doneva and Yazadjiev, 2018b), and their stability is also indicated through other methods (see Sec. IV.A.2), which makes the stability of vectorized objects even more suspect.  \nScalarized and vectorized NSs also show major differences. For example, the dependence of the strength of vectorization 15 on the NS mass has no discernible pattern, as can be seen in Fig. 24, whereas scalarization occurs in a finite NS mass interval and vanishes at the boundary mass values; see Fig. 1. It is also known that while scalarized NSs are continuously connected to those of GR as the theory parameters are smoothly changed (for example, in the limit β → 0 in Fig. 1), this is not the case in vectorization (Minamitsuji, 2020a; Ramazano˘glu, 2017). The latter is known to be an indication of instability in certain scalarization theories; see Mendes and Ortiz (2016). Thus, qualitative differences between scalarized and vectorized objects are also apparent for NS and, moreover, point to the instability of the vectorized ones.  \nThese observations suggest that the vectorized compact object solutions obtained thus far may not be stable. Hence, they are not astrophysically relevant, at least in  \nFigure 24 The strength of vectorization measured by A ( r = 0) -1 as a function the NS mass M for various values of β 0 and µ X in Eq. (37). Top row: µ X = 1 . 6 × 10 -11 eV. Middle row: µ X = 8 . 0 × 10 -12 eV. Bottom row: µ X = 4 . 8 × 10 -12 eV. Left column: β 0 = -4. Middle column: β 0 = -5. Right column: β 0 = -6. The dashed and dot-dashed lines, respectively, indicate solutions where the energy density ε and Φ = -n µ X µ (with n µ the normal vector field to spatial hypersurfaces) do not monotonically decrease with radius. From Ramazano˘glu (2017).  \n<!-- image -->  \nthe case of spherically symmetric spacetimes. A close look at the instability mechanism of vectorization reveals the underlying reasons for this in Sec. V.B.', 'B. Vectorization and ghosts': "Despite the fact that the vectorization phenomena that we have presented are directly inspired by scalarization, we have seen noteworthy differences between the two mechanisms. Even though it was puzzling for a time, this is now strongly suspected to be related to the fact that vectorization models suffer from ghost and/or gradient instabilities in addition to tachyonic ones, as demonstrated by Garcia-Saenz et al. (2021) and Silva et al. (2022). Here we first follow the latter and later outline the methodology of the former. 16  \nWe begin by reviewing the basic aspects of different types of instabilities (Demirbo˘ga et al. , 2022). Recalling our discussion of the tachyon in Secs. II.A.1 and II.A.2, we consider the linearized scalar-field equation (5) in 1+1 dimensions,  \ng tt ∂ 2 t δφ + g xx ∂ 2 x δφ = µ 2 φ, (40)  \nwhere we assume for simplicity that the metric is diagonal with constant components. The common case without instabilities is when g tt < 0, g xx > 0, and µ 2 > 0. What happens when each of these three terms changes its sign while the other two stay the same?  \nFor a plane wave mode δφ ( t, x ) ∝ exp[i ( ωt -kx )], the dispersion relation is given by  \nω ( k ) = √ ( µ 2 + g xx k 2 ) / ( -g tt ) . (41)  \nWe saw that for µ 2 < 0 the mode behaves as a tachyon, where ω ( k ) becomes imaginary for sufficiently small | k | . This leads to exponential growth as discussed, and the fastest growing mode behaves as ∼ exp[ √ µ 2 /g tt t ]; i.e., the growth rate has an upper limit.  \nA ghost occurs when g tt > 0. There is also exponential growth in this case; however, the growth rate diverges with increasing wave number as ∼ exp[ √ g xx /g tt kt ] and hence has no upper bound. The case g xx < 0 also gives the same asymptotic growth rate as the ghost and is called a gradient instability. 17 In summary, ghost and gradient instabilities as we presented them here have qualitative differences from tachyons, a with main one being the arbitrarily fast growth rates. This means that a welldefined time evolution can be problematic in such theories, and it is more difficult to devise nonlinear quenching mechanisms.  \nThe spontaneous vectorization model of Eq. (36) modifies the 'mass term' in the field equation for the vector field. Hence, the first impression might be that it leads to a tachyonic instability as in scalarization. However, this approach is naive, and a thorough analysis reveals that theories of vectorization necessarily carry ghost instabilities in addition to tachyonic ones. Following Silva et al. (2022), one way to understand this is to use the Stueckelberg trick (Ruegg and Ruiz-Altaba, 2004) where we introduce a scalar field ψ into the action (36) using the substitution  \nX α → X α + µ -1 X ∇ α ψ. (42)  \nThis leads to the scalar equation of motion  \n□ ψ = -¯ g µν [ µ X ∇ µ X ν + 1 2 ( ∇ µ log ˆ z )( ∇ ν ψ + µ X X ν ) ] , (43)  \nwhere the covariant derivative ∇ is compatible with the new metric  \n¯ g µν = ˆ z g µν , ˆ z ≡ µ 2 eff µ 2 X = 1 -8 πG ∗ A 4 Λ ˜ T µ 2 X . (44)  \nThe signature of ¯ g αβ changes with the sign of ˆ z , which controls the change of sign of the effective mass of X µ in Eq. (37). In other words, when we try to instigate a tachyon on X µ by having negative ˆ z in parts of the spacetime, ψ carries ghost and gradient instabilities; cf. Eq. (40). This is a degree of freedom whose dynamics is governed by the effective metric ¯ g µν , not the spacetime metric g µν .  \nIt is curious to see that the Stueckelberg trick reveals the ghost, even though it is not apparent in the vector field equation (37). However, a more careful analysis can clarify this issue (Silva et al. , 2022). Note first that the vector field equation (37) is not manifestly hyperbolic; i.e., its principal part, the term with the highest derivatives, is not the wave operator, since  \n∇ µ F µν = ∇ µ ∇ µ X ν -∇ µ ∇ ν X µ . (45)  \nThe second term on the right-hand side of Eq. (37) vanishes in Proca theory ( A = 1) making the principal part the wave operator. This is not the case in general due to the so-called (generalized) Lorenz condition  \n∇ µ (ˆ zX µ ) = 0 , (46)  \nwhich is obtained by acting on both sides of Eq. (37) with ∇ µ and recalling the antisymmetry of F ρµ . We can use Eq. (46) to obtain a manifestly hyperbolic form of the linearized field equations on a fixed background of vanishing vector fields as  \n∇ µ ∇ µ X ν +( ∇ µ ln ˆ z ) ∇ ν X µ = M νµ X µ , (47)  \nwhere we defined the mass-squared tensor  \nM νµ = ˆ zµ 2 X g νµ + R νµ -∇ ν ∇ µ ln | ˆ z | . (48)  \nIn Eq. (47) and (48) the metric and ˆ z have their fixed background GR values. This is a generalized massive wave equation.  \nSince ˆ z = 1 in the absence of matter and it becomes negative in a continuous manner inside matter for an instability to exist at all, ˆ z necessarily vanishes at some points of a NS spacetime. The ∇ ν ∇ µ ln ˆ z term introduced by the vectorization contains powers of ˆ z -1 , which means that the effective mass-squared tensor diverges at such points. Alternatively, if we move the ˆ z -1 terms to the other side of the equation, the principle part becomes ˆ z □ X ν ; i.e., it changes its sign with ˆ z the same way as the ψ field in Eq. (43). Hence, if we properly rewrite the vector field equation into a manifestly hyperbolic form, the ghost is apparent.  \nEquations (47) and (48) show that calling µ 2 eff = ˆ zµ 2 X the effective mass squared of the theory is misleading, unlike in the case of scalarization. Even though µ 2 eff replaces the mass term in the vector field equation of a minimally coupled Proca field, it is not actually the effective mass of physically propagating degrees of freedom. The physical mass is determined by Eq. (48), which diverges at points where ˆ z = 0.  \nWe have presented the appearance of the ghost in two separate ways, the Stueckelberg mechanism of Eq. (43) and the direct vector field formulation of Eq. (47). The degrees of freedom in the two pictures are related to each other in a nontrivial manner through Eq. (42). Hence the presence of the ghost manifests differently. Nevertheless, the criterion for the appearance of the ghost ˆ z = 0 is the same in both cases. The Stueckelberg picture has an advantage in that it can straightforwardly identify the ghost as a scalar degree of freedom ψ at decoupling, which is well known in the case of the Proca field.  \nOne might at first think that the ghost and gradient instabilities can be regularized by nonlinearities. After all, the tachyon of scalarization is also an instability, but it is eventually quenched by nonlinear effects as explained in Sec. II.A. In such a scenario, one would expect the vectorized objects to be free of instabilities even though the GR solutions are unstable. However, current evidence points away from this. Even though detailed mathematical analyses of the partial differential equations have not been performed to ensure nonlinear instability, Demirbo˘ga et al. (2022) showed that spherically symmetric vectorized NSs in the theory (36), which could be potential stable endpoints of vectorization, are also unstable to ghosts and gradient instabilities at the perturbative level.  \nEven if stable solutions can be found in other vectorization models, there is a more fundamental problem due to the modification to the principal part of the vector equations in equation (43) or (47). The coefficient of the wave operator vanishes at certain points in spacetime or, equivalently, the physical effective mass diverges if we move these coefficients to the other side of the equation. In either case, this leads to divergent time derivatives for arbitrarily small perturbative vector field values, which renders the initial assumptions of a perturbative approach invalid. Overall, it is not clear whether one can even define a well-posed time evolution, let alone one that can somehow suppress the exponential growth of the vector field (Silva et al. , 2022). Further elucidation of these issues requires a more rigorous mathematical investigation, which has not yet been attempted.  \nThus far we have investigated the ghost and gradient instabilities of the specific vectorization theory in the action (36), but similar results are known to exist for all theories of vectorization (Garcia-Saenz et al. , 2021; Silva et al. , 2022) due to similar mechanisms. To summarize, replacing the scalar of scalarization with a vector field results in theories where the central issue is not the nature  \nof the field the instability lives on, but rather the type of instabilities that arise. We cannot construct vectorization models that exhibit only the tachyonic instability commonly associated with the scalarization mechanism, which can be benign and tamed by nonlinear effects. Rather, vectorization exhibits far more threatening instabilities.  \nWe have closely followed Silva et al. (2022), but GarciaSaenz et al. (2021) reached similar conclusions via an alternative approach. They considered the most general action that contained a vector field and the metric (Heisenberg, 2014; Tasinato, 2014) and truncated its action at the quadratic order around GR to obtain the linearized field equations. The main qualitative difference from the scalar minimal action of Sec. II.B.1 is that one can demonstrate that if any of the physical degrees of freedom have a tachyonic instability, then there are also degrees of freedom with ghost or gradient instabilities, which are manifested by a change of sign of coefficients in the dispersion relationship. The conclusion is the same: vectorization based on tachyonic instabilities necessarily has ghost or gradient instabilities as well.  \nRecent work has shown that the issues with vectorization also happen in simpler theories, with an example being minimally coupled self-interacting vector field theories. For instance, for A ( x ) = 1 the action (36) becomes the minimally coupled Proca theory, which is well posed. However, generalizations of the potential beyond a mass term, that is, 2 µ 2 X x → 4 V ( x ) for some generic function V , results in ill-posedness (Clough et al. , 2022; Coates and Ramazano˘glu, 2022; Mou and Zhang, 2022). The underlying reason for this is the fact that there is an analog of the generalized Lorenz condition (46) in this case as well, which ultimately leads to the same result: there is a degree of freedom in the theory whose dynamics is governed by an effective metric, an analog of ¯ g µν in Eq. (44). This new metric depends on the vector field itself and can lose its Lorentzian nature depending on how X µ evolves in time (Coates and Ramazano˘glu, 2022).  \nEsposito-Far'ese et al. (2010) provided some of the earliest examples of the ill-posedness of self-interacting vector field theories that we have mentioned in a cosmological context, and it was recently demonstrated that initially healthy configurations of self-interacting vectors naturally evolve to a point where hyperbolic evolution becomes impossible (Clough et al. , 2022; Coates and Ramazano˘glu, 2022, 2023; Mou and Zhang, 2022). This means that nonlinear extensions of the Proca theory are more delicate than their scalar counterparts in terms of providing physical theories. It seems that the underlying reason for the 'failure' of vectorization is a much more general phenomenon in theories with conditions such as Eq. (46), at least for the models conceived thus far. In this sense, exploration of vectorization has been leading to a deeper understanding of vector fields in general, which has implications beyond theories of gravitation.  \nWe close this section with a pertinent comment on ghost  \ninstabilities. We encountered them as an unappealing artifact in vectorization theories, but one could consider making them the driver of a spontaneous growth mechanism. This option was indeed considered prior to their discovery in vectorization models. The simplest example is another analog of the DEF model where the conformal scaling function A depends on derivative terms (Ramazano˘glu, 2018a)  \nS = 1 16 πG ∗ ∫ d 4 x √ -g [ R -2 g µν ∂ µ φ∂ ν φ -2 µ 2 φ 2 ] + S m [Ψ m ; A 2 ( K ) g µν ] , (49)  \nwith K = -g µν ∂ µ φ∂ ν φ/ 2. The scalar equation of motion is  \n∇ µ (ˆ z ∇ µ φ ) = µ 2 φ, (50)  \nwhere ˆ z ≡ 1 + 4 πG ∗ ˜ T A 4 (d ln A / d K ). Unlike for the tachyonic instability where the mass term is modified, the derivative terms are modified in this case.  \nThe nature of Eq. (50) can be better understood when we consider the linearized equation for perturbative values of the scalar around a GR background and highlight the principal part as  \nˆ z □ δφ + · · · = 0 . (51)  \nFor an appropriate choice such as A = exp( β 0 K ) with large enough β 0 > 0, the sign of ˆ z can become negative in the presence of matter. In all spontaneous scalarization theories considered thus far, the source of the instability was the change of sign of the effective mass term, which resulted in a tachyon. In Eq. (51) the overall sign of the wave operator changes; hence, we have both ghost and gradient instabilities, recalling Eq. (40) and the subsequent discussion. Spontaneous scalarization that might arise from this new mechanism is called ghost-based spontaneous scalarization (Ramazano˘glu, 2018a). It has similar problems as vectorization in that it is hard to tame the ghost and other instabilities; hence, we do not end up with a theory with sensible dynamics.", 'C. Spontaneous tensorization': 'The basic mechanism of spontaneous scalarization can be generalized beyond vector fields, resulting in spontaneous tensorization . However, all such known theories also feature ghosts and suffer from similar problems with vectorization, aside from a single possible example of spinor-tensor theory.  \nThe simplest mathematical generalization of scalarization beyond vectors occurs for p -form fields, i.e., totally antisymmetric tensors X µ 1 ...µ p of rank (0 , p ). This is due to the fact that a vector field X µ is a 1-form field, and  \nthe actions that we encountered in spontaneous vectorization can be readily generalized to this class of fields as (Ramazano˘glu, 2019c)  \nS = 1 16 π ∫ d 4 x √ -g [ R -F µ 1 ...µ p +1 F µ 1 ...µ p +1 ] + S m [ Ψ m ; A 2 ( X 2 ) g µν ] , (52)  \nwhere F µ 1 ...µ p +1 = ( p + 1) ∇ [ µ 1 X µ 2 ...µ p +1 ] , with the square brackets denoting antisymmetrization and X 2 = X µ 1 ...µ p X µ 1 ...µ p . It is then straightforward to show that the analog of µ 2 eff becomes negative inside NSs for an appropriate choice of A ( x ). However, these theories carry ghost instabilities similar to those of vectors when we try to establish spontaneous growth from tachyonic instabilities since they possess an analog of the Lorenz condition (46) (Silva et al. , 2022).  \nSpontaneous tensorization of a symmetric rank (0 , 2) tensor field f µν , which is spin 2, is a natural avenue to investigate after spin-0 scalars and spin-1 vectors. Ghosts appear in this case once more since the metric in our gravity theories g µν is also a spin-2 field. Note that spontaneous tensorization theories include terms where the field that tensorizes, in this case f µν , is coupled to the metric g µν . However, interacting spin-2 fields are known to generically lead to ghost degrees of freedom, rendering most such theories unphysical (Boulware and Deser, 1972; de Rham, 2014) and have led to ghost-free theories being discovered more recently in the form of massive gravity and bigravity (Hassan and Rosen, 2012; de Rham, 2014; de Rham et al. , 2011). Matter can couple to one of the spin-2 fields in the novel ghost-free theories, but coupling to both metrics generically invokes a ghost again (Yamashita et al. , 2014). For example, trying to mimic the DEF model by having a conformal scaling function that depends on f µν , i.e., changing the matter action as g µν →A ( f µν , g µν ) g µν , leads to ghosts. Overall, there is no known form of bigravity theory that features an analog of spontaneous scalarization.  \nSpinor fields present the most interesting case of generalizing scalarization, perhaps aside from vectors. 18 The classical Lagrangian for a Dirac spinor is  \nL ψ = 1 2 [ ¯ ψγ µ ( ∇ µ ψ ) -∇ µ ˆ γ 5 γ µ ψ ] -µ ¯ ψψ, (53)  \nwhere the conventions for gamma matrices and the effect of covariant derivatives on spinors are as given by Ramazano˘glu (2018c). This action provides the usual dispersion relation ω 2 = k i k i + µ 2 for a plane wave of the form ψ ∝ exp[i( ωt -k i x i )]. Since the mass term appears linearly, not quadratically, changing its sign results in the  \nsame dispersion relation. However, a tachyonic spinor action is still possible in the following form:  \nL 5 ψ = 1 2 [ ¯ ψ ˆ γ 5 γ µ ( ∇ µ ψ ) -( ∇ µ ¯ ψ )ˆ γ 5 γ µ ψ ] -µ ¯ ψψ, (54)  \nwith the field equation  \nγ µ ∇ µ ψ -ˆ γ 5 µψ = 0 . (55)  \nIn Eq. (55) ω 2 = k i k i -µ 2 ; that is, the dispersion relation is tachyonic (Chodos et al. , 1985; Jentschura and Wundt, 2012).  \nRamazano˘glu (2018c) used Eq. (54) to obtain the first spontaneous spinorization theory where spinor fields are unstable to growth around GR backgrounds of NSs. Even though the meaning of a ghost as opposed to a tachyon is a subtle issue for spinors, the equation of motion in this theory has divergent coefficients as in the vector field case (47). Hence, this form of spinorization suffers from similar problems to vectorization.  \nMinamitsuji (2020b) proposed an alternative model of spinorization given by the action  \nS = 1 16 πG ∗ ∫ d 4 x √ -g { R + 1 2 [ ¯ ψ ˆ γ 5 γ µ ( ∇ µ ψ ) -( ∇ µ ¯ ψ )ˆ γ 5 γ µ ψ ]} + S m [ Ψ m , A 2 ( ¯ ψψ ) g µν ] . (56)  \nThis model has an unusual feature in that the derivative part of the spinor action is not the canonical one in Eq. (53), but rather the tachyonic one in Eq. (54). The resulting tachyonic equation of motion is  \nγ µ ∇ µ ψ -ˆ γ 5 (4 π A 4 Λ ˜ T ) ψ = 0 , (57)  \nwhere Λ = d ln A / d( ¯ ψψ ). The interpretation is straightforward: we have a tachyonic field with an effective mass µ eff = 4 π A 4 Λ ˜ T ; cf. Eq. (55). The most important aspect of this theory is the fact that the effective mass term and all coefficients of the field equation are regular everywhere; hence, the dynamics does not have any signs of the ill-posedness that we have found in all generalizations of scalarization thus far. In other words, the action (57) is the only known analog of spontaneous scalarization for nonscalar fields that does not suffer from ghost or gradient instabilities.  \nOur results on spontaneous tensorization suggest a nogo theorem for generalizing scalarization beyond scalar fields, at least if we want to avoid ghosts. It is curious that the only exception to this trend is the relatively exotic case of spinorization, which invites studies of the deeper reasons that make spontaneous scalarization difficult to replicate for other types of fields.', 'VI. OPEN PROBLEMS AND FUTURE PERSPECTIVES': 'We chose to review the literature by starting from the theoretical underpinnings of the scalarization mechanism,  \nmoving on to discussing NSs and BHs separately, and then returning to model building in order to discuss generalization of the mechanism to other fields. Most but not all open problems have already been mentioned in one or more of the previous sections and also discussed to some extent. Nevertheless, we revisit them in this section and cover any additional areas that require further development, opting for a concise summary of future perspectives.', 'A. Scalarization and cosmology': "As discussed in Sec. II.B.4, one of the key challenges for scalarization is understanding whether it is compatible with cosmology. Recall that the main premise of scalarization is that there is a constant value of the scalar field φ 0 for which spacetimes of GR become admissible solutions to the field equations, and it is these solutions that describe stationary objects that we expect (from an observational perspective) to carry no scalar charge. For this to be true, cosmic evolution needs to comply and drive φ to φ 0 in the late Universe; otherwise, the entire Universe will in effect be 'scalarized.' It was pointed out early on by Damour and Nordtvedt (1993) and more recently by Anderson et al. (2016) and Franchini and Sotiriou (2020) that reaching this preferred value for the scalar fields in the late Universe is not generic for simple models of scalarization and instead requires severe fine-tuning of initial conditions.  \nThis could well be an artifact of not having the complete theory, and ideally one would expect the need for finetuned initial conditions to disappear by addressing further corrections to the model. A first step in this direction was made by Antoniou et al. (2021a) [who were inspired by Damour and Nordtvedt (1993)], who showed that the mixed model of scalarization discussed in Sec. II.B.4 has GR with a constant scalar as a late-time cosmic attractor for the right sign of the coupling between the scalar and the Ricci scalar. This demonstrated that in principle scalarization models can be made compatible with late-time cosmology and provided a recipe for doing so: include in the action terms that will dominate in late cosmology (for example, Ricci coupling) over the terms that control the onset of scalarization (for example, Gauss-Bonnet invariant) and hence impose the desired cosmological behavior. However, the specific model discussed as an example by Antoniou et al. (2021a) is by no means unique, and there have been many other attempts to address similar concerns (Chen et al. , 2015; Erices et al. , 2022; Minamitsuji and Tsujikawa, 2023). Perhaps more importantly, understanding the efficiency of the attractor behavior in reproducing the behavior of GR (with a cosmological constant) quantitatively, down to evolution of perturbations, structure formations, etc., certainly merits further investigation and could lead to  \nconstraints on realistic models of scalarization.  \nA second point of friction between scalarization and cosmology relates to the early Universe. Broadly speaking, scalarization is controlled by the coupling between the scalar-field and curvature invariants, so it is inevitable that these couplings will become increasingly important as one runs cosmic evolution backward and move to higher and higher curvatures. In particular, they will tend to dominate when the size of the Universe becomes comparable to the size of the compact objects that we want scalarization to affect today (Antoniou et al. , 2021a). It was pointed out by Anson et al. (2019a,b) that quantum fluctuations could then seed a scalarizationlike instability during inflation, and that it would be hard to prevent this linear instability by adding corrections to the action. It should be stressed, however, that one does not need the scalar field to necessarily remain constant or have its current value through the evolution of the Universe. Hence, one might not need to prevent such an instability, but might instead just quench it nonlinearly, exactly as it happens in scalarization itself. The scalar field could then smoothly evolve away from the vacuum that leads to this instability before it becomes relevant as one moves backward in cosmic time. The key question here is whether there is an EFT applicable to the early Universe that contains the scalarization models as late Universe limits but that can also host an inflationary scenario compatible with observations.  \nAs a final remark in relation to cosmology, we emphasize that we have discussed here only the cosmological implications of known scalarization models and only under the assumption that the scalar field is cosmologically subdominant at late times. There are a plethora of generalized scalar-tensor theories that have been used in the context of inflation or dark energy and that exhibit couplings between the scalar and curvature invariants similar to those employed for scalarization models. In most cases, however, the nature of the couplings differs significantly from that of scalarization models. Reviewing such attempts goes well beyond the scope of this work. Exploring whether a scalar field that exhibits scalarization around compact objects could also play a role in inflation or account for dark energy would certainly be interesting, but putting together such a model (with a single scalar) seems particularly challenging. Note that in this context evading constraints on the coupling between a scalar and the Gauss-Bonnet invariant coming from the speed of GWs requires the scalar to be cosmologically subdominant (Franchini and Sotiriou, 2020).", 'B. Dynamical evolution': "The study of the dynamical nonlinear regime of scalarization requires nonperturbative time-evolution schemes, whose availability varies among models. Finding formula-  \ntions of gravity theories that are amenable to numerical time evolution, i.e., numerical relativity, is nontrivial even in GR (Baker et al. , 2006; Campanelli et al. , 2006; Pretorius, 2005), and the issue can become even more complicated for alternative gravity theories.  \nTime-domain, nonlinear evolutions of DEF-like models have been made since the early work by Novak (1998a,b); and Novak and Ibanez (2000), and binary inspirals have been performed for nearly a decade (Barausse et al. , 2013). This was possible thanks to the relative simplicity of the field equations, which allows one to establish their wellposedness (Salgado, 2006; Salgado et al. , 2008). However, detailed and fully nonlinear numerical results in relation to GWs, a primary tool to test scalarization, are still lacking, aside from the study by Shibata et al. (2014).  \nThe picture is much different for BH scalarization, where the coupling to the Gauss-Bonnet term can result in a complicated causal structure of the spacetime in dynamical situations. It was demonstrated at both linear (Bl'azquez-Salcedo et al. , 2020a,b, 2018) and nonlinear (East and Ripley, 2021a,b; Ripley and Pretorius, 2020) levels that part, but not all, of the parameter space of spontaneous scalarization where scalarized BHs can be found loses hyperbolicity when the theory is evolved in time in a certain set of gauges, which means that no predictions can be obtained. The well-posedness of the initial-value problem in broader classes of modified gravity theories, including scalar-Gauss-Bonnet gravity, was studied by Kov'acs (2019); Kov'acs and Reall (2020a,b); and Papallo (2017). At the moment there is no consensus on whether this loss of hyperbolicity can be cured with a better gauge choice or if it is intrinsic to the evolution equations, which calls for further work on the issue. Furthermore, there are no detailed studies of binary evolution and merger even for the part of the parameter space for which hyperbolic time evolution has been shown to exist. The dynamics of more general models where various coupling terms are present as in Eq. (14) are not available either. The lack of results on all of these fronts presents important future research directions.  \nIt can be the case that an EFT, obtained as a certain limit or truncation of a more fundamental theory, is ill posed while the latter is or is expected to be well posed. A typical example is relativistic hydrodynamics once viscosity has been taken into account. Indeed, it has been suggested to introduce techniques used in hydrodynamics to gravity in order to render the problematic time evolution of some theories hyperbolic by suitably modifying the equations while keeping the end point of evolution the same (Cayuso et al. , 2017). This approach has seen increasing use in alternatives to GR, including the specific case of scalar-Gauss-Bonnet theories under certain symmetries (Franchini et al. , 2022), and it is another avenue to explore the dynamics of scalarization in cases where current methods are inadequate. The applicability and power of such methods to general cases, for instance, fully  \nnonlinear 3+1 evolution of a highly dynamical system, are yet to be confirmed.", 'C. Model building in and beyond scalarization': "We have covered several models of spontaneous scalarization and their phenomenology for NSs and BHs. What they have in common is that scalarization is triggered by a tachyonic instability, whose threshold is controlled by fewer than a handful of terms in the action (Andreou et al. , 2019). Yet, the properties of the final scalarized configuration depend on the nonlinearities and coupling terms that contribute beyond the linear level. The specific choices of these nonlinear interactions have been shown to crucially affect the stability (Antoniou et al. , 2022; Bl'azquez-Salcedo et al. , 2018; Silva et al. , 2019) and the scalar charge (Antoniou et al. , 2021b; Doneva et al. , 2018a; Macedo et al. , 2019; Silva et al. , 2019) of scalarized solutions, as well as the cosmological behavior of the corresponding theory (Anson et al. , 2019a; Antoniou et al. , 2021a; Erices et al. , 2022). Thus, such choices affect the viability and observability of the models. As a result, although certain models of scalarization, such as the original DEF model, have been studied fairly exhaustively, the exploration of the broader class of theories that exhibit scalarization, and their phenomenology, has only started.  \nFurthermore, one could consider introducing new classes of scalarization based on instabilities that are not tachyonic in nature. As discussed in Sec. V, success in this direction has thus far been limited, but it is nevertheless an interesting prospect. Alternatively, one can abandon the idea of a linear instability quenched by nonlinearities altogether, as discussed in Sec. II.C.5. Doneva and Yazadjiev (2022) demostrated that a theory featuring a coupling to the Gauss-Bonnet invariant with a leading-order expansion with respect to the scalar field that is proportional to φ 4 can still have scalarized BHs below a certain mass threshold, although the Kerr solution is linearly stable with respect to massless scalar perturbations (Bl'azquez-Salcedo et al. , 2022). This observation remains true in more general theories, such as Gauss-Bonnet gravity with multiple scalar fields (Staykov and Doneva, 2022). This strongly suggests the existence of a class of theories in which scalarization is a purely nonlinear phenomenon.  \nScalarization models that have been considered thus far typically assume that the scalar field does not couple to matter directly in a suitable choice of variables, usually called the Jordan frame, which suffices for the theory to satisfy the WEP. However, as pointed out by Coates et al. (2017), this assumption might not be necessary because scalarization itself forces the scalar field into a trivial configuration away from a specific compact object. Therefore, scalarization models can be indistinguishable  \nfrom GR in regard to the WEP for all current observations. It was further argued by Coates et al. (2017) and Franchini et al. (2018) that particular couplings to matter that do not disrupt scalarization could introduce a Higgs-like mechanism in gravity: the scalar field changes value only near compact objects, and its coupling to matter changes the properties of the standard model in these regions. This is a largely unexplored possibility.  \nAs discussed in Sec. V, another important open question is whether one can generalize the scalarization mechanism to other fields. Doing so in the context of a theory that is free of pathologies has thus far proven to be challenging, except perhaps in the case of spinors (Minamitsuji, 2020b). This is in part because controlling the dynamics of the extra degree of freedom is notoriously difficult in general in gravity theories with additional vector or tensor fields (de Rham, 2014), and vectorization or tensorization models are no exception.  \nWe stress that the work on scalarization and its extensions has focused on isolated compact objects. In fact, theoretical considerations for putting together scalarization models have also been heavily influence by studies of isolated BHs and NSs. However, as previously discussed, scalarization could happen dynamically in a binary (Barausse et al. , 2013; Doneva et al. , 2022; Elley et al. , 2022; Palenzuela et al. , 2014; Shibata et al. , 2014; Silva et al. , 2021b; Taniguchi et al. , 2015). Apart from finding ways to model this effect (Khalil et al. , 2022) and potentially constraining it with GW observations, it would be particularly interesting to further explore which properties of a binary control dynamical scalarization, and how this differs among the various scalarization models.", 'D. Observational prospects': 'Perhaps the most attractive feature of scalarization is that the tachyonic instability typically leads to large scalar-field amplitudes before it is quenched, which results in nonperturbative deviations from GR in high curvature regions. At the same time, the theory mimics GR closely for weak gravity, easily satisfying existing tests. This suggests that new fundamental scalar fields that exhibit this behavior could be discovered with strong-gravity observations.  \nCurrent observations put stringent bounds on scalar radiation from pulsars that have already almost completely ruled out the original (massless) DEF model (Zhao et al. , 2022). Electromagnetic observations of compact objects in binaries or surrounded by some form of matter also have the potential to extend such results to more general models of scalarization. For example, x-ray data from NS surfaces have recently been considered as a possible way to detect scalarized objects (Silva and Yunes, 2019a,b). The masses, radii, and moment of inertia of NSs can be inferred through such observations (Bogdanov et al. ,  \n2019a,b, 2021). Future observation could be used to distinguish scalarized NSs from nonscalarized ones by, for instance, using EOS-independent relations. The astrophysical signatures of most extensions of the DEF model, especially those that allow BH scalarization, have not yet been studied in detail, because they are recent.  \nThe most interesting current development in strong gravity is the advent of GW astronomy, which is an arena for exploring spontaneous scalarization. The effects of scalarization on GWs of coalescing NSs and/or BHs have been investigated to some extent, mostly in regard to the DEF model (Khalil et al. , 2022, 2019; Niu et al. , 2021; Sennett and Buonanno, 2016) or scalar-GaussBonnet gravity (Wong et al. , 2022). The most pressing issue on this front is the development of accurate waveform models that ideally would cover the inspiral, merger, and ringdown stages of coalescing binaries. Developing such waveform models in GR has been (and remains) a formidable task that combines tools from post-Newtonian theory, numerical relativity, black-hole perturbation theory, and gravitational self-force. Progress on all these fronts is necessary for one to develop accurate scalarized waveform models. Stars that undergo core collapse are another interesting source of GWs. Albeit expected to be rare, such events may produce large bursts of scalar radiation (Kuan et al. , 2022; Rosca-Mead et al. , 2019; Sperhake et al. , 2017) and may contribute to the stochastic GW background, which could in principle be probed by pulsar-timing arrays. In addition, some models of scalarization (but not exclusively) predict the existence of extra polarizations in GWs, which provides another way to test them.  \nLast, future spaceborne GW detectors such as LISA will bring another channel for testing modified gravity. Especially interesting with respect to scalarization is the case of extreme mass-ratio inspirals when the secondary object is scalarized (Maselli et al. , 2020, 2022). Data from future GW detectors such as LISA or third generation ground-based detectors such as the Einstein Telescope and the Cosmic Explorer, combined with more precise electromagnetic observations such as those from improved x-ray observatories, will be crucial for detecting or ruling out scalar fields that exhibit scalarization.', 'ACKNOWLEDGMENTS': 'D.D.D. acknowledges financial support via the Emmy Noether Research Group, which was funded by the German Research Foundation (DFG) under Grant No. DO 1771/1-1. S.S.Y. thanks the University of Tubingen for the financial support and acknowledges partial support from Bulgarian NSF Grant No. KP-06-H28/7. T.P.S. acknowledges partial support from STFC Consolidated Grants No. ST/T000732/1 and No. ST/V005596/1. F.M.R acknowledges support from 2020 Bilim Akademisi GEB ˙ IP  \naward and T UB ˙ ITAK Project No. 122F097. H.O.S. acknowledges funding from Deutsche Forschungsgemeinschaft (DFG) Project No. 386119226. We also acknowledge the networking support of COST Actions No. CA16104 and No. CA15117.', 'REFERENCES': "- Abbott, B. P., et al. (LIGO Scientific, Virgo) (2019), Phys. Rev. X 9 (3), 031040, arXiv:1811.12907 [astro-ph.HE].\n- Abbott, R., et al. (LIGO Scientific, Virgo) (2021a), Phys. Rev. X 11 , 021053, arXiv:2010.14527 [gr-qc].\n- Abbott, R., et al. (LIGO Scientific, VIRGO, KAGRA) (2021b), 'GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run,' arXiv:2111.03606 [gr-qc].\n- Abuter, R., et al. (GRAVITY) (2018), Astron. Astrophys. 615 , L15, arXiv:1807.09409 [astro-ph.GA].\n- Abuter, R., et al. (GRAVITY) (2020), Astron. Astrophys. 636 , L5, arXiv:2004.07187 [astro-ph.GA].\n- Akiyama, K., et al. (Event Horizon Telescope) (2019), Astrophys. J. Lett. 875 , L1, arXiv:1906.11238 [astro-ph.GA].\n- Akmal, A., V. R. Pandharipande, and D. G. Ravenhall (1998), Phys. Rev. C 58 , 1804, arXiv:nucl-th/9804027.\n- Alcubierre, M., J. C. Degollado, D. Nunez, M. Ruiz, and M. Salgado (2010), Phys. Rev. D 81 , 124018, arXiv:1003.4767 [gr-qc].\n- Alexander, S., and N. Yunes (2009), Phys. Rept. 480 , 1, arXiv:0907.2562 [hep-th].\n- Alsing, J., E. Berti, C. M. Will, and H. Zaglauer (2012), Phys. Rev. D 85 , 064041, arXiv:1112.4903 [gr-qc].\n- Altaha Motahar, Z., J. L. Bl'azquez-Salcedo, D. D. Doneva, J. Kunz, and S. S. Yazadjiev (2019), Phys. Rev. D 99 (10), 104006, arXiv:1902.01277 [gr-qc].\n- Altaha Motahar, Z., J. L. Bl'azquez-Salcedo, B. Kleihaus, and J. Kunz (2017), Phys. Rev. D 96 (6), 064046, arXiv:1707.05280 [gr-qc].\n- Altaha Motahar, Z., J. L. Bl'azquez-Salcedo, B. Kleihaus, and J. Kunz (2018), Phys. Rev. D 98 (4), 044032, arXiv:1807.02598 [gr-qc].\n- Anderson, D., and N. Yunes (2019), Class. Quant. Grav. 36 (16), 165003, arXiv:1901.00937 [gr-qc].\n- Anderson, D., N. Yunes, and E. Barausse (2016), Phys. Rev. D 94 (10), 104064, arXiv:1607.08888 [gr-qc].\n- Andersson, N., and K. D. Kokkotas (1998), Mon. Not. Roy. Astron. Soc. 299 , 1059, arXiv:gr-qc/9711088.\n- Andreou, N., N. Franchini, G. Ventagli, and T. P. Sotiriou (2019), Phys. Rev. D 99 (12), 124022, [Erratum: Phys.Rev.D 101, 109903 (2020)], arXiv:1904.06365 [gr-qc].\n- Annulli, L., V. Cardoso, and L. Gualtieri (2019), Phys. Rev. D 99 (4), 044038, arXiv:1901.02461 [gr-qc].\n- Annulli, L., C. A. R. Herdeiro, and E. Radu (2022), Phys. Lett. B 832 , 137227, arXiv:2203.13267 [gr-qc].\n- Anson, T., E. Babichev, C. Charmousis, and S. Ramazanov (2019a), JCAP 06 , 023, arXiv:1903.02399 [gr-qc].\n- Anson, T., E. Babichev, and S. Ramazanov (2019b), Phys. Rev. D 100 (10), 104051, arXiv:1905.10393 [gr-qc].\n- Antoniadis, J., et al. (2013), Science 340 , 6131, arXiv:1304.6875 [astro-ph.HE].\n- Antoniou, G., A. Bakopoulos, and P. Kanti (2018), Phys. Rev. Lett. 120 (13), 131102, arXiv:1711.03390 [hep-th].  \n- Antoniou, G., L. Bordin, and T. P. Sotiriou (2021a), Phys. Rev. D 103 (2), 024012, arXiv:2004.14985 [gr-qc].\n- Antoniou, G., A. Leh'ebel, G. Ventagli, and T. P. Sotiriou (2021b), Phys. Rev. D 104 (4), 044002, arXiv:2105.04479 [gr-qc].\n- Antoniou, G., C. F. B. Macedo, R. McManus, and T. P. Sotiriou (2022), Phys. Rev. D 106 (2), 024029, arXiv:2204.01684 [gr-qc].\n- Arun, K. G., et al. (LISA) (2022), 'New Horizons for Fundamental Physics with LISA,' arXiv:2205.01597 [gr-qc].\n- Arzoumanian, Z., et al. (2014), in Space Telescopes and Instrumentation 2014: Ultraviolet to Gamma Ray , Proc. SPIE, Vol. 9144, p. 914420.\n- Astefanesei, D., C. Herdeiro, A. Pombo, and E. Radu (2019), JHEP 10 , 078, arXiv:1905.08304 [hep-th].\n- Azri, H., and S. Nasri (2021), Phys. Rev. D 103 (2), 024035, arXiv:2012.04694 [gr-qc].\n- Baker, J. G., J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter (2006), Phys. Rev. Lett. 96 , 111102, arXiv:grqc/0511103.\n- Barack, L., et al. (2019), Class. Quant. Grav. 36 (14), 143001, arXiv:1806.05195 [gr-qc].\n- Barausse, E., C. Palenzuela, M. Ponce, and L. Lehner (2013), Phys. Rev. D 87 , 081506, arXiv:1212.5053 [gr-qc].\n- Barausse, E., et al. (2020), Gen. Rel. Grav. 52 (8), 81, arXiv:2001.09793 [gr-qc].\n- Barton, S., B. Hartmann, B. Kleihaus, and J. Kunz (2021), Phys. Lett. B 817 , 136336, arXiv:2103.01651 [gr-qc].\n- Bauswein, A., T. W. Baumgarte, and H. T. Janka (2013), Phys.\n- Rev. Lett. 111 (13), 131101, arXiv:1307.5191 [astro-ph.SR]. Bauswein, A., and H. T. Janka (2012), Phys. Rev. Lett. 108 , 011101, arXiv:1106.1616 [astro-ph.SR].\n- Bauswein, A., H. T. Janka, K. Hebeler, and A. Schwenk (2012), Phys. Rev. D 86 , 063001, arXiv:1204.1888 [astro-ph.SR].\n- Bauswein, A., and N. Stergioulas (2017), Mon. Not. Roy. Astron. Soc. 471 (4), 4956, arXiv:1702.02567 [astro-ph.HE]. Bauswein, A., N. Stergioulas, and H. T. Janka (2014), Phys. Rev. D 90 (2), 023002, arXiv:1403.5301 [astro-ph.SR].\n- Baym, G., T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka (2018), Rept. Prog. Phys. 81 (5), 056902, arXiv:1707.04966 [astro-ph.HE].\n- Bekenstein, J. D. (1993), Phys. Rev. D 48 , 3641, arXiv:grqc/9211017.\n- Beltran Jimenez, J., A. L. Delvas Froes, and D. F. Mota (2013), Phys. Lett. B 725 , 212, arXiv:1212.1923 [astro-ph.CO].\n- Berti, E., L. G. Collodel, B. Kleihaus, and J. Kunz (2021),\n- Phys. Rev. Lett. 126 (1), 011104, arXiv:2009.03905 [gr-qc]. Berti, E., et al. (2015), Class. Quant. Grav. 32 , 243001, arXiv:1501.07274 [gr-qc].\n- Bettoni, D., and S. Liberati (2013), Phys. Rev. D 88 , 084020, arXiv:1306.6724 [gr-qc].\n- Bl'azquez-Salcedo, J. L., D. D. Doneva, S. Kahlen, J. Kunz, P. Nedkova, and S. S. Yazadjiev (2020a), Phys. Rev. D 101 (10), 104006, arXiv:2003.02862 [gr-qc].\n- Bl'azquez-Salcedo, J. L., D. D. Doneva, S. Kahlen, J. Kunz, P. Nedkova, and S. S. Yazadjiev (2020b), Phys. Rev. D 102 (2), 024086, arXiv:2006.06006 [gr-qc].\n- Bl'azquez-Salcedo, J. L., D. D. Doneva, J. Kunz, and S. S. Yazadjiev (2018), Phys. Rev. D 98 (8), 084011, arXiv:1805.05755 [gr-qc].\n- Bl'azquez-Salcedo, J. L., D. D. Doneva, J. Kunz, and S. S. Yazadjiev (2022), Phys. Rev. D 105 (12), 124005, arXiv:2203.00709 [gr-qc].\n- Bl'azquez-Salcedo, J. L., C. A. R. Herdeiro, S. Kahlen, J. Kunz,\n- A. M. Pombo, and E. Radu (2021), Eur. Phys. J. C 81 (2), 155, arXiv:2008.11744 [gr-qc].\n- Bl'azquez-Salcedo, J. L., C. A. R. Herdeiro, J. Kunz, A. M. Pombo, and E. Radu (2020c), Phys. Lett. B 806 , 135493, arXiv:2002.00963 [gr-qc].\n- Bl'azquez-Salcedo, J. L., F. Scen Khoo, and J. Kunz (2020d), EPL 130 (5), 50002, arXiv:2001.09117 [gr-qc].\n- Bogdanov, S., C. O. Heinke, F. Ozel, and T. Guver (2016), Astrophys. J. 831 (2), 184, arXiv:1603.01630 [astro-ph.HE]. (1), L25,\n- Bogdanov, S., et al. (2019a), Astrophys. J. Lett. 887 arXiv:1912.05706 [astro-ph.HE].\n- Bogdanov, S., et al. (2019b), Astrophys. J. Lett. 887 (1), L26, arXiv:1912.05707 [astro-ph.HE].\n- Bogdanov, S., et al. (2021), Astrophys. J. Lett. 914 (1), L15, arXiv:2104.06928 [astro-ph.HE].\n- Boulware, D. G., and S. Deser (1972), Phys. Rev. D 6 , 3368. Breu, C., and L. Rezzolla (2016), Mon. Not. Roy. Astron. Soc. 459 (1), 646, arXiv:1601.06083 [gr-qc].\n- Brihaye, Y., and B. Hartmann (2019), Phys. Lett. B 792 , 244, arXiv:1902.05760 [gr-qc].\n- Brihaye, Y., B. Hartmann, B. Kleihaus, and J. Kunz (2022), Phys. Rev. D 105 (4), 044050, arXiv:2109.12345 [gr-qc].\n- Brihaye, Y., C. Herdeiro, and E. Radu (2019), Phys. Lett. B 788 , 295, arXiv:1810.09560 [gr-qc].  \nBrihaye, Y., and Y. Verbin (2020), Phys. Rev. D  \n102  \n, 124021,  \narXiv:2004.01681 [gr-qc].  \n- Brito, R., V. Cardoso, and P. Pani (2015), Lect. Notes Phys. 906 , pp.1, arXiv:1501.06570 [gr-qc].\n- Bucciantini, N., and J. Soldateschi (2020), Mon. Not. Roy. Astron. Soc. 495 (1), L56, arXiv:2004.00322 [astro-ph.HE]. Campanelli, M., C. O. Lousto, P. Marronetti, and Y. Zlochower (2006), Phys. Rev. Lett. 96 , 111101, arXiv:gr-qc/0511048.\n- Cardoso, V., I. P. Carucci, P. Pani, and T. P. Sotiriou (2013a), Phys. Rev. Lett. 111 , 111101, arXiv:1308.6587 [gr-qc].\n- Cardoso, V., I. P. Carucci, P. Pani, and T. P. Sotiriou (2013b), Phys. Rev. D 88 , 044056, arXiv:1305.6936 [gr-qc].\n- Cardoso, V., A. Foschi, and M. Zilhao (2020), Phys. Rev. Lett. 124 (22), 221104, arXiv:2005.12284 [gr-qc].\n- Cayuso, J., N. Ortiz, and L. Lehner (2017), Phys. Rev. D 96 (8), 084043, arXiv:1706.07421 [gr-qc].\n- Chandrasekhar, S. (1970), Phys. Rev. Lett. 24 , 611.\n- Chen, P., T. Suyama, and J. Yokoyama (2015), Phys. Rev. D 92 , 124016, arXiv:1508.01384 [gr-qc].\n- Cheong, P. C.-K., and T. G. F. Li (2019), Phys. Rev. D 100 (2), 024027, arXiv:1812.04835 [gr-qc].\n- Cherubini, C., D. Bini, S. Capozziello, and R. Ruffini (2002), Int. J. Mod. Phys. D 11 , 827, arXiv:gr-qc/0302095.\n- Chiba, T. (2022), PTEP 2022 , 013E01, arXiv:2104.11362 [gr-qc].\n- Chodos, A., A. I. Hauser, and V. A. Kostelecky (1985), Phys. Lett. 150B , 431.\n- Choptuik, M. W. (1993), Phys. Rev. Lett. 70 , 9.\n- Clark, J., A. Bauswein, L. Cadonati, H. T. Janka, C. Pankow, and N. Stergioulas (2014), Phys. Rev. D 90 (6), 062004, arXiv:1406.5444 [astro-ph.HE].\n- Clifton, T., P. G. Ferreira, A. Padilla, and C. Skordis (2012), Phys. Rept. 513 , 1, arXiv:1106.2476 [astro-ph.CO].\n- Clough, K., T. Helfer, H. Witek, and E. Berti (2022), Phys. Rev. Lett. 129 (15), 151102, arXiv:2204.10868 [gr-qc].\n- Coates, A., M. W. Horbartsch, and T. P. Sotiriou (2017), Phys. Rev. D 95 (8), 084003, arXiv:1606.03981 [gr-qc].\n- Coates, A., and F. M. Ramazano˘glu (2022), Phys. Rev. Lett. 129 (15), 151103, arXiv:2205.07784 [gr-qc].  \n- Coates, A., and F. M. Ramazano˘glu (2023), Phys. Rev. Lett. 130 (2), 021401, arXiv:2211.08027 [gr-qc].\n- Collodel, L. G., B. Kleihaus, J. Kunz, and E. Berti (2020), Class. Quant. Grav. 37 (7), 075018, arXiv:1912.05382 [grqc].\n- Cowling, T. G. (1941), Mon. Not. Roy. Astron. Soc. 101 , 367. Cunha, P. V. P., C. A. R. Herdeiro, and E. Radu (2019), Phys. Rev. Lett. 123 (1), 011101, arXiv:1904.09997 [gr-qc].\n- Damour, T. (2015), Class. Quant. Grav. 32 (12), 124009, arXiv:1411.3930 [gr-qc].\n- Damour, T., and G. Esposito-Far'ese (1992), Class. Quant. Grav. 9 , 2093.\n- Damour, T., and G. Esposito-Far'ese (1993), Phys. Rev. Lett. 70 , 2220.\n- Damour, T., and G. Esposito-Far'ese (1996a), Phys. Rev. D 54 , 1474, arXiv:gr-qc/9602056.\n- Damour, T., and G. Esposito-Far'ese (1996b), Phys. Rev. D 53 , 5541, arXiv:gr-qc/9506063.\n- Damour, T., and G. Esposito-Far'ese (1998), Phys. Rev. D 58 , 042001, arXiv:gr-qc/9803031.\n- Damour, T., and K. Nordtvedt (1993), Phys. Rev. Lett. 70 , 2217.\n- Damour, T., and J. H. Taylor (1992), Phys. Rev. D 45 , 1840. Danchev, V. I., and D. D. Doneva (2021), Phys. Rev. D 103 (2), 024049, arXiv:2010.07392 [gr-qc].\n- Danchev, V. I., D. D. Doneva, and S. S. Yazadjiev (2022), Phys. Rev. D 106 (12), 124001, arXiv:2112.03869 [gr-qc]. DeDeo, S., and D. Psaltis (2003), Phys. Rev. Lett. 90 , 141101, arXiv:astro-ph/0302095.\n- DeDeo, S., and D. Psaltis (2004), 'Testing strong-field gravity with quasiperiodic oscillations,' arXiv:astro-ph/0405067.\n- Deffayet, C., S. Deser, and G. Esposito-Far'ese (2009), Phys. Rev. D 80 , 064015, arXiv:0906.1967 [gr-qc].\n- Degollado, J. C., M. Salgado, and M. Alcubierre (2020), Phys. Lett. B 808 , 135666, arXiv:2008.10683 [gr-qc].\n- Delsate, T., D. Hilditch, and H. Witek (2015), Phys. Rev. D 91 (2), 024027, arXiv:1407.6727 [gr-qc].\n- Demirbo˘ga, E. S., A. Coates, and F. M. Ramazano˘glu (2022),\n- Phys. Rev. D 105 (2), 024057, arXiv:2112.04269 [gr-qc]. Diaz Alonso, J., and J. M. Ibanez Cabanell (1985), Astrophys.\n- J. 291 , 308.\n- Dima, A., E. Barausse, N. Franchini, and T. P. Sotiriou (2020),\n- Phys. Rev. Lett. 125 (23), 231101, arXiv:2006.03095 [gr-qc]. Do, T., et al. (2019), Science 365 (6454), 664, arXiv:1907.10731 [astro-ph.GA].\n- Dolan, S. R. (2013), Phys. Rev. D 87 (12), 124026, arXiv:1212.1477 [gr-qc].\n- Doneva, D. D., L. G. Collodel, C. J. Kruger, and S. S. Yazadjiev (2020a), Phys. Rev. D 102 (10), 104027, arXiv:2008.07391 [gr-qc].\n- Doneva, D. D., L. G. Collodel, C. J. Kruger, and S. S. Yazadjiev (2020b), Eur. Phys. J. C 80 (12), 1205, arXiv:2009.03774 [gr-qc].\n- Doneva, D. D., S. Kiorpelidi, P. G. Nedkova, E. Papantonopoulos, and S. S. Yazadjiev (2018a), Phys. Rev. D 98 (10), 104056, arXiv:1809.00844 [gr-qc].\n- Doneva, D. D., and G. Pappas (2018), Astrophys. Space Sci. Libr. 457 , 737, arXiv:1709.08046 [gr-qc].\n- Doneva, D. D., K. V. Staykov, and S. S. Yazadjiev (2019), Phys. Rev. D 99 (10), 104045, arXiv:1903.08119 [gr-qc].\n- Doneva, D. D., K. V. Staykov, S. S. Yazadjiev, and R. Z. Zheleva (2020c), Phys. Rev. D 102 (6), 064042, arXiv:2006.11515 [gr-qc].\n- Doneva, D. D., A. Va˜n'o Vi˜nuales, and S. S. Yazadjiev (2022),  \nPhys. Rev. D 106 (6), L061502, arXiv:2204.05333 [gr-qc]. Doneva, D. D., and S. S. Yazadjiev (2016), JCAP 11 , 019, arXiv:1607.03299 [gr-qc].  \nDoneva, D. D., and S. S. Yazadjiev (2018a), JCAP 04 , 011, arXiv:1712.03715 [gr-qc].  \nDoneva, D. D., and S. S. Yazadjiev (2018b), Phys. Rev. Lett. 120 (13), 131103, arXiv:1711.01187 [gr-qc].  \n- Doneva, D. D., and S. S. Yazadjiev (2020a), Phys. Rev. D 101 (10), 104010, arXiv:2004.03956 [gr-qc].\n- Doneva, D. D., and S. S. Yazadjiev (2020b), Phys. Rev. D 101 (6), 064072, arXiv:1911.06908 [gr-qc].\n- Doneva, D. D., and S. S. Yazadjiev (2021a), Phys. Rev. D 103 (6), 064024, arXiv:2101.03514 [gr-qc].\n- Doneva, D. D., and S. S. Yazadjiev (2021b), Phys. Rev. D 103 (8), 083007, arXiv:2102.03940 [gr-qc].\n- Doneva, D. D., and S. S. Yazadjiev (2022), Phys. Rev. D 105 (4), L041502, arXiv:2107.01738 [gr-qc].\n- Doneva, D. D., S. S. Yazadjiev, K. D. Kokkotas, and I. Z. Stefanov (2010), Phys. Rev. D 82 , 064030, arXiv:1007.1767 [gr-qc].  \nDoneva, D. D., S. S. Yazadjiev, K. V. Staykov, and K. D. Kokkotas (2014a), Phys. Rev. D 90 (10), 104021, arXiv:1408.1641 [gr-qc].  \nDoneva, D. D., S. S. Yazadjiev, N. Stergioulas, and K. D. Kokkotas (2013), Phys. Rev. D 88 (8), 084060, arXiv:1309.0605 [gr-qc].  \nDoneva, D. D., S. S. Yazadjiev, N. Stergioulas, and K. D. Kokkotas (2018b), Phys. Rev. D 98 (10), 104039, arXiv:1807.05449 [gr-qc].  \nDoneva, D. D., S. S. Yazadjiev, N. Stergioulas, K. D. Kokkotas, and T. M. Athanasiadis (2014b), Phys. Rev. D 90 (4), 044004, arXiv:1405.6976 [astro-ph.HE].  \nDouchin, F., and P. Haensel (2001), Astron. Astrophys. 380 , 151, arXiv:astro-ph/0111092.  \nEast, W. E., and J. L. Ripley (2021a), Phys. Rev. Lett. 127 (10), 101102, arXiv:2105.08571 [gr-qc].  \n- East, W. E., and J. L. Ripley (2021b), Phys. Rev. D 103 (4), 044040, arXiv:2011.03547 [gr-qc].  \nElley, M., H. O. Silva, H. Witek, and N. Yunes (2022), Phys. Rev. D 106 (4), 044018, arXiv:2205.06240 [gr-qc].  \n- Erices, C., S. Riquelme, and N. Zalaquett (2022), Phys. Rev. D 106 (4), 044046, arXiv:2203.06030 [gr-qc].\n- Esposito-Far'ese, G. (2004), AIP Conf. Proc. 736 (1), 35, arXiv:gr-qc/0409081.\n- Esposito-Far'ese, G., C. Pitrou, and J.-P. Uzan (2010), Phys. Rev. D 81 , 063519, arXiv:0912.0481 [gr-qc].\n- Fernandes, P. G. S., C. A. R. Herdeiro, A. M. Pombo, E. Radu, and N. Sanchis-Gual (2019a), Phys. Rev. D 100 (8), 084045, arXiv:1908.00037 [gr-qc].  \nFernandes, P. G. S., C. A. R. Herdeiro, A. M. Pombo, E. Radu, and N. Sanchis-Gual (2019b), Class. Quant. Grav. 36 (13), 134002, [Erratum: Class.Quant.Grav. 37, 049501 (2020)], arXiv:1902.05079 [gr-qc].  \nFranchini, N., M. Bezares, E. Barausse, and L. Lehner (2022), Phys. Rev. D 106 (6), 064061, arXiv:2206.00014 [gr-qc].  \n- Franchini, N., A. Coates, and T. P. Sotiriou (2018), Phys. Rev. D 97 (6), 064013, arXiv:1708.02113 [gr-qc].\n- Franchini, N., and T. P. Sotiriou (2020), Phys. Rev. D 101 (6), 064068, arXiv:1903.05427 [gr-qc].\n- Freire, P. C. C. (2022), 'Pulsar mass measurements and tests of general relativity,' https://www3.mpifr-bonn.mpg.de/ staff/pfreire/NS\\_masses.html , accessed: 2021-12-28.  \nFreire, P. C. C., N. Wex, G. Esposito-Far'ese, J. P. W. Verbiest, M. Bailes, B. A. Jacoby, M. Kramer, I. H. Stairs, J. Anto-  \n- niadis, and G. H. Janssen (2012), Mon. Not. Roy. Astron. Soc. 423 , 3328, arXiv:1205.1450 [astro-ph.GA].\n- Friedman, J. L., and B. F. Schutz (1978), Astrophys. J. 222 , 281.\n- Gao, Y.-X., Y. Huang, and D.-J. Liu (2019), Phys. Rev. D 99 (4), 044020, arXiv:1808.01433 [gr-qc].\n- Garcia-Saenz, S., A. Held, and J. Zhang (2021), Phys. Rev. Lett. 127 (13), 131104, arXiv:2104.08049 [gr-qc].\n- Garcia-Saenz, S., A. Held, and J. Zhang (2022), JHEP 05 , 139, arXiv:2202.07131 [gr-qc].\n- Gendreau, K., and Z. Arzoumanian (2017), Nature Astronomy 1 , 895.\n- Gendreau, K. C., Z. Arzoumanian, and T. Okajima (2012), in Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray , Proc. SPIE, Vol. 8443, p. 844313.\n- Geng, C.-Q., H.-J. Kuan, and L.-W. Luo (2020), Eur. Phys. J. C 80 (8), 780, arXiv:2005.11629 [gr-qc].\n- Geroch, R. P. (1970), J. Math. Phys. 11 , 2580.\n- Gerosa, D., U. Sperhake, and C. D. Ott (2016), Class. Quant. Grav. 33 (13), 135002, arXiv:1602.06952 [gr-qc].\n- Gourgoulhon, E., P. Grandclement, K. Taniguchi, J.-A. Marck, and S. Bonazzola (2001), Phys. Rev. D 63 , 064029, arXiv:grqc/0007028.\n- Guo, M., J. Zhao, and L. Shao (2021), Phys. Rev. D 104 (10), 104065, arXiv:2106.01622 [gr-qc].\n- Hansen, R. O. (1974), J. Math. Phys. 15 , 46.\n- Harada, T. (1997), Prog. Theor. Phys. 98 , 359, arXiv:grqc/9706014.\n- Harada, T. (1998), Phys. Rev. D 57 , 4802, arXiv:grqc/9801049.\n- Hartle, J. B. (1967), Astrophys. J. 150 , 1005.\n- Hartle, J. B. (1978), Physics Reports 46 (6), 201.\n- Hartle, J. B., and K. S. Thorne (1968), Astrophys. J. 153 , 807.\n- Hassan, S. F., and R. A. Rosen (2012), JHEP 02 , 126, arXiv:1109.3515 [hep-th].\n- Hawking, S. W. (1972), Commun. Math. Phys. 25 , 167.\n- Heisenberg, L. (2014), JCAP 05 , 015, arXiv:1402.7026 [hepth].\n- Hellings, R. W., and K. Nordtvedt (1973), Phys. Rev. D 7 , 3593.\n- Herdeiro, C. A. R., T. Ikeda, M. Minamitsuji, T. Nakamura, and E. Radu (2021a), Phys. Rev. D 103 (4), 044019, arXiv:2009.06971 [gr-qc].\n- Herdeiro, C. A. R., and J. a. M. S. Oliveira (2020), JHEP 07 , 130, arXiv:2005.05354 [gr-qc].\n- Herdeiro, C. A. R., and E. Radu (2015), Int. J. Mod. Phys. D 24 (09), 1542014, arXiv:1504.08209 [gr-qc].\n- Herdeiro, C. A. R., and E. Radu (2019), Phys. Rev. D 99 (8), 084039, arXiv:1901.02953 [gr-qc].\n- Herdeiro, C. A. R., E. Radu, N. Sanchis-Gual, and J. A. Font (2018), Phys. Rev. Lett. 121 (10), 101102, arXiv:1806.05190 [gr-qc].\n- Herdeiro, C. A. R., E. Radu, H. O. Silva, T. P. Sotiriou, and N. Yunes (2021b), Phys. Rev. Lett. 126 (1), 011103, arXiv:2009.03904 [gr-qc].\n- Herrera, L., and N. O. Santos (1997), Phys. Rept. 286 , 53. Hilditch, D. (2013), Int. J. Mod. Phys. A 28 , 1340015, arXiv:1309.2012 [gr-qc].\n- Hod, S. (2020), Phys. Rev. D 102 (8), 084060, arXiv:2006.09399 [gr-qc].\n- Hod, S. (2022), Phys. Rev. D 105 (2), 024074.\n- Horbatsch, M., H. O. Silva, D. Gerosa, P. Pani, E. Berti, L. Gualtieri, and U. Sperhake (2015), Class. Quant. Grav.  \n32 (20), 204001, arXiv:1505.07462 [gr-qc].  \n- Horbatsch, M. W., and C. P. Burgess (2011), JCAP 08 , 027, arXiv:1006.4411 [gr-qc].\n- Horbatsch, M. W., and C. P. Burgess (2012), Class. Quant. Grav. 29 , 245004, arXiv:1107.3585 [gr-qc].  \nHorndeski, G. W. (1974), Int. J. Theor. Phys.  \n10  \n, 363.  \n- Hotokezaka, K., K. Kiuchi, K. Kyutoku, T. Muranushi, Y.-i. Sekiguchi, M. Shibata, and K. Taniguchi (2013), Phys. Rev. D 88 , 044026, arXiv:1307.5888 [astro-ph.HE].\n- Hu, Z., Y. Gao, R. Xu, and L. Shao (2021), Phys. Rev. D 104 (10), 104014, arXiv:2109.13453 [gr-qc].\n- Hulse, R. A., and J. H. Taylor (1975), Astrophys. J. Lett. 195 , L51.\n- Ikeda, T., T. Nakamura, and M. Minamitsuji (2019), Phys. Rev. D 100 (10), 104014, arXiv:1908.09394 [gr-qc].\n- Israel, G., T. Belloni, L. Stella, Y. Rephaeli, D. Gruber, P. G. Casella, S. Dall'Osso, N. Rea, M. Persic, and R. Rothschild (2005), Astrophys. J. Lett. 628 , L53, arXiv:astroph/0505255.\n- Jackiw, R., and S. Y. Pi (2003), Phys. Rev. D 68 , 104012, arXiv:gr-qc/0308071.\n- Jentschura, U. D., and B. J. Wundt (2012), J. Phys. A45 (44), 444017, arXiv:1110.4171 [hep-ph].\n- Juli'e, F.-L., and E. Berti (2019), Phys. Rev. D 100 (10), 104061, arXiv:1909.05258 [gr-qc].\n- Juli'e, F.-L., H. O. Silva, E. Berti, and N. Yunes (2022), Phys. Rev. D 105 (12), 124031, arXiv:2202.01329 [gr-qc].  \nKalogera, V., et al. (2021), 'The Next Generation Global Gravitational Wave Observatory: The Science Book,' arXiv:2111.06990 [gr-qc].  \nKanti, P., N. E. Mavromatos, J. Rizos, K. Tamvakis, and E. Winstanley (1996), Phys. Rev. D 54 , 5049, arXiv:hepth/9511071.  \n- Kase, R., M. Minamitsuji, and S. Tsujikawa (2020), Phys. Rev. D 102 (2), 024067, arXiv:2001.10701 [gr-qc].\n- Khalil, M., R. F. P. Mendes, N. Ortiz, and J. Steinhoff (2022), Phys. Rev. D 106 (10), 104016, arXiv:2206.13233 [gr-qc].  \nKhalil, M., N. Sennett, J. Steinhoff, and A. Buonanno (2019),  \nPhys. Rev. D  \n100  \n(12), 124013, arXiv:1906.08161 [gr-qc].  \n- Kleihaus, B., J. Kunz, S. Mojica, and M. Zagermann (2016), Phys. Rev. D 93 (6), 064077, arXiv:1601.05583 [gr-qc].\n- Kobayashi, T. (2019), Rept. Prog. Phys. 82 (8), 086901, arXiv:1901.07183 [gr-qc].\n- Kobayashi, T., M. Yamaguchi, and J. Yokoyama (2011), Prog. Theor. Phys. 126 , 511, arXiv:1105.5723 [hep-th].\n- Kokkotas, K. D., and B. G. Schmidt (1999), Living Rev. Rel. 2 , 2, arXiv:gr-qc/9909058.\n- Konno, K., T. Matsuyama, and S. Tanda (2009), Prog. Theor. Phys. 122 , 561, arXiv:0902.4767 [gr-qc].\n- Kov'acs, A. D. (2019), Phys. Rev. D 100 (2), 024005, arXiv:1904.00963 [gr-qc].\n- Kov'acs, A. D., and H. S. Reall (2020a), Phys. Rev. D 101 (12), 124003, arXiv:2003.08398 [gr-qc].\n- Kov'acs, A. D., and H. S. Reall (2020b), Phys. Rev. Lett. 124 (22), 221101, arXiv:2003.04327 [gr-qc].\n- Krall, V., A. Coates, and K. D. Kokkotas (2020), Phys. Rev. D 102 (12), 124065, arXiv:2012.03710 [gr-qc].\n- Kramer, M., et al. (2021), Phys. Rev. X 11 (4), 041050, arXiv:2112.06795 [astro-ph.HE].\n- Kruger, C. J., and D. D. Doneva (2021), Phys. Rev. D 103 (12), 124034, arXiv:2102.11698 [gr-qc].\n- Kruger, C. J., and K. D. Kokkotas (2020a), Phys. Rev. D 102 (6), 064026, arXiv:2008.04127 [gr-qc].\n- Kruger, C. J., and K. D. Kokkotas (2020b), Phys. Rev. Lett.  \n- 125 (11), 111106, arXiv:1910.08370 [gr-qc].\n- Kuan, H.-J., D. D. Doneva, and S. S. Yazadjiev (2021a), Phys. Rev. Lett. 127 (16), 161103, arXiv:2103.11999 [gr-qc].\n- Kuan, H.-J., J. Singh, D. D. Doneva, S. S. Yazadjiev, and K. D. Kokkotas (2021b), Phys. Rev. D 104 (12), 124013, arXiv:2105.08543 [gr-qc].\n- Kuan, H.-J., A. G. Suvorov, D. D. Doneva, and S. S. Yazadjiev (2022), Phys. Rev. Lett. 129 (12), 121104, arXiv:2203.03672 [gr-qc].\n- Landulfo, A. G. S., W. C. C. Lima, G. E. A. Matsas, and D. A. T. Vanzella (2012), Phys. Rev. D 86 , 104025, arXiv:1204.3654 [gr-qc].\n- Landulfo, A. G. S., W. C. C. Lima, G. E. A. Matsas, and D. A. T. Vanzella (2015), Phys. Rev. D 91 (2), 024011, arXiv:1410.2274 [gr-qc].\n- Lang, R. N. (2014), Phys. Rev. D 89 (8), 084014, arXiv:1310.3320 [gr-qc].  \nLanglois, D., K. Noui, and H. Roussille (2022), JCAP  \n09  \n, 019,  \narXiv:2204.04107 [gr-qc].  \n- Lattimer, J. M., and M. Prakash (2016), Phys. Rept. 621 , 127, arXiv:1512.07820 [astro-ph.SR].  \nLattimer, J. M., and B. F. Schutz (2005), Astrophys. J.  \n629  \n,  \n979, arXiv:astro-ph/0411470.  \n- Lee, D. L. (1974), Phys. Rev. D 10 , 2374.\n- Liebling, S. L., and C. Palenzuela (2023), Living Rev. Rel. 26 , 1, arXiv:1202.5809 [gr-qc].\n- Lima, W. C. C., G. E. A. Matsas, and D. A. T. Vanzella (2010), Phys. Rev. Lett. 105 , 151102, arXiv:1009.1771 [gr-qc].\n- Lima, W. C. C., R. F. P. Mendes, G. E. A. Matsas, and D. A. T. Vanzella (2013), Phys. Rev. D 87 (10), 104039, arXiv:1304.0582 [gr-qc].\n- Lima, W. C. C., and D. A. T. Vanzella (2010), Phys. Rev. Lett. 104 , 161102, arXiv:1003.3421 [gr-qc].\n- Lindblom, L. (2010), Phys. Rev. D 82 , 103011, arXiv:1009.0738 [astro-ph.HE].\n- Liu, Y., C.-Y. Zhang, Q. Chen, Z. Cao, Y. Tian, and B. Wang (2023a), Sci. China Phys. Mech. Astron. 66 (10), 100412, arXiv:2208.07548 [gr-qc].\n- Liu, Y., C.-Y. Zhang, W.-L. Qian, K. Lin, and B. Wang (2023b), JHEP 01 , 074, arXiv:2206.05012 [gr-qc].\n- Luo, W.-K., C.-Y. Zhang, P. Liu, C. Niu, and B. Wang (2022), Phys. Rev. D 106 (6), 064036, arXiv:2206.05690 [gr-qc].\n- Macedo, C. F. B. (2020), Int. J. Mod. Phys. D 29 (11), 2041006, arXiv:2002.12719 [gr-qc].\n- Macedo, C. F. B., J. Sakstein, E. Berti, L. Gualtieri, H. O. Silva, and T. P. Sotiriou (2019), Phys. Rev. D 99 (10), 104041, arXiv:1903.06784 [gr-qc].\n- Maione, F., R. De Pietri, A. Feo, and F. Loffler (2016), Class. Quant. Grav. 33 (17), 175009, arXiv:1605.03424 [gr-qc].\n- Maselli, A., N. Franchini, L. Gualtieri, and T. P. Sotiriou (2020), Phys. Rev. Lett. 125 (14), 141101, arXiv:2004.11895 [gr-qc].\n- Maselli, A., N. Franchini, L. Gualtieri, T. P. Sotiriou, S. Barsanti, and P. Pani (2022), Nature Astron. 6 (4), 464, arXiv:2106.11325 [gr-qc].\n- Mayo, A. E., and J. D. Bekenstein (1996), Phys. Rev. D 54 , 5059, arXiv:gr-qc/9602057.\n- McDermott, P. N., H. M. van Horn, and J. F. Scholl (1983), ApJ 268 , 837.\n- Mendes, R. F. (2015), Phys. Rev. D 91 (6), 064024, arXiv:1412.6789 [gr-qc].\n- Mendes, R. F., and N. Ortiz (2016), Phys. Rev. D 93 (12), 124035, arXiv:1604.04175 [gr-qc].\n- Mendes, R. F., and N. Ortiz (2018), Phys. Rev. Lett. 120 (20),  \n201104, arXiv:1802.07847 [gr-qc].  \n- Mendes, R. F., and T. Ottoni (2019), Phys. Rev. D 99 (12), 124003, arXiv:1903.11638 [gr-qc].\n- Mendes, R. F. P., G. E. A. Matsas, and D. A. T. Vanzella (2014a), Phys. Rev. D 90 (4), 044053, arXiv:1407.6405 [grqc].\n- Mendes, R. F. P., G. E. A. Matsas, and D. A. T. Vanzella (2014b), Phys. Rev. D 89 (4), 047503, arXiv:1310.2185 [grqc].\n- Miller, M. C., et al. (2019), Astrophys. J. Lett. 887 (1), L24, arXiv:1912.05705 [astro-ph.HE].\n- Minamitsuji, M. (2020a), Phys. Rev. D 101 (10), 104044, arXiv:2003.11885 [gr-qc].\n- Minamitsuji, M. (2020b), Phys. Rev. D 102 (4), 044048, arXiv:2008.12758 [gr-qc].\n- Minamitsuji, M. (2021), Phys. Rev. D 103 (8), 084002, arXiv:2104.03660 [gr-qc].\n- Minamitsuji, M., and T. Ikeda (2019a), Phys. Rev. D 99 (4), 044017, arXiv:1812.03551 [gr-qc].\n- Minamitsuji, M., and T. Ikeda (2019b), Phys. Rev. D 99 (10), 104069, arXiv:1904.06572 [gr-qc].\n- Minamitsuji, M., and H. O. Silva (2016), Phys. Rev. D 93 (12), 124041, arXiv:1604.07742 [gr-qc].\n- Minamitsuji, M., and S. Tsujikawa (2023), Phys. Lett. B 840 , 137869, arXiv:2208.08107 [gr-qc].\n- Mirshekari, S., and C. M. Will (2013), Phys. Rev. D 87 (8), 084070, arXiv:1301.4680 [gr-qc].\n- Morisaki, S., and T. Suyama (2017), Phys. Rev. D 96 (8), 084026, arXiv:1707.02809 [gr-qc].  \nMotohashi, H., and M. Minamitsuji (2018), Phys. Lett. B  \n781  \n,  \n728, arXiv:1804.01731 [gr-qc].  \n- Motohashi, H., and T. Suyama (2012), Phys. Rev. D 85 , 044054, arXiv:1110.6241 [gr-qc].\n- Mou, Z.-G., and H.-Y. Zhang (2022), Phys. Rev. Lett. 129 (15), 151101, arXiv:2204.11324 [hep-th].\n- Muther, H., M. Prakash, and T. L. Ainsworth (1987), Phys. Lett. B 199 , 469.\n- Myung, Y. S., and D.-C. Zou (2018), Phys. Rev. D 98 (2), 024030, arXiv:1805.05023 [gr-qc].\n- Myung, Y. S., and D.-C. Zou (2019a), Int. J. Mod. Phys. D 28 (09), 1950114, arXiv:1903.08312 [gr-qc].\n- Myung, Y. S., and D.-C. Zou (2019b), Eur. Phys. J. C 79 (3), 273, arXiv:1808.02609 [gr-qc].\n- Myung, Y. S., and D.-C. Zou (2019c), Eur. Phys. J. C 79 (8), 641, arXiv:1904.09864 [gr-qc].\n- Myung, Y. S., and D.-C. Zou (2021), Phys. Lett. B 814 , 136081, arXiv:2012.02375 [gr-qc].\n- Nelmes, S. G., and B. M. A. G. Piette (2012), Phys. Rev. D 85 , 123004, arXiv:1204.0910 [astro-ph.SR].\n- Niu, C., W. Xiong, P. Liu, C.-Y. Zhang, and B. Wang (2022), arXiv:2209.12117 [gr-qc].\n- Niu, R., X. Zhang, B. Wang, and W. Zhao (2021), Astrophys. J. 921 (2), 149, arXiv:2105.13644 [gr-qc].\n- Novak, J. (1998a), Phys. Rev. D 58 , 064019, arXiv:grqc/9806022.\n- Novak, J. (1998b), Phys. Rev. D 57 , 4789, arXiv:grqc/9707041.\n- Novak, J., and J. M. Ibanez (2000), Astrophys. J. 533 , 392, arXiv:astro-ph/9911298.\n- Ofengeim, D. D. (2020), Phys. Rev. D 101 (10), 103029, arXiv:2005.03549 [astro-ph.HE].\n- Oliveira, J. a. M. S., and A. M. Pombo (2021), Phys. Rev. D 103 (4), 044004, arXiv:2012.07869 [gr-qc].\n- Ozel, F., D. Psaltis, T. Guver, G. Baym, C. Heinke, and S. Guil-  \n- lot (2016), Astrophys. J. 820 (1), 28, arXiv:1505.05155 [astro-ph.HE].\n- Palenzuela, C., E. Barausse, M. Ponce, and L. Lehner (2014), Phys. Rev. D 89 (4), 044024, arXiv:1310.4481 [gr-qc].\n- Palenzuela, C., and S. L. Liebling (2016), Phys. Rev. D 93 (4), 044009, arXiv:1510.03471 [gr-qc].\n- Pani, P., and E. Berti (2014), Phys. Rev. D 90 (2), 024025, arXiv:1405.4547 [gr-qc].\n- Pani, P., E. Berti, V. Cardoso, and J. Read (2011a), Phys. Rev. D 84 , 104035, arXiv:1109.0928 [gr-qc].\n- Pani, P., V. Cardoso, E. Berti, J. Read, and M. Salgado (2011b), Phys. Rev. D 83 , 081501, arXiv:1012.1343 [gr-qc]. Papallo, G. (2017), Phys. Rev. D 96 (12), 124036, arXiv:1710.10155 [gr-qc].\n- Pappas, G., and T. A. Apostolatos (2014), Phys. Rev. Lett. 112 , 121101, arXiv:1311.5508 [gr-qc].\n- Pappas, G., D. D. Doneva, T. P. Sotiriou, S. S. Yazadjiev, and K. D. Kokkotas (2019), Phys. Rev. D 99 (10), 104014, arXiv:1812.01117 [gr-qc].\n- Pappas, G., and T. P. Sotiriou (2015), Phys. Rev. D 91 (4), 044011, arXiv:1412.3494 [gr-qc].\n- Paschalidis, V., and N. Stergioulas (2017), Living Rev. Rel. 20 (1), 7, arXiv:1612.03050 [astro-ph.HE].\n- Podkowka, D. M., R. F. P. Mendes, and E. Poisson (2018),\n- Phys. Rev. D 98 (6), 064057, arXiv:1807.01565 [gr-qc].\n- Ponce, M., C. Palenzuela, E. Barausse, and L. Lehner (2015), Phys. Rev. D 91 (8), 084038, arXiv:1410.0638 [gr-qc].\n- Popchev, D. (2015), Bifurcation of neutron star solutions in scalar-tensor theories of gravity , Master's thesis (University of Sofia).\n- Popchev, D., K. V. Staykov, D. D. Doneva, and S. S. Yazadjiev (2019), Eur. Phys. J. C 79 (2), 178, arXiv:1812.00347 [gr-qc]. Pretorius, F. (2005), Phys. Rev. Lett. 95 , 121101, arXiv:grqc/0507014.\n- Proca, A. (1936), J. Phys. Radium 7 , 347.\n- Ramazano˘glu, F. M. (2017), Phys. Rev. D 96 (6), 064009, arXiv:1706.01056 [gr-qc].\n- Ramazano˘glu, F. M. (2018a), Phys. Rev. D 97 (2), 024008, arXiv:1710.00863 [gr-qc].\n- Ramazano˘glu, F. M. (2018b), Phys. Rev. D 98 (4), 044013, arXiv:1804.03158 [gr-qc].\n- Ramazano˘glu, F. M. (2018c), Phys. Rev. D 98 (4), 044011, arXiv:1804.00594 [gr-qc].\n- Ramazano˘glu, F. M. (2019a), Phys. Rev. D 99 (4), 044003, arXiv:1901.00194 [gr-qc].\n- Ramazano˘glu, F. M. (2019b), Phys. Rev. D 99 (8), 084015, arXiv:1901.10009 [gr-qc].\n- Ramazano˘glu, F. M. (2019c), Turk. J. Phys. 43 (6), 586.\n- Ramazano˘glu, F. M., and F. Pretorius (2016), Phys. Rev. D 93 (6), 064005, arXiv:1601.07475 [gr-qc].\n- Ramazano˘glu, F. M., and K. I. Unluturk (2019), Phys. Rev. D 100 (8), 084026, arXiv:1910.02801 [gr-qc].\n- Regge, T., and J. A. Wheeler (1957), Phys. Rev. 108 , 1063. (12),\n- Rezzolla, L., and K. Takami (2016), Phys. Rev. D 93 124051, arXiv:1604.00246 [gr-qc].\n- de Rham, C. (2014), Living Rev. Rel. 17 , 7, arXiv:1401.4173 [hep-th].\n- de Rham, C., G. Gabadadze, and A. J. Tolley (2011), Phys. Rev. Lett. 106 , 231101, arXiv:1011.1232 [hep-th].\n- Rhoades, C. E., Jr., and R. Ruffini (1974), Phys. Rev. Lett. 32 , 324.\n- Ribeiro, C. C. H., and D. A. T. Vanzella (2020), Phys. Rev. Res. 2 (1), 013281, arXiv:1912.01971 [gr-qc].  \nRiley, T. E., et al. (2019), Astrophys. J. Lett. 887 (1), L21, arXiv:1912.05702 [astro-ph.HE].  \n- Ripley, J. L. (2022), Int. J. Mod. Phys. D 31 (13), 2230017, arXiv:2207.13074 [gr-qc].\n- Ripley, J. L., and F. Pretorius (2020), Class. Quant. Grav. 37 (15), 155003, arXiv:2005.05417 [gr-qc].\n- Rosca-Mead, R., C. J. Moore, M. Agathos, and U. Sperhake (2019), Class. Quant. Grav. 36 (13), 134003, arXiv:1903.09704 [gr-qc].\n- Rosca-Mead, R., C. J. Moore, U. Sperhake, M. Agathos, and D. Gerosa (2020a), Symmetry 12 (9), 1384, arXiv:2007.14429 [gr-qc].  \nRosca-Mead, R., U. Sperhake, C. J. Moore, M. Agathos, D. Gerosa, and C. D. Ott (2020b), Phys. Rev. D 102 (4), 044010, arXiv:2005.09728 [gr-qc].  \n- Ruegg, H., and M. Ruiz-Altaba (2004), Int. J. Mod. Phys. A 19 , 3265, arXiv:hep-th/0304245.\n- Ruiz, M., J. C. Degollado, M. Alcubierre, D. Nunez, and M. Salgado (2012), Phys. Rev. D 86 , 104044, arXiv:1207.6142 [gr-qc].\n- Saffer, A., H. O. Silva, and N. Yunes (2019), Phys. Rev. D 100 (4), 044030, arXiv:1903.07779 [gr-qc].\n- Salgado, M. (2006), Class. Quant. Grav. 23 , 4719, arXiv:grqc/0509001.\n- Salgado, M., D. Martinez-del Rio, M. Alcubierre, and D. Nunez (2008), Phys. Rev. D 77 , 104010, arXiv:0801.2372 [gr-qc].\n- Salgado, M., D. Sudarsky, and U. Nucamendi (1998), Phys. Rev. D 58 , 124003, arXiv:gr-qc/9806070.\n- Salgado, M., D. Sudarsky, and U. Nucamendi (2004), Phys. Rev. D 70 , 084027, arXiv:gr-qc/0402126.\n- Sampson, L., N. Yunes, N. Cornish, M. Ponce, E. Barausse, A. Klein, C. Palenzuela, and L. Lehner (2014), Phys. Rev. D 90 (12), 124091, arXiv:1407.7038 [gr-qc].\n- Sanchis-Gual, N., J. C. Degollado, C. Herdeiro, J. A. Font, and P. J. Montero (2016a), Phys. Rev. D 94 (4), 044061, arXiv:1607.06304 [gr-qc].\n- Sanchis-Gual, N., J. C. Degollado, P. J. Montero, J. A. Font, and C. Herdeiro (2016b), Phys. Rev. Lett. 116 (14), 141101, arXiv:1512.05358 [gr-qc].\n- Santiago, J., A. G. S. Landulfo, W. C. C. Lima, G. E. A. Matsas, R. F. P. Mendes, and D. A. T. Vanzella (2016), Phys. Rev. D 93 (2), 024043, arXiv:1512.02120 [gr-qc].\n- Sathyaprakash, B. S., et al. (2019), 'Extreme Gravity and Fundamental Physics,' arXiv:1903.09221 [astro-ph.HE].\n- Scheel, M. A., S. L. Shapiro, and S. A. Teukolsky (1995a), Phys. Rev. D 51 , 4208, arXiv:gr-qc/9411025.\n- Scheel, M. A., S. L. Shapiro, and S. A. Teukolsky (1995b), Phys. Rev. D 51 , 4236, arXiv:gr-qc/9411026.\n- Sen, S., and N. Banerjee (2001), Pramana 56 , 487, arXiv:grqc/9809064.\n- Sennett, N., and A. Buonanno (2016), Phys. Rev. D 93 (12), 124004, arXiv:1603.03300 [gr-qc].\n- Sennett, N., L. Shao, and J. Steinhoff (2017), Phys. Rev. D 96 (8), 084019, arXiv:1708.08285 [gr-qc].\n- Shao, L., N. Sennett, A. Buonanno, M. Kramer, and N. Wex (2017), Phys. Rev. X 7 (4), 041025, arXiv:1704.07561 [gr-qc]. Shibata, M., K. Taniguchi, H. Okawa, and A. Buonanno (2014), Phys. Rev. D 89 (8), 084005, arXiv:1310.0627 [gr-qc].\n- Shiralilou, B., T. Hinderer, S. Nissanke, N. Ortiz, and H. Witek (2021), Phys. Rev. D 103 (12), L121503, arXiv:2012.09162 [gr-qc].\n- Shiralilou, B., T. Hinderer, S. M. Nissanke, N. Ortiz, and H. Witek (2022), Class. Quant. Grav. 39 (3), 035002, arXiv:2105.13972 [gr-qc].  \n- Silva, H. O., A. Coates, F. M. Ramazano˘glu, and T. P. Sotiriou (2022), Phys. Rev. D 105 (2), 024046, arXiv:2110.04594 [grqc].\n- Silva, H. O., A. M. Holgado, A. C'ardenas-Avenda˜no, and N. Yunes (2021a), Phys. Rev. Lett. 126 (18), 181101, arXiv:2004.01253 [gr-qc].\n- Silva, H. O., C. F. Macedo, E. Berti, and L. C. Crispino (2015), Class. Quant. Grav. 32 , 145008, arXiv:1411.6286 [gr-qc].\n- Silva, H. O., C. F. B. Macedo, T. P. Sotiriou, L. Gualtieri, J. Sakstein, and E. Berti (2019), Phys. Rev. D 99 (6), 064011, arXiv:1812.05590 [gr-qc].\n- Silva, H. O., J. Sakstein, L. Gualtieri, T. P. Sotiriou, and E. Berti (2018), Phys. Rev. Lett. 120 (13), 131104, arXiv:1711.02080 [gr-qc].\n- Silva, H. O., H. Sotani, E. Berti, and M. Horbatsch (2014), Phys. Rev. D 90 (12), 124044, arXiv:1410.2511 [gr-qc].\n- Silva, H. O., H. Witek, M. Elley, and N. Yunes (2021b), Phys. Rev. Lett. 127 (3), 031101, arXiv:2012.10436 [gr-qc].\n- Silva, H. O., and N. Yunes (2019a), Class. Quant. Grav. 36 (17), 17LT01, arXiv:1902.10269 [gr-qc].\n- Silva, H. O., and N. Yunes (2019b), Phys. Rev. D 99 (4), 044034, arXiv:1808.04391 [gr-qc].\n- Soldateschi, J., N. Bucciantini, and L. Del Zanna (2020), Astron. Astrophys. 640 , A44, arXiv:2005.12758 [astro-ph.HE]. Soldateschi, J., N. Bucciantini, and L. Del Zanna (2021), Astron. Astrophys. 645 , A39, arXiv:2010.14833 [astro-ph.HE]. Sotani, H. (2012), Phys. Rev. D 86 , 124036, arXiv:1211.6986 [astro-ph.HE].\n- Sotani, H. (2014), Phys. Rev. D 89 (6), 064031, arXiv:1402.5699 [astro-ph.HE].\n- Sotani, H. (2017), Phys. Rev. D 96 (10), 104010, arXiv:1710.10596 [astro-ph.HE].\n- Sotani, H., and K. D. Kokkotas (2004), Phys. Rev. D 70 , 084026, arXiv:gr-qc/0409066.\n- Sotani, H., and K. D. Kokkotas (2005), Phys. Rev. D 71 , 124038, arXiv:gr-qc/0506060.\n- Sotani, H., and K. D. Kokkotas (2017), Phys. Rev. D 95 (4), 044032, arXiv:1702.00874 [gr-qc].\n- Sotiriou, T. P., and V. Faraoni (2012), Phys. Rev. Lett. 108 , 081103, arXiv:1109.6324 [gr-qc].\n- Sotiriou, T. P., and S.-Y. Zhou (2014a), Phys. Rev. Lett. 112 , 251102, arXiv:1312.3622 [gr-qc].\n- Sotiriou, T. P., and S.-Y. Zhou (2014b), Phys. Rev. D 90 , 124063, arXiv:1408.1698 [gr-qc].\n- Sperhake, U., C. J. Moore, R. Rosca, M. Agathos, D. Gerosa, and C. D. Ott (2017), Phys. Rev. Lett. 119 (20), 201103, arXiv:1708.03651 [gr-qc].\n- Staykov, K. V., and D. D. Doneva (2022), Phys. Rev. D 106 (10), 104064, arXiv:2209.01038 [gr-qc].\n- Staykov, K. V., D. D. Doneva, and S. S. Yazadjiev (2019), Astrophys. Space Sci. 364 (10), 178, arXiv:1902.09208 [grqc].\n- Staykov, K. V., D. D. Doneva, S. S. Yazadjiev, and K. D. Kokkotas (2015), Phys. Rev. D 92 (4), 043009, arXiv:1503.04711 [gr-qc].\n- Staykov, K. V., D. Popchev, D. D. Doneva, and S. S. Yazadjiev (2018), Eur. Phys. J. C 78 (7), 586, arXiv:1805.07818 [gr-qc]. Stefanov, I. Z., S. S. Yazadjiev, and M. D. Todorov (2007a), Phys. Rev. D 75 , 084036, arXiv:0704.3784 [gr-qc].\n- Stefanov, I. Z., S. S. Yazadjiev, and M. D. Todorov (2007b), Mod. Phys. Lett. A 22 , 1217, arXiv:0708.3203 [gr-qc].  \nStefanov, I. Z., S. S. Yazadjiev, and M. D. Todorov (2008),  \n- Mod. Phys. Lett. A 23 , 2915, arXiv:0708.4141 [gr-qc].\n- Strohmayer, T. E., and A. L. Watts (2005), Astrophys. J. Lett.\n- 632 , L111, arXiv:astro-ph/0508206.\n- Strohmayer, T. E., and A. L. Watts (2006), Astrophys. J. 653 , 593, arXiv:astro-ph/0608463.\n- Takami, K., L. Rezzolla, and L. Baiotti (2014), Phys. Rev. Lett. 113 (9), 091104, arXiv:1403.5672 [gr-qc].\n- Takami, K., L. Rezzolla, and L. Baiotti (2015), Phys. Rev. D 91 (6), 064001, arXiv:1412.3240 [gr-qc].\n- Taniguchi, K., and E. Gourgoulhon (2003), Phys. Rev. D 68 , 124025, arXiv:gr-qc/0309045.\n- Taniguchi, K., and M. Shibata (2010), Astrophys. J. Suppl. 188 , 187, arXiv:1005.0958 [astro-ph.SR].\n- Taniguchi, K., M. Shibata, and A. Buonanno (2015), Phys. Rev. D 91 (2), 024033, arXiv:1410.0738 [gr-qc].\n- Tasinato, G. (2014), JHEP 04 , 067, arXiv:1402.6450 [hep-th]. Taylor, J. H., and J. M. Weisberg (1982), Astrophys. J. 253 , 908.\n- Taylor, J. N., A. Wolszczan, and T. Damour (1993), Nature 355 , 132.\n- Thorne, K. S., and A. Campolattaro (1967), Astrophysical Journal , 591.\n- Tuna, S., K. I. Unluturk, and F. M. Ramazano˘glu (2022), Phys. Rev. D 105 (12), 124070, arXiv:2204.02138 [gr-qc]. Ventagli, G., G. Antoniou, A. Leh'ebel, and T. P. Sotiriou (2021), Phys. Rev. D 104 (12), 124078, arXiv:2111.03644 [gr-qc].\n- Ventagli, G., A. Leh'ebel, and T. P. Sotiriou (2020), Phys. Rev. D 102 (2), 024050, arXiv:2006.01153 [gr-qc].\n- Voisin, G., I. Cognard, P. C. C. Freire, N. Wex, L. Guillemot, G. Desvignes, M. Kramer, and G. Theureau (2020), Astron. Astrophys. 638 , A24, arXiv:2005.01388 [gr-qc].\n- Wald, R. M. (1993), Phys. Rev. D 48 (8), R3427, arXiv:grqc/9307038.\n- Watts, A. L., et al. (2016), Rev. Mod. Phys. 88 (2), 021001, arXiv:1602.01081 [astro-ph.HE].\n- Weih, L. R., E. R. Most, and L. Rezzolla (2018), Mon. Not. Roy. Astron. Soc. 473 (1), L126, arXiv:1709.06058 [gr-qc]. Wex, N. (2014), 'Testing Relativistic Gravity with Radio Pulsars,' arXiv:1402.5594 [gr-qc].\n- Whinnett, A. (2000), Phys. Rev. D 61 , 124014, arXiv:grqc/9911052.\n- Whinnett, A. W. (1999), Class. Quant. Grav. 16 , 2797.\n- Whinnett, A. W., and D. F. Torres (2004), Astrophys. J. Lett. 603 , L133, arXiv:astro-ph/0401521.\n- Will, C. M. (2014), Living Rev. Rel. 17 , 4, arXiv:1403.7377 [gr-qc].\n- Will, C. M. (2018), Theory and Experiment in Gravitational Physics (Cambridge University Press).\n- Witek, H., L. Gualtieri, P. Pani, and T. P. Sotiriou (2019),  \n- Phys. Rev. D 99 (6), 064035, arXiv:1810.05177 [gr-qc].\n- Wong, L. K., C. A. R. Herdeiro, and E. Radu (2022), Phys.\n- Rev. D 106 (2), 024008, arXiv:2204.09038 [gr-qc].\n- Xiong, W., P. Liu, C. Niu, C.-Y. Zhang, and B. Wang (2022), Chin. Phys. C 46 (9), 095103, arXiv:2205.07538 [gr-qc].\n- Xu, R., Y. Gao, and L. Shao (2020), Phys. Rev. D 102 (6), 064057, arXiv:2007.10080 [gr-qc].\n- Yagi, K., K. Kyutoku, G. Pappas, N. Yunes, and T. A. Apostolatos (2014), Phys. Rev. D 89 (12), 124013, arXiv:1403.6243 [gr-qc].\n- Yagi, K., L. C. Stein, N. Yunes, and T. Tanaka (2012), Phys. Rev. D 85 , 064022, [Erratum: Phys.Rev.D 93, 029902 (2016)], arXiv:1110.5950 [gr-qc].\n- Yagi, K., and M. Stepniczka (2021), Phys. Rev. D 104 (4), 044017, arXiv:2105.01614 [gr-qc].\n- Yagi, K., and N. Yunes (2013a), Science 341 , 365, arXiv:1302.4499 [gr-qc].\n- Yagi, K., and N. Yunes (2013b), Phys. Rev. D 88 (2), 023009, arXiv:1303.1528 [gr-qc].\n- Yagi, K., and N. Yunes (2017), Phys. Rept. 681 , 1, arXiv:1608.02582 [gr-qc].\n- Yamashita, Y., A. De Felice, and T. Tanaka (2014), Int. J. Mod. Phys. D 23 , 1443003, arXiv:1408.0487 [hep-th]. Yazadjiev, S. (1999), Class. Quant. Grav. 16 , L63, arXiv:grqc/9906038.\n- Yazadjiev, S. S., D. D. Doneva, and K. D. Kokkotas (2017), Phys. Rev. D 96 (6), 064002, arXiv:1705.06984 [gr-qc].\n- Yazadjiev, S. S., D. D. Doneva, and D. Popchev (2016), Phys. Rev. D 93 (8), 084038, arXiv:1602.04766 [gr-qc].\n- Yunes, N., and F. Pretorius (2009), Phys. Rev. D 79 , 084043, arXiv:0902.4669 [gr-qc].\n- Yunes, N., and X. Siemens (2013), Living Rev. Rel. 16 , 9, arXiv:1304.3473 [gr-qc].\n- Yunes, N., and L. C. Stein (2011), Phys. Rev. D 83 , 104002, arXiv:1101.2921 [gr-qc].\n- Zenginoglu, A. (2008), Class. Quant. Grav. 25 , 145002, arXiv:0712.4333 [gr-qc].\n- Zhang, C.-Y., Q. Chen, Y. Liu, W.-K. Luo, Y. Tian, and B. Wang (2022), Phys. Rev. D 106 (6), L061501, arXiv:2204.09260 [gr-qc].\n- Zhang, S.-J. (2021), Eur. Phys. J. C 81 (5), 441, arXiv:2102.10479 [gr-qc].\n- Zhang, S.-J., B. Wang, A. Wang, and J. F. Saavedra (2020),\n- Phys. Rev. D 102 (12), 124056, arXiv:2010.05092 [gr-qc].\n- Zhao, J., P. C. C. Freire, M. Kramer, L. Shao, and N. Wex (2022), Class. Quant. Grav. 39 (11), 11LT01, arXiv:2201.03771 [astro-ph.HE].\n- Zumalac'arregui, M., and J. Garc'ıa-Bellido (2014), Phys. Rev. D 89 , 064046, arXiv:1308.4685 [gr-qc]."}