aainstitutetext: Faculty of Science and Technology, University of Stavanger,
4036 Stavanger, Norway

# Quantum tachyonic preheating, revisited

Anders Tranberg a Gerhard Ungersbäck<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In certain models of inflation, the postinflationary reheating of the Universe
is not primarily due to perturbative decay of the inflaton field into
particles, but proceeds through a tachyonic instability. In the process, long-
wavelength modes of an unstable field, which is often distinct from the
inflaton itself, acquire very large occupation numbers, which are subsequently
redistributed into a thermal equilibrium state. We investigate this process
numerically through quantum real-time lattice simulations of the Kadanoff-Baym
equation, using a $1/N$-NLO truncation of the 2PI-effective action. We
identify the early-time maximum occupation number, the ”classical” momentum
range, the validity of the classical approximation and the effective IR
temperature, and study the kinetic equilibration of the system and the
equation of state.

###### Keywords:

##  1 Introduction

At the end of inflation PhysRevD.23.347 ; Liddle:2000cg , most of the energy
in the Universe resided in the potential energy of an inflaton field. This
energy is expected to have been transferred to the Standard Model degrees of
freedom through one or a combination of processes collectively referred to as
reheating PhysRevLett.73.3195 . One of these is tachyonic preheating
PhysRevLett.87.011601 , which occurs if at the end of inflation a field (the
inflaton itself, but more commonly another field coupled to it) becomes
unstable and undergoes a spinodal transition. In the process, the unstable
field acquires very large occupation numbers in the long-wavelength (IR)
modes.

The most common realisation of this scenario is a two-field model
PhysRevD.49.748 , where one field $\sigma$ plays the role of the inflaton and
a second field111The ”waterfall” field. $\phi$ undergoes spontaneous a
symmetry breaking transition at zero temperature triggered by the
(slow-)rolling of the inflaton. A generic model for this follows from the
action

$\displaystyle S=\int
d^{4}x\left[\frac{1}{2}\partial_{\mu}\sigma\partial^{\mu}\sigma+V(\sigma)+\frac{1}{2}\partial_{\mu}\phi\,\partial^{\mu}\phi+\frac{1}{2}(g^{2}\sigma^{2}-\mu^{2})\phi^{2}+\frac{1}{24}\lambda\phi^{4}+V_{0}\right].$
(1)

where $\sigma$ is real-valued and $\phi$ is here written as a real singlet,
but could have any number of components. The inflaton potential $V(\sigma)$ is
not specified but must have the property that the potential minimum is at
$\sigma=0$ and allow for slow-roll inflation at values of $\sigma$ larger than
$\mu/g$.

As the inflation rolls towards $\sigma=0$, the effective mass of the second
field $\phi$ changes sign

$\displaystyle m^{2}(t)=g^{2}\sigma^{2}(t)-\mu^{2},$ (2)

and $\phi$ undergoes a symmetry breaking transition, as the potential acquires
new minima at non-zero field value

$\displaystyle\phi^{2}\rightarrow v^{2}=6\mu^{2}/\lambda,$ (3)

with a speed determined by the evolution of $\sigma$,

$\displaystyle
u=\frac{1}{\mu^{3}}\frac{dm^{2}(t)}{dt}_{\sigma=\mu/g}=\frac{2g\dot{\sigma}}{\mu^{2}}.$
(4)

We may fix the constant $V_{0}=3\mu^{4}/2\lambda$ in such a way that the
potential is zero at $\phi=v$, $\sigma=0$.

In hybrid inflation models PhysRevD.49.748 , the end of inflation is itself
triggered by the symmetry breaking transition, as the increase of $\phi^{2}$
spoils the flatness of the $\sigma$-potential. Either way, the potential
initially stored in $V_{0}$ is transferred to all the degrees of freedom, and
the system eventually equilibrates to some reheating temperature $T_{\rm reh}$
near the potential minimum.

The process of (p)reheating is often quite complicated and model-dependent,
since in addition to triggering symmetry breaking, the $\sigma$ field itself
may also oscillate around $\sigma=0$ so that the sign of $m^{2}(t)$ flips
multiple times. Such oscillations may also lead to resonant preheating
PhysRevLett.73.3195 , where certain momentum modes in resonance with $\sigma$
grow exponentially. A substantial body of analytic and in particular numerical
work has been done over the past 25 years, exploring many model of preheating
after inflation Garcia-Bellido:2002fsq ; Smit:2002yg ; Arrizabalaga:2004iw ;
Aarts:1997kp ; Dux:2022kuk ; Garcia:2023eol ; Mahbub:2023faw . A review can be
found in Amin:2014eta .

Several aspects of a tachyonic reheating transition deserve detailed scrutiny.
One is the transition itself, where the field dynamics may be solved
numerically in real-time and the spectrum of fluctuations computed. Because
occupation numbers are large, at least for some modes, for some of the time,
most numerical treatments rely on the classical approximation to the dynamics
Aarts:1997kp ; Aarts:2001yn . Yet, the initial state of the system is the
quantum vacuum, and the unstable modes are seeded by the vacuum fluctuations.
Also, at late times classical dynamics is notoriously incomplete, in that the
classical equilibrium is badly defined and does not match the quantum
equilibrium state. To quantify the validity of the classical approximation and
follow the thermalisation process using quantum dynamics seems worthwhile
Arrizabalaga:2004iw .

The approach to equilibrium (kinetic and chemical) also has a set of
characteristic timescales and an effective equation of state which enters in
the cosmological evolution. An out-of-equilibrium regime with very large
occupation numbers may be suitable for other interesting processes to take
place. These include the creation of heavy exotic particles, primordial black
holes, gravitational waves, topological defects and a baryon asymmetry. When
model building in these contexts, it is useful to known how large occupation
numbers become, in which momentum range and over what timescales.

In the simplest analysis, one might assume that the potential energy is
instantaneously redistributed onto $g^{*}$ non-interacting relativistic
degrees of freedom, so that the final reheating temperature is simply

$V_{0}=\frac{\pi^{2}}{30}g^{*}T_{\rm reh}^{4}\rightarrow\frac{T_{\rm
reh}}{\mu}=\left(\frac{45}{\pi^{2}g^{*}\lambda}\right)^{1/4}.$ (5)

In the Standard Model below $T_{\rm reh}=100$ GeV, one expects222With the
appropriate normalisation of $V_{0}$ in the SM, $T_{\rm reh}\simeq 42$ GeV.
$g^{*}\simeq 107$, while for the O(N) model to be discussed below, the
expectation is $g^{*}=N-1$. We will see that at intermediate times, this naive
estimate has limitations.

We may quantify under what conditions we are allowed to ignore the
cosmological expansion. The longest timescale in the problem is the chemical
equilibration time $\tau_{ch}$, and so we can ignore expansion as long as
$H\tau_{ch}<1$. Since

$\displaystyle H^{2}=\frac{V_{0}}{3M_{\rm
pl}^{2}}\rightarrow\frac{H}{\mu}=\frac{1}{\sqrt{2\lambda}}\frac{\mu}{M_{\rm
pl}},$ (6)

as long as $M_{\rm pl}/\mu\gg\mu\tau_{\rm ch}$ and the coupling is not too
small, we can neglect expansion. If not, the evolution equations have to be
restated and solved in a Friedmann-Robertson-Walker metric as in
Tranberg:2008ae .

In the present work, we will ignore expansion and compute the out-of-
equilibrium quantum dynamics of a O(N) scalar field system under a mass
quench, solving the Kadanoff-Baym equations truncated at NLO in a 1/N-2PI
expansion Aarts:2001qa ; Berges:2000ur ; Aarts:2002dj . Some of the work can
be seen as a refinement and continuation of Arrizabalaga:2004iw , and shares
much of the notation.

The paper is organised as follows: In section 2 we set up our simplified model
of spinodal symmetry breaking. In section 3 we introduce the 2PI formalism,
the quantum evolution equations and their classical limit. We present our
numerical results in sections 4, 5 and 6, where the first section focuses on
classical aspects of tachyonic preheating and the applicability of classical
approximations, and the second addresses questions outside the reach of
classical approximations, where quantum equations of motions are essential. We
briefly consider the early time dynamics in finite time mass quenches in
section 6. We conclude in section 7. Some technical aspects of the simulations
are placed in an appendix A.

##  2 Symmetry and symmetry breaking in the O(N) scalar field model

We will consider a N-component scalar field with an $O(N)$-symmetric
potential, described by the following action (sum over components $a=1,...,N$
implied)

$\displaystyle S=\int
d^{4}x\left[\frac{1}{2}\partial_{\mu}\phi_{a}\partial^{\mu}\phi_{a}-\frac{m^{2}(t)}{2}\phi_{a}\phi_{a}-\frac{\lambda}{24N}(\phi_{a}\phi_{a})^{2}.\right]$
(7)

This is a generalisation of the single real field $\phi$ in the preceding
section. For $N=4$ we have a model reminiscent of the Standard Model complex
Higgs doublet. Keeping $N$ general allows us to mimic an arbitrary number of
spinodally unstable degrees of freedom to be reheated at the end of inflation.
We will also use $1/N$ as the expansion parameter to organise our 2PI diagram
series below.

When $m^{2}(t)=-\mu^{2}$, the potential has minima at
$\phi_{a}\phi_{a}=6N\mu^{2}/\lambda\equiv v^{2}$. A common prescription is to
choose a global rotation to fix this in the $\phi_{1}$-direction,
$\phi_{a}=\delta_{a1}v$. Then the mass spectrum becomes

$\displaystyle
M^{2}_{ab}=\left(m^{2}+\frac{\lambda}{6N}\phi_{c}\phi_{c}\right)\delta_{ab}+\frac{\lambda}{3N}\phi_{a}\phi_{b}=\textrm{diag}(2\mu^{2},0,0,...),$
(8)

with a massive Higgs mode with mass $m_{H}^{2}=2\mu^{2}$, corresponding to
fluctuations in the $\phi_{1}$ direction around $v$, and $N-1$ massless
Goldstone modes, corresponding to fluctuations perpendicular to $\phi_{1}$.
These two modes are often explicitly encoded in a longitudinal (along the
1-direction) and perpendicular (to the 1-direction) propagator.

On the other hand, the system has O(N) symmetry, and so the equations of
motion will conserve $\langle\phi_{a}\rangle=0$ for all $a$. We may therefore
choose to study the dynamics keeping this symmetry manifest, and consider only
one ”compound” correlator encoding both massive and massless modes
Arrizabalaga:2004iw . We expect the compound propagator to contain $N-1$
Goldstone modes and one Higgs mode in equilibrium, and that at equal times it
may be decomposed as333This is in contrast to the common procedure of defining
the connected correlator
$G_{c}=\langle\phi\phi\rangle-\langle\phi\rangle^{2}$, and from there the
transverse and longitudinal components of this correlator. Here, symmetry
ensures $\langle\phi\rangle=0$, but the correlator still solves
$dV/d\phi=0=(\lambda\langle\phi^{2}\rangle-\mu^{2})\phi$.

$\displaystyle\langle\phi({\bf x},t)\phi({\bf y},t)\rangle=G({\bf
x-y},t)=\frac{v^{2}}{N}+\frac{N-1}{N}F^{G}({\bf x-y},t)+\frac{1}{N}F^{H}({\bf
x-y},t).$ (9)

It follows that expectation value then enters as a contribution to the zero
momentum mode

$\frac{|\phi_{{\bf k}=0}|^{2}(t)}{L^{3}}\simeq\frac{v^{2}}{N}.$ (10)

We will return to this point below.

###  2.1 Triggering the symmetry breaking transition

The transition is triggered as the function $m^{2}(t)$ changes sign. This may
be parametrised in terms of a finite quench time $\tau_{Q}$, as

$\displaystyle m^{2}(t)$ $\displaystyle=\mu_{0}^{2},\quad t<0,$ (11)
$\displaystyle=\mu_{0}^{2}\left[1-\left(1+\frac{\mu^{2}}{\mu_{0}^{2}}\right)\frac{t}{\tau_{Q}}\right],\quad
0<t<\tau_{Q},$ (12) $\displaystyle=-\mu^{2},\quad t>\tau_{Q}.$ (13)

This corresponds to a quench speed
$u=-\frac{1}{\mu\tau_{Q}}(1+\frac{\mu_{0}^{2}}{\mu^{2}})$. An instantaneous
quench amounts to $\tau_{Q}=0$, cutting out the transition region altogether

$\displaystyle m^{2}(t)$ $\displaystyle=\mu_{0}^{2},\quad t<0,$ (14)
$\displaystyle m^{2}(t)$ $\displaystyle=-\mu^{2},\quad t>0.$ (15)

For simplicity, we will consider ”symmetric” quenches where
$\mu^{2}=\mu_{0}^{2}$. In the general case of a dynamical inflaton
$m^{2}(t)=g^{2}\sigma^{2}-\mu^{2}$, $\mu_{0}^{2}$ corresponds to some initial
value of $\sigma$. In the example of an instantaneous quench, modes with
momentum $k<\mu$ evolve as

$\displaystyle\phi_{k}(t)\propto e^{\pm i\omega_{k}t}\rightarrow
e^{\pm|\mu^{2}-k^{2}|t},$ (16)

which is the tachyonic (or spinodal) instability leading to exponential growth
of the long-wavelength modes. Performing a more detailed computation, using as
initial condition the quantum vacuum prior to the quench
($\pi_{k}=\partial_{t}\phi_{k}$)

$\displaystyle\langle\phi_{k}^{\dagger}\phi_{k}\rangle=\frac{n_{k}+1/2}{\omega_{k}^{+}},\quad\langle\pi_{k}^{\dagger}\pi_{k}\rangle=\omega_{k}^{+}(n_{k}+1/2),$
(17)

one finds

$\displaystyle\langle\phi_{k}^{\dagger}(t)\phi_{k}(t^{\prime})\rangle_{t=t^{\prime}}=F_{k}(t,t^{\prime})|_{t=t^{\prime}}=\frac{1}{2\omega^{+}_{k}}\Biggl{[}1-\biggl{(}\frac{\omega_{k}^{+\,2}}{\omega_{k}^{-\,2}}-1\biggr{)}\sinh^{2}(\omega_{k}^{-}t)\Biggr{]},$
(18)
$\displaystyle\langle\pi_{k}^{\dagger}(t)\pi_{k}(t^{\prime})\rangle_{t=t^{\prime}}=\partial_{t}\partial_{t^{\prime}}F_{k}(t,t^{\prime})|_{t=t^{\prime}}=\frac{\omega^{-\,2}_{k}}{2\omega^{+}_{k}}\Biggl{[}1+\biggl{(}\frac{\omega_{k}^{+\,2}}{\omega_{k}^{-\,2}}-1\biggr{)}\cosh^{2}(\omega_{k}^{-}t)\Biggr{]},$
(19)

where $\omega_{k}^{\pm 2}=\pm\mu^{2}+k^{2}$. For a linear quench, the
corresponding expression may be found in Garcia-Bellido:2002fsq . The
exponential growth continues until the self-interactions become important (see
below).

##  3 Quantum evolution equations

###  3.1 The 2PI effective action

$\Phi^{\rm LO}[G]=\begin{gathered}\leavevmode\hbox to2.8pt{\vbox
to2.8pt{\pgfpicture\makeatletter\hbox{\hskip 1.4pt\lower-1.4pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{ } \feynman \vertex(v);
\vertex[above=1.0cm of v](t); \vertex[below=1.0cm of v](b); \diagram*{ (v) --
[out=180,in=180, looseness=1.75] (t) -- [ out=0,in=0, looseness=1.75] (v) \--
[out=-180,in=180, looseness=1.75] (b) --[ out=0,in=0, looseness=1.75] (v) };
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{gathered}$

$,\quad\Phi^{\rm NLO}[G]=\begin{gathered}\leavevmode\hbox to2.8pt{\vbox
to2.8pt{\pgfpicture\makeatletter\hbox{\hskip 1.4pt\lower-1.4pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{ } \feynman \vertex(v);
\vertex[above=1.0cm of v](t); \vertex[below=1.0cm of v](b); \diagram*{ (v) --
[ out=180,in=180, looseness=1.75] (t) --[ out=0,in=0, looseness=1.75] (v) \--
[ out=-180,in=180, looseness=1.75] (b) --[ out=0,in=0, looseness=1.75] (v) };
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{gathered}+\begin{gathered}\leavevmode\hbox
to2.8pt{\vbox to2.8pt{\pgfpicture\makeatletter\hbox{\hskip
1.4pt\lower-1.4pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{ } \feynman \vertex(v1);
\vertex[right=1.5cm of v1] (v3); \diagram* { (v1) -- [ half left,
looseness=1.7] (v3) -- [ half left, looseness=1.7] (v1) }; \diagram* { (v1) --
[ half left, looseness=0.7] (v3) -- [ half left, looseness=0.7] (v1) };
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{gathered}+\begin{gathered}\leavevmode\hbox
to2.8pt{\vbox to2.8pt{\pgfpicture\makeatletter\hbox{\hskip
1.4pt\lower-1.4pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{ } \feynman \vertex(v1);
\vertex[right=1.5cm of v1] (v3); \vertex[above right=1.06cm of v1] (v2);
\diagram* { (v1) -- [ half left, looseness=1.7] (v3) -- [ half left,
looseness=1.7] (v1) }; \diagram* { (v1) -- [ half right, looseness=0.5] (v2)
-- [ half right, looseness=0.5] (v3) -- [ half left, looseness=0.7] (v1) };
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{gathered}+\begin{gathered}\leavevmode\hbox
to2.8pt{\vbox to2.8pt{\pgfpicture\makeatletter\hbox{\hskip
1.4pt\lower-1.4pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{ } \feynman \vertex(v1);
\vertex[right=1.5cm of v1] (v3); \vertex[above right=1.06cm of v1] (v2);
\vertex[below right=1.06cm of v1] (v4); \diagram* { (v1) -- [ half left,
looseness=1.7] (v3) -- [ half left, looseness=1.7] (v1) }; \diagram* { (v1) --
[ half right, looseness=0.5] (v2) -- [ half right, looseness=0.5] (v3) -- [
half right, looseness=0.5] (v4) -- [ half right, looseness=0.5] (v1) };
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{1.0pt}{0.0pt}\pgfsys@curveto{1.0pt}{0.55229pt}{0.55229pt}{1.0pt}{0.0pt}{1.0pt}\pgfsys@curveto{-0.55229pt}{1.0pt}{-1.0pt}{0.55229pt}{-1.0pt}{0.0pt}\pgfsys@curveto{-1.0pt}{-0.55229pt}{-0.55229pt}{-1.0pt}{0.0pt}{-1.0pt}\pgfsys@curveto{0.55229pt}{-1.0pt}{1.0pt}{-0.55229pt}{1.0pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{gathered}+\cdots$

Figure 1: Feynman diagrams corresponding to LO (left) and NLO (right) in an
1/N expansion. In contrast to standard perturbation theory, the diagrams are
built from self-consistently dressed propagators, resumming an infinite set of
perturbative diagrams.

For a system with large occupation numbers, classical dynamics may be a good
approximation to the full quantum evolution. In that case, variation of the
action provides classical equations of motion for the field variables. Solving
these equations straightforwardly in real time and averaging observables over
a suitable ensemble of initial realisations constitutes the classical
statistical approximation to quantum field theory (see, for instance
Tranberg:2022noe ; Berges:2013lsa ; Aarts:2001yn ). By construction, all
quantum effects are then discarded.444In some applications, one chooses to
initialise the classical simulations with ”quantum” vacuum initial
fluctuations (sometimes dubbed ”the half”). For a short enough time, the
classical state will remain quantum-like, and computing for instance effective
masses will yield result mimicking the quantum computation. This does however
not imply quantum dynamics, and for late times the classical dynamics will
relax the initial state to a classical equilibrium.

Simulating quantum dynamics in real time is highly nontrivial. One very
successful approach is to derive evolution equations555Kadanoff-Baym equations
or real-time Schwinger-Dyson equations for the propagator and the mean field
from the 2PI effective action Berges:2000ur . We define

$\displaystyle\bar{\phi}_{a}(t)=\langle\phi_{a}({\bf x},t)\rangle,\qquad
G_{ab}({\bf x-y},t,t^{\prime})=\langle\phi_{a}({\bf x},t)\phi_{b}({\bf
y},t)\rangle,$ (20)

where we have explicitly assumed a homogeneous state. The 2PI effective action
may then be written in the form Cornwall:1974vz :

$\Gamma[\bar{\phi},G]=S[\bar{\phi}]-\frac{i}{2}\rm Tr\ln{G}+\frac{i}{2}\rm
Tr\,G_{0}^{-1}[\bar{\phi}]G+\Phi[\bar{\phi},G],$ (21)

where for our model,

$\displaystyle
iG^{-1}_{0,\,ab}(x,y)=\frac{\delta^{2}\,S[\bar{\phi}]}{\delta\bar{\phi}_{a}(x)\,\bar{\phi}_{b}(y)}=\left(\partial^{2}_{x}\,\delta_{ab}+\mu^{2}\,\delta_{ab}-\frac{\lambda}{6N}\,\bar{\phi}_{c}\bar{\phi}_{c}\,\delta_{ab}-\frac{\lambda}{3N}\,\bar{\phi}_{a}\bar{\phi}_{b}\right)\delta(x-y).$
(22)

The functional $\Phi[\bar{\phi},G]$ may be written as an expansion in terms of
2PI skeleton diagrams with full propagators and mean fields and bare vertices.
The equations of motion are obtained from the stationarity conditions

$\displaystyle\frac{\delta\Gamma[\bar{\phi},G]}{\delta
G_{ab}(x,y)}=0,\quad\text{and}\quad\frac{\delta\Gamma[\bar{\phi},G]}{\delta\bar{\phi}_{a}(x)}=0.$
(23)

In what follows we consider a system initially in the symmetric phase
$\bar{\phi}=0$. Because of O(N) symmetry, if $\bar{\phi}=0$ at the initial
time it will remain zero for all time. The first equation is formally solved
by

$\displaystyle\frac{\delta\Gamma[\bar{\phi},G]}{\delta G_{ab}(x,y)}=0\to
G^{-1}_{ab}(x,y)=G_{0,ab}^{-1}(x,y)+i\Sigma_{ab}(x,y),$ (24)

in terms of the self-energy

$\displaystyle\Sigma_{ab}(x,y)=-2\frac{\delta\Phi[\bar{\phi},G]}{\delta
G_{ab}(x,y)}.$ (25)

Multiplying (24) by $G_{bc}(y,z)$ turns it into an integro-differential
equation,

$\displaystyle\delta_{ac}\delta(x,z)=\int d^{4}{\bf
y}\,\left[G^{-1}_{0,\,ab}(x,y)+i\Sigma_{ab}(x,y)\right]\,G_{bc}(y,z).$ (26)

At leading order (LO) in a $1/N$ expansion we have Aarts:2002dj

$\displaystyle\Phi^{\rm LO}[G]$ $\displaystyle=-\frac{\lambda}{4!N}\int
d^{4}x\,G_{aa}(x,x)G_{bb}(x,x),$ (27)

Since $G_{aa}\propto N$, the entire contribution is $\mathcal{O}(N)$ (see
left-most diagram in Figure 1). At NLO we have

$\displaystyle\Phi^{\rm NLO}[G]$ $\displaystyle=\frac{i}{2}\int
d^{4}x\,\ln\textrm{B}(G)(x,x),$ (28)

where

$\displaystyle\textrm{B}(x,y;G)=\delta_{\mathcal{C}}(x-y)+i\frac{\lambda}{6N}G_{ab}(x,y)G_{ab}(x,y).$
(29)

The expression (28) resums an infinite set of diagrams of $\mathcal{O}(N^{0})$
(see Figure 1, Aarts:2002dj ). We note that the figure-8 diagram appears at
both LO and NLO, providing local contributions to the self-energy at both
orders. We will work mostly at NLO, sometimes comparing to LO. Including
diagrams at NNLO is numerically highly challenging Aarts:2008wz , and we will
not do so here.

At this stage, it is useful to further decompose the propagator in terms of
the statistical ($F$) and spectral ($\rho$) components as

$\displaystyle
G_{ab}(x,y)=F_{ab}(x,y)-\frac{i}{2}\rho_{ab}(x,y)\textrm{sign}_{\mathcal{C}}(x^{0}-y^{0}),$
(30)

where

$\displaystyle
F_{ab}(x,y)=\frac{1}{2}\langle\\{\phi_{a}(x),\phi_{b}(y)\\}\rangle,\qquad\rho_{ab}(x,y)=i\langle[\phi_{a}(x),\phi_{b}(y)]\rangle,$
(31)

and the sign-function refers to ordering along the Keldysh contour in the
complex time plane666This decomposition is different from, but completely
equivalent to, the +/- formalism, often stated to be doubling the number of
degrees of freedom.. It follows that $F$($\rho$) is symmetric (antisymmetric)
with respect to space-time indices $(x,y)$. In particular,
$G_{ab}(x=y)=F_{ab}(x=y)$. An $O(N)$ symmetric state ($\bar{\phi_{a}}=0$) and
spatial homogeneity allows to further simplify $F_{ab}(x,y)=\delta_{ab}F({\bf
x-y},t,t^{\prime})$ and $\rho_{ab}(x,y)=\delta_{ab}\rho({\bf
x-y},t,t^{\prime})$. We may similarly decompose the self-energy into two
components

$\displaystyle\Sigma_{ab}(x,y)=\Sigma_{F,ab}(x,y)-\frac{i}{2}\Sigma_{\rho,ab}(x,y)\textrm{sign}_{\mathcal{C}}(x^{0}-y^{0}),$
(32)

Inserting this decomposition into eq. 26, one obtains the real-time equations
of motion for $F$ and $\rho$:

$\displaystyle\Big{(}\partial_{t}^{2}-\partial_{x}^{2}+M^{2}(x)\Big{)}F(x,y)=$
$\displaystyle-\int_{0}^{x^{0}}dz^{0}\int d^{3}z\Sigma^{\rho}(x,z)F(z,y)$
$\displaystyle+\int_{0}^{y^{0}}dz^{0}\int d^{3}z\Sigma^{F}(x,z)\rho(z,y),$
(33)
$\displaystyle\Big{(}\partial_{t}^{2}-\partial_{x}^{2}+M^{2}(x)\Big{)}\rho(x,y)=$
$\displaystyle-\int_{y^{0}}^{x^{0}}dz^{0}\int
d^{3}z\Sigma^{\rho}(x,z)\rho(z,y).$ (34)

We note that the evolution equations are explicit, and depend on the entire
past history of the evolution through the time integrals over the self-energy
and propagator components. The local parts of the self-energies are accounted
for in the mass term

$\displaystyle M_{\rm LO}^{2}(t)=m^{2}+\frac{\lambda}{6}F({\bf 0},t,t),\quad
M_{\rm NLO}^{2}(t)=M_{\rm LO}^{2}(t)+\frac{\lambda}{3N}F({\bf 0},t,t).$ (35)

At NLO the non-local self energies are given by

$\displaystyle\Sigma^{\rho}(x,y)$
$\displaystyle=-\frac{\lambda}{3N}\Big{(}F(x,y)I_{\rho}(x,y)+\rho(x,y)I_{F}(x,y)\Big{)},$
(36) $\displaystyle\Sigma^{F}(x,y)$
$\displaystyle=-\frac{\lambda}{3N}\Big{(}F(x,y)I_{F}(x,y)-\frac{1}{4}\rho(x,y)I_{\rho}(x,y)\Big{)},$
(37)

where the $1/N$ resummation is performed through the objects $I_{F,\rho}$,
that in turn satisfy777We note in passing that not iterating the $I_{F,\rho}$,
but simply inserting the local (no time integration) part of the expression
into (36, 37) reduces the evolution to NLO in a coupling expansion
Berges:2000ur ; Arrizabalaga:2005tf .

$\displaystyle I_{\rho}(x,y)=$
$\displaystyle-\frac{\lambda}{3}\int_{y^{0}}^{x^{0}}d^{3}zI_{\rho}(x,z)F(z,y)\rho(z,y)+\frac{\lambda}{3}F(x,y)\rho(x,y),$
(38) $\displaystyle I_{F}(x,y)=$
$\displaystyle-\frac{\lambda}{6}\int_{0}^{x^{0}}d^{3}zI_{\rho}(x,z)\big{(}F^{2}(z,y)-\frac{1}{4}\rho^{2}(z,y)\big{)}+\frac{\lambda}{3}\int_{0}^{y^{0}}d^{3}zI_{F}(x,z)F(z,y)\rho(z,y)$
$\displaystyle+\frac{\lambda}{6}\Big{(}F^{2}(x,y)-\frac{1}{4}\rho^{2}(x,y)\Big{)}.$
(39)

###  3.2 Numerical implementation

The non-linear integro-differential equations of motion may be subjected to
numerical evaluation, in terms of a discrete set of degrees of freedom. This
may be done either through discretising the system on a space-time lattice at
the level of the action, or through discretising momentum space and time at
the level of the equations of motion. We will here opt for the lattice
implementation.

####  3.2.1 Discretization and observables

We introduce a cubic lattice of $N_{x}^{3}=32^{3}$ sites with lattice spacing
$a\mu=0.7$, and a time discretization $a_{t}=adt$, $dt=0.1$. Although in
principle, the memory integrals over the non-local self-energies should
stretch all the way to the initial time, in practice we keep only the last
$n_{t}$ timesteps for a memory range $t_{K}=n_{t}a_{t}$. Further discussion on
this point may be found in Appendix A. The lattice momentum operator is

$\displaystyle k_{\textrm{lat}}^{2}=\sum_{i=1}^{3}\left(2-2\cos
k_{i}\right),\quad k_{i}=\frac{2\pi}{N_{x}}n_{i},\quad
n_{i}\in\Big{(}-\frac{N_{x}}{2}+1,\frac{N_{x}}{2}\Big{)}.$ (40)

The initial conditions for the spectral propagator follow from the basic
commutation relations

$\displaystyle\rho_{k}(t,t)=0,\quad\partial_{t}\rho_{k}(t,t^{\prime})|_{t=t^{\prime}}=1,$
(41)

while for the statistical propagator we choose as the initial state the free-
field quantum vacuum prior to the quench

$\displaystyle F_{k}(t,t^{\prime})|_{t=t^{\prime}=0}$
$\displaystyle=\frac{1}{\omega_{k}}\Big{(}n_{k}+\frac{1}{2}\Big{)},$ (42)
$\displaystyle\partial_{t}F_{k}(t,t^{\prime})|_{t=t^{\prime}=0}$
$\displaystyle=0,$ (43)
$\displaystyle\partial_{t}\partial_{t^{\prime}}F_{k}(t,t^{\prime})|_{t=t^{\prime}=0}$
$\displaystyle=\omega_{k}\Big{(}n_{k}+\frac{1}{2}\Big{)},$ (44)

where $n_{k}=0$ and $\omega_{k}=k_{\textrm{lat}}^{2}+m^{2}(t<0)$. Conversely,
for all later times, we may extract a particle number and a dispersion
relation from

$\displaystyle\omega_{k}(t)$
$\displaystyle=\sqrt{\frac{\partial_{t}\partial_{t^{\prime}}F_{k}(t,t^{\prime})|_{t^{\prime}=t}}{F_{k}(t,t)}},$
(45) $\displaystyle n_{k}(t)+\frac{1}{2}$
$\displaystyle=c_{k}(t)\sqrt{\partial_{t}\partial_{t^{\prime}}F_{k}(t,t^{\prime})|_{t^{\prime}=t}F_{k}(t,t)},$
(46)

where $c_{k}(t)=\sqrt{1-\frac{1}{4}dt^{2}\omega_{k}^{2}(t)}$ is a lattice
correction. Very far from equilibrium, these quantities make little sense, but
once the system settles down enough that the field excitations are quasi-
particle like, they are good representations of the occupation and frequency
of a given mode. Real time simulations derived from the 2PI effective action
have the advantage of conserving an energy functional. We monitor the energy
density and pressure density during simulations, from which we get the
equation of state parameter $\omega=\frac{P}{\rho}$. The expressions are given
by

$\displaystyle\frac{\rho(t)}{N}=$ $\displaystyle\frac{1}{2}\int
d^{3}k\partial_{t}\partial_{t^{\prime}}F_{k}(t,t)+\frac{1}{2}\int
d^{3}k\Big{(}k^{2}+m^{2}+\frac{\lambda}{2}\frac{N+2}{6N}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}F_{k^{\prime}}(t,t)\Big{)}F_{k}(t,t)$
(47) $\displaystyle+\frac{1}{4}\int_{0}^{t}dt^{\prime}\int
d^{3}k\Big{(}\Sigma^{\rho}_{k}(t,t^{\prime})F_{k}(t^{\prime},t)-\Sigma^{F}_{k}(t,t^{\prime})\rho_{k}(t^{\prime},t)\Big{)}+V_{0},$
(48)

and

$\displaystyle\frac{P(t)}{N}=$ $\displaystyle\frac{1}{2}\int
d^{3}k\partial_{t}\partial_{t^{\prime}}F_{k}(t,t)-\frac{1}{2}\int
d^{3}k\Big{(}k^{2}+m^{2}+\frac{\lambda}{2}\frac{N+2}{6N}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}F_{k^{\prime}}(t,t)\Big{)}F_{k}(t,t)$
(49) $\displaystyle-\frac{1}{4}\int_{0}^{t}dt^{\prime}\int
d^{3}k\Big{(}\Sigma^{\rho}_{k}(t,t^{\prime})F_{k}(t^{\prime},t)-\Sigma^{F}_{k}(t,t^{\prime})\rho_{k}(t^{\prime},t)\Big{)}+\frac{1}{3}\int
d^{3}k\;k^{2}F_{k}(t,t)-V_{0}.$ (50)

Initially, the system is in a free-field vacuum state with energy
$\rho_{0}/N=V_{0}$ and pressure $P_{0}/N=-V_{0}$, where
$V_{0}=\frac{3\mu^{4}}{2\lambda}$. Hence the initial equation of state
parameter is $\omega=-1$. Energy density and pressure are divergent quantities
and subject to renormalisation as described in the following.

####  3.2.2 Renormalisation

The evolution equations are UV divergent, and although the lattice provides a
regularization, a substantial cut-off dependence is introduced in particular
through the local self-energies, calling for renormalisation. Hence given an
input renormalised $a\mu=0.7$, we must perform the simulation with a bare mass
parameter given by

$\displaystyle-(a\mu_{0})^{2}=-(a\mu)^{2}-\delta(a\mu)^{2},$ (51)

where the counterterm on the lattice is given by:

$\displaystyle\delta(a\mu)^{2}=\lambda\frac{N+2}{6N}\frac{1}{N_{x}^{3}}\sum_{k_{\rm
lat}}\frac{1}{2\sqrt{(ak_{\rm
lat})^{2}+(a\mu)^{2}}}\simeq\Big{|}_{a\mu=0.7}\lambda\frac{N+2}{6N}0.212.$
(52)

Note that at LO the prefactor changes $\frac{N+2}{6N}\rightarrow\frac{1}{6}$
to account for only the LO local self-energy eq. 35. At higher loop order more
divergencies appear that may be renormalized by appropriately choosing the
coupling (see for instance Berges:2004hn ; Berges:2005hc ). As discussed in
Arrizabalaga:2004iw  the effect is less than $10\%$ for the couplings and
cutoff under consideration, and since we are not intending to take the
continuum limit, we neglect this complication.

The energy density and pressure density (47), (49) are also divergent
quantities. After implementing the mass renormalisation, we can remove the
remaining divergences by adding an overall constant counterterm (one for
energy, one for pressure). As renormalisation condition, we choose for the
initial renormalised energy (pressure) density to be $+V_{0}$ and $-V_{0}$,
respectively.

####  3.2.3 The classical approximation

The classical limit is intuitively when the occupation numbers are large,
$n_{k}\gg 1$, and indeed one may extend classical perturbation
theory888Omitting certain vertices in the Keldysh field basis at the level of
the action. to this resummed non-equilibrium framework Aarts:1997kp ;
Aarts:2001yn , in which case the classical limit amounts to

$\displaystyle F^{2}\gg\rho^{2}.$ (53)

Since in the quasi-particle limit, $F\propto n_{k}$ and $\rho\propto 1$, this
requirement is more general that $n_{k}\gg 1$, but the two are consistent when
both apply.

The 2PI-classical approximation consequently follows from discarding from the
evolution equations (36, 37, 38, 3.1) all instances of $\rho^{2}$ (and $\rho
I_{\rho}$) when added to $F^{2}$ (or $FI_{F}$) For the system and truncation
considered here, we find

$\displaystyle\Sigma_{\rm cl}^{F}(x,y)$
$\displaystyle=-\frac{\lambda}{3N}F(x,y)I_{F}(x,y),$ (54) $\displaystyle
I^{\rm cl}_{F}(x,y)$ $\displaystyle=\frac{\lambda}{6}F^{2}(x,y)$
$\displaystyle-\frac{\lambda}{6}\int_{0}^{x^{0}}d^{3}zI_{\rho}(x,z)F^{2}(z,y)+\frac{\lambda}{3}\int_{0}^{y^{0}}d^{3}zI_{F}(x,z)F(z,y)\rho(z,y).$
(55)

we note that at LO, the classical and the quantum evolution are the
same999This is true for Gaussian approximations, although beware of subtleties
such as Millington:2020vkg ..

##  4 Early time dynamics

###  4.1 Linear regime

Figure 2: Left: Time evolution of the statistical propagator $F(t,t,x=0)$ in
the free theory, LO and NLO approximation for $N=4$ and $\lambda=1$. The
dashed (dotted) horizontal line corresponds to the minimum ($\alpha=1$) and
inflection point ($\alpha=1/3$) of the potential, respectively. Right:
Comparing the free field to the LO evolution for different $\lambda$. Note
that the result is independent of the number of fields $N$. Also note the
different normalisation of $F$ in the left vs. righthand plot.

The initial state for the transition is the quantum vacuum prior to the mass
quench. Since at this stage, occupation numbers are clearly not large, the
classical approximation does not apply. However, the quantum and classical
evolution of the propagators coincide for free fields and in the LO
approximation. At small enough coupling, it is therefore possible that
classical dynamics reproduce the quantum dynamics well in this regime, after
which the classical approximation becomes robust by virtue of large occupation
numbers.

When the field self-interaction may be neglected, $\lambda=0$, the evolution
equation becomes linear

$\displaystyle\Big{(}\partial_{t}^{2}+\omega_{k}^{2}(t)\Big{)}F_{k}(t,t^{\prime})$
$\displaystyle=0,\quad\text{with}\quad\omega_{k}^{2}(t)=m^{2}(t)+k^{2},$ (56)

and closed-form solutions exist for the instantaneous (eqs. 18, 19,
Smit:2002yg ) and linear quench Garcia-Bellido:2002fsq .

Since $\lambda=0$, there is no notion yet of $v$ (the bottom of the potential)
and also no separation into Higgs and Goldstone modes. The free field
approximation to the interacting dynamics is reliable until the non-linear
term (order $\lambda$, at LO or NLO) becomes comparable to the mass parameter,

$F(x,x)=\alpha\frac{6\mu^{2}}{\lambda}.$ (57)

for some number $\alpha$. Often, the criterion is taken to be either
$\alpha=1$ (the minimum of potential is reached) or $\alpha=1/3$, (the
correlator $F$ passes the inflection point of the potential). Figure 2 (left)
shows the equal-time statistical propagator in the free field approximation,
and with LO and NLO corrections. The dashed line is the minimum of the
potential and the dotted line the inflection point. We see that the latter is
perhaps a better criterion for the onset of non-linearities, as the LO (and
NLO) depart from the pure exponential growth around $\mu t=4$. The non-linear
evolution does still overshoot the potential minimum before settling down to
an equilibrium state. We will discuss the differences between LO and NLO
below.

Figure 3: Left: Time at which free theory approximation deviates from the LO
result by 20 percent ($\Delta F(t,t,x=0)=0.2$ ). Right: Particle spectrum in
the LO approximation for $N=4$ at the time when the free field approximation
becomes unreliable ($\Delta F(t,t,x=0)=0.2$) for different values of
$\lambda$. Even relatively high values of $\lambda$ produce highly occupied
modes at this very early stage. The spectra are calculated with a lattice of
size $N_{x}^{3}=128^{3}$.

Furthermore, we can analyse the influence of the interaction strength
$\lambda$ on the free field approximation. In figure 2 (right) we see that the
inflection point of the potential (equivalently the minimum of the potential)
is reached earlier with higher $\lambda$ in the LO approximation, thus
limiting the applicability for the free field approximation to earlier times.
We can make this statement more precise by defining a time at which the equal
time statistical propagator in the free field approximation differs from the
LO approximation by 20 percent:

$\Delta F(t)=\Big{|}\frac{F_{free}-F_{\lambda}}{F_{free}}\Big{|}(t)=0.2.$ (58)

Figure 3 (left) shows the relation between this time and $\lambda$. Based on a
simple fit, we conclude that free theory is a good approximation until a time

$\mu t^{*}=a-b\ln\lambda\quad\text{where}\quad a=4.13\pm 0.20,\quad b=1.00\pm
0.03.$ (59)

The standard error of the mean is obtained by choosing 5 values of $\Delta
F\in(0.1,0.3)$. We can examine the particle spectrum at the time when the free
theory approximation breaks down. Figure 3 (right) shows the occupation
numbers in the IR for various values of $\lambda$. Although the free field
regime only lasts for a time $\mu t\simeq 2-4$ we see that the unstable modes
are already highly populated $n_{k}\gg 1$, even for large values of $\lambda$.
At least for these modes, classical evolution should be robust.

###  4.2 The classical regime

As we have seen, in the language of correlators, classicality follows in the
limit where $F^{2}\gg\rho^{2}$ (see also Aarts:1997kp  for a discussion of the
classical-statistical approximation). Often, $n_{k}>1$ is applied as a minimal
requirement, and often it is stated that as long as this is fulfilled for the
modes directly relevant to the phenomenon of interest, it is sensible to
proceed classically even though classical modes interact with non-classical
modes. We have seen that the requirement is fulfilled during a spinodal
transition, at least in the range $k_{L}/\mu<0.8$. We further test the
classical nature of the tachyonic transition by computing the evolution of
$F(t,t,x=0)$ in the LO, the quantum-NLO and the classical-NLO approximation,
shown in figure 4 (left). We see that the classical limit and the quantum
evolution agree very well, even though only unstable modes are highly
occupied. Figure 4 (right) shows the same agreement in terms of the total
particle number.

Figure 4: Left: Time evolution of the statistical correlator $F(t,t,x=0)$ at
LO, quantum-NLO and classical-NLO. Apparently, occupation numbers are large
enough for the system to act effectively classically. Right: The time
evolution of the total particle number, summed over all momentum modes. $N=4$,
$\lambda=6$.

Since self-interactions stop the instability when eq. 57 is satisfied, we can
expect the individual occupation numbers to scale as $n_{k}\propto 1/\lambda$.
During the evolution, the total particle number (summed over all the modes)
grows to a maximum value and then decreases. In figure 5 we show this maximal
value as a function of $\lambda$, at LO and NLO. We see that the LO dynamics
describes the early time particle production process very well, and we find
that for all $\lambda$, the maximal particle number occurs at times $\mu
t<15$.

Figure 5: Highest value for total particle numbers $n_{tot}$ for the LO and
NLO truncation. Even for large coupling, $n_{\rm tot}$ reaches large values.
$N=4$.

##  5 Intermediate time dynamics

###  5.1 The classical momentum range

While the early-time dynamics of the tachyonic transition is well described by
first a free field and then the LO approximation, at times $\mu t>6$, the non-
local dynamics first entering at NLO becomes important.

Figure 6: Left: The particle spectrum at the early stage of NLO evolution.
$N=4$, $\lambda=6$. Right: The largest momentum wih $n_{k}>1$ in time, for
different values of $\lambda$. $N=4$. Figure 7: Time evolution of $F(t,t,x=0)$
at NLO, in the quantum system and the classical approximation.

Figure 6 (left) shows the particle spectrum at several times during the
evolution. The dynamics at NLO allows for exchange of energy and particles
between the modes. We see that for $\mu t\lesssim 40$ the classical range of
momenta with $n_{k}>1$ grows, as the high occupancy in the IR is redistributed
towards the UV. This redistribution continues, and we show (right) the value
of the largest momentum with $n_{k}>1$ in time. After the early growth of the
classical range, at intermediate times this range shrinks again. Judging from
the local correlator $F(t,t)$ only (figure 7), the classical approximation for
the whole system is valid for time-scales $\mu t\simeq\mathcal{O}(500)$.

Several out-of-equilibrium processes may have taken place during the initial
violent stage of tachyonic preheating. Exotic heavy particles, topological
defects, sphalerons and other relics such as primordial black holes may have
been created in this environment of large IR field fluctuations. We can
attempt to describe this out-of-equilibrium stage through an effective
infrared temperature $\overline{T}_{\rm IR}$, which we define to be the
average temperature in the classical range, had the spectrum been classical
equilibrium

$\overline{T}_{\rm IR}=\frac{\sum_{n_{k}>1}\omega_{k}n_{k}}{\sum_{n_{k}>1}1}.$
(60)

This is shown in figure 8 as a function of time. We see that although the
asymptotic equilibrium temperature is $T_{\rm reh}\simeq\mu$ or less, at
intermediate times the effective temperature may be twice as large or more.
Provided the process of interest takes place over a time-scale $\mu
t\simeq\mathcal{O}(100-1000)$, this is likely the relevant temperature to
enter estimates of particle production and sphaleron rates.

Figure 8: Time evolution of the effective temperature in the IR according to
eq. 60. The inset shows the values of the temperature at the end of the
evolution. The discontinuities are caused by modes falling out of the
classicality region and are therefore a discretisation effect. $N=4$.

###  5.2 Kinetic equilibration

Figure 9: Left: Particle spectra for various stages of equilibration and the
corresponding fits to Bose-Einstein distributions (straight lines). Right: The
$R^{2}$ value of the fit quantifies how close the particle spectrum is to a
thermal distribution. $N=4$, $\lambda=6$.

The (quantum-)NLO evolution equations provide correct thermalization, in the
sense that at late times, the occupation numbers approach a quantum thermal
equilibrium state. There are two thermalization stages: kinetic equilibration,
whereby particles are redistributed into a thermal distribution with a non-
zero chemical potential; followed by chemical equilibration, where the overall
number of particles changes to eventually relax the system to zero chemical
potential and its final temperature. These processes may be identified with
certain constituent perturbative diagrams (2-to-2 scattering, 2-to-4
scattering, …), and so for small coupling, the timescales are expected to be
well separated. We can quantifiy the equilibration process by performing a fit
of the particle spectrum with a Bose-Einstein distribution:101010We trust that
the standard notation using $\mu$ as the chemical potential is not confused
with the mass parameter $\mu^{2}$.

$n_{k}=\frac{1}{e^{(\omega_{k}-\mu)/T}-1},$ (61)

where $\omega_{k}$ is the derived dispersion relation (45).

Figure 10: Time scale of kinetic equilibration: Time when $R^{2}$ of Bose-
Einstein distribution fit reaches $0.97$ (left) and time scale defined by eq.
64.

We perform a least square fit with

$\ln\big{(}1+1/n_{k}\big{)}=(\omega_{k}-\mu)/T,$ (62)

and quantify kinetic equilibration by calculating the $R^{2}$ value of this
fit:

$R^{2}=1-\frac{\sum\big{(}y_{i}-f_{i}\big{)}^{2}}{\sum\big{(}y_{i}-\bar{y}\big{)}^{2}},$
(63)

where $y_{i}$ refers to the data points, $\bar{y}$ to the mean and $f_{i}$ the
predictions of the model with best fit parameters. A $R^{2}$ value of 1
corresponds to a perfect fit where all data points are matched by the
model.111111The best fit parameters $T$ and $\mu$ in 62 can be associated with
an uncertainty from the fit procedure, quantifying a range of parameters in
agreement with data. In the following these uncertainties are propagated to
observables as $T_{\rm reh}$ and time scales. We stress that they are
systematic rather than statistical.

In figure 9 we show the process of kinetic equilibration, for $N=4$,
$\lambda=6$. On the left we show the particle spectrum at various stages of
equilibration and the corresponding fit to a Bose-Einstein distribution. We
focus on the IR defined by $\omega_{k}/\mu<1.5$, indicated by the dashed
vertical line. We perform the fit in this region and show that the UV modes
take longer to organize themselves into a thermal like distribution. In figure
9 (right) we show how the $R^{2}$ gradually approaches 1. We see that after
the initial growth of particle numbers at early times the spectrum is far from
thermal, but at intermediate times $\mu t\approx 150$ the particles are
reorganised into a thermal-like distribution. At this stage the distribution
is well approximated by eq. 61. We define kinetic equilibration to be
completed when $R^{2}$ reaches $0.97$, which happens at $\mu t=150$ in figure
9. Furthermore we can extract a kinetic equilibration time scale by fitting
the $R^{2}$ evolution to the form

$R^{2}(t)=R^{2}_{\rm final}+(R^{2}_{\rm init}-R^{2}_{\rm
final})e^{-t/\tau_{\rm kin}}.$ (64)

where we discard $R^{2}$ values $<0.97$. The fit is also shown in figure 9
(right).

Figure 10 shows the time when $R^{2}$ reaches $0.97$ (left) and the kinetic
equilibration time defined by eq. 64. The two are consistent, and reveal a
linear dependence on $N$ and that kinetic equilibration is faster for larger
coupling. At the onset of kinetic equilibration we extract a temperature,
shown in figure 11, with a chemical potential in the range of $0.2-0.33$. In
the following stage of chemical equilibration the temperature approaches its
final value and the chemical potential goes to zero.

Figure 11: Temperature at the time when kinetic equilibrium is reached.

Figure 12: Left: The particle spectrum at $\mu t=550$ for $N=4$, $\lambda=6$
at NLO and a fit to a Bose-Einstein distribution, which in this variable
corresponds to a straight line $(\omega_{k}-\mu)/T$. The red dots show the
spectrum in the classical approximation. Right: Comparison of UV part of the
particle spectrum at $\mu t=550$ for quantum evolution at NLO and in the
classical approximation. $N=4$, $\lambda=6$.

In figure 12 (left), we show an advanced stage of kinetic equilibration, for
$N=4$, $\lambda=6$. The distribution is very well approximated by eq. 61
($R^{2}>0.99$). For comparison, we show the equivalent spectrum in a
classical-NLO simulation, which does not agree very well, even for modes with
$n_{k}>1$ (shown by the dashed horizontal line). The classical spectrum is
lower than the quantum spectrum, indicating that the occupation numbers are
larger. This demonstrates that one should be wary of approximating
intermediate to late time quantum dynamics by classical dynamics. In figure 12
(right), we show the UV part of the spectrum121212But without performing the
transformation where a Bose-Einstein distribution becomes a straight line. at
this same time of $\mu t=550$. The quantum dynamics correctly conserves the
zero-point fluctuations ($n_{k}+1/2\geq 1/2$, $n_{k}\geq 0$), as shown by the
green dots. Zero point fluctuations are there to scatter on, but it is not
physical to extract net particles or energy from the vacuum. The classical
simulation however has no such constraint, and does not distinguish between
the $n_{k}$ and the $1/2$, as shown by the red dots. As such, it allows to
deplete the UV vacuum so that $n_{k}+1/2<1/2$, to incorrectly increase the
occupation in the IR and ultimately to thermalise to a classical equilibrium.
In a manifestly classical simulation, one might consider not initialising with
the zero-point fluctuations, or only in the unstable modes Arrizabalaga:2004iw
. This has some advantages, although one would lose some of the physics of
scattering on the UV fluctuations, and the late time equilibrium would still
be off. In any case, it is clear that replacing quantum dynamics by classical
evolution with all zero-point fluctuations introduces a cut-off dependent
energy reservoir, which the classical dynamics may (unphysically) tap into.

###  5.3 Chemical equilibration

Figure 13: Evolution of the Bose-Einstein fit parameters $T$ and $\mu$ in
time, and the fit to eq. 65 and eq. 66. The dashed line indicates the
beginning of the fit region. $N=4$, $\lambda=6$.

The fit parameters $\mu$ and $T$ are time-dependent, and are expected to
asymptote to zero and some finite reheating temperature $T_{\rm reh}$,
respectively. Figure 13 shows $T$ and $\mu$ in time. We may further extract a
final temperature $T_{\rm reh}$ and a chemical equilibration timescale
$\tau_{\rm ch}$ by a fit to

$T(t)=T_{\rm reh}+(T_{\rm init}-T_{\rm reh})e^{-t/\tau_{\rm ch}}.$ (65)

and

$\mu(t)=\mu_{\rm final}+(\mu_{\rm init}-\mu_{\rm final})e^{-t/\tau_{\rm ch}}.$
(66)

The fits give consistent equilibration time scales. We must note here that we
are somewhat limited by the available computational resources. The asymptotic
temperature $T_{\rm reh}$ and chemical equilibration time scale $\mu\tau_{\rm
ch}$ depend on the available time-range. Based on simulations of a smaller
lattice size $16^{3}$ to much late times ($\mu t=2800$), we find that data
until $\mu t=700$ overestimates $T_{\rm reh}$ by about $20\%$. Furthermore,
the chemical equilibration time is underestimated by about $80\%$ (see
Appendix B).

Figure 14: Left: Extrapolated equilibrium temperatures for various values of
$N$ and $\lambda$. Shown are values obtained from the simulation (dots) and a
simple Higgs-Goldstone model (solid line) and for Goldstone modes only (dotted
line). Right: Chemical equilibration time for various values of $N$ and
$\lambda$.

Figure 15: Left: Coupling strength dependence of $T_{\rm reh}$ and power fit
$\propto\lambda^{-x}$ for which $x=0.37$. Right: Coupling strength dependence
of chemical equilibration time for $N=4$.

With this caveat in mind, we show in figure 14 the reheating temperature
$T_{\rm reh}$ for various $N$ and $\lambda$ values. We also show for
comparison temperatures from a simple, non-interacting model in which all
available energy is distributed among a massive Higgs mode with
$m_{H}=\sqrt{2}\mu$ and $N-1$ massless Goldstone modes. We obtain the
temperature in this model by solving

$\displaystyle\frac{3N\mu^{4}}{2\lambda}=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\Big{(}\sqrt{k^{2}-m_{H}^{2}}\,n_{k}(m_{H},T)+(N-1)\,k\,n_{k}(0,T)\Big{)},$
(67)

for $T$, where we assume equilibrated constituents
$n_{k}(m,T)=(e^{\sqrt{k^{2}-m^{2}}/T}-1)^{-1}$. The temperatures are almost
independent of $N$. An even simpler model is to discard the Higgs mode, in
which case one returns to (5) (with $N-1$ degrees of freedom). In that case
the temperature scales with $T\propto\Big{(}\frac{N}{N-1}\Big{)}^{1/4}$, with
$T/\mu=0.93$ for $\lambda=6$. Taking into account the overestimation of
$T_{\rm reh}$ due to the time range limitation, we see that the fit
temperatures end up between the Higgs-Goldstone and Goldstone-only models.
Figure 14 (right) shows the kinetic equilibration time $\mu\tau_{\rm kin}$ and
we identify an approximate linear $N$ dependence. Furthermore, we find as for
kinetic equilibration that larger couplings lead to faster equilibration. The
systematic errors estimated from the fitting procedure are however
substantial. In addition, we note that for these rather large couplings, the
chemical equilibration times are only about twice as large as the kinetic
equilibration time scale, so that it is hard to disentangle the two processes.
In figure 15 we show the coupling strength dependence on $T_{\rm reh}$ and
$\tau_{\rm ch}$.

###  5.4 The compound propagator

Figure 16: Left: Dispersion relation for $N=4$ and $\lambda=6$ and the fit to
the result of the model of the compound propagator. Right: The field
expectation value (squared) extracted from the propagator by eq. 10.

As we discussed in section 2, the O(N) symmetric correlator is a compound of a
Higgs mode and $N-1$ Goldstone modes. Sofar, we have been fitting the particle
number with a single Bose-Einstein distribution, in effect assuming that the
$N-1$ light modes dominate. In figure 16 (left), we fit the decomposition eq.
9 to the dispersion relation, while inserting the temperatures as obtained in
the previous chapter. The field expectation value can be recovered from the
zero mode of the propagator eq. 10, which we show in fig. 16 (right). We see
that although it asymptotes to a finite value, this is not the zero-
temperature (vacuum) expectation value, but is the minimum of the finite-
temperature effective potential. Also, settling into this minimum happens
gradually on a time-scale similar to kinetic equilibration.

###  5.5 The equation of state

Figure 17: Time evolution of the equation of state for $\lambda=6$. Left:
Equation of state derived from the energy density and pressure density
functionals. Right: Equation of state from a quasi particle ansatz.

Figure 18: As figure 17 but with $\lambda=1$.

For cosmological purposes, complete equilibration is often less important than
the effective equation of state of the system. In figure 17 we show this
quantity in time for $\lambda=6$. One might have expected that the state was
an equipartitioned mixture of $g^{*}$ ”radiation” Goldstone components with
$P/\rho\simeq 1/3$ and $(N-g^{*})$ ”matter” components with $P/\rho\simeq 0$.
In that case, the equation of state should be

$\displaystyle\omega\simeq\frac{g^{*}}{3N}.$ (68)

At low temperature, one would expect $g^{*}=N-1$, so that for $N=4$,
$\omega=1/4$. But the situation is more complicated.

On the left of figure 17 we show how the equation of state evolves from the
initial $\omega=-1$ to its final value for $\lambda=6$. Already after a time
$\mu t=100-200$, the final equation of state is established. The dependence on
$N$ and the overall magnitude is not as naively expected. This is because the
renormalisation procedure calibrates to a potential minimum of $V_{0}$, while
the finite temperature effective potential has much shallower minima at
smaller values of the field expectation value/zero mode. If we instead compute
the equation of state using a quasi-particle ansatz

$\displaystyle\rho=\int d^{3}k\,n_{k}\omega_{k},\qquad P=\frac{1}{3}\int
d^{3}k\;k^{2}\frac{n_{k}}{\omega_{k}},$ (69)

we find in figure 17 (left) that the equation of state indeed approaches its
radiation value of $1/3$ for larger values of Goldstone boson components. In
figure 18 we repeat our simulation, but for $\lambda=1$, where the expectation
value is closer to the zero temperature vev. We see that both for the direct
(left) and quasi-particle (right) equation of state, the system displays the
correct behaviour.

##  6 Finite quench time

Figure 19: Energy loss for finite time quenches. Note that the post-quench,
gradual energy loss is due to a truncation of the memory kernel (see appendix
A). $N=4$, $\lambda=6$.

A natural generalisation of the instantaneous quench it to consider finite-
time linear quenches of the form (13). Implementing this as a by-hand mass
function has the drawback that total energy is no longer conserved. We
straightforwardly have that

$\displaystyle\dot{\rho}=\frac{1}{2}\dot{\mu^{2}(t)}\phi_{a}^{2}=\dot{\mu}\mu\phi_{a}^{2},$
(70)

so that the total energy loss during the transition is

$\displaystyle\Delta\rho=\int_{0}^{t}dt^{\prime}\left[\dot{\mu}(t^{\prime})\mu(t^{\prime})\phi_{a}^{2}(t^{\prime})\right].$
(71)

For a very fast quench, $\phi_{a}$ is not able to react to the quench of the
mass parameter, and since we start out at $\phi_{a}=0$, the energy loss is
negligible. Conversely, if the flip is very slow (adiabatic), the field is
able to track the bottom of the potential
$\phi_{a}^{2}=\frac{6N\mu^{2}}{\lambda}$ in which case

$\displaystyle\Delta\rho=-\frac{6N}{\lambda}\int_{0}^{t}dt^{\prime}\dot{\mu}\mu^{3}=-\frac{3\mu^{4}(t)}{2\lambda}.$
(72)

In this limit, all the potential energy is lost, as the field is slowly
deposited in the new potential minimum. In figure 19 we show the energy loss
for different quench times. Of course, in a real dynamical realisation of the
$\sigma-\phi$ coupling, the energy is strictly conserved and instead
transferred to the inflaton $\sigma$. The whole system of $N+1$ degrees of
freedom eventually thermalise to some reheating temperature. We postpone this
highly model-dependent complication to future work, and restrict ourselves to
the early time dynamics

Figure 20: Left: Highest value for total particle numbers $n_{tot}$ for
various quench times. Right: Particle spectrum in the LO approximation ($N=4$,
$\lambda=1$) at the time when the linear approximation breaks down ($\Delta
F(t,t,x=0)=0.2$). Shown are spectra of unstable modes ($k_{L}/\mu<1$) for
different values of quench times $\mu\tau_{Q}$. The spectra are calculated
with a lattice of size $N_{x}^{3}=128^{3}$.

Figure 20 (left) shows the relation between quench time and the maximum total
particle number that is achieved during evolution. Unsurprisingly, a fast
quench leads to more overall particle production. As for the instantaneous
quench we show the particle spectrum at the time when the free field
description starts to fail. This is shown in figure 20 (right). We find that
for slower quenches a smaller range of modes become classical ($n_{k}>1$) at
this time, and that for quenches slower than $\mu\tau_{Q}=100$, it is
essentially only the zero mode that grows large during the transition. The
spectrum is also not completely monotonic in the quench time, although the
total number of particles created is.

Figure 21: Time evolution of the largest momenta that has occupation number
larger than one for different values of quench time. $N=4$, $\lambda=6$.

Finally, in figure 21, we show the evolution of the classical range for
different quench times. We see a separation into ”fast” quenches
$\mu\tau_{Q}<20$, where the classical range is large at first before it
adjusts, and ”slow” quenches $\mu\tau_{Q}>50$ where the classical range is
smaller, established later and changes on a longer timescale.

##  7 Conclusion

Tachyonic preheating provides a mechanism to reheat the universe after
inflation. During the process, a transient out-of-equilibrium state with very
high occupation numbers in the IR is generated, which subsequently relaxes to
an equilibrium reheating temperature $T_{\rm reh}$. Although the large
occupation numbers allow an effective description in terms of classical field
dynamics, quantum evolution equations are necessary to reach the correct
equilibrium state and on the correct time scale131313When not having access to
full NLO dynamics, a semi-quantitative exponential relaxation from the LO
post-transition spectrum towards equilibrium may be a useful alternative
Kainulainen:2022lzp ; Kainulainen:2021eki .. We studied several aspects of
this transition in a $O(N)$ scalar field theory with a time-dependent mass,
using the 2PI-1/N formalism at NLO.

We found that even for moderate couplings, the initial instability leads to IR
occupation numbers well in the classical domain $n_{k}>1$ before non-
linearities become important around $\mu t\simeq 2-4$, and direct comparison
of quantum and classical 2PI approximations confirm that the early stages of
the transition (say, $\mu t<50$) may be treated classically (LO dynamics
agrees qualitatively, but occupation numbers tend to overshoot). Although this
is as may be expected, it does imply that on this timescale the effect of all
the unstable modes with low occupation is either negligible or at least does
not introduce quantum corrections to the IR mode dynamics. During the initial
stages of the evolution, we argued that the IR environment may be quantified
through a classical effective temperature, which may be twice as high as the
overall temperature. This could have implications for IR-dominated out-of-
equilibrium processes taking place during this time.

As the system equilibrates kinetically on time scales $\mu\tau_{\rm
kin}\simeq\mathcal{O}(100-300)$, however, classical evolution is no longer
reliable. The quantum and classical evolution is substantially different, as
energy from the initial zero-point fluctuations begins to leak from the UV
into the IR, making classical occupation numbers too large. This is true also
for ”classical” modes with $n_{k}>1$. Hence from this time onwards we must
rely on quantum 2PI-NLO evolution (or higher truncation, when available).

We proceeded to study the kinetic equilibration process in some detail, noting
that full chemical equilibration is achieved at even later times somewhat
beyond our current reach (see, however Arrizabalaga:2005tf ; inprep ). The
equilibration time-scale for moderate couplings is linear in $N$ and of order
$\mu t_{\rm kin}\simeq 100-300$. In particular we pointed out that assuming an
instantaneous redistribution of potential energy onto (in this case)
$g^{*}=N-1$ massless non-interacting degrees of freedom fails to account for
quantitatively significant effect of the equilibration stage, the massive
mode(s), the interactions and the effective chemical potential, which lingers
for $\mathcal{O}(>1000)$ after the transition.

For application to the cosmological evolution, the equation of state of the
system is established very early (around $\mu t=150$), but must be computed
carefully. The system indeed behaves as a mixture of light Goldstone modes and
a heavy Higgs mode.

2PI-simulations remain numerically challenging. Because tachyonic preheating
involves very large occupation numbers, we made conservative choices in terms
of lattice size, discretization and the length of the memory kernel.
Technology exists for adding multiple scalar fields Aarts:2007ye ;
Tranberg:2013cka , fermions Berges:2002wr ; Shen:2020jya ; Berges:2009bx ,
Hubble expansion Tranberg:2008ae  and going to NNLO Aarts:2006cv ;
Aarts:2008wz  with robust renormalisation Berges:2004hn ; Berges:2005hc , and
applications in cosmology are diverse Berges:2002cz ; Gautier:2015pca ;
Serreau:2011fu . Implementing all or most of these improvements for realistic
models of reheating for larger volumes and the entire thermalisation process
is daunting, but within reach. This would be a welcome supplement to the
information provided by classical-statistical simulations of reheating, and
the post-inflationary evolution of the Universe.

## Acknowledgements

We thank Alexander Rothkopf for useful discussions and comments on the present
manuscript.

##  Appendix A Truncating the memory kernels

Figure 22: The effect of using reduced memory integrals of size $\mu
t_{k}=38.5$ and $200$ on the zero mode of $\Sigma_{F}$ (left), and the
particle spectrum at $\mu t=700$ (right). $N_{x}=16$, $N=4$, $\lambda=6$.

Throughout this work we used lattices with $N_{x}^{3}=32^{3}$ sites and a
resolution $a\mu=0.7$. Giving a physical lattice size $(L\mu)^{3}$ of
$(22.4)^{3}$. This discretisation provides 833 distinct momentum values
$(ak_{L})^{2}$, of which 14 are unstable $(ak_{L})<a\mu$. At our choice of
time resolution, $dt=\frac{a_{t}}{a}=0.1$ a timescale of $\mu t=700$
corresponds to 10.000 timesteps. As mentioned in the main text, the equations
of motions eq. 34 are integro-differential equations that are naturally solved
by a simple Euler scheme. The computation time is heavily dependent on the
lattice size and the time extent. In principle, the calculation of the memory
integrals requires the entire history of $F(t,t^{\prime})$,
$\rho(t,t^{\prime})$ which therefore needs to be kept in memory. On the other
hand, the magnitude of contributions to the memory integrals become smaller
far-in-the-past. In the interest of numerical efficiency, we may therefore put
a cut-off on the history to reduce the required memory. Since much of the work
is done computing Fast Fourier Transforms in the memory integrals, the runtime
scales as the $N_{x}^{3}\log N_{x}\times t_{K}^{2}\times t_{f}$, where $t_{f}$
is the runtime and $t_{K}=a_{t}n_{t}$ is the length of the memory kernel.

In the following we show the effect of this procedure. We choose a lattice
size of $N_{x}^{3}=16^{3}$ and compare simulations with a memory of $\mu
t_{K}=200$ and $\mu t_{K}=38.5$ as done in Arrizabalaga:2004iw . Figure 22
shows the zero mode of the real part of the self-energy $\Sigma_{F}$ at $\mu
t=700$ where, to a good approximation, kinetic equilibrium has been reached.
We see that the simulation with memory length of $\mu t_{K}=200$ shows the
expected damping over time, whereas the simulation with $\mu t_{K}=38.5$ shows
a cut-off effect. We also found that if the memory kernel is too short, the
simulations become unstable. Although figure 22 might raise concerns about
this approximation we found that the final states are the same for both, thus
confirming previously obtained results Arrizabalaga:2004iw . The final
distributions at $\mu t=700$ are shown in figure 22, and they are almost
indistinguishable, showing that neglecting far-in-the-past quantities provides
correct thermalization temperatures. The interaction strength and occupation
numbers impact the required kernel length, as large contributes to the memory
integral need to be kept for longer. Keeping only a limited history causes an
energy loss inversely proportional to the interaction strength (shown in
figure 23), because smaller coupling leads to large occupation numbers. The
studied combinations of parameters required substantially larger physical
memory length to complete correctly. All simulations have been performed with
kernel length of $\mu t_{K}=224$ and required $150$ GB of memory and about 50
hours of computation time on a computing cluster with 30 processing units.

Figure 23: Energy loss for simulations caused by reduced memory kernels of
size $\mu t_{K}=38.5$.

##  Appendix B Chemical equilibration at late times on small lattices

Figure 24: Left: Particle spectrum comparison after $\mu t=700$ and $\mu
t=2685$ and corresponding fits to Bose-Einstein distributions. The latter
shows not only temperatures and chemical potential closer to equilibrium but
also that particles for $\omega_{k}/\mu>1.5$ are also thermal. Right: $R^{2}$
value, quantifying how close the state is to a thermal one, calculated for
times up to $t\mu=2685$. We show fits, taking data into account for up to $\mu
t=700$ and $\mu t=2685$.

Figure 25: Evolution of the Bose-Einstein fit parameters $T$ and $\mu_{\rm
ch}$ in time, and the corresponding fits. The dashed line indicates the fit
region.

The complexity of the algorithm allows us to calculate the dynamics only for
times up to $\mu t=700-1050$, on a lattice size of $N_{x}^{3}=32^{3}$. As
shown in the main text this allows to quantify equilibration properties. In
this appendix we investigate the dependence of these estimates on having data
for even longer times. Simulations of until $\mu t=2685$ are accessible for
lattice sizes of $N_{x}^{3}=16^{3}$.

In figure 24 (left) we show the particle distribution at $\mu t=700$ and $\mu
t=2685$ and the corresponding fits to Bose-Einstein distributions. Note the
reduced number of data points compared to simulations with $N_{x}^{3}=32^{3}$
lattice, e.g. figure 9. As mentioned in the main text, particles in the IR
equilibrate quicker than in the UV. In the figure we see that at $\mu t=2685$,
particles for $\omega_{k}/\mu>1.5$ also align to the straight line, although
the data points are still not included in determining the fit parameters. On
the right of figure 24 we show the $R^{2}$ value that quantifies how closely
the state has equilibrated (similar to figure 9 only that now equilibration
fit values for times up to $\mu t=2685$ are available). We perform the same
analysis procedure as for $N_{x}^{3}=32^{3}$ lattices, namely discarding
values for which $R^{2}<0.97$ and fitting the rest with eq. 64. When analysing
data for up to $\mu t=2685$ we discard the first $\mu t=300$ data points. The
figure shows both fits, taking data for up to $\mu t=700$ and up to $2685$.
Taking data up to $\mu t=2685$ into account results in a kinetical
equilibration time scale of $\mu\tau_{kin}=393$, compared to $272$. Therefore,
taking only data up to $\mu t=700$ gives an equilibration scale that is about
$50\%$ smaller than the actual value. We note that the quoted error bars,
which are not of statistical nature but quantify a range of parameters that
are in agreement with the data drop from $7.6\%$ to $0.6\%$. Figure 25 shows
the time evolution of the Bose-Einstein fit parameters $T$ (left) and $\mu$
(right) with the corresponding fits. For large values of $\mu t$ the fit that
takes values up to $\mu t=700$ into account clearly deviates from the data.
The reheating temperature extracted from the left plot is determined to be
$0.86$, whereas the short simulation predicts $1.1$. Thus data up to $\mu
t=700$ gives a temperature that is about $30\%$ too large. Again, errorbars
drop from $6.1\%$ to $0.6\%$. The chemical equilibration time scale, taken
from the right plot in figure 25 is $\mu\tau_{ch}=660$, whereas the short
simulation gives $333$. Thus, the short simulation gives only $50\%$ of the
actual value. The errorbars drop from $82\%$ to $13\%$.

## References

  * (1) Alan H. Guth.  Inflationary universe: A possible solution to the horizon and flatness problems.  Phys. Rev. D, 23:347–356, Jan 1981.
  * (2) Andrew R. Liddle and D. H. Lyth.  Cosmological inflation and large scale structure.  2000\.
  * (3) Lev Kofman, Andrei Linde, and Alexei A. Starobinsky.  Reheating after inflation.  Phys. Rev. Lett., 73:3195–3198, Dec 1994.
  * (4) Gary Felder, Juan García-Bellido, Patrick B. Greene, Lev Kofman, Andrei Linde, and Igor Tkachev.  Dynamics of symmetry breaking and tachyonic preheating.  Phys. Rev. Lett., 87:011601, Jun 2001.
  * (5) Andrei Linde.  Hybrid inflation.  Phys. Rev. D, 49:748–754, Jan 1994.
  * (6) Juan Garcia-Bellido, Margarita Garcia Perez, and Antonio Gonzalez-Arroyo.  Symmetry breaking and false vacuum decay after hybrid inflation.  Phys. Rev. D, 67:103501, 2003.
  * (7) Jan Smit and Anders Tranberg.  Chern-Simons number asymmetry from CP violation at electroweak tachyonic preheating.  JHEP, 12:020, 2002.
  * (8) Alejandro Arrizabalaga, Jan Smit, and Anders Tranberg.  Tachyonic preheating using 2PI-1/N dynamics and the classical approximation.  JHEP, 10:017, 2004.
  * (9) Gert Aarts and Jan Smit.  Classical approximation for time dependent quantum field theory: Diagrammatic analysis for hot scalar fields.  Nucl. Phys. B, 511:451–478, 1998.
  * (10) Frédéric Dux, Adrien Florio, Juraj Klarić, Andrey Shkerin, and Inar Timiryasov.  Preheating in Palatini Higgs inflation on the lattice.  JCAP, 09:015, 2022.
  * (11) Marcos A. G. Garcia and Mathias Pierre.  Reheating after Inflaton Fragmentation.  6 2023.
  * (12) Rafid Mahbub and Swagat S. Mishra.  Oscillon formation from preheating in asymmetric inflationary potentials.  Phys. Rev. D, 108(6):063524, 2023.
  * (13) Mustafa A. Amin, Mark P. Hertzberg, David I. Kaiser, and Johanna Karouby.  Nonperturbative Dynamics Of Reheating After Inflation: A Review.  Int. J. Mod. Phys. D, 24:1530003, 2014.
  * (14) Gert Aarts and Juergen Berges.  Classical aspects of quantum fields far from equilibrium.  Phys. Rev. Lett., 88:041603, 2002.
  * (15) Anders Tranberg.  Quantum field thermalization in expanding backgrounds.  JHEP, 11:037, 2008.
  * (16) Gert Aarts and Juergen Berges.  Nonequilibrium time evolution of the spectral function in quantum field theory.  Phys. Rev. D, 64:105010, 2001.
  * (17) Juergen Berges and Jurgen Cox.  Thermalization of quantum fields from time reversal invariant evolution equations.  Phys. Lett. B, 517:369–374, 2001.
  * (18) Gert Aarts, Daria Ahrensmeier, Rudolf Baier, Juergen Berges, and Julien Serreau.  Far from equilibrium dynamics with broken symmetries from the 2PI - 1/N expansion.  Phys. Rev. D, 66:045008, 2002.
  * (19) Anders Tranberg and Gerhard Ungersbaeck.  Bubble nucleation and quantum initial conditions in classical statistical simulations.  JHEP, 09:206, 2022.
  * (20) J. Berges, K. Boguslavski, S. Schlichting, and R. Venugopalan.  Basin of attraction for turbulent thermalization and the range of validity of classical-statistical simulations.  JHEP, 05:054, 2014.
  * (21) John M. Cornwall, R. Jackiw, and E. Tomboulis.  Effective Action for Composite Operators.  Phys. Rev. D, 10:2428–2445, 1974.
  * (22) Gert Aarts, Nathan Laurie, and Anders Tranberg.  Effective convergence of the 2PI-1/N expansion for nonequilibrium quantum fields.  Phys. Rev. D, 78:125028, 2008.
  * (23) Alejandro Arrizabalaga, Jan Smit, and Anders Tranberg.  Equilibration in phi**4 theory in 3+1 dimensions.  Phys. Rev. D, 72:025014, 2005.
  * (24) J. Berges, Sz. Borsanyi, U. Reinosa, and J. Serreau.  Renormalized thermodynamics from the 2PI effective action.  Phys. Rev. D, 71:105004, 2005.
  * (25) Juergen Berges, Szabolcs Borsanyi, Urko Reinosa, and Julien Serreau.  Nonperturbative renormalization for 2PI effective action techniques.  Annals Phys., 320:344–398, 2005.
  * (26) Peter Millington, Zong-Gang Mou, Paul M. Saffin, and Anders Tranberg.  Statistics on Lefschetz thimbles: Bell/Leggett-Garg inequalities and the classical-statistical approximation.  JHEP, 03:077, 2021.
  * (27) Kimmo Kainulainen, Olli Koskivaara, and Sami Nurmi.  Tachyonic production of dark relics: a non-perturbative quantum study.  JHEP, 04:043, 2023.
  * (28) Kimmo Kainulainen and Olli Koskivaara.  Non-equilibrium dynamics of a scalar field with quantum backreaction.  JHEP, 12:190, 2021.
  * (29) Anders Tranberg and Gerhard Ungersbaeck.  In preparation.
  * (30) Gert Aarts and Anders Tranberg.  Thermal effects on slow-roll dynamics.  Phys. Rev. D, 77:123521, 2008.
  * (31) Anders Tranberg and David J. Weir.  On the quantum stability of Q-balls.  JHEP, 04:184, 2014.
  * (32) Juergen Berges, Szabolcs Borsanyi, and Julien Serreau.  Thermalization of fermionic quantum fields.  Nucl. Phys. B, 660:51–80, 2003.
  * (33) Linda Shen, Jürgen Berges, Jan M. Pawlowski, and Alexander Rothkopf.  Thermalization and dynamical spectral properties in the quark-meson model.  Phys. Rev. D, 102(1):016012, 2020.
  * (34) Juergen Berges, Jens Pruschke, and Alexander Rothkopf.  Instability-induced fermion production in quantum field theory.  Phys. Rev. D, 80:023522, 2009.
  * (35) Gert Aarts and Anders Tranberg.  Nonequilibrium dynamics in the O(N) model to next-to-next-to-leading order in the 1/N expansion.  Phys. Rev. D, 74:025004, 2006.
  * (36) Juergen Berges and Julien Serreau.  Parametric resonance in quantum field theory.  Phys. Rev. Lett., 91:111601, 2003.
  * (37) Florian Gautier and Julien Serreau.  Scalar field correlator in de Sitter space at next-to-leading order in a 1/N expansion.  Phys. Rev. D, 92(10):105035, 2015.
  * (38) Julien Serreau.  Effective potential for quantum scalar fields on a de Sitter geometry.  Phys. Rev. Lett., 107:191103, 2011.