Completely dark photons from gravitational particle production during the inflationary era

Starting with the de Broglie-Proca Lagrangian for a massive vector field, we calculate the number density of particles resulting from gravitational particle production (GPP) during inflation, with detailed consideration to the evolution of the number density through the reheating. We find plausible scenarios for the production of dark-photon dark matter of mass in a wide range, as low as a micro-electron volt to 1014 GeV. Gravitational particle production does not depend on any coupling of the dark photon to standard-model particles.


Introduction
An abundance of observational evidence indicates that our Universe is filled with a mysterious, invisible substance that we call dark matter [1]. Assuming that the dark matter is a weakly-interacting collection of as-yet-unidentified elementary particles, the viable theory space is vast. Apart from weak constraints on the dark matter particle's mass, 10 −22 eV m 10 19 GeV where the lower limit is from the requirement that the de Broglie wavelength of the particle is less than the size of dark-matter dominated objects and the upper limit is the requirement that the particle is not a black hole, the cosmological and astrophysical data provide no solid additional information about the dark matter's other properties (e.g., spin), and little information about the dark matter's interactions, apart from the fact that it must couple extremely weakly to visible matter. In fact, it's useful to bear in mind that the data is consistent with a model of dark matter that only interacts gravitationally with visible-sector matter. But if the dark-matter particle has only gravitational interactions with visible matter, the question arises: "How was the dark matter produced in the early universe?" A natural answer is that the origin of the dark matter must be through its gravitational interactions. That is the explanation we pursue. In this work we assume that the dark matter is a massive spin-1 particle, and we study JHEP03(2021)283 JHEP03(2021)283 Figure 1. This diagram illustrates the comoving Hubble scale (aH) −1 and the comoving Compton wavelength (am) −1 of the vector field for the three reheating scenarios discussed in the text, and it highlights the comoving wavenumbers, k = 2π/λ, that are important for understanding the spectrum of gravitationally-produced spin-1 dark matter. Inflation ends at a = a e when H(a e ) = H e , and we assume m H e . In the de Sitter phase H ≈ const., before reheating H ∝ a −3/2 , and after reheating H ∝ a −2 . We illustrate three possibilities: Late Reheating: If reheating has not completed by the time when m ≈ H(a), then the spectrum is well-approximated by a power law with a single break at k −1 ≈ (m/H e ) −1/3 , which corresponds to the special mode that reenters the Hubble radius at the same time when m = H. Early Reheating: If reheating completes before m ≈ H(a) then the spectrum is a power law with two breaks. Immediate Reheating: In the limit where the duration of reheating goes to zero, the early reheating scenario is approximated by the immediate reheating scenario in which the spectrum is a power law with its break at k −1 ≈ (m/H e ) −1/2 . The total relic abundance (integral of the spectrum) is shown to be relatively insensitive to the reheating history for ultralight dark photons.
Our study is closely related to the work that's presented in ref. [44]. The authors of that article have also studied the gravitational production of spin-1 dark matter while accounting for the finite duration of reheating. We follow a similar analysis here, but with two notable differences in our assumptions. First, for our analytical calculations, we restrict our attention to models of reheating with equation of state w = 0, whereas the work in ref. [44] allows for a more general range −1/3 < w < 1. Our assumption is motivated by models of inflation with a quadratic inflaton potential near the minimum, which predict w = 0. Our results generally agree with the w = 0 case in ref. [44]; e.g., compare our table 1 with their eq. (3.31). Second, for our numerical calculations, we study a quadratic inflaton potential, V (φ) = m 2 φ φ 2 /2, and we solve the inflaton's equation of motion to determine the is derived in the discussion before eq. (3.3). The lower limit to T RH is due to the requirement that the RD universe is able to produce the neutrino background [47].
background spacetime, i.e. a(t), H(t), and R(t). By contrast, ref. [44] assumes an exact de Sitter phase of inflation followed by an immediate transition into reheating with equation of state −1/3 < w < 1. As we show in figure 6, accounting for the evolution of H(t) during inflation, as we have done here, can lead to an O (10) change in the predicted dark matter relic abundance.
The reader should also compare our work with ref. [40], in which the authors present a systematic study of gravitational particle production for vector dark matter (and also spin-1/2 fermions). The authors of ref. [40] recognize that gravitational particle production can be efficient even for particles with mass above the inflationary Hubble scale but below the inflaton mass scale, H inf < m < m φ . Our work focuses instead on light dark matter with m H inf ∼ m φ . Our analytic results generally agree with the light vector boson case of ref. [40], which is also in agreement with the earlier ref. [43].
For those interested in the final answer, the result for the contribution to the present mass density of dark matter, parameterized by Ωh 2 /0.12, is shown in table 1. In the table T RH is the reheat temperature (and T MAX RH is its maximum possible value) discussed in section 3, a RH /a e is the ratio of the scale factor at reheating to the scale factor at the end of inflation, m is the mass of the dark photon, and H e is the expansion rate of the Universe at the end of inflation.
The remainder of this article is organized as follows. We present the massive vector model in section 2 for a Minkowski spacetime background, and we extend it to an inflationary background in section 3 before reviewing the phenomenon of gravitational particle production in section 4. Our main results appear in sections 5 and 6, where we solve the vector field's mode equations -both analytically and numerically -to calculate the spectrum and relic abundance of gravitationally-produced spin-1 dark-matter particles. We summarize and conclude in section 7.

The massive vector model
We will only consider spin-1 fields with non-zero mass because, as we shall see, massless spin-1 fields (e.g., electrodynamics) are conformally coupled to gravity and will not be produced by the expansion of the universe. Since we are interested in GPP of massive vectors as a source of dark matter, our considerations will not apply to the massive spin-1

JHEP03(2021)283
particles of the standard model (W ± and Z). The vector field must transform as the 1 2 , 1 2 representation of the Lorentz group. It contains components with helicity 1 and 0.
For the analysis of massive spin-1 bosons, we start with the de Broglie-Proca action in Minkowski space [48][49][50]: Here η µν = diag(1, −1, −1, −1) is the Minkowski metric and F µν is the field strength tensor. We will see that in the massive (as in the massless) theory, A 0 is not dynamical. This Lagrangian is the unique renormalizable Lorentz-invariant Lagrangian for a massive spin-1 field. 2 Unlike the familiar electroweak theory, the de Broglie-Proca Lagrangian does not describe a gauge theory because the mass term explicitly breaks gauge invariance, i.e., invariance under the local transformation A µ (x) → A µ (x) + ∂ µ α(x). However, we can view the action of eq. (2.1) as the effective low-energy theory of a gauge theory, namely the Abelian-Higgs model with a complex scalar field Φ which obtains a vacuum expectation value v.
where A µ is a massless gauge field, after symmetry breaking and integrating out the massive scalar, the effective theory is equivalent to the de Broglie-Proca theory with the mass of the vector field m = gv. In this approach the de Broglie-Proca Lagrangian is the effective low-energy theory of an Abelian-Higgs model in the limit v → ∞, g → 0, and gv → const. We could relax the v → ∞, g → 0 limit and just assume an Abelian Higgs model where the mass of the Higgs (of course this is not the electroweak Higgs) is larger than H during inflation while the mass of the vector is of order or smaller than the expansion rate during inflation. Or perhaps the Higgs is produced during inflation and then decays. It would presumably decay to the massive vector, so there would be two sources of remnant vectors: GPP of the massive field during inflation, and production of the massive vector through Higgs decay. 3 The antisymmetric field-strength tensor in terms of the vector field A µ is given by The classical equation of motion is the so-called Proca equation Note that since ∂ ν ∂ µ F µν = 0, from the Proca equation we find the Lorenz gauge condition ∂ ν A ν = 0. This condition, usually set by gauge fixing in the massless theory, is a consequence of the equation of motion of the massive theory. From the Proca equation, the gauge field satisfies four copies of the Klein-Gordon equation for the four components of A µ : In this section we follow Weinberg [51]. 3 If the massive spin-1 dark photon arises from an Abelian Higgs model in the UV, then the theory predicts an additional spin-0 Higgs boson. We would have to assume that its mass is larger than 2m so that it is unstable and decays pairwise into dark photons. Additionally, we would have to assume that its mass is larger than O(few × m inflaton ) so that its gravitational production is suppressed.

JHEP03(2021)283
The conjugate momenta to A µ are π µ = ∂L/∂Ȧ µ = F 0µ . Unlike the massless vector case, the fact that π 0 = 0 will not be a problem because A 0 will be an auxiliary field. Unlike electrodynamics, which has two physical (transverse) degrees of freedom, for the massive theory there are three degrees of freedom, namely two transverse degrees of freedom, which will be denoted by A T , and one longitudinal degree of freedom, which will be denoted by A L . The m → 0 limit is tricky. The longitudinal mode survives in the m → 0 limit, but it is decoupled from the other degrees of freedom and behaves like a scalar degree of freedom (the Goldstone boson equivalence theorem).
In component form the action is 4 where √ −η = 1 is the determinant of the Minkowski metric. Note that A t does not have a kinetic term; it is an auxiliary field. The field equations in component form are Since A µ satisfies the Klein-Gordon equation we can again expand it as In terms of the normal modes the action (2.5) becomes where in the interest of notational simplicity we have suppressed the k label on A i and A t inside the integral. In order to solve for the temporal component of the field we rewrite eq. (2.8) as Now that A t is isolated it is clear it is nondynamical and we can solve for it: (2.10)

JHEP03(2021)283
After integrating out A t the action becomes Now it is useful to further decompose the spatial components of the vector field into transverse and longitudinal polarization modes. This is accomplished by first writing A i = A ki where A k (x 0 , k) is a complex 3-vector. Note that the mapping from 4-vector to 3vector is performed using the covariant 4-vector with a lowered index. We then decompose the 3-vector as where A T 1 k , A T 2 k , and A L k are complex mode functions for the two transverse and the single longitudinal polarization mode, and ε T 1 (k), ε T 2 (k), and ε L (k) are the polarization vectors, which satisfy Then the action can be broken into two terms, where Note also that the long-wavelength modes for which |k| → 0 behave identically for the two transverse polarizations and the longitudinal polarization. Although the kinetic term for A T b k is canonically normalized, the kinetic term for A L k is not. Therefore we define the field φ L via In terms of φ L k , the action for the longitudinal mode is

JHEP03(2021)283
After much manipulation we ended up with the action for two scalars, A T and φ L (A T has two degrees of freedom). Although in Minkowski space the action ends up being just the action for scalars, in a curved spacetime the result won't be quite so simple.
Using the Belifante-Rosenfeld stress-energy tensor and eq. (2.1) for L, we find This yields ρ = T 00 as

de Broglie-Proca in a Friedmann-Robertson-Walker background
Before proceeding we have to specify a background geometry. We will consider the action of eq. (2.1) in a particular curved space, namely the Friedmann-Robertson-Walker (FRW) spacetime. Since we are concerned with the early-universe evolution we are justified in taking the spatially-flat FRW metric ds 2 = dt 2 − a 2 (t)dx 2 . 5 In conformal time η the metric is simply ds 2 = a 2 (η) dη 2 − dx 2 . We will assume an initial inflationary epoch terminating at a = a e , followed by a matter-dominated (MD) era that ends with reheating at a = a RH . 6 We choose a to have dimension of length (hence, coordinates η and x are dimensionless). In the spatially-flat case we are free to scale a. We define a e to be the scale factor at the end of inflation. Since we can set the scale, a convenient choice is a e = H −1 e where H e is the expansion rate at the end of inflation. Thus, a e H e = 1 . (3.1) Since only dη is significant, we are free to add or subtract anything to η. A convenient choice is η = 0 at the end of inflation. Thus, −∞ < η < +∞, with η = 0 at the end of inflation.
We define the wavenumber of a Fourier mode, k, to be dimensionless. The physical wavenumber with units of length −1 is k/a. The physical wavenumber at the end of inflation is k/a e . Equating k/a e and H e : k/a e = H e gives k = 1 for the wavenumber crossing the Hubble radius at the end of inflation (since a e H e = 1).
Finally, it is useful to define dimensionless parameters (3.2) 5 We adopt the Landau-Lifshitz timelike conventions [52] for the signature of the metric (sign[η00] = +1 where ηµν is the Minkowski metric), the Riemann curvature tensor (R ρ σµν = +∂µΓ ρ νσ · · · ), and the sign of the Einstein tensor Gµν = +8πGN Tµν . To translate these conventions to other conventions, see the introductory material in Misner, Thorne, and Wheeler [53]. Our sign conventions correspond to (−, +, +) in their table. 6 When we refer to the values of quantities at reheating, we mean the values when the universe becomes radiation dominated after inflation.
The dependence of the scale factor, the expansion rate, and the scalar curvature on conformal time η. We assume that the de Sitter era is followed by a matter-dominated era until reheating, which initiates a radiation-dominated era. Dimensionless variables are defined by eqs. At the end of inflation and the beginning of the matter-dominated era, α = 1 and h = 1. At the end of the MD era and beginning of the RD era, α = α RH and h = h RH .
For analytic work we will assume an initial exact de Sitter (dS) phase, followed by an immediate transition to a Matter-Dominated (MD) phase at a = a e , followed by another immediate transition to a Radiation-Dominated (RD) phase at a = a RH . It will prove useful to collect the dependence of α, h, and R/6H 2 e on η for the dS, MD, and RD eras together in a single place: table 2.
For numerical results we will assume a chaotic inflation model. In chaotic inflation the dynamics of inflation is determined by the dynamics of a scalar field known as the inflaton. The inflaton potential is taken to be V = 1 2 m 2 φ φ 2 , where m φ is the inflaton mass. To be sure, this model is observationally challenged by precision CMB observations (see, e.g., ref. [54]), but it should serve our purposes and represent a large (but not exhaustive) class of slow-roll inflation models. The end of inflation for this model occurred when φ 2.5 × 10 18 GeV, or roughly the reduced Planck mass. The expansion rate at the end of inflation is H e m φ /2. 7 In the simple single-field model of inflation the expansion 7 We note that in the chaotic model of inflation the inflaton mass and He are approximately the same, but in general they can differ. For example, hybrid or hilltop models allow H inf m φ . If they are very different, then the exponential suppression in GPP for m > He can be avoided for H inf mspectator JHEP03(2021)283 rate during inflation is related to the amplitude of gravitational waves produced during inflation. The present limit on the gravitational wave contribution to the CMB limits H to be H 7.5 × 10 13 GeV. This is the limit on H approximately 30-60 e-folds in a before the end of inflation. Thirty e-folds in scale factor before the end of inflation in this model corresponds to about 4 times H e . Therefore, the limit on H e is approximately H e 3 × 10 14 GeV. We will display the dependence on H e .
After the end of inflation in the chaotic model the field reaches the minimum of the potential and commences oscillations about the minimum of the potential. During this period of oscillation about the minimum of the potential the amplitude of oscillations decreases due to the −3Hφ term in the equation of motion. In the oscillatory phasė ρ φ + 3Hφ 2 = 0. Since φ rapidly (compared to H) oscillates about the minimum of the potential,φ can be replaced by its average over an oscillation cycle, φ 2 cycle = ρ φ , anḋ ρ φ + 3Hρ φ = 0, exactly the behavior of a matter-dominated universe. Of course the oscillatory phase cannot continue indefinitely. The φ field must eventually decay into radiation. This can be modeled by including in the equation of motion a decay term Γ φφ . If Γ φ H e , the additional term will only be important during the oscillatory phase. Because of the Γ φ term the coherent energy in the φ oscillations are converted to light degrees of freedom (radiation) and the universe "reheats". 8 The temperature of the universe when it becomes radiation dominated is known as the reheat temperature, T RH . We will display the dependence on T RH .
Not much is known about the reheat temperature. Clearly the universe was radiation dominated during big-bang nucleosynthesis, so a reasonable lower bound on T RH might be a few MeV. In order to thermalize the neutrino background (as detected in the CMB) the reheat temperature must be greater than T RH > 4.7 MeV [47]. If all of the inflaton energy density is immediately converted to radiation at reheating, then (π 2 /30)g * RH T 4 RH = 3H 2 RH M 2 Pl . Here g * RH counts the effective number of degrees of freedom in the radiation at a temperature of T RH . We will set g * RH = 106.75, the value counting the number of effective degrees of freedom in the standard model. Since there are orders on magnitude uncertainty in H e and T RH we will not bother carrying the dependence on g * RH . For immediate reheating, H RH = H e , and T MAX RH /10 9 GeV = 8.4 × 10 5 (H e /10 12 GeV) 1/2 . This is the upper bound on T RH . Using the fact that during the matter-dominated phase H ∝ a −3/2 , then H RH /H e = (a e /a RH ) 3/2 , and we can relate T RH , H e , and a RH /a e : We will show below that the requirement for reheating to affect the final number density is that α RH < µ −2/3 . So we will have to take evolution through the radiation-dominated m inflaton [55,56]. 8 "Reheat" is somewhat of a misnomer since T RH is not the maximum temperature reached after inflation: see e.g., ref. [57].

JHEP03(2021)283
era into account in calculating the final value of the number density if Promoting the action of eq. (2.1) to a general spacetime with metric g µν (x) yields The tensor structure of the vector field admits two different forms of dimension-4 operators describing non-minimal interactions of the vector field with the gravitational field, here proportional to the two constants ξ 1 and ξ 2 . In our analysis eventually we will only consider minimal coupling (ξ 1 = ξ 2 = 0), but we will carry the nonminimal terms to serve as a reference for possible future investigations. The field strength tensor is The stress-energy tensor is The F 2 terms are the familiar stress-energy tensor for the Einstein-Maxwell theory (massless electromagnetism). The first square brackets are terms arising from the vector field's mass and nonminimal coupling to gravity. The second square brackets only contains terms from nonminimal gravitational involving derivative terms. Note that We calculate the trace to be (3.9)

JHEP03(2021)283
Note that T µ µ = 0 for a massless (m = 0) and minimally-coupled (ξ 1 = ξ 2 = 0) vector field. This calculation reveals that a massless vector with a minimal coupling to gravity is conformally coupled to gravity, so particle production will depend upon m and/or (ξ 1 , ξ 2 ). Finally, we calculate the energy density, ρ = g µ0 g ν0 T µν where g µν = (1, −a 2 , −a 2 , −a 2 ). We write where Now we specialize to the FRW geometry. Just as was done for massive spin-1 fields in Minkowski space, it is convenient to remove the auxiliary field and decompose the vector field into transverse and longitudinal mode functions.
In component form the action of (3.5) assuming the FRW metric is [cf. eq. (2.5)] where we have defined In the expression for ρ we have grouped the terms based on the number of derivatives of the field. In the limit m 2 eff,x = m 2 eff,t → m 2 and a → 1 = const., we recover the Minkowski result (2.19). To gain some intuition it is useful to consider static field configurations and to take A 0 = 0, which causes the energy density to reduce to If ξ 1 = ξ 2 = 0, for a relativistic vector field we find ρ ∝ a −4 from the gradient terms, and once the field becomes non-relativistic we have ρ ∝ a −2 from the non-gradient terms, which is notably different from the behavior of a non-relativistic scalar field for which ρ ∝ a 0 . The origin of the difference is that g µν appears in the mass term for vectors. Expanding the field in terms of mode functions (2.7), the action becomes (again Here we have performed the integrals over k and x to leave only the integral over k. Setting a = 1 and m 2 eff,x = m 2 eff,t = m 2 we recover the Minkowski result (2.8). Again, for notational simplicity we have suppressed the k label on A i and A t .
In order to solve for the temporal component of the field we rewrite eq. (3.17) as [cf. eq. (2.9)] (3.18)

JHEP03(2021)283
Now that A t is isolated it is clear it is nondynamical and we can solve for it [cf. eq. (2.10)]: Then, integrating out A t , the action becomes Using again the orthonormal set of basis vectors of eq. (2.13) and the mode functions, Before proceeding further we express S T and S L in conformal time η: Here we see that the action for the individual transverse modes is precisely the action for a scalar field with m 2 eff,x defined in eq. (3.13b) (recall for scalars m 2 eff = m 2 + (1 − 6ξ)R/6). We also see that we have to have a field redefinition to have a proper action for the longitudinal action. We also note that if we keep nonminimal terms in the action m 2 eff,t can be negative and the kinetic term could be negative, leading to a ghost-like action [5]; we will consider only minimal gravitational interactions. As to the rationale for only considering minimal gravitational interactions, we note that the addition of the nonminimal terms in eq. (3.5) breaks gauge symmetry (as does the mass term, but that might arise from the Stuckelberg trick).
To have a correct kinetic term for the longitudinal mode we define χ L k as [cf. eq. (2.16)] (3.23) To simplify notation we will drop the superscript L and the subscript k on χ L k and suppress the k subscript on κ k , with the understanding that χ represents the Fourier mode for the JHEP03(2021)283 longitudinal component, and that κ (not to be confused with κ = 8πG) is a function of k and η. For future use we note With the field redefinition (2.16) the kinetic term is where the second equality is the result of an integration by parts. This leads to an action for the longitudinal component of (3.26) To summarize, the transverse and longitudinal components are independent, with actions where we have defined the squared natural frequencies to be (3.28b) The mode functions A T b k and φ L k satisfy the mode equations The frequencies in general are rather complicated, but they simplify if we consider ξ 1 = ξ 2 = 0. With that choice m 2 eff,x = m 2 eff,t = m 2 , where m 2 is a constant, leading to (3.30b) Thus, the transverse mode behaves as a conformally-coupled scalar field (ξ = 1/6) with two degrees of freedom. In the limit am k, ω 2 T = k 2 is time-independent and the mode JHEP03(2021)283 will not be populated by expansion. That is not true for the longitudinal mode. In the limit am k, ω 2 L (η) = k 2 + a 2 R/6, and the longitudinal component appears as a massless, minimally-coupled scalar field (ξ = 0), which will be populated in expansion. In the latetime limit k am, R m 2 , and H 2 m 2 the frequency of both modes have the expected form, ω 2 = a 2 m 2 .

Gravitational particle production (GPP) during inflation
Now we turn to the phenomenon of gravitational particle production. The idea that the expansion of the universe may result in particle production goes back at least as far as a 1939 paper by Erwin Schrödinger [58]. Its modern field-theory incarnation started with the early work of Parker (see, e.g. [59]). Quantum field theory in curved spacetime has been well developed (see e.g., [23]), and in the context of inflation it has been studied with an eye towards producing dark matter; first studied assuming the spectator field was a fermion or scalar [30,60], and more recently assuming the spectator field is a massive vector [43]. In this paper we focus on the massive vector case.
The basic idea behind GPP is that unless the terms in the Lagrangian involving the field are invariant under conformal (Weyl) transformations (operationally this means that the trace of the stress-energy tensor must not vanish) a rapid expansion of the universe will "pull" particles from the vacuum to propagate as real particles.
It is convenient to calculate GPP by calculating the Bogoliubov coefficient relating the early-time and late-time vacua. which implies that the late-time observer detects particles: The mode equations (3.29) are solved subject to initial conditions (4.2) to obtain χ k (η). The modulus of the second Bogoliubov coefficient is extracted from the solution to the mode equations: where φ k stands for either A T b k or χ and ω k is the corresponding value of ω T or ω L . The factor of φ k ∂ η φ * k − φ * k ∂ η φ k = i as demanded by the commutation relations. We will define JHEP03(2021)283 the spectrum of the mode function, n k , and the comoving number density of particles, na 3 , as Now for initial conditions. For GPP in the inflationary era the early-time limit (η → −∞) corresponds to a → 0 and a 2 R → 0, which implies ω 2 k → k 2 . As η → −∞ the modes are deep within the Hubble radius and their mode equation is approximately that of Minkowski space. Thus, the natural initial condition, the so-called Bunch-Davies initial condition, on the mode functions (4.6) (The factor of 1/ √ 2k ensures the commutation relations are properly normalized.) Gravitational particle production of the transverse mode is exactly the same as the well-studied case of GPP of conformally-coupled scalars. Conformal symmetry is exact in the limit m → 0, so the result for the transverse component must vanish as m → 0. That is not true for the longitudinal mode.
An example of the evolution in a of |χ k | 2 for the longitudinal component of a vector field and for a minimally-coupled scalar (which is identical to the transverse component of the vector field) is illustrated in figure 2 for a particular choice of m/H e = 10 −2 and k = 10 −3 . Note the region starting at a = k/m where |χ| 2 is constant for the longitudinal component of a vector, but grows as a 2 for a minimally-coupled scalar.
A numerical solution for the massive-vector spectrum with m/H e = 0.1 is shown in figure 3. Notice that although the frequency of the longitudinal mode resembles a minimally-coupled scalar as modes cross the Hubble radius during inflation, the spectrum at small k does not resemble the scalar spectrum, which grows at small k. The integration of the mode functions is also shown in figure 3. For this figure we assumed that the mode evolves to become nonrelativistic (k/a < m) with H/a < m. In that region |χ k | 2 oscillates and decays as a −1 , i.e., it behaves as nonrelativistic matter.
From the numerical results we see the expected result that for m > H e the mode function is exponentially damped as e −πm/He , 9 and that the modes are exponentially damped for k > 1.

Analytic approximation to the comoving number density
The goal in this section is to obtain an analytic approximation for na 3 for the longitudinal component where we consider the possibility that reheating to a radiation-dominated phase occurs before the modes have reached the point where a particle description is appropriate.  In terms of dimensionless quantities α, µ, h, the frequency for the longitudinal component, eq. (3.30b), becomes terms in ω 2 L in four regions of α and µ, denoted by I dS -IV dS : In MD, h = α −3/2 , and R/6H 2 e = − 1 2 α −3 , and again there are four possible dominant terms in ω 2 L ; they are in regions I MD -IV MD : 10 10 When considering transitions between different regions in the MD era we will be cavalier about numerical factors of order unity. We will use " " to indicate equations where we have dropped order unity numerical factors.

JHEP03(2021)283
In RD, h = α 1/2 RH /α 2 , R/6H 2 e = 0, and the dominant terms in ω 2 L are 11 We are interested in the scaling of |χ k | 2 with α when various terms dominate ω 2 L . In order to solve the wave equation for |χ k | 2 in various regions, we have to convert the α-dependence of ω 2 L to a η-dependence using table 2. The wave equations for the various regions are given in table 3. Let's take each region in turn:

I dS : The wave equation and solution in this relativistic sub-Hubble region is
where we have used the fact that the Bunch-Davies boundary condition yields c 2 = 0.

II dS : The wave equation and solution in this relativistic super-Hubble region is
where we have only kept the growing mode (η −1 = α). Together the k 2 and α 2 R terms give a mode equation that resembles the Mukhanov-Sasaki evolution equation for curvature perturbations [26][27][28].
3. III dS : In this nonrelativistic super Hubble radius region Since k < αµ in this region implies kη < µ, and the magnitude of the argument of the exponentials is small and expansion is justified. Here we have kept the growing mode.

JHEP03(2021)283
where we have only taken the growing mode and have used µ < 1. This region is also relativistic super-Hubble.
The results for dS are also summarized in table 3. Various regions and the scaling with α for dS are indicated in figure 4. Also indicated in the figure is the physical significance of various regions: relativistic for k/a > m, nonrelativistic for k/a < m, super-Hubbleradius for k/a < H, and sub-Hubble radius for k/a > H. In dS, in the relativistic-sub-Hubble region |χ k | 2 ∝ α 0 , in the relativistic super-Hubble region |χ k | 2 ∝ α 2 , and in the nonrelativistic region |χ k | 2 ∝ α 0 . For all regions we are assuming m < H.

I MD : In this relativistic sub-Hubble region the wave equation and solution is
Here, some explanation is required. Since k > α −1/2 in this region, the argument of the trigonometric functions, kη 2k √ α, is much larger than unity, and χ k will oscillate with constant amplitude, hence |χ k | 2 ∝ α 0 .
6. II MD : The wave equation and solution in II MD (relativistic super-Hubble) is where we have only kept the growing mode (η 2 = α in MD).

III MD : For this penultimate region (nonrelativistic H > m) in MD, the wave equation and solution is
where here x = 5/2 (k/αµ). In this region k < αµ, and the expansion of the solution for small x yields χ k = c 1 2/x + c 2 . The mode enters this region at the end of inflation with χ k ∝ α 0 . This boundary condition implies c 1 = 0 and |χ k | 2 ∝ α 0 .

JHEP03(2021)283
8. IV MD : Here, in the final region (nonrelativistic H < m) the scaling of the solution to the wave equation with α depends upon whether α is larger or smaller than µ −2/3 .
The results for MD are also summarized in table 3. Various regions and the scaling with α for MD are indicated on the top panel of figure 5. Also indicated in the figure is the physical significance of various region: relativistic for k/a > m, nonrelativistic for k/a < m, super-Hubble-radius for k/a < H, and sub-Hubble radius for k/a > H. In words: in the relativistic super-Hubble region |χ k | 2 ∝ α 2 ; in the relativistic sub-Hubble region |χ k | 2 ∝ α 0 ; in the nonrelativistic region |χ k | 2 ∝ α 0 if H > m and |χ k | 2 ∝ α −1 if H < m. JHEP03(2021)283 Figure 5. Regions in the α-k plane for evolution in the MD phase (upper) and the RD phase (lower). In different regions we indicate the scalings of |χ k | 2 with α. Note that k = µ 1/2 α 1/4 RH will be larger than unity (hence off the graph) if α RH > µ −2 . In the RD phase, α RH will be smaller (larger) than α 1/4 RH µ −1/2 if α RH is smaller (larger) than µ −2/3 . For both MD and RD, the line denoted k = αµ is the line delineating the nonrelativistic (k < am) and the relativistic (k > am) regions. In MD the line k = α −1/2 and in RD the line k = α −1 α 1/2 RH is the line denoting k = aH; above the lines the mode is sub-Hubble-radius (k > aH), while below the line the mode is super-Hubble-radius (k < aH). The values of α = µ −2/3 for MD and α = α Note from table 2 that in RD kη kα −1/2 RH α. Thus, for α < α 1/2 RH k −1 the argument of the trigonometric functions is much less than unity and upon expansion yields for the growing mode χ k c 1 α; hence, |χ k | 2 ∝ α 2 . If α > α 1/2 RH k −1 the solution will be an oscillation in α with frequency kα −1/2 RH and constant amplitude. 10. III RD : In this relativistic super-Hubble region we find (5.17) By way of explanation, the argument of the trigonometric functions is approximately k/α 1/2 RH µα. In IV RD , k < αµ and α RH > 1, so the argument of the trigonometric functions are small, and |χ k | ∝ η ∝ α 1 .
The results for RD are also summarized in table 3. Various regions and the scaling with α for RD are indicated on the bottom panel of figure 5. Also indicated in the figure is the physical signifiance of various region: relativistic for k/a > m, nonrelativistic for k/a < m, super-Hubble-radius for k/a < H, and sub-Hubble radius for k/a > H. In words: in the relativistic super-Hubble region |χ k | 2 ∝ α 2 ; in the relativistic sub-Hubble region |χ k | 2 ∝ α 0 ; in the nonrelativistic region The evolution of the modes in dS, MD, and RD are the same in the various physical regions. There are however some differences between MD and RD. Firstly An example of the evolution of |χ k | 2 with a is shown in figure 2, where it is compared to the evolution of a minimal scalar for the same values of k and µ. The important difference between minimal-scalar and vector evolution is in the region k/m < a H 1/3 e m 2/3 . In this region |χ| 2 grows as a 2 for a minimal scalar and is constant for a vector. Thus, the final result will be a factor of (k/m)/(H

Evolution of the modes
Now we will start with a k-mode deep in the de Sitter era in the Bunch-Davies vacuum and follow |χ k | 2 until it reaches the nonrelativistic region with m > H. In this region |χ k | 2 oscillates with amplitude decreasing as α −1 . In the region the evolution is adiabatic and one can sensibly defining a number density of particles resulting from GPP. This will be the asymptotic behavior of |χ k | 2 , and thereafter α|χ k | 2 will remain constant.

de Sitter evolution
We first consider the evolution of |χ k | 2 in the de Sitter era. As α → 0, we will assume Bunch-Davies vacuum and take as initial conditions Starting with those initial conditions we can follow the evolution of |χ k | 2 easily by referring to figure 4. Consider two cases for the evolution of |χ k | 2 in the dS era: 1. 1 > k > µ. The mode begins in the Bunch-Davies vacuum and remains constant until α = k/ √ 2. Then it grows as α 2 in II dS until the end of inflation. So at α = 1, 2. µ > k > 0: Again, the mode begins in the Bunch-Davies vacuum and remains constant until α = k/ √ 2. Then it grows as α 2 until it crosses α = k/µ, after which it remains constant until the end of inflation. At α = 1, In conclusion, the mode amplitudes (squared) at the end of inflation are This result agrees with calculations by other authors such as refs. [40,44].

JHEP03(2021)283
One caveat is that we have assumed R is constant in dS. In a slow-roll model typically R grows as a logarithm in α as α → 0. We will discuss a correction for this later.

Matter-dominated evolution
Now consider the evolution of |χ k | 2 in the matter-dominated era until the evolution to the nonrelativistic, sub-Hubble-radius region. This amounts to following the evolution past α = µ −2/3 (see figure 5) assuming the mode reaches the nonrelativistic, H < m region before reheating. We will describe this possibility as the late-reheating case. We will denote this asymptotic value of |χ k | 2 as |χ k (α → ∞)| 2 .
There will be three cases, depending on the value of k. Again, with the help of the upper panel in figure 5 we can follow the evolution through the MD era. Our goal is to find the value of |χ k (∞)| 2 , which will be used to calculate n k .
2. µ 1/3 > k > µ: In this range of k the mode again enters MD in the relativistic-super-Hubble region and evolves as α 2 . Then, when α = k/µ it enters the nonrelativistic, H > m region, after which it remains constant until it crosses α = µ 2/3 , then it damps as α −1 . Gluing together the pieces of evolution, 3. µ > k > 0: For this final case the mode enters MD through the nonrelativistic H > m region and remains constant until it enters crosses into the H < m region at α = µ −2/3 when it begins damped oscillations. Thus, Note than now we have used |χ k (α = 1)| 2 = (kµ 2 ) −1 for µ > k as in eq. At low k they differ by a factor of ∼ 10, which is explained in the text.
The conclusion is that for α RH > µ −2/3 , peaked around k = µ −1/3 . We can find the total number density by integrating eq. (5.26): Modes of higher k are damped and won't contribute significantly to na 3 . Also, recall that we are only considering µ < 1 since higher-mass modes are also damped.
In the infrared the k-dependence of the analytic results is k 2 , while the numerical results for the chaotic model are better fit by a dependence of k 1.8 . This discrepancy is due to the fact that the scalar curvature R is not constant in the chaotic model, but increases as one goes further back in inflation. The numerical result k 1.8 is for the chaotic model, and a different inflation model may give a different scaling. The value of the integrated spectrum will not depend much on the exact infrared behavior so long as n k → 0 as k → 0. This is problematic for a minimally-coupled scalar, but no problem for the vector, which has a blue spectrum. The infrared dependence will have a larger effect on the isocurvature component as well as nongaussianities. So long as the spectrum decreases in the infrared faster than k, isocurvature issues should not arise [43].
Since the scaling of n k with k in the infrared leads to a convergent result for na 3 = n k d ln k, the infrared behavior does not much affect the total number density; rather, the total number density depends on the value of n k around the peak at k µ 1/3 . From eq. (5.26), the peak value scales as n k (k = µ 1/3 ) ∝ µ −1 = H e /m. This implies that the contribution to the mass density, proportional to m n a 3 , is roughly independent of m! This will be discussed in the next section.

Radiation-dominated evolution
Now consider the effects of reheating, which is important if reheating occurs before the mode reaches the nonrelativistic, sub-Hubble-radius region. We will call this the earlyreheating case. In the early reheating case we must consider RD evolution.
The evolution of |χ k | 2 in the RD era is shown in figure 5. It is useful to refer to the figure when discussing the evolution through reheating. First, we establish a hierarchy of inequalities for α RH < µ −2/3 : (5.28) We will again study the evolution for various ranges of k.
RH : In this range the mode enters the MD region as relativistic, super-Hubble and evolves as α 2 until it crosses into the relativistic sub-Hubble region at α = k −2 . Then it evolves as a constant, reheating occurs in this region and the mode continues to evolve as a constant until α = k/µ when it enters the nonrelativistic H < m region and thereafter damps as α −1 . The final value of n k will be
2. α −1/2 RH > k > µ 1/3 : In this region the evolution in MD begins as above, but reheating at α = α RH occurs in the relativistic super-Hubble region before the mode crosses α = k −2 . The mode then continues to grow as α 2 in the relativistic super-Hubble region of RD until α = α 1/2 RH /k. Then it evolves through the relativistic sub-Hubble region as a constant until α = k/µ and damped oscillations commence. This results in RH µ 1/2 , k will satisfy k > α RH µ. This implies reheating will occur while the mode is in the MD relativistic super-Hubble region before it crosses k = αµ. After reheating the mode will continue to grow as α 2 in the RD relativistic super-Hubble region until α = α 1/2 RH /k. It then remains constant until α = k/µ and begins damped oscillations. Thus, the final result will be The mode enters MD in the relativistic super-Hubble region scaling as α 2 as previously. It reheats before becoming nonrelativistic and continues to evolve in RD as α 2 until α = k/µ when it enters the nonrelativistic H > m region as remains constant until it crosses α = α RH µ −1/2 and starts damped oscillations. This leads to the result 5. α RH µ > k > µ: Now the mode scales as α 2 until it becomes nonrelativistic in MD at α = k/µ. Then it is constant in the nonrelativistic H > m region before and after reheating until it becomes nonrelativistic and commences damped oscillation. Therefore,

JHEP03(2021)283
Immediate Reheating Early Reheating Late Reheating equations. Notice that the analytic approximations have underestimated the spectrum by a factor of ∼ 10 at low k. This can be understood from the evolution of the Hubble parameter during inflation. For the numerical work, we assume a chaotic model of inflation with a quadratic inflaton potential. In this model the Hubble parameter decreases by a factor of ∼ 10 between the time of CMB mode generation and the end of inflation. Since the modes with smaller k leave the horizon earlier, they probe the larger H > H e , which leads to a larger n k relative to the analytic approximations that assume H = H e throughout inflation. Nevertheless, we see also from figure 6 that the analytic approximation works well for the modes where the spectrum is peaked, which means that the total abundance na 3 can be calculated reliably from the analytic approximations while only introducing an O(1) error. For other models of inflation in which the inflation potential is shallower and the H inf /H e is not much larger than 1, such as the α-attractor class of models [61] including Starobinsky's R 2 inflation [62], we expect that our analytic treatment will provide an even better approximation of the spectrum.

JHEP03(2021)283
dark matter, but it is not applicable when the dark photon mass becomes larger, m (1 GeV)(T RH /10 9 GeV) 2 . Here we generalize and extend the analysis by allowing for a finite duration of reheating, which is assumed to be a matter-dominated phase; see also refs. [40,44] that present closely related analyses. We calculate the vector field's mode functions during inflation and reheating both numerically (assuming a quadratic inflaton potential V ∝ φ 2 ) and analytically, finding excellent agreement between these two approaches. For the analytic calculation, we systematically decompose the mode equations into various regimes, depending on which term dominates in the dispersion relation, ω 2 k (η). This approach has a broad applicability, beyond simply the spin-1 dark matter calculation that we have performed here. As a result, we find that the finite duration of reheating causes the spectrum of gravitationally-produced spin-1 particles to develop two breaks, associated with the scales that reenter the Hubble radius at the time when reheating ends and at the time when m = H; these results are summarized in figure 6. Assuming that the spin-1 particles are stable and their comoving number density is conserved until today, we also calculate their relic abundance, which is shown in figure 7. For example, if H e ∼ 10 14 GeV then the observed dark matter relic abundance is obtained if m ∼ 10 −6 eV and 50 GeV T RH 10 16 GeV or if T RH ∼ 50 GeV and 10 −6 eV m 10 14 GeV. To avoid producing too much dark matter, the parameters H e , T RH , and m are constrained, as shown in figure 8.
In this work we have focused on understanding the gravitational production of vector dark matter during inflation and reheating. If this dark-matter candidate also has non-gravitational interactions, which simply did not play a role in its production, then a variety of observational probes become available, including direct detection in the lab. On the other hand, if the dark matter only interacts with itself and visible matter through gravity, then observational prospects are clearly more challenging, but nevertheless several detection channels could be available. Terrestrial probes, such as gravitational direct detection [64], are most sensitive to larger dark photon masses; although, even for masses as large as m ∼ H e ∼ 10 14 GeV, this signal would be very challenging to see. Cosmological probes of spectator fields include isocurvature (between the dark matter and curvature perturbations) and non-Gaussianity (of the curvature perturbations). Since the dark matter power spectrum is blue-tilted (falling toward smaller k) the isocurvature on CMB scales is predicted to be negligibly small [43]. On the other hand, in the quasi-single-field regime (m ∼ H inf ) the vector spectator may induce a detectable non-Gaussianity in the curvature perturbations [65,66] if it couples directly to the inflaton field. Finally the blue-titled spectrum enhances the small-scale power in the dark matter perturbations, which may lead to the formation of primordial black holes [67] and provide additional astrophysical probes of this scenario.