Basin of attraction for turbulent thermalization and the range of validity of classical-statistical simulations

Different thermalization scenarios for systems with large fields have been proposed in the literature based on classical-statistical lattice simulations approximating the underlying quantum dynamics. We investigate the range of validity of these simulations for condensate driven as well as fluctuation dominated initial conditions for the example of a single component scalar field theory. We show that they lead to the same phenomenon of turbulent thermalization for the whole range of (weak) couplings where the classical-statistical approach is valid. In the turbulent regime we establish the existence of a dual cascade characterized by universal scaling exponents and scaling functions. This complements previous investigations where only the direct energy cascade has been studied for the single component theory. A proposed alternative thermalization scenario for stronger couplings is shown to be beyond the range of validity of classical-statistical simulations.

Different thermalization scenarios for systems with large fields have been proposed in the literature based on classical-statistical lattice simulations approximating the underlying quantum dynamics. We investigate the range of validity of these simulations for condensate driven as well as fluctuation dominated initial conditions for the example of a single component scalar field theory. We show that they lead to the same phenomenon of turbulent thermalization for the whole range of (weak) couplings where the classical-statistical approach is valid. In the turbulent regime we establish the existence of a dual cascade characterized by universal scaling exponents and scaling functions. This complements previous investigations where only the direct energy cascade has been studied for the single component theory. A proposed alternative thermalization scenario for stronger couplings is shown to be beyond the range of validity of classical-statistical simulations.

I. INTRODUCTION
A quantum many-body system in thermal equilibrium is independent of its history in time and characterized by a few conserved quantities only. Therefore any thermalization process starting from a nonequilibrium initial state requires an effective loss of details about the initial conditions at sufficiently long times. It was pointed out previously that an effective partial memory loss of initial condition details can already be observed at earlier stages for the nonequilibrium unitary time evolution in quantum field theory [1].
An extreme case occurs if at some early stage the nonequilibrium dynamics becomes self-similar. This amounts to an enormous reduction of the sensitivity to details of the underlying theory and initial conditions. The time evolution in this self-similar regime is described in terms of universal scaling exponents and scaling functions. While normalizations of the latter depend on model parameters such as couplings, masses and initial conditions, their functional form is universal. This has the important practical consequence that the nonequilibrium dynamics for whole ranges of different model parameters and initial conditions can be grouped into universality classes. The self-similar time-evolution within a given universality class of different models can then be mapped onto each other by simple rescalings. The universal properties are described in terms of nonthermal fixed points [2] in renormalization group theory [3]. While the notion of thermal fixed points describing different universality classes is well established for systems in thermal equilibrium, the corresponding classification for nonequilibrium scaling phenomena is a rapidly progressing research topic. * K.Boguslavski@thphys.uni-heidelberg.de Universal behavior far from equilibrium has been predicted for systems ranging from early-universe inflaton dynamics to table-top experiments with cold atoms. In these examples, attractor solutions with self-similar scaling behavior are associated with the phenomenon of Kolmogorov wave turbulence [4,5] and strong turbulence with fluid-like behavior [2,[6][7][8] which leads to Bose condensation far from equilibrium [9][10][11]. A turbulent thermalization mechanism has also been predicted for non-Abelian gauge theories using classical-statistical simulations in a fixed box, where different groups obtain consistent results [12][13][14].
Recently, the existence of a nonthermal fixed point in longitudinally expanding non-Abelian plasmas was demonstrated using classical-statistical lattice simulations in the weak coupling limit [15,16]. It was argued that heavy-ion collision experiments at sufficiently high energies can provide the relevant initial conditions in the basin of attraction of this nonthermal fixed point. At sufficiently high energies, the colliding nuclei are described in the Color Glass Condensate framework [17]. The dynamics of the nonequilibrium Glasma [18] created in such a collision is that of gluon fields with typical momentum Q and a weak gauge coupling α s (Q). Since the characteristic occupancies ∼ 1/α s (Q) are large, the gauge fields are strongly correlated even for small gauge coupling and the Glasma exhibits classical-statistical dynamics.
However, in Ref. [26] it has recently been argued that by exceeding a rather small value of the gauge coupling α = g 2 /4π with g ∼ 0.5, the dynamics changes dramatically compared to previous weak-coupling estimates, us-ing small initial Gaussian distributed fluctuations superimposed to the large classical macroscopic field ∼ O(1/g) as derived in Ref. [27]. Moreover, similar findings of pressure isotropization for condensate driven initial conditions have been reported in nonequilibrium classicalstatistical scalar field theory with expansion [28].
Even without expansion, a different thermalization scenario in the presence of large initial condensates was discussed using classical-statistical simulations once the self-coupling λ reaches a certain strength [29]. In contrast, in the weak-coupling limit, the same initial conditions are known to lead to high occupancies ∼ 1/λ after a time t ∼ Q −1 log(λ −1 ) where the characteristic scale Q is determined by the energy density of the system [30,31]. The subsequent evolution exhibits the phenomenon of turbulent thermalization [4] reminiscent of the classical-statistical non-Abelian gauge theory results in the weak coupling limit [12][13][14]. The important advantage for scalars is that powerful resummation techniques exist for the description of the entire evolution in quantum field theory. It has been demonstrated that classical-statistical descriptions have a well understood range of validity for which they accurately describe the dynamics of the underlying quantum theory [2,[32][33][34].
In this work, we revisit the question of nonequilibrium dynamics from condensate driven initial conditions for the single component theory employed in Ref. [29]. We first investigate the dynamics in the weak coupling limit. In this limit, we observe turbulent thermalization, where we point out that 1. the inverse particle cascade, which has so far only been established for O(N ) symmetric theories with N ≥ 2 [2,8,35], also occurs for the single component theory.
This finding directly complements the investigations of Refs. [4,5], where the direct energy cascade has been studied for the single component theory while the phenomenon of a dual cascade was not considered. We then explore the basin of attraction for the nonthermal fixed point by comparing two generic classes of initial conditions: 2. Condensate driven initial conditions with large field and small (vacuum) fluctuations, 3. fluctuation dominated initial conditions with small (zero) field and large fluctuations, and confirm [4,5] that they lead to the same universal behavior beyond a time scale ∼ Q −1 log(λ −1 ) for the whole range of couplings where the classical-statistical approach is valid. Thus the observed differences in thermalization scenarios that have been proposed are not caused by differences in the initial conditions. We finally investigate 4. the range of validity of classical-statistical simulations and explain why the results of Ref. [29] do not lie in this range.

II. SCALAR MODEL AND INITIAL CONDITIONS
There is a restricted class of problems where the dynamics of bosonic quantum fields can be accurately mapped onto a classical-statistical problem. The most intuitive criteria for this can be formulated in situations where a kinetic description in terms of quasi-particle excitations is applicable. The system exhibits classical dynamics whenever the typical occupation numbers per mode are much larger than unity, If occupation numbers fall below unity, quantum processes will dominate the dynamics. This is clearly seen in a Boltzmann transport framework where classical scattering processes are sub-leading to quantum ones for occupation numbers smaller than unity [36,37]. It is also possible to formulate more general criteria, which do not rely on a quasi-particle picture, in the Schwinger-Keldysh formalism of nonequilibrium quantum field theory [33,[38][39][40]. This 'classicality condition' is met whenever anti-commutator expectation values for typical bosonic field modes are much larger than the corresponding commutators [33,39]. Stated differently, this concerns the large field or large occupancy limit, which is relevant for important phenomena such as nonequilibrium instabilities or wave turbulence encountered in our study. The classicality condition has been verified in detail by comparing quantum to classical-statistical results in the context of scalar quantum field dynamics [2,32,33,39] and coupled to fermions [41,42].
We emphasize that the condition for a system to exhibit classical dynamics is in general time dependent. In particular, the approach to complete thermal equilibrium is not accessible within the classical-statistical framework. As a consequence of the Rayleigh-Jeans divergence, the classical thermal state is only well defined for an ultraviolet cutoff Λ, which may be implemented by a lattice regularization. In its range of validity, observables computed from classical-statistical simulations are insensitive to the Rayleigh-Jeans divergence. Thermal equilibrium is a genuine quantum state which cannot be reached within classical-statistical field theory. Nevertheless, the classical-statistical regime may extend over large times such that the quantitative properties of the nonequilibrium quantum evolution are accurately described within the classical-statistical approach.
The time scale t quant for entering the quantum regime is not a universal quantity and depends in general on the properties of the initial state as well as the dynamics in the classical regime. For the considered cases of large initial fields or large initial fluctuations at weak coupling, this time scale is parametrically given by t quant ∼ Q −1 λ −5/4 in our case [5]. The range of validity of classical-statistical techniques in time is thus naturally confined to weak couplings. In general the use of classical-statistical methods at large couplings requires great care since genuine quantum effects may dominate the dynamics already at rather early times. This will be investigated in detail below after we present the physical weak coupling results.
We study the dynamics of a massless real scalar field theory with quartic interaction. The classical action is given by summing over µ = 0, 1, 2, 3 in Minkowski space-time.
The classical equation of motion δS/δϕ = 0 for the singlecomponent field reads Following Ref. [29], we consider Gaussian initial conditions. They are fully determined by the initial oneand two-point correlation functions of the fields ϕ(x) and π(x) = ∂ϕ(x)/∂x 0 at time t = x 0 = 0. For a spatially homogeneous system the most general initial Gaussian conditions can be characterized by The values for these initial correlation functions are taken to coincide with those from the corresponding quantum initial value problem. In the classical-statistical theory the correlation functions are obtained by performing a phase-space average over the initial configurations for a given observable with respect to the distribution function of initial fields W (ϕ(0), π(0)). The latter is chosen to give the prescribed correlation functions (4-7). Accordingly, we sample the initial field configurations until convergence to the prescribed initial correlations is achieved. Each initial configuration is separately evolved in time and the time evolution of correlation functions is obtained from the ensemble average. Using the corresponding definitions also for t > 0, we may extract the time evolution of occupation numbers f (t, p) in spatial Fourier space with p = |p| from [38] f (t, p) 1 The antisymmetric combination ∼ ϕ(x)π(y)−π(x)ϕ(y) is fixed by the corresponding field commutation relation in the quantum theory (or Poisson bracket in the classical theory) [38].
as well as the dispersion x) denotes the Fourier transform. With these definitions, the above initial conditions specify the initial field amplitude φ(t = 0) and distribution function f (t = 0, p). The initial frequency reads ω(t = 0, p) = p 2 + m 2 eff , where the effective mass term is given by In this setup the initial correlator (7) vanishes and we will also employφ(t = 0) = 0 in the following. The counter-term δm 2 Λ in Eq. (11) can be chosen to cancel the leading quadratic Λ-dependence of the threedimensional momentum integral over the initial F (t = 0, p). The same counter term is then used to cancel the associated divergence in the equation of motion (3). If not stated otherwise, we will perform such renormalized simulations. In addition, we will also present results for δm 2 Λ = 0 as employed in Ref. [29].
In the classical equation of motion (3) all model dependence on the coupling constant λ could be scaled out by the reparametrization ϕ → ϕ/ √ λ. The corresponding information is then entirely encoded in the initial conditions. The coupling drops out everywhere except for the 'quantum half' in the initial conditions, which would become 'λ/2'. Therefore, the coupling 'controls' the size of quantum corrections in the initial conditions [43]. This also reflects the fact that the classical-statistical theory will accurately describe the quantum physics as long as the coupling is sufficiently small. We do not rescale the field in the following for the purpose of our presentation.
The conserved total energy density is used to calculate the characteristic scale For the initial conditions, which we consider in the following, this scale is independent of the coupling constant. The first set of initial conditions is characterized by a large macroscopic field with σ 0 of order one. For this initial condition the mode occupancies are zero such that all modes are initialized with the vacuum 'quantum half' up to the lattice ultraviolet cutoff as in Ref. [29]. The second class of initial conditions is taken to be fluctuation dominated for comparison, with an overoccupied distribution up to some initial momentum Q 0 , (Fluctuation IC) and the initial amplitude of the particle distribution n 0 with Q ∼ 4 √ n 0 Q 0 . For both types of initial conditions one finds that the energy density scales with the coupling constant as ∼ 1/λ.

III. TURBULENT THERMALIZATION
The condensate driven initial conditions (13) lead to the well known phenomenon of parametric resonance. In the context of inflationary cosmology, the corresponding preheating dynamics [30] and the process of turbulent thermalization has been studied in great detail both using classical-statistical simulations [2,4,5,8,9,43,44] as well as directly in quantum field theory using resummation techniques based on the two-particle irreducible (2PI) effective action [2,31,32]. In the cosmological context, the preheating dynamics has also been investigated in the presence of fermion matter [45,46] and their consequences deduced for fermion production and the phenomenon of turbulent thermalization [41,42].
Initially, the energy density is stored in the large coherent field φ(t) for small coupling λ. Since φ(t) is rapidly oscillating, we will also consider the behavior of the envelope of the maximum field amplitude φ 0 (t) as a function of time. For the subsequent evolution one can identify three characteristic ranges in time, which are parametrically given as follows.
1. Instability for 0 t Q −1 log(λ −1 ): In this parametric resonance regime, the field φ(t) is a periodic oscillating function with frequency characterized by the initial rescaled field amplitude σ 0 [47,48]. This initial amplitude also determines the initial resonance band in momenta for which an exponential growth of f (t, p) can be observed [48]. The rapid growth of fluctuations leads to strong nonlinearities which broaden the resonance band and produce enhanced growth rates [31]. As a consequence, a wide range of growing modes lead to a prethermalization of the equation of state at the end of this early stage [1,49,50] while the distribution function itself is still far from equilibrium.
: In this regime, the field amplitude approaches a power-law behavior and the distribution function becomes selfsimilar [5]. Elastic scattering processes dominate and a dual cascade forms in distinct momentum ranges: A direct energy cascade towards the ultraviolet modes [5] and an inverse particle cascade towards the infrared develop [2]. These cascades are separated in momentum space by the characteristic momentum of the dominant resonance peak in the distribution function.
The direct energy cascade drives the thermalization process by energy transport towards higher momenta. The inverse particle cascade leads to the phenomenon of Bose condensation in this far from equilibrium regime [9].
The rescaled field amplitude √ λ φ0(t) for λ = 10 −4 . We employ two different cutoffs in order to demonstrate the insensitivity of the results on the lattice regularization. Also shown is the predicted power law ∼ t −1/3 (gray dashed line), where one observes very good agreement with data at later times. In the inset, the same data is shown to better resolve early times.
3. Thermalization for t Q −1 λ −5/4 : In this regime elastic and inelastic processes lead to a Bose-Einstein distribution. This regime is beyond the range of validity of classical-statistical simulations. At sufficiently late times the classical evolution will always end up showing a classical thermal distribution with a temperature parameter T Λ depending on the Rayleigh-Jeans cutoff Λ. The inability of the classical approach to describe this late-time thermalization stage has been studied in detail in the literature [33,38].
The characteristic time scales are given as weak coupling parametric estimates for our purposes. In addition, some of the prefactors are also available that can be significantly different from one for the case of the single component scalar field theory. For instance, the time scale for the end of the instability regime is more accurately given by t (2γ 0 ) −1 log(λ −1 ), where γ 0 0.033Q [48] is the largest growth rate of the primary instabilities.
The following numerical results are obtained using standard lattice discretization techniques [38]. In Fig. 1 we show the time-dependent rescaled field amplitude √ λ φ 0 (t) for a weak coupling λ = 10 −4 and σ 0 = 2.21Q. Two different data sets are shown corresponding to two different lattice momentum cutoffs Λ/Q = 10.2 and Λ/Q = 15.3. 2 The comparison confirms that the weak coupling results are insensitive to the Rayleigh-Jeans cutoff in classical-statistical field theory, which is a crucial ingredient for its ability to describe the corresponding quantum field dynamics [33].
The oscillating behavior with constant maximum amplitude at early times is visible from the inset of Fig. 1. Subsequently, for t ∼ Q −1 log(λ −1 ) the corrections from classical-statistical fluctuations change this behavior dramatically and trigger a transient rapid field decay. This happens when the size of fluctuations has grown such that their contribution to the energy density becomes comparable to that from the macroscopic field. At this stage the evolution becomes strongly non-linear which even leads to a temporary field growth [31]. The nonlinear dynamics finally leads to a power-law decay of the field amplitude with the predicted exponent δ = 1/3 [5] given by the gray dashed line in Fig. 1.
As shown in Ref. [5], during the turbulent stage between Q −1 log(λ −1 ) t Q −1 λ −5/4 the distribution function approaches a self-similar behavior, with the universal scaling exponents α, β and the stationary scaling function f S . The dynamical scaling exponents describe the evolution of typical occupation numbers and momenta. For hard modes, which dominate the energy density, β = −1/5 determines the evolution of characteristic momenta and α = −4/5 of their occupancy. The latter exponent determines the parametric time t quant ∼ Q −1 λ −5/4 at which the typical occupancies become order one and the classicality condition (1) is no longer fulfilled. Wave turbulence [51] describes the transport of conserved quantities such as energy density. Accordingly, in momentum space power-law cascades form. This can be either a 'direct cascade', for transport towards higher momenta, or an 'inverse cascade' into the infrared. In the inertial range of momenta, where the distribution function is described as a momentum power-law, one can write with a universal scaling exponent κ. The proportionality factor depends in general on time for isolated systems as in our case, i.e. without applied sources or sinks that could lead to stationary cascades. If there is more than one conserved quantity, the system may accommodate this by forming distinct cascades in different momentum regimes. It has been shown that for characteristic momenta p Q the dynamics is dominated by elastic scatterings [2]. This leads to an additional conserved quantity that emerges during the nonequilibrium time evolution even though total particle number is not conserved in the relativistic scalar field theory. As a consequence, apart from a direct energy cascade towards the ultraviolet, an inverse particle cascade towards the infrared can be observed. For the single component theory, this analysis was previously only done for the direct cascade [4,5]; the presence of the inverse cascade was not established.
In Fig. 2 we show our results for the evolution of the rescaled distribution function λf (t, p) for different fixed times. At times Qt = 160 and 240 one observes the resonance behavior of the instability stage. The resonance peak occurs first around p * /Q ≈ 1.1 and its location is proportional to the amplitude of the condensate at early times, p * = 0.52σ 0 [48]. Since the macroscopic field amplitude decreases with time, the peak is shifted to softer momenta [5].
At later times, the peak serves as a source for energy and particles, leading to momentum-scale invariant energy and particle fluxes towards the ultraviolet and the infrared, respectively. Accordingly, different momentum power-laws emerge to the left and to the right of the still visible dominant resonance peak in Fig. 2. On the double logarithmic plot the power-laws are well described by straight lines with different slopes, corresponding to different values for the exponent of (17) in distinct momentum ranges. The spectrum at Qt = 7100 is compared to three power laws with the exponents The hard scale exponent κ H = 3/2 describes the energy cascade towards higher momenta and results from effective 2 ↔ (1 + soft) processes involving the condensate [4,5]. The inverse particle cascade towards lower momentum modes shows two distinct momentum regimes depending on the size of the occupation numbers f (p) in each regime. For λf (p) O(1) the weak wave turbulence exponent κ M = 4/3 describes the transport of particles. In the nonperturbative regime of ultrasoft mo-menta where λf (p) O(1) the strong turbulence exponent κ S governs the dynamics [2].
We emphasize that the entire inverse particle cascade, which is described by κ M and κ S , can be understood in terms of elastic processes only. The enhanced scaling exponent κ S > κ M is a consequence of an emergent effective scattering matrix element for 2 ↔ 2 processes, which can be described in terms of a momentum dependent effective coupling λ eff (p) ∼ p 8 [2,8,9]. This analytical prediction is based on a systematic large-N expansion to next to leading order of the 2PI effective action for the N -component scalar quantum field theory, which gives κ S = d + 1 in d space dimensions [52]. Since we have N = 1, and in view of some reported deviations for N = 2 [35], it is remarkable that we find the predicted exponent with the observed accuracy.
In order to display the whole inverse particle cascade with both the weak and the strong turbulence regimes, we used rather large lattices up to 768 3 . Moreover, to cover the entire range of momenta shown in Fig. 2, we actually combined the data from two simulations with different lattice cutoffs Λ/Q = 10.2 and Λ/Q = 15.3. In the overlapping momentum regions the two simulations agree to very good accuracy. We emphasize that this is only done for presentational purposes of the physical weak coupling results. In particular, all results presented in the upcoming sections come from simulations including the entire range of displayed momenta.

IV. BASIN OF ATTRACTION FOR THE NONTHERMAL FIXED POINT
A crucial property of the turbulent regime is its strict independence of model parameters such as the value of the coupling constant λ or the initial conditions such as the initial field value σ 0 . This is a very powerful consequence of universality, which finds its manifestation in the self-similar behavior (16). The latter represents an enormous reduction of the possible dependence of the dynamics on variations in time and momenta, since it states that (Qt) −α f (t, p) only depends on the product (Qt) β p instead of separately depending on time and momenta. In general, renormalization group theory tells us that the fixed point distribution f S appearing in (16) will depend on all 'relevant' parameters of the system. For instance, if the initial field value σ 0 would represent such a relevant parameter then f S would in addition depend on the product (Qt) ζ σ 0 with some new exponent ζ. Plotting (Qt) −α f (t, p) only as a function of (Qt) β p/Q would then fail to describe the data. Therefore (16) represents a very strong statement about the loss of information about the parameters of the underlying system already at this transient stage of the nonequilibrium time evolution.
Since self-similarity has been extensively discussed for our theory already in the literature [5], we only address here those aspects of universality that are relevant for clarifying the conflicting statements mentioned in the in- The inset shows the corresponding occupation number distributions λf (t, p) at Qt = 6100. For the two fluctuation dominated initial conditions the inset shows in addition the initial spectra (gray). One observes that the evolution becomes independent of the initial conditions for the condensate driven as well as for the fluctuation dominated initializations.
troduction. In particular, we want to confirm that both the class of condensate driven initial conditions given by (13) as well as fluctuation dominated initial conditions (14) belong to the basin of attraction for this nonthermal fixed point. So far we have used the initial conditions characterized by a large field and small (vacuum) fluctuations. The second class of initial conditions we now consider are described in terms of highly occupied modes with initial amplitude n 0 /λ up to the momentum Q 0 . We consider n 0 = 39 as well as n 0 = 7.5. 3 The characteristic scale Q given by (12) is taken to be always the same, thus fixing Q 0 for given n 0 . We employ a weak coupling λ = 10 −4 .
Using both types of initial conditions, the evolution of the rescaled field amplitude is shown in Fig. 3. One observes that for each of the fluctuation dominated initial conditions a condensate builds up 4 , which is a consequence of the inverse particle cascade as described above [9]. Vice versa, the condensate driven initial conditions lead to fluctuations by parametric resonance. At late times the curves for all of the different initial conditions fall on top of each other to very good accuracy. This is further illustrated by showing the corresponding dis-tribution functions in the inset of the figure. While the single particle spectra are very different at initial time (Qt = 0), at late times the curves for all three distributions fall on top of each other. The fact that the evolution becomes independent of the details of the initial conditions already at this transient stage is a striking example of the phenomenon of universality far from equilibrium.

V. LIMITATIONS OF CLASSICAL-STATISTICAL SIMULATIONS
We now come back to the condensate driven initial conditions, as employed in Ref. [29], and explore the range of validity of classical-statistical simulations that provide a reliable description of the underlying quantum dynamics. Classical-statistical descriptions are restricted to the weak coupling regime as explained above. Therefore, it is important to verify up to which size of the coupling these methods can be applied. This is a time-dependent question since for sufficiently late times all classical-statistical descriptions break down once typical occupancies become order one -as is the case for thermal equilibrium 5 -and genuine quantum processes dominate. 6 For weak couplings, the time at which the instability regime ends depends logarithmically on the inverse coupling constant as explained above. The subsequent evolution becomes universal and thus independent of the coupling. The system exits the turbulent regime and enters the quantum one around the parametric time t quant ∼ Q −1 λ −5/4 . Increasing the coupling therefore means that the range of validity in time shrinks. Equivalently, for some fixed time one should observe deviations from universal behavior as the coupling is increased. The onset of any coupling dependence also signals the breakdown of the classical-statistical method. This is illustrated in Fig. 4 where we show the rescaled distribution function λf for different values of the coupling constant λ at fixed time Qt = 1200. One observes that for weak couplings λ 0.1 the spectra lie on top of each other to very good accuracy. 7 In contrast, for λ = 0.5 deviations are clearly visible and also the shape becomes qualitatively different, not showing the infrared enhancement or power law behavior. The presented results for λ ≤ 0.1 are insensitive to the employed Rayleigh-Jeans cutoff. The λ = 0.5 curve is obtained for Λ/Q = 15.3.
It is instructive to analyze the dynamics in this weakcoupling range in more detail. The classicality condition (1) implies that modes with low occupancies of f (p) 1 should not play a dominant role for the dynamics. Their contribution to the energy density may be estimated by This 'quantum' part should be compared to the total energy density contained in quasi-particle fluctuations part (t) = p ω(t, p)f (t, p) and in the macroscopic field, which we approximate by φ (t) = λφ 4 0 (t)/24. In the range of applicability of classical-statistical descriptions one expects that the ratio quant /( φ + part ) 1. Of course, such a simple separation into fluctuation and field parts The occupation number distribution f (t, p) for λ = 1 at different times for the same parameters as in Ref. [29]. Our results (lines) agree to very good accuracy with the previous study (points). The black dashed curve gives the corresponding result with mass renormalization to show that no significant differences can be observed at Qt = 6800. The orange dashed curve is a fit to the spectrum as employed in Ref. [29] to extract temperature T , mass m and chemical potential µ parameters.
is only applicable if there is a good quasi-particle description for weak couplings and such an analysis should always be taken with great care. For instance, we will see below that for λ = 1 very strong cutoff dependence occurs and the occupation number distribution becomes negative such that this analysis is inapplicable. However, we find that for the range of couplings λ 0.2 no strong cutoff dependence occurs and the distribution function fulfills f (p) ≥ 0 in the classical-statistical regime. We have seen above that quantum corrections become more important at later times. Therefore, one might hope that for short enough times one could extend the range of validity of classical-statistical simulations to stronger couplings. However, since the coupling controls the relative size of the occupation number per mode as compared to the 'quantum half' and the latter are cut off by the Rayleigh-Jeans regulator Λ, there is the danger that the sizeable cutoff dependence leads to a breakdown of classical-statistical simulations already at relatively early times in this case.
In Ref. [29] it is argued that at short enough times and strong enough couplings one does not enter the weak coupling attractor for turbulent thermalization and an alternative thermalization scenario can be observed within the framework of classical-statistical simulations. The authors find that for λ = 1 the classical-statistical thermalization dynamics is very different than what has been previously derived for weaker couplings. In particular, it is claimed that no turbulent cascades form and Bose condensation occurs as a consequence of the formation of a transient chemical potential.
In the following we reconsider the calculations of Ref. [29]. First we also use their employed 20 3 lattices and cutoff Λ/Q = 5.1 for the condensate driven initial conditions (13) in the absence of mass renormalization (δm 2 Λ = 0). In Fig. 6 our results are plotted by lines along with the data from the referenced study given by points of the same colors. Both results agree to very good accuracy. We also give the corresponding curve with mass renormalization (black dashed curve) to show that no significant differences can be observed. The orange dashed curve is a fit to the spectrum as employed in Ref. [29]. One finds that f (p) + 1/2 can be described by a classical Rayleigh-Jeans distribution T /(ω(p) − µ) with temperature parameter T , frequency ω(p) = p 2 + m 2 , (squared) mass parameter m 2 = λ 2 ϕ 2 +δm 2 Λ and chemical potential µ.
We then extend these studies to larger lattices up to the size 512 3 in order to also resolve the infrared physics. In particular, this will enable us to vary the Rayleigh-Jeans cutoff while being insensitive to finite volume effects. Our results with the same lattice cutoff but larger volume and mass renormalization are summarized in Fig. 7. The spectra start with the occupied quantum 1/2. The blue curves denote earlier times when the distribution function grows with time. For Qt = 130, one finds the resonance structure while at Qt = 6800 hard momenta can already be fitted by a classical thermal distribution, as indicated in Fig. 6. The entire spectrum can be described by a thermal function starting from the time Qt = 27000. Going to even later times, we find that the classical thermal fitting function changes its parameters with time since the effective mass decreases. Therefore, the temperature and the chemical potential of the fitting function become time dependent, where T grows while µ decreases. The difference between µ and m grows and the chemical potential decreases faster than the effective mass. With classical equilibration in the total occupation number, the ultraviolet tail of the distribution function f (p) becomes negative. Stated differently, the vacuum 'quantum half' decays [34]. This failure of the classicalstatistical approximation to describe the otherwise stable vacuum of the quantum theory is illustrated in Fig. 8, where we zoom into the hard momentum region. While for λ = 1 this unphysical decay of the 'quantum half' can already be observed at rather early times Qt ∼ 100, we emphasize that this does not happen in the weak coupling regime and clearly indicates the breakdown of the classical-statistical approach for large couplings.
Since the range of vacuum modes included in this classical-statistical setup is controlled by the Rayleigh-Jeans cutoff, one expects the results to become strongly cutoff dependent once the unphysical decay of quantum modes sets in. We will now compare the dynamics for various values of the cutoff Λ, with the mass renormalization as outlined above. Fig. 9 shows the field amplitude as a function of time. One finds that the condensate becomes zero at a finite time while the decay time itself can be varied to practically any value in a vast range by changing Λ. Since the decay happens faster the larger the cutoff, the prediction for Bose condensation of Ref. [29] has to be considered as an artifact of the employed regularization.
In Fig. 10 we show the spectra f (t, p) + 1/2 for four different lattice cutoffs. Here a similar picture emerges. The spectra at four different times ranging from Qt = 10 to Qt = 5500 exhibit cutoff dependencies. The table with the extracted fit parameters at Qt = 5500 reveals that no sensible prediction independent of the cutoff can be made.

VII. CONCLUSION
Classical-statistical simulations provide an important tool for ab initio descriptions of nonequilibrium quantum dynamics. However, they have a well defined range of validity which is restricted to weak enough couplings. Since the coupling controls the relative size of the occupation number per mode as compared to the 'quantum half' and the latter are cut off by the Rayleigh-Jeans regulator, the sizeable cutoff dependence indicates a breakdown of the classical-statistical approach for larger couplings. For the setup considered in the single component scalar field theory, we find that quantitative results can be obtained for λ 0.2 for parametric times t Q −1 λ −5/4 . In particular, our analysis shows that the results of Refs. [29,50] for λ = 1 are based on the application of the classicalstatistical approximation beyond its range of validity.
We have demonstrated that both the condensate driven as well as fluctuation dominated initial conditions belong to the basin of attraction of the same nonthermal fixed point within the range of validity of the classicalstatistical description. This leads to the phenomenon of turbulent thermalization, where the evolution is char-acterized by universality. Our results show the existence of a dual cascade for the single component scalar field theory. While this has been intensively studied for O(N ) symmetric N -component field theories, with a detailed analytic understanding using large-N techniques, our analysis points out that these provide also a remarkably accurate description of the universal properties for N = 1.
It would be extremely valuable to have a similar analysis for the range of validity of the classical-statistical approach describing the nonequilibrium quantum dynamics of longitudinally expanding systems. The latter are crucial for our understanding of ultrarelativistic collision experiments of heavy nuclei in the laboratory. Simulations in the theoretically clean weak coupling limit demonstrate the existence of a nonthermal fixed point in the space-time evolution of non-Abelian plasmas [15,16]. Alternative thermalization scenarios with the gauge coupling exceeding a certain strength [26] can be analyzed along the lines of the present work.