A singularity problem for interacting massive vectors

Interacting massive spin-1 fields have been widely used in cosmology and particle physics. We obtain a new condition on the validity of the classical limit of these theories related to the non-trivial constraints that exist for vector field components. A violation of this consistency condition causes a singularity in the time derivative of the auxiliary component and could impact, for example, the field's cosmic history and superradiance around black holes. Such a condition is expected to exist generically in many other non-trivially constrained systems.

Introduction.-The lack of direct evidence for Weakly Interacting Massive Particles (WIMP) has driven people to explore different dark matter candidates, among which light massive spin-1 (Proca) fields, the so-called "dark photons" or vector dark matter, have been drawing more attention in recent years [1][2][3]. In general Proca fields could have a variety of non-gravitational interactions and thus very rich dynamics. For example, the coupling to an axion field may allow for a significant energy transfer from axions to dark photons, and make the latter the dominant component of dark matter in the present-day universe [4,5]. If the Proca field possesses a nonlinear self-interacting potential or a non-minimal coupling to gravity, it may drive the cosmic inflation in the early universe [6,7], and support coherently oscillating localized solitonic field configurations called vector oscillons [8]. The existence of strong self-interactions would also weaken the superradiance bounds on ultralight vectors [9,10]. Moreover, the theory of vector Galileons, whose effective action contains self-interactions with higherorder derivatives, has been constructed by requiring that the equation of motion has 2nd-order time derivatives and yields three healthy propagating degrees of freedom [11]. As an IR modification of gravity, it has been shown that these self-interacting Proca fields can lead to a viable cosmic expansion history and even alleviate the Hubble tension without sabotaging the success of General Relativity on scales of the solar system [12][13][14].
Given these manifold applications, it is necessary to examine the consistency of the interacting Proca theory carefully. One guiding principle that often comes into play is the validity of an effective description of the interaction, which may arise from a low-energy approximation of coupling to other fields or non-minimal coupling to gravity. Another standard lore is that theories with ghosts or energies unbounded from below are usually unstable and problematic [15][16][17] and the initial conditions must be restricted in "islands of stability" if possible [18,19], although there may be some exceptions [20]. Regarding massive vectors, it is pointed out that if they are non-minimally coupled to gravity, the longitudinal mode may exhibit ghost instabilities and one can not discuss the vector field dynamics in a healthy way [21,22]. In practice, one performs as many sanity checks as possible to determine the scope of application for the theory in hand.
In this letter, we will discuss another type of bound that can arise in the classical limit of the theory by demanding the absence of a singularity problem forȦ 0 , due to the fact that the interacting Proca theory is a nontrivially constrained system, where the auxiliary component A 0 can not always be uniquely solved in terms of the canonical fields. A similar type of constraint is also expected for interacting massive spin-2 fields.
In what follows, we will first clarify three consistency conditions and introduce the singularity problem, then a specific example is carried out in detail both analytically and numerically. Implications of our results are discussed finally. We adopt the mostly plus convention for the metric.
The singularity problem.-For definiteness, let us consider a real-valued massive vector field A µ = (A 0 , A) with the Lagrangian 1 where F µν = ∂ µ A ν − ∂ ν A µ and there is no gauge invariance thanks to the potential V , which includes a mass term along with self-interactions. A concrete example is the Abelian-Higgs model, where a quartic selfinteraction is induced by Higgs exchange in the lowenergy limit. By varying the action S = d 4 xL with respect to the field A ν , we find the Euler-Lagrange equation ∂ µ F µν − 2V (A µ A µ )A ν = 0. In vector notation, it becomes where the prime denotes the derivative of the potential in terms of A µ A µ and we have defined the conjugate field 1 The analysis that will be carried out may be straightforwardly generalized for complex fields and for other interactions, see final discussions for more details.
One more useful equation can be obtained by noting that F µν is antisymmetric, so that ∂ µ (V A µ ) = 0. That is In the language of Hamiltonian mechanics, Π 0 = 0 is a primary constraint, equation (2) is a secondary constraint obtained by requiringΠ 0 = δH/δA 0 = 0, and equation (5) is a tertiary constraint obtained by requiring the secondary constraint to be preserved in time [23]. The foundation for these derivations is the stationary action principle, in which we have implicitly assumed that the field A µ is continuous otherwise the infinitesimal variation δA µ is ill-defined. By applying the above formalism, therefore, we require a consistent classical system to satisfy at least three conditions everywhere: (ii) The field A µ (t, x) is continuous; (iii) The second-class constraints, e.g. (2) and (5), are respected.
These conditions are not trivial, and we may gain some insights about them by using equations (3)-(5) and numerically evolving Π, A and A 0 . Given appropriate initial conditions, suppose that V − 2V A 2 0 never becomes 0, then the infinitesimal variations δΠ, δA and δA 0 are always well defined in a infinitesimal time interval δt, and the field A µ will remain smooth and unique all the time. This is indeed the case for free massive fields, where V (A µ A µ ) = m 2 A µ A µ /2. For theories with self-interactions, however, a singularity is encountered if V − 2V A 2 0 becomes 0 at some spacetime point unless −2A 0 V (A ·Ȧ) + ∇ · (V A) also vanishes in an appropriate way to ensure a finiteȦ 0 , which otherwise causes a discontinuity in A 0 and violate at least one of the consistency conditions. Thus maintaining the continuity of A 0 at this point needs an over-constraint and requires fine tuning of initial conditions. It is seen that any plausible interacting Proca theories should ensure that such a problem is avoided in its validity, to wit, the field value should never cross the boundary in field space {|A 0 |, |A|} specified by One may think that the problem identified here can be easily avoided if we use equation (2) instead of (5) to obtain A 0 . As will be shown shortly, this difficulty is actually independent of whether and how we evolve the system numerically.  (12), for repulsive (λ > 0, left) and attractive (λ < 0, right) self-interactions. There are one, two or three different real roots of A0 in equation (8) when ∆ is >, = or < 0.
A concrete model.-In order to understand the significance of this singularity bound, it is illuminating to consider the simplest possibility of a self-interaction where A 0 can be solved in closed form. We are going to show that if we stick with condition (i) and (iii), then condition (ii) will necessarily be violated if the field system hits the boundary (6) during its evolution. The secondary constraint (2) in this case becomes where c 1 = −m 2 /λ − A 2 and c 2 = −∇ · Π/λ. The general solution of this cubic equation can be given by the Cardano's formula, 0 is a real root and the other two are complex conjugate if ∆ > 0. All three are real roots with A The value of ∆ in terms of ∇ · Π and A is shown in figure  1.
A few subtleties need to be clarified when we apply the Cardano's formula (9)- (11). First, we have defined the square root of any number by its principal value. Second, we have defined the cube root of any number by its principal value when c 2 ≤ 0 and its anti-principal value when c 2 > 0. The (anti-)principal cube root returns the real cube root for a real number, and the root with the (smallest) greatest real part for a complex number. These conventions are adopted such that A (1) 0 is always real and all three roots are continuous everywhere except at ∇ · Π = 0, where A 0 can actually remain continuous by switching roots. If a field system crosses the boundary ∆ = 0 (and ∇ · Π = 0) during its evolution, then the real roots A (2,3) 0 will be annihilated or created depending on which region in figure 1 the system is in before the crossing. It is easy to see that the discontinuity of A 0 is an inevitable 0 , and if the system hits ∆ = 0 from the white region where ∆ < 0. Now we will show that A 0 can not remain continuous if the system hits the boundary specified by equation (6). To do this, we can judiciously rewrite the discriminant ∆ in terms of |A 0 | and |A| by using the secondary constraint (2). The value of ∆ in terms of |A 0 | and |A| is shown in figure 2. At ∆ = 0, the three roots (9)-(11) become |A and A 0 = 2A 0,crit and A 0 = A 0,crit are visualized as the gray dashed and solid black curves respectively in figure 2. Note that only A 0 = A 0,crit (solid black line) corresponds to the boundary V − 2V A 2 0 = 0, and the adjacent regions separated by this line both have ∆ < 0, which justifies the foregoing claim.
On the other hand, it is always safe to cross the gray dashed line, since in this case A 0 follows the root A is real and continuous. But there is no guarantee that the evolution will be healthy if the root A (1) 0 is chosen for A 0 initially, because A 0 switches roots at ∇ · Π = 0. In order to avoid the singularity problem and also allowing field values to be small, we conclude that the field evolution should be restricted in the non-meshed region in figure 2.
A minimal model of (7) is carried out numerically in 1 + 1-dimensional spacetime to support the above analysis. We present field-space trajectories for repulsive self-interactions in figure 3. The case of attractive selfinteractions is similar, and thus only shown in the appendix A, where numerical details are also provided. 3 We are only interested in the case where A 0,crit is real.  2. The value of the discriminant ∆ in terms of A0 and A for repulsive (λ > 0, left) and attractive (λ < 0, right) self-interactions. The colored and white regions correspond to ∆ > 0 and ∆ < 0 as in figure 1, and the gray dashed and black solid curves represent A0 = 2A0,crit and A0 = A0,crit at which ∆ = 0. A consistent classical system should never cross the black solid curve, which is exactly the one specified by equation (6) (see the texts for proof). Allowing field values to be small, the system during the evolution should never enter into the meshed region.  (6). The colored and white regions, and the gray dashed and black solid curves have the same meaning as in figure 1 and 2. The blue and red trajectories represent the time evolution of fields at two adjacent spatial locations starting from the solid point. As shown by the blue trajectory, when the system meets the black solid boundary, the value of A0 can no longer remain continuous and suddenly jumps to the gray dashed line, which violates the consistency conditions. Numerical details and an example for attractive self-interactions are provided in the appendix A.
Up to this point, it looks like that the "discontinuity problem" is a more appropriate name inasmuch as the temporal component A 0 can not be continuous when the field system hits the boundary specified by (6). In fact, the discontinuity is just an artificial phenomenon because in principle we can stick with the condition (i) and (ii) instead, and then we will reach a conclusion that the second-class constraints can not be obeyed. The real problem is that at least one of the consistency conditions will be violated if the boundary (6) is hit, which is closely related to a singularity inȦ 0 .
Conclusions and discussions.-We have demonstrated that there exists a generic constraint on field values for interacting Proca theory. Such a constraint can be obtained by observing that a classical massive spin-1 field A µ (t, x) should be both real-valued and continuous, and the second-class constraints of the system be preserved everywhere. Ensuring these conditions is crucial to get realistic results in some contexts, where nongravitationally interacting classical fields may play a pivotal role, such as density perturbations in dark photon production [5], vector oscillons 4 [8] and black hole superradiance of vectors [9,10]. Taking non-derivative selfinteractions as a concrete example, we have shown that the field system during its evolution should never meet the singularity bound specified by equation (6), otherwise the effective description becomes inconsistent.
Loosely speaking, trajectories in phase space (if we could ever visualize them for PDEs) would intersect at the singularity bound (6), indicating that A 0 can no longer be solved uniquely. This situation is usually avoided in physics equations because of the Picard-Lindelof theorem (also called the existence and uniqueness theorem), which states that the existence and uniqueness of solutions are guaranteed if the derivative of the variable is continuously differentiable. To fully appreciate this point, in appendix B we provide a toy model with second-class constraints where the phase portrait can be shown explicitly. We also note in there the existence of a basin of attraction towards the singularity bound. If the same goes for interacting massive vectors, the allowed field space is further restricted.
As another example, gauge-invariant interactions (such as those only involving F µν ) do not yield the singularity problem, because the gauge-invariant part of the action must satisfy ∂ µ (δS GI /δA µ ) = 0 while we have ∂ µ (δS/δA µ ) = 0 if the equation of motion is satisfied. 5 Thus the tertiary constraint like (5) can be obtained solely from gauge symmetry breaking terms. The analysis in this letter may be straightforwardly generalized for more general interactions for Proca fields, and even for interacting massive spin-2 fields.
We note that the singularity bound equally applies for complex fields. To see this, we may separate the real and imaginary part of a complex field A µ = R µ + iI µ , then the theory of A µ becomes a theory of two interacting real fields R µ and I µ , where each field is constrained by demanding the absence of the singularity. 4 Our constraint for the self-interaction considered in [8] implies that vector oscillon fields there are bounded by equation (6) and are unlikely to have very large amplitudes, in constrast to their scalar counterpart with a φ 6 self-interaction [24,25]. 5 By gauge invariance here we mean that the action is invariant under the gauge transformation Aµ → Aµ+∂µΛ. More generally, a combined gauge transformation with another scalar field can be invented for massive vectors by using the Stueckelberg trick [26]. In this case one must generalize ∂µ(δS GI /δAµ) = 0 to include the scalar field.
The existence of the singularity problem in a theory indicates that the theory can not be the complete story. Learning from perturbative unitarity [27,28], the standard solution would be to introduce new particles or to look for a UV completion above the scale where the singularity bound is met. For example, we can introduce a Higgs boson to rescue the quartic theory (7). We leave such considerations for future work.
Note: In the final stage of preparation of this work, we became aware of another related forthcoming work by Clough, Helfer, Witek and Berti, which considers ghost instabilities in massive self-interacting vector fields. We will compare it with our work in the future.

ACKNOWLEDGMENTS
We would like to thank Mustafa Amin, Ray Hagimoto, Soichiro Hashiba, Mudit Jain, Siyang Ling and Andrew Long for fruitful discussions. We would especially like to thank Mustafa Amin and Andrew Long for a careful reading of the manuscript and suggestions for its improvement. HYZ is partly supported by DOE-0000250746.

Appendix A: Numerical details
As proved in the main body of the letter, if we demand that the field is real-valued and satisfies all constraint equations, the value of A 0 can not remain continuous when the boundary of the singularity is met. This conclusion can be verified with numerical simulations.
For simplicity, we may consider a Proca field in 1+1dimensional spacetime A minimal model after discretization consists of two spatial sites x 1 and x 2 with periodic boundary conditions. For notation convenience, we define fictitious points x 3 ≡ x 1 and x 0 ≡ x 2 . Then the task of solving PDEs numerically can be simplified into a problem involving a set of ODEs by defining fields at each site where i = 0, 1, 2, 3 (while those fields with index 0 and 3 are fictitious), and the conjugate field Z i (t) ≡ F 01 (t, x i ) becomes where we have used forward difference to approximate spatial derivatives and dx = x 2 − x 1 . 6 As a result, the Lagrangian (A1) becomes where we have used the trapezoidal rule to approximate the integral, and thus the governing equations corresponding to (2)-(5) are Since we may meet a numerical singularity in equation (A6) if the denominator on the RHS becomes 0, equations (A2)-(A5) are used to evolve the system. The equation (A4) may yield one, two or three real roots for X i . Initially we may choose arbitrary one, but during the evolution we always take the closest value to X i (t − dt) as the solution of X i (t) to ensure its continuity to the most extent. The value of dx is not essential to this problem, and we simply set it unity. An example for repulsive self-interactions is given in figure 3 in the main body of the letter. Initially we set Y 1 = 0.4, Y 2 = 0.2, Z 1 = −0.2 and Z 2 = −0.4, and there are three real roots for each X i since the system stays in the region where ∆ < 0. The explicit solution shown in figure 3 corresponds to the root A (2) 0 for X 1 and A (1) 0 for X 2 , i.e. X 1 = −0.977 and X 2 = −1.105. The blue and 6 The numerical system is usually unstable if spatial derivatives are discretized by forward/backward difference method, but the numerical stability is actually irrelevant here since we do not solve the PDEs directly. red curve shows the evolution at x 1 and x 2 respectively. This is consistent with our analytical analysis: only A 0 for X 1 and X 2 respectively, i.e. X 1 = 1.076 and X 2 = 1.028.

Appendix B: A toy model
In this appendix, we study the phase space evolution of a toy model that is constrained by second-class constraints just like a Proca system. Consider a Lagrangian where x and y are real-valued variables, and the potential is This model is invented such that x and y mimic A 0 and A in Proca theory, and since it is equipped with ODEs we can draw a complete phase portrait for it. By carrying out the same analysis, it is easy to find the singularity bound where the prime denotes derivative of the potential with respect to −x 2 + y 2 . This corresponds to the bound (6) for self-interacting vectors. The canonical conjugate variable for this model is z ≡ dL/dẏ =ẏ − x for y and 0 for x. Thanks to the constraint, z can be solved in terms of x and y and the phase space is actually 2-dimensional, which is shown in figure 5.
Similar to the quartic model in the main text, the discriminant of the model is ∆ ≡ 1 4 z 4b 2 − 1 27 a 2b + y 2 3 , where there are one, two or three different real roots for x if ∆ is >, =, < 0. In figure 5, both the black solid and gray dashed lines correspond to ∆ = 0, but only the former is the bound specified by equation (B3). From the trajectories we see that there is actually a basin of attraction for the singularity bound, where different trajectories intersect. It is easy to see that finite time is needed to hit the singularity bound if a point in phase space is attracted towards it. The critical exponent behaves like |x − x 0 | ∼ |t − t 0 | 1/2 as the point approaches the bound. If this is also true for interacting massive vectors, the allowed field space is further restricted.