Effective Action for Cosmological Scalar Fields at Finite Temperature

Scalar fields appear in many theories beyond the Standard Model of particle physics. In the early universe, they are exposed to extreme conditions, including high temperature and rapid cosmic expansion. Understanding their behavior in this environment is crucial to understand the implications for cosmology. We calculate the finite temperature effective action for the field expectation value in two particularly important cases, for damped oscillations near the ground state and for scalar fields with a flat potential. We find that the behavior in both cases can in good approximation be described by a complex valued effective potential that yields Markovian equations of motion. Near the potential minimum, we recover the solution to the well-known Langevin equation. For large field values we find a very different behavior, and our result for the damping coefficient differs from the expressions frequently used in the literature. We illustrate our results in a simple scalar model, for which we give analytic approximations for the effective potential and damping coefficient. We also provide various expressions for loop integrals at finite temperature that are useful for future calculations in other models.

In cosmology involving scalar fields, it is crucial to understand their time evolution in the early universe, where they are exposed to extreme conditions including high temperature, large energy density and rapid cosmic expansion. Their evolution in a time-dependent background provided by the primordial plasma and the cosmic expansion is a nonequilibrium process. The initial conditions for this process depend on the model with the specific scalar fields, and the mechanism that set the initial conditions is unknown in almost all models. Many observational hints point towards cosmic inflation, but alternative models such as bouncing universe scenarios (see [21,22] for review) have been also studied, though there are still controversies [23,24]. The bouncing models also invoke scalar fields, typically inspired by string theory, to resolve the cosmic singularity by a bouncing phase in the very early universe [25,26]. In general, the field value at initial time is far away from the potential minimum. If the effective potential is sufficiently steep, the field will relax to its minimum very quickly. However, there are numerous examples for scalar fields with a flat potential that evolve slowly compared to the time scale related to the propagation and interactions of individual particles, including the moduli, inflaton and axions. With the present work, focusing on the finite-temperature effects in a thermal bath, we aim to make progress towards a quantitative understanding of the nonequilibrium dynamics of scalar fields in the nontrivial backgrounds of the early universe.

Motivations and assumptions
Let us consider a real scalar quantum field φ with a symmetry under φ → −φ that couples to other degrees of freedom, which we collectively refer to as X i . At this stage we do not make any assumptions about the spin or interactions of the fields X i . If we restrict ourselves to renormalizable operators, then the most general Lagrangian reads with Here O[X i ] is a sum of operators that depend on combinations of the fields X i only and L X is a Lagrangian that specifies their masses and interactions amongst each other. The generalization to non-renormalizable operators is straightforward. If we neglect the interactions with X i for a moment, then the zero mode of φ in a Friedmann-Robertson-Walker spacetime would classically follow the equation of motion where H is the Hubble constant. The full quantum field φ can be decomposed into Here ϕ is the expectation value of the field φ, while the field η contains quantum fluctuations, such as single particle excitations. If we switch on the interactions with the background medium and take into account radiative corrections, one could expect ϕ to follow an equation of motion of the form Here V(ϕ) is an effective potential for ϕ that includes radiative and thermal corrections, while Γ ϕ is the rate at which energy is dissipated from ϕ to other degrees of freedom. The equation (5) looks physically intuitive, but should be rigorously justified when being applied. It cannot always be valid because the nonequilibrium dynamics of interacting quantum fields involves non-Markovian memory kernels. It is the goal of this work to derive an effective equation of motion of the form (5) from first principles, and to determine V(ϕ) and Γ ϕ explicitly in an illustrative simple model. Together with the effective kinetic equation for the η-propagators derived in [27], this allows to understand the dynamics of the scalar field in terms of Markovian equations. Our derivation of (5) uses a somewhat different method than previous approaches in [28][29][30][31][32][33][34], but yields a result that is consistent with previous works in the limit of small field values, where a comparison can be made. The explicit expression for Γ ϕ at large field values, however, is considerably different from those found in the literature. In order to obtain a Markovian equation of motion for ϕ, we have to make several approximations. As usual in the derivation of effective kinetic equations [35], these rely on a separation of time scales. In closed systems, such a separation is automatically given if one assumes a sufficiently weak coupling. In this case the rates at which particles scatter are much smaller than their effective masses. This implies that macroscopic properties of the system change much slower than the time scale associated with individual scatterings. It also ensures that the system can be described as a dilute gas of well-defined quasiparticles with effective masses M η , M X i and widths Γ η , Γ X i with Γ η M η , Γ X i M X i . Here M η is the pole mass of the resummed η-propagator (the same for M X i ) and in general effective mass and width both depend on momentum 1 . 1 In thermal equilibrium the dispersion relations and widths of quasiparticles can be read off from the pole structure of the real time propagators [36]. This definition can be generalized to situations far from equilibrium, see e.g. [27]. The propagators in a medium can in principle have a complicated pole structure due to the appearance of collective excitations, such as luons [37]. We ignore this subtlety here.
In the early universe, Hubble expansion brings in an additional time scale. In the present work we make two basic assumptions 1) The effective masses of all particles that couple to ϕ are larger than the rate of Hubble expansion, i.e. H < M η , M X i .
2) The interactions in the primordial plasma are sufficiently weak to guarantee Γ η , Γ ϕ M η , M ϕ (and analogous for all other fields X i ). Here M 2 ϕ is the local curvature of V(ϕ). For practical purposes we make two additional assumptions that are not required for the derivation of a Markovian equation of motion, but lead to considerable simplifications 3) φ is the only field with a non-zero expectation value, i.e. X i = 0.
4) The thermodynamical state of the constituents X i of the primordial plasma can be characterized by a single parameter, an effective temperature T . This is usually justified if Γ X i > H.
Assumption 1) allows to use Minkowski-space propagators in loop diagrams. Although H is absolutely crucial in the kinetic equation for ϕ on macroscopic time scales, it can be neglected when determining the transport coefficients in this equation from microphysics. Without this assumption the computation of loops is technically very difficult [38][39][40][41].
Assumption 2) is simply necessary for any perturbative treatment to be valid. It also ensures that the effective masses (which include contributions that depend on ϕ and T , which evolve with time) are approximately constant on the microscopic time scale that is relevant in loop integrals. Assumption 3) is for simplicity only, and it is straightforward to relax it. Assumption 4) implies that X i have reached some degree of kinetic equilibration. Though thermalization in general is a complicated problem [42][43][44][45][46][47][48][49][50][51][52], this assumption seems at least qualitatively justified in the case under investigation here. If φ has a flat effective potential, then this means that it must have rather weak interactions (weaker than the interactions of the X i amongst each other). This implies that the time scales 1/Γ ϕ , 1/H on which ϕ evolves are much longer than the time scales 1/Γ X i on which the X i relax to local thermal equilibrium. Making the same assumption for η is far less justified. For simplicity, we will nevertheless assume that the occupation numbers of η modes can be characterized by the same effective temperature as the rest of the plasma. This assumption, which has been made in almost all past studies (often without mentioning), allows to use equilibrium propagators for all fields. The setup defined by the assumptions 1)-4) is not sufficient for a complete understanding of scalar fields in the early universe. In particular assumptions 1) and 4) considerably constrain its applicability. Moreover, the properties of the primordial plasma in realistic models of particle physics is far more complicated than the toy model presented in section 4. However, the analytic expressions we find are derived more systematically and consistently than any comparable results in the literature we are aware of. Moreover, in the course of their derivation we find analytic estimates for various integrals in finite-temperature field theory that will be very useful for calculations in more realistic models. Finally, in spite of the comments above, there are at least two cosmological problems in which the assumptions 1)-4) are widely believed to be justified: warm inflation [53] and the fate of moduli at late times [54].

Overview of this article
This paper is organized as follows. In section 2 we briefly review the elements of nonequilibrium quantum field theory to set up the theoretical framework and notations for the following sections. In section 3 we derive the approximate effective action and the equation of motion for the expectation value ϕ in a thermal bath. These results are illustrated in a simple scalar model in section 4, and we present approximate analytic estimates for effective potentials and dissipation coefficients. In section 5 we summarize the results of this work and discuss some related issues in comparison with previous studies. Section 6 is devoted to conclusions and outlook. In the appendices we calculate the spectral selfenergies from setting-sun diagrams, which are needed for the computation of the dissipation coefficients in section 4.

Elements of nonequilibrium quantum field theory
The standard methods to calculate S-matrix elements in particle physics are not suitable to describe systems far from equilibrium at large density because there is no well-defined notion of asymptotic states, the properties of elementary excitations may significantly differ from those of particles in vacuum and classical particle number in general is not a suitable quantity to characterize the system. However, all observables can be expressed in terms of time-dependent correlation functions of the quantum fields without reference to asymptotic states or free particles. In most practical applications, all relevant information is contained in the one-and two-point functions. The expectation value ϕ can be identified with the one-point function, That is, the average . . . = Tr( . . .) is taken over quantum as well as statistical fluctuations encoded in the density operator . We will in the following always assume that ϕ = 0 in the thermodynamic ground state and study how the system relaxes to this state if ϕ = 0 at initial time. There are two independent two-point functions. Common choices for these are the connected Wightmanfunctions or their linear combinations These have a clear physical interpretation. The Fourier transform of ∆ − is the spectral density ρ η , the poles of which in p 0 determine the spectrum of quasiparticles in a medium. The Fourier transform of ∆ + , on the other hand, characterizes the occupation numbers of different field modes and can be interpreted as a field theoretical generalization of the classical particle distribution function, see e.g. discussion in [27]. In most practical applications, ϕ, ∆ + and ∆ − contain all relevant informations.
In this formalism, correlation functions are defined with time arguments on the contour shown in figure 1, which starts from an initial time t i + i , runs parallel to the real axis to t f + i , turns around to t f − i and returns to t i − i . We consider the limit t f → ∞, → 0 and t i → −∞, where boundary conditions at finite time t 0 can be imposed by external sources localized at t 0 [27]. In the following we define the quantities required for the present analysis. A more detailed introduction to the CTP-formalism can be found in [59,60].
The standard strategy to calculate correlation functions is to split the field φ on the contour into φ + and φ − , which denote the field with time argument on the "forward" and "backward" part of C, see e.g. [60]. One can define a Feynman propagator T C φ(x 1 )φ(x 2 ) , where T C indicates time ordering along the contour path C. It can be decomposed into the correlation functions for these fields as where T is the usual time-ordering andT is anti-time-ordering (because the backwards part of the contour runs from positive to negative infinity). Considering that time arguments on the forward branch are always "earlier" along the contour than those on the backward branch, one can easily identify The same decomposition can be performed for the self-energies, see appendix of [61]. Then the combination φ c ≡ (φ + + φ − )/2 is identified with a physical field and φ ∆ ≡ φ + − φ − is treated as a response field. After integrating out the bath fields X i and φ ∆ , one obtains an effective action for φ c , from which one can obtain its equation of motion. It can be expanded in powers of φ c . At lowest non-trivial order one finds the Langevin type equation (45), including higher powers of φ c in the effective action leads to nonlinear interactions between the different field modes ("multiplicative noise"), see e.g. [30,33]. By taking quantum and thermal averages over products of φ c , one can then calculate all correlation functions for φ.

Equation of motion for the background field ϕ
In the present context we are, however, mainly interested in the one-point function ϕ. We can therefore simplify the procedure and directly write down an effective action for ϕ (instead of φ) on the closed time path C [62][63][64][65], Here the label C indicates that the time integral is to be taken along C. The Γ (n) (x 1 , . . . , x n ) are the amputated vertex functionals. The expansion in (17) is around the minimum at ϕ = 0, and the coefficients Γ (n) are taken at ϕ = 0 in the sense that they contain the field η defined in (4). However, for a consistent expansion in the small couplings requires the use of full (resumed) propagators in the loops because the T -and ϕ-dependent corrections to effective masses bring powers of the couplings into the denominator of the propagators (similar to the situation of gauge theories at high T [66]). This is done in the 2PI formalism of non-equilibrium quantum field theory, see e.g. [59]. Hence, the Γ (n) implicitly depend on ϕ via the effective masses. With (2) in mind, we expand (17) up to fourth power. Γ (1) is eliminated by the stationarity condition for vanishing external sources, Γ (2) is simply the inverse dressed propagator −iΓ (2) (3) vanishes due to the symmetry φ ↔ −φ in (2) and Γ (4) is the amputated four-point function. Due to energy-momentum conservation we can replace i.e. express the effective action in terms of an effective potential V[ϕ] ≡ −δ (4) (0)Γ [ϕ]. The functionΠ is, for example, given by the diagram in figure 5. We now perform a spatial Fourier transform to the mixed representation with arguments (t, p) and decompose ϕ into modes p, which in general are coupled to each other by the integrals in (17). We are interested in the case, in which only p = 0 mode of ϕ significantly differs from zero, and we can thus perform the spatial integrals analytically. Using the stationarity condition (18) we can find an equation of motion for ϕ on C, Here we have suppressed the momentum index and all functions are to be evaluated at p = 0. Note that this equation does not contain a Brownian noise term ξ on the RHS, as one might have expected. The reason is that this term comes from thermal fluctuations, and the expectation (7) that is taken in the definition (6) of ϕ includes an average over thermal fluctuations. A noise therm ξ does appear in the (Langevin type) equation of motion (45) for φ if one only traces over the bath degrees of freedom [28][29][30][31][32][33]. While these fluctuations cancel out in the equation for the one-point function ϕ, they crucially affect n-point correlation functions of φ with n > 1.

Equations of motion for the fluctuations η and perturbation theory
In order to calculate the coefficients Π andΠ in (20) one has to evaluate loop integrals with η and X i running in the loop. This requires knowledge of the propagators. Their equations of motion are the Kadanoff-Baym equations (see [67]), which can be derived in a similar manner as (20), see e.g. [32], where m η is a vacuum mass of η-quanta 2 and the self-energies are defined, equivalently to (10), as The usual retarded self-energy can be identified with Approximate Markovian equations -Correlation functions and self-energies for the bath fields X i can be defined analogously. In general, the coupled second order integrodifferential equations (21) and (22) are difficult to solve. However, in a weakly coupled theory, and in a background that changes slowly compared to the time scale of individual particle scatterings, the system can be described as a gas of quasiparticles. In a homogeneous medium, the spectral density defined in (11) only depends on time, ρ η (p; x) = ρ η (p; t).
LetΩ p (t) be a pole of ρ η (p; t), with Both are in general time-dependent due to the interaction with the medium. In weakly coupled theories one observes the hierarchy In that case ρ η (p; t) at fixed t features peaks of width ∼ Γ η (t) at energies p 0 ±Ω η (t), which can be interpreted as quasiparticle-resonances,and loop integrals are typically dominated by the regions near these poles. Here (and in the following) we omit the spatial momentum index. In this situation, one approximately finds [27] Here f (t) is a generalized distribution functions, which follows the Markovian equation of motion with the "gain" and "loss" terms 3 see also [68]. 2 Note that m η is in general different from m φ since it obtains ϕ-dependence due to self-interaction of φ: For example, m 2 3 Here we ignore subtleties about finite initial time in definitions.
Special case -thermal equilibrium -If the background medium is in thermal equilibrium, the functions ∆ − η and Π − η are independent of the coordinate x 1 + x 2 [32]. Moreover, the Fourier-transforms of the self-energies in the coordinate x 1 − x 2 are related by the Kubo-Martin-Schwinger (KMS) relation or equivalently where f B (p 0 ) = (e p 0 /T − 1) −1 is the Bose-Einstein distribution. In this case, ρ η reads [32] where is an infinitesimal parameter. Note that the self-energies Ω η and Γ η all depend on T . It is clear from (33) that, if condition (26) is fulfilled, the quasiparticle dispersion relation (or "mass shell") Ω η is essentially fixed by ReΠ R η (p) via the condition while ImΠ R η (p) gives the thermal width, which is just the difference of gain and loss term. ReΠ R η (p) contains a zero temperature divergence that can be absorbed into the physical mass in the same way as in vacuum [31,32]. In the following we interpret m φ as the physical mass in vacuum, with ReΠ R η (p) being finite. The effective mass M η can be defined as Ω η for vanishing spacial momentum p = 0. It is very useful to fit a Breit-Wigner function to (33) near the pole, which is parametrized by Γ η and Ω η , as Here ρ cont η (p) is the continuous part of ρ η (p). To obtain the correct residue, we introduced the parameter For weak coupling, it is often sufficient to use the "zero width limit", Other useful relations that follow directly from the definitions of the correlators are and equivalently for the self-energies of other fields. Note that all the above functions depend on T , which we have not made explicit here to simplify the notation. Correlators and self-energies for all other fields are defined analogously.
Perturbation theory -The perturbative expansion for correlation functions in a medium can be performed in terms of Feynman diagrams, analogous to the vacuum case. One difference is that one has to use the propagators ∆ ≷ or, often more conveniently, their transforms in Wigner space [68]. Knowledge of ∆ ≷ allows to determine all other correlation functions, in particular Another difference from the vacuum case is the fact that φ + and φ − formally have to be treated as two different fields. There are no vertices that directly connect a φ + line to a φ − line, i.e. the vertices are always either of "+" type or "-" type, but the propagators (43) and (44) mix the fields and can connect a "-" vertex to a "+" vertex and vice versa, see e.g. [36] for a detailed discussion.

Effective action and equation of motion for ϕ
The known equations (27), (28) and (41)- (44) allow to describe the dynamics of the field fluctuations in terms of the Markovian equations, which can be physically interpreted as effective Boltzmann equations for quasiparticles in the plasma. In the following sections we derive an approximate Markovian equation for the field ϕ, so that the system can be completely described in terms of Markovian dynamics.

Small field values: Brownian motion and Langevin dynamics
We consider the case in which all degrees of freedom except ϕ are in thermal equilibrium.
If φ at initial time is very close to its thermodynamic ground state, then the effect that this deviation has on the thermal bath ("back reaction") is negligible. This in particular requires that ϕ be small enough that we can use a quadratic approximation to V(ϕ), and the effective masses in the plasma are dominated by vacuum and thermal masses (rather than coupling to the background ϕ). Then the self-energies have the properties (32). In this case the behavior of φ can be understood in the terminology of Brownian motion. Such systems have been studied by a number of authors [28][29][30][31][32][33][69][70][71], and it was found that the zero mode in this case follow a Langevin equation [31,32], Here the integral leads to dissipation, and ξ is a Gaussian noise term with The relations (32) and (46) ensure that the fluctuation-dissipation theorem holds in (45). Strictly speaking the field in (45) is not the original field φ on the contour C, but the combination φ c = (φ + + φ − )/2, where φ + and φ − denote the field with time argument on the "forward" and "backward" part of the contour, as indicated in figure 1. We ignore this subtlety here, as both notions give the same equation of motion for ϕ. The Langevin equation (45) can be solved analytically by use of a Laplace transformation [31,32], where can be analytically obtained by Fourier-transforming (33) in the Breit-Wigner approximation as By taking the expectation value (cf. (6)) of φ in (47), we find (neglecting the p 0 -dependence of ReΠ R η , setting Ω η = M η for the zero mode and neglecting terms suppressed by Γ η /Ω η ) The solution (50) describes damped oscillations of ϕ with the plasma frequency M η near the potential minimum. In spite of the integral in (45), non-Markovian effects do not play any role for ϕ because Γ η and M η are in good approximation independent of ϕ, and we have assumed that the temperature remains constant on time scales ∼ 1/Γ η . This is realistic if Γ η > H, and if the initial ϕ is sufficiently small that the energy release into the plasma due to its dissipation is not significant. In this case an equation of the form (5) holds, with That is, V(ϕ) can be approximated by a parabola and is obtained from the free classical Lagrangian if one simply replaces the mass m φ by the effective mass M η . The damping rate of ϕ is given by the thermal quasiparticle width of η; it is obtained by evaluating ImΠ R η (p) at the quasiparticle pole p 0 = M η . Finally, each mode of φ couples to the bath, but there are no interactions of the different modes with each other. This approximation only holds for small deviations from the ground state, i.e. small field values ϕ. In cosmology this situation can be realized in the late phase of cosmic reheating after inflation. In this context, the temperature dependence of the rate Γ ϕ has been studied in detail in [72], and it was found in [73] that this temperature dependence can have a drastic effect on the thermal history of the universe during reheating.
For larger deviations from the ground state, the problem becomes much more involved because the quadratic approximation of the effective action near ϕ = 0 is not sufficient. If it is still possible to describe the dynamics in terms of an effective potential V(ϕ) and (5), then V(ϕ) certainly must involve higher powers of ϕ. This leads to "multiplicative noise terms" in the generalization of (45) [30,33], and the effective masses generally depend on ϕ. Finally, we do not expect the simple relation (52) to hold, though it is frequently used in the literature without justification. In section 3.2 we show that one can nevertheless find an equation of motion of the form (5) if ϕ changes slowly compared to the microphysical scales, but the expressions for V(ϕ) and Γ ϕ are in general more complicated.

Large field values: Slow rolling
We continue to assume that all degrees of freedom except ϕ are in equilibrium (or their thermodynamic state can at least effectively be characterized by an effective temperature), but we allow ϕ to take large values, where the quadratic approximation (51) of the effective potential is expected to break down. If φ has only feeble interactions (which is required if the effective potential should be flat), then it is justified to assume a separation of scales This implies that ϕ changes adiabatically, and we can approximate The change of the functions Π C (t, t ) andΠ C (t, t ) with respect to the time coordinate t + t occurs on the same macroscopic time scale asφ/ϕ, and we can expand them in this coordinate analogously to the expansion (54). For the present purpose it is sufficient to keep the zeroth term in this expansion, i.e. to assume that there is no explicit dependence . This will be justified later after (57). We can split the contour integral in the equation of motion (20) into With (54), the only quantities under the integrals that depend on t are Π C (t − t ) and Π C (t − t ). Since we are interested in physical times, t always lies on the "forward" part of the contour C. If t also lies on the forward part (first term in (55)), then Π C (t − t ) is to be identified with the usual time-ordered self-energy (often referred to as Π ++ (t − t ) in thermal field theory). If t lies on the backward part (second term in (55)), then Π C (t − t ) is to be identified with a "Wightman type" self-energy without time ordering Π < (t − t ), as t is then always "later" along the contour than t, see e.g. [36,60]. Their difference can be identified with the usual retarded self-energy, i.e. Π This allows to combine the two terms in (55) into a single integral and likewise forΠ C (t − t ). The RHS of (56) is nothing but the Fourier transform of Π R (t − t ) evaluated at energy ω = 0. Some care should be taken at this point. The limit ω → 0 is obtained due to the linear approximation (54), which is controlled by the separation (53) between macroscopic and microscopic time scales. It is probably best to think of ω as a parameter that characterizes this separation, i.e. ω ∼φ/ϕ, whereφ/ϕ should be identified with the macroscopic scales H or Γ. With these considerations we finally obtain a Markovian approximation to the equation of motion (20), Here we have given the retarded self-energies Π R ϕ andΠ R ϕ an index ϕ to distinguish them for the self-energies for X i . The parameter ω should be thought of as much smaller than all microphysical scales and can be set to zero as long as ω < H < Γ η < T , which always holds under the previous assumptions. Since the real part of the retarded self-energy is always symmetric and the imaginary part antisymmetric, we can easily identify and likewise forΠ R ϕ . With this limit and by comparing (57) with (5), we can express the effective potential and dissipation coefficient in terms of Π R ϕ andΠ R ϕ . The first line in (57) describes the evolution of ϕ in an effective potential with the quantum corrections to the classical potential given by Π R ϕ andΠ R ϕ : The second line describes dissipative effects, i.e. the damping of ϕ due to the interactions with the bath of X i -particles: These results allow to justify why we could neglect their dependence on t + t . The terms containing derivatives ∂ t+t acting on Π R ϕ andΠ R ϕ are suppressed with respect to terms without derivatives due to the separation of macroscopic and microscopic scales. For instance, ∂ t+t Π R ϕ (t, t ) is suppressed by three small parameters, two powers of a coupling constant and a time derivative. All terms in the first line of (57), which give the effective potential (60), contain two or less small parameters, hence ∂ t+t Π R ϕ (t, t ) can be regarded as a higher order correction. In spite of this, we have to keep the terms in the second line in (57), which in principle are also suppressed by a time derivative and two powers of a coupling constant. The reason is that these form the leading order contribution to the dissipation coefficient (61), while in the effective potential (60) they would just act as higher order correction. In this sense, (57) is a consistent expansion in gradients and the small coupling constants.
These results (60) and (61) should be compared to (51) and (52). One obvious (and expected) difference is the appearance of higher powers of ϕ in V(ϕ). Another crucial difference is that the dissipation coefficient (61) is evaluated at ω = 0 and not at ω = M η , as in (52). This supports the interpretation of ω as the rate at which ϕ changes: During oscillations around the minimum, this rate is given by the plasma frequency M η ; if the field slowly rolls down a flat potential, then this rate is much smaller than all other scales (practically ω 0 within loop diagrams).
Indeed, if we drop the terms ∝ ϕ 3 in (20) and use (55) and (56), we find This equation is formally the same as (45) with ξ = 0, so we can read off the solution from (47), Again using the Breit-Wigner approximation for (33), we can bring this into the explicit form (50). Hence, we have recovered the behavior of Brownian motion for small ϕ from the more general expression (20).

A simple scalar model
In the following we illustrate our results in a simple scalar model, in which the "bath" consists only of one other real scalar field χ,  In order to obtain explicit expressions for V(ϕ) and Γ ϕ from (60) and (61), we have to evaluate the loop diagrams given in figure 2, 3 and 5.
Here we focus on the case of adiabatically changing large ϕ discussed in section 3.2. For the case of oscillations near the potential minimum, Γ ϕ has been calculated in detail in [72].

The effective potential V(ϕ)
From (58) it is clear that V(ϕ) is obtained from the real part of retarded self-energies. To establish a consistent perturbation theory, these have to be evaluated with resummed thermal propagators. For the present purpose it is sufficient to evaluate these in the pole approximation of (38) for η and χ. The potential in (63) in terms of η and χ reads The real part of Π R ϕ | ω=0 is dominated by the well-known local diagrams due to the vertices λ φ η 4 /4! and hη 2 χ 2 /4, see figure 2. It reads with where we have used dimensional regularization. The divergent part (N terms) can be eliminated by renormalisation of m φ and coupling constants at zero temperature, using MS scheme. 4 In the following we assume that this has been done and that vacuum masses and coupling constants are the renormalized parameters at T = 0. In this work we are mainly interested in thermal corrections, which are most relevant if the temperature is larger than the effective masses (otherwise the thermal corrections are negligible since they are Boltzmann suppressed due to the appearance of the distribution functions f B in the thermal propagators). For M η , M χ T , the finite part (i.e. the T(M 2 , µ 2 ) terms), remaining after subtraction of the divergences, can be ignored since it is subdominant compared to the f B -part (the third line of (65)). 5 Due to the momentum-independence there is no need for a finite T wave function renormalisation. For M η , M χ T , one can perform the integrals analytically, for example, The same applies to the diagram with η-loop. We thus obtain for M η , M χ T There are also contributions from the setting sun diagrams in figure 3, but they are of higher order in h and λ φ .

The vertex functionalΠ R ϕ
We now turn to the vertex functional, which is given by the "fish diagrams" a) in figure  5. The f B -independent part of these diagrams is again divergent; for vanishing external momentum it is given by (including the analogous fish diagram with η-loop) Since we are only interested in the case of vanishing external momentum, we can simply absorb the N terms into the renormalisation of the vacuum couplings h and λ φ . The remaining finite part (log terms) can be ignored for M η , M χ T in comparison with f B -part given below, in the similar way to the case of Π R ϕ . 6 The f B -part from the fish diagrams is given by the the standard expression (see e.g. section 4.2.2 in [36]). Let us focus on the χ-loop with T M χ , The above integral can not be done analytically, but it is easy to see that it is strongly dominated by small Ω when M η , M χ T . Therefore we can approximate f B (Ω) T /Ω by expanding the exponential to get The contribution from the η-loop can be obtained from the above expression by replacing χ → η and h → λ φ , and altogether we obtain for M η , M χ T

The full potential
If the temperature is much smaller than the effective masses, then thermal corrections are negligible. In this case radiative corrections are mainly interesting if some mass squares are negative, leading to symmetry breaking; otherwise they just lead to small modifications of the potential. If T is comparable to the effective masses, then the loop integral in general cannot be solved analytically. For the interesting case M η , M χ T , we can use the above results in (60) and find In order to check for which range of field values the assumption M η , M χ T can be justified, we can estimate the effective masses as This leads to the conditions These conditions imply that the ϕ-dependence of the effective masses in (73) is not necessarily negligible. Taking this into account, we integrate (73) and obtain where Figure 4 shows the effective potential (77) as a function of ϕ for various temperatures.

The damping coefficient Γ ϕ
Damping coefficient (61) depends on ϕ through the effective masses and has to be evaluated at ω = 0.

Contribution from the vertexΠ R ϕ
Let us focus on diagram a) in figure 5, which gives dominant contribution 8 . The contribution from the analogous diagram with η-loop can be obtained by proper replacements of coupling constants and masses. In the zero-width approximation (38), the imaginary part ImΠ R ϕ vanishes for any finite positive ω < 2M χ and diverges for ω = 0. 9 One can understand this from kinematic reasons when recalling that the imaginary part can be interpreted in terms of microphysical processes when applying the optical theorem at finite T [72], see figure 8. A non-zero ∂ ωΠ R ϕ | ω→0 can be obtained by including the finite width of the χ-propagators in the loop, which parameterizes the effect of scattering processes with more vertices [61,72], see figure 9. This aspect has been pointed out in [74,75], but we find different results from the ones given there, as discussed in section 5.2.2. The loop in the "fish diagram" a) in figure 5 can be expressed as [61] ImΠ For notational simplicity we continue to suppress the dependence on spatial momentum; the external spatial momentum is vanishing, and for all quantities under the integral it is given by the loop momentum. Here we have first used (39) and the KMS relation (32) to express ImΠ R ϕ in terms ofΠ < ϕ . Using the Feynman rules sketched after equation (44), which are derived in detail in [36], this function can be calculated from an integral over ∆ ≷ χ alone. With (40), this can finally be expressed as an integral over a product of spectral densities, which is particularly easy to be evaluated in the approximation (38). In the analogous expressions (93) and (122) for the setting sun diagrams one can indeed use the approximation (38). In (81), on the other hand, we cannot neglect the finite width and have to use the pole approximation to (33), see e.g. [72] for a discussion, For ω Γ χ the quasiparticle resonances at p 0 Ω χ and p 0 Ω χ −ω of (82) in (81) overlap in the entire integration region, and the integral is strongly dominated order in this case.
8 See section 4.2.4 for the argument that the vertex diagram b) in figure 5 is subdominant. 9 These two facts can be seen from figure 10 b) and c), which depict the supports of the spectral densities in the expresion (81) for ImΠ R ϕ , as discussed in section 5.2.2.
The calculation of the momentum-dependent thermal width Γ χ from the two loop diagram in figure 6 is generally difficult. The full Γ χ is sum of contributions from diagram a) and b) in figure 6, i.e. Γ χ = Γ The main momentum dependence is due to the time dilatation; the factor 1/Ω χ is due to the extended lifetime of a relativistic χ-particle in the bath rest frame. The contribution from the diagram b) in figure 6 should have essentially the same behavior as that of a) as long as M χ , M φ T because the loop is dominated by hard momenta M χ , M φ in this case. Therefore one can obtain the contribution Γ The integral then reads . This integral is strongly dominated by the region Ω χ T . We expand the integrand in Ω χ /T up to second order and integrate from Ω χ = M χ up to Ω χ = T to obtain 10 This is for the fish diagram a) in figure 5 with χ-loop. The contribution from the analogous fish diagram with η-loop can be obtained from above expression by the replacement M χ → M η and h 2 /(λ 2 χ + 3h 2 ) → λ 2 φ /(λ 2 φ + 3h 2 ).

The full damping coefficient
Using above results in (61), we obtain the full damping coefficient for which is plotted in figure 7 as a function of ϕ for various temperatures.

Physical interpretation of the coefficients
The parametric dependence of (80) can easily be understood physically; the rate for 2 → 2 scatterings is proportional to h 2 and grows ∝ T 2 due to the higher density of scattering 10 It should be pointed out that for 0 < ω < Γ χ and with the previous assumptions, the integrand in (83) formally exhibits a sharp peak at Ω χ ≈ γ 2/3 In this region the integrand can be approximated by the Breit-Wigner curve. In the zero width limit the contribution from this region can be estimated as Though this contribution smoothly vanishes in the limit ω → 0, it is much bigger than (90) for most finite 0 < ω < Γ χ even if ω Γ χ . This appears to be at odds with our intuitive argument that the result of (83) should be independent of ω if |ω| Γ χ . The contribution (89) is not physical and we artificially introduced by the approximation (84), which is only valid in the region Ω χ Γ 2 χ /ω.  Figure a) shows the full damping coefficient Γ ϕ . In figure b) we plot the ratio Γ 1 /Γ 2 with Γ 1 and Γ 2 being contributions from Π R ϕ (setting-sun diagram) andΠ R ϕ (fish diagram), respectively (i.e. the first and the second term in (91)). This plot shows that Γ 1 is dominant only for very small ϕ and otherwise Γ 2 is the main contribution to the damping. Γ 1 and Γ 2 are separately plotted in figure c) and d).
partners, just as a quasiparticle damping rate would. The coefficient (90) shows a more unusual behavior: most notably it decreases with temperature and with λ χ , both of which may seem counter-intuitive if one has the scattering interpretation in mind. It should, however, be pointed out that neither the limit T → 0 nor λ χ → 0 can be applied to the approximate analytic formula (90) because it is only valid under the assumptions M χ < T andφ/ϕ < Γ χ ∝ λ 2 χ , in which case the limit ω → 0 is only justified. In addition, the limit λ χ → 0, which makes Γ χ ∝ λ 2 χ vanish, is not allowed for following reasons. In the approximation Γ χ = 0, the cut through the diagram in figure 8 at finite T includes various processes, which are obviously kinematically impossible for any finite ω < 2M χ as explained in the caption. Moreover, for ω = 0 the rate has an unphysical divergence, which can be traced back to the infrared divergence of the Bose-Einstein distribution for the quasiparticles. 11 The divergence in the limit ω → 0 is regularized by the finite width Γ χ = 0 of the intermediate χ-particle in the contributing processes, see figure 9. Therefore Γ χ needs to be kept non-vanishing.
Though not infinite, the cross section is resonantly enhanced with finite Γ χ as 1/Γ χ , and it becomes larger for smaller Γ χ ∝ λ 2 χ . This explains the λ 2 χ in the denominator of (90). The resonant enhancement of the "self-energy correction" in figure 5 a) is not present in the "vertex correction" in figure 5 b). This allows us to neglect the vertex correction, which is of the same loop order as the self-energy correction. Quantitatively there is not much difference between the rate at finite ω Γ χ and ω = 0, as one can see from figure 10. The resonant enhancement would be weaker if we were dealing with fermions with gauge interactions in the loop (instead of χ). The reason is that the damping Γ χ in the scalar model only appears at two-loop order (i.e. from the setting sun diagram in figure 3), while the mass correction is of order λ χ (from the tadpole in figure 2). This makes Γ χ /Ω χ very small for the scalar model.
Finally, in section II of [70] it has been claimed that all dissipative effects vanish if the light field χ is in a vacuum state. This appears to be confirmed by our result for Π R ϕ . Physically this conclusion seems, however, counter-intuitive. The basic laws of statistical mechanics imply that the energy in an interacting physical system should relax to a state of equilibrium, which implies equipartition of the energy amongst all degrees of freedom. If dissipation were absent for T = 0, then all energy would forever remain trapped in ϕ if the initial state of the system contains no χ-particles. The reason why we find vanishing dissipation rates for T → 0 in the present calculation lies in the approximation ω = 0. Note that ImΠ R ϕ (ω) in (83) gives a non-zero contribution for finite ω and T = 0. The precise determination of the physical damping rate at T = 0 would, however, require some extra work and a more consistent treatment of renormalization. 11 If we were dealing with real quasiparticles in the initial and final states (rather than those with vanishing external four-momentum in the condensate), then this divergence would be regularized by the thermal mass; there would be no such case as ω = 0 because the self-energies are evaluated at the quasiparticle mass shell. In a plasma the quasiparticles at rest are massive due to the thermal mass even if their vacuum mass is zero, such as for photons. In the limit T → 0 the thermal mass disappears and only the vacuum mass remains (which can be chosen to be zero), but the divergence from the Bose-Einstein distribution also disappears. a) Figure 8: Cut through the fish diagram with zero width for the propagators. Solid lines correspond to φ and dashed lines to χ. Due to the optical theorem at finite T , the imaginary part of the fish diagram can be interpreted in terms of physical processes in which the external legs and cut propagators appear as initial and final states [76][77][78][79][80]. The diagrams on the right show those processes in which initial and final state both contain particles, as only these can contribute in the present calculation. Note that all χ-particles here are on shell. Diagram a) corresponds to the process φφ → χχ, in which two φ-quanta out of ϕ-condensate become two χ-quanta. If the total four-momenta of the two φ-quanta is (ω, 0) with ω < 2M χ , this process is kinematically not allowed for on-shell χ-particles. Diagram b) represents the scattering χφ → χφ between a χ-quantum and a single φ-quantum in the condensate (in both, initial and final state). If the sum of four-momenta of the φ-quanta is (ω, 0) with any ω = 0, then it is clear that such process is also forbidden for on-shell χ-particle due to the energy-momentum conservation. Diagram c) represents a decay process χ → χφφ. Here the initial state contains a single χ-quasiparticle with non-vanishing four-momentum. In the final state, there is also a χ-quasiparticle and two φ-quanta with combined four-momentum (ω, 0) have been added to the ϕ-condensate. This process is kinematically forbidden on-shell for any ω = 0, but the cross section is divergent for ω = 0. Similarly, in diagram d) an φ-quantum in the ϕ-condensate "decays" into an φ and two χ-quasiparticles. The resulting total cross section is vanishing for any finite ω < 2M χ and diverges at ω = 0. These intuitive interpretations hold at the level of individual quanta: The ϕ-condensate is a superposition of infinitely many φ-quantum states with different particle numbers. However, if one decomposes the condensate into eigenstates of the free Hamiltonian (with well-defined particle number), then the same argument can be drawn for any multi-particle state.   Figure a), b) and c) correspond to the cases of ω > 2M χ Γ χ , 2M χ > ω > Γ χ and 2M χ Γ χ > ω, respectively. In figure a), the region where two mass-shell curves intersect gives rise to the non-vanishing contribution from kinematically allowed processes of the on-shell χ-particles, such as decay and scattering.  Figure 11: Spectral densities ρ χ (p 0 ) (blue line) and ρ χ (p 0 − ω) (red line) as a function of p 0 for fixed |p| (e.g. for |p| = 0). Figure a), b) and c) represent the same cases as in figure 10.

Main results
Effective equation of motion for the field ϕ -In section 3 we have shown from first principles that the expectation value ϕ of a scalar field in a medium approximately follows the Markovian equation of motion (57), which is of the type (5), if the coupling to the medium is weak and the properties of the medium change sufficiently slowly. That is, the system can be characterized by a complex valued effective potential. Its real part V(ϕ) is the usual finite temperature effective potential, and the imaginary part Γ ϕ is the dissipation coefficient. They depends on temperature and the field expectation value ϕ. The equation of motion (57) and the general expressions for V(ϕ) in (60) and Γ ϕ in (61) are amongst the main results of this work. The expression (61) extends our calculation of Γ ϕ in [72], which is valid near the potential minimum, to the case of large field values in the slow-roll phase. They appear to be consistent with the equation of motion for the field φ found in the literature (see e.g. [30,33] and references therein). However, our derivation is much shorter, as we directly seek an equation of motion for ϕ. Together with the expressions (27) and (28) for the propagators, the equation of motion (5) with V(ϕ) in (60) and Γ ϕ in (61) allows to describe the nonequilibrium dynamics of a slow-rolling scalar field entirely in terms of Markovian equations.
For small values of ϕ, and if the medium is composed of a sufficiently large thermal bath, it is well-known that the field φ is exposed to Brownian motion, while its expectation value ϕ performs damped oscillations and relaxes to its minimum on a time scale that is given by the thermal quasiparticle width in the plasma. We have recovered this behavior as a limiting case in (62). Far away from the potential minimum the behavior is very different.
Effective potential and dissipation coefficient in a scalar theory -In section 4 we have applied the formulae derived in section 3.2 to a simple scalar model. Analogue results for the small field case discussed in section 3.1 are given in [72]. The finite temperature effective potential and damping coefficient are approximately given by the analytic expressions (77) (equivalently (73)) and (91). Our results appear to be in some tension with the previous literature, which we discuss in section 5.2.
Loop integrals at finite temperature -In the appendices we have provided general expressions and analytic estimates for nontrivial two-loop integrals for setting-sun diagrams at finite temperature. This includes the setting-sun diagram at vanishing external fourmomentum (appendix A.2) and thermally on-shell cases (appendix B.2). In appendix C we have given analytic results for the angle integrals in setting-sun diagrams. These formulae will be very useful for further calculations in finite-temperature scalar field theory.

Tension with previous results
Our results for the damping coefficient Γ ϕ differ from those obtained in the previous literature in two ways.

ω = 0 versus ω = M η
A main difference between the damping coefficient (52) during oscillations near the potential minimum and the expression (61) for large field values in a slow-roll phase is that the diagrams in the former are evaluated on the quasiparticle mass shell, while in the latter they have to be evaluated at vanishing external energy (i.e. at vanishing external fourmomentum, which is off shell). In spite of this, the coefficient (52) is frequently used in the literature in all regimes, which is clearly incorrect. This point has previously been realized by some authors, see e.g. [30,74,81]. In spite of this realization, the authors of [30,81] used the on-shell damping rate (52) in their calculations, probably for simplicity. The authors of [74], on the other hand, calculated Γ ϕ with vanishing external four-momentum, but their results differ from ours. We clarify the reason in section 5.2.2.
The fact that the damping coefficient in (61) is evaluated at vanishing external fourmomentum raises the questions to which degree our results are sensitive to various known infrared problems of thermal field theory [82][83][84], and whether further resummations and the inclusion of vertex (and ladder) diagrams are necessary to obtain physically consistent results. For example, in the calculation ofΠ R ϕ we neglected the vertex diagram shown in figure 5 b). This can be justified in the present case because it is not subject to the resonant enhancement 1/Γ χ of the fish diagram in figure 5 a), cf. (85). In the calculation of Π R ϕ , the main argument in favor of our approach is that the thermal masses in (74) -(75) already appear at order λ χ , λ φ and √ h, while the thermal width is of higher order in the couplings. This tends to suppress all sorts of finite width and vertex corrections. However, this handwaving argument is of course not a strict proof, and we postpone the clarification of these issues to future work.

Evaluation of the fish diagram
The dissipation coefficient Γ ϕ is given by the imaginary part of the retarded self-energy from the fish diagram, ImΠ R ϕ (ω). The fact that it should be evaluated at vanishing external four-momenta (i.e. ω → 0 in (81)) was previously pointed out in [74,75]. However, the results obtained there differ from ours. In the following we try to understand the origin of the discrepancy. Again recalling the optical theorem, we can identify the region in the integration volume of (81) where both spectral densities are thermally on-shell (i.e. at the intersection of two mass-shell curves in figure 10 a)) as the non-vanishing contribution from kinematically allowed processes (e.g. decay and scattering) of the on-shell χ-quasiparticles. From the kinematic considerations given in figure 8, it is clear that such a region exists for ω > 2M χ . For Γ χ Ω χ it strongly dominates the integral (81), see e.g. [61,72] for a detailed discussion of this integral. In this case one can use the zero-width approximation ρ χ (p 0 ) 2πsign(p 0 )δ(p 2 0 − Ω 2 χ ) to evaluate the integral. For Γ χ < ω < 2M χ (see figure 10 b)), the use of this approximation is not allowed. In this case the integral is dominated by the regions in which one of the ρ χ is on-shell (see figure 11 b)). Within these regions, the part, where the distribution functions f B have their maxima, gives the biggest contribution. This led the authors of [74,75] to the conclusion that the integral in (81) for ω → 0 would be dominated by the p 0 0 region and that one could use the replacement to evaluate the integral. 12 However, for the case ω Γ χ under consideration in this work (and in their papers), the pole regions overlap in the entire integration volume, as depicted in figures 10 c) and figure 11 c), and the integral is always dominated by the on-shell region p 0 Ω χ . Hence, the replacement (92) cannot be justified in this case. Therefore we have used the full expression (82) for the spectral densities and obtained the more involved results (83) and (85).

Conclusions and outlook
We have studied the effective action and equation of motion for the expectation value ϕ of a scalar field φ in a dense medium from first principles of nonequilibrium quantum field theory. We focused on two cosmologically important processes, damped oscillations near the ground state and a slow-roll phase in a flat potential. In a series of controlled approximations, we showed that these processes can be described in terms of Markovian effective equations of motion for both, the occupation numbers in the plasma and the field expectation value. This allows to describe the system in terms of a complex valued effective potential. The real part V(ϕ) is the usual effective potential at finite temperature, which includes radiative and thermodynamic corrections to the scalar potential. The imaginary part Γ ϕ is the dissipation coefficient that leads to damping and particle production.
When applying our results to the damped oscillations near the potential minimum, we found that our method reproduces the well-known results that have been obtained by studying the Brownian motion of φ near the ground state. As expected, the effective potential and damping coefficient for ϕ coincide with the thermal mass and width of quasiparticles in the plasma. However, far away from the minimum we found a very different behavior. In particular, during a slow-roll phase, the loop integrals that determine the coefficients in the effective action have to be evaluated at vanishing external four-momentum and in resummed perturbation theory. This shows that the common practice, which uses the thermal mass and width of quasiparticles as the order-of-magnitude estimates for the thermal corrections to the effective potential and dissipation coefficient, is clearly incorrect in this situation.
We illustrated our results in a simple scalar model. Our results for the dissipation coefficient significantly differ from those obtained in [75]. We identified the origin of the difference and argued that an analytic estimate for the damping coefficient has been applied beyond the range of its validity in the previous works. In this context, we provided explicit expressions and analytic estimates for the setting-sun and fish diagrams at finite temperature. These will be very useful for future computations in scalar field theories at finite temperature.
Our results mark an important step forward towards a quantitative understanding of the evolution of scalar fields in the early universe. For instance, they can be applied to study the fate of moduli, curvatons, axions and other scalars with a flat potential during and after reheating. They may also be used to models of warm inflation. Finally, the scalar field may also represent an order parameter in applications outside the domain of cosmology. Our analysis is, to the best of our knowledge, the most systematic and comprehensive one of the subject to date. However, it still relies on several restrictive assumptions. We have used thermal propagators in Minkowski space to evaluate loop integrals. This can only be justified if the particles in the loop have reached a local kinetic equilibrium, and when the thermal masses in the plasma are larger than the rate of Hubble expansion. Both of these assumptions should be relaxed in future work. Finally, the model in which we illustrate our results is a simple toy model in which the primordial plasma is composed only of another real scalar field. In most realistic applications, the scalar field in question couples directly or indirectly to fermions with non-Abelian gauge interactions. Some progress towards an inclusion of medium effects on the effective potential has recently been made e.g. in [33,34,50,51,69,70,72,[85][86][87][88][89], and phenomenological implications have been discussed in [73,[90][91][92][93][94][95][96][97][98], but it is still a long way to go to a complete quantitative understanding of scalar fields in the early universe in realistic models.
A The self-energy Π − ϕ from setting-sun diagram The main goal of this appendix is to evaluate the self-energy in (59) from the setting-sun diagram shown in figure 12, using the zero-width spectral density (38) neglecting ρ cont . We consider the spectral self-energy Π − ϕ , which is related to the imaginary part of the retarded self-energy.

A.1 General expression
The general expression of Π − ϕp (p 0 ) = 2iImΠ R ϕp (p 0 ) is given by which can be derived analogously to (81). After performing an integral over q for ηpropagator, we obtain (suppressing subscript ϕ for Π − ) with Here The D terms represent decay and inverse decay η ↔ ηχχ. Note that the integration limits in (95) imply that D p (p 0 ) vanishes for p 0 < 2M χ +M η , as expected from energy-momentum conservation in the processes η ↔ ηχχ that are impossible for on-shell particles. S p1 (p 0 ) and S p2 (p 0 ) correspond to the damping by scattering processes ("Landau damping"), ηη ↔ χχ and ηχ ↔ ηχ, respectively. The variables x and y are the cosines of the two nontrivial angles, i.e. x = p·k |p||k| , y = p·l |p||l| , so that I(a, b, c) represents integral over angles. The integrals cannot be performed analytically in general. For the special case p → 0 under consideration, we can, however, find analytic approximations, see appendix A.2 .
A.2 Analytic approximation for |p| = 0 and ω → 0 Here we derive approximate estimates of Π − 0 (ω) for zero spatial momentum and arbitrarily small energy.
In (94), if the energy p 0 = ω is smaller than 2M χ + M η , then the decay term D doesn't contribute as the decay process is not allowed by the energy conservation. Therefore let us focus on Landau damping terms (96) and (97). The integral over Ω k does not vanish, only if there exists an overlap between the integral region of Ω k and the region allowing a 2 j < 1 (with j = 1, 2) due to the step function θ(1 − a 2 j ) from the angle integral in (175) (now with a = a j ).
Let us find the range of Ω k that allows a 2 j < 1. Using the definition of a j given in (108) (with p = 0) and introducing dimensionless variables as the range of t for a 2 j < 1 reads with t j± being solutions of a 2 j = 1, given (up to the first order in infinitesimalω) by First of all, in order to have real t j± , it should be that u ≥ 2 (i.e. M η ≥ 2M χ ). If u < 2 (i.e. M η < 2M χ ), then it follows that a 2 j > 1 and integrals vanish. When u ≥ 2, one can see that t 2± < 0 for arbitrarily smallω since (u − 1)s > u(u − 2)(s 2 − 1), so there is no overlap between the integral range of Ω k and the range allowing a 2 2 < 1. S 2 thus vanishes due to the angle integral 13 .
For S 1 corresponding to Landau damping ηη ↔ χχ, which is the only potential nonvanishing contribution, one can show that M χ t 1± > Ω + 1 for arbitrarily smallω, so that the integral over Ω k becomes In the second and third equality we have analytically performed the integral and taken limit ω → 0 to obtain a result up to the first order in ω. Here which depend on Ω l and M η through s and u.
Let us now consider the Ω l -integration of (115), i.e.
This integral does not seem to allow an exact analytical result, so we try to get an approximate estimate of this integral. First, when u is close to 2 we use the fact that the integrand has a peak 14 around s = 2 (i.e. Ω l = 2M χ ) as shown in figure 13 and fix the logarithm at the peak to pull it out of the integral. Then we can perform the integral (with lower limit M 0 instead of M χ ) to get where we have put an overall factor 2 to compensate the change of the lower limit of the integral. Indeed numerical analysis shows that the above expression with M 0 = 2M χ is a good approximation of the integral for 2 ≤ M η /M χ ≤ 3.
On the other hand, if u becomes larger, then Ω 1+ and Ω 1− approach the original integration boundary (see figure 14). Therefore for sufficiently big u, one can take Ω 1+ → ∞ and Ω 1− → Ω + 1 = max M η − Ω l , M χ and then it is possible to analytically perform the integral giving We numerically checked that the above result agrees with the exact integral to a factor of Altogether we obtain 15 an approximation of Π − 0 (ω) for infinitesimally small ω (to a 14 To see that the integrand has a peak, note that the logarithm is zero at s = 1, at which Ω 1+ = Ω 1− , d Ω l falls with s ≥ 1. Furthermore, we numerically verified that the position of peak 1 < s 2 little changes with temperature T . 15 So far we have considered only S 1 (ω) in (94) and have to include S 1 (−ω) to get Π − . Inclusion of this term gives rise to an overall factor 2 since S 1 (ω) = −S 1 (−ω) if 3M χ M η (120 and zero otherwise (i.e. if M η ≤ 2M χ ). When M χ M η T , the above expression (i.e. the second one) can be approximated, to leading order, as B The self-energy Π − χ from setting-sun diagram In this appendix we calculate the spectral self-energy Π − for the self-interacting scalar particle χ from the setting-sun diagram shown in figure 15. In particular we find approximate estimates of Π − for the on-shell case (i.e. when the external momentum is on-shell) in appendix B.2. This is needed to determine Γ χ in (85) in section 4.2.2. For p = 0, Π − has already been calculated in [72] and here we generalize the results to arbitrary p = 0.

B.1 General expression
The spectral self-energy Π − χ from the setting-sun diagram in figure 15 has a general expression which can be derived in the same way as (93). After integrations over q and some angles of momenta k and l, we get with Here The notations here are the same as in appendix A.1. The integration limits imply that D p (p 0 ) vanishes for p 0 < 3M χ , due to the energy-momentum conservation in the decay and inverse decay processes χ ↔ χχχ, which are not allowed for on-shell particles.

B.2 Approximate estimates for on-shell quasiparticles
Here we derive approximate estimates of the self-energy (123) for the on-shell particle with arbitrary spatial momentum. In the on-shell case (i.e. p 0 = Ω p ), we have 16 16 Note that D terms and S  Using the result on the angle integral given in (198) in appendix C.2.2, the above integral becomes In the second line, symmetry with respect to Ω k ↔ Ω l has been taken into account, resulting in overall factor 2. Here I, II and III refer to the corresponding integral regions in Ω k -Ω l plane as shown in figure 16.
The first line (integral over region I) can be analytically performed, giving When Ω p = M χ , F [I] s = F s and (138) corresponds to the result of the entire F s for zero mode, so that Now let's turn to integrals over region II and III. Concerning II, at most one integral can be performed, giving The above integral and the one over region III, can not be done analytically. Therefore we try to obtain approximations of the integrals over II and III by considering cases of Ω p < T and Ω p > T separately, under the assumption M χ < T .
Here we have introduced variables t = Ω l /Ω p , s = Ω k /Ω p and x = M χ /Ω p . From above expressions for the G functions it follows that We have numerically checked that for x < 1/2 (i.e. Ω p > 2M χ ) the G functions in (146)-(148) change rather slow with x and that for x < 1/3 (i.e. Ω p > 3M χ ) they are little different from the values at x = 0. Thus for 3M χ Ω p < T we can replace the G functions in (144)-(145) by their values at zero, giving Together with (138), the entire result for the self-energy (134) reads -Another way: The function F s (Ω k , Ω l , A), given in (127), has a maximum at Ω l = M χ and Ω k = Ω p (also at Ω k = M χ and Ω l = Ω p , by symmetry Ω l ↔ Ω p ). Around this point, F s is dominated by the term f B (Ω l ) 2 for Ω l < T . Then one can use f B (Ω l ) ≈ T /Ω l to see that |l| f B (Ω l ) 2 has a peak at |l| ∼ M χ . Therefore the integrand in (141) has a maximum at |l| ∼ M χ and in a similar way one can show that the integrand in (143) has a peak when A 2 − M 2 χ ∼ M χ . Using these facts, we can obtain approximations to (141) and (143) by replacement |l| → M χ in (141) and A 2 − M 2 χ → M χ in (143). Then the integrals can be analytically performed, giving and These expressions have more involved forms than (152)-(153), which can be seen as approximations to the above expressions.

B.2.2 Ω p T
In this case the integrands in (141) and (143) are maximal at |l| ∼ T and A 2 − M 2 χ ∼ T , respectively. Therefore we can obtain approximations by replacing |l| → T in (141) and Here we have replaced M χ in the low limits of the integral by some constant M 0 , which we choose M 0 ∼ T /2. For Ω p T , from above we get and this is the result of the whole F s since F [I] s 0 for Ω p T (see (139)). To sum up, when M χ T the spectral self-energy for on-shell χ-particle with any momentum p can be approximated as

C The angle integral in setting-sun diagrams
In this appendix we estimate the angle integrals for setting-sun diagrams I(a, b, c), see (105) Here A is a function of p 0 , Ω k and Ω l , see e.g. (129). First we consider the cases with zero external momentum p = 0 (on-shell or off-shell) and then turn to the on-shell cases with non-zero momentum p = 0.

C.1 Zero mode (p = 0)
In this case we have b = 0 = c and z = a, so that If a 2 < 1 and x 2 < 1, there exists −1 < y ± < 1 that solve I(x, y) = 0, so that one can write (see figure 17) I(x, y, z) = −(y − y − )(y − y + ) with One can easily see that |y ± | < 1 for a 2 < 1 and x 2 < 1. Then assuming a 2 < 1 the angle integral can be performed, giving Note that in the second line the y-integral gives π irrespectivly of x, though y ± depends on x. If a 2 > 1, there exists no real y ± for x 2 < 1. This implies I(x, y) < 0 for x 2 < 1 and y 2 < 1, giving vanishing result for the integral due to θ(I(x, y)) in the integrand. The final result for zero mode is thus Note that a depends on p 0 , |k|, |l| and M χ by (167), so that the step-function θ(1 − a 2 ) in (175) restricts the integral region on Ω k -Ω l plane to the one allowing a 2 < 1.
C.2 Non-zero mode (p = 0) Here we generalize the zero-mode result obtained above to the cases of non-zero mode. First we consider generic cases, which include off-shell as well as on-shell, and then focus on the latter, which is needed to evaluate Γ χ in (85) in section (4.2.2).
The range of x satisfying this condition is where x ± are solutions of (a + c x) 2 = (1 + b 2 − 2b x), and they are real when b 2 + 2abc + c 2 + b 2 c 2 > 0.
Furthermore it follows from (167) that The right hand side is positive over the integral region in (124) and (125) since it has been chosen in such a way that A 2 > M 2 χ holds, by using Ω ± for the integral limits. Thus x ± are always real on the integral region.
If X = ∅, then the angle integral vanishes. If X = ∅, we define x − and x + as a lower and upper bound of X, respectively and then the angle integral can be performed as Here x − = max[−1, x − ] and x + = min [1, x + ] provided that X = ∅. Since x ± are functions of a, b and c, the condition X = ∅ can be reduced in terms of these parameters, which are functions of energy, momenta and masses of particles through (167). In the following we consider on-shell cases and obtain analytical results for the angle integral.
Using the second property above in (182), we get and one can show that Therefore Then in (188) we have In the second line on the right hand side we have used the fact that x ± solve 1 + b 2 − 2b x = (a + c x) 2 . Noting that a + cx − < 0 and using properties for c > 1, we obtain the final expression for the angle integral (188), In the last line we have used (184) and assumed 18 A ≥ M χ (i.e. Ω k + Ω l > Ω p + M χ in (125)). In terms of energies and momenta, the result above reads if Ω k , Ω l < Ω p and M χ ≤ Ω k + Ω l − Ω p 0 otherwise. (198)