Vector dark matter, inflation, and non-minimal couplings with gravity

We propose a cosmological dark matter production mechanism in the form of a longitudinal massive vector boson. We build upon the work [1] including non-minimal couplings of the massive vector with gravity, developing a well motivated set-up from an effective field theory perspective. We carefully track the dynamics of vector field in passing from inflation to radiation dominated universe to show that the late time abundance of longitudinal modes — excited initially by the quantum fluctuations during inflation — can provide the observed dark matter abundance for sufficiently weak non-minimal coupling and wide range of vector masses 5 × 10-7 ≲ m [eV] ≲ 5 × 103. The final abundance of dark matter depends on two parameter, the vector mass and its non-minimal coupling with gravity. We discuss experimental venues to probe this framework, including the production of a stochastic gravitational wave background. The latter is especially interesting, as the same mechanism that generates dark matter can potentially lead to the production of gravitational waves in the LISA frequency band, through the second-order effects of large dark matter iso-curvature perturbations at small scales. We take a first step in this direction, identifying the potential information that gravitational wave experiments can provide on the parameter space of dark matter within this scenario.


Introduction
The nature of dark matter is one of the most pressing open problems in physics.Current observations provide convincing evidence for gravitational interactions between dark matter belonging to a hidden sector, and the visible sector of the Standard Model [2][3][4].If gravity is the only mediator between the hidden and visible sectors, we need some efficient mechanism to produce dark matter in the early universe.Cosmic inflation is able to do so by means of the phenomenon of particle production in curved space.See e.g.[5,6] for classic examples exploiting this possibility.In this work we consider the recent, compelling scenario [1], in which dark matter corresponds to the longitudinal components of a massive vector field created during inflation.The longitudinal vector fluctuations have a distinctive kinetic term structure, and specific cosmological behaviour during inflation and radiation domination.Their spectrum grows towards small scales up to a peak, whose position and amplitude depends on the vector mass.From a cosmological viewpoint, the longitudinal vector fields correspond to iso-curvature fluctuations.Their size is small at large cosmic microwave background (CMB) scales, and satisfy CMB constraints.Besides its elegance, the scenario is compelling since its predictions depend on two parameters only: the present-day dark matter abundance depends only on the vector mass and on the energy scale of inflation.Various developments and further studies of the original scenario [7][8][9][10][11][12][13][14] have been considered.Additional early cosmological phases different from radiation enlarge the possible ranges for the vector mass [9,10], as well as additional couplings with particles in the dark sector [14].The drawback of these extensions is that often the dark matter abundance depends on additional free parameters, making the set-up less predictive.
We propose a generalization of [1] with an action including a non-minimal coupling of the vector field with gravity.The coupling we consider is quadratic in the vector field and it has dimension four as the vector field strength.It is the only non-minimal coupling with these characteristics, and which avoids Ostrogradsky instabilities.Hence, from an effective field theory point of view, once the Abelian symmetry is broken by a vector mass, our scenario is as minimal and well motivated as the original [1].Its predictions on the dark matter amplitude again depend on two parameters only: the vector mass and the strength of the non-minimal vector coupling to gravity.The cosmological dynamics of vector longitudinal fluctuations is more complex than [1] for the presence of gradient instabilities generated by the non-minimal couplings.We discuss in detail how to tame such instabilities, and obtain a working scenario with testable predictions.
One interesting aspect of dark matter constituted by longitudinal massive vector fluctuationsa dark photon -is that it can be tested in a variety of ways.From a particle physics perspective, the dark photon can be tested through its milli-charged couplings with the Standard Model of particle physics (see e.g.[15][16][17] for recent comprehensive reviews).The dark photon cosmological evolution has distinctive properties, see e.g.[18][19][20][21][22].It can have testable effects when surrounding black holes, or forming cosmic strings, see e.g.[23][24][25].It can be probed with accelerometers [26] and gravitational wave experiments, see e.g.[27][28][29][30][31].For the range of vector masses allowed by our set-up, we discuss for the first time the possibility of testing this proposal by means of gravitational waves in the LISA frequency band, induced at second order in perturbations by the vector longitudinal modes (see e.g.[32] for a general review on the subject).
We start our discussion with sections 2, 3, and 4 elaborating on our set-up in general terms, and the dynamics of longitudinal vector fluctuations during inflation and radiation domination.In section 5 we show that the abundance of dark matter in our scenario depends on two parameters only.Section 6 discusses phenomenological implications of vector dark matter in our set-up.Four appendixes cover technical aspects of our results.2 The set-up: Cosmology with non-minimally coupled massive vector field We study the dynamics of linearized vector fluctuations propagating through a cosmological space-time, controlled by the quadratic vector action [33,34] where A µ is a spin-1 vector field, F µν = ∇ µ A ν − ∇ ν A µ its field strength, and is the Einstein tensor.The vector A µ belongs to a dark sector independent from the Standard Model of particle physics.In our viewpoint, we consider A µ a dark photon (see e.g.[16] for a review on the subject), and we assume that the dark photon (and the dark sector in general) has only gravity-induced interactions with Standard Model matter.We use eq.( 2.1) to study the generation of longitudinal vector dark matter, extending the analysis of [1] by including the effects of the α 2 contributions.The additional term enriches the dynamics of vector fluctuations during radiation domination.It enlarges the parameter space of scenarios leading to the correct dark matter abundance, allowing for a wider range of vector masses with respect to [1].In order to do so, there is no need of introducing non-standard cosmological eras, or to make any hypothesis on the reheating process occurring between the end of inflation and the onset of radiation domination [9,10].Importantly, the final dark matter abundance depends only on the parameters m and α appearing in eq.(2.1), with no dependence on the underlying cosmology (not even the scale of inflation).This makes a dark matter scenario based on the theory (2.1) very minimal, and hopefully easier to test since there are only two parameters to constrain.
The Abelian symmetry A µ → A µ + ∂ µ λ(x) of the field strength is broken by the vector mass term proportional to m, and by the non-minimal coupling to gravity proportional to α 2 (α being a dimensionless quantity).From an effective field theory perspective, once the Abelian symmetry is broken by the vector mass, the contribution proportional to α 2 in action (2.1) is unavoidable.In fact, it is normally generated by loop corrections, since it respects the (lack of) symmetries of the system.It has dimension four as the remaining terms in the Lagrangian density, hence its contributions are in principle as important as the ones of the other terms, since it is not suppressed by a high-energy scale.Furthermore, it is the only operator quadratic in the vector field that non-minimally couples the vector system to gravity, without introducing Ostrogradsky instabilities [33,34].Hence, our set-up is theoretically well motivated, and we can consider the structure of action (2.1) to be as minimal as the one of the original scenario [1].We intend to quantitatively determine the physically interesting ranges for the parameter α which allow us to realize dark matter in this set-up1 .
To study the cosmological evolution of vector fluctuations, we expand the action (2.1) in terms of the vector components A µ = (A 0 , A i ).The resulting action is quadratic in these quantities.The time component A 0 is non-dynamical, and can be integrated out.We are left with three dynamical degrees of freedom: two transverse vector modes A λ (λ = ±), and one longitudinal mode A L .The dynamics of transverse and longitudinal modes decouple around a Friedmann-Robertson-Walker (FRW) background, which is sourced by a perfect fluid parametrized by a constant equation of state P = wρ.
The resulting vector equations of motion in Fourier space read (see Appendix A): Each mode depends on conformal time τ and co-moving momentum k = |k|, H is the Hubble rate evolving as H ′ = −3(1 + w) H 2 a/2.We introduce an effective time dependent mass squared as In writing equations (2.2), we indicate with primes derivatives along conformal time.In the course of this work, depending on what is more convenient, we work with conformal time τ or cosmic time t, related by dτ = dt/a(t).In the limit of vanishing non-minimal coupling α → 0, the system of equations (2.2) reduces to minimally coupled Proca set-up studied in [1].But, for α ̸ = 0, the non-minimal couplings with gravity can affect the vector longitudinal mode dynamics considerably.
In what follows, we neglect the back-reaction of vector fields on the geometry to analyze the vector evolution during inflation (w ≃ −1) and the subsequent radiation dominated universe (RDU) (w = 1/3).We assume that RDU starts after an efficient, nearly instantaneous reheating process.During inflation the longitudinal vector dynamics is qualitatively very similar to the original scenario [1], leading to an identical scale dependence of the longitudinal power spectrum spectrum that however exhibits a quantitatively different amplitude due to non-minimal coupling.On the other hand, a novel vector dynamics can appear at sub-horizon scales during RDU.In fact, due to the contribution proportional to α 2 which multiplies the k 2 in the first of equation of (2.2), the longitudinal mode develops a gradient instability, lasting a short period of time.Although brief, this phase is sufficient to potentially affect the vector evolution, and needs to be handled with care.We will do so in what comes next (see section 4.3), finding conditions to avoid dangerous instabilities.
The resulting cosmological dynamics is rich, and depends on the cosmological epoch one considers, as well as on the size of the vector co-moving wave-number k with respect to the remaining parameters of the system.We visually represent the various cosmological phases in Fig 1 .Each epoch is denoted by a capital letter in boldface.It has different consequences for the longitudinal vector evolution, which we study in the next sections.We make use both numerical production in the early universe through misalignment mechanisms.and analytical tools, inspired from the important papers [1,[8][9][10].The results will then be used to investigate the properties of dark matter consisting of longitudinal vector modes associated with the theory of eq.(2.1).

Dynamics of a light vector field during inflation
As we discuss below, the vector dynamics during inflation is qualitatively similar to the results in [1].However, for the range of parameters on which we are interested, we will find that the amplitude of the vector power spectrum is considerably affected by the non-minimal coupling to gravity.To study the dynamics of A L and A ± during inflation, we neglect sub-leading slow-roll corrections.We work with a constant Hubble rate H ≃ H I , and an equation of state parameter w = −1.The scale factor is given by a(τ ) = 1/(−H I τ ) in term of the (negative) conformal time −∞ < τ ≤ 0 during inflation.The equations of motion (EoM) satisfied by longitudinal and transverse modes are The effective mass reads and is constant during inflation.The structure of these equations is identical to the one studied in [1].Only the value (3.2) of the effective mass m eff,I changes, since it receives a contribution proportional to α 2 due to the non-minimal coupling with gravity.
We focus on a region in parameter space where the effects of the non-minimal coupling with gravity is sizeable.This will make a difference with respect to [1].For this reason, we satisfy the inequalities The first inequality ensures that the non-minimal couplings to gravity, controlled by α, plays an important role in our discussion.Furthermore, the second inequality places us in a light-vector regime during inflation, characterized by m eff,I ≪ H I .Hence the co-moving Hubble horizon (aH I ) −1 is always much smaller than the Compton co-moving horizon (am eff,I ) −1 of the vector field during inflation.In this regime, individual modes of the vector are guaranteed to satisfy the relativistic condition k ≫ am eff,I during inflation, at least until the time of horizon crossing where k = aH I .This condition is particularly important for understanding the time evolution of the longitudinal mode A L .

Dynamics of the transverse modes
We first focus on the transverse vector modes.We can be brief since their dynamics is the standard one of vector fields during inflation.In a de Sitter limit, the solution to the mode functions reducing to the standard Bunch Davies vacuum deep inside the horizon gives where H (1) ν is the Hankel function of first kind and we defined x = −kτ .In the light vector field regime m eff,I /H I → 0 we have ν ≃ 1/2, H 1/2 (x) ∝ e ix / √ x and so the transverse modes remain in their vacuum configuration as they are swept outside the horizon during inflation.This is expected since their equation of motion (EoM) (3.1) obeys the same equation of a scalar field conformally coupled to gravity in the small effective mass limit.Hence a sizeable contribution to the transverse mode population, if any, can only be due to some dynamics after the end of inflation.We address this topic in Appendix D, where we show that the post-inflationary dynamics does not manage to appreciably increase the amplitude of transverse vector fields.

Dynamics of the longitudinal modes
The dynamics of longitudinal modes is definitely more interesting.To understand the evolution of A L during inflation, we first notice that -due to the hierarchy m eff,I ≪ H I -all modes start their evolution in the relativistic regime k ≫ a m eff,I , where they remain until after horizon crossing.The dynamics of the longitudinal mode is then equivalent to a massless scalar field, up to an overall (important!) re-scaling.In fact, we define the canonical variable . Time evolution of two longitudinal modes for m eff,I /H I = 10 −2 from 4-folds before the horizon exit until the end of inflation where a/a end = 1.The dashed vertical lines denote the time of horizon crossing and the dotted vertical line refers to the time when the corresponding mode becomes non-relativistic.Note that the small scale mode labeled by red color stays relativistic throughout its evolution.
We rewrite the EoM (3.1) as in a k ≫ am eff,I limit.During inflation a ′′ /a = 2/τ 2 ; consequently, the solution of eq (3.6) matching the Bunch Davies vacuum deep inside the horizon (−kτ ≫ 1) is Inverting the field redefinition (3.5) by writing A L = (k/am eff,I )Q L , the solution for the longitudinal mode in the relativistic regime is where we defined the k-dependent super-horizon amplitude as Utilizing (3.8) and (3.9), we numerically solve the full mode equation (3.1) during inflation, by initializing each mode 4 e-folds before they exit the horizon.An example of evolution for two different modes initialized in the relativistic regime is shown in Figure 2. Shortly after horizon exit, all modes settles into a constant amplitude given by (3.9), irrespective of a given mode becoming non-relativistic or not.We represent this behaviour in the red curve of Figure 2, which shows the dynamics of a mode that never becomes non-relativistic during inflation.Hence, all longitudinal modes freeze out quickly on super-horizon scales, and settle to a constant amplitude given by The behaviour of A L we described can be also understood analytically, by studying the mode evolution in a piece-wise manner during inflation [1,9,10].To do so, it is convenient to pass from conformal to standard time, dτ = dt/a(t).Recall that m 2 eff = m 2 eff,I ≃ const., ∂ t m eff = 0 and w = −1 during inflation.The EoM (A.24) of the longitudinal modes, when expressed in terms of cosmic time, reads We now focus on the super-horizon relativistic and non-relativistic regimes, by taking the corresponding limits of (3.11).We follow the analytic methods developed in [1] and in [9,10].From now on, capital letters in boldface refer to different evolution phases as represented in Fig 1.

Phase (A): Super-horizon relativistic regime,
In the relativistic regime the solution to the longitudinal mode is given by (3.8).Shortly after horizon exit this solution can be approximated as which leads to a function containing a (growing) constant mode, plus a decaying mode ∝ a −2 : 3.2.2Phase (B): Super-horizon non-relativistic regime, H I > m eff,I ≫ k/a In this limit, we can approximate the EoM as (Hdt = da/a): We distinguish two situations.In the case of ultralight fields, the mass contribution in the right hand side can be neglected, m 2 eff,I /H 2 I = 0.The solution is given by For not-so-light vector fields however, the right hand side of eq.(3.14) plays a role.We generate a particular solution of (3.14) using iteratively the growing constant solution derived from the homogeneous equation.This solution has a secular growth suppressed by the square of the ratio m eff /H I , which leads to the final solution At first sight the secularly growing term might appear as an issue.However, since all the modes that become non-relativistic evolve through the relativistic regime during inflation, the contribution proportional to ln a never becomes dominant against the constant term.This conclusion is confirmed by the behavior of the blue curve in Figure 2 where the mode enters into the non-relativistic regime (at a = k/m eff,I ) at an intermediate time denoted by vertical dotted lines.Therefore, the only influence of the secular term is on the derivative of longitudinal mode, which changes its behavior from A ′ L (a) ∝ a −3 (see e.g.eq.(3.13)) to A ′ L (a) ∝ a −1 (see eq. (3.16)) as the mode evolve from relativistic to non-relativistic regime.However, in the ultra-light mass limit, eq.(3.15) is still valid in the non-relativistic regime and the dynamics of A L is captured by eqs.(3.13) to (3.15) as the modes evolve from relativistic to non-relativistic regime.(There are also modes that never becomes non-relativistic during inflation, the dynamics of them being described by the solution (3.13).)

The power spectrum of longitudinal modes at the end of inflation
Collecting the previous results, we conclude that super-horizon longitudinal vector modes acquire a power spectrum given by (3.17) We learn, as [1], that the power spectrum at the end of inflation is suppressed at large scales (k → 0) while it increases as k 2 towards large k (small scales).Such a characteristic slope of the longitudinal spectrum is one of the main findings of [1].Given that longitudinal vector perturbations act as iso-curvature modes, we then find negligible contributions from iso-curvature modes at large cosmic microwave background (CMB) scales.Hence, this scenario is safe from the point of view of CMB constraints on iso-curvature fluctuations.Recall that, given the inequality (3.3), we are interested on the regime of small bare mass m 2 ≪ α 2 H 2 I with respect to the contribution induced by non-minimal couplings to gravity.Consequently, during inflation, m 2 eff,I ≃ α 2 H 2 I , and the longitudinal vector spectrum is given by Notice that the amplitude of the power spectrum is independent from the vector mass and from the Hubble parameter during inflation.This remarkable feature, which is associated with the non-minimal coupling with gravity, makes a difference from the Proca case [1], and will be important in what follows when computing the dark matter abundance.Starting from these initial conditions, the behaviour of the power spectrum at later times depends on the post-inflationary evolution of the fluctuations, as they evolve from super-horizon to sub-horizon scales after inflation        ends.This is the argument of our next discussion, where we explore the consequences of the non-minimal coupling to gravity for the abundance of vector dark matter.

Dynamics of the longitudinal vector mode during radiation domination
During RDU the vector dynamics at small scales is considerably influenced by the non-minimal coupling with gravity.As we will learn, the latter induces a gradient instability in the longitudinal vector sector, which must be handled with care to maintain the system under control.As a consequence, to avoid catastrophic instabilities we will find constraints on the available parameter space.
In order to understand the effects of the non-minimal coupling α of the vector with gravity, we first study the system analytically.We make simplifying assumptions on the time-dependent profiles of the quantities involved, borrowing from previous studies.Our analytical findings will then be corroborated by a numerical analysis.Since part of the dynamics is identical to the Proca case, we can be brief where there is strong overlap with [1].After our analytical considerations, when turning to numerics, we focus most of our attention on the epoch (G), which is specific of our set-up, and where the dangerous gradient instability might develop.See Fig 3.

Analytical considerations
We set w = 1/3, m 2 eff = m 2 + α 2 H 2 .We recall that the Hubble parameter during RDU scales as where H I is the constant Hubble parameter during inflation, and a end the scale factor at the end of inflation.In writing eq. ( 4.1), we make the hypothesis of instantaneous transition from a (quasi) de Sitter inflationary stage to RDU.The evolution equation for the vector longitudinal mode, when expressed in terms of cosmic time t, is given in eq (A.24) and reads Recall that we are interested in the regime where non-minimal coupling dominates during inflation, see the inequality (3.3).The hierarchy αH I ≫ m implies that time dependent effective mass is decreasing in the initial stages of RDU, until it reaches its asymptotic value m at a critical time denoted by a c : 3) The critical scale factor a c sets an important time for the cosmological evolution of longitudinal modes and is defined by the relation During a fraction of the cosmic evolution a end ≤ a ≤ a c , corresponding to phase (G) in Fig 3, the coefficient of k 2 in eq. ( 4.2) becomes negative, signalling a gradient instability.This is a new effect of the non-minimal coupling with gravity, that we need to track with care.Notice that the duration of this phase depends on the size non-minimal coupling or more precisely on the ratio αH I /m whose amplitude will be important to keep the fluctuations that re-enter the horizon at this phase under control (see section 4.3).
Given the time dependence of the effective mass in (4.3), we parametrize its time evolution in a piece-wise manner as which constitutes a useful approximation for analytically handling the friction term in eq (4.2).At the light of the discussion above, we now study the analytic evolution of the longitudinal modes both for super and sub-horizon scales, during the different epochs represented in Fig 1.

Phases (E, F): super-horizon relativistic regime, H ≫ k/a ≫ m eff
In this regime, the EoM (4.2) is reduced to In writing the second line we assume that the factor multiplying the mass term in eq. ( 4.2) quickly reaches an order-one value: This limit can be justified by re-writing the second term inside the parentheses above as α 2 H 2 /m 2 eff ≃ (a c /a) 4 for a > a c (recall the scaling (4.1)).The solutions of eqs (4.6) in the super-horizon relativistic regime of RDU are then given by a end ≤ a ≤ a c , (E) Notice that there is a growing mode in epoch (E).However, the inflationary stage drives the time derivative of the longitudinal mode to extremely small values, and singles out the constant mode to be the dominant one on super-horizon scales.Therefore in the subsequent evolution during RDU, the solution in epoch (E) remains constant on super-horizon scales.We numerically confirmed this behaviour within the parameter space of interest.
Similarly, in the sub-horizon regime, the EoM of longitudinal mode can be reduced to (we switch back to conformal time for convenience) The solutions to these equations can be expressed as where we defined c s = 1/ √ 3 since the first line of (4.9) resembles the equation for a driven harmonic oscillator with an imaginary sound speed.The first exponentially growing solution in (4.10) is due to the aforementioned gradient instability for the short scale modes within the horizon (phase (G)) associated with the non-minimal coupling with gravity.See eq (4.2) and Fig 3 .Given its dangerous exponential amplification, we need to be sure to keep these short-scale modes under control.In fact, soon we will learn that this condition imposes limitations on the available parameter space.Notice also that after crossing phase (G) the modes enter in phase (H) (see Fig. 1) where their amplitude decreases.This behaviour will be important for the numerical considerations of Section 4.3.In this regime we neglect terms proportional to k in the damping contributions, as well as the mass term in (4.2) as compared to the damping term.We obtain the evolution equation with the following solution Similar to the analysis we made for epoch (F), the growing mode in eq.(4.12) does not take over, since the initial conditions inherited from the inflationary stage singles out the constant mode as the dominant one.In this regime, we can again neglect k 2 contributions in the damping term and the mass term, see (4.2): The solution of this equation is As the modes enter the non-relativistic regime, it starts oscillating with a decaying envelope a −1/2 (see Fig. 1 and 3).As first noticed in [1], for light vector fields, the energy density contained in the longitudinal mode starts to act like a dark matter in this phase.

Numerical Analysis
The analytical considerations of section 4.1 show that during RDU the longitudinal mode dynamics is quite rich and diverse.We now need a more careful analysis and a numerical treatment to follow the evolution of modes from inflation deep into RDU, taking into account the transitions among different epochs, and examining the role of the exponential gradient instability we anticipated in section 4.1.2.We intend to explore a region in parameter space where the gradient instability can be tamed, and where at the same time non-minimal couplings to gravity lead to sizeable, interesting effects for the phenomenology of dark matter.
For this purpose, we follow the cosmological evolution of a dimensionless quantity T L playing the role of a transfer function: where the initial condition A (0) L is given in (3.10) in terms of the value of the longitudinal mode at the onset of inflation: see section 3.2.To determine accurate initial conditions for the transfer function at the beginning of RDU, we evolve T L for a given mode k, from the time of its horizon crossing until the end of inflation.We use the numerical procedure described in Appendix B. This method allows us to determine the initial conditions at the beginning of RDU.We use them to follow the mode evolution during RDU, to finally find the longitudinal vector power spectrum at late times.We express the latter as in terms of the transfer function T L .We introduce the quantity k * corresponding to the wavenumber that re-enters the horizon at the time when the co-moving Hubble horizon is equal to the Compton horizon, as set by the bare mass (see Figs. 1 and 3):

.17)
As the modes evolve from super-horizon to sub-horizon scales, we study the evolution of T L during RDU.We find convenient for our numerical implementation to introduce a new time variable y, and a quantity x * related with the wave-number k as In terms of these variables, we use the EoM (2.2) for the longitudinal modes to obtain the evolution where H is the Hubble rate during RDU (4.1): H(y) = H I (y end /y) 2 with denoting the epoch signalling end of inflation, and the onset of RDU.Notice from the terms proportional to α in (4.19) that, the role of the non-minimal coupling with gravity can be important at early times during RDU due to the hierarchy αH I /m, but its role quickly becomes negligible as the universe expands and y ≫ y end .
Using accurate initial conditions we derive on the inflation side (see Appendix B), we numerically evolve (4.19) to obtain the power spectrum (4.16) evaluated at late times corresponding to y f = e 5 (e.g. 5 e-folds after the epoch H = m) and y f ≫ e 5 .We show in Figs. 4 and 5 the resulting (re-scaled) power spectrum y f x * |T L | 2 for representative parameter choices.The power spectrum exhibits a peak amplitude at an intermediate momentum k * = a * m corresponding to the scale that re-enters the horizon when H(a * ) = m.The k 2 rise of the spectrum from large towards small scales is inherited from the inflationary initial conditions (see section 3.3).It is due to the constant behaviour of the transfer function T L for long wavelength modes that mostly evolve at super-horizon scales.On the other hand, smaller scale modes with wave-numbers k ≳ k * start to experience sub-horizon evolution with decaying amplitudes (as in region (H)).At late times, this property causes their power amplitude to decrease inversely proportional to the wave-number, k −1 .It is interesting to obtain with no fine-tunings a spectrum of fluctuations with a peak structure, a task that is not easy in the context of adiabatic fluctuations and primordial black hole dark matter (see e.g.[38] for a recent review).
The profile for the longitudinal mode spectrum shown in Figs 4 and 5 shares the same qualitative behaviour of the Proca model (with α = 0) as introduced in [1].In fact, for not so large choices of αH I /m we adopt in Figures 4 and 5, the super-horizon dynamics of the vector modes around the peak k * is not much influenced from by the non-minimal coupling (see e.g.Figs. 1 and 2), as they go through phases (A)-(E)-(F) or (B)-(C).However, for larger values of the coupling α we increase the duration of phase (G) (see eq. (4.4)), hence the range of modes that enter the horizon during this phase.What is the fate of these modes?Can they contribute to the power spectrum in a dangerous way, especially for small scales that are not shown in Figure 4? What are the conditions we need to impose on the size of α in order to avoid catastrophic gradient instabilities?

4.3
The fate of the modes that re-enter the horizon during phase (G) We now discuss how to tame the instabilities of short-scale modes going through phase (G) of Fig 1 .(Recall the analytical discussion in section 4.1.2.)These modes are dangerous because they might considerably enhance the spectrum of vector modes at small scales, generating a second pronounced peak besides the one we met in Fig 4 .In this section we discuss how to avoid this fact to occur.
First of all, to understand the behavior of these modes, we compute the range of vector wave-numbers affected by this epoch.We evaluate the range of scales that re-enter the horizon between the end of inflation, a end , and the end of the dangerous phase (G), a c .This range is k end ≤ k ≤ k c .In terms of the coordinate y, the quantity a c is where we use eqs (4.4) and (4.20).Therefore, the duration of the phase (G) is longer for larger non-minimal coupling.In terms of the quantity y c , the mode k c that re-enters the horizon at a c can be expressed as Hence, the range of wave-numbers that re-enters the horizon in phase (G) shifts towards larger and larger scales with increasing non-minimal coupling.Introducing the wave-number that leaves the horizon right at then end of inflation as we find that the range of wave-numbers that re-enter the horizon during the phase (G) can be Figure 6.Time evolution of the maximally enhanced longitudinal modes (with colors) that re-enters the horizon in phase (G), corresponding to very small scales x * ≃ 2.8 × 10 7 (top-red) and x * ≃ 8.5 × 10 9 (bottom-orange) that exits the horizon just before the end of inflation, x * ≲ x end * .The part of the graphs highlighted by light blue indicates the time a mode spends in the unstable region as it enters inside the horizon.The scenario shown in the bottom panel has a larger αH I /m corresponding a longer (G) phase (see eq. (4.4)) such that a given short mode spends more time in the unstable region with an increasing amplitude.For α = 0 (black curves), the unstable era is not present and these short scale modes starts to decay as soon as they re-enter the horizon (the left most vertical dotted line) [1].
summarized in terms of the model parameters as This expression indicates that, choosing smaller and smaller values for m/H I , unstable modes will be pushed more towards UV, at a fixed non-minimal coupling.On the other hand at fixed vector field mass, the unstable modes have a tendency to shift towards IR by increasing the non-minimal coupling with gravity.
In Fig 6 we represent two different scenarios to show the evolution of modes that leave the horizon right before the end of inflation.These modes spend the most time within the unstable phase (G): therefore the gradient instability maximally enhances their size with respect to modes at larger scales.The unstable short-scale modes enhance their size in proportion to the duration of phase (G), as determined by eq.(4.4).
The danger -as mentioned above -is that the gradient instability acts on the small scale modes within the interval (4.24), leading to a second peak in the spectrum more pronounced than the first peak analyzed in Fig 4 .We wish to avoid this.Fortunately, the mode growth stops as Table 1.The peak amplitude of the re-scaled power spectrum for the maximally enhanced mode x * that re-enters the horizon during the gradient instability phase (G) for different parameter choices parametrized by the (α2 , m/H I ) pairs.See main text for additional explanations.
soon as they step inside phase (H), characterized by y > y c .(See section 4.1.2.)It is precisely this property that allows us to tame the instabilities.At this stage, indeed, they start to decay |T L | ∝ a −2 , spending a large amount of in the (H) era, until they become non-relativistic at very late times: The smaller the bare mass of the vector field is, the larger the amount of time these small scale modes spend inside the horizon in the decaying oscillatory phase (H).This can be visually realised by comparing the overall scales in the horizontal axis of Figure 6, and by noticing that for smaller bare mass,i.e.m/H I , individual modes tend to become non-relativistic later, see e.g.Fig Collecting the previous considerations, we can now estimate whether the size of the aforementioned second instability peak can be set under control.The parameter choice should ensure that the instability interval a end < a < a c remains small, and at the same time, we have a prolonged phase (H) following the unstable phase (G).We consider the examples studied in Figures 4 and 5.We extrapolate the behaviour of the transfer function to late times by fitting a decaying envelope to the maxima of oscillations of |T L | 2 , as shown in Fig. 6.In this way, we compute the re-scaled power spectrum y x 2 * |T L | 2 at y nr (defined in eq.(4.25)) for each parameter choices shown in Figures 4 and 5.We present these results in Table 1, which informs us that the second peak generated at small scales by the gradient instability is well smaller than the first peak at k * .Similar considerations and numerical checks can be done for other pairs of parameters (α, m/H I ), so to circumscribe the region of parameters which avoid large effects from the gradient instability 2 .We will return to this point in the next Section, when discussing the final abundance of longitudinal vector modes.

Abundance of longitudinal vector modes
We now collect the previous results, and compute the final abundance of longitudinal vector modes.We show that they can constitute the entirety of dark matter, for a broad range of vector masses.The final dark matter abundance depends on two parameters only, m and α, making this scenario very economical.
First of all, we compute the longitudinal-mode contribution to the universe energy density.The presence of the non-minimal coupling coupling complicates the computation of the the energy density contained in the longitudinal modes.However, as we show in Appendix C, the contributions due to the non-minimal coupling in the energy-momentum tensor (EMT) becomes negligible when all the modes associated with the peak in the power spectrum becomes nonrelativistic.In particular, the corrections induced by the non-minimal interactions appear at order α 2 H 2 O(A ′2 , A 2 ).Their size is then negligible compared to terms of order m 2 O(A ′2 , A 2 ) for a ≫ a c , as can be realized by noticing that α 2 H 2 is comparable to m 2 at around a ∼ a c , see eq (4.4).This situation can also be interpreted as a consequence of the null energy condition, which ensures that Hubble rate is always decreasing.The effects introduced by to non-minimal coupling are maximal during inflation (where H decreases very slowly), but subsequently their size diminishes as the Hubble rate decays fast during the post-inflationary era.Therefore, at late times the energy density contained in the fluctuations of the longitudinal modes is [1], When all the modes we consider are inside the horizon, the fluctuations are virialized so that all quantities contribute equally to the energy density in (5.1).Parametrizing these contributions in terms of the power spectrum, and using the time variable y = a/a * adopted in the previous sections, we can therefore re-write ρ A L in terms of the transfer function (4.15) of the longitudinal modes as As can be inferred from the evolution of longitudinal modes deep in phase (D) (see section 4.1.4),the integrand in this expression is time-invariant, as far as we evaluate the power spectrum at epochs when all the modes x * are in the non-relativistic regime.Therefore, considering the peaked profile for the power spectrum (as in Fig. 4) in (5.2), and using the energy density in the longitudinal modes today gives where we assume m 2 eff,I ≃ α 2 H 2 I in the αH I ≫ m regime we focus on.Taking into account the observed dark matter abundance ρ dm (t 0 ) = 1.26 × 10 −6 GeV/cm 3 , we adopt the central value for the red-shift at equality z eq = 3402 [2].We finally obtain the fraction of longitudinal modes in dark matter density today as (5.5) of transverse modes remain subdominant).This scenario is very economical; nevertheless, the expression for the dark matter abundance (5.5) leaves some degeneracy because the parameters α and m can be chosen arbitrarily, at least within the mass interval (5.6).In the next section we discuss further phenomenological consequences of our framework, aimed to distinguish the effects of the two parameters (α, m) on observable quantities, which might allow us to obtain independent measurements of their size.

Phenomenological considerations
We learned in the previous section that our version of the scenario introduced in [1] realizes dark matter in the form of longitudinal modes of a massive vector non-minimally coupled with gravity.
The dark matter abundance depends on two parameters only, making the model very economical.
In order to provide the totality of dark matter, the non-minimal coupling α -as appearing in action (2.1) -and the vector mass m are related by (see eq.Hence, for relatively light vector fields (recall the allowed mass range of eq.(5.6)), the nonminimal coupling of vector to gravity of in action (2.1) acquires very small values3 .At the same time, the vector mass (for high-scale inflation) can be chosen within the interval (5.6), leaving some degeneracy in parameter space for realizing scenarios where the totality of dark matter is constituted by longitudinal modes.
Besides cosmological consequences of dark photon dark matter, we now ask whether there are any further experimental avenues to independently probe the parameter m appearing in action (2.1).Let us consider scenarios whose tree level actions do not contain gauge interactions between visible and dark sectors, besides gravity.Nevertheless, loop effects inevitably induce portal couplings with the Standard Model, for example through a kinetic mixing with the SM photon A (γ) µ [42][43][44][45] F µν (γ) is the field strength for the SM photon field, while F µν the dark vector field strength.This gauge-invariant dimension-four operator is controlled by the dimensionless parameter ε, which induces milli-charged couplings of dark photons to ordinary matters.Constraints exist on ε, which depend on the dark photon mass (see e.g.[16] for a comprehensive review).Even for set-up where there are no tree-level gauge couplings between visible and dark sectors, it has been demonstrated (see [46]) that loop effects mediated by gravity alone generate kinetic mixing of the form (6.2).The resulting values of the parameter ε are very tiny: ε ≃ 10 −13 ; nevertheless such small values can be probed for dark vector masses in the eV − keV range [47] (see also [48,49]), which is within the mass interval (5.6) we found.Hence, direct experiments testing milli-charged couplings with the dark sector can independently probe the range of vector masses which can be realized through non-minimal couplings of the dark vector to gravity.
Another option is to make use of gravitational wave experiments.As mentioned in the introduction, it is conceivable to use gravitational wave detectors to directly probe effects of spacetime deformations associated with vector degrees of freedom.See e.g.[28] who uses a vectorial version of the so-called Khmelnitsky-Rubakov effect [50] in the context of pulsar timing-array experiments.
Alternatively, one can study primordial black holes [51,52] or induced stochastic gravitational wave backgrounds (SGWB) [53] produced at second order in perturbations by large-size isocurvature modes 4 .This is particularly interesting for vector mass ranges in the eV − keV interval, as we are going to discuss.The subject of scalar-induced SGWB by adiabatic fluctuations, first investigated in [54][55][56][57][58][59][60][61][62][63], is very well developed by now (see the review [32]).Much less studied are SGWB induced by iso-curvature modes -the case relevant for us since we study vector longitudinal modes corresponding to iso-curvature fluctuations.A notable exception is [53], where the authors carry on a detailed analysis of how an enhanced spectrum of iso-curvature modes at super-horizon scales can induce a stochastic GW background after re-entering the horizon in a radiation-dominated universe.This is possible because the early iso-curvature modes can get converted into curvature fluctuations upon horizon crossing, which in turn source GWs.Therefore, similar to the well studied adiabatic case an enhanced, peaked spectrum of longitudinal iso-curvature modes (with a maximum at scale k * ) could lead to a stochastic GW background peaked at scales of the same order of magnitude.Since in our case the scale of the scalar peak is proportional to the vector mass k * = a * m (see eq. ( 5.3)), converting this quantity to the frequency domain (f = k/(2π)), we expect to generate a SGWB enhanced at a peak frequency Hz. (6.3) A choice of vector dark matter mass around the eV − keV range -as can be realized in our scenario (5.6) -can then enhance the SGWB signal in the milli-Hertz band, making it potentially detectable by the LISA interferometer [64].We realize that the formulas developed in [53] can not be directly used in our context, without some extra work.The derivations of [53] take into due account the back-reaction of iso-curvature fluctuations to gravity, while in this work we do not consider the back-reaction of vector fields in the geometry.Moreover, the formulas of [53] assume that the primordial iso-curvature spectrum is evaluated in a super-horizon regime, while the peaked shape of the spectrum presented in Figs 4 and 5 include also contributions from modes evolving at sub-horizon scales.Nevertheless, we expect that a future, more careful analysis of energy fluctuations associated with iso-curvature modes (including back-reaction in the geometry) again will lead to a spectrum exhibiting a peaked structure, with a peak around the scale k * = a * m found in the previous section.The analysis of iso-curvature density fluctuations carried on in [1] indicates that the corresponding spectrum grows as (k/k * ) 3 from large towards small scales, up to a maximum at around k ≃ k * , to then decrease as (k/k * ) −1 for smaller scales k ≥ k * .We take such a spectrum profile as representative The GW density spectrum induced at second order by longitudinal vector iso-curvature fluctuations.See the main text, in particular the discussion after eq. ( 6.3) for a description of the set-up.
Ansatz of longitudinal vector density fluctuations, and denote A s as the maximum of the spectrum at k = k * .We then numerically compute the integrals of [53], using their kernels and obtain the profile for the GW density Ω gw as a function of f /f * as represented in Fig 9 .The spectrum has a peaked structure, with a relatively broad maximum at scales corresponding to the frequency region around the f * of eq. ( 6.3) detectable with LISA.In the IR scales shown, it has a relatively mild slope as compared to GWs induced by adiabatic fluctuations (Ω gw ∼ f 3 ) with a f <2 behavior, while in the UV it decays as fast as f −6 .It would be interesting to perform a more careful analysis on the shape and amplitude A s in a physically realistic set-up as it could potentially provide us information about the couplings and parameters in our proposal (2.1).We leave this important task to future studies.

Conclusions and Outlook
We proposed a scenario of dark matter in the form of longitudinal massive vector bosons, extending the work [1] by including non-minimal couplings of the vector with gravity.Our set-up, besides being well motivated from an effective field theory perspective, is very economical since the final dark matter abundance depends on two parameters only (see eq. (5.5)): the vector mass and its non-minimal coupling with gravity.We determined and analysed in detail the rich dynamics of vector fluctuations during inflation and radiation domination.In particular, we showed that during inflation, quantum fluctuations of longitudinal modes obtain a large super-horizon power, which depends on the inverse strength of the non-minimal coupling with gravity.The latter induces gradient instabilities on vector field dynamics during radiation dominated era, which can nevertheless be tamed.For weak non-minimal coupling with gravity, the inflationary era thus sets suitable initial conditions for the longitudinal component of vector field.The resulting post-inflationary vector energy density can account for the totality of dark matter, for a wide range of vector masses (5.6).We emphasize that the set-up we described is also minimal in the sense that it assumes standard post-inflationary history following the inflationary epoch.It would also be interesting to explore variations of our proposal on dark matter production by utilizing non-trivial post-inflationary histories beyond the instantaneous reheating approximation we undertake in this work, as done in [9,65].
Due to non-minimal coupling, the resulting abundance of dark matter leads to a wider range of allowed vector masses with respect to [1].Interestingly, such mass range can be experimentally probed by means of induced gravitational waves in the LISA frequency band.We touched this argument in section 6 in a phenomenological way: it certainly deserves further investigations.
Given the absence of direct evidence of dark matter besides its gravitational interactions with visible matter, it is perhaps a good idea to focus on gravitational experiments as dark matter detectors.Gravitational wave experiments can provide new ways to probe light and feebly interacting dark matter (see e.g.[66] for a review).It will important in the near future to further explore the opportunities they can provide.

A Equations of motion and mode functions of the gauge field
In this appendix, we derive the equations of motion of the gauge fields following two different methods: i) first by expanding the action (2.1) and decomposing the vector field into its temporal and spatial components ii) by following a covariant approach we will discuss below.We begin with the former by expanding the kinetic term of the gauge fields as where we decomposed the non-vanishing component of the field strength tensor as 1), we re-write the action (2.1) as (A.2) Notice that the action does not contain time derivatives of the temporal component of the gauge field which informs us that A 0 is a Lagrange multiplier whose equation of motion indicates a constraint of the system.To derive the latter, we require the knowledge of the Einstein tensor at Using (A.24), we present an analytical understanding of mode evolution during the cosmic history in sections 3 and 4.

B Mode evolution of A L during inflation: Numerical Procedure
In order to accurately capture the cosmological evolution of the individual A L modes in the post-inflationary universe, we need to determine their initial conditions at the beginning of RDU when they are outside the co-moving horizon.For this purpose, we solve for the transfer function T L (4.15) using the EoM (3.1) during inflation.Working with time variable y = a/a * and noting x ≡ −kτ = −kτ * (τ /τ * ) = x * y −1 and x * ≡ k/k * , the dynamics of the transfer function during inflation can be characterized by the following equation: Focusing on sizeable non-minimal coupling regime αH I /m ≫ 1, we work with sufficiently light longitudinal modes where m 2 eff,I = m 2 + α 2 H 2 I ≃ α 2 H 2 I ≪ H 2 I implying the strength of the non-minimal coupling to be small, α ≪ 1.In this case, we can initialize the modes at horizon exit by safely assuming they are in the relativistic regime parametrized by the solution (3.8).In terms of the new time variable y, we therefore adopt the following initial conditions during inflation, where for each mode we have y in = x * y 2 end with y end = m/H I denoting the end of inflation/beginning of RDU.To determine accurate final conditions that can be used for the post-inflationary evolution of the longitudinal modes, we solve (B.1) for a grid of x * = k/k * values using (B.2) for various α and m/H I ≪ 1 choices.As we mentioned before, we do so in a regime where effects introduced by the non-minimal coupling is sufficiently large during inflation, i.e.
where the final inequality is required to keep the longitudinal mode light during inflation.This procedure provides us with a pair of {T L (y end , x * ), T ′ L (y end , x * )} for all x * = k/k * values that can be used as an input for the post-inflationary evolution we discuss in Section 4.2.

C Energy density of the longitudinal modes
In this appendix, we focus our attention to the energy density and pressure of the vector fields.For this purpose, we note that the energy-momentum tensor (EMT) of the action (2.1) can be written as: where in passing from second to third and from third to last line we have used the symmetry properties of the Riemann tensor.
To derive the contribution to T + g µν g ρσ − δ ρ µ δ σ ν g αβ + δ α ν δ σ µ g βρ + δ β µ δ σ ν g αρ − δ α ν δ β µ g σρ − g µν g αρ g βσ ∇ ρ ∇ σ (A α A β ) , where Noting that T 00 = ρ, we first expand the non-derivative terms explicitly using the components of Einstein and Ricci tensor in the first line of (C.9) together with the higher derivative terms in the second line in terms of the temporal and spatial component of the vector field as (C.11) Utilizing the definition of the Riemann tensor (A.13), the second term above read as (C.12) Finally unpacking the last term in (C.11) according to the definition (C.10) and after a bit of algebra, we compile our findings to obtain the energy density due to non-minimal coupling as On the other hand, energy density in the absence of non-minimal coupling can be obtained from (C.2) as To obtain the total energy density of the longitudinal modes from (C.13) and (C.14), we first notice that the curly of the vector field ⃗ ∇ × A does not contain the longitudinal mode which is a direct consequence of the property of the polarization vector ϵ (L) i noted in the second line of (A.8).On the other hand going into the Fourier space using the eqs.(A.21)-(A.23),one can show that the combination of the last four terms in the first line, as well as the first term in the last line of (C.13) vanishes.Furthermore, as far as the longitudinal mode is concerned, the envelope of the last term in the non-minimal EMT behaves as ∼ H Ȧk A k ∼ H 2 A k A k at late times, as can be inferred from eq. (4.14).Therefore, in the late time limit the energy density of the longitudinal modes that stem from the non-minimal coupling behave as ρ (α) Comparing the overall scaling of these terms with the last terms in (C.14), we conclude that the corrections to the energy density due to non-minimal interactions are much less important in the late time limit at which we are interested to evaluate the abundance of the longitudinal modes.In particular, for a ≫ a * ≫ a c , we have m ≫ αH because αH is already comparable to m at the end of the phase (G) in the RDU (see Fig. 3).Therefore, we can safely adopt the standard energy density (C.

D Abundance of transverse vector modes
In the main text we ignored the contribution to the abundance of the vector field from the transverse modes, relying on the fact that they stay in their vacuum configuration during inflation in the small mass limit, m/H I → 0. However, post-inflationary evolution of the transverse modes could potentially provide significant contribution to the abundance since they obtain a tachyonic mass following inflationary phase in the RDU era for the sizeable non-minimal coupling regime αH I /m ≫ 1 we are operating.This could be seen from eq. (2.2), where the mass of the transverse mode in the post-inflationary era behave as where the end of tachyonic regime with respect to the end of inflation can be identified as Notice that the end of the tachyonic era almost coincides with the end of phase (G) (4.4) during which the longitudinal modes have gradient type instability.This similarity notwithstanding, the transverse modes that are affected by the tachyonic instability are different as compared to the former.To understand this, we rewrite the the mode equation of the transverse modes in terms of the variables we are familiar from the the previous sections as where we again note H(y) = H I (y end /y) 2 during RDU.From the effective dispersion relation of the expression above, we infer that the shortest modes that could be affected by the presence of the tachyonic regime satisfy Therefore, as compared to the longitudinal modes, relatively large scale transverse modes could be potentially amplified during the tachyonic regime a end < a < a tc .
To understand if the transverse modes can be amplified in the tachyonic era, we followed a similar approach we carry in Appendix B to first generate initial conditions for the modes in the post-inflationary era by evolving A T during inflation from deep inside the horizon until the end of inflation.We then use these initial condition to numerically evolve the mode equation of the transverse modes (D.3) during RDU until sufficiently late times.In this way, we found that transverse modes that could potentially exhibit instability do not behave any interesting way during the tachyonic era as compared with the vanishing non-minimal coupling case α → 0 for which the tachyonic phase is absent.We will not present these null results here, however we refer the interested to the Mathematica notebook that lead us draw this conclusion.As can be inferred < l a t e x i t s h a 1 _ b a s e 6 4 = " o d y D 1 b l j D I 6 0 e y t i I y x w s T Y r E o 2 B H / 5 5 V X S q l X 9 i + r 5 f a 1 S v 8 7 j K M I x n M A Z + H A J d b i F B j S B Q A L P 8 A p v z p P z 4 r w 7 H 4 v W g p P P H M E f O J 8 / 4 P m T O Q = = < / l a t e x i t > Inflation < l a t e x i t s h a 1 _ b a s e 6 4 = " L F 3 / N 0 d c + 2 2 b 4 h J M F M I x r N w J I k I T H 5 I z U S J 1 w c k 8 e y T N 5 C R 6 C p + A 1 e J t G C 0 E + s 0 V + I H j / A v S G n D M = < / l a t e x i t > comoving scales < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 f 9 I O J L F t W p m k f D O r p m L j 6 o 7 x h A = " > A A A B 8 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I x g Q 2 S 5 i d n S R D 5 r H M z A p h y W d 4 8 a C w D T W 4 h z r 4 Q I H D M 7 z C m y O d F + f d + Z i 3 F p x 8 5 h D + w P n 8 A a W r j p g = < / l a t e x i t > a ⇤ < l a t e x i t s h a 1 _ b a s e 6 4 = " y H i j U 1 / D l C i 3 J R 5 K S 4 + b E p i W X z A = " > A A A B 8 X i c d V D L S g M x F M 3 U V 6 2 v q k s 3 w S L U h S V T t L a 7 o h u X F e w D 2 7 F k 0 k w b m m S G J C O U o X / h x o U i b v 0 b d / 6 N m b a C i h 6 4 c D j n X u 6 9 x 4 8 4 0 w a h D y e z t L y y u p Z d z 2 1 s b p n 1 6 e F + s U i j i w 4 A I e g C F x w D u r g C j R A E x A g w Q N 4 A s + O d h 6 d F + d 1 3 p p x F j P 7 4 A e c t 0 9 b 7 5 A Y < / l a t e x i t > (am) 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " a 5 c D L J e h c 3 g w 7 h F d 4 8 7 b 1 4 7 9 7 H v L X g 5 T O H 8 A f e 5 w + o J J B U < / l a t e x i t > a c < l a t e x i t s h a 1 _ b a s e 6 4 = " + a l c l Z c V m W o 5 u l I V z w V t x z N 0 5 o 4 T w t e n 8 H / S y D t u y S l e F r L V s 3 k c a b A P D k A O u O A U V M E F q I E 6 I O A O P I A n 8 G z d W 4 / W i / U 6 a 0 1 Z 8 5 l d 8 A P W 2 y d 1 O p S f < / l a t e x i t > (am e↵ ) 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " u / O V M N M F J M H Y r R 1 V L L Y n o 1 d 4 4 P g = " > A A A B 8 X i c b V D L S g N B E O z 1 G e M r 6 t H L Y B A 8 h V 3 x d Q x 6 8 R j B P D B Z w u y k N x k y O 7 v O z A p h y V 9 4 8 a C I V / / G m 3 / j J N m D J h Y 0 F F X d d H c F i e D a u O 6 3 s 7 S 8 s r q 2 X t g o b m 5 t 7 + y W 9 v Y b O k 4 V w z q L R a x a A d U o u M S 6 4 U Z g K 1 F I o 0 B g M x j e T P z m E y r N Y 3 l v R g n 6 E e 1 L H n J G j Z U e a D f r q I j g 4 7 h b K r s V d w q y S L y c l j a g i j J j Q y r a E L z 5 l x d J 4 7 T i X V T O 78 7 K 1 e s 8 j g I c w h G c g A e X U I V b q E E d G E h 4 h l d 4 c 7 T z 4 r w 7 H 7 P W J S e f O Y A / c D 5 / A H 7 h k N E = < / l a t e x i t > a eq < l a t e x i t s h a 1 _ b a s e 6 4 = " p I N c O Z b Z W / M S D S o s x r 1 x 7 f F B D r U = " > A A A B 6 H i c d V B N S w M x E J 2 t X 7 V + V T 1 6 C R b B U 8 k W r e 2 t K I j H F u w H t E v J p t k 2 N p t d k q x Q S n + B F w + K e P U n e f P f m G 0 r q O i D g c d 7 M 8 z M 8 2 P B t c H 4 w 8 m s r K 6 t b 2 Q 3 c 1 v b O 7 t 7 + f 2 D l o 4 S R V m T R i J S H Z 9 o J r h k T c O N Y J 1 Y M R L 6 g r X 9 8 V X q t + + Z 0 j y S t 2 Y S M y 8 k Q 8 k D T o m x U u O 6 n y / g I r Y o l 1 F K 3 A p 2 L a l W K 6 V S F b l z C + M C L F H v 5 9 9 7 g 4 g m I 2 t i I 6 I I t T Y b H I 2 h K 9 P 0 f + k V S q 6 5 e J 5 4 6 x Q u 1 z G k Y U j O I Z T c O E C a n A D d W g C B Q Y P 8 A T P z p 3 z 6 L w 4 r 4 v W j L O c O Y Q f c N 4 + A f T 2 j R E = < / l a t e x i t > F < l a t e x i t s h a 1 _ b a s e 6 4 = " 1 I + l b x h L T Q 5 / D v q y e E P q v y + O l 7 s = " > A A A B 6 H i c d V B N S w M x E J 2 t X 7 V + V T 1 6 C R b B U 8 k W r e 2 t 6 s V j C / Y D 2 q V k 0 2 w b m 8 0 u S V Y o p b / A i w d F v P q T v P l v z L Y V V P T B w O O 9 G W b m + b H g 2 m D 8 4 W R W V t f W N 7 K b u a 3 t n d 2 9 / P 5 B S 0 e J o q x J I x G p j k 8 0 E 1 y y p u F G s E 6 s G A l 9 w d r + + D r 1 2 / d M a R 7 J W z O J m R e S o e Q B p 8 R Y q X H Z z x d w E V u U y y g l b g W 7 l l S r l V K p i t y 5 h X E B l q j 3 8 + + 9 Q U S T k E l D B d G 6 6 + L Y e F O i D K e C z X K 9 R L O Y 0 D E Z s q 6 l k o R M e 9 P 5 o T N 0 Y p U B C i J l S x o 0 V 7 9 P T E m o 9 S T 0 b W d I z E j / 9 l L x L 6 + b m K D i T b m M E 8 M k X S w K E o F M h N K v 0 Y A r R o 2 Y W E K o 4 v Z W R E d E E W p s N j k b w t e n 6 H / S K h X d c v G 8 c V a o X S 3 j y M I R H M M p u H A B N b i B O j S B A o M H e I J n 5 8 5 5 d F 6 c 1 0 V r x l n O H M I P O G + f 7 W K N D A = = < / l a t e x i t > A < l a t e x i t s h a 1 _ b a s e 6 4 = " k Q l F x + e W B x n d B G a 0 L / r + D l j j W R 4 = " > A A A B 6 H i c d V B N S w M x E J 2 t X 7 V + V T 1 6 C R b B U 8 k W r e 2 t K I L H F u w H t E v J p t k 2 N p t d k q x Q S n + B F w + K e P U n e f P f m G 0 r q O i D g c d 7 M 8 z M 8 2 P B t c H 4 w 8 m s r K 6 t b 2 Q 3 c 1 v b O 7 t 7 + f 2 D l o 4 S R V m T R i J S H Z 9 o J r h k T c O N Y J 1 Y M R L 6 g r X 9 8 V X q t + + Z 0 j y S t 2 Y S M y 8 k Q 8 k D T o m x U u O 6 n y / g I r Y o l 1 F K 3 A p 2 L a l W K 6 V S F b l z C + M C L F H v 5 9 9 7 g 4 g m I 2 t i I 6 I I t T Y b H I 2 h K 9 P 0 f + k V S q 6 5 e J 5 4 6 x Q u 1 z G k Y U j O I Z T c O E C a n A D d W g C B Q Y P 8 A T P z p 3 z 6 L w 4 r 4 v W j L O c O Y Q f c N 4 + A f N y j R A = < / l a t e x i t > E < l a t e x i t s h a 1 _ b a s e 6 4 = " j a p 0 b 1 w d o O m 7 / v c K R s h 4 k 6 g + q y

Figure 1 .
Figure 1.A schematic diagram that summarizes the evolution of the longitudinal mode A L at different stages of the cosmic history.In the main text we will repeatedly refer to this figure, explaining the dynamics in each phase denoted by a capital letter.
< l a t e x i t s h a 1 _ b a s e 6 4 = " o d y D 1 b l j D I 6 0 x a C 0 4 + c w x / I H z + Q P V E 4 6 n < / l a t e x i t > ln(a)< l a t e x i t s h a 1 _ b a s e 6 4 = " i c 8 S q A 2 G X gV 2 K s f s + Y N Z h 3 4 a b T k = " > A A A B + H i c b V A 9 S w N B E J 2 L X z F + 5 N T S Z j E I V u E u i F p Y B G y 0 i 2 A + I A l h b 7 O X L N n d O 3 b 3 h H j k l 9 h Y K G L r T 7 H z 3 7 i 5 X K G J D w Y e 7 8 0 w M y + I O d P G 8 7 6 d w t r 6 x u Z W c b u 0 s 7 u 3 X 3 Y P D l s 6 S h S h T R L x S H U C r C l n k j Y N M 5 x 2 Y k W x C D h t B 5 O b u d 9 + p E q z S D 6 Y a U z 7 A o 8 k C x n B x k o D t 5 z 2 l E B 3 M u S Z M B u 4 F a / q Z U C r x M 9 J B X I 0 B u 5 X b x i R R F B p C M d a d 3 0 v N v 0 U K 8 M I p 7 N S L 9 E 0 x m S C R 7 R r q c S C 6 n 6 a H T 5 D p 1 Y Z o j B S t q R B m f p 7 I s V C 6 6 k I b K f A Z q y X v b n 4 n 9 d N T H j V T 5 m M E 0 M l W S w K E 4 5 M h O Y p o C F T l B g + t Q Q T xe y t i I y x w s T Y r E o 2 B H / 5 5 V X S q l X 9 i + r 5 f a 1 S v 8 7 j K M I x n M A Z + H A J d b i F B j S B Q A L P 8 A p v z p P z 4 r w 7 H 4 v W g p P P H M E f O J 8 / 4 P m T O Q = = < / l a t e x i t > Inflation < l a t e x i t s h a 1 _ b a s e 6 4 = " L F 3 / N 0 d c + 2 2 b 4 h J M F M I x r N w J I k I = " > A A A B 7 n i c b V B N T w I x E J 3 F L 8 Q v 1 K O X R m L i i e w a g h 6 J e v C I x g U S 2 J B u 6 U J D 2 9 2 0 X R O y 4 U d 4 8 a A x X v 0 9 3 v w 3 F t i D g i + Z 5 O W 9 m c z M C x P O t H H d b 6 e w t r 6 x u V X c L u 3 s 7 u 0 f l A z s W g t O P n M M f y B 8 / k D i m + P D w = = < / l a t e x i t > RDU < l a t e x i t s h a 1 _ b a s e 6 4 = " 8 u o n F h P 8 K n H k F e m B H I d t q 6 H 2 z p 8 = " > AA A B / H i c d V D L S g M x F M 3 U V 6 2 v 0 S 7 d B I t Q Q Y e M t L X u i m 6 6 r G A f 0 N a S y W T a 0 M z D J C O U o f 6 K G x e K u P V D 3 P k 3 p u 0 I K n o u F w 7 n 3 E t u j h N x J h V C H 0 Z m a X l l d S 2 7 n t v Y 3 N r e M X f 3 W j K M B a F N E v J Q d B w s K W c B b S q m O O 1 E g m L f 4 b T t j C 9 n f v u O C s n C 4 F p N I t r 3 8 T B g H i N Y a W l g 5 n v H u m 5 j 7 B Y x r B / d J C f 2 d G A W k F U + 1 6 h C Z F U Q Q n Y 5 J a g M b Q v N U Q A p G g P z ve e G J P Z p o A j H U n Z t F K l + g o V i h N N p r h d L G m E y x k P a 1 T T A P p X 9 Z H 7 8 F B 5 q x Y V e K H Q H C s 7 V 7 x s J 9 q W c + I 6 e 9 L E a y d / e T P z L 6 8 b K q w b N w b j 8 a L 8 b o Y z R j p T h 7 8 g P H 2 C d N N k 5 8 = < / l a t e x i t > (aH) 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " n F 2 P T M W j L f h H 6 X d r E 5 y r z 5 4 + a D 6 a P d N T x V o J F L 5 l w r C h P s Z M y i 4 B J G x X b q I G H 8 h v W h 5 a l m C l w n m y w z o r t e 6 d K e s f 5 o p B P 1 + 0T G l H N D F f u k Y n j t f n t j 8 T + v l W L v q J M J n a Q I m k 8 f 6 q W S o q H j Z m h X W O A o h 5 4 w b o X / K + X X z D K O v r + i L y H 6 v f J f 0 q h W o o P K / n m 1 f H y S 1 7 F I t s k O 2 S MR O S T H 5 I z U S J 1 w c k 8 e y T N 5 C R 6 C p + A 1 e J t G C 0 E + s 0 V + I H j / A v S G n D M = < / l a t e x i t > comoving scales < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 f 9 I O J L F t W p m k f D O r p m L j 6 o 7 x h A = " > A A A B 8 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I x g Q 2 S 5 i d n S R D 5 r H M z A p h y W d 4 8 a C 1 p y E z 9 P Z F h p P U 4 C m 1 n h G a o F 7 2 p + J / X T k 3 / J s i 4 T F L D J J 0 v 6 q e C m J h M P y c 9 r h g1 Y m w J U s X t r Y Q O U S E 1 N p + S D c F b f H m Z N M 6 r 3 l X 1 8 u G i U r v N 4 y j C E R z D K X h w D T W 4 h z r 4 Q I H D M 7 z C m y O d F + f d + Z i 3 F p x8 5 h D + w P n 8 A a W r j p g = < / l a t e x i t > a ⇤ < l a t e x i t s h a 1 _ b a s e 6 4 = " y H i j U 1 / D l C i 3 J R 5 K S 4 + b E p i W X z A = " > A A A B 8 X i c d V D L S g M x F M 3 U V 6 2 v q k s 3 w S L U h S V T t L a 7 o h u X F e w D 2 7 F k 0 k w b m m S G J C O U o X / h x o U i b v 0 b d / 6 N m b a C i h 6 4 c D j n X u 6 9 x 4 8 4 0 w a h D y e z t L y y u p Z d z 2 1 s b m 3 v 5 H f 3 W j q M F a F N E v J Q d X y s K W e S N g 0 z n H Y i R b H w O W 3 7 4 8 v U b 9 9 T p V k o b 8 w k o p 7 A Q 8 k C R r C x 0 m 0 R Q 3 F 8 l 5 y 4 0 3 6 + g E r I o p n 1 6 e F + s U i j i w 4 A I e g C F x w D u r g C j R A E x A g w Q N 4 A s + O d h 6 d F + d 1 3 p p x F j P 7 4 A e c t 0 9 b 7 5 A Y < / l a t e x i t > (am) 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " a 5 c D L J e h c 3 g w c / R A y Z v H y g T A R + o = " > A A A B 8 H i c bV D J S g N B E K 2 J W 4 x b 1 K O X x i B 4 C j P i d g x 6 8 R j B L J I M o a f T k z T p Z e j u E c K Q r / D i Q R G v f o 4 3 / 8 Z O M g d N f F D w e K + K q n p R w p m x v v / t F V Z W 1 9 Y 3 i p u l r e 2 d 3 b 3 y / k H T q F Q T 2 i C K K 9 2 O s K G c S d q w z H L a T j T F I u K 0 F Y 1 u p 3 7 r i W r D l H y w 4 4 S G A g 8 k i x n B 1 k m P u J d 1 t U B k 0 i t X / K o / A 1 o m Q U 4 q k K P e K 3 9 1 + 4 q k g k p L O D a m E / i J D T O s L S O c T k r d 1 N A E k x E e 0 I 6 j E g t q w m x 2 8 A S d O K W P Y q V d S Y t m 6 u + J D A t j x i J y n Q L b o V n 0 p u J / X i e 1 8 X W Y M Z m k l k o y X x S n H F m F p t + j P t O U W D 5 2 B B P N 3 K 2 I D L H G x L q M S i 6 E Y P H l Z d I 8 q w a X 1 Y v 7 8 0 r t J o + j C E d w D K c Q w B X U 4 A 7 q 0 A A C A p 7 h F d 4 8 7 b 1 4 7 9 7 H v L X g 5 T O H 8 A f e 5 w + o J J B U < / l a t e x i t > a c < l a t e x i t s h a 1 _ b a s e 6 4 = " + a l c l Z c V m W o 5 u l I V z w V t x z N 0 5 o 4 = " > A A A B / X i c d V D L S g M x F M 3 U V 6 2 v 8 b F z E y x C X T h k a l v b X d G N y w q 2 F t p a M m m m D U 1 m h i Q j 1 K H 4 K 2 5 c K O L W / 3 D n 3 5 g + B B U 9 c O F w z r 3 c e 4 8 X c a Y 0 Q h 9 W a m F x a X k l v Z p Z W 9 / Y 3 L K 3 d x o q j C W h d R L y U D Y 9 r C h n A a 1 r p j l t R p J i 4 X F 6 7 Q 3 P J / 7 1 L Z W K h c G V H k W 0 I 3 A / Y D 4 j W B u p a + / l M B T d p C 0 F p L 4 / P r p J j t 1 x 1 8 4 i p 1 J B B b c C k V N E K F 8 u G Y J O 8 u V i E b o O m i I L 5 q h 1 7 f d 2 L y S x o I E m H C v V c l G k O w m W m h F O x 5 l 2 r G i E y R D 3 a c v Q A A u q O s n 0 + j E 8 N E o P + q E 0 F W g 4 V b 9 P J F g o N R K e 6 R R Y D 9 R v b y L + 5 b V i 7 Z c 7 C Q u i W N O A z B b 5 M Y c 6 h J M o Y I 9 J S j Q f G Y K J Z O Z W S A Z Y Y q J N Y B kT w t e n 8 H / S y D t u y S l e F r L V s 3 k c a b A P D k A O u O A U V M E F q I E 6 I O A O P I A n 8 G z d W 4 / W i / U 6 a 0 1 Z 8 5 l d 8 A P W 2 y d 1 O p S f < / l a t e x i t > (am e↵ ) 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " u /O V M N M F J M H Y r R 1 V L L Y n o 1 d 4 4 P g = " > A A A B 8 X i c b V D L S g N B E O z 1 G e M r 6 t H L Y B A 8 h V 3 x d Q x 6 8 R j B P D B Z w u y k N x k y O 7 v O z A p h y V 9 4 8 a C I V / / G m 3 / j J N m D J h Y 0 F F X d d H c F i e D a u O 6 3 s 7 S 8 s r q 2 X t g o b m 5 t 7 + y W 9 v Y b O k 4 V w z q L R a x a A d U o u M S 6 4 U Z g K 1 F I o 0 B g M x j e T P z m E y r N Y 3 l v R g n 6 E e 1 L H n J G j Z U e a D f r q I j g 4 7 h b K r s V d w q y S L y c l C F H r V v 6 6 v R i l k Y o D R N U 6 7 b n J s b P q D K c C R w X O 6 n G h L I h 7 W P b U k k j 1 H 4 2 v X h M j q 3 S I 2 G s b E l D p u r v i Y x G W o + i w H Z G 1 A z 0 v D c R / / P a q Q m v / I z L J D U o 2 W x R m A p i Y j J5 n / S 4 Q m b E y B L K F L e 3 E j a g i j J j Q y r a E L z 5 l x d J 4 7 T i X V T O 7 8 7 K 1 e s 8 j g I c w h G c g A e X U I V b q E E d G E h 4 h l d 4 c 7 T z 4 r w 7 H 7 P W J S e f O Y A / c D 5 / A H 7 h k N E = < / l a t e x i t > a eq < l a t e x i t s h a 1 _ b a s e 6 4 = " S v 7 n G 4 M A f e Z 6 B V B r b o B I h 8 G n S S 8 = " > A A A B 6 H i c d Z D L S g M x F I Y z 9 V b r r e r S T b A I r o a M t q P u i i 5 0 2 Y K 9 Q D u U T H q m j c 1 c S D J C G f o E b l w o 4 t Z H c u f b m F 4 E F f 0 h c P j + c 8 g 5 v 5 8 I r j Q h H 1 Z u a X l l d S 2 / X t j Y 3 N r e K e 7 u N V W c S g Y N F o t Y t n 2 q Q P A I G p p r A e 1 E A g 1 9 A S 1 / d D X 1 W / c g F Y + j W z 1 O w A v p I O I B Z 1 Q b

Figure 3 .
Figure 3.The evolution of the individual modes as they evolve from sub-horizon during inflation to super-horizon then back into the horizon in the post-inflationary era.

Figure 4 .
Figure 4. Re-scaled power spectrum y f x 2 * |T L | 2 (see e.g.(4.16)) of the longitudinal modes evaluated for two different parameter choices (left and right panels) evaluated at 5 e-folds after H(a * ) = m, i.e. when ln(y f ) = ln(a f /a * ) = ∆N = 5 (top) and much later times y f ≫ e 5 (bottom).In the top panels, modes that reside on the right hand side of the dashed vertical line are still relativistic at the time a f of the evaluation satisfying k > a f m → k/k * = e ∆N = e 5 .At later times, the re-scaled power spectrum preserves its amplitude, however more modes shown in the UV tail becomes non-relativistic, adopting the k −1 behavior of the modes around the peak.

Figure 5 .
Figure 5. Same as Figure 4, but with different parameter choices.