Perturbations of ultralight vector field dark matter

We study the dynamics of cosmological perturbations in models of dark matter based on ultralight coherent vector fields. Very much as for scalar field dark matter, we find two different regimes in the evolution: for modes with $k^2\ll {\cal H}ma$, we have a particle-like behaviour indistinguishable from cold dark matter, whereas for modes with $k^2\gg {\cal H}ma$, we get a wave-like behaviour in which the sound speed is non-vanishing and of order $c_s^2\simeq k^2/m^2a^2$. This implies that, also in these models, structure formation could be suppressed on small scales. However, unlike the scalar case, the fact that the background evolution contains a non-vanishing homogeneous vector field implies that, in general, the evolution of the three kinds of perturbations (scalar, vector and tensor) can no longer be decoupled at the linear level. More specifically, in the particle regime, the three types of perturbations are actually decoupled, whereas in the wave regime, the three vector field perturbations generate one scalar-tensor and two vector-tensor perturbations in the metric. Also in the wave regime, we find that a non-vanishing anisotropic stress is present in the perturbed energy-momentum tensor giving rise to a gravitational slip of order $(\Phi-\Psi)/\Phi\sim c_s^2$. Moreover in this regime the amplitude of the tensor to scalar ratio of the scalar-tensor modes is also $h/\Phi\sim c_s^2$. This implies that small-scale density perturbations are necessarily associated to the presence of gravity waves in this model. We compare their spectrum with the sensitivity of present and future gravity waves detectors.

though the field evolution is generically anisotropic, the average energy-momentum tensor is not. In particular, for massive fields it is straightforward to show that the average energy density scales as a −3 . This opens the possibility of extending the wave DM scenario to higher spin fields.
In this work, we will consider the case of massive abelian vector fields. The interest in homogeneous vector fields as cosmological fluids has been growing in the last years, see [57][58][59][60][61][62] for dark energy examples and [63,64] for inflation models based on vector fields. The possibility that a condensate of very light vector particles could play the role of DM was explored in [65] and a wide phenomenological study of this model was made in [66]. Such a condensate could be produced during inflation and its small mass could be generated by the Stuckelberg mechanism. A small kinetic mixing with the photon could make this dark photon detectable [67][68][69].
Here, in particular, we will concentrate in the dynamics of cosmological perturbations in such vector DM models. The main difficulty compared to the scalar case is the presence of a non-vanishing vector field already at the background level. This implies that the usual decoupling of the evolution of scalar, vector and tensor perturbations at the linear level no longer holds in this case. However, this fact provides a potential way of discriminating vector models from scalar ones. In particular, we will show that although in the particle regime with k 2 Hma the model is indistinguishable from CDM, in the wave regime k 2 Hma the scalar and vector modes are coupled to the tensors. This implies that unlike scalar field models, density perturbations generate a specific gravity wave spectrum together with a non-vanishing anisotropic stress.
The work is organized as follows. Firstly, in Section II, we will review the time averaging procedure in cosmology. Then, in Section III, we will consider the anisotropy problem for homogeneous vector fields. In Section IV we obtain the basic equations for perturbations of massive vectors and in Section V we write them for scalar and vector modes. Sections VI and VII are devoted to the results for scalar and vector perturbations where the adiabatic solutions of the perturbations equations are obtained in the different regimes. In Section VIII we concentrate on the generation of gravity waves and Section IX in the possibility of detection. Finally Section X includes the main conclusions of the work.

II. TIME AVERAGING IN COSMOLOGY
Many of the results we will obtain in this work are based on the assumption that in the presence of rapidly oscillating fields, it is possible to time average the energy-momentum tensor so that the resulting solutions of Einstein equation are a good approximation to the exact ones. In order to determine when this is the case, we will consider a simple example, which will help us to understand the key aspects of this procedure.
Let us consider a homogeneous scalar field oscillating in a power-law potential with n and even integer. In a flat FLRW metric in proper time, the equation of motion can be written asφ where the dot represents the t derivative. Making the change φ =φ a − 6 n+2 and dt = a where is the derivative with respect to the new time variableη. Let us assume that the frequency of the oscillations ω is large compared to the rate of expansion of the universe i.e. ω H withH = a /a. Thus we can define the small parameter ≡H/ω. Accordingly, the terms proportional toφ in (4) will be suppressed by O( 2 ) compared to the other ones and we can writeφ Thus, the solution can be written in terms of the field φ as where F = a − 6 n+2 is a slowly evolving fuction ofη with F /F ∼H and P is a periodic fast oscillating function with period 2π/ω, i.e. P /P ∼ ω.
Let us now try to obtain the scale factor a(η) from Einstein equations. The system formed by the Friedmann and conservation equations read from which we can obtain,Ḣ so that integrating twice in time we get: The φ terms in the integrand are dominated by the derivatives of the rapidly oscillating function so that we can approximate Since P 2 (η) is also a periodic function we can Fourier expand it as: Let us now perform the first time integration Integrating by parts the m > 0 terms we get, Notice thatF 2 (η 1 ) is proportional to the first derivative of I 0 which in general is expected to be, Thus we see that compared to the I 0 term, the amplitude of the oscillating I m>0 contributions are generically suppressed by: Moreover, the second integration in (10), reduces in another O( ) factor the oscillatory contributions. Notice also that the periodic factor of the O( ) correction term in (11) can be expressed as a total time derivative, P (η)P (η) = ∂ηP 2 (η), which does not contribute to the zero mode of the Fourier expansion, c 0 . Thus, in general, we can expand the scale factor as where a m=0 (η) is the contribution from the c 0 term whereas a m>0 (η) are the oscillatory contributions. We can conclude that, up to O( 2 ), it is a good approximation to neglect the oscillatory terms c m>0 in the source of Einstein equations provided the solution involves two time integrations. Thus we will denote by the operation of extracting the m = 0 mode of the Fourier expansion, i.e.: a(η) = a m=0 (η). (18) Notice that this operation is equivalent to time averaging up to O( ) terms. Indeed where in the last step we have considered that both uncertainties are of the same order. Consequently, in general if we consider the average Einstein equations the corresponding solutions for g µν would differ from the exact ones in O( 2 ) terms.

III. MASSIVE VECTOR COSMOLOGY
Let us consider a massive abelian vector field in an expanding universe [56]. The corresponding action reads with The equations of motion are given by: We will first consider the dynamics of the homogeneous background fields. For simplicity we will work with linearly polarized fields Assuming that the energy-momentum tensor is dominated by the vector field, the background geometry can be represented through a Bianchi I metric, The µ = 0 component of the equation of motion reads so that the temporal component identically vanish, whereas the µ = i equations implÿ where dot represents derivative respect to the conformal time η. On the other hand, from the exact Einstein equations Let us now assume that the field A z is oscillating rapidly around the minimum of the potential, i.e. we will consider that ma H, with H =ȧ/a the comoving Hubble parameter and ma ḃ , then the equation of motion (24) can be solved in the WKB approximation as, with = {H/(ma),ḃ/(ma)}. Introducing this solution in the system and averaging (extracting the zero mode as discussed in the previous section), the average Einstein equations readä a +ȧ 2 a 2 = 2πG The second equation shows that there is no source for anisotropy in the average equations, so that Thus, if the initial conditions are isotropic then b(η) = 0 at all times and the third equation readṡ with solution to leading order in Thus, as shown in [56] the average geometry generated by a rapidly oscillating massive abelian vector field is isotropic and evolves as in a matter dominated universe. If the vector field is responsible for all the DM contribution, its amplitude will be given by with Ω c the CDM density parameter and ρ c = 3H 2 0 /(8πG) the critical density.

IV. PERTURBATIONS OF MASSIVE VECTORS.
In the previous section we have shown that despite the anisotropic evolution of the background vector field, the average geometry can be described by an isotropic FLRW metric. Thus, we will consider the most general form of the perturbations around the Robertson-Walker geometry where Q is a solenoidal vector field and h ij a symmetric traceless transverse tensor. From (23), Fourier transforming the spatial dependence, the equation of motion for δA i results, and δA 0 satisfies the constraint The first order perturbations of the energy-momentum tensor can be written in Fourier space as and the corresponding perturbations of the average Einstein equations (33) read Thus, we can define the average energy density and pressure as: Unlike the scalar field case, the perturbations of a vector field can source the three kinds of perturbations. This means that the standard separation in the evolution can be more involved in this case. In this work we will proceed as follows: we will first consider the dynamics of scalar and vector modes neglecting the contributions from gravity waves. We will then analyze the generation of gravity waves and will find that they are generically suppressed compared to the scalar and vector modes, thus proving that our initial assumption was correct.

V. SCALAR AND VECTOR PERTURBATIONS: BASIC FORMULAE AND PRELIMINARIES
From equations (45)(46)(47) setting h ij = 0 we can obtain the following set of equations, which together with the equation of motion (42) will be the starting point of our analysis: withk the unitary vector in the wavenumber direction.
• From the longitudinal part of G 0 i we geṫ • From the solenoidal part of G 0 i we can write Before analysing the modes in the different regimes, we would like to make the following preliminary considerations: 1. For simplicity we will consider a linearly polarized background field, which can be written without loss of generality as with A B (η) a slowly varying amplitude.
We will work in the matter dominated era, assuming that all the DM is generated by the vector field and ignoring for simplicity the small baryon contribution. Thus, from the Friedmann equation, with ρ the average energy density: we get 2. As we are dealing with vector equations, it is very helpful to adopt the orthonormal basis {û a ,û pk ,û p } ≡ û A , (k ×û A )/ sin θ, k − cos θû A / sin θ , for modes withk ∦û A ; whereû A is the normalized vector in A direction and cos θ ≡k ·û A . On the other hand, in the degenerate case withk û A , we can use a new orthonormal basis {û a ,û pk1 ,û pk2 }, where {û pk1 ,û pk2 } span the orthogonal plane toû a . It can be seen that the perturbations in those directions are purely vector with the same dynamics as theû pk , whereas theû a components generate purely scalar modes also with the same behaviour as in the non-parallel case. Finally, no tensor modes are sourced fork û A .
3. Theû a component of (53) gives us an algebraic equation for Q a . Once it is solved, it is straightforward to write the other two components of Q as a function of the scalar perturbations of the metric and the vector field. After that, combining equations (50), (51) and (52), we reach an algebraic system from which we obtain Ψ and Φ depending on A and δA.
4. We have three independent (comoving) scales in the problem, namely, ma, H and k. The main assumption of this work is that ma {H, k}, so that we can define two small parameters = H/(ma) and k/(ma). Depending on the relation between these two ratios, the evolution of the perturbations will behave differently. Thus, as we will show, the case k/(ma) ∼ , i.e. k ∼ H will lead to the standard CDM behaviour, whereas the k/(ma) ∼ 1/2 case, i.e. k ∼ (Hma) 1/2 will correspond to the wave DM behaviour as commented above. We will perform an expansion in and obtain only the leading order term. The sub-leading correction will in general receive contributions from the oscillating terms (m > 0) mentioned in Section II and are beyond the scope of this work. 5. Provided k ma, as mentioned before, we will take for the perturbations an adiabatic ansatz similar to that of the background: with δ A (s,c) slowly evolving amplitudes. On the other hand, in the regime with k ∼ ma, the perturbed field oscillates with a different frequency and as a result in the averaging procedure all the perturbed quantities vanish, consequently a cut off in perturbations is expected in the high wavenumber region.
6. In the very-low wavenumber regime with k/(ma) ∼ 3/2 i.e. k ∼ (H 3 /(ma)) 1/2 the leading order equations get contributions from the oscillating terms that cannot be neglected. Thus, our perturbative approach does not allow to explore this region. Fortunately, for the masses usually considered [66] this range is out of the cosmologically observable band.
In Fig.1 we show the evolution of the different comoving scales involved in the adiabatic expansion, namely, (H 3 /(ma)) 1/2 , H, √ Hma and ma as a function of a from matter-radiation equality to the present time. We see that modes in the wave regime can cross into the particle regime, but this is not possible in the opposite way.
FIG. 1: Evolution of comoving scales from matter-radiation equality for m = 10 −22 eV. The blue line sets the limit between the particle and wave regimes. The yellow region corresponds to super-Hubble modes. The green one corresponds to sub-Hubble modes. The blue area is the wave regime and the pink one is the cutoff region. The orange region on the left corresponds to the region where the perturbative approach breaks down.

VI. SCALAR PERTURBATIONS: RESULTS
As mentioned above, we will concentrate in two different regimes, namely, k/(ma) ∼ and k/(ma) ∼ 1/2 . We will present the results of the perturbative analysis to the leading order in the adiabatic expansion, i.e. up to relative corrections of order . For the components orthogonal to the background vector i = p, pk the solution is straightforward, with C i,(s,c) constants. These components do not contribute to the scalar perturbations Φ and Ψ nor to the density and pressure perturbations. The equations for the i = a component read: Solving we obtain, The behaviour is the same as that of standard CDM. Notice that the non-decaying mode of the scalar perturbation Φ is constant independently of the mode and the gravitational slip vanishes since Φ = Ψ. Moreover, the perturbed energy density is controlled by k 2 η 2 , making the density contrast δ constant for super-Hubble modes and growing as δ ∼ a for sub-Hubble modes as expected. This case corresponds to modes whose wavelength is comparable to the de Broglie wavelength of a comoving DM particle. In this regime the wave properties of DM could have important effects.
As in the previous case, we study the evolution of the different components. For the i = p, pk components, we get to the leading order and with C i1 ( k) and C i2 ( k) constants. Again these components do not contribute to the scalar perturbations of the metric. In order to solve for theû a component, we will use equation (52) and average A· (42) obtaining δA a,c + 2Hδ A a,c + k 4 4m 2 a 2 − H 2 δA a,c = 0 .
By solving (72) we get, The expressions of the metric scalar perturbations are trivially deduced from (69), (70) and the leading order solution (73). The perturbed energy density and pressure can be written as δp(η, k) = − 3 8πG whereas for the scalar metric perturbations we get: The evolution of the scalar perturbation potential is shown in Fig. 2. We see that for the i = p component it is not possible to define a sound speed. For the i = a component the effective sound speed takes a very simple form, however, this expression can become negative. As a matter of fact, since as we will show below, there is a non-vanishing gravitational slip, this is not going to be the characteristic propagation velocity of scalar perturbations. Even though to the leading order we get Φ = Ψ, it is possible to derive the sub-leading contribution from (50) Thus for the i = a components we can write for the gravitational slip: which, as commented before, is O( ) in the adiabatic expansion. As expected, the previous expressions smoothly tend to the standard CDM behaviour discussed in the previous section for k 2 /(maH) 1 (see Fig. 2), indeed In the opposite limit k 2 /(maH) 1 expressions (73)- (75) imply that the perturbations δρ(η, k) 3 8πG are all oscillating and decaying (see Fig. 2). This is completely analogous to the scalar field DM case, and this is the reason why on scales with k 2 maH we expect a suppression in the matter power spectrum as compared to the standard CDM. See Fig. 3 for the modification on the linear transfer function Φ k (a 0 )/Φ k (a eq ) induced on small scales. Notice however that the possibility of generating a gravitational slip is absent in the scalar field case.

VII. VECTOR PERTURBATIONS: RESULTS
From equation (53) , we get to the leading order, the following results in the two regimes in which we are interested We see that to the leading order, only the i = p, pk components contribute to the vector modes. Such modes did not contribute to the scalar perturbations in this regime, and accordingly they are purely vector like. Thus we get Q(η, k) = 8πG 3 3i √ 2H ka 3/2 sin(θ)C p,s ( k)û a − cos(θ)C p,s ( k)û p + cos(θ)C pk,s ( k)û pk . (87) All the components decay as a −2 in the matter dominated era. This is the same behaviour expected for vector modes in standard CDM. Notice also that only the sine components of the vector perturbations δ A s actually contribute to the vector modes. This can be understood since being the background A a cosine function, the first term in (86), which is the only one contributing to the leading order, contains a˙ A factor which in the average procedure is only non-vanishing for sine perturbations. In this regime also only the i = p, pk components contribute to the vector modes We see that once again all the components decay as a −2 but with an oscillating behaviour. Also in this case only the δ A s actually contribute in the average. In both regimes, even though we have a source for the vector modes, they actually decay in the same fashion as in standard CDM, so that we do not expect large contributions at late times, unless they were produced with very large initial amplitudes.
To summarize this section, we have seen that in both regimes the i = a component contributes to the scalar but not to the vector perturbations, whereas for the i = p, pk components the situation is the other way around, contributing to the vector perturbations only. In addition, in the k ∼ H regime perturbations behave exactly as in standard CDM, whereas in the k ∼ (Hma) 1/2 regime we find a different behaviour implying that all the scalar perturbations decay with expansion in the same way as in scalar field DM, but unlike the scalar case, a small but non-vanishing gravitational slip is generated for the i = a component.

VIII. TENSOR PERTURBATIONS
So far we have neglected the tensor perturbations in the average equations. In order to extract the equations for such modes we use the projector Thus, contracting with Einstein equations we obtain We will calculate the components of the tensor perturbation in the orthonormal basis defined by: û 1 =û pk ,û 2 = cos θû p − sin θû a ,û 3 =k = sin θû p + cos θû a . In this basis the tensor perturbation takes the standard form, From the projection of Einstein equations we reacḧ A. Particle regime (k ∼ H, i.e. k/(ma) ∼ ) In this regime the sources vanish to the leading order so that the generation of gravity waves will be negligible. In this regime, the average sources read Thus we see that the i = a modes associated to the scalar perturbations are not purely scalar, but rather scalar-tensor modes whose tensor components have only + polarization. The vector perturbations with i = p are actually vectortensor modes also with + polarization and finally the i = pk component generates vector-tensor perturbations with × polarization. Redefining the field as h (+,×) = a −1h (+,×) , we can obtain solutions in terms of the Green's functions: with solution, Let us consider for example h × and assume C pk1 = 0 In this regime k H so that sin(k(η − η )) in the integrand oscillates rapidly whereas S (+,×) (η ) evolves slowly with time, so that can use partial integration as in (14) so that the leading term will be (100) Thus, we obtain waves with an amplitude that decays as a −1 and propagate at the speed of light. A completely analogous result can be obtained for the h + polarization.

IX. GRAVITATIONAL WAVE DETECTION
As shown above, scalar perturbations given by the C a, (1,2) components generate a gravity wave background with h + polarization in the k ∼ (Hma) 1/2 regime. If all the cosmological DM is generated by the vector field, it is possible to estimate the spectrum of gravity waves associated to such components and compare with the sensitivity of present and future detectors.
Note that for the mentioned components in this regime, the source S + in (95) can be written in terms of the scalar perturbation Φ given in (77) as Thus, from (98), integrating by parts as in the previous section we get: where in order to simplify the calculation we have assumed an instantaneous change to matter domination at equality, Ω m (a eq ) = Ω rad (a eq ), for both background and perturbations. Thus, we set the initial amplitude of the gravity waves to zero at equality. The term in brackets oscillates with an amplitude that decays as a −2 so that it is dominated by the lower integration time. Evaluating it at η = η eq , we obtain h + (η, k) = k 2 2m 2 a(η) sin 2 (θ) 1 − 4 cos 2 (θ) cos (k (η − η eq )) a −1 eq Φ(η eq , k).
We can see at equality when the amplitude of the gravity wave is largest, we have which is O( ). This is the reason why we could neglect the contribution of tensor modes in the evolution of scalar and vector perturbations in Section V. In this regime, with k k eq where k eq = H eq = 0.073 Ω m h 2 Mpc −1 it is possible to obtain Φ(η eq , k) directly from the linear scalar transfer function [70] Φ(η eq , k) = 9 10 Φ prim T (k) 9 10 Φ prim 12k 2 eq k 2 ln with Φ prim the primordial amplitude of perturbations generated during inflation, with a spectrum where the values of the parameters, ln(10 10 A s ) = 3.089 ± 0.036, n s = 0.9655 ± 0.0062 and the pivot scale k 0 = 0.05 Mpc −1 correspond to Planck observations [71].  [74,75]. The black dashed line show the upper bound limit of the massive vector gravitational waves production, as it can be seen its detection is unlikely with the future detectors.
Finally, the energy density today of gravitational waves with + polarization per energy interval and solid angle unit reads [72,73], Integrating over the whole solid angle we obtain the spectral energy density, In Fig. 4, we compare the prediction of the vector field DM model with the sensitivity of present and future gravity wave experiments. The best sensitivity at low frequencies correspond to the CMB data, but unfortunately the spectral range only reaches k = k eq just in the limit of the wave regime. On the other hand, at higher frequencies, SKA pulsar timing limit is twelve orders of magnitude above the production prediction.
If we integrate over k in the gravitational wave production band we get where we have approximated n s = 1 for simplicity. In Fig. 5 we plot Ω GW (η 0 ) as a function of m. The expected sensitivity for the combined analysis of the future COrE and Euclid missions [76] on the effective number of relativistic degrees of freedom can be translated into a limit on gravity wave abundance As can be seen in Fig.5 the maximum production corresponds to masses m ∼ 10 −27 eV, but still it is a few orders of magnitude below the mentioned limit.

X. CONCLUSIONS
Ultralight bosonic fields are natural DM candidates which can avoid some of the small-scale problems of the standard CDM model. Most of the work developed so far in this field has focused on the simplest implementation of this scenario based on scalar fields. In this work we have considered the case of ultralight vector fields.
The first difficulty in the higher spin case already appears at the background level, since such fields generically break isotropy. Fortunately, a general result [56] shows that for massive fields without self interactions and masses much larger than the expansion rate, the average energy-momentum tensor is isotropic and behaves as CDM.
At the perturbation level, we have considered an adiabatic expansion in two different regimes, the so called particle regime with k/(ma) ∼ and the wave regime with k/(ma) ∼ 1/2 . Very much as in the scalar case, we find that for vectors, the particle regime is indistinguishable from CDM. However, in the wave regime important differences with respect to the scalar case arise. Thus, perturbations in the vector field support three kinds of metric perturbations. On one hand, we have two scalar modes with a small but non-vanishing sound speed c 2 s ∼ k 2 /(m 2 a 2 ) which suppresses structure formation for k > √ Hma. Such modes generate a non-vanishing gravitation slip of order (Φ − Ψ)/Φ ∼ k 2 /(m 2 a 2 ) which is a specific feature of this ultralight vector field DM model. In addition, the scalar modes source tensor modes with a small amplitude h/Φ ∼ k 2 /(m 2 a 2 ) and a characteristic spectrum which peaks around k max ∼ ma eq . The amplitude of this gravity wave spectrum is however below the sensitivity of present and future detectors. Nevertheless, the calculation done in this work has been conservative in the sense that we have focused only in the potential generation of gravity waves in the matter dominated era. A complete study would require to consider also perturbations in the radiation era. On the other hand, we have four vector modes which are also sources of gravity waves. The vector modes decay as a −2 , i.e. in a similar fashion as in standard CDM cosmologies, so that we expect a negligible amount of vector modes at late times also in this model. In Fig.6 a summary of the perturbations behaviour in the different regions of the spectrum is shown for massive vector and scalar models. In the standard CDM scenario, DM can be the source of scalar and vector perturbations. In this case there are only two relevant regimes, namely, sub-and super-Hubble modes. Coherent scalar DM can only source scalar perturbations, but their scaling depends not only on the wavenumber but also on the mass of the field, defining a total of four different regimes. For coherent vectors, we have the same regimes as for the scalar, but in this case vector and tensor perturbations can also be sourced, as well as a gravitational slip. Both massive scalar and vector fields mimic CDM in yellow and green regions.