Thermalization and dynamical spectral properties in the quark-meson model

We investigate the nonequilibrium evolution of the quark-meson model using 2PI effective action techniques. Our numerical simulations, which include the full dynamics of the order parameter of chiral symmetry, show how the model thermalizes into different regions of its phase diagram. In particular, by studying quark and meson spectral functions, we shed light on the real-time dynamics approaching the crossover transition, revealing e.g. the emergence of light effective fermionic degrees of freedom in the infrared. At late times in the evolution, the fluctuation-dissipation relation emerges naturally among both meson and quark degrees of freedom, confirming that the simulation successfully reaches the universal thermal fixed point.


I. INTRODUCTION
The quest to discover the conjectured critical point of the QCD phase diagram is a central motivation of modern heavy-ion collision (HIC) experiments at collider facilities, such as the Large Hadron Collider (LHC) at CERN and the Relativistic Heavy-Ion Collider (RHIC) at Brookhaven National Laboratory. In the beam energy scan currently executed at RHIC, the phase diagram of QCD is explored over a wide range of temperatures and baryon densities by depositing different amounts of energy in the initial collision volume. As the fireball expands and cools, the efficient exchange of energy and momentum among quarks and gluons leads to local thermalization over time. The question to answer is: if a critical point exist and some of the volume of the fireball evolves close to it, does the dynamical build-up of long range fluctuations leave any discernible mark on the yields of measurable particles?
Understanding the out-of-equilibrium dynamics of heavy-ion collisions thus remains one of the most pressing theory challenges in heavy-ion physics. So far, genuinely nonperturbative ab-initio calculations of the equilibration process of the quark-gluon plasma and the dynamics close to the phase transition remain out of reach. In order to make progress, we therefore set out to shed light onto pertinent aspects of the physics of dynamical thermalization in heavy-ion collisions by deploying a low-energy effective theory of QCD, the two-flavor quark-meson model. This model incorporates the offshell dynamics of the lowest mass states in QCD, the pseudoscalar pions and the scalar sigma-mode as well as the light up and down quarks. Further degrees of freedom, and in particular the gluons, heavier quark flavors as well as higher mass hadronic resonances carry masses 500 MeV and are neglected here. This low energy effective theory reflects the central and physically relevant feature of low-energy QCD: chiral symmetry breaking in vacuum and its restoration at finite temperature and density. At its critical endpoint, the model is expected to lie in the same universality class as QCD and hence con-stitutes a viable low-energy effective theory to explore dynamical critical phenomena in QCD at finite temperature and density at scales 500 MeV.
In the present work we consider the real-time dynamics of the two-flavor quark-meson model with small current quark masses in a nonexpanding scenario, for progress on the out-of-equilibrium quark-meson model see [1][2][3][4]. In the presence of such an explicit chiral symmetry breaking the equilibrium chiral transition at finite temperature is a crossover as confirmed for QCD at vanishing and small density, for recent results see [5][6][7]. By the help of different initial conditions defined via the initial occupations of sigma and quark fields, we map out the thermalization dynamics for different regions of the phase diagram. This allows, for the first time, to study full the thermalization dynamics including that of order parameter of chiral symmetry. An extension of the present study to the scenario of an expanding fireball should give access to the freeze-out physics of heavy-ion collisions.
The evolution towards thermal equilibrium is viewed through the lens of the one-and two-point functions of the theory, which are computed with the two-particle irreducible (2PI) approach by means of their quantum equations of motion. These correlation functions not only provide complementary order parameters for the study of chiral symmetry restoration but also give direct access to the spectral properties, including the quasiparticle content of the system. Being genuine nonequilibrium quantities they map out the whole time evolution of the system including the physics of the crossover transition in the late-time limit.
This paper is organized as follows. In Section II we briefly review the quark-meson model and give an overview over our nonequilibrium and nonperturbative treatment. The numerical setup for the time evolution starting from free-field initial conditions quenched to a highly nonequilibrium environment is described. In Section III we discuss the spectral functions of the bosonic and fermionic degrees of freedom, which provide information about the masses as well as the lifetimes of the dynamical degrees of freedom. We investigate the late-time limit of our simulations, which reveals the dy-namical emergence of the fluctuation-dissipation relation and hence allows us to define a thermalization temperature. Finally, Section IV covers the results for the sigma field describing the order parameter of the quark-meson model. We further discuss the behavior of different order parameters in equilibrium which lead to a consistent pseudocritical temperature. In Section V we conclude with a summary. Appendix A provides details about the evolution equations of the model including the relevant expressions for the deployed approximation scheme.

II. THE QUARK-MESON MODEL
QCD evolves from a theory of dynamical quarks and gluons at large momentum scales, the fundamental degrees of freedom, to a theory of dynamical hadrons at low momentum scales. This transition of the dynamical degrees of freedom is related to the mass gaps of the respective fields. It is by now well understood that the gluon degrees of freedom start to decouple at about 1 GeV, that is above the chiral symmetry breaking scale k χ of about 400 MeV. Most of the hadron resonances are too heavy for taking part in the offshell dynamics and we are left with the up, down and to some extend the strange quarks, as well as the pions and the scalar sigma mode, for details see [8,9]. Indeed, low-energy effective theories emerge naturally at low momentum scales from first principle QCD, and their systematic embedding leads us to the quark-meson model and its Polyakov loop enhanced version as QCD-assisted low-energy effec-tive theories. While its quantitative validity has been proven for momentum scales k with k 300 MeV [10], it reproduces qualitative QCD features up to k 700 MeV. It is this natural QCD embedding as well as its robust QCD-type chiral properties that has triggered a plethora of works with the quark-meson model on the QCD phase structure with functional methods, see e.g. [11][12][13][14][15][16]. More recently also real-time correlation functions in equilibrium have been investigated in e.g. [17][18][19][20][21][22][23][24][25][26][27].
(Pre)-Thermalization has been studied in the O(N = 4) symmetric scalar model coupled to fermions using a 1/N expansion to next-to-leading order of the 2PI effective action in [1,28]. The model was studied extensively in Refs. [2,3] in the context of inflaton dynamics to describe nonequilibrium instabilities with fermion production from inflaton decay. In [4] the model was investigated for highly occupied bosonic fields, where the predictions were shown to agree well with lattice simulation results in the classical-statistical regime. Further results for spectral functions in and out-of-equilibrium with 2PI effective action techniques can be found in [29], and with classical-statistical simulations in [30,31] for scalar theories, and in [32] for Yang-Mills theory.
The π mesons play the role of the light Goldstone bosons in the chirally broken phase whereas the σ meson represents the heavy mode. Assigning these roles to the components of the scalar field is achieved by choosing a coordinate system in field space where the field expec-tation value has a single component which defines the σ direction, i.e. φ a (x) = ϕ(x) = { σ(x) , 0, 0, 0}.
The quasiparticle excitation spectrum of the quarkmeson model is encoded in the spectral functions of the respective fields. For the bosonic and fermionic fields, the spectral function is defined as the expectation value of the commutator and anticommutator, respectively, where a, b = 1, . . . , N denote field space and A, B = 1, . . . , 4 correspond to Dirac spinor indices. Fermion flavor indices are omitted and the operator nature of the quantum fields is implied. We consider systems with spatial isotropy and homogeneity such that the spectral functions depend on times and relative spatial coordinates, i.e. ρ(t, t , |x − y|) or in momentum space ρ(t, t , |p|), while the field expectation value only depends on time, i.e. σ(t) . Due to the remaining O(N − 1) symmetry of the chirally broken model, the bosonic spectral function can be written as ρ φ ab = diag(ρ σ , ρ π , ρ π , ρ π ) where the components ρ i with i = σ, π describe the respective mesons. The fermionic spectral function can be decomposed into Lorentz components according to with σ µν = i 2 [γ µ , γ ν ] and γ 5 = iγ 0 γ 1 γ 2 γ 3 . The corresponding Lorentz components are given by where the trace acts in Dirac space. In spatially homogeneous and isotropic systems with parity and CP invariance, the only nonvanishing components are the scalar, vector and 0i-tensor components. Rotational invariance allows us to write where we refer to the two-point functions ρ S , ρ 0 , ρ V and ρ T on the right hand sides as the scalar, vector, vectorzero and tensor components. The relevant contributions to the quark spectral function are the scalar, vector-zero and vector components, where the vector-zero component represents the quark excitations of the system [25,33]. For chiral symmetric theories with m ψ = 0 the scalar and tensor components vanish. The spectral functions also encode the equal-time commutation and anticommutation relations of the quantum theory, implying that while all other fermion components vanish at equal time.
In addition to the spectral functions we may consider the so-called statistical functions. These are the anticommutator and commutator expectation values where field space, Dirac and flavor indices are suppressed. The statistical functions carry information about the particle density of the system, i.e. the occupation of the available modes in the system. Together, the spectral and statistical functions fully describe the time-ordered connected two-point correlation function, commonly denoted as G(x, y) = T ϕ(x)ϕ(y) − ϕ(x) ϕ(y) for the bosonic and ∆(x, y) = T ψ(x)ψ(y) for the fermionic sector. Note that in nonequilibrium settings, the timeordering occurs along the closed time path also known as Schwinger-Keldysh contour.
A. 2PI effective action real-time formalism at NLO One can derive closed and nonsecular evolution equations for the one-and two-point functions of the quarkmeson model out of equilibrium. These equations follow from the 2PI effective action Γ[φ, G, ∆], the quantum counterpart of the classical action S[ψ, ψ, σ, π], via a variational principle (see e.g. [34]). The relevant expressions for the one-and two-point functions have the form (explicit expression can be found in Appendix A): with shorthand notation t2 t1 dz ≡ t2 t1 dz 0 d 3 z and the dependence of the self-energies Σ i on φ, G, ∆ is implied. Similar expressions hold for the statistical functions. On the left hand side the Klein-Gordon or Dirac operators act on the corresponding expectation value. Thereby, effective masses take into account local quantum corrections. On the right, the effects of quantum fluctuations appear in so-called memory integrals that encode the generally non-Markovian effects of fluctuations in the past. The source term J φ in the field equation arises in the chirally broken case and describes the fermion backreaction on the field. It pushes the field to nonzero field expectation values even in the case where φ(t) = ∂ t φ(t) = 0 at initial time.
In order to carry out explicit computations, the selfenergies Σ i need to be approximated. Here we deploy an expansion to next-to-leading order (NLO) in 1/N for the bosons, where N is the number of scalar field components, and a NLO expansion in g, the Yukawa coupling. The large N expansion provides a controlled nonperturbative approximation scheme, which at NLO includes scattering as well as off-shell and memory effects, capable of handling relatively large couplings [35]. The loop-expansion in g to NLO contributes with a fermionboson loop originally discussed in [1]. The explicit equa-tions of motions are presented in Appendix A, where also the self-energy expressions for the given approximation scheme are provided. To study the time evolution of the system, we iteratively solve the equations of motion without further approximations.

B. Initial conditions
The derivation of the nonequilibrium 2PI effective action and the equations of motions following from it rely on the assumption of a Gaussian initial state. This corresponds to a system initially exhibiting the characteristics of a noninteracting theory. However, higher order correlation functions build up during the subsequent time evolution. While this appears at first sight to correspond to a very limited choice of initial conditions, it still allows for a wide variety of different configurations through which we can determine for instance the energy density ε init at the beginning of our computation. In particular, the Gaussian initial state represents a genuine nonequilibrium state in the fully interacting nonequilibrium system, in which the time evolution takes place.
We allow for spontaneous symmetry breaking by using a negative mesonic bare mass squared m 2 < 0 in the classical potential of the system. Since the initial state is determined by a free theory with m 2 = m 2 init > 0, the sign flip of m 2 leads to a quench of the classical potential from positive to negative curvature in the first time step. At initial time, the classical potential is minimal at vanishing field expectation value while the minimum at t > 0 becomes nonzero by taking m 2 < 0.
A Gaussian initial state can be fully specified in terms of the one-point and two-point functions. Since the field evolution equation involves second order time derivatives, one has to specify both the sigma field value and its initial time derivative. We choose the latter to vanish and refer to the initial field expectation value as σ 0 , where σ(t) now denotes the expectation value of the sigma field. As pointed out above, due to the presence of a finite bare quark mass m ψ the field can move away from σ 0 = 0 due to the backreaction with the fluctuations of the theory. We specify the initial conditions for the two-point functions in terms of the spectral and statistical components. The initial conditions for the bosonic (fermionic) spectral functions are fully determined by the equal-time commutation and (anti)commutation relations (6). For the remaining statistical functions we employ free-field expressions with a given initial particle number. The bosonic statistical function then reads with i = σ, π and where at initial time t = t = 0 the dispersion is set to ω i (0, |p|) = |p| 2 + m 2 init with initial mass squared m 2 init > 0 and the particle distribution given by n i (0, |p|) = 0. For the fermions the free statistical function can be written as where we choose the initial dispersion to be ω ψ (0, |p|) = |p| 2 + m 2 ψ and the initial particle distribution to be constant, i.e. n ψ (0, |p|) = n 0 .
The energy contained in the initial state via ε init determines the temperature at which the system thermalizes. By preparing different initial conditions, we can study the thermalization process towards different temperatures and hence phases of the model as sketched in Figure 1.

C. Numerical implementation
As is customary in the context of the 2PI effective action, we discretize the system on the level of the equations of motion (8). The explicit form of the equations allows us to deploy a leap-frog scheme, where in particular the fermionic two-point functions are discretized in a temporally staggered fashion. The two-point functions, as the name suggests, carry an explicit dependence on two temporal coordinates. Since the memory integrals contain the full time history, the required memory grows quadratically with the number of time steps. In order to keep the computation manageable we reduce the memory burden by exploiting isotropy and homogeneity, which reduces the effective spatial dimensions to one. A modified Fourier transform based on Hankel functions allows us to evaluate the self-energy contributions in coordinate space and to simplify the convolutions in the memory integrals in momentum space. For this project we extended the code used in Ref. [2] to include the additional nonvanishing fermionic two-point functions present in our setup. (The source code for this project is publicly accessible via the Zenodo repository under [36]. ) We choose the parameters of the theory and the UV cutoff such that the phenomenology of low-energy effective theory is reproduced. The numerical time evolution is computed using a spatial grid with N x = 200 lattice points and a lattice spacing of a x = 0.2. The time step size is chosen to be a t = 0.05 a x guaranteeing energy conservation at the level of a few percent for the times analyzed. In the following all dimensionful quantities will be given in units of the pseudocritical temperature T pc , which has the value T pc = 1.3 a −1 x determined according to the procedure described in Section IV B, see Figure 18. Subsequently, we use m to denote the dimensionless ratio m/T pc and likewise for all other dimensionful quantities.
Interactions between the macroscopic field, the bosonic and the fermionic propagators lead to an exchange of en-T crossover region χs broken phase χ symmetric phase n 0 = 0, 0 > 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " + L v r T g X R D W + J U 7 E b g i B X g z c R L I s = " > A A A G P n i c h V T d j t t E F H Y L g S Z Q 2 M J V x c 2 U 3 Z U Q s l z b i f N T K V X V R a I X R S p s d 1 t p E 0 X j 8 d g Z M j / W z L j b Y E U 8 D b f w C r w G L 8 A d 4 p Z L z j g p 1 J t F T B T r + J z v f H P + f N K S M 2 P D 8 L c b N 9 9 5 t / P e + 7 e 6 v Q 8 + v P 3 R x w d 3 P j k 3 q t K E n h H F l X 6 Z Y k M 5 k / T M M s v p y 1 J T L F J O X 6 S r E 2 d / 8 Y p q w 5 R 8 b t c l n Q t c S J Y z g i 2 o F g d 3 j + Q i R F M U H v n o a G Z Y I T C 8 P 4 T 3 x c F h G I T N Q f t C t B M O v d 1 5 t r j T u T 3 L F K k E l Z Z w b M x F F J Z 2 X m N t G e F 0 0 5 t V h p a Y r H B B L 0 C U W F A z r 5 s c N u g Y N B n K l Y a / t K j R v u 1 R Y 2 H M W q S A F N g u z V W b U 1 5 n u 6 h s P p 7 X T J a V p Z J s L 8 o r j q x C r i A o Y 5 o S y 9 c g Y K I Z x I r I E m t M L J S t d U u h c b l k 5 H U r k 9 q 8 K t o K y 1 Y / b D V O 4 i z V W K 9 r s 8 Q l N U G B L T y 5 K h g J z k 5 9 q y k 1 f q k M c + 1 g s t h 3 z C h R u u m W C U r I T i g N U x M G 0 q l W 0 z w + q k N v g K j C Q N o j S 3 I z J O P U N M U M W 6 P z N I H 1 + a G p 2 r 2 M 7 e R 4 3 R R 7 Y p w d H / g R G j g u F k 4 W g h h h F F I 6 w / 6 R Q v s b j j / S M e u 0 F n T c p I g E r H Z q + M w n q D 7 G T A / g 2 t r d S 2 b i S A 6 A R r i c E + D u s 4 n g o F w Y D 7 z Q r s F A 6 7 P o x I r J 2 D 2 G f y m z J B S Y c U X k C q U g E R q a r C Q c K e K d Q P S Y Y e d f o C K 0 Z e X e J s y 1 y S A f B 0 4 D p t J r 3 X z z 1 h O L Z g J v l F f j q / v y 4 i C R e i 4 5 H l 9 V Z E T j K N q 3 3 A p O e l K K c 6 n 1 h B y P R h m d F s J m o 6 4 Q a 2 b 0 x E 8 V U K y + 0 N / 1 K / G 5 2 u c w m V X I z K c H r y U K M h u D J s n H Y + h / d 0 L 9 + h y r 7 S S p e 4 1 u 2 z S z x r 7 K 0 1 2 + v N t e 2 1 l c 2 n k w V f C G 4 H d 4 K 7 Q T t 4 G G w G L 4 K t Y D d g w d f g V / A 7 + F P 7 V P t S + 1 b 7 P g 6 9 M D f h 3 A q m T u 3 H X 9 D Y t l Y = < / l a t e x i t > FIG. 1. Sketch of the setup deployed in this study. We consider the real-time evolution from nonequilibrium initial states characterized by an energy density sourced either through a finite σ field expectation value (blue circle) or a nonzero occupancy of fermionic modes (orange triangle). Depending on the initial energy contained in the system, one of three discernible final states, the chiral broken phase, the crossover regime or the (almost) symmetric phase is approached.
ergy between the different sectors. To observe an efficient energy exchange and equilibration process at computationally accessible times, it is necessary to study large couplings. We choose the quartic self-coupling λ = 90.0, the Yukawa coupling g = 5.0, the bare mass squared m 2 = −0.0047 and the bare fermion mass m ψ = 0.15. These parameters not only allow us to observe the equilibration of the system on timescales accessible computationally but also lead to reasonable values for the observables when compared to the phenomenological values known at T = 0, where the pion decay constant is f π 93.5 MeV, the meson masses are m σ 400 MeV and m π 135 MeV, and the constituent quark mass is m q = 350 MeV [37]. The above choices are close to that used in equilibrium computations of the quark-meson model with functional methods and a physical ultraviolet cutoff Λ UV ≈ 1 GeV. Moreover, a large value of the mesonic self interaction λ is also present in functional QCD computations [9,38]. There it can be shown that the self-interaction is of subleading relevance for the fluctuation dynamics despite the large size of the coupling. Moreover, these studies as well as a comparison of the quark-meson model to QCD, (see e.g. [10]), reveals, that a one-to-one correspondence of the low-energy limits of both theories in quantitative approximations to the full dynamics in the quark-meson model either requires a far smaller UV cutoff for the latter or a systematic improvement of the model towards QCD-assisted low-energy effective theories [8,9]. In the present work we restrict ourselves to studying the qualitative properties of the nonequilibrium dynamics as a first step.
When identifying the sigma field expectation as pion decay constant, we can reproduce f π < m π < m q < m σ at low temperatures. At the lowest temperatures considered in this work, we find f π /m π 0.65, which is very close to the phenomenologically known value of approximately 0.69, the meson mass ratio m σ /m π 1.75, smaller than the vacuum value of around 2.9 but expected to increase when going to lower temperatures, and m q /m π = 1.45, being on the order of magnitude with zero-temperature value of 2.6. Hence we expect our findings to qualitatively reproduce the QCD dynamics. Note however, that the meson mass ratio m σ /m π < 2 leads to another order of the thresholds for scattering processes, and hence respective difference in the spectral functions.
For the bosonic sector, we use vacuum initial conditions, i.e. n φ (t = 0, |p|) = 0. The initial mass is fixed at m 2 init = 0.008. The fermion initial distribution is chosen to be constant n ψ (t = 0, |p|) = n 0 . We study simulations with fluctuation dominated initial conditions where the fermion number n 0 is varied between 0 and 1 while the initial field value is σ 0 = 0. Furthermore, the field dominated initial conditions with a nonvanishing field value of σ(t = 0) = σ 0 between 0 and 2.0 with vanishing fermion number n 0 = 0 are investigated. Unless otherwise specified, plots are shown for the case n 0 = 0 and σ 0 = 0. For plots showing spectral and statistical functions in frequency space a cubic spline interpolation of the data points is employed.

III. SPECTRAL FUNCTIONS
In this section we explore the nonequilibrium evolution of the quark-meson model from the point of view of its quark and meson spectral functions. As these quantities are derived from the two-point correlation functions, they provide insight on the (quasi)particle content of the theory, the dispersion relation of propagating modes and their decay widths, providing insight into the modification of the system due to the presence of a (non)equilibrium medium. Our numerical simulations find clear indications for quasiparticles in both the IR and UV, revealing the presence of additional light propagating fermion modes for temperatures above the pseudocritical temperature.
It is convenient to analyze the spectral functions in the Wigner representation where the Fourier transformed spectral function can be interpreted as the density of states such that its structure provides information about the quasiparticle states of the system. Therefore, the temporal dependence of the unequal-time two-point correlation functions on the two times t and t is rephrased in terms of Wigner coordinates: the central time τ = (t + t )/2 and the relative time ∆t = t − t . The dynamics in ∆t describes microscopic properties of the system while the evolution in τ describes macroscopic properties governed by nonequilibrium characteristics of the system. In order to study the frequency spectrum of the spectral functions, we then apply a Wigner transformation to the propagators. This corresponds to a finite range Fourier transformation of the propagators with respect to the relative time ∆t, which is constrained by ±2τ in initial value problems where t, t ≥ 0. As a result, we obtain the frequency space spectral function with analogous expressions for all statistical functions. For a real and antisymmetric spectral function (as in the bosonic case and for the fermionic scalar, vector and tensor components) as well as for an imaginary and symmetric spectral function (as for the fermionic vector-zero component) the Wigner transform ρ(τ, ω, |p|) is imaginary. Due to symmetry, it is sufficient to present the Wigner transformed spectral functions for positive frequencies ω. Since the frequency space spectral functions are imaginary in our definition, we always plot −iρ in the subsequent sections, thereby omitting the −i in the plot labels to ease notation. The commutation and anticommutation relations (6) can be rephrased in frequency space, where they are referred to as sum rules. In our numerical computations, the bosonic and fermionic sum rules are satisfied at the level of O(10 −2 ) and O(10 −6 ), respectively.

A. Establishing thermal equilibrium at late times
Before embarking on a detailed study of the dynamical approach to thermal equilibrium, we first ascertain that our simulations of the quark-meson model reach the unique equilibrium fixed point at late times. We do so by observing the dynamic emergence of the fluctuationdissipation theorem. One needs to keep in mind that as discussed in [1], the idealized thermal equilibrium state cannot be reached in principle due to the timereversibility of the evolution equations. The simulation approaches the state more and more closely over time and at some point becomes indistinguishable from it for a given resolution. Hence we expect the computation to approach a steady state.
The fluctuation-dissipation theorem is reflected in a particular property of the spectral and statistical functions in thermal equilibrium: they are not independent of each other. In four-dimensional Fourier space it reads with n BE (ω) = (e βω − 1) −1 being the Bose-Einstein and n FD (ω) = (e βω + 1) −1 the Fermi-Dirac distribution. In (14) the frequency ω is the Fourier conjugate to the relative time ∆t = t − t as the time-dependence of F eq and ρ eq can be fully described in terms of ∆t due to the time-translation invariance of thermal equilibrium. Out of equilibrium the independence of F and ρ manifest itself in the fact that the ratio F/ρ in general carries a momentum dependence. The equilibrium relation (14) on the other hand allows us to define the generalized particle distribution function [34] with a negative (positive) sign for bosonic (fermionic) components. This kind of distribution function has been studied in the context of nonthermal fixed points in relativistic as well as nonrelativistic scalar field theories [39]. Considering (15) the approach of thermal equilibrium in a general nonequilibrium time evolution setup is characterized by n i (τ, ω, |p|) → n BE/FD (ω). In Figure 2 we show the time evolution of the particle distribution defined in (15) for low and high momenta (left and right column). One can see that at late times (red curves) the same shape is approached for small and large momenta, whereas at early times the distribution functions differ from each other. This loss of momentum dependence is required for the thermalization process and reflects the emergence of the fluctuation-dissipation relation in the equilibrium state. From the late-time distributions shown in Figure 2 one can already guess that thermal distribution functions are reached.
We also observe that the evolution of the effective particle number is different for fermions and bosons. The bosonic distribution functions n σ and n π show strong oscillations along frequencies at low momenta whereas oscillations at high momenta are weak. Since the particle distributions are computed by taking the ratio of the statistical and spectral functions, n i plotted against ω essentially describes how similar the peaks shapes of F and ρ are. In the high-momentum range we find that the quasiparticle peaks of the bosonic statistical and spectral functions resemble one another from early times on, while in the low-momentum range more time is required for the peak shapes become aligned. In contrast, the  (15) for bosonic and fermionic components (rows) and two different momenta (left and right column). At late times (red curve) the effective partcle number becomes timeand momentum-independent and approaches the shape of a Bose-Einstein and Fermi-Dirac distribution, respectively. The shown data is interpolated using a cubic spline. Dimensionful quantities are given in units of the pseudocritical temperature Tpc (cf. Figure 18). quarks show an opposite behavior. Their distributions have much stronger frequency oscillations for large momenta than for small momenta, i.e. it takes longer for the high-momentum modes to approach a thermal distribution.
Putting the pieces together, we can see that a redistribution of the occupancies in fermionic and bosonic degrees of freedom occurs during the nonequilibrium time evolution. While the timescales to converge to thermal distribution functions depend on the particle species and the momentum modes we find that the distribution functions all become stationary for times τ 100, reflecting the time-translation invariant property of thermal equilibrium.
Although Figure 2 already indicates the approach of thermal distribution functions, we still need to prove whether our final state actually fulfils the fluctuationdissipation theorem. For a quantitative analysis, we compute the generalized Boltzmann exponents with positive (negative) sign for bosonic (fermionic) components. In thermal equilibrium, the fluctuationdissipation theorem (14) requires these exponents to suffice A i (τ, ω, |p|) = βω, implying in particular that they become independent of momentum |p| and time τ , where the latter is fulfilled by our late-time states. A linear fit of our simulation data for the generalized Boltzmann exponents to βω yields the thermalization temperature T p = β −1 p , which can in general be τ -  (16) shown as a function of frequency ω at a given momentum |p| for bosonic and fermionic components. For better visibility, only every 39th data point is shown. Using a linear fit one can determine the slope β and hence the temperature T for each component. The temperature T indicated in the plot is averaged over all momenta and the three components. Right: The relative deviation from the thermalization temperature ∆ = (Ti −T )/T shown for all three components as a function of momentum. The results for the bosonic and fermionic sectors agree very well.
In both plots dimensionful quantities are given in units of the pseudocritical temperature Tpc (cf. Figure 18).
dependent. An example for such a fit is presented on the left side of Figure 3. The plot shows that the Boltzmann exponent of all three components i = σ, π, V nicely fits to the same line with slope β. We compute the temperature averaged over all momenta to obtain T i for each component. The system temperature denoted by T is taken to be the mean over all three components.
For every simulation, we compute the temperatures at each momentum |p| and study the momentum dependence of the obtained temperature T p . As pointed out in [28], thermodynamic relations can become valid before real thermal equilibrium is attained, a phenomenon known as prethermalization. Thermal equilibrium is characterized by T p being equal to some equilibrium temperature for all modes |p|. On the right side of Figure 3, the deviations from the mean thermalization temperature T are plotted. As can be seen, the deviations are very small. Hence, the Boltzmann exponents at late times τ become momentum-independent and the late-time states are thermal in the sense that they fulfill the fluctuation-dissipation theorem. The thermalization temperatures for all simulations in this work have been determined at time τ = 130. For the example shown in Figure 3 it was checked that the thermalization temperatures found in the time range between τ = 100 and τ = 160 are constant at the level of O(10 −3 ). We have checked for all simulations in this work that the temperature has reached a stationary value at time τ = 130.
Having clarified the successful approach to quantum thermal equilibrium in our system, we are now able to study the differences during the out-of-equilibrium evolution leading to the thermal states in detail.

B. Nonequilibrium time evolution of the spectral and statistical functions
In this section we study the dynamics of the thermalization process, starting from fluctuation or field dominated initial conditions. We investigate the time evolution of the spectral and statistical functions and consider derived quantities such as particle masses and widths. While the initial conditions strongly influence the nonequilibrium dynamics taking place, the final states are universal and characterized by the initial energy density ε init that translates into a unique temperature.
The time evolution leads to the emergence of quasiparticle peaks in the spectral functions of both quark and mesons. The value of the particle mass and its decay width are a consequence of the interactions taking place among the microscopic degrees of freedom. While the initial states correspond to free particles, which would have a spectrum given by a δ-distribution located at the mass parameters of the classical action, the scattering effects included in the nonequilibrium evolution lead to peaks with finite widths in the spectrum.
In Figure 4 we present a representative set of fermionic spectral functions from the vector-zero channel, which describes the quark excitation spectrum [33,40]. The three columns correspond to three different field dominated initial conditions of increasing initial energy density, as sketched by the blue dots in Figure 1. The top row shows the Wigner space spectral function at the lowest available momentum (IR), the bottom row at the highest momentum (UV). We can identify several characteristic properties of these spectral functions from a simple inspection by eye.
In the UV a single quasiparticle structure is present  Figure 18). at all times and at all energy densities. With increasing energy density in the initial state, corresponding to an increasing final temperature, the position of the peak as well as its width increase. This is consistent with the expectation that a fermion in an energetic medium will be imbued with an in-medium mass (to lowest order in perturbation theory it would be proportional to the temperature). Higher energy densities go hand in hand with an increased chance of scattering between the fermion and the other medium constituents, which also leads to a larger in-medium width. In the UV no qualitative difference exists between the broken, crossover or symmetric phase behavior.
On the other hand, in the IR a clear distinction between the crossover region and all other energy density regimes is visible. While we also find a single quasiparticle structure at low and high initial energy densities, in the crossover region at early times no well-defined peaks are present at all. Instead as times passes two structures emerge. One dominant peak is located where one would expect the usual quasiparticle excitation to reside, another peak sits close to the frequency origin, denoting a significantly lighter additional propagating mode.
In general we find that also for the other fermionic and bosonic spectral functions the approach of the equilibrium state depends on the initial conditions. In the presence of a nonzero initial field value σ 0 , the spectral functions evolve differently than in the case where σ 0 = 0 but the fermion occupation is finite, i.e. when the initial state contains more energy in terms of fermion occupations. As pointed out in Figure 4, the most interesting dynamical features can be seen in the low-momentum area, which we therefore focus on during the following analysis.

High energy densities
Here, we study the quark-meson model at high enough initial energy densities such that the late-time evolution thermalizes in the high-temperature phase, where chiral symmetry is restored. For our analysis we compare two simulations starting from different initial conditions characterized by almost indistinguishable energy densities. One is dominated by the field σ 0 = 1.36 and n 0 = 0, while the other is dominated by fermion fluctuations σ 0 = 0 and n 0 = 0.8. The final states feature similar thermalization temperatures of T = 3.15 and T = 3.18, respec-tively. However, since the initial states are very different from each other, the evolution towards thermal equilibrium takes significantly different paths.
For such high initial energy densities, the differences in the time evolution are most apparent in the bosonic sector. This can be studied by looking at the bosonic spectral and statistical functions. Numerical results are shown in Figure 5, where only the pion spectral and statistical functions are presented since the behavior of the sigma meson is analogous. The final states of both simulations (red curve) are characterized by the same peak shapes for both spectral and statistical functions. However, the functions at intermediate times exhibit a completely different behavior.
For field dominated initial conditions (left column in Figure 5) the peak position of the spectral function moves towards smaller frequencies with time, which means that the mass of the quasiparticle state decreases during the time evolution. In addition, the nonzero initial field leads to large amplitudes in the pion statistical function at early times (lower left plot in Figure 5) which corresponds to relatively high occupancies in the bosonic sector compared to the final thermal distribution. These occupancies have to redistribute to other bosonic momentum modes |p| and the fermionic sector to let the system equilibrate.
This behavior can be readily understood from the microscopic evolution equations of the system. The finitevalued initial field drives the fluctuations in the bosonic sector because it contributes to the bosonic self-energy at initial time t = 0. Since the nonequilibrium time evolution takes into account the full time history since t = 0, these initial fluctuations not only play a role at initial time but also at intermediate times. Only at late times, the system looses memory about the details of the initial state. Since the macroscopic field only couples to the bosons directly but not to the fermions, the energy provided by the initial field is first turned into bosonic fluctuations before being transferred to fermionic modes. As a consequence, the thermalization of an initial state with nonzero initial field value shows rich dynamics in the bosonic spectra.
In contrast, for fluctuation dominated initial conditions (right column in Figure 5) one observes a continuous increase of the amplitudes of both spectral and statistical functions until the maximum is reached in the thermal state. If the initial energy density is provided via fermionic fluctuations, the thermal final state is found to be realized already at intermediate times.
The spectral functions can be used to deduce the dispersion relation and lifetimes of the corresponding quasiparticle species. Following [41] we assume for the moment that the spectral function decays exponentially and can be approximated as ρ(t, t , |p|) = e −γp|t−t | ω −1 p sin[ω p (t − t )] with a dispersion ω p and a damping rate γ p , which are both allowed to be τdependent. The corresponding Wigner transform is given by ρ(τ, ω, |p|) = ρ BW (τ, ω, |p|) + δρ(τ, ω, |p|) where ρ BW denotes the relativistic Breit-Wigner function ρ BW (τ, ω, |p|) = 2ωΓ(τ, |p|) which describes a peak with width Γ(τ, |p|) = 2γ p (τ ) at position ω = ω(τ, |p|). The term δρ ∼ exp(−2τ γ p ) describes boundary effects due to the finite integration range in (12). Since δρ decreases exponentially with τ γ p , this term is negligible for sufficiently large damping ratios and/or sufficiently late times [41]. Otherwise, the frequency space spectral function suffers under severe noise coming from boundary effects. For all times shown in this work, we find that boundary effects are irrelevant. We observe that peak shapes of the bosonic spectral functions can be well approximated by the Breit-Wigner function (17). At some given time τ , performing Breit-Wigner fits of the spectral function at all momenta |p| yields the dispersion relation ω(τ, |p|) and the momentum-dependent width Γ(τ, |p|). For initial states with high energy densities, such as considered in this section, the spectral and statistical functions exhibit quasiparticle peak structures already at early times (see Figure 5). Consequently, it is possible to fit a Breit-Wigner function to the spectral functions at any stage such that the time evolution of the dispersion relation ω i (τ, |p|) and momentum-dependent width Γ i (τ, |p|) for i = σ, π can be mapped out.
In the left plot of Figure 6 we show the dispersion relation of the pion at different times τ encoded in the color scheme. A fit of ω(τ, |p|) to the relativistic dispersion relation Z |p| 2 + m 2 at various times τ yields the quasiparticle masses m(τ ), which are shown in the inset. In the following, the stationary late-time value is denoted as m. We note that the mass corresponds to the dispersion relation in the limit of vanishing momentum, i.e. m = ω(τ, |p| → 0). The right plot of Figure 6 displays the momentum-dependent width of the pion extracted from the Breit-Wigner fits. We find a plateau in the IR and a maximum in the UV. In analogy to the dispersion, where the quasiparticle mass describes the zero-momentum limit, we can extract the asymptotic value of the width in the limit of vanishing momentum, Γ = Γ(τ, |p| → 0). Since Γ corresponds to the width of the spectral function that is peaked at the quasiparticle mass, we can be viewed as the width of the quasiparticle. As the right plot in Figure 6 indicates, Γ is increasing with time.
We can now work out the differences observed in Figure 5 in a quantitative fashion. There is an apparent difference in the approach of the late-time values of the mass m π and the width Γ π when comparing the time evolution starting from the two different initial conditions. The results are shown in Figure 7, where again only the pion data is shown because the sigma meson behaves accordingly.
For field dominated initial conditions the effective mass of the pion meson decreases during the time evolution,  Figure 18).
whereas for fluctuation dominated initial conditions it grows, albeit only slightly. This is in accordance to the previous observation of the shifting peak position for field dominated initial conditions. It is important to note that the mass of the quasiparticles is not contained in the initial state, since m init is much smaller than the particle masses of the thermal state, but generated dynamically during the time evolution. The quasiparticle masses build up from the fluctuations contained in the self-energies. Since the nonzero initial field value leads to large bosonic self-energy contributions in the beginning of the time evolution, at early times the masses are larger than in the case of vanishing initial field. The time-dependence of the spectral width shown in the right plot of Figure 7 can be understood in terms of the sum rule (13) according to which the bosonic spectral functions are normalized. Due to the additional factor of ω in the integrand, which arises from the time-derivative on one of the fields in the boson commutation relation, a larger mass automatically implies smaller widths. Consequently, the behavior of mass and width in the time evolution must be converse to each other.
After discussing the dynamics of the meson spectral and statistical functions at high initial energy densities, we now turn to the quark sector. After decomposing the Dirac structure of fermionic two-point functions and imposing symmetries, we are dealing with four components for the quark spectral and statistical functions, the scalar, vector-zero, vector and tensor components as introduced in (5). Of these four components, the vector-zero component contains information about which states can be occupied [33,40]. Since it is normalized to unity according to the sum rule (13), the vector-zero component quark spectral function can be interpreted as the density of states for the quarks.
We note that in a chiral symmetric theory with vanishing fermion bare mass one finds ρ S = ρ P = ρ µν T = 0 since only components in (4) that anticommute with γ 5 are allowed. Here, we consider a setup where chiral symmetry restoration takes place. For initial conditions with high energy densities and the corresponding final states in the high-temperature chiral symmetric regime, the quark dynamics can be studied in terms of the vector-zero and vector component.
As was shown in Figure 4, for high energy densities there is not much dynamics taking place in the excitation spectrum of the quarks. More insight can be gained by looking at the vector component which is presented   Figure 18).
in Figure 8 for the same field or fluctuation dominated initial conditions as discussed before for the bosons. The interesting case is again the evolution starting from field dominated initial conditions. The corresponding vector spectral function (upper left plot) shows that the peak position moves towards smaller frequencies, just as in the bosonic case. It indicates that the energy of both meson and quark quasiparticles decreases during the time evolution. However, it is important to note that -in contrast to the mesons -the amplitude of the fermion statistical function is increasing during the time evolution. As discussed before, the nonzero initial field leads to strong fluctuations and hence occupancies in the bosonic sector. It takes time for these fluctuations to be transferred to the fermionic sector, which is why we observe that the fermion occupation grows slowly during the time evolution.
For the fluctuation dominated initial conditions we again observe that the spectral and statistical function approach their late-time behavior very quickly. We conclude that the available states and their occupation quickly approach their thermal final state if energy is provided in terms of particles rather than the field in the initial state.

Intermediate energy densities
From Figure 4 we can see that the most interesting dynamics is taking place for systems thermalizing in the  Figure 18). crossover region. Thus, we aim to study the evolution of the vector-zero quark spectral function for two simulations thermalizing in the cross-over region.
Again we compare two simulations employing field or fluctuation dominated initial conditions, respectively, but in this case we are able to probe initial conditions that lead to the same late-time state. When comparing the late-time field expectation value σ, the mass ratio m σ /m π , and the temperature T of the final state of these two simulations, we find that the respective quantities differ by less than 0.5 %. Also, the shape of the spectral and statistical functions in frequency space are the same for both bosonic as well as fermionic components. Quantitatively, we find that |ρ 1 − ρ 2 |/ max (ρ 1 ) is smaller than O(10 −2 ) for all frequencies ω and momenta |p|, where the indices 1 and 2 denote the two simulations compared and max (ρ) the maximal amplitude of ρ. Larger deviations are observed for the vector-zero component statistical function and for the tensor component spectral and statistical functions, where the amplitudes are of order O(10 −7 ) such that numerical inaccuracies come into play. In conclusion, we consider the late-time state of the two simulations to be the same thermal state, universal in the sense that the dependence on the initial conditions is lost. It is characterized solely by a temperature of T = 1.04, a mass ratio of m σ /m π = 1.46 and a field expectation value of σ = 0.33. As we will see later, this corresponds to a state in the crossover region.
The regime of intermediate energy densities distinguishes itself from high and low energy density initial conditions by showing a double-peak structure in the quark spectral functions. Our findings in a nonperturbative real-time setting corroborate previous observations of such double peak structures with perturbative computations or spectral reconstructions reported e.g. in [33,[42][43][44][45][46][47][48].
First, consider the vector-zero component describing the excitation spectrum of the quarks. In Figure 9 we show the time evolution of both spectral and statistical functions. As before, for fluctuation dominated initial conditions (right column) the system quickly approaches the shape of the late-time two-point functions. However, in the case of field dominated initial conditions, the double-peak structure of the spectral function only emerges at later times. At early times, the spectral function reveals a single broad structure.
We further point out that the statistical function F 0 decays to zero during the time evolution, implying that the fermion occupation is not contained in the vectorzero component but in other components. This agrees well with the effective quasiparticle number that has been  Figure 18). employed previously [3,40], with effective mass M ψ (t) = m ψ + hσ(t). This definition of an effective particle number only provides a good description of the quark content in the system if the occupations in the vector-zero and tensor component are negligible. In our computations we find that F 0 and F T are of the order O(10 −7 ) and hence irrelevant for the quark particle number. In order to study the particle content, we take into account the vector component which is shown in Figure 10. We can see that the double-peak structure observed in the vector-zero component is also visible in the vector component, in particular in both spectral and statistical functions. From this we learn that the additional light degrees of freedom, provided in the low-frequency peak of the quark spectral density, is actually occupied in terms of the vector component quark statistical function. Hence, for states thermalizing in the crossover temperature regime, there is an additional light mode with finite occupation in the quark sector available to participate in the dynamics.
We further observe that for fixed momentum |p| the energy of the light mode increases with rising temperature. At sufficiently high temperatures this additional mode reaches energies comparable with the main quasiparticle mode such that the two peaks merge into the single peak persistent in the high-temperature regime.
We conclude this section with a comment on the dynamics found for initial states with low energy-densities. In contrast to the cases of intermediate and high energy densities, we find well-defined quasiparticle peaks for both quarks and mesons. The smaller energy density leads to lower thermalization temperatures and a stronger chiral symmetry breaking, reflected by a mass difference between the σ and π mesons. After discussing the nonequilibrium time evolution of the spectral functions, we now turn to the equilibrium properties.  temperatures reflect the crossover transition of the quark-meson model from the chiral broken to a chiral symmetric phase. We find that the shapes of the final states are universal in the sense that they only depend on the temperature and not on the details of the initial state.

Mesons
Information about the different phases of the model can be obtained from the temperature-dependence of the late-time thermal spectral functions of the mesons. We find that the shape of the bosonic spectal functions is described by a Breit-Wigner function for all considered temperatures. Thereby, the width and the position of the Breit-Wigner peak only depend on the temperature but not on the initial conditions chosen.
As discussed in Section III B 1, the momentumdependent width and the dispersion relation are obtained by applying Breit-Wigner fits to the spectral functions. Although the Breit-Wigner function (17) has two parameters, the width Γ(τ, |p|) and the peak position given by ω(τ, |p|), there is only one free parameter since the normalization condition given by the sum rule (13) must be satisfied.
In the right plot of Figure 6 we already saw that there is a characteristic momentum mode |p| at which the momentum-dependent width becomes maximal. This corresponds to the momentum at which the decay is strongest and can be considered as the main decay mode, in the following denoted by Q. In the left plot of Figure 11 we show the main decay mode Q as a function of temperature for both meson species. At low temperatures, the strongest decays are found in the IR, whereas at high temperatures the strongest decays occur in the UV. There is an abrupt change at some critical temperature, above which Q > 0 meaning that the momentumdependent width has a maximum at a nonzero momentum, as shown by the upper line in the inset. Comparing the momentum-dependent width at low T and high T , we can see that the transition from the chiral broken to the chiral symmetric phase is characterized by new decay modes in the UV. Thereby, the main decay mode is suddenly shifted from the IR to the UV.
Another prominent signature for the crossover transition is provided by the quasiparticle masses of the σ and π mesons. The two meson species are distinguishable in the chiral broken phase, where they have different masses, while they become identical in the chiral sym- Right: Temperature-dependence of quasiparticle masses. Restoration of chiral symmetry is reflected in identical masses of the σ and π mesons at high temperatures. The quark q quasiparticle mass is obtained from the dominant peak of the vector-zero component quark spectral function. We also plot the "plasmino" branch p obtained from the quark spectral function. In both plots gray lines show cubic spline interpolations of the data points. Dimensionful quantities are given in units of the pseudocritical temperature Tpc (cf. Figure 18). metric phase. When plotting the meson masses as a function of temperature, as shown in the right plot of Figure 11, we can nicely visualize the restoration of chiral symmetry, manifest in the quasiparticle masses of σ and π becoming identical (pink and cyan data points). We observe a softening of the masses at intermediate temperatures, i.e. the quasiparticle masses are minimal in the temperature region where the crossover phase transition occurs. Decreasing masses indicate growing correlation lengths. In the limit of a second order phase transition, which is characterized by diverging correlation lengths, the masses would vanish at the transition point. In the high-temperature range, masses grow with rising temperatures. This reflects that the quasiparticle masses can be considered as thermal masses in the sense that they contain self-energy contributions and are generated by quantum fluctuations, which increase with temperature.
We further note that one could also study the temperature dependence of the width Γ = Γ(τ, |p| → 0) instead of m = ω(τ, |p| → 0). However, the information is equivalent due to the normalization of the spectral functions, as pointed out above. Consequently, the behavior of Γ is converse to the behavior of m and not presented here explicitly. The with Γ is small at low temperatures, strongly grows towards intermediate temperatures where it reaches a maximum value in the crossover temperature regime, and then decays slowly when going to higher temperatures.

Quarks
Let us now consider the thermal spectral functions for the quark sector. Several aspects of the different components invite discussion. Let us begin with a recap of the findings shown in the vector-zero component of the quark spectral function. As presented in Figure 4 the spectral density has different shapes at low, intermediate and high temperatures. In particular, the intermediate temperature range of the crossover transition is characterized by a double-peak structure. The temperature-dependence of the fermionic quasiparticle masses is depicted in Figure 11. The mass of the low-frequency mode (plasmino branch, denoted by p) grows continuously with rising T until it merges with the main peak (denoted by q), forming the wide quasi-particle peak found for initial states with large energy densities. For related studies with perturbative computations or spectral reconstructions see e.g. [33,[42][43][44][45][46][47][48]. Note also that this double-peak structure is only visible in the small momentum regime. It can be studied by considering the dispersion relation obtained from the vector-zero quark spectral function.
For temperatures below some critical temperature in the crossover regime, the vector-zero spectral function reveals the shape of a nonrelativistic Breit-Wigner function, also known as the Lorentz function, where A is a normalization constant, Γ(τ, |p|) the width q is shown by the black dashed line. For higher temperatures the behavior at small and large momenta differs as the additional low-frequency peak and the main peak merge into one peak. We perform separate fits at low and high momenta, shown by the dashed and dotted black lines. Dimensionful quantities are given in units of the pseudocritical temperature Tpc (cf. Figure 18).  13. The vector-zero quark spectral function as a function of frequency ω shown for a range of spatial momenta |p|. The three plots correspond to the same three temperatures as in Figure 12. The purple line indicates peak position of the spectral function in the |p|-ω-plane and is therefore equivalent to the dispersion relation shown in Figure 12. The spectral function reveals a narrow quasiparticle peak at low temperatures. As the temperature is increased the light mode interferes with the low-momentum spectral function, leading to a broad peak at small momenta. At high momenta, the quasiparticle-peak remains narrow. Dimensionful quantities are given in units of the pseudocritical temperature Tpc (cf. Figure 18).
and ω(τ, |p|) the dispersion. When temperature is increased, the vector-zero quark spectral function ceases to be described in terms of (19) as the low-frequency mode arises and grows in amplitude. Due to appearance of the additional peak, it is not possible to perform a Lorentz fit at all temperatures. As a consequence, we choose to compute the dispersion relation of the quarks by determining the peak position of the main peak of ρ 0 . The obtained dispersion relation is shown for three temperatures in Figure 12.
For low temperatures, where no additional peak is present, the quark dispersion is well-described by a relativistic dispersion relation, see left plot of Figure 12. When going to intermediate temperatures, the additional light mode leads to a double-peak structure. As long as the two peaks are distinguishable, one can determine the dispersion relation of the main peak, which yields the same shape as in the low-temperature regime. However, when the main peak and the side peak merge into a single peak, the dispersion relation obtained from the overlap of the two peaks has a dispersion relation of the form shown by the middle plot of Figure 12. There is a clearly visible dip in the dispersion, showing that for small momenta the peak position is determined by the light mode, while for large momenta the peak position is determined by the main peak. We can fit the low-momentum and high-momentum areas separately to a relativistic dispersion relation, as shown by the dashed and dotted lines in Figure 12. When considering higher temperatures, the position of the dip moves towards larger frequencies and is not visible by eye anymore. However, we find that the dispersion relation cannot be described by the relativistic dispersion relation Z |p| 2 + m 2 over the whole momentum range but still distinguishes between highmomentum and low-momentum regimes. We conclude that the single peak of ρ 0 at large temperatures is still the result of an overlap of a small low-frequency peak with the main peak. More insight is gained by considering the momentum-dependence of the corresponding spectral functions, which is shown in Figure 13. The spectral function ρ 0 (τ, ω, |p|) is shown at some late time τ where the system has approached thermal equilibrium. The peak position of the spectral functions corresponds to the dispersion relations shown in Figure 12. At low temperatures, we find a single narrow quasiparticle peak. For higher temperatures, however, an additional light mode interferes with the main peak. At intermediate temperatures, where a softening of the mass occurs, the light mode and the main peak have comparable frequencies in the infrared. The superposition of the main peak and the light mode leads to a broad peak at small momenta, whereas the peak remains narrow at high momenta. As temperature increases, the light mode is only visible at higher momenta. An example is shown by the right plot in Figure 13, where one can see a small enhancement of the spectral function at low frequencies for intermediate momenta. This observation indicates that the quark spectral function harbors additional degrees of freedom at high temperatures, as compared to the lowtemperature regime.
From the dispersion relation of the vector-zero quark spectral function, we determine the constituent quark mass by taking the asymptotic value at vanishing momentum, i.e. m q = ω 0 (τ, |p| → 0). The constituent quark mass behaves analogously to the bosonic masses, i.e. a softening of the mass in the crossover temperature range occurs, see violet data points in the right plot of Figure 11. At low temperatures, the constituent quark mass lies between the σ and π masses, which is in qualitative agreement with the particle masses known at T = 0. For temperature below T 2 we find that the pion is the lightest particle in the theory. This supports chiral perturbation theory as an effective theory for QCD where only pion degrees of freedom are considered. On the contrary, at high temperatures the constituent quark mass is smaller than the meson mass. As light modes are easier to excite, they dominate the dynamics in a system. Hence, our observation matches our idea that the chiral symmetric phase is dominated by quark degrees of freedom whereas the chiral broken phase is described by hadronic degrees of freedom, in particular by pions.
Finally, we shortly discuss the scalar component of the quark spectral function. In a chiral symmetric theory with vanishing fermion bare mass, the scalar component of the quark spectral function vanishes, i.e. ρ S = 0. Although chiral symmetry is broken explicitly here, we expect the system to restore chiral symmetry at high temperatures, implying that the limit of a vanishing scalar component quark spectral function is approached. In Figure 14 we present the numerical results for a range of temperatures. The clear quasiparticle peak existing at low temperatures widens and flattens with rising temperature. The amplitude of the scalar component finally decays to zero, visualizing the predicted restoration of chiral symmetry in the course of the crossover transition. We further note that the peak position of the scalar component spectral function qualitatively shows the same behavior as the vector-zero component. The peak moves towards small frequencies at intermediate temperatures, corresponding to the softening of a mass, and is shifted towards higher frequencies at low and high temperatures.

IV. THE MACROSCOPIC FIELD
In this section, we study the time evolution of the expectation value of the macroscopic field σ(t) , for the set of different initial conditions deployed also in the previous section. In addition, we study the model for different fermion bare masses in order to analyze the effects of spontaneous symmetry breaking in the model.

A. Nonequilibrium time evolution of the field
The classical potential of the sigma field is given by where the parameter choice of m 2 < 0 allows for spontaneous symmetry breaking. Thus, the potential has the shape of a double-well with minima located at σ = ± −3!N m 2 /λ. For the parameters employed in this work the minimum is located at σ ≈ 0.04. The time evolution of a classical field in this potential is described by the classical equation of motion where spatial homogeneity and isotropy is assumed. If the initial field value, or the initial field derivative, is nonzero, the field rolls down a potential hill and oscillates until it equilibrates at the minimum of the potential. Here, we go beyond the classical theory and compute the nonequilibrium time evolution including additional  Figure 18). quantum fluctuations. As discussed above, we employ an approximation that includes quantum corrections at NLO in 1/N and g. The quantum corrections lead to an effective potential and additional terms in the field equation (21). The full evolution equation at the given approximation can be found in Appendix A.
Depending on the initial conditions the time evolution of the field shows different properties. Let us first consider field dominated initial conditions, where the initial field is set to a finite value σ 0 . The time evolution for the expectation value of the field σ(t) is shown for different σ 0 in the left plot of Figure 15. One can see that the field oscillates and eventually reaches a stationary value. In contrast to the classical theory, where the field always reaches the same equilibrium value given by the position of the potential minimum, the field reaches different latetime values. The reason is that the field itself generates quantum fluctuations as it rolls down a potential hill. These dynamically emerging fluctuations again influence the effective potential in which the nonequilibrium time evolution takes place. As the initial field value effects the amount of quantum fluctuations in the system and hence the shape of effective potential, different values of σ 0 lead to different late-time values for σ(t) . Before we come to a more detailed discussion of the plots in Figure 15, we provide some intuition for the influence of fluctuations on the effective potential.
Quantum fluctuations can be represented as loop corrections of the effective action. The effective potential is obtained when evaluating this effective action at a constant field. For a nonvanishing fermion bare mass, i.e. m ψ = 0, the chiral symmetry breaking tilts the effective potential towards negative values. Thereby, larger m ψ cause stronger tilts. On the other hand, bosonic fluctuations provide positive contributions, pushing the potential towards a symmetric shape. Together, this leads to a tilted Mexican hat potential with a minimum at some finite field expectation value. The position of the minimum of the effective action can easily be much larger than the position of the minimum of the classical potential.
The influence of these quantum corrections to the effective potential can be visualized by looking at the energy density of the system, which we compute from the energy-momentum tensor. We distinguish classical, bosonic and fermionic contributions to the energy density, with the relevant expressions presented in Appendix A 3. Quantum fluctuations are taken into account for the fermionic and bosonic parts of the energy density. Hence, the energy density reflects the amount of fluctuations in the system.
In order to study the influence of the initial field, we consider the energy density computed at initial time ε init . In Figure 16 we show the contributions from the field, the bosons and the fermions separately. The blue and pink lines show, how the field and the bosonic energy densities exhibit a positive curvature. In contrast, the fermionic contribution shown in violet leads to a tilt towards negative curvature, which is a consequence of the explicit chiral symmetry breaking. Summing the three parts together, one obtains the total energy density that has a minimum at a nonzero field value, as represented by the gray curve.
It is important to note that the energy density ε init computed at time t = 0 does not include the quantum fluctuations that are generated dynamically by the field. As the field generates further fluctuations, the effective potential is pushed towards a more symmetric form with its minimum moving towards smaller field expectation values. In order to see this, we also look at the energy density computed at late-times, where the system is thermalized. The result is shown by the black line in Figure 16. It can be seen that the energy density indeed becomes steeper and the minimum moves towards smaller field values. Thus, the energy density provides a useful quantity in order to study the impact of quantum fluctuations on the effective potential, although we emphasize that the energy density and the effective potential are two different quantities.
Having this qualitative picture of the effective quan-tum potential in mind, we can understand the behavior of the three curves shown in the left plot of Figure 15. If the initial field sits close to the minimum of the effective potential, it barely oscillates and hence almost no additional fluctuations are created dynamically. Accordingly, the shape of the potential does not change with the time evolution such that the position of the minimum stays the same. This is shown by the black line. In contrast, the field can be placed on a point away from the potential minimum. As it starts moving towards the potential minimum, the field dynamically generates fluctuations. These fluctuations change the shape of the potential, thereby altering the position of the minimum. The further away the field is from the potential minimum in the beginning, the more fluctuations are generated and the stronger the potential deforms. As we increase the distance of the initial field from the potential minimum at time t = 0, the minimum of the potential at late times moves towards zero. Examples of this behavior are depicted by the green and red curve in the left plot Figure 15.
As discussed in the previous section, energy cannot only be provided in terms of a nonzero initial field value (and the fluctuations this field generates), but also in terms of occupancies. Hence, the same late-time field value can be approached for different initial conditions. In the right plot of Figure 15 the time evolution of σ(t) is shown for the two simulations discussed in Section III B 2. The blue line displays the time evolution of the field starting from field dominated initial conditions, while the orange line shows the time evolution starting from fluctuation dominated initial conditions. For both initial conditions the quantum potential has the same minimum, characterized by a late-time stationary field value of σ = 0.33.
Although we commonly say that initial conditions with the same energy density lead to the same thermal state, there is a caveat. Two initial states thermalizing at the same late-time state usually do not have the same energy density at time t = 0 because the energy density computed at initial time does not include dynamically generated fluctuations. What one means is that different initial conditions provide the same amount of fluctuations to the system. The way they are provided depends on the initial state and partly they are generated dynamically. However, for the quantum thermal equilibrium state that is approached at late times only the amount of fluctuations introduced to the system is relevant.
B. Thermal equilibrium

Field expectation value
After discussing the time evolution of the field expectation value we now turn to its late-time properties. We denote the stationary value of the field at late times by σ. As discussed in III A, at these times the fluctuation- The energy density at initial time εinit = ε(t = 0) and at late times ε th = ε as a function of the initial field value. We present the classical, bosonic and fermionic contributions to the initial energy density separately. Together, they form a bounded shape with minimum at a nonzero initial field value (gray curve). At late times, the energy density reaches the constant shape shown by ε therm in black. The minimum of the energy density at late times corresponds to the maximal field values found. The initial field is given in units of the pseudocritical temperature Tpc (cf. Figure 18). dissipation theorem is satisfied and the system state is considered to be thermal. Thus, we consider σ to be the thermal field expectation value.
The late-time field values σ are determined by the average of field values over a time range [t * , t * + ∆t] with t * being a time at which the field is sufficiently stationary. For the results shown in this work we use t * = 130 and ∆t = 130, such that the standard deviation of the mean is O(10 −4 ) to O(10 −11 ) depending on the initial conditions used.
First, let us look at the time evolution of the macroscopic field for different initial field values σ 0 . Naively, one might expect that larger field values automatically imply increasing energy densities in the initial state, hence a higher thermalization temperature and smaller field value. However, as the discussion above already pointed out, this is not the case. In the left plot of Figure 17 we show how the late-time field value σ depends on the initial field σ 0 . With increasing σ 0 , the thermal field σ first grows and then decays to zero. The maximal value for σ is expected, when the least amount of fluctuations is generated dynamically, as these fluctuations would push the minimum of the potential and thus σ towards zero. We indeed find the largest late-time field values for σ 0 ≈ σ, which in indicated by the gray dashed line in the plot.
Secondly, we consider fluctuation dominated initial conditions where the field value is set to σ 0 = 0 while the initial fermion occupation is taken to be constant, i.e. n ψ (t = 0, |p|) = n 0 , and varied between zero and one. In the right plot of Figure 17 we can see that increasing the fermion occupation number n 0 goes along with smaller thermal field values σ. Thus, for larger n 0 higher temperatures are reached, emphasizing again that larger fermion occupation numbers lead to a rise of the fluctuations that make the effective potential more symmetric.

Crossover phase transition
Our results regarding the thermal state of the system can be summarized in an analysis of the crossover transition between the chiral broken and the chiral symmetric phase of the quark-meson model. When a system becomes thermal, the thermodynamic concept of a phase diagram can be applied. The conjectured phase diagram of the quark-meson model contains important features of the QCD phase diagram. It exhibits a chiral symmetric phase with vanishing field expectation value at high temperature T , as well as a chiral broken phase with nonzero field corresponds at low T .
In order to study the phase transition and the transition temperature, we employ two different order parameters, one deduced from the one-point function and one from the two-point functions. The first one is the field  Figure 18). expectation value of the thermalized field σ. It is nonzero in the chirally broken phase and zero in the chirally symmetric phase. Often this field value is identified as the pion decay constant f π . The second one is the mass ratio m σ /m π , where m σ and m π are the masses of the σ-meson and the pion, respectively. The masses are determined from the bosonic spectral functions as discussed in Section III C 1. In the chiral limit, the mass ratio is expected to go to unity. Starting from the different initial states analyzed, we find that the system thermalizes at different temperatures. Thereby, the dependence of an order parameter on the temperature provides insight into the nature of the phase transition. In Figure 18 we show our numerical results for the temperature-dependence of the two order parameters, σ in the upper and m σ /m π in the lower plot. Every point in the diagram corresponds to a simulation with a different initial state. As indicated in the legend, we are considering initial states of various fermion occupations, described by n 0 , and initial field values, described by σ 0 . It is reassuring to see that the order parameters obtained from field or fluctuation dominated initial conditions align themselves on a single curve, which is characteristic for a smooth crossover transition. This is yet another way of seeing that the thermal states are independent of the details of the initial conditions.
As chiral symmetry is restored with rising temperature, the field value decays to zero while the mass ratio goes down to one. The field expectation value σ is often considered as a first approximation for the pion decay constant f π . As can be seen, in the limit T → 0 some value σ O(1) is approached. At the lowest temperature considered we find σ/m π 0.65, matching the phenomenological value f π /m π 0.69 [37]. Further, we can see from the lower plot in Figure 18 that the mass ratio is only m σ /m π 1.8 at the lowest temperatures available, which is smaller than the expectation from the known values of the masses. However, the mass ratio is expected to further increase with decreasing temperature.
We perform a cubic spline fit to the data points and identify the inflection point of the field σ as the pseudocritical temperature of the crossover T pc . We indicate the inflection point of both the field and the mass ratio in the plots of Figure 18. It can be seen that the temperatures deduced from the two different order parameters are comparable with each other. We find that the pseudocritical temperature is of the order of the pion mass. This is in agreement with the expectation of the QCD phase transition being at around 150 MeV for vanishing baryon density.

C. Spontaneous symmetry breaking
We have seen that the explicit chiral symmetry breaking in the system leads to nonzero field expectation values. Here, we analyze the limit of vanishing explicit symmetry breaking, i.e. m ψ → 0, with spontaneous symmetry breaking still present.
If the fermion bare mass vanishes, i.e. m ψ = 0, the action of the quark-meson model (1) is invariant under chiral SU L (2) × SU R (2) ∼ O(4) transformations and therefore symmetric under chiral symmetry. Still, a nonzero field expectation value can break this symmetry spontaneously. For nonzero fermion bare masses, chiral symmetry is explicitly broken and the minimum of the potential is located at some nonzero field value. If the field expectation value stays nonzero for m ψ → 0, we expect to observe spontaneous symmetry breaking.
We compare simulations with different fermion bare masses m ψ while all other parameters of the theory are kept fixed. The system is studied for initial conditions with σ 0 = 0.62 and vanishing fermion and boson occupations, i.e. n 0 = 0. In the left plot of Figure 19, we show the time evolution of the field expectation value σ(t) , for three examples with different fermion bare masses. As before, the field oscillates before it equilibrates to the thermal late-time value σ.
Since the fermion bare mass m ψ governs the strength of the chiral symmetry breaking and thus the deformation of the potential, increasing fermion bare masses yield larger values for σ. At the same time m ψ determines the fermion backreaction on the field, i.e. how strong the field is pushed away from its current value. The field only reaches a stationary value, if the backraction from the fermions on the field and the bosonic interactions with the field balance out.
Here, fermion bare masses with values from O(10 −4 ) ranging to O(1) are considered. We find that the field approaches the asymptotic value σ = 0.48 for m ψ → 0, which is shown in the right plot of Figure 19. This anal-ysis shows that our numerical simulations of the quarkmeson model reproduce the expected spontaneous symmetry breaking in the limit of vanishing fermion bare mass.

V. SUMMARY AND CONCLUSION
Motivated by current experimental studies of the QCD phase diagram in heavy-ion collisions, we investigated the dynamical approach of the quark-meson model to thermal equilibrium using a range of different initial conditions dominated by either the sigma field or fermionic fluctuations. The time evolution of one-and two-point functions was computed numerically using closed equations of motion derived from the 2PI effective action at NLO in 1/N and the Yukawa coupling.
We show that our simulations correctly capture the universal thermal fixed point of the dynamics, which depends only on the energy density of the initial condition. The crossover phase transition from the chiral broken phase at low temperatures to the chiral symmetric phase at high temperatures is reproduced by the late-time equi-  Figure 18).
librium states. Thermalization in the chiral broken phase is characterised by a finite field expectation value, a mass difference between the sigma meson and the pions as well as narrow quasiparticle peaks in the spectrum. The restoration of chiral symmetry in the high-temperature regime expresses itself in the field expectation value decreasing to zero, the mass ratio of σ and π mesons going to unity and the scalar component of the quark spectral functions decaying to zero. Our investigation focused in detail on the dynamical thermalization revealing differences in the time evolution depending on the initial state employed. We not only studied the time evolution of the field expectation value but also probed the dynamical properties of the two-point functions, expressed in terms of the spectral and statistical functions, which carry information about the available quasiparticle states and their occupation in the system, respectively. For initial states with vanishing initial field but energy supplied by fermion occupation, the spectral and statistical functions of both quarks and mesons approach their late-time thermal shapes already at early times. In contrast, if the energy density is predominantly provided by the nonzero initial field value, the redistribution of energy from the field first to the bosonic sector and subsequently to the fermionic sector leads to high occupancies of the mesons at intermediate stages. This is also reflected by the different behavior found in the time evolution of the quasiparticle masses depending on the initial conditions.
The deployed nonequilibrium setup of the quark-meson model captures important features of the low-energy behavior of QCD. By studying the temperature dependence of the quasiparticle masses, we find that the lightest degrees of freedom are given by the pions at temperatures below and by quarks above the phase transition. This implies that quarks are the relevant degrees of freedom at high temperatures while pions dominate below the critical temperature. Furthermore, we learn from the width that at high temperatures the more energetic highmomentum decay modes are more pronounced than for low temperatures.
The nonvanishing expectation value of the sigma field describes the order parameter of the chiral phase transition. Its dynamics depends on the initial state. If the initial field value is close to the minimum of the effective potential, the field remains almost constant. Otherwise, the field rolls down a potential hill and starts oscillating, thereby dynamically generating fluctuations.
Having shown that the dynamics of the thermalization process reveal interesting features before approaching the final thermal state, we lay the foundation for future investigations of the quark-meson model with nonzero baryon chemical potential. In particular, the possibility of probing the dynamical thermalization of systems surpassing the critical point of the chiral phase transition is of outermost interest.

Approximation Scheme
At NLO in the large N expansion, the nonlocal interaction terms in the field equation (A1) are given by dt p I ρ (t, t , |p|)F σ (t, t , |p|) + I F (t, t , |p|)ρ σ (t, t , |p|) σ(t ) , where I F and I ρ are the spectral and statistical components of the summation functions presented below in (A12).