Oscillating scalar dissipating in a medium

We study how oscillations of a scalar field condensate are damped due to dissipative effects in a thermal medium. Our starting point is a non-linear and non-local condensate equation of motion descending from a 2PI-resummed effective action derived in the Schwinger-Keldysh formalism appropriate for non-equilibrium quantum field theory. We solve this non-local equation by means of multiple-scale perturbation theory appropriate for time-dependent systems, obtaining approximate analytic solutions valid for very long times. The non-linear effects lead to power-law damping of oscillations, that at late times transition to exponentially damped ones characteristic for linear systems. These solutions describe the evolution very well, as we demonstrate numerically in a number of examples. We then approximate the non-local equation of motion by a Markovianised one, resolving the ambiguities appearing in the process, and solve it utilizing the same methods to find the very same leading approximate solution. This comparison justifies the use of Markovian equations at leading order. The standard time-dependent perturbation theory in comparison is not capable of describing the non-linear condensate evolution beyond the early time regime of negligible damping. The macroscopic evolution of the condensate is interpreted in terms of microphysical particle processes. Our results have implications for the quantitative description of the decay of cosmological scalar fields in the early Universe, and may also be applied to other physical systems.


Introduction
Scalar fields may have played a key role in shaping the observable Universe. Cosmic inflation [1][2][3][4] provides the (probably) simplest explanation for the overall geometry of the cosmos as well as the small perturbations that seeded the formation of galaxies. In a large number of models [5] the quasi-exponential expansion is driven by the potential energy of one or several scalar fields dubbed inflaton(s). The transition from the inflationary phase to the radiation-dominated epoch occurs when the inflaton dissipates its energy into other degrees of freedom through particle production that fills the Universe with a dense plasma [6][7][8][9][10][11][12]. Even though this process of cosmic reheating is of uttermost importance for a complete understanding of the cosmic history as it sets the initial conditions of the hot big bang, it has been studied far less than inflation itself, cf. [13][14][15] for reviews. A quantitative understanding of the reheating phase is highly desirable because it affects the mapping between cosmological observables and inflationary model parameters [16][17][18][19], and fundamental physics parameters [20][21][22][23]. Moreover, cosmological relics such as gravitational waves [24,25] or Dark Matter [26][27][28][29][30][31] can be produced during this period.
In many scenarios reheating occurs during coherent oscillations of the inflaton condensate around the minimum of its effective potential. For the most part this process takes place in a plasma of produced particles which practically constitute a thermodynamic ensemble described by a von Neumann density operator ̺. In the present work we are interested in the background evolution, in particular the evolution of the inflaton condensate defined as an average over both quantum and statistical fluctuations, Even though the underlying laws of physics are local, i.e. equations governing the evolution of the statistical ensemble are local (as is the von Neumann equation for ̺), when expressed in terms of a finite number of momenta of the statistical ensemble they generally exhibit non-Markovian (non-local) character. Such are the Kadanoff-Baym equations supplemented by the effective condensate equation that arise in the two-particle-irreducible (2PI) effective action formalism that we utilize here. In this work we shall specifically be analysing the non-Markovian equation, describing the relaxation of ϕ as it oscillates around the minimum, where for simplicity we neglect the expansion of the Universe. This equation occurs for a broad class of interactions with both bosons and fermions, and the information about specific interactions is encoded in the time translation-invariant integral kernels π R (t−t ′ ) and v R (t−t ′ ), which are closely related to the retarded self-energies and proper fourvertices in a given particle physics model. Equation (1.2) arises in the small field expansion, and the non-local terms can be seen as due to perturbative particle interactions and decays. Therefore, the applicability of this equation is limited to the mildly non-linear regime, where resonant particle production is inefficient, but where non-linearities cannot be neglected. Non-Markovian equations, such as (1.2) for the one-point function or similar equations for two-point functions, are analytically tractable in a limited number of simple cases, such as in the linear regime when all non-linearities can be neglected [32][33][34][35][36] and the equation can be solved exactly. However, this is not satisfactory for applications when condensate oscillations extend beyond the linear regime, which they typically do, and thus a better understanding of the non-linear evolution is highly desirable. In simplified models the problem can be addressed by resorting directly to lattice computations [37]. However, having a semi-analytic understanding of how the inflaton relaxes during reheating, and how this is affected by different types of interactions and non-linear evolution is clearly advantageous. This work is devoted to developing both qualitative and quantitative understanding of this problem.
The first part of the paper is dedicated to developing a systematic approximation scheme for solving non-Markovian equations of the type (1.2) in the mildly non-linear regime, based on applying the multiple-scale perturbation theory [38,39]. This approximation method is designed to describe systems where physical processes happen on different time scales. The system at hand exhibits such behaviour when the coupling constants are perturbatively small and when the kernels of non-local terms in Eq. (1.2) provide an effective window spanning no more than few oscillations: dissipative and dispersive effects at play are locally small, but build up over long times and many oscillations to appreciably alter the evolution. Such effects cannot be captured by the naive time-dependent perturbation theory, which is plagued by secular artefacts that grow without bound. Multiple-scale perturbation theory is designed to avoid such artefacts by effectively resumming naive perturbation theory, thus yielding uniform approximations valid for very long times. Even the leading order approximation captures the evolution with great precision, as we confirmed by comparing to direct numerical solutions of Eq. (1.2) for a number of cases. The solutions exhibit interesting features: (i) a power-law damping behaviour arises from the cubic non-local term in Eq. (1.2), which competes with exponential damping arising from the linear non-local term, and (ii) there are in general time-dependent frequency corrections to the oscillation amplitude.
The second part of the paper is devoted to considering other approximation schemes, and examining how they fare when compared to directly solving the nonlocal equation (1.2). More complicated systems might be described by a complex system of coupled non-local equations, which might be difficult to solve, even in the approximation scheme developed here. Instead of solving such equations directly, a much more popular and intuitive approach to non-Markovian equations is to derive their Markovian counterparts by making system-specific local approximations of the non-local terms. 1 These are Markovian quantum kinetic equations (effective Boltzmann type equations) for quantities that describe important properties of the system under consideration, such as on-shell single particle distribution functions, cf. e.g. [45][46][47][48][49][50][51][52]. In cosmology this approach has been heavily used in the treatment of baryogenesis (c.f. e.g. [53][54][55][56]), and the case of a scalar field we examine here [40,42,57,58], where typically the following Markovian equation is considered, where V(ϕ) is the effective potential, and Γ(ϕ) is the field-dependent damping coefficient, both containing thermal and quantum corrections. This equation is much easier to understand than the non-local one in (1.2), as it is a generalization of the Duffing's equation essentially describing a damped an-harmonic oscillator. While justified for adiabatically evolving condensate, it is not exactly the form of the equation one finds when carefully applying the Markovainisation procedure appropriate for oscillating systems [40], which also produces corrections to the kinetic term encoded by K(ϕ), 2 K(ϕ)φ + 1 2 K ′ (ϕ)φ 2 + Γ(ϕ)φ + V ′ (ϕ) = 0 . (1.4) Neglecting quantum and thermal corrections to the kinetic term effectively neglects some corrections to the inflaton oscillation frequency that might be relevant when analysing the spectrum of produced particles. The equation of motion (1.4) correctly reproduces the leading order solutions when solved using multiple-scale perturbation theory. It is worth noting that the Markovian equation cannot correctly account for the subleading corrections to the evolution, but it can nevertheless be used if one is not interested in that level of precision. We examine in detail the Markovianisation of Eq. (1.2) and point out ambiguities in its construction, and its limitations by comparing the solutions obtained using multiple-scale perturbation theory to the ones obtained by solving the non-local equation (1.2) directly. Lastly, we examine to which extent the time-dependent perturbation theory can describe the evolution of damped inflaton oscillations, and demonstrate that it can only be applied to early times when damping is approximately linear in time. It is in general not possible to extract how the oscillations dampen past this point, by using this method only. An exception to this is the linear regime, where one can correctly guess that the damping has to be exponential, and correctly resum the early time behaviour. This is not really an advantage because in the linear regime we can solve for the evolution exactly. It is important to note that in general we cannot resum the early time evolution in this way, as we do not know what functional form the damping must take.
In Sec. 2 we present an explicit derivation of (1.2) in a specific toy model starting from the 2PI effective action in the Schwinger-Keldysh formalism of nonequilibrium quantum field theory, and making a small field expansion. Explicit expressions for π R and v R in other models can be obtained by straightforward generalisations of the steps given there. In Sec. 3 we obtain analytic approximations for solutions of (1.2) by means of multiple-scale perturbation theory. In Sec. 4 we apply a Markovianisation procedure to (1.2) to obtain a local equation of the form (1.4). We solve this 2 For more general non-renormalizable effective interactions quantum and thermal corrections will produce Markovian equations with a more general structure, where X =φ 2 . Such equations can be seen to arise in Galileon theories [59] and Horndeski theories [60], hinting at a possibility that cosmologically interesting theories with non-canonical kinetic terms, such as e.g. k-essence [4] or kinetic gravity braiding [61], could arise effectively as quantumcorrected canonical theories.
Markovian equation using the very same approximation methods from Sec. 3, allowing us to establish limitations of Markovian equations. Lastly, in Sec. 5 we consider to which degree standard time-dependent perturbation theory is able to describe the relaxation of the oscillating scalar condensate, confirming it is able to do so only at early times. The concluding section is devoted to discussing and comparing different approaches considered in this work, with respect to their range of validity, their capability to compute subleading corrections, and the possibility to systematically include quantum statistical effects, and also to interpreting the relaxation of the condensate in terms of microphysical on-shell processes.
It is worth remarking that the methods considered and results obtained here, though motivated by a cosmological problem, are relevant for other fields, whenever a system is modelled as an oscillating dissipating scalar field. Such descriptions arise for example in studies of many-body theory, e.g. [62][63][64].

2PI-resummed effective action
Effective equations for the condensate evolution, such as Eq. (1.2), arise generically in a broad class of models. Here we consider a simple one, with two scalar fields Φ and χ, to demonstrate how to derive Eq. (1.2) in the Schwinger-Keldysh (in-in, closed-time-path) formalism [65][66][67][68], by constructing the 2PI-resummed effective action in the small field expansion. Our starting point is the 2PI effective action [67,69], dependent on the condensate ϕ(x) = Φ(x) (the expectation value of χ is set to vanish from the beginning), and the two-point functions ∆ ϕ and ∆ χ , where the kinetic operators are, and where Γ 2 contains all the 2PI Feynman diagrams with two and more loops. For simplicity we consider its two-loop truncation only, All the ± polarity indices carried by lines and vertices of the Schwinger-Keldysh diagrammatic expansion are implied and summed over when not explicitly denoted. Double lines stand for the connected two-point functions, , crossed circles stand for the ϕ insertion, and vertices stand for, As evident from definitions above, a given vertex is determined by how many solid and dashed lines connect to it, regardless of whether they are double lines, or single lines to be defined in Sec. 2.1 below. We are interested in writing down the effective theory for just the condensate. This is accomplished by first solving on-shell for the 2-point functions as functionals of the condensate, and then plugging these solutions back into the 2PI action, which produces the 2PIresummed effective action for the condensate, This is the form of the effective action that we derive in the following subsection, and is our starting point for studying the dynamics of the oscillating condensate, whose equation of motion then follows from the variation of (2.8), We are interested in describing the evolution of the condensate in the regime where non-linearities are small, but not negligible. A scheme well adapted to such a regime is the small field expansion, where the 2PI-resummed effective action is organized in powers of the condensate. In particular, we will be constructing Γ res. 2PI [ϕ] up to quartic order in ϕ.

Solving for two-point functions
Constructing the 2PI-resummed effective action requires us to solve for the two-point functions up to quadratic order in the condensate. We do this we diagrammatically, starting from the equations of motion (2.7) for the two-point functions, , (2.10) . (2.11) Note that vertices labeled by coordinates and Schwinger-Keldysh polarities are not integrated over, while the unlabelled ones are both integrated and summed over.
Since we are interested in the small field expansion of the 2PI-resummed effective action, it suffices to solve for the two-point functions above up to quadratic order. The leading contributions to the two-point functions in the small field expansion, which we denote by single solid and dashed lines, are independent of the condensate, (2.12) They are determined from Eqs. (2.10) and (2.11) where the condensate-dependent parts are dropped, , (2.14) with the kinetic operators on the left hand sides defined in (2.3a)-(2.3b). The solutions to these equations, (2.15) are free thermal equilibrium propagators, with the remaining components being complex conjugates, distribution, and where the thermal masses in the large temperature limit are, Note that the coincidence limits of equilibrium thermal propagators are independent of space-time coordinates, hence M φ and M χ are constants. Including the widths of the propagator would require going to two loops in Eqs. (2.13) and (2.14), generated by the three-loop effective action, which should be a straightforward extension. Solving for the condensate-dependent corrections now gives, . (2.20) Note that we have solved for the propagators up to (i) quartic order in the condensate at zero loops, and (ii) quadratic order at one loop. This is sufficient for the order to which we compute the 2PI-resummed effective action, which we obtain by plugging the solutions (2.19) and (2.20) back into the 2PI effective action (2.2), where we neglect writing constant terms which do not contribute to the dynamics. The reasoning behind the above truncation of the 2PI-resummed effective action is as follows. Different classes of diagrams induce qualitatively rather different behaviour for the condensate. The first big difference is the one between local and non-local diagrams. The local ones cannot induce dissipative effects, while non-local ones can. Local diagrams can only be quadratic in the fields, and we truncate them to the leading one-loop order in (2.21b). The non-local diagrams are ones where condensate fields appear at non-equal space-time points. The kinds of effects that non-local diagrams can induce depend on the power of condensate fields appearing in them. Here we keep lowest order quadratic non-local diagrams, which are the twoloop ones in (2.21c), and the lowest order quartic diagrams, which are the one-loop ones in (2.21d). Thus, in each qualitatively different class we keep only the diagrams of lowest loop order. The truncated 2PI-resumed effective action (2.21) is for our purposes conveniently written as, where M φ was given in (2.18a), and where, is the two-loop non-local self-energy and, , (2.24) is the one-loop non-local proper 4-vertex function, with the lines in the diagrams being thermal propagators (2.15). Note that, seemingly, at this order of truncation there is no dependence on the self-coupling λ χ of the χ-field. This is, however, not the case as the thermal masse (2.18b) of the corresponding resummed propagator depends on this coupling constant.

Homogeneous isotropic condensate
When the 2PI-resummed effective action (2.22) is specialized to a homogeneous and isotropic scalar condensate that we are interested in here, it reads, is an irrelevant spatial volume factor, and where, 26) are, respectively, the self-energy and the proper four-vertex in the mixed (two-time) representation evaluated at the vanishing spatial momentum. In the following we study the dynamics of the condensate encoded by this effective action, described by the equation of motion, 0 = δΓ res.

2PI
δϕ The condensate equation of motion takes the explicit form (1.2), where t 0 = 0 is the initial time at which initial conditions ϕ(0) andφ(0) are specified, and where the retarded self-energy and the retarded proper four-vertex function are, respectively, Important quantities that appear in the following sections are Fourier transforms of the retarded self-energy and the proper four-vertex, Note that we just as well could have written t as the upper limit of integration due to the step functions in (2.29) and (2.30). The diagrams making up the quantities above appear in different contexts in related studies, e.g. [42,44,58,[70][71][72].
For the two scalar model at hand (2.1), and at the level of truncation of the 2PI-resummed effective action (2.21), the imaginary part of the quantity v R (ω) can be computed analytically [73], The real part of the proper four-vertex has been studied in e.g. [74]. Even though this contribution is often neglected as it does not contribute to dissipation, we do take it into account here, and find that it does produce an interesting effect of timedependent corrections to the oscillation frequency of the condensate. For the comparison with time-dependent perturbation theory in Sec. 5 we shall require the T → 0 limit of expression (2.32) above, The leading contribution to π R (ω) comes from the "setting sun" two-loop diagrams in Eq. (2.23). We are not aware of any closed-form solutions for these diagrams. For the self-interaction piece there are estimates for its imaginary part in the following regimes [75,76], (2.34) For the gΦ 2 χ 2 /4 interaction one can expect a similar behaviour in the high-temperature regime.
The results for diagrams we quote in this section pertain to the simple two-field model from (2.1). If, for instance, χ happens to couple to other degrees of freedom, the collision terms in (2.28) in the high-temperature regime can look much more complicated, cf. e.g. Refs. [58,76,77].

Solving non-local condensate equations
In this section we consider a non-local (integro-differential) equation, corresponding to the quantum corrected equation of motion for the condensate (2.28) derived from the 2PI effective action, where M = M φ and λ = Λ φ , and t 0 = 0 is the initial time. We systematically construct approximate solutions of this equation in four cases of increasing complexity. The approximation scheme relies on a couple of assumptions. The equation of motion (3.1) was derived in the preceding section in the small field expansion, assuming that nonlinear terms can be considered small. In addition we make two more assumptions: (i) that the coupling constants appearing in kernels are perturbatively small so that the linear non-local term in (3.1) can also be considered small, and (ii) that the kernels of the non-local terms are such that they provide an effective window that spans no longer than several oscillation periods, which in thermal theories is typically true. Given that the last three terms in Eq. (3.1) are small by assumption, it might be tempting to approach the problem using standard perturbation theory. However, this approach fails completely for time-dependent problems, even when applied to familiar local equations (cf. [67] for an illustrative example). The reason behind this is a feature of secular growth of subleading corrections in perturbation theory applied to nonequilibrium problems. This secular growth causes the corrections to grow without bound, eventually becoming comparable to the leading approximation, thus invalidating the perturbative expansion. Such features are ubiquitous, and appear even for systems with a bounded evolution, and were explored in [78] for the case of an oscillating scalar condensate we consider here. These issues are artefacts of the perturbative expansion, where the solution to the previous order equation provides a resonant driving force for the equation for the following order, thus leading to secular growth. For these reasons a more sophisticated approximation scheme needs to be applied to Eq. (3.1).
Our method of choice is the multiple-scale perturbation theory (multiple-scale analysis) [38,39] (for a particularly lucid and insightful conceptual introduction to the method cf. Ref. [79]), which we adapt to non-local equations. This method is perfectly suited for treating time-dependent problems consisting of processes happening on different time scales, which is typically the case when some parts of the dynamical equation can be treated locally as perturbations, such as in the problem at hand. Even though locally small, these can build up over longer times to appreciably alter the evolution. It is this physical observation, together with realising the technical limitations of standard perturbation theory, that motivate the biggest conceptual leap of multiple-scale perturbation theory: the solution is assumed to depend on two times, ϕ(t) = F (t, τ = εt), and the two times t and τ are treated as formally independent, the former accounting for processes taking place on shorter time scales, and the latter for ones evolving on longer time scales. Mathematically, this turns our eqaution of motion (3.1) from an ordinary differential equation in single time into a partial differential equation in two times, and thus turns a well-posed initial value problem into an underdetermined one, since no additional initial conditions are being specified. It is this introduced ambiguity that is crucial when applying perturbation theory to the evolution in two times, as it is fixed by requiring that spurious secular growth be absent from subleading corrections, and for perturbative expansion to be valid uniformly for late times. This step effectively resums an infinite subclass of naive perturbative corrections.
In the interest of clarity, we treat all cases in an algorithmic form. All steps are explained in the simplest case example of the linear equation. For all other cases the algorithm is applied without further explanation. The approximate solutions for the condensate evolution in the cases we consider all take the same form, where Ω(t) is the local frequency of condensate oscillations, and the amplitude A(t) is a monotonically decreasing function of time describing damping of the condensate oscillations, conveniently expressed in terms of a local damping rate, The evolution of Ω(t) and Υ(t) depends crucially on which interactions we consider, and on the coupling governing those interactions. The physical phenomena exhibited in individual cases are commented at the end of each subsection.
Approximate solutions constructed using multiple-scale perturbation theory are checked versus direct numerical solutions of the non-local equation (3.1). Since the approximation scheme does not depend on the specificities of the two integral kernels (self-energy and proper four-vertex) in non-local terms, apart from them providing a finite window, for numerical checks we make use of simplified mock kernels instead, which satisfy the basic property of being anti-symmetric, The two specific kernels we use are the exponential and the Lorentzian, illustrated in Fig. 1. They provide a window of a typical inverse width a 1,2 > 0, and c 1,2 > 0 play the role of the coupling constants in the sense that kernels are normalized to it, What eventually appears in the analytic approximate solutions of Eq. (3.1) are the Fourier transforms of integral kernels; for the mock kernels they read, where the hyperbolic sine and cosine integrals are, respectively, and γ E is the Euler-Mascheroni constant.

Case 1: Linear equation
The first case we consider is the simplest, where in (2.28) we only keep terms linear in the condensate,φ The approximate solution is constructed in a series of steps described below.
(i) Identify small terms. The terms that can be treated as small are multiplied by a formal parameter ε for bookkeeping purposes. The solution is then constructed as an asymptotic series in ε, which is set back to one at the very end. It is natural to treat the non-local term as a small correction, (ii) Assume two physical time scales. When there are hierarchies between typical physical scales of a given time-dependent problem, it is reasonable to expect processes taking place on different time scales. For small contributions in the equation of motion it will typically take a long time until their influence has built up appreciably. This observation is formally implemented by assuming the system evolves in two independent timesfast time t and slow time τ = εt -so that the solution takes the following form, This assumption allows us to write the equation of motion (3.11) as, where we have used the chain rule to write out the derivatives. Note that inside of the integral we still have to consider that τ ′ = εt ′ .
(iii) Look for perturbative solution. Typically non-local equation such as (3.12) are not analytically tractable, and we have to resort to perturbation theory in the formally small parameter ε. The crucial step is to not treat the εdependence of the slow time τ as a perturbative parameter, but only the εdependence outside of it. Thus, we look for a power series solution of the form, (3.15) where the coefficient functions are, In a sense, we have assumed a partially resummed ansatz, where the τ dependence of the solution eventually accounts for the resummation of naive perturbation theory, extending the validity of perturbative solution to late times.
(iv) Assume limited extent of non-locality. If the non-locality effectively extends over a window of time which the solution varies only slightly in τ , it is appropriate to Taylor-expand the τ -dependence of the integrand, This effectively localizes the τ -dependence of the equation of motion, by recasting the non-local term as a power series expansion in ε which is introduced as an expansion parameter in (3.15), This step is novel compared to the usual multiple-scale perturbation theory applied to local equations. In physical terms, it corresponds to assuming the amplitude and the frequency of the condensate oscillations change only slightly during the window provided by the kernel of the non-local term.
(v) Organize equation of motion in powers of ε. Introducing the expansion (3.18) for the non-local term, and the perturbative expansion of the solution (3.15) into the equation of motion (3.14) allows to collect the terms of the same order in ε, which defines the perturbative tower of equations that have to be satisfied independently, and are in general easier to solve. The resulting equations are solved order by order, where the solution of the previous one provides a source for the following one. We need only the leading and first subleading order, (vi) Solve leading order equation. This is a harmonic oscillator equation, with a well known solution, Note that the leading equation is blind to the τ -dependence, resulting in R(τ ) being left arbitrary as a constant of integration. This is a crucial feature of multiple-scale analysis that plays an important role at the next step.
(vii) Simplify sources in subleading equation. The leading solution (3.21) now provides a source for the subleading equation, Bearing in mind the assumption of a finite extent of non-locality, here we make the approximation of neglecting any early time transient effects due to initial conditions. This amounts to eliminating explicit t 0 dependence of the source appearing in the lower limit of integration by taking t 0 → −∞, which allows us to recognize the integrals as Fourier transforms defined in (2.31).
It is convenient to define separately real and imaginary parts of the Fourier transform, (3.24) as they will be responsible for different physical effects.
(viii) Remove spurious resonances from subleading correction. The equation of motion for the subleading correction now reads, Note that the τ dependence appears only parametrically in this equation, which can be considered as an ordinary differential equation in time t. It describes a harmonic oscillator F 1 with a natural frequency M, being driven by an oscillating force of the same frequency M. It is well known that such driven oscillator exhibits resonant behaviour where its amplitude grows without bound. However, this resonance is not a physical effect, but an artefact of organizing the computation in a perturbative expansion. Physically, the amplitude of oscillations can only decrease in this system. This is where the assumption of two times comes into play. In order to avoid the spurious resonances that invalidate the perturbative expansion at late times and give unphysical answers for the evolution, we require that R(τ ), which is thus far unspecified, to be such that resonant driving forces are absent from Eq. (3.25), implying, This simultaneously fixes the ambiguity from the leading solution (3.21) and removes unphysical secular behaviour from the subleading order.
When solving for the complex function R(τ ), it is convenient to represent it in the polar basis in terms of two real functions, R(τ ) = A(τ )e −if (τ ) , which then breaks the complex equation (3.26) into two real ones, which are readily solved to give, where A 0 and f 0 are real constants of integration. Removing the spurious resonances from the subleading correction thus fully fixes the leading order approximation as it determines R(τ ). Note that here we did not have to solve for the subleading correction to accomplish that.
(ix) Obtain leading solution and determine initial conditions. Removing spurious resonances from the subleading equation in the previous step fixes completely the leading approximation, where A 0 and f 0 are constants of integration related to initial conditions ϕ(t 0 ) = ϕ 0 andφ(t 0 ) =φ 0 , This connection is obtained by evaluating the approximate solution (3.29) and its time derivative at the initial time t 0 = 0.
Comparison between the approximate solution (3.29) and a direct numerical solution of Eq. (3.11) is given in Fig. 2 demonstrating excellent agreement. This is not surprising for this example, since the linear equation can be solved exactly yielding the functional form captured by the approximation in (3.29). The effect that the linear non-local correction in the equation of motion (3.11) has on the oscillating dynamics of the condensate is a constant shift of the oscillation frequency, (3.31) and the damping of the oscillation amplitude which takes exponential form with the rate γ/2, giving a constant local damping rate, The approximate solution (3.29) has the precise form of the solution one obtains for a linear damped harmonic oscillator described by a local equation. This is not accidental, and is related to the assumption on the kernel providing a finite time window. There exists a local equation approximation to the problem which captures the leading solution we found here, which we discuss in detail in Sec. 4.

Case 2: Quartic potential
Taking into account the non-linear local potential in addition to the linear non-local correction, the equation of motion for the condensate reads, We find the approximate solution to this equation applying the algorithm developed in Sec. 3.1.
(i) Consider all non-linear and non-local terms as small corrections, (ii) Assume the solution depends on the two time scales, t and τ = εt, for which the equation of motion is, (iii) Assume a perturbative series solution of the form, The τ dependence of the integrand in the non-local term is expanded as, The equation of motion is organized in powers of ε as, The solution of the leading equation is, The subleading equation is, We neglect early time transient effects by taking t 0 → −∞, recognizing the Fourier transform from (2.31) so that the subleading equation reads, where we have defined, It is best to adopt a polar representation for the solution, R(τ ) = A(τ )e −if (τ ) , upon which the complex equation above breaks into two real ones, Integrating first the left equation, and then the right one gives, where A 0 and f 0 are constants of integration.
(ix) The leading approximation in this case is therefore, where A 0 and f 0 are related to initial conditions via, . (3.50) Comparison between the approximate solution (3.48) and direct numerical solutions of Eq. (3.33) is given in Fig. 3, showing excellent agreement. The local damping rate in this case is constant, just as in the case of the linear equation from the preceding section. This is not surprising since only the linear non-local term contributes to dissipation, while the local quartic potential does not. In fact, in the absence of the non-local term, it is well known that the quartic potential only shifts the oscillation frequency (cf. Duffing's equation in [38]), which is clearly seen in the early time limit, These oscillations do not continue at this frequency forever. As time progresses the amplitude is reduced due to dissipation from the non-local term. As the amplitude is reduced, so is the influence of the quartic potential, which leads to a time-varying effective frequency, which at late times becomes the frequency of the linear oscillator from Sec. 3.1. In fact, at late times the only information about the quartic potential that survives is the constant phase shift,

Case 3: Cubic non-locality
Consider now the case where the linear non-local correction in (2.28) is either absent or negligible, and where we keep only the cubic non-local correction,  (i) Consider the non-local term to be a small correction, (ii) Assume the solution depends on the two time scales t and τ = εt, so that the equation of motion is, (iii) Assume a perturbative power series solution of the form, The τ -dependence of the integrand in the non-local term is expanded in ε as, (v) The equation of motion organized in powers of ε reads, The solution of the leading equation is, The subleading equation is, We neglect early time transient effects by taking t 0 → −∞, where we recognized the Fourier transform from (2.31), so that the subleading equation reads, where we have defined, , (3.66) and α = α 0 +α 2 . Note that the imaginary part of the physical proper four-vertex is antisymmetric Im v R (−ω) = −Im v R (ω) , so that Im v R (0) = 0.
(viii) The terms in the first four lines on the right hand side cause spurious resonances in the subleading correction if allowed to appear, and we remove them by requiring, Adopting a polar representation for the solution, R(τ ) = A(τ )e −if (τ ) , breaks up the complex equation above into two real ones, Integrating these equations gives, where A 0 and f 0 are constants of integration.
(ix) The leading approximation for the condensate evolution in this case is therefore, where the constants of integration are related to the initial conditions via, Comparison between the approximate solution (3.48) and direct numerical solutions of Eq. (3.33) is given in Fig. 4, showing excellent agreement. The solution in (3.70) exhibits behaviour rather different than the solution of the linear equation (3.29). Firstly, it oscillates with a continuously changing frequency, Secondly, and more strikingly, the oscillations are damped, but the attenuation is a power-law, giving a decreasing local damping rate, This signals that the condensate dissipates less and less as time progresses, which is not strange if one notes that the non-local correction in the condensate equation (3.55) gets smaller as the condensate relaxes to smaller amplitudes. We emphasize that the exponential damping (constant local damping rate) is not a general feature of a dissipating system, but a feature of a linearized dissipative system, i.e. of linear response theory. For non-linear systems the form of the damping depends on particular interactions giving rise to dissipation.

Case 4: Full equation
Lastly, we consider the full condensate equation of motion (2.28), and construct the approximate solution by applying the algorithm outlined in Sec. 3.1.
(i) Treat all non-local and non-linear terms as small perturbations, (ii) Assume the solution depends on the two time scales t and τ = εt, so that the equation of motion is, (3.78) (iii) Assume a perturbative series solution of the form, (iv) The τ dependence in the integrands of the non-local terms is expanded in ε as, (v) The equation of motion organized in powers of ε, We neglect early time transient effects by taking t 0 → −∞,
(viii) Spurious resonances are removed from the subleading correction by requiring, Adopting a polar representation for the solution, R(τ ) = A(τ )e −if (τ ) , breaks up this complex equation into two real ones, These are readily integrated to give, where A 0 and f 0 are constants of integration.
(ix) The leading approximation for the solution of Eq. (3.75) is, where A 0 and f 0 are related to initial conditions for the condensate, This approximate solution captures very well the behaviour of the system, as it shows excellent agreement with numerical solutions given in Fig. 5. The solution (3.92) incorporates properties from all three previous cases. Most notably, the condensate exhibits both the exponential and the power-law damping, and the interplay of the two depends on the ratio between the two inverse time scales γ/(σA 2 0 ), as is seen from the local damping rate, At early times the damping is dominated by the cubic non-locality, +α ln 1+σA 2 0 t , (3.96) which matches the result (3.70) from the preceding section, with the effect of the local quartic potential included. Eventually the power-law damping is shut off as the amplitude of oscillations becomes smaller, and only the exponential damping is left. This is exhibited by the late-time limit, different regimes in case of a hierarchy between the dissipation coefficients, σϕ 2 0 /γ ∼ 10.
An important check of the approximate solutions we found in the four cases are the two assumptions mentioned at the very beginning of the section, • For all the non-local and non-linear effects to be small enough to consider them as perturbations we must have that the amplitude and the frequency of oscillations do not change much during one oscillation, which implies, γ, µ, σϕ 2 0 , αϕ 2 0 ≪ M.
• The window of the kernels has to be small enough so that slow processes can be neglected within it (steps (v) in the algorithm), meaning that the amplitude of the solution has to vary only slightly within the 1/a 1,2 extend of mock kernels. This is true if γ, µ, σϕ 2 0 , αϕ 2 0 ≪ a 1,2 .
Both conditions are satisfied for all the plots shown.

Comparison to Markovianised equations
In this section we consider a different approach to the non-local equation (1.2), which first calls for approximating it by a Markovian equation. In the first part of this section we examine in detail the Markovianisation method appropriate for an oscillating scalar system, and find it comes with certain intrinsic limitations stemming from the inability to formulate it as a systematic expansion in some parameter, and with an ambiguity related to perturbative manipulations of the local equation of motion. The latter issues can be removed by other considerations. The second part of the section is devoted to solving the Markovianised equation using the method of multiple-scale perturbation theory developed in the preceding section. We find that the leading approximate solution matches exactly the one found by solving the non-local equation directly in Sec. 3.4. This suggests Markovian equations can be used if we limit ourselves to precision where subleading corrections are neglected.

Markovainising non-local terms
Markovainising (localising) the non-local terms in the condensate equation (1.2) is performed utilizing the approach from [40], which relies on assuming that the kernels of non-local terms, π R (t−t ′ ) and v R (t−t ′ ), provide effective windows of some typical width t, beyond which they are suppressed, the precise details depending on the particular model and interactions (cf. e.g. [36,44,80] for a discussion). This approach was applied in [58,81] to the equation (1.2) we study here, but with some differences in the resulting Markovian equation which we comment on below. Adopting the assumptions given at the beginning of Sec. 3 -that field elongations are restricted to the mildly non-linear regime, that coupling constants are perturbatively small, and that kernel windows span no more than several oscillations -then, within a given window being integrated over, the condensate can be approximated by the tree-level solution, where ϕ 0 andφ 0 are initial conditions at the beginning of the window. Starting from this assumption, it is a simple matter of making use of addition formulas for trigonometric functions to derive the relation, which effectively Markovianises the non-local terms. For the linear term this is rather straightforward.
where the two parameters are, Markovianising the non-linear term superficially seems just as straightforward, as it involves only a square of the Markovianising expression (4.2), which Markovianises the non-local term, (4.6) where the three parameters are with the Fourier transform defined in (2.31), and where we have neglected early time the Markovian Eq. (4.10) when solving it in the following section. In [58,81] the equation (4.10) above is derived, but with keeping only γ and σ coefficients and neglecting coefficients µ, α 0 , and α 2 . The former are responsible for dissipative behaviour as they are not time-reversal invariant, while the latter are responsible for time-dependent frequency corrections to the condensate oscillations, as can be seen in Sec. 3 and as we demonstrate below. Even though the physical behaviour induced by these coefficients is different, they do arise as corrections of the same order to the evolution equation for the condensate.

Solving Markovian equation
We found in Sec 3.4 that the approximate solution obtained using multiple-scale perturbation theory captures the evolution of the system very well up to very late times. Here we apply the same algorithm to solve the Markovian equation (4.10).
(i) Identify small terms by introducing ε in equation of motion (4.10), (ii) Assume two time scales t and τ = εt , (4.12) and expand derivatives using the chain rule, (4.14) (iii) Assume a perturbative series solution of the form, (iv) The step of approximating the non-locality is already accomplished by the Markovianisation relation (4.9).
(v) Organize the equation in powers of ε. (viii) The subleading equation is, Note that the dependence on the arbitrary parameter ξ completely disappears at this level after solving the leading equation. Removing spurious resonances from the subleading equation implies, This is the exact equation for A(τ ) we derived in Sec. 3.4 when solving the non-local equation directly. In fact, the entire subleading equation is exactly the same. Therefore, the leading order approximate solution is the same as in (3.92), × cos f 0 + (M +µ)t + 1 σ We see that, even though the Markovianized condensate equation (4.10) contains ambiguities, and even though it is unclear how to write its subleading corrections, it does capture well the behaviour of the system when solved using the methods from Sec. 3. The ambiguity is fixed by requiring the equation to follow from the action principle, in which case it takes the form, which is precisely equation (1.4) advertised in the introductory section, with the three field-dependent functions, The issues with not accounting for subleading corrections is tied to the Markovianising relation (4.9), which we do not know how to generalize past leading order. Accounting for subleading corrections would require finding a more sophisticated Markovianising relation, which likely would require an analysis closely related to directly solving the non-local equation using multiple-scale perturbation theory, where at each level of perturbation the equations are local. Nevertheless, with this restriction in mind, Markovianised equations can be used to capture the behaviour of the full non-local system.

Comparison to time-dependent perturbation theory
Particle production from an oscillating scalar is often described by the standard timedependent perturbation theory. This approach has e.g. been used to study reheating in [84]. While this method does have its range of applicability, it misses several important effects if applied beyond this range. Namely, this method applies only to early times, as it neither captures the thermal effects nor the backreactions on the condensate oscillations. Therefore, in general it is not possible to infer how the condensate dissipates.
In this section we illustrate explicitly the limitations of the time-dependent perturbation theory. To this end we consider a simplified system in (2.1), where the self-interaction coupling constants are set to zero, λ φ = λ χ = 0, and only the coupling constant g of the bi-quadratic interaction is kept. We briefly go through the derivation of how one infers the dissipation of the condensate by considering the perturbative particle production as a power series in the coupling constant g. The computation is done at zero temperature for the sake of simplicity. Considering thermal effects within this approach is cumbersome at best, and they are often neglected within it. However, reformulating the approach in terms of particle distributions offers a way to incorporate thermal effects (cf. e.g. [85,86]), even though one has to pay attention to the early time restriction. Going beyond these requires additional assumptions, which start resembling more and more the approach we set up in sections 2 and 3. The Φ field in the action (5.1) is assumed to have a homogeneous condensate, which is very close to the classical condensate ϕ cl satisfying the classical equation of motion,φ So it is expanded as Φ(x) = ϕ cl (t) + φ(x), where φ captures its deviation from the classical behaviour. The χ field is assumed to be a spectator at the classical level, meaning it has no classical condensate (nor can it acquire it due to quantum corrections because of the Z 2 symmetry). The resulting action reads, The assumption of the time-dependent perturbation theory is that φ and χ can be treated as small perturbations. This is already where the restriction of describing the system at early times comes in. Quantum effects will introduce dissipation into the evolution of the condensate, which inevitably leads to large deviations from its classical behaviour. What we will be able to describe is the early-time damping of the condensate, where its amplitude is attenuated linearly in time. This is a generic feature of early time dissipation, as is evident from Taylor-expanding all the solutions from Sec. 3 for early times. Unfortunately, it is not possible to extract how the condensate dissipates at later times from this limited information, as it is not possible to resum the first two terms of the Taylor series into the full function.
Often it is assumed that the condensate amplitude dissipates exponentially, which is only valid for linear systems, but not in general. This assumption is based on our intuition on the damped harmonic oscillator, whose evolution is linear. We should not rely on this intuition to guide our conclusions in non-linear cases. The evolution of the system in (5.3) is solved for as a power series in g. The leading order then corresponds to treating φ and χ as free quantum fields, that are described in the standard canonical quantization, and their space of states given in the Fock basis of particle excitations. The free χ field operator takes the form, where ω k = k 2 +m 2 χ , and whereâ k andâ † k are annihilation and creation operators satifying canonical commutation relations, â k ,â † q = δ 3 (k−q), and which define the Fock space of particle states in the standard way. Analogous decomposition applies to the φ field operator.
We proceed by assuming that the initial state is the vacuum of the free theory, and ask the question of what is the probability that the evolution of the system excites particles from the vacuum, i.e. for the vacuum state to decay. This will be caused by the interaction terms in the second line of (5.3). In fact, at leading order in g only the first of these interaction terms will contribute, as that is the only one associated with a kinematically allowed process. As this term depends on the classical oscillating condensate, we will interpret it as the classical oscillations induced decay of the vacuum into particle states. The basic quantity needed to describe this decay is the transition amplitude for the initial vacuum state at t 0 = 0 to get excited to some other state at a later time t, where Dyson's time-evolution operator is, At leading order in the coupling g only the first term above contributes, and limits the possible out states to the two-χ-particle subspace. The only non-vanishing transition amplitude at leading order is computed straightforwardly making use of Wick contractions, where δ 3 (k + q) signifies spatial momentum conservation and we have substituted t 0 = 0. One of the usual quantities of interest is the probability density rate for the vacuum to decay to any other particle state, d dt where the 3-dimensional volume is V = (2π) 3 δ 3 (0). As already mentioned, at leading order the sum over the out states is restricted to the 2-χ-particles subspace, on which the decomposition of the identity operator reads, In fact, having cosmological applications in mind, it is better to consider the rate of change of the expectation value of energy density of produced χ particles, which is then given at leading order by, (5.11) These quantities are proportional to each other in the limit of Fermi's golden rule.
Computing the quantity in (5.11) is now done by assuming the time t to be much larger than one classical oscillation period, t ≫ 2π/m φ , which is the limit of Fermi's golden rule, and implies that we are considering quantities averaged over several oscillation periods. In that limit most of the terms in (5.11) average to zero, and the only non-vanishing contribution comes from the following term in the integrand, where the limiting definition for the δ-function has been used, sin(zx)/(πz) x→∞ − −− → δ(z). Making use of this, the rate of change of energy density is, where the averaging over classical oscillations is defined as This result is, in essence, Fermi's golden rule, in a slightly different form. What this result tells us is that a condensate oscillating with frequency m φ can produce two χ particles each of energy m φ = m 2 χ +k 2 , and with spatial momenta ±k. This can only happen if the kinematic conditions allow, which implies m φ > m χ , as indicated by the step function in the result above. Note that the energy density of the produced χ particles grows linearly in time, as suggested by the result in (5.13). Clearly this cannot proceed forever, and the result is limited to the regime where the energy density of produced particles is only a fraction of the energy density ρ ϕ stored in the oscillating condensate, ρ χ ≪ ρ ϕ . Moreover, this condition has to be satisfied over many oscillation periods for the approximations going into the result (5.13) to be valid. Our next task is to infer what is the correction to the evolution of the scalar condensate in this regime, so we can compare it to the results of Sec. 3.3. This can be done by appealing to energy conservation, which tells us that,ρ ϕ = −ρ χ . (5.15) Assuming that the energy density of the condensate ϕ is given by the classical functional, 16) and that condensate evolution at leading order in g takes the form, where δϕ(t) is of order g 2 , we get by linearizing the energy conservation equation in δϕ, where the left hand side is averaged over several oscillations in the same sense as the result in (5.13), and where we have used that the purely classical contribution to the energy density is conserved by itself, and that the classical contribution to the condensate satisfies the classical equation of motion. One can see that the averaged equation (5.18) is solved by the following ansatz, that vanishes at initial time. This collapses the equation into an algebraic one, providing a solution for one of the constants, Note that due to averaging over oscillations it is not possible to determine constant b from ansatz (5.19). Therefore, the evolution of the condensate we are able to infer using the method of time-dependent perturbation theory is Yet again it is clear that this result applies to early times only, in particular it applies in the regime where It is worth noting that the second term b is often not considered at all, in which case the first term is readily identified as damped classical oscillations. Comparing this result to the zero temperature limit of the result (3.70) from Sec. 3.3, expanded at early times, we see that (i) the early time correction that we were not able to determine is nonvanishing and corresponds to frequency correction of the condensate oscillations, and (ii) the early time correction we did determine indeed corresponds to damping provided that σ from Sec. 3.3 matches the appropriate coefficient computed here, which indeed is the case when the zero temperature limit of the imaginary part of the proper 4-vertex diagram (2.33) is plugged into the definition of σ. This illustrates that the method of time-dependent perturbation theory can capture the behaviour at early times. However, it neglects any kind of feedback effects, including the decrease in the condensate amplitude due to dissipation, the time-dependent frequency shifts due to the non-linearities in the system, and the impact of the produced particles on the damping rate. The restriction to early times comes from assuming that the condensate evolution is approximately the classical one at leading order. Inferring the evolution at late times, when the condensate has dissipated appreciably, would require resumming an infinite subclass of subleading corrections in the Dyson series (5.6). This could possibly be circumvented by making more sophisticated assumptions. Namely, that the condensate evolution has to take the form ϕ(t) = A(t)ϕ cl (t), and that the amplitude correction can be effectively treated as approximately constant locally, i.e., inside of the integral in (5.8). This assumption is very close to the assumption of step (iv) in the algorithm of Sec. 3, but it still neglects the frequency correction. On top of this assumption, one would still need to construct an equation of motion for the condensate amplitude A(t), which would yield the late-time evolution. Whatever this procedure might be, if correctly done, this equation has to take the form of equation for A(τ ) in (3.68) which follows from multiple-scale perturbation theory. It becomes clear that accomplishing this would be far more cumbersome than simply setting up the problem in the non-equilibrium 2PI formalism of Sec. 2, and solving it using methods of Sec. 3 or Sec. 4. Moreover, the time-dependent perturbation theory as presented here assumes an initial state without particles in (5.5). This cannot account for important effects that modify the time evolution of ϕ, including quantum statistical effects, screening, and scatterings of particles off the field ϕ. Implementing these here would be cumbersome, but the 2PI approach systematically incorporates them from the beginning through the quantum statistical averages introduced in (1.1).
By construction the result (5.22) has a simple physical interpretation as a decay of the "vacuum" 0; t 0 to a two-particle state 1 χ k 1 χ q ; t , induced by the external field ϕ cl (t). However, the split into the classical condensate and fluctuations, Φ(x) = ϕ cl (t)+φ(x), is not a physical one since φ = 0. This suggests that we view ϕ cl as a leading contribution to the full condensate ϕ = Φ . The classical condensate is then seen as quantum coherent state built out of a large number of zero mode quanta φ, with an average number of quanta per unit volume being, Now, the production of χ particles induced by the external field can better be seen as the zero mode quanta φ annihilating into χ particles, φφ → χχ. This process depletes the coherent state making up the condensate, whose evolution then necessarily must deviate from ϕ cl after some time. The zero mode quanta can also be seen as particles at rest with their energy equalling their mass m φ (but they also have to be seen as completely delocalised). Thus the energy of the process where two φ quanta annihilate is 2m φ , which explains the kinematic threshold in (5.8). The coefficient σ from (5.25) carries the information about the damping of the oscillating condensate amplitude at early times, and we have just argued it is constructed from an amplitude for φφ → χχ decay. This is a special case of a general result of cutting rules that follow from the optical theorem. The coefficient σ is proportional to the imaginary part of the second diagram in (2.21d), which precisely corresponds to this on-shell decay.

Discussion and conclusion
Summary. In this work we studied the damping of scalar condensate oscillations in the mildly non-linear regime in the formalism that naturally accommodates quantum statistical effects of the thermalised medium where dissipation takes place. Thermal corrections to the damping rate are known to be important for applications in cosmology, in particular for the thermal history of the reheating epoch [20,22,26,71,76,[87][88][89][90][91], or the production of axions [92,93] and sterile neutrinos [94]. Our starting point was the non-local equation of motion (1.2) following from the 2PI-resummed effective action (2.21), derived as a loop expansion and a small field expansion in the 2PI effective action formalism of non-equilibrium quantum field theory. Even though in Sec. 2 we derived Eq. (1.2) for a particular two-scalar model given by action (2.1), the form of the equation is rather generic, and captures a broad class of interactions with both bosons and fermions. The information about particular interactions is encoded in the integral kernels π R (t−t ′ ) and v R (t−t ′ ), which are closely related to the retarded self-energies and proper four-vertex functions in a given particle physics model. In Sec. 3 we have developed a systematic approximation scheme to directly solve the non-local equation (1.2) describing a damped oscillating scalar condensate. We accomplished this by applying multiple-scale perturbation theory, and were able to find a leading-order analytic approximation in (3.92) that captures the full solution very well, as demonstrated by comparing to exact numerical solutions in Fig. 5. However, solving non-local equations directly is notoriously difficult, and it is often advantageous to approximate them by local equations. For that reason in Sec. 4 we considered the Markovianisation procedure [40] applied to our non-linear equation (1.2). We found an ambiguity in the resulting equation that we were able to fix by requiring the equation to follow from an action principle. Solving the resulting Markovian equation using the same methods of multiple-scale perturbation theory from Sec. 3 produces the same leading order solution (3.92), justifying the use of Markovian equations. Examining the more common time-dependent perturbation theory approach in Sec. 5, we found it is neither capable of describing the evolution beyond early times when feedback effects are very small, nor does it naturally incorporate thermal statistical effects as the non-equilibrium 2PI formalism of Sec. 2 does.
Analytic approximate solution. The main feature of the approximate solution (3.92) of the non-linear non-local equation (1.2) exhibits are the two competing regimes of damping -a well known exponential one dominated by linear evolution at late times, and a power-law one dominated by non-linear evolution at early and intermediate times -conveniently expressed in terms of a local damping rate, defined in (3.3) and computed in (3.95), Another noteworthy feature is the time-dependent frequency of oscillations, given in Eq. (3.98). Depending on the parameters of the model the power-law damping regime can persist for a significant amount of time, before the condensate amplitude attenuates enough, and oscillations inevitably transition into the exponential damping regime. Thus, the power-law damping regime can leave imprints on the late time linear evolution in the form of amplitude correction and frequency shift. An example of such transition between damping regimes in given in Fig. 6. We emphasise that damping of the oscillations beyond the linear regime in general differs from the exponential damping found for linear oscillators, and depends on the particular dissipative interactions.
Thermal effects and microscopic interpretation. The damping of the condensate (6.1) is a macroscopic effect, but it derives from microscopic dissipative processes. At the end of Sec. 5 a picture of the condensate made up of a large number of zero mode quanta φ decaying into χ-particles was presented. A similar, but more intricate picture applies to the condensate in a thermal medium of non-vanishing particle numbers that is characterised by a statistical density operator ̺. In such a medium the physical excitations are quasi-particles that have modified properties compared to particles in vacuum; in Sec. 2 we have recounted how the first such important modifications arise -thermally generated quasi-particle masses, that can significantly modify kinematics of various processes. Moreover, new channels of dissipation open up in a thermal medium due to the possibility that zero-mode quanta φ making up the condensate interact with quasi-particles in the medium. The connection between the macroscopic behaviour of the condensate and the microphysical dissipative processes is given by the dissipation coefficients γ and σ appearing in (6.1), which are related to imaginary parts of the self-energy and the proper four-vertex as in (3.87). The physical interpretation of these quantities is systematised by the cutting rules at finite temperature [95][96][97][98][99][100]. 5 In short, these cutting rules applied to diagrams in (2.21c) related to the self energy π R , and the ones in (2.21d) related to the proper four-vertex v R suggest that processes contributing to the γ dissipation coefficient consist of φ ⇄ χχφ and φ ⇄ φφφ and rearrangements thereof (e.g. scattering χφ → χφ), and that processes contributing to the σ dissipation coefficient consist of φφ ⇄ χχ, φφ ⇄ φφ and permutations thereof. The use of the 2PI formalism for non-equilibrium quantum field theory ensures all the processes contribute with appropriate weights determined by the Bose-Einstein distribution. The number of zero mode quanta participating in the microscopic processes can also be seen from the arguments of π R and v R in (3.87), where it is the multiple of the tree-level frequency M that determines the number of quanta. Even though this discussion is based on the simple two-scalar model in (2.1), it is easily generalised to many other models and interactions since all the microphysical processes are encoded by the two integral kernels π R and v R appearing in the condensate equation of motion (1.2). The microphysical interpretation of these coefficients in terms of cutting rules holds irrespectively of the particle content.
Naturally, the result for damping (6.1) we found is an approximate one, and that is the reason why only processes involving one and two condensate quanta contribute to it. The evolution of the system is non-linear, and will necessarily generate contributions where many more condensate quanta decay into two or three φ and χ particles, (e.g. processes such as φφφφ → χχ or φφφ → χχφ), that would lead to the appearance of higher dissipation coefficients associated with higher number of quanta participating in the microphysical process. One way to see how higher dissipation coefficients arise is to continue with the algorithm in Sec. 3 beyond leading order and work out subleading corrections. It should be clear after some inspection that higher harmonics are generated in the condensate's evolution and that this leads to non-local terms generating higher dissipation coefficients at subleading orders. In addition, were we to include higher loops in the effective action (2.21) we would also account for higher order processes where one or two zero mode quanta decay into four or more φ and χ quasi-particles, producing further corrections to dissipation coefficients. However, the contributions of higher processes, either as corrections to dissipation coefficients or by generating higher dissipation coefficients, are expected to be subleading compared to the leading order approximation (6.1), and likely negligible in most applications. It is nonetheless important to note that they are present in principle, because it is due to them that a self-interacting scalar oscillating in the mildly non-linear regime will relax to the minimum via dissipating into fluctuations.
The Markovianised equation (4.21), on the other hand, will not capture the higher dissipation coefficients. By construction, it contains only γ and σ coefficients inside of the damping coefficient in (4.22b), and solving it to subleading orders will never generate any other dissipation coefficients apart from the products of these two. Going past this would require generalising the main Markovianisation relation (4.2) to include subleading corrections, which is not clear how it should be accomplished. Nevertheless, Markovianised equation captures the behaviour of the system very well, and as long as extreme precision is not required they are more than adequate to describe the system, as long as we fix the ambiguities arising by requiring that it follows from an action principle.
Damping coefficient vs. dissipation rate. It is very important to note that the local damping rate Υ(t) given in (6.1) is not the same thing as the damping coefficient Γ(ϕ) in (6.2) above, that appears in the Markovianised equation, and is of key interest in deriving effective quantum kinetic equations that generalise the standard Boltzmann equations. The two are sometimes conflated since in the case of linear condensate evolution, that is most often studied, the two are just constants of the same dimension, and thus proportional to each other, Γ(ϕ) σ=0 = 2Υ(t) σ=0 = γ. However, in general the local damping rate is a property of the solution and therefore is a function of time, while the damping coefficient is a property of the Markovian equation of motion and therefore is a function of the condensate. Whenever a diagram that contributes to dissipation descends from an interaction term containing more than one condensate field it is crucial to draw this distinction. 6 While solving the non-local equation directly and solving the Markovianised equation both produce the same local damping rate (6.1) and the same local frequency (3.98), the time-dependent perturbation theory we examined in Sec. 5 cannot describe the evolution of the condensate beyond very early times where damping is almost negligible and only linear in time. This is evident from the fact that the method produces a constant local damping rate, which is a zero temperature and early time limit of (6.1). In addition to not incorporating thermal statistical effects (except if formulated in terms of occupation numbers), this method also cannot account for feedback effects and damping. Both shortcomings make the application of time-dependent perturbation theory for the computation of transport coefficients in quantum kinetic equations questionable. Nevertheless, this method is sometimes used to work out the damping coefficients in Markovianised equations. While this works in cases where condensate evolution is linear (most often cubic interactions Φχ 2 are considered), in non-linear cases such as one considered in Sec. 5 the heuristic prescription Γ(ϕ) = 2Υ ϕ 0 →ϕ = σϕ 2 definitively fails, as it is off by a factor of 4 compared to the correct result (6.2).

Limits of applicability.
Having discussed the physical interpretation of our results, and different methods and approaches of obtaining them, it is important to discuss the limitations of the results -the restriction to the mildly non-linear regime, and perturbative processes, implemented via the small field expansion introduced in Sec. 2. By construction, this expansion neglects the non-perturbative effect of resonant particle production [9,10,12], which puts a restriction on the physical systems we can describe. Even though when present the resonant particle production typically is the dominant mechanism of dissipation in oscillating systems, there are potentially interesting regimes where this effect is subdominant, and where dissipation is dominated by perturbative processes. Such regimes are well described by the formalism outlined in Sec. 2 -the small field expansion within the 2PI formulation of the Schwinger-Keldysh formalism for non-equilibrium quantum field theory. The conditions restricting a given system to such a regime are illustrated on the example of the simple two-scalar model (2.1) we considered here, but are generalised to more complex models in a straightforward manner. There are two similar, but in principle distinct conditions which determine the applicability of our description. To understand both it is useful to note that the full effective masses, M 2 φ,χ = M 2 φ,χ +∆m 2 φ,χ , of φ and χ quasi-particles do not consist solely of constant tree-level and thermal contributions M 2 φ,χ from (2.18), but also receive time-dependent contributions ∆m 2 φ,χ on the account of coupling to the oscillating condensate ϕ(t) at the linear level. These can be read off from the inverse tree-level propagators in (2.3), The small field expansion used throughout this paper assumes these time-dependent contributions are much smaller than the respective tree-level and thermal masses, These restrictions justify the solutions for two-point functions in (2.19) and (2.20) being given as a small field expansion, which is controlled by the small parameters in (6.5) or (6.6). Since for simplicity we had presented the formalism in Sec. 2 in the high temperature limit, where thermal masses dominate over the tree-level ones, we quote the restrictions in the same limit of the two-scalar model, where the small field conditions above translate into for the particular two-scalar model, and are expressed as restrictions on the amplitude of oscillations in units of the temperature. The range of values that the initial condensate amplitude ϕ 0 is allowed to take crucially depends on the coupling constants and possible hierarchy between them.
In addition to assuming the small field expansion, we are neglecting the nonperturbative effect of resonant particle production, which is absent provided the effective quasi-particle masses change only adiabatically, For the simple two-scalar model this implies where we also made use of the already assumed small field restriction (6.5). Note that the left condition above coming from φ quasi-particles is the same as the respective condition from the small field expansion (6.6), while the one coming from χ quasiparticles differs. The reason behind this is that the frequency at which the effective masses oscillate are set by the thermal mass of the condensate, which equals the thermal mass of φ quasi-particles, but differs in general from the one of χ quasiparticles.
It is interesting to ask whether the transition between damping regimes seen in Fig. 6, which requires σϕ 2 0 /γ ≫ 1, could be supported by the simple two-scalar model (2.1). In answering this question we are limited by the knowledge of twoloop self-energy from (2.23). We only have an estimate for the part coming from self-interactions governed by the coupling constant λ φ , in the upper case of (2.34), and that only in the high temperature limit. That is the dominant contribution to self-energy if we assume a hierarchy λ φ ≫ g. In that case conditions (6.6) and (6.8) imply ϕ 2 0 /T 2 ≪ g/λ φ , and the relevant ratio is σϕ 2 0 /γ ∼ [ϕ 2 0 /T 2 ]×[g 2 /λ 2 φ ]×λ −1/2 φ , which could possibly be larger than one if the coupling λ φ is small enough to overcome the other two small ratios, confirming the relevance of our result (3.92) even in this simple model. On the other hand, assuming a different hierarchy λ χ ≫ g ≫ λ φ would allow the for ϕ 2 0 /T 2 1, and, assuming that the self-energy contribution governed by the coupling constant g would take the same form as the one in (2.34), we get for the relevant ratio to be of the order σϕ 2 0 /γ ∼ [ϕ 2 0 /T 2 ]×λ −1/2 φ , which could indeed be very large. However, the kinematic condition from (2.32) in this case is not met, and because of it σ dissipation coefficient actually vanishes. Nevertheless, assuming this hierarchy for a system in vacuum (at vanishing temperature) might lead to different conclusions, as the restrictions (6.6) and (6.8) would be modified and the kinematic condition in (2.33) could be met by tree-level masses, though some of the approximations we make to derive the results might not apply in that case [78]. More complex models should allow for interesting regimes more readily. Our work provides a description of the evolution of such systems with an oscillating scalar condensate, once the appropriate self-energies and proper four-vertices are computed in a given particle physics model. Detailed studies of phenomenological applications of the results presented here are left for future work.

Conclusion.
Non-linear effects in the evolution of an scalar condensate oscillating in a medium can be relevant even in the regime where non-perturbative resonant particle production is negligible. This regime, which we refer to as the mildly nonlinear regime, has not been investigated nearly as much as the linear regime or the non-perturbative regime, and this work helps to fill this gap (see also [58] for a recent study). The non-equilibrium and the non-linear nature of the problem necessitates a judicious choice of formalism and approximation methods for describing the evolution. The 2PI formalism in the Schwinger-Keldysh formulation we utilize is appropriate for describing quantum and statistical effects, and yields an effective non-local equation of motion for the scalar condensate (1.2). This macroscopic equation is connected to elementary microscopic processes via the self-energy and the proper four-vertex calculable in a given particle physics model. We examined solutions of Eq. (1.2) and methods of obtaining them. The most reliable, though most difficult, approach is solving the non-local equation directly. This we accomplished using multiple-scale perturbation theory, appropriate for time-dependent problems, to find the effects of non-linear evolution in solution (3.92) -power-law damping (6.1) and time-dependent frequency shifts -in addition to the well known effects of exponential damping and constant frequency shifts known from linear evolution. A considerably simpler and far more used approach to solving equations such as (1.2) is approximating them by Markovian counterparts. We found this approach to describe the system at leading order equally well as solving the non-local equation directly, though some care is necessary when attempting more precise descriptions. Standard time-dependent perturbation theory fails to describe systems for long times, and is able to describe the oscillating condensate only during early times while its damping is negligible. Our results provide a guideline for the appropriate choice of methods when describing dissipative effects and particle production in the early Universe, e.g. during reheating, Dark Matter production or moduli decay. Since they were obtained with minimal asumptions, they may be applied to a wide range of other systems that can effectively be described by an oscillating scalar field.