Thermalization of Quantum Fields from Time-Reversal Invariant Evolution Equations

We study the time evolution of correlation functions in closed quantum systems for nonequilibrium ensembles of initial conditions. For a scalar quantum field theory we show that generic time-reversal invariant evolutions approach equilibrium at large times. The calculation provides a first principles justification of Boltzmann's conjecture that the large-time behavior of isolated macroscopic systems can be described by thermal ensemble averages.


I. INTRODUCTION
One of the fundamental questions in statistical physics concerns the observation of macroscopic irreversible or dissipative processes arising from microscopic reversible dynamics. A classical approach to this question employs the separation of a system into observed and unobserved or environmental degrees of freedom and an averaging procedure over the latter [1]. The reduced system may lose the detailed information about initial nonequilibrium conditions -a necessary condition to reach thermal equilibrium. In contrast, in a closed system with a unitary time evolution no information is lost and it cannot exhibit dissipation or reach thermal equilibrium at a fundamental level.
In this Letter we study the time evolution of closed quantum systems: for a real, scalar quantum field theory we show that generic time-reversal invariant evolutions of nonequilibrium quantum fields approach thermal equilibrium. We compute time-reversal invariant evolution equations for correlation functions Tr {ρ(t 0 ) Φ(t 1 ) . . . Φ(t n )}, i.e. expectation values of products of Heisenberg fields Φ for given initial density matrix ρ(t 0 ). We explicitly calculate the large-time limit numerically for a set of different initial conditions with given average energy density. We observe that the correlation functions approach the thermal distribution asymptotically without reaching equilibrium at accessible finite times. Typical evolutions of correlation functions for spatially homogeneous fields in 1+1 dimensions are shown in fig. 2. Our calculation provides for a quantum field theory a first principles justification of Boltzmann's conjecture: the large-time behavior of isolated macroscopic systems can effectively be described by equilibrium ensemble averages. Apart from its fundamental meaning the treatment of closed quantum systems is of considerable practical importance. In many situations a clear separation into reduced system and environment is not obvious. This applies in particular to strongly coupled theories or systems far away from equilibrium where a strict separation of scales often does not exist. The nonequilibrium description of many-body quantum systems or quantum field theories has most diverse applications ranging from mesoscopic quantum devices to relativistic heavy-ion collision experiments.
Nonequilibrium quantum field theory has attracted much interest in recent years. Nonperturbative approaches based on mean field or large-N approximations as in [2][3][4] or the use of evolution equations for generating functionals [5] can give first principle insights into the dynamics of quantum fields. However, the large-time behavior of nonequilibrium quantum field theory posed an unresolved problem [3][4][5]. We study the time evolution of correlation functions in a mean field approximation, and in the Born collision approximation first suggested by Kadanoff and Baym for nonrelativistic realtime thermal Green's functions [6]. Both descriptions lead to causal and manifestly time-reversal invariant evolution equations. While the mean field approximation fails to describe the long-time behavior [2][3][4][5]7], we show that thermalization can be described in the collision approximation without further assumptions.

II. METHOD: LOOP EXPANSION OF THE 2P I EFFECTIVE ACTION
We investigate a quantum field theory for a real, scalar field ϕ with classical action The time evolution of n-point correlation functions for a given initial density matrix ρ(t 0 ) can be obtained from the generating functional by taking derivatives of Z with respect to the external sources J or K. In our notation C,x ≡ d d x C dx 0 and T c denotes time ordering along a closed time contour C starting at t 0 and chosen to support J, K at times of interest [8]. The introduction of the external bilocal source term in (2) is used to construct the generating functional for 2P I Green's functions [9]. The evaluation of the 2P I generating functional for a closed time path has been presented by Chou et al. [10] and by Calzetta & Hu [11], and we refer the reader to their work for details.
exploits the fact that the information contained in the density matrix ρ(t 0 ) can alternatively be described by specifying all initial correlations. Here the source terms α i are non-vanishing at initial times t = t 0 only. (An explicit representation of the sources in terms of the density matrix is given in [10], section 6.) The 2P I effective action is defined as the Legendre transform of Here we have absorbed α 1 into J and α 2 into K. The effective action (4) is parametrized by the macroscopic field φ(x) ≡ −i δ(ln Z)/δJ(x) = ϕ(x) and the exact connected propagator G(x, y) = ϕ(x)ϕ(y) − ϕ(x) ϕ(y) in the presence of the initial-time source terms α 3 , α 4 etc. Solutions for φ and G require where the RHS is zero for vanishing external sources J, K (t > t 0 ). Eq. (5) yields time evolution equations for the macroscopic fields φ and G. We note that these equations are obtained from an action functional by a variational principle which guarantees a unitary time evolution. We restrict the discussion to a quartic initial density matrix, i.e. α i = 0 for i ≥ 5, which is no approximation but constrains the initial state. We approximate the effective action (4) by a series expansion in orders of or, equivalently, an expansion in the number of loops for two-particle irreducible graphs [9,11,10]. For simplicity we consider in this section the three-loop result for the effective action for a vanishing field expectation value φ.
Deviations from φ = 0 can be included in a straightforward way and results are described below. For the Z 2symmetric theory this implies α 3 ≡ 0 and the effective action Γ[G; with the four-point function III. NONEQUILIBRIUM TIME-REVERSAL INVARIANT DYNAMICS The time evolution equation for the time ordered propagator is obtained from (5) and the equal-time commutation re- with the four-point function Here L 0 is the initial-time four-point function which is determined by (7) for given initial α 4 and G. Note that G * > (x, y) = G > (y, x). The evolution for G > (x, y) is time-reflection invariant and causal since the "memory integral" in (10) only depends on times smaller than max(x 0 , y 0 ). The evolution equation for G > is a nonlinear, integraldifferential equation which can be solved numerically. We use a standard lattice discretization for a spatial volume V with periodic boundary conditions and study the large volume limit to remove finite size effects. The time discretization respects time-reversal invariance. Numerical results are calculated for 1+1 dimensions. All quantities will be expressed in units of appropriate powers of m with ≡ 1. We consider spatially homogeneous fields such that G > (x, y) = G > (x 0 , y 0 ; x − y) and Fourier transform with respect to the spatial variables. The initial conditions are specified in terms of the momentum modes of the propagator G > (0, 0; p), its first derivative with respect to timeĠ, and the four-point function L 0 at initial time.

A. Two-loop (mean field) approximation
The solution of (9) for the propagator to order 0 corresponds to the free field solution and the propagator modes G > (t, 0; p) oscillate with frequency p 2 + m 2 /2π. The two-loop contribution to the effective action (6) adds a time dependent mass shift ∆m(t) = λ 2 , t; q) to the free field evolution equation. As a consequence the equation (9) for G becomes nonlinear. The mass shift term in the evolution equation is the same for all Fourier modes G > (t, t ′ ; p) and each mode propagates collisionlessly with a time dependent effective mass. Fig. 1 shows a typical time evolution of the mass shift ∆m(t) to order . As initial conditions we use the mean field thermal solution for G(0, 0; p) with inverse temperature β = 0.5 andĠ(0, 0; p) chosen to deviate from the thermal solution. For ∆m(t) one observes an initial effective damping of oscillations in the mean field approximation termed dephasing [2]. The damping is more efficient in higher dimensions [2][3][4][5]. Since the evolution equation is time-reversal invariant dissipation is absent and the effective damping is due to the superposition of oscillatory functions with continuous frequency spectrum. This absence of dissipation can be easily demonstrated for finite systems with a discrete frequency spectrum. In this case the time-reversal invariant evolution equation, which is local in time to order , leads to characteristic recurrence times after which the effective damping is lost. Since we use a lattice discretization for a spatial volume V we have a discrete frequency spectrum for finite volumes. In 1+1 dimensions we explicitly verify that the observed recurrence times for ∆m(t) scale with V or the number of degrees of freedom to infinity.
The lower curve in fig. 1 shows a typical evolution of the macroscopic field φ(t) in the mean field approximation with nonzero initial φ. The evolution of the field, but also the evolution of Fourier modes of correlation functions, are typically undamped in 1+1 dimensions. The mean field approximation does not describe the approach to a thermal equilibrium distribution [2][3][4][5]7].

B. Three-loop (collision) approximation: thermalization
A crucial improvement for the description of the longtime behavior comes from the three-loop contribution to the effective action (6). The resulting effective four-point function (10) in the evolution equation for G > introduces scattering. In fig. 2 we show the time dependence of the equal-time propagator G > (t, t; p) for three Fourier modes |p| = 0, 3, 5 and three different initial conditions with the same energy density. For the solid line the initial conditions are close to a mean field thermal solution with β = 0.1, the initial mode distribution for the dashed and the dashed-dotted lines deviate more and more from thermal equilibrium. It is striking to observe that propagator modes with very different initial values but with the same momentum p approach the same large-time value. The asymptotic behavior of the two-point function modes are uniquely determined by the initial energy density which is characteristic for thermal equilibrium. We explicitly verify that energy is conserved and that the higher n-point functions approach asymptotic values uniquely determined by the energy. Note that all higher correlation functions are expressed in terms of G and their effective thermalization results from the behavior of the two-point function. We find that small deviations from φ = 0 approach zero in the long-time limit for the initial conditions considered in fig. 2.
The evolution of correlation function modes exhibits an effective damping of oscillations. The decrease of the maximum amplitude quickly approaches an exponential behavior, but the oscillations never damp out completely. This can be seen from the logarithmic plot of the propagator zero mode G > (t, t ′ ; 0) in fig. 3. The upper curve shows |G > (t, 0; 0)| and exhibits a strong effective reduction of correlations of the field at time t with the initial field. The lower curve presents the equal-time propagator G > (t, t; 0), subtracted by its time average over 20 time steps G t and limited below to make it suitable for a logarithmic plot. We emphasize that time-reversal invariance is manifest, i.e. for each evolution towards the thermal solution there exists the reversed evolution away from equilibrium. However, if one follows these reversed evolutions generically the system approaches equilibrium again after large enough times. We note that recurrence times can be infinite in the collision approximation even for small volumes V with discrete Fourier momentum modes. This is possible because the evolution equation for the ensemble average (9),(10) is nonlocal in time. The memory effects quickly lead to contributions which can be described by the superposition of oscillatory functions with a continuous frequency spectrum. For the initial conditions employed in figs. 2 simulations for rather small volumes of order tens the correlation length are found to describe the large-time behavior well.
The efficient damping observed in figs. 2 is due to a large coupling constant λ = 10. The damping time τ , defined as the time for which the envelope amplitude of G > (t, t; 0) is reduced by a factor e, is presented in fig. 4 as a function of λ. We find that the critical coupling for which the amplitudes are undamped is zero. The qualitative thermalization behavior does not differ in the weak and the strong coupling regime in the present approximation, which neglects higher-order scattering effects. These are relevant, in particular, for the computation of the bulk viscosity or the description of critical phenomena [12,11].
Of course, there are examples, like an energy eigenstate as initial condition, where thermalization can not occur. For few body systems like a set of coupled anharmonic oscillators it is known that a transition to non-ergodic behavior can occur for small but nonzero couplings [13]. For a finite volume V and ultraviolet cutoff for the momentum modes our model can be mapped onto a finite set of coupled anharmonic oscillators. We indeed find that for small enough volumes and couplings below a certain critical coupling λ c the long-time evolution can be undamped. The values for the critical coupling strongly depend on the initial conditions and the volume or number of degrees of freedom.