{'fraction_non_alphanumeric': 0.04302373928366365, 'fraction_numerical': 0.012037935335696042, 'mean_word_length': 4.472156013001084, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thermalization of fermionic quantum fields

We solve the nonequilibrium dynamics of a 3+1 dimensional theory with Dirac fermions coupled to scalars via a chirally invariant Yukawa interaction. The results are obtained from a systematic coupling expansion of the 2PI effective action to lowest non-trivial order, which includes scattering as well as memory and off-shell effects. The dynamics is solved numerically without further approximation, for different far-from-equilibrium initial conditions. The late-time behavior is demonstrated to be insensitive to the details of the initial conditions and to be uniquely determined by the initial energy density. Moreover, we show that at late time the system is very well characterized by a thermal ensemble. In particular, we are able to observe the emergence of Fermi--Dirac and Bose--Einstein distributions from the nonequilibrium dynamics.


Introduction
The abundance of experimental data on matter in extreme conditions from relativistic heavy-ion collision experiments, as well as applications in astrophysics and cosmology has lead to a strong increase of interest in the dynamics of quantum fields out of equilibrium. Experimental indications of thermalization in collision experiments and the justification of current predictions based on equilibrium thermodynamics, local equilibrium or (non-)linear response pose major open questions for our theoretical understanding [1,2]. One way to resolve these questions is to try to understand quantitatively the far-from-equilibrium dynamics of quantum fields, without relying on the assumption of small departures from equilibrium, or on a possible separation of scales that forms the basis of effective kinetic descriptions [3]. In contrast to close-to-equilibrium approaches, the quantum-statistical fluctuations of the fields are not assumed to be described by a thermal ensemble. Moreover, one goes beyond the range of applicability of the usual gradient expansion and dilute-gas approximation.
In contrast to thermal equilibrium, which keeps no information about the past, nonequilibrium dynamics poses an initial value problem: timetranslation invariance is explicitly broken by the presence of the initial time, where the system has been prepared. The question of thermalization investigates how the system effectively looses the dependence on the details of the initial condition, and becomes approximately time-translation invariant at late times. According to basic principles of equilibrium thermodynamics, the thermal solution is universal in the sense that it is independent of the details of the initial condition and is uniquely determined by the values of the (conserved) energy density and of possible conserved charges. 1 There are two distinct classes of universal behavior, corresponding to Bose-Einstein and Fermi-Dirac statistics respectively.
The description of the effective loss of initial conditions and subsequent approach to thermal equilibrium in quantum field theory requires calculations beyond so-called "Gaussian" (leading-order large-N, Hartree or mean-field type) approximations [4,5,6]. Similar to the free-field theory limit, these approximations typically exhibit an infinite number of additional conserved quantities which are not present in the underlying interacting theory [7,8].
These spurious constants of motion constrain the time evolution and lead to a non-universal late-time behavior [9]. It has recently been demonstrated in the context of scalar field theories, that the approach to quantum thermal equilibrium can be described by going beyond these approximations [10,9]. In particular, this has been achieved by using a systematic coupling expansion [10] or a 1/N expansion to next-to-leading-order [9,11] of the two-particle irreducible generating functional for Green's functions, the so-called 2PI effective action [12,13,14,15]. In this context, the corresponding equations of motion have been shown to lead to a universal late time behavior, in the sense mentioned above, without to have recourse to any kind of coarsegraining.
In this work, we study the far-from-equilibrium time evolution of relativistic fermionic fields and their subsequent approach to thermal equilibrium. Nonequilibrium behavior of fermionic fields has been previously addressed in several scenarios [5,6,16,17,18,19], including the full dynamical problem with fermions coupled to inhomogeneous classical bosonic fields [20]. Here we go beyond these approximations and compute the nonequilibrium evolution in a 3+1 dimensional quantum field theory of Dirac fermions coupled to scalars in a chirally invariant way. As a consequence, we are able to study the approach to quantum thermal equilibrium. For the considered nonequilibrium initial conditions we find that the late-time behavior is universal and characterized by Fermi-Dirac and Bose-Einstein distributions, respectively. The results are obtained from a systematic coupling expansion of the 2PI effective action to lowest nontrivial (twoloop) order, which includes scattering as well as memory and off-shell effects. The nonequilibrium dynamics is solved numerically without further approximations. We emphasize that, given the limitations of a weak-coupling expansion, this is a first-principle calculation with no other input than the dynamics dictated by the considered quantum field theory for given initial conditions.
In Sect. 2, we review the 2PI effective action for fermions, which we use to derive exact time evolution equations for the spectral function and the statistical two-point function in Sect. 3. We discuss the Lorentz structure of our equations in Sect. 4 and exploit some symmetries in Sect. 5. In Sect. 6, we specify to a chiral quark model for which we solve the nonequilibrium dynamics. The initial conditions are discussed in Sect. 7 and some details concerning the numerical implementation are given in Sect. 8. The numerical results are presented and discussed in Sect. 9 and Sect. 10. We attach two appendices discussing in detail the quasiparticle picture we used to interprete some of our results.

2PI effective action for fermions
We consider first a purely fermionic quantum field theory with classical action for i = 1, . . . , N f "flavors" of Dirac fermions ψ i , a mass parameter m f and an interaction term V (ψ, ψ) to be specified below. Here ∂ / ≡ γ µ ∂ µ , with Dirac matrices γ µ (µ = 0, . . . , 3). Summation over repeated indices and contraction in Dirac space is implied. All correlation functions of the quantum theory can be obtained from the corresponding two-particle-irreducible (2PI) effective action Γ. For the relevant case of a vanishing fermionic "background" field the 2PI effective action can be written as [13] Γ The exact expression for the functional Γ 2 [D] contains all 2PI diagrams with vertices described by V (ψ, ψ) and propagator lines associated to the full connected two-point function D. In coordinate space the trace Tr includes an integration over a closed time path C along the real axis [21], as well as integration over spatial coordinates and summation over flavor and Dirac indices. The free inverse propagator is given by 3) The equation of motion for D in absence of external sources is obtained by extremizing the effective action [13] δΓ According to (2.2) one can write (2.4) as an equation for the exact inverse propagator with the proper self-energy .
(2.6) Equation (2.5) can be rewritten in a form, which is suitable for initial value problems by convoluting it with D. One obtains the following time evolution equation for the propagator: where we employed the shorthand notation z = C dz 0 d 3 z. Note that to keep the notation clear, we omit Dirac indices and reserve the Latin indices i, j, k, . . . to denote flavor.
3 Exact evolution equations for the spectral and statistical components of the two-point function To simplify physical interpretation we rewrite (2.7) in terms of equivalent equations for the spectral function, which contains the information about the spectrum of the theory, and the statistical two-point function. The latter will, in particular, provide an effective description of occupation numbers. For this we write for the time-ordered two-point function D ij (x, y) and proper self-energy Σ ij (x, y; D): 2 Note that for convenience we omit the explicit D-dependence in the notation for Σ >,< ij . Inserting the above decompositions into the evolution equation (2.7), we obtain evolution equations for the functions D > ij (x, y) and D < ij (x, y) as well as the identity

3)
2 If there is a local contribution to the proper self-energy, we write and the decomposition (3.2) is taken for Σ (nonlocal) (x, y; D). In this case the local contribution gives rise to an effective space-time dependent fermion mass term ∼ m f + Σ (local) (x; D).
which corresponds to the anticommutation relation for fermionic field operators. For later use we also note the hermiticity property and equivalently for D < ij (x, y). Note that here the hermitean conjugation denotes complex conjugation and taking the transpose in Dirac space only.
Similar to the discussion for scalar fields in Refs. [9,22] we introduce the spectral function ρ ij (x, y) and the statistical propagator F ij (x, y) defined as 3 The corresponding components of the self-energy are given by 4 and equivalently for A ij (x, y) and C ij (x, y).
Using the above notations the evolution equation (2.7) written for ρ ij (x, y) and F ij (x, y) are given by 3 Equivalently, one can decompose where we have taken the initial time to be zero, and For known self-energies the equations (3.11) and (3.12) are exact. We note that the form of their RHS is identical to the one for scalar fields [9,22]. To solve the evolution equations one has to specify initial conditions for the two-point functions, which is equivalent to specifying a Gaussian initial density matrix 5 . We note that the fermion anticommutation relation or (3.3) uniquely specifies the initial condition for the spectral function: (3.13) Suitable nonequilibrium initial conditions for the statistical propagator F ij (x, y) will be discussed below.

Lorentz decomposition
It is very useful to decompose the fields ρ ij (x, y) and F ij (x, y) into terms that have definite transformation properties under Lorentz transformation. We will see below that, depending on the symmetry properties of the initial state and interaction, a number of these terms remain identically zero under the exact time evolution, which can dramatically simplify the analysis. Using a standard basis and suppressing flavor indices we write where σ µν = i 2 [γ µ , γ ν ] and γ 5 = iγ 0 γ 1 γ 2 γ 3 . For given flavor indices the 16 (pseudo-)scalar, (pseudo-)vector and tensor components 5 We emphasize that a Gaussian initial density matrix only restricts the initial conditions or the "experimental setup" and represents no approximation for the time evolution. More general initial conditions can be discussed using additional source terms in defining the generating functional for Green's functions.
are complex two-point functions. Here we have definedtr ≡ 1 4 tr where the trace acts in Dirac space. Equivalently, there are 16 complex components for F ij , A ij and C ij for given flavor indices i, j. Using (3.9) and (3.10), one sees that they obey where Γ = {S, P, V, A, T }. Inserting the above decomposition into the evolution equations (3.11) and (3.12) one obtains the respective equations for the various components displayed in Eq. (4.2). For a more detailed discussion, we first consider the LHS of the evolution equations (3.11) and (3.12). In fact, the approximation of a vanishing RHS corresponds to the standard mean-field or Hartree-type approaches frequently discussed in the literature [4,5]. However, to discuss thermalization we have to go beyond such a "Gaussian" approximation: it is crucial to include direct scattering which is described by the nonvanishing contributions from the RHS of the evolution equations. Starting with the LHS of (3.11) one finds, omitting flavor indices (see also Ref. [16]): The corresponding expressions for the LHS of the evolution equation (3.12) for F follow from (4.4) with the replacement ρ → F . Considering now the various component (4.2) of the integrand on the RHS of Eq. (3.11), we find With the above expressions one obtains the evolution equations for the various Lorentz components in a straightforward way using (3.11). We note that the convolutions appearing on the RHS of the evolution equation (3.12) for F are of the same form than those computed above for ρ. The respective RHS can be read off Eqs. (4.5)-(4.9) by replacing ρ → F for the first term and A → C for the second term under the integrals of Eq. (3.12). We have now all the relevant building blocks to discuss the most general case of nonequilibrium fermionic fields. However, this is often not necessary in practice due to the presence of symmetries, which require certain components to vanish identically.

Symmetries
In the following, we will exploit symmetries of the action (2.1) and of the initial conditions in order to simplify the fermionic evolution equations derived in the previous section.
Spatial translation invariance and isotropy: We will consider spatially homogeneous and isotropic initial conditions. In this case it is convenient to work in Fourier space and we write and similarly for the other two-point functions. Moreover, isotropy implies a reduction of the number of independent two-point functions: e.g. the vector components of the spectral function can be written as where p ≡ | p| and v = p/p.

Parity:
The vector components ρ 0 V (x 0 , y 0 ; p) and ρ V (x 0 , y 0 ; p) are unchanged under a parity transformation, whereas the corresponding axialvector components get a minus sign. Therefore, parity together with rotational invariance imply that The same is true for the axial-vector components of F , A and C. Parity also implies the pseudo-scalar components of the various two-point functions to vanish.
CP -invariance: For instance, under combined charge conjugation and parity transformation the vector component of ρ transforms as and similarly for A 0 V and A V . The F -components transform as and similarly for C 0 V and C V . Combining this with the hermiticity relations (4.3), one obtains for these components that for all times x 0 and y 0 and all individual momentum modes.
We stress that a nonequilibrium ensemble respecting a particular symmetry does not imply that the individual ensemble members exhibit the same symmetry. For instance, a spatially homogeneous ensemble can be build out of inhomogeneous ensemble members and clearly includes the associated physics. The most convenient choice of an ensemble is mainly dictated by the physical problem to be investigated. Below we will study a "chiral quark model" with Dirac fermions coupled to scalars in a chirally invariant way. Since in this paper we will restrict the discussion of this model to the phase without spontaneous breaking of chiral symmetry, it is useful to exploit this symmetry as well.
Chiral symmetry: The only components of the decomposition (4.1) allowed by chiral symmetry are those which anticommute with γ 5 . We therefore have and similarly for the corresponding components of F , A and C. In particular, chiral symmetry forbids a mass term for fermions and we have m f ≡ 0.

Equations of motion
In conclusion, for the above symmetry properties we are left with only four independent propagators: the two spectral functions ρ 0 V and ρ V and the two corresponding statistical functions F 0 V and F V . They are either purely real or imaginary and have definite symmetry properties under the exchange of their time arguments x 0 ↔ y 0 . These properties as well as the corresponding ones for the various components of the self-energy are summarized below: The exact evolution equations for the spectral functions read (cf. Eq. (3.11)): 6 Similarly, for the statistical two-point functions we obtain (cf. Eq. (3.12)): The above equations are employed below to calculate the nonequilibrium fermion dynamics in a chiral quark-meson model.

Chiral quark-meson model
As an application we consider a quantum field theory involving two fermion flavors ("quarks") coupled in a chirally invariant way to a scalar σ-field and a triplet of pseudoscalar "pions" π a (a = 1, 2, 3). The classical action reads where π 2 ≡ π a π a and where τ a denote the standard Pauli matrices. The above action is invariant under chiral SU L (2) × SU R (2) transformations. For a quartic scalar self-interaction ∼ (σ 2 + π a π a ) 2 , this model corresponds to the well known linear σ-model [30], which has been extensively studied in thermal equilibrium in the literature using various approximations [31]. For simplicity we consider here a purely quadratic scalar potential: which is sufficient to study thermalization in this model. We note that this theory has the same universal properties than the corresponding linear σmodel. Extending our study to take into account quartic self-interaction is straightforward. It gives additional contributions to the scalar self-energies, Eqs. (6.13) and (6.14) below. These contributions can be found in Refs. [9,11] and we will point out the respective changes below.

Equations of motion for the scalar field
The 2PI effective action for nonequilibrium scalar fields has been extensively studied in the literature [14,15,10,9,11]. Here, we briefly recall the main features of the scalar sector and stress those aspects which are relevant for the present paper (for details see e.g. Refs. [9,11]). For the model considered here, Eq. (6.1), the 2PI effective action is a functional of fermionic propagators as well as scalar propagators. 7 The scalar fields form an O(4)vector φ A (x) ≡ (σ(x), π(x) ) and we denote the full scalar propagator by G AB (x, y) with A, B = 0, . . . , 3. The 2PI effective action (2.2), augmented by the scalar sector, reads 3) with the free scalar inverse propagator Similar to the fermionic case discussed above, the equation of motion for the scalar two-point function is obtained by minimizing the 2PI effective action with respect to G, and the proper self-energy is given by [13,9,11] rotation. Without loss of generality, because of chiral symmetry the effective action Γ[G, D] and its functional derivatives can be evaluated for G taken to be the unit matrix in O(4)-space: The same holds for the corresponding self-energy (6.5). Similarly, in the fermionic sector, because of chiral symmetry the most general fermion twopoint function can be taken to be proportional to unity in flavor space: Similar to the discussion in Sect. 3, the scalar spectral and statistical twopoint functions are defined as [9,22] and equivalently for the spectral and statistical self-energies Σ ρ φ and Σ F φ . These are all real functions and F -like components are symmetric under the exchange of x and y, whereas the ρ-like components are antisymmetric. The equal-time commutation relation of two scalar field operators implies [9,22] which uniquely specifies the initial conditions for the spectral function. Initial conditions for the statistical two-point function will be discussed below. Finally, the equations of motion for the scalar propagators read [9,22]: (6.10) where we have assumed spatially homogeneous initial conditions and Fourier transformed with respect to spatial coordinates. For known self-energies these are the exact equations for the theory described by the classical action (6.1) together with (6.2). In the presence of scalar self-interactions, the self-energy receives in particular a local contribution −iΣ , which amounts to a shift of the bare mass squared appearing in the above equations. In that case, the exact equations of motion have the same form as above, with the replacement [9,11].

2PI coupling expansion
We consider a systematic coupling expansion of the 2PI effective action (6.3) to lowest non-trivial order, which includes scattering as well as memory and off-shell effects, such that thermalization can be studied. The first nontrivial order in a coupling expansion corresponds to the two-loop contribution depicted in Fig. 1 (cf. also the discussion in Ref. [25]). 8 Using (6.6) and (6.7), one can express the two-loop contribution to Γ 2 directly in terms of G φ and D.
We obtain for the chirally symmetric theory: (6.12) where N f = 2 is the number of fermion flavors and N s = 4 is the number of scalar components. From there, it is straightforward to compute the spectral and statistical components of the two-loop self-energies 9 . We obtain for the scalar self-energies , (6.13) 8 We note that this approximation can also be related to a nonperturbative 1/N f expansion at next-to-leading order. 9 The relevant self-energies can be obtained with 2i and similarly for the fermionic self energies Σ ij = Σ δ ij : and for the fermion self-energies Finally, we note that for the case of a non-vanishing scalar self-interaction, the only additional two-loop contribution gives rise to a local mass shift as described above. In particular, at this order there is no additional contribution to the direct scattering part, i.e. to the RHS of Eqs. (6.10) and (6.11) relevant for thermalization.

Initial conditions
The time evolution for the fermions is described by first-order (integro-)differential equations for F and ρ: Eqs. (5.5)-(5.8). As pointed out above, the initial fermion spectral function is completely determined by the equal-time anticommutation relation of fermionic field operators (cf. Eq. (3.13)). To uniquely specify the time evolution for F we have to set the initial conditions. The most general (Gaussian) initial conditions for F respecting spatial homogeneity, isotropy, parity, charge conjugation and chiral symmetry can be written as Here n f 0 (p) denotes the initial particle number distribution, whose values can range between 0 and 1 (the definition of effective particle number distribution in terms of equal-time two-point functions is detailed in Appendix A). At late times, when thermal equilibrium is approached, this will lead to a canonical description with zero chemical potential. 10 The evolution equations (6.11) and (6.10) for the scalar correlators are second-order in time and one needs to specify initial conditions for the propagators and their time derivatives. As for the fermions, the initial conditions for the scalar spectral function is completely specified by the field commutation relations (cf. Eq. (6.9)). For the scalar two-point function F φ we consider (cf. also Refs. [9,22,26]) with an initial particle number distribution n 0 (p) and initial mode energy ǫ 0 (p) (see Appendix A). It is instructive to consider for a moment the solution of the free field equations, which can be obtained from Eqs. (5.5)-(5.6) and (5.7)-(5.8) by neglecting the memory integrals on their RHS. The solution of the fermionic free field equations with the above initial conditions reads and for the spectral functions one obtains One observes that each mode of the equal-time correlator F V (t, t, p) is strictly conserved in the absence of the memory integrals. Since this correlator is directly related to particle number (see Appendix A), this means that the latter is conserved mode by mode in this approximation. Although this is expected in the free field limit, this is of course not the case in the fully interacting theory. We emphasize that such additional conservation laws do not only appear in the free field limit, but are a property of mean-field-type approximations, which include local corrections to the bare mass, but neglect the scattering contributions described by the memory integrals on the RHS of Eqs. (5.5)-(5.8). In these approximations, the existence of this infinite number of spurious conserved quantities prevents the system to approach the thermal equilibrium limit at late times. This aspect has been discussed in detail in the context of scalar field theories in Refs. [9,32]. It is therefore crucial to go beyond such "Gaussian" approximations in order to correctly describe in particular the late-time evolution of the system in the interacting theory.

Numerical implementation
We numerically solve the evolution equations (5.5)-(5.8) and (6.10)-(6.11), together with the self-energies Eqs. (6.13)-(6.16). The structure of the fermionic equations is reminiscent of the form of classical canonical equations. In this analogy, F V (t, t ′ ) plays the role of the canonical coordinate and F 0 V (t, t ′ ) is analogous to the canonical momentum. This suggests to discretize F V (t, t ′ ) and ρ V (t, t ′ ) at t − t ′ = 2na t (even) and F 0 V (t, t ′ ) and ρ 0 V (t, t ′ ) at t − t ′ = (2n + 1)a t (odd) time-like lattice sites with spacing a t . This is a generalization of the "leap-frog" prescription for temporally inhomogeneous two-point functions. This implies in particular that the discretization in the time direction is coarser for the fermionic two-point functions than for the bosonic ones. This "leap-frog" prescription may be easily extended to the memory integrals on the RHS of Eqs. (5.5)-(5.8) as well.
We emphasize that the discretization does not suffer from the problem of so-called fermion doublers [24]. The spatial doublers do not appear since (5.5)-(5.8) are effectively second order in x-space. Writing the equations for F V (t, t ′ , x) and ρ V (t, t ′ , x) starting from (5.5)-(5.8) one realizes that instead of first order spatial derivatives there is a Laplacian appearing. Hence we have the same Brillouin zone for the fermions and scalars. Moreover, timelike doublers are easily avoided by using a sufficiently small stepsize in time a t /a s .
The fact that Eqs. (5.5)-(5.8) and (6.10)-(6.11) contain memory integrals makes numerical implementations expensive. Within a given numerical precision it is typically not necessary to keep all the past of the two-point functions in the memory. A single PIII desktop workstation with 2GB memory allows us to use a memory array with 470 timesteps (with 2 temporal dimensions: t and t ′ ). We have checked for the presented runs that a 30% change in the memory interval length did not alter the results. For a typical run 1-2 CPU-days were necessary.
The shown plots are calculated on a 470 × 470 × 32 3 lattice. (The dimensions refer to the t and t ′ memory arrays and the momentum-space discretization, respectively.) By exploiting the spatial symmetries described in Section 5 the memory need could be reduced by a factor of 30. We have checked that the infrared cutoff is well below any other mass scales and that the UV cutoff is greater than the mass scales at least by a factor of three.
To extract physical quantities we follow the time evolution of the system for a given lattice cutoff up to late times and measure the renormalized scalar mass m which is then used to set the scale. In the evolution equations we analytically subtract only the respective quadratically divergent terms obtained from a standard perturbative analysis. We emphasize that for the results presented below we use the late-time thermal mass to set the scale, and not the vacuum mass for convenience.
We made runs for a range of couplings g 2 = 0.49 -1 which show very similar qualitative behavior. Below, we present plots corresponding to g = 1, for which the time needed to closely approach thermal equilibrium is the shortest. This allows us to obtain an accurate thermalization with the lowest numerical cost.

Far-from-equilibrium dynamics
In Fig. 2 we present the time evolution of the fermion equal-time two-point function F V (t, t; p) for three momenta p. Results are given for two very different initial particle number distributions, which are displayed in the insets (see also Eqs. (7.1)-(7.3)). The (conserved) energy density is taken to be the same for both runs. In this case, since thermal equilibrium is uniquely specified by the value of the energy density, the correlator modes should approach universal values at late times if thermalization occurs.
It is striking to observe from Fig. 2 that after a comparably short time, much before the correlation modes reach their late-time values, the dynamics becomes rather insensitive to the details of the initial conditions: for a given momentum, the curves corresponding to the two different runs come close to each other rather quickly. During the slow subsequent evolution the system is still far away from equilibrium before the approach to the latetime values sets in. From both runs one observes that the characteristic time needed to effectively loose the information about the details of the initial conditions is much shorter than the time needed to approach the late-time    Figure 3: The same as in Fig. 2 but for the bosonic two-point function F φ (t, t; p) for three different momenta.
The initial particle number distributions for the two runs denote by "A" and "B" are displayed in the insets. For the employed parameters one observes that, in the scalar sector, the time needed to become effectively insensitive to the initial distributions is comparable to the time scale describing the approach to the universal late-time value.
result. Moreover, the late-time values are found to be universal in the sense that the different runs agree with each other to very good precision. Fig. 3 shows the corresponding behavior of the scalar correlator modes F φ (t, t; p). The respective initial particle number distributions in the scalar sector for the two runs are given in the inset. For the two different runs the modes having the same momenta approach each other rather slowly as compared to the fermionic sector. However, they reach their final values on a time scale which is comparable to that observed for the fermions in Fig 2. To characterize these time scales in more detail, we consider in Fig. 4 the unequal-time two-point function F V (t, t ′ ; p). As expected, if the system is to become insensitive to the details of the initial conditions, we observe that the correlation between some time t and another time t ′ is suppressed for sufficiently large t − t ′ . We note that the oscillation envelope of F V (t, t ′ , p) One observes a very good agreement with an exponential behavior. The rate is well described by the width of the corresponding spectral function in frequency space as shown in the inset (cf. the text for details).
can be well described in terms of an exponential for sufficiently late times. The corresponding damping rate approaches a constant value. To estimate the asymptotic rate we show in Fig. 4 the unequal-time two-point function F V (t, t ′ ; p) as a function of t − t ′ for the late time mt = 600. The fit to an exponential yields the damping rate γ . We emphasize that the rate γ (damp) f (p) can be related to the width of the Fourier transform of the spectral function with respect to the time difference t − t ′ . In principle, the latter involves an integration over an infinite time interval. However, since we consider an initial-value problem for finite times we know ρ V (t, t ′ , p) only on a finite interval. To overcome this problem, we fit the data for ρ V (t, t ′ , p) by a 7-parameter formula that is capable to account for the observed oscillations and damping, but which is more general than the usual 2-parameter Breit-Wigner formula. We perform the Fourier transformation on the extrapolated data. The resulting function ρ V (ω, p) is displayed as a function of frequency ω in the inset of Fig. 4. One clearly observes a nonzero width of the spectral function, the value of which may be obtained from a fit to a Breit-Wigner formula. By doing so, we obtain a very good agreement of the damping rate inferred from the width of the spectral function on the one hand and from the linear fit on the log-plot for F V (t, t ′ ; p) on the other hand.
We can use the fermion damping rate to quantify the time scale characterizing the effective loss of the details of the initial conditions: Comparing with Fig. 2, we observe that the inverse fermion damping rate at p ≃ 0 1/γ (damp) f = 15(1)m −1 characterizes well the time for which the dynamics becomes rather insensitive to the initial distributions. In contrast, we find that this time scale does not characterize the late-time behavior. For the latter, we observe to very good approximation an exponential relaxation of each mode F V (t, t; p) to its universal late-time value. Carrying out the measurement for different modes, we observe that the corresponding rate is almost independent of momentum and is given by 1/γ One observes that although the damping rates for fermions and bosons are very different, the respective thermalization rates are rather similar. For the scalars the thermalization rate is larger than the damping rate, as is observed above for the fermions. However, the difference is much less pronounced for the scalars. Similar studies in 1 + 1 dimensional quantum [9] and classical [8] scalar theories typically find a substantial difference between damping and thermalization rates. However, this may be a consequence of the stringent phase-space restrictions and therefore particular to 1 + 1 dimensional systems.
We finally note that approximate rates describing the early-time behavior may be determined by an exponential fit to the functions F φ (t, 0), F 0 V (t, 0) and F V (t, 0) in a finite time interval. These rates depend on time and approach the late-time values given above. At early times we observe an approximate exponential damping with a rate about twice as big as the latetime value for the fermions, and about half the late-time rate for the scalars.
In the previous section, we have seen that the out-of-equilibrium evolution of the system leads to a universal late-time behavior, uniquely characterized by the initial energy density. We now analyze in detail if quantum thermal equilibrium, characterized by Bose-Einstein and Fermi-Dirac statistics, is approached 11 . In thermal equilibrium, the spectral function and the statistical two-point function are not independent of each other, but are related by the fluctuation-dissipation relation. The latter is an exact relation, which can be stated in 4-dimensional Fourier space as [9,22] for the scalar correlators, where n BE (ω) = 1/[exp(ω/T )−1] denotes the Bose-Einstein distribution function. The frequency ω is the Fourier conjugate of t − t ′ (in thermal equilibrium, time-translation invariance implies that twopoint functions only depend on t − t ′ ). The value of the temperature T is determined by the energy density of the system. The fluctuation-dissipation relation for the fermion correlators is given by the equivalent expression with the replacement (n BE (ω) + 1 2 ) → (n F D (ω) − 1 2 ) in Eq. (10.1) and the Fermi-Dirac distribution n F D = 1/[exp(ω/T ) + 1]. The same relations hold as well for the statistical and spectral components of the self-energies in thermal equilibrium. They provide an unambiguous way to extract the distribution functions, which are specific for quantum thermal equilibrium. Out of equilibrium, the spectral and statistical two-point functions are completely independent in general. However, if thermal equilibrium is approached at late times, they become related by the fluctuation-dissipation relations.
To extract the answer about the late-time distributions, without relying on any assumptions, we consider the statistical and spectral correlators in Wigner coordinates. For this we express F φ (t, t ′ ; p) and ρ φ (t, t ′ ; p), as well as the corresponding fermion two-point functions, in terms of the center n BE (ω) Figure 5: The late-time ratio of the statistical two-point function and the spectral function in frequency space, both for fermions (F V /ρ V ) and for scalars (F φ /ρ φ ). In thermal equilibrium the quotient corresponds to the Bose-Einstein (BE) distribution function for scalars and to the Fermi-Dirac (FD) distribution for fermions -independently of any assumption on a quasiparticle picture (see Eq. (10.1)). The BE/FD distributions are displayed by the continuous curves parametrized by the same temperature. The value of the latter, T = 0.94m, is actually not fitted but has been taken from the inverse slope of Fig. 8. This shows the correspondence with the quasi-particle picture described in the text. coordinate X 0 = (t + t ′ )/2 and the relative coordinate s 0 = t − t ′ and write and equivalently for the other correlators. Since we consider an initial-value problem, the time integral over s 0 is bounded by ±2X 0 (cf. also the detailed discussion in Ref. [22]). If thermal equilibrium is approached for sufficiently large X 0 , then the correlators do no longer depend on X 0 and a Fourier transform with respect to t−t ′ can be performed to very good approximation (cf. the discussion in Sect. 9). The distribution functions may then be  Figure 6: The time-dependent fermion quasi-particle distribution n f (t, p) as a function of mode energy p for various times t. We have plotted the inverse slope function log(1/n f − 1), which reduces to a straight line intersecting the origin when n f (t, p) approaches a Fermi-Dirac distribution. This plot shows the data for run "A" of Fig. 2.
extracted from the quotient of the Wigner transformed two-point functions.
In Fig. 5 we show the respective ratios at late times. Both functions are in good agreement with the equilibrium distributions n F D (ω) and n BE (ω), respectively (cf. Eq. (10.1)). We emphasize that the displayed continuous curves are no separate fits. They are the Bose-Einstein and Fermi-Dirac distribution functions parametrized by the same temperature. In particular, the value for the temperature is not fitted but has been extracted from the inverse slope parameter as explained below. We stress that the above procedure to extract the distribution functions is independent of any assumption about a quasi-particle picture. However, for many practical purposes it is very convenient to have an effective description of particle number and mode energy directly in real time -without the need of a Fourier transform. An efficient description is elaborated in Appendix A. The value for the temperature in the thermal equilibrium distributions of Fig. 5 has been actually measured based on this quasi-particle picture. In   Figure 7: The time-dependent boson quasi-particle distribution n(t, p) as a function of mode energy ǫ(t, p) (see text) for various times. In this case the inverse slope function is log(1/n + 1), which reduces to a straight line in case of a Bose-Einstein distribution. This plot shows the data for run "A" of Fig. 3.
Appendix A we define the effective particle number and energy to be used, both for the fermionic and bosonic fields. For scalar field theories, this quasiparticle picture has been successfully employed previously to investigate thermalization [9,28,22]. In Figs. 6 and 7 we show the effective quasi-particle number distributions defined as for fermions and for scalars, where ε(t, p) is the quasi-particle mode energy as discussed in Appendix A. 12 The curves correspond to the initial conditions of run "A" shown in Figs. 2 and 3. One observes how the effective fermion and boson particle numbers change with time, the former approaching a Fermi-Dirac and the latter a Bose-Einstein distribution. To emphasize this point, we plot the corresponding "inverse slope functions" log(1/n f − 1) and log(1/n+1), which reduce to straight lines for Fermi-Dirac and Bose-Einstein distributions respectively. The associated inverse slopes correspond to the temperature of the thermal equilibrium distributions. We see in Figs. 6 and 7 that both inverse slope functions approach straight lines at late times. The associated temperatures for fermions and bosons are independent of time and agree very well with each other, as shown in Fig. 8. For comparison, we display in Fig. 8 the fermion inverse slope function evaluated with the bosonic effective particle number, and vice versa. This illustrates the degree of sensitivity of the inverse slope functions to the different statistics and in turn the degree of precision with which we are able to probe thermalization.

Conclusions
In this paper we have discussed the far-from-equilibrium dynamics and subsequent thermalization of a system of coupled fermionic and bosonic quantum fields. We solved the nonequilibrium dynamics beyond mean-field type approximations by calculating the complete lowest non-trivial order in a systematic coupling expansion of the 2PI effective action, which includes direct scattering as well as memory and off-shell effects. To our knowledge this is the first time that such a calculation is performed without further approximations. As a result, we show that, for various far-from-equilibrium initial conditions, the late-time behavior is universal and uniquely determined by the value of the initial energy density. Moreover, we are able to probe the approach to quantum thermal equilibrium, characterized by the emergence of Fermi-Dirac and Bose-Einstein distribution functions, with high accuracy. We emphasize that in the present calculation, besides the limitations of a coupling expansion, there is no other input than the dynamics dictated by the considered quantum field theory for given nonequilibrium initial conditions. This work can be extended in many directions. Most of the equations we derive are also valid in the phase with spontaneous symmetry breaking. Combined with earlier work on scalar theories [9,11,26], this provides a description of the nonequilibrium dynamics of the linear σ-model for QCD with two quark flavors. The model has served for many years as a valuable testing ground for ideas on the equilibrium phase structure of low energy QCD at nonzero temperature and density. With the present techniques a quantitative understanding of the out-of-equilibrium physics of this model is within reach. the case of neutral scalar fields, we recover the particle number definition used in previous studies to exhibit thermalization in purely scalar field theories [22,9].

Fermions
We start by considering the case of charged fields and construct the effective particle number from the conserved current generated by the U(1) symmetry. 14 The associated 4-current for each given flavor is ∼ψγ µ ψ. Fourier transformed with respect to spatial momenta, the expectation value of the latter can be written as J µ f (t, p) = tr[γ µ D < (t, t, p)], where the subscript stands for fermions (cf. Eq. (3.1)). In terms of the equal-time statistical two point function, its temporal and spatial components read: 15 A nice property is that the above expressions in terms of the full correlators do not contain any explicit dependence on the interaction part. In order to obtain an effective particle number, we want to identify these expressions with the corresponding ones in a quasi-particle description with free-field expressions. These are given by: 16 where Q f (t, p) = n f −n f is the difference between particle and anti-particle effective number densities and N f (t, p) = (n f +n f )/2 is their half-sum. The physical content of these expressions is simple: the temporal component J 0 directly represents the net-charge density per mode Q f (t, p), whereas the spatial part J is the net current density per mode and is therefore sensitive to the sum of particle and anti-particle number densities. Identifying the 14 As we shall see below for scalar fields, the expression for the effective particle number one obtains in this way can be directly applied to the case of neutral fields as well. 15 The constant factor in the temporal component comes from the fact that we define the current without the standard normal ordering. 16 Here, we have explicitly used our assumptions of parity symmetry and rotational invariance (cf. Sect. 5), which imply that the different spin states contribute the same, therefore the factor of two. above expressions, we define Of course, these expressions are only meaningful for the case that the interacting theory is well-described by a quasi-particle picture. The important point is that the above procedure allows one to construct an effective particle number density without knowing a priori which part of the interaction is to be considered as the "dressing" of the quasi-particles and which part describe their "residual" interactions. We note that the equal-time two-point functions on the RHS of these equations are real by definition (see Eq. (5.3)). Moreover, using the anticommutation relation for the fermion fields, it is shown in Appendix B that these definitions always satisfy 0 ≤ N f (t, p) ≤ 1 and −1 ≤ Q f (t, p) ≤ 1. These properties are important for the above definitions to be physically meaningful.
Introducing a non-vanishing net charge density per mode was useful for the above general construction. However, in the present paper, we consider only CP -invariant systems, which imply that the latter should vanish. Indeed, we see from Eqs. (5.3) that the requirement of CP -invariance imply that our above definition of net charge density per mode vanishes identically for all times and for all modes. Therefore, the effective particle and antiparticle numbers are equal and we have where n f (t, p) is the effective particle number.

Bosons
Following the same lines as above, let us consider for a moment the case of a single charged scalar field φ. The 4-current associated the corresponding and, as for the case of fermions, its expectation value has a simple expression in terms of the equal-time statistical two-point function of the charged field. In momentum space, it reads: where, as before, the constant contribution comes from the fact that we define the current without the usual normal ordering. Here, the statistical two-point function for the charged scalar field is defined as the expectation value of the anticommutator of two fields operators: It has the symmetry property F φ (t, t ′ ; p) = F * φ (t ′ , t; p), so that the above expressions are real. The corresponding quasi-particle expressions read: where Q b (t, p) and N b (t, p) have the same meaning as before in terms of effective particle and anti-particle number densities and where ǫ(t, p) is the quasi-particle energy. We therefore define: Note that the RHS of both expressions are real quantities, as they should if the LHS are to be interpreted as charge and quasi-particle number densities respectively. It remains to define the effective quasi-particle energy ǫ(t, p). 17 For this purpose, we use the free-field like expression for the average energy per mode, which in the case of a single charged scalar field can be written as ∂ t ∂ t ′ F φ (t, t ′ ; p)| t ′ =t + ǫ 2 (t, p)F φ (t, t; p) ≡ ǫ(t, p) 2n(t, p) + 1 in terms of the statistical two-point function. Therefore, one obtains for the effective quasi-particle energy: It is straightforward to show that the combination of equal-time correlators appearing on the RHS of (A.5) is indeed positive. Moreover, using commutation relations for the scalar field, one can show that the effective particle number N b (t, p), as given by Eqs. (A.4) and (A.5), is a positive quantity (see Appendix B). We note that for the case of CP -invariant systems, the effective charge density per mode (A.3) vanishes, as it should. 18 When dealing with neutral scalar fields, as it is the case in the present paper, we can use the same formula derived above. Notice that this is consistent since in that case one has F φ (t, t ′ ; p) = F φ (t ′ , t; p), which directly implies that Q b (t, p) = 0. Therefore, in the present paper, we use the definition 1 2 + n(t, p) = ǫ(t, p) F φ (t, t; p) (A. 6) together with (A.5) for the boson effective particle number n(t, p) and mode energy ǫ(t, p) for each individual scalar species [22,9].

Appendix B
Here we show that the combinations of equal-time two-point functions used in this paper to define effective particle number densities are always positive and, for the fermionic case, smaller than one. It is demonstrated to be a simple consequence of the (anti-)commutation relations of field operators.

Fermions
As a first exercise, we show that for any operators ϕ and ϕ † satisfying the anticommutation relations [ϕ, ϕ † ] + = 1 and [ϕ, ϕ] + = [ϕ † , ϕ † ] + = 0 , where the brackets denote an average with respect to any density matrix. The left inequality above is trivially obtained for example by inserting a complete sum of states between the operators ϕ † and ϕ. One obtains a sum of positive quantities which is of course positive. To show the second inequality, we introduce the hermitian operator N = ϕ † ϕ. Using the anticommutation relations above, it is easy to see that N 2 = N, from which it follows that It is clear that ∆n 2 ≥ 0 and we have just shown that N ≥ 0. One therefore conclude from the above equation that N ≤ 1, as announced.
We now come to our effective particle number densities. Using Eqs. (A.1) and (A.2) and recalling that Q f ≡ n f −n f and N f ≡ (n f +n f )/2, we get for the effective particle and anti-particle number densities: In terms of the fermionic field operators ψ(t, p) andψ(t, p), which satisfy the anticommutation relations Therefore, one can rewrite (from now on, we drop the explicit t and p dependence):

Now we introduce the operator
in term of which, where we explicitely wrote the sum over Dirac indices. For a given Dirac index, it is straightforward to check that and [ϕ α , ϕ α ] + = [ϕ † α , ϕ † α ] + = 0 , where we specialized to the Dirac basis to write the last equality of Eq. (B.3). 19 Using (B.1) for each individual α, we conclude from (B.2) that 0 ≤ n f (t, p) ≤ 1 for any time t and any momentum p. It is straightforward to repeat the above arguments to show that 0 ≤n f (t, p) ≤ 1 .
Recalling the definition of the statistical propagator in terms of field operators: F φ (t, t; p) = φ † (t, p) φ(t, p) , ∂ t ∂ t ′ F φ (t, t ′ ; p)] t ′ =t = Π † (t, p) Π(t, p) , where Π(t, p) and Π † (t, p) are the canonical momenta conjugate to φ † (t, p) and φ(t, p) respectively (e.g. Π(t, p) =φ † (t, p)), we see that where we dropped in the notation the explicit time and momentum dependence of field operators. The second inequality is nothing but Heisenberg's uncertainty relation, a direct consequence of the equal-time canonical commutation relations of field operators [33].