The ghost of vector fields in compact stars

Spontaneous scalarization is a mechanism that allows a scalar field to go undetected in weak gravity environments and yet develop a nontrivial configuration in strongly gravitating systems. At the perturbative level it manifests as a tachyonic instability around spacetimes that solve Einstein's equations. The endpoint of this instability is a nontrivial scalar field configuration that can significantly modify a compact object's structure and can produce observational signatures of the scalar field's presence. Does such a mechanism exists for vector fields? Here we revisit the model that constitutes the most straightforward generalization of the original scalarization model to a vector field and perform a perturbative analysis. We show that a ghost appears as soon as the square of the naive effective mass squared becomes negative anywhere. This result poses a serious obstacle in generalizing spontaneous scalarization to vector fields.


I. INTRODUCTION
The first gravitational wave (GW) signal from a compact binary coalescence detected by the LIGO-Virgo collaboration [1] in 2015 opened a new vista into the nonlinear and highly dynamical regime of gravity. Moreover, and perhaps more excitingly, GWs now allow us to probe (or constrain) new physics beyond GR and the Standard Model [2][3][4][5][6][7]. This had so far been limited to astronomical probes either in the weak gravitational field and slow velocity in our Solar System or in the strong gravitational field, but small velocity and large separation regime of binary pulsars [8,9].
In this context, a particularly appealing new physics scenario is one where new fundamental fields lie "dormant" in weak-gravity environments and yet manage to have significant effects in strongly-gravitating bodies and systems. The prototypical theory that achieves this was first introduced by Damour and Esposito-Farèse [10,11] and involves a massless scalar field ϕ. The theory is described by the action where g is the metric determinant, R is the Ricci scalar, S m is the action of matter fields Ψ m , which couple to Ω 2 DEF (ϕ) g µν , with Ω DEF = exp βϕ 2 /2 ; β being a dimensionless constant.
The scalar field satisfies the field equation whereT is the trace of the matter field's energymomentum tensor. Equation (2) clearly admits a vanishing scalar field as a solution. However this is not the only solution for a given matter configuration. Linearized scalar field perturbations δϕ on the background of a neutron star can be shown to obey a wave equation where µ 2 eff is a position dependent effective mass squared. For a neutron star described by a perfect fluid,T = 3p −ε (wherep is the pressure andε the total energy density). TypicallyT < 0 and thus these perturbations can become tachyonic when β < 0 [12,13], with only a weak dependence on the equation of state [14][15][16]. Numerical simulations show that this linear instability is ultimately nonlinearly quenched and thus the star becomes spontaneously scalarized. Due to Eq. (2), these scalarized stars coexist with the GR solutions [defined as stars with ϕ = 0] and, importantly, are energetically favored: thus they can form dynamically from stellar collapse [14,[17][18][19][20] or in neutron star binaries [21][22][23][24][25][26].
The Damour-Esposito-Farèse scalarization model cannot lead to black holes scalarization unless the latter is induced by surrounding matter [27,28]. However, more general models which fashion couplings with the Gauss-Bonnet invariant, have been shown to lead to black hole scalarization, controlled by the mass [29,30] or by the spin of the black hole [31][32][33] and can take place in stellar collapse [34]. Black hole scalarization can also have potentially observable effects in binary black hole binaries [35][36][37] and be induced by other curvature scalars, such as the Pontryagin invariant [38,39]. The instability leading to scalarization can also be understood from a quantum field theory perspective, see e.g., Refs. [40][41][42].
A different type of generalization of the Damour-Esposito-Farèse mechanism that has been explored is to extend it to vector fields [43]. Inspired by [10,11], arXiv:2110.04594v3 [gr-qc] 20 Jan 2022 Ref. [43] studied the action where F αβ = ∇ α A β −∇ β A α is the antisymmetric Faraday tensor and A α is a vector field with bare mass µ v . 1 In analogy with Ω DEF , the conformal factor is chosen as where β is a free parameter of the theory. The field equation of A µ is This equation promotes the bare mass µ v of the Proca field to what appears to be an effective mass squared for linearized vector field perturbations. This effective mass squared can become negative in the presence of dense matter as in the theory (1). This property is not specific to the theory (4), and is shared with other vectortensor theories with curvature coupling terms [45][46][47][48] or disformal couplings [49,50]. Based on the similarity between the field equations (2) and (5) it is natural to expect that in the theory (4), A µ could also become tachyonically unstable around sufficiently compact neutron stars and a spontaneous vectorization mechanism exists. Although nonlinear vectorized neutron star solutions have indeed been shown to exist in [43], the perturbative manifestation of vectorization has not been explored yet. This leaves a number of open questions unanswered. In particular, a massive vector A µ is known to propagate an additional longitudinal degree of freedom. What is its role in this process? Could vectorization be scalarization in disguise to some extent? More generally, can it be understood intuitively, as is the case for scalarization, as a tachyonic instability quenched by nonlinearities? Answering these questions is important from a model-building perspective, but also from a phenomenological perspective. They become even more pressing once one observes that, intuitively speaking, the aforementioned longitudinal mode gets contributions in its kinetic term from the A µ A µ terms in the action. That kinetic term will therefore have a nontrivial structure, which in turn raises doubts about whether this mode is well behaved.
Motivated by these questions, here we revisit the model of Ref. [43] from a perturbative perspective and indeed uncover a ghost instability. Therefore vectorization appears to be fundamentally different from scalarization. It also strongly suggests that the time-evolution problem of a star undergoing vectorization is potentially ill-posed, casting serious doubts on the viability of this theory and other related ones. Combined with the work [51] which also found ghost (and gradient) instabilities in generalized Proca theories in compact object backgrounds, our work raises serious questions about the possibility to generalize the original mechanism of Damour and Esposito-Farèse beyond scalars since all proposed vectorization theories feature at least ghost instabilities.
The remainder of this paper explains how we arrived at these conclusions. In Secs. II and III we review the model introduced in [43], restore gauge invariance by performing the Stuckelberg trick and analyse the resulting field equations. In Sec. IV we linearize the theory's action in the background of a nonrotating, spherically symmetric star and show how ghost instability appears. In Sec. V we lift the assumption of linearized gauge field perturbations and consider the complete set of field equations. We show how ghosts, which first went unnoticed in [43], arise. In Sec. VI we summarize our main results.

II. A MODEL FOR SPONTANEOUS VECTORIZATION WITH GAUGE SYMMETRY
Action (4) has been constructed in analogy with (1), but a caveat of the resulting tensor-vector theory is absence of gauge invariance under A α → A α + ∂ α λ (λ being a scalar function) due to the mass term µ 2 v A µ A µ . To restore gauge invariance, and at the same time more easily investigate the different degrees of freedom in the vector field, we can apply the Stueckelberg trick [52]. It consists of introducing a scalar field ψ (the Stueckelberg field) through the substitution which results in a scalar-vector-tensor theory, with conformal factor The theory is now gauge invariant under the simultaneous transformations: We see that ψ can be set to zero by a suitable choice of λ and thus the action (4) is a gauge-fixed version of action (8).
Indeed, for β = 0 the conformal factor Ω becomes unity and we recover the Stueckelberg theory minimally coupled to gravity (see e.g., [53]). If we fix a gauge where ψ = 0 (we call this the "Proca gauge"), we obtain the nonminimally coupled Einstein-Proca theory of Ref. [43]. If we instead take µ v → 0 we obtain the Einstein-Maxwell theory with the addition of a scalar field. In Proca theory the µ v → 0 limit has an apparent discontinuity of the longitudinal polarization mode of A α . In the "Stueckelberged" version of the same theory, the µ v → 0 limit is manifestly continuous and corresponds to the decoupling between ψ and A α (the latter associated with the usual Maxwell theory). Note that, when β = 0, maintaining regularity of Ω requires that β approaches zero at least as fast as µ 2 v when taking the limit µ v → 0. The action (8) is written in the Einstein frame (thus we call g αβ the Einstein frame metric). We will refer tõ g αβ = Ω 2 g αβ as the Jordan frame metric. We will use tildes to denote objects in the Jordan frame, some of which, asT , already appeared in the Introduction.

III. THE FIELD EQUATIONS
The field equations of the theory can be obtained by varying the action (8) with respect to ψ, A α and g αβ : where, and we defined the individual energy-momentum contributions from "pure electromagnetic" theory T e µν and from the "Stueckelberg contribution" to the action T s µν , The Jordan frame energy-momentum tensor of matter fields and its trace are defined as We also have by construction: Going back to Eq. (12) and due to ∇ µ ∇ ν F µν = 0, it is convenient to define a current j α as: which is conserved In terms of j α we have: In the absence of matter (T = 0) and in the Proca gauge (ψ = 0), Eq. (21) becomes the Lorenz constraint on A α of Proca theory. Thus, the field equation for ψ [cf. Eq. (11)] and the Lorenz constraint on A µ are tightly connected.
We can see the first sign of the ghost by introducing a third metric, in terms of which the scalar field equation becomes This third metric can, in principle, have a signature change in some parts of the spacetime due to theẑ −1 term. If this happens, the field will be a ghost in at least some region compared to any field which is coupled to a fixed signature metric. Another potential problem is the fact that this metric changes sign by diverging, rather than crossing zero, in a similar vein discussed in [54,55]. It is unclear whether there is a rectification for such a problem, or, worse, whether the theory can evolve from a state where this metric has a fixed signature to another where the signature changes. It is also instructive to consider the limit µ v = 0, with β → 0 as fast as µ 2 v . In this limit, the Stueckelberg field ψ is no longer affected by gauge transformations, so A µ becomes a gauge field. Then A µ smoothly decouples from ψ and the matter fields. However, there is still coupling to gravity and ψ continues to be coupled to matter. In particular, Eq. (23) becomes, note thatẑ does not depend on A α here. So, ψ will become a ghost when theḡ µν metric changes signature and, as it is coupled to gravity and matter, its ghostly nature is physical. This same procedure is used in the Stueckelberg picture of Proca theory to show that ψ and A µ decouple and hence there is no discontinuity as µ v = 0 (i.e., no degree of freedom disappears). In this setting, one has β = 0, flat spacetime, and no matter. One may object that ψ can be completely removed by a gauge choice such as the Proca gauge ψ = 0, and thus the ghost can be exorcised. For this reason we will use the rest of the paper to assuage any doubts. We will begin by examining the quadratic Lagrangian for scalar-vector perturbations around a neutron star GR solution. Doing so we will find there exists a gauge invariant scalar field that suffers the same problems.

A. Background spacetime and overview of the calculation
In this section we explore the test-field limit of our theory, where we study the dynamics of ψ and A α in a background corresponding to a stellar solution of Einstein's field equations, i.e., a solution of the Tolman-Oppenheimer-Volkoff (TOV) equations [56,57] whose line element we write as where ν (lapse) and µ (mass function) are functions of the radial coordinate r only. In Sec. IV B we will linearize the field equations for small field perturbation δψ and δA µ at the level of the field equations (11)- (12), and show how the ghost arises in this background. Then, in Sec. IV C, we reach the same conclusion by directly perturbing the Lagrangian, by expanding it to second order in the fields on the same background.

B. Linearized field equations
We are interested in studying the dynamics of δA α and δψ propagating on the background line element (25). To proceed we decompose δψ and δA α in scalar and vector harmonics respectively. This is the convenient basis to expand scalar and vector fields on the unit two-sphere and, thus, in problems with spherical symmetry. We follow closely the presentation by Rosa and Dolan [58], but with a slightly different normalization. More specifically, we write δψ as where Y m = Y m (θ, φ) are the spherical harmonics with = 0, 1, 2 . . . , and |m| . For the vector perturbations, we decompose δA α as where c 1 = c 2 = 1, c 3 = c 4 = 1/ ( + 1), and Z (i) m α are the vector harmonics given by These functions are orthonormal when integrated on the unit two-sphere, according to the inner product, Under parity inversion x → −x (or equivalently, in spherical coordinates, θ → π − θ and φ → φ + π), the first three harmonics (i = 1, 2, 3) pick a factor of (−1) , while the fourth (i = 4) picks a factor of (−1) +1 . We follow the literature convention and call the former "even parity" modes and the latter "odd parity" modes. The scalar perturbation δψ is of even parity.
At this point it will be useful to follow a similar procedure to [51]. We expand the Stueckelberged action (8) around a GR solution to second order in the test field approximation and find, where, which is unity outside the star, whereT = 0. Note that we could have arrived at an action of this form by using the Stueckelberg trick in the Proca Lagrangian with a "dressed mass" zµ 2 v . Therefore, the results of this section apply to any theory whose quadratic Lagrangian can be put in this form, i.e., where one would naively expect just a screened Proca field prone to develop a tachyonic instability. Substituting the decompositions of δA α and ψ in harmonics, results in a Lagrangian, with even and odd-parity sector decoupled from one another. We look at each sector next.

C. Monopolar even-parity quadratic Lagrangian
We first focus on the monopole perturbations ( = m = 0), which have the lowest instability threshold and belong to the even-parity sector. Since Y 00 = constant, only the i = 1, 2 vector harmonics are defined [58]. This means that we would need to work with three variables σ 00 , u where, to shorten the notation, we use σ = σ 00 , u 1 = u 00 , and u 2 = u (2) 00 hereafter.
Inserting Eqs. (35) in the action (33) and integrating over the angular coordinates leaves us with, where we defined (·) = ∂ r (·) and( ·) = ∂ t (·). It can be readily verified that under the gauge transformation (10) with λ = l/(2 √ πr) that, and that the action (36) is invariant under this transformation. In fact, it can be verified that, the combination, is itself gauge invariant (proportional to the = 0 component of the electric field). If we introduce the auxiliary field φ such that, on shell, φ = r 2 Φ, we can rewrite Eq. (36) as, In this formulation, φ, u 1 and u 2 are all nondynamical: their equations of motion can be solved algebraically in the form, e.g., We can then replace this solution directly into the action, "integrating out" whichever field. Integrating out u 1 and u 2 one arrives at an action that is a functional of φ alone (all terms involving σ cancel). This transfers all of the dynamics from σ to φ. The resulting action has the form, where We see immediately that the sign of the kinetic contribution changes if z does (and also diverges when z crosses 0). That is, we have shown that, in this situation, there is a gauge invariant statement of the problems discussed in Sec. II, arising from Eq. (23).

D. Odd-parity quadratic Lagrangian
Having identified the presence of a ghost in the evenparity sector, it is natural to ask whether such ghosts also arise in the odd-parity sector, which contains a single degree of freedom u 4 , with multipole 1. We find, after integration over the angular coordinates, where we defined u 4 = u (4) 0 and set m = 0 due to the background's spherical symmetry.
Hence, we see that u 4 is prone to a tachyonic instability controlled by the same effective mass squared zµ 2 v also responsible for inducing a ghost instability in the evenparity sector. Indeed, the term between square brackets is the effective potential for massive vector axial perturbations found in [58], Eq. (13), for z = 1. We then conclude that the axial sector can become tachyonic unstable, but the dominant effect occurs at lower multipole: the ghost instability in the even-sector.

V. UNVEILING THE GHOST IN THE PROCA GAUGE
We have identified a ghost instability in the scalar sector of our theory, however no ghosts were reported in the spontaneous vectorization theory introduced in Ref. [43], or related theories investigated in Refs. [45][46][47][48]. In this section and related appendices, we will demonstrate that these theories contain divergent terms in their field equations irrespective of whether one uses the Stueckelberg trick to restore gauge symmetry or not.
Recall that the Proca gauge (ψ = 0) is equivalent to the spontaneous vectorization theory of Ref. [43]. Effectively, this gauge undoes the Stueckelberg trick (7) and we only need to consider Eq. (5). Since there is no separate equation for ψ in this picture and there are no divergent terms in this field equation, it is unclear where the ghost lurks. This is elucidated by considering the constraint equation.
Since ∇ µ ∇ ν F µν = 0 still holds due to the antisymmetry of F µν in Eq. (5), we obtain This is the generalized version of the ∇ µ A µ = 0 constraint for a minimally coupled Proca field. The puzzling aspect of Eq. (5) is that it does not have any explicit indication of a ghost, however we now know from our discussion in Sec. III that the constraint (46) given in the form of a conserved current in Eq. (20) is also crucial to understand the time evolution. Indeed, the constraint imposes a time evolution for A 0 that will reveal the ghost. 2 Let us rewrite the constraint in terms ofẑ [cf. Eq. (6)], We can convert the covariant derivatives to partial derivatives to obtain where i runs over the spatial coordinates. We see that this time-evolution equation has divergent terms due to the behavior ofẑ, even if all fields other than A 0 are regular. Outside any matter distributionẑ = 1 and we requireẑ < 0 in some part of spacetime if we want an astrophysical object to vectorize. Sinceẑ is continuous, it has to vanish at some point. There is no symmetry to ensure that √ −gẑA 0 vanishes whereẑ vanishes sinceẑ and its derivatives do not vanish at the same spacetime points in general. This means, A 0 will generically diverge even if √ −gẑA 0 stays regular. Alternatively, we can move theẑ term outside the derivative on the left-hand side, which means that now the coefficient of the leading time derivative of A 0 vanishes at certain points. This means that the divergent terms we observed in the ghost instabilities of ψ manifest themselves not directly in the field equation (5), but in the constraint equation (46), or equivalently, in Eq. (48).
The dynamics of A 0 implied by Eq. (48) is first order in time, thus not strictly of the same nature of the wave equation obeyed by ψ. Nonetheless, the change of sign in the time derivative leads to an analogous pathology. This can be understood by recasting the field equation (5) into an explicitly hyperbolic form.
Let us start by rearranging the constraint (47) as Next, we manipulate Eq. (5) as followŝ 2 Note that A 0 is not a dynamical degree of freedom in the standard Hamiltonian sense [59]. where we related the commutator of two covariant derivatives to the Riemann tensor in the third line, and used the constraint equation (49) in the fourth line. We finally obtain where we defined the mass-squared tensor We should be cautious about the fact thatẑ contains A α terms [inside the conformal factor; cf. Eq. (6)], which, strictly speaking, means that ∇ α ∇ β ln |ẑ| also belongs to the principal part of the differential equation. However, for perturbative values of A α , such as in a fixed background calculation of Sec. IV, this dependence can be ignored to leading order and M αβ becomes a proper mass-square tensor. Hence, Eq. (51) can be viewed as a generalized massive wave equation. Equation (51) has a divergent mass term due to various factors ofẑ −1 on its right-hand side. We have the option of moving these factors to the left-hand side, which means the principal part becomesẑ A µ . This is a field equation prone to a ghost instability sinceẑ changes sign as we discussed before in Eq. (23). One can also analyze the equation of motion for each vector harmonic, which likewise leads to divergent effective mass terms.
The behavior ofẑ is slightly modified for a vector field with no intrinsic mass, µ v = 0. In this case Eq. (6), and correspondingly Eq. (47), are modified aŝ where we assume the neutron star matter to behave as a perfect fluid with Jordan frame total energy densityε and pressurep as before. We see thatẑ vanishes outside the star and is generally negative within it; thus it never crosses zero. However, there are still divergences. The first case of the divergence in the field equations for µ v = 0 occurs at the surface of the neutron star. The relevant part of the TOV equations for a spherically symmetric star is [56,57], In the outer layers of the star one hasp ε and 4πpr 3 µ [60,61], which allows us to approximate Eq. (54) as where we approximated the total energy-density as equal to the rest mass density (ε ≈ρ), introduced the proper radial length [related to the coordinate radius r as d /dr = (1 − 2µ/r) −1/2 ], and defined the "local gravitational acceleration" g = (µ/r 2 )(1 − 2µ/r) −1/2 [60]. Focusing on the outer envelope of the star [62], we can approximate the spacetime as being Schwarzschild, i.e., µ ≈ M and ν = ln(1 − 2M/r s ) in Eq. (25), where M is the mass and r s the radius of the star. We can further introduce a local proper depth ÿ = (R−r)(1−2M/r s ) −1/2 , in terms of which we can recast Eq. (55) as, dp dÿ = g sρ , i.e., the equation of a plane-parallel atmosphere with a relativistic-corrected surface gravity g s = g(r s ) = (M/r s )(1 − 2M/r s ) −1/2 . In the outermost stellar layers, the main contribution to the pressure is due to a nonrelativistic degenerate electron gas, for which Eq. (56) can be solved exactly (see Ref. [60], Sec. 6.9), yielding the scalingρ ∝ ÿ 3/2 and, within the same assumptions, T ∝ −4πβ ÿ 3/2 . This, in turn, means that ∇ µ ln |ẑ| and the effective mass diverge on the surface, completing our argument. The same reasoning can in principle be applied to other systems which have an interface of vacuum and matter, suggesting that any such interfaces would lead to a divergence in the vector field equations in general. These divergences at the surface of the star are not exclusive to the vector-tensor model considered here, but are known to also arise, albeit with a different origin, in Palatini f (R) [63][64][65] and in Eddington-inspired Born-Infeld [66] theories. See also [67,68]. The second case of divergence in the field equations for µ v = 0 is related to massive neutron stars. AlthoughT is negative in general, it can switch sign and become positive in the core of such stars for some equations of state (see e.g., [69][70][71][72]). This means thatẑ vanishes somewhere inside the star [cf. Eq. (53)], where our previous results for the µ v = 0 case directly apply.
Overall, the above discussion provides a heuristic tool to identify ghosts in spontaneous vectorization theories. If the spacetime dependent µ 2 eff vanishes in nonvacuum regions in a theory with field equation ∇ µ F µα = µ 2 eff A α , this generically leads to divergent terms in the explicitly hyperbolic field equations. In other words, despite the appearances and the naming we used, µ eff is not the effective mass of all physical degrees of freedom. A careful analysis reveals that the true effective mass diverges as in Eq. (51), which was overlooked in the original spontaneous vectorization theory of Ref. [43] and other similar theories. We work this out explicitly in Appendix A (for the Hellings-Nordtvedt vector-tensor theory [73,74] studied in [45]) and in Appendix B (for the vector-Gauss-Bonnet theory of [46,48]).

VI. CONCLUSIONS
We revisited the tensor-vector gravity model proposed in Ref. [43] and explored the vectorization process using perturbation theory. This was done by working with a gauge invariant, Stueckelberg version of the theory and complemented with an analysis of the Lorenz constraint in the Proca gauge. In analogy with scalarization, one would expect to see the vector field develops a tachyonic instability, which is then quenched nonlinearly, and this process gives rise to the vectorized configurations found in previous work. Instead, we have uncovered a ghost instability. This results demonstrates quite clearly that the strong resemblance of this model of vectorization to the Damour-Esposito-Farèse model of scalarization is in fact rather misleading and a phase transition process that is physically similar to scalarization does not take place.
A potential way out may exist if one can tame the ghost instability nonlinearly, similar to the quenching of the tachyonic instability in scalarization. Indeed, "ghost-based spontaneous tensorization" has been investigated [75]. In the vectorization model studied here, ghosts appear inadvertently, and there is no explicit derivative coupling before the introduction of the Stueckelberg mechanism. Yet, if a nonlinear quenching mechanism exists, it could, in principle, suppress the ghost. Note that theẑ term in Eq. (47) that controls the instability approaches its GR value ofẑ = 1 when A µ A µ → ∞ (for β < 0). Hence, a solution with large vector field values can lead to a case whereẑ > 1 everywhere. This possibility was recently investigated for action (4) in Ref. [76], and all computed static and spherically symmetric vectorized neutron stars were shown to still carry ghost or gradient instabilities. Hence, there is no sign of a quenching of the instabilities so far.
The main issue however with the ghost instabilities we investigated is that it is not known whether their time evolution can be done. Even if a vector field growing to large values might quench the ghost, it is not clear if the very time evolution of the vector field that leads to growth can be formulated as a well-posed initial value problem due to the divergent terms such as those in Eq. (51). The resolution of this issue requires a mathematical analysis of the partial differential equations we have, which is beyond the scope of this work. We remark that these are not problems in the Proca limit of our model and in the absence of matter, in which numerical relativity simulations have been performed, e.g., in Refs. [77][78][79].
Spontaneous vectorization theories with restored gauge symmetry were also conceived using the Higgs mechanism rather than the Stueckelberg mechanism [80], inspired by the gravitational Higgs mechanism [81][82][83]. However, this theory [80] also has divergent terms in its field equations akin to Eq. (50), hence, it is susceptible to the same ill-posedness problems we discussed here.
We worked on the specific theory of Eq. (5), but other spontaneous vectorization models in the literature have similar field equations where ∇ µ F µα directly appears as the principal part [44][45][46][47][48][49]. Hence, a constraint can be obtained the same way as we did, which leads to divergent terms using the arguments in Sec. V or related ones, as we show in Appendices A and B.
Lastly, we stress that our results are relevant for most known extensions of spontaneous scalarization to other fields, not just the vectors, and our study can be considered as a first step to obtain a no-go theorem for extending spontaneous scalarization to other fields. For vector fields, Garcia-Saenz et al. [51] has identified the presence of ghost and gradient instabilities in the background of compact objects in a broad class of generalized Proca theories [84,85]. Similar concerns were also raised in the context of cosmology in Ref. [86]. Going beyond vector fields, all known formulations of nonminimally coupled spin-2 fields that could spontaneously grow are known to lead to ghost instabilities as well [75]. Likewise, pform fields also have the same constraint structure we discussed in Sec. V, hence they suffer from similar divergent terms [87]. Spontaneous growth of spinor fields as it was introduced in Ref. [88] also contains divergent terms.
The only potential exception to our long list of problematic theories is a second form of spontaneous spinorization theory proposed in Ref. [89], whose equations of motion are not known to feature divergences. It remains to be seen if other well-posed theories exist. If this is the case, understanding what distinguishes these theories at a fundamental level from the problematic ones may lead to a proper no-go theorem for arbitrary generalizations of spontaneous scalarization. The authors also acknowledge networking support by the GWverse COST Action CA16104, "Black holes, gravitational waves and fundamental physics". Some of our calculations were performed with the Mathematica packages xPert [90] and xCoba, parts of the xAct/xTensor suite [91,92].

Appendix A: The Hellings-Nordtvedt theory
In this Appendix we apply the approach of Sec. V to examine the field equations in the Hellings-Nordtvedt [73,74] vector-tensor theory studied in Ref. [45] as a vectorization model.
In this theory, the vector field obeys the field equation, where ω and η are dimensionless coupling constants. Let us first obtain a generalized Lorenz constraint satisfied by A α by taking a covariant derivative of Eq. (A1) and using ∇ ν ∇ µ F νµ = 0, We can expand this equation and replace the Ricci tensor with the Einstein tensor and the Ricci scalar. The resulting constraint equation is, We can now return to Eq. (A1), write F αβ in terms of A α , follow the same steps that lead to Eq. (50), and find: (A4) At last, using Eq. (A3) we obtain, where M αβ = 1 2 ωRg αβ + 1 + η 2 R αβ − ∇ α ∇ β ln |ωR|, (A6) which should be compared against Eq. (52). Note that in Eq. (A5) the last term in the first line is also second order, hence, it contributes to the principal part of the differential equation in addition to the wave operator. Hence, this equations is not in an explicitly hyperbolic form, and we cannot immediately identify M αβ as a squared-mass tensor whose eigenvalues are related to the effective masses of the individual degrees of freedom. However, such identification is possible in the special case η = 0 in which the problematic term vanishes and then: M (η=0) αβ = (ω/2)Rg αβ + R αβ − ∇ α ∇ β ln |ωR|. (A7) We see, by comparing with Eqs. (A7) and (52), that ωR plays the role ofẑ. We then conclude that a ghost arises for the same reasons discussed in Sec. V. For the general case η = 0 it is more convenient to analyse the constraint (A2) which we write as, where Ξ α β = η G α β + (ω + η/2) R δ α β .
Let us focus on the perturbative regime where the background metric is fixed and the Einstein equations hold, i.e., G αβ = 8πT αβ [45]. For a static, spherically symmetric perfect fluid star with energy density ε and pressure p, which is diagonal. The constraint can then be written as where we wrote the summation over the spatial coordinates k explicitly to avoid confusion. This means the diagonal elements have the role of a generalizedẑ in the massless case in Eq. (53). We see that ∂ 0 A 0 has a contribution in the form of The behavior of this term is given by the dependence of the energy density and the pressure on the radial coordinate at the surface of the star. We normally encounter power law dependence in stars due to the TOV equations as we mentioned in relation to Eq. (53). Hence, ∂ 0 A 0 diverges for generic configurations of A µ . We conclude by noticing that the constraint equations of disformally coupled vector-tensor theories of Ref. [49,50] have a similar structure to Eq. (A8), which would lead to similar results in terms of divergences.