format: text

text:
# Monte Carlo simulations of the classical two-dimensional discrete frustrated ϕ⁴ model

## 1 Introduction

One of the microscopic origins of the spatially modulated structures is given by the competing interactions between the particles. The Frenkel-Kontorova model, the axial next-nearest-neighbour Ising (ANNNI) model, and the discrete frustrated $`\varphi ^4`$ (DIFFOUR) model provide basic microscopic Hamiltonians for the systems of this class. These models have been studied extensively, specially the first two ANNNI1 ; ANNNI2 . Competing interactions lead to the appearance of incommensurate and commensurate regions in the phase diagrams of these models. The order parameter in these systems is a vector that determines the amplitude and the phase of modulation. For the incommensurate state, the free energy is degenerate in the phase of modulation. The two-dimensional case is of special interest here, since there is no long-range order in such planar systems if the order parameter has more than one component PP . A nonmonotonic algebraic decay of correlations appears instead of the usual incommensurate state ANNNI2 .

The correlation function in two-dimensional systems with competing interactions has at least three possible types of behaviour, as found in the spin models ANNNI2 . The disordered phase shows an exponential monotonic or modulated ($`q0`$) fall-off of the correlation function:

$$C(r)=<S(0)S(r)>\mathrm{exp}(r/\xi )\mathrm{cos}(qr).$$
(1)
The case $`q0`$ is usually called floating-fluid phase. An algebraical decay of correlations at $`r\mathrm{}`$ is a property of the floating-incommensurate (FIC) or floating-solid phase:

$$C(r)r^\eta \mathrm{cos}(qr+\phi ).$$
(2)
Finally, a commensurate phase with locked-in wavevector $`q_0`$ is not destroyed by the fluctuations in two-dimensional case. Therefore the correlation function in this case takes the form

$$C(r)\mathrm{cos}(q_0r).$$
(3)
The commensurate case includes the ferro-phase also ($`q_0=0`$).

One should note that in the two-dimensional case the phase transition from the para-phase to the FIC-phase is believed to be related with the formation and dissociation of the dislocations or vortices ANNNI2 . This scenario of a phase transition can be described by the Kosterlitz-Thouless theory which is basic for the two-dimensional $`XY`$ universality class KT ; Pelissetto ; Gupta . The power law dependence of $`C(r)`$ for this universality class appears everywhere below the critical temperature, down to $`T=0`$ or to another phase transition (to a commensurate state, for example). The critical index $`\eta `$ is equal to $`0.25`$ in the point of the phase transition ($`T_c`$) and decreases for decreasing temperature as $`\eta (T)=0.25C\sqrt{T_cT}`$, where $`C`$ is a constant ($`\eta T`$ for $`T0`$). The thermal dependence of the correlation length above the transition point is of the form $`\xi \mathrm{exp}(b/\sqrt{TT_c})`$, where $`b`$ is a constant.

One of the well studied spin models with competing interactions is the ANNNI model. The ANNNI model is a spin lattice model with $`\pm 1`$ classical spins, ferromagnetic nearest-neighbours interaction $`J_1`$ in all directions and antiferromagnetic next-nearest-neighbour interaction $`J_2`$ in a single direction. The ANNNI model in two dimensions has been studied by numerous methods ANNNI1 ; ANNNI2 . The ground state of this model at zero temperature is a ferromagnetic phase for $`J_2/J_1<0.5`$ and a $`<2>`$ phase (two spins up and two spins down $``$) for $`J_2/J_1>0.5`$. Obviously, an FIC-phase can appear only at finite temperature between the paramagnetic phase and the $`<2>`$ phase for $`J_2/J_1>0.5`$. The behaviour of the model was studied by transfer-matrix calculations and finite-size-scaling analysis Beale , Monte Carlo simulations Selke ; Sato , dynamical Monte Carlo method Barber , fermion approximation Villain  and its modification Grynberg , the interface free energy method Muller . Some of these results seem to indicate a Kosterlitz-Thouless KT  type of the phase transition from the paramagnetic to the FIC-phase and a Pokrovsky-Talapov PT  type of phase transition from this phase to $`<2>`$. The latter scenario implies a square-root power law change of the wavevector with respect to temperature near the transition to the $`<2>`$ phase. One should point out that the one-dimensional quantum ANNNI model studied in Rieger  also reveals all phases and types of the phase transitions mentioned above. The approximations give a spread of the phase transition temperature over a rather wide region. On the other hand, other calculations predict an FIC-phase with a width comparable to that spread. The recent results on a cluster heat bath simulation Sato  of systems of about $`64\times 128`$ atoms demonstrate that the FIC region is at least 2 times smaller than the previous papers predicted Beale . Further, the nonequilibrium relaxation method Shirahata  for $`6400\times 6400`$ systems does not predict any area for the FIC-phase at all, so that even the presence of this phase was questioned.

In this paper we consider the classical two-dimensional DIFFOUR, which has continuous variables and competing interactions in one direction. The DIFFOUR model was introduced in TJ  as a displacement translationally invariant model for the incommensurate phases. Instead of $`\pm 1`$ spins, this model considers double-well anharmonic 2-4 oscillators at the nodes of a lattice. Similarly, there is a harmonic nearest-neighbour and next-nearest-neighbour coupling. The potential energy can be written as

$`V`$ $`=`$ $`{\displaystyle \frac{a}{2}}{\displaystyle \underset{i}{}}u_i^2+{\displaystyle \frac{a}{4}}{\displaystyle \underset{i}{}}u_i^4+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}(u_iu_j)^2\sigma _{ij}`$ (4)
$`+`$ $`{\displaystyle \frac{d}{2}}{\displaystyle \underset{i,j}{}}(u_iu_j)^2\sigma _{ij}^{}.`$
Here $`\sigma _{ij}=1`$ for the nearest neighbours and $`\sigma _{ij}=0`$ otherwise, while $`\sigma _{ij}^{}=1`$ for the next-nearest-neighbours in the direction $`z`$ (see Fig. 1a). The phase diagram for the analogous three-dimensional model is presented in AR . The system has commensurate and incommensurate regions in the phase diagram for $`d<0`$ i.e. with competing interactions.

The ANNNI model can be considered as a limiting case of the DIFFOUR model ($`a+\mathrm{}`$). It corresponds to the order-disorder limit of the latter model. In the other limit ($`a+0`$) of the DIFFOUR model the ground state may be an incommensurate phase which is reached via a displacive transformation TJ1 . Furthermore, the greater freedom in the DIFFOUR model could also possibly lead to the existence of a broader the FIC-phase than in the ANNNI model. In view of the unclear situation in the latter we want to study the existence of FIC-phase in the DIFFOUR model. We shall do this using Monte Carlo simulations. The mean-field phase diagram can easily be determined and is presented in Fig. 1b.

Correlation functions of the model are analyzed for certain parameters; their behaviour is always one of the three scenarios mentioned above (1,2,3). An FIC-phase has been found. Indications of the Kosterlitz-Thouless scenario of the phase transition from the floating-fluid to the FIC-phase are discussed. The transition to the ferro-phase occurs via the floating-fluid phase. Surprisingly, the period of this phase increases as temperature decreases.

## 2 Algorithm

We study the two-dimensional DIFFOUR model using the Monte Carlo method. The difficulty of the simulation of the models with competing interactions is in the complexity of their configuration space which usually has many deep almost degenerate local minima. Then, it is difficult to overcome a potential barrier and to change the period of the modulation. This situation can lead to an unacceptably long computation time using conventional Monte Carlo schemes, particularly for large systems. Special algorithms should be used to achieve a reasonable result Sato ; Shirahata .

In this paper, we use the ”local-heating” algorithm which worked well for the three-dimensional DIFFOUR model AR . The idea is to modify the potential energy in certain regions of the system to allow an easier switch between local minima. Actually, the temperature of the system is not changed but the constants of the potential energy are modified. However, for the largest part of the system its potential energy is still determined by (4). The regions with a modified potential energy are not included in the calculation of the physical quantities (order parameter, correlation functions), while the usual Metropolis algorithm for Monte Carlo sampling is applied to the whole system.

We use periodic boundary conditions and the following modification of the potential energy:

$`V`$ $`=`$ $`{\displaystyle \frac{a}{2}}{\displaystyle \underset{i}{}}\theta _iu_i^2+{\displaystyle \frac{a}{4}}{\displaystyle \underset{i}{}}\theta _iu_i^4+{\displaystyle \frac{1}{4}}{\displaystyle \underset{i,j}{}}(\theta _i+\theta _j)(u_iu_j)^2\sigma _{ij}`$
$`+`$ $`{\displaystyle \frac{d}{4}}{\displaystyle \underset{i,j}{}}\theta _{(i+j)/2}(u_iu_j)^2\sigma _{ij}^{},`$
$`\theta _i`$ $`=`$ $`(1+l_1(\mathrm{exp}[(x_i^2+z_i^2)/l_2^2]`$ (5)
$`+`$ $`\mathrm{exp}[((N_xx_i)^2+(N_zz_i)^2)/l_2^2]))^1,l_1>0.`$
$`N_x`$ and $`N_z`$ are the sizes of the lattice along the $`x`$ and $`z`$-directions, respectively. The new coefficients $`\theta _i`$ (5) mean that their values in the vicinity of the lines $`x_i=0`$ and $`z_i=0`$ are smoothly decreased. Parameters $`l_1,l_2`$ determine the height and the width of this decrease, respectively. The change of the parameters mimics a local heating of the corresponding regions. The idea of this procedure is to put these parts of the system in a para-phase and to destroy the local modulation phase. If the profile $`\theta (x,z)`$ is smooth enough, the phases of modulation at the two sides of the heated line are in an ideal case independent and can randomly fluctuate. The whole system may change the period of the modulation now. Use of ”local heating” leads to faster thermalization of the system.

Additionally we realised that it is useful to simulate the systems with a larger size along the $`x`$ direction. The possible explanation is that once the period is changed at some $`x`$, it is relatively easy for this change to spread to the whole system. The probability for this initial change is proportional to the $`x`$-size.

We calculate the correlation function as follows:

$$C(\stackrel{}{r})=<u(\stackrel{}{r}^{})u(\stackrel{}{r}^{}+\stackrel{}{r})>|_\stackrel{}{r}^{}.$$
(6)
Here $`<>|_\stackrel{}{r}^{}`$ means an ensemble average. We consider two functions $`C_x=C(\stackrel{}{r}=r_x)`$ (along $`x`$ axis) and $`C_z=C(\stackrel{}{r}=r_z)`$ (along $`z`$ axis) hereafter. The presence of the heated regions is taken into account as follows. We add data to the point $`\stackrel{}{r}`$ of correlation function (6) if there is no ”local heated” region between $`\stackrel{}{r}`$ and $`\stackrel{}{r}^{}`$ and add data to the point $`\stackrel{}{r}\stackrel{}{r}^{}`$ otherwise. If one of the points $`\stackrel{}{r}`$ or $`\stackrel{}{r}^{}`$ appears to be in the ”local heated” region the data are not included in the calculation of the correlation function. This procedure reduces the influence of the finite-size effects. The data for the correlation function appear to be reliable even for distances larger than half the system size.

We check whether the approach of the ”local heating” algorithm by changing the size of the system and by starting with various initial distributions of atom’s displacements is reasonable. One should use large enough sizes of the system since it is difficult to distinguish exponential and power law decay of correlations. The situation becomes worse when we deal with modulated structures with large modulation period. We found that an optimal size for the present algorithm is about $`200\times 100`$. The total number of Monte Carlo samplings for the given parameters of the model $`a,d`$ and temperature $`t`$ is about $`10^{10}÷10^{11}`$. Parameters of ”local heating” are varied in a region $`2l_1,l_24`$. We have checked that a two-times increase of the slab size in the $`x`$ or $`z`$-direction does almost not change the exponents of power or exponential decaying or the period of the modulation within the numerical accuracy. The result does not depend on the initial configuration of the system. Thus, we believe that the presented algorithm allows us to perform accurate numerical simulations.

## 3 Results and their analysis

We present here detailed results for two values of the parameter $`d`$ at $`a=5`$: $`d=0.45`$ and $`d=0.35`$. We choose these values of the parameter $`d`$ since at least the mean-field approximation predicts various ground states for zero temperature (commensurate phase with $`N=6`$ for $`d=0.45`$ and ferro-phase for $`d=0.35`$, see Fig. 1b). The parameter value $`a=5`$ of the DIFFOUR model (4) is used since the amplitude of modulation then has a value comparable with the value of $`C(0)`$. It makes the analysis of the data easier. We have found that the modulation amplitude is much less than the value $`C(0)`$ for larger values of $`a`$ (for example, at $`a20`$). The mean-field approach for the DIFFOUR model confirms these observations.

The correlation functions $`C_z`$ and $`C_x`$ are calculated at several values of the temperature. The results for $`C_z`$ are presented in Fig. 2 and Fig. 3 for $`d=0.45`$ and $`d=0.35`$, respectively. Typical dependencies of $`C_x`$ for both values of the parameter $`d`$ are shown in Fig. 4. We use a log-log scale for the $`C_x`$ functions. The most informative dependencies here are $`C_z`$ (Fig. 2,3). Qualitatively the behaviour of the functions $`C_x`$ and $`C_z`$ is different: non-monotonic and monotonic dependence on the distance, respectively. We believe that the influence of the finite-size effects starts at $`x100`$ and at $`z70`$ and we do not consider these regions in the analysis of the data. The error bar for $`C_x`$ and $`C_z`$ at moderate values of temperature for $`d=0.45`$ and $`d=0.35`$ is about a few percent and becomes larger at lower temperatures.

We analyze the obtained results by fitting the data for $`C_z`$ (Fig. 2,3) using one of the functions (1,2). We also try to find a power law decay for $`C_x`$ (Fig. 4) if it is possible.

### 3.1 Case $`d=0.45`$

We start with the case $`d=0.45`$ (Fig. 2). We look for the possibility of a phase transition from the para-phase to the FIC-phase. The dependencies of the correlation functions $`C_z`$ at $`t1.2`$ can be fitted by equation (1), i.e. by non-monotonic exponential decay. There is no power law decay for $`C_x`$ either for these values of temperature. Thus, we conclude that it is floating-fluid phase at $`t>1.2`$ and $`d=0.45`$.

One can fit the dependencies $`C_z`$ in the vicinity of the $`t1.2`$ by nonmonotonic power law decay (2). The results of this fitting are presented in Fig.2 by solid lines. It is possible also to find a power law decay for $`C_x`$. This region is believed to be close to the possible point of the phase transition from the floating-fluid phase to the FIC-phase. The value of the index $`\eta `$ is approximately equal to $`0.25`$ at the temperature $`t_11.2`$. This coincides with the condition of the Kosterlitz-Thouless phase transition. One can find, that data for $`C_z`$ (Fig. 2) for the region $`0.55<t<1.2`$ can also be fitted by a non-monotonic power law dependence (2) with index $`\eta `$ depending on temperature. We plot the dependence of the $`(0.25\eta (t))^2`$ on $`t_1t`$ for the $`t<1.2`$ (Fig.5a). One sees a linear dependence, which is a feature of the Kosterlitz-Thouless type behaviour KT ; Gupta . Thus, the phase transition from the floating-fluid phase to the FIC-phase at $`t_11.2`$ is believed to be of Kosterlitz-Thouless type. The data for $`C_x`$ at $`0.55<t<1.2`$ can also be fitted by a monotonic power law decay. Typical dependence of $`C_x`$ for this temperature region is presented in Fig.4 ($`t=1.0`$). In our calculations, the indices $`\eta `$ for $`C_x`$ and $`C_z`$ do not coincide (the difference is about $`1020\%`$). Although, we don’t know whether this difference is of physical or numerical origin, we rather attribute it to finite-size effects in simulations. Nevertheless, we use here the index $`\eta `$ obtained from the data of fitting in the $`z`$ direction, i.e. by equation (2) for analysis of the data.

We also show in Fig. 5b how the wavevector changes as a function of temperature. We find that the modulation period is equal to $`N=6`$ at $`t0.55`$. The appearance of the domain walls (discommensurations - small regions with pseudo-period $`N6`$) can be seen at $`t0.55`$. It prevents the system to get a commensurate modulation period ($`N=6`$). The point $`t_20.55`$ is supposed to be the point of the phase transition from the FIC-phase to the commensurate phase with period $`N=6`$ (see Fig. 5b). This picture is similar to that for the two-dimensional ANNNI model Sato  for the transition from the FIC-phase to the $`<2>`$ phase. One should note that the accuracy of our results is not sufficient to determine whether the phase transition from FIC-phase to the $`<2>`$ phase is of Pokrovsky-Talapov type or not (Fig. 5). Present results are in agreement with previous conclusions by Schulz and Haldane et al that only commensurate phase with sufficiently small periodicity can transform into a FIC-phase Schulz ; Haldane . It is also possible that narrow commensurate regions exist in the temperature range $`0.55<t<1.2`$. Study of this problem demands more accurate methods of numerical simulations and analysis of the data.

### 3.2 Case $`d=0.35`$

Before presenting the numerical results, let us stress several features of the $`d=0.35`$ situation. The mean-field predicts two transitions: from the para-phase to an incommensurate one, and than to the ferro-phase. The three-dimensional DIFFOUR model indeed behaves so at certain parameters AR . On the other hand, the ANNNI model does not show two transitions; there is only a single transition from the floating-fluid (disordered) to the ferro-phase. For the present system, because of the relatively small value $`a=5`$, the modulation is much more pronounced and falls off slower than in the ANNNI model. As well, the mean-field region for the incommensurate phase is quite wide. Therefore one could expect a transition from a floating-fluid to an FIC-phase for $`d=0.35`$. Surprisingly, we did not found evidence for this; qualitatively the result coincides with that of the ANNNI model.

Results for the case $`d=0.35`$ for $`C_z`$ are presented in Fig. 3. Data for $`t1.0`$ can be fitted by non-monotonic exponential decay (1). Results of fitting $`C_z`$ are shown in Fig. 3 by dashed lines. We have performed a detailed studies of the region $`1.0<t<0.8`$ (part of them is shown in Fig. 3). Since the period of the modulation is more than $`14`$ in this region it is difficult to distinguish whether there is a non-monotonic power law decay or not. Besides, the fluctuations in this region are increased and our algorithm does not allow us to judge the type of decay of the correlations. We try to fit these dependencies by exponential and power law decay. The results for $`t=0.95`$ are presented in Fig. 3 by dashed and solid lines, respectively. Both dependencies show agreement with Monte Carlo results within the numerical error. The index $`\eta `$ for power law decay is found to be about $`0.5`$. Thus, we believe that there is no Kosterlitz-Thouless phase transition for $`d=0.35`$ since there is no agreement with a power law decay (2) with the exponent value $`\eta =0.25`$ in this temperature region. One can see that the period of the modulation is increased, while the amplitude of modulation is decreased. This situation is opposite to the case $`d=0.45`$. A typical result for the $`C_x`$ is presented in Fig. 4 ($`t=0.95`$). It is difficult to find power law decay in the region $`1.0<t<0.8`$ for $`C_x`$. Certainly, there is no index $`\eta =0.25`$ in this case either. It is difficult to obtain clear dependence of correlations decay in the vicinity of the phase transition point due to increase of fluctuations and strong influence of the finite-size effects (see Fig. 3, $`t=0.925`$). Thus, we conjecture that there is only one second order phase transition for $`d=0.35`$ from the floating-fluid phase into the ferro-phase. We estimate the temperature of this phase transition to be $`t_c0.91`$. We would like to note that one can suggest a transition from floating-fluid phase to a ferro-phase via a phase with monotonic exponential decay of correlations with $`q=0`$ (1) and the existence of a disorder line by analogy with the ANNNI model Beale . This phase may exist in a very narrow temperature region at given parameters of the model and data for $`t=0.925`$ confirm this. The question of the presence of a disorder line in the phase diagram of the two-dimensional DIFFOUR model needs an additional study.

## 4 Conclusions

In conclusion, we have studied the DIFFOUR model in two dimensions by Monte Carlo simulations. The classical DIFFOUR model has two parameters $`a`$ and $`d`$, where the latter one is the effective constant of the competition of the interactions. Two sets of parameters of the model have been considered here in detail: $`a=5,d=0.45`$ and $`a=5,d=0.35`$. A ”local heating” algorithm has been used for numerical simulations. This algorithm uses modifications of the constants and parameters of the potential energy of the DIFFOUR model. Correlation functions in both directions have been calculated for analysis of the behaviour of the model.

We have found that correlation functions for the $`a=5,d=0.45`$ can be fitted by a non-monotonic exponential decay at $`t>1.2`$ and by non-monotonic power law decay at temperatures $`0.55<t<1.2`$. The phase transition at $`t_11.2`$ of Kosterlitz-Thouless type is believed to take place with $`\eta =0.25`$ from the floating-fluid phase to the FIC-phase. The dependence of the index $`\eta (t)`$ below the point of the phase transition has a square root of $`t_1t`$ behaviour in agreement with the Kosterlitz-Thouless theory. The ground state of the model at $`t<0.55`$ ($`a=5,d=0.45`$) is a commensurate phase with modulation period $`N=6`$. The second phase transition at these parameters of the model from FIC-phase to the commensurate one occurs at $`t_20.55`$.

The correlation functions for the case $`a=5,d=0.35`$ can be fitted by a non-monotonic exponential decay for $`t>0.95`$. Dependencies at $`t0.95`$ can be fitted by non-monotonic power law decay also but with $`\eta >0.25`$. Numerical data do not allow us to distinguish exponential from power law decay for this case since the period of modulation is too large. Nevertheless, we believe that there is no FIC-phase in this case and only one phase transition from the floating-fluid phase to the ferro-phase, at $`t_c0.91`$. The ground state of the model for $`t<0.9`$ is a ferro-phase.

We would like to point out that the results of the Monte Carlo simulations (Fig. 2,3) and the mean-field data (Fig. 1) are different. Mean-field approximation predicts an incommensurate phase, while Monte Carlo results show a FIC-phase with power law decay of correlations. Mean-field theory predicts two phase transitions for $`d=0.35`$ and a region of the incommensurate phase, while Monte Carlo simulations show only one phase transition from the floating-fluid phase to the ferro-phase. Nevertheless, the low-temperature phase (ferro-phase for $`d=0.35`$ and commensurate phase with period $`N=6`$ for $`d=0.45`$) is predicted correctly in the mean-field approximation, as in the limit $`t=0`$ fluctuations are absent and the mean-field result is exact.

###### Acknowledgements.

This work was supported by INTAS project for young scientists (grant YSF 2001/1 - 135).


text:
# Thermalization of fermionic quantum fields

## 1 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 . 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 . 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: time-translation 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.222Here we consider closed systems without coupling to a heat bath or external fields, which could provide sources or sinks of energy. 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 . 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 . These spurious constants of motion constrain the time evolution and lead to a non-universal late-time behavior . 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 . In particular, this has been achieved by using a systematic coupling expansion  or a $`1/N`$ expansion to next-to-leading-order  of the two-particle irreducible generating functional for Green’s functions, the so-called 2PI effective action . 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 coarse-graining.

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 , including the full dynamical problem with fermions coupled to inhomogeneous classical  bosonic fields . 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 (two-loop) 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.

## 2 2PI effective action for fermions

We consider first a purely fermionic quantum field theory with classical action

$$S=\mathrm{d}^4x(\overline{\psi }_i(x)[i/m_f]\psi _i(x)+V(\overline{\psi },\psi ))$$
(2.1)
for $`i=1,\mathrm{},N_f`$ “flavors” of Dirac fermions $`\psi _i`$, a mass parameter $`m_f`$ and an interaction term $`V(\overline{\psi },\psi )`$ to be specified below. Here $`/\gamma ^\mu _\mu `$, with Dirac matrices $`\gamma _\mu `$ ($`\mu =0,\mathrm{},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 $`\mathrm{\Gamma }`$. For the relevant case of a vanishing fermionic “background” field the 2PI effective action can be written as

$$\mathrm{\Gamma }[D]=i\text{Tr}\mathrm{ln}D^1i\text{Tr}D_0^1D+\mathrm{\Gamma }_2[D]+\mathrm{const}.$$
(2.2)
The exact expression for the functional $`\mathrm{\Gamma }_2[D]`$ contains all 2PI diagrams with vertices described by $`V(\overline{\psi },\psi )`$ 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 $`𝒞`$ along the real axis , as well as integration over spatial coordinates and summation over flavor and Dirac indices. The free inverse propagator is given by

$$iD_{0,ij}^1(x,y)=(i/m_f)\delta _𝒞^4(xy)\delta _{ij}.$$
(2.3)
The equation of motion for $`D`$ in absence of external sources is obtained by extremizing the effective action

$$\frac{\delta \mathrm{\Gamma }[D]}{\delta D_{ij}(x,y)}=0.$$
(2.4)
According to (2.2) one can write (2.4) as an equation for the exact inverse propagator

$$D_{ij}^1(x,y)=D_{0,ij}^1(x,y)\mathrm{\Sigma }_{ij}(x,y;D),$$
(2.5)
with the proper self-energy

$$\mathrm{\Sigma }_{ij}(x,y;D)i\frac{\delta \mathrm{\Gamma }_2[D]}{\delta D_{ji}(y,x)}.$$
(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:

$$(i/_xm_f)D_{ij}(x,y)i_z\mathrm{\Sigma }_{ik}(x,z;D)D_{kj}(z,y)=i\delta _𝒞^4(xy)\delta _{ij},$$
(2.7)
where we employed the shorthand notation $`_z=_𝒞dz^0\mathrm{d}^3z`$. Note that to keep the notation clear, we omit Dirac indices and reserve the Latin indices $`i,j,k,\mathrm{}`$ 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 $`\mathrm{\Sigma }_{ij}(x,y;D)`$:333 If there is a local contribution to the proper self-energy, we write    $`\mathrm{\Sigma }_{ij}(x,y;D)=i\mathrm{\Sigma }^{(\mathrm{local})}(x;D)\delta _𝒞^4(xy)\delta _{ij}+\mathrm{\Sigma }_{ij}^{(\mathrm{nonlocal})}(x,y;D),`$  and the decomposition (3.9) is taken for $`\mathrm{\Sigma }^{(\mathrm{nonlocal})}(x,y;D)`$. In this case the local contribution gives rise to an effective space-time dependent fermion mass term $`m_f+\mathrm{\Sigma }^{(\mathrm{local})}(x;D)`$.

$`D_{ij}(x,y)`$ $`=`$ $`\mathrm{\Theta }_𝒞(x^0y^0)D_{ij}^>(x,y)\mathrm{\Theta }_𝒞(y^0x^0)D_{ij}^<(x,y),`$ (3.8)
$`\mathrm{\Sigma }_{ij}(x,y;D)`$ $`=`$ $`\mathrm{\Theta }_𝒞(x^0y^0)\mathrm{\Sigma }_{ij}^>(x,y)\mathrm{\Theta }_𝒞(y^0x^0)\mathrm{\Sigma }_{ij}^<(x,y).`$ (3.9)
Note that for convenience we omit the explicit $`D`$–dependence in the notation for $`\mathrm{\Sigma }_{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

$$\gamma ^0\left(D_{ij}^>(x,y)+D_{ij}^<(x,y)\right)|_{x^0=y^0}=\delta (\stackrel{}{x}\stackrel{}{y})\delta _{ij},$$
(3.10)
which corresponds to the anticommutation relation for fermionic field operators. For later use we also note the hermiticity property

$$\left(D_{ji}^>(y,x)\right)^{}=\gamma ^0D_{ij}^>(x,y)\gamma ^0,$$
(3.11)
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.  we introduce the spectral function $`\rho _{ij}(x,y)`$ and the statistical propagator $`F_{ij}(x,y)`$ defined as444Equivalently, one can decompose
$$D_{ij}(x,y)=F_{ij}(x,y)\frac{i}{2}\rho _{ij}(x,y)\left[\mathrm{\Theta }_𝒞(x^0y^0)\mathrm{\Theta }_𝒞(y^0x^0)\right].$$

$`\rho _{ij}(x,y)`$ $`=`$ $`i\left(D_{ij}^>(x,y)+D_{ij}^<(x,y)\right),`$ (3.12)
$`F_{ij}(x,y)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(D_{ij}^>(x,y)D_{ij}^<(x,y)\right).`$ (3.13)
The corresponding components of the self-energy are given by555 Besides the dynamical field degrees of freedom $`D_{ij}`$ we introduce quantities which are functions of these fields. These functions are denoted by either boldface or Greek letters.

$`𝐀_{ij}(x,y)`$ $`=`$ $`i\left(\mathrm{\Sigma }_{ij}^>(x,y)+\mathrm{\Sigma }_{ij}^<(x,y)\right),`$ (3.14)
$`𝐂_{ij}(x,y)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\mathrm{\Sigma }_{ij}^>(x,y)\mathrm{\Sigma }_{ij}^<(x,y)\right).`$ (3.15)
With (3.11) the two-point functions have the properties

$`\left(\rho _{ji}(y,x)\right)^{}`$ $`=`$ $`\gamma ^0\rho _{ij}(x,y)\gamma ^0,`$ (3.16)
$`\left(F_{ji}(y,x)\right)^{}`$ $`=`$ $`\gamma ^0F_{ij}(x,y)\gamma ^0,`$ (3.17)
and equivalently for $`𝐀_{ij}(x,y)`$ and $`𝐂_{ij}(x,y)`$.

Using the above notations the evolution equation (2.7) written for $`\rho _{ij}(x,y)`$ and $`F_{ij}(x,y)`$ are given by

$`(i/_xm_f)\rho _{ij}(x,y)`$ $`=`$ $`{\displaystyle _{y^0}^{x^0}}dz𝐀_{ik}(x,z)\rho _{kj}(z,y),`$ (3.18)
$`(i/_xm_f)F_{ij}(x,y)`$ $`=`$ $`{\displaystyle _0^{x^0}}dz𝐀_{ik}(x,z)F_{kj}(z,y)`$ (3.19)
$``$ $`{\displaystyle _0^{y^0}}dz𝐂_{ik}(x,z)\rho _{kj}(z,y),`$
where we have taken the initial time to be zero, and $`_{y^0}^{x^0}dz_{y^0}^{x^0}dz^0\mathrm{d}^3z`$. For known self-energies the equations (3.18) and (3.19) are exact. We note that the form of their RHS is identical to the one for scalar fields . 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 matrix666We 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.. We note that the fermion anticommutation relation or (3.10) uniquely specifies the initial condition for the spectral function:

$$\gamma ^0\rho _{ij}(x,y)|_{x^0=y^0}=i\delta (\stackrel{}{x}\stackrel{}{y})\delta _{ij}.$$
(3.20)
Suitable nonequilibrium initial conditions for the statistical propagator $`F_{ij}(x,y)`$ will be discussed below.

## 4 Lorentz decomposition

It is very useful to decompose the fields $`\rho _{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

$$\rho =\rho _S+i\gamma _5\rho _P+\gamma _\mu \rho _V^\mu +\gamma _\mu \gamma _5\rho _A^\mu +\frac{1}{2}\sigma _{\mu \nu }\rho _T^{\mu \nu },$$
(4.21)
where $`\sigma _{\mu \nu }=\frac{i}{2}[\gamma _\mu ,\gamma _\nu ]`$ and $`\gamma _5=i\gamma ^0\gamma ^1\gamma ^2\gamma ^3`$. For given flavor indices the 16 (pseudo-)scalar, (pseudo-)vector and tensor components

$`\rho _S`$ $`=`$ $`\stackrel{~}{\text{tr}}\rho ,`$
$`\rho _P`$ $`=`$ $`i\stackrel{~}{\text{tr}}\gamma _5\rho ,`$
$`\rho _V^\mu `$ $`=`$ $`\stackrel{~}{\text{tr}}\gamma ^\mu \rho ,`$ (4.22)
$`\rho _A^\mu `$ $`=`$ $`\stackrel{~}{\text{tr}}\gamma _5\gamma ^\mu \rho ,`$
$`\rho _T^{\mu \nu }`$ $`=`$ $`\stackrel{~}{\text{tr}}\sigma ^{\mu \nu }\rho ,`$
are complex two-point functions. Here we have defined $`\stackrel{~}{\text{tr}}\frac{1}{4}\text{tr}`$ where the trace acts in Dirac space. Equivalently, there are 16 complex components for $`F_{ij}`$, $`𝐀_{ij}`$ and $`𝐂_{ij}`$ for given flavor indices $`i,j`$. Using (3.16) and (3.17), one sees that they obey

$$\rho _{ij}^{(\mathrm{\Gamma })}(x,y)=\left(\rho _{ji}^{(\mathrm{\Gamma })}(y,x)\right)^{},F_{ij}^{(\mathrm{\Gamma })}(x,y)=\left(F_{ji}^{(\mathrm{\Gamma })}(y,x)\right)^{},$$
(4.23)
where $`\mathrm{\Gamma }=\{S,P,V,A,T\}`$. Inserting the above decomposition into the evolution equations (3.18) and (3.19) one obtains the respective equations for the various components displayed in Eq. (4.22).

For a more detailed discussion, we first consider the LHS of the evolution equations (3.18) and (3.19). In fact, the approximation of a vanishing RHS corresponds to the standard mean–field or Hartree–type approaches frequently discussed in the literature . 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.18) one finds, omitting flavor indices (see also Ref. ):

$`\stackrel{~}{\text{tr}}\left[(i/m_f)\rho \right]`$ $`=`$ $`\left(i_\mu \rho _V^\mu \right)m_f\rho _S,`$
$`i\stackrel{~}{\text{tr}}\left[\gamma _5(i/m_f)\rho \right]`$ $`=`$ $`i\left(i_\mu \rho _A^\mu \right)m_f\rho _P,`$
$`\stackrel{~}{\text{tr}}\left[\gamma ^\mu (i/m_f)\rho \right]`$ $`=`$ $`\left(i^\mu \rho _S\right)+i\left(i_\nu \rho _T^{\nu \mu }\right)m_f\rho _V^\mu ,`$ (4.24)
$`\stackrel{~}{\text{tr}}\left[\gamma _5\gamma ^\mu (i/m_f)\rho \right]`$ $`=`$ $`i\left(i^\mu \rho _P\right)+{\displaystyle \frac{1}{2}}ϵ^{\mu \nu \gamma \delta }\left(i_\nu \rho _{T,\gamma \delta }\right)m_f\rho _A^\mu ,`$
$`\stackrel{~}{\text{tr}}\left[\sigma ^{\mu \nu }(i/m_f)\rho \right]`$ $`=`$ $`i\left(i^\mu \rho _V^\nu i^\nu \rho _V^\mu \right)+ϵ^{\mu \nu \gamma \delta }\left(i_\gamma \rho _{A,\delta }\right)m_f\rho _T^{\mu \nu }.`$
The corresponding expressions for the LHS of the evolution equation (3.19) for $`F`$ follow from (4.24) with the replacement $`\rho F`$. Considering now the various component (4.22) of the integrand on the RHS of Eq. (3.18), we find

$`\stackrel{~}{\text{tr}}\left[𝐀\rho \right]`$ $`=`$ $`𝐀_S\rho _S𝐀_P\rho _P+𝐀_V^\mu \rho _{V,\mu }𝐀_A^\mu \rho _{A,\mu }`$ (4.25)
$`+{\displaystyle \frac{1}{2}}𝐀_T^{\mu \nu }\rho _{T,\mu \nu },`$
$`i\stackrel{~}{\text{tr}}\left[\gamma _5𝐀\rho \right]`$ $`=`$ $`𝐀_S\rho _P+𝐀_P\rho _Si𝐀_V^\mu \rho _{A,\mu }+i𝐀_A^\mu \rho _{V,\mu }`$ (4.26)
$`+{\displaystyle \frac{1}{4}}ϵ^{\mu \nu \gamma \delta }𝐀_{T,\mu \nu }\rho _{T,\gamma \delta },`$
$`\stackrel{~}{\text{tr}}\left[\gamma ^\mu 𝐀\rho \right]`$ $`=`$ $`𝐀_S\rho _V^\mu +𝐀_V^\mu \rho _Si𝐀_P\rho _A^\mu +i𝐀_A^\mu \rho _P+i𝐀_{V,\nu }\rho _T^{\nu \mu }`$ (4.27)
$`+i𝐀_T^{\mu \nu }\rho _{V,\nu }+{\displaystyle \frac{1}{2}}ϵ^{\mu \nu \gamma \delta }\left(𝐀_{A,\nu }\rho _{T,\gamma \delta }+𝐀_{T,\nu \gamma }\rho _{A,\delta }\right),`$
$`\stackrel{~}{\text{tr}}\left[\gamma _5\gamma ^\mu 𝐀\rho \right]`$ $`=`$ $`𝐀_S\rho _A^\mu +𝐀_A^\mu \rho _Si𝐀_P\rho _V^\mu +i𝐀_V^\mu \rho _P+i𝐀_{A,\nu }\rho _T^{\nu \mu }`$ (4.28)
$`+i𝐀_T^{\mu \nu }\rho _{A,\nu }+{\displaystyle \frac{1}{2}}ϵ^{\mu \nu \gamma \delta }\left(𝐀_{V,\nu }\rho _{T,\gamma \delta }+𝐀_{T,\nu \gamma }\rho _{V,\delta }\right),`$
$`\stackrel{~}{\text{tr}}\left[\sigma ^{\mu \nu }𝐀\rho \right]`$ $`=`$ $`𝐀_S\rho _T^{\mu \nu }+𝐀_T^{\mu \nu }\rho _S{\displaystyle \frac{1}{2}}ϵ^{\mu \nu \gamma \delta }\left(𝐀_P\rho _{T,\gamma \delta }+𝐀_{T,\gamma \delta }\rho _P\right)`$ (4.29)
$`i\left(𝐀_V^\mu \rho _V^\nu 𝐀_V^\nu \rho _V^\mu \right)+ϵ^{\mu \nu \gamma \delta }\left(𝐀_{V,\gamma }\rho _{A,\delta }𝐀_{A,\gamma }\rho _{V,\delta }\right)`$
$`+i\left(𝐀_A^\mu \rho _A^\nu 𝐀_A^\nu \rho _A^\mu \right)+i\left(𝐀_T^{\mu \gamma }\rho _{T,\gamma }^{}{}_{}{}^{\nu }𝐀_T^{\nu \gamma }\rho _{T,\gamma }^{}{}_{}{}^{\mu }\right).`$
With the above expressions one obtains the evolution equations for the various Lorentz components in a straightforward way using (3.18). We note that the convolutions appearing on the RHS of the evolution equation (3.19) for $`F`$ are of the same form than those computed above for $`\rho `$. The respective RHS can be read off Eqs. (4.25)–(4.29) by replacing $`\rho F`$ for the first term and $`𝐀𝐂`$ for the second term under the integrals of Eq. (3.19). 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.

## 5 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

$$\rho (x,y)\rho (x^0,y^0;\stackrel{}{x}\stackrel{}{y})=\frac{\mathrm{d}^3p}{(2\pi )^3}e^{i\stackrel{}{p}(\stackrel{}{x}\stackrel{}{y})}\rho (x^0,y^0;\stackrel{}{p}),$$
(5.30)
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

$`\rho _V^0(x^0,y^0;\stackrel{}{p})`$ $`=`$ $`\rho _V^0(x^0,y^0;p),`$
$`\stackrel{}{\rho }_V(x^0,y^0;\stackrel{}{p})`$ $`=`$ $`\stackrel{}{v}\rho _V(x^0,y^0;p),`$
where $`p|\stackrel{}{p}|`$ and $`\stackrel{}{v}=\stackrel{}{p}/p`$.

Parity: The vector components $`\rho _V^0(x^0,y^0;p)`$ and $`\rho _V(x^0,y^0;p)`$ are unchanged under a parity transformation, whereas the corresponding axial-vector components get a minus sign. Therefore, parity together with rotational invariance imply that

$$\rho _A^0(x^0,y^0;p)=\rho _A(x^0,y^0;p)=0.$$
(5.31)
The same is true for the axial-vector components of $`F`$, $`𝐀`$ and $`𝐂`$. 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 $`\rho `$ transforms as

$`\rho _V^0(x^0,y^0;p)`$ $``$ $`\rho _V^0(y^0,x^0;p),`$
$`\rho _V(x^0,y^0;p)`$ $``$ $`\rho _V(y^0,x^0;p),`$
and similarly for $`𝐀_V^0`$ and $`𝐀_V`$. The $`F`$–components transform as

$`F_V^0(x^0,y^0;p)`$ $``$ $`F_V^0(y^0,x^0;p),`$
$`F_V(x^0,y^0;p)`$ $``$ $`F_V(y^0,x^0;p),`$
and similarly for $`𝐂_V^0`$ and $`𝐂_V`$. Combining this with the hermiticity relations (4.23), one obtains for these components that

$$\begin{array}{ccc}\mathrm{Re}\rho _V^0(x^0,y^0;p)\hfill & =\hfill & \mathrm{Im}\rho _V(x^0,y^0;p)=0,\hfill \\ \mathrm{Re}F_V^0(x^0,y^0;p)\hfill & =\hfill & \mathrm{Im}F_V(x^0,y^0;p)=0,\hfill \\ \mathrm{Re}𝐀_V^0(x^0,y^0;p)\hfill & =\hfill & \mathrm{Im}𝐀_V(x^0,y^0;p)=0,\hfill \\ \mathrm{Re}𝐂_V^0(x^0,y^0;p)\hfill & =\hfill & \mathrm{Im}𝐂_V(x^0,y^0;p)=0,\hfill \end{array}$$
(5.32)
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.21) allowed by chiral symmetry are those which anticommute with $`\gamma _5`$. We therefore have

$$\rho _S(x^0,y^0;\stackrel{}{p})=\rho _P(x^0,y^0;\stackrel{}{p})=\rho _T^{\mu \nu }(x^0,y^0;\stackrel{}{p})=0,$$
(5.33)
and similarly for the corresponding components of $`F`$, $`𝐀`$ and $`𝐂`$. In particular, chiral symmetry forbids a mass term for fermions and we have $`m_f0`$.

### 5.1 Equations of motion

In conclusion, for the above symmetry properties we are left with only four independent propagators: the two spectral functions $`\rho _V^0`$ and $`\rho _V`$ and the two corresponding statistical functions $`F_V^0`$ and $`F_V`$. They are either purely real or imaginary and have definite symmetry properties under the exchange of their time arguments $`x^0y^0`$. These properties as well as the corresponding ones for the various components of the self-energy are summarized below:

| $`\rho _V^0`$, $`𝐀_V^0`$: | imaginary, symmetric; |
| --- | --- |
| $`\rho _V`$, $`𝐀_V`$: | real, antisymmetric; |
| $`F_V^0`$, $`𝐂_V^0`$: | imaginary, antisymmetric; |
| $`F_V`$, $`𝐂_V`$: | real, symmetric. |

The exact evolution equations for the spectral functions read (cf. Eq. (3.18)):777We note that the following equations do not rely on the restrictions (5.32) imposed by $`CP`$–invariance: they have the very same form for the case that all two-point functions are complex.

$`i{\displaystyle \frac{}{x^0}}\rho _V^0(x^0,y^0;p)=p\rho _V(x^0,y^0;p)`$
$`+`$ $`{\displaystyle _{y^0}^{x^0}}dz^0\left[𝐀_V^0(x^0,z^0;p)\rho _V^0(z^0,y^0;p)𝐀_V(x^0,z^0;p)\rho _V(z^0,y^0;p)\right],`$
$`i{\displaystyle \frac{}{x^0}}\rho _V(x^0,y^0;p)=p\rho _V^0(x^0,y^0;p)`$
$`+`$ $`{\displaystyle _{y^0}^{x^0}}dz^0\left[𝐀_V^0(x^0,z^0;p)\rho _V(z^0,y^0;p)𝐀_V(x^0,z^0;p)\rho _V^0(z^0,y^0;p)\right].`$
Similarly, for the statistical two-point functions we obtain (cf. Eq. (3.19)):

$`i{\displaystyle \frac{}{x^0}}F_V^0(x^0,y^0;p)=pF_V(x^0,y^0;p)`$
$`+`$ $`{\displaystyle _0^{x^0}}dz^0\left[𝐀_V^0(x^0,z^0;p)F_V^0(z^0,y^0;p)𝐀_V(x^0,z^0;p)F_V(z^0,y^0;p)\right]`$
$``$ $`{\displaystyle _0^{y^0}}dz^0\left[𝐂_V^0(x^0,z^0;p)\rho _V^0(z^0,y^0;p)𝐂_V(x^0,z^0;p)\rho _V(z^0,y^0;p)\right],`$
$`i{\displaystyle \frac{}{x^0}}F_V(x^0,y^0;p)=pF_V^0(x^0,y^0;p)`$
$`+`$ $`{\displaystyle _0^{x^0}}dz^0\left[𝐀_V^0(x^0,z^0;p)F_V(z^0,y^0;p)𝐀_V(x^0,z^0;p)F_V^0(z^0,y^0;p)\right]`$
$``$ $`{\displaystyle _0^{y^0}}dz^0\left[𝐂_V^0(x^0,z^0;p)\rho _V(z^0,y^0;p)𝐂_V(x^0,z^0;p)\rho _V^0(z^0,y^0;p)\right].`$
The above equations are employed below to calculate the nonequilibrium fermion dynamics in a chiral quark-meson model.

## 6 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 $`\sigma `$–field and a triplet of pseudoscalar “pions” $`\pi ^a`$ ($`a=1,2,3`$). The classical action reads

$`S`$ $`=`$ $`{\displaystyle }\mathrm{d}^4x\{\overline{\psi }i/\psi +{\displaystyle \frac{1}{2}}[_\mu \sigma ^\mu \sigma +_\mu \pi ^a^\mu \pi ^a]`$ (6.38)
$`+g\overline{\psi }[\sigma +i\gamma _5\tau ^a\pi ^a]\psi V(\sigma ^2+\pi ^2)\},`$
where $`\pi ^2\pi ^a\pi ^a`$ and where $`\tau ^a`$ denote the standard Pauli matrices. The above action is invariant under chiral $`SU_L(2)\times SU_R(2)`$ transformations. For a quartic scalar self-interaction $`\left(\sigma ^2+\pi ^a\pi ^a\right)^2`$, this model corresponds to the well known linear $`\sigma `$–model , which has been extensively studied in thermal equilibrium in the literature using various approximations . For simplicity we consider here a purely quadratic scalar potential:

$$V=\frac{1}{2}m_0^2(\sigma ^2+\pi ^2),$$
(6.39)
which is sufficient to study thermalization in this model. We note that this theory has the same universal properties than the corresponding linear $`\sigma `$-model. Extending our study to take into account quartic self-interaction is straightforward. It gives additional contributions to the scalar self-energies, Eqs. (6.50) and (6.2) below. These contributions can be found in Refs.  and we will point out the respective changes below.

### 6.1 Equations of motion for the scalar field

The 2PI effective action for nonequilibrium scalar fields has been extensively studied in the literature . 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. ). For the model considered here, Eq. (6.38), the 2PI effective action is a functional of fermionic propagators as well as scalar propagators.888As emphasized in the previous section, we do not consider the possibility of a broken symmetry in this paper. Therefore we restrict the discussion to a vanishing scalar field expectation value. The scalar fields form an $`O(4)`$-vector $`\varphi _A(x)(\sigma (x),\stackrel{}{\pi }(x))`$ and we denote the full scalar propagator by $`𝒢_{AB}(x,y)`$ with $`A,B=0,\mathrm{},3`$. The 2PI effective action (2.2), augmented by the scalar sector, reads

$$\mathrm{\Gamma }[𝒢,D]=\frac{i}{2}\text{Tr}\mathrm{ln}𝒢^1+\frac{i}{2}\text{Tr}𝒢_0^1𝒢i\text{Tr}\mathrm{ln}D^1i\text{Tr}D_0^1D+\mathrm{\Gamma }_2[𝒢,D]+\mathrm{const},$$
(6.40)
with the free scalar inverse propagator

$$i𝒢_{0,AB}^1(x,y)=(\mathrm{}_x+m_0^2)\delta _𝒞^4(xy)\delta _{AB}.$$
(6.41)
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

$$\mathrm{\Sigma }_{AB}(x,y)=2i\frac{\delta \mathrm{\Gamma }_2[𝒢,D]}{\delta 𝒢_{BA}(y,x)}.$$
(6.42)
Under chiral transformations, the matrix $`𝒢𝒢^{}`$, where $``$ is an $`O(4)`$ rotation. Without loss of generality, because of chiral symmetry the effective action $`\mathrm{\Gamma }[𝒢,D]`$ and its functional derivatives can be evaluated for $`𝒢`$ taken to be the unit matrix in $`O(4)`$–space:

$$𝒢_{AB}(x,y)=G_\varphi (x,y)\delta _{AB}.$$
(6.43)
The same holds for the corresponding self-energy (6.42). Similarly, in the fermionic sector, because of chiral symmetry the most general fermion two-point function can be taken to be proportional to unity in flavor space:

$$D_{ij}(x,y)=D(x,y)\delta _{ij}.$$
(6.44)
Similar to the discussion in Sect. 3, the scalar spectral and statistical two-point functions are defined as

$$G_\varphi (x,y)=F_\varphi (x,y)\frac{i}{2}\rho _\varphi (x,y)\left[\mathrm{\Theta }_𝒞(x^0y^0)\mathrm{\Theta }_𝒞(y^0x^0)\right],$$
(6.45)
and equivalently for the spectral and statistical self-energies $`\mathrm{\Sigma }_\varphi ^\rho `$ and $`\mathrm{\Sigma }_\varphi ^F`$. These are all real functions and $`F`$–like components are symmetric under the exchange of $`x`$ and $`y`$, whereas the $`\rho `$–like components are antisymmetric. The equal-time commutation relation of two scalar field operators implies

$$\rho _\varphi (x,y)|_{x^0=y^0}=0,_{x^0}\rho _\varphi (x,y)|_{x^0=y^0}=\delta (\stackrel{}{x}\stackrel{}{y}),$$
(6.46)
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 :

$`\left[_{x^0}^2+\stackrel{}{p}^2+m_0^2\right]\rho _\varphi (x^0,y^0;\stackrel{}{p})={\displaystyle _{y^0}^{x^0}}dz^0\mathrm{\Sigma }_\varphi ^\rho (x^0,z^0;\stackrel{}{p})\rho _\varphi (z^0,y^0;\stackrel{}{p}).`$
(6.47)
$`\left[_{x^0}^2+\stackrel{}{p}^2+m_0^2\right]F_\varphi (x^0,y^0;\stackrel{}{p})={\displaystyle _0^{x^0}}dz^0\mathrm{\Sigma }_\varphi ^\rho (x^0,z^0;\stackrel{}{p})F_\varphi (z^0,y^0;\stackrel{}{p})`$
$`+{\displaystyle _0^{y^0}}dz^0\mathrm{\Sigma }_\varphi ^F(x^0,z^0;\stackrel{}{p})\rho _\varphi (z^0,y^0;\stackrel{}{p}),`$
(6.48)
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.38) together with (6.39). In the presence of scalar self-interactions, the self-energy receives in particular a local contribution $`i\mathrm{\Sigma }_\varphi ^{(\mathrm{local})}(x)\delta _𝒞^4(xy)`$, 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 $`m_0^2M^2(x)=m_0^2+\mathrm{\Sigma }_\varphi ^{(\mathrm{local})}(x)`$ .

### 6.2 2PI coupling expansion

We consider a systematic coupling expansion of the 2PI effective action (6.40) to lowest non-trivial order, which includes scattering as well as memory and off-shell effects, such that thermalization can be studied. The first non-trivial order in a coupling expansion corresponds to the two-loop contribution depicted in Fig. 1 (cf. also the discussion in Ref. ).999We note that this approximation can also be related to a nonperturbative $`1/N_f`$ expansion at next-to-leading order. Using (6.43) and (6.44), one can express the two-loop contribution to $`\mathrm{\Gamma }_2`$ directly in terms of $`G_\varphi `$ and $`D`$. We obtain for the chirally symmetric theory:

$$\mathrm{\Gamma }_2^{(2\mathrm{loop})}[𝒢,D]=ig^2\frac{N_fN_s}{2}_𝒞\mathrm{d}^4x\mathrm{d}^4y\text{tr}[D(x,y)D(y,x)]G_\varphi (x,y),$$
(6.49)
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-energies101010The relevant self-energies can be obtained with
$$2i\frac{\delta \mathrm{\Gamma }_2}{\delta G_\varphi }=2i\frac{\delta \mathrm{\Gamma }_2}{\delta 𝒢_{AB}}\frac{\delta 𝒢_{AB}}{\delta G_\varphi }=\delta _{AB}\mathrm{\Sigma }_{BA}=N_s\mathrm{\Sigma }_\varphi ,$$
and similarly for the fermionic self energies $`\mathrm{\Sigma }_{ij}=\mathrm{\Sigma }\delta _{ij}`$:
$$i\frac{\delta \mathrm{\Gamma }_2}{\delta D}=N_f\mathrm{\Sigma }.$$
. We obtain for the scalar self-energies

$`\mathrm{\Sigma }_\varphi ^\rho (x^0,y^0;\stackrel{}{p})=8g^2N_f{\displaystyle \frac{\mathrm{d}^3q}{(2\pi )^3}\rho _V^\mu (x^0,y^0;\stackrel{}{q})F_{V,\mu }(x^0,y^0;\stackrel{}{p}\stackrel{}{q})},`$ (6.50)
$`\mathrm{\Sigma }_\varphi ^F(x^0,y^0;\stackrel{}{p})=4g^2N_f{\displaystyle }{\displaystyle \frac{\mathrm{d}^3q}{(2\pi )^3}}[F_V^\mu (x^0,y^0;\stackrel{}{q})F_{V,\mu }(x^0,y^0;\stackrel{}{p}\stackrel{}{q})`$
$`{\displaystyle \frac{1}{4}}\rho _V^\mu (x^0,y^0;\stackrel{}{q})\rho _{V,\mu }(x^0,y^0;\stackrel{}{p}\stackrel{}{q})],`$ (6.51)
and for the fermion self-energies

$`𝐀_V^\mu (x^0,y^0;\stackrel{}{p})=g^2N_s{\displaystyle }{\displaystyle \frac{\mathrm{d}^3q}{(2\pi )^3}}[F_V^\mu (x^0,y^0;\stackrel{}{q})\rho _\varphi (x^0,y^0;\stackrel{}{p}\stackrel{}{q})`$
$`+\rho _V^\mu (x^0,y^0;\stackrel{}{q})F_\varphi (x^0,y^0;\stackrel{}{p}\stackrel{}{q})],`$ (6.52)
$`𝐂_V^\mu (x^0,y^0;\stackrel{}{p})=g^2N_s{\displaystyle }{\displaystyle \frac{\mathrm{d}^3q}{(2\pi )^3}}[F_V^\mu (x^0,y^0;\stackrel{}{q})F_\varphi (x^0,y^0;\stackrel{}{p}\stackrel{}{q})`$
$`{\displaystyle \frac{1}{4}}\rho _V^\mu (x^0,y^0;\stackrel{}{q})\rho _\varphi (x^0,y^0;\stackrel{}{p}\stackrel{}{q})].`$ (6.53)
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.47) and (6.48) relevant for thermalization.

## 7 Initial conditions

The time evolution for the fermions is described by first-order (integro-)differential equations for $`F`$ and $`\rho `$: Eqs. (5.1)–(5.1). 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.20)). 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

$`F_V(t,t^{},p)|_{t=t^{}=0}`$ $`=`$ $`{\displaystyle \frac{1}{2}}n_0^f(p),`$ (7.54)
$`F_V^0(t,t^{},p)|_{t=t^{}=0}`$ $`=`$ $`0.`$ (7.55)
Here $`n_0^f(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.111111In particular, for nonzero chemical potential or net charge density the BCS mechanism can lead to the condensation of Cooper pairs of fermions, which will be discussed elsewhere.

The evolution equations (6.48) and (6.47) 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.46)). For the scalar two-point function $`F_\varphi `$ we consider (cf. also Refs. )

$`F_\varphi (t,t^{},p)|_{t=t^{}=0}`$ $`=`$ $`{\displaystyle \frac{1}{ϵ_0(p)}}\left[n_0(p)+{\displaystyle \frac{1}{2}}\right],`$
$`_tF_\varphi (t,t^{},p)|_{t=t^{}=0}`$ $`=`$ $`0,`$ (7.56)
$`_t_t^{}F_\varphi (t,t^{},p)|_{t=t^{}=0}`$ $`=`$ $`ϵ_0(p)\left[n_0(p)+{\displaystyle \frac{1}{2}}\right],`$
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.1)–(5.1) and (5.1)–(5.1) by neglecting the memory integrals on their RHS. The solution of the fermionic free field equations with the above initial conditions reads

$`F_V(t,t^{},p)`$ $`=`$ $`\left({\displaystyle \frac{1}{2}}n_0^f(p)\right)\mathrm{cos}[p(tt^{})],`$
$`F_V^0(t,t^{},p)`$ $`=`$ $`i\left({\displaystyle \frac{1}{2}}n_0^f(p)\right)\mathrm{sin}[p(tt^{})]`$
and for the spectral functions one obtains

$$\rho _V(t,t^{},p)=\mathrm{sin}[p(tt^{})],\rho _V^0(t,t^{},p)=i\mathrm{cos}[p(tt^{})].$$

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.1)-(5.1). 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. . 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.

## 8 Numerical implementation

We numerically solve the evolution equations (5.1)–(5.1) and (6.47)–(6.48), together with the self-energies Eqs. (6.50)–(6.53). 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_V^0(t,t^{})`$ is analogous to the canonical momentum. This suggests to discretize $`F_V(t,t^{})`$ and $`\rho _V(t,t^{})`$ at $`tt^{}=2na_t`$ (even) and $`F_V^0(t,t^{})`$ and $`\rho _V^0(t,t^{})`$ at $`tt^{}=(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.1)–(5.1) as well.

We emphasize that the discretization does not suffer from the problem of so–called fermion doublers . The spatial doublers do not appear since (5.1)–(5.1) are effectively second order in $`\stackrel{}{x}`$-space. Writing the equations for $`\stackrel{}{F}_V(t,t^{},\stackrel{}{x})`$ and $`\stackrel{}{\rho }_V(t,t^{},\stackrel{}{x})`$ starting from (5.1)–(5.1) 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, time-like doublers are easily avoided by using a sufficiently small stepsize in time $`a_t/a_s`$.

The fact that Eqs. (5.1)–(5.1) and (6.47)–(6.48) 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\times 470\times 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.

## 9 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.54)–(7.56)). 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 late-time 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 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_\varphi (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 $`tt^{}`$. We note that the oscillation envelope of $`F_V(t,t^{},p)`$ 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 $`tt^{}`$ for the late time $`mt=600`$. The fit to an exponential yields the damping rate $`\gamma _f^{(\mathrm{damp})}(p=0.78m)=0.03(1)m`$. We find a moderate momentum dependence of this rate with $`\gamma _f^{(\mathrm{damp})}(0)=0.067(1)m`$ and $`\gamma _f^{(\mathrm{damp})}(0)\gamma _f^{(\mathrm{damp})}(p>0)`$.

We emphasize that the rate $`\gamma _f^{(\mathrm{damp})}(p)`$ can be related to the width of the Fourier transform of the spectral function with respect to the time difference $`tt^{}`$. 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 $`\rho _V(t,t^{},p)`$ only on a finite interval. To overcome this problem, we fit the data for $`\rho _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 $`\rho _V(\omega ,p)`$ is displayed as a function of frequency $`\omega `$ 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 $`p0`$ $`\left(1/\gamma _f^{(\mathrm{damp})}=15(1)m^1\right)`$ 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/\gamma _f^{(\mathrm{therm})}=95(5)m^1`$. The corresponding time scale $`1/\gamma _f^{(\mathrm{therm})}`$ is therefore much larger than the characteristic damping time $`1/\gamma _f^{(\mathrm{damp})}`$.

A similar analysis can be performed for the scalar sector. Here we find at $`p0`$ the corresponding values $`1/\gamma _\varphi ^{(\mathrm{damp})}=50(5)m^1`$ and $`1/\gamma _\varphi ^{(\mathrm{therm})}=90(5)m^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  and classical  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_\varphi (t,0)`$, $`F_V^0(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 late-time value for the fermions, and about half the late-time rate for the scalars.

## 10 Approach to quantum thermal equilibrium

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 approached121212We emphasize that thermal equilibrium cannot be reached on a fundamental level from time-reversal invariant evolution equations at any finite time. The results demonstrate that thermal equilibrium can be approached very closely at sufficiently late time, without again deviating from it for practically accessible times. For a more detailed discussion of this aspect see Ref. .. 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

$$F_\varphi ^{(\mathrm{eq})}(\omega ,\stackrel{}{p})=i\left(n^{BE}(\omega )+\frac{1}{2}\right)\rho _\varphi ^{(\mathrm{eq})}(\omega ,\stackrel{}{p}),$$
(10.57)
for the scalar correlators, where $`n^{BE}(\omega )=1/[\mathrm{exp}(\omega /T)1]`$ denotes the Bose–Einstein distribution function. The frequency $`\omega `$ is the Fourier conjugate of $`tt^{}`$ (in thermal equilibrium, time-translation invariance implies that two-point functions only depend on $`tt^{}`$). 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}(\omega )+\frac{1}{2})(n^{FD}(\omega )\frac{1}{2})`$ in Eq. (10.57) and the Fermi–Dirac distribution $`n^{FD}=1/[\mathrm{exp}(\omega /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_\varphi (t,t^{};p)`$ and $`\rho _\varphi (t,t^{};p)`$, as well as the corresponding fermion two-point functions, in terms of the center coordinate $`X^0=(t+t^{})/2`$ and the relative coordinate $`s^0=tt^{}`$ and write

$$\rho _\varphi (X^0;\omega ,p)=_{2X^0}^{2X^0}ds^0e^{i\omega s^0}\rho _\varphi (X^0+s^0/2,X^0s^0/2;p),$$
(10.58)
and equivalently for the other correlators. Since we consider an initial-value problem, the time integral over $`s^0`$ is bounded by $`\pm 2X^0`$ (cf. also the detailed discussion in Ref. ). 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 $`tt^{}`$ can be performed to very good approximation (cf. the discussion in Sect. 9). The distribution functions may then be 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^{FD}(\omega )`$ and $`n^{BE}(\omega )`$, respectively (cf. Eq. (10.57)). 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 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 quasi-particle picture has been successfully employed previously to investigate thermalization .

In Figs. 6 and 7 we show the effective quasi-particle number distributions defined as

$$\frac{1}{2}n_f(t,p)=F_V(t,t;p)$$
(10.59)
for fermions and

$`{\displaystyle \frac{1}{2}}+n(t,p)`$ $`=`$ $`ϵ(t,p)F_\varphi (t,t;p),`$ (10.60)
$`ϵ(t,p)`$ $`=`$ $`\left({\displaystyle \frac{_t_t^{}F_\varphi (t,t^{};p)}{F_\varphi (t,t^{};p)}}\right)_{t=t^{}}^{1/2}`$ (10.61)
for scalars, where $`\epsilon (t,p)`$ is the quasi-particle mode energy as discussed in Appendix A.131313Note that because of chiral symmetry there is no mass term present for the fermions and the corresponding quasi-particle mode energy is simply $`p`$. 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” $`\mathrm{log}(1/n_f1)`$ and $`\mathrm{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.

## 11  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 , this provides a description of the nonequilibrium dynamics of the linear $`\sigma `$-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.

### Acknowledgment

We thank C. Wetterich for many discussions and collaboration on related work. We are also grateful to W. Wetzel for his continuous support with computers. Sz. B. acknowledges the hospitality of the Institute für Theoretische Physik, Heidelberg. His work was supported by the short-term scholarship program of the DAAD.

## Appendix A Effective particle number and energy

A particle number can only be strictly defined in the presence of a conserved charge. However, for physical interpretation it is often convenient to define an effective particle number even when there is no conserved charge. In particular, besides a total particle number it is often useful to have a definition of an effective particle number per (momentum) mode. The latter is typically not conserved in an interacting theory, and in the context of thermalization one would like to find a time-dependent particle number distribution, which allows one to observe a Bose-Einstein or Fermi-Dirac distribution at sufficiently late times141414We emphasize, however, that the definition of an effective particle number distribution is not necessary to analyse the approach to quantum thermal equilibrium, as is discussed in Sect. 10.. Of course, the notion of an effective particle number is typically only meaningful if the relevant degrees of freedom or “quasi-particles” are weakly interacting.

There are many different ways of how one can introduce the notion of an effective particle number (for a recent discussion see e.g. Ref.  and references therein). For example, one can define it as the average energy per mode divided by the energy of the corresponding mode. In an interacting theory, this procedure can be ambiguous because the expression of the total energy receives contributions from interactions and the average energy per mode is not uniquely defined. In this appendix, we discuss an elegant way to circumvent this difficulty. More importantly for our purposes, the procedure leads to a definition of an effective particle number density, which indeed allows one to directly observe the emergence of a Fermi-Dirac or Bose-Einstein distribution from the nonequilibrium dynamics for the theory considered in this paper. In particular, it can be explicitely shown that the effective particle numbers are always positive and, for the fermions, smaller than one (see Appendix B). Applying the present construction to the case of neutral scalar fields, we recover the particle number definition used in previous studies to exhibit thermalization in purely scalar field theories .

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.151515As 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. The associated $`4`$–current for each given flavor is $`\overline{\psi }\gamma ^\mu \psi `$. Fourier transformed with respect to spatial momenta, the expectation value of the latter can be written as $`J_f^\mu (t,p)=\text{tr}[\gamma ^\mu D^<(t,t,p)]`$, where the subscript stands for fermions (cf. Eq. (3.8)). In terms of the equal-time statistical two point function, its temporal and spatial components read:161616The constant factor in the temporal component comes from the fact that we define the current without the standard normal ordering.

$`J_f^0(t,p)`$ $`=`$ $`2[12F_V^0(t,t;p)],`$
$`\stackrel{}{J}_f(t,p)`$ $`=`$ $`4\stackrel{}{v}F_V(t,t;p).`$
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:171717Here, 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.

$`J_f^{0(\mathrm{QP})}(t,p)`$ $`=`$ $`2[1+Q_f(t,p)],`$
$`\stackrel{}{J}_f^{(\mathrm{QP})}(t,p)`$ $`=`$ $`2\stackrel{}{v}[12N_f(t,p)],`$
where $`Q_f(t,p)=n_f\overline{n}_f`$ is the difference between particle and anti-particle effective number densities and $`N_f(t,p)=(n_f+\overline{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 $`\stackrel{}{J}`$ is the net current density per mode and is therefore sensitive to the sum of particle and anti-particle number densities. Identifying the above expressions, we define

$`{\displaystyle \frac{1}{2}}Q_f(t,p)`$ $`=`$ $`F_V^0(t,t;p),`$ (A.62)
$`{\displaystyle \frac{1}{2}}N_f(t,p)`$ $`=`$ $`F_V(t,t;p).`$ (A.63)
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.32)). Moreover, using the anticommutation relation for the fermion fields, it is shown in Appendix B that these definitions always satisfy $`0N_f(t,p)1`$ and $`1Q_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.32) 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 anti-particle numbers are equal and we have

$`Q_f(t,p)`$ $`=`$ $`0,`$
$`{\displaystyle \frac{1}{2}}n_f(t,p)`$ $`=`$ $`F_V(t,t;p).`$
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 $`\varphi `$. The $`4`$-current associated the corresponding $`U(1)`$ symmetry is $`i[\varphi ^{}(^\mu \varphi )(^\mu \varphi ^{})\varphi ]`$ 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:

$`J_b^0(t,p)`$ $`=`$ $`i(_t_t^{})F_\varphi (t,t^{};p)|_{t^{}=t}1,`$
$`\stackrel{}{J}_b(t,p)`$ $`=`$ $`2\stackrel{}{p}F_\varphi (t,t;p),`$
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: $`F_\varphi (t,t^{};p)=\frac{1}{2}[\varphi (t,\stackrel{}{p}),\varphi ^{}(t^{},\stackrel{}{p})]_+`$. It has the symmetry property $`F_\varphi (t,t^{};p)=F_\varphi ^{}(t^{},t;p)`$, so that the above expressions are real. The corresponding quasi-particle expressions read:

$`J_b^{0(\mathrm{QP})}(t,p)`$ $`=`$ $`Q_b(t,p)1,`$
$`\stackrel{}{J}_b^{(\mathrm{QP})}(t,p)`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{p}}{ϵ(t,p)}}[1+2N_b(t,p)],`$
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:

$`Q_b(t,p)`$ $`=`$ $`i(_t_t^{})F_\varphi (t,t^{};p)|_{t^{}=t},`$ (A.64)
$`1+2N_b(t,p)`$ $`=`$ $`2ϵ(t,p)F_\varphi (t,t;p).`$ (A.65)
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)`$.181818Note that this was not necessary for the fermionic case, because of our assumption of chiral symmetry, which prevents an effective mass term. In fact, one can repeat the following argument to obtain for the fermion effective quasi-particle energy $`ϵ_f(t,p)=p`$, as employed above. 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_\varphi (t,t^{};p)|_{t^{}=t}+ϵ^2(t,p)F_\varphi (t,t;p)ϵ(t,p)\left[2n(t,p)+1\right]$$

in terms of the statistical two-point function. Therefore, one obtains for the effective quasi-particle energy:

$$ϵ^2(t,p)=\left(\frac{_t_t^{}F_\varphi (t,t^{};p)}{F_\varphi (t,t^{};p)}\right)_{t=t^{}}.$$
(A.66)
It is straightforward to show that the combination of equal-time correlators appearing on the RHS of (A.66) 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.65) and (A.66), is a positive quantity (see Appendix B). We note that for the case of $`CP`$–invariant systems, the effective charge density per mode (A.64) vanishes, as it should.191919This immediately follows from the behavior of the statistical propagator under $`CP`$–transformation: $`F_\varphi (t,t^{};p)F_\varphi (t^{},t;p)`$.

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_\varphi (t,t^{};p)=F_\varphi (t^{},t;p)`$, which directly implies that $`Q_b(t,p)=0`$. Therefore, in the present paper, we use the definition

$$\frac{1}{2}+n(t,p)=ϵ(t,p)F_\varphi (t,t;p)$$
(A.67)
together with (A.66) for the boson effective particle number $`n(t,p)`$ and mode energy $`ϵ(t,p)`$ for each individual scalar species .

## 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 $`\phi `$ and $`\phi ^{}`$ satisfying the anticommutation relations

$$[\phi ,\phi ^{}]___+=1$$

and

$$[\phi ,\phi ]__+=[\phi ^{},\phi ^{}]__+=0,$$

one has

$$0\phi ^{}\phi \mathrm{\hspace{0.17em}1},$$
(B.68)
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 $`\phi ^{}`$ and $`\phi `$. One obtains a sum of positive quantities which is of course positive. To show the second inequality, we introduce the hermitian operator $`N=\phi ^{}\phi `$. Using the anticommutation relations above, it is easy to see that $`N^2=N`$, from which it follows that

$$\mathrm{\Delta }n^2\left(NN\right)^2=N^2N^2=N\left(1N\right).$$

It is clear that $`\mathrm{\Delta }n^20`$ and we have just shown that $`N0`$. One therefore conclude from the above equation that $`N1`$, as announced.

We now come to our effective particle number densities. Using Eqs. (A.62) and (A.63) and recalling that $`Q_fn_f\overline{n}_f`$ and $`N_f(n_f+\overline{n}_f)/2`$, we get for the effective particle and anti-particle number densities:

$`{\displaystyle \frac{1}{2}}n_f(t,p)`$ $`=`$ $`F_V(t,t,p)+F_V^0(t,t,p)={\displaystyle \frac{1}{4}}\text{tr}\left[(\gamma _0+\stackrel{}{v}\stackrel{}{\gamma })F(t,t,p)\right],`$
$`{\displaystyle \frac{1}{2}}\overline{n}_f(t,p)`$ $`=`$ $`F_V(t,t,p)F_V^0(t,t,p)={\displaystyle \frac{1}{4}}\text{tr}\left[(\gamma _0\stackrel{}{v}\stackrel{}{\gamma })F(t,t,p)\right].`$
In terms of the fermionic field operators $`\psi (t,\stackrel{}{p})`$ and $`\overline{\psi }(t,\stackrel{}{p})`$, which satisfy the anticommutation relations

$$[\psi (t,\stackrel{}{p}),\overline{\psi }(t,\stackrel{}{p})]__+=\gamma _0$$

and

$$[\psi (t,\stackrel{}{p}),\psi (t,\stackrel{}{p})]__+=[\overline{\psi }(t,\stackrel{}{p}),\overline{\psi }(t,\stackrel{}{p})]__+=0,$$

one has

$$F(t,t,p)=\frac{1}{2}[\psi (t,\stackrel{}{p}),\overline{\psi }(t,\stackrel{}{p})]_{_{}}.$$

Therefore, one can rewrite (from now on, we drop the explicit $`t`$ and $`\stackrel{}{p}`$ dependence):

$`n_f`$ $`=`$ $`{\displaystyle \frac{1}{4}}\overline{\psi }(\gamma _0+\stackrel{}{v}\stackrel{}{\gamma })\psi ,`$
$`1\overline{n}_f`$ $`=`$ $`{\displaystyle \frac{1}{4}}\overline{\psi }(\gamma _0\stackrel{}{v}\stackrel{}{\gamma })\psi .`$
Now we introduce the operator

$$\phi =\frac{\gamma _0+\stackrel{}{v}\stackrel{}{\gamma }}{\sqrt{2}}\psi ,$$

in term of which,

$$n_f=\frac{1}{4}\underset{\alpha =1}{\overset{4}{}}\phi _\alpha ^{}\phi _\alpha ,$$
(B.69)
where we explicitely wrote the sum over Dirac indices. For a given Dirac index, it is straightforward to check that

$$[\phi _\alpha ,\phi _\alpha ^{}]___+=\left(1\gamma _0\stackrel{}{v}\stackrel{}{\gamma }\right)_{\alpha \alpha }=1$$
(B.70)
and

$$[\phi _\alpha ,\phi _\alpha ]__+=[\phi _\alpha ^{},\phi _\alpha ^{}]__+=0,$$

where we specialized to the Dirac basis to write the last equality of Eq. (B.70).202020This is more convenient, but not necessary. It is simple to adapt the argument to any basis. Using (B.68) for each individual $`\alpha `$, we conclude from (B.69) that

$$0n_f(t,p)1$$

for any time $`t`$ and any momentum $`p`$. It is straightforward to repeat the above arguments to show that

$$0\overline{n}_f(t,p)1.$$

Bosons

For scalars, we first use Eqs. (A.66) and (A.66) to rewrite

$$\frac{1}{2}+n(t,p)=\sqrt{F_\varphi (t,t;p)[_t_t^{}F_\varphi (t,t^{};p)]_{t^{}=t}}.$$

Recalling the definition of the statistical propagator in terms of field operators:

$`F_\varphi (t,t;p)`$ $`=`$ $`\varphi ^{}(t,\stackrel{}{p})\varphi (t,\stackrel{}{p}),`$
$`_t_t^{}F_\varphi (t,t^{};p)]_{t^{}=t}`$ $`=`$ $`\mathrm{\Pi }^{}(t,\stackrel{}{p})\mathrm{\Pi }(t,\stackrel{}{p}),`$
where $`\mathrm{\Pi }(t,\stackrel{}{p})`$ and $`\mathrm{\Pi }^{}(t,\stackrel{}{p})`$ are the canonical momenta conjugate to $`\varphi ^{}(t,\stackrel{}{p})`$ and $`\varphi (t,\stackrel{}{p})`$ respectively (e.g. $`\mathrm{\Pi }(t,\stackrel{}{p})=\dot{\varphi }^{}(t,\stackrel{}{p})`$), we see that

$$n(t,p)0<=>\varphi ^{}\varphi \mathrm{\Pi }^{}\mathrm{\Pi }\frac{1}{2},$$

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 .