Numerical study of the Schwinger effect
in axion inflation

Yann Cado, Mariano Quirós

Institut de Física d’Altes Energies (IFAE) and
The Barcelona Institute of Science and Technology (BIST),
Campus UAB, 08193 Bellaterra, Barcelona, Spain

Abstract

> Previous studies demonstrate that the inflaton, when coupled to the
> hypercharge Chern-Simons density, can source an explosive production of
> helical hypermagnetic fields. Then, in the absence of fermion production,
> those fields have the capability of preheating the Universe after inflation
> and triggering a successful baryogenesis mechanism at the electroweak phase
> transition. In the presence of fermion production however, we expect a
> strong damping of the gauge fields production from the fermion backreaction,
> a phenomenon called Schwinger effect, thus jeopardizing their original
> capabilities. Using numerical methods we study the backreaction on the
> generated gauge fields and revisit the processes of gauge preheating and
> baryogenesis in the presence of the Schwinger effect. We have found that
> gauge preheating is very unlikely, while still having a sizable window in
> the parameter space to achieve the baryon asymmetry of the Universe at the
> electroweak phase transition.

###### Contents

  1. 1 Introduction
  2. 2 The model
    1. 2.1 Equations of motion
    2. 2.2 Observable quantities
    3. 2.3 The gauge vacuum
  3. 3 Slow Roll Analysis
    1. 3.1 Absence of Schwinger effect
    2. 3.2 Presence of Schwinger effect
    3. 3.3 Numerical results at the end of inflation
      1. 3.3.1 Constant $\Delta$ and $\xi$ approximation
      2. 3.3.2 Variable $\Delta$ and $\xi$
  4. 4 Full Analysis
    1. 4.1 Full numerical results at the end of inflation
    2. 4.2 Inflationary models
      1. 4.2.1 $\alpha$-attractor models
      2. 4.2.2 Hilltop quartic models
    3. 4.3 Numerical results beyond the end of inflation
    4. 4.4 End of reheating
    5. 4.5 Baryon asymmetry
  5. 5 Conclusions
  6. A Numerical method: the slow roll case
  7. B Numerical method: full analysis

##  1 Introduction

Cosmological inflation Guth:1980zm ; Linde:1981mu ; Albrecht:1982wi  is
nowadays a well-established paradigm to solve the classical (flatness,
horizon, …) problems of the Standard Cosmological Model, and to generate the
primordial density perturbations giving rise to the present Universe
structure. The achievements of cosmological inflation usually require the
presence of one (or several) scalar field –the inflaton– giving rise to
physics beyond the Standard Model (SM) of Particle Physics (BSM). In this way,
along with the classical problems of the SM (hierarchy problem, baryogenesis,
strong CP-problem, dark matter,…), cosmological inflation provides yet another
motivation for BSM physics.

Although the existence of a period of cosmological inflation is pretty well
stablished by observational cosmological data Planck:2018jri , there is no
consensus on a detailed model. An interesting candidate for the inflaton is a
pseudoscalar field $\phi$, denoted in this paper as axionlike particle111With
an (obvious) abuse of language, we are identifying in this paper axions with
axionlike (or pseudoscalar) particles, which allow a wider choice for
inflationary potentials Adshead:2019lbr ; Adshead:2019igv ., which can then
couple to the Chern-Simons density $F_{\mu\nu}\tilde{F}^{\mu\nu}$ of a $U(1)$
gauge field. In this case, and depending on the size of the coupling of the
inflaton to the Chern-Simons term, there can be an explosive production of
helical gauge fields at the end of inflation Anber:2009ua ; Anber:2015yca ;
Cado:2016kdp . This exponential production can dominate the energy density of
the Universe during the coherent oscillations of the inflaton around its
minimum, a phenomenon dubbed as gauge preheating Adshead:2015pva ;
Sfakianakis:2018lzf ; Cuissa:2018oiw , and lead to a rapid production of
inhomogeneities sourcing a significant gravitational wave background, leading
to strong constraints on the inflaton Chern-Simons coupling from the Planck
(and future CMB-S4) limits on the net energy density in gravitational waves
Adshead:2019lbr ; Adshead:2019igv .

When the gauge field is identified with the SM hypercharge, $Y_{\mu}$, with
strength $Y_{\mu\nu}$, the inflaton coupling to the Chern-Simons density
$Y_{\mu\nu}\tilde{Y}^{\mu\nu}$ gives rise to the production of helical
hypermagnetic fields which can then survive until the electroweak phase
transition (crossover), and trigger the baryon asymmetry of the Universe (BAU)
Kamada:2016eeb ; Kamada:2016cnb ; Jimenez:2017cdr ; Domcke:2019mnd ;
Cado:2021bia ; Cado:2022evn . However, in the presence of strong gauge fields,
light fermions charged under the gauge group are produced by the backreaction
of gauge fields that source the fermion equations of motion (EoM)
Domcke:2018eki ; Kitamoto:2021wzl . The corresponding currents can then, in
turn, backreact on the produced gauge fields, a phenomenon called Schwinger
effect (see e.g. Ref. Cohen:2008wz ). The backreaction of fermion currents on
the produced gauge fields acts as a damping force in the explosive production
of helical gauge fields, and many of the conclusions from the gauge field
production should be revised in the presence of the Schwinger effect222One
possible way out is if there are no light charged fields when gauge fields are
produced, a condition that is not fulfilled here., in particular those
concerning the preheating capabilities and the baryogenesis mechanism.

In this paper, we will study the effect of the Schwinger particle production
on the helical hypermagnetic fields produced at the end of inflation, and in
particular its influence on the gauge preheating efficiency and baryogenesis
capability. In order to consider the backreaction of the produced gauge fields
on the inflationary equations of motion and that of the Schwinger effect on
the gauge field production, we will use numerical methods, in particular, the
fourth order Runge-Kutta (RK4) algorithm. Our numerical results are validated
as they overlap with some recent semianalytical methods and the gradient
expansion formalism of Refs. Sobol:2019xls ; Gorbar:2021rlt ; Gorbar:2021zlr ;
Fujita:2022fwc . Our general finding is that the gauge field production is
much less explosive than in the absence of the Schwinger effect, which will
jeopardize the conclusions concerning the possibility of gauge preheating,
although they leave an open window for baryogenesis.

The contents of this paper are as follows. In Sec. 2 we present the general
lines of the model and the methods we will consider, including the relevant
equations of motion in momentum space and the observable quantities we will
compute. In Sec. 3 we will present numerical results for the gauge sector
assuming the slow roll conditions in the inflaton equations of motion. In
order to check the validity of our approach, numerical results will be
compared with some estimates from the literature, namely the backreactionless
solution, the Schwinger equilibrium and maximal estimates as well as the
gradient expansion formalism where dynamical results are obtained analytically
and numerically in configuration space. Some details about the numerical
methods will be explained in App. A. In Sec. 4 we will perform the full
numerical calculation for two kinds of models that predict cosmological
observables in agreement with the observed values: the $\alpha$-attractor
models and the quartic hilltop models. In all cases the gauge preheating
efficiency does not seem good enough to ensure complete reheating, which has
to be completed by other perturbative or nonperturbative mechanisms. Moreover
we have reanalyzed the baryogenesis predictions in the presence of the
Schwinger effect and found a sizable window where the BAU is correctly
predicted. Again, some details about the numerical methods we used are
described in App. B. Finally, our conclusions are presented in Sec. 5.

##  2 The model

The model action is given by

$\displaystyle S=$ $\displaystyle\int
d^{4}x\left[\sqrt{-g}\left(\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-\dfrac{1}{4}Y_{\mu\nu}Y^{\mu\nu}-V(\phi)\right)-\dfrac{\phi}{4f_{\phi}}\;Y_{\mu\nu}\tilde{Y}^{\mu\nu}\right]+\int
d^{4}x\;\sqrt{-g}\;i\,\bar{\psi}\cancel{\mathcal{D}}\psi,$ (2.1)

where $\phi$ is the pseudoscalar inflaton, $V$ the inflaton potential, and
$f_{\phi}$ provides the inverse coupling of the inflaton to the Chern-Simons
term. $Y^{\mu\nu}$ is the field strength of the hypercharge gauge field
$Y^{\mu}$ and
$\widetilde{Y}^{\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}Y_{\rho\sigma}$
its dual tensor. We also have included the interaction of fermionic currents,
corresponding to hypercharge $Q_{Y}$ fermions, with the hypercharge fields
(encoded in the covariant derivative
$\mathcal{D}_{\mu}\equiv\partial_{\mu}-g^{\prime}Q_{Y}A_{\mu}$). All gauge
field quantities are $U(1)$ hypercharge fields, i.e. $\bm{A}_{Y}$,
$\bm{E}_{Y}$, $\bm{B}_{Y}$, etc. To make the notation lighter, we drop the
index $Y$ as there will be no ordinary electromagnetic fields in this work.

###  2.1 Equations of motion

Variation of the action with respect to $\phi$ and the hypercharge gauge field
$A_{\mu}=(A_{0},\bm{A})$ leads to the gauge equations of motion in the
radiation gauge, $A_{0}=0$ and $\bm{\nabla}\cdot\bm{A}=0$,

$\displaystyle\ddot{\phi}+3H\dot{\phi}+V^{\prime}(\phi)=\dfrac{\bm{E}\cdot\bm{B}}{a^{4}f_{\phi}},$
(2.2a)
$\displaystyle\left(\dfrac{\partial^{2}}{\partial\tau^{2}}-\nabla^{2}-\dfrac{a\,\dot{\phi}}{f_{\phi}}\;\bm{\nabla}\times\right)\bm{A}=\bm{J},$
(2.2b)

where we have used $Y_{\mu\nu}\tilde{Y}^{\mu\nu}=-4\,\bm{E}\cdot\bm{B}$ and
$J^{\mu}=(\rho_{c},\bm{J})=ig^{\prime}Q_{Y}\bar{\psi}\gamma^{\mu}\psi$. We
assume that initially the Universe does not contain any asymmetry of charged
particles and that these ones are produced only later in particle-antiparticle
pairs. Therefore, we set the charge density to zero, $\rho_{c}=0$. Finally,
the current $\bm{J}$ is given by the Ohm’s law

$\bm{J}=\sigma\bm{E}=-\sigma\dfrac{\partial\bm{A}}{\partial\tau},$ (2.3)

where $\sigma$ is the generalized conductivity, which will be defined later.

As it can be seen from the above system, we use cosmic time $t$ for the
inflaton dependence and the conformal time $\tau$, defined by
$g_{\mu\nu}=a^{2}(\tau)\,\eta_{\mu\nu}$, for the gauge field dependence. We
denote the derivative with respect to conformal time $\tau$ with a prime and
the derivative with respect to the cosmic time $t$ with a dot, e.g.
$a^{\prime}=da/d\tau$ and $\dot{a}=da/dt$. The Hubble parameter is defined as
$H=\dot{a}(t)/a(t)$ where $a$ is the scale factor. We assume a homogeneous
inflaton with only zero mode, $\phi(t,\bm{x})=\phi(t)$.

We now quantize the gauge field $\bm{A}$ in momentum space

$\bm{A}(\tau,\bm{x})\,=\,\sum_{\lambda=\pm}\int\frac{d^{3}k}{(2\pi)^{3}}\,\left[\bm{\epsilon}_{\lambda}(\bm{k})\,a_{\lambda}(\bm{k})\,A_{\lambda}(\tau,\bm{k})\,e^{i\bm{k}\cdot\bm{x}}+\,\text{h.c.}\right],$
(2.4)

where $\lambda$ is the photon polarization, and $a_{\lambda}(\bm{k})$
($a_{\lambda}^{\dagger}(\bm{k})$) are annihilation (creation) operators that
fulfill the canonical commutation relations

$[a_{\lambda}(\bm{k}),a_{\lambda^{\prime}}^{\dagger}(\bm{k}^{\prime})]=(2\pi)^{3}\delta_{\lambda\lambda^{\prime}}\delta^{(3)}(\bm{k}-\bm{k^{\prime}})\,.$
(2.5)

The polarization vectors $\bm{\epsilon}_{\lambda}(\bm{k})$ satisfy the
conditions333A simple realization can be given in terms of a real basis with
the orthonormal vectors $(\bm{k}/|\bm{k}|,\bm{e}_{i})$, ($i=1,2$), such that
$\bm{k}\cdot\bm{e}_{i}=\bm{e}_{1}\cdot\bm{e}_{2}=0$ and
$\bm{e}_{i}\cdot\bm{e}_{i}=1$, with
$\bm{\epsilon}_{\lambda}\equiv(\bm{e}_{1}+i\lambda\bm{e}_{2})/\sqrt{2}$, from
where identities (2.6) follow.

$\displaystyle\bm{k}\cdot\bm{\epsilon}_{\lambda}(\bm{k})$
$\displaystyle=0\,,\hskip 56.9055pt$
$\displaystyle\bm{k}\times\bm{\epsilon}_{\lambda}(\bm{k})$
$\displaystyle=-i\lambda k\,\bm{\epsilon}_{\lambda}(\bm{k})\,,$ (2.6)
$\displaystyle\bm{\epsilon}^{*}_{\lambda^{\prime}}(\bm{k})\cdot\bm{\epsilon}_{\lambda}(\bm{k})$
$\displaystyle=\delta_{\lambda\lambda^{\prime}}\,,$
$\displaystyle\bm{\epsilon}^{*}_{\lambda}(\bm{k})$
$\displaystyle=\bm{\epsilon}_{\lambda}(-\bm{k})\,,$

where $k\equiv|\bm{k}|$. Therefore, the equation of motion for the gauge modes
yields

$A^{\prime\prime}_{\lambda}+\sigma
A^{\prime}_{\lambda}+k\left(k-\lambda\,\dfrac{a\,\dot{\phi}}{f_{\phi}}\right)A_{\lambda}=0\,.$
(2.7)

In some special cases ($\sigma=0$ and slow-roll inflation), this equation can
be solved analytically, and we will do it in Sec. 3.1. In the general case, we
will solve it using numerical methods.

###  2.2 Observable quantities

Once we obtain a solution to the modes $A_{\lambda}$, we can compute the
(hyper)electromagnetic energy densities as

$\displaystyle\rho_{E}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{4}}\int^{k_{c}}_{k_{\rm
min}}dk\,\frac{k^{2}}{4\pi^{2}}\left(|A^{\prime}_{+}|^{2}+|A^{\prime}_{-}|^{2}\right),$
(2.8a) $\displaystyle\rho_{B}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{4}}\int^{k_{c}}_{k_{\rm
min}}dk\,\frac{k^{4}}{4\pi^{2}}\left(|A_{+}|^{2}+|A_{-}|^{2}\right).$ (2.8b)

The upper integration limit comes because subhorizon modes have an oscillatory
behavior and should be regarded as quantum fluctuations. Therefore, such modes
do not contribute to the above classical observables and are excluded from the
integration. More details and precise value of $k_{c}$ will be given in Sec.
2.3. For the lower integration limit $k_{\rm min}$, see Eq. (A.5).

In this work, we will also make use of the (hyper)magnetic helicity and its
derivative, defined as

$\displaystyle\mathcal{H}$ $\displaystyle\equiv$
$\displaystyle\lim_{V\to\infty}\dfrac{1}{V}\int_{V}d^{3}x\;\frac{\langle\bm{A}\cdot\bm{B}\rangle}{a^{3}}=\frac{1}{a^{3}}\int^{k_{c}}_{k_{\rm
min}}dk\,\frac{k^{3}}{2\pi^{2}}\left(|A_{+}|^{2}-|A_{-}|^{2}\right),$ (2.9a)
$\displaystyle\mathcal{G}$ $\displaystyle\equiv$
$\displaystyle\dfrac{1}{2a}\dfrac{d\mathcal{H}}{d\tau}=-\lim_{V\to\infty}\dfrac{1}{V}\int_{V}d^{3}x\;\frac{\langle\bm{E}\cdot\bm{B}\rangle}{a^{4}}.$
(2.9b)

In the case of one Dirac fermion with mass $m$ and hypercharge $Q_{Y}$, the
conductivity can be written as444As the conductivity $\sigma$ relates $\bm{J}$
and $\bm{E}$ in (2.3), it is a comoving quantity, i.e. it scales with the
Universe expansion. Our definition differs from the one in Gorbar:2021rlt ;
Gorbar:2021zlr  where the authors used a physical conductivity that we will
denote $\hat{\sigma}$ in this paper, their relation being
$\sigma=a\,\hat{\sigma}$. Domcke:2018eki

$\sigma=\frac{|g^{\prime}Q_{Y}|^{3}}{6\pi^{2}}\frac{a}{H}\;\sqrt{2\rho_{B}}\;\coth\left(\pi\sqrt{\frac{\rho_{B}}{\rho_{E}}}\right)\exp{\left\\{-\frac{\pi
m^{2}}{\sqrt{2\rho_{E}}\,|g^{\prime}Q_{Y}|}\right\\}},$ (2.10)

where $g^{\prime}\simeq 0.4$ is computed at the characteristic scale
$\mu\simeq(\langle\bm{E}\rangle^{2}+\langle\bm{B}\rangle^{2})^{1/4}$ where the
Schwinger effect takes place Gorbar:2021zlr . This estimation is valid in the
case of collinear electric and magnetic fields, an assumption that we have
numerically checked by verifying that

$\cos\theta\equiv\frac{|\mathcal{G}|}{2\sqrt{\rho_{E}\rho_{B}}}\simeq 1,$
(2.11)

where $\theta$ is the spatial angle between $\bm{E}$ and $\bm{B}$.

Moreover, the massless hypercharged fermions that are continuously produced
during inflation have an energy density given by

$\rho_{\psi}=\lim_{V\to\infty}\dfrac{\sigma}{V}\int_{V}d^{3}x\;\frac{\langle\bm{A}\cdot\bm{E}\rangle}{a^{4}}=\frac{\sigma}{a^{4}}\int^{k_{c}}_{k_{\rm
min}}dk\,\frac{k^{2}}{2\pi^{2}}\dfrac{d}{d\tau}\left(|A_{+}|^{2}+|A_{-}|^{2}\right)$
(2.12)

Notice that the observable quantities $\rho_{E}$, $\rho_{B}$, $\rho_{\psi}$,
$\mathcal{H}$ and $\mathcal{G}$ are physical555They relate to the comoving
ones $\rho_{E}^{c}$, $\rho_{B}^{c}$, $\mathcal{H}^{c}$, and $\mathcal{G}^{c}$
by the relations $\rho_{B,E}^{c}=a^{4}\rho_{B,E}$,
$\mathcal{H}^{c}=a^{3}\mathcal{H}$, $\mathcal{G}^{c}=a^{4}\mathcal{G}$., while
the fields $\bm{A}$, $\bm{E}$ and $\bm{B}$ as well as the conductivity
$\sigma$ and current $\bm{J}$ are comoving.

Concerning the Higgs vacuum expectation value, there are two possibilities
during the inflationary period:

i)

The first possibility, which we will consider throughout this paper, is that
$\langle h\rangle=0$, and so the electroweak symmetry is unbroken during the
inflationary period. In order to ensure unbroken electroweak symmetry and
hence massless SM fermions, which all contribute to the conductivity (2.10),
we assume that the SM Higgs field $h$ remains stabilized at the origin in
field space by a large mass term throughout the inflationary period. Such a
large mass can, e.g., be induced by a nonminimal coupling to the Ricci
curvature scalar as $\mathcal{L}=\frac{1}{2}\xi h^{2}R$ with $\xi>3/16$ (see
e.g. Ref. Espinosa:2015qea ). Hence, we get

$\sigma\simeq\frac{41\,g^{\prime
3}}{72\,\pi^{2}}\frac{a}{H}\;\sqrt{2\rho_{B}}\;\coth\left(\pi\sqrt{\frac{\rho_{B}}{\rho_{E}}}\right).$
(2.13)

ii)

The second possibility is that the electroweak symmetry is broken during the
inflationary period. In this case after $\Delta N$ $e$-folds of inflation,
there is a Gaussian distribution of values of the Higgs field with zero mean
and variance $\langle h^{2}\rangle=H^{2}\Delta N/(4\pi^{2})$ with probability
$P(h,\Delta N)\propto\exp(-\frac{1}{2}\frac{h^{2}}{\langle h^{2}\rangle})$
dominated by the values $h\lesssim\sqrt{\langle h^{2}\rangle}$, see Ref.
Espinosa:2015qea 666The SM Higgs potential is still unstable at a value of the
Higgs field $h=h_{I}\simeq 10^{11}$ GeV and the condition for $P(h_{I},\Delta
N)<e^{-3\Delta N}$ (so that it is unlikely to find the Higgs away from its EW
vacuum in any of the $e^{3\Delta N}$ causally disconnected regions formed
during inflation) implies $H_{E}<\sqrt{2/3}\,\pi h_{I}/\Delta N$, a condition
that is not fulfilled by any of the models of inflation we have considered.
Therefore this possibility would require stabilization of the Higgs potential
by some new physics.. In this case, the electroweak symmetry is broken and the
hypercharge field strength $Y_{\mu\nu}$ in Eq. (2.1) is projected onto the
electromagnetic field strength $F_{\mu\nu}$ with a coupling to the inflaton
given by $f_{\phi}/\cos^{2}\theta_{W}$ where $\theta_{W}$ is the electroweak
angle. Now the conductivity for the hypermagnetic field in Eq. (2.10) should
be replaced by a similar expression for the magnetic field, with the
replacement $|g^{\prime}Q_{Y}|\to|eQ|$, where
$e=gg^{\prime}/\sqrt{g^{2}+g^{\prime 2}}$ and $Q$ is the fermion electric
charge. The condition for a fermion $f$ to contribute to the magnetic
conductivity $\pi m_{f}^{2}<\sqrt{2\rho_{E}}\,|eQ_{f}|$ translates into the
condition, for the fermion Yukawa coupling,

$Y_{f}\lesssim 0.45\left(\frac{\rho_{E}}{H^{4}}\right)^{1/4}\sqrt{|Q_{f}|},$
(2.14)

and we have computed all couplings at the characteristic scale
$\mu\simeq(\langle\bm{E}\rangle^{2}+\langle\bm{B}\rangle^{2})^{1/4}$ where the
Schwinger effect takes place. If the three generations of fermions satisfy the
above condition then the conductivity for the magnetic field is given by Eq.
(2.13) with the replacement $\frac{41g^{\prime
3}}{72\pi^{2}}\to\frac{e^{3}}{\pi^{2}}$. We have checked that, in this case,
the results for $f_{\phi}\lesssim 0.2\,M_{\text{pl}}$ are consistent with all
three generation fermions contributing to the magnetic conductivity. For
$f_{\phi}\gtrsim 0.2\,M_{\text{pl}}$ only the top quark does not contribute.
Given that $41g^{\prime 3}/72\simeq 0.37$ while $e^{3}\simeq 0.36$, at the
scales where the Schwinger effect takes place, we have found that the results
in this second case are qualitatively similar to those for the previous case,
which will be worked out in detail in this paper.

Considering then the case i) above, the conductivity (2.13) yields a
nontrivial integro-differential system as the damping term grows with the
magnetic energy and hence backreacts on the amount of produced
electric/magnetic fields. We aim to solve this setup of the Schwinger effect
numerically. In the next sections we will consider specific cases where this
system can be further simplified.

###  2.3 The gauge vacuum

At very early times, when $|a\dot{\phi}|\ll kf_{\phi}$, the modes are in their
Bunch-Davies (BD) vacuum, hence

$A_{\lambda}(\tau,k)=\dfrac{1}{\sqrt{2k}}\;{\rm e}^{-ik\tau}\hskip
42.67912pt(\tau\to-\infty).$ (2.15)

Initially, we can consider all the modes in the BD vacuum (which would be
possible by initializing the numerical simulation such that $a_{0}\ll
k_{0}/H_{0}$). In that case, since $|A_{+}|=|A_{-}|$, the fields $\bm{E}$ and
$\bm{B}$ are plane waves perpendicular to each other, as $\mathcal{G}=0$ in
(2.11) yields $\cos\theta=0$. Therefore, there is no Schwinger effect and
$\sigma=0$.

It has recently been shown that in the presence of the conductivity $\sigma$,
the BD vacuum amplitude of the modes that are still in the vacuum get damped
by the ones that left it Gorbar:2021rlt . Indeed, consider we are at a time
$a_{\ast}$ where modes $k>k_{\ast}$ are still in the BD vacuum, while modes
$k<k_{\ast}$ were amplified by both tachyonic and parametric instabilities
from Eq. (2.7). Then, the equation of motion for modes such that
$|a_{\ast}\dot{\phi}(\tau_{\ast})|\ll kf_{\phi}$ does not reduce to a plane
wave in the presence of a non-zero $\sigma$, but instead to
$A^{\prime\prime}_{\lambda}+\sigma A^{\prime}_{\lambda}+k^{2}A_{\lambda}=0$,
and Eq. (2.15) is not a solution anymore. To derive the generalized BD vacuum,
we write the gauge equation of motion (2.7) in cosmic time:

$\ddot{A_{\lambda}}+\left(\hat{\sigma}+H\right)\dot{A_{\lambda}}+\frac{k}{a}\left(\frac{k}{a}-\lambda\frac{\dot{\phi}}{f_{\phi}}\right)A_{\lambda}=0,$
(2.16)

where we used the identity
$a^{-2}A^{\prime\prime}_{\lambda}=\ddot{A}_{\lambda}+H\dot{A}_{\lambda}$, and
perform the transformation $A_{\lambda}=\sqrt{\Delta}\,\mathcal{A}_{\lambda}$
with Gorbar:2021rlt

$\Delta(t)=\exp{\left\\{-\int_{-\infty}^{t}\hat{\sigma}(t^{\prime})\,dt^{\prime}\right\\}}.$
(2.17)

We recall that we have defined $\hat{\sigma}=\sigma/a$ as the physical
conductivity in footnote 4. The above equation hence becomes

$\mathcal{A}^{\prime\prime}_{\lambda}+\left[\frac{k}{a}\left(\frac{k}{a}-\lambda\frac{\dot{\phi}}{f_{\phi}}\right)-\frac{\dot{\hat{\sigma}}}{2}-\frac{\hat{\sigma}^{2}}{4}-\frac{H\hat{\sigma}}{2}\right]a^{2}\mathcal{A}=0,$
(2.18)

where we used the fact that $\dot{\Delta}(t)=-\hat{\sigma}(t)\Delta(t)$. A
mode crosses the horizon when the expression in the square brackets vanishes
for the first time at least for one polarization, at $k=k_{c}$. The modes in
the vacuum are then characterized by $k\gg k_{c}$. This yields the momentum of
the mode that crosses the horizon at time $t$, namely the cutoff of the
integrals:

$k_{c}=\left|\frac{a\dot{\phi}}{2f_{\phi}}\right|+\sqrt{\left(\frac{a\dot{\phi}}{2f_{\phi}}\right)^{2}+\frac{a^{2}}{2}\left[\dot{\hat{\sigma}}+\hat{\sigma}\left(\frac{\hat{\sigma}}{2}+H\right)\right]}\,.$
(2.19)

Deep inside the horizon, when the first term in square brackets of (2.18)
dominates, the solution must satisfy the BD condition (2.15). As we have seen,
in the presence of finite conductivity, this equation does not fully describe
the gauge-field mode function inside the horizon, as the damped BD condition
includes an exponential damping factor

$A_{\lambda}(\tau,k)=\sqrt{\dfrac{\Delta(t)}{2k}}e^{-ik\tau}\hskip
42.67912pt(\tau\to-\infty).$ (2.20)

The bottom line of this section is that the modes still in their BD vacuum see
their amplitudes damped because of the effect of the modes that left their
vacuum earlier and participate in the equations of motion (2.2a) and (2.7).
The parameter $\Delta$ was first introduced in the context of the gradient
expansion formalism in Ref. Gorbar:2021rlt , where it was dynamically solved,
while in Ref. Gorbar:2021zlr  it was also considered as a free parameter and
validated the corresponding procedure by numerical calculations. In order to
compare with results from the gradient expansion formalism in configuration
space, we will also both compute $\Delta$ numerically and consider it as a
free parameter, although our final results will be based upon the dynamical
calculation of $\Delta$.

##  3 Slow Roll Analysis

The slow-roll inflation paradigm has been used by many authors to compute the
amount of electromagnetic energy density Anber:2006xt ; Anber:2009ua ;
Durrer:2013pga , or baryogenesis through helicity Anber:2015yca ; Cado:2016kdp
; Domcke:2019mnd ; Cado:2022evn  at the end of inflation, with or without
taking into account the Schwinger effect. Here we aim to validate our
numerical results by comparison with the known analytical results at the end
of inflation.

In this section we will take $\phi$ as a slowly rolling inflaton field such
that $\ddot{\phi}\simeq 0$, $3H\dot{\phi}\simeq-V^{\prime}(\phi)$, and so we
can consider $\dot{\phi}$ and $H=H_{E}$ as constant. Doing so, we are
neglecting the gauge field backreaction in the right-hand side of Eq. (2.2a),
a hypothesis that we have consistently checked a posteriori. The results of
this section will be model independent, within the hypothesis of the slow roll
approximation.

###  3.1 Absence of Schwinger effect

Here we are assuming there is no Schwinger effect777This condition should be
considered as being fulfilled by some physical systems, as e.g. systems with
no massless fermions, more than as an approximation to the full (more
realistic) case., i.e. $\sigma=0$, hence we can rewrite (2.7) as

$A^{\prime\prime}_{\lambda}+k\left(k+\lambda\,\dfrac{2\xi}{\tau}\right)A_{\lambda}=0,$
(3.1)

where, following the slow roll equations,

$\xi=-\dfrac{\dot{\phi}}{2H_{E}f_{\phi}}$ (3.2)

is a constant. Since we are in de Sitter space, we can use the scale factor
definition $a=-(H\tau)^{-1}$ and solve (3.1) asymptotically. At early time,
when $|k\tau|\gg 2\xi$, the modes are in their BD vacuum given by (2.15), as
here $\Delta=1$. When $|k\tau|\sim 2\xi$, one of the modes develops both
parametric and tachyonic instabilities leading to exponential growth while the
other stays in the vacuum. During the last $e$-folds of inflation, i.e.
$|k\tau|\ll 2\xi$, the growing mode has the solution Anber:2006xt ;
Anber:2015yca

$A_{\lambda}\simeq\dfrac{1}{\sqrt{2k}}\left(\dfrac{k}{2\xi
a_{E}H_{E}}\right)^{\frac{1}{4}}\exp{\left\\{\pi\xi-2\sqrt{\dfrac{2\xi
k}{a_{E}H_{E}}}\right\\}},$ (3.3)

where $a_{E}$ and $H_{E}$ are, respectively, the scale factor and the Hubble
parameter at the end of inflation. Here, as we assume a slow roll regime, we
consider $H_{E}$ constant and we take the convention $a_{E}=1$.

Using (2.8) and (2.9) we can compute all electromagnetic quantities:

$\rho_{E}\simeq\dfrac{63}{2^{16}}\;\dfrac{H_{E}^{4}}{\pi^{2}\xi^{3}}\;{\rm
e}^{2\pi\xi},\quad\rho_{B}\simeq\dfrac{315}{2^{18}}\;\dfrac{H_{E}^{4}}{\pi^{2}\xi^{5}}\;{\rm
e}^{2\pi\xi},\quad\mathcal{H}\simeq\dfrac{45}{2^{15}}\;\dfrac{H_{E}^{3}}{\pi^{2}\xi^{4}}\;{\rm
e}^{2\pi\xi},\quad\mathcal{G}\simeq\dfrac{135}{2^{16}}\dfrac{H_{E}^{4}}{\pi^{2}\xi^{4}}\;{\rm
e}^{2\pi\xi}.$ (3.4)

These results are only valid when the absence of backreaction on the inflaton
equation of motion (2.2a) is guaranteed, hence when
$\left|\mathcal{G}/V^{\prime}(\phi)\right|\ll f_{\phi}$. This model-dependent
condition puts a lower bound on the parameter $f_{\phi}$ or, equivalently, a
higher bound on $\xi$. Using the slow roll equations and the definition of the
slow roll parameters, this parameter can be written as

$\xi=\dfrac{M_{\text{pl}}}{f_{\phi}}\sqrt{\dfrac{\epsilon}{2}},$ (3.5)

where $\epsilon=(M_{\text{pl}}^{2}/2)\,(V^{\prime}/V)^{2}$. Therefore, at the
end of inflation, where by definition $\epsilon=1$, one has
$\xi=M_{\text{pl}}/\sqrt{2}f_{\phi}$, and the no backreaction condition in Eq.
(2.2) provides the bound $\xi<5.73$ (or equivalently
$f_{\phi}>0.12\,M_{\text{pl}}$). In Fig. 1 we show, with orange lines, the
quantities $\rho_{E}$, $\rho_{B}$, $\mathcal{H}$ and $\mathcal{G}$ evaluated
at the end of inflation obtained from the analytical backreactionless
solutions from Eqs. (3.4), while the blue dots are the numerical solutions,
which correspond to the case $\sigma=0$ (no Schwinger effect) and
correspondingly $\Delta=1$. We have used a Runge-Kutta method which is
explained in App. A.

Figure 1: Electric $\rho_{E}$ and magnetic $\rho_{B}$ energy densities, and
the helicity $\mathcal{H}$ and its derivative $\mathcal{G}$, at the end of
inflation (i.e. for $\epsilon(a_{E})=1$), in units of $H_{E}$, as functions of
the coupling $f_{\phi}$ assuming $\Delta$ constant. We see the plots confirm
the result from Fig. 1 of Gorbar:2021zlr . Here, we also assumed $\xi$
constant.

###  3.2 Presence of Schwinger effect

The Schwinger effect is taken into account by means of the conductivity
$\sigma$ in Eq. (2.7), as given by Eq. (2.10) Domcke:2018eki . The growth of
$\sigma$ with time then backreacts on the gauge field, as the damping term
grows in its differential equation. We will compare our numerical calculations
with three analytical (or semianalytical) results: the Schwinger maximal and
equilibrium estimates Domcke:2018eki ; Gorbar:2021zlr , as well as the
gradient expansion formalism Sobol:2019xls ; Gorbar:2021rlt ; Gorbar:2021zlr .
From the numerical point of view however, we aim to solve Eq. (2.7) with
$\sigma$ computed at each time step using (2.13). The details about the
numerics will be displayed in Sec. 3.3.

#### Schwinger equilibrium estimate

In this case, the backreaction of the chiral fermions on the gauge fields is
taken into account by just replacing the parameter $\xi$ with the effective
one Domcke:2018eki

$\xi_{\rm eff}=\xi-\frac{41\,g^{\prime
3}}{144\,\pi^{2}}\coth\left(\pi\sqrt{\frac{\rho_{B}}{\rho_{E}}}\right)\frac{\sqrt{2\rho_{E}}}{H_{E}^{2}}\,,$
(3.6)

in the backreactionless solutions (3.4). This amounts to solving

$\,\dfrac{63}{2^{15}\pi^{2}}\;\dfrac{{\rm e}^{2\pi\xi_{\rm eq}}}{\xi_{\rm
eq}^{3}}=\left(\frac{144\,\pi^{2}}{41\,g^{\prime 3}}\right)^{2}(\xi-\xi_{\rm
eq})^{2}\;\tanh^{2}\left(\sqrt{\frac{5}{4}}\,\frac{\pi}{\xi_{\rm eq}}\right),$
(3.7)

which provides the function $\xi_{\rm eq}=\xi_{\rm eq}(\xi)$ that we plug in
(3.4) instead of the bare $\xi$ to obtain the quantities $\rho_{\rm eq}^{E}$,
$\rho_{\rm eq}^{B}$, $\mathcal{H}_{\rm eq}$ and $\mathcal{G}_{\rm eq}$. These
equilibrium estimates are shown with a purple line in the plots of Fig. 1.

#### Schwinger maximal estimate

In this case, we assume the exponential behaviors of the backreactionless
solutions to be valid until they saturate the maximal value that we will
display hereafter. We numerically determine the value of crossing, which
happens for $\xi\simeq 4.4$-$4.7$ depending on each quantity.

The maximum helicity density can be estimated as the solution of
Domcke:2018eki

$|E|^{2}+|B|^{2}=\xi_{\rm eff}|E|\,|B|.$ (3.8)

This replacement yields an equation relating the $|E|$ and $|B|$ fields that
can be solved analytically. We then choose, as definition of our maximal
estimate, the solution $(|E|,|B|)$ of (3.8) that maximizes the product
$|E|\cdot|B|$ 888Notice that our definition of maximal solution departs from
that given in Refs. Domcke:2018eki ; Gorbar:2021zlr , where the fields $|E|$
and $|B|$ are separately maximized, while we are maximizing the product
$|E|\cdot|B|$, the relevant quantity for the baryon asymmetry generation..
This yields for $\xi\gg 1$

$\displaystyle\rho^{E}_{\rm max}$ $\displaystyle\simeq$
$\displaystyle\frac{8}{9}\left(\frac{72\,\pi^{2}}{41\,g^{\prime
3}}\right)^{2}\;\xi^{2}H_{E}^{4},$ (3.9a) $\displaystyle\rho^{B}_{\rm max}$
$\displaystyle\simeq$
$\displaystyle\frac{8}{81}\left(\frac{72\,\pi^{2}}{41\,g^{\prime
3}}\right)^{2}\;\xi^{4}H_{E}^{4},$ (3.9b) $\displaystyle\mathcal{H}_{\rm max}$
$\displaystyle=$ $\displaystyle\frac{2\,\mathcal{G}_{\rm
max}}{3H_{E}}\simeq\frac{32}{81}\left(\frac{72\,\pi^{2}}{41\,g^{\prime
3}}\right)^{2}\;\xi^{3}H_{E}^{3}.$ (3.9c)

The maximal estimates for the quantities $\rho^{E}_{\rm max}$, $\rho^{B}_{\rm
max}$, $\mathcal{H}_{\rm max}$ and $\mathcal{G}_{\rm max}$ are shown with a
pink line in the plots of Fig. 1.

#### Gradient expansion formalism

This method was introduced in Refs. Sobol:2019xls ; Gorbar:2021rlt ;
Gorbar:2021zlr  and transforms the EoM for the vector field $\bm{A}$ into EoM
for observable quantities, in particular the electric $\bm{E}$ and magnetic
$\bm{B}$ fields. As the spatial gradients in the EoM do always appear as ${\rm
rot}\,\bm{E}$ and ${\rm rot}\,\bm{B}$, the EoM can be written as an infinite
series in terms of the bilinears
$\mathcal{E}^{(n)}=\langle\bm{E}\cdot\textrm{rot}^{n}\bm{E}\rangle/a^{n}$,
$\mathcal{G}^{(n)}=\langle\bm{E}\cdot\textrm{rot}^{n}\bm{B}\rangle/a^{n}$ and
$\mathcal{B}^{(n)}=\langle\bm{B}\cdot\textrm{rot}^{n}\bm{B}\rangle/a^{n}$,
with $n=0,1,\dots$. In this way the coupled system of EoM for the fields
$\bm{E}$ and $\bm{B}$ transforms into a system of coupled differential
equations for the quantities $\mathcal{E}^{(n)}$, $\mathcal{B}^{(n)}$ and
$\mathcal{G}^{(n)}$. This system is not block diagonal in the space of the $n$
index so that the system has to be truncated to find solutions.

Moreover, the parameter $\Delta(t)$ in Eq. (2.17), which suppresses the gauge-
field amplitude on small scales depends on the conductivity at all times
$t^{\prime}<t$. So, a precise determination of $\Delta(t)$ would require a
complete analytical solution of the infinite-dimensional system of equations.
While $\Delta$ was dynamically computed in Ref. Gorbar:2021rlt , for the sake
of simplicity and generality, it was considered as a free parameter in Ref.
Gorbar:2021zlr  and fixed to the values
$\Delta=1,\,10^{-2},\,10^{-4},\,10^{-6}$. In our numerical approach we will
consider $\Delta$ as a function of the conductivity $\sigma$, as the initial
condition for $\bm{E}$ and $\bm{B}$ are plane waves, such that
$\bm{E}\cdot\bm{B}=0$ and therefore initially $\sigma=0$ and so $\Delta=1$.
However, as time is evolving $\bm{E}$ and $\bm{B}$ will become collinear, and
a nonvanishing conductivity will appear, as well as the function
$\Delta(t)<1$. In order to compare our numerical results with those from Ref.
Gorbar:2021zlr , we also will eventually enforce $\Delta$ to be a constant in
our code. Upon considering a constant value of $\Delta$, our results will
agree pretty well with those obtained in the gradient expansion formalism, see
Fig. 1. In the more realistic cases where we just compute the value of
$\Delta(t)$, we will see that at the beginning $t=t_{0}$, just very deep
inside the inflationary period, $\Delta(t_{0})=1$, while the value of $\Delta$
will decrease very fast and at the end of inflation $t=t_{E}$,
$\Delta(t_{E})\ll 1$.

###  3.3 Numerical results at the end of inflation

We will find it more convenient to change the variable from the time $t$ to
the scale factor $a$. The gauge field equation of motion (2.7) then becomes

$\frac{\partial^{2}A_{\lambda}}{\partial
a^{2}}+\frac{1}{a}\left(2+\frac{\sigma}{aH_{E}}\right)\frac{\partial
A_{\lambda}}{\partial
a}+\frac{k}{a^{3}H_{E}}\left(\frac{k}{aH_{E}}-2\lambda\xi\right)A_{\lambda}=0.$
(3.10)

We recall that, as we are considering the slow roll regime in this section, we
do not need to solve the equation of motion for $\phi$.

The Bunch-Davis solutions can now be written as

$\begin{aligned} A_{\lambda}(a,k)&=\sqrt{\frac{\Delta(a)}{2k}}\;{\rm
e}^{ik/aH_{E}}\\\ \frac{\partial A_{\lambda}}{\partial
a}(a,k)&=\frac{\sqrt{\Delta(a)}}{a^{2}H_{E}}\left(-i\,\sqrt{\dfrac{k}{2}}-\frac{\sigma}{2}\dfrac{1}{\sqrt{2k}}\right){\rm
e}^{ik/aH_{E}}\end{aligned}\hskip 34.14322pt(a\to 0),$ (3.11)

with

$\Delta(a)=\exp{\left\\{-\int_{a_{0}}^{a}\frac{\sigma(a^{\prime})}{a^{\prime
2}H_{E}}da^{\prime}\right\\}}.$ (3.12)

The technical details of the numerical simulations for solving Eq. (3.10),
subject to the boundary conditions (3.11), can be found in App. A. We display
in Fig. 2 the spectra of all the observable quantities in order to see how the
BD vacuum is dominating the spectra for large $k$ and how the cutoff
$k_{c}(a)$, given by (2.19), efficiently removes that part of the integration.
The difference between the BD vacuum and the damped BD vacuum is also clear,
as the first goes like $k^{3}$ whereas the second goes like $\Delta(a)k^{3}$
with $\Delta$ decreasing with time. Hence the asymptotic behaviors are not
superimposed since $\Delta$ changes. Finally, we also see explicitly how the
growth of $\rho_{E}$ and $\rho_{B}$ with the scale factor $a$ is due to the
increase in amplitude of the spectrum hump and its shift to larger values of
$k$. For this illustrative purpose we used a constant value of $\xi$. Here we
have fixed $f_{\phi}=0.1\,M_{\text{pl}}$, while for other values of this
parameter the plots are similar.

Before moving to the full numerical results, we will compare our slow roll
based inflaton numerical results with the recent literature on the subject.

Figure 2: Spectra of the magnetic energy (top left), electric energy (bottom
left), helicity (top right) and its derivative (bottom right), i.e. the
integrands of (2.8) and (2.9), for different values of $a$ during inflation
simulation. Here we used variable $\sigma(a)$ and $\Delta(a)$ with constant
$\xi$. The color matching dashed vertical lines show the cutoff values
$k_{c}(a)$ computed from (2.19).

####  3.3.1 Constant $\Delta$ and $\xi$ approximation

We will first assume that the parameters $\Delta$ and $\xi$ are constants. As
we already mentioned, the parameter $\Delta$ was fixed to constant values in
Ref. Gorbar:2021zlr  while $\xi$, as defined in Eq. (3.2), is often considered
as a constant in the slow roll approximation. In Fig. 1, we displayed several
results already present in the literature that we successfully reproduced with
our numerical method. First the backreactionless case, where there is no
conductivity, by simply enforcing $\sigma=0$ (therefore $\Delta=1$) in the
code. The data set are displayed in blue and match the corresponding
analytical value given by Eq. (3.4). Then, in order to reproduce results from
Gorbar:2021zlr , we considered a non-zero conductivity given by (2.13) while
assuming $\Delta$ constant during inflation, thus making it a free parameter.
In Fig. 1, we plot the quantities $\rho_{B}$, $\rho_{E}$, $\mathcal{H}$ and
$\mathcal{G}$ at the end of inflation for chosen values of $\Delta$. We can
see that the results agree well with those using the gradient expansion
formalism in Ref. Gorbar:2021zlr .

####  3.3.2 Variable $\Delta$ and $\xi$

The benefit of the slow roll approximation is that the results look “model
independent”. However, the tradeoff comes with the need of having a constant
parameter $\xi$ as the slow roll regime implies an approximately constant
$\dot{\phi}$. Besides, we know that this parameter can also be expressed in
terms of the slow roll parameter $\epsilon$ (see Eq. (3.5)), which is indeed
small and constant during inflation but then quickly becomes unity during the
last $e$-folds. We also know that the modes produced during the last $e$-folds
are the ones that contribute the most to the integrals (2.8) and (2.9), as all
the modes previously generated get washed out by the Universe expansion.

All these observations lead us to conclude that the most important
contribution to the quantities $\rho_{E}$, $\rho_{B}$, $\mathcal{H}$ and
$\mathcal{G}$ is taking place during an epoch when the constant $\xi$
approximation loses its relevance. Hence, in this section, we will instead
specify an inflation model, namely the Starobinsky potential, and make its
study in the slow roll regime with a function $\xi(a)$ that can be obtained
from the model. We have chosen in this section the Starobinsky potential as it
provides a realistic model of inflation, and will be a particular case of a
more general class of models we will consider to make predictions using the
full solution of the system. The purpose of this section will thus be to
assess the goodness of the slow roll approximation when computing the full
solution to the system (2.2).

The Starobinsky potential is given by

$V=\Lambda^{4}\left[1-\exp{\left\\{-\sqrt{\dfrac{2}{3}}\dfrac{|\phi|}{M_{\text{pl}}}\right\\}}\right]^{2}.$
(3.13)

Using the slow roll regime, the inflaton field $\phi$ is given by

$\sqrt{\frac{2}{3}}\,\frac{\phi(a)}{M_{\text{pl}}}=-\log\left(\dfrac{a_{E}}{a}\right)^{\frac{4}{3}}-\mathcal{W}_{-1}\left[-\beta
e^{-\beta}\left(\dfrac{a_{E}}{a}\right)^{-\frac{4}{3}}\right]-\beta+\log\beta,\quad\beta=1+\dfrac{2}{\sqrt{3}}$
(3.14)

where $\mathcal{W}_{n}$ is the $n$th branch of the Lambert function. The value
of the function $\xi$ is then given by

$\xi(a)=\sqrt{\frac{2}{3}}\frac{M_{\text{pl}}}{f_{\phi}}\frac{1}{\exp\left[\sqrt{\frac{2}{3}}\frac{\phi(a)}{M_{\text{pl}}}\right]-1}\,.$
(3.15)

In Fig. 3, we display in blue results for the Starobinsky model, for various
values of $\epsilon$, when $\sigma$ and $\Delta$ vary dynamically. Although
the slow roll approximation loses its relevance for values of $\epsilon$
closer to 1 (an issue we address in the next section), we already see a
difference with the plots in Fig. 1. This is because, no matter the value of
the initial time, the function $\Delta(a)$ rapidly goes to extremely small
values, thus killing the BD modes that would have been amplified at the very
end of inflation and that would have contributed the most to the integrals
(2.8) and (2.9). With a constant $\Delta$, this suppression is less effective
and the tachyonic amplification yields higher energy densities and helicity.

Figure 3: Comparison between the slow roll approximation and the full solution
for the Starobinsky model. The analytical estimates are given for
$\epsilon=1$. As expected, the slow roll computation diverges from the full
solution as inflation is nearing the end, since the slow roll approximation is
only valid in the regime $\epsilon\ll 1$. Hence the slow roll computation
overshoots the value of all quantities, closer to the value given by the
Schwinger equilibrium estimate for $f_{\phi}\lesssim 0.05\,M_{\text{pl}}$. As
expected, we also have compared both analysis, slow roll and full solution,
for values of $a$ such that $\epsilon(a)\ll 1$ (in particular
$\epsilon=10^{-1},\,10^{-2},\,10^{-3}$) and found good agreement.

##  4 Full Analysis

In this section, we are not using the slow roll hypothesis for the inflaton
equation of motion and consider the full solution to the system (2.2) in
specific models of inflation. We will choose a set of inflationary models that
are well known to be in agreement with all cosmological constraints. Also, we
do not assume any peculiar geometry of the Universe.

The equations to be solved during inflation are the system (2.2) written in
terms of the variable $a$. Unlike in the previous section, the current change
of variables must take into account that the Hubble parameter is not constant,
but moreover we have $\frac{da}{dt}=\dot{a}=aH$, and we will define the
auxiliary quantity $\mathcal{F}$ as

$\mathcal{F}=-\frac{a}{H}\frac{dH}{da}=-\frac{a}{2H^{2}}\frac{dH^{2}}{da}.$
(4.1)

We will relate it to the Friedmann equations

$\displaystyle H^{2}$ $\displaystyle=$
$\displaystyle\frac{\rho}{3M_{\text{pl}}^{2}},$ (4.2a)
$\displaystyle\frac{\ddot{a}}{a}$ $\displaystyle=$
$\displaystyle-\frac{3p+\rho}{6M_{p}^{2}},$ (4.2b)

which combine themselves into

$\frac{a}{2}\frac{dH^{2}}{da}=\frac{dH}{dt}=\frac{\ddot{a}}{a}-H^{2}=-\frac{p+\rho}{2M_{p}^{2}},$
(4.3)

where the total energy density and pressure are

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{1}{2}\dot{\phi}^{2}+V+\rho_{\rm EM}+\rho_{\psi},$ (4.4a)
$\displaystyle p$ $\displaystyle=$
$\displaystyle\frac{1}{2}\dot{\phi}^{2}-V+\frac{\rho_{\rm
EM}}{3}+\frac{\rho_{\psi}}{3}.$ (4.4b)

Hence we have

$H^{2}\mathcal{F}=-\frac{a}{2}\frac{dH^{2}}{da}=\frac{1}{M_{\text{pl}}^{2}}\left(\frac{1}{2}\dot{\phi}^{2}+\frac{2}{3}\rho_{\rm
EM}+\frac{2}{3}\rho_{\psi}\right).$ (4.5)

and the system (2.2) becomes

$\displaystyle\frac{d^{2}\phi}{da^{2}}+\frac{4-\mathcal{F}}{a}\frac{d\phi}{da}$
$\displaystyle+$
$\displaystyle\frac{V^{\prime}(\phi)}{a^{2}H^{2}}+\frac{\mathcal{G}}{a^{2}H^{2}f_{\phi}}=0$
(4.6a)
$\displaystyle\frac{d^{2}A_{\lambda}}{da^{2}}+\frac{1}{a}\left(2-\mathcal{F}+\frac{\sigma}{aH}\right)\frac{dA_{\lambda}}{da}$
$\displaystyle+$
$\displaystyle\frac{k}{a^{2}H}\left(\frac{k}{a^{2}H}-\frac{\lambda}{f_{\phi}}\frac{d\phi}{da}\right)A_{\lambda}=0.$
(4.6b)

The Hubble parameter can be computed from the Friedmann equation (4.2a), where
$\rho$ is given by (4.4a). This way, we can compute the value of $H$ and
$\mathcal{F}$ at each time step recursively to feed the equations of motion,
like we already do for $\sigma$ and $\mathcal{G}$. The BD vacuum modes are
identical to the previous case, see Eqs. (3.11). Finally, for comparison
purposes, we can define a generalized time dependent instability parameter
$\xi(a)$ as

$\xi(a)=-\frac{a}{2f_{\phi}}\frac{d\phi}{da}$ (4.7)

such that it corresponds to the definition (3.2). The simulations show that
this parameter, obtained from full solution computation, significantly differs
from the slow roll one at the very end of inflation.

###  4.1 Full numerical results at the end of inflation

In this subsection we will compare our results at the end of inflation, where
we are making a full numerical analysis of the EoM, with those obtained using
the slow roll approximation for the inflationary potential. For the sake of
comparison we will concentrate on the Starobinsky model given by (3.13). In
this current framework, we see in Fig. 3 that the four studied quantities,
namely $\rho_{B}$, $\rho_{E}$, $\mathcal{H}$ and $\mathcal{G}$, are much
closer to the Schwinger equilibrium estimate at the end of inflation.

We present in Fig. 3 the values of the physical observables evaluated at
various stages of inflation, i.e. various values of the scale factor $a$, from
$\epsilon(a)=10^{-3}$ to $\epsilon(a)=1$, as a function of the coupling
$f_{\phi}$ for the Starobinsky model. We superimpose the analytical results
from Secs. 3.1 and 3.2, and hence the backreactionless solution as well as the
Schwinger maximal and equilibrium estimates. From the plots we see that for
$f_{\phi}\lesssim 0.05\,M_{\text{pl}}$ the equilibrium estimate is a good
approximation, especially for $\rho_{E}$ where the predictions of maximal and
equilibrium estimates merge. We also verify that $\cos\theta\simeq 1$ hence
satisfying the assumption on parallel electric and magnetic fields leading to
the conductivity definition (2.10).

In this setup, our numerical code is computing a value of the conductivity
$\sigma$ and $\Delta$ for each time step, hence we got the functions
$\sigma(a)$ and $\Delta(a)$. The variation and presence of $\Delta(a)$ is not
without effect on the final results. Indeed, the smallest ($k\gg H_{E}$) modes
are the ones that most contribute to the integrals (2.8) and (2.9). Without
the Schwinger effect, these modes are produced last, just at the end of
inflation, and only briefly leave the horizon. They therefore should have a
significant impact on preheating. When the Schwinger effect prevents their
generation, by reducing them by a $\ll 1$ factor, while they are still in the
BD vacuum, we can ask ourselves about the effectiveness of gauge preheating.
It was shown in previous studies of gauge preheating Cuissa:2018oiw  that its
efficiency mainly depends on the electromagnetic energy fraction available at
the end of inflation $\rho_{\rm EM}/\rho_{\rm tot}$. To shed light on the last
point, we will extend, in the next section, our numerical results beyond the
end of inflation when the inflaton is coherently oscillationg around its
potential minimum. We will do that in a set of particularly interesting
phenomenological models that we describe in the next section.

###  4.2 Inflationary models

We will here introduce two classes of models that all satisfy the cosmological
constraints. They should be considered as a sample of possible models, and
they are just chosen for illustrative purposes, as they do not exhaust by any
means the allowed inflationary models.

####  4.2.1 $\alpha$-attractor models

The $\alpha$-attractor potential is given by Kallosh:2013yoa

$V_{\alpha}(\phi)=\Lambda_{\alpha}^{4}\left[1-\exp{\left\\{-\sqrt{\dfrac{2}{3\alpha}}\dfrac{|\phi|}{M_{\text{pl}}}\right\\}}\right]^{2}.$
(4.8)

Setting $\alpha=1$ yields the $R^{2}$ model or Starobinsky potential (3.13).
To make the comparison interesting, we choose to have $1\leq\alpha\leq 100$,
where cosmological observables are correctly reproduced. In the slow roll
approximation, the field value at the end of inflation is

$\phi_{E}=\sqrt{\dfrac{3\alpha}{2}}\;M_{\text{pl}}\,\log{\left(1+\dfrac{2}{\sqrt{3\alpha}}\right)}.$
(4.9)

We can readily compute $\phi_{\ast}$, and evaluate the slow roll parameters
$N_{\ast}=60$ $e$-folds before the end of inflation. The slow roll parameters
and the cosmic observables are in agreement with the cosmological contraints
for the range

$1\leqslant\alpha\lesssim 100.$ (4.10)

In particular, for $\alpha=1\,(100)$ we get

$\displaystyle\epsilon_{\ast}$ $\displaystyle\simeq 0.00019\ (0.00387),\hskip
56.9055pt\eta_{\ast}\simeq-0.0159\ (-0.00331)$ (4.11) $\displaystyle n_{s}$
$\displaystyle\simeq 0.967\ (0.97),\hskip 28.45274ptr_{\ast}\simeq 0.003\
(0.062),\hskip 28.45274ptH_{E}\simeq 0.82\ (1.13)\,\cdot 10^{13}\ \text{GeV}.$

in agreement with the observed values Planck:2018jri

$n_{s}^{\rm obs}\simeq 0.9649\pm 0.0042,\quad r_{\ast}^{\rm obs}\lesssim
0.06,\quad H_{\ast}^{\rm obs}\lesssim 6\cdot 10^{13}\,\textrm{GeV}\
(95\%\,\textrm{CL}).$ (4.12)

Using the observed value of $A_{s}$ from Ref. Planck:2018jri , $A_{s}^{\rm
obs}=2.2\cdot 10^{-9}$, we fix the vacuum energy. The result depends on
$\alpha$ and is approximately given by $\Lambda_{\alpha}\simeq 3.4\cdot
10^{-3}\,\alpha^{1/5}\,M_{\text{pl}}$. We then obtain the values
$\Lambda_{1}=3.152\cdot 10^{-3}\,M_{\text{pl}}$ and $\Lambda_{100}=8.313\cdot
10^{-3}\,M_{\text{pl}}$.

Figure 4: Inflaton kinetic and potential energy density, as well as electric,
magnetic and fermion energy density ratios to the initial total energy density
of the Universe for the $\alpha$-attractor models with $\alpha=1$ (upper
panels) and $\alpha=100$ (lower panels). The vertical gray lines display the
value $a$ for which $\epsilon(a)=1$ and the dashed line shows the expected
scaling of the dominant sector.

####  4.2.2 Hilltop quartic models

The hilltop model potential is given by Boubekeur:2005zm

$V_{\rm h}(\phi)=\Lambda_{\rm
h}^{4}\left[1-\left(\dfrac{\phi}{\mu}\right)^{p}\right]^{2}.$ (4.13)

The case $p=4$ can be compatible with the Planck measurements. There are two
ways for the field to relax to the minimum at $\phi=\mu$, with different
initial conditions:

  1. 1.

$\phi_{\ast}>\phi_{E}$: In this case the field $\phi>\mu$ is relaxing in a
potential region that can be approximated by $V_{\rm h}\sim\phi^{8}$, and
thus, the slow roll conditions are not met, as chaotic inflation is ruled out.

  2. 2.

$\phi_{\ast}<\phi_{E}$: In this case the field $\phi<\mu$ is relaxing in a
flat potential region and the model predicts correct inflationary observables
for a large range of the parameter. In this work, we will study this option.

The slow roll parameters and the cosmic observables are in agreement with the
contraints for the range

$10M_{\text{pl}}\lesssim\mu\lesssim 50M_{\text{pl}}.$ (4.14)

We fix the vacuum energy from the constraint on the amplitude of scalar
fluctuations. The result depends on $\mu$ and is approximately $\Lambda_{\rm
h}\simeq 6\cdot 10^{-4}\,\mu^{2/3}\,M_{\text{pl}}^{1/3}$. We then have the
values $\Lambda_{\rm h}=3.243\cdot 10^{-3}\,M_{\text{pl}}$ for
$\mu=10\,M_{\text{pl}}$ and $\Lambda_{\rm h}=8.081\cdot
10^{-3}\,M_{\text{pl}}$ for $\mu=50\,M_{\text{pl}}$.

In particular, for $\mu=10\ (50)\,M_{\text{pl}}$ we get

$\displaystyle\epsilon_{\ast}$ $\displaystyle\simeq 0.00021\ (0.0041)\hskip
68.28644pt\eta_{\ast}\simeq-0.0207\ (-0.00328)$ (4.15) $\displaystyle n_{s}$
$\displaystyle\simeq 0.957\ (0.97)\hskip 36.98866ptr_{\ast}\simeq 0.00335\
(0.0654),\hskip 36.98866ptH_{E}\simeq 0.64\ (1.1)\,\cdot 10^{13}\ \text{GeV}.$

Figure 5: Inflaton kinetic and potential energy density, as well as electric,
magnetic and fermion energy density ratios to the initial total energy density
of the Universe for the hilltop models with $\mu=10\,M_{\text{pl}}$ (upper
panels) and $\mu=50\,M_{\text{pl}}$ (lower panels). The vertical gray lines
display the value $a$ for which $\epsilon(a)=1$ and the dashed line shows the
expected scaling of the dominant sector.

###  4.3 Numerical results beyond the end of inflation

Now that we have established a method to numerically compute the quantities
$\rho_{E}$, $\rho_{B}$, $\rho_{\psi}$, $\mathcal{H}$ and $\mathcal{G}$, we aim
to study the system evolution past $\epsilon=1$, and the onset of reheating.
Indeed, the system (4.6) describes the most general interaction of the zero
mode of both hypercharge gauge and inflaton fields. In particular, no
assumption was made on the Universe geometry, hence there is no specific
reason to stop its numerical computation at the end of inflation. We will also
find it convenient to present some numerical results using as the variable the
number of $e$-folds before the end of inflation $N$, instead of the scale
factor $a$, and related to it by

$N=-\log{\frac{a_{E}}{a}}$ (4.16)

such that $N=0$ corresponds to the time $a_{E}$ when $\epsilon(a_{E})=1$.

Figure 6: Time evolution of the electromagnetic to total energy density
fraction, during and after inflation for various values of the coupling
$f_{\phi}$. The upper panels correspond to the $\alpha$-attractor model with
$\alpha=1$ (top left) and $\alpha=100$ (top right) and the lower panels to the
hilltop model with $\mu=10\,M_{\text{pl}}$ (bottom left) and
$\mu=50\,M_{\text{pl}}$ (bottom right).

We show the postinflationary energy breakdown for selected value of
$f_{\phi}$, for the $\alpha$-attractor models in Fig. 4, $\alpha=1$ (upper
panels) and $\alpha=100$ (lower panels), and the hilltop models of inflation
in Fig. 5, with $\mu=10\,M_{\text{pl}}$ (upper panels) and
$\mu=50\,M_{\text{pl}}$ (lower panels). From the inflaton behavior, we see
that the Universe enters a matter domination era as $\rho_{\phi}\sim a^{-3}$.
For high enough values of $f_{\phi}$, i.e. $f_{\phi}\gtrsim
0.1\,M_{\text{pl}}$, we reproduce the results shown in Ref. Cuissa:2018oiw ,
whereas for $f_{\phi}\lesssim 0.1\,M_{\text{pl}}$ the electric and magnetic
fields exhibit a different behavior: the former decays faster than the latter
while oscillating. This is due to the fact, already mentioned in Ref.
Gorbar:2021rlt , that the energy density for the electric component
$\bm{E}=-\bm{A}^{\prime}$ is much more sensitive to the Schwinger effect than
the magnetic component $\bm{B}$, because it directly couples to the
conductivity in the gauge field equation of motion (2.7). On the other hand,
the magnetic component reflects spatial effects, as it is defined by
$\bm{B}=\bm{\nabla}\wedge\bm{A}$. In this work, we do not consider the
inflaton spatial effects, $\bm{\nabla}\phi$, because this would require one to
implement real fermion interactions in a lattice simulation. Hence, for low
values of $f_{\phi}$, when the Schwinger effect is strongly affecting the
system, the behavior of $\rho_{B}$ is expected to be subject to changes when
the spatial effects are enabled; namely we expect to see a faster decay, like
that of $\rho_{E}$. As also observed in Ref. Gorbar:2021rlt , the electric
field, which is dominant during inflation, becomes subdominant afterwards.
Finally, we can see that for low values of $f_{\phi}$ the fermion energy
density dominates the radiation energy density at the end of inflation as
already highlighted in Ref. Gorbar:2021rlt .

The authors of Ref. Cuissa:2018oiw  quote a sufficient criterion for gauge
preheating to happen, namely that at least an 80% fraction of the total energy
density of the Universe is electromagnetic energy. In the absence of the
Schwinger effect, they found that this criterion is satisfied for values
$f_{\phi}\lesssim 0.1\,M_{\text{pl}}$. However, as expected, the Schwinger
effect significantly reduces the share of electromagnetic energy, as shown on
Fig. 6 for the considered models, which displays the ratio $\rho_{\rm
EM}/\rho_{\rm total}$ for the four previous considered cases. We can see that
the maximum is attained with a value $\sim 10^{-3}$, which precludes any gauge
preheating, at least for $f_{\phi}\gtrsim 0.01\,M_{\text{pl}}$. Another
conclusion from Ref. Cuissa:2018oiw  is that the spatial effects of the
inflaton become relevant for sufficiently low values of $f_{\phi}$ and
contribute to preheating. Since we are neglecting them in our simplified
calculation, any negative statement concerning the possibility of gauge
preheating due to the lack of enough electromagnetic energy should be a
conservative one.

Figure 7: Maximum value of the electromagnetic to total energy fraction as a
function of $f_{\phi}$ for the four considered models: $\alpha$-attractor
models, with $\alpha=1,\,100$, and hilltop models,
$\mu=10,\,50\,M_{\text{pl}}$. Preheating seems unlikely to occur.

The final results from our analysis can be summarized in Fig. 7, where we plot
the maximum value of the electromagnetic to total energy fraction as a
function of $f_{\phi}$ (preheating efficiency) for the Starobinsky model, the
$\alpha$-attractor model with $\alpha=100$ and the hilltop models with
$\mu/M_{\text{pl}}=10,50$. For $f_{\phi}\gtrsim 0.01\,M_{\text{pl}}$, we
obtain

$\frac{\rho_{\rm EM}}{\rho_{\rm tot}}\lesssim 0.01,$ (4.17)

which seems to prevent gauge preheating as its efficiency is far from the
value of $\sim 0.8$ established in the numerical analysis of Ref.
Cuissa:2018oiw .

###  4.4 End of reheating

If gauge preheating does not occur, the inflaton will eventually decay by
perturbative processes which depend on the inflaton total decay width
$\Gamma_{\phi}$. Therefore at the time $t_{\rm rh}\sim 1/\Gamma_{\phi}$, the
inflaton has completely decayed and the radiation domination era starts.

Results from last sections have shown that shortly after inflation ends, the
Universe is dominated by matter, hence we can approximate the Hubble parameter
by

$H\simeq\left(\frac{a_{E}}{a}\right)^{\frac{3}{2}}H_{E},\hskip
28.45274ptH\simeq\frac{2}{3t},$ (4.18)

where $H_{E}\equiv H(a_{E})$, such that

$a_{\rm rh}\simeq
a_{E}\left(\dfrac{3H_{E}}{2\Gamma_{\phi}}\right)^{\frac{2}{3}}$ (4.19)

is the end value after reheating by inflaton perturbative decays. Of course
$a_{\rm rh}$ is a model-dependent quantity, which depends on the value of
$\Gamma_{\phi}$, which in turn, depends on the couplings of the inflaton to
the matter.

In particular, the coupling $1/f_{\phi}$ of the inflaton to the hypercharge
Chern-Simons density provides a channel for the perturbative decay of the
inflaton into a pair of hyperphotons $A$, as $\phi\to AA$. This decay has a
width given by Adshead:2015pva

$\Gamma(\phi\to AA)\simeq\frac{m^{3}_{\phi}}{64\pi f^{2}_{\phi}}.$ (4.20)

where $m_{\phi}$ is the inflaton mass given by

$m^{2}_{\phi}=\left.\dfrac{\partial^{2}V}{\partial\phi^{2}}\right|_{\phi=\phi_{\rm
min}}.$ (4.21)

For the $\alpha$-attractor (hilltop quartic) model, we have $\phi_{\rm
min,\alpha}=0$ ($\phi_{{\rm min},\,h}=\mu$) and

$m^{2}_{\phi,\,\alpha}=\dfrac{4\Lambda_{\alpha}^{4}}{3\alpha
M_{\text{pl}}^{2}},\hskip
56.9055ptm^{2}_{\phi,\,h}=\dfrac{32\Lambda_{h}^{4}}{\mu^{2}}.$ (4.22)

In the simplest case where the inflaton is only coupled to the hypercharge
gauge bosons through the Chern-Simons density, the total width is
$\Gamma_{\phi}=\Gamma(\phi\to AA)$. Using the masses found above we have that

$\Gamma_{\phi}\simeq 12\ (3.0)\cdot
10^{-18}\cdot\dfrac{M_{\text{pl}}^{3}}{f_{\phi}^{2}},$ (4.23)

for $\alpha=1\ (100)$ in the $\alpha$-attractor models, and

$\Gamma_{\phi}\simeq 4.2\ (21)\cdot
10^{-19}\cdot\dfrac{M_{\text{pl}}^{3}}{f_{\phi}^{2}},$ (4.24)

for $\mu=10\ (50)\,M_{\text{pl}}$ in the hilltop models.
The value of the scale factor and the temperature at reheating, $a_{\rm rh}$
and $T_{\rm rh}$, are given by

$\frac{a_{\rm rh}}{a_{E}}\simeq 0.4\left(\frac{T_{\rm rh}}{T_{\rm rh}^{\rm
ins}}\right)^{-4/3},\hskip 56.9055pt\frac{T_{\rm rh}}{T_{\rm rh}^{\rm
ins}}\simeq\sqrt{\frac{\Gamma_{\phi}}{H_{E}}}.$ (4.25)

Consequently we can express $a_{\rm rh}$ and $T_{\rm rh}$ as functions of all
the involved parameters, namely $f_{\phi}$, and $\alpha$ ($\mu$) for
$\alpha$-attractor (hilltop quartic) model. In particular, the relevant
parameter for baryogenesis is the ratio $T_{\rm rh}/T_{\rm rh}^{\rm ins}$
given by

$\frac{T_{\rm rh}}{T_{\rm rh}^{\rm ins}}\simeq 1.9\ (0.8)\cdot
10^{-4}\,\left(\dfrac{0.01}{f_{\phi}/M_{\text{pl}}}\right),$ (4.26)

for $\alpha=1\ (100)$ in the $\alpha$-attractor models, and

$\frac{T_{\rm rh}}{T_{\rm rh}^{\rm ins}}\simeq 0.4\ (0.7)\cdot
10^{-4}\left(\dfrac{0.01}{f_{\phi}/M_{\text{pl}}}\right),$ (4.27)

for $\mu=10\ (50)\,M_{\text{pl}}$ in the hilltop models. As we will see in the
next section the obtained values of the ratio $T_{\rm rh}/T_{\rm rh}^{\rm
ins}$ are fully consistent with the general baryogenesis results, see Fig. 8,
provided that $f_{\phi}\lesssim 0.03\,M_{\text{pl}}$.
In the presence of extra couplings of the inflaton to matter, the predictions
for the inflaton decay width, Eqs. (4.23) and (4.24), and the reheating
temperature, Eqs. (4.26) and (4.27), will change in a model-dependent way, as
well as the model predictions concerning the generation of the baryon
asymmetry.
Of course in the hypothetical case where the explosive production of gauge
fields should have prevailed over the perturbative inflaton decays, gauge
preheating would have taken place over a few $e$-folds after the end of
inflation. As we see from the previous results, this is never the case and
gauge preheating is never strong enough to reheat the Universe after the
period of cosmological inflation. This result does not preclude that, in the
presence of a strong coupling $\lambda$ of the inflaton with some other field,
e.g. a scalar (or a fermion), there could exist an explosive production of
that scalar (or fermion), triggering preheating of the Universe after
inflation Cosme:2022htl .

###  4.5 Baryon asymmetry

Before concluding this paper we wish to make a small comment on the
baryogenesis issue at the electroweak phase transition. In Ref. Cado:2022evn ,
we presented a model of inflation that leads to a successful BAU. The
effective potential for the inflaton, labeled therein as $\chi$, was the
Starobinsky potential999In fact, we used in Ref. Cado:2022evn  an scalar field
$\phi$ non-minimally coupled with gravity as
$\mathcal{L}=-\frac{1}{2}g\phi^{2}R+\dots$, which yields for the canonically
normalized field $\chi$ in the Einstein frame an $\alpha$-attractor potential
with $\alpha=1+\tfrac{1}{6g}\in[4.3,17.6]$, where the lower bound was coming
from imposing the naive unitarity bound $g\phi^{2}<M_{\text{pl}}^{2}$. As the
dependence in $\alpha$ (hence in $g$) is tiny, we choose to show in the
present paper the result for $\alpha=1$, hence for the Starobinsky potential
(which would correspond in Ref. Cado:2022evn  to the limit $g\gg 1$)., and we
did consider the Schwinger equilibrium and maximal estimates. Hence it is
straightforward, using our numerical analysis in this paper, to make an update
of the final results for the BAU for inflation driven by the
$\alpha$-attractor models with $\alpha=1$.

As all details are explained in Secs. 6 and 7 of Ref. Cado:2022evn , we skip
them here and go straight to the final result. First of all we show in Fig. 8
the analogous plot to Fig. 9 of Ref. Cado:2022evn , namely the parameter space
that provides a successful BAU.

Figure 8: The baryogenesis window in the parameter space $(f_{\phi},T_{\rm
rh}/T_{\rm rh}^{\rm ins})$ for the Starobinsky potential ($\alpha$-attractor
model with $\alpha=1$). The dashed line corresponds to Eq. (4.26).

In particular, we display in blue the region where the asymmetry parameter
meets its observational value given by

$\eta_{B}\simeq 4\cdot
10^{-12}\,f_{\theta_{W}}\frac{\mathcal{H}}{H_{E}^{3}}\left(\frac{H_{E}}{10^{13}\,\text{GeV}}\right)^{\frac{3}{2}}\left(\frac{T_{\rm
rh}}{T_{\rm rh}^{\rm ins}}\right)\,\,\simeq\,9\cdot 10^{-11},$ (4.28)

where we have imposed the observed value Zyla:2020zbs  in the right-hand side.
Following Refs. Domcke:2019mnd ; Cado:2021bia  we define the parameter
$f_{\theta_{W}}$, which encodes all the details of the EW phase transition and
its uncertainties, as

$f_{\theta_{W}}=-\sin(2\theta_{W})\,\dfrac{d\theta_{W}}{d\ln
T}\bigg{\rvert}_{T=135\text{ GeV}},\quad 5.6\cdot 10^{-4}\lesssim
f_{\theta_{W}}\lesssim 0.32.$ (4.29)

In addition to their dependence on the gauge sector observables, the
quantities used in this section vary according to the ratio of the reheating
over the instant reheating temperature. This parameter hence adds to
$f_{\phi}$ in the parameter space. The reheating temperature is computed as

$T_{\rm
rh}=\left(\frac{90}{\pi^{2}g_{\ast}}\right)^{\frac{1}{4}}\sqrt{\Gamma_{\phi}M_{\text{pl}}},$
(4.30)

where $g_{\ast}=106.75$ is the SM number of relativistic degrees of freedom,
and we define $T_{\rm rh}^{\rm ins}$ as a reference temperature given by the
above equation with $\Gamma_{\phi}\simeq H_{E}$, which is obtained from the
simulation. It would correspond to the reheating temperature for instant
reheating, and takes the value $T_{\rm rh}^{\rm ins}\simeq 2.87\cdot 10^{15}$
GeV. Using Eq. (4.26) it is possible to link the reheat temperature to the
parameter $f_{\phi}$. The corresponding plot is shown in Fig. 8 which shows
that it provides a wide window for baryogenesis.

Second, we display in orange the region where the magnetic Reynold’s number at
reheating $\mathcal{R}_{m}^{\rm rh}$ is bigger than one, hence ensuring that
the required magnetohydrodynamical conditions are fulfilled for the
(hyper)magnetic fields to survive until the electroweak crossover. As we are
in the viscous regime, it can be computed as Cado:2022evn

$\mathcal{R}_{m}^{\rm rh}\approx 5.9\cdot
10^{-6}\;\frac{\rho_{B}\ell_{B}^{2}}{H_{E}^{2}}\left(\frac{H_{E}}{10^{13}\,\text{GeV}}\right)\left(\frac{T_{\rm
rh}}{T_{\rm rh}^{\rm ins}}\right)^{\frac{2}{3}},$ (4.31)

where $\ell_{B}$ is the physical correlation length of the magnetic field
given by

$\ell_{B}=\frac{2\pi}{\rho_{B}\,a^{3}}\int^{k_{c}}dk\,\frac{k^{3}}{4\pi^{2}}\left(|A_{+}|^{2}+|A_{-}|^{2}\right),$
(4.32)

which can be numerically computed during the simulation in the same way as the
other observables.

Third, and last101010Besides, we checked that the generation of baryon
isocurvature perturbation provides no constraint., we show in green the
condition on the chiral plasma instability (CPI) temperature, ensuring that
the CPI time scale is long enough to allow all right-handed fermionic states
to come into chemical equilibrium with the left-handed ones via Yukawa
coupling interactions (so that sphalerons can erase their corresponding
asymmetries in particle number densities) before CPI can happen. The estimated
temperature at which CPI takes place is

$T_{\rm CPI}/\textrm{GeV}\approx 4\cdot
10^{-7}\,\;\frac{\mathcal{H}^{2}}{H_{E}^{6}}\,\left(\frac{H_{E}}{10^{13}\,\text{GeV}}\right)^{3}\left(\frac{T_{\rm
rh}}{T_{\rm rh}^{\rm ins}}\right)^{2}\,.$ (4.33)

The constraint $T_{\rm CPI}\lesssim 10^{5}$ GeV (the temperature at which
$e_{R}$ comes into chemical equilibrium) guarantees that the CPI cannot occur
before the smallest Yukawa coupling reaches equilibrium and all particle
number density asymmetries are erased, preventing thus the cancellation of the
helicity generated at the reheating temperature.

Therefore, as we can see from Fig. 8, the resulting baryogenesis window for
the Starobinsky potential is close to the Schwinger equilibrium estimate for
$f_{\phi}\lesssim 0.06\,M_{\text{pl}}$, just as the corresponding results on
helicity and magnetic energy density suggest (see the green dots of Fig. 3).
However, for $f_{\phi}\gtrsim 0.06\,M_{\text{pl}}$ there is no space for the
BAU, as the production of gauge fields is too weak, unlike in the previous
results from Ref. Cado:2022evn . In addition to this we have seen that the
reheating temperature is constrained by the model, see Eq. (4.26), as we can
see from Fig. 8 and compatibility of the model reheating temperature with the
baryogenesis results translates into the baryogenesis region on the parameter
$f_{\phi}$

$f_{\phi}\lesssim 0.03\,.M_{\text{pl}}$ (4.34)

Finally, one of the results of this paper is then that baryogenesis at the
electroweak phase transition is favored by low reheating temperatures, in the
range $10^{-6}\,T_{\rm rh}^{\rm ins}\lesssim T_{\rm rh}\lesssim
10^{-3}\,T_{\rm rh}^{\rm ins}$.

##  5 Conclusions

In this paper, we have studied by means of numerical computations the effect
of the Schwinger particle production on the helical hypermagnetic fields
produced at the end of inflation. The inflaton field $\phi$ can decay, through
its coupling to the Chern-Simons density
$\tfrac{\phi}{4f_{\phi}}Y_{\mu\nu}\tilde{Y}^{\mu\nu}$, into helical
hypermagnetic fields in a nonperturbative process. When exiting the vacuum,
the gauge modes are strong enough to create particle/antiparticle pairs of
light fermions, which contribute to the electrical conductivity of the plasma.
The backreaction of fermion currents on the produced gauge fields acts as a
damping force in the explosive production of helical gauge fields. This
effect, called Schwinger effect, was already considered in numerous studies of
inflation and/or baryogenesis, where some analytical and numerical estimates
were computed, mainly in configuration space while our calculation is done in
momentum space.

The equations of motion are in fact a nontrivial integro-differential system.
It was solved numerically by using a fourth order Runge-Kutta method, with
details being displayed in the Appendices. The computed observables of
interest are the electric and magnetic energy density, the helicity as well as
the helicity time derivative. We assumed a homogeneous inflaton with only zero
mode, hence we did not treat any spatial effects. Besides, we also ensured the
convergence of the algorithm and its invariance to the initial conditions.

First of all we have checked that we recover previous results in the slow roll
inflation regime by making the same approximations required by an analytical
resolution. In this way, we validate our code, i.e. we verify that our code
produces the right results in known cases such as the backreactionless case,
where the Schwinger effect is turned off, and the gradient expansion
formalism, where the Bunch-Davies parameter $\Delta$ was first introduced.

In a second step, still in the slow roll regime, we considered a specific
model of inflation, namely the Starobinsky potential, in order to account for
the instability parameter as a function, $\xi(a)$, instead of the constant
imposed by the analytical approximations. That way, we could also implement
the effects of a function $\Delta(a)$ obtained from the plasma evolution on
the gauge production itself.

We then simulated, in a third step, the full system, where neither the slow
roll conditions nor the Universe geometry (e.g. de Sitter) are imposed. In
order words, the inflaton equation of motion was computed alongside with the
gauge one, taking the backreaction of the latter to the former into account
along with the Schwinger effect. We compare our result to the previous setup
and found perfect agreement as long as the slow roll conditions are met. When
inflation is near its end, the full solution diverges from the slow roll
results and produces, as expected, less energy density and helicity.

Finally we will comment on the implication about two related topics: gauge
preheating and baryogenesis. As our code is free from any geometrical issues,
and only requires a model of inflation, we let the simulations run until the
onset of reheating to compute the electromagnetic to total energy density
ratio. We choose two well-known classes of models that satisfy the
cosmological constraints as illustrative examples. Previous studies have
quoted a sufficient criterion for gauge preheating to happen, namely that this
fraction should be at least $\gtrsim$ 80% Cuissa:2018oiw . However, our
numerical estimates suggest that the Schwinger effect significantly reduces
the share of electromagnetic energy for the considered models and preheating
is unlikely to occur. Moreover, since we are neglecting all spatial effects,
any negative statement concerning the possibility of gauge preheating due to
the lack of electromagnetic energy should be a conservative one. On the other
hand our results do apply to the considered class of inflationary models. They
show a certain degree of model dependence, so we cannot exclude a
qualitatively different result for models of inflation other than the
considered ones.

On the other hand, as a successful baryogenesis does depend on a delicate
equilibrium between the amount of helicity, magnetic energy density, and
magnetic correlation length, damped fields do not necessarily mean no baryon
asymmetry in the late Universe. Actually, as a result of our numerical
calculation, we have found there is still a window in the parameter space for
baryogenesis to happen as long as $f_{\phi}\lesssim 0.05\,M_{\text{pl}}$,
while consistency from the perturbative decay channel of the inflaton into
hypergauge bosons implies the bound $f_{\phi}\lesssim 0.03\,M_{\text{pl}}$.
Moreover, baryogenesis is favored for low enough values of the reheating
temperature $T_{\rm rh}\lesssim 10^{-3}\,T_{\rm rh}^{\rm ins}$. Of course, the
baryogenesis predictions should, to some extent, depend on the model of
inflation. In this way, our result here is restricted to the Starobinsky model
and should be considered just as a “proof of existence” for baryogenesis in
the presence of the Schwinger effect.

These two comments should be viewed as hints for future studies that address
the production of gauge fields at the end of inflation. Of course, a full
lattice simulation of the Schwinger effect involving fermions remains to be
done.

## Acknowledgments

YC would like to thank Adrien Florio, for interesting discussions at the early
stage of this work. MQ is indebted to Daniel Figueroa for interesting
discussions, full of new ideas, on the subject. This work is supported by the
Departament d’Empresa i Coneixement, Generalitat de Catalunya Grant No.
2017SGR1069, by the Ministerio de Economía y Competitividad Grant No.
FPA2017-88915-P. IFAE is partially funded by Centres de Recerca de Catalunya.
YC is supported by the European Union’s Horizon 2020 research and innovation
programme under the Marie Sklodowska-Curie Actions No. 754558.

##  Appendix A Numerical method: the slow roll case

We provide here the technical details for the solution of Eq. (3.10) subject
to the initial condition (3.11). For convenience, we implement the numerical
computation in units of $H_{E}$. Writing

$x_{\lambda}(a)=A_{\lambda}(a),\hskip
56.9055pty_{\lambda}(a)=\frac{dA_{\lambda}}{da}(a),$ (A.1)

Eq. (3.10) becomes the following system:

$\frac{d}{da}\left(\begin{array}[]{c}x_{\lambda}\\\
y_{\lambda}\end{array}\right)=\left(\begin{array}[]{cc}0&1\\\
\frac{k}{a^{3}}\left(2\lambda\xi-\frac{k}{a}\right)&-\frac{1}{a}\left(\frac{\sigma}{a}+2\right)\end{array}\right)\left(\begin{array}[]{c}x_{\lambda}\\\
y_{\lambda}\end{array}\right)\quad\Leftrightarrow\quad\frac{d\bm{x}}{da}=\bm{f}(a,\bm{x}).$
(A.2)

To perform each time step $\Delta a$, we use the fourth order Runge-Kutta
(RK4) algorithm:

$\displaystyle\bm{\lambda}_{1}$ $\displaystyle=$
$\displaystyle\bm{f}(a_{i},\bm{x}_{i})$ (A.3a) $\displaystyle\bm{\lambda}_{2}$
$\displaystyle=$ $\displaystyle\bm{f}\left(a_{i}+\tfrac{1}{2}\Delta
a,\bm{x}_{i}+\tfrac{1}{2}\Delta a\bm{\lambda}_{1}\right)$ (A.3b)
$\displaystyle\bm{\lambda}_{3}$ $\displaystyle=$
$\displaystyle\bm{f}\left(a_{i}+\tfrac{1}{2}\Delta
a,\bm{x}_{i}+\tfrac{1}{2}\Delta a\bm{\lambda}_{2}\right)$ (A.3c)
$\displaystyle\bm{\lambda}_{4}$ $\displaystyle=$
$\displaystyle\bm{f}\left(a_{i}+\Delta a,\bm{x}_{i}+\Delta
a\bm{\lambda}_{3}\right)$ (A.3d) $\displaystyle a_{i+1}$ $\displaystyle=$
$\displaystyle a_{i}+\Delta a$ (A.3e) $\displaystyle\bm{x}_{i+1}$
$\displaystyle=$ $\displaystyle\bm{x}_{i}+\tfrac{1}{6}\Delta
a(\bm{\lambda}_{1}+2\bm{\lambda}_{2}+2\bm{\lambda}_{3}+\bm{\lambda}_{4})$
(A.3f)

Note that $\bm{x}$ is complex, hence we solve the above system for both real
and imaginary parts but with their specific initial conditions. These are
mode-dependent as it takes longer for modes with bigger wave number to leave
the BD vacuum. Therefore, we choose as initial condition for each mode

$a_{k,0}=\frac{k}{x_{\rm BD}},$ (A.4)

where we choose the factor $x_{\rm BD}$ in order to make sure that we
initialize the gauge field sufficiently deep inside the Hubble radius. Its
exact value is subject to analysis and is discussed later. As we can see from
(2.8) and (2.9), high values of $k$ are dominating the integral hence large
modes are negligible compared to small ones. This makes us to choose a lower
bound on the $k$ range such that the initial time of the simulation is

$a_{0}=\frac{k_{\rm min}}{x_{\rm BD}}.$ (A.5)

In that way, at $a_{0}$ we make sure that all the modes are in their
respective vacua, which implies $\sigma=0$ as explained above.

In practice, this means that the modes with $k>x_{\rm BD}a$ are given by the
following relations

$\displaystyle{\rm Re}(x_{\lambda,i}^{\rm BD})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\Delta_{i}}{2k}}\cos\frac{k}{a_{i}},$ (A.6a)
$\displaystyle{\rm Im}(x_{\lambda,i}^{\rm BD})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\Delta_{i}}{2k}}\sin\frac{k}{a_{i}},$ (A.6b)
$\displaystyle{\rm Re}(y_{\lambda,i}^{\rm BD})$ $\displaystyle=$
$\displaystyle\frac{1}{a_{i}^{2}}\sqrt{\frac{\Delta_{i}}{2}}\left(\sqrt{k}\sin\frac{k}{a_{i}}-\frac{\sigma_{i}}{2\sqrt{k}}\cos\frac{k}{a_{i}}\right),$
(A.6c) $\displaystyle{\rm Im}(y_{\lambda,i}^{\rm BD})$ $\displaystyle=$
$\displaystyle\frac{1}{a_{i}^{2}}\sqrt{\frac{\Delta_{i}}{2}}\left(-\sqrt{k}\cos\frac{k}{a_{i}}-\frac{\sigma_{i}}{2\sqrt{k}}\sin\frac{k}{a_{i}}\right)$
(A.6d)

while the others are evolving with the RK4 algorithm.

The time steps are distributed on a logarithmic scale

$\log{a_{i}}-\log{a_{i-1}}=\log{a_{i+1}}-\log{a_{i}},$ (A.7)

so that the discretization is the same for each order of magnitude. This means
$\Delta a$ grows exponentially with $a$. The advantage of this method is that
there is a refinement of the grid for small values of $a$, at the beginning of
inflation. The same is done for the discretization in $k$.

We explored the numerical convergence of the solution, both in the number of
$a_{i}$’s, labeled as $N_{a}$, and in the number of $k_{j}$’s, labeled as
$N_{k}$. Provided that $N_{a}>2000$ and $N_{k}>200$, the simulations are very
stable and the output does not depend on the discretization. For big values of
$f_{\phi}$, $f_{\phi}\gtrsim 0.1$, we can even lower the number of time steps
needed.

Besides, we must choose the BD penetration factor $x_{\rm BD}$ such that it
produces trustable results. We have done a numerical analysis and conclude
that depending on the value of $N_{a}$, a range $20<x_{\rm BD}<50$ yields
trustable results. We hence choose throughout this work the following values

$x_{\rm BD}=20,\hskip 28.45274ptN_{a}=500,\,1000,\,2000,\hskip
28.45274ptN_{k}=300.$ (A.8)

At each time step, we compute the electric and magnetic energy density as

$\displaystyle\rho_{E}^{i}$ $\displaystyle=$ $\displaystyle\int_{k_{\rm
min}}^{k_{c}^{i}}dk\,\frac{k^{2}}{4\pi^{2}}\left(|y_{i}^{+}(k)|^{2}+|y_{i}^{-}(k)|^{2}\right),$
(A.9a) $\displaystyle\rho_{B}^{i}$ $\displaystyle=$
$\displaystyle\frac{1}{a_{i}^{4}}\int_{k_{\rm
min}}^{k_{c}^{i}}dk\,\frac{k^{4}}{4\pi^{2}}\left(|x_{i}^{+}(k)|^{2}+|x_{i}^{-}(k)|^{2}\right),$
(A.9b)

where we choose

$k_{c}^{i}=a_{i}\xi+\sqrt{\left(a_{i}\xi\right)^{2}+\frac{a_{i}^{2}}{2}\left[\frac{\sigma_{i}-\sigma_{i-1}}{a_{i}-a_{i-1}}+\frac{\sigma_{i}}{a_{i}}\left(\frac{\sigma_{i}}{2a_{i}}+1\right)\right]}.$
(A.10)

such that we cut off the spectra to retain only modes outside the horizon. The
helicity (2.9a) and its derivative (2.9b) become

$\displaystyle\mathcal{H}_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{a_{i}^{3}}\int_{k_{\rm
min}}^{k_{c}^{i}}dk\,\frac{k^{3}}{2\pi^{2}}\left(|x_{i}^{+}(k)|^{2}-|x_{i}^{-}(k)|^{2}\right),$
(A.11a) $\displaystyle\mathcal{G}_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{a_{i}^{2}}\int_{k_{\rm
min}}^{k_{c}^{i}}dk\,\frac{k^{3}}{2\pi^{2}}\left(|x_{i}^{+}(k)y_{i}^{+}(k)|-|x_{i}^{-}(k)y_{i}^{-}(k)|\right).$
(A.11b)

In the numerics, these integrals are performed numerically over the range of
$k$ that takes $N_{k}$ discrete values. If the Schwinger effect is taken into
account, we turn on the possibility of having $\sigma$ computed at each time
step $a_{i}$ of the numerical computation with

$\sigma_{i+1}=\frac{41\,g^{\prime
3}}{72\,\pi^{2}}\,a_{i}\,\sqrt{2\rho_{B}^{i}}\;\coth\left(\pi\sqrt{\frac{\rho_{B}^{i}}{\rho_{E}^{i}}}\right)$
(A.12)

and injected into the calculation of the next step. Otherwise, we keep it
zero. Last, the fermion energy density is computed as

$\rho_{\psi}^{i}=\frac{\sigma_{i}}{a_{i}^{2}}\int_{k_{\rm
min}}^{k_{c}^{i}}dk\,\frac{k^{2}}{\pi^{2}}\sum_{\lambda=\pm}\left[{\rm
Re}(x_{i}^{\lambda}){\rm Re}(y_{i}^{\lambda})+{\rm Im}(x_{i}^{\lambda}){\rm
Im}(y_{i}^{\lambda})\right].$ (A.13)

Finally, we stop the simulation at $a=a_{E}$. Quantities at that time are
compared to the known analytical results. The color matching dashed vertical
lines in Fig. 2 show the cutoff values $k_{c}^{i}$ computed from (A.10). They
agree perfectly with the point where the BD vacuum modes become dominant for
large $k$.

##  Appendix B Numerical method: full analysis

The numerical implementation follows from the previous case. Defining the
variables

$w=\phi,\qquad x=\frac{d\phi}{da},\qquad y_{\lambda}=A_{\lambda},\qquad
z_{\lambda}=\frac{dA_{\lambda}}{da}$ (B.1)

we transform the above coupled system of differential equations (4.6) into the
system

$\displaystyle\frac{dw}{da}$ $\displaystyle=$ $\displaystyle x$ (B.2a)
$\displaystyle\frac{dx}{da}$ $\displaystyle=$
$\displaystyle-\dfrac{\mathcal{G}}{a^{2}H^{2}f_{\phi}}-\frac{4-\mathcal{F}}{a}\;x-\frac{V^{\prime}(w)}{a^{2}H^{2}}$
(B.2b) $\displaystyle\frac{dy_{\lambda}}{da}$ $\displaystyle=$ $\displaystyle
z_{\lambda}$ (B.2c) $\displaystyle\frac{dz_{\lambda}}{da}$ $\displaystyle=$
$\displaystyle\frac{k}{a^{2}H}\left(\frac{\lambda}{f_{\phi}}\;x-\frac{k}{a^{2}H}\right)\,y_{\lambda}-\frac{1}{a}\left(2-\mathcal{F}+\frac{\sigma}{aH}\right)\;z_{\lambda}$
(B.2d)

which is equivalent to writing

$\frac{d\bm{x}}{da}=\bm{f}(a,\bm{x}).$ (B.3)

We recall that $w,\,x\in\mathbb{R}$ and
$y_{\lambda},\,z_{\lambda}\in\mathbb{C}$. Similarly to the previous
calculation with the slow roll approximation, we use the RK4 algorithm (A.3)
with the values of $H$, $\sigma$, $\mathcal{F}$ and $\mathcal{G}$ computed at
each time step.

Inflaton initial condition could be set to

$w_{0}=\phi_{\ast},\hskip 56.9055ptx_{0}=0.$ (B.4)

However, the number of $e$-folds sets the initial time as $a_{0}={\rm
e}^{-|N_{\ast}|}\sim 10^{-26}$, which is too small a number for the numerical
implementation. We then proceed as follows. For $a\lesssim k_{\rm min}/x_{\rm
BD}$, and sufficiently low $k_{\rm min}$, the gauge field modes stay in their
vacuum and the total contribution to $\Delta(a)$ is negligible. Hence we do
not need to perform the numerical simulation before that time, as the inflaton
is the main player, so we can solve its equation of motion analytically.
Instead, we fix the start of the simulation like before, at $a_{0}=k_{\rm
min}/x_{\rm BD}$ and we compute the corresponding number of $e$-folds $N$
which leads us to the corresponding value of $\phi(N)$. Therefore, the initial
condition must be set to $w_{0}$ such that

$\int_{\phi_{E}}^{w_{0}}\dfrac{V(\phi)}{V^{\prime}(\phi)}\,d\phi=-M_{\text{pl}}^{2}\log{a_{0}}$
(B.5)

and, using $\dot{\phi}\simeq-\frac{V^{\prime}(\phi)}{3H}$ which is valid at
the early stages of inflation,

$x_{0}=-\dfrac{V^{\prime}(w_{0})}{3a_{0}H_{0}^{2}}.$ (B.6)

As for the gauge field, initial conditions are set in the same way as in the
slow roll approximation, see App. A.

## References

  * (1) A. H. Guth, _The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems_ , _Phys. Rev. D_ 23 (1981) 347–356.
  * (2) A. D. Linde, _A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems_ , _Phys. Lett. B_ 108 (1982) 389–393.
  * (3) A. Albrecht and P. J. Steinhardt, _Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking_ , _Phys. Rev. Lett._ 48 (1982) 1220–1223.
  * (4) Planck collaboration, Y. Akrami et al., _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10, [1807.06211].
  * (5) P. Adshead, J. T. Giblin, M. Pieroni and Z. J. Weiner, _Constraining axion inflation with gravitational waves from preheating_ , _Phys. Rev. D_ 101 (2020) 083534, [1909.12842].
  * (6) P. Adshead, J. T. Giblin, M. Pieroni and Z. J. Weiner, _Constraining Axion Inflation with Gravitational Waves across 29 Decades in Frequency_ , _Phys. Rev. Lett._ 124 (2020) 171301, [1909.12843].
  * (7) M. M. Anber and L. Sorbo, _Naturally inflating on steep potentials through electromagnetic dissipation_ , _Phys. Rev. D_ 81 (2010) 043534, [0908.4089].
  * (8) M. M. Anber and E. Sabancilar, _Hypermagnetic Fields and Baryon Asymmetry from Pseudoscalar Inflation_ , _Phys. Rev. D_ 92 (2015) 101501, [1507.00744].
  * (9) Y. Cado and E. Sabancilar, _Asymmetric Dark Matter and Baryogenesis from Pseudoscalar Inflation_ , _JCAP_ 04 (2017) 047, [1611.02293].
  * (10) P. Adshead, J. T. Giblin, T. R. Scully and E. I. Sfakianakis, _Gauge-preheating and the end of axion inflation_ , _JCAP_ 12 (2015) 034, [1502.06506].
  * (11) E. I. Sfakianakis and J. van de Vis, _Preheating after Higgs Inflation: Self-Resonance and Gauge boson production_ , _Phys. Rev. D_ 99 (2019) 083519, [1810.01304].
  * (12) J. R. C. Cuissa and D. G. Figueroa, _Lattice formulation of axion inflation. Application to preheating_ , _JCAP_ 06 (2019) 002, [1812.03132].
  * (13) K. Kamada and A. J. Long, _Baryogenesis from decaying magnetic helicity_ , _Phys. Rev. D_ 94 (2016) 063501, [1606.08891].
  * (14) K. Kamada and A. J. Long, _Evolution of the Baryon Asymmetry through the Electroweak Crossover in the Presence of a Helical Magnetic Field_ , _Phys. Rev. D_ 94 (2016) 123509, [1610.03074].
  * (15) D. Jiménez, K. Kamada, K. Schmitz and X.-J. Xu, _Baryon asymmetry and gravitational waves from pseudoscalar inflation_ , _JCAP_ 12 (2017) 011, [1707.07943].
  * (16) V. Domcke, B. von Harling, E. Morgante and K. Mukaida, _Baryogenesis from axion inflation_ , _JCAP_ 10 (2019) 032, [1905.13318].
  * (17) Y. Cado, B. von Harling, E. Massó and M. Quirós, _Baryogenesis via gauge field production from a relaxing Higgs_ , _JCAP_ 07 (2021) 049, [2102.13650].
  * (18) Y. Cado and M. Quirós, _Baryogenesis from combined Higgs – scalar field inflation_ , _Phys. Rev. D_ 106 (2022) 055018, [2201.06422].
  * (19) V. Domcke and K. Mukaida, _Gauge Field and Fermion Production during Axion Inflation_ , _JCAP_ 11 (2018) 020, [1806.08769].
  * (20) H. Kitamoto and M. Yamada, _Semiclassical analysis of axion-assisted and axion-driven pair production_ , _JHEP_ 06 (2022) 103, [2109.14782].
  * (21) T. D. Cohen and D. A. McGady, _The Schwinger mechanism revisited_ , _Phys. Rev. D_ 78 (2008) 036008, [0807.1117].
  * (22) O. O. Sobol, E. V. Gorbar and S. I. Vilchinskii, _Backreaction of electromagnetic fields and the Schwinger effect in pseudoscalar inflation magnetogenesis_ , _Phys. Rev. D_ 100 (2019) 063523, [1907.10443].
  * (23) E. V. Gorbar, K. Schmitz, O. O. Sobol and S. I. Vilchinskii, _Gauge-field production during axion inflation in the gradient expansion formalism_ , _Phys. Rev. D_ 104 (2021) 123504, [2109.01651].
  * (24) E. V. Gorbar, K. Schmitz, O. O. Sobol and S. I. Vilchinskii, _Hypermagnetogenesis from axion inflation: Model-independent estimates_ , _Phys. Rev. D_ 105 (2022) 043530, [2111.04712].
  * (25) T. Fujita, J. Kume, K. Mukaida and Y. Tada, _Effective treatment of U(1) gauge field and charged particles in axion inflation_ , 2204.01180.
  * (26) J. R. Espinosa, G. F. Giudice, E. Morgante, A. Riotto, L. Senatore, A. Strumia et al., _The cosmological Higgstory of the vacuum instability_ , _JHEP_ 09 (2015) 174, [1505.04825].
  * (27) M. M. Anber and L. Sorbo, _N-flationary magnetic fields_ , _JCAP_ 10 (2006) 018, [astro-ph/0606534].
  * (28) R. Durrer and A. Neronov, _Cosmological Magnetic Fields: Their Generation, Evolution and Observation_ , _Astron. Astrophys. Rev._ 21 (2013) 62, [1303.7121].
  * (29) R. Kallosh, A. Linde and D. Roest, _Superconformal Inflationary $\alpha$-Attractors_, _JHEP_ 11 (2013) 198, [1311.0472].
  * (30) L. Boubekeur and D. H. Lyth, _Hilltop inflation_ , _JCAP_ 07 (2005) 010, [hep-ph/0502047].
  * (31) C. Cosme, D. G. Figueroa and N. Loayza, _Gravitational wave production from preheating with trilinear interactions_ , 2206.14721.
  * (32) Particle Data Group collaboration, P. Zyla et al., _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.