UG-FT 333-21, CAFPE 203-21

# Dark photon dark matter from a rolling inflaton

Mar Bastero-Gil<EMAIL_ADDRESS>Departamento de Física Teórica y del Cosmos and
CAFPE, Universidad de Granada, Campus de Fuentenueva, E-18071 Granada, Spain
Jose Santiago<EMAIL_ADDRESS>Departamento de Física Teórica y del Cosmos and
CAFPE, Universidad de Granada, Campus de Fuentenueva, E-18071 Granada, Spain
Lorenzo Ubaldi<EMAIL_ADDRESS>SISSA and INFN Sezione di Trieste,Via
Bonomea 265, 34136 Trieste, Italy Institute for Fundamental Physics of the
Universe, Via Beirut 2, 34014 Trieste, Italy Roberto Vega-Morales
<EMAIL_ADDRESS>Departamento de Física Teórica y del Cosmos and CAFPE,
Universidad de Granada, Campus de Fuentenueva, E-18071 Granada, Spain

###### Abstract

We study in detail a recently proposed mechanism for producing non-thermal
dark photon dark matter at the end of inflation in the mass range $\mu\,{\rm
eV}\lesssim m\lesssim 10\,{\rm TeV}$. A tachyonic instability induced by a
rolling inflaton leads to the _coherent_ production of dark (abelian) gauge
bosons with a peak in the power spectrum corresponding to the Hubble scale at
the end of inflation. As the Universe expands after inflation the dark photons
redshift and, at some point in their cosmic evolution, they obtain a mass. We
focus in particular on the case where the dark photons are _relativistic_ at
the time their mass is generated and examine the associated cosmic evolution
to compute the relic abundance today. We also examine the late time power
spectrum demonstrating explicitly that it preserves the peak generated at the
end of inflation. We show that the peak corresponds to small physical scales
today, $\ell_{\rm today}\sim{\rm cm}-100\,{\rm km}$, with large density
fluctuations at $\ell_{\rm today}$ leading to a clumpy nature for the dark
photon dark matter. We also discuss potential phenomenology and future
directions, briefly commenting on the non-relativistic case.

###### Contents

  1. I Introduction
  2. II Vector dark matter production from end of inflation
    1. II.1 Tachyonic production in an expanding universe
    2. II.2 Energy density at the end of inflation
    3. II.3 Evolution of total energy density
    4. II.4 Vector dark matter parameter space
  3. III Energy Density Spectrum and Clumpy Dark Matter
    1. III.1 Energy density spectrum at the end of inflation
    2. III.2 Cosmological evolution of modes
    3. III.3 Late time energy density spectrum
    4. III.4 Isocurvature and density contrast power spectrum
    5. III.5 Clumping of vector dark matter
  4. IV Summary and Outlook
  5. IV.1 Derivation of equations of motion, energy density, and pressure
  6. IV.2 Analytic study of the energy density spectrum during inflation
  7. IV.3 Power spectrum of energy density fluctuations
  8. IV.4 Numerical procedure for solving equations of motion

##  I Introduction

The nature and production mechanism of dark matter still remains a mystery.
While most of the experimental effort in the past was aimed at detecting
weakly interacting massive particles (WIMPs), the lack of observation has
necessitated new theoretical ideas and proposed search strategies to cover as
many alternatives as possible Bertone:2018xtm . In this context, it is
important to explore new production mechanisms for dark matter and assess
their impact on search strategies. As the WIMP paradigm comes into increasing
tension, candidates for non-thermal dark matter have gained renewed interest
Bertone:2018xtm . Recently, one class of non-thermal dark matter candidates
receiving considerable attention is that of vector (or “dark photon” 111We use
‘dark photon’ and ‘dark vector’ interchangeably throughout to refer to a
neutral spin-one massive vector boson associated with a broken dark $U(1)_{D}$
gauge symmetry.) dark matter for which several production mechanisms have been
proposed. These include production mechanisms associated with inflation
Nelson:2011sf ; Arias:2012az ; Graham:2015rva ; Bastero-Gil:2018uel ;
Ema:2019yrd ; Nakayama:2020rka ; Nakai:2020cfw ; Ahmed:2020fhc ; Kolb:2020fwh
; Salehian:2020asa ; Firouzjahi:2020whk  as well as oscillating scalars after
inflation Agrawal:2018vin ; Co:2018lka ; Dror:2018pdh . The common feature
among them is a coupling between a dark abelian gauge boson and a separate
sector which induces a time dependence into the dispersion relation of the
dark vector field. This separate sector can come in the form of a scalar like
the inflaton Bertone:2018xtm ; Salehian:2020asa  or an axion Agrawal:2018vin ;
Co:2018lka  or simply gravity Graham:2015rva ; Ahmed:2020fhc ; Kolb:2020fwh .
Over large regions of parameter space this leads to exponential dark vector
production and can reproduce the observed dark matter relic abundance over
many orders of magnitude of dark matter mass.

Motivated by, but not limited to, scenarios of axion inflation, in Bastero-
Gil:2018uel  it was shown that a dark abelian gauge field coupled to an
inflaton via a $\phi F\tilde{F}$ coupling can be produced by a tachyonic
instability and generate the observed dark matter relic abundance in the mass
range $\mu\,{\rm eV}\lesssim m\lesssim 10\,{\rm TeV}$. More specifically, the
time dependance induced by the rolling inflaton leads to a tachyonic
enhancement and exponential production of _one transverse polarization_ of the
dark photon. As the Universe expands after inflation the dark gauge bosons
redshift and, at some point in their cosmic evolution, they obtain a mass and
become non-relativistc. As in Graham:2015rva  where the longitudinal mode is
produced by inflationary fluctuations, there is a peaked structure in the
power spectrum. However, in this case the peak is not due to redshifting, but
instead to the time dependence of the inflaton as it rolls down its potential.
As we examine in detail below, the dark photon production is exponentially
sensitive to the inflaton velocity. This leads to the maximum production just
at the end of inflation as the inflaton exits slow roll. This in turn gives
rise to a peak in the dark photon power spectrum at scales corresponding to
the Hubble scale at the end of inflation.

The goal of this work is to examine this mechanism in detail and to compute
the late time energy density spectrum. We focus in particular on the case
where the dark vector is _relativistic_ at the time its mass is generated and
examine the associated cosmic evolution. We first review the production
mechanism and track the total energy density to estimate the parameter space
which can reproduce the observed dark matter relic abundance. We then examine
the energy density spectrum at the end of inflation as well as its evolution
to late times once the dark vector has become non-relativistic. We obtain the
late time power spectrum demonstrating explicitly that the peak generated at
the end of inflation is preserved after redshifting. We then show that the
peak corresponds to small physical scales today, $\ell_{\rm today}\sim{\rm
cm}-100\,{\rm km}$, with potentially large density fluctuations at $\ell_{\rm
today}$ leading to a clumpy nature for the dark photon matter. We also discuss
potential phenomenology and future directions, briefly commenting on the non-
relativistic case.

##  II Vector dark matter production from end of inflation

Here we expand on the discussion presented in Bastero-Gil:2018uel  to show
that a dark vector field coupled to a slow rolling inflaton via a $\phi
F\tilde{F}$ coupling can be produced by a tachyonic instability and generate
the observed dark matter relic abundance. To do this we first derive the
equations of motion and obtain (approximate) analytic solutions for the
tachyonic modes. We then compute the total energy density at the end of
inflation and track its cosmic evolution.

Since we work in the weak field regime where backreaction effects can be
neglected, the evolution of the dark vector modes is linear. Thus, while the
vector field begins as quantum fluctuations during inflation that turn
classical only later on, the mode functions of the creation and destruction
operators can be obtained using the classical equations of motion. The
evolution of the energy density and power spectrum can then be directly
extracted from the mode functions obeying classical evolution.

###  II.1 Tachyonic production in an expanding universe

Our starting point is the action for an inflaton with potential $V(\phi)$
coupled to a spin-1 vector boson which is neutral under the Standard Model
gauge group,

$\displaystyle S$ $\displaystyle=$ $\displaystyle-\int
d^{4}x\,\sqrt{-g}\Big{[}\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi+V(\phi)+\leavevmode\nobreak\
\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\frac{1}{2}m^{2}A_{\mu}A^{\mu}+\frac{\alpha}{4f}\phi
F_{\mu\nu}\tilde{F}^{\mu\nu}\Big{]}\,,$ (1)

where $\phi$ is the inflaton field that drives inflation, $A_{\mu}$ is the
dark vector field, $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is
the field strength, and
$\tilde{F}^{\mu\nu}=\epsilon^{\mu\nu\alpha\beta}F_{\alpha\beta}/2$ with
$\epsilon^{\mu\nu\alpha\beta}$ the completely antisymmetric tensor. We use the
Friedmann-Robertson-Walker metric with $ds^{2}=-dt^{2}+a^{2}(t)d\vec{x}^{2}$
and the convention $\epsilon^{0123}=1/\sqrt{-g}$. The vector boson mass $m$
can be zero or non-zero during inflation and can be of a Stueckelberg or
‘Higgs-ed’ type with an associated symmetry breaking phase transition. As long
as it is smaller than the Hubble scale at the end of inflation the vector mass
has negligible effects on the tachyonic production mechanism and only becomes
relevant when considering cosmological evolution and the final relic
abundance. There could in principle be a kinetic mixing between the visible
and dark photons, but it does not spoil the production mechanism as long as it
is small enough to prevent thermalization of the dark vector with the visible
sector and decay to Standard Model particles for masses above the electron
threshold as well as satisfy other experimental constraints Gherghetta:2019coi
. We also do not specify the inflaton potential $V(\phi)$ as its precise form
is not crucial for the production mechanism. However, the form of the
inflationary potential can affect the shape of the dark photon energy density
spectrum as we examine in more detail in Sec. III.

The only crucial ingredients needed for the production mechanism are: $(i)$
the scale of the inflaton potential which sets the Hubble scale during
inflation,

$\displaystyle H=\frac{\sqrt{V(\phi)}}{\sqrt{3}M_{\rm Pl}}\,,$ (2)

and $(ii)$ the coupling of the inflaton to $F\tilde{F}$ responsible for
exponential production of only one polarization of the transverse modes. Such
a coupling is generically present in models of natural inflation, where $\phi$
is a pseudoscalar (odd under parity) axion-like field subject to a shift
symmetry. For this reason, this class of models provides a well motivated
theoretical framework for the mechanism presented here. However, since the
dynamics of the mechanism do not depend on whether $\phi$ is an axion or not,
$\phi$ can be a generic scalar (or function of $\phi$ Barnaby:2011qe ) as long
as it is rolling towards its minimum and couples to $F\tilde{F}$ 222These
ingredients are also present in some models aimed at the relaxation of the
electroweak scale Tangarife:2017rgl  or of the cosmological constant
Graham:2019bfu .. Note that we do not need to impose that the Lagrangian
respects parity so $\phi$ could also be parity even and, in particular, is not
necessarily a pseudo Goldstone boson.

We quantize the vector field by expanding in the helicity basis in terms of
creation and annihilation operators and their mode functions as follows,

$\displaystyle\hat{\vec{A}}(\vec{x},t)$ $\displaystyle=$
$\displaystyle\sum_{\lambda=\pm,L}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}}\
\vec{\epsilon}_{\lambda}(\vec{k})\times[A_{\lambda}(k,t)a_{\lambda}(\vec{k})+A_{\lambda}(k,t)^{\ast}a_{\lambda}^{\dagger}(-\vec{k})],$
(3)

where we include transverse and longitudinal polarisation in the sum and the
creation and annihilation operators satisfy the commutation relation,

$\displaystyle\Big{[}a_{\lambda}(\vec{k}),\,a_{\lambda}^{\dagger}(\vec{k}^{\prime})\Big{]}=(2\pi)^{3}\delta_{\lambda\lambda^{\prime}}\,\delta^{3}(\vec{k}-\vec{k}^{\prime}).$
(4)

The mode functions obey the equations of motion derived in the Appendix
starting from the action in Eq.(1) and which in Fourier space read
Anber:2009ua ; Barnaby:2010vf ,

$\displaystyle\ddot{\phi}$ $\displaystyle+$ $\displaystyle
3H\dot{\phi}+V^{\prime}=\frac{\alpha}{4f}F\tilde{F},$ (5)
$\displaystyle\ddot{A}_{\pm}$ $\displaystyle+$ $\displaystyle
H\dot{A}_{\pm}+\left(\frac{k^{2}}{a^{2}}\pm\frac{k}{a}\frac{\alpha\dot{\phi}}{f}+m^{2}\right)A_{\pm}=0\,,$
(6) $\displaystyle\ddot{A}_{L}$ $\displaystyle+$
$\displaystyle\frac{3k^{2}+a^{2}m^{2}}{k^{2}+a^{2}m^{2}}H\dot{A}_{L}+\left(\frac{k^{2}}{a^{2}}+m^{2}\right)A_{L}=0\,,$
(7)

where we have also included the inflaton equation of motion. The overdots
denote derivatives with respect to physical time $t$ and $k\equiv|\vec{k}|$ is
the magnitude of the comoving momentum. We consider only the spatially
homogeneous zero momentum mode ($k=0$) of the inflaton. We have separated the
three degrees of freedom of the vector into transverse and longitudinal
components, $\vec{A}_{T}$ and $A_{L}$ respectively, where
$\vec{k}\cdot\vec{A}=kA_{L}$ and $\vec{k}\cdot\vec{A}_{T}=0$, and we have
written the transverse component in terms of the two helicities,
$\vec{A}_{T}=\vec{\epsilon}_{+}A_{+}+\vec{\epsilon}_{-}A_{-}$. We see
explicitly the $\phi F\tilde{F}$ coupling only enters into the equations of
motion for the transverse modes. The equation of motion for $A_{L}$ then
corresponds to the one derived in Graham:2015rva  and thus, as demonstrated in
Graham:2015rva , if $m\neq 0$ during inflation the longitudinal mode is
produced via inflationary fluctuations and will also contribute to the dark
vector energy density.

In what follows we concentrate on the equation of motion of the transverse
vector modes in Eq.(6). It is convenient to introduce the dimensionless
“instability parameter’,

$\displaystyle\xi\equiv\frac{\alpha\dot{\phi}}{2Hf}=\sqrt{\frac{\epsilon}{2}}\frac{\alpha}{f}M_{\rm
Pl}\,,$ (8)

where $\epsilon\equiv-\dot{H}/H^{2}$ and for single field inflation we have,

$\displaystyle|\dot{\phi}|\approx V^{\prime}/3H,\leavevmode\nobreak\
\leavevmode\nobreak\ \epsilon=\frac{\dot{\phi}^{2}}{2H^{2}M_{Pl}^{2}}.$ (9)

We then rewrite the equation of motion in terms of conformal time $\tau$
defined as $ad\tau=dt$,

$\displaystyle\Big{[}\frac{\partial^{2}}{\partial\tau^{2}}+k^{2}\pm
2\,k\,\frac{\xi}{\tau}+\frac{\bar{m}^{2}}{\tau^{2}}\Big{]}A_{\pm}(k,\tau)\equiv\Big{[}\frac{\partial^{2}}{\partial\tau^{2}}+\omega^{2}(k,\tau)\Big{]}A_{\pm}(k,\tau)=0\,,$
(10)

where we have defined the dimensionless ratio,

$\displaystyle\bar{m}\equiv\frac{m}{H}.$ (11)

Without loss of generality we use the convention $\dot{\phi}>0$ which gives
$\xi>0$, implying that _only_ the mode $A_{+}$ experiences a tachyonic
instability when,

$\displaystyle\omega^{2}(k,\tau)=k^{2}+2k\frac{\xi}{\tau}+\frac{\bar{m}^{2}}{\tau^{2}}=k^{2}-2k\xi
aH+\bar{m}^{2}(aH)^{2}<0\,,$ (12)

where we have used the fact that during inflation ($\tau<0$) we have
$\tau\simeq-\frac{1}{aH}$. On the other hand, the opposite polarization
$A_{-}$ does not have tachyonic modes and is therefore neglected. The
condition in Eq.(12) then leads to the tachyonic conditions on the vector mass
(which could be zero during inflation) and physical momentum ($q$),

$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ q\equiv\frac{k}{a}<\xi H+\xi
H\sqrt{1-\bar{m}^{2}/\xi^{2}}\leavevmode\nobreak\ \leavevmode\nobreak\
\rm{\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
(tachyonic\leavevmode\nobreak\ condition)},$ (13)

which requires $\bar{m}<\xi$ to avoid the mass term quenching the tachyonic
production Meerburg:2012id . When $\bar{m}\ll\xi$, we also see the tachyonic
condition on the physical momentum becomes $q<2\xi H$. As we will see below,
we are interested in the (weak coupling) regime where $\xi$ is
$\mathcal{O}(1)$ implying that the modes become tachyonic as their physical
wavelength is stretched to be of order the horizon, $\lambda\equiv q^{-1}\sim
H^{-1}$. Since $-\omega^{2}$ is maximal at $q\simeq\xi H$, this implies the
vector field power spectrum should also have a peak at scales around the co-
moving horizon. Within the horizon these modes add up coherently and have very
large occupation number due to the exponential enhancement from the tachyonic
instability. Thus, these dark vectors are well described by a classical dark
‘electromagnetic’ field and since only one helicity is enhanced exponentially,
we are left at the end of inflation with a maximally _helical dark
electromagnetic field_. This is in analogy with magnetogenesis scenarios
Adshead:2016iae  constructed to solve the puzzle of the origin of primordial
magnetic fields.

Treating $\xi$ as a constant, a good approximation early on during inflation
as the inflaton slow-rolls, one can solve Eq.(10) analytically Meerburg:2012id
in terms of the Whittaker functions. The overall normalization is determined
by the requirement that the gauge field is initially (in the sub-horizon limit
$-k\tau\to\infty$) in the Bunch-Davies (BD) vacuum 333This choice of
normalization makes the classical mode solutions consistent with the quantum
ones.,

$\displaystyle\lim_{-k\tau\to\infty}A_{\pm}(k,\tau)=\frac{e^{-ik\tau}}{\sqrt{2k}}\equiv
A_{\rm{BD}}\,.$ (14)

Neglecting the dark photon mass (with $\xi$ as constant) the full analytic
solution to Eq.(10) can be found in terms of Coulomb functions Meerburg:2012id
; Adshead:2015pva ; Domcke:2018eki  which, in the tachyonic regime we are
interested in $-k\tau<2\xi$ ($k<2\xi aH$), are very well approximated by
Barnaby:2011vw ,

$\displaystyle A_{+}(k,\tau)$ $\displaystyle\simeq$ $\displaystyle
e^{\pi\xi}\sqrt{\frac{-2\tau}{\pi}}K_{1}\left[2\sqrt{-2\xi k\tau}\right]\,,$
(15)

where $K_{1}$ is a modified Bessel function of the second kind. In this form
we see the exponential dependence on $\xi$ (or $\dot{\phi}$) explicitly. In
the super-horizon (SH) limit ($-k\tau\to 0$) Eq.(15) gives,

$\displaystyle\lim_{-k\tau\to 0}A_{+}(k,\tau)=\frac{e^{\pi\xi}}{2\sqrt{2\pi
k\xi}}\equiv A_{\rm{SH}}\,.$ (16)

Useful analytic solutions can be obtained via the WKB approximation
Anber:2009ua ; Barnaby:2011vw ; Tangarife:2017rgl ,

$\displaystyle A_{+}(k,\tau)_{\rm WKB}$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\sqrt{2k}}\left(\frac{-k\tau}{2\xi}\right)^{1/4}e^{\pi\xi-2\sqrt{-2\xi
k\tau}}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ (\frac{1}{8\xi}<-k\tau<2\xi)\,,$ (17)

where the regime of validity is dictated by the adiabatic condition
$|\Omega^{\prime}/\Omega^{2}|\ll 1$ with $\Omega^{\prime}\equiv
d\Omega/d\tau$. For modes with $k$ in the range $aH/(8\xi)<k<2\xi aH$,
$A_{+}(k,\tau)_{\rm WKB}$ approximates very well the solution obtained in
Eq.(15) and gives us intuition into the behavior of the modes around horizon
crossing as they become exponentially enhanced. This also allows us to use
analytic solutions to obtain an estimate of the total vector dark matter relic
abundance and the viable regions of parameter space for the mechanism.
Eventually, we will need to compute the power spectrum at the end of inflation
to use as input when tracking the cosmological evolution of the energy density
spectrum. However, as we emphasize in Sec. III, a more precise calculation of
the dark matter energy density and shape of the power spectrum requires
accounting for the time dependence of $\xi$ which necessitates solving the
system of equations in Eq.(5) and Eq.(6) numerically.

###  II.2 Energy density at the end of inflation

To eventually obtain the final relic abundance for the vector dark matter we
need to track the evolution of the energy density starting from the time of
its production at the end of inflation. Thus we need to first compute the
total vector energy density at the end of inflation. As we review in the
Appendix, starting from the action in Eq.(1) we can obtain the total energy
density for the transverse component of the dark vector field ($\rho_{D}$) in
terms of the tachyonic mode amplitude and its (conformal) time derivative,

$\displaystyle\rho_{D}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi^{2}a^{4}}\int_{0}^{\infty}dk\,k^{2}\Big{(}|\partial_{\tau}A_{+}(k,\tau)|^{2}+\left(k^{2}+a^{2}m^{2}\right)|A_{+}(k,\tau)|^{2}\Big{)}$
(18) $\displaystyle=$ $\displaystyle\frac{1}{2a^{4}}\int d\ln
k\Big{(}\mathcal{P}_{\partial_{\tau}A_{+}}(k,\tau)+\left(k^{2}+a^{2}m^{2}\right)\mathcal{P}_{A_{+}}(k,\tau)\Big{)}$
$\displaystyle=$
$\displaystyle\langle(\partial_{\tau}A_{+})^{2}\rangle+\langle
A_{+}^{2}\rangle\,=\frac{1}{2}\langle\vec{E}^{2}+\vec{B}^{2}\rangle.$

where $\rho_{D}\equiv\langle\rho_{D}\rangle$ represents the spatial average as
defined in the Appendix (see Eq.(IV.1)) and we identify the ‘magnetic’ ($B$)
and ‘electric’ ($E$) components respectively. We have used $dt=ad\tau$ as well
as defined the field and (time) derivative power spectra,

$\displaystyle\mathcal{P}_{X}(k,\tau)$ $\displaystyle=$
$\displaystyle\frac{k^{3}}{2\pi^{2}}|X(k,\tau)|^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ X=A_{+}\ {\rm or}\ \
\partial_{\tau}A_{+}\,,$ (19)

which allows us to define the electric and magnetic energy density spectra
respectively,

$\displaystyle\frac{d\rho_{E}}{d\,{\rm ln}\,k}$ $\displaystyle=$
$\displaystyle\frac{1}{2a^{4}}\mathcal{P}_{\partial_{\tau}A_{+}}(k,\tau),\leavevmode\nobreak\
\leavevmode\nobreak\ \frac{d\rho_{B}}{d\,{\rm
ln}\,k}=\frac{1}{2a^{4}}\left(k^{2}+a^{2}m^{2}\right)\mathcal{P}_{A_{+}}(k,\tau)\,.$
(20)

Note as part of the magnetic component we have included the mass term which of
course is not present in the case of the visible electromagnetic field. Since
only one transverse mode is exponentially enhanced (which we take as $A_{+}$)
by the tachyonic instability, we can safely neglect the contribution from
$A_{-}$ to Eq.(18). For the case where the dark vector already has a mass
during inflation, we also compute the energy density contained in the
longitudinal mode and obtain the same result as in Graham:2015rva .

Taking $\bar{m}\ll 1$ during inflation and neglecting the time dependence of
$\xi$, we can use the WKB solution for $A_{+}$ in Eq.(17) to compute
analytically the dark electric field contribution given by
$\mathcal{P}_{\partial_{\tau}A_{+}}(k,\tau)$ as well as the dark magnetic
field contribution $k^{2}\mathcal{P}_{A_{+}}(k,\tau)$, where the former always
dominates over the latter during inflation. Integrating over momenta we can
then estimate the energy density contained in the dark vector during
inflation,

$\displaystyle\rho_{D}\approx 10^{-4}\frac{e^{2\pi\xi}}{\xi^{3}}\,H^{4}.$ (21)

In reality of course, both the Hubble parameter and $\xi$ depend on time. The
parameter $\xi$, which controls the dark photon production and grows with
$\dot{\phi}$ (or the slow roll parameter $\epsilon$), is largest at the end of
inflation. Thus, the largest contribution to the dark electromagnetic energy
density comes from the end of inflation, as confirmed in our numerical
analysis discussed below. From Eq.(21) we can estimate the energy density at
the end of inflation as,

$\displaystyle\rho_{D}(a_{\rm{end}})\equiv\rho^{\rm{end}}_{D}\approx
10^{-4}\frac{e^{2\pi\xi_{\rm end}}}{\xi_{\rm end}^{3}}H_{\rm
end}^{4}\,=10^{-4}\frac{e^{2\pi\xi_{\rm end}}}{\xi_{\rm
end}^{3}}\,\epsilon_{H}^{4}H^{4},$ (22)

with $\xi_{\rm end}$ the value of $\xi$ at the end of inflation and given by
(using Eq.(8)),

$\displaystyle\xi_{\rm end}=\frac{\alpha}{\sqrt{2}}\frac{M_{\rm Pl}}{f},$ (23)

where we have assumed $\epsilon=1$ at the end of inflation444Note in hybrid
inflation models Copeland:1994vg ; Dvali:1994ms , one can have values
different from $\epsilon\approx 1$ at the end of inflation.. In the last
equality in Eq.(22) we account for a decreasing Hubble parameter during
inflation and parametrize it as,

$\displaystyle H_{\rm end}=\epsilon_{H}H\,,$ (24)

with $\epsilon_{H}$ a dimensionless parameter that can be calculated in a
particular model of inflation. Typically the slow-roll parameter
$\epsilon=-\dot{H}/H^{2}$ is $\mathcal{O}(10^{-2}-10^{-3})$ during inflation
and $\approx 1$ at the end of inflation. This translates into values of
$\epsilon_{H}$ that are model dependent, but for most models of inflation we
expect $\epsilon_{H}$ to be in the range $10^{-3}<\epsilon_{H}<10^{-1}$. In
axion inflation models involving tachyonic production of vector fields, which
have roughly 60 e-folds of inflation, $\xi$ during the first few e-folds is
constrained by CMB measurements Barnaby:2010vf ; Barnaby:2011vw  to be less
than $\xi_{\rm CMB}\lesssim 2.5$, but $\xi_{\rm end}$ is allowed to be
significantly larger. However, if $\xi$ is too large this will induce back
reaction effects Barnaby:2011qe ; Peloso:2016gqs  which must be accounted for
and which could potentially destroy the production mechanism. Thus, in order
to neglect these back reaction effects we limit ourselves to $\xi_{\rm
end}<\mathcal{O}(10)$ in the following. On the other hand, to obtain
sufficient dark vector production requires $\xi_{\rm end}\gtrsim 1$ leading us
to consider the range $1\lesssim\xi_{\rm end}\lesssim 10$ when examining the
viable vector dark matter parameter space below.

During inflation, when $\phi$ dominates the energy density, we have for the
inflaton,

$\displaystyle\rho_{I}=V(\phi)=3H^{2}M_{\rm Pl}^{2}\,,$ (25)

while at the end of inflation, with the inflaton still dominating, we have
$H=H_{\rm end}$ and,

$\displaystyle\rho^{\rm{end}}_{I}=3H_{\rm end}^{2}M_{\rm
Pl}^{2}=3\,\epsilon_{H}^{2}H^{2}M_{\rm Pl}^{2}\,.$ (26)

A fraction of $\rho^{\rm{end}}_{I}$ is transferred into the dark vector while
another fraction goes into visible radiation to reheat the Universe whose
energy density we can write as,

$\displaystyle\rho_{R}(T_{\rm RH})=\frac{\pi^{2}}{30}g_{*}(T_{\rm RH})T_{\rm
RH}^{4}\,\equiv\epsilon_{R}^{4}\,\rho_{I}=3H_{\rm RH}^{2}M_{\rm Pl}^{2}.$ (27)

Combining with Eq.(25), this allows us to define the reheating temperature,

$\displaystyle T_{\rm RH}=\epsilon_{R}\,\left(\frac{90}{\pi^{2}g_{*}(T_{\rm
RH})}\right)^{1/4}\sqrt{HM_{\rm Pl}}\,,$ (28)

as well as the Hubble scale at reheating in terms of the Hubble scale during
inflation,

$\displaystyle H_{\rm RH}=\epsilon_{R}^{2}H\,.$ (29)

The dimensionless parameter $\epsilon_{R}<1$ parametrizes the fraction of the
inflaton energy that goes into visible radiation. We take $g_{*}(T_{\rm RH})$
to denote the number of relativistic degrees of freedom which we fix to
$g_{*}(T_{\rm RH})\sim 100$ and restrict ourselves to reheating temperatures
above the electroweak scale. We work in the approximation of instantaneous
reheating assuming it takes place as soon as the inflaton exits slow-roll.
Thus, we assume the Universe has a temperature $T_{\rm RH}$ at $a=a_{\rm end}$
(see Fig. 3) implying $\rho^{\rm{end}}_{D}=\rho_{D}(T_{\rm RH})$. Note
however, we do not assume $H_{\rm RH}=H_{\rm end}$ which will be important to
account for when we discuss constraints on the vector dark matter parameter
space below.

###  II.3 Evolution of total energy density

In this section we track the redshift of the dark vector energy density after
production at the end of inflation and estimate the present day relic
abundance. Here we only consider the case where the dark vector is
relativistic at the time its mass is generated. We also assume the dark vector
mass is already present during inflation and remains non-zero throughout its
cosmic evolution. This scenario can be applied to either a Stueckelberg or
dark Higgs mechanism for generating the dark vector mass. We study the non-
relativistic case and other possibilities for cosmic evolution when the dark
vector mass is generated via a dark Higgs mechanism in Higgsfollowup  where we
find a different result from Eq.(37) for the relic density. For purposes of
tracking the evolution and estimating the viable dark matter parameter space,
it is enough to track the redshift of the energy density considering only
modes around the peak of the power spectrum for which we assume they all
redshift together.

As discussed in Sec. II.1, for a given scale factor, the power spectrum will
be peaked around scales the size of the co-moving horizon. Thus, at the end of
inflation the modes which give the largest contribution to
$\rho^{\rm{end}}_{D}$ have physical momentum,

$\displaystyle q(T_{\rm RH})\equiv\frac{k}{a_{\rm end}}\simeq H_{\rm end}\gg
m\,,$ (30)

where again we have assumed instantaneous reheating allowing us to define
$T_{\rm RH}$ and track the redshift using temperature instead of scale factor.
At reheating the dark vectors are relativistic with the physical momentum then
redshifting as,

$\displaystyle q(T)=q(T_{\rm RH})\frac{T}{T_{\rm RH}}\,.$ (31)

The dark photons become non relativistic at a temperature $T=\bar{T}$ defined
by the condition,

$\displaystyle q(\bar{T})=m\,,$ (32)

which combining with Eq.(24), Eq.(28) and Eq.(30) allows us to solve for
$\bar{T}$ as,

$\displaystyle\bar{T}=\frac{m}{H_{\rm end}}\,T_{\rm
RH}=m\left(\frac{90}{\pi^{2}g_{*}(T_{\rm
RH})}\right)^{1/4}\frac{\epsilon_{R}}{\epsilon_{H}}\left(\frac{M_{\rm
Pl}}{H}\right)^{1/2}\,.$ (33)

Above $\bar{T}$ (or before scale factor $a=\bar{a}$) the vector energy density
redshifts like radiation,

$\displaystyle\rho_{D}(T)=\rho_{D}(T_{\rm RH})\left(\frac{T}{T_{\rm
RH}}\right)^{4}\,,$ (34)

while below $\bar{T}$ (after $a=\bar{a}$) it redshifts like matter giving,

$\displaystyle\rho_{D}(T)=\rho_{D}(T_{0})\left(\frac{T}{T_{0}}\right)^{3}\,,$
(35)

where $T_{0}\approx 10^{-13}$ GeV is today’s CMB temperature. The cosmic
evolution described above and in the previous section can be summarized with
Fig. 1 which shows the evolution of the various energy densities defined in
Eq.(22), Eq.(25), and Eq.(27).

Figure 1: Schematic representation of the cosmic evolution of the vector dark
matter energy density ($\rho_{D}$) as well as the energy density in visible
radiation ($\rho_{R}$) and the inflaton ($\rho_{I}$).

Equating Eq.(34) and Eq.(35) at $T=\bar{T}$ and combining with Eq.(22),
Eq.(28), and Eq.(33), we obtain the energy density today in terms of the
energy density at the end of inflation,

$\displaystyle\rho_{D}(T_{0})$ $\displaystyle=$ $\displaystyle
m\,T_{0}^{3}\left(\frac{90}{\pi^{2}g_{*}(T_{\rm
RH})}\right)^{-3/4}\left(\frac{\epsilon_{H}}{\epsilon_{R}}\right)^{3}\left(\frac{H}{M_{\rm
Pl}}\right)^{3/2}\left(\frac{\rho^{\rm{end}}_{D}}{H_{\rm end}^{4}}\right)\,.$
(36)

This expression assumes the modes around the peak all redshift together which
is strictly speaking not true as modes of different momenta become non-
relativistic at different times, but for present purposes the approximation in
Eq.(36) is sufficient. The dark vector energy density at the end of inflation
$\rho^{\rm{end}}_{D}$ (normalized to $H_{\rm end}^{4}$) can be obtained
numerically, but as found in Eq.(22) a useful analytic result can be obtained
with the WKB approximation.

Using the WKB approximation and taking the observed energy density of cold
dark matter today Ade:2015xua  to be $\rho_{\rm CDM}=9.6\times 10^{-48}\ {\rm
GeV}^{4}$, the contribution from the transverse dark vector mode produced via
tachyonic instability can then be written as,

$\displaystyle\frac{\rho_{D}(T_{0})}{\rho_{\rm
CDM}}=\frac{\Omega_{T}}{\Omega_{\rm CDM}}\simeq 2\times
10^{-4}\cdot\left(\frac{m}{\rm MeV}\right)\left(\frac{H}{10^{14}\ {\rm
GeV}}\right)^{3/2}\left(\frac{e^{2\pi\xi_{\rm end}}}{\xi_{\rm
end}^{3}}\right)\,\left(\frac{\epsilon_{H}}{\epsilon_{R}}\right)^{3}.$ (37)

We see in this case the final relic abundance depends on five parameters. The
first two are the dark vector mass $m$ and the Hubble scale during inflation
$H$, as in the case of the longitudinal mode Graham:2015rva . The third
parameter $\xi_{\rm end}$ parametrizes the strength of the inflaton - dark
vector coupling (see Eq.(8)) and depends weakly on the precise shape of the
inflaton potential. The final two parameters parametrize our ignorance of the
inflaton sector, $\epsilon_{H}$ which parametrizes how much $H$ has decreased
by the end of inflation (see Eq.(24)) and $\epsilon_{R}$ which parametrizes
the fraction of energy transferred from the inflaton to the reheating sector
(see Eq.(27)). Both can in principle be calculated in a particular model of
inflation.

If the dark vector has a mass already during inflation, whether of
Stueckelberg or Higgs-ed type (with the $U(1)_{D}$ _not_ restored during
reheating), which persists throughout the entirety of its cosmic evolution
until today, there is also a contribution to the relic density from the
longitudinal mode produced via inflationary quantum fluctuations
Graham:2015rva ,

$\displaystyle\frac{\Omega_{L}}{\Omega_{\rm CDM}}=\left(\frac{m}{6\times
10^{-15}\ {\rm GeV}}\right)^{1/2}\left(\frac{H}{10^{14}\ {\rm
GeV}}\right)^{2}\,.$ (38)

We see this only depends on the vector mass and Hubble scale during inflation.
In this case both the transverse and longitudinal modes could in principle
contribute appreciably to the dark matter relic abundance. Thus, for these
cases of dark vector mass generation, we will also include the longitudinal
mode when exploring the viable parameter space below.

###  II.4 Vector dark matter parameter space

Having followed the cosmic evolution of the the dark vector energy density and
found an estimate for the relic abundance today, we can go on to estimate the
regions of parameter space for a viable dark matter candidate. To find these
regions, we must ensure that not only the observed cold dark matter relic
abundance is reproduced, but also that a number of constraints on the
parameters in Eq.(37) are satisfied for consistency of the mechanism. As
discussed in the previous section, we will focus only on scenarios where the
dark vector is relativistic at the time its mass is generated and assume in
the Higgs-ed case that once broken, the dark $U(1)_{D}$ is not restored by
reheating or any other phase transition Higgsfollowup .

_Constraints on model parameters_

  1. 1.

For the constraints on the dark vector mass the upper bound results from
requirng efficient tachyonic production at the end of inflation and constrains
only the case of $m\neq 0$ during inflation. At the same time $m$ is bounded
from below by the condition that the dark vector becomes non-relativistic,
thus behaving like cold dark matter, before matter-radiation equality. We
write this condition as $\bar{T}>T_{\rm CMB}$ (see Eq.(33)), with $T_{\rm
CMB}\simeq 10^{-9}$ GeV. Thus $m$ is constrained to lie in the window,

$\displaystyle\left(\frac{\pi g_{*}(T_{\rm
RH})}{90}\right)^{1/4}\frac{\epsilon_{H}}{\epsilon_{R}}\sqrt{\frac{H}{M_{\rm
Pl}}}\,T_{\rm CMB}<m$ $\displaystyle<$ $\displaystyle\epsilon_{H}H\,.$ (39)

This window allows for dark vector masses spanning many orders of magnitude.

  2. 2.

The energy density in the inflaton at the end of inflation must be larger than
in the radiation which reheats the Universe,
$\rho^{\rm{end}}_{I}>\rho_{R}(T_{\rm RH})$ implying $H_{\rm end}>H_{\rm RH}$
and,

$\displaystyle\frac{\epsilon_{H}}{\epsilon_{R}^{2}}>1\,,$ (40)

to ensure the Universe does not reheat to energies greater than
$\rho^{\rm{end}}_{I}$.

  3. 3.

At $a=a_{\rm RH}=a_{\rm end}$, the energy density of radiation in the
reheating sector must be greater than the energy density of the dark vector,
$\rho_{R}(T_{\rm RH})>\rho^{\rm{end}}_{D}$ which leads to,

$\displaystyle\frac{\epsilon_{H}}{\epsilon_{R}}\ll 10\ \xi_{\rm
end}^{3/4}e^{-\pi\xi_{\rm end}/2}\left(\frac{M_{\rm Pl}}{H}\right)^{1/2}\,.$
(41)

Otherwise, the Universe would become (dark) matter dominated at a temperature
above $T_{\rm CMB}$, thus violating matter-radiation equality at $T_{\rm
CMB}$.

  4. 4.

If the inflaton - dark vector coupling is too large, the dark vector can
thermalize with the inflaton which must also couple to Standard Model
particles to reheat the Universe. This would lead to thermalization of the
dark vector with the visible sector and spoil our dark matter production
mechanism. Ensuring the inflaton and dark vector do not thermalize
Ferreira:2017lnd ; Ferreira:2017wlx  puts an upper bound on $\xi$,

$\displaystyle\xi<0.44\ln\frac{f}{\alpha H}+3.4\,.$ (42)

Using Eq.(8) and taking the slow-roll parameter to be $\epsilon=1$ at $a_{\rm
end}$ then gives,

$\displaystyle\xi_{\rm end}<0.44\ln\left(\frac{1}{\sqrt{2}\xi_{\rm
end}}\frac{M_{\rm Pl}}{H}\right)+3.4\,.$ (43)

  5. 5.

We also assume that back-reaction effects on the inflaton dynamics are
negligible which leads to two conditions. The first is $3H\dot{\phi}\simeq
V^{\prime}\gg\alpha/f\langle\vec{E}\cdot\vec{B}\rangle$ meaning that the
$\langle F\tilde{F}\rangle=\langle\vec{E}\cdot\vec{B}\rangle$ term is
negligible in the inflaton equation of motion in Eq.(5). The second condition
is that $3H^{2}M_{\rm Pl}^{2}\gg\langle\vec{E}^{2}\rangle/2$, ensuring that
the inflaton dominates the energy density during inflation rather than the
dark vector. Both conditions are satisfied as long as $\xi$ is not too large
Barnaby:2011qe ; Peloso:2016gqs . Requiring that they hold all the way to the
end of inflation results in the constraint on $\xi_{\rm end}$,

$\displaystyle\frac{\epsilon_{H}H}{M_{\rm Pl}}\ll 10^{2}\,\xi_{\rm
end}^{3/2}\,e^{-\pi\xi_{\rm end}}\,.$ (44)

  6. 6.

When the inflaton exits the slow-roll regime, it starts oscillating about the
minimum of its potential and reheats the Universe. If the coupling $\alpha/f$
is moderately large, roughly $\alpha/f>35M_{\rm Pl}^{-1}$, the production of
dark vectors during these oscillations can be important. This phenomenon is
referred to as gauge-preheating and has been studied in Adshead:2015pva .
However, in order to satisfy the constraints listed above, here we consider
the range $1\lesssim\xi_{\rm{end}}\lesssim 10$ which results in the window for
$\alpha/f$,

$\displaystyle\frac{\sqrt{2}}{M_{\rm{Pl}}}\lesssim\frac{\alpha}{f}\lesssim\frac{10\sqrt{2}}{M_{\rm{Pl}}}\,,$
(45)

where we have used Eq.(23). In this range of $\alpha/f$, preheating into dark
vectors is not efficient Adshead:2015pva  implying that during the inflaton
oscillations only a negligible fraction of its energy density is transferred
to the dark vector. We then assume that reheating proceeds via the
perturbative decay of the inflaton into visible radiation.

Finally, we also implicitly assume there are no other light scalar or fermion
fields in the dark sector which couple to the dark vector. If such light
fields were present, they would be produced by the strong dark electromagnetic
field via the Schwinger effect Tangarife:2017vnd ; Tangarife:2017rgl  and
would in principle contribute to the dark matter abundance today. Here we
assume they are heavy enough to avoid this interesting possibility which we
leave to forthcoming work DSfollowup .

_Final viable parameter space_

In Fig. 2 we show the relic abundance given in Eq.(37) as a function of $m$
and $H$, imposing the constraints listed above, for different values of the
parameters $\xi_{\rm end}$, $\epsilon_{R}$, $\epsilon_{H}$. In practice,
requiring the dark vector is not over abundant along with the first two
constraints automatically ensures the remaining constraints are satisfied, so
only these are shown. Along the contours labeled “Transverse” for different
values of $\xi_{\rm end}$ which separate the purple shaded regions we obtain
the observed relic abundance with the transverse mode making up the entirety
of the dark matter. In the colored regions to the right of these lines the
dark matter is overabundant. We see that the transverse mode of the dark
vector can make a viable dark matter candidate over a wide range of parameter
space: $\mu{\rm eV}\lesssim m\lesssim{\rm TeV}$, $100\ {\rm GeV}\lesssim
H\lesssim 10^{14}\ {\rm GeV}$ for $\xi_{\rm end}\sim\mathcal{O}(1-10)$ (see
Eq.(23)).

Figure 2: Parameter space in the dark photon mass versus Hubble scale $m-H$
plane for values of the parameters $\xi_{\rm end}$, $\epsilon_{R}$, and
$\epsilon_{H}$ as indicated above each plot and described in the text. Along
the lines labeled “Transverse” for different values of $\xi_{\rm end}$, we
obtain the observed relic abundance for the transverse mode while in the
colored regions to the right of these lines the dark matter is overabundant.
The region in the grey band at large masses is excluded by requiring efficient
tachyonic production during inflation while the region in gray at low masses
is excluded by requiring the dark photons are non-relativistic by the time of
CMB formation (see Eq.(39)). For the case when the dark vector has a mass
during inflation, we also show along the black line labeled “Longitudinal” the
contour where the longitudinal mode makes up all of the observed dark matter.

For comparison, we also plot the relic abundance of the longitudinal mode in
Eq.(38) in the case the dark vector has a mass during inflation and is thus
also produced via inflationary fluctuations Graham:2015rva . In the regions
where the line labelled ‘Transverse’ is to the left of the ones labelled
‘Longitudinal’, the transverse mode gives the dominant contribution to the
relic density. We see large regions of parameter space where this is the case.
On the left in Fig. 2 we also see a region of parameter space where the
longitudinal and transverse modes give comparable contributions to the dark
matter relic abundance which raises interesting possibilities to be discussed
more below. Since for the parameters shown on the right in Fig. 2 the
transverse mode always dominates we do not show the contour for the
longitudinal mode. Note however the relic abundance for the longitudinal mode
only depends on $m$ and $H$ so it has the same contour regardless of the other
parameters.

Finally, for illustration purposes we consider a specific benchmark point,

$\displaystyle\xi_{\rm end}=6\,,\ \ \epsilon_{R}=10^{-1}\,,\ \
\epsilon_{H}=10^{-1}\,,$ (46) $\displaystyle H=10^{9}\ {\rm GeV}\,,\ \ m=1.3\
{\rm keV}\,.$

This leads to a reheating temperature $T_{\rm RH}=2.7\times 10^{12}$ GeV and
an initial radiation energy density $\rho_{R}(T_{\rm RH})=1.7\times 10^{51}\
{\rm GeV}^{4}$, several orders of magnitude larger than the initial energy
density in the dark electromagnetic field $\rho_{D}(T_{\rm RH})=10^{42}\ {\rm
GeV}^{4}$. The dark photons become non-relativistic at $\bar{T}=36$ MeV, then
redshift like matter for some time before matching the energy density of
radiation at $T_{\rm CMB}$. Note that the momentum of the dark photon has a
long time to redshift from $\bar{T}$ to $T_{\rm CMB}$ so it is very ‘cold’ by
the time of matter-radiation equality. Note also that at the time of Big Bang
Nucleosynthesis (BBN), $T_{\rm BBN}\sim 1$ MeV, the dark photon is already
non-relativistic and still constitutes a small fraction of the total energy
density. Therefore, bounds on extra relativistic species ($N_{\rm eff}$) are
easily avoided.

##  III Energy Density Spectrum and Clumpy Dark Matter

Here we examine in detail the cosmological evolution of the dark vector energy
density starting from the end of inflation. In particular, we examine the
energy density spectrum and show explicitly that the peak in the spectrum at
the end of inflation survives cosmic evolution until late times. We then
examine the spectrum of _fluctuations_ in the energy density around the time
of matter radiation equality. We confirm that power at large scales remains
highly suppressed allowing for a vector produced in this manner to evade
constraints on isocurvature from measurements of the CMB Akrami:2018odb , thus
making it a viable dark matter candidate. The power spectrum of density
fluctuations also serves as the starting point for studying implications on
structure formation Kolb:1994fi ; Graham:2015rva ; Alonso-Alvarez:2018tus ;
Berges:2019dgr  which we leave for future work.

###  III.1 Energy density spectrum at the end of inflation

To gain further intuition for the tachyonic production mechanism and evolution
of the modes during inflation as well as facilitate the discussion below, we
can consider the analytic solutions in Eq.(14)-Eq.(17) together with the
conformal diagram in Fig. 3. Here we show co-moving length scales versus scale
factor (or conformal time) with the co-moving horizon indicated by the contour
(solid black) at $k^{-1}=(aH)^{-1}$. The last mode to exit the horizon during
inflation is indicated by the black line labeled
$k_{\rm{end}}^{-1}=(a_{\rm{end}}H_{\rm end})^{-1}$. The Compton wavelength
contour at $\lambda=q^{-1}=(k/a)^{-1}=m^{-1}$ applies in the case that the
dark vector already has a mass during inflation and defines the time when it
becomes non-relativistic. We also show the last mode to cross this contour
(and become non-relativistic) during inflation labeled $k_{m}$ as well as the
maximum momentum tachyonic mode $k_{\rm max}\simeq 2\xi k_{\rm max}$.
Utilizing Fig. 3 and Eq.(13) we find the ratios of these scales,

Figure 3: Conformal diagram during inflation zoomed in around the horizon
region (green) showing co-moving scales versus scale factor (see text for more
information).

$\displaystyle\frac{k_{m}}{k_{\rm{end}}}$ $\displaystyle=$
$\displaystyle\frac{m}{H_{\rm end}},\leavevmode\nobreak\ \leavevmode\nobreak\
\frac{k_{\rm{max}}}{k_{\rm{end}}}=2\xi_{\rm
end}+\mathcal{O}(\frac{m^{2}}{H_{\rm end}^{2}}),$ (47)

where $H_{\rm end}$ indicates the Hubble scale at the end of inflation. In the
super horizon limit ($-k\tau\to 0$) shown in blue the amplitudes approach the
asymptotic solution in Eq.(16). In the orange region deep inside the horizon
($-k\tau\to\infty$) the modes are in the Bunch-Davies vacuum given by the
solution in Eq.(14). As the modes enter the horizon region shown in green and
Eq.(13) is satisfied, the tachyonic instability leads to exponential growth of
the amplitude for one of the transverse modes (which we take to be $A_{+}$).
As we’ll see in Sec. III.1, for a fixed point in time (or scale factor), the
power in the dark electromagnetic field is dominated by modes contained within
this region leading to a peak in the power spectrum at $k\sim aH$ when the
modes have a wavelength of order the horizon. This leads to a power spectrum
at the end of inflation that is peaked at co-moving scales around
$k_{\rm{end}}^{-1}$.

The result for the total dark vector energy density at the end of inflation
obtained in Eq.(22) relied on the WKB approximation for the amplitude in
Eq.(17) which assumed the time dependence of $\xi$ could be neglected. This
was sufficient for obtaining an estimate of the relic density and viable dark
matter parameter space in Fig. 2. However, the instability parameter
$\xi\propto\dot{\phi}\propto\sqrt{\epsilon}$ is not only largest towards the
end of inflation, but also experiences the largest _growth_ just as inflation
is ending and the slow roll parameter approaches $\epsilon\approx 1$ (in
single field inflation scenarios). Since the energy density depends
exponentially on $\xi$ (see Eq.(22)), to obtainin an accurate density spectrum
it is crucial to account for this time dependence. However, to account for the
time dependence of $\xi$ and the Hubble parameter as well as the breakdown of
the slow roll approximation, it is necessary to numerically solve the
equations of motion in Eq.(5) and Eq.(6). This requires a robust integration
procedure which we describe in the Appendix.

Having performed this numerical integration, we show in Fig. 4 a comparison
between the energy density spectrum obtained numerically (solid) versus
analytically in Eq.(15) (dashed) for both the electric (blue) and magnetic
(red) components (see Eq.(20)). We consider a $\phi^{4}$ type inflaton
potential and $m/H_{\rm end}\ll 1$. On the left we show the spectra at early
times during inflation as the CMB modes leave the horizon where we require
$\xi$ to satisfy $\xi_{\rm{CMB}}<2.5$ Barnaby:2010vf ; Barnaby:2011vw . On the
right we show the same spectra, but now just at the end of inflation for
$\xi_{\rm end}=9$. As we can see, at early times during inflation when
$\xi_{\rm end}\approx 1$ the analytic and numerical solutions are in good
agreement. However, at the end of inflation we see the shape, location, and
exact height of the peak depends not just on $\xi_{\rm end}$ (which is the
same in both cases), but on how $\xi$ changes with time. We also see that the
time dependence in $\xi$ leads to stronger suppression of power at large
scales than the analytic case. Note this way of suppressing power at large
scales is distinct from other dark matter production mechanisms connected to
inflation which also lead to a peaked power spectrum Graham:2015rva ; Alonso-
Alvarez:2018tus ; Berges:2019dgr .

Figure 4: Energy density spectrum for the electric (E) and magnetic (B)
components in Eq.(20) obtained numerically (solid) versus analytically
(dashed) at both early times during inflation (left) as the CMB modes leave
the horizon and just at the end of inflation (right) for $\xi_{\rm end}=9$ and
$m=0$. The spectrum for the Bunch-Davies vacuum modes is also shown (dotted).

With the numerical solutions in hand, we obtain the energy density spectrum at
the end of inflation for various inflationary scenarios. We first examine the
effects of the dark vector having a mass already during inflation555In this
case there may be a constraint from the Swampland conjecture which requires
$m>60$ eV Bastero-Gil:2018uel ; Reece:2018zvv . which can arise through either
a Stueckelberg or Higgs mechanism. In Fig. 5 we show the energy density
spectrum for the electric (solid) and magnetic (dashed) component at the end
of inflation as a function of $k/k_{\rm{end}}$ for $m/H_{\rm
end}=0\leavevmode\nobreak\ (\rm{black}),\leavevmode\nobreak\ 3.15\cdot
10^{-5}\leavevmode\nobreak\ (\rm{red}),\leavevmode\nobreak\ 3.15\cdot
10^{-1}\leavevmode\nobreak\ (\rm{green})$ and $\xi_{\rm end}=9$. The spectrum
for the Bunch-Davies vacuum (black dotted) is also shown. We see all the
spectra are above the vacuum for modes with $k/a_{end}<2\xi_{\rm end}H_{\rm
end}$ showing the particle production effect. Higher momentum modes stay in
the Bunch-Davies vacuum and have a spectrum below the vacuum one once a proper
subtraction scheme has been implemented (see Appendix). This signals the
absence of the tachyonic instability and particle production effects for these
modes.

Figure 5: Left: Energy density spectrum for the electric (solid) and magnetic
(dashed) component at the end of inflation for $m/H_{\rm
end}=0\leavevmode\nobreak\ (\rm{black}),\leavevmode\nobreak\ 3.15\cdot
10^{-5}\leavevmode\nobreak\ (\rm{red}),\leavevmode\nobreak\ 3.15\cdot
10^{-1}\leavevmode\nobreak\ (\rm{green})$ and $\xi_{\rm end}=9$. The red and
green dotted vertical lines indicate $k_{m}/k_{\rm{end}}=m/H_{\rm end}$ (see
Fig. 3 and Eq.(47)). The spectrum for the Bunch-Davies vacuum (black dotted)
is also shown. Right: Same energy density spectra as left, but zoomed in
around the peak for $\xi_{\rm end}=3,6,9$ and $m/H_{\rm end}\ll 1$.

We also see a number of features which arise when the dark vector has a mass
already during inflation. For instance, we see at around
$k_{m}/k_{\rm{end}}=m/H_{\rm end}$ (see Fig. 3 and Eq.(47)) indicated by the
red and green dotted vertical lines, the slope in the spectrum for the
magnetic component changes from decreasing like $k^{-4}$ to one decreasing
like $k^{-2}$ as we go to larger scales. This change occurs when $q=m$ and the
modes at large co-moving scale (top in Fig. 3) cross the Compton wavelength
contour in Fig. 3 causing them to become non-relativistic and damp more slowly
with expansion than the still relativistic modes at smaller scales. As
discussed above, the last mode for which this occurs is at $k_{m}$ so modes at
scales larger than $k_{m}^{-1}$ will see an enhancement relative to the
massless case as seen in Fig. 5. We see also that at the end of inflation this
leads to domination by the magnetic component at scales $k^{-1}>k_{m}^{-1}$
where $m>>q$. At even larger scales $k^{-1}\gg k_{m}^{-1}$, we see the
electric component also changes slope from one decreasing like $k^{-4}$ to one
decreasing like $k^{-2}$ after going through a kink when the field time
derivative $\partial_{\tau}A_{+}$ changes sign. However, we see it still
remains subdominant to the magnetic component at these scales. This is in
contrast to the massless case in which the electric component dominates at all
scales666We discuss these mass effects on the spectrum at large scales in more
detail in the Appendix.. Finally, we see that around the peak the mass effects
are negligible and, in particular, the massive and massless cases have the
same spectrum for modes $k>k_{m}$ which will always contain the majority of
the peak. This is the case unless $H_{\rm end}\lesssim m<\xi H_{\rm end}$
which we do not consider since the tachyhonic production begins to be
suppressed. However, this super heavy mass case could be an interesting
possibility.

We next examine the dependence of the energy density spectrum on $\xi_{\rm
end}$. Since around the peak the electric component always dominates and
largely does not depend on the vector mass, we can focus on the region around
the peak and consider different values for $\xi_{\rm end}$. On the right in
Fig. 5 we show the same spectra as on the left, but only around the peak for
$\xi_{\rm end}=3,6,9$ and $m/H_{\rm end}\ll 1$. Here the exponential
sensitivity to $\xi_{\rm end}$ becomes clear as well as the domination of the
electric component of the energy density. Note also that only for the
$\xi_{\rm end}=9$ case does the energy density in the dark vector begin to
approach the energy density of the inflaton ($\frac{d\rho_{E}}{d\,{\rm
ln}\,k}\approx\rho^{\rm{end}}_{I}$) so we can safely neglect back reaction
effects.

Figure 6: Left:  Instability parameter $\xi$ as a function of number of
e-folds for various models of inflation. Right: Energy density spectrum at the
end of inflation for the same models with $\xi_{\rm end}=9$.

Finally, we examine how the shape of the inflaton potential affects the dark
vector energy density spectrum. In Fig. 6 we show on the left how the
instability parameter $\xi$ changes as a function of the number of e-folds for
various models of inflation while the energy density spectrum at the end of
inflation is shown on the right for the same models. We see that regardless of
the behaviour of $\xi$ early on during inflation, the largest growth occurs
just at the very end of inflation as $\xi$ approaches $\xi_{\rm end}$. We also
see that different inflationary potentials can lead to different spectra for
the same $\xi_{\rm end}$. This opens the possibility that by precisely
measuring the dark vector power spectrum, we can potentially infer properties
of the inflaton potential, but we leave an exploration of this possibility to
future work. These energy density spectra will serve as the input needed to
track the cosmic evolution of the dark vector energy density. Below we study
the cosmic evolution to obtain the late time energy density spectra and its
fluctuations.

###  III.2 Cosmological evolution of modes

Here we track the evolution of the tachyonic modes starting from the end of
inflation through to matter radiation equality following a similar analysis to
that found in Graham:2015rva ; Alonso-Alvarez:2018tus , keeping track of both
the (dark) magnetic and electric components as the Universe continues to
expand. As we saw above, when produced at the end of inflation, the electric
component dominates over the magnetic. However, as we will see below, the
magnetic component quickly grows relative to the electric after inflation and
‘catches up’ by the time the Hubble scale becomes comparable to the dark
matter mass. After this point they redshift together like radiation and then
eventually like matter by the time of matter radiation equality.

The cosmological evolution of the field as a function of scale factor or time
(either physical or conformal) can be parametrized in the same way as,

$\displaystyle\mathcal{P}_{X}(k,x)$ $\displaystyle=$
$\displaystyle\mathcal{P}_{X}(k,x_{\rm{end}})\frac{|X(k,x)|^{2}}{|X(k,x_{\rm{end}})|^{2}},\leavevmode\nobreak\
\leavevmode\nobreak\ X=A_{+},\,\partial_{\tau}A_{+},\leavevmode\nobreak\
\leavevmode\nobreak\ x=a,t,\tau,$ (48)

where the power spectrum is defined in Eq.(19) (with $\tau\to x$) and
$x_{\rm{end}}$ indicates the end of inflation. The input amplitude (plus
derivative) and power spectrum are taken at $x_{\rm end}$ and obtained by
numerically solving the equations of motion in Eq.(10) as discussed in Sec.
III.1 and the Appendix. The late time modes can similarly be obtain
numerically, but it becomes computationally intensive to evolve them to late
times. An approximate analytic solution for the late time amplitude and its
derivative can be found in different limits of the equations of motion in
Eq.(10). These different limits can be depicted geometrically via the
conformal diagram shown in Fig. 7 where each limit corresponds to one of the
five colored regions as indicated. As we discuss in more detail below, taking
the amplitude (plus derivative) and power spectrum at the end of inflation as
input, we can obtain approximations to $X(k,x)$ in the various limits and then
‘glue’ them together to construct the full late time amplitude as well as the
mean energy density and fluctuation spectrum Graham:2015rva ; Alonso-
Alvarez:2018tus .

Figure 7: Conformal diagram showing co-moving scales versus scale factor from
early during inflation until matter radiation equality (see text for more
information).

Using Fig. 7 we derive a second set of ratios of scales (in addition to
Eq.(47)),

$\displaystyle\frac{k_{\ast}}{k_{\rm{end}}}=\frac{a_{\rm{end}}}{a_{\ast}}$
$\displaystyle=$ $\displaystyle\sqrt{\bar{m}},\leavevmode\nobreak\
\leavevmode\nobreak\ \frac{a_{\rm{end}}}{\bar{a}}=\bar{m},\leavevmode\nobreak\
\leavevmode\nobreak\
\frac{a_{\rm{end}}}{a_{\rm{max}}}=\frac{\bar{m}}{2\xi_{\rm end}},$ (49)

where $a_{\ast}$ is defined as the scale factor when $H=m$ and the $k_{\ast}$
mode777As we discuss below, the peak in the power spectrum of the longitudinal
mode occurs at $k_{\ast}.$ becomes non-relativistic (crosses the Compton
wavelength contour in Fig. 7) while again $a_{\rm end}$ is the scale factor at
the end of inflation. The scale factors $\bar{a}$ and $a_{\rm max}$ (or
e-folds) indicate respectively the time when the modes $k_{\rm end}$ and
$k_{\rm max}$ become non-relativistic. Since we are taking the end of
inflation as the initial condition for cosmological evolution, unless
otherwise stated we define $\bar{m}$ as the ratio of the dark vector mass to
Hubble scale at the end of inflation,

$\displaystyle\bar{m}=\frac{m}{H_{\rm end}}.$ (50)

For the following discussion it will also be useful to rewrite the electric
and magnetic energy density spectra respectively as a function of scale factor
$a$,

$\displaystyle\frac{d\rho_{E}}{d\,{\rm ln}\,k}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{a_{\rm{end}}}{a}\right)^{4}H_{\rm{end}}^{2}\mathcal{P}_{\partial_{a}A_{+}}(k,a),\leavevmode\nobreak\
\leavevmode\nobreak\ \frac{d\rho_{B}}{d\,{\rm
ln}\,k}=\frac{1}{2a^{2}}\left(\frac{k^{2}}{a^{2}}+m^{2}\right)\mathcal{P}_{A_{+}}(k,a)\,,$
(51)

where we have used $da=a^{2}Hd\tau$ and $H=(a_{\rm{end}}/a)^{2}H_{\rm{end}}$
for a radiation dominated era.

  * •

Region I - End of Inflation: $H=H_{\rm end}$

During inflation (region I in Fig. 7) quantum fluctuations of the dark vector
boson are amplified by the expansion of the Universe leading to the tachyonic
instability and exponential production of one polarization. As discussed, this
production is maximal just at the end of inflation leading to a power spectrum
peaked at $k_{\rm end}$ that serves as input for cosmological evolution after
inflation (see Eq.(48)). Explicitly, at scale factor $a=a_{\rm end}$, we take
the tachyonically enhanced amplitude and its derivative to be,

$\displaystyle A_{\rm{I}}$ $\displaystyle=$ $\displaystyle
A_{+}(k,a_{\rm{end}}),$ (52) $\displaystyle\partial_{a}A_{\rm{I}}$
$\displaystyle=$ $\displaystyle\partial_{a}A_{+}(k,a)|_{a=a_{\rm{end}}},$

where $A_{+}$ and $\partial_{a}A_{+}$ are obtained numerically as discussed in
Sec. III.1. These will serve as the initial conditions for the cosmological
evolution of the tachyonic modes.

  * •

Region II - Super-horizon radiation era relativistic: $H\gg q\gg m$

Just after inflation we are in the super horizon regime $H\gg q,m$ for modes
with $am<k<k_{\rm end}$ (region II in Fig. 7). At this point, the inflaton
energy has been converted into reheating and thus, $\xi=0$ in the equation of
motion in Eq.(6). For these relativistic Hubble damped modes with $q\gg m$
this gives an equation of motion,

$\displaystyle(\partial_{t}^{2}+H\partial_{t})\,A_{+}\simeq
0,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\Longleftrightarrow\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ (\partial_{a}aH\partial_{a}+H\partial_{a})\,A_{+}\simeq
0,$ (53)

where we have used $\dot{a}=aH$ and $da=a^{2}Hd\tau$. Defining $A_{\rm{II}}$
as the approximate solution for the amplitude in this regime, during the
radiation dominated era when the Hubble parameter scales as $H\propto a^{-2}$
we have,

$\displaystyle A_{+}\simeq A_{\rm{II}}$ $\displaystyle=$ $\displaystyle
c^{\rm{II}}_{1}+c^{\rm{II}}_{2}a,$ (54) $\displaystyle\partial_{a}A_{\rm{II}}$
$\displaystyle=$ $\displaystyle c^{\rm{II}}_{2},$

where $c^{\rm{II}}_{1,2}$ are constants in time, but functions of scale
$k^{-1}$, and we see a term in the amplitude that grows linearly with scale
factor. Using Eq.(51), we see the electric and magnetic energy density scale
as,

$\displaystyle\frac{d\rho_{E}}{d\,{\rm ln}\,k}\,\propto
a^{-4},\leavevmode\nobreak\ \leavevmode\nobreak\ \frac{d\rho_{B}}{d\,{\rm
ln}\,k}\,\propto\frac{q^{2}}{a^{2}}(c^{\rm{II}}_{1}+c^{\rm{II}}_{2}a)^{2}\propto
a^{-(2-4)},$ (55)

where again $q=k/a$ is the physical momentum. However, imposing continuity of
the amplitude and derivative at $a_{\rm{end}}$ we also have,

$\displaystyle c^{\rm{II}}_{1}$ $\displaystyle=$ $\displaystyle
A_{\rm{I}}-a_{\rm end}(\partial_{a}A_{\rm{I}})$ $\displaystyle
c^{\rm{II}}_{2}$ $\displaystyle=$ $\displaystyle\partial_{a}A_{\rm{I}},$ (56)

where we have used Eq.(52). So we see the size of the term that grows linearly
with $a$ depends on the the size of the input amplitude derivative at the end
of inflation. Since this depends on $k$ and has a peak at $\sim k_{\rm{end}}$,
modes at small scales around $k_{\rm end}$ have a large $c^{\rm{II}}_{2}$.
However, they spend less time growing linearly before re-entering the horizon
while the converse is true for modes at larger scales (see Fig. 7). The net
effect is a brief period of damping like $a^{-2}$ which is much slower than
the $a^{-4}$ damping of the electric component. As we show in Fig. 8 and
discuss more below, this brief period of enhancement during super-horizon
evolution allows the magnetic component of the energy density to quickly
‘catch up’ to the electric component after which they redshift together. Note
that in the case of the longitudinal mode the factor of 3 in front of the
$\partial_{t}H$ Hubble damping term leads instead to a solution of the form
$A_{L}\simeq c_{1}+c_{2}a^{-1}$, which quickly leads to a constant as the
Universe expands.

  * •

Region III - Super-horizon radiation era non-relativistic: $H\gg m\gg q$

If the dark vector has a mass already during inflation or one is generated
while some modes are still super-horizon during radiation era, we have the
possibility of non-relativistic Hubble damped evolution (region III in Fig.
7). Here we have the same equation of motion as in II leading to the same
solutions, including the coefficients ($c^{\rm{III}}_{i}=c^{\rm{II}}_{i}$)
which at these scales $k\ll k_{\rm end}$ are very small. Note, this solution
has the same form as the longitudinal case which has the same equation of
motion in this region with differing coefficients due to the different input
spectrum at the end of inflation. Unlike in region II, the mass term in
Eq.(51) now dominates over the momentum term for the magnetic component of the
dark vector energy density spectrum. The electric and magnetic energy
densities then damp as,

$\displaystyle\frac{d\rho_{E}}{d\,{\rm ln}\,k}\,\propto
a^{-4},\leavevmode\nobreak\ \leavevmode\nobreak\ \frac{d\rho_{B}}{d\,{\rm
ln}\,k}\,\propto\frac{m^{2}}{a^{2}}(c^{\rm{II}}_{1}+c^{\rm{II}}_{2}a)^{2}\propto
a^{-(0-2)}.$ (57)

Since we are far from the peak, these contributions to the energy density are
negligible unless $m\sim H_{\rm end}$ which we do not consider as we assume
$m\ll H_{\rm end}$.

  * •

Region IV - Sub-horizon radiation era relativistic: $q\gg m,H$

This is the region in Fig. 7 containing the modes around the peak in the dark
vector energy density spectrum. In this regime we have for the equation of
motion,

$\displaystyle\Big{(}\partial_{t}^{2}+H\partial_{t}+\frac{k^{2}}{a^{2}}\Big{)}\,A_{+}\simeq
0\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\Longleftrightarrow\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \Big{(}\partial_{\tau}^{2}+k^{2}\Big{)}\,A_{+}\simeq 0,$
(58)

where we have used $\partial_{\tau}a^{-1}=-H$. Note this differs from the
equation of motion for the longitudinal component which in conformal time has
a $2aH\partial_{\tau}$ ‘damping’ term Graham:2015rva . This is due to a factor
of 3 in front of the $H\partial_{t}$ Hubble damping term in the physical time
equation of motion. Thus Eq.(58) has solution for the transverse mode,

$\displaystyle A_{+}\simeq A_{IV}$ $\displaystyle=$
$\displaystyle(c^{\rm{IV}}_{1}e^{ik\tau}+c^{\rm{IV}}_{2}e^{-ik\tau})$
$\displaystyle\partial_{a}A_{IV}$ $\displaystyle=$
$\displaystyle\frac{ik}{a_{\rm{end}}^{2}H_{\rm{end}}}(c^{\rm{IV}}_{1}e^{ik\tau}-c^{\rm{IV}}_{2}e^{-ik\tau}),$
(59)

which we see has no overall damping. The longitudinal mode amplitude on the
other hand has an overall $a^{-1}$ suppression and the derivative damps like
$a^{-2}$. Note that the factor of $a_{\rm end}^{2}H_{\rm end}$ comes from the
Jacobian in going from $\partial_{\tau}$ to $\partial_{a}$. For the transverse
mode the electric and magnetic energy density then scale as,

$\displaystyle\frac{d\rho_{E}}{d\,{\rm ln}\,k}\,\propto
a^{-4},\leavevmode\nobreak\ \leavevmode\nobreak\ \frac{d\rho_{B}}{d\,{\rm
ln}\,k}\propto a^{-2}q^{2}\propto a^{-4}.$ (60)

We see that the magnetic and electric components of the energy density
redshift the same with scale factor in this regime. Note that while the
amplitude scales differently, the damping of the energy density is the same as
for the longitudinal mode Graham:2015rva .

  * •

Region V - Non-relativistic massive regime: $m\gg H,q$

In the non-relativistic massive regime we have $m\gg H,q$ which gives for the
equations of motion (now in terms of physical time, see Eq.(6)),

$\displaystyle\Big{(}\partial_{t}^{2}+H\partial_{t}+m^{2}\Big{)}\,A_{+}\simeq
0.$ (61)

This is the same as for the longitudinal component Graham:2015rva  and has
solution,

$\displaystyle A_{+}\simeq A_{V}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{a}}(c^{V}_{1}e^{imt}+c^{V}_{2}e^{-imt})$
$\displaystyle\partial_{a}A_{V}$ $\displaystyle=$
$\displaystyle\frac{im\sqrt{a}}{a_{\rm{end}}^{2}H_{\rm{end}}}(c^{V}_{1}e^{imt}-c^{V}_{2}e^{-imt}),$
(62)

where we have used the change of variables from physical time to scale factor,

$\displaystyle
t\to\frac{1}{2H_{\rm{end}}}(\frac{a^{2}}{a_{\rm{end}}^{2}}-1)+t_{\rm{end}}.$
(63)

We see the amplitude damps like $a^{-1/2}$ while the derivative grows like
$a^{1/2}$ leading to a scaling for the electric and magnetic energy densities
respectively,

$\displaystyle\frac{d\rho_{B}}{d\,{\rm ln}\,k}$ $\displaystyle\propto$
$\displaystyle\frac{m^{2}}{a^{3}}\propto a^{-3},\leavevmode\nobreak\
\leavevmode\nobreak\ \frac{d\rho_{E}}{d\,{\rm ln}\,k}\propto\frac{a_{\rm
end}^{4}H_{\rm end}^{2}}{a^{3}}\,\propto a^{-3}.$ (64)

Thus both the magnetic and electric components of the energy density redshift
like matter at late times. In this regime the longitudinal component has the
same equation of motion and therefore same solution as in Eq.(• ‣ III.2). So
again we have the same damping behaviour with scale factor Graham:2015rva ,
but with different coefficients in Eq.(• ‣ III.2).

###  III.3 Late time energy density spectrum

Figure 8: On the left and right we show the electric and magnetic energy
density spectra respectively (near the peak) for various numbers of e-folds
after inflation, $N_{e}=0\,({\rm dashed\leavevmode\nobreak\
black}),N_{\ast}/10\,({\rm blue}),N_{\ast}\,({\rm orange}),N_{\rm{max}}\,({\rm
light\leavevmode\nobreak\ brown})$ where $a=a_{\rm end}e^{N_{e}}$ defines the
scale factor (see Fig. 7). We have normalized the spectrum to the energy
density of the inflaton at the end of inflation for $\xi_{\rm{end}}\simeq
9,\bar{m}\simeq 10^{-5}$ with a $\phi^{4}$ inflaton potential and factored out
an $a^{4}$.

In Fig. 8 and Fig. 9 we summarize the evolution of the energy density spectrum
after inflation focusing on modes with co-moving momentum $k_{\ast}\lesssim
k\lesssim k_{\rm max}$ which contain the vast majority of the power. In Fig. 8
on the left and right we show the electric and magnetic energy density spectra
respectively while in Fig. 9 we show their sum for various e-folds after
inflation with $\xi_{\rm end}\simeq 9,\bar{m}\simeq 10^{-5}$ and a
$V\propto\phi^{4}$ inflaton potential. The oscillatory behaviour seen in Fig.
8 arises from the oscillatory solutions for the modes (see Eq.(59) and Eq.(• ‣
III.2)) in regions IV and V of Fig. 7. For modes $k_{\ast}<k<k_{\rm end}$ this
occurs once they re-enter the horizon while modes with $k_{\rm end}<k<k_{\rm
max}$, which contain the peak of the power spectrum, approach but never exit
the horizon during inflation. Modes with $k<k_{\ast}$ at $N>N_{\ast}$ remain
outside the horizon for longer than modes with $k>k_{\ast}$ (see Fig. 7). Thus
they have not had enough time to begin oscillating upon reentering the horizon
so we see a still smooth spectrum in this regime. For the input values
$\xi_{\rm end}=9$ and $\bar{m}=10^{-5}$ this gives in terms of e-folds
$(N_{e})$ after inflation $N_{e}=N_{\ast}\simeq 5.65,\,\bar{N}\simeq
11.30,\,N_{\rm{max}}\simeq 14.19$ where $a=a_{\rm end}e^{N_{e}}$. There are a
number of features which are evident that reflect the behavior of the modes in
the different limits of the equations of motion discussed in Sec. III.2 and
can be understood with the help of Fig. 7 and Eq.(47), Eq.(49).

Figure 9: Total energy density spectra (electric plus magnetic shown in Fig.
8) for $N_{e}=0\,({\rm dashed\leavevmode\nobreak\ black}),N_{\ast}\,({\rm
red}),\bar{N}\,({\rm blue}),N_{\rm{max}}\,({\rm
orange}),1.5\,N_{\rm{max}}\,({\rm grey})$ where $a=a_{\rm end}e^{N_{e}}$ and
again normalized to $\rho^{\rm{end}}_{I}$ with $a^{4}$ factored out. See text
for more information.

Looking first at the electric energy density spectrum on the left in Fig. 8,
we normalize to the energy density of the inflaton at the end of inflation and
divide by $a^{-4}$ in order to compare to radiation like damping. We see that
once the input electric energy density spectrum is set at the end of inflation
(black dashed) and which dominates over the magnetic component, the electric
energy density then redshifts like radiation until $a_{\ast}$ (orange). At
this point scales larger than $k^{-1}_{\ast}$ begin redshifting like matter
while modes at smaller scales around the peak continue redshifting like
radiation. Modes around the peak start becoming non-relativistic, and thus
relatively enhanced as they begin redshifting like matter, once $q\leq m$
which happens at progressively later times for modes with larger and larger
momenta as can be understood geometrically in Fig. 7. The last mode to exit
the horizon during inflation, $k_{\rm end}=a_{\rm end}H_{\rm end}$, becomes
non-relativistic at $\bar{a}$ while the highest momentum tachyonic mode,
$k_{\rm max}=2\xi\,k_{\rm end}$, becomes non-relativistic at $a_{\rm{max}}$
(light brown). After this point, all of the modes redshift together like
matter with $\propto a^{-3}$ damping and the shape of the electric energy
density spectrum no longer changes.

On the right of Fig. 8 we examine the magnetic energy density spectrum which
has a more interesting evolution. We again normalize to the energy density of
the inflaton at the end of inflation and factor out an $a^{4}$. We see that
after initially being sub dominant to the electric energy density at the end
of inflation (black dashed), there is a brief period where modes around the
peak damp like $\propto a^{-(2-3)}$ (blue) as compared to $\propto a^{-4}$
damping for the electric component. As discussed in Sec. III.2, this brief
period of relative enhancement is due to the linear growth that the mode
amplitudes in region II experience while they are super horizon (see Eq.(54))
combined with continuity of the amplitude and derivative at the various
boundaries. After this brief period, the modes around the peak quickly begin
redshifting like radiation with the usual $a^{-4}$ damping well before the
time they reach $a_{\ast}$ (orange). It is this very brief period of growth
relative to matter, typically lasting around an e-fold or less after
inflation, that allows the magnetic energy density to “catch up” to the
electric by the time $H=m$ at $a_{\ast}$ after which they redshift together
like radiation and then eventually like matter (see Fig. 7).

In Fig. 9 we show the total dark vector energy density spectrum given by the
sum of the electric and magnetic components. We see that when summed the
oscillations of the electric and magnetic component cancel one another to give
a smooth spectrum as can be shown analytically using Eq.(59) or Eq.(• ‣ III.2)
and plugging it into Eq.(51). As we can also see, the total dark vector energy
density spectrum essentially evolves like (the envelope of) the electric
component. This is because the electric component dominates at the end of
inflation when the dark vector is produced and then redshifts along with the
magnetic component once the magnetic component catches up as discussed above.
Thus to a good approximation the total energy density evolution has the same
qualitative behavior as the electric component. In particular, after inflation
ends modes around the peak redshift like radiation until $q=m$ after which
they redshift like matter (see Fig. 1) as the mass term dominates the dark
vector energy density.

As we see from the energy density spectrum in Fig. 8, the vast majority of the
power is located at scales $\sim k_{\rm end}^{-1}$ around the peak which are
vastly smaller than the scales probed by CMB measurements. To see this
explicitly we write for modes around the peak $\sim k_{\rm end}^{-1}$,

$\displaystyle 1/k_{\rm end}$ $\displaystyle=$ $\displaystyle(a_{\rm
end}H_{\rm end})^{-1}=(a_{\rm RH}\epsilon_{H}H)^{-1}=(\frac{T_{0}}{T_{\rm
RH}}\epsilon_{H}H)^{-1}$ (65) $\displaystyle\approx$
$\displaystyle\frac{10^{-1}}{T_{0}}\,\frac{\epsilon_{R}}{\epsilon_{H}}\left(\frac{M_{\rm
Pl}}{H}\right)^{1/2}\approx 10\ {\rm km}\
\frac{\epsilon_{R}}{\epsilon_{H}}\left(\frac{100\ {\rm
GeV}}{H}\right)^{1/2}\,,$

where we have assumed $a_{\rm end}=a_{\rm RH}$ and used Eq.(28) as well as
$a_{\rm RH}/a_{0}=T_{0}/T_{\rm RH}$ with the scale factor today set to
$a_{0}=1$. In the absence of an extreme hierarchy between $\epsilon_{R}$ and
$\epsilon_{H}$, we see that the typical co-moving scale associated with the
peak is $\lesssim 10\,\rm{km}$ which is tiny on cosmological scales and thus
we expect isocurvature to be negligible on the large scales relevant for the
CMB. However, since it is relevant for how the dark matter is distributed
spatially after matter radiation equality and the evolution of density
perturbations, below we compute the power spectrum of isocurvature
perturbations. If the vector has a mass already during inflation, there could
be regions of parameter space where the longitudinal and transverse components
give comparable contributions to the dark matter energy density. As discussed
below, in this case there would be a double peaked structure in the energy
density spectrum with one peak corresponding to the transverse component
located at co-moving momenta $\sim k_{\rm{end}}$ and a second one
corresponding to the longitudinal mode at $k_{\ast}$ Graham:2015rva .

###  III.4 Isocurvature and density contrast power spectrum

Up until now, we have considered only the mean energy density and spectra
which can be written in terms of the dark vector field power spectrum as given
in Eq.(18). However, the dark vector energy density is subject to fluctuations
which can be of the same order as the mean energy density. These fluctuations
will be independent of the inflaton ones which set the curvature (or
adiabatic) perturbations that are imprinted on the metric. The fluctuations in
the dark vector energy density will therefore contribute to isocurvature
perturbations and can have implications for the CMB as well as structure
formation and clumping Kolb:1994fi  of the dark matter. Since measurements of
the CMB Akrami:2018odb  severely constrain the amplitude of isocurvature
perturbations, we must ensure that they are suppressed on long length scales.
Following Graham:2015rva ; Alonso-Alvarez:2018tus , here we compute these
isocurvature perturbations and demonstrate explicitly that they are highly
suppressed at the large scales relevant for the CMB.

The starting point is the density contrast field $\delta(\vec{x})$ which
describes deviations from the mean dark vector energy density
$\langle\rho\rangle$ and is defined as,

$\displaystyle\rho(\vec{x})=\langle\rho\rangle(1+\delta(\vec{x})).$ (66)

Near the time of matter radiation equality and once the vector dark matter
begins redshifting like matter, we can describe the energy density via the
mass term in the magnetic component. As we derive in the Appendix, the Fourier
transform of $\delta(\vec{x})$ can be written in terms of products of the
power spectrum for the tachyonically enhanced transverse mode as,

$\displaystyle{\cal P}_{\delta}(k,t)$
$\displaystyle=\frac{k^{2}}{\left[\int_{0}^{\infty}\frac{dk^{\prime}}{k^{\prime}}{\cal
P}_{A_{+}}(k^{\prime},t)\right]^{2}}\int_{0}^{\infty}dq\int_{|q-k|<p<q+k}dp\
\frac{1}{q^{2}p^{2}}{\cal P}_{A_{+}}(p,t){\cal P}_{A_{+}}(q,t)\,,$ (67)

where we have taken the mass term in Eq.(18) to represent the energy density
at late times and defined the power spectrum in terms of the 2-point function
of a random variable $X$,

$\displaystyle\langle
X(\vec{k})X(\vec{k}^{\prime})\rangle=(2\pi)^{3}\delta^{3}(\vec{k}+\vec{k}^{\prime})\frac{2\pi^{2}}{k^{3}}{\cal
P}_{X}(k)\,,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
X\equiv A_{+},\delta.$ (68)

We see in Eq.(67) that the power spectrum for $\delta$ corresponding to the
transverse vector mode 888We note that Eq.(67) agrees with the result found in
Alonso-Alvarez:2018tus  for _scalar_ dark matter produced during inflation.
differs from the one corresponding to the longitudinal mode Graham:2015rva .

Performing the integral in Eq.(67) numerically we show the density contrast
power spectrum in Fig. 10 for $\bar{m}\ll 1$ and $\xi_{\rm end}\approx 9$. We
see explicitly that the spectrum falls off sharply at large scales. We also
see that the location of the peak does not shift noticeably from the peak in
the mean energy density spectrum. At scales relevant for the CMB we see the
power is completely negligible though for $m\sim H_{\rm end}$ and
$\xi_{\rm{CMB}}\sim 2.5$ perhaps these isocurvature perturbations may have
observable consequences, but we do not explore this possibility here. To
determine the matter distribution and scale of clumping of the vector dark
matter today, one needs to follow the evolution of these density perturbations
through the non-linear regime until today. The power spectrum in Fig. 10
serves as the input for this evolution Berges:2019dgr , beginning at matter
radiation equality, but further investigation is left to future work.

Figure 10: The late time density contrast power spectrum for $\bar{m}\ll 1$
and $\xi_{\rm end}\approx 9$.

Finally, we comment that in the mechanism presented here, isocurvature for the
transverse mode is suppressed due to the fact that the tachyonic modes
experience maximal growth when they are of order the horizon combined with the
time dependence of the tachyonic instability parameter $\xi$ _during_
inflation. This is in contrast to the mechanism for suppressing isocurvature
in the energy density of the longitudinal mode 999Note this is also distinct
from how isocurvature is suppressed on long length scales for scalar dark
matter production mechanisms Alonso-Alvarez:2018tus ; Berges:2019dgr ;
Markkanen:2018gcw  which are also connected to inflation. which is due to how
the different modes redshift _after_ inflation Graham:2015rva .

###  III.5 Clumping of vector dark matter

As discussed in other ‘clumpy’ dark matter scenarios Graham:2015rva ; Alonso-
Alvarez:2018tus ; Berges:2019dgr , the scale corresponding to the peak of the
power spectrum of energy density fluctuations also has implications for the
scale on which the dark matter ‘clumps’. In one of these scenarios, a dark
vector has a mass during inflation and the longitudinal mode is necessarily
produced by inflationary fluctuations Graham:2015rva . In this case, due to
how the modes redshift after inflation, a peak is produced in the energy
density spectrum at $k_{\ast}$ (see Fig. 7) where,

$\displaystyle 1/k_{\ast}\sim 10^{10}\,{\rm km}\times\sqrt{\frac{10^{-5}\,{\rm
eV}}{m}}.$ (69)

We see that for the longitudinal mode the location of the peak only depends on
the dark vector mass which is constrained to be $m\geq 10^{-5}\,{\rm eV}$ if
it makes up the entirety of the dark matter Graham:2015rva . Considering the
range of Hubble scales during inflation $10^{2}<H/{\rm GeV}<10^{14}$ this
leads to the range of ‘clumping’ scales for the longitudinal component,

$\displaystyle 10^{-1}\,{\rm km}\ <k_{\ast}^{-1}<10^{10}\,{\rm km}\,.$ (70)

In the case of the transverse component produced via tachyonic instability, we
saw in Eq.(65) that the location of the peak does not depend on the dark
vector mass, but instead on the Hubble scale during inflation as well as
$\epsilon_{R}$ and $\epsilon_{H}$. For the same range of Hubble scales during
inflation $10^{2}<H/{\rm GeV}<10^{14}$ and taking
$\epsilon_{H}/\epsilon_{R}\sim 1$ we find the transverse component instead
clumps on scales in the range,

$\displaystyle{\rm cm}<k_{\rm end}^{-1}<10\,{\rm km}.$ (71)

The locations of the two peaks are related by the Hubble scale and dark vector
mass,

$\displaystyle\frac{k_{\rm
end}^{-1}}{k_{*}^{-1}}=\frac{\epsilon_{R}}{\epsilon_{H}}\left(\frac{m}{H}\right)^{1/2}\sim\left(\frac{m}{H}\right)^{1/2}\,.$
(72)

We see in general the transverse mode clumps on much smaller scales than the
longitudinal mode though when $m\sim H$ the scales can be comparable. As
discussed, this also opens the possibility of a double peak in the power
spectrum when both the longitudinal and transverse components contribute
appreciably to the relic abundance and would imply dark matter clumping on two
different scales. This ‘doubly clumpy’ possibility can of course only occur if
the dark vector has a mass already during inflation and thus would constitute
a striking signal of vector dark matter with an inflationary origin.

##  IV Summary and Outlook

In this study we have examined in detail the recently proposed mechanism
Bastero-Gil:2018uel  for producing non-thermal dark photon dark matter at the
end of inflation. This mechanism can generate the observed dark matter relic
abundance for dark vector masses in the range $\mu\,{\rm eV}\lesssim m\lesssim
10\,{\rm TeV}$ and Hubble scales during inflation in the range $100\,{\rm
GeV}\lesssim H\lesssim 10^{14}\,{\rm GeV}$. We have focused in particular on
the case where the dark vectors are relativistic at the time their mass is
generated and examined the associated cosmic evolution to compute the relic
abundance today. We have also examined the power spectrum and cosmic evolution
of the dark vector modes demonstrating explicitly that the late time spectrum
preserves the peak generated at the end of inflation. We have shown that the
peak corresponds to small physical scales today, $\ell_{\rm today}\sim{\rm
cm}-100\,{\rm km}$, with large density fluctuations at $\ell_{\rm today}$
implying a clumpy nature for the vector dark matter. The case of a non-
relativistic dark vector at the time its mass is generated has been left to
forthcoming work.

There are a number of interesting avenues to explore the phenomenology
associated with the dark photon dark matter production mechanism presented
here. If there are other dark charged particles present during inflation they
can potentially be produced via a dark Schwinger production mechanism. In this
case dark charged fermions and scalars may also contribute to the final dark
matter relic abundance. If the dark vector obtains its mass via a dark Higgs
mechanism, this can lead to an alternative cosmic evolution to the one
explored in this study and there may also be a possibility of generating
gravitational wave signals associated with the dark Higgs phase transition
Breitbach:2018ddu . Allowing for non-zero (but small) kinetic mixing between
the dark and visible photons can lead to interesting dark matter phenomenology
as has been thoroughly explored in many studies. Whether any of the initial
polarization of the dark vector survives cosmic evolution and has observable
effects is also worth investigating. Explorations of these various interesting
possibilities are ongoing.

_Acknowledgements_ : We thank Prateek Agrawal, Diego Blas, Adam Falkowski,
Bohdan Grzadkowski, Takeshi Kobayashi, Eric Madge, Manuel Masip, Gilad Perez,
Maxim Pospelov, Jennifer Schober, Pedro Schwaller, Javi Serra, Anna Socha,
Walter Tangarife, Tomer Volansky, and Tien-Tien Yu for useful comments and
discussions. This work has been partially supported by MINECO grants
PID2019-106087GB-C22, including ERDF (J.S., R.V.M.), PID2019-105943GB-I00
(M.B.G.), Junta de Andalucía Projects FQM-101, A-FQM-211-UGR18, P18-FR-4314,
(fondos FEDER), and SOMM17/6104/UGR (M.B.G., J.S., R.V.M.). L.U. acknowledges
support from the PRIN project “Search for the Fundamental Laws and
Constituents” (2015P5SBHT_002). R.V.M. would also like to acknowledge the
Mainz Institute for Theoretical Physics (MITP) of the Cluster of Excellence
PRISMA+ (Project ID 39083149) for their hospitality and partial support as
well as participants of _The Mysterious Universe_ workshop for useful
discussions and stimulating atmosphere while part of this work was completed.

## Appendix

###  IV.1 Derivation of equations of motion, energy density, and pressure

Here we derive the equations of motion as well as the energy and pressure
densities.

_Conventions_

Before presenting the derivation we define the conventions used. For the
metric we have,

$\displaystyle
ds^{2}=-dt^{2}+a^{2}(t)d\vec{x}^{2}=a^{2}(\tau)(-d\tau^{2}+d\vec{x}^{2})\,,$
(73) $\displaystyle g_{\mu\nu}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-1&0\\\
0&a^{2}(t)\delta_{ij}\end{pmatrix}=a^{2}(\tau)\begin{pmatrix}-1&0\\\
0&\delta_{ij}\end{pmatrix}\,,$ $\displaystyle g^{\mu\nu}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-1&0\\\
0&a^{-2}(t)\delta^{ij}\end{pmatrix}=a^{-2}(\tau)\begin{pmatrix}-1&0\\\
0&\delta^{ij}\end{pmatrix}\,,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $ (74)

with $t$ the cosmic time and $\tau$ the conformal time. The Levi-Civita tensor
is given by,

$\displaystyle\epsilon^{\mu\nu\rho\sigma}=\frac{\tilde{\epsilon}^{\mu\nu\rho\sigma}}{\sqrt{-g}}\,,\quad\epsilon_{\mu\nu\rho\sigma}=\sqrt{-g}\
\tilde{\epsilon}_{\mu\nu\rho\sigma}\,,$ (75)

with the following conventions for the anti-symmetric Levi-Civita symbol and
metric,

$\displaystyle\tilde{\epsilon}^{0123}=+1\,,\quad\tilde{\epsilon}_{0123}=-1\,,$
$\displaystyle\sqrt{-g}\equiv\sqrt{-\det(g_{\mu\nu})}=a^{3}(t)=a^{4}(\tau)\,.$
(76)

The three-dimensional Levi-Civita symbol is related to the four-dimensional
one as

$\displaystyle\tilde{\epsilon}_{ijk}=\tilde{\epsilon}^{ijk}=\tilde{\epsilon}^{0ijk}\,.$
(77)

For the Hubble parameter we have in terms of the scale factors,

$\displaystyle
H\equiv\frac{1}{a(t)}\frac{da(t)}{dt}=\frac{\dot{a}(t)}{a(t)}\,,\qquad{\cal
H}\equiv\frac{1}{a(\tau)}\frac{da(\tau)}{d\tau}=\frac{a^{\prime}(\tau)}{a(\tau)}\,,$
(78)

where we use an overdot for the cosmic time derivative and a prime for the
conformal time derivative. We work in comoving momentum space where for the
classical field $A_{\mu}$ we have:

$\displaystyle A_{0}(\vec{x},\tau)$
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}A_{0}(\vec{k},\tau)e^{i\vec{k}\cdot\vec{x}}\,,$
(79) $\displaystyle A_{i}(\vec{x},\tau)$
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}A_{i}(\vec{k},\tau)e^{i\vec{k}\cdot\vec{x}}\,.$
(80)

As the gauge field is real, $A_{\mu}(\vec{x},\tau)=A^{*}_{\mu}(\vec{x},\tau)$
which implies $A^{*}_{\mu}(\vec{k},\tau)=A_{\mu}(-\vec{k},\tau)$. We decompose
$\vec{A}(\vec{k},\tau)$ along its transverse ($\vec{A}_{T}$) and longitudinal
components ($A_{L}$),

$\displaystyle\vec{k}\cdot\vec{A}_{T}$ $\displaystyle=0$ (81)
$\displaystyle\vec{k}\cdot\vec{A}$ $\displaystyle=|\vec{k}|A_{L}\equiv
kA_{L}\,,$ (82)

and further decompose the transverse modes into the usual two polarizations

$\vec{A}_{T}(\vec{k},\tau)=\sum_{\lambda=\pm}\vec{\epsilon}_{\lambda}(\vec{k})\left[A_{\lambda}(k,\tau)+A_{\lambda}(k,\tau)^{*}\right]\,,$
(83)

where the polarization vectors satisfy,

$\displaystyle\vec{k}\cdot\vec{\epsilon}_{\pm}(\vec{k})=0\,,\quad\vec{k}\times\vec{\epsilon}_{\pm}(\vec{k})=\mp
ik\vec{\epsilon}_{\pm}(\vec{k})\,,$ (84)
$\displaystyle\vec{k}\cdot\vec{\epsilon}_{L}(\vec{k})=k\,,\quad\vec{k}\times\vec{\epsilon}_{L}(\vec{k})=0\,,$
(85)
$\displaystyle\vec{\epsilon}_{\lambda}(\vec{k})^{*}=\vec{\epsilon}_{\lambda}(-\vec{k})\,,\quad\vec{\epsilon}_{\lambda}(\vec{k})\cdot\vec{\epsilon}_{\lambda^{\prime}}(-\vec{k})=\delta_{\lambda\lambda^{\prime}}\,.$
(86)

The power spectra associated to the two-point correlation functions of the
classical field is,

$\displaystyle\langle\vec{A}(\vec{k},t)\cdot\vec{A}(\vec{k}^{\prime},t)\rangle$
$\displaystyle=(2\pi)^{3}\delta^{3}(\vec{k}+\vec{k}^{\prime})\frac{2\pi^{2}}{k^{3}}{\cal
P}_{A}(k,t)\,,$
$\displaystyle\langle\partial_{t}\vec{A}(\vec{k},t)\cdot\partial_{t}\vec{A}(\vec{k}^{\prime},t)\rangle$
$\displaystyle=(2\pi)^{3}\delta^{3}(\vec{k}+\vec{k}^{\prime})\frac{2\pi^{2}}{k^{3}}{\cal
P}_{\partial_{t}A}(k,t)\,,$ (87)

where we can then write for the two point function in position space,

$\displaystyle\langle\vec{A}(\vec{x},t)^{2}\rangle=\int_{0}^{\infty}\frac{dk}{k}{\cal
P}_{A}(k,t)=\int\frac{d^{3}k}{(2\pi)^{3}}\frac{2\pi^{2}}{k^{3}}{\cal
P}_{A}(k,t)\,.$ (88)

For the quantum field we expand in terms of creation and annihilation
operators,

$\displaystyle\hat{\vec{A}}(\vec{x},t)=\sum_{\lambda=\pm,L}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}}\
\vec{\epsilon}_{\lambda}(\vec{k})\
[A_{\lambda}(k,t)a_{\lambda}(\vec{k})+A_{\lambda}(k,t)^{*}a_{\lambda}^{\dagger}(-\vec{k})]\,.$
(89)

where the creation and annihilation operators satisfy,

$\displaystyle a_{\lambda}(\vec{k})|0\rangle=0\,,\quad\langle
0|a_{\lambda}^{\dagger}(\vec{k})=0\,,\quad\left[a_{\lambda}(\vec{k}),a^{\dagger}_{\lambda^{\prime}}(\vec{k}^{\prime})\right]=(2\pi)^{3}\delta_{\lambda\lambda^{\prime}}\delta^{3}(\vec{k}-\vec{k}^{\prime})\,.$
(90)

Two point correlation functions with the quantum field are then obtained by
sandwiching the field operators between vacuum states,

$\displaystyle\langle\hat{\vec{A}}(\vec{x},t)^{2}\rangle$
$\displaystyle\equiv$ $\displaystyle\langle
0|\hat{\vec{A}}(\vec{x},t)^{2}|0\rangle$ (91) $\displaystyle=$
$\displaystyle\sum_{\lambda,\lambda^{\prime}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{d^{3}k^{\prime}}{(2\pi)^{3}}e^{i(\vec{k}+\vec{k}^{\prime})\cdot\vec{x}}\vec{\epsilon}_{\lambda}(\vec{k})\cdot\vec{\epsilon}_{\lambda^{\prime}}(\vec{k}^{\prime})$
$\displaystyle\times$ $\displaystyle\langle
0|(A_{\lambda}(k,t)a_{\lambda}(\vec{k})A_{\lambda^{\prime}}(k^{\prime},t)^{*}a^{\dagger}_{\lambda^{\prime}}(-\vec{k}^{\prime}))|0\rangle$
$\displaystyle=$
$\displaystyle\sum_{\lambda}\int\frac{d^{3}k}{(2\pi)^{3}}|A_{\lambda}(k,t)|^{2}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $

where to go from second to third line we used the commutation relations.
Comparing to the definition of the power spectrum in Eq.(88) we find in the
quantum case,

$\displaystyle{\cal
P}_{\hat{A}}(k,t)=\frac{k^{3}}{2\pi^{2}}\sum_{\lambda}|A_{\lambda}(k,t)|^{2}\,.$
(92)

Similarly we can compute the power spectrum for the field time derivative
$\langle(\partial_{t}\hat{\vec{A}}(\vec{x},t))^{2}\rangle$,

$\displaystyle\langle(\partial_{t}\hat{\vec{A}}(\vec{x},t))^{2}\rangle$
$\displaystyle=$
$\displaystyle\sum_{\lambda}\int\frac{d^{3}k}{(2\pi)^{3}}|\partial_{t}A_{\lambda}(k,t)|^{2}\,,$
$\displaystyle{\cal P}_{\partial_{t}\hat{A}}(k,t)$ $\displaystyle=$
$\displaystyle\frac{k^{3}}{2\pi^{2}}\sum_{\lambda}|\partial_{t}A_{\lambda}(k,t)|^{2}\,.$
(93)

Note that the time dependence in the Fourier expansion is only in the function
$A_{\lambda}(k,t)$.

_Equations of motion_

To derive the equations of motion we consider the following action,

$\displaystyle S=\int d^{4}x\,L$ $\displaystyle=$ $\displaystyle\int
d^{4}x\sqrt{-g}\Big{[}-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-V(\phi)-\frac{1}{4}g^{\mu\nu}g^{\rho\sigma}F_{\mu\rho}F_{\nu\sigma}$
(94) $\displaystyle-$
$\displaystyle\frac{\alpha}{4f}\phi\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}-\frac{1}{2}m^{2}g^{\mu\nu}A_{\mu}A_{\nu}\Big{]}\,,$

where the field strength is defined in the usual way as,

$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}\,.$ (95)

Note that we have inserted the mass term for the vector field à la Proca so
this Lagrangian describes a model with a scalar field and a massive vector
field and is strictly speaking not a gauge theory, but this distinction is not
important for present purposes. Note also that,

$\displaystyle\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}$
$\displaystyle=$
$\displaystyle\frac{4}{\sqrt{-g}}\left[(\partial_{\tau}\vec{A}-\nabla
A_{0})\cdot(\nabla\times\vec{A})\right]$ $\displaystyle\rightarrow$
$\displaystyle\frac{4}{\sqrt{-g}}\int\frac{d^{3}k\
d^{3}k^{\prime}}{(2\pi)^{6}}e^{i(\vec{k}+\vec{k}^{\prime})\cdot\vec{x}}(\partial_{\tau}\vec{A}(\vec{k},\tau)-i\vec{k}A_{0}(\vec{k},\tau))\cdot(i\vec{k}^{\prime}\times\vec{A}(\vec{k}^{\prime},\tau))\,,$

where we see that this vanishes for the longitudinal mode for which
$\vec{k}^{\prime}\times\vec{A}(\vec{k}^{\prime},\tau)=0$. Thus, the term with
$F\tilde{F}$ only affects the transverse modes as we will see explicitly
below.

The equations of motion are obtained from the lagrangian density via,

$\displaystyle\partial_{\alpha}\frac{\delta
L}{\delta(\partial_{\alpha}\phi)}-\frac{\delta L}{\delta\phi}$
$\displaystyle=0\,,$ (97) $\displaystyle\partial_{\alpha}\frac{\delta
L}{\delta(\partial_{\alpha}A_{\beta})}-\frac{\delta L}{\delta A_{\beta}}$
$\displaystyle=0\,.$ (98)

From the first we equation obtain in cosmic time,

$\displaystyle\ddot{\phi}-\frac{1}{a^{2}(t)}\partial_{i}^{2}\phi+3H\dot{\phi}+\frac{\partial
V}{\partial\phi}+\frac{\alpha}{4f}\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}=0\,,$
(99)

or in conformal time,

$\displaystyle\partial_{\tau}^{2}\phi-\partial_{i}^{2}\phi+2{\cal
H}\partial_{\tau}\phi+a^{2}(\tau)\frac{\partial
V}{\partial\phi}+a^{2}(\tau)\frac{\alpha}{4f}\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}=0\,.$
(100)

Next, we turn to the equations of motion for the vector field which we derive
in conformal time. From here on every time we write $a$ for the scale factor
we mean implicitly $a(\tau)$. After some algebra Eq.(98) can be cast into the
form,

$\displaystyle
g^{\alpha\nu}g^{\beta\sigma}\partial_{\alpha}F_{\nu\sigma}+\frac{\alpha}{2f}(\partial_{\alpha}\phi)\epsilon^{\alpha\beta\rho\sigma}F_{\rho\sigma}-m^{2}g^{\beta\mu}A_{\mu}=0\,.$
(101)

This gives four equations of motion for the massive vector field where the
$\beta=0$ and $\beta=l$ components can be written respectively as,

$\displaystyle(\partial_{i}^{2}A_{0}-\partial_{0}\partial_{i}A_{i})+\frac{\alpha}{f}\tilde{\epsilon}^{ijk}(\partial_{i}\phi)\partial_{j}A_{k}-a^{2}m^{2}A_{0}$
$\displaystyle=$ $\displaystyle 0\,,$ (102) $\displaystyle
g^{00}g^{li}\partial_{0}F_{0i}+g^{kj}g^{li}\partial_{k}F_{ji}+\frac{\alpha}{2f}\left[(\partial_{0}\phi)\epsilon^{0ljk}F_{jk}+(\partial_{i}\phi)\epsilon^{ilk0}F_{k0}\right]-m^{2}g^{li}A_{i}$
$\displaystyle=$ $\displaystyle 0\,.$ (103)

After some algebra the $\beta=l$ component can be massaged into the form,

$\displaystyle\partial_{0}\partial_{0}A_{l}$ $\displaystyle-$
$\displaystyle\partial_{i}\partial_{i}A_{l}-\partial_{l}\partial_{0}A_{0}+\partial_{l}\partial_{i}A_{i}$
(104) $\displaystyle-$
$\displaystyle\frac{\alpha}{f}\tilde{\epsilon}^{ljk}\left[(\partial_{0}\phi)\partial_{j}A_{k}+\frac{1}{2}(\partial_{j}\phi)(\partial_{k}A_{0}-\partial_{0}A_{k})\right]+a^{2}m^{2}A_{l}=0\,.$

As a massive vector has three degrees of freedom, we need to apply one
constraint which can be obtained by acting on Eq.(101) with
$\partial_{\beta}$,

$\displaystyle
g_{\beta\gamma}\partial^{\gamma}\left[g^{\alpha\nu}g^{\beta\sigma}\partial_{\alpha}F_{\nu\sigma}+\frac{\alpha}{2f}(\partial_{\alpha}\phi)\epsilon^{\alpha\beta\rho\sigma}F_{\rho\sigma}-m^{2}g^{\beta\mu}A_{\mu}\right]=0\,,$
(105)

which after some algebra leads to,

$\frac{\alpha}{f}\tilde{\epsilon}^{ijk}(\partial_{i}\phi)\partial_{j}A_{k}=\frac{a^{2}m^{2}}{4{\cal
H}}(2{\cal H}A_{0}+\partial^{0}A_{0}-\partial^{i}A_{i})\,.$ (106)

We now proceed to write more explicitly these equations, first projecting
along the longitudinal mode and then along the transverse modes.

We first consider the longitudinal mode for which we have,

$\displaystyle\tilde{\epsilon}^{ijk}(\partial_{i}\phi)\partial_{j}A_{k}\rightarrow(\nabla\phi(\vec{x},\tau))\cdot(i\vec{k}\times\vec{A}(\vec{k},\tau))=0\,.$
(107)

This then gives for Eq.(102),

$\displaystyle(\partial_{i}^{2}A_{0}-\partial_{0}\partial_{i}A_{i})-a^{2}m^{2}A_{0}=0\,,$
(108)

which in Fourier space has the solution,

$\displaystyle
A_{0}(\vec{k},\tau)=\frac{-i\vec{k}\cdot\partial_{\tau}\vec{A}(\vec{k},\tau)}{k^{2}+a^{2}m^{2}}\,.$
(109)

We next consider the terms in squared parentheses in Eq.(104),

$\displaystyle\tilde{\epsilon}^{ljk}\left[(\partial_{0}\phi)\partial_{j}A_{k}+\frac{1}{2}(\partial_{j}\phi)(\partial_{k}A_{0}-\partial_{0}A_{k})\right]$
$\displaystyle\rightarrow(\partial_{\tau}\phi)\left(i\vec{k}\times\vec{A}_{L}(\vec{k},\tau)\right)+\frac{1}{2}(\nabla\phi)\times\left(i\vec{k}A_{0}(\vec{k},\tau)-\partial_{\tau}\vec{A}_{L}(\vec{k},\tau)\right)\,,$
(110)

This results in vectors which are orthogonal to the direction of $\vec{k}$ and
so do not contribute to the equation of the longitudinal mode. In Fourier
space Eq.(104) then simplifies to,

$\displaystyle\partial_{\tau}^{2}A_{L}(\vec{k},\tau)-ik\partial_{\tau}A_{0}(\vec{k},\tau)+a^{2}m^{2}A_{L}(\vec{k},\tau)=0\,.$
(111)

Using the solution for $A_{0}$ in Eq.(109) this becomes,

$\displaystyle\partial_{\tau}^{2}A_{L}(\vec{k},\tau)+\frac{2k^{2}}{k^{2}+a^{2}m^{2}}{\cal
H}\partial_{\tau}A_{L}(\vec{k},\tau)+(k^{2}+a^{2}m^{2})A_{L}(\vec{k},\tau)=0\,.$
(112)

Switching to cosmic time the equation of motion for the longitudinal mode
becomes,

$\displaystyle\ddot{A}_{L}(\vec{k},t)+\frac{3k^{2}+a^{2}(t)m^{2}}{k^{2}+a^{2}(t)m^{2}}H\dot{A}_{L}(\vec{k},t)+\left(\frac{k^{2}}{a^{2}(t)}+m^{2}\right)A_{L}(\vec{k},t)=0\,,$
(113)

which is in agreement with Graham:2015rva .

Turning to the transverse modes for which $\vec{k}\cdot\vec{A}=0$, we apply
the constraint of Eq.(106) to Eq.(102) which leads to,

$\displaystyle\left(\partial_{i}^{2}-\frac{a^{2}m^{2}}{2}+\frac{a^{2}m^{2}}{4{\cal
H}}\partial^{0}\right)A_{0}=0\,.$ (114)

Integrated from some initial conformal time $\tau_{0}$, the solution to this
differential equation in Fourier space is given by,

$\displaystyle
A_{0}(\vec{k},\tau)=\frac{a^{2}(\tau)}{a^{2}(\tau_{0})}A_{0}(\tau_{0})e^{-2\left(\frac{k^{2}}{a^{2}(\tau)m^{2}}-\frac{k^{2}}{a^{2}(\tau_{0})m^{2}}\right)}\,.$
(115)

Here we are not free to choose the initial condition $A_{0}(\tau_{0})$ and the
only possible choice is the one which makes the solution for $A_{0}$
consistent with Eq.(109), that is $A_{0}(\tau_{0})=0$. Note that the
constraint in Eq.(106) forces,

$\tilde{\epsilon}^{ijk}(\partial_{i}\phi)\partial_{j}A_{k}=0\,,$ (116)

consistently for both longitudinal and transverse modes. Thus, we conclude
that,

$\displaystyle A_{0}(\vec{k},\tau)$
$\displaystyle=\frac{-ik\cdot\partial_{\tau}A_{L}(\vec{k},\tau)}{k^{2}+a^{2}m^{2}}\,,$
$\displaystyle{\rm longitudinal}\,,$ (117) $\displaystyle A_{0}(\vec{k},\tau)$
$\displaystyle=0\,,$ $\displaystyle{\rm transverse}\,.$ (118)

Note in particular these solutions imply that the time component of the vector
field does _not_ mix the transverse and longitudinal components. With these
solutions Eq.(104) for the transverse modes in Fourier space becomes,

$\displaystyle(\partial_{\tau}^{2}+k^{2}+a^{2}m^{2})\vec{A}_{T}(\vec{k},\tau)-\frac{\alpha}{f}\left[(\partial_{\tau}\phi(\vec{x},\tau))(i\vec{k}\times\vec{A}_{T}(\vec{k},\tau))-\frac{1}{2}(\nabla\phi(\vec{x},\tau))\times(\partial_{\tau}\vec{A}_{T}(\vec{k},\tau))\right]=0\,.$

Rewriting this in cosmic time and switching back to coordinate space we
have101010Note this differs from Eq. (4) in Agrawal:2018vin  where we have a
factor of two difference in our last term and we do not have a term
proportional to $\nabla\phi\times\nabla A_{0}$.,

$\displaystyle\ddot{\vec{A}}_{T}+H\dot{\vec{A}}_{T}-\frac{\nabla^{2}}{a^{2}}\vec{A}_{T}+m^{2}\vec{A}_{T}-\frac{1}{a}\frac{\alpha}{f}\left[\dot{\phi}(\nabla\times\vec{A}_{T})-\frac{1}{2}\nabla\phi\times\dot{\vec{A}}_{T}\right]=0\,.$
(120)

Next we consider the case of a homogeneous scalar field where we approximate
$\nabla\phi\approx 0$ so that $\phi(\vec{x},t)\approx\phi(t)$, as appropriate
for $\phi$ during inflation. Using Eq.(83) for the transverse modes Eq.(IV.1)
then simplifies to,

$\displaystyle\partial_{\tau}^{2}A_{\pm}(k,\tau)+\left[k^{2}\mp\frac{\alpha}{f}(\partial_{\tau}\phi)k+a^{2}m^{2}\right]A_{\pm}(k,\tau)=0\,.$
(121)

Using the fact that the conformal time during inflation is
$\tau\approx-(aH)^{-1}$ we can rewrite,

$\displaystyle\partial_{\tau}\phi=a\frac{\partial\phi}{\partial
t}=aH\frac{\dot{\phi}}{H}\approx-\frac{1}{\tau}\frac{\dot{\phi}}{H}\,.$ (122)

The equation of motion then reads in conformal time,

$\displaystyle\partial_{\tau}^{2}A_{\pm}(k,\tau)+\left[k^{2}\pm\frac{\alpha}{f}\frac{\dot{\phi}}{H}\frac{k}{\tau}+a^{2}m^{2}\right]A_{\pm}(k,\tau)=0\,.$
(123)

Using $\partial_{\tau}^{2}A=a^{2}(t)H\dot{A}+a^{2}(t)\ddot{A}$ this can be
rewritten in cosmic time as,

$\displaystyle\ddot{A}_{\pm}(k,t)+H\dot{A}_{\pm}(k,t)+\left[\frac{k^{2}}{a^{2}(t)}\mp\frac{\alpha}{f}\dot{\phi}\frac{k}{a(t)}+m^{2}\right]A_{\pm}(k,t)=0\,.$
(124)

_Energy and pressure densities_

Starting from the action, we can compute the stress-energy tensor via,

$\displaystyle T_{\alpha\beta}=-\frac{2}{\sqrt{-g}}\frac{\delta S}{\delta
g^{\alpha\beta}}\,,$ (125)

where the following relations will be useful for the calculation,

$\displaystyle\frac{\delta}{\delta
g^{\alpha\beta}}\sqrt{-g}=-\frac{1}{2}\sqrt{-g}\
g_{\alpha\beta}\,,\qquad\frac{\delta}{\delta g^{\alpha\beta}}(-g)=g\
g_{\alpha\beta}\,.$ (126)

One can verify that after computing $\frac{\delta S}{\delta g^{\alpha\beta}}$,
there are two terms proportional to $\phi F\tilde{F}$ which cancel exactly.
Thus, we see that the operator responsible for inducing the tachyonic
instability does not contribute to the energy density. Explicitly we find,

$\displaystyle T_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\partial_{\alpha}\phi\partial_{\beta}\phi+g^{\mu\nu}F_{\mu\alpha}F_{\nu\beta}+m^{2}A_{\alpha}A_{\beta}$
(127) $\displaystyle-$ $\displaystyle
g_{\alpha\beta}\left(\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi+V(\phi)+\frac{1}{4}g^{\mu\nu}g^{\rho\sigma}F_{\mu\rho}F_{\nu\sigma}+\frac{1}{2}m^{2}g^{\mu\nu}A_{\mu}A_{\nu}\right)\,.$

In cosmic time $t$ the energy density is then given by,

$\displaystyle\rho=T_{00}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\dot{\phi}^{2}+\frac{1}{2a^{2}}(\partial_{i}\phi)^{2}+V(\phi)$
$\displaystyle+$
$\displaystyle\frac{1}{2a^{2}}(\partial_{t}A_{i}-\partial_{i}A_{0})^{2}+\frac{1}{4a^{4}}(\partial_{i}A_{j}-\partial_{j}A_{i})^{2}+\frac{1}{2}m^{2}A_{0}^{2}+\frac{1}{2a^{2}}m^{2}A_{i}^{2}$
$\displaystyle=$ $\displaystyle\rho_{\phi}+\rho_{A}\,,$

where $\rho_{\phi}$ and $\rho_{A}$ denote the contributions from the inflaton
and dark vector respectively. The pressure density can then be computed as,

$\displaystyle p$ $\displaystyle=$
$\displaystyle\frac{1}{3}(g^{\alpha\beta}T_{\alpha\beta}+\rho)=\frac{1}{2}\dot{\phi}^{2}-\frac{1}{6a^{2}}(\partial_{i}\phi)^{2}-V(\phi)$
$\displaystyle+$
$\displaystyle\frac{1}{6a^{2}}(\partial_{t}A_{i}-\partial_{i}A_{0})^{2}+\frac{1}{12a^{4}}(\partial_{i}A_{j}-\partial_{j}A_{i})^{2}+\frac{1}{2}m^{2}A_{0}^{2}-\frac{1}{6a^{2}}m^{2}A_{i}^{2}$
$\displaystyle=$ $\displaystyle p_{\phi}+p_{A}\,,$

where $p_{\phi}$ and $p_{A}$ denote the contributions from the inflaton and
dark vector respectively.

In Fourier space, we can write the energy density of the vector field
(Eq.(IV.1)) as,

$\displaystyle\rho_{A}(\vec{x},t)$ $\displaystyle=$
$\displaystyle\frac{1}{2a^{2}}\int\frac{d^{3}k\
d^{3}k^{\prime}}{(2\pi)^{6}}e^{i(\vec{k}+\vec{k}^{\prime})\cdot\vec{x}}\Big{\\{}[\partial_{t}\vec{A}(\vec{k},t)\cdot\partial_{t}\vec{A}(\vec{k}^{\prime},t)]+m^{2}[\vec{A}(\vec{k},t)\cdot\vec{A}(\vec{k}^{\prime},t)]$
$\displaystyle-$
$\displaystyle\frac{1}{a^{2}}\left([\vec{k}\cdot\vec{k}^{\prime}][\vec{A}(\vec{k},t)\cdot\vec{A}(\vec{k}^{\prime},t)]-[\vec{k}\cdot\vec{A}(\vec{k}^{\prime},t)][\vec{k}^{\prime}\cdot\vec{A}(\vec{k},t)]\right)$
$\displaystyle-$ $\displaystyle
iA_{0}(\vec{k},t)[\vec{k}\cdot\partial_{t}\vec{A}(\vec{k}^{\prime},t)]-[\vec{k}\cdot\vec{k}^{\prime}][A_{0}(\vec{k},t)A_{0}(\vec{k}^{\prime},t)]$
$\displaystyle-$ $\displaystyle
iA_{0}(\vec{k}^{\prime},t)[\vec{k}^{\prime}\cdot\partial_{t}\vec{A}(\vec{k},t)]+a^{2}m^{2}A_{0}(\vec{k},t)A_{0}(\vec{k}^{\prime},t)\Big{\\}}\,.$

We can “bracket” the expression above and use the definition of the power
spectrum Eq.(87) to obtain the spatial average $\langle\rho_{A}(t)\rangle$.
Separating the transverse and longitudinal modes, we have $A_{0}=0$ for the
former so last two lines in Eq.(IV.1) vanish. After some algebra we have,

$\displaystyle\langle\rho_{A}^{T}(t)\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2a^{2}}\int_{0}^{\infty}\frac{dk}{k}\left[{\cal
P}_{\partial_{t}A_{T}}(k,t)+\left(\frac{k^{2}}{a^{2}}+m^{2}\right){\cal
P}_{A_{T}}(k,t)\right],\,$ $\displaystyle\langle\rho_{A}^{T}(\tau)\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{2a^{4}}\int_{0}^{\infty}\frac{dk}{k}\left[{\cal
P}_{\partial_{\tau}A_{T}}(k,\tau)+\left(k^{2}+a^{2}m^{2}\right){\cal
P}_{A_{T}}(k,\tau)\right]\,.$ (131)

in cosmic and conformal time respectively.

For the longitudinal mode, we use Eq.(109) and its analog in cosmic time,

$\displaystyle
A_{0}(\vec{k},t)=\frac{-i\vec{k}\cdot\partial_{t}\vec{A}(\vec{k},t)}{k^{2}+a^{2}m^{2}}\,,$
(132)

and substitute them into Eq.(IV.1). After some algebra this gives,

$\displaystyle\langle\rho_{A}^{L}(t)\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2a^{2}}\int_{0}^{\infty}\frac{dk}{k}\left[\frac{a^{2}m^{2}}{k^{2}+a^{2}m^{2}}{\cal
P}_{\partial_{t}A_{L}}(k,t)+m^{2}{\cal P}_{A_{L}}(k,t)\right]\,,$
$\displaystyle\langle\rho_{A}^{L}(\tau)\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2a^{4}}\int_{0}^{\infty}\frac{dk}{k}\left[\frac{a^{2}m^{2}}{k^{2}+a^{2}m^{2}}{\cal
P}_{\partial_{\tau}A_{L}}(k,\tau)+a^{2}m^{2}{\cal
P}_{A_{L}}(k,\tau)\right]\,,$ (133)

in agreement with Graham:2015rva . Note that the expressions in Eq.(IV.1) and
Eq.(IV.1) are valid both for the classical and the quantum gauge field upon
using the corresponding definitions of the power spectra discussed in the
previous section.

###  IV.2 Analytic study of the energy density spectrum during inflation

In this section we utilize the (approximate) analytic solutions of the
equations of motion to examine the energy density spectrum of the dark vector
at the end of inflation. Starting from Eq.(123) and using
$\tau\approx-(aH)^{-1}$ during inflation we define,

$\displaystyle\xi\equiv\frac{\alpha\dot{\phi}}{2Hf}\,,\qquad\bar{m}\equiv\frac{m}{H}\,,$
(134)

which allows us to write the equation of motion for the transverse modes as,

$\displaystyle\partial_{\tau}^{2}A_{\pm}(k,\tau)+\left[k^{2}\pm
2\xi\frac{k}{\tau}+\frac{\bar{m}^{2}}{\tau^{2}}\right]A_{\pm}(k,\tau)=0\,.$
(135)

Introducing the dimensionless variable,

$\displaystyle x=-k\tau\approx\frac{k}{aH}\,,$ (136)

the equation of motion then becomes,

$\displaystyle\partial_{x}^{2}A_{\pm}(x)+\left[1\mp\frac{2\xi}{x}+\frac{\bar{m}^{2}}{x^{2}}\right]A_{\pm}(x)=0\,.$
(137)

Neglecting the time dependence in $\xi$ and the Hubble parameter, this
equation of motion can be solved analytically. Noting that $x>0$ and using the
convention $\xi>0$, the mode that gets exponentially enhanced is $A_{+}$ and
we can neglect $A_{-}$ in what follows. The solution to Eq.(137), once
properly normalized, is given in terms of the Whittaker function
Meerburg:2012id ,

$\displaystyle
A_{+}(x)=\frac{e^{\pi\xi/2}}{\sqrt{2k}}W_{-i\xi,\mu}(-2ix)\,,\qquad\mu=\sqrt{1/4-\bar{m}^{2}}\,.$
(138)

From this we can obtain the power spectra as,

$\displaystyle{\cal
P}_{\partial_{\tau}A_{+}}=\frac{k^{3}}{2\pi^{2}}|\partial_{\tau}A_{+}|^{2}\,,\qquad{\cal
P}_{A_{+}}=\frac{k^{3}}{2\pi^{2}}|A_{+}|^{2}\,.$ (139)

Starting with the spatially averaged energy density as a function of $\tau$ in
Eq.(IV.1) and the definition of $x$ in Eq.(136), we can write the energy
density as,

$\displaystyle\langle\rho_{A}^{T}\rangle$
$\displaystyle=\frac{H^{4}}{8\pi^{2}}\int\frac{dx}{x}\ 2k\
x^{4}\left[|\partial_{x}A_{+}(x)|^{2}+\left(1+\frac{\bar{m}^{2}}{x^{2}}\right)|A_{+}(x)|^{2}\right]\,,$
(140) $\displaystyle\frac{d\langle\rho_{A}^{T}\rangle}{d\ln x}$
$\displaystyle=\frac{H^{4}}{8\pi^{2}}2k\
x^{4}\left[|\partial_{x}A_{+}(x)|^{2}+\left(1+\frac{\bar{m}^{2}}{x^{2}}\right)|A_{+}(x)|^{2}\right]$
(141)
$\displaystyle=\frac{H^{4}}{8\pi^{2}}e^{\pi\xi}x^{4}\left[|\partial_{x}W_{-i\xi,\mu}(-2ix)|^{2}+\left(1+\frac{\bar{m}^{2}}{x^{2}}\right)|W_{-i\xi,\mu}(-2ix)|^{2}\right]\,.$
(142)

Note the above is only a function of $x$ under the assumptions that
$\tau\approx-(aH)^{-1}$ and $\xi$ is constant, both of which are good
approximations during slow roll inflation.

We next separate the contributions to the energy density spectrum into the
electric and magnetic components (and drop the brackets for $\rho_{A}$),

$\displaystyle\frac{d\rho_{E}}{d\ln x}$ $\displaystyle=$
$\displaystyle\frac{H^{4}}{8\pi^{2}}e^{\pi\xi}x^{4}|\partial_{x}W_{-i\xi,\mu}(-2ix)|^{2}\,,$
$\displaystyle\frac{d\rho_{B}}{d\ln x}$ $\displaystyle=$
$\displaystyle\frac{H^{4}}{8\pi^{2}}e^{\pi\xi}(x^{4}+x^{2}\bar{m}^{2})|W_{-i\xi,\mu}(-2ix)|^{2}\,.$
(143)

In Fig. 11 we plot $\frac{d\rho_{E}}{d\ln x}$ and $\frac{d\rho_{B}}{d\ln x}$
as a function of $x$ for the listed choice of parameters $\bar{m}$, $\xi$, and
$H$ during inflation, $H_{I}$. Both components show a peak at $x\sim 0.1-1$
which indicates the point of maximal tachyonic enhancement. Moving to lower
values of $x$, the magnetic component drops as $x^{4}$ until it reaches
$x=\bar{m}$ and then decreases as $x^{2}$. This behavior can be understood
directly from Eq.(IV.2) where for $x<\bar{m}$ the term $x^{2}\bar{m}^{2}$
dominates over $x^{4}$. This can be traced back easily to $a^{2}m^{2}$
dominating over $k^{2}$ in the last two terms of Eq.(IV.1). The electric
component, dominant for $\bar{m}<x<0.1$, continues decreasing as $x^{4}$ until
values of $x$ smaller than $\bar{m}^{2}/(2\xi)$ after which point it goes
through a kink and then decreases as $x^{2}$. It is also clear from the plot
that at small $x$ the energy density is dominated by the mass term which we
have included in the magnetic component.

Figure 11: Electric (solid) and magnetic (dashed) components of the energy
density spectra normalized to the energy density of the inflaton
$\rho_{I}=3M_{Pl}^{2}H_{I}^{2}$.

Looking more closely at the electric component to understand the kink and the
change of slope, the location of the kink occurs where $\partial_{x}A_{+}(x)$
changes sign and the electric energy density goes to zero. It is located at
$x<\bar{m}^{2}/(2\xi)$, but a compact analytic formula for the exact $x$ is
difficult to obtain due to the complicated form of the solution in terms of
the Whittaker function. We can however easily understand the shape to the
right and left of the kink by inspection of the equation of motion. For
$\bar{m}^{2}/(2\xi)\ll x<2\xi$, we can neglect the mass term in Eq.(137). Then
the solution for the tachyonic mode becomes,

$\displaystyle
A_{+}(x)=\frac{1}{2k}\sqrt{\frac{4x}{\pi}}e^{\pi\xi}K_{1}(2\sqrt{2x\xi})\,,$
(144)

where $K_{1}$ is the Bessel function. We plug this it into Eq.(141) and expand
for $x\ll 1$ to find,

$\displaystyle\frac{d\rho_{E}}{d\ln x}=\frac{H^{4}\xi
e^{2\pi\xi}}{4\pi^{3}}\left(2\gamma+\ln(2x\xi)\right)^{2}x^{4}+{\cal
O}(x^{9/2})\,,\quad\bar{m}^{2}/(2\xi)<x\ll 1\,,$ (145)

where $\gamma=0.577$ is the Euler’s constant and we see in Fig. 12 this is an
excellent approximation to the full solution in the regime
$\bar{m}^{2}/(2\xi)\ll x<2\xi$. We also see from the approximate solution in
Eq.(145) that the slope to the right of the kink is $\sim x^{4}$.

Figure 12: Electric component of the energy density spectrum. For
$x>\bar{m}^{2}/(2\xi)$ the full solution (solid red) to the equation of motion
is well approximated by the massless one (dashed gray) given in Eq.(144). For
$x\ll\bar{m}^{2}/(2\xi)$ the mass term dominates and the full solution is
approximated by Eq.(147) (dashed black). The slope changes from $x^{4}$ on the
right of the kink to $x^{2}$ on the left.

For $x\ll\bar{m}^{2}/(2\xi)$, the mass term dominates in the equation of
motion,

$\displaystyle\partial_{x}^{2}A_{+}(x)+\frac{\bar{m}^{2}}{x^{2}}A_{+}(x)=0\,,$
(146)

which then leads to the approximate solution,

$\displaystyle A_{+}(x)=c_{1}x^{\bar{m}^{2}}+c_{2}x^{1-\bar{m}^{2}}\,,\qquad
x\ll\bar{m}^{2}/(2\xi)\,.$ (147)

Fixing $c_{1}$ and $c_{2}$ to match the normalization of the full solution and
plugging into Eq.(141),

$\displaystyle\frac{d\rho_{E}}{d\ln
x}=\frac{H^{4}}{8\pi^{2}}c_{1}^{2}\bar{m}^{4}x^{2}+{\cal O}(x^{3})\,,\qquad
x\ll\bar{m}^{2}/(2\xi)\,,\quad\bar{m}\ll 1\,.$ (148)

The solution in this regime is shown in dashed black in Fig. 12 where we see
it is an excellent approximation to the full solution for
$x\ll\bar{m}^{2}/(2\xi)$ and the slope goes like $x^{2}$.

###  IV.3 Power spectrum of energy density fluctuations

Here we derive the power spectrum of the energy density fluctuations. At late
times once the dark vector is non-relativistc the energy density is well
approximate by the mass term in the lagrangian, $\rho\sim m^{2}A^{2}$. We can
then define the energy density contrast $\delta(\vec{x})$ which measures the
deviation from the mean energy density,

$\displaystyle\rho(\vec{x})=\langle\rho\rangle(1+\delta(\vec{x}))=\frac{1}{2}m^{2}\left(\langle\hat{\vec{A}}(\vec{x})^{2}\rangle+\hat{\vec{A}}(\vec{x})^{2}\right)\,,$
(149)

where we have implicitly dropped the cosmic time variable $t$. From this we
identify,

$\displaystyle\langle\rho\rangle$
$\displaystyle=\frac{1}{2}m^{2}\langle\vec{A}(\vec{x})^{2}\rangle,$ (150)
$\displaystyle\delta(\vec{x})$
$\displaystyle=\frac{\vec{A}(\vec{x})^{2}}{\langle\vec{A}(\vec{x})^{2}\rangle}=\frac{1}{\langle\vec{A}(\vec{x})^{2}\rangle}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{d^{3}q}{(2\pi)^{3}}e^{i(\vec{p}+\vec{q})\cdot
x}\vec{A}(\vec{p})\cdot\vec{A}(\vec{q})$
$\displaystyle=\frac{1}{\langle\vec{A}(\vec{x})^{2}\rangle}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}}\int\frac{d^{3}q}{(2\pi)^{3}}\vec{A}(\vec{k}-\vec{q})\cdot\vec{A}(\vec{q})$
(151)
$\displaystyle\equiv\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}}\delta(\vec{k})\,,$

where we have defined the momentum shift $\vec{k}=\vec{p}+\vec{q}$. From this
we can then read off the Fourier transform of the density contrast,

$\displaystyle\delta(\vec{k})=\frac{1}{\langle\vec{A}(\vec{x})^{2}\rangle}\int\frac{d^{3}q}{(2\pi)^{3}}\vec{A}(\vec{k}-\vec{q})\cdot\vec{A}(\vec{q})\,.$
(152)

Even though the exponentially enhanced tachyonic modes are highly classical,
here we can work explicitly with the quantum field defined in Eq.(89). To
compute the $\hat{\vec{A}}(\vec{x})^{2}$ we have,

$\displaystyle\hat{\vec{A}}(\vec{x})^{2}$ $\displaystyle=$
$\displaystyle\sum_{\lambda,\lambda^{\prime}}\int\frac{d^{3}pd^{3}q}{(2\pi)^{6}}e^{i(\vec{p}+\vec{q})\cdot
x}\vec{\epsilon}_{\lambda}(\vec{p})\cdot\vec{\epsilon}_{\lambda^{\prime}}(\vec{q})$
$\displaystyle\times$
$\displaystyle\left(A_{\lambda}(p)a_{\lambda}(\vec{p})+A_{\lambda}(p)^{*}a^{\dagger}_{\lambda}(-\vec{p})\right)\left(A_{\lambda^{\prime}}(q)a_{\lambda^{\prime}}(\vec{q})+A_{\lambda^{\prime}}(q)^{*}a^{\dagger}_{\lambda^{\prime}}(-\vec{q})\right)$
$\displaystyle=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}}\sum_{\lambda,\lambda^{\prime}}\int\frac{d^{3}q}{(2\pi)^{3}}\vec{\epsilon}_{\lambda}(\vec{k}-\vec{q})\cdot\vec{\epsilon}_{\lambda^{\prime}}(\vec{q})$
$\displaystyle\times$
$\displaystyle\left(A_{\lambda}(k-q)a_{\lambda}(\vec{k}-\vec{q})+A_{\lambda}(k-q)^{*}a^{\dagger}_{\lambda}(\vec{q}-\vec{k})\right)\left(A_{\lambda^{\prime}}(q)a_{\lambda^{\prime}}(\vec{q})+A_{\lambda^{\prime}}(q)^{*}a^{\dagger}_{\lambda^{\prime}}(-\vec{q})\right)\,.$

where again we have defined the momentum shift $\vec{k}=\vec{p}+\vec{q}$. From
Eq.(152) we then have,

$\displaystyle\delta(\vec{k})$ $\displaystyle=$
$\displaystyle\frac{1}{\langle\hat{\vec{A}}(\vec{x})^{2}\rangle}\sum_{\lambda,\lambda^{\prime}}\int\frac{d^{3}q}{(2\pi)^{3}}\vec{\epsilon}_{\lambda}(\vec{k}-\vec{q})\cdot\vec{\epsilon}_{\lambda^{\prime}}(\vec{q})$
$\displaystyle\times$
$\displaystyle\left(A_{\lambda}(k-q)a_{\lambda}(\vec{k}-\vec{q})+A_{\lambda}(k-q)^{*}a^{\dagger}_{\lambda}(\vec{q}-\vec{k})\right)\left(A_{\lambda^{\prime}}(q)a_{\lambda^{\prime}}(\vec{q})+A_{\lambda^{\prime}}(q)^{*}a^{\dagger}_{\lambda^{\prime}}(-\vec{q})\right)\,.$

The space-dependent part of the energy density contrast is obtained from the
two-point correlation function,

$\displaystyle\langle\delta(\vec{k})\delta(\vec{k}^{\prime})\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{\langle\hat{\vec{A}}(\vec{x})^{2}\rangle^{2}}\sum_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}\int\frac{d^{3}q\
d^{3}q^{\prime}}{(2\pi)^{6}}\vec{\epsilon}_{\lambda_{1}}(\vec{k}-\vec{q})\cdot\vec{\epsilon}_{\lambda_{2}}(\vec{q})\
\vec{\epsilon}_{\lambda_{3}}(\vec{k}^{\prime}-\vec{q}^{\
\prime})\cdot\vec{\epsilon}_{\lambda_{4}}(\vec{q}^{\ \prime})$
$\displaystyle\times$ $\displaystyle\langle
0|\left(A_{\lambda_{1}}(k-q)a_{\lambda_{1}}(\vec{k}-\vec{q})\right)\left(A_{\lambda_{2}}(q)a_{\lambda_{2}}(\vec{q})+A_{\lambda_{2}}(q)^{*}a^{\dagger}_{\lambda_{2}}(-\vec{q})\right)$
$\displaystyle\times$
$\displaystyle\left(A_{\lambda_{3}}(k^{\prime}-q^{\prime})a_{\lambda_{3}}(\vec{k}^{\prime}-\vec{q}^{\
\prime})+A_{\lambda_{3}}(k^{\prime}-q^{\prime})^{*}a^{\dagger}_{\lambda_{3}}(\vec{q}^{\
\prime}-\vec{k}^{\prime})\right)\left(A_{\lambda_{4}}(q^{\prime})^{*}a^{\dagger}_{\lambda_{4}}(-\vec{q}^{\
\prime})\right)|0\rangle\,.$

Focusing on the last two lines, there is only one combination of creation and
annihilation operators that leads to a space-dependent result,

$\displaystyle\langle
0|a_{\lambda_{1}}a_{\lambda_{2}}a^{\dagger}_{\lambda_{3}}a^{\dagger}_{\lambda_{4}}|0\rangle$
$\displaystyle=$
$\displaystyle(2\pi)^{6}\Big{[}\delta_{\lambda_{2}\lambda_{3}}\delta_{\lambda_{1}\lambda_{4}}\delta^{3}(\vec{k}^{\prime}-\vec{q}^{\
\prime}+\vec{q})\delta^{3}(-\vec{q}^{\ \prime}-\vec{k}+\vec{q})$ (156)
$\displaystyle+$
$\displaystyle\delta_{\lambda_{1}\lambda_{3}}\delta_{\lambda_{2}\lambda_{4}}\delta^{3}(\vec{k}-\vec{q}-\vec{q}^{\
\prime}+\vec{k}^{\prime})\delta^{3}(\vec{q}+\vec{q}^{\ \prime})\Big{]}\,.$

Plugging this into Eq.(IV.3) and keeping only the tachyonic transverse mode
$A_{+}$ we find,

$\displaystyle\langle\delta(\vec{k})\delta(\vec{k}^{\prime})\rangle$
$\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{3}(\vec{k}+\vec{k}^{\prime})\frac{1}{\langle\hat{A}_{+}(\vec{x})^{2}\rangle^{2}}\int\frac{d^{3}q}{(2\pi)^{3}}2|A_{+}(k-q)|^{2}|A_{+}(q)|^{2}$
$\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{3}(\vec{k}+\vec{k}^{\prime})\frac{2}{\langle\hat{A}_{+}(\vec{x})^{2}\rangle^{2}}\frac{2\pi}{(2\pi)^{3}}\int
q^{2}dqd\cos\theta\frac{2\pi^{2}}{(k-q)^{3}}{\cal
P}_{A_{+}}(k-q)\frac{2\pi^{2}}{q^{3}}{\cal P}_{A_{+}}(q)$ $\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{3}(\vec{k}+\vec{k}^{\prime})\frac{2\pi^{2}}{\langle\hat{A}_{+}(\vec{x})^{2}\rangle^{2}}\int
dq\ dp\ q^{2}\frac{p}{kq}\frac{1}{p^{3}q^{3}}{\cal P}_{A_{+}}(p){\cal
P}_{A_{+}}(q)\,,$

where we have used Eq.(92) and the change of variables from $q,\cos\theta$ to
$q,p$ using,

$\displaystyle\vec{p}=\vec{k}-\vec{q}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ p^{2}=k^{2}+q^{2}-2kq\cos\theta\,,\leavevmode\nobreak\
\leavevmode\nobreak\ d\cos\theta=-\frac{p}{kq}dp\,,\leavevmode\nobreak\
\leavevmode\nobreak\ dqd\cos\theta\rightarrow\frac{p}{kq}dqdp\,,$ (158)

as well as trivially performed the $d\phi$ integral. Defining the power
spectrum for the energy density contrast in the quantum case (and for
transverse modes) as,

$\displaystyle\langle\delta(\vec{k})\delta(\vec{k}^{\prime})\rangle$
$\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{3}(\vec{k}+\vec{k}^{\prime})\frac{2\pi^{2}}{k^{3}}{\cal
P}_{\delta}(k)\,,$ (159)

we arrive finally at the expression given in Eq.(67),

$\displaystyle{\cal P}_{\delta}(k)$ $\displaystyle=$
$\displaystyle\frac{k^{2}}{\langle\hat{A}_{+}(\vec{x})^{2}\rangle^{2}}\int_{|q-k|<p<q+k}dq\
dp\ \frac{1}{q^{2}p^{2}}{\cal P}_{A_{+}}(p){\cal P}_{A_{+}}(q)\,,$ (160)
$\displaystyle\langle\hat{A}_{+}(\vec{x})^{2}\rangle^{2}=\left[\int_{0}^{\infty}\frac{dk}{k}{\cal
P}_{A_{+}}(k)\right]^{2}\,.$

This is the power spectrum of the energy density contrast for $\rho\sim
m^{2}A^{2}$. For completeness we have also computed the power spectrum for the
energy density contrast corresponding to the kinetic term in the energy
density for which we find the same expression as Eq.(160) with
$A_{+}\to\partial_{t}A_{+}$. For the longitudinal mode we find the same result
as in Graham:2015rva .

###  IV.4 Numerical procedure for solving equations of motion

Here we sketch the numerical solutions to the equations of motion and how the
input spectra are obtained. We have integrated the equations of motion for the
longitudinal and transverse vector perturbations at linear order in cosmic
time111111We have used the SUNDIALS package, “SUite of Nonlinear and
DIfferential/ALgebraic equation Solvers” for the numerical integration.,

$\displaystyle\ddot{A}_{L}(\vec{k},t)+\frac{3k^{2}+a^{2}(t)m^{2}}{k^{2}+a^{2}(t)m^{2}}H\dot{A}_{L}(\vec{k},t)+\left(\frac{k^{2}}{a^{2}(t)}+m^{2}\right)A_{L}(\vec{k},t)$
$\displaystyle=$ $\displaystyle 0\,,$ (161)
$\displaystyle\ddot{A}_{\pm}(k,t)+H\dot{A}_{\pm}(k,t)+\left[\frac{k^{2}}{a^{2}(t)}\mp\frac{\alpha}{f}\dot{\phi}\frac{k}{a(t)}+m^{2}\right]A_{\pm}(k,t)$
$\displaystyle=$ $\displaystyle 0\,,$ (162)

together with the background inflaton equation of motion,

$\ddot{\phi}+3H\dot{\phi}+V^{\prime}=0\,,$ (163)

where $V^{\prime}=dV(\phi)/d\phi$. Working in the regime of no dark vector
backreaction during single-field inflation, the Hubble parameter is given by,

$H^{2}=\frac{1}{3M_{PL}^{2}}\left(\frac{\dot{\phi}^{2}}{2}+V(\phi)\right)\,,$
(164)

and the number of e-folds is obtained by integrating,

$\frac{d\ln a}{dt}=H\,.$ (165)

We start the integration in the slow-roll regime with initial field velocity,

$\dot{\phi}=\frac{V^{\prime}}{3H}\,,$ (166)

and the initial value of the inflaton field for any given inflationary
potential $V(\phi)$ will set the total number of e-folds $N_{e}$ up to the end
of inflation at $\epsilon_{end}=\dot{\phi}^{2}/(2H^{2}M_{PL}^{2})=1$. For
example, for the models considered in Fig. 6, $\phi(0)=24M_{PL}$ gives
$N_{e}\simeq 72$ for a chaotic quartic potential, while $\phi(0)=17M_{PL}$
gives $N_{e}\simeq 73$ for a chaotic quadratic potential, $\phi(0)=M_{PL}$
gives $N_{e}\simeq 73$ for the hilltop quadratic potential with $v=6M_{PL}$,
and $\phi(0)=60M_{PL}$ gives $N_{e}\simeq 62$ for the axion-like model with
$\Lambda=24M_{PL}$.

For each $k$ mode, we start the integration with the vector fluctuations
initially in the Bunch-Davies vacuum defined as,

$\displaystyle A_{\lambda}^{(R)}(k,0)$
$\displaystyle=\frac{1}{\sqrt{2\omega_{k}}}\,,$ $\displaystyle
A_{\lambda}^{(I)}(k,0)$ $\displaystyle=0\,,$ (167)
$\displaystyle\dot{A}_{\lambda}^{(R)}(k,0)$ $\displaystyle=0\,,$
$\displaystyle\dot{A}_{\lambda}^{(I)}(k,0)$
$\displaystyle=-\frac{\omega_{k}}{\sqrt{2}}\,,$ (168)

where $\lambda$ refers to either longitudinal or transverse mode,
$\omega_{k}^{2}=k^{2}/a_{0}^{2}+m^{2}$, $a_{0}$ is the initial value of the
scale factor, and $R$($I$) refers to the real (imaginary) component of the
perturbation. Therefore, the initial vacuum power spectrum is given by,

$\displaystyle{\cal P}_{A_{\lambda}}(k,0)$ $\displaystyle=$
$\displaystyle\frac{k^{3}}{2\pi^{2}a_{0}^{3}}|A_{\lambda}(k,0)|^{2}=\frac{k^{3}}{4\pi^{2}\omega_{k}a_{0}^{3}}\,,$
(169) $\displaystyle{\cal P}_{\dot{A}_{\lambda}}(k,0)$ $\displaystyle=$
$\displaystyle\frac{k^{3}}{2\pi^{2}a_{0}^{3}}|\dot{A}_{\lambda}(k,0)|^{2}=\frac{k^{3}\omega_{k}}{4\pi^{2}a_{0}^{3}}\,.$
(170)

Two-point functions and energy densities for the dark vector perturbations
obtained directly from the momentum integration of the power spectrum are
clearly UV divergent and must be regularized. In an expanding universe this
can be accomplished using the adiabatic regularization method Parker:2009uva
which is based on a WKB-type expansion in powers of the time derivatives of
the scale factor and frequency modes. For our purposes, given that we are
interested in particle production effects driven by the tachyonic instability,
it is enough to regularize the expressions at zero-order, i.e. we only
subtract the vacuum contribution as in Minkowsky space,

$\displaystyle{\cal P}^{reg}_{A_{\lambda}}(k,t)$ $\displaystyle=$
$\displaystyle{\cal P}_{A_{\lambda}}(k,t)-{\cal P}_{A_{\lambda}}(k,0)\,,$
(171) $\displaystyle{\cal P}^{reg}_{\dot{A}_{\lambda}}(k,t)$ $\displaystyle=$
$\displaystyle{\cal P}_{\dot{A}_{\lambda}}(k,t)-{\cal
P}_{\dot{A}_{\lambda}}(k,0)\,.$ (172)

On the left hand side in Fig. 5 the dotted line shows the subtracted vacuum
contribution. We see all the spectra are above this line up to modes
$k/a_{end}\sim O(10)H_{\rm end}$ showing the particle production effect.
Higher momentum modes stay in the Bunch-Davies vacuum and the partial
subtraction we have performed leads to a spectra below the vacuum one. This
signals the absence of the tachyonic instability and particle production
effects for these modes.

## References