# Argument for the radiation-dominated behavior of matter fields in the
preinflationary era

Taylor M. Ordines  Eric D. Carlson  Department of Physics, Wake Forest
University, 1834 Wake Forest Road, Winston-Salem, North Carolina 27109, USA

###### Abstract

We investigate the leading-order behavior of matter fields in the
preinflationary era using the semiclassical approximation. Many inflationary
models assume without supporting arguments that the Universe was radiation
dominated prior to inflation, leading to modifications of cosmological
observables, such as the Cosmic Microwave Background power spectrum. In
previous work, we demonstrated that conformally coupled scalar fields do have
a radiation-like contribution to the stress-energy tensor at sufficiently
early times. In this work, we extend these arguments to apply to massless
spin-1 fields and massive or massless spin-$\frac{1}{2}$ fields. We find
massless spin-1 fields always have a radiation-like contribution. For
spin-$\frac{1}{2}$ fields, we find the contribution at early times is
radiation-like assuming this is the dominant contribution to the stress-energy
tensor.

##  I Introduction

The inflationary paradigm explains many characteristics of our observed
Universe. Many models of inflation make use of a radiation-dominated behavior
of matter fields in the preinflationary era in order to obtain modifications
to the standard predictions of inflation and better explain observed phenomena
(see [1] and references therein for a review of models using a radiation-
dominated preinflationary era). For instance, various models involving the
transition from a radiation-dominated era to inflation lead to modifications
of the Cosmic Microwave Background anisotropy spectrum, such as a lowering of
the quadrupole moment which appears to be anomalously suppressed [2]. The
assumption that fields are radiation-dominated prior to inflation is central
to such models, though it is non-trivial to demonstrate that such behavior was
the case for our Universe, and it would be hopeless to try to detect particles
due to these fields today given the effects of inflation. It is therefore
interesting to analyze the preinflationary era and the behavior of matter
fields in it.

The primary objective of this work is to investigate the behavior of quantum
fields in the preinflationary era. A complete analysis of the the early
universe would require a theory of quantum gravity, which for now is out of
reach, but one could anticipate that after the Planck era some span of the
preinflationary era would have curvature well below the Planck scale, in which
case the semiclassical approximation should hold. In this paper, we will
assume the semiclassical approximation to hold in some portion of the
preinflationary era, and we will use this framework to analyze quantum fields
in a Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime and obtain their
corresponding contributions to the energy density.

The case of massive scalar fields was investigated using the semiclassical
approximation in [3]. There it was shown that under a set of conditions on the
scale factor the fields did result in radiation-dominated contributions to the
energy density. The conditions on the scale factor were also argued to likely
be valid for our Universe. In the present work, we will investigate the
corresponding behavior of other fields, namely massless spin-1 and massive or
massless spin-$\frac{1}{2}$ fields.

The analysis for spin-1 and spin-$\frac{1}{2}$ fields is generally more
complicated than for scalar fields. For spin-$\frac{1}{2}$ fields, the mass
appears explicity in the expression for the counterterms at higher than zeroth
adiabatic order, which leads to additional complications. Performing adiabatic
regularization with a zeroth-order parameterization of the states does not
explicitly appear to produce finite energy densities until one considers the
higher order contributions buried in the parameterization. For an overview of
the adiabatic regularization procedure for spin-$\frac{1}{2}$, see [4, 5, 6].
For spin-1 fields one must consider the massive and massless cases separately.
Massive vector fields do appear in the Standard Model, but they are ultimately
due to interactions with the Higgs field and are out of the scope of this
paper. For the massless case one can show that the analysis decomposes into
four decoupled copies of a scalar field, as was done in Ref. [7]. There the
analysis was performed to obtain the trace anomaly, but obtaining the
renormalized energy density from this groundwork is non-trivial.

In the following, we work in Planck units with $c=G=\hbar=1$ and use the
$(-,+,+,+)$ signature for the metric. We assume a spatially flat universe
described by the FLRW metric

$ds^{2}=a^{2}\left(-d\eta^{2}+d\vec{x}^{2}\right)\;,$ (1)

where $a\equiv a(\eta)$ is the scale factor and $\eta$ is conformal time. We
will employ a prime, such as $a^{\prime}\equiv\partial_{\eta}a$, to denote
differentiation with respect to conformal time. We are interested in taking
$a\rightarrow 0$ at early times and determining whether one can assume
radiation-dominated behavior. However, we do not want to consider times in the
Planck era, during which the semiclassical approximation is not assumed to be
valied, so we denote by $\eta_{0}$ the earliest conformal time of
consideration and implicitly assume it to be past the Planck era. In Planck
units, this corresponds to a value of the Hubble parameter
$H\equiv\frac{a^{\prime}}{a^{2}}\ll 1$. We therefore aim to determine if the
renormalized energy density $\rho_{r}$ for each field is appropriately
radiation dominated during a range of conformal time $\eta_{0}<\eta<\eta_{1}$.
We will allow for the possibility of a non-zero mass $m$ for fermions, in
which case we will insist that $m\ll 1$ in Planck units. If there are fermion
fields with $m\gtrsim 1$, these will have mass comparable to or greater than
the Planck mass, and then we will assume the fields have no contribution to
the total energy density when $\eta\gtrsim\eta_{0}$.

The body of this article is split into three sections, one for each of the two
spins of fields and one for a discussion of the results obtained. Section II
contains the analysis for spin-1 vector fields, beginning with a summary of
the adiabatic regularization procedure and concluding with the renormalization
of the energy density and our main result for the spin-1 case. Section III
contains the analysis for spin-$\frac{1}{2}$ fields and begins with a summary
of the modified adiabatic regularization procedure, as described in [6]. The
subsequent renormalization and analysis is more complicated than that for the
vector fields and is split into subsections: in Sec. III.1 we give the
renormalization counterterms and preview the assumptions built into our
analysis; in Sec. III.2, we derive bounds on the renormalized energy density
by splitting the contributions into high and low energy terms and analyzing
each in turn; and in Sec. III.3 we obtain the leading order behavior of the
renormalized energy density and our main result for the spin-$\frac{1}{2}$
case. We close in Sec. IV with a discussion of our main results and final
remarks.

##  II Energy Density for Vector Fields

As discussed in the introduction, massive and massless vector fields must be
treated separately. In the Lagrangian description, the massive field is
described by the Proca action [8], which in flat spacetime is a gauge-fixed
theory involving the Higgs mechanism. Working with interactions and the curved
spacetime form of the action is beyond the scope of this paper. We will focus
instead on the massless case.

According to an argument in [9], any conformally invariant theory in a flat
FLRW spacetime will have a stress-energy tensor that contains two terms, one
of which is radiation-like and the other of which is the anomalous term.
Renormalizing the electromagnatic field, however, requires that masses are
introduced for the photon and ghost fields, which break the conformal
invariance. These masses are then taken to zero to obtain physical results. We
therefore feel it is worth working through the details of this procedure to
confirm that the argument in [9] works.

The massless vector field in curved spacetime is given by the massless limit
of the theory described by Lagrangian [7]

$\mathcal{L}=\sqrt{-g}\left[-\frac{1}{4}g^{\mu\rho}g^{\nu\sigma}F_{\mu\nu}F_{\rho\sigma}-\frac{1}{2\xi}\nabla^{\mu}A_{\mu}\nabla^{\nu}A_{\nu}-\frac{1}{2}m^{2}g^{\mu\nu}A_{\mu}A_{\nu}+ig^{\mu\nu}\partial_{\mu}\bar{\chi}\partial_{\nu}\chi+im_{\chi}\bar{\chi}\chi\right]\;,$
(2)

where $A_{\mu}$ is the four-vector, $\xi$ is the gauge fixing parameter,
$\chi$ is the (complex) ghost field used to maintain gauge invariance,
$m_{\chi}$ is the mass of the ghost field, and
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$. Masses are included
in the Lagrangian (2) in order to properly renormalize the theory. These
masses can trivially be taken to zero before computing any unrenormalized
observables, but they are crucial in obtaining the appropriate renormalization
counterterms. We will make clear when they are to be taken to zero in the
following. The energy-momentum tensor from Eq. (2) is

$\displaystyle T_{\mu\nu}=$ $\displaystyle
T_{\mu\nu}^{\text{Maxwell}}+T_{\mu\nu}^{\text{G}}+T_{\mu\nu}^{\text{ghost}}\;,$
(3a) $\displaystyle T_{\mu\nu}^{\text{Maxwell}}\equiv$
$\displaystyle-\frac{1}{4}g_{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}F_{\alpha\rho}F_{\beta\sigma}+g^{\alpha\beta}F_{\alpha\mu}F_{\beta\nu}-\frac{1}{2}g_{\mu\nu}m^{2}g^{\alpha\beta}A_{\alpha}A_{\beta}+m^{2}A_{\mu}A_{\nu}\;,$
(3b) $\displaystyle T_{\mu\nu}^{\text{G}}\equiv$
$\displaystyle\frac{1}{\xi}\bigg{[}-\frac{1}{2}g_{\mu\nu}\left(g^{\alpha\beta}\nabla_{\alpha}A_{\beta}\right)^{2}+\left(g_{\mu\nu}g^{\rho\sigma}A_{\sigma}\nabla_{\rho}-A_{\nu}\nabla_{\mu}-A_{\mu}\nabla_{\nu}\right)\left(g^{\alpha\beta}\nabla_{\alpha}A_{\beta}\right)\bigg{]}\;,$
(3c) $\displaystyle T_{\mu\nu}^{\text{ghost}}\equiv$ $\displaystyle
ig_{\mu\nu}g^{\rho\sigma}\partial_{\rho}\bar{\chi}\partial_{\sigma}\chi-i\left(\partial_{\mu}\bar{\chi}\partial_{\nu}\chi+\partial_{\nu}\bar{\chi}\partial_{\mu}\chi\right)+ig_{\mu\nu}m_{\chi}^{2}\bar{\chi}\chi\;.$
(3d)

Following the procedure in Ref. [7], one can define the components of the
four-vector $A_{\mu}$ as a combination of temporal, transverse, and
longitudinal parts,

$A_{\mu}\equiv\left(A_{0},B_{i}+\partial_{i}A\right)\;.$ (4)

These components and the ghost field can be expanded in terms of mode
functions $Y_{a}$, for $a=0,L,T,\chi$, where $L$ represents longitudinal and
$T$ transverse contributions:

$A_{0}=\frac{1}{m^{2}a^{2}}\int\frac{d^{3}\vec{k}}{\left(2\pi\right)^{3}}\left[a_{\vec{k}}^{(0)}\left(\partial_{\eta}-\frac{2a^{\prime}}{a}\right)\left(maY_{0}\right)e^{i\vec{k}\cdot\vec{x}}+a_{\vec{k}}^{(3)}kmaY_{L}e^{i\vec{k}\cdot\vec{x}}+\text{H.c.}\right]\;,$
(5a)
$B_{i}=\int\frac{d^{3}\vec{k}}{\left(2\pi\right)^{3}}\sum_{p=1,2}\left(\epsilon_{i}^{p}a_{\vec{k}}^{(p)}Y_{T}e^{i\vec{k}\cdot\vec{x}}+\text{H.c.}\right)\;,$
(5b)
$A=\frac{1}{m^{2}a^{2}}\int\frac{d^{3}\vec{k}}{\left(2\pi\right)^{3}}\left(a_{\vec{k}}^{(0)}maY_{0}e^{i\vec{k}\cdot\vec{x}}-a_{\vec{k}}^{(3)}\partial_{\eta}\left(\frac{ma}{k}Y_{L}\right)e^{i\vec{k}\cdot\vec{x}}+\text{H.c.}\right)\;,$
(5c)
$\chi=\int\frac{d^{3}\vec{k}}{\left(2\pi\right)^{3}}\left(\frac{b_{\vec{k}}Y_{\chi}}{a}e^{i\vec{k}\cdot\vec{x}}+\frac{b^{\dagger}_{\vec{k}}Y_{\chi}^{*}}{a}e^{-i\vec{k}\cdot\vec{x}}\right)\;,$
(5d)
$\bar{\chi}=\int\frac{d^{3}\vec{k}}{\left(2\pi\right)^{3}}\left(\frac{\bar{b}_{\vec{k}}Y_{\chi}}{a}e^{i\vec{k}\cdot\vec{x}}+\frac{\bar{b}^{\dagger}_{\vec{k}}Y_{\chi}^{*}}{a}e^{-i\vec{k}\cdot\vec{x}}\right)\;,$
(5e)

where $a_{\vec{k}}^{(\mu)}$, $b_{\vec{k}}$, and $\bar{b}_{\vec{k}}$ and their
hermitian conjugates are annihilation and creation operators,
$\epsilon_{i}^{p}$ are the two polarization vectors of the transverse modes,
and H.c. represents the hermitian conjugate of all preceding terms. The
creation and annihilation operators satisfy the commutation and
anticommutation relations

$\displaystyle\left[a_{\vec{k}}^{(a)},a_{\vec{k}^{\prime}}^{(b)\dagger}\right]$
$\displaystyle=\eta^{ab}\left(2\pi\right)^{3}\delta^{(3)}(\vec{k}-\vec{k}^{\prime})\;,$
(6a)
$\displaystyle\left\\{b_{\vec{k}},\bar{b}_{\vec{k}^{\prime}}^{\dagger}\right\\}$
$\displaystyle=-\left\\{\bar{b}_{\vec{k}},b_{\vec{k}^{\prime}}^{\dagger}\right\\}=i\left(2\pi\right)^{3}\delta^{(3)}(\vec{k}-\vec{k}^{\prime})\;,$
(6b)

where $\eta^{ab}=\text{diag}(-1,1,1,1)$. The polarization vectors satisfy

$\sum_{i}k_{i}\epsilon_{i}^{p}=0\;,$ (7a)
$\sum_{i}\epsilon_{i}^{p}\epsilon_{i}^{p^{\prime}}=\delta^{pp^{\prime}}\;,$
(7b)
$\sum_{p=1,2}\epsilon_{i}^{p}\epsilon_{j}^{p}=\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\;,$
(7c)

with $p$ the polarization index. The modes $Y_{a}$ satisfy the decoupled set
of differential equation

$\left(\partial_{\eta}^{2}+\Omega_{a}^{2}\right)Y_{a}=0\;,$ (8)

where

$\Omega_{a}^{2}\equiv\omega_{a}^{2}+\zeta_{a}$ (9)

is defined for each contribution by

$\displaystyle\omega_{0}^{2}$ $\displaystyle\equiv k^{2}+\xi m^{2}a^{2}\;,$
(10a) $\displaystyle\omega_{L,T}^{2}$ $\displaystyle\equiv
k^{2}+m^{2}a^{2}\;,$ (10b) $\displaystyle\omega_{\chi}^{2}$
$\displaystyle\equiv k^{2}+m_{\chi}^{2}a^{2}\;,$ (10c)

and

$\displaystyle\zeta_{0,\chi}$
$\displaystyle\equiv-\frac{a^{\prime\prime}}{a}\;,$ (11a)
$\displaystyle\zeta_{L}$
$\displaystyle\equiv\frac{a^{\prime\prime}}{a}-\frac{2a^{\prime 2}}{a^{2}}\;,$
(11b) $\displaystyle\zeta_{T}$ $\displaystyle\equiv 0\;.$ (11c)

The modes satisfy the standard normalization conditions for a scalar field,

$Y_{a}Y_{a}^{*\prime}-Y_{a}^{\prime}Y_{a}^{*}=i\;.$ (12)

Note that in the massless limit $Y_{0}$ and $Y_{\chi}$ satisfy the same
differential equations, so we can choose

$Y_{0}=Y_{\chi}\quad\text{for}\quad m=m_{\chi}=0\;.$ (13)

Using this mode decomposition, one can find contributions to the energy
density $\rho\equiv\langle T_{\hat{0}\hat{0}}\rangle$ from Eqs. (3b)–(3d) in
terms of $Y_{a}$:

$\displaystyle\rho^{\text{Maxwell}}=$ $\displaystyle\lim\limits_{m\to
0}\frac{1}{4\pi^{2}a^{4}}\int_{0}^{\Lambda}dk\,k^{2}\bigg{[}\omega_{T}^{2}\left|Y_{T}\right|^{2}+\left|Y_{T}^{\prime}\right|^{2}-\left(k^{2}+\frac{a^{\prime
2}}{a^{2}}\right)\left|Y_{0}\right|^{2}+\frac{a^{\prime}}{a}\partial_{\eta}\left(\left|Y_{0}\right|^{2}\right)-\left|Y_{0}^{\prime}\right|^{2}$
$\displaystyle+\left(\omega_{L}^{2}+\frac{a^{\prime
2}}{a^{2}}\right)\left|Y_{L}\right|^{2}+\frac{a^{\prime}}{a}\partial_{\eta}\left(\left|Y_{L}\right|^{2}\right)+\left|Y_{L}^{\prime}\right|^{2}\bigg{]}$
(14a) $\displaystyle\rho^{\text{G}}=$ $\displaystyle\lim\limits_{m\to
0}\frac{1}{2\pi^{2}a^{4}}\int_{0}^{\Lambda}dk\,k^{2}\left[\left(\omega_{0}^{2}-\frac{1}{2}\xi
m^{2}a^{2}+\frac{a^{\prime
2}}{a^{2}}\right)\left|Y_{0}\right|^{2}-\frac{a^{\prime}}{a}\partial_{\eta}\left(\left|Y_{0}\right|^{2}\right)+\left|Y_{0}^{\prime}\right|^{2}\right]$
(14b) $\displaystyle\rho^{\text{ghost}}=$
$\displaystyle\lim\limits_{m_{\chi}\to
0}\frac{1}{2\pi^{2}a^{4}}\int_{0}^{\Lambda}dk\,k^{2}\left[-\left(\omega_{\chi}^{2}+\frac{a^{\prime
2}}{a^{2}}\right)\left|Y_{\chi}\right|^{2}+\frac{a^{\prime}}{a}\partial_{\eta}\left(\left|Y_{\chi}\right|^{2}\right)-\left|Y_{\chi}^{\prime}\right|^{2}\right]\;,$
(14c)

where $\Lambda$ is a cutoff regulator which we will later demonstrate can be
taken to infinity in the massless limit.

One may adiabatically renormalize the energy density by writing the vacuum
states with the standard WKB ansatz

$\displaystyle
Y_{a}=\frac{1}{\sqrt{2W_{a}}}e^{-i\int^{\eta}d\bar{\eta}W_{a}(\bar{\eta})}\;,$
(15)

where

$\left(W_{a}\right)^{2}=\Omega_{a}^{2}-\left[\frac{W_{a}^{\prime\prime}}{2W_{a}}-\frac{3}{4}\left(\frac{W_{a}^{\prime}}{W_{a}}\right)^{2}\right]\;.$
(16)

Solutions to Eq. (16) can be approximated using $W_{a}^{(0)}=\Omega_{a}$ as
the lowest order and iterating to higher orders, keeping to the appropriate
adiabatic order, given by the number of time derivatives on the scale factor,
at each iteration. Substitution of $W_{a}$ to some adiabatic order $A$ into
Eq. (15) would then require expanding the square root only to terms of
adiabatic order $A$.

In order to renormalize the energy density, one must take Eq. (15) to the
appropriate adiabatic order and substitute into Eqs. (14a)–(14c) to produce
the renormalization counterterms. On dimensional grounds, one would need to
keep to fourth adiabatic order to renormalize $\rho$. The fourth order
counterterms are

$\displaystyle\rho_{c}^{(4)}=$ $\displaystyle\lim\limits_{m\to
0}\frac{1}{4a^{2}}\int\frac{d^{3}k}{(2\pi)^{3}}\bigg{\\{}W_{0}+W_{L}+2W_{T}-2W_{\chi}+\frac{\omega_{0}^{2}}{W_{0}}+\frac{\omega_{L}^{2}}{W_{L}}+\frac{2\omega_{T}^{2}}{W_{T}}-\frac{2\omega_{\chi}^{2}}{W_{\chi}}+\frac{a^{\prime
2}}{a^{2}W_{0}}$ $\displaystyle+\frac{a^{\prime
2}}{a^{2}W_{L}}-\frac{2a^{\prime
2}}{a^{2}W_{\chi}}+\frac{a^{\prime}W_{0}^{\prime}}{aW_{0}^{2}}-\frac{a^{\prime}W_{L}^{\prime}}{aW_{L}^{2}}-\frac{2a^{\prime}W_{\chi}^{\prime}}{aW_{\chi}^{2}}+\frac{W_{0}^{\prime
2}}{4W_{0}^{3}}+\frac{W_{L}^{\prime 2}}{4W_{L}^{3}}+\frac{W_{T}^{\prime
2}}{2W_{T}^{3}}-\frac{W_{\chi}^{\prime 2}}{2W_{\chi}^{3}}\bigg{\\}}^{(4)}\;,$
(17)

where $\\{\dots\\}^{(4)}$ implies that all $W_{a}$ are taken to fourth order.

In order to analyze the early-time behavior of the energy density, we will
parameterize the mode functions in terms of zeroth-order adiabatic states

$Y_{a}=\alpha_{k,a}Y_{a}^{(0)}+\beta_{k,a}Y_{a}^{(0)*}\;,$ (18)
$Y_{a}^{\prime}=\alpha_{k,a}Y_{a}^{(0)\prime}+\beta_{k,a}Y_{a}^{(0)*\prime}\;.$
(19)

From the normalization condition (12), one has

$\left|\alpha_{k,a}\right|^{2}-\left|\beta_{k,a}\right|^{2}=1\;.$ (20)

One could instead parameterize the mode functions in terms of higher-order
adiabatic states, in which case the Bogoliubov coefficients $\alpha_{k,a}$ and
$\beta_{k,a}$ would be constant to the given adiabatic order, but zeroth order
will be sufficient to properly renormalize the theory. Substituting Eqs. (18)
and (19) into Eq. (8), one obtains differential equations for the
coefficients,

$\displaystyle\alpha_{k,a}^{\prime}$
$\displaystyle=\frac{\Omega_{a}^{\prime}}{2\Omega_{a}}\beta_{k,a}e^{2i\theta_{a}}$
(21a) $\displaystyle\beta_{k,a}^{\prime}$
$\displaystyle=\frac{\Omega_{a}^{\prime}}{2\Omega_{a}}\alpha_{k,a}e^{-2i\theta_{a}}\;,$
(21b)

where $\theta_{a}\equiv\int^{\eta}d\bar{\eta}W_{a}(\bar{\eta})$.

One can then obtain the renormalized energy density in terms of $\alpha_{k,a}$
and $\beta_{k,a}$ by subtracting these renormalization counterterms from the
unrenormalized energy density given in Eqs. (14a)–(14c). Using a zeroth order
parameterization, one finds that the renormalized energy density $\rho_{r}$
separates into analytic and mode terms

$\rho_{r}=\rho_{\text{an}}+\rho_{\alpha\beta}\;,$ (22)

where the analytic terms $\rho_{\text{an}}$ are finite higher-order terms,
independent of the cutoff regulator $\Lambda$, coming from the counterterms
(17),

$\displaystyle\rho_{\text{an}}=$
$\displaystyle\frac{1}{2880\pi^{2}}\left[62\,{}^{(3)}H_{00}+\left(3+\frac{5}{2}\ln\xi\right)\,{}^{(1)}H_{00}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{2880\pi^{2}}\bigg{[}\frac{186a^{\prime
4}}{a^{6}}+\frac{216a^{\prime
2}a^{\prime\prime}}{a^{5}}+\frac{54a^{\prime\prime
2}}{a^{4}}-\frac{108a^{\prime}a^{\prime\prime\prime}}{a^{4}}$
$\displaystyle+\ln\xi\left(\frac{180a^{\prime
2}a^{\prime\prime}}{a^{5}}+\frac{45a^{\prime\prime
2}}{a^{4}}-\frac{90a^{\prime}a^{\prime\prime\prime}}{a^{4}}\right)\bigg{]}\;,$
(23)

where ${}^{(1)}H_{00}$ and ${}^{(3)}H_{00}$ are higher order corrections to
the Einstein field equations [10], and the mode terms $\rho_{\alpha\beta}$ are
those coming from the zeroth order parameterization of Eq. (18),

$\displaystyle\rho_{\alpha\beta}\equiv$
$\displaystyle\lim\limits_{m,m_{\chi}\to
0}\frac{1}{a^{2}}\int_{0}^{\Lambda}\frac{d^{3}k}{(2\pi)^{3}}\bigg{[}\frac{k^{2}+\omega_{T}^{2}}{a^{2}\omega_{T}}\left|\beta_{k,T}\right|^{2}+\frac{m^{2}}{2\omega_{L}}\left|\beta_{k,L}\right|^{2}+\frac{1}{\omega_{0}}\left(\frac{3\omega_{0}^{2}+k^{2}}{2a^{2}}+\frac{a^{\prime
2}}{a^{4}}\right)\left|\beta_{k,0}\right|^{2}$
$\displaystyle-\frac{1}{\omega_{\chi}}\left(\frac{\omega_{\chi}^{2}+k^{2}}{a^{2}}+\frac{a^{\prime
2}}{a^{4}}\right)\left|\beta_{k,\chi}\right|^{2}-\frac{m^{2}}{\omega_{T}}\text{Re}\\!\left(\alpha_{k,T}\beta_{k,T}^{*}e^{-2i\theta_{T}}\right)+\frac{m^{2}}{2\omega_{L}}\text{Re}\\!\left(\alpha_{k,L}\beta_{k,L}^{*}e^{-2i\theta_{L}}\right)$
$\displaystyle+\frac{1}{\omega_{0}}\left(\frac{a^{\prime 4}}{a^{2}}-\frac{\xi
m^{2}}{2a^{2}}\right)\text{Re}\\!\left(\alpha_{k,0}\beta_{k,0}^{*}e^{-2i\theta_{0}}\right)+\frac{1}{\omega_{\chi}}\left(m_{\chi}^{2}-\frac{a^{\prime
2}}{a^{4}}\right)\text{Re}\\!\left(\alpha_{k,\chi}\beta_{k,\chi}^{*}e^{-2i\theta_{\chi}}\right)$
$\displaystyle-\frac{2a^{\prime}}{a^{3}}\text{Im}\\!\left(\alpha_{k,0}\beta_{k,0}^{*}e^{-2i\theta_{0}}\right)+\frac{2a^{\prime}}{a^{3}}\text{Im}\\!\left(\alpha_{k,\chi}\beta_{k,\chi}^{*}e^{-2i\theta_{\chi}}\right)\bigg{]}\;.$
(24)

Assuming $\beta_{k,a}$ falls faster than $k^{-2}$, integrating the terms in
Eq. (24) will yield finite results, even if $\Lambda\to\infty$, and hence the
massless limit can be freely taken inside the integral. One finds from
substitution of Eqs. (15), (18), and (19) into Eq. (8) that
$\alpha_{k,0}^{\prime}=\alpha_{k,\chi}^{\prime}$ and
$\beta_{k,0}^{\prime}=\beta_{k,\chi}^{\prime}$ after taking the massless
limit, which when combined with Eq. (13) allows one to choose the coefficients
for the $0$ and $\chi$ contributions to be identical. The mode term
contribution to the energy density therefore drastically simplifies to

$\rho_{\alpha\beta}=\frac{1}{\pi^{2}a^{4}}\int_{0}^{\Lambda}dk\,k^{3}\left|\beta_{k,T}\right|^{2}\;.$
(25)

The ultraviolet cutoff can now be removed, and we take the limit
$\Lambda\to\infty$. The renormalized energy density ultimately only depends on
the transverse mode functions, which are the only physical modes of the
theory, and higher order dependencies on the background curvature:

$\displaystyle\rho_{r}=$
$\displaystyle\frac{1}{\pi^{2}a^{4}}\int_{0}^{\infty}dk\;k^{3}\left|\beta_{k,T}\right|^{2}+\frac{1}{2880\pi^{2}}\left[62\,{}^{(3)}H_{00}+\left(3+\frac{5}{2}\ln\xi\right)\,{}^{(1)}H_{00}\right]\;.$
(26)

Assuming the higher order corrections are subdominant in the semiclassical
approximation below the Planck scale, $\rho_{r}$ for massless vector fields
does have the expected radiation-dominated behavior. This result agrees with
the prediction in [9] that other than the anomalous term all of the
contributions in a flat FLRW metric of a conformally invariant field will act
like radiation.

The higher-order terms in (26) are of the same form as those found for the
trace anomaly in [11, 7], in which the ${}^{(1)}H_{00}$ term has a gauge-
dependent coefficient. This coefficient corresponds to a $\Box R$ term
appearing in the trace anomaly. The exact value of this coefficient is
dependent on the regularization scheme used, unlike for scalar and
spin-$\frac{1}{2}$ fields which respectively have the same coefficient
regardless of regularization scheme.

##  III Energy Density for Dirac Fields

We now turn our attention to Dirac spinor fields in the preinflationary era of
an FLRW universe. We follow the work and notation given in Ref. [6]. There,
expressions are in terms of cosmic time $t$, related to conformal time by
$a\,d\eta=dt$. We summarize the procedure for obtaining the unrenormalized
energy density here.

Consider Dirac spinor fields $\Psi(x)$ that obey the Dirac equation in curved
spacetime,

$\left(i\gamma^{a}e\indices{{}_{a}^{\mu}}\nabla_{\mu}-m\right)\Psi=0\;,$ (27)

where $e\indices{{}_{a}^{\mu}}$ is the vierbein, $\gamma^{a}$ are the flat
spacetime Dirac matrices satisfying $\\{\gamma^{a},\gamma^{b}\\}=2\eta^{ab}$,
and $\nabla_{\mu}\equiv\partial_{\mu}+\Gamma_{\mu}$ is the covariant
derivative associated with the spin connection $\Gamma_{\mu}$. For the metric
(1), the Dirac equation (27) becomes

$\left[\gamma^{0}\left(\partial_{\eta}+\frac{3a^{\prime}}{2a}\right)+\gamma^{i}\partial_{i}+ima\right]\Psi=0\;.$
(28)

The field can be written in terms of creation and annihilation operators
$D^{\dagger}_{\vec{k}\lambda}(\eta)$ and $B_{\vec{k}\lambda}(\eta)$ as

$\Psi=\sum_{\lambda=\pm\frac{1}{2}}\int
d^{3}k\left(B_{\vec{k}\lambda}\psi_{\vec{k}\lambda}+D^{\dagger}_{\vec{k}\lambda}C\bar{\psi}_{\vec{k}\lambda}^{\intercal}\right)\;,$
(29)

where $C$ is the charge conjugation matrix, $\lambda=\pm 1/2$ represents the
helicity eigenvalue, and

$\left\\{B_{\vec{k},\lambda},B_{\vec{k},\lambda}^{\dagger}\right\\}=\delta_{\lambda\lambda^{\prime}}\delta^{(3)}(\vec{k}-\vec{k}^{\prime})$
(30)

and similarly for $D_{\vec{k}\lambda}$ and $D_{\vec{k}\lambda}^{\dagger}$,
with all other anticommutators vanishing. Working in the Dirac-Pauli
representation for $\gamma^{a}$, the modes $\psi_{\vec{k}\lambda}$ can be
written as

$\psi_{\vec{k}\lambda}(\eta,\vec{x})=\frac{e^{i\vec{k}\cdot\vec{x}}}{\sqrt{8\pi^{3}a^{3}}}\begin{pmatrix}h^{I}_{k}(\eta)\xi_{\lambda}(\vec{k})\;,\\\
h^{II}_{k}(\eta)\hat{k}\cdot\vec{\sigma}\xi_{\lambda}(\vec{k})\end{pmatrix}\;,$
(31)

where $\xi_{\lambda}(\vec{k})$ are two-component spinors and are eigenvectors
of the spin component along the $\vec{k}$ direction, so that
$\frac{1}{2}(\hat{k}\cdot\vec{\sigma})\xi_{\lambda}(\vec{k})=\lambda\xi_{\lambda}(\vec{k})$,
with normalization $\xi^{\dagger}_{\lambda}\xi_{\lambda}=1$, and $h^{I}_{k}$
and $h^{II}_{k}$ are scalar functions that satisfy coupled first-order
differential equations

$\displaystyle\partial_{\eta}h^{I}_{k}(\eta)$
$\displaystyle=-ikh^{II}_{k}(\eta)-ima(\eta)h^{I}_{k}(\eta)\;,$ (32a)
$\displaystyle\partial_{\eta}h^{II}_{k}(\eta)$
$\displaystyle=-ikh^{I}_{k}(\eta)+ima(\eta)h^{II}_{k}(\eta)\;,$ (32b)

and have the normalization

$\left|h^{I}_{k}(\eta)\right|^{2}+\left|h^{II}_{k}(\eta)\right|^{2}=1\;.$ (33)

The energy density for the Dirac field in terms of the mode functions
$h_{k}^{I,II}$ can be written as

$\rho=\frac{1}{\pi^{2}a^{4}}\int_{0}^{\infty}dkk^{2}\left[ma\left(\left|h^{II}_{k}\right|^{2}-\left|h^{I}_{k}\right|^{2}\right)-k\left(h^{I}_{k}h^{II*}_{k}+h^{I*}_{k}h^{II}_{k}\right)\right]\;.$
(34)

At this point, we will diverge from this procedure coming from [6], who
themselves proceeded to obtain counterterms to a generic unrenormalized energy
in an FLRW universe. These counterterms were then used to prove conservation
of the energy density and were applied to a de Sitter spacetime and a
radiation-dominated universe. As demonstrated there, one can make use of
adiabatic regularization to renormalize the energy density, but the ansatz
used for the WKB approximation of adiabatic states must be a modified form of
the standard ansatz. This renormalization procedure was first demonstrated in
[4] and has been applied in other cases [5, 12, 13, 14]. Here, we will use the
energy density (34) as well as the modified WKB ansatz to obtain
renormalization counterterms, but we will use a different process to obtain an
explicit form of the unrenormalized energy density. Namely, we will expand the
mode functions $h_{k}^{I,II}$ in terms of adiabatic modes $g_{k}^{I,II}$ via a
Bogoliubov-like expansion in order to analyze early time behavior using the
Bogoliubov coefficients.

First, one expands $h_{k}^{I,II}$ in terms of adiabatic modes $g^{I,II}_{k}$,

$\displaystyle h^{I}_{k}=\alpha_{k}g^{I}_{k}-\beta_{k}g_{k}^{II*}\;,$ (35a)
$\displaystyle h^{II}_{k}=\alpha_{k}g^{II}_{k}+\beta_{k}g_{k}^{I*}\;,$ (35b)

where $g^{I,II}_{k}$ are given to adiabatic order $A$,

$g^{I,II}_{k}=g^{I,II(0)}_{k}+g^{I,II(1)}_{k}+\dots+g^{I,II(A)}_{k}\;,$ (36)

with adiabatic order understood to be the number of conformal time
derivatives, and satisfy the differential equations (32a) and (32b), and
$\alpha_{k}$ and $\beta_{k}$ are the time-dependent Bogoliubov coefficients
and are constant to order $A$. Coupled first order differential equations for
$\alpha_{k}$ and $\beta_{k}$ can be obtained from Eqs. (32a), (32b), and (33).
The unrenormalized energy density in terms of the adiabatic states is

$\displaystyle\rho=$
$\displaystyle\frac{1}{\pi^{2}a^{4}}\int_{0}^{\infty}dkk^{2}\bigg{\\{}2\left|\beta_{k}\right|^{2}\left[ma\left(\left|g^{I}_{k}\right|^{2}-\left|g^{II}_{k}\right|^{2}\right)+2k\text{Re}\left(g^{I}_{k}g^{II*}_{k}\right)\right]+4ma\text{Re}\left(\alpha_{k}\beta_{k}^{*}g^{I}_{k}g^{II}_{k}\right)$
$\displaystyle+2k\text{Re}\left[\alpha_{k}\beta_{k}^{*}\left(\left(g^{II}_{k}\right)^{2}-\left(g^{I}_{k}\right)^{2}\right)\right]+ma\left(\left|g^{II}_{k}\right|^{2}-\left|g^{I}_{k}\right|^{2}\right)-2k\text{Re}\left(g^{I}_{k}g^{II*}_{k}\right)\bigg{\\}}\;,$
(37)

and the adiabatic renormalization counterterms are

$\displaystyle\rho_{c}=$
$\displaystyle\frac{1}{\pi^{2}a^{4}}\int_{0}^{\infty}dkk^{2}\left[ma\left(\left|g^{II}_{k}\right|^{2}-\left|g^{I}_{k}\right|^{2}\right)-2k\text{Re}\left(g^{I}_{k}g^{II*}_{k}\right)\right]^{(4)}\;,$
(38)

where $[\dots]^{(4)}$ indicates that $g^{I}_{k}$, $g^{II}_{k}$ are fourth
order states. In order for the energy density to be renormalized and all
divergences eliminated, the adiabatic states in Eq. (III) must be at least of
the adiabatic order at which the counterterms in Eq. (38) are divergent.

In order to obtain forms for the adiabatic states, one can obtain uncoupled
second order equations from Eqs. (32a) and (32b),

$\displaystyle\left(\partial_{\eta}^{2}+\frac{a^{\prime}}{a}\partial_{\eta}-ima^{\prime}+\omega^{2}\right)g^{I}_{k}=0\;,$
(39)
$\displaystyle\left(\partial_{\eta}^{2}+\frac{a^{\prime}}{a}\partial_{\eta}+ima^{\prime}+\omega^{2}\right)g^{II}_{k}=0\;,$
(40)

where

$\omega^{2}\equiv k^{2}+m^{2}a^{2}\;,$ (41)

and assume formal WKB series solutions, truncating at the desired adiabatic
order. However, Ref. [4] and later Refs. [5, 6] pointed out that the usual WKB
ansatz does not satisfy the normalization condition (33), so one must use a
modified WKB ansatz of the form

$\displaystyle g^{I}_{k}$
$\displaystyle=\sqrt{\frac{\omega+ma}{2\omega}}Fe^{-i\theta_{k}}\;,$ (42a)
$\displaystyle g^{II}_{k}$ $\displaystyle=\sqrt{\frac{\omega-
ma}{2\omega}}Ge^{-i\theta_{k}}\;,$ (42b)

and the functions

$\displaystyle F$ $\displaystyle=1+F^{(1)}+\dots+F^{(A)}\;,$ (43a)
$\displaystyle G$ $\displaystyle=1+G^{(1)}+\dots+G^{(A)}\;,$ (43b)
$\displaystyle\theta_{k}$
$\displaystyle=\int^{\eta}d\tilde{\eta}\left(\omega+\omega^{(1)}+\dots+\omega^{(A)}\right)$
(43c)

are determined by repeated substitution of Eqs. (42a) and (42b) into Eqs.
(32a), (32b), and (33). There is an ambiguity in the exact forms following
this method, but all local observables are independent of the ambiguity [4],
so one may fix the ambiguity by choosing $F^{(n)}(-m)=G^{(n)}(m)$ for each
order $n\geq 1$. However, we are only interested in using zeroth order states,
for which one obtains

$\displaystyle g^{I}_{k}$
$\displaystyle=\sqrt{\frac{\omega+ma}{2\omega}}e^{-i\theta_{k}}\;,$ (44a)
$\displaystyle g^{II}_{k}$ $\displaystyle=\sqrt{\frac{\omega-
ma}{2\omega}}e^{-i\theta_{k}}\;,$ (44b)

which have a normalization from Eq. (33) of

$\left|g^{I}_{k}\right|^{2}+\left|g^{II}_{k}\right|^{2}=1\;.$ (45)

Note that we will continue writing $\theta_{k}$ like we have in Eqs. (44a) and
(44b) for simplicity and assume it to be understood that only the zeroth order
term is kept. Substituting Eqs. (35a) and (35b) into the differential
equations (32a) and (32b), one obtains differential equations for $\alpha_{k}$
and $\beta_{k}$,

$\displaystyle\alpha_{k}^{\prime}$
$\displaystyle=\frac{-kma^{\prime}}{2\omega^{2}}\beta_{k}e^{2i\theta_{k}}\;,$
(46) $\displaystyle\beta_{k}^{\prime}$
$\displaystyle=\frac{kma^{\prime}}{2\omega^{2}}\alpha_{k}e^{-2i\theta_{k}}\;,$
(47)

and substituting them into the normalization condition (33), one finds

$\left|\alpha_{k}\right|^{2}+\left|\beta_{k}\right|^{2}=1\;.$ (48)

###  III.1 Renormalized Energy Density

In order to obtain finite results so that we may inspect the behavior of the
energy density, one must subtract counterterms up to fourth order from Eq.
(III). Order by order, these counterterms $\rho^{(n)}$ are

$\displaystyle\rho^{(0)}_{c}=$ $\displaystyle-\omega\;,$ (49)
$\displaystyle\rho^{(2)}_{c}=$ $\displaystyle\frac{k^{2}\omega^{\prime
2}}{8m^{2}a^{2}\omega^{3}}\;,=\frac{k^{2}m^{2}a^{\prime 2}}{8\omega^{5}}$ (50)
$\displaystyle\rho^{(4)}_{c}=$
$\displaystyle\mathcal{O}\\!\left(k^{-5}\right)\;.$ (51)

The fourth order counterterms produce finite contributions to the energy
density. Using the zeroth order adiabatic states (44a) and (44b) in Eq. (34)
and subtracting the counterterms (49)-(51), one obtains the renormalized
energy density

$\rho_{r}=\frac{1}{\pi^{2}a^{4}}\int_{0}^{\infty}dk\,k^{2}\left[2\omega\left|\beta_{k}\right|^{2}-\frac{k^{2}m^{2}a^{\prime
2}}{8\omega^{5}}\right]+\frac{2}{2880\pi^{2}}\left[-\frac{1}{2}\tensor{*}[^{(1)}]{H}{{}_{00}}+\frac{11}{2}\tensor{*}[^{(3)}]{H}{{}_{00}}\right]\;,$
(52)

where the finite renormalization terms coming from the fourth order
counterterms [10],

$\displaystyle\frac{2}{2880\pi^{2}}\left[-\frac{1}{2}\tensor{*}[^{(1)}]{H}{{}_{00}}+\frac{11}{2}\tensor{*}[^{(3)}]{H}{{}_{00}}\right]$
$\displaystyle=\frac{2}{2880\pi^{2}a^{4}}\left(\frac{33a^{\prime
4}}{2a^{4}}+\frac{18a^{\prime}a^{\prime\prime\prime}}{a^{2}}-\frac{9a^{\prime\prime
2}}{a^{2}}-\frac{36a^{\prime 2}a^{\prime\prime}}{a^{3}}\right)\;,$ (53)

are assumed small beyond the Planck era.

We ultimately will attempt to solve the Friedmann equation

$H^{2}=\frac{8\pi}{3}\rho\;,$ (54)

where

$H\equiv\frac{a^{\prime}}{a^{2}}\;,$ (55)

and $\rho$ contains $\rho_{r}$ and may also include other terms such as a
cosmological constant or other classical contributions. Because we are working
with a semiclassical approximation, we do not assume our analysis to be valid
during the Planck era. Hence we will work starting at an initial time
$\eta_{0}$, which is assumed to be after the Planck era and corresponds to a
scale factor $a_{0}$ that is above the Planck scale, and demonstrate that
$\rho_{r}\sim a^{-4}$ for some region $\eta_{0}<\eta<\eta_{1}$.

Above the Planck scale, the Hubble parameter $H\lesssim 1$, so given Eqs. (55)
and (53), we will insist that the following set of inequalities of derivatives
of the scale factor must hold:

$\displaystyle a^{\prime}$ $\displaystyle\lesssim a^{2}\;,$ (56)
$\displaystyle a^{\prime\prime}$ $\displaystyle\lesssim a^{3}\;,$ (57)
$\displaystyle a^{\prime}a^{\prime\prime\prime}$ $\displaystyle\lesssim
a^{6}\;.$ (58)

We will use these inequalities in order to compute the integral in Eq. (52).

###  III.2 Bounds on the Renormalized Energy Density

We will begin by splitting the integral into infrared and ultraviolet regions
$I_{1}$ and $I_{2}$ by a cutoff $k_{c}$. The infrared contribution is

$\displaystyle I_{1}$
$\displaystyle=\int_{0}^{k_{c}}dk\,k^{2}\left[2\omega\left|\beta_{k}\right|^{2}-\frac{k^{2}m^{2}a^{\prime
2}}{8\omega^{5}}\right]$
$\displaystyle=\int_{0}^{k_{c}}dk\,k^{2}\left[2k\left|\beta_{k}\right|^{2}+2\left(\omega-k\right)\left|\beta_{k}\right|^{2}-\frac{k^{2}m^{2}a^{\prime
2}}{8\omega^{5}}\right]\;,$ (59)

and the ultraviolet contribution is

$\displaystyle I_{2}$
$\displaystyle=\int_{k_{c}}^{\infty}dk\,k^{2}\left[2\omega\left|\beta_{k}\right|^{2}-\frac{k^{2}m^{2}a^{\prime
2}}{8\omega^{5}}\right]\;.$ (60)

At time $\eta_{0}$, we define

$B_{0}\equiv\int_{0}^{k_{c}}dk\,2k^{3}\left|\beta_{k}(\eta_{0})\right|^{2}\;.$
(61)

If $B_{0}$ is the dominant contribution to the renormalized energy density at
time $\eta$, then $\rho_{r}\propto a^{-4}$ as desired. However, as one may
anticipate given the apparent logarithmically divergent term in Eq. (52), this
may not always be the case. We will investigate this in the following
sections. Given that $\left|\beta_{k}(\eta_{0})\right|\leq 1$ from Eq. (48),
one finds from Eq. (61) a lower bound on $k_{c}$ of

$k_{c}\gtrsim B_{0}^{1/4}\;.$ (62)

####  III.2.1 Ultraviolet Region

One is tempted to assume $\beta_{k}\rightarrow 0$ sufficiently quickly at high
$k$, as is often done with scalar fields [3]. However, because the final term
in Eq. (60) produces a logarithm divergence, it is evident that doing so will
introduce a divergent contribution to the energy density. This situation is
occuring because until now we have been working with a zeroth order
parameterization of the states, but because the logarithmic divergence comes
in at higher than zeroth adiabatic order, one would expect to need to work
with at least a second order parameterization to eliminate this higher-order
divergence. These problems would indeed disappear working with a higher-order
parameterization, but it becomes much more difficult to analytically obtain
generic bounds on $I_{2}$ doing so.

One can instead continue to work with a zeroth-order parameterization using
$\alpha_{k}$ and $\beta_{k}$. We can understand what the appropriate higher
order states look like at large $k$ by integrating Eq. (47) by parts, which
shows that the asymptotic behavior of $\beta_{k}$ will take the form

$\beta_{k}\xrightarrow[k\to\infty]{}-\frac{ikma^{\prime}}{4\omega^{3}}\alpha_{k}e^{-2i\theta_{k}}\;.$
(63)

This motivates us to use the parameterization $\bar{\beta}_{k}$, defined as

$\bar{\beta_{k}}\equiv\beta_{k}+\frac{ikma^{\prime}}{4\omega^{3}}\alpha_{k}e^{-2i\theta_{k}}\;,$
(64)

as the appropriate description to use in calculating the energy density at
large $k$. We then anticipate that $\bar{\beta}_{k}$ will fall faster than the
leading order term at large $k$, so

$\left|\bar{\beta}_{k}(\eta_{0})\right|<\frac{A_{0}}{k^{2}}\left(\frac{k_{c}}{k}\right)^{b_{0}}$
(65)

for $k>k_{c}$, some $b_{0}>0$, and $A_{0}$ independent of $k$. The ultraviolet
integral (60) written in terms of $\bar{\beta}_{k}$ is

$\displaystyle I_{2}$
$\displaystyle=\int_{k_{c}}^{\infty}dk\,k^{2}\left[\left(1-\frac{k^{2}m^{2}a^{\prime
2}}{16\omega^{6}}\right)\left(2\omega\left|\bar{\beta}_{k}\right|^{2}+\frac{kma^{\prime}}{\omega^{2}}\text{Im}\\!\left(\alpha_{k}\bar{\beta}_{k}^{*}e^{-2i\theta_{k}}\right)\right)-\frac{k^{4}m^{4}a^{\prime
4}}{128\omega^{11}}\left|\alpha_{k}\right|^{2}\right]\;.$ (66)

Because $\bar{\beta}_{k}$ encodes the cancellation of the divergent term in
the renormalized energy density, one expects that the contributions to the
energy density from Eq. (66) will converge.

In order for the energy density to be radiation-dominated the ultraviolet
contribution (66) must either be the dominant contribution and itself
radiation-dominated or be subdominant to the radiation-dominated part of the
infrared contribution (59). As we will demonstrate, every term in Eq. (66) is
in fact subdominant to $B_{0}$ (61), which itself produces a radiation-
dominated term in the energy density, and hence the latter is true.

To show this, we will first simplify the first factor in Eq. (66) by using
Eqs. (62) and (56) and $\omega>k>k_{c}$ to show that $k^{2}m^{2}a^{\prime
2}\ll 16\omega^{6}$ and therefore $1-\frac{k^{2}m^{2}a^{\prime
2}}{16\omega^{6}}\approx 1$, provided the condition $a^{2}\ll
B_{0}^{1/2}m^{-1}$ is satisfied, so we will need

$a\ll B_{0}^{1/4}m^{-1}\;.$ (67)

This is the first of several conditions we will need. We will consider the
complete set of conditions collectively later.

With this simplification, one finds a bound on $I_{2}$ of

$\displaystyle\left|I_{2}\right|$
$\displaystyle\lesssim\int_{k_{c}}^{\infty}dk\left[2k^{2}\omega\left|\bar{\beta}_{k}\right|^{2}+\frac{1}{2}kma^{\prime}\left|\bar{\beta}_{k}\right|\right]+\frac{m^{4}a^{\prime
4}}{512k_{c}^{4}}$
$\displaystyle<\int_{k_{c}}^{\infty}dk\left[2k^{3}\left|\bar{\beta}_{k}\right|^{2}+km^{2}a^{2}\left|\bar{\beta}_{k}\right|^{2}+\frac{1}{2}kma^{\prime}\left|\bar{\beta}_{k}\right|\right]+\frac{m^{4}a^{\prime
4}}{512k_{c}^{4}}\;,$ (68)

where we have bounded the integral on the final term using the normalization
(48) to obtain $\left|\alpha_{k}\right|<1$. For $B_{0}$ to dominate, the
contributions from $I_{2}$ must remain less than $B_{0}$, and therefore the
full contribution from $I_{2}(\eta)=I_{2}(\eta_{0})+\Delta I_{2}$ must be less
than $B_{0}$. We will find the conditions under which
$\left|I_{2}(\eta_{0})\right|$ and $\left|\Delta I_{2}\right|$, and therefore
$\left|I_{2}(\eta)\right|$, are subdominant to $B_{0}$.

Given Eq. (65), one finds from Eq. (68) that the contributions from the
integrand of $\left|I_{2}(\eta_{0})\right|$ are subdominant to $B_{0}$
provided that

$\displaystyle A_{0}^{2}$ $\displaystyle\ll B_{0}\;,$ (69a) $\displaystyle
m^{2}a^{2}A_{0}^{2}$ $\displaystyle\ll B_{0}^{3/2}\;,$ (69b) $\displaystyle
ma^{\prime}A_{0}$ $\displaystyle\ll B_{0}\;.$ (69c)

To satisfy (65), one may increase $k_{c}$ which allows for decreasing $A_{0}$
and ensuring (69a) can be satisfied. Furthermore, provided Eq. (69a) is
satisfied and using Eq. (56), one can show Eqs. (69b) and (69c) become

$\displaystyle a$ $\displaystyle\ll B_{0}^{1/4}m^{-1}\;,$ (70) $\displaystyle
a$ $\displaystyle\ll B_{0}^{1/4}m^{-1/2}\;.$ (71)

The final term in (68) is also less than $B_{0}$ provided (71) is satisfied.

For the bound on the integral in $\Delta I_{2}$ to converge,
$\Delta\bar{\beta}_{k}$ must fall faster than $k^{-2}$. We will assume that
$\Delta\bar{\beta}_{k}$ falls at least as fast as $k^{-2-b_{\Delta}}$, and
then we will need to show

$\left|\Delta\bar{\beta}_{k}\right|<\frac{A_{\Delta}}{k^{2}}\left(\frac{k_{c}}{k}\right)^{b_{\Delta}}\;,$
(72)

with $A_{\Delta}$ independent of $k$ and $b_{\Delta}>0$ chosen appropriately
for each term contributing to $\Delta\bar{\beta}_{k}$. One finds from Eq. (68)
that the contributions from $\Delta I_{2}$ are less than $B_{0}$ provided that

$\displaystyle A_{\Delta}^{2}$ $\displaystyle\ll B_{0}\;,$ (73) $\displaystyle
m^{2}a^{2}A_{\Delta}^{2}$ $\displaystyle\ll B_{0}^{3/2}\;,$ (74)
$\displaystyle ma^{\prime}A_{\Delta}$ $\displaystyle\ll B_{0}\;.$ (75)

One can obtain bounds on the contributions to $\Delta\bar{\beta}_{k}$, and
hence the conditions under which Eqs. (73)–(75) are satisfied, using the
differential equation for $\bar{\beta}_{k}$. From Eqs. (47) and (64), the
differential equation is

$\displaystyle\bar{\beta}_{k}^{\prime}$
$\displaystyle=-\frac{ikma^{\prime\prime}}{4\omega^{3}}\alpha_{k}e^{-2i\theta_{k}}+\frac{ik^{2}m^{2}a^{\prime
2}}{8\omega^{5}}\beta_{k}+\frac{3ikm^{3}aa^{\prime
2}}{4\omega^{5}}\alpha_{k}e^{-2i\theta_{k}}\;.$ (76)

Writing
$\bar{\beta}_{k}(\eta)=\bar{\beta}_{k}(\eta_{0})+\Delta\bar{\beta}_{k}$, one
then finds by integrating by parts on the phase

$\displaystyle\Delta\bar{\beta}_{k}=$
$\displaystyle\int_{\eta_{0}}^{\eta}dx\left[\frac{-ikma^{\prime\prime}}{4\omega^{3}}\alpha_{k}e^{-2i\theta_{k}}+\frac{ik^{2}m^{2}a^{\prime
2}}{8\omega^{5}}\beta_{k}+\frac{3ikm^{3}aa^{\prime
2}}{4\omega^{5}}\alpha_{k}e^{-2i\theta_{k}}\right]$ $\displaystyle=$
$\displaystyle\left.\frac{kma^{\prime\prime}}{8\omega^{4}}\alpha_{k}e^{-2i\theta_{k}}\right|_{\eta_{0}}^{\eta}+\int_{\eta_{0}}^{\eta}dx\,\bigg{[}\frac{-kma^{\prime\prime\prime}}{8\omega^{4}}\alpha_{k}e^{-2i\theta_{k}}+\frac{3km^{3}aa^{\prime}a^{\prime\prime}}{8\omega^{6}}\alpha_{k}e^{-2i\theta_{k}}+\frac{k^{2}m^{2}a^{\prime}a^{\prime\prime}}{16\omega^{6}}\beta_{k}$
$\displaystyle+\frac{ik^{2}m^{2}a^{\prime
2}}{8\omega^{5}}\beta_{k}+\frac{3ikm^{3}aa^{\prime
2}}{4\omega^{5}}\alpha_{k}e^{-2i\theta_{k}}\bigg{]}\;,$ (77)

and therefore

$\displaystyle\left|\Delta\bar{\beta}_{k}\right|<$
$\displaystyle\frac{m\left(\left|a^{\prime\prime}(\eta_{0})\right|+\left|a^{\prime\prime}(\eta)\right|\right)}{8k^{3}}+\frac{1}{8k^{3}}\int_{\eta_{0}}^{\eta}dx\left(m\left|a^{\prime\prime\prime}\right|+m^{2}a^{\prime
2}\right)+\frac{1}{16k^{4}}\int_{\eta_{0}}^{\eta}dx\,m^{2}\left|a^{\prime}\right|\left|a^{\prime\prime}\right|$
$\displaystyle+\frac{3}{4k^{4}}\int_{\eta_{0}}^{\eta}dx\,m^{3}aa^{\prime
2}+\frac{3}{8k^{5}}\int_{\eta_{0}}^{\eta}dx\,m^{3}a\left|a^{\prime}\right|\left|a^{\prime\prime}\right|\;.$
(78)

Each term in Eq. (78) can be written in the form of Eq. (72) to obtain
conditions under which each will satisfy Eq. (73). Term by term, using the
inequalities in Eqs. (56) and (57), we can obtain conditions on the scale
factor under which each term is subdominant to the infrared contribution. In
order to obtain these conditions, we will assume $a^{\prime}$,
$a^{\prime\prime}$, and $a^{\prime\prime\prime}$ have definite signs; that is,
each of them is either always positive or always negative. The conditions are

$\displaystyle\frac{ma^{\prime\prime}}{k^{3}}=\frac{ma^{\prime\prime}}{k^{2}k_{c}}\left(\frac{k_{c}}{k}\right)\implies$
$\displaystyle B_{0}\gg\left(\frac{ma^{\prime\prime}}{k_{c}}\right)^{2}$
$\displaystyle a\ll B_{0}^{1/4}m^{-1/3}\;,$ (79a)
$\displaystyle\frac{1}{k^{3}}\int dx\,m^{2}a^{\prime 2}<\frac{1}{k^{3}}\int
dx\,m^{2}a^{2}a^{\prime}=\frac{m^{2}a^{3}}{k^{3}}$
$\displaystyle=\frac{m^{2}a^{3}}{k^{2}k_{c}}\left(\frac{k_{c}}{k}\right)$
$\displaystyle\implies$ $\displaystyle
B_{0}\gg\left(\frac{m^{2}a^{3}}{k_{c}}\right)$ $\displaystyle a\ll
B_{0}^{1/4}m^{-2/3}\;,$ (79b) $\displaystyle\frac{1}{k^{4}}\int
dx\,m^{2}a^{\prime}a^{\prime\prime}=\frac{m^{2}a^{\prime 2}}{k^{4}}$
$\displaystyle=\frac{m^{2}a^{\prime
2}}{k^{2}k_{c}^{2}}\left(\frac{k_{c}}{k}\right)^{2}$ $\displaystyle\implies$
$\displaystyle B_{0}\gg\left(\frac{m^{2}a^{\prime 2}}{k_{c}^{2}}\right)^{2}$
$\displaystyle a\ll B_{0}^{1/4}m^{-1/2}\;,$ (79c)
$\displaystyle\frac{1}{k^{4}}\int dx\,m^{3}aa^{\prime 2}<\frac{1}{k^{4}}\int
dx\,m^{3}a^{3}a^{\prime}=\frac{m^{3}a^{4}}{k^{4}}$
$\displaystyle=\frac{m^{3}a^{4}}{k^{2}k_{c}^{2}}\left(\frac{k_{c}}{k}\right)^{2}$
$\displaystyle\implies$ $\displaystyle
B_{0}\gg\left(\frac{m^{3}a^{4}}{k_{c}^{3}}\right)^{2}$ $\displaystyle a\ll
B_{0}^{1/4}m^{-3/4}\;,$ (79d) $\displaystyle\frac{1}{k^{5}}\int
dx\,m^{3}aa^{\prime}a^{\prime\prime}<\frac{1}{k^{5}}\int
dx\,m^{3}a^{4}a^{\prime}=\frac{m^{3}a^{5}}{k^{5}}$
$\displaystyle=\frac{m^{3}a^{5}}{k^{2}k_{c}^{3}}\left(\frac{k_{c}}{k}\right)^{3}$
$\displaystyle\implies$ $\displaystyle
B_{0}\gg\left(\frac{m^{3}a^{5}}{k_{c}^{3}}\right)$ $\displaystyle a\ll
B_{0}^{1/4}m^{-3/5}\;.$ (79e)

Even if one of the three quantities $a^{\prime}$, $a^{\prime\prime}$, or
$a^{\prime\prime\prime}$ is not of definite sign, one can subdivide the
integrals appearing in (78) into regions in which it is of definite sign, and
assuming the number of regions is not too large the sum of these regions can
be similarly bounded if the inequalities (79a)–(79e) are satisfied.

We are working at times beyond the Planck era during which the mass is small
compared to the Planck mass, and hence $m\ll 1$, so the strongest restriction
among Eq. (67) and Eqs. (79a)–(79e) for which the ultraviolet contributions
are subdominant to that from $B_{0}$ is

$a\ll B_{0}^{1/4}m^{-1/3}\;.$ (80)

Equation (80) is stronger than the conditions in Eqs. (70)–(71) and those in
Eqs. (74)–(75). One therefore has a complete set of conditions under which the
ultraviolet contribution is subdominant to that from $B_{0}$: $k_{c}$ must be
large enough such that Eq. (69a) is true, and $a$ must be small enough to
satisfy Eq. (80).

####  III.2.2 Infrared Region

It remains to be shown that $B_{0}$ is indeed the dominant contribution to the
energy density among all the terms in the infrared contribution (59). Given
Eq. (61), the infrared contribution (59) can be written

$\displaystyle I_{1}$
$\displaystyle=B_{0}+\left(B-B_{0}\right)+\int_{0}^{k_{c}}dk\,\left[2k^{2}\left(\omega-k\right)\left|\beta_{k}\right|^{2}-\frac{k^{4}m^{2}a^{\prime
2}}{4\omega^{5}}\right]\;,$ (81)

where we have defined

$B\equiv\int_{0}^{k_{c}}dk\,2k^{3}\left|\beta_{k}\right|^{2}\;.$ (82)

The $\left(B-B_{0}\right)$ term can be written

$\displaystyle B-B_{0}$
$\displaystyle=\int_{0}^{k_{c}}dk\,2k^{3}\left(\left|\beta_{k}\right|^{2}-\left|\beta_{k}(\eta_{0})\right|^{2}\right)$
$\displaystyle\leq\int_{0}^{k_{c}}dk\,2k^{3}\left(\left|\Delta\beta_{k}\right|^{2}+2\left|\beta_{k}(\eta_{0})\right|\left|\Delta\beta_{k}\right|\right)\;,$
(83)

where $\Delta\beta_{k}=\beta_{k}-\beta_{k}(\eta_{0})$. From Eq. (47),

$\displaystyle\left|\Delta\beta_{k}\right|<\int_{\eta_{0}}^{\eta}dx\,\frac{kma^{\prime}}{2\omega^{2}}=\frac{1}{2}\int_{\eta_{0}}^{\eta}dx\,\partial_{x}\tan^{-1}\left(\frac{ma}{k}\right)<\frac{ma}{2k}\;,$
(84)

which implies

$\displaystyle\int_{0}^{k_{c}}dk\,2k^{3}\left|\Delta\beta_{k}\right|^{2}<\frac{1}{4}k_{c}^{2}m^{2}a^{2}\;$
(85)

and

$\int_{0}^{k_{c}}dk\;4k^{3}\left|\beta_{k}(\eta_{0})\right|\left|\Delta\beta_{k}\right|<\frac{2}{3}mak_{c}^{3}\;,$
(86)

so, using Eq. (62), one finds $\left(B-B_{0}\right)$ is dominated by $B_{0}$
if

$a\ll B_{0}^{1/4}m^{-1}\;.$ (87)

Similarly, using the normalization condition (48) to bound
$\left|\beta_{k}\right|^{2}\leq 1$, the next term in Eq. (81) is

$\displaystyle\int_{0}^{k_{c}}dk\,2k^{2}\left(\omega-k\right)\left|\beta_{k}\right|^{2}\leq\frac{1}{2}k_{c}^{2}m^{2}a^{2}\;,$
(88)

which is also less than $B_{0}$ if $a\ll B_{0}^{1/4}m^{-1}$. This condition is
weaker than Eq. (80) and hence will be satisfied if the set of conditions
under which the ultraviolet contribution is subdominant to $B_{0}$ is
satisfied.

The last term in $I_{1}$ is

$\displaystyle\int_{0}^{k_{c}}dk\,\frac{-k^{4}m^{2}a^{\prime 2}}{4\omega^{5}}$
$\displaystyle=\frac{m^{2}a^{\prime
2}}{4}\left(\frac{k_{c}^{3}}{3\omega_{c}^{3}}+\frac{k_{c}}{\omega_{c}}\right)+\frac{1}{4}m^{2}a^{\prime
2}\ln\left(\frac{ma}{\omega_{c}+k_{c}}\right)$
$\displaystyle\approxeq\frac{1}{3}m^{2}a^{\prime 2}+\frac{1}{4}m^{2}a^{\prime
2}\ln\left(\frac{ma}{2k_{c}}\right)\;,$ (89)

where we have used $a\ll B_{0}^{1/4}m^{-1}<k_{c}$ to simplify
$\omega_{c}^{2}\equiv k_{c}^{2}+m^{2}a^{2}\approx k_{c}^{2}$.

###  III.3 Friedmann Equation

The energy density is given by the sum of the contribution from the infrared
and ultraviolet regions. As shown in Secs. III.2.1–III.2.2, provided one is
looking early enough such that Eq. (87) and hence Eq. (80) are satisfied, the
energy density is dominated by the contributions from $B_{0}$ and the
logarithm in Eq. (89):

$\displaystyle\rho_{r}\approxeq\frac{1}{\pi^{2}a^{4}}\left[B_{0}+\frac{1}{4}m^{2}a^{\prime
2}\ln\left(\frac{ma}{2k_{c}}\right)\right]\;.$ (90)

Using the Friedmann equation

$\left(\frac{a^{\prime}}{a^{2}}\right)^{2}=\frac{8\pi}{3}\rho_{r}$ (91)

with the assumption that the contribution from massless spin-half fields
dominates the energy density, one finds

$a^{\prime
2}\approxeq\frac{8B_{0}}{3\pi}\left[1-\frac{2}{3\pi}m^{2}\ln\left(\frac{ma}{2k_{c}}\right)\right]^{-1}\;.$
(92)

The renormalized energy density can be written using (92) as

$\rho_{r}\approxeq\frac{B_{0}}{\pi^{2}a^{4}}\left\\{1+\frac{2}{3\pi}m^{2}\left[\ln\left(\frac{2k_{c}}{ma_{0}}\right)-\ln\left(\frac{a}{a_{0}}\right)\right]\right\\}^{-1}\;.$
(93)

At this point, we will choose $a_{0}\approx B_{0}^{1/4}$ to ensure we are past
the Planck era, and given the condition (80) we have $a_{1}\approx
B_{0}^{1/4}m^{-1/3}$, so $m^{2}\ln(a/a_{0})<m^{2}\ln(m^{-1/3})\ll 1$ is
irrelevant compared to the $1$ term. Hence, in the range
$B_{0}^{1/4}<a<B_{0}^{1/4}m^{-1/3}$, the energy density will indeed be
radiation dominated.

Note that $a^{\prime}$ (92) is approximately constant in the range of
interest, so $a^{\prime\prime}$ will be strongly suppressed. This allows one
to relax some of the accumulated constraints, for example Eq. (79a), implying
$\rho\propto a^{4}$ over a larger range of $a$. However, limits such as (87)
will generally not relax, but this is not surprising because at
$a=B_{0}^{1/4}m^{-1}$ one expects fermions to become non-relativistic. Of
course, realistic models would not only have fermion fields but also an
inflaton field, which would look like a cosmological constant, but, until the
mass term becomes important at $a\sim B_{0}^{1/4}m^{-1}$ or inflation takes
over, things will still be radiation dominated.

##  IV Discussion and Conclusion

In this paper we have analyzed the early-Universe, preinflationary behavior of
massless vector fields of spin-1 and massive or massless fermion fields of
spin-$\frac{1}{2}$ in the semiclassical approximation. We showed for a range
of conformal time after the Planck era that both types of fields have
radiation-dominated behavior. Along with the same result obtained for scalar
fields in [3], we have demonstrated that all matter fields which one might
anticipate to play a role in the preinflationary era do in fact produce a
radiation-dominated energy density that is typically assumed in inflationary
models.

In Sec. II, we summarized the adiabatic renormalization procedure for a
massless vector field in a spatially flat FLRW universe, following the
groundwork laid out in [7]. We then used this procedure to renormalize the
energy density contribution for such a field following a parameterization of
the mode functions in terms of adiabatic states. We found the renormalized
energy density to have a radiation-dominated form similar to that for the
scalar field in [3].

In Sec. III, we summarized a modified version of the adiabatic renormalization
procedure for a massive or massless fermion field, using the modified WKB
ansatz given in [4]. We found that parameterizing the mode functions in terms
of adiabatic states required higher than zeroth order contributions in order
to properly renormalize the energy density at high energies. We then made use
of the leading second adiabatic-order contributions to the parameterization
coefficients to obtain the leading order behavior of such high-energy
contributions and used this to show the high-energy contributions are in fact
subdominant to the radiation-like term in the remaining energy density
contributions given a set of constraints on the scale factor. These
constraints are more stringent than those required for the scalar field [3],
for which only a constraint on $(a^{2})^{\prime\prime}$ was necessary. We then
found the next-to-leading order behavior of the energy density to be a
logarithmic contribution and showed that it is subdominant compared to the
radiation-dominated energy density for a the range of conformal time beyond
the Planck era, where the constraint on the scale factor coming from the
analysis of the high-energy contribution served as the upper bound on the
range.

Our anaylsis of the fermion field assumed that fermionic matter was the
dominant contribution to the energy density in order to make use of the
Friedmann equation to obtain the approximate leading order behavior. This
result for fermions is different than that for scalars [3] and vectors, which
work in any case. Our argument for fermions also only applies for the range
$a_{0}<a<a_{1}$ described in Sec. III, though we expect this range to be in
the preinflationary era prior to the fields becoming matter dominated. This
result is weaker than those for the other two fields, but in application it is
not. One would not be able to observe these matter fields during this era
directly but rather through their effects on other phenomena, such as in the
Cosmic Microwave Background, so the results are not necessarily weaker in
application.

Having demonstrated in the semiclassical approximation that the matter content
in the post-Planck, preinflationary era will indeed be radiation-dominated,
one can proceed with the procedure for the inflaton field in [3]. There, it
was assumed that spin-$\frac{1}{2}$ and spin-$1$ massive fields could be
modeled with conformally coupled massive scalar fields, which were argued to
be radiation-dominated themselves in the preinflationary era. Our results that
spin-$\frac{1}{2}$ fermion fields and massless spin-$1$ vector fields are
radiation-dominated, themselves, supports the procedure in [3], and the
analysis for obtaining a renormalized energy density for a universe with a
mixture of cosmological constant and classical radiation, and hence for
obtaining the power spectrum of the Cosmic Microwave Background, is identical.

###### Acknowledgements.

We would like to thank Paul Anderson for helpful insights and discussion.

## References

  * Das  _et al._ [2015] S. Das, G. Goswami, J. Prasad, and R. Rangarajan, Revisiting a pre-inflationary radiation era and its effect on the CMB power spectrum, J. Cosmol. Astropart. Phys. 2015 (06), 001.
  * Bennett  _et al._ [2003] C. L. Bennett, M. Halpern, G. Hinshaw, N. Jarosik, A. Kogut, M. Limon, S. S. Meyer, L. Page, D. N. Spergel, G. S. Tucker, E. Wollack, E. L. Wright, C. Barnes, M. R. Greason, R. S. Hill, E. Komatsu, M. R. Nolta, N. Odegard, H. V. Peiris, L. Verde, and J. L. Weiland, First-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Preliminary Maps and Basic Results, Astrophys. J. Suppl. 148, 1 (2003).
  * Anderson  _et al._ [2020] P. R. Anderson, E. D. Carlson, T. M. Ordines, and B. Hicks, Semiclassical predictions regarding a preinflationary era and its effects on the power spectrum, Phys. Rev. D 102, 063528 (2020).
  * Landete  _et al._ [2013] A. Landete, J. Navarro-Salas, and F. Torrentí, Adiabatic regularization for spin-$1/2$ fields, Phys. Rev. D 88, 061501(R) (2013).
  * Landete  _et al._ [2014] A. Landete, J. Navarro-Salas, and F. Torrentí, Adiabatic regularization and particle creation for spin one-half fields, Phys. Rev. D 89, 044030 (2014).
  * del Rio  _et al._ [2014] A. del Rio, J. Navarro-Salas, and F. Torrentí, Renormalized stress-energy tensor for spin-$1/2$ fields in expanding universes, Phys. Rev. D 90, 084017 (2014).
  * Chu and Koyama [2017] C.-S. Chu and Y. Koyama, Adiabatic regularization for gauge fields and the conformal anomaly, Phys. Rev. D 95, 065025 (2017).
  * Schambach and Sanders [2018] M. Schambach and K. Sanders, The proca field in curved spacetimes and its zero mass limit, Reports on Mathematical Physics 82, 203 (2018).
  * Wald [1978] R. M. Wald, Axiomatic renormalization of the stress tensor of a conformally invariant field in conformally flat spacetimes, Annals of Physics 110, 472 (1978).
  * Birrell and Davies [1982] N. Birrell and P. Davies,  _Quantum Fields in Curved Space_ , edited by P. Landshoff, W. McCrea, D. Sciama, and S. Weinberg (Cambridge University Press, 1982).
  * Endo [1984] R. Endo, Gauge Dependence of the Grabvitational Conformal Anomaly for the Electromagnetic Field, Progress of Theoretical Physics 71, 1366 (1984), https://academic.oup.com/ptp/article-pdf/71/6/1366/5313387/71-6-1366.pdf .
  * Torrentí [2015] F. Torrentí, Adiabatic renormalization of the stress-energy tensor for spin-1/2 fields., Journal of Physics: Conference Series 600, 012029 (2015).
  * Barbero G.  _et al._ [2018] J. F. Barbero G., A. Ferreiro, J. Navarro-Salas, and E. J. S. Villaseñor, Adiabatic expansions for dirac fields, renormalization, and anomalies, Phys. Rev. D 98, 025016 (2018).
  * Beltrán-Palau  _et al._ [2020] P. Beltrán-Palau, J. Navarro-Salas, and S. Pla, Adiabatic regularization for dirac fields in time-varying electric backgrounds, Phys. Rev. D 101, 105014 (2020).