The role of vectors in reheating

We explore various aspects concerning the role of vector bosons during the reheating process. Generally, reheating occurs during the period of oscillations of the inflaton condensate and the evolution of the radiation bath depends on the inflaton equation of state. For oscillations about a quadratic minimum, the equation of state parameter, w = p/ρ = 0, and the evolution of the temperature, T(a) with respect to the scale factor is independent of the spin of the inflaton decay products. However, for cases when w > 0, there is a dependence on the spin, and here we consider the evolution when the inflaton decays or scatters to vector bosons. We also investigate the gravitational production of vector bosons as potential dark matter candidates. Gravitational production predominantly occurs through the longitudinal mode. We compare these results to the gravitational production of scalars.


Introduction
While the inflationary paradigm [1] is well-embedded into standard Big Bang cosmology, it lies beyond the Standard Model (SM) framework.Consequently, there is considerable model dependence regarding how the inflationary sector couples to the Standard Model-couplings that are necessary for establishing an early period of radiation domination.It is often assumed that the inflaton must decay or scatter to reheat the universe [2,3].This transfer of energy typically occurs after the period of accelerated expansion has ended, and the inflaton begins to oscillate around its potential minimum.The thermalization of these decay or scattering products is responsible for generating the thermal background that eventually comes to dominate the energy density.The process of establishing this thermal background is not instantaneous [4][5][6], and the reheating temperature, along with the scaling of the radiation density during reheating, depends on the spin of the final state particle and the shape of the potential near the minimum that governs inflaton oscillations [7][8][9].
Indeed, when the inflaton potential around its minimum is dominated by a quadratic term, V (ϕ) ∝ ϕ2 , the energy density of the inflaton, ρ ϕ , scales as ρ ϕ ∝ a −3 , where a represents the cosmological scale factor.As oscillations begin, and assuming that the inflaton decays, a bath of relativistic particles is produced which subsequently thermalizes, 1 and the energy density of radiation, ρ R , scales as ρ R ∝ a −3/2 independent of the spin of the final state products [8].This difference in scaling allows for the process of reheating to occur, defined as the moment when ρ ϕ (a RH ) = ρ R (a RH ).Furthermore, if the inflaton primarily decays to vector bosons, this also remains true.However, when the potential minimum is approximated by λϕ k for k > 2, the situation changes.The evolution of temperature with respect to the scale factor and the temperature at reheating is highly dependent on whether the decay products are bosons or fermions.In this work, we find that decay to vectors results in a temperature scaling T ∝ a −1 , which holds regardless of the value of k when k ≥ 4, and is distinctly different from decays to fermions or scalars [8].Moreover, for k = 2, the scattering of inflaton condensate to scalars cannot reheat the Universe since the energy density of radiation ρ R ∝ a −4 redshifts faster than ρ ϕ ∝ a −3 [8].For k ≥ 4, the evolution is different, and in this case, two-to-two scatterings of inflatons to scalars can also contribute to reheating.We demonstrate that the same is true for scatterings to vectors as well.
One of our motivations for considering vector decay products stems from models of inflation that are consistent with no-scale supergravity [18].As the natural framework for a low-energy field theory derived from string theory [19], no-scale supergravity serves as an excellent starting point for the construction of a theory of inflation [20].In fact, it is straightforward to construct models of inflation [20][21][22][23][24][25][26][27] that lead to Starobinsky inflation [28].However, if there are no direct superpotential couplings between the inflaton and Standard Model fields, there are no decay channels directly to matter scalars or fermions [25,29]. 2 On the other hand, the gauge kinetic function must be fielddependent if gaugino masses are generated at the tree level when supersymmetry is broken.If the gauge kinetic function depends on the inflaton, it naturally follows that the inflaton could decay into gauge bosons (and gauginos) [20,25,29,[32][33][34].In many models of inflation derived from no-scale supergravity this is the dominant mode for inflaton decay and hence decays to vectors will play a significant role in the reheating process.In this work, we will derive the evolution of the radiation bath when vectors are the primary products of either inflaton decay or scattering.
In addition to the role of vectors in the reheating process, massive vectors may also be produced through scattering, involving either a single graviton exchange or a gravitational portal [35][36][37][38][39][40][41][42][43][44].Related mechanisms for the gravitational production of particles during reheating were also considered in .In this work, to compute the production of massive spin-1 bosons, we adopt the methods previously used for the production of spin-0 and spin- 1  2 particles [36,39], as well as those for spin- 3  2 particles [44].For related work on the gravitational production of vectors, see [38,47,50,53].Vectors can be produced gravitationally either through the scatterings of SM particles in the newly created thermal bath or directly through inflaton scatterings.We will demonstrate that the latter is by far dominant.Furthermore, we will show that vector bosons produced at the earliest stage of reheating can account for the dark matter if the product of the vector mass and reheating temperature, m A T RH ≃ 10 17 GeV 2 , when k = 2.For higher values of k, we find significantly stronger constraints on the vector mass for a given reheating temperature.We compare these results to the production of scalars through gravitational interactions.
The paper is organized as follows: In Section 2, we explore the possibility that reheating occurs primarily through the decay of gauge bosons.We compute the evolution of the radiation bath, which originates in the form of vectors, and compute the reheating temperature in terms of k, the parameter that determines the shape of the potential about its minimum.In Section 3, we examine the gravitational production of vectors, considering both scattering within the thermal bath and direct production from the inflaton condensate.Additionally, we discuss two mechanisms that generate vector masses while preserving gauge invariance, namely, the Stueckelberg and Higgs mechanisms.We summarize our results in Section 4, and provide some additional details in Appendices A, B, and C.

The Inflaton Potential
Reheating is an essential component of any inflationary model, necessitating a coupling between the inflaton sector and the SM.Although gravitational interactions always mediate interactions between the inflaton and the SM, it was shown in [43,70] that minimal gravitational couplings alone are insufficient for achieving a radiation-dominated universe with a sufficiently high reheating temperature.Consequently, some form of inflaton decay, or alternatively, an inflaton scattering mechanism, is necessary.In this section, we explore the possibility for non-gravitational decays and scatterings to massive gauge bosons.We first discuss the motivation for this possibility in the context of no-scale supergravity and subsequently examine the effects of an inflaton-vector coupling (regardless of its origin) on the evolution of the radiation energy density and reheating temperature.As a phenomenologically acceptable example of an inflaton potential V (ϕ), we begin by considering the Starobinsky model of inflation [28] V where m ϕ denotes the inflaton mass, the sole independent parameter in the Starobinsky model.Here, the reduced Planck mass is defined as M P = 1/ √ 8πG N .The amplitude of the cosmic microwave background (CMB) anisotropies fixes the inflaton mass to m ϕ ≃ 3 × 10 13 GeV.This scalar potential ensures sufficient inflation, and the resulting slow-roll parameters ϵ and η will give a tilt to the CMB anisotropy spectrum.Additionally, it predicts a scalar tilt n s and the tensor-to-scalar ratio r [21,71], both of which are in good agreement with current Planck observations [72].Importantly, this potential can be derived from an ultraviolet (UV) perspective in the context of no-scale supergravity.

Motivations from No-Scale Supergravity
The Starobinsky model can be straightforwardly derived from no-scale supergravity.We define a Kähler potential of the following form (in Planck units) where T is associated with the volume modulus, φ represents the non-canonical inflaton field, 3 and the X i are associated with matter fields.The superpotential in the inflationary sector is given by (also expressed in Planck units) leads to the Starobinsky potential once the volume modulus, T , is fixed (or stabilized) [21], and the inflaton field φ is canonically normalized.Although this superpotential for the Starobinsky potential is not unique [22], various superpotentials that yield the same scalar potential can be related by the underlying SU(2,1)/SU(2)×U(1) no-scale symmetry [23].
If there are no superpotential couplings between φ and other SM fields, inflaton decay becomes impossible [29].As previously mentioned, in this case (k = 2 for the Starobinsky potential), scatterings do not lead to reheating.An efficient way to introduce a coupling of the inflaton to SM fields is through the gauge kinetic function [25,29,32,34].In supergravity, the kinetic terms for gauge fields can be written as where represent the field strength and its dual, respectively.The inflaton couplings to gauge bosons can arise from the gauge kinetic term if the gauge kinetic function f αβ depends on the inflaton.
In a supergravity model, an inflaton-dependent modification of the gauge kinetic function leads not only to a coupling of the inflaton to gauge bosons but also to gauginos, which will also be produced through inflaton decay.The derivative coupling between the gauge boson and the inflaton is given by (2.4) and the relevant gaugino-inflaton coupling arises from Here G ≡ K + ln W + ln W * is the Kähler function, with G i ≡ ∂G/∂ϕ i and G j ≡ ∂G/∂ϕ * j , and λ α denotes the gaugino fields.By an appropriate choice of the gauge kinetic function a coupling proportional to m ϕ ϕλλ emerges, with the dependence on the inflaton mass originating from the superpotential (2.3).If we do not associate supersymmetry breaking (at low energies) with the inflationary sector, we can expect that at the minimum, ϕ min (ϕ min = 0 for the Starobinsky potential), h(ϕ min ) = 1, and we restore canonical kinetic terms for the gauge bosons and gauginos. 4or the moment, we will focus on the coupling of the inflaton to gauge bosons.Later in this section, we will also explore the implications of decays to gauginos.While we do not specify the function h(ϕ), it is clear that any linear or quadratic term in ϕ, naturally provides inflaton decay or annihilation channels, respectively, through derivative couplings to the vector bosons.To maintain generality, in the following subsection, we will limit our discussion to couplings of the form5 : where we normalize the coupling constants using the reduced Planck mass, so that g, g, κ, κ are dimensionless.Note that the coupling of the inflaton to gauge bosons may also induce the nonperturbative production of gauge bosons during inflation, i.e., superhorizon modes, which is on the other hand strongly dependent on the form of h(ϕ) for ϕ ≫ M P .We here assume that the couplings (2.7) are valid only after the end of inflation and will focus on the gauge boson production during reheating.

Reheating Through Decays to Vector Bosons
The most direct reheating process involves the decay of the inflaton field.However, computing particle production through the decay of a scalar condensate such as the inflaton can differ from calculating the decay width of a quantum field.For oscillations about a quadratic potential, the results are identical.In more general cases, when the potential can be expanded around its minimum as the derived rates depend on an effective coupling, averaged over a single oscillation, which, in turn, depends on the shape of the potential and hence the parameter k.In general, these effective couplings must be computed numerically [8,74,75], and details of this computation are given in Appendix A.
The Starobinsky potential, when expanded about the minimum, gives k = 2 and does not easily generalize to larger values of k.However, the related α-attractor T-models [76] are also phenomenologically viable [71].The scalar potential for these models can be expressed as which reduces to Eq. (2.8) when expanded about the minimum at ϕ = 0. We note that this class of potentials can also be generated from no-scale supergravity models by choosing the following superpotential form [7] 3(k + 6) . (2.10) As in the Starobinsky model, the lone parameter λ is determined from the normalization of the CMB anisotropies [72].The normalization of the potential for different values of k is given by [8] λ ≃ where N * is the number of e-folds of inflation and A S * ≃ 2.1 × 10 −9 is the amplitude of the curvature power spectrum.For this potential with k = 2, we find λ = 2.05 × 10 −11 and m ϕ = √ 2λM P = 1.5 × 10 13 GeV for N * = 55 e-folds.
From Eq. (2.7), the inflaton decay rate to vector bosons can be computed as where α 2 = (g 2 eff +g 2 eff )/(64π) and the effective couplings g eff and geff arising from the Lagrangian (2.7) are given by Eqs.(A.5) and (A.9), respectively.The Lagrangian coupling g (g) is defined by Eq. (A.5) (Eq.(2.8)) in Appendix A. Note that the couplings to F F and F F do not interfere due to the antisymmetry of the Levi-Civita tensor.The inflaton mass is given by For the last equality, we use ρ ϕ = V (ϕ 0 ), given by Eq. (2.8), where ϕ 0 is the time-dependent amplitude of the oscillation of ϕ, as discussed in more detail in Appendix B.
Following the procedure described in [8], it is useful to rewrite the decay rate in terms of ρ ϕ so that and The dependence of the decay rate on the energy density in this case differs from that of decays to fermions or scalars.In the case of decays to fermions, ϕ → f f , the power is given by l = 1 2 − 1 k , whereas for decays to scalar bosons ϕ → bb, l = 1 k − 1 2 [8].As expected, they follow the redshift behavior Γ ϕ ∝ m ϕ for the fermion final state, and Γ ϕ ∝ 1/m ϕ for boson final states.In the case we are considering here, the form of the couplings (2.7) imposed by gauge invariance implies a width which undergoes a much more pronounced redshift, with Γ ϕ ∝ m 3 ϕ .This dependence is very important, because it indicates that the behaviour of vector reheating is closer to that of fermions (a width that decreases with time) than to scalar boson decay (a width that increases with time).Furthermore, since the decay rate decreases over time, the energy transfer is most effective at the beginning of the process.
The decay rate determines the evolution of the radiation density through the coupled equations of motion: (2.17) where H = ȧ a is the Hubble parameter, a is the cosmological scale factor, and the equation of state parameter w ϕ = P ϕ ρ ϕ is given by Using d dt = aH d da , we can integrate Eq. (2.17) to obtain valid when γ ϕ ≪ H.Here ρ end = ρ ϕ (a end ), where a end is the value of the scale factor when ϕ(a end ) = ϕ end and ϕ end corresponds to the value of ϕ when the inflationary expansion ends, namely when ä = 0 which implies that ρ end = 3 2 V (ϕ end ).An approximate solution for the inflaton field value at the end of inflation as a function of k is given by [8] For k = 2, we find ρ 1/4 end = 5.2 × 10 15 GeV, and it slightly decreases with increasing k [8,39].Using the solution for ρ ϕ , we integrate Eq. (2.18) and obtain so at late times, when a ≫ a end , and for 8 + k − 6kl = 26 − 8k > 0, we obtain where g ρ = 915/4 is the number of relativistic degrees of freedom which we take from the minimal supersymmetric standard model (MSSM).When k = 2 and l = 0, T ∝ a − 3 8 which is the same behavior as other decay channels, ϕ → f f and ϕ → bb because the width Γ ϕ is constant for a quadratic potential.However, for k ≥ 4, we have 8 + k − 6kl < 0, and the radiation density in Eq. (2.22) is dominated by the second term, and the temperature evolves as ρ R ∝ a −4 .
This result illustrates what we discussed above: reheating primarily occurs at the very beginning of the process, subsequently, the density then follows the classical isentropic redshift law, where the decay rate ∝ ρ ϕ .In addition, the inflaton energy density from Eq. (2.20) scales ρ ϕ ∝ a − 6k k+2 , which gives a −3 , a −4 , a − 9 2 for k = 2, 4, 6, respectively.Consequently, reheating is not possible for k = 4 since ρ ϕ and ρ R scale in the same way.However, reheating is possible for k = 2 and k > 4. Nonetheless, as we discuss below, k > 10 is necessary for reheating temperatures higher than T RH > 100 GeV.This behavior can be seen from the analytic expressions for the reheating temperature.Since we define the reheat temperature when we can determine a RH /a end to find ρ RH = ρ R (a RH ), which for 8 + k − 6kl > 0 gives For 8 + k − 6kl < 0 and k > 4, which is clearly seen to be proportional to ρ end for large k.
In the absence of a term linear in ϕ in h(ϕ), a quadratic term may also produce vector bosons through annihilation, ϕϕ → A µ A µ , driven by the last two terms of Eq. (2.7).These two couplings will give rise to the same dependence of the rate on ρ ϕ and m ϕ .Defining β 2 ≡ (κ 2 eff + κ2 eff )/(4π), where κ eff and κeff are given as a function of κ and κ by Eqs.(A.13) and (A.17), respectively, we have from which we can extract (2.30) As 8 + k − 6kl = 14 − 8k is always negative when k ≥ 2, the second term in the brackets in Eq. (2.22) dominates and the temperature decreases as T ∝ a −1 at late times.This is because Γ ϕ ∝ a − 9k−6 k+2 whereas H ∝ a − 3k k+2 .For k ≥ 2, the scattering rate decreases faster than the Hubble rate, and reheating occurs at the very end of inflation.Very quickly, the rate of energy transfer from the condensate to the thermal bath is unable to compensate for the expansion rate.However, as ρ ϕ ∝ a − 6k k+2 and ρ R ∝ a −4 , reheating can only take place if k ≥ 6.The reheating temperature in this case is given by This situation is very different from scalar boson scattering [8].In this case, the radiation density evolves as ρ R ∝ a − 18 k+2 , which means that reheating is impossible for k = 2, but feasible for k ≥ 4.Moreover, in the case of scalar boson scattering, for k ≥ 4 reheating is most efficient at the end of the process, not at the beginning.The behavior of the temperature during reheating is given in Table 1 for decays into fermions, scalars, and vectors as well as scattering to scalars and vectors.
Table 1: Dependence of the temperature T as function of the scale factor a for the different cases that we analyze in this work.The 'generic' result assumes the validity of Eq. (2.23).This is not true when 8 + k − 6kl < 0, in which case the scaling is simply T ∝ a −1 .
We show in Figs.1-3 the evolution of ρ ϕ and ρ R for k = 2, 4, and 6 in the case of inflaton decay and scattering for α = 10 −2 (α = 1 for k = 6) and β = 10 −2 , respectively, in the left panels.The temperature of the radiation bath is shown in the right panels and the reheating temperature is indicated by a star.For k = 2, shown in Fig. 1, we see that reheating through the decay is possible, and quite efficient, giving a reheating temperature of T RH ∼ 10 9 GeV for the adopted value of α.On the other hand, as we predicted, reheating is not possible via scattering in this case, since ρ R ∝ a −4 dilutes more rapidly than ρ ϕ ∝ a −3 .For k = 4, Fig. 2 shows that neither inflaton decay nor scattering to vectors is efficient enough to achieve reheating.Indeed, in both cases, the inflaton density and the energy density of the radiation decrease as a −4 .However, as demonstrated in Fig. 3, reheating is possible for both decays and scatterings when k ≥ 6.For this value of k, with α = 1 and β = 10 −2 , we see that reheating from decays occurs very late resulting in a very low reheating temperature, as seen in the right panel.Reheating from scatterings also occurs but for this value of β, it occurs even later, beyond the range shown in the plot.Though it is hard to see, the slopes for ρ R from scattering is shallower than ρ ϕ .
As one can infer from the k dependence in Eqs.(2.28) and (2.31), the reheating temperature for k > 4 increases with increasing k.This is demonstrated in Fig. 4, where we show the resulting reheating temperature as a function of k.Here we have assumed α = 1 and β = 1 for decays and scatterings, respectively.As already noted, the reheating temperature is too small for k = 6, but is considerably higher when k ≥ 8.The analysis presented above is perturbative, and in scenarios when reheating is significantly delayed (e.g. when k = 6), non-perturbative effects such as the fragmentation of the condensate become important [77][78][79][80][81][82].For the above perturbative analysis to hold, reheating must occur before the fragmentation of the condensate.The fragmentation of the condensate for the T-models considered here has been recently examined for k ≥ 6 [82].As one might expect, the relative value of the scale factor, a frag /a end increases substantially with k, whereas a RH /a end decreases with k.We find that the a RH ∼ a frag for k = 10 and α = 1.Note that a RH /a end ∝ α −(k+2)/(k−4) (a frag /a end is independent of α).For k = 12, reheating occurs before fragmentation for α ≳ O(0.1), and the limit on α weakens at higher k.Thus, for values of k sufficiently large to provide a reasonable reheating temperature, the effects of fragmentation can be safely ignored.Although similar constraints could be imposed on β, as evident from Fig. 4, achieving the same reheating temperature would require larger values of β, and the fragmentation limits discussed above would also correspond to increased values of β.

Production of gauginos
Before concluding this section, we recall that in a supersymmetric model, generally the decays of the inflaton to gauge bosons will be accompanied by decays to gauginos (so long as the channel is kinematically accessible).The results discussed above and in [8] assume that there is a dominant    decay channel with definite spin.However, it was shown in [25] that the partial width to gauginos is the same as the partial width to gauge bosons where y ′ is a dimensionless factor encoding the number of final states and the form of the gauge kinetic function.The above equality implies that the gauge boson and the gaugino result in exactly the same evolution of the reheating temperature, namely, T ∝ a − 3 8 for k = 2 and T ∝ a −1 for k ≥ 4. Compared to the case with only gauge bosons in the final state, the effect of a gauge boson-gaugino mixture in the final states amounts simply to doubling the decay rate in the ρ R and ρ ϕ evolution equations.In Fig. 5, we show the numerical solution of the reheating temperature when decays to gauge bosons and gauginos are included.As one can see, there are no changes in the evolutionary slopes.
It is interesting to note that for a generic Yukawa coupling yϕ f f where y is independent of the inflaton mass, the inflaton decay rate differs from Γ ϕ→λλ .The former is given by [8] Γ ϕ→ f f = y 2 8π m ϕ .
(2.33)   In this case, the temperature scales as T ∝ a −3/8 , a −3/4 , a −15/16 , respectively, for k = 2, 4, 6, as seen in the Table, and T ∝ a −1 for k ≥ 8.As a result, the qualitative behavior of the temperature is different for a gauge boson-gaugino mixture and for a gauge boson-matter fermion mixture in the range 4 ≤ k ≤ 6.In Fig. 6, we show the numerical solution for T in the case of gauge boson-matter fermion mixing for k = 2 (left) and k = 4 (right).As one can see, for k = 2, since all decay channels lead to the same scaling (T ∝ a −3/8 ), the lines are parallel, though for k = 4, they are not.
Note that the above discussion neglects the effects of the kinematic suppression due to the effective masses of the fermions and inflaton [8].When there is a non-zero field value for the inflaton, the coupling leading to the decay of the condensate also provides an effective mass to fermions (and  bosons, though not to gauge bosons). 6This effect leads to a suppression in the decay rate to fermions proportional to R −1/2 with R ∝ m 2 eff /m 2 ϕ .For example, for the generic Yukawa coupling, m eff = yϕ, and , and R ∝ (ϕ/M P ) (4−k) and can be quite large for k ≥ 6.Thus, for large k, we expect that our results for gauge boson final states are unaffected by the inclusion of fermions.

Gravitational production
In the previous section, we have considered the possible role for vector bosons in the reheating process.That is, given a coupling between the inflaton and some vector boson produced either by decay or scattering, the evolution of the radiation bath differs from the cases where the inflaton is coupled to fermions or scalars.Gauge invariance dictates a derivative coupling leading to a Planck mass suppression (in the case motivated by supergravity in Eq. (2.7)).A non-suppressed coupling to vectors would lead to an evolution identical to that in the case of scalars.However, within the Standard Model, there is an unavoidable coupling between a massive vector and the inflaton, mediated by gravity.Indeed, the scattering through the exchange of an s-channel graviton, as shown in Fig. 7, is a process that must be present, and therefore provides a minimal (unavoidable) production of vector bosons (and potentially a dark matter candidate) during the reheating process.In this section, we calculate the abundance of vectors produced solely through gravitational interactions.

The amplitude
To compute the gravitational production rate of vector bosons, we begin by first expanding the space-time metric around a flat Minkowski background.Using the relation g µν ≃ η µν + h µν being the canonically normalized graviton field, the interaction Lagrangian can be expressed as [39,41,83 where T µν SM , T µν ϕ , and T µν Aµ represent the energy-momentum tensors for the Standard Model particles, the inflaton, and the massive vector respectively.We will consider a massive vector field described by the Proca Lagrangian: and the canonical energy-momentum tensors for the massive fields with spins i = 0, 1/2, 1 are given by where V (S) is the potential of either the inflaton field or the SM Higgs boson, with S = ϕ, H, and The production of spin-3 2 particles was recently considered in [44].For the gravitational processes we parametrize the scattering amplitudes for the production rate of massive vectors as where i = 0, 1/2, 1 denote the initial state spins involved in the scattering process.Here, Π µνρσ is the graviton propagator for the canonical field h µν with momentum The partial amplitudes, M i µν , can be expressed by [39] M Here we keep our discussion completely general and argue that any inflationary model that agrees with the slow-roll parameters determined by Planck data [72] can be used, as long as it has a welldefined minimum and the potential can be expanded about its minimum as given by Eq. (2.8).Therefore, the currently favored Starobinsky model [28] (for k = 2) and α-attractor models [76] (for arbitrary k) are perfectly sufficient [71].

Gravitational production from the inflaton
The gravitational production of massive vector dark matter processes is illustrated by the Feynman diagram in Fig. 7.We consider first the process of vector dark matter production from the inflaton condensate ϕ+ϕ → A µ +A µ .For the case of a quadratic potential, the inflaton behaves like a massive particle at rest with four-momentum p 1,2 . 8We compute the partial amplitude using Eq.(3.9) with V (ϕ) given by Eq. (2.8) and use the inflaton condensate zero mode.For more details, see [41].The dominant particle production typically occurs right after the end of inflation when the amplitude of the oscillation is maximal.The gravitational particle production is Planck-suppressed, but the vector dark matter production from the inflaton is still significant since the energy density of the inflaton is huge at the beginning of reheating and a significant amount of it is transferred during the beginning of the reheating process.

Quadratic potential minimum
To begin, we first focus on the simplest case with a quadratic potential minimum given by V (ϕ) To compute the rate, we first find the square of the matrix element in Eq. (3.7) and use M 0 ρσ for the incoming state of the inflaton.For the inflaton condensate, we assume that the incoming 3-momentum of the inflaton vanishes, with ⃗ p 1,2 = 0, and proceed to compute the amplitude element squared |M 01 | 2 by summing over all polarizations of the outgoing vector states.We find that the matrix element squared is given by:9 where in the second equality we used the Mandelstam variable s = 4m 2 ϕ and we introduced the definition τ ≡ m 2 A /m 2 ϕ .It is worth noting that the massless limit of the result in Eq. (3.12) is finite.While there is no gravitational production of massless fermions due to helicity conservation, or to massless gauge bosons due to conformal invariance [36,39,45], the production of spin- 3  2 particles (raritrons) diverges in the massless limit [44] due to the pathology of spins greater than one.As we will see, the massless limit of (3.12) corresponds to the amplitude for producing a real scalar.We will return to this question when we consider the Stueckelberg mechanism for massive vectors in Section 3.5.Note that the above expression for |M| 2 is summed over both transverse and longitudinal modes.We can distinguish between the production of transverse (±1) and longitudinal (0) modes by introducing the following forms of the polarization vectors The relative contributions to |M| 2 break down as: 2τ 2 for the transverse mode and (2 − τ ) 2 for the longitudinal mode, each with the same prefactor of (m ϕ /2M P ) 4 .Thus we see that for massless gauge bosons, the transverse contribution vanishes due to the conformal structure of the gauge kinetic term.The remaining contribution in the massless limit found the longitudinal mode corresponds to the scalar degree of freedom which is produced in the massless limit [39,48].
Using Eq. (3.12), we find that the production rate, R ϕ k , for a quadratic minimum with k = 2, can be expressed as [84] where the outgoing momentum is given by A , and combining it with the matrix element squared (3.12), one obtains [56] Here we included the factor of 2 in the numerator to explicitly account for the 2 produced massive vectors per annihilation.10

General Potentials
As noted earlier, the gravitational production of the massive vector bosons from the inflaton condensate depends on the shape of the potential minimum.We discuss the derivation of the Boltzmann equation for a decaying inflaton condensate in Appendix B. Therefore, in this subsection, we consider more general potentials of the form (2.8).From Eq. (B.11), we can express the potential as V (ϕ) = V (ϕ 0 ) • P(t) k .When we Fourier expand P(t) k as in Eq. (B.12), we can write the potential in terms of its corresponding Fourier modes [8,75,86] where ⟨ρ ϕ ⟩ is the mean energy density averaged over the oscillations and ω is the oscillation frequency of ϕ, given by [8] with m 2 ϕ = ∂ 2 V /∂ϕ 2 | ϕ 0 .To calculate the scattering rate of the inflaton condensate, we follow the approach outlined in Appendix B. We find that the massive vector production rate is given by with where E n = nω is the n-th mode energy of the inflaton oscillation.For the k = 2 case, with ω = m ϕ (see Eq. (3.19)) and P(t , only the second mode in the Fourier expansion contributes to the sum.Using E 2 = 2m ϕ , we find that the production rate (3.20) reduces to Eq. (3.17).For other values of k, when m A ≪ m ϕ , Σ k 1 ≃ |P k,n | 2 .For example, for k = 4(6), Σ k 1 ≃ 0.063 (0.056) [39].The massive vector boson production rate (3.20) can then be decomposed as where the transverse spin ±1 contribution can be expressed as and the longitudinal spin 0 contribution is We note that the sum of transverse and longitudinal components satisfy the expression (3.21), with Σ k 1 = Σ k 1,1 +Σ k 1,0 .We also find in Eqs.(3.23) and (3.24) the behaviour we noted above, namely that in the limit m A → 0, the production of transverse modes vanishes, and only the longitudinal component is produced gravitationally.This result is illustrated in Fig. 8, which shows the evolution of the transverse and longitudinal production rates as a function of the ratio τ in the case of the quadratic potential, k = 2. What is striking is that production of the longitudinal mode greatly dominates for τ < 1.Already for τ = 1/2, 90% of the production of the vector is through the longitudinal mode.This agrees with the results found in [46,47,52,53] using very different techniques.Also note that the absolute value of the rates depends on the Fourier coefficients P n,k , which themselves become very similar for large values of k, hence we do not expect big differences for larger values of k.

Relic Abundance
To compute the relic abundance of massive vector fields, we start with the Boltzmann equation where n is the dark vector number density.If we introduce the yield Y ≡ a 3 n, we can rewrite the Boltzmann equation as a function of the scale factor, Assuming that the inflaton energy density dominates the total density, we can write and ρ ϕ (a) given by Eq.(2.20), the Boltzmann equation (3.26) can be expressed as (3.28) We first focus on the simplest case k = 2, which leads to ρ ϕ ∼ a −3 , ρ R ∼ T 4 ∼ a −3/2 , and the inflaton mass m 2 ϕ = 2λM 2 P .Using R ϕ 2 (a) given by Eq. (3.17), Eq. (3.28) is easily integrated and we obtain where we assumed that a RH ≫ a end and g RH is the number degrees of freedom at reheating.We can then compute the relic abundance [84] Ωh where we used g 0 = 43/11 and g RH = 915/4, together with the assumption m A ≪ m ϕ .Here, we normalized the values of m ϕ and ρ end for an α-attractor model of inflation with k = 2, though it weakly depends on T RH , see Refs.[7,39,71] for a detailed discussion.
We can compare the result obtained in (3.31) with that obtained for the relic density of scalars [36,39], fermions [36,39], or raritrons [44].For scalar dark matter, the relic density is approximately the same as the expression in Eq. (3.29) for vectorial dark matter (one need only replace (1 − τ + 3 4 τ 2 ) → (1 + τ ) 2 ) because it is the longitudinal mode which is mainly produced when T RH ≲ m ϕ .They are identical when in the limit τ → 0. In contrast, the fermion (χ) relic density is suppressed by a factor proportional to m 2 χ /m 2 ϕ (in this case, one need only replace (1 ).Thus fermion masses 10 4 times larger are required before saturating the matter content of the Universe.On the other hand, the relic density of a raritron of mass 10 7 GeV leads to Ω 3/2 h 2 ≳ 10 11 [44] due to the very efficient production of longitudinal modes with R ϕ 2 3/2 ∝ m 2 ϕ /m 2 3/2 .From Eqs. (3.22) -(3.24) we can repeat the exercise of separating the production of transverse and longitudinal modes.For the transverse contribution for k = 2, we find and  As expected and discussed previously, the gravitational production of the transverse mode is completely negligible.Similarly, we can compute the longitudinal contribution Finally, we can also generalize this result to cases when k ̸ = 2.We find that the number density can be expressed as which is identical to the result for a scalar with the replacement of Σ k 1 → Σ k 0 which are in fact equal in the limit m A → 0 [39].The dark matter abundance becomes As expected, for k = 2 this expression reduces to Eq. (3.31).
We illustrate these results in Fig. 9 where we show the lines of constant Ω 1 h 2 on the (m A , T RH ) plane. 11The line for k = 2, corresponds roughly to m A T RH ≲ 10 17 GeV 2 .For k = 4, there is no dependence on the reheating temperature, and there is an upper limit to the vector mass, m A ≲ 120 GeV.For k = 6, we have m A /T 1/3 RH ≲ .0017GeV 2/3 .We also note that our analysis only captures the post-inflationary production, corresponding to the UV modes in the particle number abundance spectrum.For a more detailed discussion that includes the IR particle power spectrum contribution arising during inflation, see Refs.[47,50,53].

Thermal Production of Massive Vector Dark Matter
A second production mechanism for producing massive vector dark matter arises from the Standard Model thermal background, SM + SM → A µ + A µ , also depicted in Fig. 7, and uses the partial amplitudes (3.9)- (3.11) for SM particles on the right-hand side of Eq. (3.7).The scattering processes of SM particles include the Higgs scalars, gauge bosons, and fermions in the initial state.Given our assumption of the MSSM, massless superpartners are also included in these initial states.Since the initial particle momenta p 1 and p 2 are large (of order m ϕ ) at the beginning of reheating, dominating < l a t e x i t s h a 1 _ b a s e 6 4 = " V y o q w h o H f u h Y d j J y I X 4 l u 7 K x e t A = " > A A A B 6 3 i c d V B N S 8 N A E J 3 4 W e t X 1 a O X x S J 4 C k l b 0 l y E o h e P F e w H t K F s t t t 2 6 e 4 m 7 G 6 J p e V o q D U t m x P b 9 e 9 T 3 k 2 K 7 r e L V q T q q e 7 z r I t Z 0 F y r B C c 1 B 6 7 w 8 j k g g q D e F Y 6 5 7 r < l a t e x i t s h a 1 _ b a s e 6 4 = " w A x p T 2 / m P V J q q M m 6 L K over electroweak scale quantities, we can assume these initial particle states to be effectively massless.We provide a detailed computation of the thermal rates in Appendix C.
For scalar initial states, including Higgses and sfermions, we use Eq.(C.3) with the association ϕ → s and set m ϕ = 0 (i.e., we neglect the scalar masses, and all MSSM masses relative to the reheating temperature in the thermal bath).Similarly, we use Eqs.(C.4) and (C.5) for the initial massless fermions and gauge bosons.To calculate the production rate, we integrate the amplitude squared over the incoming states [84,87,88], where we assume that s ≫ 4m 2 A , with s = (p 1 +p 2 ) 2 .Here the factor of two accounts for two massive vectors produced per scattering, E i denotes the energies of the initial and final state particles, θ 13 and θ 12 are the angles formed by momenta p 1,3 (in the center-of-mass frame) and p 1,2 (in the laboratory frame), respectively, and dΩ 13 = 2π d cos θ 13 .We assume the following thermal distributions for the incoming particles and obtain the expression for |M| 2 by summing over the three amplitudes associated with the different spins of the initial states where N b , N f , and N V represent the number of bosonic, fermionic, and vector fields in the model, respectively.For SM initial states, these values are N b = 4, N f = 45, and N V = 12.For MSSM initial states, which we focus on here, the corresponding values are N b = 98, N f = 61, and N V = 12.
The matrix elements and their computation is detailed in Appendix C. The gravitational production rate can then be expressed as with the coefficients β i given by Eqs.(C.9-C.11) for SM and Eqs.(C.12-C.14)for MSSM.
Since the thermal production rate is strongly sub-dominant compared to the production from the inflaton condensate, we consider only k = 2.To compute the thermally-produced vector number density, we replace the rate in Eq. (3.28) by the thermal production rate (3.41).The temperature can be expressed as function of the scale factor by solving We find that the thermally-produced massive vector number density is given by As before, we assumed that 4m 2 A ≪ s, which approximately corresponds to m A ≲ T RH , and we integrated Eq. (3.28) between a end and a RH .Since the β 1 terms dominates the thermal production when m A ≪ T RH , we only keep this term, and using Eq.(3.30) for the relic abundance, we obtain Comparing Eq. (3.44) with Eq. (3.31), we conclude that thermal production is negligible compared to the production from the inflaton condensate for k = 2.The result is also valid for larger values of k.Consequently, the relic density generated by the thermal source is negligible compared with that generated by inflaton oscillations across the entire parameter space.This result remains valid for k > 2.

The Vector Mass from the Stueckelberg and Higgs Mechanisms
Our discussion so far focused on a massive vector field described by the Proca Lagrangian (3.2), with a generic mass term.Despite its simplicity, the Proca Lagrangian has several disadvantages.Contrary to the massless case, the gauge invariance is explicitly broken by the mass term and one can not smoothly apply the massless limit as one loses one degree of freedom in the massless limit.On the other hand, the mismatch of the degrees of freedom can also be seen when we consider the squared amplitude for the production of A µ , e.g., Eq. (3.12), when the vector mass is sent to 0. An alternative description of a massive vector field is given by Stueckelberg [89,90].In fact, the Proca Lagrangian can be made gauge invariant through the field redefinition: so we get The new Lagrangian is gauge invariant under the transformations We easily recover the Proca Lagrangian by setting the unitary gauge σ = 0. On the other hand, the massless limit now gives rise to a massless vector field along with a scalar, hence the number of degrees of freedom is conserved.Indeed, the new scalar corresponds to the longitudinal mode of the massive vector, that drops out when m A → 0.
To quantize the theory, we add to the Stueckelberg Lagrangian (3.46) a gauge-fixing term: from which we recover a massive vector field and a decoupled massive scalar.The Feynman gauge choice amounts to ξ = 1 and the unitary gauge corresponds to ξ → ∞.We still have three degrees of freedom because, in addition to the subsidiary condition ∂ µ A µ = −ξm A σ, the gauge invariance (3.47) also holds, but now the gauge parameter must satisfy ∂ 2 + ξm 2 A Λ = 0.The Stueckelberg mechanism can allow for, e.g., a massive abelian gauge boson whose mass is not provided by a Higgs boson through spontaneous symmetry breaking.Examples of the Stueckelberg extension to the SM or supersymmetry can be found in [91][92][93].
Assuming minimal coupling to gravity, it is possible to generalize this Lagrangian to the curved space [94] and linearly perturb the metric around Minkowski, as we have done in the beginning of this section, then the canonically-normalized graviton field h µν will not only couple to A µ , but also to σ.For simplicity, we choose the Feynman gauge ξ = 1, so the action becomes The Stueckelberg scalar is decoupled and has the same mass as the vector boson.Note also that in this gauge, the polarization sum for A µ reads λ ϵ λ,µ ϵ * λ,ν = −g µν , as if A µ was "massless".In the Minkowski limit, the energy-momentum tensor is: with (3.52) One can easily check the conservation ∂ ν T µν = 0 upon using the Klein-Gordon equations.The energy-momentum tensor for σ is well-known and is identical to that in Eq. (3.3) when In contrast, the energy-momentum tensor for A µ is different from that derived from the Proca Lagrangian, Eq. (3.5).In the following, we compute the production rate of Stueckelberg vector field A µ from the process in Fig. 7, and compare that to the production of the Proca vector field A µ .
The partial amplitude for A µ is now given by After some involved algebra, one can verify the simple relation between squared amplitudes: which holds for any massive or massless initial states of spin 0, 1/2, 1.The above relation indicates that the production of the Stueckelberg A µ differs from the Proca A µ by a scalar contribution.To understand this point, recall that the Stueckelberg Lagrangian is obtained from the Proca Lagrangian by a redefinition of A µ and adding a gauge fixing term The two Lagrangians are equivalent if one requires for the former that the physical states be those for which (see [90]): Therefore, removing from the amplitude |M| 2 Stueckelberg, A the unphysical scalar state, which is of mass m A , we obtain |M| 2 Proca, A .Once the unphysical state is removed, the production rate for the vector with a Stueckelberg induced mass is the same as that discussed in Section 3.2 for the Proca Lagrangian.
Of course, another way to generate a gauge boson mass is through the Higgs mechanism.The relevant Lagrangian is where H denotes the complex Higgs scalar.In this case, A µ acquires a mass m A = ev when the Higgs is set to its vacuum expectation value v/ √ 2. Going to the unitary gauge and ignoring higher point interactions, one has a gravitational pair production channel for A µ and another for h (the physical Higgs).In this case, the squared amplitude for A µ production is just |M| 2 Proca, A .Interestingly, the sum of A µ and h productions will coincide with |M| 2  Stueckelberg, A only if the vector and scalar fields have the same mass.Note also that, when h is heavier than A µ , a secondary production of A µ may occur by the decay or annihilation of the produced h, which is absent in the Stueckelberg case.

Summary
Inflation is defined by the accelerated expansion of the Universe.To some extent, models of inflation can be phenomenologically triaged using the slow-roll parameters, which are determined during the final stages of the exponential expansion.A model of inflation embedded in a UV theory, which includes the SM (or the MSSM), can be further triaged by its ability to successfully reheat and produce an early stage of radiation domination, and perhaps participate in the process of dark matter production.When examined in detail, the reheating process is not instantaneous and depends on the couplings of the inflaton to matter.Reheating through inflaton decays is most efficient.However the evolution of the radiation bath may depend on the spin of the final state decay products [8].
Here we have concentrated on the reheating properties of so-called T-models of inflation [76], which yield predictions for the inflationary observables similar to that predicted by the Starobinsky model and are almost independent of k [8].For k = 2 the radiation density scales as a −3/2 , independent of the spin of the decay product.Reheating to gauge bosons then is relatively efficient [25].However for larger k, there can be significant differences, as is easily seen in Table 1.For k = 4, reheating through the decay to gauge bosons does not occur, though it does for decays to scalars and fermions.Note that decays to gauginos also do not allow for reheating for k = 4.For k > 6, decays to gauge bosons can lead to reheating, though its efficiency is poor and the reheating temperature is low unless k is relatively large as we have shown in Fig. 4.
We have also considered the mechanism of the gravitational production of vectors from the inflaton condensate as well as the thermal bath.The latter was found to be highly sub-dominant.The production of the vectors from the condensate is saturated by the production of the longitudinal mode and is equivalent to the production of a scalar in the limit the vector mass vanishes.We have also discussed the distinction in the production process depending on whether the vector mass is generated by either the Stueckelberg or Higgs mechanism.
In this case, the decay rate is given by where the effective coupling can be expressed as

B Boltzmann Equation for a Decaying Inflaton Condensate
In this appendix, we describe the derivation of the evolution equation for the energy density of a decaying inflaton condensate.Assuming the inflaton decay is a perturbative process, the inflaton condensate, ϕ, is spatially homogeneous, and one can express its phase space distribution (PSD) as f ϕ (k, t) = (2π) 3 n ϕ (t)δ (3) (k), where n ϕ represents the instantaneous inflaton number density.When neglecting the effects of Bose enhancement or Pauli blocking for the decay products of ϕ, the integrated Boltzmann equation for the number density can be expressed as follows [95]: Here A, B denote the inflaton decay products, dΨ ϕ,A,B represents the phase space measure for the products A, B and the condensate ϕ, and M denotes the transition amplitude.We note that there are no backreaction effects that would produce the inflaton particles in the condensate.We can express the right-hand side of Eq. (B.1) as where M n represents the transition amplitude for each oscillating field mode of ϕ during one oscillation, transitioning from the coherent state |ϕ⟩ to a two-particle final state |A, B⟩.Next, we introduce the following inflaton condensate normalization where p n = (E n , 0) and E n represents the n-th oscillation mode energy.
The production rate can now be readily obtained from the energy transfer rate from the inflationary to the decay product sector.We introduce the equation of state parameter w ϕ = p ϕ /ρ ϕ for the inflaton field, and the energy density evolution of ρ ϕ can be expressed as where S is a symmetry factor for the final states.
We express the oscillating inflaton field as ϕ(t) ≃ ϕ 0 (t) • P(t) , (B.11 where P is the oscillatory contribution and ϕ 0 is a slowly-evolving envelope that redshifts with time.One can approximately treat ϕ 0 as a constant, given its relatively slow evolution over a single oscillation.Next, we decompose the oscillatory contribution as where ω is the inflaton oscillation frequency.For an oscillating inflaton field, one can show that the equation of state parameter is given by [8] w We use these expressions in following appendix to compute the decay and scattering rates for different inflaton and vector boson couplings, as well as different values of k.

C Amplitudes and Thermal Rates
In this appendix, we calculate the thermal production rate of massive vector fields resulting from the scattering of SM or MSSM states, SM/MSSM+SM/MSSM → A µ +A µ .These states are assumed to be massless since the initial particle momenta p 1 and p 2 are typically of order of the inflaton mass m ϕ at the beginning of reheating.We first determine the squared amplitudes for different spins of the initial state particles.Following that, the thermal rate is derived using the general expression (3.38).We assume that 4m 2 A ≪ s and include a factor of 2 to account for the production of 2 particles per scattering.
The squared amplitudes are expressed in terms of the Mandelstam variables s and t, with The general squared amplitude for thermal processes involving SM initial states is given by Eq. (3.40), where we include 4 degrees for 1 complex Higgs doublet, 12 degrees for 8 gluons and 4 electroweak bosons, and 45 degrees for 6 (anti)quarks with 3 colors, 3 (anti)charged leptons and 3 neutrinos.
For MSSM, we include 98 degrees for squarks, sleptons, and Higgs bosons, 61 degrees for quarks, leptons, gauginos and higgsinos, and 12 degrees for gluons and electroweak bosons.We note that the squared amplitudes include the symmetry factors of both the initial and final states, and this is indicated with an overbar.
When summing over all polarizations, the total squared amplitude of the gravity-mediated production of massive vector from massless scalar is given by We parametrize β i , with i = 1, 2, 3, as the expansion rate H ∝ ρ 1 2

Figure 1 :
Figure 1: The evolution of the energy density of inflaton, ρ ϕ (solid red), and radiation, ρ R (blue dashed for decay and green dot-dashed for scattering) as a function of the scale factor a/a end for k = 2, α = 10 −2 and β = 10 −2 (left panels).The evolution of the temperature of the radiation is shown in the right panels.The star indicates the moment when ρ R (T RH ) = ρ ϕ (T RH ).

Figure 4 :
Figure 4: The reheating temperature defined when ρ ϕ = ρ R as a function k for k ≥ 6 in the case of decays (solid blue) with α = 1 and scatterings (solid orange) with β = 1.

Figure 5 :
Figure 5:The evolution of the temperature as a function of a/a end for decays to either gauge bosons or gauginos (dashed lines) and the sum when both are included (solid lines) for k = 2 (left) and k = 4 (right).Here we have chosen α = 10 −2 and g ρ = 915/4.The star indicates the temperature and scale factor at reheating.

Figure 7 :
Figure 7: Feynman diagram for the production of massive vector particles through the gravitational scattering of the inflaton condensate or the Standard Model particle bath.

Figure 8 :
Figure 8: Longitudinal (solid) and transverse (dashed) gravitational production rates of massive vector bosons for k = 2 in units of R × M 4 P /ρ 2 ϕ as a function of τ = m 2 A /m 2 ϕ .

. 34 )
which for m A ≪ m ϕ , gives the result in Eq. (3.29) since the production of massive vector field is completely dominated by the longitudinal component and Ω 0 h 2 ≃ 0.12 × T RH 10 10 GeV ρ end (5.2 × 10 15 GeV)