#  When does the Schwinger Preheating Occur?

So Okano  and Tomohiro Fujita

###### Abstract

When the inflaton couples to photons and amplifies electric fields, charged
particles produced via the Schwinger effect can dominate the universe after
inflation, which is dubbed as the Schwinger preheating. Using the hydrodynamic
approach for the Boltzmann equation, we numerically study two cases, the
Starobinsky inflation model with the kinetic coupling and the anisotropic
inflation model. The Schwinger preheating is not observed in the latter model
but occurs for a sufficiently large inflaton-photon coupling in the first
model. We analytically address its condition and derive a general attractor
solution of the electric fields. The occurrence of the Schwinger preheating in
the first model is determined by whether the electric fields enter the
attractor solution during inflation or not.

##  1 Introduction

The Schwinger effect, which is the non-perturbative pair production in the
strong electromagnetic fields, has been theoretically predicted [1, 2, 3]. The
reaction rate in the case of homogeneous and constant electric fields, $E$, is
proportional to $\exp\left(-\pi\frac{m^{2}}{eE}\right)$, where $m$ and $e$
represent the mass and the coupling constant of the produced particles,
respectively. The critical strength of electric fields is given by $E_{\rm
cr}=\pi m^{2}/e\sim 10^{18}\,{\rm V/m}$ for the electrons and the positrons,
and it has not yet been reached by experiments. However, the European Extreme
Light Infrastructure (ELI) project and the X-ray free electron laser (XFEL)
project aim to achieve $E\sim 10^{-1}E_{\rm cr}$ [4, 5], and therefore the
Schwinger effect would be experimentally verified in the near future.

Several physical situations which may cause the Schwinger effect are
discussed, such as relativistic heavy ion collisions [6], neutron stars [7,
8], charged black holes [9, 10], and inflation with the amplification of the
electromagnetic fields [16, 11, 17, 18, 12, 19, 21, 13, 20, 14, 15, 22, 23,
24]. In particular, such inflation models are intensively studied in the
context of magnetogenesis and baryogenesis. Since the production of the
charged particles affects the evolution of the electromagnetic fields, the
Schwinger effect should be taken into account in the analysis. However, it is
a complicated task and various approaches have been discussed [25, 26, 27, 28,
29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. The difficulties of the analysis
arise from the time dependence of the electric fields and the curved
spacetime. In fact, even in the flat spacetime, the Schwinger effect in
general dynamical electric fields is difficult to analyze because we cannot
analytically obtain the mode functions of the charged particles, and thus we
need some prescriptions. One of the well studied prescriptions for the uniform
and dynamical electromagnetic fields is the quantum kinetic approach [39, 40,
41, 42, 43, 44, 45, 46], where the time evolution of the momentum distribution
function of the charged particles is described by the quantum Boltzmann
equation.

Ref. [46] approximately applied the quantum kinetic approach to the kinetic
coupling inflation model without magnetic fields. In this model, the inflaton
is coupled to the electromagnetic fields through the gauge kinetic function,
$f(\phi)^{2}F^{\mu\nu}F_{\mu\nu}/4$ [11]. Ref. [46] found that the charged
particles produced by the Schwinger effect temporarily dominate the universe
right after inflation, which they call “Schwinger preheating”. This is a new
fascinating phenomenon and is worth a closer investigation. In addition, the
oscillation of the electric fields after inflation was found. With a large
amount of charged particles, the electric fields are naively expected to
exponentially decay due to the electric conductivity. In their analysis,
however, the electric fields do not quickly decay but persistently oscillate
with the charged particles. These findings indicate the importance of the
Schwinger effect in the inflation models with the electromagnetic fields.
Nonetheless, they studied only one inflation model and did not give the
condition for the Schwinger preheating. The calculation in Ref. [46] also
assumed the isotropic background spacetime, which is inconsistent with the
existence of the background electric field. When the Schwinger preheating
occurs, the energy density of the electric fields become comparable to the
total energy density at the end of the inflation, and the spacetime anisotropy
may not be negligible. Thus, the generalities of the Schwinger preheating as
well as the electric field oscillation are still unclear.

In this paper, we numerically scrutinize two inflaton models. One is the same
model as Ref. [46] and the other is the anisotropic inflation model [48]. With
analytic arguments, we derive the attractor solution of the electric field and
find the general condition for the Schwinger preheating. Furthermore, we take
into account the anisotropy of the spacetime and estimate its impact on the
dynamics. Our analysis implies that the oscillation of the electric fields can
occur independently from the model of inflation, but the Schwinger preheating
does not always occur.

This paper is organized as follows. In section 2, we introduce the formulation
of the kinetic approach to the inflationary universe and present some
numerical results. In section 3, we derive the attractor solution for the
electric fields and give the condition for the Schwinger preheating. The
influence of the spacetime anisotropy is discussed in section 4. Finally,
section 5 is devoted to the summary and discussion.

##  2 Kinetic approach to the Schwinger effect in the early universe

In this section, we study the dynamics of the charged particles produced via
the Schwinger effect by using Boltzmann equation for one particle distribution
function in the inflationary background and the subsequent reheating era. We
first review the inflation dynamics in an anisotropic background induced by
electric fields. Then we introduce the Boltzmann equation in the anisotropic
background and give the simple representations of the source term and the
collision term by following the procedure of Ref. [46]. We introduce the
moment of the distribution function and reduce the Boltzmann equation into the
hydrodynamic equations. Finally we show the numerical results of two models in
the case of the isotropic background. In Sec. 4, we will discuss the influence
of background anisotropy on the numerical results which we neglect in this
section.

###  2.1 Inflation in the anisotropic background

Here we briefly review the inflation dynamics in the anistropic background and
derive the background equations for the spacetime, inflaton, and electric
fields. Our Lagrangian for the gravity, inflaton, electromagnetic fields, and
charged particles are given by

$\displaystyle\mathcal{L}=-\frac{M_{\rm
Pl}^{2}}{2}\mathcal{R}+\frac{1}{2}\left(\partial_{\mu}\phi\right)^{2}-V(\phi)-\frac{1}{4}f(\phi)^{2}F_{\mu\nu}F^{\mu\nu}+\mathcal{L}_{\psi}.$
(2.1)

Here
$\mathcal{R},~{}\phi,~{}F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$
are the Ricci scalar, the inflaton, the field strength of the $U(1)$ gauge
field $A_{\mu}$, respectively. $V(\phi)$ is the potential of the inflaton, and
$f(\phi)$ is the kinetic coupling function of the inflaton to the $U(1)$ gauge
fields, which will be specified later. $M_{\rm Pl}$ is the reduced Planck
mass. $\mathcal{L}_{\psi}$ is the lagrangian of charged particles interacting
with electromagnetic fields. To be concrete, we consider the Dirac fermions
$\mathcal{L}_{\psi}=\bar{\psi}(i\not{D}-m)\psi$, while the extension to
charged scalars is straightforward. Without loss of generality, we can choose
the $x$-axis as the direction of the electric fields. For simplicity, we
assume that the direction of the electric fields does not change in time. In
the homogeneous universe, we consider $\phi=\phi(t)$ and
$A_{\mu}=\left(0,A(t),0,0\right)$, where we used the gauge fixing condition
$A_{0}=0$ due to the gauge invariance. Note that the magnetic fields are
neglected in this treatment. The consistent background spacetime metric for
this field configuration is known as the Bianchi-I spacetime,

$\displaystyle
ds^{2}=dt^{2}-e^{2\alpha(t)}\left(e^{-4\sigma(t)}dx^{2}+e^{2\sigma(t)}(dy^{2}+dz^{2})\right).$
(2.2)

Here $\alpha(t)$ and $\sigma(t)$ represent the e-foldings of the isotropic
scale factor and deviation from the isotropy respectively. We solve the
following evolution euqations of the spacetime $\alpha(t),\sigma(t)$, the
inflaton $\phi(t)$ and the electric energy density $\rho_{E}(t)$:

$\displaystyle\dot{\alpha}^{2}-\dot{\sigma}^{2}$
$\displaystyle=\frac{1}{3M_{\rm
Pl}^{2}}\left(\rho_{E}+\rho_{\psi}+V(\phi)+\frac{1}{2}\dot{\phi}^{2}\right),$
(2.3) $\displaystyle\ddot{\sigma}+3\dot{\alpha}\dot{\sigma}$
$\displaystyle=\frac{2\rho_{E}}{3M_{\rm Pl}^{2}},$ (2.4)
$\displaystyle\ddot{\phi}+3\dot{\alpha}\dot{\phi}+V^{\prime}(\phi)$
$\displaystyle=2\frac{f^{\prime}}{f}\rho_{E},$ (2.5)
$\displaystyle\dot{\rho}_{E}+4\left(\dot{\alpha}+\dot{\sigma}\right)\rho_{E}$
$\displaystyle=-2\frac{f^{\prime}}{f}\dot{\phi}\rho_{E}-e^{\alpha}\tilde{j}_{\rm
tot}\mathcal{E}+\frac{\dot{\alpha}^{3}}{4\pi^{2}}\left[\dot{\alpha}^{2}+\left(\frac{\dot{f}}{f}\right)^{2}\right],$
(2.6)

where dot denotes the cosmic time derivative and prime denotes the inflaton
derivative, e.g. $f^{\prime}\equiv\partial_{\phi}f(\phi)$. The physical
electric fields are defined as $E(t)=-e^{-\alpha+2\sigma}\dot{A}$ and its
energy density is $\rho_{E}=\frac{1}{2}f^{2}E^{2}$. $\rho_{\psi}$ represents
the energy density of the charged particles. $\tilde{\bm{j}}_{\rm tot}$ is the
total electric current of the charged particles, which flows along the
electric fields $\bm{E}$. We used $e\tilde{\bm{j}}_{\rm
tot}\cdot\bm{E}=\tilde{j}_{\rm tot}\mathcal{E}$ by introducing
$\mathcal{E}\equiv eE(t)$, where $e$ is the coupling constant of the fermion
with the electric fields. The last term in the right hand side of Eq. (2.6) is
phenomenologically added. With the varying kinetic function $f$, the
fluctuation of the gauge field is produced in the inflationary background and
its contribution to the background electric energy is roughly estimated as
[47]

$\displaystyle(\dot{\rho}_{E})_{H}=\left.\frac{d\rho_{E}}{dk}\right|_{k=k_{H}}\cdot\frac{dk_{H}}{dt}\simeq\frac{\dot{\alpha}^{3}}{4\pi^{2}}\left[\dot{\alpha}^{2}+\left(\frac{\dot{f}}{f}\right)^{2}\right].$
(2.7)

This contribution gives a finite initial amplitude for the electric field
which is subsequently amplified by the inflaton. Since the electric field is
eventually stabilized at its attractor value, the initial amplitude is not
very important and we adopt the same expression for $(\dot{\rho}_{E})_{H}$ as
Ref. [47]. The drawback of Eq. (2.6) is that we cannot track the direction of
the electric field. Therefore, when we numerically solve the set of the
equations, we use not Eq. (2.6) but the EoM for electric fields
$\mathcal{E}=ef^{-1}\sqrt{2\rho_{E}}$,111In our numerical calculation, we
multiply the last term by the Heaviside step function
$\theta(\mathcal{E}-10^{-100}\mu^{2})$ to avoid the singularity at
$\mathcal{E}=0$.

$\displaystyle\dot{\mathcal{E}}+2(\dot{\alpha}+\dot{\sigma})\mathcal{E}+2\frac{\dot{f}}{{f}}\mathcal{E}=-\frac{e^{2}}{f^{2}}j_{\rm
tot}+\frac{e^{2}}{f^{2}}\frac{\dot{\alpha}^{3}}{4\pi^{2}\mathcal{E}}\left[\dot{\alpha}^{2}+\left(\frac{\dot{f}}{f}\right)^{2}\right],$
(2.8)

where we introduce $j_{\rm tot}\equiv e^{\alpha}\tilde{j}_{\rm tot}$, which is
the electric current density for the comoving observer.

A well-known example of inflation accompanied by significant electric fields
is the anisotropic inflation model [48]. In this model, the inflaton potential
and the kinetic coupling function are given by

$\displaystyle V_{W}(\phi)=\frac{1}{2}\mu^{2}\phi^{2},\qquad
f_{W}(\phi)=\exp\left(\frac{c\phi^{2}}{2M_{\rm Pl}^{2}}\right),$ (2.9)

where $\mu$ is the mass of inflaton, and $c$ is a dimensionless coupling
constant. When we ignore the charged particles, this system has an attractor
solution of the electric field [48],

$\displaystyle\rho_{E}^{\rm
att}=\frac{3}{2}\frac{c-1}{c^{2}}\epsilon_{V}M_{\rm Pl}^{2}\dot{\alpha}^{2},$
(2.10)

with $\epsilon_{V}=\frac{M_{\rm
Pl}^{2}}{2}\left(\frac{V^{\prime}}{V}\right)^{2}$. Note that we should set
$c>1$ to have a non-vanishing value of the electric field in the attractor
solution.

###  2.2 Boltzmann equation in the anisotropic background

We employ the Boltzmann equation in the curved spacetime with the
electromagnetic fields to describe the dynamics of the charged particles,

$\displaystyle\left[P^{\mu}\partial_{\mu}-\Gamma_{\nu\lambda}^{\mu}P^{\nu}P^{\lambda}\frac{\partial}{\partial
P^{\mu}}-eP^{\mu}F_{\mu j}\frac{\partial}{\partial
P_{j}}\right]\mathcal{F}(x,P)=u_{\mu}P^{\mu}\left(\mathcal{S}[\mathcal{F}]+C[\mathcal{F}]\right),$
(2.11)

where $u,P,\mathcal{S},\mathcal{C}$ and $\mathcal{F}(x,P)$ represent the four-
velocity of the local flow, the four-momentum of the particle, the Schwinger
production term, the collision term, and the distribution function,
respectively. $\Gamma^{\mu}_{\nu\lambda}$ denotes the Levi-Civita
connetiction. When we consider the anisotropic background, Eq. (2.2), and the
electric fields in the $x$-direction, the Boltzmann equation can be rewritten
as

$\displaystyle\frac{\partial\mathcal{F}_{p}}{\partial
t}+\left(\mathcal{E}(t)+3\dot{\sigma}p_{\parallel}\right)\frac{\partial\mathcal{F}_{p}}{\partial
p_{\parallel}}-(\dot{\alpha}+\dot{\sigma})\bm{p}\cdot\frac{\partial\mathcal{F}_{p}}{\partial\bm{p}}=\mathcal{S}[\mathcal{F}_{p}]+\mathcal{C}[F_{p}],$
(2.12)

where we introduced the physical momentum $\bm{p}$ with
$P^{\mu}=(\epsilon_{\bm{p}},\ e^{-\alpha+2\sigma}\bm{p}_{\parallel},\
e^{-\alpha-\sigma}\bm{p}_{\perp})$ and
$\bm{p}^{2}=\bm{p}_{\parallel}^{2}+\bm{p}_{\perp}^{2}.$ One can reproduce the
equation in the isotropic background by taking the limit
$\dot{\sigma}\rightarrow 0$.

To solve the Boltzmann equation, we need to specify the collision term and the
Schwinger source term. In this part, we basically follow the procedure used in
the previous work [46]. We assume that the collision term can be approximated
by a restoring force towards the thermal equilibrium with a relaxation time
$\tau$,

$\displaystyle\mathcal{C}[\mathcal{F}]=-\frac{\mathcal{F}_{\bm{p}}-\mathcal{F}_{\bm{p}}^{\rm
eq}}{\tau}$ (2.13)

where $\mathcal{F}_{p}^{\rm eq}$ represents the Fermi-Dirac distribution
function without chemical potential $\mathcal{F}_{p}^{\rm
eq}=\left(\exp\left[(u^{\nu}P_{\nu}/T)\right]+1\right)^{-1}$ with
$T=[60\rho_{\psi}/(7\pi^{2})]^{1/4}$. It is known that the relaxation time of
the electron in the flat spacetime is $\tau=\frac{c_{0}(4\pi)^{2}}{e^{4}T\ln
e^{-1}}$, where $c_{0}$ is an order 1 parameter [49]. We adjust this
expression to our setup. $e$ should be replaced by the effective coupling
$e_{{\rm eff}}=e/f$. The temperature $T$ in $\tau$ is approximated by the
average kinetic energy of the electron, $\rho_{\psi}/n$. Finally we chose
$c_{0}=1$. Then we obtain the relaxation time as

$\displaystyle\tau=\frac{(4\pi)^{2}}{\left(\frac{\rho_{\psi}(t)}{n(t)}\right)e_{{\rm
eff}}^{4}\ln e_{{\rm eff}}^{-1}}.$ (2.14)

The Schwinger source term in the kinetic approach in the inflationary universe
was estimated as [46],

$\displaystyle\mathcal{S}[f]=\left(1-2\mathcal{F}(t,\bm{p})\right)\sqrt{|\mathcal{E}(t)|}\exp\left(-\pi\frac{m^{2}+\bm{p}^{2}}{|\mathcal{E}(t)|}\right).$
(2.15)

The Pauli blocking factor $(1-\mathcal{F})^{2}$ is expanded into
$\left(1-2\mathcal{F}\right)$ in the leading order of $\mathcal{F}$. One can
show that Eq. (2.15) reproduces the pair production rate for the constant
electric fields,

$\displaystyle\Gamma=\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\mathcal{S}[\mathcal{F}=0]=\frac{|\mathcal{E}|^{2}}{(2\pi)^{3}}\exp\left(-\frac{\pi
m^{2}}{|\mathcal{E}|}\right),$ (2.16)

where we ignored the Pauli blocking factor. We also find the polarization
current density, $j_{\rm pol}$, of the produced charged particles which is not
yet accelerated by the electric fields. In the strong field limit,
$|\mathcal{E}|\gg H^{2},m^{2}$, which is applicable to our case, $j_{\rm pol}$
is computed as

$\displaystyle j_{\rm
pol}=\frac{2g}{\mathcal{E}}\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\sqrt{m^{2}+\bm{p}^{2}}\mathcal{S}[\mathcal{F}=0]\simeq
g\,{\rm sign}(\mathcal{E})\frac{|\mathcal{E}|^{3/2}}{(2\pi)^{3}},$ (2.17)

where $g=2s+1$ is the number of spin projection and $s=1/2$ in our case. Here
we ignored $\mathcal{F}$ again. We distinguish $j_{\rm pol}$ from the
conductive current, $j_{\rm cond}$, which is the current density of the
charged particles accelerated by the electric fields after their production.
When we consider the interaction between the electric fields and charged
particles, we should include the both current densities, $j_{{\rm tot}}=j_{\rm
pol}+j_{\rm cond}$.

Our numerical calculations use Eqs. (2.16) and (2.17) which assume that the
number of the charged particles per unit volume in the phase space is
negligible, $\mathcal{F}\approx 0$. The approximation of the vanishing
$\mathcal{F}$ may be justified, since the rapid expansion of the spacetime and
the acceleration by the electric field change the momenta of the charged
particles and empty the phase space that were filled up by the previous
particle production. Note that Ref. [46] directly solved the Boltzmann
equation without using this approximation $\mathcal{F}\approx 0$ and confirmed
that the hydrodynamic approach reproduced the qualitatively equivalent results
with this approximation. Thus we adopt the same approximation.

Let us introduce a $(k,r)$ moment to apply the hydrodynamic approach for Eq.
(2.12) as,

$\displaystyle\mathcal{J}_{i_{1},\cdots,i_{k}}^{(k,r)}(t)=2g\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}\mathcal{F}(t,\bm{p}),$ (2.18)

where the factor 2 appears from the contribution of the electron and positron.
Multiplying the Boltzmann equation by $\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}$ and integrating over the momentum
space, we obtain the following chain of ordinary differential equations:

$\displaystyle\left[\frac{d}{dt}+\frac{1}{\tau}-(r+k-4)(\dot{\alpha}+\dot{\sigma})-3\dot{\sigma}\right]\mathcal{J}_{i_{1}\cdots
i_{k}}^{(k,r)}-3\dot{\sigma}\sum_{\ell=1}^{k}\mathcal{J}_{i_{1}\cdots\hat{i_{\ell}}\cdots
i_{k}x}^{(k,r)}+3\dot{\sigma}(r+2k-1)\mathcal{J}_{i_{1}\cdots
i_{k}xx}^{(k+2,r-2)}$
$\displaystyle+\mathcal{E}\left[(r+2k-1)\mathcal{J}_{i_{1}\cdots
i_{k},x}^{(k+1,r)}-\sum_{\ell=1}^{k}\delta_{\ell,x}\mathcal{J}_{i_{1}\cdots\hat{i_{\ell}}\cdots
i_{k}}^{(k-1,r+2)}\right]+(r+2k-1)(\dot{\alpha}+\dot{\sigma})m^{2}\mathcal{J}_{i_{1}\cdots
i_{k}}^{(k,r+2)}$
$\displaystyle=2g\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{p_{i_{1}}p_{i_{2}}\dots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}\mathcal{S}(\bm{p})+2g\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{p_{i_{1}}p_{i_{2}}\dots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}\mathcal{F}_{\rm eq}(\bm{p},T).$ (2.19)

The detailed derivation is described in Appendix. A. Here $\hat{i}_{l}$ means
that the index $i_{l}$ is absent in the sequence. Since we consider that the
charged particles are relativistic and $\mathcal{E}\gg m^{2}$, we ignore the
contribution from the mass term in the above equation. We cannot solve the
infinite chains of equations and hence we need closure conditions to cut off
the chains. Now we extract the three pieces of the chains of equations for the
$(0,1),(0,0),(1,0)$ moments, which are important to describe the Schwinger
pair production. We specify each moments as the energy density of the charged
particles, $\mathcal{J}^{(0,0)}\equiv\rho_{\psi}$, their number density,
$\mathcal{J}^{(0,1)}\equiv n$, and their conductive current density along
electric fields, $\mathcal{J}_{x}^{(1,0)}=j_{\rm cond}$. The evolution
equations for them are written as,

$\displaystyle\frac{dn}{dt}+\left(\frac{1}{\tau}+3\dot{\alpha}\right)n=2g\Gamma+\frac{n_{\rm
eq}}{\tau},$ (2.20)
$\displaystyle\frac{d\rho_{\psi}}{dt}+\left(4\dot{\alpha}+\dot{\sigma}\right)\rho_{\psi}-3\dot{\sigma}\mathcal{J}_{xx}^{(2,-2)}=\mathcal{E}j_{\rm
pol}+\mathcal{E}j_{\rm cond},$ (2.21) $\displaystyle\frac{dj_{\rm
cond}}{dt}+\left(\frac{1}{\tau}+3\dot{\alpha}\right)j_{\rm
cond}-3\dot{\sigma}j_{\rm
cond}+3\dot{\sigma}\mathcal{J}_{xxx}^{(3,-2)}=\mathcal{E}\mathcal{J}_{\perp\perp}^{(2,0)},$
(2.22)

where $n_{\rm eq}=2g\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\mathcal{F}_{\rm eq}$.
The reaction rate $\Gamma$ and the polarization current $j_{\rm pol}$ are
given by Eqs. (2.16) and (2.17), respectively. The variables whose expressions
have not yet been specified are $\mathcal{J}_{xx}^{(2,-2)}$,
$\mathcal{J}_{xxx}^{(3,-2)}$ and $\mathcal{J}_{\perp\perp}^{(2,0)}$.

Here, for simplicity, we ignore the anistropy for a while. Taking the
isotropic limit, $\dot{\sigma}\rightarrow 0$, one finds that
$\mathcal{J}_{xx}^{(2,-2)}$ and $\mathcal{J}_{xxx}^{(3,-2)}$ become irrelevant
to the evolution of the system. $\mathcal{J}_{\perp\perp}^{(2,0)}$ is the only
variable remaining relevant but unspecified. Note that we will restore the
anisotropy and discuss its influence in Sec. 4.

Let us identify $\mathcal{J}_{\perp\perp}^{(2,0)}$ through a microscopic
arguments. Assuming the conductive current is approximated by the bulk flow of
the charged particles in the $x$-direction, $j_{\rm cond}=nv_{x}$, we obtain
the EoM for the charged particle by combining Eq. (2.20) and Eq. (2.22) with
$\dot{\sigma}\rightarrow 0$ as,

$\displaystyle\frac{dv_{x}}{dt}=-\frac{v_{x}}{n}\left(2g\Gamma+\frac{n_{\rm
eq}}{\tau}\right)+\frac{\mathcal{E}}{n}\mathcal{J}_{\perp\perp}^{(2,0)}.$
(2.23)

The first term in the R.H.S of the above equation describes the slowing down
of the particle by creating new particles whose averaged velocity is zero and
scattering with the other particles. The second term represents the Lorentz
force from the electric fields. The EoM for the relativistic charged particle
in the isotropic background,
$ds^{2}=dt^{2}-e^{2\alpha}\left(dx^{2}+dy^{2}+dz^{2}\right)$ under the
electric field in the $x$-direction is written as,

$\displaystyle\frac{dv_{x}}{dt}=-Cv_{x}+\mathcal{E}\frac{n}{\rho_{\psi}}(1-v_{x}^{2}),$
(2.24)

where $C$ is the decay constant and we replaced the typical momentum of the
particle $p$ by $\rho_{\psi}/n$ in the second term [23]. Comparing Eq. (2.23)
with Eq. (2.24), we obtain the closure condition as,

$\displaystyle\mathcal{J}_{\perp\perp}^{(2,0)}$
$\displaystyle=\frac{n^{2}-j_{\rm cond}^{2}}{\rho_{\psi}}.$ (2.25)

###  2.3 Time evolution of the system

In this section, we will present the numerical results of the time evolution
of the system. We solve the following EoMs for the inflaton, the spacetime,
and electric fields in the isotropic background:

$\displaystyle\dot{\alpha}^{2}$ $\displaystyle=\frac{1}{3M_{\rm
Pl}^{2}}\left(\rho_{E}+\rho_{\psi}+V(\phi)+\frac{1}{2}\dot{\phi}^{2}\right),$
(2.26) $\displaystyle\ddot{\phi}$
$\displaystyle=-3\dot{\alpha}\dot{\phi}-V^{\prime}(\phi)+2\frac{f^{\prime}}{f}\rho_{E},$
$\displaystyle\dot{\mathcal{E}}+2\dot{\alpha}\mathcal{E}+2\frac{\dot{f}}{{f}}\mathcal{E}$
$\displaystyle=-\frac{e^{2}}{f^{2}}j_{\rm
tot}+\frac{e^{2}}{f^{2}}\frac{\dot{\alpha}^{3}}{4\pi^{2}\mathcal{E}}\left[\dot{\alpha}^{2}+\left(\frac{\dot{f}}{f}\right)^{2}\right].$
(2.27)

We simultaneously solve the evolution of the hydrodynamic variables introduced
in the previous section. The hydrodynamic equations are given by

$\displaystyle\frac{dn}{dt}+\left(\frac{1}{\tau}+3\dot{\alpha}\right)n$
$\displaystyle=2g\Gamma+\frac{n_{\rm eq}}{\tau},$ (2.28)
$\displaystyle\frac{d\rho_{\psi}}{dt}+4\dot{\alpha}\rho_{\psi}$
$\displaystyle=\mathcal{E}j_{\rm pol}+\mathcal{E}j_{\rm cond},$ (2.29)
$\displaystyle\frac{dj_{\rm
cond}}{dt}+\left(\frac{1}{\tau}+3\dot{\alpha}\right)j_{\rm cond}$
$\displaystyle=\mathcal{E}\frac{n^{2}-j_{\rm cond}^{2}}{\rho_{\psi}}.$ (2.30)

To fix its mass, we identify the fermions as electrons and positrons. Their
mass is given by the Yukawa coupling and the Higgs field vev, $h(t)$, which
obtains a large root mean square value due to the quantum fluctuations during
inflation. The mass of the electron is estimated as

$\displaystyle
m(t)=\frac{y_{e}}{\sqrt{2}}h(t)\simeq\frac{y_{e}}{\sqrt{2}}\lambda^{-1/4}\dot{\alpha}(t),$
(2.31)

where $y_{e}=3\times 10^{-6}$ and $\lambda=0.26$. Note that when the particle
production via the Schwinger effect is the most efficient, this mass is
negligibly smaller than the electric field, $m\ll\mathcal{E}$ and thus its
precise value does not affect our final results.

####  2.3.1 Starobinsky model with Ratra coupling

Figure 1:  The numerical result in case of the Starobinsky model $V_{S}$ and
the Ratra coupling $f_{R}$ with the coupling constant $\beta=10$. We show the
time evolution of the energy densities of the electric field $\rho_{E}$ (red
solid), the charged fermions $\rho_{\psi}$ (blue dashed), the inflaton
$\rho_{\phi}$ (orange solid) and their total $\rho_{\rm tot}$ (green dot-
dashed). The calculation starts at $\mu t=0$ corresponding to $N\simeq 60$ and
inflation ends at $\mu t=136$ (vertical black dashed line) where
$\epsilon_{H}\equiv-\ddot{\alpha}/\dot{\alpha}^{2}=1$. The produced fermions
dominate the universe after inflation, $\rho_{\psi}>\rho_{\phi},\rho_{E}$,
which is called the Schwinger preheating.

The first model that we study is the Starobinsky inflation model [50] with the
Ratra coupling function [11],

$\displaystyle V_{\rm S}(\phi)=\frac{3\mu^{2}M_{\rm
Pl}^{2}}{4}\left[1-\exp\left(-\sqrt{\frac{2}{3}}\frac{\phi}{M_{\rm
Pl}}\right)\right]^{2},\qquad f_{R}(\phi)=\exp\left(\beta\frac{\phi}{M_{\rm
Pl}}\right),$ (2.32)

where $\mu\simeq 1.3\times 10^{-5}M_{\rm Pl}$ is the mass of the inflaton
around $\phi=0$. To perform the numerical calculation, we set the initial
condition at $N=60$ based on

$\displaystyle\phi(0)=\sqrt{\frac{3}{2}}M_{\rm
Pl}\ln\left(\frac{4}{3}N\right),\qquad\phi^{\prime}(0)=-\frac{\mu M_{\rm
Pl}}{2N}\sqrt{\frac{3}{2}},$ (2.33)

where $N$ is the e-folding number.

The time evolution of the energy densities obtained by the numerical
calculation is presented in Fig. 1. For most of the inflationary phase, the
electric fields sustained by the quantum fluctuation are not amplified and the
fermion production is negligible. A few e-folds before the inflation end,
however, the electric fields rapidly grow and enter the attractor (see also
Fig. 4), and the Schwinger effect becomes active. After the inflation end,
$\rho_{\psi}$ overwhelms $\rho_{\phi}$ and $\rho_{E}$. Therefore, in this
model, one observes the Schwinger preheating in which the fermions produced by
the Schwinger effect dominate the universe after inflation. The oscillation of
the electric fields is also seen as reported in Ref. [46].

The Schwinger preheating is very fascinating, since it is a novel mechanism of
preheating, it may provide some inflation models with a new channel to reheat
the universe, and it might change the thermal history. Nevertheless, we do not
know the generality of this phenomenon. We investigate if the Schwinger
preheating occurs in another model in the next section.

####  2.3.2 The anisotropic inflation model

Figure 2: The numerical result in the case of the anisotropic inflation model
$V_{W}$ and $f_{W}$ with the coupling constant $c=2$. The plot scheme is the
same as Fig. 1. $\mu t=0$ corresponds to $N\simeq 56$ and inflation ends at
$\mu t=23.5$. The Schwinger preheating does not occur in this case.

As the second model, we investigate the anisotropic inflation model which was
originally proposed in Ref. [48]. Its potential and kinetic coupling have been
given in Eq. (2.9). The advantage of this model is that its dynamics has been
well studied and some analytic solutions are available (e.g. Eq. (2.10)).
Setting the initial condition as

$\displaystyle\phi(0)=M_{\rm Pl}\sqrt{2N},\qquad\phi^{\prime}(0)=0,$ (2.34)

we run the numerical calculation in the same way as the previous model. The
result is shown in Fig. 2. Although the electric fields enter the attractor
much earlier than the previous model and the Schwinger effect lets
$\rho_{\psi}$ exceed $\rho_{E}$ at around the inflation end, the produced
fermions never dominate the universe. Thus, we find a counter-example of the
Schwinger preheating. At least, the anistropic inflation model does not
exhibit it for the coupling constant $c=2$.

Figure 3:  The maximum value of the energy fraction of the charged particles
against the coupling constant. The left and right panels adopt the Starobinsky
model with the Ratra coupling and the anisotropic inflation model,
respectively. The black dashed line denotes the energy fraction equal to 1/2,
and the Schwinger preheating is observed above this line. The orange dotted
line in the right panel represents $(c-1)/2c$ which asymptotes to 1/2. The
Schwinger preheating occurs only in the former model with the sufficiently
strong coupling constant, $\beta\gtrsim 8$.

It is naturally expected that the occurrence of the Schwinger preheating
depends on the coupling constant. We show the dependence in the both models in
Fig. 3. When the maximum value of $\rho_{\psi}/\rho_{\rm tot}$ exceeds $1/2$,
we observe the Schwinger preheating in the numerical calculations and vice
versa. The anisotropic inflation model does not exhibit the Schwinger
preheating, at least for $1<c\leq 10$. Here, although one may think that even
$\rho_{\psi}/\rho_{\rm tot}\approx 0.4$ is practically eligible for the
Schwinger preheating, we adopt $1/2$ as the criteria. On the other hand, the
Schwinger preheating occurs in the first model, if the coupling constant is
large enough, $\beta\gtrsim 8$. Therefore, we find that the occurrence of the
Schwinger preheating depends on models and coupling constants. We will further
explore the dependence on the coupling constant in the next section.

##  3 Analytic Conditions for the Schwinger Preheating

In this section, we analytically derive the conditions for the Schwinger
preheating. Since the Schwinger preheating occurs when the charged particles
dominate the universe and $\rho_{E}$ supplies the energy to $\rho_{\psi}$ via
the Schwinger effect, the necessary condition for the Schwinger preheating is
$\rho_{E}\sim\rho_{{\rm tot}}$. Then we can obtain a model independent
requirement for the Schwinger preheating by using the analytic expression for
$\rho_{E}$, which is a function of $V(\phi)$ and $f(\phi)$.

Fortunately, this system is expected to have an attractor solution of the
electric field for the following reason. The electric fields are amplified by
the inflaton, because the inflaton’s kinetic energy is transfer to the
electric fields through the kinetic function. The amplification might be
weaker than the dilution effect of the spacetime expansion in the early stage
of the inflationary phase. Towards the end of inflation, however, the inflaton
velocity $\dot{\phi}$ is accelerated by the potential force, the amplification
effect is enhanced and eventually overcomes the dilution effect, and the
electric fields begin to grow. When the electric fields are sufficiently
amplified and its backreaction on the inflaton becomes significant, the growth
stops and $\rho_{E}$ is stabilized. This is because, if the electric fields
are further amplified, the inflaton is slowed down by the backreaction, its
kinetic energy pumped to the electric fields is reduced, and they decay. If
the electric fields are smaller than the stabilized value, the inflaton is
accelerated and gives more energy to the electric fields. Therefore, the
inflaton and the electric fields constitute a negative feedback system and
they should have an attractor solution.

In what follows, we first derive the analytic expression for $\rho_{E}$ in the
attractor phase as a function of $V$ and $f$. Second, we show that in the
Starobinsky model with the Ratra coupling, once $\rho_{E}$ reaches the
attractor, $\rho_{E}$ always dominates $\rho_{{\rm tot}}$ at the inflation
end, regardless of the value of $\beta$. Finally we derive the condition for
$\rho_{E}\sim\rho_{\rm tot}$.

###  3.1 The attractor solution for $\rho_{E}$

Here we find the general attractor solution for the energy density of the
electric fields, $\rho_{E}^{\rm att}$. As we discussed above, when the
electric field reaches the attractor, its energy density is stabilized.
Solving Eq. (2.6) without the fermions and the perturbations, we obtain

$e^{4\alpha}f^{2}\rho_{E}={\rm const}.$ (3.1)

Thus, in the attractor phase, the kinetic function is expected to satisfy
$f\propto e^{-2\alpha}$. We rewrite this relation as

$\displaystyle\frac{df}{dt}=-2\dot{\alpha}f\quad\Longrightarrow\quad\dot{\phi}=-2\frac{\dot{\alpha}}{(\ln
f)^{\prime}}.$ (3.2)

Substituting it into Eq. (2.5) and ignoring $\ddot{\phi}$ as a slow-roll
correction term, we obtain

$\displaystyle\rho_{E}^{\rm att}=-\frac{V}{M_{\rm Pl}^{2}(\ln f)^{\prime
2}}+\frac{V^{\prime}}{2(\ln f)^{\prime}},$ (3.3)

where we used the slow-roll equaiton, $3M_{\rm Pl}^{2}\dot{\alpha}^{2}\simeq
V$.

Now we confirm that the above solution, Eq. (3.3), is the attractor solution
by studying the behavior of $\rho_{E}$ around $\rho_{E}^{\rm att}$.
Considering a small deviation from the attractor solution,
$\rho_{E}=\rho_{E}^{\rm att}+\delta\rho_{E}$, we rewrite the EoM for the
inflaton, Eq. (2.5), as

$\displaystyle\delta\rho_{E}=\frac{3\dot{\alpha}}{2(\ln f)^{\prime
2}}\left[2\dot{\alpha}+(\ln f)^{\prime}\dot{\phi}\right].$ (3.4)

The EoM for $\rho_{E}$, Eq. (2.6), reads $\dot{\rho}_{E}+2[2\dot{\alpha}+(\ln
f)^{\prime}\dot{\phi}]\rho_{E}=0$ in the isotropic background without the
fermions. Substituting $\rho_{E}=\rho_{E}^{\rm att}+\delta\rho_{E}$ and
ignoring $\dot{\rho}_{E}^{\rm att}$ as a slow-roll correction as well as the
second order perturbation terms with $\delta\rho_{E}^{2}$, we obtain the
evolution equation for $\delta\rho_{E}$ as

$\displaystyle\delta\dot{\rho}_{E}$ $\displaystyle=-\frac{4}{3}\frac{(\ln
f)^{\prime 2}\rho_{E}^{\rm att}}{\dot{\alpha}}\delta\rho_{E}.$ (3.5)

This equation implies that $\delta\rho_{E}$ exponentially decays with a time
scale, $\tau^{\rm att}=\frac{3\dot{\alpha}}{4(\ln f)^{\prime 2}\rho_{E}^{\rm
att}}$. Therefore, $\rho_{E}^{\rm att}$ is an attractor solution, when
$\tau^{\rm att}$ is sufficiently short. The comparison between $\tau^{\rm
att}$ and the Hubble time $\dot{\alpha}^{-1}$ is given by

$\displaystyle\dot{\alpha}\tau^{\rm att}=\frac{\rho_{{\rm tot}}/\rho_{E}^{\rm
att}}{4M_{\rm Pl}^{2}(\ln f)^{\prime 2}}=2.5\times
10^{-2}\times\begin{cases}\left(\frac{\beta}{10}\right)^{-2}\left(\frac{\rho_{E}^{\rm
att}/\rho_{{\rm
tot}}}{10^{-1}}\right)^{-1}\,&\big{(}f=e^{\beta\frac{\phi}{M_{\rm
Pl}}}\big{)}\\\ c^{-2}\left(\frac{\phi}{10M_{\rm
Pl}}\right)^{-2}\left(\frac{\rho_{E}^{\rm att}/\rho_{{\rm
tot}}}{10^{-1}}\right)^{-1}\,&\big{(}f=e^{\frac{c}{2}\left(\frac{\phi}{M_{\rm
Pl}}\right)^{2}}\big{)}.\end{cases}$ (3.6)

Since $\dot{\alpha}\tau^{\rm att}$ is proportional to $\rho_{\rm
tot}/\rho_{E}^{\rm att}$, this attractor is very weak for negligibly small
electric fields, $\rho_{E}\ll\rho_{\rm tot}$. However, the attractor becomes
stronger, as the electric fields are amplified and significantly backreact on
the inflaton, as expected.

In the case of the anistropic inflation model (i.e. $V=\mu^{2}\phi^{2}/2$ in
the second line of Eq. (3.6)), we can further reduce the time scale as
$\dot{\alpha}\tau^{\rm att}=\frac{1}{4(c-1)}$, which is consistent with the
known solution of $\rho_{E}$. Thus, in the anisotropic inflation model,
$\rho_{E}^{\rm att}$ is always the attractor solution irrespective of $\phi$,
unless $c-1$ is negative or tiny.

Figure 4: The energy density of the electric fields $\rho_{E}$ (blue) and that
in the attractor phase $\rho_{E}^{\rm att}$ given in Eq. (3.3) (orange), which
are normalized by the total energy density $\rho_{\rm tot}$. These numerical
calculations do not include the fermions. The left and right panels adopt the
Starobinsky model with the Ratra coupling for $\beta=10$ and the anisotropic
inflation model for $c=2$, respectively. The horizontal axis denotes the
e-folding number which is zero at the end of inflation ($\epsilon_{H}=1$) and
takes negative values during inflation. In the left panel, $\rho_{E}$ does not
reach the attractor until the last few e-folds of inflation, while it is well
approximated by $\rho_{E}^{\rm att}$ even after inflation ends. In the right
panel, on the other hand, $\rho_{E}$ stays at the attractor solution during
inflation and deviates from it at around the end of inflation. When we include
the fermions in the numerical calculation, the right panel does not change
much, but in the left panel $\rho_{E}$ abruptly drops after the inflation end.

We find that our attractor solution, Eq. (3.3), matches the numerical results
very well. Fig. 4 presents the comparisons between the numerically computed
$\rho_{E}$ and our analytic expression $\rho_{E}^{\rm att}$. Note that these
numerical calculations do not include the fermion for illustrative purposes.
The right panel shows that the electric fields stay at the attractor for most
of the inflationary phase in the second model. However, the left panel shows
that the electric fields enter the attractor only a few e-folds before the
inflation end in the first model for $\beta=10$. Since $\dot{\alpha}\tau^{\rm
att}\propto\beta^{-2}$ in the Ratra coupling case, for a smaller coupling
constant, the attractor becomes weaker and the electric fields may not enter
the attractor phase during inflation.

In Fig. 5, we show the evolution of the energy densities in the first model
for the coupling constant $\beta=7.7$, which is slightly smaller than the
threshold value for the Schwinger preheating $\beta\approx 8$. In this case,
the electric fields never enter the attractor and the Schwinger preheating
does not occur. It infers that whether the electric fields reach the attractor
solution is crucial to the Schwinger preheating. We will study this issue in
more detail in Sec. 3.3. Before that, let us first investigate the cases where
the electric fields enter the attractor.

Figure 5:  The same plot as Fig. 1 for a smaller coupling constant,
$\beta=7.7$. The calculation starts at $\mu t=0$ corresponding to $N\simeq 60$
and inflation ends at $\mu t=120$ (vertical black dashed line). The attractor
solution $\rho_{E}^{\rm att}$ given in Eq.(3.3) is shown as a black solid line
until inflation ends. The electric fields do not enter the attractor phase
during inflation. The maximum value of $(\rho_{E}+\rho_{\psi})/\rho_{\phi}$ is
roughly equal to 1, and the Schwinger preheating is marginally absent.

###  3.2 The inflation end with the attractor

We have found that the Schwinger preheating occurs in the Starobinsky model
with the Ratra coupling, if the coupling constant is large enough,
$\beta\gtrsim 8$. In the same model, the electric fields enter the attractor
during inflation for the same range of the coupling constant, $\beta\gtrsim
8$. Now we investigate the inflation end in the case that the electric fields
are in the attractor phase. One finds that the slow-roll parameter is
rewritten as

$\displaystyle\epsilon_{H}\equiv-\frac{\ddot{\alpha}}{\dot{\alpha}^{2}}=\frac{\dot{\phi}^{2}/2+2\rho_{E}/3}{\dot{\alpha}^{2}M_{\rm
Pl}^{2}}=\frac{2}{(\ln f)^{\prime 2}M_{\rm Pl}^{2}}+2\Omega_{E},\qquad({\rm
attractor\ phase})$ (3.7)

where $\dot{\phi}$ was evaluated with Eq. (3.2) which holds in the attractor
phase, and $\Omega_{E}\equiv\rho_{E}/\rho_{\rm tot}$. The first and the second
term in the R.H.S. represent the contributions from the inflaton kinetic
energy and the electric energy density, respectively. As they grow, the slow-
roll parameter $\epsilon_{H}$ increases and inflation eventually ends at
$\epsilon_{H}=1$.

In the case of the Ratra coupling $f_{R}$, we have $(\ln
f)^{\prime}=\beta/M_{\rm Pl}$ and obtain

$\displaystyle\epsilon_{H}=\frac{2}{\beta^{2}}+2\Omega_{E}.\qquad({\rm
attractor\ phase\ \&\ Ratra\ coupling})$ (3.8)

Since the first term is constant, $\Omega_{E}$ must be larger than one-half to
end inflation. $\Omega_{E}$ does not have to stop increasing at $1/2$, and may
become even larger after inflation, until the electric fields are
significantly consumed to produce and accelerate the fermions. This is a
general property of inflationary models with the Ratra coupling in the
attractor phase of the electric fields.

The above argument apparently suggests that the models with Ratra coupling
always have a large $\Omega_{E}$ at the inflation end and lead to the
Schwinger preheating. However, as shown in Fig. 3, the Starobinsky model with
the Ratra coupling does not cause the Schwinger preheating for $\beta\lesssim
8$. This is because, for a smaller coupling constant, the electric fields do
not reach the attractor solution, as seen in Fig. 5.

###  3.3 Condition to enter the attractor

Entering the attractor phase of the electric fields is a requirement to cause
the Schwinger preheating in the Starobinsky model with the Ratra coupling. In
this section, we will derive the condition to reach the attractor solution.

Let us introduce the inflaton value $\phi^{\ast}$ at which the electric fields
begin to grow by overcoming the cosmic dilution. At this moment, the kinetic
function transiently satisfies $\dot{\phi}^{\ast}=-2\dot{\alpha}/(\ln
f)^{\prime}=-2\dot{\alpha}M_{\rm Pl}/\beta$ for the same reason as Eq. (3.2).
Using the slow-roll approximation, we obtain $\phi^{\ast}$ as

$\displaystyle\phi^{\ast}$
$\displaystyle=\sqrt{\frac{3}{2}}\log\left(1+\frac{\sqrt{6}}{3}\beta\right)M_{\rm
Pl}.$ (3.9)

At $\phi=\phi^{\ast}$, the energy density of the electric field sourced by the
perturbation is also obtained as

$\dot{\rho}_{E}+8\dot{\alpha}\rho_{E}=\frac{5\dot{\alpha}^{5}}{4\pi^{2}},\quad\Longrightarrow\quad\rho_{E}(\phi^{\ast})=\frac{5\dot{\alpha}^{4}_{\ast}}{32\pi^{2}}=\frac{5V^{2}(\phi^{\ast})}{288\pi^{2}M_{\rm
Pl}^{4}}.$ (3.10)

Here we ignored the anisotropy and the fermions, and used
$\dot{f}/f=2\dot{\alpha}$ in Eq. (2.6). After $\phi=\phi^{\ast}$, from Eq.
(3.1), the electric energy density is given by

$\displaystyle\rho_{E}(\phi)=\frac{5V^{2}(\phi^{\ast})}{288\pi^{2}M_{\rm
Pl}^{4}}\exp\left[2\beta(\phi^{\ast}-\phi)+4N\right],$ (3.11)

where we introduced the e-folding number which is computed as,

$\displaystyle
N=\int_{t^{\ast}}^{t}Hdt=\int^{\phi^{\ast}}_{\phi}\frac{V}{V^{\prime}}d\phi=\left[\frac{3}{4}e^{\sqrt{\frac{2}{3}}\phi}-\sqrt{6}\phi\right]^{\phi^{\ast}}_{\phi}.$
(3.12)

Now we can calculate when $\rho_{E}(\phi)$ reaches $\rho_{E}^{\rm att}$.
Numerically solving

$\displaystyle\rho_{E}(\phi^{\rm att})=\rho^{\rm att}_{E}(\phi^{\rm att}),$
(3.13)

we obtain the inflaton value $\phi^{\rm att}(\beta)$, at which the electric
fields enter the attractor phase, as a function of $\beta$. Finally, we impose
a condition that the electric fields enter the attractor before the end of
inflation at $\phi=\phi_{\rm end}$,

$\displaystyle\phi^{\rm att}(\beta)>\phi_{\rm end}.$ (3.14)

Figure 6:  The inflaton value $\phi^{\rm att}(\beta)$ at which the electric
fields enter the attractor phase and which solves Eq. (3.13). For the electric
fields to enter the attractor during inflation, $\phi^{\rm att}(\beta)$ must
be larger than $\phi_{\rm end}$, which sets the lower bound on $\beta$. The
dashed lines represent the values of $\phi$ satisfying
$\epsilon_{V}(\phi)=1,0.3$ and $0.1$ from bottom to top, and they cross
$\phi^{\rm att}(\beta)$ at $\beta=6.6,7.6$ and $9.4$

In Fig. 6, we plot $\phi^{\rm att}(\beta)$ and compare it with the inflaton
values satisfying $\epsilon_{V}(\phi)=0.1,0.3$ and 1. If one adopt
$\epsilon_{V}(\phi_{\rm end})=1$ to characterize the inflation end, one
obtains $\beta=6.5$ as the threshold value, which is a slightly milder
condition than the numerical result shown in Fig. 3. Note that since we
extrapolated the dynamics around the end of inflation by using the expressions
derived under the slow-roll approximation, we may well have a residual error.
In fact, when one uses $\epsilon_{V}(\phi_{\rm end})=0.3$ instead, the
threasold would be $\beta=7.6$ which is closer to the numerical result.

Consequently, we found that the occurrence of the Schwinger preheating in the
Starobinsky model with the Ratra coupling is determined by whether the
electric fields enter the attractor solution during inflation or not. For a
weaker coupling than $\beta=8$, the amplification of the electric fields is
too slow to reach the attractor by the end of inflation, and the charged
fermions cannot receive a sufficiently large energy to dominate the universe.

##  4 Influence of the Anisotropy

We have neglected the spacetime anisotropy since Eq. (2.22). However, the
isotropic spacetime is inconsistent with the background electric field, when
it is strong enough to significantly contribute to the spacetime evolution. In
this case, we should work in the anisotropic background, Eq. (2.2). In this
section, we restore the anisotropy and investigate its influence on the
evolution of the matter fields.

In Eqs. (2.20)-(2.22), we wrote the hydrodynamic equations in the anisotropic
background, while we left $\mathcal{J}_{xx}^{(2,-2)}$ and
$\mathcal{J}_{xxx}^{(3,-2)}$ unspecified, which became irrelevant in the
isotropic limit $\dot{\sigma}\to 0$. Here let us determine their expressions
through phenomenological arguments in similar ways to
$\mathcal{J}_{xx}^{(2,0)}$. Assuming $j_{\rm cond}=nv_{x}$ and combining Eq.
(2.20) and (2.22), we obtain the EoM for the charged particles as

$\displaystyle\frac{dv_{x}}{dt}=-\frac{v_{x}}{n}\left(2g\Gamma+\frac{n_{\rm
eq}}{\tau}\right)+\frac{\mathcal{E}}{n}\mathcal{J}_{\perp\perp}^{(2,0)}+3\dot{\sigma}\left(v_{x}-\frac{\mathcal{J}_{xxx}^{(3,-2)}}{n}\right).$
(4.1)

Compared to the istropic case, Eq. (2.23), the anistropy induces the last term
in the R.H.S.222Since we took the massless limit, the contribution from the
isotropic expansion vanishes. The EoM for the relativistic charged particles
in the anisotropic background (2.2) under the electric fields in the
$x$-direction is written as

$\displaystyle\frac{dv_{x}}{dt}=-Cv_{x}+\mathcal{E}\frac{n}{\rho_{\psi}}(1-v_{x}^{2})+3\dot{\sigma}(v_{x}-v_{x}^{3}).$
(4.2)

Comparing Eq. (4.1) with Eq. (4.2), we obtain an additional closure condition
as

$\displaystyle\mathcal{J}_{xxx}^{(3,-2)}$ $\displaystyle=j_{\rm
cond}\left(\frac{j_{\rm cond}}{n}\right)^{2}.$ (4.3)

In the meanwhile, the same expressions for $\mathcal{J}_{xxx}^{(3,-2)}$ as
well as $\mathcal{J}_{\perp\perp}^{(2,0)}$ can be obtained by approximately
manipulating their definitions,

$\displaystyle\mathcal{J}_{\perp\perp}^{(2,0)}$ $\displaystyle\equiv
n\left<\frac{p_{\perp}^{2}}{\epsilon_{\bm{p}}^{3}}\right>\simeq
n\frac{1-\braket{\frac{p_{\parallel}}{\epsilon_{\bm{p}}}}^{2}}{\left<\epsilon_{\bm{p}}\right>}=\frac{n^{2}-j_{\rm
cond}^{2}}{\rho_{\psi}},$ $\displaystyle\mathcal{J}_{xxx}^{(3,-2)}$
$\displaystyle\equiv
n\left<\left(\frac{p_{\parallel}}{\epsilon(\bm{p})}\right)^{3}\right>\simeq
n\left<\frac{p_{\parallel}}{\epsilon(\bm{p})}\right>^{3}=j_{\rm
cond}\left(\frac{j_{\rm cond}}{n}\right)^{2},$ (4.4)

where we used $v_{x}=\braket{\frac{p_{\parallel}}{\epsilon_{\bm{p}}}}$,
$\rho_{c}=n\braket{\epsilon_{\bm{p}}}$, and the bracket $\braket{\cdots}$
denotes the average over the phase space weighted by $\mathcal{F}(t,\bm{p})$.
Although the phenomenological argument to identify $\mathcal{J}_{xx}^{(2,-2)}$
is unknown, we can apply a similar approximation to find its expression as

$\displaystyle\mathcal{J}_{xx}^{(2,-2)}\equiv
n\left<\frac{p_{\parallel}^{2}}{\epsilon}\right>\simeq
n\left<\epsilon\right>\left<\frac{p_{\parallel}}{\epsilon}\right>^{2}=\rho_{\psi}\left(\frac{j_{\rm
cond}}{n}\right)^{2}.$ (4.5)

Without the verification from the microscopic EoM, one may be concerned with
the accuracy of the closure condition for $\mathcal{J}_{xx}^{(2,-2)}$.
Nevertheless, since the term with $\mathcal{J}_{xx}^{(2,-2)}$ evaluated by the
above expression is always subdominant in Eq. (2.21) in our numerical
calculations, a small error in the expression, if any, would not affect our
final results.

Figure 7: The time evolution of the energy fraction
$(\rho_{E}+\rho_{\psi})/\rho_{\phi}$ computed in the anisotropic (blue) and
the isotropic background (orange). The left and right panels adopt the
Starobinsky model with the Ratra coupling for $\beta=10$ and the anisotropic
inflation model for $c=2$, respectively. The background anisotropy makes the
peak height and its timing of the energy fraction slightly higher and earlier
in the left panel, respectively, while it is opposite (i.e. slightly lower and
later) in the right panel.

We obtain the full hydrodynamic equations for this system, and we can perform
numerical calculations in the anisotropic background. We consider the same
models as Sec. 3, which are the Starobinsky model with the Ratra coupling, Eq,
(2.32), and the anisotropic inflation model, Eq (2.9). The numerical results
in the anisotropic and isotropic background are compared in Fig. 7. The
anisotropy does not have a dramatic impact on the evolution of the energy
densities. One observes that $(\rho_{E}+\rho_{\psi})/\rho_{\phi}$ is slightly
shifted by the anisotropy, while the overall behavior stays almost the same in
the both models.

##  5 Summary and Discussion

In this paper, we applied the kinetic hydrodynamic approach to study the
Schwinger effect in the Starobinsky inflation model with the Ratra coupling
and the anisotropic inflation model. The former model was analyzed in the
previous work [46], and they concluded that the Schwinger preheating occurred
when the coefficient of the Ratra coupling satisfies $\beta\geq 8$ without any
analytical discussion. While we confirmed the results of the previous work, we
have shown that the Schwinger preheating does not occur in the latter model.

We also analytically discussed the condition for the Schwinger preheating. The
Schwinger preheating occurs when the energy density of the charged particles
dominates the universe and the energy density of the charged particle is
supplied from electric fields. Thus the necessary conditions for the Schwinger
preheating is $\rho_{E}/\rho_{{\rm tot}}\sim 1/2$, which leads to the
significant backreaction from the electric fields on the inflaton. We found
the general attractor solution of $\rho_{E}$ by using the slow-roll
approximation, and investigated the stability of the attractor solutions. Our
result is consistent with the conventional analysis of the anisotropic
inflation model. Using the attractor solution, we evaluated the electric field
energy density at the inflaton end. In the Starobinsky model with the Ratra
coupling, the analytic evaluation apparently indicates that $\rho_{E}$ always
dominates the universe at the end of inflation, independent of $\beta$.
However, we found that $\rho_{E}$ does not reach the attractor, and hence it
remains much smaller than $\rho_{\phi}$ in the case of a small $\beta$. The
threshold to enter the attractor is estimated as $\beta=7.5$ when the
inflation end is evaluated as $\epsilon_{H}=0.3$, while the numerical results
indicate the threshold is $\beta\approx 8$.

We have also addressed the impact of the spacetime anisotropy. Two more
closure conditions to cut the chain equation are needed, and we specified
them. However we still do not know the physical meaning of one of the two
closure conditions. Although the energy fraction of the electric fields
approaches 1/2 and a significant effect of the anisotropy on the background
evolution could be expected, the numerical results showed that the dynamics is
almost same as the isotropic case except for the slight shifts of the timing
of the Schwinger preheating and the maximum value of
$\left(\rho_{E}+\rho_{\psi}\right)/\rho_{\phi}$. We found that the anisotropy
hastens the preheating in Starobinsly and Ratra model but delays it in
anisotropic inflation. The origin of these difference is beyond the scope of
this paper.

The Schwinger preheating may give a new reheating channel to inflation models
and change the thermal history of the universe. Hence it is meaningful to
investigate its dynamics and find its condition. In our analysis, we do not
consider a comprehensive model including inflaton decay processes. If the
inflaton itself or its decay products are coupled to the charged fermions
produced by the Schwinger effect, for instance, the Schwinger preheating might
trigger the reheating. It would be interesting to introduce these couplings
and study the dynamics after the Schwinger preheating. We hope to come back to
these issues in the future work.

The quantum kinetic approach is well studied approach to the Schwinger effect
in the case of the spatially uniform dynamical electric fields. According to
Ref. [52], the general form of the Schwinger source term includes the
integration of the distribution function in time except for the case of low
number density of fermions, which is the same situation as our case. Since the
deviation from the thermal equilibrium is not always small in the low number
density, the low number density approximation is possibly in conflict with the
ansatz where the collision term is described by the restoring force. In order
to the more general analysis, we should relax either one of the two
assumptions. We will investigate the other prescription to deal the kinetic
approach in the future work.

## Acknowledgments

We would like to thank Keigo Shimada, Shin Kobayashi, Shintaro Sato and
Teruaki Suyama for useful discussions and comments. The work of T.F. was
supported by JSPS KAKENHI No. 17J09103 and No. 18K13537.

##  Appendix A Derivation of the hydrodynamic Boltzmann equation

Let us obtain the Eq. (2.2) from the Boltzmann equation in the Bianchi-I
spacetime, Eq. (2.12).

### The first term in Eq. (2.12)

Since the moment of the distribution function is independent from the
momentum, the partial time derivative for $\mathcal{F}_{\bm{p}}$ becomes a
total derivative as,

$\displaystyle 2g\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}\frac{\partial\mathcal{F}_{\bm{p}}}{\partial
t}=\frac{d}{dt}\mathcal{J}^{(k,r)}_{i_{1},\cdots,i_{k}}.$ (A.1)

### The second term in Eq. (2.12)

By using the integration by part, the term proportional to $\mathcal{E}(t)$ is
written as

$\displaystyle 2g\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}\frac{\partial\mathcal{F}_{\bm{p}}}{\partial
p_{\parallel}}=(r+2k-1)\mathcal{J}^{(k+1,r)}_{i_{1},\cdots
i_{k},x}-\sum_{\ell=1}^{k}\delta_{i_{\ell},x}\mathcal{J}^{(k+1,r+2)}_{i_{1},\cdots,\hat{i}_{i_{\ell}},\cdots,x},$
(A.2)

where we used the relation,

$\displaystyle\epsilon(\bm{p})$ $\displaystyle=\sqrt{\bm{p}^{2}+m^{2}},$ (A.3)
$\displaystyle\frac{\partial}{\partial
p_{\parallel}}\left(\epsilon(\bm{p})^{-r-2k+1}\right)$
$\displaystyle=-(r+2k-1)\epsilon(\bm{p})^{-r-2k-1}p_{\parallel}.$ (A.4)

In the same way as $\mathcal{E}(t)$ term, the term proportional to
$3\dot{\sigma}$ is written as,

$\displaystyle 2g\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}p_{\parallel}\frac{\partial\mathcal{F}_{\bm{p}}}{\partial
p_{\parallel}}$
$\displaystyle=-2g\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{\partial}{\partial
p_{\parallel}}\left(\frac{p_{i_{1}}\cdots
p_{i_{k}}p_{\parallel}}{\epsilon(\bm{p})^{r+2k-1}}\right)\mathcal{F}_{p},$
$\displaystyle=(r+2k-1)\mathcal{J}_{i_{1}\cdots
i_{k}xx}^{(k+2,r-2)}-\mathcal{J}_{i_{1}\cdots
i_{k}}^{(k,r)}-\sum_{\ell=1}^{k}\mathcal{J}_{i_{1}\cdots\hat{i_{\ell}}\cdots
i_{k}x}^{(k,r)},$ (A.5)

where we used the following relation,

$\displaystyle\frac{\partial}{\partial
p_{\parallel}}\left(\frac{p_{i_{1}}\cdots
p_{i_{k}}p_{\parallel}}{\epsilon(\bm{p})^{r+2k-1}}\right)=-(r+2k-1)\frac{p_{i_{1}}\cdots
p_{i_{k}}p_{\parallel}^{2}}{\epsilon(\bm{p})^{r+2k+1}}+\sum_{\ell=1}^{k}\delta_{\ell,x}\frac{p_{i_{1}}\cdots\hat{p}_{i_{\ell}}\cdots
p_{i_{k}}p_{\parallel}}{\epsilon(\bm{p})^{r+2k-1}}+\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}.$ (A.6)

### The third term in Eq. (2.12)

The term proportional to the $(\dot{\alpha}+\dot{\sigma})$ is written as,

$\displaystyle
2g(\dot{\alpha}+\dot{\sigma})\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{p_{i_{1}},\cdots,p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}\bm{p}\cdot\frac{\partial\mathcal{F}_{p}}{\partial\bm{p}}=-(r+2k-1)(\dot{\alpha}+\dot{\sigma})m^{2}\mathcal{J}^{(k,r+2)}_{i_{1},\cdots,i_{k}}-(4-k-r)(\dot{\alpha}+\dot{\sigma})\mathcal{J}^{(k,r)}_{i_{1},\cdots,i_{k}},$
(A.7)

where we used the relations,

$\displaystyle
2g(\dot{\alpha}+\dot{\sigma})\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{p_{i_{1}},\cdots,p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}\bm{p}\cdot\frac{\partial\mathcal{F}_{p}}{\partial\bm{p}}=-2g(\dot{\alpha}+\dot{\sigma})\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{\partial}{\partial\bm{p}}\left(\frac{p_{i_{1}}\cdots
p_{i_{k}}\bm{p}}{\epsilon(\bm{p})^{r+2k-1}}\right)\mathcal{F}_{p},$ (A.8)

and

$\displaystyle\frac{\partial}{\partial\bm{p}}\left(\frac{p_{i_{1}}\cdots
p_{i_{k}}\bm{p}}{\epsilon(\bm{p})^{r+2k-1}}\right)$
$\displaystyle=(k+3)\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}-(r+2k-1)|\bm{p}|^{2}\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k+1}},$
$\displaystyle=(k+3)\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}-(r+2k-1)(\epsilon^{2}-m^{2})\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k+1}},$
$\displaystyle=-(r+k-4)\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k-1}}+(r+2k-1)m^{2}\frac{p_{i_{1}}\cdots
p_{i_{k}}}{\epsilon(\bm{p})^{r+2k+1}}.$ (A.9)

Thus we obtain Eq. (2.2) by using these discussions.

## References

  * [1] F. Sauter, Über das Verhalten eines Elektrons im homogenen elektrischen Feld nach der relativistischen Theorie Diracs (On the behavior of an electron in the homogeneous electric field according to the relativistic theory of Dirac), Z. Phys. 69, 742 (1931).
  * [2] W. Heisenberg and H. Euler, Folgerungen aus der Diracschen Theorie des Positrons (Conclusions from Dirac’s theory of the positron), Z. Phys. 98, 714 (1936).
  * [3] J. S. Schwinger, Phys. Rev. 82, 664-679 (1951) doi:10.1103/PhysRev.82.664
  * [4] G. V. Dunne, Eur. Phys. J. D 55, 327-340 (2009) doi:10.1140/epjd/e2009-00022-0 [arXiv:0812.3163 [hep-th]].
  * [5] A. Ringwald, Phys. Lett. B 510, 107-116 (2001) doi:10.1016/S0370-2693(01)00496-8 [arXiv:hep-ph/0103185 [hep-ph]].
  * [6] A. Casher, H. Neuberger and S. Nussinov, Phys. Rev. D 20, 179-188 (1979) doi:10.1103/PhysRevD.20.179
  * [7] R. Ruffini, L. Vitagliano and S. S. Xue, Phys. Lett. B 559, 12-19 (2003) doi:10.1016/S0370-2693(03)00299-5 [arXiv:astro-ph/0302549 [astro-ph]].
  * [8] A. Yatabe and S. Yamada, doi:10.3847/1538-4357/aa97d8 [arXiv:1712.03960 [astro-ph.HE]].
  * [9] G. W. Gibbons, Commun. Math. Phys. 44, 245-264 (1975) doi:10.1007/BF01609829
  * [10] C. M. Chen, S. P. Kim, I. C. Lin, J. R. Sun and M. F. Wu, Phys. Rev. D 85, 124041 (2012) doi:10.1103/PhysRevD.85.124041 [arXiv:1202.3224 [hep-th]].
  * [11] B. Ratra, Astrophys. J. Lett. 391, L1-L4 (1992) doi:10.1086/186384
  * [12] K. Bamba and J. Yokoyama, Phys. Rev. D 69, 043507 (2004) doi:10.1103/PhysRevD.69.043507 [arXiv:astro-ph/0310824 [astro-ph]].
  * [13] T. Kobayashi, JCAP 05, 040 (2014) doi:10.1088/1475-7516/2014/05/040 [arXiv:1403.5168 [astro-ph.CO]].
  * [14] T. Fujita and R. Namba, Phys. Rev. D 94, no.4, 043523 (2016) doi:10.1103/PhysRevD.94.043523 [arXiv:1602.05673 [astro-ph.CO]].
  * [15] S. Vilchinskii, O. Sobol, E. Gorbar and I. Rudenok, Phys. Rev. D 95, no.8, 083509 (2017) doi:10.1103/PhysRevD.95.083509 [arXiv:1702.02774 [astro-ph.CO]].
  * [16] M. S. Turner and L. M. Widrow, Phys. Rev. D 37, 2743 (1988) doi:10.1103/PhysRevD.37.2743
  * [17] W. D. Garretson, G. B. Field and S. M. Carroll, Phys. Rev. D 46, 5346-5351 (1992) doi:10.1103/PhysRevD.46.5346 [arXiv:hep-ph/9209238 [hep-ph]].
  * [18] G. B. Field and S. M. Carroll, Phys. Rev. D 62, 103008 (2000) doi:10.1103/PhysRevD.62.103008 [arXiv:astro-ph/9811206 [astro-ph]].
  * [19] M. M. Anber and L. Sorbo, JCAP 10, 018 (2006) doi:10.1088/1475-7516/2006/10/018 [arXiv:astro-ph/0606534 [astro-ph]].
  * [20] T. Fujita, R. Namba, Y. Tada, N. Takeda and H. Tashiro, JCAP 05, 054 (2015) doi:10.1088/1475-7516/2015/05/054 [arXiv:1503.05802 [astro-ph.CO]].
  * [21] C. Caprini and L. Sorbo, JCAP 10, 056 (2014) doi:10.1088/1475-7516/2014/10/056 [arXiv:1407.2809 [astro-ph.CO]].
  * [22] C. Caprini, M. C. Guzzetti and L. Sorbo, Class. Quant. Grav. 35, no.12, 124003 (2018) doi:10.1088/1361-6382/aac143 [arXiv:1707.09750 [astro-ph.CO]].
  * [23] T. Fujita and R. Durrer, JCAP 09, 008 (2019) doi:10.1088/1475-7516/2019/09/008 [arXiv:1904.11428 [astro-ph.CO]].
  * [24] S. Okano and T. Fujita, JCAP 03, 026 (2021) doi:10.1088/1475-7516/2021/03/026 [arXiv:2005.13833 [astro-ph.CO]].
  * [25] J. Garriga, Phys. Rev. D 49, 6343-6346 (1994) doi:10.1103/PhysRevD.49.6343
  * [26] J. Martin, Lect. Notes Phys. 738, 193-241 (2008) doi:10.1007/978-3-540-74353-8_6 [arXiv:0704.3540 [hep-th]].
  * [27] M. B. Fröb, J. Garriga, S. Kanno, M. Sasaki, J. Soda, T. Tanaka and A. Vilenkin, JCAP 04, 009 (2014) doi:10.1088/1475-7516/2014/04/009 [arXiv:1401.4137 [hep-th]].
  * [28] R. G. Cai and S. P. Kim, JHEP 09, 072 (2014) doi:10.1007/JHEP09(2014)072 [arXiv:1407.4569 [hep-th]].
  * [29] T. Kobayashi and N. Afshordi, JHEP 10, 166 (2014) doi:10.1007/JHEP10(2014)166 [arXiv:1408.4141 [hep-th]].
  * [30] C. Stahl, E. Strobel and S. S. Xue, Phys. Rev. D 93, no.2, 025004 (2016) doi:10.1103/PhysRevD.93.025004 [arXiv:1507.01686 [gr-qc]].
  * [31] E. Bavarsad, C. Stahl and S. S. Xue, Phys. Rev. D 94, no.10, 104011 (2016) doi:10.1103/PhysRevD.94.104011 [arXiv:1602.06556 [hep-th]].
  * [32] T. Hayashinaka, T. Fujita and J. Yokoyama, JCAP 07, 010 (2016) doi:10.1088/1475-7516/2016/07/010 [arXiv:1603.04165 [hep-th]].
  * [33] V. Domcke, Y. Ema, K. Mukaida and R. Sato, JHEP 03, 111 (2019) doi:10.1007/JHEP03(2019)111 [arXiv:1812.08021 [hep-ph]].
  * [34] H. Kitamoto, Phys. Rev. D 98, no.10, 103512 (2018) doi:10.1103/PhysRevD.98.103512 [arXiv:1807.03753 [hep-th]].
  * [35] V. Domcke, Y. Ema and K. Mukaida, JHEP 02, 055 (2020) doi:10.1007/JHEP02(2020)055 [arXiv:1910.01205 [hep-ph]].
  * [36] O. O. Sobol, E. V. Gorbar, A. I. Momot and S. I. Vilchinskii, Phys. Rev. D 102, no.2, 023506 (2020) doi:10.1103/PhysRevD.102.023506 [arXiv:2004.12664 [gr-qc]].
  * [37] H. Kitamoto, Phys. Rev. D 103, no.6, 063521 (2021) doi:10.1103/PhysRevD.103.063521 [arXiv:2010.10388 [hep-th]].
  * [38] V. Domcke, Y. Ema and K. Mukaida, JHEP 05, 001 (2021) doi:10.1007/JHEP05(2021)001 [arXiv:2101.05192 [hep-ph]].
  * [39] J. Rau, Phys. Rev. D 50, 6911-6916 (1994) doi:10.1103/PhysRevD.50.6911 [arXiv:hep-ph/9402256 [hep-ph]].
  * [40] S. M. Schmidt, D. Blaschke, G. Ropke, S. A. Smolyansky, A. V. Prozorkevich and V. D. Toneev, Int. J. Mod. Phys. E 7, 709-722 (1998) doi:10.1142/S0218301398000403 [arXiv:hep-ph/9809227 [hep-ph]].
  * [41] Y. Kluger, E. Mottola and J. M. Eisenberg, Phys. Rev. D 58, 125015 (1998) doi:10.1103/PhysRevD.58.125015 [arXiv:hep-ph/9803372 [hep-ph]].
  * [42] J. C. R. Bloch, V. A. Mizerny, A. V. Prozorkevich, C. D. Roberts, S. M. Schmidt, S. A. Smolyansky and D. V. Vinnik, Phys. Rev. D 60, 116011 (1999) doi:10.1103/PhysRevD.60.116011 [arXiv:nucl-th/9907027 [nucl-th]].
  * [43] C. K. Dumlu, Phys. Rev. D 79, 065027 (2009) doi:10.1103/PhysRevD.79.065027 [arXiv:0901.2972 [hep-th]].
  * [44] S. P. Kim and C. Schubert, Phys. Rev. D 84, 125028 (2011) doi:10.1103/PhysRevD.84.125028 [arXiv:1110.0900 [hep-th]].
  * [45] A. Huet, S. P. Kim and C. Schubert, Phys. Rev. D 90, no.12, 125033 (2014) doi:10.1103/PhysRevD.90.125033 [arXiv:1411.3074 [hep-th]].
  * [46] E. V. Gorbar, A. I. Momot, O. O. Sobol and S. I. Vilchinskii, Phys. Rev. D 100, no.12, 123502 (2019) doi:10.1103/PhysRevD.100.123502 [arXiv:1909.10332 [gr-qc]].
  * [47] O. O. Sobol, E. V. Gorbar, M. Kamarpour and S. I. Vilchinskii, Phys. Rev. D 98, no.6, 063534 (2018) doi:10.1103/PhysRevD.98.063534 [arXiv:1807.09851 [hep-ph]].
  * [48] M. a. Watanabe, S. Kanno and J. Soda, Phys. Rev. Lett. 102, 191302 (2009) doi:10.1103/PhysRevLett.102.191302 [arXiv:0902.2833 [hep-th]].
  * [49] M. H. Thoma, Rev. Mod. Phys. 81, 959-968 (2009) doi:10.1103/RevModPhys.81.959 [arXiv:0801.0956 [physics.plasm-ph]].
  * [50] A. A. Starobinsky, Phys. Lett. B 91, 99-102 (1980) doi:10.1016/0370-2693(80)90670-X
  * [51] T. Fujita and I. Obata, JCAP 01, 049 (2018) doi:10.1088/1475-7516/2018/01/049 [arXiv:1711.11539 [astro-ph.CO]].
  * [52] S. M. Schmidt, D. Blaschke, G. Ropke, A. V. Prozorkevich, S. A. Smolyansky and V. D. Toneev, Phys. Rev. D 59, 094005 (1999) doi:10.1103/PhysRevD.59.094005 [arXiv:hep-ph/9810452 [hep-ph]].


# CGAR: Critic Guided Action Redistribution in Reinforcement Leaning

Tairan Huang  SCSE
Beihang University
Beijing, China
<EMAIL_ADDRESS>Xu Li  Cognitive Computing Lab
Baidu Research
Beijing, China
<EMAIL_ADDRESS>Hao Li  Mingming Sun  ECE
Peking University
Beijing, China
<EMAIL_ADDRESS>Cognitive Computing Lab
Baidu Research
Beijing, China
<EMAIL_ADDRESS>Ping Li  Cognitive Computing Lab
Baidu Research
Seattle, USA
<EMAIL_ADDRESS>
###### Abstract

Training a game-playing reinforcement learning agent requires multiple
interactions with the environment. Ignorant random exploration may cause a
waste of time and resources. It’s essential to alleviate such waste. As
discussed in this paper, under the settings of the off-policy actor critic
algorithms, we demonstrate that the critic can bring more expected discounted
rewards than or at least equal to the actor. Thus, the Q value predicted by
the critic is a better signal to redistribute the action originally sampled
from the policy distribution predicted by the actor. This paper introduces the
novel Critic Guided Action Redistribution (CGAR) algorithm and tests it on the
OpenAI MuJoCo tasks. The experimental results demonstrate that our method
improves the sample efficiency and achieves state-of-the-art performance. Our
code can be found at https://github.com/tairanhuang/CGAR.

###### Index Terms:

Reinforcement Learning, Soft actor-critic, MuJoCo tasks

##  I Introduction

In recent years, reinforcement learning has been widely used in games and has
made excellent progress in Atari, StarCraft, Dota2, Honor of Kings, and other
games [1, 2]. However, training a reinforcement learning model is time-
consuming due to the massive interactions between the learning and the
environment [3, 4, 5]. Also, interactions with naive exploration strategies
slow down the model’s learning speed and waste resources during model training
[6, 7]. These deficiencies restrict the applications of reinforcement learning
in games.

Off-policy algorithms store the experience in the buffer and reuse them to
reduce interaction costs. Actor critic algorithms use an actor to select
actions and use a critic to estimate the value function to reduce the variance
of policy gradient and accelerate the convergence. Deep deterministic policy
gradient (DDPG) [8] and Soft Actor-Critic (SAC) algorithm [9] combine the
advantages of off-policy algorithms and actor critic algorithms. DDPG learns a
critic and an actor at the same time. It uses the Bellman function to optimize
the critic and uses the critic to optimize the actor. SAC [9] maximizes both
expected return and entropy [10, 11] to balance exploration and exploitation.
The actor of SAC predicts the action distribution and then samples action in
the training stage and uses the mean value of the distribution as the action
in the evaluation stage. It further improves the sample efficiency and
stability of reinforcement learning.

Following the same purpose, we propose a novel Critic Guided Action
Redistribution (CGAR) mechanism and show that SAC with CGAR achieves state-of-
the-art performance. In this paper, we first give a theoretical analysis of
actor critic algorithms and demonstrate that the critic can bring more
expected discounted rewards than or at least equal to the actor in off-policy
actor critic algorithms. To utilize such advancement of critic, we
redistribute the action distribution predicted by the actor through the Q
value predicted by the critic.

Figure 1: CGAR: Critic Guided Action Redistribution in Reinforcement Leaning

In the proposed CGAR model, after getting an action distribution predicted by
the actor, we sample $K$ actions from that action distribution. Then, the
critic model predicts the Q values for these $K$ actions conditioned on the
state. We reset the selection probability of these actions based on their
corresponding Q values. Since the critic model is optimized to predict the
expected return, with the distribution positively correlated to the Q value,
the algorithm tends to select the action with a higher expected return. We
deploy our method on SAC and conduct experiments on the OpenAI MuJoCo tasks
[12]. The experimental results demonstrate that our approach is effective and
achieves state-of-the-art performance. Here we summarize our main
contributions:

  * •

We demonstrate that the critic can bring more expected discounted rewards than
or at least equal to the actor in the off-policy actor critic algorithm.

  * •

We propose a novel Critic Guided Action Redistribution (CGAR), which uses the
Q value predicted by the critic to resample action from the action
distribution predicted by the actor.

  * •

We apply our method to SAC and achieve state-of-the-art performance on OpenAI
MuJoCo tasks.

##  II Algorithm

###  II-A Motivation

Under the classic setting of the off-policy actor critic learning procedure,
at each interaction with the environment (environment step), the agent
collects data from the environment under the policy network with the action
distribution of ${\pi_{\phi_{a}(s)}}$. Then, given the updated data buffer
$\mathcal{D}$, the critic is optimized to estimate the future reward. Then,
the actor is updated to maximize the estimated expected future reward given
the currently learned critic. The procedure forms a dependency circle and
loops at each environment step during training.

At a certain environment step, suppose that the critic ${Q_{\phi_{c}}}$ is
optimized to ${Q_{\phi_{c}^{i}}}$ after $i$ gradient step under the currently
collected data $\mathcal{D}$. The loss function to optimize actor is generally
written as,

$J_{\pi}(\phi_{a})=\mathbb{E}_{s_{t}\sim\mathcal{D}}\mathbb{E}_{a_{t}\sim\mathcal{\pi}_{\phi_{a}}(a_{t}|s_{t})}(-Q_{\phi_{c}^{i}}(s_{t},a_{t})).$
(1)

When the loss function is minimized, the action with the highest probability
predicted by the actor is the same as the action leading to the highest
${Q_{\phi_{c}^{i}}}$ predicted by the critic. Then the performance of the
actor and the critic is the same. However, most of the time, the actor is less
optimized to the optimal distribution of action given the current critic,
which leads to poor performance compared with the current critic.

This learning procedure is similar to the knowledge distillation or teacher-
student method [13]. In our situation, we use the Q value output by the critic
to calculate the loss function to train the actor. Hence, critic and actor
correspond to teacher and student, respectively. The performance gap between
the student and teacher has been demonstrated in many previous works, and the
gap exists even though the student network has the same size as the teacher
network [14, 15, 16]. In Section II-B, we also give an empirical demonstration
of the motivation. As the performance of the critic is better than or at least
equal to the actor in each environment step, we can expect that the critic can
bring more expected discounted rewards than or at least equal to the actor
during the RL training procedure. This paper proposes a CGAR algorithm to
ameliorate the performance gap between the actor and critic.

###  II-B Empirical Demonstration of the Motivation

We conduct experiments under supervised learning settings to empirically
demonstrate the above motivation. In detail, We use MNIST [17] as a dataset
and maintain two identical multilayer perceptron models, $M_{1}$ and $M_{2}$.
We let $M_{1}$ fit the dataset and let $M_{2}$ fit $M_{1}$. Given input $x$
and its label $y$, the loss function of $M_{1}$ is:
${\mathcal{L}_{1}=CrossEntropy(M_{1}(x),y)}$, and the loss function of $M_{2}$
is: ${\mathcal{L}_{2}=CrossEntropy(M_{2}(x),M_{1}(x))}$. Under this setting,
$M_{1}$ is the signal provided to $M_{2}$, and $y$ is the signal provided to
$M_{1}$. Besides, $M_{1}$ and $M_{2}$ are updated one after another
iteratively. This setting is quite similar to the off-policy actor critic
algorithm in Section II-A, where the reward is provided to calculate the
critic’s target, and the Q is the signal provided to the actor. To simulate
the loss function of critic, we further design another loss function,
${\mathcal{L}_{1}=MSE(M_{1}(x),\bar{y})}$, where $\bar{y}$ is the one-hot
representation of $y$, and ${\mathcal{L}_{2}}$ remains unchanged. We test the
evaluation accuracy of $M_{1}$ and $M_{2}$ every epoch over five seeds and
report the results in Fig. 2. We can see that under both loss functions,
$M_{2}$ learns more slowly than $M_{1}$, demonstrating our motivation’s
correctness.

Figure 2: Comparison results of $M_{1}$ and $M_{2}$.

###  II-C Critic Guided Action Redistribution

We propose our Critic Guided Action Redistribution (CGAR) based on the above
motivation and empirical demonstration. In actor critic algorithms, given
state ${s_{t}}$ at environment step $t$, the actor ${\pi_{\phi_{a}}}$ predicts
the action distribution ${\mathcal{P}_{t}}$. After that, it samples action
${a_{t}}$ from ${\mathcal{P}_{t}}$ to interact with the environment.

$a_{t}\sim\mathcal{P}_{t}=\pi_{\phi_{a}}(s_{t}).$ (2)

In our algorithm, we first sample $K$ actions
${\\{a_{t}^{0},a_{t}^{1},...,a_{t}^{K-1}\\}}$ from $\mathcal{P}_{t}$. These
actions construct the actions set ${\mathbb{L}_{a}^{t}}$.

$\mathbb{L}_{a}^{t}:=\\{a_{t}^{i}\sim\mathcal{P}_{t}|i\in[0,K-1]\\}.$ (3)

Then we use the critic to predict the Q value for every action in
${\mathbb{L}_{a}^{t}}$. The Q value ${Q_{t}^{i}}$ is calculated on every
action ${a_{t}^{i}}$ in ${\mathbb{L}_{a}^{t}}$ conditioned with state $s_{t}$.
These Q values ${\\{Q_{t}^{0},Q_{t}^{1},...,Q_{t}^{K-1}\\}}$ construct the Q
value set ${\mathbb{L}_{Q}^{t}}$.

$\mathbb{L}_{Q}^{t}:=\\{Q_{t}^{i}=Q_{\phi_{c}}(s_{t},a_{t}^{i})|i\in[0,K-1]\\}.$
(4)

We use the Q value ${Q_{t}^{i}}$ to calculate the probability ${p_{t}^{i}}$ of
${a_{t}^{i}}$. The probability set results from the Q value set calculated by
the Softmax function. In this way, we could get the action probability
distribution $\mathcal{P}_{t}^{\prime}$.

$\displaystyle\mathcal{P}_{t}^{\prime}$
$\displaystyle=\text{Softmax}(\mathbb{L}_{Q}^{t})$ (5)
$\displaystyle=\\{\frac{e^{Q_{t}^{i}}}{\sum_{j=0}^{K-1}e^{Q_{t}^{j}}}|i\in[0,K-1]\\}.$

We use the Softmax function to make actions with large Q values more likely to
be sampled while maintaining a certain degree of exploration. Finally, we
select the action ${a_{t}}$ from the new distribution. Fig. 1 is our model’s
diagram.

$a_{t}\sim\mathcal{P}^{\prime}_{t}.$ (6)

###  II-D Implementation

We apply CGAR to SAC. SAC is an off-policy algorithm based on the Maximum
Entropy Principle. Its optimization goal is simultaneously maximizing both the
expected return and the entropy. It learns a policy ${\pi_{\phi_{a}}}$, a Q
value function ${Q_{\phi_{c}}}$, and a temperature coefficient ${\alpha}$,
with parameters ${\phi_{a}}$, ${\phi_{c}=(\phi_{c}^{1},\phi_{c}^{2})}$, and
${\alpha}$ separately. The loss functions of SAC are defined below, which are
introduced in the SAC paper.

$\displaystyle J_{Q}(\phi_{c})=$
$\displaystyle\mathbb{E}_{s_{t},a_{t}\sim\mathcal{D}}\frac{1}{2}(Q_{\phi_{c}}(s_{t},a_{t})$
(7) $\displaystyle-$ $\displaystyle(r_{t}+\gamma\mathbb{E}_{s_{t+1}\sim
p}(V_{\bar{\phi}_{c}}(s_{t+1}))))^{2},$ $\displaystyle
V_{\bar{\phi}_{c}}(s_{t})=$
$\displaystyle\mathbb{E}_{a_{t}\sim\pi_{\phi_{a}}(a_{t}|s_{t})}Q_{\bar{\phi}_{c}}(s_{t},a_{t})$
(8) $\displaystyle-$ $\displaystyle\alpha\log\pi_{\phi_{a}}(a_{t}|s_{t}),$
$\displaystyle J_{\pi}(\phi_{a})=$
$\displaystyle\mathbb{E}_{s_{t}\sim\mathcal{D}}\mathbb{E}_{a_{t}\sim\mathcal{\pi}_{\phi_{a}}(a_{t}|s_{t})}$
(9)
$\displaystyle(\alpha\log\pi_{\phi_{a}}(a_{t}|s_{t})-Q_{\phi_{c}}(s_{t},a_{t})).$
$J_{\alpha}(\alpha)=\mathbb{E}_{a_{t}\sim\mathcal{\pi}_{\phi_{a}}(a_{t}|s_{t})}-\alpha\log\pi_{\phi_{a}}(a_{t}|s_{t})-\alpha\mathcal{H},$
(10)

Note that different from Eq. (1), there’s an entropy term in Eq. (9), which
doesn’t affect the demonstration in Section II-A. Our complete algorithm is
shown in Algorithm 1. Our unique operations are marked with red. And the
deleted operations in SAC are marked with blue and strikethrough.

Algorithm 1 CGAR applied to SAC

1: Input: env

2: Output: $\phi_{c}$, $\phi_{a}$ and $\alpha$

3: Initialize parameters $\phi_{c}$, $\phi_{a}$ and $\alpha$

4: Set $\bar{\phi}_{c}=\phi_{c}$, $\mathcal{D}=\emptyset$

5: for $k=1,\dots,N_{\rm init}$ do

6: Sample random action: ${a_{t}\sim\pi_{random}(s_{t})}$

7: Execute action: ${s_{t+1},r_{t+1},\rm done\sim\textsc{env}(a_{t})}$

8: Collect data:
${\mathcal{D}\leftarrow\mathcal{D}\cup\\{s_{t},a_{t},s_{t+1},r_{t+1}\\}}$

9: end for

10: for $k=1,\dots,N_{\rm train}$ do

11: Predict action distribution: ${\mathcal{P}_{t}=\pi(s_{t})}$

12: Sample action: ${a_{t}\sim\mathcal{P}_{t}}$

13: Sample ${K}$ actions from ${\mathcal{P}_{t}}$:
${\mathbb{L}_{a}^{t}=\\{a_{t}^{i}\sim\mathcal{P}_{t}|i\in[0,K-1]\\}}$

14: Calculate the Q value
set:${\mathbb{L}_{Q}^{t}=\\{Q_{t}^{i}=Q_{\phi_{c}}(s_{t},a_{t}^{i})|i\in[0,K-1]\\}}$

15: Get the new distribution:
${\mathcal{P}^{\prime}_{t}=\text{Softmax}(\mathbb{L}_{Q}^{t})}$

16: Sample action: ${a_{t}\sim\mathcal{P}^{\prime}_{t}}$

17: Execute action: ${s_{t+1},r_{t+1},\rm done\sim\textsc{env}(a_{t})}$

18: Collect data:
${\mathcal{D}\leftarrow\mathcal{D}\cup\\{s_{t},a_{t},s_{t+1},r_{t+1}\\}}$

19: Train critic:
${\phi_{c}\leftarrow\phi_{c}-\lambda_{Q}\nabla_{\phi_{c}}J_{Q}(\phi_{c})}$

20: Train actor:
${\phi_{a}\leftarrow\phi_{a}-\lambda_{\pi}\nabla_{\phi_{a}}J_{\pi}(\phi_{a})}$

21: Train alpha:
${\alpha\leftarrow\alpha-\lambda_{\alpha}\nabla_{\alpha}J_{\alpha}(\alpha)}$

22: Update the target critic:
${\bar{\phi}_{c}\leftarrow\tau\phi_{c}+(1-\tau)\bar{\phi}_{c}}$

23: end for

##  III Experiment

###  III-A Implementation Details

We implement CGAR on the SAC algorithm noted as CGAR-SAC. The implementation
of the SAC comes from [18]. We evaluate our method on OpenAI MuJoCo tasks,
including tasks such as standing, walking, and running. The state of the agent
consists of parameters such as positions and velocities. Action is a real-
valued vector that represents the control of the agent’s joints. The purpose
of model learning is to maximize the expected discounted rewards. We test the
agent performance every 10,000 environment steps. We compute the mean episode
returns an agent obtains over ten episodes for every evaluation. All results
are over five different seeds, and we keep the minimum, maximum, and mean
values over these seeds.

###  III-B Comparison between CGAR-SAC with SAC

TABLE I: Mean value of the average return for the whole training process Task | BCC | CR | RE | WW | WS | FS
---|---|---|---|---|---|---
CGAR-SAC | ${\bf{843_{\pm 73}}}$ | ${\bf{593_{\pm 51}}}$ | ${\bf{821_{\pm 81}}}$ | ${\bf{761_{\pm 71}}}$ | ${{{\bf{857}}_{\pm 32}}}$ | ${\bf{793_{\pm 47}}}$
SAC | ${{841_{\pm 89}}}$ | ${{574_{\pm 59}}}$ | ${{817_{\pm 82}}}$ | ${{746_{\pm 110}}}$ | ${{846_{\pm\bf{29}}}}$ | ${{760_{\pm 49}}}$

We compare the performance of CGAR-SAC with SAC. We counted the mean, maximum
and minimum values of the average return obtained overall seeds during each
evaluation and plotted them in Fig. 3. The curve represents the average
return, and the shading represents the range between the maximum and minimum.
From Fig. 3 we can see that in most tasks, CGAR-SAC converges faster and
achieves better final performance than SAC, especially Cheetah Run, Walker
Walk, and Finger Spin. And CGAR-SAC is not weaker than SAC in other tasks. We
also counted the average return obtained overall seeds during each evaluation
and calculated their mean over the entire training process. In Table I, we use
the initials of the first letter of every environment as the table title, and
the full names can be found in Fig. 3. The value of each item in the table
represents the mean value of the average return, and the subscript denotes the
standard deviation. From Table I, we can see that the mean value of the
average return of our method during the training process is higher than SAC in
every task, and the standard deviation is lower than SAC in most tasks. We can
conclude that CGAR improves the sample efficiency of SAC.

Figure 3: Comparison results between CGAR-SAC and SAC.

##  IV Conclusion

This paper proposes a novel action redistribution algorithm, Critic Guided
Action Redistribution for game playing. We demonstrate that the critic can
bring more expected discounted rewards than or at least equal to the actor in
the off-policy actor critic algorithm. Based on the demonstration, we use the
Q value predicted by the critic to redistribute the actions probability
distribution generated by the actor. Then we sample actions from the new
distribution to interact with the environment. We implement our algorithm on
SAC and test it on the OpenAI MuJoCo tasks. The experimental results
demonstrate that our method improves the sample efficiency and achieves state-
of-the-art performance. Future research can be done by applying CGAR to other
games or analyzing the distribution of Q value on the multimodal distribution.

## References

  * [1] D. Ye, G. Chen, W. Zhang, S. Chen, B. Yuan, B. Liu, J. Chen, Z. Liu, F. Qiu, H. Yu, Y. Yin, B. Shi, L. Wang, T. Shi, Q. Fu, W. Yang, L. Huang, and W. Liu, “Towards playing full MOBA games with deep reinforcement learning,” in _NeurIPS 2020, December 6-12, 2020, virtual_.
  * [2] N. Brown, A. Bakhtin, A. Lerer, and Q. Gong, “Combining deep reinforcement learning and search for imperfect-information games,” in _NeurIPS 2020, December 6-12, 2020, virtual_.
  * [3] R. S. Sutton and A. G. Barto, “Reinforcement learning: An introduction,” 2011\.
  * [4] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski _et al._ , “Human-level control through deep reinforcement learning,” _nature_ , vol. 518, no. 7540, pp. 529–533, 2015.
  * [5] T. Malloy, C. R. Sims, T. Klinger, M. Liu, M. Riemer, and G. Tesauro, “Capacity-limited decentralized actor-critic for multi-agent games,” in _COG 2021, Copenhagen, Denmark, August 17-20, 2021_. IEEE, 2021, pp. 1–8.
  * [6] R. Liang, Y. Zhu, Z. Tang, M. Yang, and X. Zhu, “Proximal policy optimization with elo-based opponent selection and combination with enhanced rolling horizon evolution algorithm,” in _COG 2021, Copenhagen, Denmark, August 17-20, 2021_. IEEE, 2021, pp. 1–4.
  * [7] J. Büttner and S. von Mammen, “Training a reinforcement learning agent based on XCS in a competitive snake environment,” in _COG 2021, Copenhagen, Denmark, August 17-20, 2021_. IEEE, 2021, pp. 1–5.
  * [8] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” in _ICLR_ , 2016.
  * [9] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _ICML_ , 2018.
  * [10] T. Haarnoja, H. Tang, P. Abbeel, and S. Levine, “Reinforcement learning with deep energy-based policies,” in _ICML 2017, Sydney, NSW, Australia, 6-11 August 2017_. PMLR, 2017, pp. 1352–1361.
  * [11] B. D. Ziebart, A. L. Maas, J. A. Bagnell, and A. K. Dey, “Maximum entropy inverse reinforcement learning,” in _AAAI 2008, Chicago, Illinois, USA, July 13-17, 2008_. AAAI Press, pp. 1433–1438.
  * [12] Y. Tassa, Y. Doron, A. Muldal, T. Erez, Y. Li, D. de Las Casas, D. Budden, A. Abdolmaleki, J. Merel, A. Lefrancq, T. P. Lillicrap, and M. A. Riedmiller, “Deepmind control suite,” _CoRR_ , vol. abs/1801.00690, 2018.
  * [13] G. Hinton, O. Vinyals, and J. Dean, “Distilling the knowledge in a neural network,” in _NIPS Deep Learning and Representation Learning Workshop_ , 2015\.
  * [14] S. I. Mirzadeh, M. Farajtabar, A. Li, N. Levine, A. Matsukawa, and H. Ghasemzadeh, “Improved knowledge distillation via teacher assistant,” in _AAAI 2020_ , vol. 34, no. 04, 2020, pp. 5191–5198.
  * [15] J. H. Cho and B. Hariharan, “On the efficacy of knowledge distillation,” in _ICCV 2019_.
  * [16] X. Deng and Z. Zhang, “Can students outperform teachers in knowledge distillation based model compression?”
  * [17] Y. LeCun and C. Cortes, “MNIST handwritten digit database,” 2010.
  * [18] D. Yarats and I. Kostrikov, “Soft actor-critic (sac) implementation in pytorch,” https://github.com/denisyarats/pytorch_sac.