#  Dark photon dark matter from charged inflaton

Hassan Firouzjahi  Mohammad Ali Gorji  Shinji Mukohyama  Borna Salehian

###### Abstract

We present a scenario of vector dark matter production during inflation
containing a complex inflaton field which is charged under a dark gauge field
and which has a symmetry breaking potential. As the inflaton field rolls
towards the global minimum of the potential the dark photons become massive
with a mass which can be larger than the Hubble scale during inflation. The
accumulated energy of the quantum fluctuations of the produced dark photons
gives the observed relic density of the dark matter for a wide range of
parameters. Depending on the parameters, either the transverse modes or the
longitudinal mode or their combination can generate the observed dark matter
relic energy density.

YITP-20-139, IPMU-0115

##  1 Introduction

There are numerous compelling evidences for the existence of dark matter but
understanding its nature remains one of the important questions in cosmology.
Dark photon dark matter, alternatively known as “vector dark matter”, is one
of the candidate for dark matter which is widely studied in recent years [1,
2, 3, 4, 5, 6, 7]. It can be produced before the time of matter and radiation
equality through their direct interactions with other particles in the dark
sector or with the inflaton field and then playing the role of dark matter
after the time of matter and radiation equality. It should have negligible
direct interactions with the Standard Model (SM) particles in order not to
destroy the predictions of the SM [8, 9]. However, dark photon dark matter can
be indirectly probed by means of the gravitational waves observations through
its inevitable minimal coupling with gravity [10, 11, 12, 13, 14, 15, 16].

The vector dark matter models are usually constructed with a mechanism for the
production of dark photons after inflation but before the time of matter and
radiation equality [4, 17, 6, 7]. However, it is noticed that they can be
produced even during inflation from the longitudinal mode of a massive gauge
boson minimally coupled to the inflaton field [3, 18]. For the transverse
modes, while their energy density usually decays exponentially, such a decay
can be avoided and they can instead be produced if the gauge field is coupled
non-minimally to the inflaton field [19, 20]. In this way, the dark gauge
field acquires energy from the inflaton field and the transverse modes can be
produced efficiently during inflation to play the roles of dark matter after
the time of matter and radiation equality [21, 22, 23, 24]. For this purpose,
we need a large enough mass for dark photons so that they become non-
relativistic before the time of matter and radiation equality. On the other
hand, having a large mass compared to the Hubble scale causes an inefficient
particle production during inflation. That is why most of the previous studies
assumed a very small mass for dark photons. In Ref. [25] we have circumvented
these issues by generating a dynamical mass for the dark photon at the end of
inflation through the symmetry breaking i.e. the Higgs mechanism. In [25], a
real scalar field is the inflaton field while a massive complex field plays
the roles of symmetry breaking and terminating inflation, much like the
waterfall field of hybrid inflation. As the heavy complex field rolls rapidly
to the global minimum of its potential a mass is generated for the vector
field via the Higgs mechanism. The induced mass can be large such that the
excitations of the transverse modes can furnish the current dark matter relic
energy density. An important effect in [25] was that the longitudinal mode can
not contribute to the dark matter relic density as it is non-dynamical during
inflation while being too heavy to be produced efficiently at the end of
inflation.

In this work we consider a model of inflation in which the inflaton itself is
charged under the dark photon. The potential near the minimum is similar to a
simple chaotic model but with a non-zero vacuum expectation value (vev). In
order to amplify the dark photon excitations, as in [25], we allow a coupling
between the radial component of the inflaton field and the dark photon through
the gauge kinetic function. As the inflaton rolls towards its minimum with a
non-zero initial value, the dark photons acquire mass due to the Higgs
mechanism. Consequently, unlike [25], the longitudinal mode becomes dynamical
during inflation and can contribute to the dark matter relic density. In a
sense, this setup may be viewed as a hybrid of [25] where dark matter relic
comes from the energy density of transverse modes and [3] where only the
longitudinal mode is the origin of dark matter.

The rest of the paper is organized as follows. In section 2 we present our
setup and in section 3 we specify the form of the gauge kinetic function and
solve the mode functions. In section 4 we calculate the relic energy density
followed by a summary and discussions in section 5. Some technicalities
associated with the longitudinal mode are relegated to the Appendix.

##  2 Inflationary model with Higgs mechanism

We consider an inflationary model in which the inflaton is a complex scalar
field $\phi$ charged under a $U(1)$ gauge field $A_{\mu}$. We assume no
coupling between the gauge field and the SM particles, namely it is a dark
photon. The theory is described by the following action

$\displaystyle S=\int d^{4}x\sqrt{-g}\bigg{[}\frac{M_{\rm
Pl}^{2}}{2}R-\frac{1}{2}D_{\mu}\phi\overline{D^{\mu}\phi}-V(|\phi|)-\frac{1}{4}f^{2}(|\phi|)F_{\mu\nu}F^{\mu\nu}\bigg{]}\,,$
(2.1)

where $M_{\rm Pl}=1/\sqrt{8\pi G}$ is the reduced Planck mass and $G$ is the
Newton constant. The components of the covariant derivative acting on the
inflaton field are given by

$\displaystyle D_{\mu}\phi={\partial}_{\mu}\phi+i{\bf{e}}\,\phi\,A_{\mu}\,,$
(2.2)

in which ${\bf{e}}$ is the charge of the inflaton field which is a
dimensionless quantity. The components of the field strength tensor associated
to the gauge field is given by
$F_{\mu\nu}={\partial}_{\mu}A_{\nu}-{\partial}_{\nu}A_{\mu}$. Furthermore, we
have introduced the gauge kinetic function $f(|\phi|)$ so that the energy can
be transferred from the inflaton to the gauge field. This form of coupling is
extensively used in the context of anisotropic inflation [20, 26, 27, 28, 29,
30] and also primordial magnetogenesis [31, 32, 19, 33, 34, 35, 36, 37, 38,
39, 40, 41].

By assuming the $U(1)$ gauge symmetry, the potential and the gauge kinetic
function depends on only the modulus of the inflaton field $|\phi|$. We
decompose the inflaton into the radial and angular components as follows

$\displaystyle\phi\equiv\varphi e^{i\theta}\,,$ (2.3)

where $\theta$ determines the argument of the complex inflaton field $\phi$
while $\varphi$ characterizes the amplitude $\varphi=|\phi|$. By using the
local gauge redundancy of the theory we choose $\theta=0$, i.e. the unitary
gauge, and deal with a real and positive valued inflaton field $\varphi$. The
gauge kinetic function and the potential are then only a function of
$\varphi$. We will not assume any specific form for $f(\varphi)$ in this
section. However, we assume a power law form for it as a function of conformal
time in the next sections.

For the potential we do not assume any specific functional form other than a
couple of general features. First, the potential must satisfy the slow-roll
conditions to provide enough number of e-folds of accelerated expansion.
Second, we assume that the inflaton field settles finally into its global
minimum with $V=0$ with a a nonzero vev. This is to generate a mass for the
dark photon, as in Higgs mechanism, which is related to the inflaton field as
[25]

$\displaystyle m_{\tiny A}^{2}={\bf{e}}^{2}\varphi^{2}\,.$ (2.4)

Since the dark photon is meant to form the observed dark matter content of the
universe, it has to be massive after inflation. As a result, we need to have
$\varphi\neq 0$ at the end of inflation. For later usage, we define the ratio
of the gauge field mass to the Hubble scale during inflation as

$R\equiv\frac{m_{\tiny
A}^{2}}{H^{2}}=\frac{{\bf{e}}^{2}\varphi^{2}}{H^{2}}\,.$ (2.5)

Note that one can tune parameters of the theory such as ${\bf{e}}$ so that
dark photons can be heavier ($R>1$) or lighter ($R<1$) than the Hubble scale
during inflation. Further, the parameter $R$ generically depends on time. More
specifically we have

$\frac{\dot{R}}{HR}=\left(\sqrt{2\epsilon}+\frac{2M_{\rm
Pl}}{\varphi}\right)\sqrt{2\epsilon}\,,$ (2.6)

where $\epsilon\equiv-\dot{H}/H^{2}$ is the slow roll parameter and a dot
denotes the derivative with respect to the cosmic time. If we assume
$\varphi\sim\order{M_{\rm Pl}}$, which is in fact the case in many
inflationary models, then Eq. (2.6) shows that the fractional change of the
parameter $R$ in one e-fold is slow-roll suppressed. As we neglect slow-roll
corrections in the following analysis, we assume that $R$ does not change with
time. Including the effects of the running of the parameter $R$ will add
subleading corrections to our results.

One can easily build a potential with the above mentioned features. In the
following subsection we consider a simple example of such inflationary model.

###  2.1 Inflaton dynamics

As a simple explicit realization of the above setup, we consider the standard
symmetry breaking potential of the form

$V(\varphi)=\frac{\lambda}{4}\left(\varphi^{2}-v^{2}\right)^{2}\,,$ (2.7)

in which $\lambda$ is a dimensionless coupling and $v$ is the position of the
global minimum of the potential. The above potential is well motivated in the
context of symmetry breaking when a complex scalar field is charged under a
$U(1)$ field, such as in the Higgs symmetry breaking mechanism. The potential
is depicted in Fig. 1. We can imagine two types of inflationary scenarios for
$\varphi<v$ and $\varphi>v$. We will review these two types in the following.
Here, we do not consider them as realistic inflationary models as their
predictions are disfavoured by the Planck observation [42]. Instead, since the
computations in the next sections do not depend on the details of inflationary
models, we use these simple models just to show that generically $R$ can be
large.

Figure 1: A schematic view of the potential $V(\varphi)$. One can imagine two
inflationary scenarios depending on $\varphi<v$ or $\varphi>v$.

Hilltop case. First consider the case where inflation proceeds as in hilltop
model in the region $0<\varphi<v$. In this region, the background expansion is
given by

$\displaystyle 3M_{\rm Pl}^{2}H^{2}\simeq\frac{\lambda v^{4}}{4},$ (2.8)

and the number of e-folds needed to solve the flatness and the horizon
problems is given by

$\displaystyle\mathcal{N}=\frac{1}{M_{\rm Pl}^{2}}\int
d\varphi\,\frac{V}{V_{\varphi}}\simeq\frac{v^{2}}{4M_{\rm
Pl}^{2}}\ln{\frac{\varphi_{\rm e}}{\varphi_{\rm i}}}\,,$ (2.9)

where $\varphi_{\rm i}$ and $\varphi_{\rm e}$ are the initial and the final
values of the inflaton field respectively. We expect that $\varphi_{\rm e}$ is
not exponentially larger than $\varphi_{\rm i}$, i.e. $\varphi_{\rm
e}\sim\varphi_{\rm i}$ up to an $\order{1}$ factor, so one typically has
$\mathcal{N}\sim v^{2}/4M_{\rm Pl}^{2}$.

The slow-roll conditions require that $v>2M_{\rm Pl}$ and $\varphi<v/\sqrt{2}$
for a successful inflation model. One can compute the change in the parameter
$R$ during inflation as

$\frac{R_{\rm e}}{R_{\rm i}}=\left(\frac{H_{\rm i}}{H_{\rm
e}}\right)^{2}\left(\frac{\varphi_{\rm e}}{\varphi_{\rm
i}}\right)^{2}\simeq\left(\frac{\varphi_{\rm e}}{\varphi_{\rm
i}}\right)^{2}\,,$ (2.10)

which as discussed above is an $\order{1}$ value. Finally, the parameter $R$
at the end of inflation can be expressed as

$\displaystyle R_{\rm e}=\frac{{\bf{e}}^{2}\varphi_{\rm e}^{2}}{H_{\rm
e}^{2}}\simeq\frac{{\bf{e}}^{2}}{4\pi^{2}{\cal{P}}}\left(\frac{v}{2M_{\rm
Pl}}\right)^{4}\sim\frac{{\bf{e}}^{2}\mathcal{N}^{2}}{4\pi^{2}{\cal{P}}}\,,$
(2.11)

where we have used the definition of the dimensionless curvature perturbation
power spectrum as

$\displaystyle{\cal{P}}=\frac{H^{2}}{8\pi^{2}\epsilon M_{\rm Pl}^{2}}\simeq
2\times 10^{-9}\,.$ (2.12)

The relation (2.11) shows that with ${\bf{e}}\sim 1$, the parameter $R$ can be
as large as $10^{12}$. By reducing ${\bf{e}}$, smaller values of $R$ can be
obtained. For example, with ${\bf{e}}\sim 10^{-3}$ we obtain $R\sim 10^{6}$.

Large field. Now let us consider the case of large field model where
$\varphi>v$ so the field rolls from the large value of $\varphi$ toward its
global minimum at $\varphi=v$. For large values of $\varphi$, the potential
looks like $V\simeq\lambda\varphi^{4}/4$ so inflation proceeds as in chaotic
models with large field excursions. Of course, these types of large field
models are disfavoured by the Planck data as they produce large tensor
perturbations [42]. Here, however, we presented this setup as a simple working
example for the purpose of dark photon and symmetry breaking while a variant
of this potential can be invoked to be consistent with the Planck constraint
on primordial gravitational waves as well.

The number of e-fold to solve the flatness and the horizon problems in this
case is given by

$\displaystyle\mathcal{N}=\frac{1}{M_{\rm Pl}^{2}}\int
d\varphi\,\frac{V}{V_{\varphi}}\simeq\frac{1}{8M_{\rm Pl}^{2}}(\varphi_{\rm
i}^{2}-\varphi_{\rm e}^{2})\,.$ (2.13)

Note that $\varphi_{\rm i}>\varphi_{\rm e}>v$. Assuming $\varphi_{\rm
i}\gg\varphi_{\rm e}$ we approximately have $\mathcal{N}\sim\varphi_{\rm
i}^{2}/8M_{\rm Pl}^{2}$ and the change in the parameter $R$ is easily computed
to be

$\frac{R_{\rm e}}{R_{\rm i}}\simeq\left(\frac{\varphi_{\rm i}}{\varphi_{\rm
e}}\right)^{2}\sim\mathcal{N}\,,$ (2.14)

which shows that $R$ changes as the logarithm of the conformal time. As we
will see, this running can be neglected compared to other terms involved in
the model in agreement with the general argument (2.6). Also, the value of $R$
at the start of inflation is obtained to be

$\displaystyle R_{\rm i}=\frac{{\bf{e}}^{2}\varphi_{\rm i}^{2}}{H_{\rm
i}^{2}}=\frac{{\bf{e}}^{2}}{8\pi^{2}{\cal{P}}\epsilon}\frac{\varphi_{\rm
i}^{2}}{M_{\rm
Pl}^{2}}\simeq\frac{{\bf{e}}^{2}\mathcal{N}^{2}}{\pi^{2}{\cal{P}}}\,.$ (2.15)

Thus, as in the case of hilltop inflation, $R_{\rm i}$ can be very large
depending on the value of ${\bf{e}}$.

We comment that the standard symmetry breaking potential Eq. (2.7) may be used
as a working example for our inflationary dynamics. However, as we mentioned
before, other potentials can be used too as long as the potential has a non-
zero vev at the end of inflation to generate mass for the dark photon. In this
way, we can avoid the strong bounds from the Planck observation on potential
(2.7) such as the bound on the tensor-to-scalar power spectra ratio. As we
shall see in next sections the details of the relic density will not be
strongly dependent on the shape of the potential at the start of inflation.
Furthermore, as discussed above, we neglect the time dependence of $R$ in the
following analysis so hereafter we write $R$ without any reference to the
corresponding time (with the index ${\rm i}$ or ${\rm e}$).

Contrary to the model discussed in [25], where the mass of the vector field is
generated by the vev of the waterfall field after hybrid inflation, in the
present model the $U(1)$ symmetry breaking does not result in the formation of
cosmic strings. This is because, in the present model, a single value of the
phase $\theta$ (up to quantum fluctuations) is chosen by the inflationary
trajectory in the observable patch of the universe. As a result, we do not
have to impose the CMB constraint on the tension of cosmic strings.

In the following subsection we focus on the gauge field sector of the theory
to study the process of dark photon production.

###  2.2 Vector field dynamics

As we explained in the previous section, due to the Higgs mechanism, the dark
gauge field acquires mass during inflation which depends on the inflaton field
as given in Eq. (2.4). Therefore, the dark photon sector of the action (2.1)
would take the following form

$\displaystyle S_{\tiny A}=\int
d^{4}x\sqrt{-g}\left[-\frac{1}{4}f^{2}(\varphi)F_{\mu\nu}F^{\mu\nu}-\frac{1}{2}m_{\tiny
A}^{2}\,A_{\mu}A^{\mu}\right]\,.$ (2.16)

The inflationary background is given by the spatially flat Friedmann-Lemaître-
Robertson-Walker (FLRW) metric as

$\displaystyle{\rm d}s^{2}=a^{2}(\tau)\big{(}-{\rm
d}\tau^{2}+\delta_{ij}~{}{\rm d}x^{i}{\rm d}x^{j}\big{)}\,,$ (2.17)

where $a(\tau)$ is the scale factor and $\tau$ is the conformal time. Taking
the advantage of the rotational symmetry of the spatial part of the metric, we
decompose the spatial components of the gauge field as

$\displaystyle A_{i}={\partial}_{i}\chi+A^{T}_{i}\,,\hskip
28.45274pt\delta^{ij}{\partial}_{i}A^{T}_{j}=0\,,$ (2.18)

where $\chi$ and $A^{T}_{i}$ characterize the longitudinal and transverse
parts of the gauge field respectively. The component $A_{0}$ is non-dynamical
and it can be solved in terms of the longitudinal mode as [25]

$A_{0}=\Big{(}\frac{f^{2}\nabla^{2}}{f^{2}\nabla^{2}-m_{\tiny
A}^{2}a^{2}}\Big{)}\chi^{\prime}\,,$ (2.19)

where a prime denotes the derivative with respect to the conformal time.

Substituting the above solution in the action, it can then be decomposed into
the transverse and longitudinal parts as

$S_{\tiny A}=S_{T}+S_{L}\,,$ (2.20)

where

$\displaystyle S_{T}$
$\displaystyle=\frac{1}{2}\int\differential[3]{x}\differential{\tau}\Big{[}f^{2}(A^{T}_{i}{}^{\prime})^{2}-f^{2}({\partial}_{i}A^{T}_{j})^{2}-m_{\tiny
A}^{2}a^{2}(A^{T}_{i})^{2}\Big{]}\,,$ (2.21) $\displaystyle S_{L}$
$\displaystyle=\frac{1}{2}\int\differential[3]{x}\differential{\tau}m_{\tiny
A}^{2}a^{2}\Big{[}\chi^{\prime}\Big{(}\frac{f^{2}\nabla^{2}}{f^{2}\nabla^{2}-m_{\tiny
A}^{2}a^{2}}\Big{)}\chi^{\prime}-({\partial}_{i}\chi)^{2}\Big{]}\,.$ (2.22)

Finally, the total energy density of the gauge field during inflation can be
decomposed as $\rho_{\tiny A}=\rho_{T}+\rho_{L}$ where $\rho_{T}$ and
$\rho_{L}$ are the energy densities for the transverse and longitudinal modes
respectively given by (see Ref. [24, 25] for the details)

$\rho_{T}=\frac{f^{2}}{2a^{4}}\Big{[}(A^{T}_{i}{}^{\prime})^{2}+({\partial}_{i}A^{T}_{j})^{2}+\frac{m_{\tiny
A}^{2}a^{2}}{f^{2}}(A^{T}_{i})^{2}\Big{]}\,,\quad\rho_{L}=\frac{m_{\tiny
A}^{2}}{2a^{2}}\Big{[}\chi^{\prime}\big{(}\frac{f^{2}\nabla^{2}}{f^{2}\nabla^{2}-m_{\tiny
A}^{2}a^{2}}\big{)}\chi^{\prime}+({\partial}_{i}\chi)^{2}\Big{]}\,.$ (2.23)

The energy density of the vector field is the accumulated energy of its
quantum fluctuations for both transverse and longitudinal modes which are
produced during inflation.

The quantization procedure follows almost exactly the standard textbooks. For
the transverse modes we define the canonical field $V_{i}\equiv fA^{T}_{i}$
and expand it in terms of the creation and annihilation operators
$a^{\dagger}_{\mathbf{k},\lambda}$ and $a_{\mathbf{k},\lambda}$ as usual

$V_{i}=\sum_{\lambda}\int\frac{\differential[3]{k}}{(2\pi)^{3/2}}\varepsilon^{\lambda}_{i}({\bf
k})\big{[}v_{k}(\tau)a_{\mathbf{k},{\lambda}}+v_{k}(\tau)^{*}a^{\dagger}_{-\mathbf{k},{\lambda}}\big{]}e^{i\mathbf{k}.\mathbf{x}}\,,\hskip
28.45274pt[a_{\mathbf{k},{\lambda}},a^{\dagger}_{\mathbf{k}^{\prime},{\lambda^{\prime}}}]=\delta_{\lambda\lambda^{\prime}}\delta(\mathbf{k}-\mathbf{k}^{\prime})\,,$
(2.24)

where the time dependence of the vector field is encoded in the mode function
$v_{k}(\tau)$. In the above relation, $\varepsilon^{\lambda}_{i}({\bf k})$ are
the polarization vectors for $\lambda=\pm$ (see appendix A of Ref. [13] for
properties of polarization vectors). The mode function satisfies the equation
of motion

$\displaystyle
v_{k}^{\prime\prime}+\Big{(}k^{2}-\frac{f^{\prime\prime}}{f}+\frac{m_{\tiny
A}^{2}a^{2}}{f^{2}}\Big{)}v_{k}=0\,,$ (2.25)

subject to the Bunch-Davies initial condition.

Similarly, for the longitudinal mode we define the canonical field

$\displaystyle U\equiv h\,\chi\,,\qquad h\equiv m_{\tiny
A}af\sqrt{\frac{\nabla^{2}}{f^{2}\nabla^{2}-m_{\tiny A}^{2}a^{2}}}\,.$ (2.26)

Then we expand the field $U$ in terms of the creation and annihilation
operators $b^{\dagger}_{\mathbf{k}}$ and $b_{\mathbf{k}}$ as follows

$U=\int\frac{\differential[3]{k}}{(2\pi)^{3/2}}\big{[}u_{k}(\tau)b_{\mathbf{k}}+u_{k}(\tau)^{*}b^{\dagger}_{-\mathbf{k}}\big{]}e^{i\mathbf{k}.\mathbf{x}}\,,\hskip
28.45274pt[b_{\mathbf{k}},b_{\mathbf{k}^{\prime}}^{\dagger}]=\delta(\mathbf{k}-\mathbf{k}^{\prime})\,,$
(2.27)

and the equation for the mode function is easily found to be

$\displaystyle
u_{k}^{\prime\prime}+\Big{(}k^{2}-\frac{h^{\prime\prime}}{h}+\frac{m_{\tiny
A}^{2}a^{2}}{f^{2}}\Big{)}u_{k}=0\,,$ (2.28)

which is also subject to the Bunch-Davies initial condition.

Finally, the energy density of the transverse modes and the longitudinal mode
can be computed from Eq. (2.23) in terms of the corresponding mode functions

$\rho_{T}=\frac{1}{a^{4}}\int\frac{\differential[3]{k}}{(2\pi)^{3}}\Big{[}\big{|}{f\big{(}\frac{v_{k}}{f}\big{)}^{\prime}}\big{|}^{2}+\big{(}k^{2}+\frac{m_{\tiny
A}^{2}a^{2}}{f^{2}}\big{)}\absolutevalue{v_{k}}^{2}\Big{]}\,,$ (2.29)
$\rho_{L}=\frac{1}{2a^{4}}\int\frac{\differential[3]{k}}{(2\pi)^{3}}\Big{[}\big{|}{h\big{(}\frac{u_{k}}{h}\big{)}^{\prime}}\big{|}^{2}+\big{(}k^{2}+\frac{m_{\tiny
A}^{2}a^{2}}{f^{2}}\big{)}\absolutevalue{u_{k}}^{2}\Big{]}\,.$ (2.30)

We will be interested in the total energy density at the end of inflation so
it is useful to parameterize the energy density at that time as follows

$\rho_{T}(\tau_{\rm e})=\mathcal{C}_{T}\,H_{\rm
e}^{4}\,,\qquad\rho_{L}(\tau_{\rm e})=\mathcal{C}_{L}\,H_{\rm e}^{4}\,,$
(2.31)

where $H_{\rm e}$ is the Hubble expansion rate at the end of inflation and
$\mathcal{C}_{T}$ and $\mathcal{C}_{L}$ are dimensionless functions which
depend on the parameters of our model. The energy density of the produced dark
photons at the end of inflation is $\rho_{\rm e}=\rho_{T}(\tau_{\rm
e})+\rho_{L}(\tau_{\rm e})$ and we similarly define $\rho_{\rm
e}\equiv\mathcal{C}\,H_{\rm e}^{4}$ with
$\mathcal{C}\equiv\mathcal{C}_{T}+\mathcal{C}_{L}$.

So far in our analysis, we did not specify the form of the gauge kinetic
function $f$. In the following section, we assume a specific form for the
gauge kinetic function and derive an explicit expression for $\mathcal{C}_{T}$
and $\mathcal{C}_{L}$. Then we show that such a model can account for the
observed dark matter energy density.

##  3 The phenomenological model

As we mentioned before, during inflation the gauge field energy density
usually decays quickly. This issue can be remedied if energy is pumped from
the inflaton field to the gauge field through the coupling $f(\varphi)$. Since
in single clock inflation the inflaton field can be exchanged by time, the
coupling $f$ can be written equivalently as a function of time. To simplify
the analysis we assume that the time-dependence is of the form of a power law
as

$\displaystyle f(\tau)=\left(\dfrac{\tau}{\tau_{\rm e}}\right)^{p}\,,$ (3.1)

where we have normalized the coupling so that $f(\tau_{\rm e})=1$. In order to
have a well-defined perturbative theory, we require $f^{-1}<1$ during
inflation which implies $p>0$. To simplify the analysis further, we assume
that after inflation the dependence of $f$ on $\varphi$ is sufficiently weak
such that we can safely ignore any particle production during (p)reheating.
Note that the gauge kinetic function in the form of Eq. (3.1) has been
employed extensively in the context of magnetogenesis [31, 32, 19, 33, 34, 35,
36, 37, 38, 39, 40, 41], anisotropic inflation [20, 26, 28, 27, 43, 44, 29,
30], pseudoscalar inflaton [19, 45], and inflation with multiple vector fields
[46, 47, 48, 49].

In the following we solve the mode functions for the transverse and
longitudinal modes which will be used to calculate their contributions to the
relic density of dark matter.

###  3.1 Transverse modes

Assuming the specific form (3.1) for the gauge kinetic function, the equation
for the transverse mode function Eq. (2.25) takes the following form

$\displaystyle
v_{k}^{\prime\prime}+\Big{(}k^{2}-\frac{p(p-1)}{\tau^{2}}+\frac{R}{\tau^{2}}\left(\frac{\tau_{\rm
e}}{\tau}\right)^{2p}\Big{)}v_{k}=0\,.$ (3.2)

As we have discussed previously the time dependence of the parameter $R$ can
be ignored up to slow roll corrections and it can be treated as a constant
throughout the period of inflation. Note that the value of $R$ can be much
larger than unity.

From Eq. (3.2) we see the competition between the last two terms in the big
bracket. The term containing $p(p-1)/\tau^{2}$ is from the gauge kinetic
function and for the particular case $p=2$ yields the usual term $2/\tau^{2}$
in the equation of motion for the mode function of a massless scalar field
which is the hallmark of the scale invariant perturbations in a fixed de
Sitter background. If $R$ is chosen large enough, then the last term can
dominate over the second term towards the final stage of inflation. This
happens at the time $\tau_{\rm c}$ when the two terms become equal [36]

$\frac{\tau_{\rm c}}{\tau_{\rm
e}}=\Big{(}\frac{R}{p(p-1)}\Big{)}^{\frac{1}{2p}}\,.$ (3.3)

If $R$ is large enough then the mass term of the dark photon plays important
roles in the evolution of mode functions. Specifically, assuming $R>p(p-1)$
hereafter, the dominance of the mass term happens sometime during inflation.
To have a sense as when the effect of mass becomes important it is better to
express Eq. (3.3) in terms of the number of e-folds. Defining the number of
e-folds from the time when the condition (3.3) is met till the time of end of
inflation by $\Delta\mathcal{N}$, we obtain

$\displaystyle\Delta\mathcal{N}\equiv\ln\frac{\tau_{\rm c}}{\tau_{\rm
e}}=\frac{1}{2p}\ln\Big(\frac{R}{p(p-1)}\Big{missing})\,.$ (3.4)

The larger $\Delta\mathcal{N}$ is, the sooner the effect of mass becomes
important. For example for $p=2$ and $R=10^{12}$ we obtain
$\Delta\mathcal{N}\simeq 7$, so only towards the last seven e-folds the mass
term becomes important. Of course, with smaller values of $R$,
$\Delta\mathcal{N}$ becomes smaller.

Technically speaking, the effect of the mass in the mode function is captured
by the combination111The vector field defined in the action (2.16) is not
canonically normalized. After defining the canonical field all couplings
receive a factor of $f^{-1}$. $m_{A}^{2}a^{2}/f^{2}$ as is evident from Eq.
(2.25). Since $f$ is very large initially and become $\order{1}$ only at the
end of inflation, the effects of the mass of dark photon are subleading during
inflation so they can be excited efficiently during most of the period of
inflation. Thus, even for dark photons much heavier than the Hubble scale with
$R\gg 1$, we can have efficient particle production during inflation. Finally,
note that $\tau_{\rm c}$ and $\Delta\mathcal{N}$ are independent of scale.
Hence, the domination of the mass term happens at the same time for all
momenta.

From the above discussion, one can ignore the last term in the brackets of Eq.
(3.2) at early times during inflation. Starting with the Bunch-Davies initial
condition, the solution during this stage is given by

$\displaystyle
v_{k}(\tau)=\frac{\sqrt{-\pi\tau}}{2}e^{i(1+2\nu)\pi/4}H_{\nu}^{(1)}(x)\,,\qquad\text{early
times, }-\tau\gg-\tau_{\rm c}\,,$ (3.5)

in which $\nu\equiv p-1/2$, $H_{\nu}^{(1)}$ is the Hankel function of the
first kind and its argument is $x\equiv-k\tau$. On the other hand, towards the
final stage of inflation the last term in Eq. (3.2) containing $R$ dominates
and the solution is given by [36]

$\displaystyle
v_{k}(\tau)=\frac{\sqrt{-\pi\tau}}{2}e^{i(1+2\nu)\pi/4}\left[c_{2}\,H_{\mu}^{(1)}(y)+c_{1}\,H_{\mu}^{(2)}(y)\right]\,,\qquad\text{late
times, }-\tau\ll-\tau_{\rm c}\,,$ (3.6)

where $\mu\equiv 1/2p$ and $c_{1}$ and $c_{2}$ are two constants of
integration and we have defined

$y\equiv\frac{\sqrt{R}}{p}\left(\frac{\tau_{\rm e}}{\tau}\right)^{p}\,.$ (3.7)

Matching the late time solution (3.6) with the early time solution (3.5) at
$\tau=\tau_{c}$ fixes the coefficients $c_{1}$ and $c_{2}$. More specifically,
demanding that both $v_{k}$ and $v_{k}^{\prime}$ to be continuous at
$\tau_{\rm c}$ yields

$\displaystyle c_{\ell}=\frac{\mp i\pi}{4p}\left[py_{\rm
c}H_{\nu}^{(1)}(x_{\rm c})\,H_{\mu}^{(\ell)}{}^{\prime}(y_{\rm c})+x_{\rm
c}\,H_{\nu}^{(1)}{}^{\prime}(x_{\rm c})\,H_{\mu}^{(\ell)}(y_{\rm
c})\right]\,,\qquad\ell=1,2$ (3.8)

in which the upper (lower) sign is for $\ell=1$ ($\ell=2$) and the primes
denote derivatives with respect to the corresponding arguments. Also from the
above definitions we have

$x_{\rm c}=-k\tau_{\rm c}\,,\qquad y_{\rm
c}=\frac{\sqrt{R}}{p}\Big{(}\frac{\tau_{\rm e}}{\tau_{\rm
c}}\Big{)}^{p}=\sqrt{\frac{p-1}{p}}\,,$ (3.9)

where in the last equality we have used Eq. (3.3). Note that the coefficients
$c_{\ell}$ depend on $k$ through the dependence on $x_{\rm c}$. Fig. 2
compares our analytic solution for the mode function with the numerical
solution of Eq. (3.2). The plot shows that our analytic expression is in good
agreement with the numerical result.

Figure 2: Comparison between the numerical solution of Eq. (3.2) and the
analytic expression from matching Eqs. (3.5) and (3.6) at $\tau=\tau_{\rm c}$.
We have set $R=10^{12}$, $p=2.2$ and $-k\tau_{\rm e}=10^{-20}$. The latter has
been chosen such that the mode lies between $k_{\rm min}$ and $k_{\rm max}$.

We are interested in the range of momenta for which dark photons can be
produced via excitation of quantum fluctuations. This happens when the
adiabatic condition for the mode function is broken. Most efficiently this
happens if the mode function undergoes a tachyonic growth. Roughly speaking,
the only chance for the mode function to become tachyonic is until the time
$\tau_{\rm c}$ when the middle term in the brackets of Eq. (3.2) is dominant
while after that the mode cannot experience tachyonic growth. As a result, we
define

$k_{\rm min}\equiv\frac{\sqrt{p(p-1)}}{-\tau_{\rm i}}\,,\hskip
28.45274ptk_{\rm max}\equiv\frac{\sqrt{p(p-1)}}{-\tau_{\rm c}}\,,$ (3.10)

corresponding to the mode that becomes tachyonic at the beginning of inflation
$\tau_{\rm i}$ and at the time of $\tau_{\rm c}$ respectively. Hereafter, we
assume that $\tau_{\rm i}<\tau_{\rm c}$ (i.e. $k_{\rm min}<k_{\rm max}$) so
that the tachyonic growth occurs during inflation. We also assume that
$R>p(p-1)$ so that $\tau_{\rm c}<\tau_{\rm e}$. Thus, dark photons with
momenta in the range between $k_{\rm min}$ and $k_{\rm max}$ can be produced
during inflation.

We need to compute the energy density of the transverse modes at the end of
inflation. We use Eq. (2.29), which for clarity is repeated as

$\rho_{T}(\tau_{\rm e})\approx\frac{1}{a_{\rm e}^{4}}\int_{k_{\rm
min}}^{k_{\rm
max}}\frac{\differential[3]{k}}{(2\pi)^{3}}\Big{[}\big{|}{\big{(}\frac{v_{k}}{f}\big{)}^{\prime}}\big{|}^{2}+m_{\tiny
A}^{2}a_{\rm e}^{2}\absolutevalue{v_{k}}^{2}\Big{]}_{\tau_{\rm e}}\,.$ (3.11)

Here, we have integrated only over the range of momenta specified in Eq.
(3.10). Note also that we have neglected a term containing $k^{2}$ compared to
$m_{\tiny A}^{2}a_{\rm e}^{2}$ which is always valid after $\tau_{\rm c}$ for
the range given by Eq. (3.10).

In the following we find an approximate closed form for the energy density at
the end of inflation. The matching time $\tau_{\rm c}$ does not depend on the
momentum $k$ and assuming it happens during inflation (i.e. $\tau_{\rm
i}<\tau_{\rm c}<\tau_{\rm e}$), we can use the analytic expression (3.6) for
the late time. In the expression (3.6), the only $k$ dependence comes from the
coefficients $c_{\ell}$ ($\ell=1,2$) given in Eq. (3.8) through the variable
$x_{\rm c}=-k\tau_{\rm c}$. From the range of momenta given in Eq. (3.10),
$x_{\rm c}$ which is the argument of $H^{(1)}_{\nu}$ in $c_{\ell}$, can take
values between

$\sqrt{p(p-1)}e^{-(\mathcal{N}-\Delta\mathcal{N})}<x_{\rm c}<\sqrt{p(p-1)}\,.$
(3.12)

Further, the argument of the Hankel functions $H^{(\ell)}_{\mu}$ in $c_{\ell}$
is $y_{\rm c}$ which is always smaller than unity. As a result, for the
purpose of analytic studies it is reasonable to use the small argument
approximation for all the Hankel functions as follows

$H^{(1)}_{\alpha}(x)\simeq-i\frac{\Gamma(\alpha)}{\pi}\left(\frac{x}{2}\right)^{-\alpha}\,,\qquad
H_{\alpha}^{(1)}{}^{\prime}(x)\simeq
i\frac{\Gamma(\alpha+1)}{2\pi}\left(\frac{x}{2}\right)^{-\alpha-1}\,,$ (3.13)

with $\alpha>0$. For $H^{(2)}_{\alpha}$ similar expressions can be obtained
noting that $H^{(2)}_{\alpha}={H^{(1)}_{\alpha}}^{*}$.

Using the above approximations and after some lines of algebra we find the
following approximate expressions for the coefficients (3.8)

$c_{1}\simeq
c_{2}\simeq-i\frac{\Gamma(\nu)\Gamma(\mu)}{4\pi}\left(\frac{x_{\rm
c}}{2}\right)^{-\nu}\left(\frac{y_{\rm c}}{2}\right)^{-\mu}\,.$ (3.14)

Plugging these values of $c_{\ell}$ in the final mode function in Eq. (3.11)
for $\rho_{T}(\tau_{\rm e})$ and using the definition Eq. (2.31) we obtain

$\mathcal{C}_{T}\simeq\frac{e^{-2\nu\Delta\mathcal{N}}}{2\pi^{2}}\Big{(}\frac{\Gamma(\nu)\Gamma(\mu)}{2^{2-\nu-\mu}\pi}\Big{)}^{2}\Big{(}\frac{p}{p-1}\Big{)}^{\frac{1}{2p}}\int_{x_{\rm
min}}^{x_{\rm max}}\differential{x}(L_{\rm e}+R\,K_{\rm e})x^{2-2\nu}\,,$
(3.15)

where as before, the index ${\rm e}$ means the function is evaluated at
$\tau_{\rm e}$. Here, we have defined $L(\tau)$ and $K(\tau)$ which are
functions of time and not momentum $k$ as follows

$\displaystyle
L(\tau)\equiv\pi\big{[}pyJ_{\mu-1}(y)+(p-1)J_{\mu}(y)\big{]}^{2}\,,\hskip
28.45274ptK(\tau)\equiv\pi J_{\mu}(y)^{2}\,,$ (3.16)

where $J_{\mu}$ is the Bessel function of the first kind, and the integration
domain is given by

$x_{\rm min}=\sqrt{p(p-1)}e^{-\mathcal{N}}\,,\qquad x_{\rm
max}=\sqrt{p(p-1)}e^{-\Delta\mathcal{N}}\,.$ (3.17)

For the arguments of the Bessel functions in Eq. (3.16) at the end of
inflation, we have $y_{\rm e}=\sqrt{R}/p>1$. Note that we are mostly
interested in the regime where $y_{\rm e}\gg 1$ (i.e. $R\gg p^{2}$) so we can
use the large argument approximation of the Bessel functions to write

$L_{\rm e}\simeq
p\sqrt{R}\Big{[}1-\sin\big(\frac{4\sqrt{R}-\pi}{2p}\big{missing})\Big{]}\,,\qquad
K_{\rm
e}\simeq\frac{p}{\sqrt{R}}\Big{[}1+\sin\big(\frac{4\sqrt{R}-\pi}{2p}\big{missing})\Big{]}\,,$
(3.18)

which yields $L_{\rm e}+R\,K_{\rm e}\simeq 2p\sqrt{R}$.

As a result, from Eq. (3.15) we obtain

$\mathcal{C}_{T}\simeq\frac{p[p(p-1)]^{5/2-p}}{\pi^{4}(4-2p)}\left(\frac{\Gamma(\nu)\Gamma(\mu)}{2^{2-\nu-\mu}}\right)^{2}\left(\frac{p}{p-1}\right)^{\frac{1}{2p}}\begin{cases}-e^{(2p-4)\mathcal{N}}e^{-(p-1)\Delta\mathcal{N}}\,,\hskip
14.22636ptp>2\\\ e^{-(3-p)\Delta\mathcal{N}}\,,\hskip
41.25641pt1<p<2\end{cases}\,,$ (3.19)

where we have used Eq. (3.4) to express the factor of $\sqrt{R}$ in terms of
$\Delta\mathcal{N}$. This is our final expression for the energy density of
the transverse modes of the gauge field at the end of inflation.

###  3.2 Longitudinal mode

Now, we compute the same quantity for the longitudinal mode. The equation of
motion for the longitudinal mode is given by Eq. (2.28) which we repeat as

$u_{k}^{\prime\prime}+\Big{(}k^{2}-\frac{h^{\prime\prime}}{h}+\frac{m_{\tiny
A}^{2}a^{2}}{f^{2}}\Big{)}u_{k}=0\,,\qquad\quad h=\frac{m_{\tiny
A}ak}{\sqrt{k^{2}+m_{\tiny A}^{2}a^{2}/f^{2}}}\,.$ (3.20)

Here, we have written $h$, already defined in (2.26), in momentum space. The
most straightforward way to find a solution for (3.20) is to study its
solutions in the two regimes $k^{2}\gg\frac{m_{\tiny A}^{2}a^{2}}{f^{2}}$ and
$k^{2}\ll\frac{m_{\tiny A}^{2}a^{2}}{f^{2}}$ separately.

For short modes $k^{2}\gg\frac{m_{\tiny A}^{2}a^{2}}{f^{2}}$, we have
$h\approx m_{\tiny A}a$ which leads to $h^{\prime\prime}/h\approx 2/\tau^{2}$
and Eq. (3.20) simplifies to

$u_{k}^{\prime\prime}+\Big{(}k^{2}-\frac{2}{\tau^{2}}\Big{)}u_{k}=0\,.$ (3.21)

Starting with a Bunch-Davies vacuum the solution is given by

$u_{k}=\frac{\sqrt{-\pi\tau}}{2}e^{i(1+2{\tilde{\nu}})\pi/4}H^{(1)}_{{\tilde{\nu}}}(x)\,,$
(3.22)

where ${\tilde{\nu}}\equiv 3/2$ and as before the argument of the Hankel
function is $x=-k\tau$.

For long modes $k^{2}\ll\frac{m_{\tiny A}^{2}a^{2}}{f^{2}}$, we have $h\approx
kf$ which gives $h^{\prime\prime}/h=f^{\prime\prime}/f\approx p(p-1)/\tau^{2}$
and Eq. (3.20) simplifies to

$u_{k}^{\prime\prime}+\Big{(}-\frac{p(p-1)}{\tau^{2}}+\frac{R}{\tau^{2}}\Big{(}\frac{\tau_{\rm
e}}{\tau}\Big{)}^{2p}\Big{)}u_{k}=0\,,$ (3.23)

which has the following solution

$u_{k}=\frac{\sqrt{-\pi\tau}}{2}e^{i(1+2{\tilde{\nu}})\pi/4}\Big{[}{\tilde{c}}_{2}H^{(1)}_{{\tilde{\mu}}}(y)+{\tilde{c}}_{1}H^{(2)}_{{\tilde{\mu}}}(y)\Big{]};\hskip
28.45274pt{\tilde{\mu}}\equiv 1-\frac{1}{2p}\,.$ (3.24)

Note that the above expression is an exact solution of Eq. (3.23) which, apart
from the $k^{2}$ term, is the same as Eq. (3.2) for the transverse modes. As a
result, similar to the transverse case, the second term in the big bracket of
(3.23) dominates over the first term around and after the time $\tau_{\rm c}$
defined in (3.3). Such behaviour is encoded in the analytic solution (3.24).

The two solutions (3.22) and (3.24) should be matched at the time
$\tilde{\tau}_{\rm c}$ defined by

$k^{2}=\frac{m_{\tiny
A}^{2}a^{2}}{f^{2}}\quad\implies\quad\frac{\tilde{\tau}_{\rm c}}{\tau_{\rm
e}}=\Big{(}\frac{\sqrt{R}}{-k\tau_{\rm e}}\Big{)}^{\frac{1}{p+1}}\,.$ (3.25)

As it is evident from the above definition, unlike $\tau_{\rm c}$, the time
$\tilde{\tau}_{\rm c}$ depends on the momentum. To go further, we identify two
types of scales. First, short modes, namely $-k\tau_{\rm e}>\sqrt{R}$, for
which the critical time $\tilde{\tau}_{\rm c}$ happens after inflation. For
these short scales the form of the mode function is well approximated by
(3.22) during inflation. Note that we are mostly interested in $R>1$. As a
result, these scales never exit the horizon and from Eq. (3.21) they are not
excited during inflation. Thus, we can safely ignore the contributions of
these scales in the energy density of the longitudinal mode. Second, we have
long modes $-k\tau_{\rm e}<\sqrt{R}$ for which the time $\tilde{\tau}_{\rm c}$
occurs during inflation. For these scales the late time behaviour is described
by Eq. (3.24). The unknown coefficients ${\tilde{c}}_{\ell}$ are obtained by
requiring the continuity of the mode function and its time derivative at
$\tilde{\tau}_{\rm c}$. Performing these matching conditions, we obtain

$\displaystyle\tilde{c}_{\ell}=\frac{\mp i\pi}{4p}\left[p\tilde{y}_{\rm
c}H_{\tilde{\nu}}^{(1)}(\tilde{x}_{\rm
c})\,H_{\tilde{\mu}}^{(\ell)}{}^{\prime}(\tilde{y}_{\rm c})+\tilde{x}_{\rm
c}\,H_{\tilde{\nu}}^{(1)}{}^{\prime}(\tilde{x}_{\rm
c})\,H_{\tilde{\mu}}^{(\ell)}(\tilde{y}_{\rm c})\right]\,,\qquad\ell=1,2$
(3.26)

where a prime denotes derivative with respect to the corresponding argument
and by definition we have

$\tilde{x}_{\rm c}\equiv-k\tilde{\tau}_{\rm c}\,,\qquad\tilde{y}_{\rm
c}\equiv\frac{\sqrt{R}}{p}\left(\frac{\tau_{\rm e}}{\tilde{\tau}_{\rm
c}}\right)^{p}=\frac{\tilde{x}_{\rm c}}{p}\,.$ (3.27)

One important difference compared to the transverse modes is that
$\tilde{y}_{\rm c}$ now depends on the momentum $k$ while the corresponding
quantity $y_{\rm c}$ for the transverse modes is independent of $k$.

Here we should comment that the analytic mode function obtained by matching
early and late time solution at $\tilde{\tau}_{\rm c}$ is not entirely an
accurate expression for the mode function. The reason is that for a time
interval around $\tilde{\tau}_{\rm c}$, the equation of motion is not in the
form given by neither (3.21) nor (3.23) but takes a more complicated form. We
have checked numerically that although the qualitative behaviour of our
analytic result is accurate, this causes an $\order{1}$ discrepancy in the
amplitude of the mode function at the end of inflation for the range of
parameters of our interest. It is possible to improve the analytic result by
solving the equation for the mode function around the time $\tilde{\tau}_{\rm
c}$ and then match the result with both early and late solutions given above.
However, the resulting expression would be much more complicated which is not
worth bothering the reader with. As a result, we continue exploiting the above
analytic expression understanding that for a high precision result one needs
to perform numerical analysis.

We are only interested in the range of momenta that experience growth through
breaking the adiabaticity condition which are given by (see appendix A for the
details)

$\displaystyle x_{\min}=\sqrt{2}\,e^{-{\cal N}}\,,\hskip
28.45274ptx_{\max}=\begin{cases}\sqrt{p(p-1)}e^{-\Delta\mathcal{N}}\hskip
78.24507ptp>2\\\
2^{\frac{p+1}{2p}}[p(p-1)]^{-\frac{1}{2p}}e^{-\Delta\mathcal{N}}\hskip
24.18501pt1<p<2\end{cases}\,,$ (3.28)

where $x_{\min}=-k_{\min}\tau_{\rm e}$ and $x_{\max}=-k_{\max}\tau_{\rm e}$ as
before and the upper bound corresponds to the momentum for which
$-k_{\max}\tilde{\tau}_{\rm c}(k_{\max})=\sqrt{2}$. Correspondingly, the
energy density of the longitudinal mode at the end of inflation (2.30) is
approximated by

$\rho_{L}(\tau_{\rm e})\approx\frac{1}{2a_{\rm e}^{4}}\int_{k_{\rm
min}}^{k_{\rm
max}}\frac{\differential[3]{k}}{(2\pi)^{3}}\Big{[}\big{|}{\big{(}\frac{u_{k}}{f}\big{)}^{\prime}}\big{|}^{2}+m_{\tiny
A}^{2}a_{\rm e}^{2}\absolutevalue{u_{k}}^{2}\Big{]}_{\tau_{\rm e}}\,,$ (3.29)

where we have used the fact $h(\tau_{\rm e})\simeq kf(\tau_{\rm e})=k$ and
also neglected the $k^{2}$ term in comparison to the mass term for the momenta
in the range given in Eq. (3.28). Note that the mode function at the end of
inflation is given by Eq. (3.24) in which the coefficients
${\tilde{c}}_{\ell}$ given in Eq. (3.26).

To proceed further, we note that from (3.28) we have

$\begin{cases}\sqrt{2}\,e^{-(\mathcal{N}-\Delta\mathcal{\tilde{N}})}<\tilde{x}_{\rm
c}<\sqrt{p(p-1)}e^{\Delta\mathcal{\tilde{N}}-\Delta\mathcal{N}}\hskip
78.24507ptp>2\\\
\sqrt{2}\,e^{-(\mathcal{N}-\Delta\mathcal{\tilde{N}})}<\tilde{x}_{\rm
c}<2^{\frac{p+1}{2p}}[p(p-1)]^{-\frac{1}{2p}}e^{\Delta\mathcal{\tilde{N}}-\Delta\mathcal{N}}\hskip
24.18501pt1<p<2\end{cases}\,,$ (3.30)

where $\Delta{\cal\tilde{N}}=\ln\frac{{\tilde{\tau}}_{\rm c}}{\tau_{\rm e}}$.
As a result, $\tilde{x}_{\rm c}$ is effectively a small quantity. More
precisely, since $\Delta{\cal\tilde{N}}-\Delta{\cal
N}=\ln\frac{{\tilde{\tau}}_{\rm c}}{\tau_{\rm c}}$, we have
$\Delta\mathcal{\tilde{N}}<\Delta\mathcal{N}$ and
$\Delta\mathcal{\tilde{N}}>\Delta\mathcal{N}$ for ${\tilde{\tau}}_{\rm
c}<\tau_{\rm c}$ and ${\tilde{\tau}}_{\rm c}>\tau_{\rm c}$ respectively.
Therefore, in the case of ${\tilde{\tau}}_{\rm c}<\tau_{\rm c}$ the
approximation $\tilde{x}_{\rm c}\ll 1$ is more accurate than the case
${\tilde{\tau}}_{\rm c}>\tau_{\rm c}$. So, similarly to the case of the
transverse modes, we use the small argument approximation of the Hankel
functions yielding

${\tilde{c}}_{1}\simeq{\tilde{c}}_{2}\simeq-i\Big{(}\frac{p+1}{p}\Big{)}\frac{\Gamma({\tilde{\nu}})\Gamma({\tilde{\mu}})}{4\pi}\left(\frac{{\tilde{x}}_{\rm
c}}{2}\right)^{-{\tilde{\nu}}}\left(\frac{{\tilde{y}}_{\rm
c}}{2}\right)^{-{\tilde{\mu}}}\,.$ (3.31)

Performing a similar analysis which led to Eq. (3.19) we obtain

$\mathcal{C}_{L}\simeq\frac{p(p+1)^{3}}{\pi^{4}(4-2p)}\Big{(}\frac{\Gamma(\tilde{\nu})\Gamma(\tilde{\mu})}{2^{2-\tilde{\nu}-\tilde{\mu}}}\Big{)}^{2}\frac{\big{(}p(p-1)\big{)}^{\frac{1}{2}-\frac{3}{2p}}}{p^{\frac{1}{p}}2^{-\frac{2p+5}{1+p}}}\begin{cases}-[p(p-1)]^{\frac{2-p}{p(p+1)}}e^{\frac{2p-4}{p+1}\mathcal{N}}e^{\frac{p^{2}-4p+1}{p+1}\Delta\mathcal{N}}\hskip
17.07182ptp>2\\\ 2^{\frac{2-p}{p(p+1)}}e^{-(3-p)\Delta{\cal N}}\hskip
85.35826pt1<p<2\end{cases},$ (3.32)

where $\mathcal{C}_{L}$ is defined in Eq. (2.31) as $\rho_{L}(\tau_{\rm
e})=\mathcal{C}_{L}H_{\rm e}^{4}$.

##  4 Dark matter relic density

In the previous section we obtained the energy density of the gauge field for
both transverse and longitudinal modes at the end of inflation. In this
section we use the above results to compute its present energy density after
getting redshifted after inflation due to cosmic expansion.

First of all, note that the energy of the produced dark photons during
inflation must be sub-dominant compared to the inflationary background.
Besides, the backreaction of the gauge field on the inflaton slow-roll
dynamics must be negligible. As discussed in Refs. [24, 25], in our
phenomenological model with the coupling (3.1), the backreaction condition on
the inflaton field’s slow-roll dynamics is stronger. This condition constrains
the value of $p$. For instance, similar to [25], for $\mathcal{N}=50$ e-folds
of inflation one must have $p\lesssim 2.2$ in order to avoid any backreaction.

In addition, there are strong constraints from the dark matter isocurvature
perturbations. In this setup the origin of dark matter is different than the
origin of baryonic matter and radiation which are generated from the decay of
inflaton at the end of inflation. Therefore, dark matter isocurvature will
inevitably be generated. This constraint was studied in [25] yielding to the
bound $p<2.01$ in our setup with the coupling (3.1). Interestingly, we see
that the constraint on $p$ from the isocurvature perturbation is stronger than
the backreaction condition.

Assuming that the above conditions are satisfied, we must compute the present
time energy density of the gauge field. Note that since we assume a large
value for the parameter $R$, the dark photons become non-relativistic right
after inflation and their energy density decrease like $a^{-3}$ in the
subsequent cosmic history. The present energy density of dark photon can be
approximated as [25]

$\Omega_{A,0}=0.5\times 10^{-4}\left(\frac{T_{\rm r}}{10^{12}{\rm
GeV}}\right)^{5}\mathcal{C}(p,{\cal N},\Delta{\cal N})\,,$ (4.1)

where $\mathcal{C}=\mathcal{C}_{T}+\mathcal{C}_{L}$ in which $\mathcal{C}_{T}$
and $\mathcal{C}_{L}$ are given by Eqs. (3.19) and (3.32) respectively and
$T_{\rm r}$ is the reheating temperature. In deriving above result we have
assumed an instantaneous reheating for simplicity.

Figure 3: Left: present energy density of dark photons given by Eq. (4.1) as a
function of $p$ for $R=10$ (dashed-dotted curves) and $R=10^{6}$ (solid
curves) for different reheat temperatures with $\mathcal{N}=50$. The solid
black horizontal line represents the observed dark matter density
$\Omega_{D}=0.27$. Right: comparison of the contributions of the transverse
and longitudinal modes to the energy density for $R=10^{6}$.

We are mostly interested in the case $R>1$ where the dark matter is non-
relativistic right after inflation. The left panel of Fig. 3 shows the relic
energy density of dark photons for two values of $R=10$ (dashed-dotted curves)
and $R=10^{6}$ (solid curves) for different reheat temperatures. We see that
the dependence of the relic energy density on $R$ is stronger for $p<2$ than
$p>2$. For a fixed value of $p$, the reheat temperature is higher for larger
values of $R$ in order for the dark photons furnish the observed dark matter
relic density. If the dark photon comprises the whole dark matter energy
density, the reheat temperature is typically in the range
$10^{12}\mathrm{GeV}<T_{\mathrm{r}}<10^{14}\mathrm{GeV}$.

In the right panel of Fig. 3 we have compared the contribution of the
longitudinal mode in the total energy density to that of the transverse modes.
As we see the contribution from the longitudinal mode for $p>2$ is negligible
in agreement with the results of [24, 25]. However, for $1<p<2$, we see that
the longitudinal mode is dominant. This is a new result and roughly in
agreement with the results of Ref. [3] where it is shown that for $p=0$, the
observed dark matter relic can be achieved from the longitudinal mode. But,
let us look more carefully to the origin of these results. Note that one main
difference between the longitudinal mode and transverse modes is that the
particle production for the longitudinal mode is governed initially by the
term $-2/\tau^{2}$ and then the term $-p(p-1)/\tau^{2}$ while for the
transverse modes this is always the term $-p(p-1)/\tau^{2}$ which is
responsible for growth. For $p>2$, the factor $-2/\tau^{2}$ has less weight
compared to the factor $-p(p-1)/\tau^{2}$ and thus the longitudinal mode
experiences less growth compared to the transverse modes. However for $1<p<2$,
the stronger growth happens for the longitudinal mode.

Throughout this paper, we always assumed that $R>1$ and we did not study the
case $R<1$ as it is already studied in the literature [22, 23, 24, 25].
However, some comments for the case $R<1$ are in order. For $p>2$ and $R<1$,
it is shown that one can always neglect the longitudinal mode [24] while this
scenario is restricted to $p<2.01$ from the bound on the isocurvature
perturbations [25]. In the case of $1<p<2$ and $R<1$ the short scale modes
with $-k\tau_{\rm e}>\sqrt{R}$ give the main contribution to the energy
density of the longitudinal mode and it can be shown that this contribution is
also larger than the contribution coming from the transverse modes. Therefore,
we can have a scenario of light vector dark matter in which the whole
contribution comes from the longitudinal mode. This scenario is qualitatively
similar to the model considered in Ref. [3] which corresponds to the special
case of $p=0$.

##  5 Summary and discussion

In this paper we considered a new mechanism of vector dark matter production
during inflation. The setup contains a complex inflaton field which is charged
under the dark $U(1)$ gauge field. As the inflaton field rolls towards the
global minimum of the potential the gauge field acquires mass due to the Higgs
mechanism. In order for the gauge field to play the roles of dark matter we
require the potential to have a global minimum with a non-zero vacuum
expectation value of the field. The accumulated energy density associated with
the excited quantum fluctuations of the transverse and the longitudinal modes
play the roles of the observed dark matter energy density. We mostly
concentrated in the cases where the vector field is massive during inflation
with $R>1$ so the dark photons become non-relativistic right after inflation.
This is motivated from simple models of inflation such as the Higgs symmetry
breaking potential which is under the category of large field inflation, i.e.
fields value larger than $M_{\rm{Pl}}$. However, there is no fundamental
reason for the potential be in the category of large field setup, so one may
consider potentials in which the field displacement is not large so a smaller
value of $R$ is also possible.

In order to get energy from inflaton and prevent the decay of the gauge field
energy density, we have added a coupling between the inflaton field and the
gauge field through the gauge kinetic function. Such a coupling not only pumps
energy from the inflaton sector to the gauge field sector but also reduces the
effective mass of the vector field so the quanta of the vector field can be
excited efficiently during most of the period of inflation. An interesting
feature of this setup is that the vector field is massive during the entire
period of inflation so the longitudinal mode also contributes to the dark
matter relic density. For a gauge kinetic function with a power law time
dependence with the free parameter $p$ we have calculated the relic energy
density associated with the transverse and the longitudinal modes. Requiring a
negligible backreaction on the background inflaton dynamics and imposing the
dark matter isocurvature constraint implies that $p<2.01$. We have shown that
the observed dark matter relic density can be generated with the reheat
temperature in the range
$10^{12}\mathrm{GeV}<T_{\mathrm{r}}<10^{14}\mathrm{GeV}$. For a fixed value of
$p$ models with larger values of $R$ require higher reheat temperature to
obtain the observed dark matter relic density.

The current analysis can be extended in various directions. First, one can
consider the setup with $R<1$ so the vector field is not non-relativistic at
the end of inflation. In this case we have to wait until the time
$t_{\mathrm{NR}}$ when the mode function becomes non-relativistic and from
then the energy density of the dark photon can be carried to the time of
matter and radiation equality as the seed of dark matter. This case was
studied in details in [25] and a similar analysis can be performed here as
well. The second direction is to look for the dark photon productions during
(p)reheating. More specifically, the oscillations of the inflaton field near
its minimum at the end of inflation induce an oscillatory coupling to the
gauge field via the gauge kinetic function $f(\varphi)$. This can excite the
quanta of dark photon via parametric resonance. This channel of particle
production is more significant for the case of $R<1$ where the dark photon is
light while such particle creation is expected to be negligible for the case
$R>1$ as the excitations of heavy dark photon may not be efficient at the end
of inflation. This is an interesting question which is beyond the scope of
this work. Finally, while we have considered the inflaton as a complex scalar
field whose vev gives mass to the dark photon, the analysis in the present
paper can be applied to more general setups in which $\phi$ is not necessarily
the inflaton as far as the $U(1)$ gauge symmetry is broken for the whole
period of inflation. Therefore it is certainly intriguing to explore other
possibilities.

Acknowledgments: The work of M.A.G. was supported by Japan Society for the
Promotion of Science Grants-in-Aid for international research fellow No.
19F19313. The work of S.M. was supported in part by Japan Society for the
Promotion of Science Grants-in-Aid for Scientific Research No. 17H02890, No.
17H06359, and by World Premier International Research Center Initiative, MEXT,
Japan.

##  Appendix A Growing momenta for longitudinal mode

As we discussed in the main text, the short modes with $-k\tau_{\rm
e}>\sqrt{R}$ do not grow and we do not consider them further here. For the
modes with $-k\tau_{\rm e}<\sqrt{R}$ the critical time $\tilde{\tau}_{\rm c}$
occurs during inflation. We define the number of e-folds from
$\tilde{\tau}_{\rm c}$ until the time of end of inflation from (3.25) as

$\displaystyle\Delta\tilde{\mathcal{N}}\equiv\ln\frac{{\tilde{\tau}}_{\rm
c}}{\tau_{\rm e}}=\frac{1}{p+1}\ln(\frac{\sqrt{R}}{-k\tau_{\rm e}})\,.$ (A.1)

Depending on the parameters, $\tilde{\tau}_{\rm c}$ can occur before or after
$\tau_{\rm c}$. To see this, we look at

$\Delta\mathcal{N}-\Delta\tilde{\mathcal{N}}=\ln\frac{\tau_{\rm
c}}{{\tilde{\tau}}_{\rm c}}=\frac{1}{p+1}\ln(\frac{-k\tau_{\rm
e}}{\sqrt{p(p-1)}e^{-\Delta\mathcal{N}}})\,.$ (A.2)

First we consider $-k\tau_{\rm e}>\sqrt{p(p-1)}e^{-\Delta\mathcal{N}}$ for
which $\tilde{\tau}_{\rm c}$ happens after $\tau_{\rm c}$. As a result, for
these scales tachyonic instability is possible only if the expression in the
bracket in (3.21) becomes negative before ${\tilde{\tau}}_{\rm c}$. Namely,
the tachyonic instability in this case requires $-k\tilde{\tau}_{\rm
c}<\sqrt{2}$, which is rewritten as $-k\tau_{\rm
e}<2^{\frac{p+1}{2p}}[p(p-1)]^{-\frac{1}{2p}}e^{-\Delta\mathcal{N}}$ and thus
is inconsistent with the our primary condition $-k\tau_{\rm
e}>\sqrt{p(p-1)}e^{-\Delta\mathcal{N}}$ for $p>2$. Therefore there is no
chance to have any growing mode for $-k\tau_{\rm
e}>\sqrt{p(p-1)}e^{-\Delta\mathcal{N}}$ if $p>2$. For $p<2$, however, we have
growing modes for

$\displaystyle\sqrt{p(p-1)}e^{-\Delta\mathcal{N}}<-k\tau_{\rm
e}<2^{\frac{p+1}{2p}}[p(p-1)]^{-\frac{1}{2p}}e^{-\Delta\mathcal{N}}\,,$ (A.3)

Second, we consider the range $-k\tau_{\rm
e}<\sqrt{p(p-1)}e^{-\Delta\mathcal{N}}$ for which $\tilde{\tau}_{\rm c}$
happens before $\tau_{\rm c}$. For these scales the equation of the mode
function can be approximated as

$u_{k}^{\prime\prime}+\Big{(}k^{2}-\frac{2}{\tau^{2}}\Big{)}u_{k}=0\,,\qquad\tau_{\rm
i}<\tau<\tilde{\tau}_{\rm c}$ (A.4)
$u_{k}^{\prime\prime}+\Big{(}-\frac{p(p-1)}{\tau^{2}}\Big{)}u_{k}=0\,,\qquad\tilde{\tau}_{\rm
c}<\tau<\tau_{\rm c}$ (A.5)
$u_{k}^{\prime\prime}+\Big{(}\frac{R}{\tau^{2}}\left(\frac{\tau_{\rm
e}}{\tau}\right)^{2p}\Big{)}u_{k}=0\,,\qquad\tau_{\rm c}<\tau<\tau_{\rm e}\,.$
(A.6)

We thus find the following range of momenta for the growing modes

$\displaystyle\sqrt{2}\,e^{-{\cal N}}<-k\tau_{\rm
e}<\sqrt{p(p-1)}e^{-\Delta\mathcal{N}}\,.$ (A.7)

From Eqs. (A.3) and (A.7) we find the range of momenta for the all growing
modes as follows

$x_{\min}=\sqrt{2}\,e^{-{\cal N}}\,,\hskip
28.45274ptx_{\max}=\begin{cases}\sqrt{p(p-1)}e^{-\Delta\mathcal{N}}\hskip
56.9055ptp>2\\\
2^{\frac{p+1}{2p}}[p(p-1)]^{-\frac{1}{2p}}e^{-\Delta\mathcal{N}}\hskip
24.18501ptp<2\end{cases}\,,$ (A.8)

which we use to find the energy density of the produced longitudinal dark
photons. Here, the lower bound corresponds to the mode for which the
expression in the bracket in (A.4) changes its sign at the beginning of
inflation.

## References

  * [1] A. E. Nelson and J. Scholtz, _Dark Light, Dark Matter and the Misalignment Mechanism_ , _Phys. Rev. D_ 84 (2011) 103501 [1105.2812].
  * [2] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo and A. Ringwald, _WISPy Cold Dark Matter_ , _JCAP_ 06 (2012) 013 [1201.5902].
  * [3] P. W. Graham, J. Mardon and S. Rajendran, _Vector Dark Matter from Inflationary Fluctuations_ , _Phys. Rev. D_ 93 (2016) 103520 [1504.02102].
  * [4] P. Agrawal, G. Marques-Tavares and W. Xue, _Opening up the QCD axion window_ , _JHEP_ 03 (2018) 049 [1708.05008].
  * [5] R. T. Co, A. Pierce, Z. Zhang and Y. Zhao, _Dark Photon Dark Matter Produced by Axion Oscillations_ , _Phys. Rev. D_ 99 (2019) 075002 [1810.07196].
  * [6] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi and F. Takahashi, _Relic Abundance of Dark Photon Dark Matter_ , _Phys. Lett. B_ 801 (2020) 135136 [1810.07188].
  * [7] J. A. Dror, K. Harigaya and V. Narayan, _Parametric Resonance Production of Ultralight Vector Dark Matter_ , _Phys. Rev. D_ 99 (2019) 035036 [1810.07195].
  * [8] M. Goodsell, J. Jaeckel, J. Redondo and A. Ringwald, _Naturally Light Hidden Photons in LARGE Volume String Compactifications_ , _JHEP_ 11 (2009) 027 [0909.0515].
  * [9] J. Alexander et al., _Dark Sectors 2016 Workshop: Community Report_ , 8, 2016, 1608.08632.
  * [10] A. Pierce, K. Riles and Y. Zhao, _Searching for Dark Photon Dark Matter with Gravitational Wave Detectors_ , _Phys. Rev. Lett._ 121 (2018) 061102 [1801.10161].
  * [11] C. S. Machado, W. Ratzinger, P. Schwaller and B. A. Stefanek, _Audible Axions_ , _JHEP_ 01 (2019) 053 [1811.01950].
  * [12] C. S. Machado, W. Ratzinger, P. Schwaller and B. A. Stefanek, _Gravitational wave probes of axionlike particles_ , _Phys. Rev. D_ 102 (2020) 075033 [1912.01007].
  * [13] B. Salehian, M. A. Gorji, S. Mukohyama and H. Firouzjahi, _Analytic study of dark photon and gravitational wave production from axion_ , 2007.08148.
  * [14] Y. Michimura, T. Fujita, S. Morisaki, H. Nakatsuka and I. Obata, _Ultralight vector dark matter search with auxiliary length channels of gravitational wave detectors_ , _Phys. Rev. D_ 102 (2020) 102001 [2008.02482].
  * [15] R. Namba and M. Suzuki, _Implications of Gravitational-wave Production from Dark Photon Resonance to NANOGrav Observations and Effective number of Relativistic Species_ , 2009.13909.
  * [16] N. Kitajima, J. Soda and Y. Urakawa, _Nano-Hz gravitational wave signature from axion dark matter_ , 2010.10990.
  * [17] R. T. Co, L. J. Hall and K. Harigaya, _QCD Axion Dark Matter with a Small Decay Constant_ , _Phys. Rev. Lett._ 120 (2018) 211602 [1711.10486].
  * [18] E. W. Kolb and A. J. Long, _Completely Dark Photons from Gravitational Particle Production During Inflation_ , 2009.03828.
  * [19] M. M. Anber and L. Sorbo, _N-flationary magnetic fields_ , _JCAP_ 10 (2006) 018 [astro-ph/0606534].
  * [20] M.-a. Watanabe, S. Kanno and J. Soda, _Inflationary Universe with Anisotropic Hair_ , _Phys. Rev. Lett._ 102 (2009) 191302 [0902.2833].
  * [21] M. Bastero-Gil, J. Santiago, L. Ubaldi and R. Vega-Morales, _Vector dark matter production at the end of inflation_ , _JCAP_ 04 (2019) 015 [1810.07208].
  * [22] K. Nakayama, _Vector Coherent Oscillation Dark Matter_ , _JCAP_ 10 (2019) 019 [1907.06243].
  * [23] K. Nakayama, _Constraint on Vector Coherent Oscillation Dark Matter with Kinetic Function_ , _JCAP_ 08 (2020) 033 [2004.10036].
  * [24] Y. Nakai, R. Namba and Z. Wang, _Light Dark Photon Dark Matter from Inflation_ , 2004.10743.
  * [25] B. Salehian, M. A. Gorji, H. Firouzjahi and S. Mukohyama, _Vector dark matter production from inflation with symmetry breaking_ , 2010.04491.
  * [26] M.-a. Watanabe, S. Kanno and J. Soda, _The Nature of Primordial Fluctuations from Anisotropic Inflation_ , _Prog. Theor. Phys._ 123 (2010) 1041 [1003.0056].
  * [27] R. Emami, H. Firouzjahi, S. Sadegh Movahed and M. Zarei, _Anisotropic Inflation from Charged Scalar Fields_ , _JCAP_ 02 (2011) 005 [1010.5495].
  * [28] N. Bartolo, S. Matarrese, M. Peloso and A. Ricciardone, _Anisotropic power spectrum and bispectrum in the $f(\phi)F^{2}$ mechanism_, _Phys. Rev. D_ 87 (2013) 023504 [1210.3257].
  * [29] R. Emami, _Anisotropic Inflation and Cosmological Observations_ , Ph.D. thesis, Hong Kong U. Sci. Tech., 2015.  1511.01683.
  * [30] S. Shakeri, M. A. Gorji and H. Firouzjahi, _Schwinger Mechanism During Inflation_ , _Phys. Rev. D_ 99 (2019) 103525 [1903.05310].
  * [31] B. Ratra, _Cosmological ’seed’ magnetic field from inflation_ , _Astrophys. J. Lett._ 391 (1992) L1.
  * [32] W. Garretson, G. B. Field and S. M. Carroll, _Primordial magnetic fields from pseudoGoldstone bosons_ , _Phys. Rev. D_ 46 (1992) 5346 [hep-ph/9209238].
  * [33] J. Martin and J. Yokoyama, _Generation of Large-Scale Magnetic Fields in Single-Field Inflation_ , _JCAP_ 01 (2008) 025 [0711.4307].
  * [34] V. Demozzi, V. Mukhanov and H. Rubinstein, _Magnetic fields from inflation?_ , _JCAP_ 08 (2009) 025 [0907.1030].
  * [35] S. Kanno, J. Soda and M.-a. Watanabe, _Cosmological Magnetic Fields from Inflation and Backreaction_ , _JCAP_ 12 (2009) 009 [0908.3509].
  * [36] R. Emami, H. Firouzjahi and M. Movahed, _Inflation from Charged Scalar and Primordial Magnetic Fields?_ , _Phys. Rev. D_ 81 (2010) 083526 [0908.4161].
  * [37] T. Fujita and S. Mukohyama, _Universal upper limit on inflation energy scale from cosmic magnetic field_ , _JCAP_ 10 (2012) 034 [1205.5031].
  * [38] C. Caprini and L. Sorbo, _Adding helicity to inflationary magnetogenesis_ , _JCAP_ 10 (2014) 056 [1407.2809].
  * [39] C. Caprini, M. C. Guzzetti and L. Sorbo, _Inflationary magnetogenesis with added helicity: constraints from non-gaussianities_ , _Class. Quant. Grav._ 35 (2018) 124003 [1707.09750].
  * [40] J. Schober, T. Fujita and R. Durrer, _Generation of chiral asymmetry via helical magnetic fields_ , _Phys. Rev. D_ 101 (2020) 103028 [2002.09501].
  * [41] A. Talebian, A. Nassiri-Rad and H. Firouzjahi, _Revisiting Magnetogenesis during Inflation_ , 2007.11066.
  * [42] Planck collaboration, Y. Akrami et al., _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * [43] R. Emami and H. Firouzjahi, _Curvature Perturbations in Anisotropic Inflation with Symmetry Breaking_ , _JCAP_ 10 (2013) 041 [1301.1219].
  * [44] A. A. Abolhasani, R. Emami, J. T. Firouzjaee and H. Firouzjahi, _$\delta N$ formalism in anisotropic inflation and large anisotropic bispectrum and trispectrum_, _JCAP_ 08 (2013) 016 [1302.6986].
  * [45] L. Sorbo, _Parity violation in the Cosmic Microwave Background from a pseudoscalar inflaton_ , _JCAP_ 06 (2011) 003 [1101.1525].
  * [46] K. Yamamoto, M.-a. Watanabe and J. Soda, _Inflation with Multi-Vector-Hair: The Fate of Anisotropy_ , _Class. Quant. Grav._ 29 (2012) 145008 [1201.5309].
  * [47] K. Yamamoto, _Primordial Fluctuations from Inflation with a Triad of Background Gauge Fields_ , _Phys. Rev. D_ 85 (2012) 123504 [1203.1071].
  * [48] H. Firouzjahi, M. A. Gorji, S. A. Hosseini Mansoori, A. Karami and T. Rostami, _Charged Vector Inflation_ , _Phys. Rev. D_ 100 (2019) 043530 [1812.07464].
  * [49] M. A. Gorji, S. A. Hosseini Mansoori and H. Firouzjahi, _Inflation with multiple vector fields and non-Gaussianities_ , 2008.08195.