Gravitational production of vector dark matter

A model of vector dark matter that communicates with the Standard Model only through gravitational interactions has been investigated. It has been shown in detail how does the canonical quantization of the vector field in varying FLRW geometry implies a tachyonic enhancement of some of its momentum modes. Approximate solutions of the mode equation have been found and verified against exact numerical ones. De Sitter geometry has been assumed during inflation while after inflation a non-standard cosmological era of reheating with a generic equation of state has been adopted which is followed by the radiation-dominated universe. It has been shown that the spectrum of dark vectors produced gravitationally is centered around a characteristic comoving momentum $k_\star$ that is determined in terms of the mass of the vector $m_X$, the Hubble parameter during inflation $H_{\rm I}$, the equation of state parameter $w$ and the efficiency of reheating $\gamma$. Regions in the parameter space consistent with the observed dark matter relic abundance have been determined, justifying the gravitational production as a viable mechanism for vector dark matter. The results obtained in this paper are applicable within various possible models of inflation/reheating with non-standard cosmology parametrized effectively by the corresponding equation of state and efficiency of reheating.


Introduction
One of the outstanding puzzles of high energy physics is to understand the nature of dark matter (DM). There is overwhelming evidence that the dominant component of matter density in the universe is due to DM. However, all the observational pieces of evidence in favor of DM have been established by its interactions with the Standard Model (SM) only through gravity [1][2][3][4]. A wide range of DM models have been proposed with different DM production mechanisms depending on the DM interaction strength with the SM and/or beyond the SM new physics, for a review see [5] and original references therein. Arguably the most popular DM models involve weakly interacting massive particles (WIMPs) whose mass ∼ O(TeV) and interaction strength ∼ O(0.1) are similar to those of the SM particles. The production mechanism for WIMPs is a thermal freeze-out scenario where it is assumed that the SM and DM were produced in the early universe during the reheating phase. The WIMP DM was in thermal equilibrium with the SM and as the universe cooled down it went out of the equilibrium to constitute the observed relic abundance. Thermal freeze-out and other DM production mechanisms [5] are theoretically appealing, however there has been no sign of DM interactions with the SM in a variety of experiments thus far, apart from its gravitational interactions.
The lack of DM signals motivates us to envision novel mechanisms of DM production where the DM and SM interactions are absent or negligible. One such mechanism is the particle production due to quantum fluctuations in a rapidly expanding universe [6,7]. An inflationary scenario in the early universe not only explains successfully some of the puzzles of modern cosmology (for a review see [8]), it can also provide an ideal phase for DM production due to quantum fluctuations [9,10]. Furthermore, a very heavy DM can also be produced during the (p)reheating period -the last stage of inflation -where the inflaton field transfers its energy density to the visible and dark sectors through coherent oscillations [11][12][13]. In particular if the reheating does not happen instantaneously the maximum temperature achieved during the reheating phase can be larger than the temperature at the end of reheating which allows supper-heavy DM production [14]. Recently, there has been a renewed interest in models of gravitational production of DM during and after inflation due to quantum fluctuations [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33].
In this work, we study a minimal model of gravitationally produced DM where the DM is a massive Abelian gauge boson X µ . In particular, we focus on vector DM production due to the quantum fluctuations during and after inflation in an early era of non-standard cosmology parameterized effectively by the equation of state w. We assume the dark vector field couples minimally to gravity and it has no other interactions with the visible sector. The vector DM mass is assumed to be generated via the Stueckelberg mechanism and it is non-zero during and after inflation.
A massive vector field X µ has three physical polarizations: two transverse X ± and one longitudinal X L , whereas X 0 is an auxiliary component. It is well known that the transverse components of a minimally coupled vector field produce scale-invariant density fluctuations during inflation [15,[34][35][36] therefore they can not constitute all the observed DM relic density. However, if the vector DM is non-minimally coupled to gravity or the inflaton field then the transverse components of a vector field may constitute all the observed DM density, see e.g. [18][19][20][21] for the latter case. On the other hand, the longitudinal modes in a minimally or non-minimally coupled scenarios do not lead to scale-invariant density fluctuations during and after inflation [15,24,27] and hence can account for all the observed DM abundance. Since we are interested in a minimally coupled vector DM scenario, therefore our focus is confined to the production of the longitudinal modes of the vector DM during and after inflation.
The main new feature of our work in comparison to the previous works on the gravitational production of vector DM is an early epoch of non-standard cosmology during the reheating period parameterized by a general equation of state w of the inflaton/reheaton field. During inflation, without specifying the inflationary dynamics, we assume a de Sitter phase with constant Hubble parameter H I . After the end of inflation a reheating phase starts where a non-standard cosmological evolution is assumed with the dominant energy density scaling as a −3(1+w) , where a is the scale factor and −1/3 < w < 1. This specific range of w during reheating period is motivated by our analytic analysis which is robustly verified against exact numerical calculations with specific values of w = {−1/4,−1/6, 0, 1/3, 2/3}. At the end of reheating phase, with non-standard cosmology, the standard cosmological evolution resumes with the radiation-dominated (RD) era followed by the matterdominated (MD) universe. We define the end of the reheating phase when the energy density of the non-standard cosmology ∝ a −3(1+w) is dominated by the SM radiation energy density ρ RD ∝ a −4 . The production of dark matter during the non-instantaneous reheating phase due to the coherent oscillations/decays of the inflaton field are studied in Refs. [37][38][39][40][41], however we neglect such effects here. Furthermore, we do not specify the dynamics of reheating mechanism, however we assume that such mechanism exists and leads to a non-standard cosmology where the dominant energy density scales as a −3(1+w) and the SM radiation energy density during this phase is subdominant. In recent years there has been many works done on DM scenarios in the presence of an early universe non-standard cosmology [42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60]. Our focus in this paper is to study, for the first time, the production of vector DM purely due to quantum fluctuations during and after inflation employing non-standard cosmological period of reheating in the early universe.
In Sec. 2 we provide a detailed description of the gravitational production of a minimally coupled vector DM in an expanding universe. Focusing only on the longitudinal modes, we give approximate analytic and exact numerical solutions for the mode functions during and after inflation for different regimes of DM mass and wavelengths. In particular, we note that for certain mass and wavelength range the frequency squared of the longitudinal modes becomes negative leading to a tachyonic enhancement of those modes. We calculate the vector DM relic abundance in Sec. 3 in two DM mass ranges, H rh ≤ m X < H I and m X < H rh , where H rh is the Hubble scale at the end of reheating with non-standard cosmology. In Sec. 4 we present our conclusions. Supplementary material including details of the quantization of a vector field on the curved background is given in Appendix A.

Gravitational vector DM production in expanding universe
During an epoch of rapidly expanding universe (inflation and partially reheating) the gravitational production dominates over other DM production mechanisms. In this work we consider renormalizable (dim-4) Lagrangian for the vector DM minimally coupled to gravity. The action for the vector DM in a background metric g µν is given by where the background metric g µν is of the FLRW form with the line element whereas, X µν ≡ ∂ µ X ν − ∂ ν X µ denotes the field strength of the vector field X µ . The mass for the dark vector boson m X is generated via the Stueckelberg mechanism 1 . Now, using the gravitational definition of the energy-momentum tensor one can find the energy density of the vector DM as where X 0 and X ≡ X i are components of the vector field, the over-dot denotes derivative w.r.t. t and ∇ ≡ ∂/(∂ x).
In Appendix A, we provide details of the quantization of a vector field in a curved background. Adopting the notation from the appendix, we write the vacuum expectation values for the longitudinal (L) and transverse (±) components of the energy density as 1 Alternatively one may assume an Abelian Higgs mechanism such that the extra neutral Higgs boson (radial mode) is very heavy. Then it is possible to show that, in an appropriate limit of constant mX , effectively the Higgs model reduces to the Stueckelberg mechanism. However, this Higgs mechanism has to be effective already during the inflationary period.
where ρ L ≡ 0|:ρ L : |0 and ρ ± ≡ 0|:ρ ± : |0 2 and Above the X L (τ, k) and X ± (τ, k) are Fourier transforms of the vector-field longitudinal X L (τ, x) and transverse X ± (τ, x) components, respectively (see Appendix A). The prime denotes derivative w.r.t. to the conformal time τ , defined as dt = a(τ )dτ . Equations of motion for the longitudinal and transverse Fourier modes, i.e. X L (t, k) and X ± (t, k), are where the time-dependent frequencies are given by Note that the transverse mode frequency ω 2 ± (τ ) is always positive, however, the longitudinal mode frequency ω 2 L (τ ) can be negative in some regions of the parameter space. Therefore, for ω 2 L (τ ) < 0, the production of longitudinal modes would receive the tachyonic enhancement. Hence the production of the longitudinal modes can be parametrically larger than that of the transverse modes, see also [15,24,27]. In the following, we only focus on the gravitational production of the longitudinal modes in the expanding universe. In particular, we will present both numeric and approximate analytic solutions for the longitudinal mode functions X L in two phases, (i) during slow-roll inflation, and (ii) after inflation during the reheating phase with non-standard cosmology followed by the RD universe.
Cosmological evolution: First, let us focus on the background dynamics. We assume that during the inflationary period the evolution of the scale factor is approximated by the de Sitter universe with an exactly exponential variation of the scale factor, a I (t) ∝ e H I t , with H I being constant Hubble parameter defined as H ≡ȧ/a = a /a 2 during inflation. In general, the evolution during inflation is governed dynamically by a homogeneous inflaton field that rolls slowly down in a properly adjusted potential interacting with gravity [8]. This more elaborate approach is not necessary here, as it would not alter the qualitative results presented in this work. Hence, neglecting corrections from the inflaton dynamics, the inflationary stage can be well approximated by the de Sitter scale factor a I (τ ) ∝ −1/(τ + const.) in the conformal coordinates. Inflation starts at conformal time τ i → −∞ and lasts up to τ e 0 + . Then, at τ = τ e the universe smoothly (in the presence of the inflaton/reheaton field φ) transfers from the de Sitter stage into the reheating period which employs non-standard cosmology. During the reheating the energy density of the universe is dominated by the energy density ρ φ of the inflaton/reheaton field φ. In this period, we assume that the inflaton equation of state is parametrized by a parameter w describing the non-standard cosmology that we will vary and investigate its relevance: where e.g. w = 0 and w = 1/3 correspond to an early MD and RD universe, respectively. Strictly speaking, we consider transfer of energy density from ρ φ to the SM (radiation) density ρ SM during the epoch of non-standard cosmology until the very end of this period when ρ SM ≈ ρ RD ∝ a −4 dominates and standard cosmological universe resumes with RD era. In particular, the evolution of non-standard cosmology is governed by the following set of Boltzmann equations for ρ φ and ρ SM , where Γ φ is the rate of transfer of inflaton energy to the SM. Assuming ρ φ ρ SM and Γ φ 3(1 + w)H, the generic form of the equation of state for the inflaton/reheaton (2.10) implies that during the reheating period the scale factor a RH evolves as a function of the conformal time τ as (2.12) The reheating period is followed, in turn, by the RD epoch when the total energy density is dominated by the SM radiation during which the scale factor varies as a RD ∝ τ . We require the scale factor a(τ ) and the Hubble rate H(τ ) to be continuous at the two transition points, i.e. where τ rh refers to the conformal time at which the transition from the reheating to the RD era happens 3 . The continuity of the scale factor and Hubble rate imply rh , τ rh < τ . Consequently, we can express the Hubble rate as a function of the scale factor a, where a e ≡ a(τ e ) and a rh ≡ a(τ rh ) define the scale factors at the end of inflation and reheating periods, respectively. The total energy density ρ(a) = ρ φ + ρ SM in these different epochs during the evolution of universe is, , a e < a ≤ a rh 0, a rh < a , where we have used Γ φ = (5−3w)H rh /2, such that H rh ≡ H(a rh ) defines the end of reheating period when ρ φ (a rh ) = ρ SM (a rh ). As mentioned in Introduction we confine w ∈ (−1/3, 1), therefore during the non-standard reheating period the SM radiation energy density ρ SM approximately scales as a − 3(1+w) 2 (the second term in the parenthesis in the middle line of (2.20) proportional to a −4 is negligible). Hence the temperature during the reheating period scales as a − 3(1+w) 8 , this is a useful result for later use. Note that the SM energy density ρ SM is zero during the inflationary period and is continuous at the end of reheating since a −4 term in the middle line of (2.20) is negligible in comparison to a −3(1+w)/2 term for −1/3 < w < 1.
Finally, it is convenient to define reheating efficiency γ as which parametrizes the duration of reheating with non-standard cosmology. By fixing γ and w, we can express a rh in terms of a e , or equivalently τ rh in terms of τ e , as  The efficiency of reheating γ is a constrained parameter. Requiring H rh < H I sets an upper limit on γ < 1. Whereas the lower limit on γ 10 −18 can be deduced by using the fact that curvature perturbations during inflation constrain H I ≤ 6.6 × 10 13 GeV at 95% C.L. [61] and Big Bang Nucleosynthesis (BBN) sets the lower limit on the reheating temperature T rh 10 MeV [62], (see Eq. (3.28)).
Evolution of the Hubble rate H(a) and the energy density ρ(a) as a function of the scale factor a are illustrated in Fig. 1. During the reheating period, we consider nonstandard cosmological evolution parameterized by w ∈ (−1/3, 1). The end of inflation a e , end of reheating a rh , and the end of RD, i.e. the matter-radiation equality (mre) a mre , periods are shown as gray dashed vertical lines. In the left-panel, the gray dotted horizontal lines represent the Hubble rate at the end of inflation H I , the Hubble rate at the end of reheating H rh and the Hubble rate at the end of RD epoch H mre . In the right-panel, we show the total energy density ρ(a) = ρ φ + ρ SM as solid curve. The inflaton energy density ρ φ is shown as the red curve, whereas the SM energy density ρ SM is shown as green curve. Note that the ρ SM is zero during the inflationary period, however, it has non-zero value proportional to a − 3 2 (1+w) during the non-standard reheating, shown as dashed green curve.

Longitudinal modes during inflation
During the slow-roll period of inflation the Hubble rate H I is almost constant, so we adopt the de Sitter solution. The relation between the scale factor and conformal time is given by (2.15), . (2.23) Consequently, during the de Sitter stage the longitudinal frequency ω 2 L (2.8) simplifies as where we have used the fact that during the slow-roll period of inflation It is seen that for vector DM mass m X H I , ω 2 L remains positive with no chance for the tachyonic enhancement. Furthermore, another reason not to consider vector DM heavier than the inflationary scale H I is the fact that we assume during inflation the inflaton energy density 3M 2 Pl H 2 I to be the dominant energy density. Therefore, hereafter we consider only vector DM masses m X H I .
In the far past during inflation (subhorizon limit) all modes were deep inside the horizon, so that a I (τ )m X a I (τ )H I k. In this limit, the frequency ω 2 L becomes constant which results in a simple harmonic oscillator equation for modes, The unique solution that minimizes the energy is known as the Bunch-Davies state, Hereafter we will take the above solution as our initial condition. For the vector DM mass m X H I , there are two other regimes relevant during inflation: In Fig respectively. Early enough both these modes were in the red region (sub-horizon) and they satisfied the initial condition provided by (2.26) known as the Bunch-Davies vacuum. Next they enter the intermediatewavelength region (purple) and before the end of inflation at a e the mode with momentum k 1 crosses the Compton wavelength (am X ) −1 at τ 2 and enters the long-wavelength regime (blue) region. We solve the mode equation in different regimes of k until the end of inflation.
After that, at a e , we match solutions found during inflation with those during the reheating phase, that will be discussed in the next subsection 2.2.
Starting with the first (purple) region during inflation, i.e. intermediate-wavelength a I (τ )m X k a I (τ )H I , the equation of motion for modes takes the form Log a ae arh Figure 2: Evolution of various cosmological distances during and after inflation for heavy vector DM i.e. H rh ≤ m X < H I (left diagram) and light vector DM m X < H rh (right diagram). The red region corresponds to modes with wavevector in range am X , aH k, purple refers to the region where am X < k < aH, blue corresponds to the condition k < am X < aH and in the green region aH, k am X . Here a c = k/H refers to the second horizon crossing, a ≡ a(τ ), k m ≡ a e m X , k = a m X , k e ≡ a e H I and k rh ≡ a rh H rh . The plot assumes −1/3 < w < 1/3 during the reheating phase.
The solution to the above equation is given by where the integration constants are obtained, by imposing the Bunch-Davies initial condition (2.26), as Hence the mode solution takes the form where the approximation in the second line above (k a I H I ) implies the modes are well inside the intermediate wavelength (purple) region. In this limit longitudinal modes grow linearly with the scale factor a I (τ ) until the mode momentum k crosses the Compton wavelength, i.e. k = a I (τ )m X line. This evolution of the longitudinal modes of the vector DM is similar to a massless scalar.
In the second case with long-wavelength k a I (τ )m X a I (τ )H I (blue region), the mode equation (2.7) reduces to where J n (x) and Y n (x) are the Bessel functions of the first and second kind, respectively. The order of Bessel functions n is Consequently, the mode X L simplifies as follows We fix coefficients C 2 by requiring, where τ 2 is defined by the condition k ≈ a I (τ 2 )m X , i.e.
hereb is an auxiliary parameter introduced to improve the agreement between numerical and analytical solutions, in the following calculationsb = 1.2 has been adopted. For the coefficients we find, Expanding the above result in the limit k a I (τ )m X a I (τ )H I we can approximate Eq. (2.36) by the following constant solution The factor (sin (b) +b cos (b))/b 2 ≈ 0.95 forb = 1.2 and could be safely neglected. Note that in this case, with k a I (τ )m X a I (τ )H I , we have ω 2 L > 0 (2.33). Hence there should be no tachyonic enhancement for k values in this regime, as it is indeed confirmed by the above mode function being constant. Note however that in the long-wavelength region (blue) in Fig. 2 the constant mode function amplitude square | X (b) L | 2 scales as k −1 which is different from that of a massive scalar field which scales as k −3/2 .
In Fig. 3 we show | X L | 2 and the frequency ω 2 L as a function of the scale factor a during inflation for m X = 10 8 GeV, H I = 10 13 GeV, and two values of the momentum k = 10 −8 GeV (left-panels) and k = 10 −3 GeV (right-panels) representing the two regimes k a e m X a e H I and a e m X k a e H I , respectively. In the upper panels the approximate analytical solutions (2.26), (2.32) and (2.41) are compared with exact numerical solutions. As it is seen the agreement is excellent. The vertical gray dashed lines show the value of the scale factor corresponding to the radius of the Hubble sphere during the inflation a = k/H I and a = k/m X . At those points we have matched the approximate analytical solutions as discussed above. In the lower panels we plot ω 2 L and show its zeros by vertical green dashed lines. Note nearly perfect correlation between enhancement (tachyonic) of | X L | 2 and regions of negative ω 2 L . This behavior has been anticipated. Our goal is to find energy density of the longitudinal modes X L during inflation, where ρ L denotes vacuum expectation value of the energy density calculated with respect to the Bunch-Davies vacuum. In the case of the sub-horizon modes a I (τ )H I k (red region in Fig. 2) energy density per log momentum is given by The first three terms are divergent in the limit k → ∞ while integrating over k. However, since we are going to use a cut-off Λ = a e H I for the integration, that will eliminate contributions of those sub-horizon modes.
For modes with a I m X k a I H I (purple region in Fig. 2) the energy density per log momentum is (2.44) The modes with k a I m X a I H I (blue region in Fig. 2) contribute to the energy density per log momentum as, (2.45) Therefore, at the end of inflation, a e ≡ a(τ e ), the energy density per ln k of the longitudinal modes with m X H I scales as, sub-horizon modes k a e H I (red), , super-horizon modes k a e H I (purple, blue). (2.46)

Longitudinal modes after inflation
In this subsection our goal is to find solutions to the equation of motion (2.7) for the longitudinal modes after inflation. In particular, we do not assume standard cosmological evolution, in which the inflationary period is followed by the RD epoch until the matterradiation equality (mre). Instead, we allow the possibility for a non-standard cosmological evolution after the end of inflation by assuming a general equation of state (2.10) such that the dominant energy density scales as a −3(1+w) . The period of non-standard cosmology, referred to as reheating, is then taken over by the standard cosmology with the RD universe. After inflation, for a > a e , the longitudinal mode frequency ω 2 L (2.8) takes the form  .17), respectively. Note that the longitudinal mode frequency ω 2 L (2.8) contains a term proportional to a , which is discontinuous across the two transition points τ e and τ rh , therefore the longitudinal mode frequency also has a jump at these two points. Technically, discontinuity disappears by use of the Friedmann equations (2.48), however, in some sense, it is still present due to the fact that we assume an instantaneous transition from the inflationary period to a non-standard cosmological epoch of reheating parameterized by w and similarly for the transition from non-standard cosmological epoch to RD universe. In a case of describing the inflation and reheating by a dynamical inflaton/reheaton field, the discontinuity would disappear as w would be a continuous function of the inflaton/reheaton field. However, we would like to emphasize that the energy density (2.42) contains only terms proportional to a, a and X L , X L . Since the scale factor and the Hubble rate are continuous across the transition points, by requiring the mode functions to match at τ e and τ rh : To find an approximate solution to the equation of motion for the longitudinal mode after inflation, we consider the following regimes in parameter space together with the corresponding frequency ω 2 L values: where we have indicated colors which, in Fig. 2, mark regions where a given formula for the frequency applies.
In the following analysis we restrict our considerations to modes that went outside the horizon during inflation since only these modes receive tachyonic enhancement, i.e. k < k e = a e H I . After the end of inflation, modes with mass m X H I can enter two regions (see Fig. 2): purple region, in which am X < k < aH, and blue region, which describes modes with k < am X < aH I . The modes with different k values evolve differently and they all eventually enter the green area in which am X k, aH is the dominant energy scale.
Our aim is to find approximate analytical solutions for the mode equation in different mass and wavelength regions after the end of inflation.

Heavy vector dark matter: H rh ≤ m X < H I
In the case of heavy vector DM we consider the following three regimes w.r.t. the wavelength of the modes.
Long-wavelength case (k k m ≡ a e m X ): After the end of inflation, modes with k < k m ≡ a e m X , for instance k −1 1 line in Fig. 2, pass through the blue region to the green area, where they eventually re-enter the horizon. At the early times, when the modes k a(τ )m X a(τ )H(τ ) are still in the blue region, the longitudinal mode frequency is negligible, i.e. ω 2 L → 0. Consequently, in this regime, we obtain the following solution to the mode equation  Since in this regime ω 2 L ≈ a 2 m 2 X , therefore the above condition implies which usually (except a short period just after τ = τ ) is satisfied in this region. Hence, we can use the following approximation to the solution of the mode equatioñ X (II) The τ dependance could be easily converted into a dependance on the scale factor a. To find C we should merge the above solution with X (I) L at τ = τ . The approximate form of these coefficients is The form of solutions in this region are oscillatory decaying as shown in the left-panel of Fig. 4 (dotted green) which match well with the exact numerical solution (solid cyan).

Intermediate-wavelength case (k m < k < k ≡ a(τ )m X ):
After the end of inflation modes with wavevector in the range k m < k < k , e.g. k −1 2 line in Fig. 2, pass through the three regions (purple, blue, green). Their post-inflationary evolution starts in the purple region, where am X k aH so that a 2 m 2 X /k 2 1. This implies that the second term in the parenthesis in Eq. (2.50b) can be neglected. Hence the mode equation reduces tõ with the solutionX We fix the coefficients C (III) 1 and C (III) 2 by merging the above solution at τ = τ e with X (a) L (τ e ) (2.32). However, as long as w < 1 the second term in (2.57) is small compared to the first one. Since we confine ourselves to the range −1/3 < w < 1, therefore, in the subsequent calculations we neglect this part of the solution. Hence the above mode solution reduces to So, in this regime (purple region) the longitudinal modes evolve linearly with the scale factor a. We show this approximate analytic solution as the dashed purple curve in the upper-middle panel of Fig. 4. Next, at τ =τ 2 = τ e (b a e m X /k) −(1+3w)/2 such that k =ba(τ 2 )m X , the mode crosses the Compton wavelength line (am X ) −1 and enter the blue area in Fig. 2 (left-panel). In this region, we have found an asymptotic analytical solution given by (2.51). In order to get coefficients parametrizing solution (2.51) for k m < k < k we merge solution (2.51) with (2.58) at τ =τ 2 . Note that for w > −1/3 the second term in Eq. (2.51) is growing. Contrary to the long-wavelength case, now the modes with k am X aH follow the purple region, in which mode functions increase in amplitude according to solution (2.58). One would naively expect that the growing term in (2.51) should dominate the solution following (2.58). However, exact numerical solution presented in upper-middle panel of  Fig. 4. The approximate solution in this region is in the form given by (2.55) with the following coefficients: . (2.60)

Short-wavelength case (k < k < k e ):
After inflation the modes with k < k < k e also start their evolution in the purple region (Fig. 2 left-panel), but unlike the previous cases they do not evolve into blue region. Instead at τ =τ 1 , they re-enter the horizon, i.e. the red area, in which k is the dominant energy scale. In this case, it is convenient to consider those two regimes collectively, i.e. in purple and red regions, we neglect in ω 2 L terms proportional to a 2 m 2 X and arrive at the following equation of motion, with the solution where the order of Bessel functions is ν ≡ −3(1−w) 2(1+3w) . Merging the above solution (2.62) with X : where k e ≡ a e H I . This analytic solution in the purple and red regions is shown in the right-panel of Fig. 4 as dashed line in their respective colors. Then at τ =τ 3 the modes enter the green region with the solution given by Eq. (2.55). Again after some straightforward but tedious calculation we find that both coefficients appearing in (2.55) are non-zero and the solution is decaying oscillatory as shown in the right-panel of Fig. 4 (dotted green).
In lower panels of Fig. 4 we present frequency ω 2 L as a function of the scale factor a in order to verify the expected correlation between regions of strong increase of mode functions and negativity of ω 2 L . The correlation is indeed confirmed in Fig. 4. Furthermore, we would like to emphasize that the solutions for mode functions and therefore for energy densities obtained here are continuous functions, that can be seen just from inspection of upper panels of Figs. 3 and 4. However a small discontinuity appears is plots of the frequency ω 2 L as a function of a, a careful reader may notice it comparing ω 2 L at a e ≈ 10 −12 in lower panels of Figs. 3 and 4.

Light vector dark matter: m X < H rh
The above approximate solutions are also valid for the lighter vector DM case. As shown in the right-panel of Fig. 2, in this case, all modes are in the purple region at the beginning of reheating. Previously, we have demonstrated that in this region modes increase in amplitude according to the solution (2.58). This growth is eventually terminated at the point when modes go outside this region. For modes with k > a rh m X it happens during the RD epoch. In this case, one can simplify the above results using the fact that in this era w = 1/3. However, for the modes with k < a rh m X the transition from the purple region to blue happens during the reheating phase and results are similar to the heavy vector DM case discussed above. The modes k < k < k rh , e.g. k 2 in the right-panel of Fig. 2, cross the horizon (transition from the purple to red region) during the RD epoch. In this case the solution (2.62) is valid, however with the fixed value of equation of state w = 1/3. On the other hand, the evolution of the mode functions with k rh < k < k e during the reheating period is same as for the heavy DM regime discussed above.

Energy density scaling
Before closing this section, it is instructive to get the approximate analytic results for the energy density per ln k in each colored region of Fig. 2. Employing the approximate analytic mode function solutions, we use Eq. (2.4) to get the energy density per log momentum after the end of inflation and its scaling with a in different regions. In the red region where k is the dominant energy scale one can approximate the mode solution (2.62) as which implies that energy density redshifts as, In the regime with am X k aH, i.e. purple region, we approximate solution to mode equation by (2.58) and the energy density scales as (2.66) Note that the same scaling would be obtained if we used (2.62) in the limit k aH. In the blue region k am X aH, using Eq. (2.51) with C (I) This scaling of energy density of the longitudinal modes ∝ a −2 is a very unique feature of the vector DM, see also [15]. For instance, the energy density of a massive scalar DM scales as constant in this blue region which leads to enhanced iso-curvature perturbations at large scales, such large iso-curvature perturbations are severely constrained by the CMB data [61]. Hence, the damping of the energy density ∝ a −2 observed here is a welcome property of the vector DM that implies suppression of unwanted iso-curvature perturbations at the large scales. Finally, in the green region where am X is the dominant energy scale we obtain the following scaling where we have neglected the oscillatory terms proportional to exp [±i dτ a(τ )m X ]. Note that is this region the energy density behave as of the matter density.
Log a ae arh The above red-shifting of the energy density and its scaling w.r.t. mode momentum k can be summarized as, where the colors represent different regions as in Fig. 5. Note that the dependance of energy density per ln k on momentum k is w-dependent for the modes with k a(τ )m X , a(τ )H(τ ), i.e. the red region. However, if these modes are during the inflationary (de Sitter) period then w = −1 and if they are in the RD epoch then w = 1/3. Whereas, during the reheating phase the equation of state parameter w takes the value in the range w ∈ (−1/3, 1).
To summarize the scaling of energy density w.r.t. the scale factor a, we present the exact numerical results for different wavelength modes discussed above in Fig. 6. We note that as the universe expands the redshift of the energy density varies for different modes and the approximate analytic scaling of the energy density (2.68) matches well with those of the exact numerical results. We see that at the end of inflation (a = a e ) the main contribution to the total energy density comes from modes with the shortest wavelength (k ∼ k e ). However, after the end of inflation, those modes receive the strongest suppression proportional to a −4 . On the other hand, modes with longer wavelength initially have a smaller contribution to the total energy density, but their energy is also redshifted to a lesser extent. The intermediate-wavelength modes contribute the largest to the energy density as noted above.

Relic abundance of vector dark matter
The present energy/number density of DM particles could be expressed in terms of energy/number density at H(a ) = m X . The number density of longitudinal modes, n L , is related to their energy density ρ L as where E L denotes single-particle energy. Consequently at H(a ) = m X the number density per ln k is computed as, In the following, we evaluate the DM number density in two mass regimes, the heavy vector DM H rh ≤ m X < H I and the light vector DM m X < H rh . Moreover, it should be stressed that we consider the evolution of number density for super-horizon modes with k < k e , i.e. the modes which were outside of the horizon at the end of inflation. Before presenting detailed numerical analysis it is instructive to estimate the number density per ln k (3.2) of the longitudinal modes at H(a ) = m X , from the energy density scaling observed in Eq. (2.68). In Fig. 5, we show the d n L (a ) /d ln k scaling with k for the longitudinal modes at a = a(τ ). Note that the dominant number density is generated by the modes with comoving momentum k represented as dashed (red) line. After exiting the horizon, these dominant modes pass through the purple region only until they reach H(a ) = m X . Focusing first on the heavy vector DM case (Fig. 5 left-panel), the modes with wavelengths longer that k −1 , i.e. k k , would have to pass through the blue region where the energy density of these modes redshifts as a −2 . Therefore in this region the energy density associated with these longer wavelength modes k k scales proportional to k 2 . Hence the long wavelength longitudinal vector modes k < k lose their power proportional to k 2 and therefore are safe from generating dangerous iso-curvature perturbations at the CMB scales. This result would remain same for the number density of modes in the blue region, i.e. d n L (a ) /d ln k ∝ k 2 . On the other hand, the shorter wavelength modes, k > k , re-enter the horizon at a < a and follow the evolution in the red region (Fig. 5 left-panel) where the energy density of these modes redshifts as radiation i.e. a −4 . Since these shorter wavelength modes k > k re-enter the horizon during the reheating phase when the Hubble rate scales are a −3(1+w)/2 , therefore the number density of the modes scale as k − 3(1−w) (3w+1) . Note that for w < 1, the number density per ln k at a has a peak structure and the dominant modes are k where most of the number density is contained. This peak structure was first noted in [15] for instantaneous reheating and RD universe after inflation. Here, we extended this observation for non-standard cosmology with equation of state w < 1.
Similarly for the light DM case, m X < H rh , we observe that the dominant modes (shown as dashed red line in the right-panel of Fig. 5) are those with k = k . The number density for longer wavelength modes (k < k ) scales as k 2 as they pass through the blue region to reach a = a(τ ). However, in the light DM case, the modes with wavenumber k rh < k < k re-enter the horizon during the RD universe (H(a) ∝ a −2 ) and their energy density also redshifts like the radiation a −4 . Therefore, the number density of the modes with k rh < k < k scales as k −1 . The modes with comoving momentum k e < k < k rh re-enter the horizon during the period of reheating, i.e. H(a) ∝ a −3(1+w)/2 and their energy density (in the red region) scales as a −4 , and the number density d n L (a ) /d ln k of these modes scales as k − 3(1−w) (3w+1) . Hence, for w < 1, the number density per log momentum has a peak structure with maximum values around k modes. In the following we present detailed numerical analysis of this qualitative discussion by calculating the number density of the vector DM.

Number density of the heavy vector DM: H rh ≤ m X < H I
For vector DM in the mass range H rh ≤ m X < H I , we consider two cases for the momentum vector k: The number density for modes with k < k < k e can be written as where a c corresponds to the horizon crossing and it is defined as . (3.4) Note that in this case the horizon crossing occurs during reheating, therefore, one can express a c as a c | k <k<ke = a e H I a c k . (3.5) Consequently, in this limit we can rewrite Eq. (3.3) as . (3.6) In the other case, k < k , the number density at a = a is given by . (3.8)

Number density of the light vector DM: m X < H rh
Now we focus on the case when the vector DM is lighter than the Hubble scale at the end of reheating, i.e. m X < H rh . In this scenario, we calculate the DM number density for the following three cases: where k e ≡ a e H I , k rh ≡ a rh H rh . The scale factor at the end of reheating a rh ≡ a(τ rh ) is given in terms of a e and γ in Eq. (2.22).
In the case (a), for large k values, i.e. k rh < k < k e , we get the number density as , (3.10) Hence the number density for the longitudinal modes for the case (a) can be rewritten as, . (3.11) In the case (b), i.e. k < k < k rh , we get the number density for longitudinal modes as where again a c corresponds to value of the scale factor at the second horizon crossing and for light vector DM this happened during the RD epoch. Therefore, (3.14) Hence the number density in the case (b) is a e k .
(3.15) Finally, for the longitudinal modes with k < k (c), we obtain the number density as (3.17) We can summarize our results for number density for the longitudinal modes at a = a in both mass regimes as follows: • For heavy vector DM mass H rh ≤ m X < H I , for k < k < k e , for k < k .
• For light vector DM mass m X < H rh , a e k , for k < k < k rh , for k < k .
In Fig. 7 we show results of exact numerical computation of the number density per ln k as a function of k/k for different values of w at a = a . Note that in the above two vector DM mass regimes, the number density per log momentum has a peak structure if and only if w ∈ (−1/3, 1). A similar peak structure was also observed in [15] with standard cosmological history assuming instantaneous reheating and radiation dominated universe after the end of inflation. In this case, d n L (a ) /d ln k is dominated by modes with k ∼ k . In Fig. 8 we compare the exact numerical results (solid curves) with the approximate analytic form [(3.18a)-(3.18b)] (dashed curves) of the number density per ln k as a function of k for w = 0, 1/3 at a = a . As seen from the plots the exact numerical results agree very well with the approximate analytic solutions.
The total number density n L (a ) for the case of heavy vector DM, H rh ≤ m X < H I , at the moment when H(a ) = m X is given by, (3.20) where in the last approximation we have assumed that w ∈ (−1/3, 1) such that H I m X 1−w 1+w 1. In this case of heavy DM, the dominant mode momentum k is given as . (3.21) Note that we have introduced a cut-off at Λ = k e , so only modes that exit the horizon during inflation and re-enter during reheating contribute to the total number/energy density. 1+3w dk , where we have assumed that w ∈ (−1/3, 1) such that in the last step we used γ 2(w−1) 1+w 1 approximation. For the case of light DM, the dominant mode momentum k is (3.23) For typical values of parameters considered in this work, i.e. γ ∈ [10 −1 , 10 −10 ], H I 6.6×10 13 GeV, a e ≈ 10 −12 , w ∈ (−1/3, 1), and m X H I , we find that k m X for both light and heavy DM cases. In other words, since a = k /m X a 0 ≈ 1, it implies that the dominant modes (i.e. those with the momentum k ) re-entered the horizon in the far past.
Finally, we are ready to calculate the present day relic abundance of DM particles where ρ c is the critical density and T 0 refers to the present temperature. The present num- only. Therefore detection of such DM in laboratory experiments is rather unlikely. However, absolute stability of DM is not necessary as long as the lifetime of DM is larger than age of the Universe. In this case, one may allow a small mixing of vector DM with the SM photon via the dark U (1) X kinetic mixing with the SM hypercharge U (1) Y , i.e. 1 2 B µν X µν , where B µν is the field strength tensor for the SM U (1) Y gauge boson B µ . For very small values of the mixing parameter the vector DM can be stable at the time scales of the age of the Universe and there is possibility of direct detection, see [15] and references therein.

Conclusions
The aim of this work was to investigate the possibility of gravitational production of an Abelian vector dark matter due to rapidly expanding early universe. In this scenario the SM is extended by a U (1) X gauge group equipped with a stabilizing Z 2 symmetry, so that the corresponding vector boson is a DM candidate. It has been assumed that the dark sector communicates with the SM only through gravitational interactions described by the General Relativity. We have focused here on the possibility of generating vacuum expectation value of energy density for the longitudinal component of vector DM in the presence of time dependent FLRW metric in the early universe. We have shown in detail how does the canonical quantization of the vector field in this varying gravitational background imply the tachyonic enhancement of some momentum modes of the field. In all cases approximate solutions of the mode equation have been found and verified against exact numerical solutions.
We have assumed the period of inflation described effectively by de Sitter geometry, however for the following period of reheating we have adopted a generic equation of state with its parameter varying in the window −1/3 < w < 1. That way we have effectively taken into account possibilities of unknown dynamics modeled by some unspecified inflation scenario followed by an extended period of reheating with non-standard early universe cosmology. It has been shown that the spectrum of dark vectors produced that way is centered around a characteristic comoving momentum k that is determined in terms of the mass of the vector, the Hubble parameter during inflation H I , the equation of state parameter w and the efficiency of reheating γ. The ultimate result of this work was to calculate the present total vector-dark-matter abundance produced purely gravitationally. Regions in the parameter space consistent with the Planck measurement of Ω DM have been determined justifying the gravitational production as a viable mechanism for vector dark matter production for a wide range of DM masses. In particular, we found non-trivial dependance of the relic abundance on the vector DM mass m X and the Hubble parameter during inflation H I . This dependance of the relic abundance on m X and H I is different for the heavy DM H rh < m X < H I and light DM m X < H rh regimes. The results obtained in this paper are applicable within various possible models of inflation/reheating with nonstandard cosmology parametrized by corresponding equation of state.
where X µν ≡ ∂ µ X ν −∂ ν X µ denotes the field strength and m X is the vector boson mass. The background metric is in the FLRW form (2.2). The above action results in the following equations of motion, that allows to rewrite Eq. (A.8) in the desired form of an oscillator equation. We obtain the equation of motion for X L in the following form X L + ω 2 L (τ ) X L = 0, (A. 10) with the frequency (A.11) Let us here recall some basic mathematical facts about time-dependent oscillator equation. Such equations have a two-dimensional space of solutions, spanned by X (2) ± , respectively. The general solutions are given by where a ± k , b ± k are complex time-independent constants and L , X ± ≡ X (1) ± . (A.14) Using Eqs. (A.4, A.12, A.13) we get where we have promoted X L(±) to the quantum field operatorsX L(±) X L (τ, x) = d 3 k (2π) 3/2 L ( k)â k X L (τ, k)e i k· x + L ( k)â † k X * L (τ, k)e −i k· x , (A.16) if the Wronskian: is time-independent and normalized as follows To solve equations of motion for the two transversely-polarized mode (A.7) and the single longitudinally-polarized mode (A.10) we impose the Bunch-Davies initial conditions: This boundary condition together with (A.21) completely fixes the mode functions X L,± .