Non-equilibrium dynamics of a scalar field with quantum backreaction

We study the dynamical evolution of coupled one- and two-point functions of a scalar field in the 2PI framework at the Hartree approximation, including backreaction from out-of-equilibrium modes. We renormalize the 2PI equations of motion in an on-shell scheme in terms of physical parameters. We present the Hartree-resummed renormalized effective potential at finite temperature and critically discuss the role of the effective potential in a non-equilibrium system. We follow the decay and thermalization of a scalar field from an initial cold state with all energy stored in the potential, into a fully thermalized system with a finite temperature. We identify the non-perturbative processes of parametric resonance and spinodal instability taking place during the reheating stage. In particular we study the unstable modes in the region where the vacuum 1PI effective action becomes complex and show that such spinodal modes can have a dramatic effect on the evolution of the one-point function. Our methods can be easily adapted to simulate reheating at the end of inflation.


Introduction
Classical scalar fields coupled to quantum matter play an important role in various settings in cosmology. They are used to study the creation of seed perturbations for structure formation, reheating processes, particle production and the creation of baryon asymmetry. Almost exclusively in these treatments it is assumed that the scalar field evolves in some classical, possibly quantum corrected but fixed effective potential. One rarely accounts for the backreaction of the non-equilibrium quanta that may be created during the dynamical process. However, such quanta may be produced copiously during out-of-equilibrium phase transitions [1,2] by parametric resonance [3][4][5][6][7] or by spinodal instability [6,[8][9][10][11][12][13], and they could significantly affect the evolution of the system [14][15][16][17][18]. In this paper we study the effects of quantum backreaction on the scalar field evolution using two-particle irreducible (2PI) effective action methods.
A crucial step in the rigorous analysis of the problem is performing a consistent renormalization of the equations of motion derived from the 2PI effective action. This is a highly non-trivial task, because in any finite truncation of the 2PI expansion, a number of auxiliary vertex and self-energy functions appear that require setting up consistent renormalization conditions [19]. Other works on the renormalization of 2PI-truncated theories include for example references [20][21][22]. In this paper we carefully go through the renormalization of our model using the method of cancellation of the sub-divergences [23][24][25][26]. We emphasize that while the renormalization counterterms are constants, the divergences that get subtracted, and hence also the vacuum state of the system, depend on the infrared physics, such as temperature, or even the shape of the non-equilibrium particle spectrum.
To be specific, we study a simple λφ 4 -model with a spontaneous symmetry breaking tree-level potential. We work in the Hartree approximation and perform the auxiliary renormalizations using the MS subtraction scheme. The renormalized equations of motion and the 2PI effective action are however scale independent and completely specified in terms of physical parameters. We present explicit results for the vacuum and finite temperature effective potentials as well as for the vacuum potential in the presence of non-equilibrium fluctuations. We stress that in the non-equilibrium case the effective potential can only be constructed a posteriori and it is not in general a useful quantity for solving the equations of motion.
With our renormalized equations we can follow in real time how the potential energy of the classical field is transferred into quantum fluctuations by the non-perturbative processes. We identify a strong parametric resonance, even though our self-coupled system is too complicated to admit a comprehensive analytical stability analysis. We also show that due to backreaction from spinodal instability the field can pass through a potential barrier even when starting with less energy than the initial barrier height. We also follow the full thermal history of a system that starts with pure potential energy, until it is fully thermalized with nearly all of its energy stored in thermal plasma. We also show that at the initial stages of reheating the quantum system is highly coherent, but the coherence is gradually erased by interactions as the system thermalizes.
This paper is organized as follows. In section 2 we review the 2PI effective action techniques and introduce our truncation scheme, the Hartree approximation. In section 3 we show how to self-consistently renormalize the 2PI equations of motion and express them in terms of physical quantities. We also study both resummed vacuum and thermal effective potentials in the Hartree case and compare them with other approximations. In section 4 we write our equations of motion in the Wigner space in terms of moment functions following references [27,28], and also complement the equations with phenomenological friction terms. Section 5 is dedicated to numerical results. We compute the evolution of various quantities, such as the classical field, particle number and coherence functions using the fully coupled 2PI equations. Finally, section 6 contains our conclusions.

2PI effective action and equations of motion
We will study the non-equilibrium dynamics of a scalar field theory with the potential V (φ) = − 1 2 µ 2 φ 2 + 1 4! λφ 4 using the two-particle irreducible (2PI) effective action technique of non-equilibrium quantum field theory [29,30]. The 2PI effective action for this theory is Similarly, the inverse classical propagator is given by Finally, Γ 2 consists of all 2PI vacuum graphs with lines corresponding to the full propagator ∆ and interactions inferred from the shifted action where φ = ϕ + φ q and φ q is the quantum field. The stationarity conditions of Γ 2PI will give the equations of motion for the one-and two-point functions: δΓ 2PI δϕ a = 0 and δΓ 2PI δ∆ ab = 0. (2.5) When the classical solution to the latter equation, parametrized in terms of ϕ, is reinserted back into the effective action, we formally recover the 1PI actionΓ In the full dynamical case the two equations are however strongly coupled and should be solved simultaneously, as we will do in our study. For the classical field ϕ + (x) = ϕ − (x) and we may drop the branch index and find: We also left the branch indices out from the local correlation function ∆(x, x), which is the same for all components of the two-point function ∆ ab (x, y). The stationarity condition for Figure 2. The first few terms contributing to Γ 2 , including their precise coupling constant dependences.
∆ ab (x, y) leads to the Schwinger-Dyson equation where summation over b is implied and the self-energy function is given by To proceed we also have to specify an approximation for the interaction term Γ 2 .

Hartree approximation
The first few terms contributing to Γ 2 , arising from the action (2.4), are shown in figure 2 (the role of the indices in the couplings is related to renormalization and will be explained in the next section). In this work we shall work in the Hartree approximation, which includes only the first term in the series, given by In this case the self-energy has only a singular or local part: while Π ab nsg (x, y) = 0. Obviously ∂Γ H 2 /∂ϕ = 0 as well, so there is no contribution to equation (2.6) in the Hartree approximation. We can now write the non-renormalized equations of motion compactly as Eventually we will move to the Wigner space defined in section 4 and solve these equations numerically in some example cases for homogeneous systems, but before we can do that, we have to address the divergences in ∆ ab and in particular in the local correlation function ∆(x, x).

Renormalization
Systematic renormalization in the context of the 2PI expansion was thoroughly discussed in reference [19]. Here we use the method introduced in reference [23], and later used in references [24,26], and we include also a connection to physical parameters. The key issue is that any finite order truncation of Γ 2 [ϕ, ∆] leads to an approximation forΓ 1PI [ϕ] that contains infinite resummations of 1PI diagrams and the associated counterterms. This gives rise to a number of auxiliary n-point functions which need independent renormalization conditions. These conditions can be defined by requiring that all sub-divergences cancel [23], but one needs to introduce a different renormalized parameter for each different operator.
To be precise, all n-point functions can be classified in terms of the number of classical fields that connect to them, and all functions that are connected also to propagator lines are auxiliary. Below we shall first renormalize the auxiliary n-point functions in the MS-scheme and show that the resulting 1PI action is independent of the renormalization scale. We start by defining the renormalized fields, propagators, couplings and masses: where the index, I = 0, 1, 2, 4, follows the power of the classical field associated with the n-point function. Written in terms of the renormalized quantities, the 2PI effective action becomes: where S[ϕ R ] is the same as in equation (2.2) with ϕ → ϕ R , µ 2 → µ 2 R(2) and λ → λ (4) R , and ∆ −1 0R is the same as (2.3) with ϕ → ϕ R , µ 2 → µ 2 R(0) and λ → λ (2) R . Moreover we defined the classical counterterm action and the inverse classical counterterm propagator where δ (I) ϕ ≡ Z (I) − 1 and the other effective counterterms are defined as: Also in the interaction term in (3.2) the renormalized couplings in the combination λ (I) R +δ (I) λ follow the power of the classical field in the interaction term (2.4), rewritten in terms of the renormalized quantities. The renormalized equations of motion can now be derived from the renormalized effective action, or more directly from (2.11a) and (2.11b), by writing the the non-renormalized quantities in terms of the renormalized ones: Here we suppressed the arguments in the local functions ϕ R (x) and ∆ R (x, x), as well as in δ (4) (x − y), for brevity. We now proceed to determine the various counterterms appearing in these equations and in the end find the renormalized equations of motion that include the effects of quantum corrections.
Auxiliary renormalization conditions. Because the operator acting on ∆ ab R in (3.6b) is independent of branch indices, we can concentrate on the time ordered component ∆ 11 R of the two-point function. We choose the mass-shell renormalization condition in the vacuum configuration ϕ R = v R , which simultaneously minimizes the effective action. That is, we These conditions imply that Z (0) = 1 in the Hartree approximation, and in our current scheme we can also set Z (2) = 1 (see footnote 2 below). The renormalization conditions (3.7) then become: where ∆ R (v R ) refers to the still divergent local correlator computed at the renormalization point. It should be noted that ∆ ab R is an auxiliary function and the parameter m 2 R is not yet related to any physical mass. Similarly, none of the couplings are yet related to observables, and there is considerable amount of choice related to their definition. We will choose the following conditions: 1 1 These choices are partly specific for the Hartree approximation, where the self-energy Π ab is proportional to the local correlation function. In any higher order 2PI truncation λ (0) R and λ (2) R would need to be renormalized separately.
Because Z (0,2) = 1 here, equation (3.9), together with (3.1) and (3.5) implies that λ (0) R = λ (2) R . Equation (3.9b) is less restrictive: it merely states that both renormalized mass terms are related to the same bare mass term. Conditions (3.9b) and (3.9c) still allow us to choose δ (2) µ and δ (4) λ such that m 2 R and λ (4) R can be matched to a physical mass parameter and a physical coupling. Using the conditions (3.9) and equation (3.8b) we can write equation (3.8a) simply as That is, we are able to keep the tree-level relation between the coupling λ (4) R , the background field v R and the mass parameter m 2 R .
Cancelling the sub-divergences. In order to proceed, we need to find out the divergence structure of the local correlation function. Using dimensional regularization we can write where ≡ 4 − d and 2/ ≡ 2/ − γ E + ln(4π) and Q is an arbitrary renormalization scale. We now separate ∆ R into divergent and finite parts as follows: In what follows we will suppress the Q-dependence of the function ∆ F0 for brevity. Next we insert the decomposition (3.12) back into equation (3.8b), use relations (3.9) and let the leading order terms define the renormalized mass independently from the terms containing divergences or counterterms. In this way we get two equations out of the equation (3.8b): Using definition (3.13) again in equation (3.14) and rearranging we get This equation has a consistent solution where the leading constant terms and the terms multiplying the combination (the sub-divergence part) v 2 R +∆ F0 cancel independently. This gives two constraint equations, from which we can finally solve the counterterms δ (2) λ and δ (2) µ : Scale dependence. The scale dependence of the auxiliary couplings and the mass parameters can now be worked out as usual by requiring that the bare parameters do not run: = 0. Using equations (3.17) one then immediately finds: The latter equation applies for both mass parameters, assuming that δ (0) µ and δ (2) µ differ by at most a finite and Q-independent term, which is the case in the Hartree approximation. Equations (3.18) can be easily integrated: Remember that as a result of equation ( On the other hand, the coupling λ (4) R does not run at all; indeed, to keep λ (4) R finite, we must have δ (4) λ = 3δ (2) λ up to finite terms according to equation (3.9c), which implies We shall see below that λ (4) R can be further fixed by some physical condition.

Renormalized equations of motion
It is essential to impose a correct treatment of the local correlation function away from the renormalization point in the equations of motion (3.6a) and (3.6b). Analogously to (3.13), we first define a leading order mass function that contains all finite terms in equation (3.6b): Here ∆ F is the finite part of the local correlation function, which must be defined analogously to equation (3.12): Note that both the finite part and the divergence contain non-trivial contributions from both the vacuum and the non-equilibrium fluctuations. Using definitions (3.21) and (3.22) we can write equation (3.6b) as follows: Then, noting that λ (4) R + δ (4) λ = −2λ (4) R + O( ), the same manipulations can be done also in equation (3.6a). This results in the final renormalized equations of motion: These equations appear deceivingly simple: when written for the Wightman function ∆ < R = ∆ +− R , equation (3.24b) takes the form of a wave equation with a time-dependent mass and, as we shall see in the next section, equation (3.24a) describes the motion of the onepoint function in a quantum corrected effective potential including backreaction from nonequilibrium modes. However, despite their apparent simplicity, the equations are strongly coupled through the local correlator in the gap equation (3.21) for the mass term.

Effective potential and physical parameters
Let us now consider the adiabatic limit of the evolution equations, where ∆ F is given purely by vacuum fluctuations with no physical excitations. In this case definition (3.21) reduces to and correspondingly Note that m 2 (ϕ R ) and ∆ R differ from definitions (3.21) and (3.22) only through a different value of the background field ϕ R . Using the equation of motion (3.6b) in the renormalized 2PI action (3.2) we can write down the 1PI effective potential in the Hartree approximation as follows: where V is the space-time volume and the tree-level effective potential is where in the last step we dropped all terms of order . Writing iTr ln ∆ R = V dm 2 ∆ R and using equation (3.26) one finds that the divergences cancel between the two last terms in equation (3.27) and the finite part of Tr ln ∆ R creates the one-loop correction to the effective potential. After a little algebra one finds the result: This is the vacuum effective potential in the Hartree approximation, found for example in reference [32]. Despite the apparent Q-dependence, V H (ϕ R ) is in fact scale-independent. Indeed, one can first show from definition (3.25), using (3.18), that ∂ Q m 2 (ϕ R ) = 0. Then by a direct differentiation and using equations (3.18) and (3.20) one finds that ∂ Q V H (ϕ R ) = 0.
Physical parameters. Differentiating (3.25) with respect to ϕ R one can first derive the identity Using (3.30) it is simple to show that the first derivative of the potential can be written as Comparing equation (3.31) with equation (3.24a) we can see that in the case of pure vacuum fluctuations the equation of motion can be written as ∂ 2 t ϕ R + ∂V H /∂ϕ R = 0. Differentiating equation (3.31) once more with respect to ϕ R one finds we now see that the on-shell mass parameter m R of the auxiliary propagator can be identified with a physical mass, 2 if we at the same time define This is the choice of parameters we shall use in the rest of this paper. So far we have defined the counterterm δ (4) λ only up to a finite constant. This, and other remaining freedom in choosing the counterterms (see footnote 2) would allow us to further match λ (4) R to some observable. Given that λ (4) R does not run, equations (3.32) and (3.33) are enough to fix the parameters of the theory. Going beyond the Hartree approximation would lead to more complicated calculations and relations between the auxiliary parameters, but would not change the derivation conceptually.

Finite temperature effective potential
In our derivation in section 3.1 we did not specify the finite part of the local correlation function ∆ F , and in what follows we will compute it numerically from the equations of motion. Before that it is useful to make one more observation concerning thermal corrections. Indeed, we can include thermal corrections by a simple generalization of equations (3.25) and (3.26): where I x = 2∂ x J x and J x is the dimensionless bosonic one-loop thermal integral Here the infinitesimal imaginary part iε defines the correct branch of the logarithm for a negative m 2 . With these definitions one can go through the analysis of the previous paragraph and show that the equation of motion of the homogeneous field is now is the thermally corrected, scale independent effective potential in the Hartree approximation: Note that in the 2PI approach also the vacuum part V H (ϕ R ) of the potential is computed with the thermally corrected mass, which is the solution to equations (3.34) and (3.35). It is worth mentioning that in each special case considered above, from the vacuum renormalization (3.12) to the quantum corrected effective action with (3.35) and without (3.26) thermal corrections, and finally to the general case (3.22), the divergence that gets removed by counterterms is different and depends on the value of the background field, the temperature and the particle distribution. This is an unavoidable feature of the 2PI equations with a finite order truncation. However, in all cases the counterterms themselves remain the same, uniquely defined constants.
Comparison to ring-resummed potentials. Equations (3.34) and (3.35) and (3.37) provide a consistent resummation of the thermal potential to super-daisy level. They can be seen as a consistent generalization of the Parwani resummation method [33]. In these approaches the thermal corrections affect all modes on equal footing, while in the usual ring resummation method [34,35] only the long wavelength modes are screened by the short wavelength modes in a thermal plasma. The advantage of the potential (3.37) is that it provides a consistently renormalized, smooth continuation between the non-relativistic and relativistic regimes. In all other ring-resummed potentials this behaviour has to be put in by hand.
To effect a fair comparison of different approximations, we write all potentials using the same renormalization conditions. To be concrete, we use the conditions With these conditions the standard one-loop corrected potential without the ring-corrections becomes where m 2 0 (ϕ R ) = −µ 2 R + 1 2 λ R ϕ 2 R and the standard one-loop vacuum potential is (this potential also satisfies the condition ∂ ϕ V 1−loop (v R ) = 0)  In the Parwani approximation [33] one replaces m 2 0 (ϕ R ) with the lowest order thermal mass m 2 0 (ϕ R , T ) = m 2 0 (ϕ R ) + 1 24 λ R T 2 in equation (3.38) and in the ring approximation [34,35], where only the zero-mode is dressed by thermal corrections, one finds: Above we wrote the Hartree potential in terms of the scale dependent variables. However, since the potential is actually scale independent, we can rewrite it at the scale Q = m R , explicitly in terms of the physical parameters: where m 2 T is the solution to the gap equation, which now becomes as can be shown by direct differentiation of equation (3.41).
In the left panel of figure 3 we show the evolution of the second ϕ-derivatives of the potentials near the critical temperature at ϕ R = 0. The sharp kinks seen in the ring (green dashed line) and Parwani (red dash-dotted line) cases at T = T 0 result from the non-analytic dependence of the resummed potentials on the thermally corrected mass term (we are using the high-temperature expansion for m 2 (ϕ, T ) in these schemes). The oneloop result (blue dotted line) does not share this feature, because there we are using the non-resummed mass term. Interestingly, the Hartree result (black line) does not show signs of similar non-analyticity. In the middle panel we show the evolution of the ratio v(T )/T , where v(T ) is the position of the asymmetric minimum as a function of T . There are significant differences between the approximations: in all resummed potentials a metastable minimum emerges, and it has the largest jump in the Hartree case. In the one-loop case the metastability does not develop, but there is a jump in v(T )/T at T ≈ 120 GeV due to the non-analytic behaviour, now of the vacuum mass term as a function of ϕ. In the right panel we show the potentials at the critical temperature (whose value for each model is given in the figure caption). The transition strength is dramatically different in the different approximations and it is by far the strongest in the Hartree case. Of course one should keep in mind that this is a very simple model, with only a single scalar field. However, when one compares the one-loop results with lattice calculations, one typically finds that both ring and Parwani approximations give weaker transitions than the lattice or 3d-perturbation theory calculations [36]. It would be interesting to see if the Hartree approximation was in better agreement with these schemes when applied in more complex models.

Wigner space equations
We now proceed to solving the general equations (3.24a) and (3.24b) for homogeneous nonequilibrium systems. Of these, equation (3.24a) is already in the desired form, when we assume that field ϕ R is homogeneous, but equation (3.24b) for the correlation function will be easier to handle in the mixed representation. Because of the homogeneity an ordinary Fourier transformation is sufficient for the spatial coordinates, but for the time variable we need the Wigner transformation: where t ≡ 1 2 (x 0 + y 0 ) and r 0 ≡ x 0 − y 0 . Because all correlation functions ∆ ab (x, y) have the same local limit, it suffices to consider the equation for the lesser Wightman function ∆ +− ≡ ∆ < . Starting from equation (3.24b), we find that in the Wigner representation it satisfies the equation, Here the index m in the derivative ∂ m t signals that the time-derivative acts only on the mass function and not on the propagator. Equation (4.2) is still equivalent to (3.24b) and highly complicated because of the infinite tower of t-and k 0 -derivatives involved. It can be recast into a simpler form by introducing a moment expansion. Following reference [27] we first introduce the moment functions: Then taking the real and imaginary parts of equation (4.2) integrated over k 0 and the imaginary part of the same equation integrated over k 0 and weighted by k 0 , we get three equations coupling the moments ρ nk with n = 0, 1, 2 to the field equation for a homogeneous field ϕ R (t): We used the shorthand m 2 eff (t) ≡ m 2 (ϕ R , ∆ R ) for the mass defined in (3.21) and (3.22) and defined ω 2 k (t) ≡ k 2 + m 2 eff (t). The gap equation (3.21), including the out-of-equilibrium (or thermal) modes, can be written explicitly as where we defined k ≡ 1 2π 2 ∞ 0 d|k||k| 2 , and the Heaviside theta-function θ ω 2 k (t) ensures that the vacuum does not contain the unstable spinodal modes. 3 If ρ 0k (t) is identified with a thermal distribution (including the vacuum part), equation (4.5) clearly reduces to (3.34). After discretizing the momentum variable, equations (4.4a), (4.4c), (4.4d) and (4.5) can be written as a closed set of coupled first order differential equations, which include backreaction from the non-equilibrium modes into the evolution of the homogeneous field ϕ R . The gap equation (4.5) must be solved first at the entry to the routine, after which the solution is advanced using (4.4a), (4.4c) and (4.4d). In practice one must introduce a UV-cutoff for the magnitude of the momentum |k|, but results should not depend on its precise value, because all non-trivial physics results from gradient terms acting in the infrared region. We have indeed shown that this is the case in our numerical examples.
Friction. Our main goal is to study the dynamical evolution of ϕ R including the backreaction from the modes excited during the zero-crossings (parametric resonance) and from the unstable modes (spinodal, or tachyonic, instability). We would also like to study dissipative interactions in our system. To do this correctly, we should go beyond the Hartree approximation. This would be in principle a straightforward but very laborious task. Some formal results can be found for example in [37]. Here we will instead add phenomenological friction terms to our equations. Following references [27] and [28] we generalize equations (4.4a)-(4.4c) as follows: where δρ nk ≡ ρ nk − ρ vac nk with ρ vac nk being the vacuum moments defined in equations (4.8) below, and the explicit forms for the equilibrium distributions ρ eq nk have to be provided externally depending on the problem. Collision integrals could be computed accurately in the context of the cQPA formalism following reference [28] (see also [38]), but here we are only interested in qualitative effects, for which the above phenomenological approach is sufficient. Even then the coefficients c i could be some momentum dependent functions, but for simplicity we will assume that they are constants. Note that ρ nk and ρ eq nk in general have different vacuum distributions due to different respective solutions to mass gap equations.
Number densities and coherence function. We can get a better understanding of the physical meaning of the moments by comparing them with the spectral cQPA solutions found in reference [27]. As explained in section 4.2 of reference [27], the moments are in a one-to-one correspondence with the cQPA mass-shell functions f m k± and the coherence function f c k . The former can be further related to the particle and antiparticle number densities n k and n k , so that one eventually finds [27,28]: The coherence functions f c± k measure the degree of quantum coherence, or squeezing, between particle-antiparticle pairs with opposite 3-momenta [39]. A non-coherent vacuum state must then be defined as a state with no squeezing in addition to having no particles. This corresponds to setting n k = n k = f c± k ≡ 0, which is equivalent to: where Θ k ≡ θ ω 2 k (t) . Because we are assuming that ϕ R is a real scalar field we also have n k = n k , which implies that ρ 1k = −1/2 at all times, so that the equation for ρ 1k is actually redundant. This is indeed a consistent solution even with the friction terms included.

Numerical results
We shall now solve the coupled dynamical equations (4.4d), (4.5) and (4.6a)-(4.6c) in a few examples chosen to illustrate the rich physics of the strongly coupled system including the quantum backreaction. We will uncover some known results and find new phenomena associated with spinodal and resonant particle production at phase transitions 4 . We will show that a strong spinodal instability can cause a quantum assisted barrier penetration without tunneling, and we emphasize the difficulty of giving any sensible definition for the effective potential in a non-equilibrium system. Eventually, we will follow the full thermalization process of a scalar field starting at rest in the vacuum potential until the end, when the energy in the field is almost completely transformed into thermal fluctuations. 5

Particle production and reheating via parametric resonance
We first consider a case where the field starts from a relatively large value and oscillates several times between positive and negative field values. Because we are also interested in the spinodal instability, we consider a tree-level potential with a negative mass term. As physical parameters we use m R = 100 GeV and λ (4) R = 1, given which, µ 2 R(2) (Q 0 ) can be solved from (3.13), while the running couplings and masses are defined in (3.19). We compute the initial value for the effective mass function m 2 (ϕ R , ∆ R ) using the vacuum Hartree approximation (3.25). We used running parameters everywhere in our calculations. This served as a useful consistency check, since all final results must be (and indeed were) scale independent. In this example we also set the friction terms to zero, c i = 0.
The essential results for this run are shown in figures 4 and 5. In the left panel of figure 4 we show the evolution of the classical field ϕ R , which here displays an orderly oscillation pattern with a slowly decaying amplitude. The middle panel of figure 4 shows the evolution of the fluctuations in the zeroth moment integrated over the 3-momentum, which is the non-equilibrium contribution to the local correlation function: k δρ 0k = k ρ 0k − ρ vac 0k ≡ δ∆ F (t, t). These results are in good agreement with reference [16], where this problem was studied earlier using the mode equation approach. The rapid increase of δ∆ F (t, t) at early times is caused by two non-perturbative processes, the spinodal instability and the parametric resonance.
Spinodal instability. The presence of a spinodal instability is manifest in the right panel of figure 4, where the effective mass term m 2 (ϕ R , ∆ R ) is seen to become periodically negative in the region t 0.25. Indeed, whenever the mass-function is negative, all k-modes satisfying k 2 + m 2 (ϕ R , ∆ R ) < 0 (5.1) are unstable and can grow exponentially. This is the spinodal or tachyonic instability. One might then be tempted to associate the growth in fluctuations in the period t 0.25 fully to the spinodal instability. If this was true, the excited modes should satisfy the condition (5.1), which here translates to |k| 60 GeV. However, from figure 5 we see that this is not the case. The fast production of modes is clearly visible in the upper panels which show the integrated particle number (left) and the integrated modulus of the coherence functions (right). But from the lower panels, showing time-momentum heat plots of the same quantities, we see that the excited modes are concentrated on a frequency band which lies entirely above the spinodal region (5.1).
Parametric resonance. While our equations are highly non-linear and strongly selfcoupled, it is apparent that the structures seen in the heat plots in figure 5 correspond to Mathieu instabilities associated with parametric resonance, familiar from the studies of inflationary reheating [4]. This problem was also studied using 2PI methods in reference [7], albeit with a different set of approximations and a different potential. If we identify the mass squared of the mode function in the Mathieu equation with our mass function m 2 (ϕ R , ∆ R ), and follow the analysis of section V in reference [4], we can (very roughly) estimate the Mathieu equation q-parameter in our case to be where ∆m 2 eff ≈ 2 × 10 4 GeV 2 is the instantaneous amplitude and ν ≈ 21 GeV is the local frequency of oscillations of the effective mass term m 2 (ϕ R , ∆ R ), shown in figure 4. The value of the q-parameter, which remains roughly the same throughout the calculation, suggests an intermediate resonance between the narrow and broad regimes. Similarly, the expected position of the first resonance band is by and large estimated to be This result, and the expected width of the resonance [4] ∆|k| ∼ |k| rb ≈ 60 GeV are in qualitative agreement with our results. In figure 5 we can even observe a second, much narrower band below the first one, which dominates the particle production at t ≈ 1.
While this is again in agreement with the qualitative expectations, its interpretation via Mathieu equation methods becomes even more tenuous. At late times t 0.3 the shape of the growth pattern fits well in the standard picture [4], but in the spinodal region the resonant production appears to be more efficient than usual: upon spinodal zero-crossings the resonant production that normally shows (as it indeed does at later times also in our example) a period of anti-correlation, is here always positively correlated. While individual growth bursts are not enhanced, this positive correlation leads to particularly strong particle production.
Because we did not include interactions in this run, the fluctuation band structure remains stable at all times. The system also remains highly coherent, as is evident both from the increase of the integrated coherence function and the stability of the heat plot of the coherence function shown in the right panels of figure 5.

Strong spinodal instability
In the above analysis we made little reference to the effective potential. Indeed, the oneparticle irreducible effective action is not a very useful quantity in an out-of-equilibrium setting and it can even be defined only after the equations of motion have been solved. Even then one cannot define it universally, but only as a quantity evaluated locally in time.
We will now study this question in the case of a very strong spinodal instability. To be specific, we still use the values m R = 100 GeV, λ (4) R = 1 and ∂ t ϕ R,in = 0, but we take ϕ R,in = 243.5 GeV and include also friction. We assume that collisions drive the system to the vacuum state, i.e. we take δρ eq nk ≡ 0, and we specify the coefficients to be c 1,2 = 0.6 GeV. 6 In this case the initial potential energy of the field is lower than the peak of the vacuum potential at ϕ R = 0. This can be seen in the right panel of figure 6, where we plot the Hartree-resummed vacuum potential (red dashed line) and indicate the initial field value by the black dot. Obviously, if the potential was held fixed, the field would simply oscillate around the positive minimum with a decaying amplitude. However, when backreaction is included, the picture changes dramatically. The actual field evolution is shown in the upper left panel of figure 6. Curiously, the field stays around the positive minimum during only one oscillation cycle, after which it apparently passes through the potential barrier, spending a rather long time near the middle of the potential with the effective mass function close to zero. Of course what happens is that in the first passage of the field into the spinodal region, an explosive creation of fluctuations takes place. This is clearly demonstrated in figure 7, which shows the integrated fluctuations in the moment functions (upper panels) and the associated heat plots in the time-momentum plane (lower panels). These fluctuations absorb a large amount of entropy, which decreases the free energy in the system and lowers the barrier between the minima allowing the field to pass to the negative side. The key issue is to not confuse the total internal energy of the system and the free energy, which may vary strongly 6 Although we gave the friction terms only in a qualitative form, we can provide an estimate for the magnitude of the ci-coefficients. From equations (4.6) it is clear that ci have the dimensions of mass. The lowest order contribution to the collision integrals arises at the second order in coupling in the 2PI expansion. Hence the naïve scale of the coefficients ci is given by λ  depending on the entropy production.
Non-equilibrium effective potentials. While the effective potential cannot be defined a priori, it is illustrative to construct it a posteriori as a time dependent potential that reproduces the equation of motion (4.4d) at all times. This potential can be constructed as the definite integral where ϕ R and ∆ R are the solutions of the equations of motion. We show this potential as the dashed black line in figure 6. After the crossing to the negative side, the shape of the potential function settles and the field oscillates around the negative minimum with a decaying amplitude. We stress that V 1PI is only useful for the visualization and interpretation of results and there is no unique definition of the effective potential in the non-equilibrium case.
As was already mentioned in section 3.3, in any finite truncation the renormalized 2PI vacuum becomes dependent on the IR-physics. Another interesting potential 7 function then is the equivalent of the vacuum Hartree potential in the presence of fluctuations. This potential is defined as where V H (ϕ R , ∆ R ) is the 2PI vacuum potential (3.29) evaluated replacing the vacuum mass function m 2 (ϕ R ) with the general mass function m 2 (ϕ R , ∆ R ). Note that the integral term over the fluctuations of the zeroth moment is a part of the vacuum Hartree potential, similarly to the case with the thermal potential (3.37). The potential (5.5) is shown with the blue solid line in the right panel of figure 6. It represents changes in the 2PI Hartree vacuum energy including the backreaction effects, and like the instantaneous V 1PI -potential, its barrier around ϕ R = 0 is temporarily lowered by the backreaction. This example demonstrates that the final stages of a phase transition may involve very complicated quantum dynamics, where classical expectations and constraints do not hold. We conclude this subsection by stressing on the difference of the fluctuation spectra in the present case, shown in the lower panels of figure 7, and in the parametric resonance case shown in figure 5. Even though we used the same mass and coupling parameters, essentially all fluctuations are here created by the spinodal instability. Indeed, they occupy a region in the phase space which is consistent with the instability constraint (5.1), continues all the way to zero momentum and lies entirely below the parametric resonance band.

Self-thermalization
As our final example we study thermalization of the scalar field energy in a self-interacting system. We use the same physical parameters and initial conditions as in section 5.1 but include collision terms with the friction coefficients c 0,1 = 0.6 GeV, and assume that the collisions drive the system to thermal equilibrium, i.e. we take δρ eq nk ≡ δρ th nk . With rigorously computed collision terms the thermal state would emerge automatically as an attractor solution, but in our phenomenological approach we need to give a definition for the instantaneous temperature. In thermal equilibrium a general moment can be written as where n BE (k 0 ) = (e k 0 /T − 1) −1 is the Bose-Einstein distribution function. In particular δρ th 0k = 1 ω k n BE (ω k ) and δρ th 2k = ω k n BE (ω k ). while δρ th 1k = 0. We define the equivalent temperature T = T (t) by requiring that the thermal state has the same energy as what is stored in the fluctuations: In all these equations ω 2 k = k 2 + m 2 (ϕ R , ∆ R ) is a function of time. The energy stored in the classical field is With our definitions of the temperature and the collision integrals the total energy H = H ϕ +H ∆ should be conserved, and we checked that this is indeed the case to a high accuracy in our calculations. For more details on this, and on the numerical setup in general, see appendix A.
Spinodal slowing. In the left panel of figure 8 we show the evolution of the classical field ϕ R . Initially ϕ R evolves as in the collisionless case, oscillating with a nearly constant frequency and a large amplitude, but around t ∼ 2 the frequency starts to decrease until it reaches a minimum around t ∼ 3. After this the field gets trapped around the positive minimum while the oscillation frequency increases again. This spinodal slowing effect was already seen in connection with the barrier crossing in section 5.2. The bearing of the spinodal modes is revealed in the inset in the left panel of figure 11, which shows that the effective mass term m 2 (ϕ R , ∆ R ) repeatedly becomes negative in this region. In the right panel of figure 8 we show the energy components H ϕ and H ∆ . Initially all energy is stored in the classical field, but the fraction of energy in the fluctuations increases until the system is reheated, with almost all of the energy contained in the fluctuations.  Figure 10. Shown are the momentum distributions k 2 2π 2 δρ 2k (left) and k 2 2π 2 f c± k (right) for three different times: t = 0.2 (solid blue lines) t = 1.3 (red dotted lines) and t = 6 (black dashed lines). Also shown in the left plot is the weighted thermal distribution k 2 2π 2 ω k n BE (ω k ) for the equivalent temperature T (t = 6) = 144.9 GeV (black dotted line).
Mode transfer and decoherence. In figure 9 we again show the evolution of the number density and coherence functions, including both the integrated quantities and the timemomentum heat plots. There are striking, but expected differences between these plots and the corresponding non-interacting results shown in figure 5. First, the number density stops growing already at t ∼ 1 and eventually starts to decrease for t 2. As is seen from figure 8, fluctuations dominate the total energy already for t 1, and the subsequent decrease of particle number results from a transfer of modes to higher energies. Thermalization process should also lead to decoherence, and this is indeed clearly visible in the upper right panel of figure 9, which shows the integrated function f c± k . From the heat plots we see that particle production gets progressively less efficient and moves to smaller frequencies, as less and less energy is left in the classical field. From the heat plot in the lower right panel we see that coherence is erased throughout the phase space at late times.
Thermalization. In figure 10 we show the |k|-distributions of δρ 2k (left panel) and the coherence function f c± k (right panel) weighted by the phase space factor, for selected times during the evolution. At a relatively early time t = 0.2 the distributions shown in solid blue still display a clear parametric resonance band structure. At a later time t = 1.3 (red dotted lines) the resonant spectrum is already much more complex, apparently with contributions from many narrow bands. Also a significant mode-transfer to the thermal region has already taken place. Indeed, from the main plot in the left panel of figure 11 we see that the equivalent temperature at t = 1.3 is roughly 140 GeV, and as the field is relatively light, m 2 eff 1/2 /T 1 with m 2 eff being the local average of the oscillating effective mass function, the expected maximum of the thermal spectrum is located at |k| ≈ 3T ≈ 400 GeV. At the end of the simulation, t = 6 (black dashed curve), the system has essentially thermalized. Almost all energy is in the fluctuations and very little particle production activity remains. The particle number in the resonance bands is small and the coherence is almost vanishing everywhere and in particular in the thermal region. Also the fluctuations in the equivalent temperature have but a small residual amplitude left. For the final time we also plotted (black dotted line in the left panel of figure 10) the equivalent thermal spectrum k 2 2π 2 ω k n BE (ω k ) with T = 144.9 GeV, corresponding to the equivalent temperature at t = 6. The close agreement between the actual and thermal distributions shows that the system has indeed thermalized to a very high accuracy. where H = H ϕ + H ∆ is the total energy and the total pressure P = P ϕ + P ∆ is similarly the sum of the pressures in the classical field and in the fluctuations. The former is given by where V H∆ was defined in (5.5). The pressure contained in the fluctuations can be computed as the spatial component of the energy-momentum tensor [27], and it can be written in terms of the moment functions as follows: It is easy to see that in the thermal limit (5.12) reduces to the negative of the thermal part of the effective potential in the Hartree approximation: P ∆ = −T 4 J m 2 T /T 2 . We plot the EOS-parameter w in the right panel of figure 11. The EOS-parameter starts from w = −1 and initially oscillates between w = −1, corresponding to total vacuum energy dominance, and w = 1, corresponding to kinetic energy dominance (kination) in the classical field sector. However, as the energy is moved out from the field and the system thermalizes, the EOS-parameter moves to the band 0 < w < 1/3 corresponding to normal matter. From the inset of the left panel we see that the average value |x eff | = |m 2 eff | 1/2 /T ≈ 0.6 at late times. This indicates that the reheated thermal plasma is almost relativistic and indeed, the EOS-parameter is asymptoting close to w = 1/3 at late times. (In a purely thermal plasma with x eff = 0.6 one would get w ≈ 0.315.) The periodic deviation below this value seen in figure 11 is due to the field contributions to energy and pressure. with a finite discretization. This problem can be avoided by a careful choice of binned variables for the vacuum distribution. That is, we replace the vacuum distribution by a coarse-grained distribution defined by an integration over each momentum bin q ∈ [q i , q i+1 ]: where q ci ≡ 1 2 (q i + q i+1 ), ∆q i ≡ q i+1 − q i and i 0 (q) ≡ 1 2 qω q − m 2 artanh q ω q . (A.5) When the bin width goes to zero, the replacement (A.4) does not make any difference. However, for a finite discretization it avoids the singularity that would occur in the spinodal region when the effective mass function coincides with one of the bin-momenta squared, m 2 (ϕ R , ∆ R ) = −q 2 ci .
Energy conservation. In figure 12 we show the relative change in the total energy δ H ≡ H/H 0 −1 in the example we studied in section 5.3. The total energy is H = H ϕ +H ∆ , where partial energies in the fluctuations H ∆ and in the classical field H ϕ were defined in equations (5.8) and (5.9). In this example the total energy should be conserved, and this is indeed true to a very high accuracy. In this run we used a discretized momentum |k| ∈ [0, 2000] GeV with 1000 grid points. As can be seen in the figure, the error is essentially negligible between the spinodal regions. Within the spinodal regions there is some residual noise at early times. This arises from the integrable singularity near m 2 (ϕ R , ∆ R ) = 0, even with the coarse grained binning, but even this error is small and can be further reduced by reducing the bin width. We conclude that numerical errors are well under control in our calculations.