Quantum tachyonic preheating, revisited

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.

A Truncating the memory kernels 26 B Chemical equilibration at late times on small lattices 27

Introduction
At the end of inflation [1,2], 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 [3].One of these is tachyonic preheating [4], 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 [5], where one field σ plays the role of the inflaton and a second field 1 ϕ 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 where σ is real-valued and ϕ is here written as a real singlet, but could have any number of components.The inflaton potential V (σ) is not specified but must have the property that the potential minimum is at σ = 0 and allows for slow-roll inflation at values of σ larger than µ/g.As the inflation rolls towards σ = 0, the effective mass of the second field ϕ changes sign and ϕ undergoes a symmetry breaking transition, as the potential acquires new minima at non-zero field value with a speed determined by the evolution of σ, u = 1 µ 3 dm 2 (t) dt σ=µ/g = 2g σ µ 2 . (1.4) We may fix the constant V 0 = 3µ 4 /2λ in such a way that the potential is zero at ϕ = v, σ = 0.
In hybrid inflation models [5], the end of inflation is itself triggered by the symmetry breaking transition, as the increase of ϕ 2 spoils the flatness of the σ-potential.In other implementations, inflation has ended well before the transition.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 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 σ field itself may also oscillate around σ = 0 so that the sign of m 2 (t) flips multiple times.Such oscillations may also lead to resonant preheating [3], where certain momentum modes in resonance with σ 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 [6][7][8][9][10][11][12].A review can be found in [13].
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 [9,14].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 [8].
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 . (1.5) In the Standard Model below T reh = 100 GeV, one expects2 g * ≃ 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 also quantify under what conditions we are allowed to ignore cosmological expansion in our simulations.The longest timescale in the problem is the chemical equilibration time τ ch , and so we can ignore expansion as long as Hτ ch < 1, where H is the Hubble rate.Since as long as M pl /µ ≫ µτ ch and the coupling is not too small, expansion is negligible.If not, the evolution equations have to be restated and solved in a Friedmann-Robertson-Walker metric as in [15].
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 [16][17][18].Early work considering such a transition include [19][20][21], with a substantial body of work based on the Hartree approximation (see for instance [22][23][24]), which in the 2PI-language is LO in a coupling expansion.Truncations at LO are very useful for the early stages of the evolution, and are numerically straightforward to implement.However, including only local (in time) selfenergies, they are unable to capture the quantum equilibration and ultimate thermalisation of the system (see however [25,26] for inhomogeneous systems).
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 and conclude in section 7. Some technical aspects of the simulations are placed in an appendix A. Some of the present work can be seen as a refinement and continuation of [8], and shares much of the notation.

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) This is a generalisation of the single real field ϕ 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) = −µ 2 , the potential has minima at ϕ a ϕ a = 6N µ 2 /λ ≡ v 2 .A common prescription is to choose a global rotation to fix this in the ϕ 1 -direction, ϕ a = δ a1 v. Then the mass spectrum becomes with a massive Higgs mode with mass m 2 H = 2µ 2 , corresponding to fluctuations in the ϕ 1 direction around v, and N − 1 massless Goldstone modes, corresponding to fluctuations perpendicular to ϕ 1 .These two modes are often explicitly encoded in a longitudinal (along the 1-direction) and transverse (to the 1-direction) propagator.
On the other hand, the system has O(N) symmetry, and so the equations of motion will conserve ⟨ϕ a ⟩ = 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 [8].Close to equilibrium, the compound propagator contains N − 1 Goldstone modes and one Higgs mode, and at equal times it may be decomposed as (2. 3) It follows that the expectation value v enters as a contribution to the zero momentum mode (2.4) We will return to this point below.

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 τ Q , as and the corresponding quench speed is given by u ).An instantaneous quench amounts to τ Q = 0, cutting out the transition region altogether (2.9) In the general case of a dynamical inflaton m 2 (t) = g 2 σ 2 − µ 2 , µ 2 0 corresponds to some initial value of σ.For simplicity, we will consider "symmetric" quenches where µ 2 = µ 2 0 .In the example of an instantaneous quench, modes with momentum k < µ evolve as which is the tachyonic (or spinodal) instability leading to exponential growth of the longwavelength modes.Performing a more detailed computation [7,19,21,22], using as initial condition the quantum vacuum prior to the quench ( one finds ) where ω ±2 k = ±µ 2 + k 2 .For a linear quench, the corresponding expression may be found in [6].The exponential growth continues until the self-interactions become important (see below).
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 selfconsistently dressed propagators, resumming an infinite set of perturbative diagrams.
3 Quantum evolution equations

The 2PI effective action
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 [14,27,28]).By construction, all quantum effects are then discarded.Simulating quantum dynamics in real time is highly nontrivial.One very successful approach is to derive evolution equations 4 for the propagator and the mean field from the 2PI effective action [17].We define where we have explicitly assumed a homogeneous state, but have not yet explicitly imposed that the mean field vanishes.The two-point function is time-ordered along the Keldysh contour C. The 2PI effective action may then be written in the form [29]: where for our model, The functional Φ[ φ, 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 δΓ[ φ, G] δG ab (x, y) = 0, and In what follows we consider a system initially in the symmetric phase φa = 0, for all a.Because of the O(N) symmetry of the equations of motion, if φa = 0 at the initial time it will remain zero for all time.The first equation is formally solved by in terms of the self-energy Multiplying (3.5) by G bc (y, z) turns it into an integro-differential equation, At leading order (LO) in a 1/N expansion we have [18] Φ Since G aa ∝ N , the entire contribution is O(N ) (see left-most diagram in Figure 1).At NLO we have where The expression (3.9) resums an infinite set of diagrams of O(N 0 ) (see Figure 1, [18]).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 [30], 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 (ρ) components as where and the sign-function again refers to ordering along the Keldysh contour in the complex time plane 5 .It follows that F (ρ) is symmetric (antisymmetric) with respect to spacetime indices (x, y).In particular, G ab (x = y) = F ab (x = y).An O(N ) symmetric state ( φa = 0) and spatial homogeneity allows to further simplify F ab (x, y) = δ ab F (x − y, t, t ′ ) and ρ ab (x, y) = δ ab ρ(x − y, t, t ′ ).We may similarly decompose the self-energy into two components Inserting this decomposition into eq.(3.7), one obtains the real-time equations of motion for F and ρ: 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 effective mass At NLO the non-local self energies are given by where the 1/N resummation is performed through the objects I F,ρ , that in turn satisfy6

Numerical implementation
The non-linear integro-differential equations of motion may be evaluated numerically, 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.

Discretization and observables
We introduce a cubic lattice of N 3 x = 32 3 sites with lattice spacing aµ = 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 The initial conditions for the spectral propagator follow from the basic commutation relations while for the statistical propagator we write Choosing n k = 0 and ω k = k 2 lat + m 2 (t < 0) amounts to selecting the free-field vacuum state prior to the mass quench.Conversely, for all later times, we may extract a particle number and a dispersion relation from7 ) Very far from equilibrium, the physical interpretation of these quantities is not clear, but once the system settles down enough that the field excitations are quasi-particle like, they are good representations of the occupation number and frequency of a given mode (see for instance [25]).
Real time evolution derived from the 2PI effective action have the advantage of conserving an energy functional.We monitor the energy and pressure densities during the simulations, from which we compute the equation of state parameter ω = P ρ .The expressions are given by and Initially, the system is in a free-field vacuum state with energy ρ 0 /N = V 0 and pressure P 0 /N = −V 0 , where V 0 = 3µ 4 2λ .Hence the initial equation of state parameter is ω = −1.Energy density and pressure are divergent quantities and subject to renormalisation as described in the following.

Renormalisation
The evolution equations are UV divergent, and although the lattice provides a regularization, a substantial cut-off dependence is present, calling for renormalisation.2PI-truncated evolution equations are renormalisable [32,33], but because of the infinite resummation of diagrams, a precise cancellation of all divergences comes at a substantial numerical cost [34].An approximate, but much simpler renormalisation procedure addresses only the local self-energies and quadratic divergences.In effect, given an input renormalised aµ = 0.7, we must perform the simulation with a bare mass parameter given by where the counterterm on the lattice at NLO is given by: At LO the prefactor changes N +2 6N → 1 6 to account for only the LO local self-energy eq.(3.16).This approximate subtraction works well8 for small-to-moderate values of the coupling, and not too close to the continuum limit (as discussed for instance in [8]).
The energy density and pressure density (3.28), (3.29) 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 and pressure density to be +V 0 and −V 0 , respectively.

The classical approximation
The classical limit is intuitively when the occupation numbers are large, n k ≫ 1, and indeed one may extend classical perturbation theory9 to this resummed non-equilibrium framework [9,14], in which case the classical limit amounts to This requirement is more general than n k ≫ 1, but since in the quasi-particle limit, F ∝ n k and ρ ∝ 1, the two are consistent when both apply.The 2PI-classical approximation consequently follows from discarding from the evolution equations (3.17, 3.18, 3.19, 3.20) all instances of ρ 2 (and ρI ρ ) when added to F 2 (or F I F ).For the system and truncation considered here, we find We note that at LO, the classical and the quantum evolution are the same, and so quantum effects only enter at NLO10 .It is also only at NLO that non-local effects such at equilibration and thermalisation enter, which is why the LO truncation (as well as the Hartree approximation) are insufficent for simulating preheating beyond the initial stages.
4 Early time dynamics To get us started, we will first consider the very early stages of the spinodal transition.The initial state for the transition is the quantum vacuum prior to the mass quench, and since occupation numbers are not large, the classical approximation would seem to not apply.However, because the quantum and classical evolution of the propagators coincide for free fields and in the LO approximation, at small coupling the classical and quantum evolution will agree for some time.
When the field self-interaction may be neglected, λ ≃ 0, the evolution equation becomes linear and closed-form solutions exist for the instantaneous (eqs.2.12, 2.13, [7]) and linear quench [6,22].For λ = 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 λ, at LO or NLO) becomes comparable to the mass parameter, with some number α. Often, the criterion is taken to be either α = 1 (the minimum of potential is reached) or α = 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 µ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.In figure 2 (right) we see that the inflection point of the potential (equivalently the minimum of the potential) is reached earlier with higher λ 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: The quoted error of the mean corresponds to varying the tolerance over an interval ∆F ∈ (0.1, 0.3).We can also 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 λ.Although the free field regime only lasts for a time µt ≃ 2 − 4 we see that the unstable modes are already highly populated n k ≫ 1, even for large values of λ.
At least for these modes, as we leave the early time, small-coupling regime, the classical approximation should be robust.

The classical regime
As we have seen, in the language of correlators, classicality follows in the limit where F 2 ≫ ρ 2 (see also [9] 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 for all modes.That requirement is fulfilled during a spinodal transition, at least in the range k L /µ < 0.8 and we may 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 and 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.Since self-interactions stop the instability when eq. (4.2) is satisfied, we can expect the individual occupation numbers to scale as n k ∝ 1/λ.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 λ, 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 λ, the maximal particle number occurs at times µt < 15.
5 Intermediate time dynamics

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 µt > 6, the non-local dynamics first entering at NLO becomes important.Figure 6 (left) shows the particle spectrum at several times during the evolution, and we see that for µt ≲ 40 the range of momenta with n k > 1 grows, as the high occupancy    in the IR is redistributed towards the UV.In figure 6 (right) we show the value of the largest momentum with n k > 1 (k max ), which at first grows and subsequently shrinks on a time-scale of µt = 500.Judging from the local correlator F (t, t) only (figure 7), the classical approximation for the whole system is valid throughout.
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 the out-of-equilibrium stage through an effective infrared temperature T IR , which we define to be the average temperature in the classical range, had the spectrum been classical equilibrium This is shown in figure 8 as a function of time, and see that although the asymptotic equilibrium temperature is T reh ≃ µ or less, at intermediate times the effective temperature may be twice as large or more.Provided the out-of-equilibrium process of interest takes place over a time-scale µt ≃ O(100 − 1000), this effective temperature should likely enter estimates of particle and relic production.

Kinetic equilibration
The (quantum-)NLO evolution equations provide correct thermalization, in the sense that at late times, the occupation numbers approach a quantum thermal equilibrium state [17,31,36].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 quantify the equilibration process by performing a fit of the particle spectrum with where ω k is the derived dispersion relation (3.26).We perform a least square fit with and quantify kinetic equilibration by calculating the R 2 value of this fit 12 In figure 9 we show the process of kinetic equilibration, for N = 4, λ = 6.On the left we see the particle spectrum at various stages of equilibration and the corresponding fit to a Bose-Einstein distribution.We restrict our fits to the IR defined by ω k /µ < 1.5, indicated by the dashed vertical line.In figure 9 (right) we see that the R 2 gradually approaches 1, and that although after the initial growth of particle numbers the spectrum is far from thermal, at intermediate times µt ≈ 150 the particles are reorganised into a thermal-like distribution.At this stage the distribution is well approximated by eq. ( 5.2).We define kinetic equilibration to be completed when R 2 reaches 0.97, which happens at µt = 150 in figure 9. Furthermore we can extract a kinetic equilibration time scale by fitting the R 2 evolution to the form 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.(5.5) (right).The two are consistent, and reveal a linear dependence on N and that larger coupling produces faster kinetic equilibration.At the onset of kinetic equilibration we extract a temperature, shown in figure 11, with a chemical potential in the range µ ch µ = 0.2 − 0.33.In the following stage of chemical equilibration the temperature approaches its final value and the chemical potential goes to zero.
In figure 12 (left), we show an advanced stage of kinetic equilibration, for N = 4, λ = 6.The distribution is very well approximated by eq.(5.2) (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 line is lower than the quantum spectrum, corresponding to larger occupation numbers.Only where yi refers to the data points, ȳ to the mean and fi 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.The best fit parameters T and µ in 5.3 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 reh and time scales.We stress that they are systematic rather than statistical.for n k > 3 do the two spectra agree.This implies that one should be wary of approximating intermediate to late time quantum evolution by classical dynamics.In figure 12 (right), we show the UV part of the spectrum13 at this same time of µt = 550.The quantum dynamics correctly conserves the zero-point fluctuations (n k + 1/2 ≥ 1/2, n k ≥ 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 an incorrect, classical equilibrium [27].In a manifestly classical simulation, one might consider not initialising with the zeropoint fluctuations, or only in the unstable modes [8].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 wrong.

Chemical equilibration
Figure 13: Evolution of the Bose-Einstein fit parameters T and µ ch in time, and the fit to eq. ( 5.6) and eq.(5.7).The dashed line indicates the beginning of the fit region.N = 4, λ = 6.
The fit parameters µ ch and T are time-dependent, and are expected to asymptote to zero and some finite reheating temperature T reh , respectively.Figure 13 shows T and µ ch in time, and we may further extract a final temperature T reh and a chemical equilibration timescale τ ch by a fit to (5.6) and µ ch (t) = µ final ch + (µ init ch − µ final ch )e −t/τ ch . (5.7) We must note here that we are somewhat limited by the available computational resources.
The asymptotic temperature T reh and chemical equilibration time scale µτ ch depend on the available time-range.Based on simulations of a lattice size 16 3 to much late times (µt = 2800), we find that data until µt = 700 overestimates T reh by about 20%.Furthermore, the chemical equilibration time is underestimated by about 80% (see Appendix B).With this caveat in mind, we show in figure 14 the reheating temperature T reh for various N and λ values.We also show for comparison temperatures from a simple toy model in which all available energy is distributed among a massive Higgs mode with m H = √ 2µ and N − 1 massless Goldstone modes.We obtain the temperature in this model by solving for T , where we assume equilibrated constituents The temperatures are almost independent of N .An even simpler model is to discard the Higgs mode, in which case one returns to (1.5) (with N − 1 degrees of freedom).In that case the temperature scales with T ∝ , with T /µ = 0.93 for λ = 6.Taking into  account the overestimation of T 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 µτ 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 reh and τ ch .

The compound propagator
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. ( 2.3) to the dispersion relation, while inserting the temperatures as obtained in the previous chapter.We indeed see that the dispersion relation is a composite of two particle modes, one with non-zero mass (Higgs, purple line) and one massless (Goldstone, green line).The field expectation value can be recovered from the zero mode of the propagator eq.(2.4), 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.

The equation of state
Figure 17: Time evolution of the equation of state for λ = 6.Left: Equation of state derived from the energy density and pressure density functionals.Right: Equation of state from a quasi particle ansatz.
For cosmological purposes, complete equilibration is often less important than the effective equation of state of the system, which often settles much faster [37].In figure 17 we show this quantity in time for λ = 6.One might have expected that the state was an equipartitioned mixture of g * "radiation" Goldstone components with P/ρ ≃ 1/3 and (N − g * ) "matter" components with P/ρ ≃ 0. In that case, the equation of state should be ω ≃ g * 3N . (5.9) At low temperature, one would expect g * = N − 1, so that for N = 4, ω = 1/4.But the situation is more complicated.On the left of figure 17 we show how the equation of state parameter evolves from the initial ω = −1 to its final value for λ = 6, within a time µt = 100−200.The dependence on N and the overall magnitude is not as naively expected.This is due to our renormalisation procedure, which calibrates to a potential minimum of V 0 , while the finite temperature effective potential has much shallower minima.If we instead compute the equation of state using a quasi-particle ansatz we find in figure 17 (left) that the equation of state indeed approaches its radiation value of 1/3 in the case of many Goldstone boson components (large N ).In figure 18 we repeat our simulation, but for λ = 1, where the expectation value is closer to the zero temperature vev.We see that both for the direct and quasi-particle (right) equation state, the system displays the correct behaviour.

Finite quench time
A natural generalisation of the instantaneous quench it to consider finite-time linear quenches of the form (2.7).Implementing this as a by-hand mass function has the drawback that total energy is no longer conserved.We straightforwardly have that 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, λ = 6.so that the total energy loss during the transition is For a very fast quench, ϕ a is not able to react to the quench of the mass parameter, and since we start out at ϕ 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 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.In a real dynamical realisation of the σ −ϕ coupling, the energy is strictly conserved and instead transferred to the inflaton σ.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 only briefly discuss the early time dynamics.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 in figure 20 (right) the particle spectrum at the time when the free field description starts to fail.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 µτ 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.Finally, in figure 21, we show the evolution of the classical range for different quench times.We see a separation into "fast" quenches µτ Q < 20, where the classical range is large at first before it adjusts, and "slow" quenches µτ Q > 50 where the classical range is smaller, established later and changes on a longer timescale.

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 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 scale 14 .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 14 When 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 [38,39].
numbers well in the classical domain n k > 1 before non-linearities become important around µt ≃ 2−4, and direct comparison of quantum and classical 2PI approximations confirm that the early stages of the transition (say, µ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 µτ kin ≃ 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) or only consider modes with n k ≫ 1.
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 [31,40]).The equilibration time-scale for moderate couplings is linear in N and of order µt kin ≃ 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 O(> 1000) after the transition.
For application to the cosmological evolution, the equation of state of the system is established very early (around µt = 150, see also [37]), 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 [41,42], fermions [36,43,44], Hubble expansion [15] and going to NNLO [30,45] with robust renormalisation [32,33], and applications in cosmology are diverse [46][47][48].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.

A Truncating the memory kernels
Throughout this work we used lattices with N 3 x = 32 3 sites and a resolution aµ = 0.7.Giving a physical lattice size (Lµ) 3 of (22.4) 3 .This discretisation provides 833 distinct momentum values (ak L ) 2 , of which 14 are unstable (ak L ) < aµ.At our choice of time resolution, dt = at a = 0.1 a timescale of µt = 700 corresponds to 10.000 timesteps.As mentioned in the main text, the equations of motions eq.(3.15) 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 ′ ), ρ(t, t ′ ) 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 3 x log N x × t 2 K × 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 3 x = 16 3 and compare simulations with a memory of µt K = 200 and µt K = 38.5 as done in [8]. Figure 22 shows the zero mode of the real part of the self-energy Σ F at µt = 700 where, to a good approximation, kinetic equilibrium has been reached.We see that the simulation with memory length of µt K = 200 shows the expected damping over time, whereas the simulation with µ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 [8].The final distributions at µ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 µt K = 224 and required 150 GB of memory and about 50 hours of computation time on a computing cluster with 30 processing units.to quantify equilibration properties.In this appendix we investigate the dependence of these estimates on having data for even longer times.Simulations of until µt = 2685 are accessible for lattice sizes of N 3 x = 16 3 .In figure 24 (left) we show the particle distribution at µt = 700 and µt = 2685 and the corresponding fits to Bose-Einstein distributions.Note the reduced number of data points compared to simulations with N 3

B Chemical equilibration at late times on small lattices
x = 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 µt = 2685, particles for ω k /µ > 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 µt = 2685 are available).We perform the same analysis procedure as for N 3 x = 32 3 lattices, namely discarding values for which R 2 < 0.97 and fitting the rest with eq.5.5.When analysing data for up to µt = 2685 we discard the first µt = 300 data points.The figure shows both fits, taking data for up to µt = 700 and up to 2685.Taking data up to µt = 2685 into account results in a kinetical equilibration time scale of µτ kin = 393, compared to 272.Therefore, taking only data up to µ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 µ (right) with the corresponding fits.For large values of µt the fit that takes values up to µ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 µ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 µτ 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%.

4. 1 Figure 2 :
Figure2: Left: Time evolution of the statistical propagator F (t, t, x = 0) in the free theory, LO and NLO approximation for N = 4 and λ = 1.The dashed (dotted) horizontal line corresponds to the minimum (α = 1) and inflection point (α = 1/3) of the potential, respectively.Right: Comparing the free field to the LO evolution for different λ.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.

Figure 3 :
Figure 3: Left: Time at which free theory approximation deviates from the LO result by 20 percent (∆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 for different values of λ.Even relatively high values of λ produce highly occupied modes at this very early stage.The lattice has the size N 3 x = 128 3 .

. 3 )Figure 3 (
Figure 3 (left) shows the relation between this time and λ.Based on a simple fit, we conclude that free theory is a good approximation until a time µt * = a − b ln λ where a = 4.13 ± 0.20, b = 1.00 ± 0.03.(4.4)

Figure 4 :
Figure4: 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, λ = 6.

Figure 5 :
Figure 5: Highest value for total particle numbers n tot for the LO and NLO truncation.Even for large coupling, n tot reaches large values.N = 4.

Figure 6 :
Figure 6: Left: The particle spectrum at the early stage of NLO evolution.N = 4, λ = 6.Right: The largest momentum wih n k > 1 in time, for different values of λ.N = 4.

Figure 7 :
Figure 7: Time evolution of F (t, t, x = 0) at NLO, in the quantum system and the classical approximation.

Figure 8 :
Figure 8: Time evolution of the effective temperature in the IR according to eq. (5.1).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.

Figure 9 :
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, λ = 6.

Figure 10 :
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.(5.5) (right).

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

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

Figure 14 :
Figure 14: Left: Extrapolated equilibrium temperatures for various values of N and λ.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 λ.

Figure 15 :
Figure 15: Left: Coupling strength dependence of T reh and power fit ∝ λ −x for which x = 0.37.Right: Coupling strength dependence of chemical equilibration time for N = 4.

Figure 16 :
Figure 16: Left: Dispersion relation for N = 4 and λ = 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.(2.4).

Figure 20 :
Figure 20: Left: Highest value for total particle numbers n tot for various quench times.Right: Particle spectrum in the LO approximation (N = 4, λ = 1) at the time when the linear approximation breaks down (∆F (t, t, x = 0) = 0.2).Shown are spectra of unstable modes (k L /µ < 1) for different values of quench times µτ Q .The spectra are calculated with a lattice of size N 3 x = 128 3 .

Figure 21 :
Figure 21: Time evolution of the maximum momentum with n k > 1 for different values of quench time.N = 4, λ = 6.

Figure 23 :
Figure 23: Energy loss for simulations caused by reduced memory kernels of size µt K = 38.5.

Figure 24 :
Figure24: Left: Particle spectrum comparison after µt = 700 and µ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 ω k /µ > 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µ = 2685.We show fits, taking data into account for up to µt = 700 and µt = 2685.

Figure 25 :
Figure 25: Evolution of the Bose-Einstein fit parameters T and µ ch in time, and the corresponding fits.The dashed line indicates the fit region.