Proca Stars: gravitating Bose-Einstein condensates of massive spin 1 particles

We establish that massive complex Abelian vector fields (mass $\mu$) can form gravitating solitons, when minimally coupled to Einstein's gravity. Such Proca stars (PSs) have a stationary, everywhere regular and asymptotically flat geometry. The Proca field, however, possesses a harmonic time dependence (frequency $w$), realizing Wheeler's concept of geons for an Abelian spin 1 field. We obtain PSs with both a spherically symmetric (static) and an axially symmetric (stationary) line element. The latter form a countable number of families labelled by an integer $m\in \mathbb{Z}^+$. PSs, like (scalar) boson stars, carry a conserved Noether charge, and are akin to the latter in many ways. In particular, both types of stars exist for a limited range of frequencies and there is a maximal ADM mass, $M_{max}$, attained for an intermediate frequency. For spherically symmetric PSs (rotating PSs with $m=1,2,3$), $M_{max}\simeq 1.058 M_{Pl}^2/\mu$ ($M_{max}\simeq 1.568,\, 2.337, \, 3.247 \, M_{Pl}^2/\mu$), slightly larger values than those for (mini-)boson stars. We establish perturbative stability for a subset of solutions in the spherical case and anticipate a similar conclusion for fundamental modes in the rotating case. The discovery of PSs opens many avenues of research, reconsidering five decades of work on (scalar) boson stars, in particular as possible dark matter candidates.

Introduction. According to the latest cosmological data [1], ∼ 26% of the Universe's energy content is dark matter (DM). The fundamental nature of DM, however, is unknown. A widespread viewpoint is that DM consists of weakly interactive massive particles (WIMPs); popular candidates are the lightest supersymmetric particles, with masses GeV, such as the neutralino [2]. Despite its success in modeling structure formation, some shortcomings of the WIMPs DM model arise in small scales, such as the "missing satellite" and the "cuspy core" problems. Different proposals, some of which have been claimed to solve these problems (see [3,4] for reviews), introduce light or ultra-light bosonic particles/fields, with masses ≪ eV, which may form, gravitationally, macroscopic Bose-Einstein condensates. Even though such proposals have essentially a phenomenological character, ultra-light particles (axions) are a natural ingredient of the Peccei-Quinn mechanism to solve the strong CP problem [5] and may be motivated at a fundamental level by string theoretical constructions -the axiverse [6].
Gravitationally-bound bosonic structures are thus relevant in the context of DM searches. The study of (scalar) Bose-Einstein condensates as DM candidates is often performed in a Newtonian limit. In the fully relativistic regime, such models yield gravitating solitons: (scalar) boson stars (SBSs). These objects, initially proposed as a (scalar) realisation of Wheeler's geon idea [7,8], have found a variety of applications, from black hole mimickers in astrophysics, to particle models in TeV gravity scenarios (see [9,10] for reviews).
In light of recent proposals advocating massive spin 1 particles as a DM ingredient [11][12][13][14], the study of the cor-responding self-gravitating structures and their dynamics is of special importance [15]. In this letter we show that much like massive spin 0 particles, massive spin 1 particles can cluster as everywhere smooth, asymptotically flat lumps of energy under their own weight, producing gravitating solitons we dub Proca stars (PSs). Moreover, we observe a close parallelism between the physical properties of PSs and SBSs. Einstein-(complex)-Proca theory. We consider one complex Proca field, with mass µ (or equivalently, two real Proca fields with the same mass). It is described by the potential 1-form A, and field strength F = dA. We denote the corresponding complex conjugates by an overbar,Ā andF . The minimal Einstein-(complex)-Proca model is described by the action: The Einstein and Proca field equations are, respectively, G αβ = 8πGT αβ , ∇ α F αβ = µ 2 A β , where the energymomentum tensor reads: The Proca equations imply the Lorentz condition (which is not a gauge choice, but a dynamical requirement): The global U (1) invariance of the action, under the transformation A µ → e iα A µ , with α constant, implies the existence of a conserved 4-current, From the Proca equation, ∇ α j α = 0. Thus a conserved Noether charge exists, Q, obtained integrating the temporal component of the 4current on a space-like slice Σ: Spherically symmetric PSs. We first consider spherically symmetric solutions, taking the following line element form where N (r) ≡ 1−2m(r)/r and the Proca potential ansatz σ(r), m(r), f (r), g(r) are all real functions of the radial coordinate only and w is a real frequency parameter. As for SBSs the harmonic time dependence of A is crucial and the complex nature of A makes it compatible with the time-independent line element. But unlike SBSs there are now two independent radial functions for the 'matter' field, making the analysis more involved (even more so in the rotating case). The Proca equations yield where ′ denotes radial derivative. These two equations imply the Lorentz condition constraint (which determines f (r) in terms of the remaining functions). The essential Einstein equations read For the ansatz (1)-(2), the Noether charge reads Q = 4πµ 2 w ∞ 0 dr r 2 g(r) 2 σ(r)N (r), and the energy density measured by a static observer is Finally, PSs satisfy the virial relation: used to test the numerical accuracy of the results. This relation can also be used to rule out non-gravitating solutions, i.e 'Proca-balls' without backreaction [16]. Asymptotic behaviours. An analysis of the field equations both near the origin and at spatial infinity reveals smooth behaviours. Close to r = 0 we find where f 0 , σ 0 are constants, the values of f (r) and σ(r) at the origin, respectively. As r → ∞, we find where M is the ADM mass and c 0 is a constant; observe that w < µ, which is a bound state condition.
Numerical results. The solutions that smoothly interpolate between the two above asymptotic behaviours are found numerically. In numerics -and in the followingwe set µ = 1, 4πG = 1, by using a scaled radial coordinate r → rµ (together with m → mµ, w → w/µ) and The equations are solved by using a standard Runge-Kutta ODE solver and implementing a shooting method in terms of the parameter f (0).
In Fig. 1 we plot the ADM mass, M , and the Noether charge, Q, as a function of the scalar field frequency, w. As w → 1, the mass and Noether charge of the solutions vanish, but Q/M → 1. Approaching this limit, PSslike SBSs -become spatially diluted with an effective size much larger than their Schwarzschild radii [17], and trivialize at w = 1. Reducing w from this maximal value, PSs become more compact; both M and Q follow a spiral, towards a central critical configuration located around w ≃ 0.891. The maximum of both M and Q occurs for w max ≃ 0.875, with M max ≃ 1.058 < Q max ≃ 1.088. We observe that M < Q from w = 1 down to almost the minimal allowed frequency. In the lower part of the spiral M > Q, and the binding energy 1 − M/Q is negative, corresponding to a region of "energy excess", where the solutions must be unstable against perturbations. All these features follow closely those of spherical SBSs. We remark that for the family of solutions exhibited in Fig. 1, f (r) has one node [18] and g(r) is nodeless. These are the fundamental solutions. There are also excited solutions, with more nodes for these functions, but we did not attempt to study them in any systematic way. Studies of SBSs suggest that excited solutions are unstable [19]. In Fig. 2 we show the behaviour of physical and profile functions for two PS solutions. In particular observe the energy concentration near r = 0.

Stability.
A positive binding energy -as that observed along the upper branch of Fig. 1 -does not guarantee linear stability. For SBSs, a perturbative stability analysis shows that the solutions are only stable from the maximal frequency w = 1 down until the frequency at which the maximal mass is attained, w max . Beyond this point an unstable mode develops [20,21]. An analogous conclusion holds for PSs, as we now show by considering their linear radial perturbations. We assume that all perturbations have a harmonic time dependence of the form e −iΩt , with Ω being the characteristic vibrational frequencies of the PS. Gauge freedom allows to write the perturbed metric as: the vector field is perturbed as where h 0 , h 1 , f 1 , f 2 , g 1 and g 2 are radial perturbations around the background solution, and ǫ is a small, bookkeeping parameter. The perturbed Einstein-Proca system can be reduced to a pair of second order ODE's. Imposing regularity of the perturbations at the origin and at infinity, the resulting system is a two dimensional eigenvalue problem for Ω and one other constant which we have chosen to be the value of h 0 at the origin. A numerical solution is then obtained by a two dimensional shooting, with the result shown in Fig. 3, where we plot Ω 2 versus the PS's ADM mass around the maximal mass. For Ω 2 > 0 the mode is purely real and corresponds to stable normal modes. This occurs for w > w max . For Ω 2 < 0 the mode is a purely positive imaginary number, indicating that these configurations are unstable (cf. Eqs. (3) and (4)). Thus, in complete agreement with the SBSs case [20][21][22][23], we can see that the star's maximum mass corresponds to a transition point between stable and unstable configurations. For unstable configurations with mass 0.1% below the threshold, the instability timescale τ is already smaller than τ 100M . Thus, we expect the unstable branch of PS to share similar properties to those observed in the scalar case: configurations which reach this branch will either quickly collapse to black holes or migrate back to the stable branch via mass ejection (the "gravitational cooling" mechanism) [22,[24][25][26].
The Proca field ansatz is given in terms of another four functions (H i , V ) which depend also on r, θ with m ∈ Z + . The Einstein-Proca equations are solved with the following boundary conditions, which we have found to be compatible with an approximate construction of the solutions on the boundary of the domain of integration: ,π = 0 (note that for m = 1 the conditions for H 2 , H 3 are different, with ∂ θ H 2 θ=0,π = ∂ θ H 3 θ=0,π = 0). All rotating PSs we have constructed so far are symmetric w.r.t. a reflection along the equatorial plane. Odd-parity composite configurations, however, are also likely to exist. As usual, the ADM mass M and angular momentum J are read off from the asymptotic expansion, g tt = −1 + 2GM /r + . . . , g ϕt = −2GJ sin 2 θ/r + . . . . In analogy with the SBSs case [27], one can show that the angular momentum is a multiple of the Noether charge, J = mQ, but in contrast to the SBSs case, the angular momentum density T t ϕ and the Noether charge density j t are not proportional any longer.
The Einstein-Proca equations for this rotating case are quite involved and shall not be presented here. They are solved numerically, subject to the above boundary conditions, by using the elliptic PDE solver FIDISOL [28] based on the Newton-Raphson procedure, employing a compactified radial coordinate x = r/(1 + r) [29]. In Fig. 4 we exhibit the mass vs. angular momentum, together with the mass and the Noether charge vs. the Proca field frequency, for rotating PSs with m = 1, 2, 3 (inset). Both panels exhibit the familiar shape from the study of SBSs. In particular M and J are always positively correlated [30]. The physical quantities for the largest mass PSs are exhibited in the following table, for m = 0 (spherical) and m = 1, 2, 3 (rotating), where they are compared with the results in the literature for the corresponding largest mass SBSs (last column); one observes that M max is always smaller for the latter.  Other physical properties of the rotating PSs mimic those of SBS: the rotating PSs solutions spatially delocalize as w → 1 (similarly to the static case) and trivialize in that limit; the energy density has an essentially toroidal distribution; also, for m = 1, neither the energy density nor the Noether charge density vanish on the symmetry axis [32].
We have not attempted to study in detail the stability of rotating PSs. For rotating SBSs, catastrophe theory arguments [33] support a similar conclusion to that obtained in the static case: the solutions are stable from the maximal frequency down to the point where the maximal mass is attained. We expect a similar conclusion to hold for rotating PSs. In the rotating case there can also occur ergoregion instabilities [34]; for SBSs ergoregions only appear in the unstable branch of solutions and an analogous situation is true for PSs [35]. Discussion and Outlook. The discovery of PSs suggests various novel directions of research, besides the detailed analysis of the solutions presented herein; let us mention two. 1) In the scalar case, a single (rather than complex) real field allows the existence of very long lived quasi-solitonsoscillatons; the same is true for the Proca case [26]. 2) Can one add a black hole inside a PS, like for other gravitating solitons? While for spherically symmetric configurations one can show the absence of black hole solutions, this is possible inside a spinning PSs, again in complete analogy with the SBSs case [36][37][38].
A final remark relating PSs with DM phenomenology. Since ordinary matter, which constitutes only ∼ 5% of the Universe's energy content, is comprised of over ten elementary particles, it is reasonable to admit that DM is composed of different kinds of fundamental entities, but which, when gravitationally clustered into macroscopic lumps, display some universality. If Bose-Einstein condensates of ultra-light scalar fields are viable DM models, the existence of PSs akin to SBSs suggests the former may be another DM component. In this context it would be interesting to study the Newtonian limit of PSs.