Spectral, statistical and vertex functions in scalar quantum field theory far from equilibrium

We compute the far-from-equilibrium dynamics of relativistic scalar quantum fields in 3+1 space-time dimensions starting from over-occupied initial conditions. We determine universal scaling exponents and functions for two-point correlators and the four-vertex in a self-similar regime in time and space or momenta. The scaling form of the momentum-dependent four-vertex exhibits a dramatic fall-off towards low momenta. Comparing spectral functions (commutators) and statistical correlations (anti-commutators) of field operators allows us to detect strong violations of the fluctuation-dissipation relation in this non-perturbative infrared regime. Based on a self-consistent expansion in the number of field components to next-to-leading order, a wide range of interaction strengths is analyzed and compared to weak-coupling estimates in effective kinetic theory and classical-statistical field theory.

We compute the far-from-equilibrium dynamics of relativistic scalar quantum fields in 3+1 spacetime dimensions starting from over-occupied initial conditions. We determine universal scaling exponents and functions for two-point correlators and the four-vertex in a self-similar regime in time and space or momenta. The scaling form of the momentum-dependent four-vertex exhibits a dramatic fall-off towards low momenta. Comparing spectral functions (commutators) and statistical correlations (anti-commutators) of field operators allows us to detect strong violations of the fluctuation-dissipation relation in this non-perturbative infrared regime. Based on a self-consistent expansion in the number of field components to next-to-leading order, a wide range of interaction strengths is analyzed and compared to weak-coupling estimates in effective kinetic theory and classical-statistical field theory.

I. INTRODUCTION
Non-equilibrium dynamics of relativistic scalar quantum field theory is an important cornerstone in our understanding of the evolution of the early universe. For a significant class of inflationary scenarios, where the scalar field describes the inflaton or even the Higgs field, the evolution traverses a far-from-equilibrium regime triggered by dynamical instabilities [1,2]. In general, these processes create highly occupied excitations at some characteristic momentum scale Q 0 , which differs significantly from the temperature of a thermally equilibrated system with the same energy density and particle number. The over-occupied system subsequently relaxes in terms of self-similar cascades, transporting particles towards low momentum scales [3,4] and energy to higher momenta [5,6]. Similar phenomena are predicted to characterize the dynamics of the highly excited quark gluon plasma during the early stages of a heavy-ion collision [7][8][9][10]. Recently, self-similar scaling phenomena have been experimentally discovered in ultra-cold quantum gases far from equilibrium, showing remarkable universal properties in the non-perturbative infrared regime at sufficiently low momenta [11][12][13]. Theoretical descriptions of the underlying non-thermal infrared fixed points typically involve effective kinetic theory [4,[14][15][16][17][18] or classical-statistical approximations [4,[19][20][21][22][23][24] in the weak coupling limit.
In this work, we compute the time evolution of relativistic scalar quantum fields in 3+1 space-time dimensions starting from over-occupied initial conditions. We consider a self-interacting N -component field theory often employed in the context of scalar inflaton models, and which can also be taken to describe the Higgs sector of the standard model of particle physics for N = 4. A self-consistent expansion in powers of 1/N to next-toleading order (NLO) provides a non-perturbative account of the dynamics, such that we may analyze the highly occupied infrared for a wide range of interaction strengths. This has been previously employed to study far-from-equilibrium dynamics of this model, focusing on the role of a symmetry breaking field expectation value for its evolution [25]. Here, we analyze for the first time the selfsimilar scaling properties of two-point correlators and the four-vertex in quantum field theory at NLO without further approximations by numerically extracting the universal dynamical exponents and scaling functions. Our results are compared to previous weak-coupling estimates for equal-time two-point correlators using effective kinetic theory or classical-statistical field theory [4].
A further important focus of our work is the computation of unequal-time correlation functions far from equilibrium, namely the commutator expectation value of two field operators and the respective anti-commutator at different times. The former gives the spectral function, which provides essential information about the nature of the excitations, such as the possible existence or absence of long-lived quasi-particles far from equilibrium. In contrast, the anti-commutator expectation value captures quantum-statistical aspects, such as the occupation number of modes. While in thermal equilibrium the spectral (commutator) and statistical (anti-commutator) correlation functions are related by the fluctuation-dissipation relation, this is violated out of equilibrium in general. The question of whether some generalized fluctuationdissipation relation may be defined, such as underlying effective kinetic descriptions of the dynamics, has been recently addressed in this context with the help of classical-statistical simulations [26,27]. Here we compute both spectral and statistical correlation functions directly based on the underlying quantum field. While at sufficiently high momenta we recover the expected quasi-particle structure with a generalized fluctuationdissipation relation, we demonstrate that significant violations occur in the non-perturbative infrared regime.
The paper is organized as follows: In Sec. II we introduce the model and specify the class of far-fromequilibrium initial conditions we consider. The results for the equal-time two-point correlation functions and the four-vertex are presented in Sec. III. We extract the selfsimilar scaling properties for a wide range of couplings, and compare the full numerical results to analytical estimates that are obtained assuming universal scaling. In Sec. IV also unequal-time correlations are computed and the spectral as well as statistical functions are analyzed. We identify three characteristic momentum regimes and investigate the role of the fluctuation-dissipation relation in each of these regimes. We conclude in Sec. V. Three appendices provide details about the error estimates for the extraction of the scaling exponents, the Wigner transformation employed in the read-out of the spectral function, and the procedure for the numerical computation of the time evolution.

II. SCALAR QUANTUM FIELD THEORY FAR FROM EQUILIBRIUM
We consider a relativistic N -component scalar quantum field theory interacting via a quartic self-coupling λ.
Its O(N )-symmetric classical action for massless fields is given by where a summation over field indices a = 1, . . . , N and Lorentz indices µ = 0, 1, 2, 3 is implied, with four-vector x = (x 0 , x). Quantum corrections are taken into account using a large-N expansion at next-to-leading order (NLO) of the two-particle irreducible (2PI) effective action [28,29]. This self-consistent expansion scheme is uniform in time, resumming secular terms such that it can be applied also to study late-time dynamics [30].
A full description of general out-of-equilibrium dynamics can be based on both commutator and anticommutator expectation values of products of field operators. Here we consider the statistical (F ) and the spectral two-point function (ρ), where {. , . } describes the anti-commutator and [. , . ] the commutator that is applied to the field operatorφ a (x). The evolution equations at NLO in a 1/N expansion have been derived in Ref. [28]. For their solution we have to supply initial conditions at time x 0 = t = 0. Here we employ initial conditions for the spatially homogeneous system with a vanishing macroscopic field, i.e. φ a (t = 0) = 0 and ∂ tφa (t = 0) = 0 for all field components. By virtue of the O(N ) symmetry, the field expectation value then remains zero at all times. This allows us to write F ab (x, y) = F (x, y)δ ab , ρ ab (x, y) = ρ(x, y)δ ab , (3) such that we only need to consider the diagonal elements F (x, y) and ρ(x, y).
The initial conditions at t = 0 for the spectral function are fixed by its anti-symmetry, ρ(t, t, x, y) = 0, and the equal-time commutation relations of the bosonic quantum theory, For the statistical function at initial time we consider Gaussian correlations, which in spatial Fourier space with momentum p are parameterized as Here f p denotes an initial particle number distribution with dispersion ω p . The latter is set at initial time by ω p = |p| 2 + m 2 0 in the limit of a vanishing initial mass m 0 → 0 + . Following Ref. [25], we consider the initial condition of a highly occupied system of particles with distribution function where the characteristic momentum Q 0 sets the initial scale, with occupancy parameter n 0 , and Heaviside function Θ. Since the initial occupancy is inversely proportional to the coupling, this represents a non-perturbative problem even for λ 1. Since the NLO approximation for the dynamics is non-Gaussian, higher correlations will build up during the time evolution.
Starting from this over-occupied initial situation, we compute the time evolution of the system by numerically solving the coupled NLO evolution equations for the statistical and spectral correlation functions [28]: . The NLO approximation for the statistical and spectral components of the self-energies, Σ F (x, z) and Σ ρ (x, z), entails a geometric series summation of the correlation functions F (x, y) and ρ(x, y), which can be conveniently expressed in terms of the summation functions I F (x, y) and I ρ (x, y): F (x, y)I ρ (x, y) + ρ(x, y)I F (x, y) .

(8b)
The summation functions contain as building blocks the 'one-loop' self-energies Π F (x, y) and Π ρ (x, y) with The latter spectral component is directly related to the retarded self-energy We emphasize that the large N approximation only involves a systematic power counting of factors of 1/N . As a consequence, the terms neglected are suppressed by an additional power of 1/N . Since this is not an expansion in powers of the coupling λ, even non-perturbative situations far from equilibrium such as the over-occupied initial conditions (6) can be addressed. Moreover, the NLO contributions are essential: at LO (N → ∞) the right hand side of the evolution equations (7) would vanish since the self-energies (8) are zero at that order.
The numerical results presented in this work employ N = 4 and occupation parameter n 0 = 100. We study a wide range of coupling parameters λ = 0.01, 0.10, 1.0, 2.0. If not stated otherwise, λ = 1.0 is used. In the following, all quantities are given in units of the characteristic initial scale Q 0 and are stated as dimensionless numbers. The equations of motion are discretized on a lattice with temporal step a t and spatial spacing a s . The results shown are obtained from computations using N s = 500 points on a spatial grid with a s = 0.75, corresponding to an ultraviolet (UV) momentum cutoff of Λ UV = 4.19 and an infrared (IR) cutoff of Λ IR = 0.0084. We have checked that all relevant results are insensitive to both IR and UV cutoffs. In order to efficiently approach late times, the time step a t is tuned to be as large as possible while checking numerical convergence to runs with smaller time steps. For the presented numerical results a time step of a t = 0.3 is used.

III. UNIVERSAL SCALING DYNAMICS OF EQUAL-TIME CORRELATIONS
We first compute the non-equilibrium evolution of the statistical correlation function (2a) at equal times, for which the spectral function (2b) vanishes because of its anti-symmetry. A similar calculation for a larger class of initial conditions (with non-zero initial field expectation value) has been done in Ref. [25] pointing out the independence of rescaled results on system parameters, such as the value of the coupling λ or details of the initial conditions, in an emergent universal regime associated to a non-thermal fixed point [3,4]. Exploiting this insensitivity to initial condition details, we restrict ourselves to the symmetric regime with initial conditions described in Sec. II. This reduces the numerical efforts and we do not have to employ an adaptive grid size as done in Ref. [25], which allows us for the first time to accurately extract the universal self-similar scaling exponents and scaling function from the complete NLO evolution equations (7) -(10) by numerical computations. A calculation of the far-from-equilibrium scaling exponents and function has so far only been done using additional assumptions, such as a quasi-particle ansatz for an effective kinetic description at the non-thermal fixed point in Refs. [4,15,16], or based on classical-statistical field theory approximations in Ref. [4] in the weak-coupling limit. We emphasize that our approach is not restricted to weak couplings and includes genuine quantum effects at NLO in the large-N expansion.

A. Particle distribution
So far, the phenomenon of scaling has mainly been discussed in terms of a particle number distribution function, whose time-dependence we extract from the twopoint correlation function as [28,30] f (t, |p|) Similarly, one may define a time-dependent effective dispersion such that we have at initial time f (t = 0, |p|) = f p and ω(t = 0, |p|) = ω p according to (5). Starting from the over-occupied initial state at t = 0, the non-equilibrium evolution leads to a redistribution of both particle number and mode energy. In the upper plot of Fig. 1, we show the time evolution of the distribution function f (t, p). The effective mass m eff , which is given by the approximately time-independent value of the dispersion at zero momentum as analyzed in Sec. IV B, is also indicated. For t 1500, the dynamics slows down considerably and we analyze in the following whether the system becomes self-similar. We concentrate on the nonperturbative behavior for sufficiently low momenta, and refer for the analysis of the perturbative high-momentum properties to Refs. [6,25]. For a self-similar time evolution the distribution obeys the scaling property with scaling exponents α, β and time-independent scaling function f S . Therefore, in the scaling regime t −α f (t, |p|) does not depend on time and momentum separately but only on the product t β |p| for a set of exponents α and β. Universality implies that the shape of the distribution function and the values of the exponents do not depend on the microscopic model parameters, such as the value of the coupling λ, which is discussed in the following. As shown in the lower plot of Fig. 1, the distribution function at different times can be rescaled such that the curves of t −α f (t, |p|) as a function of t β |p| lie on top of each other for lower momenta.
Numerically, the scaling exponents are obtained by comparing the distribution function at some reference time t ref with several earlier times t, where we perform comparisons within a time window of ∆t w = 720. For details on the method and the employed error estimation we refer to Appendix A. In Fig. 2 we show our results for the scaling exponents α and β for different couplings where the data points are binned for the plots. Both exponents approach approximately constant values given in Table I. For later comparison, we also present values for the exponent α λ as defined in the table caption.
The exponents for all couplings studied here agree within errors with each other. Furthermore, they are consistent with the results found in Ref. [4] using classicalstatistical lattice simulations, confirming that statistical fluctuations dominate over genuine quantum fluctuations in this highly occupied regime. As pointed out in Refs. [4,15,16], the values of these exponents coincide with those of the corresponding non-relativistic scalar model, since infrared momenta below m eff of the relativistic theory behave non-relativistically. Within errors, we also find α = dβ for d = 3 spatial dimensions such that p f (t, |p|) = t α−dβ q f S (|q|) is approximately conserved, reflecting a transport of particles towards lower momenta for the β > 0 observed.  I. Exponents α and β obtained from the scaling analysis of the particle distribution f (t, |p|). Using these values we also give α λ = 2(β − α) for later comparison, which is the expected scaling exponent of the four-vertex as introduced in Sec. III C. Our results are shown for different values of the coupling parameter λ.

B. Mode energy
In addition, we consider the time-dependent mode energy ε(t, |p|), which is given at NLO in the 1/N expansion by [25] where the approximation employed for the last line is only used here to analyze the quasi-particle content of the dynamics. 1 The momentum sum of the mode energy (15) is equal to the conserved total energy density at NLO, ε(t) = d 3 p/(2π) 3 ε(t, |p|). In our simulations, we checked that ε(t) is conserved at the level of 1 % accuracy for the times under consideration. The scaling analysis for the mode energy density employs where we anticipate that the scaling exponents are the same as for the distribution function, and ε S denotes the energy scaling function. The upper graph of Fig. 3 shows the time evolution of ε(t, |p|), while the lower one displays the rescaled quantity t −α ε(t, |p|) as a function of t β |p| employing the same values for the exponents as obtained from the particle number distribution. One observes that in the scaling regime all curves collapse to a single one in the infrared to very good accuracy.
In order to illustrate this agreement, we analyze the approximate quasi-particle expression (16). In Fig. 4 the mode energy at NLO (15) is compared to the quasiparticle expression (16). One observes rather good agreement, in particular in terms of the scaling properties.

C. Scaling of the effective coupling
In this section we analyze the scaling properties of the four-vertex in an approximation based on the large-N expansion to NLO. Since we are interested in a slowly evolving self-similar scaling regime, we may simplify the computation considerably by relying on a derivative expansion in time. More precisely, at lowest order in derivatives it is convenient to consider the (on-shell) effective coupling in Fourier space given by [4,14,31] This effective coupling approximates the full four-vertex, which at this order receives its momentum dependence from the resummed geometric series underlying the NLO approximation [28,29]. Moreover, the effective coupling is directly related to the summation functions (9) according to [4,14,31] λ eff (t, |p|) In the following this is used to numerically compute the time-dependent effective coupling, where I F (t, |p|) and Π F (t, |p|) are obtained from the expressions (9) and (10) in spatial Fourier space evaluated at equal times. Using the scaling property (14) of the distribution function in the expression for the self-energies (10) and (8) in the non-relativistic regime, one expects [4,15,16] with scaling function λ eff,S and coupling scaling exponent directly related to the occupation number exponents α and β in d spatial dimensions. The upper panel of Fig. 5 shows the momentum dependence of the effective coupling at different times. Remarkably, the effective coupling drops over several orders of magnitude in the non-perturbative infrared regime, where the occupation number grows larger with time. This non-equilibrium phenomenon of a dynamically reduced four-vertex counteracts the dramatic Bose enhancement from the very high occupancies in the infrared, which would otherwise lead to faster and faster dynamics for growing occupancies. The corresponding phenomenon has recently also been experimentally demonstrated in an atomic Bose gas far from equilibrium [13]. The lower graph of Fig. 5 demonstrates that the numerical data for the effective coupling collapses rather well to a common curve in the infrared momentum range when rescaled accordingly. As a check, we also determine the scaling exponents α λ and β defined by (20) directly from our numerical data without assuming the scaling relation (21), see Appendix A for details on the method. The binned data obtained for α λ in this way is plotted in Fig. 6, where we also show the results for the exponent if computed according to (21) using the scaling exponents α and β. We observe a good agreement of the data although the errors for the analysis using λ eff are larger and fluctuate stronger, which is reflected in enhanced statistical errors as discussed in Appendix A. The asymptotic values approached are presented in Table II.

IV. SPECTRAL AND STATISTICAL CORRELATIONS AT UNEQUAL TIMES
For the study of correlation functions at different times t and t , it is convenient to rephrase the time-dependence in terms of Wigner coordinates, employing the central time τ = (t + t )/2 and the relative time ∆t = t − t . For the spatially homogeneous system, we then denote the two-point functions using the new temporal coordinates as F (τ, ∆t, |p|) and ρ(τ, ∆t, |p|).
In order to study the frequency spectrum of these statistical and spectral functions, we consider a finite-range Fourier transformation of the Wigner space propagators  with respect to the relative time ∆t, where the factor of i is introduced such that both ρ(τ, ω, |p|) and F (τ, ω, |p|) are real. To ease the notation, we will neglect the tilde for the real spectral function in frequency space,ρ(τ, ω, |p|), in the following having in mind the extra factor of i in its definition. The integrals with respect to the relative times in (22) are fundamentally restricted by ±2τ for the initial value problems with t, t ≥ 0. Moreover, it is sufficient to present the propagators for positive ∆t or ω, since the statistical (spectral) function is (anti-)symmetric in ∆t and hence ω. Effects resulting from the finite-time boundary vanish in the limit τ → ∞ and are discussed further in Sec. IV A. Details on the numerical treatment are presented in Appendix B.
We emphasize that in the interacting quantum field theory the spectral function ρ(τ, ω, |p|) is, in general, not positive for ω > 0. In contrast, simple quasi-particle descriptions typically rely on positivity, such as the freefield spectral function with a positive particle peak at the frequency that is equal to the mass of the particle. We can check whether positivity is approximately realized by comparing to the corresponding expression with the absolute value of the integrand of (23), Fig. 7 compares this modified expression to the sum rule as a function of spatial momentum. Here the data for the Fourier transformed spectral function is computed using a discrete Fourier transformation, see Appendix B for details. For sufficiently large momenta, we indeed observe agreement as expected based on the validity of perturbation theory in this regime. However, at low momenta significant deviations are seen, indicating that there is no simple quasi-particle spectral function in the deep infrared.
To further analyze this behavior, we distinguish three momentum regimes which we identify with the help of the . Shown is also the particle distribution function from which we identify the time-dependent scales K(t) and Q(t).
particle distribution function defined in (12), as shown in Fig. 7. There we indicate a characteristic time-dependent momentum scale K(t), where the distribution function shows maximum positive curvature in the scaling regime, along with Q(t) defined by the scale of maximum negative curvature. These two scales are used to identify the (I) infrared "plateau" regime, |p| K(t), (II) infrared "power-law" regime, K(t) |p| Q(t), (III) high-momentum perturbative regime, |p| Q(t).
We emphasize that both momentum regimes (I) and (II) constitute the inverse particle cascade and are characterized by the same universal scaling exponents. In particular, K(t) ∼ t −β , which can be inferred from the upper graph of Fig. 8. The scale Q(t) also evolves towards the infrared as can be seen in lower plot of Fig. 8. It appears not to be described by a simple power law, as the evolution becomes faster with time.
In the following, we discuss the differences of the unequal-time two-point functions in the three momentum regimes in detail.

A. Violations of the fluctuation-dissipation relation
In contrast to our non-equilibrium situation considered, thermal equilibrium is time-translation invariant, such that correlation functions only depend on relative coordinates and there is no τ -dependence. Moreover, the statistical and spectral functions in thermal equilibrium, F (eq) and ρ (eq) , are related by the fluctuation-dissipation theorem, see e.g. [30], according to This implies that the ratio F (eq) (ω, p)/ρ (eq) (ω, p) is independent of spatial momentum p and determined by the frequency-dependent Bose-Einstein distribution f BE (ω) = (e βω − 1) −1 with inverse temperature β in the absence of conserved charges. Out of equilibrium, the statistical and spectral correlation functions are linearly independent in general, but one may hope to find some generalized fluctuationdissipation relation where the Bose-Einstein distribution is replaced by a time-dependent distribution function. Such a generalized relation provides, for instance, the basis for standard kinetic descriptions.
In order to study the non-equilibrium frequency space two-point functions, we consider the Fourier Wigner space propagators F (τ, ω, |p|) and ρ(τ, ω, |p|) for times τ at which the system has reached the scaling regime. Details on the numerical computation of the frequency- space propagators can be found in Appendix B. As discussed above, we expect the spectral function to describe a quasi-particle excitation spectrum for sufficiently large momenta of regime (III) and maybe (II). Fig. 9 shows the numerical results for (f (t = τ, |p|) + 1/2)ρ(τ, ω, |p|) and F (τ, ω, |p|) as a function of ω at some time τ = 2250 in the scaling regime for three different momenta in the ranges (I), (II), and (III) as indicated in the inset. As anticipated, from the upper graph one observes that both expressions are clearly different in the deep infrared regime (I). The middle graph shows that in regime (II) both quantities become similar in shape, however, the respective peaks are shifted relative to each other and the height is not the same. In contrast, F and (f + 1/2)ρ have the same Breit-Wigner shape for the considered momentum in regime (III) although the amplitudes do not fully agree yet. As the momenta become larger, we checked that this agreement gets more accurate. This establishes a well-defined generalized fluctuation-dissipation relation in terms of the nonequilibrium distribution function f (τ, |p|) in the highmomentum regime.
We now analyze the behavior of the spectral function in regime (I) in more detail. The upper graph of Fig. 7 reveals that the spectral function has the shape of an enveloped oscillation in this regime. We emphasize that this behavior is insensitive to the numerical discretization, i.e. to changes of the time-step size, the spatial grid step or volume size. Instead, it is related to the initial time boundaries at t ≥ 0 and t ≥ 0, which enter the integration boundaries in (22). In order to understand this effect better, it is helpful to consider the Wigner space propagators F (τ, ∆t, |p|) and ρ(τ, ∆t, |p|), i.e. without the Fourier transformation with respect to relative time ∆t. Both functions oscillate in ∆t with a frequency that in general can depend on both τ and |p|. The envelopes of these oscillations are plotted in Fig. 10, where the insets show the actual oscillations. The oscillation frequency can be used to define a dispersion relation ω(τ, |p|), which we call Wigner dispersion in order to distinguish it from the effective dispersion defined in (13). Since the oscillation frequencies observed are practically constant during the self-similar time evolution, the Wigner dispersion is quasi-stationary.
In principle, one could obtain Wigner dispersions for the statistical and spectral functions separately. This seems necessary at first sight, since the behavior of F and ρ is clearly different in the infrared, as seen in Fig. 9. However, we find that the relative difference (ω F −ω ρ )/ω F is smaller than 0.1 % such that for all practical purposes the Wigner dispersions of the statistical and spectral functions are treated as being equal.
While the oscillation frequencies of F (τ, ∆t, |p|) and ρ(τ, ∆t, |p|) are practically equal, the envelopes of the oscillations can differ significantly from each other in the infrared. Fig. 10 compares the envelopes of F and ρ, which are plotted on a logarithmic scale against the relative time ∆t, in the different momentum regimes. In (II) and (III), one observes that both F and ρ decay exponentially in ∆t, which leads to well-defined quasiparticle peaks in frequency space as seen in Fig. 9. In regime (I), however, the spectral function grows towards the initial-time boundaries whereas the statistical function is damped stronger than exponentially (upper plot of Fig. 10).
To analyze this further, it is helpful to consider the simplified case of strictly exponentially damped oscilla- tions. Omitting the τ -dependence for the moment, we approximately write where γ p is the damping constant, F 0 denotes the ampli-tude of the statistical function at ∆t = 0 and the factor ω −1 p in (26b) ensures that the spectral function suffices the commutation relation (4). For this ansatz, calculating the Fourier transform with respect to the relative time ∆t analytically yields the frequency space propagators where the latter corresponds to the relativistic Breit-Wigner function [32]. On-shell, where (ω 2 + ω 2 p ) ≈ 2ω 2 , the statistical function is also described by a Breit-Wigner shape. The damping constant γ p determines the width of the quasi-particle peak.
In regime (III), F and ρ decay with the same damping constant and consequently have the same Breit-Wigner shape in Fourier space, see Fig. 9. We can determine the parameters γ p , ω p and F 0 by fitting our data to either (26) or (27). Since the Fourier Wigner space propagators are numerically obtained by a discrete Fourier transformation of the Wigner propagators, both methods are equivalent up to numerical uncertainties.
When looking at regime (II) the statistical function has a slightly larger damping constant than the spectral function. The main difference, however, appears at large relative times ∆t, i.e. at the temporal boundaries t = 0 and t = 0, where the spectral function grows somewhat while the statistical function decays even faster. This means that spectral correlations at large time-separations are enhanced whereas statistical correlations are suppressed. As a consequence, the Fourier Wigner propagators shown in Fig. 9 reveal different shapes.
For exponentially decaying Wigner space propagators, the effects from finite integration bounds are negligible at sufficiently late times, where τ is large. However, in regime (I) the spectral function does not decay towards the initial-time boundaries whereas the statistical function drops very quickly. Consequently, only the statistical function is described by a peak in Fourier Wigner space, however, the peak cannot be described by a Breit-Wigner function since no exponential decay is involved. In contrast, the growth of the spectral function at large relative times is sharply cut off by the time boundary |∆t| < 2τ appearing in initial value problems. Hence, the Wigner transformation gives rise to fast oscillations of the spectral propagator in frequency space, as seen in the upper plot of Fig. 9. The oscillation frequency is determined by the integration range. The spectral propagator ρ(τ, ω, |p|) oscillates in ω with frequency 2τ and has an envelope that is peaked around ω = ω p .
The differences of the Wigner space propagators between the three momentum regimes can be visualized more clearly when looking at the contour plots shown in Fig. 11, where the relative time ∆t labels the horizontal axis, the central time τ the vertical axis, and the envelopes of F (τ, ∆t, |p|) and ρ(τ, ∆t, |p|) are encoded in the color scheme. The initial time bounds t, t ≥ 0, equivalent to −2τ ≤ ∆t ≤ 2τ , are marked by the gray lines. Each horizontal slice corresponds to a fixed central τ and can be Wigner-transformed with respect to ∆t in order to obtain the corresponding Fourier Wigner space propagators. The propagators shown in Fig. 9 correspond to the latest available time slices at τ = 2250.
In accordance with the behavior of the particle distribution function discussed above, the amplitude of F decreases by several orders of magnitude when going from low to high momenta (left plots of Fig. 11 from top to bottom). In contrast, the spectral function ρ does not differ much in amplitude since it is normalized according to the sum rule (23).
The contour plots in Fig. 11 visualize how the envelopes of F and ρ evolve with time τ . Due to the self-similar time evolution, going towards later times τ is equivalent to moving to larger momenta |p|. In the infrared momentum regime, the increase of the amplitude of the spectral function towards the initial time bounds declines with evolving time τ . For higher momenta, the decay rate increases as the time τ evolves. Hence, the effect of high amplitudes at large ∆t vanishes at sufficiently late times, where the time scale is larger for small mo-menta. If one considers, for instance, some fixed momentum |p * | < K(t) in region (I), then during the time evolution the characteristic momentum K(t) ∼ t −β moves towards the infrared such that at some later time t > t one finds |p * | > K(t ). Since the momentum moved from region (I) into region (II), where the Wigner space propagators decay exponentially, no boundary effects occur. In that sense, due to the self-similar time-evolution it is always possible to wait long enough to overcome the initial-time boundary effect for a given momentum.

B. Dispersion and effective mass
Motivated by the results of the last section, we may consider a "Wigner" particle distribution function f (τ, |p|) obtained from which is a priori different from the "equal-time" definition employed in (12). From the upper panel of Fig. 12 one observes that both definitions are in good agreement with each other. Small deviations in the high-momentum range can be cured in the limit a t → 0. Similarly, we introduced the effective dispersion (13) and the Wigner dispersion (see Sect IV A). Both definitions for the dispersion relations can actually be fitted rather well to a relativistic dispersion, with effective mass m eff that incorporates quantumstatistical fluctuations. Our numerical computations show that they agree well with each other, as seen from the lower graph of Fig. 12. Although both the dispersion and the effective mass are in general time-dependent, they turn out to be practically constant in time in the scaling regime. Because the extraction of the effective mass from the Wigner dispersion turns out to be numerically more stable than using the equal-time dispersion, the values cited in the text are obtained from the former.
In particular, we find from fitting the data to the relativistic dispersion relation (29) the following values for the effective mass for different values of the interaction parameter, 0.638 ± 0.005 λ = 0.01 0.641 ± 0.005 λ = 0.10 0.677 ± 0.007 λ = 1.00 0.716 ± 0.008 λ = 1.00 (30) where the error indicates that m eff is not exactly timeindependent during the self-similar evolution. The presence of an effective mass explains why the infrared exponents found in Sec. III are close to the predictions for a non-relativistic theory [4]. As can be seen in the upper plot of Fig. 1, the momenta of the whole infrared regime are much smaller than the effective mass (note the logscale) and thus effectively non-relativistic.

V. CONCLUSION
We analyzed the far-from-equilibrium scaling properties of equal-time and unequal-time correlation functions in scalar quantum field theory. The numerical results are obtained from a fully self-consistent large-N expansion to NLO. Our results for scaling exponents and scaling functions are in agreement within errors with previous weak-coupling estimates for equal-time correlations using effective kinetic theory or classical-statistical field theory. We find these universal results for a wide range of couplings even beyond the weak-coupling regime. Moreover, we have established the self-similar behavior of the strongly momentum-dependent effective coupling.
The computation of the unequal-time spectral and statistical functions allowed us to observe the validity of a generalized fluctuation-dissipation relation in the perturbative regime at high momenta, while we demonstrated that significant violations occur in the non-perturbative infrared. We identified a characteristic deep-infrared regime, where the corresponding distribution function approaches a "plateau". In this regime the spectral function does not decay as a function of relative time, which leads to an enhanced sensitivity to initial times and the absence of positivity for frequencies ω > 0. The statistical function in this regime shows a characteristic peak structure with significant deviations from a Breit-Wigner form.
Our results give unprecedented insights into the nonperturbative nature of collective excitation far from equilibrium. While many features of the system are indeed seen to become universal near a non-thermal fixed point also beyond the weak-coupling regime, unequaltime properties encoded in spectral functions can reveal intriguing properties such as an enhanced sensitivity to initial times. Since our results are based on a large-N expansion, it would be very interesting to also analyze the behavior of weakly-coupled systems with small N using real-time lattice simulations along the lines of Refs. [26,27].
which allows us to compare the distribution at a reference time t ref with several distributions at earlier times t i , i = 1, . . . , n, within a given time window ∆t w . For the correct set of scaling exponents, a perfect self-similar time evolution implies that the rescaled distribution becomes independent of time, i.e.
for times t i and t ref within the scaling regime. In practice, the numerically obtained distribution functions deviate from the perfect scaling behavior. The scaling exponents are determined by minimizing these deviations, which we quantify by We sum over the n comparisons between the reference time t ref and earlier times t i < t ref . The integration over d log p enhances the low-momentum range.
In order to compute the differences ∆ i (p), the numerical data of the distribution functions is interpolated using a cubic spline provided by the Python SciPy library. The integration range is chosen dynamically in terms of the time-dependent characteristic momentum scales K(t) and Q(t). We include 95 % of the momentum range [log Λ IR , log K(t ref )] and 92 % of the momentum range [log K(t ref ), log Q(t ref )], which ensures that the high-momentum range is excluded from the fitting procedure for the analysis of the infrared fixed point.
From the minimization of (A3) at different reference times we obtain scaling exponentsᾱ andβ at which χ 2 is minimal. Here, we use optimization routines of Python SciPy library for the minimization procedure. Thereby we observe that χ 2 (ᾱ,β) decreases with time and converges to a constant value when the system is approaching the non-thermal fixed point.
The error for the values ofᾱ andβ is estimated by the marginal likelihood functions where the normalization constants N α,β are chosen such that dα W (α) = dβ W (β) = 1. Fitting W α,β to a Gaussian distribution allows us to estimate the error of the exponents by the standard deviation of the Gaussian function. We refer to this uncertainty as the fit errors denoted by ∆α fit and ∆β fit , which determine the errorbars in the plots shown in Figs. 2 and 6. In the analysis of the distribution function f , the fit errors of both exponents α and β decrease with the time evolution and approach asymptotic values.
The asymptotic values of the exponents are obtained by averaging over the values obtained at late reference times. For the analysis of the main text we compute the mean over a time window ∆t av ≥ 600. The corresponding standard deviation is used to quantify the statistical errors ∆α stat and ∆β stat . These errors reflect how strong the exponents fluctuate during the time-evolution. For the values presented in Tables I and II, we provide the error ∆α = (∆α fit ) 2 + (∆α stat ) 2 , and accordingly for the other exponents.
In the scaling analysis presented, n = 4 comparisons within a time window of ∆t fit = 720 were used. The fit and statistical errors of our simulations are shown in Tables III and IV. While the analysis for f yields very small statistical errors, the exponents obtained for both the effective coupling λ eff as well as the mode energy ε reveal strong fluctuations represented by much larger statistical errors.
We have checked that the values for the exponents are insensitive to the numerical discretization (time step size and cutoffs). Naturally, there exists a dependence on the parameters of the method used, such as the number of comparisons n, the time windows ∆t w and ∆t fit , and the momentum range used to compute (A3). We checked that the values obtained are not sensitive to n and ∆t w .

Appendix B: Wigner transformation
In this appendix we present the methods used in order to compute the Wigner transformed spectral and statistical functions according to (22). In order to compute the Wigner transform f (ω) of a temporal signal f (t), we need to numerically evaluate integrals of the form where we use t instead of ∆t here for notational convenience. We are interested in signals oscillating with a given frequency ν that are enveloped by some function g(t) which can be written as f ± (t) = g(t) 1 2 (e iνt ± e −iνt ) , where the relative sign determines whether the signal is symmetric or antisymmetric. In the following, we present the methods that we employ in order to compute (B1) for such signals.
If g(t) decays sufficiently strong, i.e. if it becomes sufficiently small at the boundaries ±2τ , effects from the finite integration boundaries are negligible. In this case, the Wigner transform can be computed using a discrete Fourier transformation (DFT), where N t denotes the number of data points, and ω m and t m are the discretized frequency and time, respectively, with m = 0, . . . , N t − 1. In our simulations, the DFT is implemented using the FFTW library [33].
For an exponentially decaying envelope, which is described by g(t) ∼ exp(−γ|t|) with some decay constant γ, the DFT can be compared to the analytically calculated Wigner transform with finite as well as infinite integration boundaries. For the parameters relevant in our simulations, we confirmed that the DFT of the exponentially decaying Wigner functions appearing in the high momentum regime yields accurate results. In addition, we checked the applicability of the DFT for other momentum ranges by comparing the decaying behavior of the propagator functions with exponential decays. Thus, for the analysis shown in the main text we use the DFT to compute the Wigner transform of the statistical function in all momentum ranges and the spectral function for medium or high momenta.
In general, the finite integration boundaries in (B1) lead to oscillations of the Wigner transform f (ω). These are particularly relevant if the envelope function g(t) increases towards larger |t|. In this case, we determine the envelope g(t) by a polynomial fit and the oscillation frequency ν using the Lomb-Scargle periodogram implemented in the Python SciPy library. We can then analytically compute the Wigner transformation of (B2) for the relevant parameters. This method is employed for computing the spectral function at small momenta. However, when comparing the Wigner transform at different spatial momenta |p| as in Fig. 7, we employ the DFT for both F and ρ. with an ultraviolet cutoff given by the discretization.