IPARCOS-UCM-23-119 IFT-UAM/CSIC-23-128 Prepared for submission to JCAP Vector dark matter production during inflation and reheating Jose A. R. Cembranos,a Luis J. Garay,a Álvaro Parra-Lópeza and Jose M. Sánchez Velázquezb aDepartamento de Física Teórica and IPARCOS, Facultad de Ciencias Físicas, Universidad Complutense de Madrid, Ciudad Universitaria, 28040 Madrid, Spain bInstituto de Física Teórica UAM/CSIC, c/ Nicolás Cabrera 13-15, Cantoblanco, 28049, Madrid, Spain E-mail: cembra@fis.ucm.es, luisj.garay@ucm.es, alvaparr@ucm.es, jm.sanchez.velazquez@csic.es Abstract. Gravitational particle production of spectator fields due to the expansion universe during the inflationary and reheating phases of the early universe is of particular interest in the context of dark matter, since it allows to constrain the properties of the dark candidate by comparing the density of particles produced with the observed dark matter abundance. In such processes, tachyonic instabilities arise as a consequence of the coupling to the curvature, greatly enhancing mode production. In this work, we consider a massive vector field that is coupled to the curvature scalar and the Ricci tensor only, and study its gravitational production through inflation and reheating. We show how the mechanism is more efficient than in the case of a non-minimally coupled scalar field, giving rise to larger abundances. Moreover, we analyze the importance of the coupling to the Ricci tensor, which increases tachyonic instabilities in the system, and constrain the mass of the dark particle and the values of the coupling constants by comparing the corresponding abundance with observations. Keywords: Cosmology of Theories beyond the SM, Effective Field Theories, Classical Theories of Gravity ArXiv ePrint: - ar X iv :2 31 0. 07 51 5v 1 [ gr -q c] 1 1 O ct 2 02 3 mailto:cembra@fis.ucm.es mailto:luisj.garay@ucm.es mailto:alvaparr@ucm.es mailto:jm.sanchez.velazquez@csic.es https://arxiv.org/abs/- Contents 1 Introduction 1 2 Dynamics of a massive vector field in flat FLRW cosmologies 2 2.1 Dynamics of transverse modes 4 2.2 Dynamics of longitudinal modes 5 2.3 Background dynamics 6 3 Particle production of transverse modes 8 3.1 Solution to the transverse mode equation 8 3.2 Spectra of produced transverse modes 10 4 Particle production of longitudinal modes 11 4.1 Solution to the longitudinal mode equation 11 4.2 Spectra of produced longitudinal modes 13 5 Abundance of produced particles 16 6 Conclusions 18 A Parameters 20 B Slow-roll approximation for the solution to the mode equation 20 1 Introduction Although dark matter is one of the keystones of our modern explanation of astrophysical and cosmological phenomena, there are still many fundamental properties about its own nature that are unknown. Observational evidence suggests that any suitable dark matter model must interact weakly with the Standard Model fields, it must be non-relativistic in the present, and furthermore, it must have appropriate clustering properties in order to explain large scale structure formation [1–4]. However, due to the rather indirect nature of the observational evidence, fundamental properties such as the mass or the spin of the dark matter field are unknown and more theoretical work must be done in order to unveil them. Within the current scenario, it becomes imperative to explore different paths that may shed light into its fundamental properties. One of the most paved ways in the last years has been to explore how different creation mechanism for dark matter particles can constrain the hypothetical space of parameters for its mass, its coupling to gravity and other fields, or even its spin. Within all the theoretical proposals for plausible creation mechanisms, we want to draw attention to one which is inherent to quantum fields that have only gravitational interactions: gravitational particle production. Due to the dynamical nature of space-time, any field will undergo particle creation [5–7]. Moreover, this phenomenon becomes particularly important during the stages of inflation [8–10] and reheating [11–15] in the early universe, since the geometry of spacetime rapidly changes. Furthermore, if the field under consideration is free, there is no mechanism to dilute the created particles, as no interaction with other – 1 – fields is allowed, hence making it one of the most appealing mechanisms to explore for dark matter candidates. These ideas have been applied mainly to scalar fields [16–24], but also to vectors and fermions in the references [25–29]. Moreover, recent works have shown that if one wants to consider gravitational production to obtain robust constraints on the parameters of the field, a detailed characterization of spacetime dynamics is needed. Indeed, it has been analyzed in the few past years how the long-time ignored oscillations that the curvature undergoes during reheating or the specific slow-roll dynamics during inflation can change the theoretical production of particles by several orders of magnitude [30–34]. Hence it is natural to revisit the production of vector fields as dark matter candidates, using a more sophisticated description of the background geometry during the inflationary and reheating epochs of our universe. In this work, we analyze the gravitational production of a vector dark matter field which is only coupled to gravity (through the curvature scalar and the Ricci tensor) in the most generic way considering only renormalizable terms. The inflationary model considered involves a single inflaton field ϕ that slowly rolls down a quadratic potential and oscillates around its minimum in the transition to reheating. The dynamics is obtained after numerically solving the corresponding equation of motion without any approximations. Note that this model, in its simplest version, is ruled out by CMB observations, although it allows for comparison with previous literature, in particular with the recent work [34]. Nevertheless, the procedure can be applied similarly to other choices of potential. We consider that the spectator dark matter field is initially prepared in the Bunch-Davies state, and particle production takes place during the expansion of the geometry. This will generate a non-vanishing, comoving particle density that will reach a constant value once the expansion of the geometry becomes adiabatic. In order to extract this quantity, we solve the mode equation until the adiabatic regime is reached. Then we use the customary adiabatic vacuum prescription for obtaining the Bogoliubov coefficients relating the latter with the in vacuum. Additionally, we use the analytical slow-roll approximation to the mode equation derived in [34] in order to speed up the numerical calculations. Within our derivation, we can set constraints on the possible values for each of the coupling strengths and the mass of the produced quanta so that it does not lead to overproduction, by comparing the resulting abundance with observations. The remainder of this paper is organized as follows. In section 2, we describe the dynamics of a massive vector field in an expanding cosmology, and write its mode content in terms of an effective mass tensor that incorporates all couplings to the geometry. At the same time, we obtain the mode equations and characterize the dynamics of the inflaton. Particle production of transverse and longitudinal modes is obtained in sections 3 and 4, respectively, whereas we perform an analysis of the total abundance yield in section 5. Lastly, we elaborate our conclusions in section 6. Notation. We set MP = 1/ √ G, ℏ = c = kB = 1, and use the metric signature (−, +, +, +). Furthermore, greek indices µ, ν run from 0 to 3, while latin indices i, j run from 1 to 3. 2 Dynamics of a massive vector field in flat FLRW cosmologies Let us consider a massive vector field Aµ in a flat Friedmann-Lemaître-Robertson-Walker (FLRW) geometry [35–41] which is non-minimally coupled to the geometry, extending therefore the minimally coupled model presented in references [26–28]. – 2 – The interaction with gravity is taken into account through terms of the form RAµAµ and R̃µνAµAν in the action S = −1 2 ∫ d4x √ −g (1 2FµνF µν + m2AµAµ + γRAµAµ + σR̃µνAµAν ) , (2.1) where g is the determinant of the metric, m is the mass of the boson, Fµν = ∂µAν − ∂νAµ is the field strength, and γ and σ are the couplings to the Ricci scalar R and the traceless Ricci tensor R̃µν = Rµν − gµνR/4, respectively. The FLRW metric can be written, considering conformal time η and using Cartesian coordinates for the spatial sections, in a conformally flat form, gµν = a2(η)ηµν , (2.2) where cosmological time t follows from dt = a(η)dη, and ηµν is the Minkowski metric, which will be used in the following for raising and lowering indices. Introducing (2.2) in (2.1), we can explicitly write S = −1 2 ∫ d4x [1 2F µνFµν + MµνAµAν ] , (2.3) where we have defined the mass tensor Mµν as Mµν ≡ a2ηµν ( m2 + γR ) + σR̃µν . (2.4) This action leads to the equation of motion ∂µ∂µAν − ∂µ∂νAµ = MµνAµ, (2.5) which is a generalization of the Proca equation for massive bosons [42] (the latter being recovered in the case of vanishing couplings γ and σ). Homogeneity allows us to we write the field Aµ in momentum space through Fourier modes, Aµ(η, x) = ∫ d3k (2π)3/2 õ(η, k)eik·x, (2.6) with the conjugate field fulfilling õ(η, k) = Ã∗ µ(η, −k) in order to ensure that Aµ is real. From now on, we will work with õ but drop the tilde for clarity. Additionally, spatial components are divided in longitudinal and transverse modes defined as A = AT + AL k |k| , (2.7) with k · A = kAL and k · AT = 0. (2.8) Writing (2.5) in momentum space, we observe that A0 is not dynamical, but can actually be expressed in terms of the longitudinal part,[ k2 − M00 ] A0(η, k) = −ikA′ L(η, k). (2.9) Here, the prime denotes derivative with respect to conformal time. Moreover, as long as k2 − M00 is different than 0, one can write A0(η, k) = − ikA′ L(η, k) k2 − M00 . (2.10) – 3 – This condition requires that M00 be negative for all values of m, γ, and σ, which restricts the allowed parameter space of our theory. Note that this is indeed a fundamental limitation of the theory, and not an operational requirement, since positive values of M00 would lead to instability problems (see equation (2.11) below). Since the explicit value of M00 depends on the background dynamics, we will analyze the allowed region of parameters in subsection 2.3. Expanding (2.3) in Fourier modes and making use of (2.10), we arrive at S = 1 2 ∫ dη d3k (2π)3/2 [ |A′ T|2 − ( k2 + MTT ) |AT|2 − M00 k2 − M00 |A′ L|2 − MLL|AL|2 ] , (2.11) We note that the dynamics of transverse and longitudinal degrees of freedom are independent of each other, and as such can be treated separately. Moreover, the fact that a positive M00 leads to instabilities is clear in (2.11), since in this case, the kinectic term of the longitudinal mode becomes negative. 2.1 Dynamics of transverse modes The action for the transverse modes is given by ST = 1 2 ∫ dη d3k (2π)3/2 [ |A′ T|2 − ( k2 + MTT ) |AT|2 ] , (2.12) and leads to the equation of motion A′′ T(η, k) + [ k2 + MTT(η) ] AT(η, k) = 0. (2.13) The most general solution can be written as AT(η, k) = ∑ h=± [ aT,h(k)vT(η, k)ϵh + a∗ T,h(−k)v∗ T(η, k)ϵ∗ h ] , (2.14) where ϵh is a polarization vector (h = +1 or −1). Note that this expression fulfills the condition A∗ T(η, k) = AT(η, −k), as required by the reality of Aµ. In turn, this implies that vT(η, k) only depends on the norm of the 3-momentum. Therefore, the mode equation for AT boils down to v′′ T(η, k) + ω2 T(η, k)vT(η, k) = 0, (2.15) where ω2 T(η, k) = k2 + MTT(η) (2.16) and MTT = a2 ( m2 + γR ) + σR̃TT. (2.17) This has the shape of a harmonic oscillator equation with a time-dependent frequency, and for σ = 0, it has exactly the shape of the equation of the non-minimally coupled scalar field discussed in [34]. – 4 – 2.2 Dynamics of longitudinal modes On the other hand, the action for the longitudinal mode reads SL = 1 2 ∫ dηd3k (2π)3/2 [ − M00 k2 − M00 |A′ L|2 − MLL|AL|2 ] , (2.18) which can be written in a form similar to (2.12), provided we introduce the auxiliary field AL as AL(η, k) ≡ AL(η, k) f(η, k) , f(η, k) = √ −M00(η)√ k2 − M00(η) , (2.19) where the corresponding element of the mass tensor reads M00(η) = −a2(η) [ m2 + γR(η) ] + σR̃00(η). (2.20) Then, equation (2.18) can be expressed as SL = 1 2 ∫ dηd3k (2π)3/2 { |AL|′ 2 − [ MLL f2 − f ′′ f ] |AL|2 } . (2.21) As before, the field can generally be written as the linear combination AL(η, k) = aL(k)vL(η, k) + a∗ L(−k)v∗ L(η, k), (2.22) leading to the mode equation v′′ L(η, k) + ω2 L(η, k)vL(η, k) = 0, (2.23) where the frequency is given by ω2 L(η, k) = MLL(η) f(η, k)2 − f ′′(η, k) f(η, k) , MLL(η) = a2(η) [ m2 + γR(η) ] + σR̃LL(η). (2.24) We see here that, contrary to the transverse mode and the scalar field case, the shape of the longitudinal frequency is not simply k2 plus an η-dependent term, but indeed more involved, and temporal and wavenumber dependencies mix. Note that, after canonical quantization, the coefficients of the expansions of both the transverse and the longitudinal modes will be promoted to creation and annihilation operators. As long as solutions to their respective equations of motion are normalized as vkv′ ∗ k − v′ kv∗ k = i, these operators will fulfill the standard commutation relations. For any of the spatial modes, the relation between two different basis of solutions vk and uk—each of which expands the corresponding component, AT or AL, with different coefficients—is given in terms of the Bogoliubov transformation uk = αkvk + βkv∗ k, with |αk|2 − |βk|2 = 1, and similarly for the associated creation and annihilation operators, b̂k = α∗ kâk − β∗ k ↠k [43]. If the frequency is time-independent, βk = 0 and there is no particle production. However, if the geometry changes with time, two particlar solutions vk and uk will correspond, in general, to different notions of vacuum and particle [44]. As a consequence, the mean number density of b-particles in the a-vacuum will not be zero, but ⟨0a| a†,b k ab k |0a⟩ = |βk|2, (2.25) – 5 – from which the total mean density is obtained after integration over all wavenumbers. The latter will remain finite as long as |βk|2 → 0 faster than k−3 for increasing k. One can now associate each notion of particle to observers living before and after the expansion. We consider the Bunch-Davies vacuum as the initial state, so that the observer before the expansion measures no particles. Then, the abundance of dark matter resulting from the expansion of the geometry during inflation and reheating is obtained by computing the mean number of particles an inertial observer would measure after reheating in the state of the system. If expansion is adiabatic enough after this point, particle production will be negligible. As long as the dark matter field is non-interacting, one can extrapolate this abundance to the present and compare with observations. In order to solve the mode equations (2.15) and (2.23) to obtain the number of particles that are produced in each mode, it remains to specify the background dynamics, analyzed in the following subsection. 2.3 Background dynamics In order to describe the geometry during inflation and reheating we consider a chaotic inflationary model consisting of a single scalar field ϕ whose equation of motion, assuming homogeneity and isotropy, is given by 0 = ϕ̈ + 3H(t)ϕ̇ + ∂ϕV (ϕ), (2.26) where the potential is of the form V (ϕ) = 1 2m2 ϕϕ2, mϕ denotes the inflaton mass, H(t) ≡ ȧ(t)/a(t) is the Hubble parameter, and the dot is derivative with respect to the cosmological time t. The relation between cosmological and conformal time is obtained from η = η0 +∫ t 0 dt/a(t). Together with (2.26), we solve the corresponding Friedmann equation, H2 = 4π 3M2 P [ ϕ̇2 + 2V (ϕ) ] , (2.27) considering that the only contribution to the energy-momentum tensor comes from the inflaton. Therefore, the Ricci curvature scalar and the components of the Ricci tensor in terms of ϕ and η read R = 8π M2 P [ 4V (ϕ) − a−2ϕ′ 2 ] , R̃00 = 6π M2 P a−4ϕ′ 2, R̃ii = 2π M2 P a−4ϕ′ 2, (2.28) which allows us to explicitly write both (2.16) and (2.24). The dynamics is obtained by numerically solving equation (2.26), together with (2.27), up to the point in which particle production becomes negligible. Inflation starts at ηi and ends at η = t = 0, when reheating begins, leading to the characteristic oscillations of the inflaton field ϕ and, as a consequence, of the rest of the curvature-dependent quantities. This can be seen in figure 1. Note that, although they oscillate, the traceless Ricci tensor components have always positive sign. This fact will be important when discussing instabilities further below. Let us now discuss in which region of parameter space our theory is well-defined, starting with the coupling γ. Stability of initial conditions (which will be analyzed further below) requires γ ≥ 1/6. Any value of γ within this region is valid, as long as m and σ are such that M00 < 0. However, back-reaction has to be taken into account for very large values of γ (see ref. [45] for such an analysis in the context of a non-minimally coupled scalar field), and production of particles for γ ≳ 1 is qualitatively similar to that of γ = 1. Therefore, – 6 – 0 2 4 6 0 5 10 15 -1 0 1 2 3 0 1 2 -1 0 1 2 3 0 0.2 0.4 0.6 0.8 1 2 3 -0.25 0 0.25 0.5 1.5 2.0 2.5 0 1 2 1.5 2.0 2.5 0 1 2 Figure 1. Inflaton field ϕ(η) (upper-left panel), curvature scalar R(η) (upper-right panel) and components of the traceless Ricci tensor (bottom row), as functions of conformal time. The range of time corresponds to the end of inflation and the beginning of reheating. The parameters used for all figures in this article are given in Appendix A. we consider the range 1/6 ≤ γ ≤ 1. At the same time, we want to allow for the possibility of vanishing σ, so that we can compare with scalar field production, for example. This, together with the range of γ discussed above, sets a limit for the smallest mass we can consider, m = 0.5mϕ. This is clear when looking at the explicit form of M00, given in (2.20). The Ricci scalar oscillates during reheating around 0, and therefore the mass term in M00 has to compensate for this behavior so that M00 remains negative. On the other hand, the traceless Ricci tensor component is always positive or 0 (see figure 1). This means that a positive value of σ is going to contribute always against the mass term. Therefore, σ has to be negative if we want to keep the chosen range in γ while having m ≥ 0.5mϕ. Because of all the previous considerations, we restrict ourselves to the following region of parameter space: m ≥ 0.5mϕ, γ ∈ [1/6, 1], σ ≤ 0. (2.29) – 7 – One can consider smaller values of the mass while keeping the same allowed range in γ, but then the coupling σ has to become negative and σ = 0 is not allowed. At the same time, positive values of σ are possible, but in this case it is required that the mass of the test field becomes greater than 0.5mϕ. 3 Particle production of transverse modes The dynamics of the transverse modes resembles that of a scalar field, and thus we closely follow the analysis done previously in [33, 34], where now we have an additional coupling to take into account. 3.1 Solution to the transverse mode equation If we want to calculate (2.25) in order to obtain the number of gravitationally produced particles in transverse modes, we need to compare at the same time the two solutions corresponding to the (Bunch-Davies) vacuum at the beginning of inflation and to the one at the time when the process is over. In particular, making the evaluation of (2.25) at this final time requires solving equation (2.15) from the beginning of inflation, ηi, until a time ηf in the reheating era for which the expansion of the geometry is adiabatic enough so that there is almost no production of particles after this point. The transverse mode equation, however, has no analytical solution given the background described in subsection 2.3. Nevertheless, during the first stages of inflation, within the slow-roll of the inflaton field ϕ, one can use the analytical approximation vT,SR(η, k) ≃ √ πτ(η, k)/2eiπνH(1) ν (kτ(η, k)) , τ(η, k) = ∣∣∣∣∣ ωT(η, k) ωT,dS(η, k)(η−η∗)+η∗−η0 ∣∣∣∣∣, (3.1) which is compatible with Bunch-Davies initial conditions and valid until η∗, whose particular value depends on the wavenumber k, the mass of the field m, and the couplings γ and σ (see appendix B) . Here, ωT,dS is the frequency in a de Sitter geometry, ωT,dS(η, k)2 = k2 + µ2 (η − η0)2 , with µ2 = m2/H2 0 + 12γ, (3.2) where H0 = H(ηi) = 1/η0 and ν = √ 1/4 − µ2. This coincides, precisely, with the configuration of the geometry at the beginning of inflation. Note that if we write γ = ξ − 1/6 we recover the equations for the case of gravitational production of scalar particles in the same background, ξ being the coupling to R. Thus, the initial vacuum state becomes unstable for γ < 0 in the case of transverse mode production, being γ = 0 the conformal point (equivalent to ξ = 1/6 in the scalar field case). Equation (2.15) is then solved in the following two regions, η = { ηi ≤ η ≤ η∗, Slow-roll approximation, η∗ ≤ η ≤ ηf, Numerical solution, (3.3) where numerical computations is needed once the slow-roll approximation described above is no longer valid. Let us remark that the use of this analytical solution, studied in detail in [34], is crucial for alleviating the numerical workload, especially given that the space of parameters is now three-dimensional, in contrast to the non-minimally coupled scalar case. – 8 – 0.01 0.1 1 10 100 0.5 1.0 1.5 0.01 0.1 1 10 100 -5 -4 -3 -2 -1 0 Figure 2. Error squared ϵ2 k for the transverse modes as function of the wave number k and the field mass m, for γ = 1/6 (left) and γ = 1 (right), and σ = −1, for η∗ = −500mϕ. The behaviour for σ = 0 is qualitatively similar. Then, it remains only to specify the vacuum of the observer living at ηf, for which we use the customary zeroth-order adiabatic prescription [43, 44], uT(ηf, k) = 1√ ωT(ηf, k) , u′ T(ηf, k) = − 1√ ωT (ηf, k) ( iωT(ηf, k) + 1 2 ω′ T(ηf, k) ωT(ηf, k) ) . (3.4) Note that this prescription is a good notion of vacuum as long as the dynamics are adiabatic at ηf. The regime of adiabaticity is reached at a different point in time depending on the parameters k, m, γ and σ, and therefore one has to choose an end point ηf that fulfills the adiabaticity condition for all the parameter space considered. With this, we can calculate the number of produced particles, nT(m, γ, σ) = ∫ dk 2π2 k2|βT|2(k, m, γ, σ) = nT,SR + ∫ ∞ 0 dk 2π2 k2|βT,SR|2O(ϵ2 k), (3.5) with nT,SR = ∫ ∞ 0 dk 2π2 k2|βT,SR|2 being the result obtained when using the approximation vT,SR. The error in (3.5) due to the use of the approximate solution (3.1) is essentially given by ϵ2 k, with ϵk defined in appendix B, and depends on the value of k, m and ξ for a fixed η∗. This is shown in figure 2 for η∗ = −500mϕ on specific but representative values of the parameters for illustrative purposes. We note that ϵ2 k is essentially independent of the value of the mass m within this range, and is small for k > mϕ. For wavenumbers smaller than the mass of the inflaton, the error increases. However, the density of produced particles is suppressed in this range by the k2 factor in (3.5), and therefore nT(m, γ, σ) ≈ nT,SR is a good approximation. – 9 – 0.1 1 10 100 0.001 0.01 0.1 1 0.1 1 10 100 10-8 10-6 10-4 10-2 1 0.1 0.5 1 5 10 50 10-9 10-5 10-1 0.1 0.5 1 5 10 50 10-7 10-3 101 Figure 3. Spectra of produced transverse modes, in log scale, for a vanishing coupling σ. The upper panels show the spectra for several values of γ, given m = 0.5mϕ (left) or m = mϕ (right). On the other hand, the lower panels show the spectra for several values of m, given γ = 1/6 (left) or γ = 1 (right). 3.2 Spectra of produced transverse modes In this subsection we show the results of particle production of transverse modes, given the inflationary model described in subsection 2.3. In particular, we focus on the study of k2|βT|2 as it is the interesting quantity for particle density production. We first consider a vanishing value of the coupling σ, that is to say, no coupling to the traceless Ricci tensor. The resulting spectra can be found in figure 3, where the top panels – 10 – correspond to a fixed value of m, whereas the bottom panels are for a fixed value of γ. A logarithmic scale was used so that the UV and IR behaviour can be observed. Note that increasing the mass of the dark matter particle leads to less production in all scales. In particular, spectra for m = mϕ feature oscillations for large values of k, and fall off much faster. This falling is enhanced when the mass goes beyond this value. On the other hand, increasing the value of γ strengthens the coupling to the geometry, and therefore particle production increases. Interestingly, the maximum of the spectra depends on both the value of m and γ. Let us consider now a negative value of σ. In particular, resulting spectra for σ = −1 are depicted in figure 4. Note that a linear scale is used in this case for illustrating the magnitude of nT,SR. The behaviour is similar to the case σ = 0, but particle production is enhanced by the extra term in the frequency (2.16). Indeed, a negative value of σ implies that the traceless Ricci tensor term in the frequency always contributes as a negative term (see figure 1). This means that it causes tachyonic instabilities as well as enhances those coming from the oscillations of the curvature scalar around 0, which are precisely the main source of non-adiabaticity, as analyzed in [30–34]. Similarly, one also finds that a positive value of σ yields less production, since the corresponding term contributes always against tachyonic instability. Note that this is not in contradiction with the requirement M00 < 0. A large, positive value of σ has to be compensated with the values of m and γ. However, once within the validity of the theory, a negative value of σ will contribute towards the frequency becoming imaginary. 4 Particle production of longitudinal modes Production of longitudinal modes is more involved than that of transverse modes or a scalar field since the form of the frequency (2.24) is qualitatively different. However, in the limit when we approach the beginning of inflation, it can be written as a de Sitter-like frequency, which we will use to develop a slow-roll approximation similar to (3.1). 4.1 Solution to the longitudinal mode equation Let us again use the fact that the geometry approaches de Sitter at the onset of inflation, as we did in section 3. In such geometry, in which a(η) = −1/H0(η −η0) and MLL = −M00 = µ2H2 0 , the longitudinal frequency (2.24) becomes (we shift the conformal time η−η0 → η for simplicity in the notation) ω2 L,dS(η, k) = k2 + µ2 η2 − k2 ( 2η2k2 − µ2) (η2k2 + µ2)2 . (4.1) We can now define the dimensionless quantity x ≡ 1/(kη), so that ω2 L,dS(η, k) k2 = 1 + µ2x2 − x2(2 − µ2x2) (1 + µ2x2)2 . (4.2) Note that in the limit µ2x2 ≪ 1, we can expand (4.2) in powers of µ2x2, ω2 L,dS(η, k) k2 = 1 + µ2x2 − x2 ( 2 − µ2x2 ) [ 1 − 2µ2x2 + O ( µ4x4 )] , (4.3) which up to terms of order O(µ2x2) can be written as ω2 L,dS(η, k) k2 ≃ 1 + x2(µ2 − 2). (4.4) – 11 – 0 10 20 30 0 2 4 6 0 10 20 30 0 0.4 0.8 0 10 20 30 0 0.1 0.2 0.3 0 10 20 30 0 2 4 6 Figure 4. Spectra of produced transverse modes, in linear scale, for a fixed value σ = −1. The upper panels show the spectra for several values of γ, given m = 0.5mϕ (left) or m = mϕ (right). On the other hand, the lower panels show the spectra for several values of m, given γ = 1/6 (left) or γ = 1 (right). This allows us to write the frequency (4.1) in this limit as ω2 L, dS(η, k) ≃ k2 + µ2 − 2 η2 , (4.5) which is a good approximation for kη → kηi. Note that in this limit, the conformal case is, similarly to the scalar field, γ = 1/6, value below which the initial vacuum becomes unstable. – 12 – 0.01 0.1 1 10 100 0.5 1.0 1.5 0.01 0.1 1 10 100 -5 -4 -3 -2 -1 0 Figure 5. Error squared ϵ2 k for the longitudinal modes as function of the wave number k and the field mass m, for γ = 1/6 (left) and γ = 1 (right), and σ = −1, fixed the value of η∗ = −500mϕ. The behaviour for σ = 0 is qualitatively similar. This is precisely the reason why we required γ ≥ 1/6. Then, we can write an approximation to the solution during slow-roll for the longitudinal mode equation as vL,SR(η, k) ≃ √ πτ(η, k)/2eiπνH(1) ν (kτ(η, k)) , τ(η, k) = ∣∣∣∣∣ ωL(η, k) ωL,dS(η, k)(η−η∗)+η∗−η0 ∣∣∣∣∣, (4.6) where the de Sitter frequency is now given by (reintroducing η0) ωL,dS(η, k)2 = k2 + µ2 − 2 (η − η0)2 , with µ2 = m2/H2 0 + 12γ. (4.7) From this point on, the analysis follows closely what has been discussed in section 3. In particular, the error in the number density of produced longitudinal modes due to the use of (4.6) is shown in figure 5 for the same value η∗ = −500mϕ. Note that for these modes, the slow-roll approximation is slightly worse, although the behaviour is qualitatively similar to that of the transverse modes. However, the error becomes important again for low values of k, which are suppressed in the density, and the integral of the spectra is not significantly affected. 4.2 Spectra of produced longitudinal modes Let us now show the particle production spectra for longitudinal modes, according to the background dynamics discussed in subsection 2.3. As before, we analyze the quantity k2|βT|2 for different values of the parameters. – 13 – 0 20 40 0 25 50 75 100 0 20 40 0 0.5 1.0 0.1 1 10 100 10-14 10-9 10-4 101 0.1 1 10 100 10-8 10-5 10-2 101 104 Figure 6. Spectra of produced longitudinal modes for a vanishing coupling σ. The upper panels, in linear scale, show the spectra for several values of γ, given m = 0.5mϕ (left) or m = mϕ (right). On the other hand, the lower panels, in log scale, show the spectra for several values of m, given γ = 1/6 (left) or γ = 1 (right). Because the shape of the frequency is different, the spectra of produced particles are qualitatively different from those obtained in the case of transverse modes. Let us analyze the same cases in the same order. We first concentrate in the situation in which the σ coupling vanishes, represented in figure 6. Contrary to the transverse case, all spectra feature significant oscillations in k, not only for m = mϕ. As before, for fixed mass (top panels), particle production increases with the coupling γ, whereas it increases when decreasing the dark matter mass (bottom panels) given a fixed value of γ. Crucially, spectra for m = mϕ fall off faster, as it was observed before. Note that, in general, production of particles is much more important in the longitudinal mode than in the transverse ones: Spectra are broader and maxima can be several orders of magnitude larger. These effects are particularly noticeable – 14 – 0 40 80 0 50 100 150 0 20 40 0 0.1 0.2 0.3 0.1 1 10 100 10-9 10-6 10-3 1 0.1 1 10 100 10-10 10-6 10-2 102 Figure 7. Spectra of produced longitudinal modes for a fixed value σ = −1. The upper panels, in linear scale, show the spectra for several values of γ, given m = 0.5mϕ (left) or m = mϕ (right). On the other hand, the lower panels, in log scale, show the spectra for several values of m, given γ = 1/6 (left) or γ = 0.5 (right). for small masses and large values of the coupling γ, i.e., in scenarios in which production is enhanced (see top left and bottom right panels of figure 6). As anticipated, non-vanishing, negative values of σ will induce tachyonic instabilities, yielding larger production. This is illustrated in figure 7, where one can observe that oscillations in the spectra are much more violent as well, as compared to the σ = 0 case. Because spectra are wider, we resolve them up to k ≈ 250mϕ and restrict ourselves to γ ≤ 0.5 in order to – 15 – calculate the total density of particles produced. Similarly to the transverse modes, a positive value of σ would lead to less production, since instabilities would become less important. This amplification of the tachyonic behavior of the frequency through the coupling σ is absent in the non-minimally coupled scalar field case for obvious reasons, but it is an important mechanism for producing a large abundance of dark matter. We therefore devote the next section to the analysis of this quantity. 5 Abundance of produced particles As we have seen, the mechanism of gravitational particle production discussed above results in a non-vanishing density of particles of the test field Aµ. The total particle density is obtained after integrating k2|βk|2 over all values of k, as described in equation (3.5). If the dark matter field is non-interacting (i.e., it is a spectator field), and particle production becomes negligible after ηf, this comoving density will remain the same until today. In other words, it will be related to the physical density only by a factor that takes into account the dilution due to the expansion of the universe, nT,phys(ηtoday)a3(ηtoday) = nT (and similarly for the longitudinal mode). Writing the abundance today in terms of the background radiation temperature, one arrives at [33, 34] Ω(m, γ, σ) = 8π 3M2 P H2 today gtoday ∗S greh ∗S ( Ttoday Treh )3 m n(m, γ, σ) a3 reh , (5.1) where Ttoday, Treh and gtoday ∗S , greh ∗S are the radiation temperatures today and at the end of reheating, and the corresponding relativistic degrees of freedom, respectively. Note that n = nL + 2nT, in order to include in the abundance the production for the three modes of the vector field. Additionally, the value of the scale factor at the end of reheating can be obtained as function of the reheating temperature, which is a free parameter. Indeed, when radiation dominates, at ηreh, the Hubble rate can be expressed as H2 reh = 8π 3M2 P π2 30greh ∗S T 4 reh, (5.2) thus allowing to write ηreh and areh as function of Treh. Let us remark that in the whole parame- ter space considered, adiabaticity is reached way before the end of reheating (namely ηf < ηreh) as long as Treh ≲ 1013 GeV, thus providing an upper limit for the reheating temperature. We show the longitudinal and transverse abundances in figure 8, for two values of the coupling, σ = 0 and σ = −1, and fixed reheating temperature Treh = 102 GeV. As expected from the spectra analysis, the abundance of longitudinal modes is much larger, several orders of magnitude even, depending on the values of m and γ. At the same time, a negative value of the coupling σ increases both abundances by around one order of magnitude. It is also interesting to note that the transverse abundance changes slower with mass, until m ∼ mϕ is reached. More importantly, there is no qualitative change in the behavior of the longitudinal mode abundance in the region close to the conformal point γ = 1/6, σ = 0, in contrast to the scalar field studied in [34]. On the other hand, total abundances are depicted in figure 9, for σ = 0 (upper panels) and σ = −1 (bottom panels) as well as for two different reheating temperatures, Treh = 102 GeV and Treh = 104 GeV. As before, a negative value σ increases the total abundance (which is dominated by the longitudinal contribution), and so does increasing the value of the reheating – 16 – 0.5 0.6 0.7 0.8 0.9 1 1/6 0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9 1 0.5 0.6 0.7 0.8 0.9 1 1/6 0.2 0.3 0.4 0.5 0.5 0.6 0.7 0.8 0.9 1 -4 -3 -2 -1 0 1 2 3 4 Figure 8. Transverse (left panels) and longitudinal (right panels) abundances for σ = 0 (upper panels) and σ = −1 (bottom panels), and a reheating temperatures of Treh = 102 GeV. Note that ΩT includes only one transverse degree of freedom. temperature. We observe that gravitational production is able to reproduce observations for masses below the inflaton, and enhancing production (for example via a negative coupling σ) shifts the observed abundance towards larger masses. In general, gravitational production of vector fields is much more efficient than that of scalar fields, since one is able to explain observations for a dark matter candidate with a mass of the order of the inflaton mass, in contrast to the case analyzed in [34]. In particular, note that the reheating temperatures in figures 8 and 9 are much lower than standard values. If we were to reproduce the observed – 17 – 0.5 0.6 0.7 0.8 0.9 1 1/6 0.2 0.3 0.4 0.5 0.5 0.6 0.7 0.8 0.9 1 0.5 0.6 0.7 0.8 0.9 1 1/6 0.2 0.3 0.4 0.5 0.5 0.6 0.7 0.8 0.9 1 -4 -3 -2 -1 0 1 2 3 4 Figure 9. Total abundance for σ = 0 (upper panels) and σ = −1 (bottom panels), and two different reheating temperatures, Treh = 102 GeV (left panels) and Treh = 104 GeV (right panels). The dashed line corresponds to the observed dark matter abundance. abundance for a typical value of Treh (for instance, Treh ∼ 109 GeV), the mass of the dark matter candidate would have to be orders of magnitude higher than the inflaton mass. 6 Conclusions Gravitational production during inflation and reheating can be used as a mechanism to constrain the properties of spectator dark matter by comparing the corresponding abundance – 18 – with observations. Previous works have concentrated mostly on scalar spectator fields, although minimally coupled massive vector fields have also been considered. In this work, we have studied particle production of a massive vector boson that is coupled to both the Ricci scalar and the Ricci tensor, in the context of inflation and reheating, extending therefore previous works on the subject. We showed that such interactions can be rewritten in terms of a mass tensor that results in an effective, time-dependent mass in the equations of motion of the respective modes. Crucially, the fact that this effective mass has to be always possitive restricts the valid parameter space of the theory. In this way, the dynamics of a Proca field are generalized to incorporate the evolution of the relevant background. The latter is given by numerically solving the equation of motion of the inflaton field, without the need of assuming slow-roll, for a particular inflationary potential, although the procedure is completely independent of this choice. We made use of an analytical approximation to the solution of the mode equations of both transverse and longitudinal modes, valid during the slow-roll part of inflation, which allowed us to increase efficiency of calculations and to properly define initial conditions for the spectator field. As a result, we presented the spectra of produced particles together with the corresponding abundances obtained after dilution of the total density of produced particles after the adiabatic regime has been reached. The longitudinal mode production dominates over that of the transverse modes, being the leading contribution to the final abundance. As expected, production increases when decreasing the test field mass, and increasing the coupling to the curvature scalar, following a pattern similar to scalar fields. However, the coupling to the Ricci tensor can take negative values, leading to tachyonic instabilities that facilitate highly efficient particle production. This feature—absent when studying gravitational production of non-minimally coupled scalar fields—is especially relevant for the longitudinal case, which is generally more sensitive to instabilities. Overall, gravitational production of such a massive vector field is much more efficient and much more influenced by tachyonic behavior than that of scalar fields. As a consequence, the observed abundance is recovered for a heavy dark matter candidate (close to the inflaton mass) only if reheating temperatures are lower than typical values. Indeed, if one considers usual reheating temperatures, the mass of the dark matter particle would have to be orders of magnitude higher than that of the inflaton, since lighter particles would be overproduced. Acknowledgements This work was partially supported by the MICINN (Ministerio de Ciencia e Innovación, Spain) projects PID2019-107394GB-I00/AEI/10.13039/501100011033 (AEI/FEDER, UE), PID2020- 118159GBC44, and PID2022-139841NB-I00, COST (European Cooperation in Science and Technology) Actions CA21106 and CA21136. Additionally, Á. P.-L. is supported by the MIU (Ministerio de Universidades, Spain) fellowship FPU20/05603. J. M. S. V. acknowledges the support of the Spanish Agencia Estatal de Investigación through the grant “IFT Centro de Excelencia Severo Ochoa CEX2020-001007-S". Finally, J. A. R. C. acknowledges support by Institut Pascal at Université Paris-Saclay during the Paris-Saclay Astroparticle Symposium 2022, with the support of the P2IO Laboratory of Excellence (program “Investissements d’avenir” ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038), the P2I axis of the Graduate School of Physics of Université Paris-Saclay, as well as IJCLab, CEA, APPEC, IAS, OSUPS, and the IN2P3 master project UCMN. – 19 – A Parameters In this appendix, we provide all the parameters which were used in the numerical calculations. Most of the results are left in terms of the mass of the inflaton, mϕ. When necessary, we have taken mϕ = 1.2 × 1013 GeV, value that sets the scale of the problem. Therefore, the Planck mass MP has the value MP = 1.02 × 106mϕ. The inflaton field has an initial value of ϕi = 3MP . Inflation ends at t = 0, point at which, under the slow-roll approximation, the field woudl reach the value ϕ0 = 0.5MP . Therefore, one can extract the initial time of inflation, ti ≃ −15.35/mϕ. For solving the exact equation of motion (2.26), we take ϕ′(ti) = ϕ′ SR(ti), i.e., the value of the derivative as given by the slow-roll approximation with the conditions imposed above. Thus, ϕ(t = 0) slightly deviates from ϕ0. Furthermore, we make the choice a(t = 0) = a0 = 1. Lastly, we assume slow-roll is a good approximation until η∗ = −500/mϕ. Particle production is obtained at ηf = 7.58/mϕ, which is within the adibatic regime for all the parameter space explored, namely 0.5mϕ ≤ m, 1/6 ≤ γ ≤ 1 and σ ≤ 0. B Slow-roll approximation for the solution to the mode equation Let us elaborate on the approximation for the solution to the mode equations (3.1) and (4.6). In what follows, we write η − η0 as η, and drop the mode index k for clarity. The derivation, detailed in [Out Ref.], requires that the following quantities are small, given values of k, m, γ and σ. • Let us first define ϵ(m, γ, σ) = max η∈I1 ∣∣∣∣∣1 − ωSR(η; m, γ, σ) ωdS(η; m, γ, σ) ∣∣∣∣∣, with I1 = (−∞, η1), (B.1) where η1 is chosen such that ϵ ≪ 1. Then, we can define h(η; m, γ, σ) by ωSR ωdS = 1 + ϵh. (B.2) By construction, |h(η)| ≤ 1 for η ∈ I1. Moreover, h′(η) ≥ 0. • Second, we define the quantity δ(m, γ, σ) = max η∈I2 ∣∣∣h′(η; m, γ, σ)η ∣∣∣, with I2 = (−∞, η2), (B.3) where we choose η2 such that δ ≤ ϵ. Then, we introduce g(η; m, γ, σ) as h′(η) = δg(η) η , (B.4) for which again we have that |gk(η)| ≤ 1 for η ∈ I2. • Finally, ρ(m, γ, σ) = max η∈I3 ∣∣∣∣∣ω′ dS(η) ωdS(η)η ∣∣∣∣∣, with I3 = (−∞, η3), (B.5) and choose η3 such that ρ ≤ ϵ. Now, we take η∗ = min(η1, η2, η3) and I = (−∞, η∗), where I is the interval for which the three parameters ϵ, δ, ρ are small. Note that η∗ < 0 since inflation ends at η = 0. – 20 – • Additionally, we require that |η∗/η0| > 1. Provided the above is fulfilled, and translation of η to η − η0 is undone, one can write the solution to the mode equation as (3.1) and (4.6) up to terms of order O(ϵ2), where by ϵ we mean here the largest of the three small parameters defined above. References [1] Y. Sofue and V. Rubin, Rotation curves of spiral galaxies, Annual Review of Astronomy and Astrophysics 39 (2001) 137. [2] M. Bartelmann and P. Schneider, Weak gravitational lensing, Physics Reports 340 (2001) 291. [3] D. Clowe, A. Gonzalez and M. Markevitch, Weak-lensing mass reconstruction of the interacting cluster 1e 0657-558: Direct evidence for the existence of dark matter, The Astrophysical Journal 604 (2004) 596. [4] Planck collaboration, Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641 (2020) A10 [1807.06211]. [5] L. Parker, Quantized Fields and Particle Creation in Expanding Universes. I, Phys. Rev. 183 (1969) 1057. [6] L.H. Ford, Gravitational particle creation and inflation, Phys. Rev. D 35 (1987) 2955. [7] L.H. Ford, Cosmological particle production: a review, Rept. Prog. Phys. 84 (2021) [2112.02444]. [8] A. Starobinsky, A new type of isotropic cosmological models without singularity, Physics Letters B 91 (1980) 99. [9] A.H. Guth, Inflationary universe: A possible solution to the horizon and flatness problems, Phys. Rev. D 23 (1981) 347. [10] A. Linde, A new inflationary universe scenario: A possible solution of the horizon, flatness, homogeneity, isotropy and primordial monopole problems, Physics Letters B 108 (1982) 389. [11] E.W. Kolb and M.S. Turner, The Early Universe, vol. 69 (1990), 10.1201/9780429492860. [12] L. Kofman, A.D. Linde and A.A. Starobinsky, Reheating after inflation, Phys. Rev. Lett. 73 (1994) 3195 [hep-th/9405187]. [13] L. Kofman, A. Linde and A.A. Starobinsky, Towards the theory of reheating after inflation, Phys. Rev. D 56 (1997) 3258. [14] R. Allahverdi, R. Brandenberger, F.-Y. Cyr-Racine and A. Mazumdar, Reheating in inflationary cosmology: Theory and applications, Annual Review of Nuclear and Particle Science 60 (2010) 27. [15] D. Baumann and L. McAllister, Inflation and String Theory, Cambridge Monographs on Mathematical Physics, Cambridge University Press (5, 2015), 10.1017/CBO9781316105733, [1404.2601]. [16] D.J.H. Chung, E.W. Kolb and A. Riotto, Superheavy dark matter, Phys. Rev. D 59 (1998) 023501. [17] D.J.H. Chung, P. Crotty, E.W. Kolb and A. Riotto, Gravitational production of superheavy dark matter, Phys. Rev. D 64 (2001) 043503. [18] T. Markkanen, A. Rajantie and T. Tenkanen, Spectator dark matter, Phys. Rev. D 98 (2018) 123532. – 21 – https://doi.org/10.1146/annurev.astro.39.1.137 https://doi.org/10.1146/annurev.astro.39.1.137 https://doi.org/10.1016/s0370-1573(00)00082-x https://doi.org/10.1086/381970 https://doi.org/10.1086/381970 https://doi.org/10.1051/0004-6361/201833887 https://doi.org/10.1051/0004-6361/201833887 https://arxiv.org/abs/1807.06211 https://doi.org/10.1103/PhysRev.183.1057 https://doi.org/10.1103/PhysRev.183.1057 https://doi.org/10.1103/PhysRevD.35.2955 https://doi.org/10.1088/1361-6633/ac1b23 https://arxiv.org/abs/2112.02444 https://doi.org/https://doi.org/10.1016/0370-2693(80)90670-X https://doi.org/https://doi.org/10.1016/0370-2693(80)90670-X https://doi.org/10.1103/PhysRevD.23.347 https://doi.org/10.1103/PhysRevD.23.347 https://doi.org/https://doi.org/10.1016/0370-2693(82)91219-9 https://doi.org/10.1201/9780429492860 https://doi.org/10.1103/PhysRevLett.73.3195 https://doi.org/10.1103/PhysRevLett.73.3195 https://arxiv.org/abs/hep-th/9405187 https://doi.org/10.1103/PhysRevD.56.3258 https://doi.org/10.1103/PhysRevD.56.3258 https://doi.org/10.1146/annurev.nucl.012809.104511 https://doi.org/10.1146/annurev.nucl.012809.104511 https://doi.org/10.1017/CBO9781316105733 https://arxiv.org/abs/1404.2601 https://doi.org/10.1103/PhysRevD.59.023501 https://doi.org/10.1103/PhysRevD.59.023501 https://doi.org/10.1103/PhysRevD.64.043503 https://doi.org/10.1103/PhysRevD.98.123532 https://doi.org/10.1103/PhysRevD.98.123532 [19] M. Fairbairn, K. Kainulainen, T. Markkanen and S. Nurmi, Despicable dark relics: generated by gravity with unconstrained masses, Journal of Cosmology and Astroparticle Physics 2019 (2019) 005. [20] D.J.H. Chung, E.W. Kolb and A.J. Long, Gravitational production of super-hubble-mass particles: an analytic approach, Journal of High Energy Physics 2019 (2019) . [21] S. Hashiba and J. Yokoyama, Gravitational particle creation for dark matter and reheating, Phys. Rev. D 99 (2019) 043008. [22] N. Herring, D. Boyanovsky and A.R. Zentner, Nonadiabatic cosmological production of ultralight dark matter, Phys. Rev. D 101 (2020) 083516. [23] K. Kainulainen, O. Koskivaara and S. Nurmi, Tachyonic production of dark relics: a non-perturbative quantum study, Journal of High Energy Physics 2023 (2023) . [24] M.A. Garcia, M. Pierre and S. Verner, Scalar dark matter production from preheating and structure formation constraints, Physical Review D 107 (2023) . [25] P.W. Graham, J. Mardon and S. Rajendran, Vector Dark Matter from Inflationary Fluctuations, Phys. Rev. D 93 (2016) 103520 [1504.02102]. [26] Y. Ema, K. Nakayama and Y. Tang, Production of purely gravitational dark matter: the case of fermion and vector boson, Journal of High Energy Physics 2019 (2019) . [27] M. Bastero-Gil, J. Santiago, L. Ubaldi and R. Vega-Morales, Vector dark matter production at the end of inflation, Journal of Cosmology and Astroparticle Physics 2019 (2019) 015. [28] A. Ahmed, B. Grzadkowski and A. Socha, Gravitational production of vector dark matter, Journal of High Energy Physics 2020 (2020) . [29] T. Sato, F. Takahashi and M. Yamada, Gravitational production of dark photon dark matter with mass generated by the higgs mechanism, Journal of Cosmology and Astroparticle Physics 2022 (2022) 022. [30] Y. Ema, R. Jinno, K. Mukaida and K. Nakayama, Gravitational particle production in oscillating backgrounds and its cosmological implications, Phys. Rev. D 94 (2016) 063517. [31] T. Markkanen and S. Nurmi, Dark matter from gravitational particle production at reheating, Journal of Cosmology and Astroparticle Physics 2017 (2017) 008. [32] Y. Ema, K. Nakayama and Y. Tang, Production of purely gravitational dark matter, Journal of High Energy Physics 2018 (2018) . [33] J.A. Cembranos, L.J. Garay and J.M.S. Velázquez, Gravitational production of scalar dark matter, Journal of High Energy Physics 2020 (2020) . [34] J.A.R. Cembranos, L.J. Garay, A. Parra-López and J.M. Sánchez Velázquez, Late vacuum choice and slow roll approximation in gravitational particle production during reheating, JCAP 08 (2023) 060 [2301.04674]. [35] A. Friedman, Über die Krümmung des Raumes, Z. Phys. 10 (1922) 377. [36] A. Friedman, Über die Möglichkeit einer Welt mit konstanter negativer Krümmung des Raumes, Z. Phys. 21 (1924) 326. [37] A.G. Lemaître, A Homogeneous Universe of Constant Mass and Increasing Radius accounting for the Radial Velocity of Extra-galactic Nebulæ, Mon. Not. R. Astron. Soc. 91 (1931) 483. [38] H.P. Robertson, Kinematics and World-Structure, Astrophys. J. 82 (1935) 284. [39] H.P. Robertson, Kinematics and World-Structure II, Astrophys. J. 83 (1936) 187. [40] H.P. Robertson, Kinematics and World-Structure III, Astrophys. J. 83 (1936) 257. – 22 – https://doi.org/10.1088/1475-7516/2019/04/005 https://doi.org/10.1088/1475-7516/2019/04/005 https://doi.org/10.1007/jhep01(2019)189 https://doi.org/10.1103/PhysRevD.99.043008 https://doi.org/10.1103/PhysRevD.99.043008 https://doi.org/10.1103/PhysRevD.101.083516 https://doi.org/10.1007/jhep04(2023)043 https://doi.org/10.1103/physrevd.107.043530 https://doi.org/10.1103/PhysRevD.93.103520 https://arxiv.org/abs/1504.02102 https://doi.org/10.1007/jhep07(2019)060 https://doi.org/10.1088/1475-7516/2019/04/015 https://doi.org/10.1007/jhep08(2020)059 https://doi.org/10.1088/1475-7516/2022/08/022 https://doi.org/10.1088/1475-7516/2022/08/022 https://doi.org/10.1103/PhysRevD.94.063517 https://doi.org/10.1088/1475-7516/2017/02/008 https://doi.org/10.1007/jhep09(2018)135 https://doi.org/10.1007/jhep09(2018)135 https://doi.org/10.1007/jhep06(2020)084 https://doi.org/10.1088/1475-7516/2023/08/060 https://doi.org/10.1088/1475-7516/2023/08/060 https://arxiv.org/abs/2301.04674 https://doi.org/10.1007/BF01332580 https://doi.org/10.1007/BF01328280 https://doi.org/10.1093/mnras/91.5.483 https://doi.org/10.1086/143681 https://doi.org/10.1086/143716 https://doi.org/10.1086/143726 [41] R.M. Wald, General Relativity, Chicago University Press, Chicago (1984), 10.7208/chicago/9780226870373.001.0001. [42] M.E. Peskin and D.V. Schroeder, An Introduction to quantum field theory, Addison-Wesley, Reading, USA (1995). [43] N.D. Birrell and P.C.W. Davies, Quantum Fields in Curved Space, Cambridge Monographs on Mathematical Physics, Cambridge University Press, Cambridge (1982), 10.1017/CBO9780511622632. [44] V. Mukhanov and S. Winitzki, Introduction to Quantum Effects in Gravity, Cambridge University Press, Cambridge (2007), 10.1017/CBO9780511809149. [45] O. Lebedev, T. Solomko and J.-H. Yoon, Dark matter production via a non-minimal coupling to gravity, JCAP 02 (2023) 035 [2211.11773]. – 23 – https://doi.org/10.7208/chicago/9780226870373.001.0001 https://doi.org/10.1017/CBO9780511622632 https://doi.org/10.1017/CBO9780511809149 https://doi.org/10.1088/1475-7516/2023/02/035 https://arxiv.org/abs/2211.11773 Introduction Dynamics of a massive vector field in flat FLRW cosmologies Dynamics of transverse modes Dynamics of longitudinal modes Background dynamics Particle production of transverse modes Solution to the transverse mode equation Spectra of produced transverse modes Particle production of longitudinal modes Solution to the longitudinal mode equation Spectra of produced longitudinal modes Abundance of produced particles Conclusions Parameters Slow-roll approximation for the solution to the mode equation