Large-N kinetic theory for highly occupied systems

We consider an effective kinetic description for quantum many-body systems, which is not based on a weak-coupling or diluteness expansion. Instead, it employs an expansion in the number of field components N of the underlying scalar quantum field theory. Extending previous studies, we demonstrate that the large-N kinetic theory at next-to-leading order is able to describe important aspects of highly occupied systems, which are beyond standard perturbative kinetic approaches. We analyze the underlying quasiparticle dynamics by computing the effective scattering matrix elements analytically and solve numerically the large-N kinetic equation for a highly occupied system far from equilibrium. This allows us to compute the universal scaling form of the distribution function at an infrared nonthermal fixed point within a kinetic description and we compare to existing lattice field theory simulation results.


I. INTRODUCTION
A fully microscopic description of the real-time dynamics of quantum many-body systems in terms of quantum field theory can be very demanding. Often, effective theories with a well-defined range of validity at some (long) time and distance scales provide an efficient alternative description. A well-known example is kinetic theory, which describes the state of the system in terms of a classical phase-space distribution of particles, f (t, x, p), at time t with position x and momentum p [1].
Accordingly, the derivation of kinetic theory from the underlying quantum field theory involves a series of crucial assumptions [2][3][4][5][6]. An important condition is that the de Broglie wavelength ∼ 1/|p| of relevant (quasi-)particles must be small compared to the mean free path between collisions. Otherwise, a description in terms of classical particles with a well-defined position and momentum between collisions would not be valid. Likewise, quantum interference effects between successive scattering events should not spoil a description in terms of independent scatterings. The existence of quasiparticle modes with a well-defined dispersion relation ω(p) translates in the language of quantum field theory to sufficiently narrow peaks of the spectral function [7].
These conditions can be often met in the presence of a sufficiently weak coupling or small diluteness parameter controlling the strength of the scatterings. In particular, controlled perturbative kinetic descriptions exist for fermionic quantum field theories, and scalar field theories close to equilibrium where the relevant modes with momenta of the order of the temperature have occupancies of order one [2]. Likewise, perturbative descriptions exist far from equilibrium [8][9][10] if the occupancies of typical particle modes are not too high such that f 1/λ with λ representing the relevant coupling constant or diluteness parameter. Though gauge theories are more involved, perturbative kinetic descriptions dealing with the problem of quantum interference have been given [11][12][13].
Much less is known about effective kinetic descriptions for general far-from-equilibrium situations. Pressing applications concern systems where the occupancies of relevant modes are non-perturbatively large (f ∼ 1/λ) such that a perturbative power counting in terms of a small coupling parameter fails.
An important example concerns the early stages of a relativistic heavy-ion collision (for recent reviews see [14,15], and [11][12][13]16] for current perturbative kinetic descriptions). In this situation, typical gauge boson occupancies can become non-perturbatively large at low momenta below the Debye mass scale. These modes may influence the evolution of important quantities like the longitudinal pressure P L of the expanding plasma, evidence of which was found from real-time lattice simulations [17][18][19]. Recent studies have found universal scaling behavior of infrared modes [20] that may be connected to nontrivial field configurations [21]. The possible existence and influence of an enhanced low-momentum region for non-Abelian plasmas out of equilibrium has been extensively discussed in the literature [22][23][24][25][26][27][28][29], as well as methods for accessing spectral information at the Debye scale and below [30][31][32][33].
Remarkably, longitudinally expanding non-Abelian plasmas and self-interacting scalar field theories are found to share important universal aspects of their far-from-equilibrium evolution [17,34]. Similar selfsimilar scaling properties are known for a wide variety of highly occupied systems. These include relativistic scalar systems, often used in inflationary models of the early universe after a period of resonant particle production [9,[35][36][37], and non-relativistic systems such as ultracold quantum gases or other condensed matter systems after a strong quench [38][39][40][41]. Characteristic infrared properties of these highly occupied systems turn out to be quantitatively the same for both relativistic and non-relativistic models [6,42]. The universal scaling properties are associated to a nonthermal renormalization group fixed point [37,43] defining a universality class out of equilibrium, which encompasses nonrelativistic and relativistic N -component field theories [6,44], scalar systems in different geometries [17], in different spatial dimensions [45,46], and it can even be observed for attractive quartic interactions as long as arXiv:1710.11146v1 [hep-ph] 30 Oct 2017 mean interactions are repulsive [47]. Perturbative kinetic approaches [48,49] are not able to reproduce key features of this low-momentum dynamics.
In order to describe the evolution also for nonperturbatively large occupation numbers, we consider an effective kinetic description for scalar systems which is not based on a weak-coupling or diluteness expansion. Developed in [6,45], it exploits the fact that often one describes complex many-body problems with more than one particle species. In this case, alternative kinetic descriptions with an extended range of validity may be derived based on non-perturbative expansions in the number of species available. For scalar field theories with N species and quartic self-interactions, it results from a large-N expansion to next-to-leading order (NLO) based on a two-particle irreducible (2PI) resummation of self-energy diagrams [37,39,43,50,51], which translates to a vertex resummation in the kinetic framework [6,45].
So far, the large-N kinetic theory at NLO has been successfully applied to analytically compute the selfsimilarity exponents near nonthermal fixed points at low momenta, well agreeing with lattice results [6]. However, a complete characterization of the nonperturbative infrared regime involves also the scaling form of the distribution function, which has not been established from the large-N kinetic theory yet. In this work, we present the first numerical solution of the large-N kinetic equation applied to the universal low-momentum scaling regime in three spatial dimensions. Our results are found to compare rather well to available lattice simulation data of the underlying field theory, in particular, establishing a ∼ |p| −4 tail of the distribution in the regime dominated by number conservation. We analyze in detail the range of validity and quasiparticle picture of the large-N or "vertexresummed" kinetic theory, and show how it encompasses and extends standard perturbative descriptions.
The paper is organized as follows: In Sec. II we consider scalar N -component field theory. Starting from relativistic models with quartic self-interaction, we discuss the non-relativistic low-energy limit relevant, e.g., also for the description of ultracold Bose gases. Sec. III summarizes main aspects of perturbative kinetic theory, before we present the large-N kinetic description in Sec. IV. The latter has an extended range of validity based on the inclusion of vertex corrections, which is analyzed in detail in Secs. V and VI. We present a numerical solution of the large-N kinetic theory for the description of a nonthermal fixed point in Sec. VII. After concluding in Sec. VIII, we end with two appendices on calculational details of the collision integrals (App. A) and on integration boundaries (App. B).

II. RELATIVISTIC AND NON-RELATIVISTIC SCALAR FIELDS
We consider an O(N ) symmetric quantum field theory for the field components ϕ a (t, x), a = 1, . . . , N with time t and space variable x in three dimensions and quartic self-interactions. For N = 4 the relativistic model describes the Higgs sector of the Standard Model of particle physics [52]. In the context of low-energy descriptions of quantum chromodynamics, such a model encodes the three pions and the sigma resonance. Inflaton models for early universe cosmology often employ related multi-component field theories [53]. In a nonrelativistic setting, the Heisenberg magnet for N = 3 is a prominent example, and the case N = 2 can be used to describe the two real components of a complex Bose field in systems of ultracold atomic gases dominated by s-wave scattering [54].
The considered relativistic quantum theory is described, on a classical level, by the action (1) with the notation t,x ≡ dt d 3 x and the (renormalized) mass m and coupling parameter λ. Here summation over repeated Lorentz indices µ = 0, . . . , 3 and field indices a = 1, . . . , N is implied. We will always employ natural units, with the speed of light, Boltzmann's constant and the reduced Planck constant equal to unity, i.e. c = k B = = 1.
For processes with characteristic momenta below the mass scale m, one expects an effectively non-relativistic description to become relevant even for the relativistic microscopic model (1). More generally, the descriptions of ultracold quantum gases or other condensed matter systems typically employ non-relativistic field theories. To provide an explicit link between relativistic and nonrelativistic descriptions, one may have in mind the phenomenologically important case of the N = 2 component theory. It can be alternatively described by the complex field ϕ = (ϕ 1 + iϕ 2 )/ √ 2. Following standard procedures [55], the effective low-energy description for the non-relativistic limit of (1) may then be characterized by the action with non-relativistic field φ = ω/2 ϕ + i/ √ 2ω π * in terms of the original complex field ϕ(t, x) and its conjugate momentum field π(t, x) for dispersion ω. Here we also introduced the effective non-relativistic coupling which is no longer dimensionless. For dilute Bose systems, g can be related to the s-wave scattering length, given by a = mg/(4π) [55]. Similar relations can also be given for general N , which is employed in section VI. For the non-relativistic quantum theory the expectation value of the particle density n = φ * φ is conserved, and we will consider spatially homogeneous systems. The density n and scattering length a can be used to define a characteristic "coherence length" whose inverse is the momentum scale Q = √ 16π a n ∼ √ m g n .
We also define the "diluteness parameter" which provides a dimensionless expansion parameter for the non-relativistic system, similar to the dimensionless coupling λ for the relativistic system. Specifically, with (3) we obtain ζ ∼ (Q/m) λ.

III. PERTURBATIVE KINETIC THEORY
The derivation of perturbative kinetic equations from the underlying quantum-statistical field theory employs an expansion in terms of a small coupling λ 1 or small diluteness parameter ζ 1, together with a gradient expansion for not too early times [3,4]. The phase-space distribution function of particles, f (t, p), describing the state of the spatially homogeneous system at time t and momentum p, is obtained from the expectation value of two-field correlators evaluated at equal times.
More precisely, for the relativistic field theory the time derivative of the anti-commutator expectation value {ϕ a , ϕ b } ≡ ϕ a ϕ b + ϕ b ϕ a determines (in the absence of external forces) the change of the distribution function according to [45] Here the frequency ω and spatial momentum p arise from the Fourier transform with respect to the relative space-time arguments of the two fields, while the remaining time dependence describes the breaking of time-translation invariance of the spatially homogeneous system out of equilibrium. The kinetic description involves the projection onto positive frequency contributions by integrating over ω, respectively, and we exploit O(N ) symmetry assuming no spontaneous symmetry breaking such that {ϕ a , ϕ b } ∼ δ ab . The kinetic equation may then be computed perturbatively by taking into account interaction effects, which are subsumed into the "collision term" C[f ] to obtain In its range of validity, the leading contributions to C[f ] in a coupling expansion and an expansion to lowest order in gradients for the massive scalar field theory (1) leads to the well-known Boltzmann equation with the collision integral for elastic 2 ↔ 2 scatterings [55], with the notation q ≡ d 3 q/(2π) 3 and the relativistic dispersion The functional I 2↔2 [f ] contains the distribution functions f p ≡ f (t, p) describing the changes by loss or gain through scattering: In the last equation, we give the approximate expression for large occupancies that will be useful later.
Since the overall collision term C[f ] is of order λ 2 , further perturbative corrections to terms appearing in the integrand of (8) are subleading. In particular, it allows one to employ in the integrand a well-defined dispersion relation (9). Phrased in terms of the underlying field theory, this corresponds to employing a quasiparticle approximation for the spectral function given by the expectation value of the commutator of two fields [8,45]: Beyond lowest order in the coupling, the spectral function receives corrections leading to a non-zero width of the spectral function encoding "off-shell" contributions to processes. However, since they are of higher order in the perturbative power counting for the collision term, there is a well-defined quasiparticle description of the collision dynamics.
Similarly, assuming for simplicity a non-relativistic dispersion ω p = |p| 2 /2m, the collision term for the kinetic equation of the theory with action (2) becomes [10] C nr [f ](t, p) = l,q,r The assumption of a quadratic dispersion relation appears natural in this context since it corresponds to a low-momentum limit of the relativistic collision integral (8) [9,49]. However, we note that for non-relativistic theories one generally has a Bogoliubov dispersion relation with a quadratic dispersion at higher and a linear dispersion at lower momenta if a Bose-Einstein condensate exists [10]. The perturbative power counting for λ 1 leading to (8), or (12) for ζ 1, with elastic 2 ↔ 2 scatterings as in (10) assumes that the relevant occupancies f p for typical momenta are not too high. More precisely, only for f p 1/λ in the relativistic, or f p 1/ζ in the non-relativistic case, the higher-order corrections are parametrically small. Before we discuss this issue in more detail below in section VI, we will introduce in the following an alternative kinetic description based on a large-N expansion of the underlying quantum field theory.

IV. LARGE-N KINETIC THEORY
The standard Boltzmann equation described in the last section is based on a weak-coupling expansion, Scattering processes at next-to-leading order in the large-N expansion, which are mediated by an effective interaction. While the solid lines represent particles with given four-momenta, the dashed line represents the function (16) in Fourier space describing s-, t-and u-channel exchange.
which restricts the range of validity of the kinetic theory to perturbative problems. However, for the Ncomponent field theory an alternative kinetic description with an extended range of validity may be derived based on a non-perturbative expansion in N . For details about its derivation from the underlying quantum field theory we refer to Refs. [6,55]. Here we give the relevant expressions that are used below to solve the large-N kinetic equation. At large N the classical action (1) scales proportional to N employing ϕ a ϕ a ∼ N . Genuine quantum corrections due to scatterings appear at subleading orders in a large-N expansion, i.e., they are down by factors of 1/N compared to classical contributions [39,51]. In particular, the spectral function at leading order (LO) in a large-N expansion reads where the dispersion ω rel p = M 2 + p 2 for the relativistic theory contains an effective mass term M 2 that, additionally to the initial mass parameter m, involves contributions from fluctuations [6,51]. The mass term is constant at lowest order in the gradient expansion underlying kinetic descriptions [4].
We will describe in the following that at NLO in the large-N expansion, there is a well-defined effective kinetic description in terms of scatterings between quasiparticles. Similar to the previous section III, we start by considering the relativistic theory and extend the discussion to the non-relativistic case in the end.
In order to discuss subleading corrections in the 1/N expansion, it is convenient to employ the auxiliary field formulation of the same model [51]. For this purpose, we rewrite the original action (1) by introducing an auxiliary field χ(x) as Integrating out χ in the defining functional integral yields the original action, and from the Heisenberg equations of motion one sees that the auxiliary field represents the composite operator While the auxiliary field is not a dynamical degree of freedom, it can be used to conveniently express scatter-ing corrections in terms of the expectation value Since χ represents a two-point function according to (15), the function D(x, y) encodes a four-point function or vertex. Specifically, at NLO in the 1/N expansion, scatterings are mediated by (16) [51]. This is indicated in Fig. 1, where dashed lines represent the two-point function (16) in Fourier space. The modified vertex at NLO is shown for scatterings in the s-, t, and u-channel, respectively. Accordingly, the effective kinetic equation at NLO is given by the same kinetic equation (7), however, with the different collision term [6] C rel where ω rel p = M 2 + p 2 . In the derivation of the collision integral, the LO expression of the spectral function (13) is used since its subleading corrections in 1/N would result in subleading corrections of the collision integral, which are part of the large-N kinetic theory at next-to-next-to-leading order (NNLO), and are thus omitted. The time-and momentum-dependent effective coupling function incorporates the vertex corrections for the different scattering channels according to Fig. 1. The appearance of the renormalized "one-loop" retarded selfenergy in the denominator of Eq. (18) is the result of an analytical continuation of a geometric series summation of an infinite number of scattering processes at NLO in the large-N expansion [50]. We emphasize that Π rel R , and thus also λ 2 eff , is time-dependent since it depends on the evolving distribution function.
From (18) one observes that for |Π rel R | 1, which is the case for weak enough coupling (λ 1) and not too large typical occupancies (f 1/λ), the vertex corrections encoded in the momentum-dependent effective coupling become irrelevant such that λ 2 eff λ 2 . In this case the collision term (17) essentially describes standard perturbative two-to-two scatterings, however, at large N . 1 In contrast, for high characteristic occupancies with f ∼ 1/λ the collision term (17) can be strongly modified if Π rel R starts to become of order one. We discuss the corresponding behavior of the effective coupling in more detail below, for which we first discuss the non-relativistic effective kinetic equation relevant at low momenta.
Following along the lines of section (III), we consider first the case N = 2 to illustrate the effective kinetic equation for a non-relativistic complex scalar field, i.e., with two real field components. The case of general N then proceeds accordingly [56].
As above, we assume again a quadratic nonrelativistic dispersion relation ω p = |p| 2 /2m instead of the full Bogoliubov dispersion that would also contain a linear regime [10]. This is motivated by a recent study of unequal-time correlation functions, where the transport coefficient z ≈ 2 with ω p ∼ |p| z has been measured [57], and it also leads to technical simplifications. One obtains in this case [6] C nr The effective coupling in the collision integral reads with the "one-loop" retarded self-energy and the momentum difference P = p − q. As for the relativistic theory, the collision integral (20) reduces to its perturbative expression 2 for small |Π R | 1, while it encounters modifications otherwise. Moreover, when this non-relativistic theory is considered to be the low-momentum limit of a relativistic theory as in (14), then the mass m in non-relativistic kinetic expressions should be replaced by the effective mass M , implying ω p = |p| 2 /2M , ζ ∼ (Q/M ) λ, etc. Therefore, also relativistic theories with a vanishing mass parameter m = 0 have this non-relativistic limit for small momenta below the effective mass |p| < M . To keep the notation simple, we will stay with the symbol m for the mass while implying the replacement m → M whenever we compare to a relativistic theory.
For later use, it is helpful to further evaluate the expressions for C nr NLO and Π R . Using magnitudes of momenta p = |p|, and similar for q, k and P , we find for an isotropic system (details are given in App. A) and where we have dropped the labels of C nr NLO to shorten the notation. Accordingly, the effective kinetic equation (7) depends on time t and the magnitude of the momentum p.

V. BEHAVIOR OF THE LARGE-N RESUMMED EFFECTIVE VERTEX
To discuss the extended range of validity of large-N kinetic theory, we first consider the effective coupling g 2 eff appearing in the collision integral (23), which is a function of the difference in energies of the in-and outgoing particles ω = ω p − ω q = (p 2 − q 2 )/2m and of the magnitude of the momentum change P = |p − q| in a scattering event. It is beneficial to consider limiting cases of its arguments to analyze its behavior. We distinguish three typical collision scenarios. In the first case, the momentum of a particle within a collision is strongly decelerated so that p q and thus 2mω ≈ p 2 ≈ P 2 . Similarly, q p leads to the same effective coupling because g 2 eff [f ](t, ω, P ) is symmetric in ω. To ease the discussion we will therefore use ω ≥ 0 in the following. In the second case, the magnitude of the momentum of the particle undergoing the collision stays at the same order p ∼ q while it changes its direction. The nearly collinear regime is discussed separately and constitutes the third scenario.
The effective coupling in (21) can be calculated from the retarded self-energy Π R (t, ω, P ) given in (24). For the distribution function f (t, k) that enters the integrals in Π R , we assume that for l = 0, 1, 2 integrals of the form are dominated at the possibly time-dependent momentum scale K, if it lies within the integration limits. This scale K can be defined as the momentum where such that it provides the dominant contributions to the particle number density For the integrals in (25) to converge between the maximal limits of 0 and ∞, the distribution function f (t, k) should fall off faster than k −3 at large momenta k K and decrease slower than k −1 at low momenta k K, or not decrease there at all.
We start with the regime p q. Then one has P ≈ p and ω ≈ ω p = p 2 /2m ≈ P 2 /2m. The second relation states that the energy difference of in-and outgoing momenta ω follows the non-relativistic dispersion relation with momentum (difference) P , and we refer to this regime as dispersive. The "one-loop" retarded selfenergy (24) in this limit reads Its real part involves an integration over all momenta and we can use (25) in the limiting cases of P K and P K. In the first case, the integrand is dominated at low momenta k and we can approximate log(k + P ) − log |k − P | ≈ 2k/P + O((k/P ) 3 ). Similarly, the second case leads to log(k + P ) − log |k − P | ≈ 2P/k + O((P/k) 3 ). Hence, the limiting expressions are For the imaginary part, large and small ingoing momenta P K and P K lead to the expressions In (32) we used that kf (t, k) should be a growing function at low momenta to be consistent with (25).
To get the corresponding limiting expressions for the effective coupling, we first assume that for typical soft momenta K the occupation number f (t, K) is sufficiently large such that for the considered momenta P , the 1 in the denominator of g 2 eff in (21) can be neglected, and the effective coupling reads g 2 eff ≈ g 2 (Re Π R ) 2 + (Im Π R ) 2 −1 . Since P f (t, P ) is limited by K f (t, K), the real part (30) dominates at low momenta P K. On the other hand, the imaginary part (31) decreases slower than the real part at high momenta P K, and is thus larger. With this, the effective coupling is parametrically Accordingly, it is constant below K and follows the power law P 2 beyond K. We note that at even larger momenta, it becomes constant g 2 when the +1 in the denominator of its definition becomes important.
As noted above, the regime q p leads to the same expressions (33), (34), with P ≈ q. The frequency argument gets a minus sign −q 2 /2m, which does not change the values of g 2 eff because of its symmetry.
In the momentum-dominated regime, we consider the situation when in-and outgoing momenta are of the same order p ∼ q. For further simplifications, we also assume that the frequency difference is small ω = ω p − ω q P 2 /2m. This occurs for most of the scattering angles cos θ pq = pq/pq, since this assumption is equivalent to the condition cos θ pq min(q/p, p/q). The remaining case of nearly collinear collisions cos θ pq ∼ 1 will be discussed in Sec. V C.
In the considered regime, the imaginary and real parts of the "one-loop" retarded self-energy (24) become Expanding the logarithm of the real part and approximating the integral as in (25), one arrives at similar expressions as in the dispersive case A closer look on the logarithm of the original expression in (24) reveals that if 2mω 2KP is satisfied, the estimate for lower momenta (38) is even valid in the collinear regime where 2mω exceeds P 2 . Hence, the full range of validity for (38) is 2K P 2mω/2K, which may only hold if (2K) 2 2mω. Otherwise, for (2K) 2 2mω, no region with the value (38) exists.
To compute the effective coupling, we again assume large occupation numbers and thus neglect the 1 in the denominator in (21), which yields g 2 eff ≈ g 2 (Re Π R ) 2 + (Im Π R ) 2 −1 . For both small and large momenta P , the real part is larger than the imaginary part Re Π R Im Π R . This follows from 2mω P 2 and, for small momenta P/2 K, from (P/2) f (t, P/2) K f (t, K) while for larger momenta P/2 K it results from (P/2) 3 f (t, P/2) K 3 f (t, K), that are both requirements for f (t, k) and were formulated below Eq. (25). Hence, the effective coupling parametrically follows The main difference to the dispersive regime is the steep power law P 4 at large momenta. Interestingly, the transition between small and large momentum expressions proceeds at the slightly larger scale 2K.
The remaining case is when in-and outgoing momenta are nearly collinear cos θ pq ∼ 1, i.e. the case where P 2 2mω. The logarithm appearing in Re Π R in (24) can be written as where we have expanded it in the second line. Assuming that the integral is dominated at momenta k ∼ K as in the cases above, this expansion is justified for large momenta P 2K while the condition 2mω 2KP is additionally required at low momenta P 2K. With this, we can readily estimate the real and imaginary parts of Π R as where we have introduced the inverse momentum P inv = 2mω/P . Recall that for the real part, the remaining situation of 2mω 2KP for low momenta P 2K in the collinear regime has been discussed in Sec. V B, where it led to the expression (38).
From this, we can compute the effective coupling in the collinear regime. Interestingly, the real part is negative and 1 + Re Π R in the denominator of the effective coupling (21) may become zero, leading to a resonant increase of the coupling. Otherwise, if f (t, K) is sufficiently large, the 1 in the denominator of (21) can again be neglected and the real and imaginary parts of Π R can be compared in order to estimate At first, we consider P inv 2K. This condition translates to 2mω Eff. coupling: g 2 ω, P ) without 1 in its denominator as a function of momentum for the momentumdominated regime (blue line, ω = 0) and for the dispersive regime (green line, ω = P 2 /2m). The distribution function is chosen as in Eq. (45) with parameters A = 84100/(4m 2 ) and B = 0.0521×(2m). These values are taken from Ref. [6], where this form was fitted to the distribution function in a classical lattice simulation at time t = 300/(2m). Keeping κ< = 0 fixed, we show bands for 3.5 ≤ κ> ≤ 6.
2KP and corresponds to the real part as given by (42).
Similarly, the case P inv 2K translates to 2mω 2KP and the real part is then given by (38).
for low momenta, one again finds Im Π R |Re Π R |. Thus, as in the momentum-dominated regime, the real part dominates the effective coupling. At low momenta P 2K it leads to (40) for 2mω 2KP while in all other cases in the collinear regime, one has Hence, different from the momentum-dominated regime, the effective coupling decreases here as P −4 .

D. Comparison to numerical results
We now compute the effective coupling g 2 eff numerically as defined in (21). For the distribution function we use which has been suggested in Ref. [6] to approximate the distribution function at low momenta during the self-similar regime. The parameter B is related to the momentum scale K defined in Eq. (26) via K = B ((2− κ < )/(κ > −2)) 1/(κ>−κ<) , such that both are of the same order K ∼ B and, for κ < = 0 and κ > = 4, even equal K = B. The function in Eq. (45) reduces to power laws ∼ p −κ> for large p B and to ∼ p −κ< for small p B momenta. The effective coupling is shown in Fig. 2 for the momentum-dominated regime (labeled ω = 0) and for the dispersive regime (labeled ω = P 2 /2m) with κ < = 0. In order to show that g 2 eff follows the parametric estimates derived in Secs. V A and V B and is thus insensitive to the precise functional form of f (t, p), we vary 3.5 ≤ κ > ≤ 6 and pool the resulting curves into a band for each regime. Moreover, the 1 in the denominator of the effective coupling (21) has been removed to be able to follow the power laws to larger momenta. 3 One observes that while at low momenta g 2 eff stays constant, it grows as P 4 in the momentumdominated regime and as P 2 in the dispersive regime, irrespective of the detailed form of f (t, p). In addition, the transition between constant and power law behavior occurs roughly at 2B in the prior and at B in the other case. We checked that using κ < = 0.5 leads to similar results. These observations confirm the parametric expressions for the effective coupling in Eqs. (33,34,39,40). Combining these scenarios with the collinear regime, we can also understand the functional form of the effective coupling for any fixed 2mω = p 2 − q 2 as a function of P . Before showing our numerical results, we use our parametric estimates to illustrate the functional form of g 2 eff and its relevant momentum scales in Fig. 3 for sufficiently low frequency differences 2mω (2K) 2 , where this time we also include the 1 in the denominator of (21). For large momenta P 2K, the coupling g 2 eff is in the momentum-dominated regime and follows (39), growing as a power law P 4 . In contrast, at low momenta P mω/K, it is in the collinear regime and decreases as a power law P −4 due to (44). Between these regions mω/K P 2K, it takes its minimal value g 2 eff ∼ (m K f (t, K)) −2 as in Eq. (40). Additionally, it becomes g 2 eff ≈ g 2 at very large and very low momenta, parametrically given by respectively, where we used (39) and (44). The inverse coherence scale Q defined in Eq. (4) enters the expressions because of Therefore, the momentum-dominated regime with ω = 0 shown in Fig. 2 is a special case of this functional form. According to our considerations in Secs. V B and V C, there is no constant region with g 2 eff ∼ (m K f (t, K)) −2 for larger frequencies 2mω (2K) 2 , but the P −4 power law changes over to P 4 at the transition scale P 2 ≈ 2mω between collinear and momentumdominated regimes.
All of this is also seen in our numerical results of g 2 eff for different 2mω in Fig. 4. For low frequency differences 2mω (2K) 2 each of the parts visualized in Fig. 3 can be identified (blue, green and red curves) while for larger frequency differences (2K) 2 2mω (2Q) 2 , the constant middle part is absent and the transition between the power laws occurs at P 2 ≈ 2mω (cyan and violet curves), as expected from parametric estimates. At large frequency differences 2mω (2Q) 2 , the conditions in (46) are true for almost all momenta P and thus, one has g 2 eff ≈ g 2 such that no clear power laws are visible (yellow and black curves). 4 Besides these general features, we can also understand the origin of the spikes that are visible in Fig. 4. The upper spikes result from the negative Re Π R in the collinear regime approaching −1 and thus canceling the 4 The frequency of the orange curve is of the order of 2mω ∼ (2Q) 2 , which is between the described functional forms.
1 in the denominator of g 2 eff . This occurs roughly at the scale when g 2 eff ≈ g 2 becomes constant at low momenta. Because of Im Π R |Re Π R | ≈ 1, the effective coupling becomes g 2 eff ≈ g 2 /(Im Π R ) 2 g 2 at its maximum and exceeds there g 2 , which results in a peak. The lower spikes are a consequence of the real part changing its sign between the collinear and momentum-dominated regimes for P 2K. Then the real part vanishes Re Π R ≈ 0 and one again has g 2 eff ≈ g 2 /(Im Π R ) 2 . Since this occurs at P 2 ≈ 2mω, we can use the expression for the imaginary part in the dispersive regime (31), such that the effective coupling grows as P 2 in that case. For comparison, a corresponding power law curve is shown in Fig. 4, which is in good accordance with the peaks of the lower spikes.
This constitutes a detailed understanding of the functional form of g 2 eff based on parametric estimates. It also shows the importance of the inverse coherence scale Q. All nontrivial values of the effective coupling g 2 eff < g 2 occur for momenta below this scale, such that Q can be regarded as the separation scale between the non-perturbative infrared region and the perturbative regime at hard momenta.
In the language of relativistic scalar theory for large infrared occupation numbers, the inverse coherence length is related to a mass shift contained in the effective mass M . 5 Hence, one has M Q, which implies that the whole low-momentum region where the effective coupling takes nontrivial values is below the mass scale and thus, essentially non-relativistic. This has indeed been observed in lattice simulations [6,17]. Therefore, our analysis of the effective coupling should also be valid for the relativistic case.

VI. EXTENDED RANGE OF VALIDITY OF LARGE-N KINETIC THEORY
Having understood the behavior of the effective coupling g 2 eff in the previous section, we can study its consequences for the large-N kinetic theory. Different from the perturbative kinetic approach, we show here that no limitations on the typical occupancy of the distribution function exist and one can therefore consider also highly occupied systems. In particular, we demonstrate that the large-N kinetic theory coincides with the standard perturbative description at large N for sufficiently low typical occupancies or if only large momenta p Q are considered.

A. Role of the effective vertex in the collision integral
Using the results of Sec. V, we first discuss which values of g 2 eff actually occur in an evaluation of the collision 5 The effective mass M 2 includes a term that results from local self-interactions.
integral (23) for incoming momentum p. In particular, this will allow us to identify standard perturbative descriptions in the regime where the occupancies are not large.
The integration variables relevant for the effective coupling are q and P . While the q-integration proceeds over all momenta, the P integration is limited by P < = |p − q| and P > = p + q and thus, not all values of g 2 eff that were discussed above occur for given p. The relation provides a useful link between the frequency difference 2mω = p 2 − q 2 and the integration limits.
In the dispersive regime, one has P 2 ≈ 2m|ω|, with P ≈ p or P ≈ q, and the resulting effective coupling is written in (33), (34). The P -integration can then be simplified to While for q p the effective coupling is approximately g 2 eff [f ] t, p 2 /2m, p , which stays constant during the qintegration, the coupling g 2 eff [f ] t, q 2 /2m, q varies for q p. When the momenta q ∼ p are of the same order, the P -integration range becomes nontrivial. However, because of (49) there is always a region with P 2 ≈ 2m|ω| as in the dispersive regime for any p and q. The integration limits P 2 > > 2m|ω| and P 2 < < 2m|ω| are thus located around this region and correspond to the momentum-dominated and the collinear regimes, respectively. There the value of the effective coupling is given by (39), (40) and (44). We note that for the special case that both momenta p, q K, the effective coupling takes its minimal value g 2 eff ∼ (m K f (t, K)) −2 for all momenta P in the integration.
The effective coupling becomes trivial g 2 eff ≈ g 2 at sufficiently large frequency differences 2m|ω| (2Q) 2 . Because of 2mω = p 2 − q 2 , this statement depends on the integration variable q in the collision integral and may occur as part of the integration. Similarly, this happens at low and high momentum differences P as in (46).
In this context, an interesting case is when the incoming momentum exceeds the inverse coherence length p Q. Then the effective coupling stays mostly trivial, i.e., g 2 eff ≈ g 2 during the q-integration. For p q and p q this results from P Q. The qmomentum range where g 2 eff takes nontrivial values is located around p ∼ q and, because of P ∼ Q q, it is small when compared to typical q values and to the whole q-integration range that proceeds from 0 to ∞. This shows that the large-N kinetic theory reduces to the perturbative kinetic theory if dynamics at only large momenta p Q is studied. Furthermore, g 2 eff also becomes trivial when K Q, which is independent of integration parameters. This condition is equivalent to m g K f (t, K) 1, which, as − + − . . . will be shown in Eq. (53), corresponds to ζ f (t, K) 1, where ζ is the diluteness parameter. Hence, when typical occupation numbers are sufficiently small, the effective coupling becomes g 2 for all values of the integration parameters and the large-N kinetic theory reduces to the perturbative kinetic theory.

B. Interference terms and their large-N resummation
The collision integrals of perturbative and large-N kinetic theories differ in the coupling, which is constant g 2 or resummed g 2 eff , respectively. To understand the range of validity of the theories due to interference, we consider a series of diagrams that all contribute to the elastic 2 ↔ 2 channel, as depicted in Fig. 5. Since the retarded self-energy Π R as defined in Eq. (24) is originally a convolution of the distribution function and the retarded propagator, it corresponds to a "blob" connecting two vertices. With the results of Sec. V, it can be parametrically estimated as where we have included the diluteness parameter ζ ∼ m g Q with the scale Q 2 ∼ m g K 3 f (t, K), which are defined in Eqs. (5) and (4). As long as |Π R | 1, the series of chain diagrams in Fig. 5 corresponds to a (perturbative) loop expansion whose leading term in the collision integral is given by the coupling constant g 2 . Hence, is a necessary condition for the validity of the perturbative kinetic theory. 6 The equivalent condition on the right hand side restricts the occupation numbers for that case. Otherwise, when |Π R | 1, all of the considered scattering processes are of the same order and interfere, which leads to the break-down of the perturbative approach.
This problem is absent in the large-N kinetic description because of the effective coupling g 2 eff in (21). Its perturbative expansion corresponds to the series of 6 Although we consider a non-relativistic theory with a quadratic dispersion relation, it is enlightening to check how this condition translates into the language of a relativistic scalar theory. Assuming that the effective mass M is of the order of Q, one can identify λ ∼ ζ. Hence, perturbative kinetic theory is valid for small enough occupancies up to the inverse coupling f (t, K) 1/λ, which is consistent with what is commonly known (see, e.g., [3,55].
scattering processes that was depicted in Fig. 5. Hence, the effective coupling appears as a geometric series at small Π R where the series converges. More precisely, it resums all diagrams emerging at NLO in a 1/N expansion [37,51].
Hence, the interference of an infinite number of scattering processes is encoded into the vertex resummation in the large-N kinetic theory and results in an effective coupling g 2 eff . In this effective kinetic theory, characteristic occupancies are not restricted. Thus, it comprises and extends the range of validity of standard perturbative kinetic descriptions.

C. Cross section and mean free path
In Sec. IV, we have argued that the large-N kinetic theory appears as a consistent quasiparticle description in the large-N limit since the spectral function can be approximated by an on-shell form in the collision integral. Moreover, this effective theory extends perturbative kinetic approaches and can also be applied to highly occupied systems, as we argued in this section. In addition to these general arguments, we give here a more illustrative explanation of the validity of the large-N kinetic theory in terms of (perturbative) concepts as cross sections and mean free paths.
We first get an expression for the cross section σ for elastic scattering in the presence of the vertex corrections encoded in the effective coupling. We start with a general Lorentz invariant expression of the differential cross section for 2 ↔ 2 scatterings in the center-of-mass system (CMS) (see, e.g., [58]) In the considered non-relativistic limit, the Mandelstam variable s simplifies to s ≈ 4m 2 . The matrix element for the perturbative kinetic theory can be expressed as |M| 2 ∼ m 4 g 2 /N , where we included the number of field components N explicitly. For the large-N kinetic theory, the coupling constant is replaced by the effective coupling g 2 → g 2 eff [f ] (t, ω, P ). In the CMS, in-and outgoing momenta are equal p = q and thus, one may set ω = 0.
The total cross section can be obtained from (54) by integration. For this, we substitute the integration variables from angular variables to P 2 = 2p 2 − 2p 2 cos θ pq with dΩ = 2π dcos θ pq = −(π/p 2 ) dP 2 , and arrive at Hence, for constant coupling g 2 as in the standard (perturbative) kinetic theory, one has the familiar expression 7 where we included the relativistic coupling λ for comparison. The same expression σ eff ≈ σ pert is found for the effective coupling whenever it reduces to g 2 , which occurs at sufficiently low typical occupancies ζ f (t, K) 1 or at large momenta p Q (Secs. VI A and V D). Otherwise g 2 eff [f ] (t, ω = 0, P ) is in the momentum-dominated regime that was discussed in Sec. V B. One then has Let us now proceed with the mean free path L of quasiparticles between two collisions for typical momenta K. A condition for a valid quasiparticle picture is that the mean free path L of quasiparticles is larger than their spatial extent, which can be estimated by their de Broglie wavelength λ DB ∼ 1/K [5], such that The mean free path of classical particles is usually estimated by means of the cross section and of the particle number density as L class ∼ 1/(σn). Since the frequency of collisions that a typical particle encounters may be enhanced by the Bose factor (1 + f (t, K)) [16], an alternative estimate of the mean free path is Because of L class ≥ L, the quasiparticle relation (59) is valid for L class whenever this is the case for L. Therefore, we will use L in Eq. (60) for our estimates.
On general grounds, we can already compare the large-N effective mean free path (referred to as L eff ) and the one resulting from standard (perturbative) kinetic theory L pert . Because of g 2 eff ≤ g 2 , one has σ eff ≤ σ pert ⇔ L eff ≥ L pert .
Hence, the effective coupling leads to a reduced cross section and an increased mean free path, counteracting in this way even high quasiparticle densities.
To be more specific, for the perturbative case of constant g 2 interactions, the mean free path can be estimated as The same expression is encountered for the large-N kinetic theory in the cases when the effective coupling g 2 eff reduces to g 2 . Otherwise, one has ζ f (t, K) 1 and limit. For an (ultra-)relativistic theory, one would have s ≈ 4p 2 , with p = |p|, which would lead to σ rel pert ∼ λ 2 /(p 2 N ). In equilibrium, typical momenta are of the order of the temperature scale p ∼ T , which would lead to familiar expressions for the cross section and the mean free path [5].
For perturbative kinetic theory, the mean free path (62) and the imposed restriction on occupation numbers m g K f (t, K) 1 in (53) imply L pert 1/K and condition (59) is satisfied.
In the case of the large-N kinetic theory, no restrictions on the occupation numbers have been posed. Whenever the effective coupling reduces to the coupling constant at momenta of the order K at vanishing ω, one has m g K f (t, K) 1 and for large N , this leads to L pert N/K 1/K. Otherwise, we can use (63) with L eff ∼ N/K 1/K, which is again consistent with the underlying assumption of large N and satisfies the condition (59).
In summary, based on the interference argument in Sec. VI B, typical occupation numbers are restricted to ζ f (t, K) 1 for perturbative kinetic theory. Since f (t, K) is a dynamical variable, this condition is time dependent. On the other hand, no restrictions on occupation numbers exist for large-N kinetic theory. In addition, Eq. (59) is satisfied in both cases, which contributes to a consistent quasiparticle interpretation.

VII. NONTHERMAL FIXED POINT FROM LARGE-N KINETIC THEORY
In this section, we determine the distribution function in the highly occupied self-similar regime at low momenta of scalar systems [6]. For the first time, we compute the universal scaling function f S (p) within the large-N kinetic theory, performed in three spatial dimensions. This is done numerically by solving a rescaled version of the kinetic equation that leads to the fixed point solution f S (p). We present our numerical approach and discuss the solution. Although we focus here on the self-similar regime, our numerical approach can be applied to more general situations of the thermalization process of scalar systems.
If not stated otherwise, all dimensionful quantities in this section are provided in suitable powers of 2m, and occupation numbers are rescaled by 4m 2 g. Hence, time, momenta and occupation numbers are obtained by the rescalings t → 2m t, p → p/(2m) and f → f /(4m 2 g).

A. Review of lattice simulation results
It was found in classical-statistical lattice simulations that non-relativistic and relativistic O(N ) symmetric scalar field theories for different number of components N approach a universal scaling regime at low momenta where the distribution function is large ζ f (t, K) 1 and follows a self-similar evolution [6,44,57] The scaling exponents α and β as well as the scaling function f S (p) are the same for a wide range of initial conditions once the evolution is sufficiently close to its nonthermal fixed point [37]. Since these properties are universal in the sense that they are approximately the same even in different theories, they form a far-fromequilibrium universality class [6,34]. In three spatial dimensions, the scaling exponents have been measured to be close to the values The scaling function f S (p) was observed to be approximately described by the functional form (45). This form can be approximated by two power laws for large and small momenta around a typical (stationary) infrared scale K S being the momentum where the particle number density is dominated according to (26) f S (p) The larger spectral exponent κ > has been investigated in different studies, often related to strong wave turbulence, and lies between 4 and 5 [6,37,39,40,44,45,59].
The lower exponent has been observed to lie within the range 0 ≤ κ < ≤ 1/2 [6,44]. Simulations within the 2PI framework confirm the existence of this nonthermal fixed point and show that it even survives at moderate couplings [37,60].
With the observed values for the scaling exponents α and β in (65), the dynamics can be understood as an inverse particle cascade that occupies low momenta and leads to the formation of a Bose-Einstein condensate far from equilibrium [61]. The condensate formation time grows with volume and the condensate growth follows a power law behavior in time [6,57].

B. Self-similar solution
Because of its large occupation numbers ζ f (t, K) 1, this scaling region cannot be described by perturbative kinetic theory, the range of validity of which is restricted by Eq. (53). Instead, the large-N kinetic theory can be applied. Indeed, it was shown in Ref. [6] that it correctly reproduces the observed scaling exponents (65). In the following, we extend this calculation by also obtaining the scaling function f S , and thus, we provide a complete scaling solution of the large-N kinetic theory in the low-momentum region.
To find the scaling form of the distribution function in the infrared, we first follow Refs. [6,9] and plug the self-similarity ansatz (64) into the effective kinetic equation (7). With the rescaled momentump = t β p, the left hand side can be expressed as Similarly, one can rescale the collision integral 8 Using ω p = p 2 /2m, one finds with the expressions for the collision integral (20), (21) and (22) Π R (t, ω p − ω q , P ) = t α−β Π R ωp − ωq,P , (69) which leads to the scaling exponent of the collision integral [6] Setting α −1 = µ, we get rid of the explicit time factors on both sides of the kinetic equation and we arrive at β = 1/2. This value of β is even the same in different spatial dimensions d [6]. Moreover, particle number density conservation leads to the additional relation α = 3β for the considered case of d = 3. Hence, the computed values agree with the observations (65). The remaining fixed point equation then reads We solve this equation numerically for the scaling function f S (p). While details of our numerical approach are postponed to the following subsections, we present here its solution in Fig. 6. The scaling function has the following properties: with the transition ("bending") scale at K S = 1, where we used its definition in Eq. (26). The exponents κ i of the distribution at larger and lower momenta compare well with former studies (see Sec. VII A). The functional form of the here obtained f S is so similar to Eq. (45), which was used in Ref. [6] to approximate f S on the lattice, that it is difficult to visually distinguish between them, as will be clearer in Fig. 7. We yet note that there are small deviations from that simple form in the bending region at the scale K S . Apart from the scaling form, we can also compare non-universal quantities that depend on parameters of our system like mass or particle number density. These are the amplitude f S (K S ) and the typical momentum K S . To be more specific, we compare these quantities to the classical-statistical lattice results from nonrelativistic scalar theory shown in Fig. 3 of Ref. [6]. There, the rescaled distribution function within the self-similar regime f S = (t/t ref ) −α f is plotted as a function of the rescaled momentum p = (t/t ref ) β p, where we used the overline notation to distinguish from our conventions for the scaling function and rescaled momentum. The two notations are related by The reference time was set to t ref = 300 in the same units as used here. With this, we can transform the typical momentum K S and the amplitude f S (K S ) obtained here within large-N kinetic theory to the overline notation, as and f S (K S ) = t α ref f S (K S ), and compare them to the corresponding quantities in Ref. [6]. Since in our prescription of the large-N kinetic theory for very high occupation numbers (see footnote 8) the amplitude drops out of the kinetic equation, it can thus be adjusted arbitrarily even after the simulation, depending on the particle number density in the system n ∼ Hence, it is not suitable for a comparison with the lattice results. On the other hand, being independent of n, the scale K S (or K S ) is fixed by the mass parameter m (or by both m and the reference time t ref ) and enables a quantitative comparison between large-N and lattice simulation results. In Fig. 3 of Ref. [6] the transition scale K S is located within the momentum range 0.05 K S 0.08. This is consistent with our result from the large-N kinetic theory (75).
We note that the small deviations between large-N kinetic and lattice results may have different reasons. First of all, we use the large-N kinetic theory at NLO and omit higher orders in 1/N . Moreover, the functional form measured on the lattice in Ref. [6] may suffer from finite-time effects and may slightly change at later times beyond the simulation times shown there. 9 However, we have seen that the large-N kinetic theory provides an even quantitatively good description of classical-statistical lattice data. This confirms its applicability to systems with very high occupation numbers, extending perturbative kinetic frameworks. The scaling function in Fig. 6 and its properties are the main results of this section.

C. Numerical setup
Here we discuss the numerical setup that led to the scaling solution in Fig. 6. In order to solve the fixed point equation (72), we start again with the full kinetic equation in (7) with the collision integral C[f ] as given by (23). Our strategy is to rescale the kinetic equation such that it relaxes to the fixed point equation with time. With this, we follow Ref. [26] where this strategy was used for the self-similar region at hard momenta in non-Abelian gauge theory. Therefore, instead of using a self-similarity ansatz, we rescale the distribution function and momenta according to with the values for the scaling exponents from (65). Comparing (76) to the self-similar evolution (64), one finds that the scaling function is the stationary limit of this rescaled distributionf (t,p) → f S (p). Plugging (76) into the kinetic equation (7), one arrives at the rescaled kinetic equation since the explicit time factors cancel. Note that Eq. (77) is equivalent to the original kinetic equation (7) but reduces to the fixed point form (72) whenf becomes time independent. Hence, a stationary solution of this equation corresponds to the scaling function f S (p), as has been noted above. In this sense, and because of the logarithmic time derivative ∂f /∂ log t, it can be regarded as a relaxation algorithm for f S (p) in time.
Moreover, the overall amplitude off drops out of the kinetic equation (77) because we neglect the 1 in the denominator of the effective coupling (21) (see also footnote (8)). This corresponds to the high-occupancy limitf (t,K) → ∞ withf (t,p)/f (t,K) kept fixed for each momentum. HereK is the typical (rescaled) infrared momentum scale where particle number density n is dominated with respect to the distributionf . Similarly, the coupling constant g also drops out of the kinetic equation.
To numerically solve (77), we discretize time logarithmically with constant ∆log t = log t k+1 − log t k = const between successive times t k+1 and t k . The time can then be calculated as t k = e k ∆log t . Moreover, we choose a logarithmically spaced momentum grid for the distribution functionf (t,p), its derivative and the collision integral, in order to resolve the distribution function at very low momenta. The grid involves N p momenta between Λ IR and Λ UV such that the ratio between successive momenta is constantp k+1 /p k = const. We employed ∆log t = 0.1, Λ IR = 0.017, Λ UV = 17 -52 and N p = 100 for the plots of this section.
For the computation of the collision integral and for the interpolation of the distribution function, we use methods 10 from the GNU Scientific Library [62]. Interpolation is required since the integration methods need continuous functions for the integrand, and it is performed based on the sampling points (p k ,f n,k ≡ f (t n ,p k )) at time t n . For momenta outside the momentum grid [Λ IR , Λ UV ], we set the distribution function to zero. Therefore, some terms in the functional I 2↔2 [f ] in (10) become zero when one of the momenta is outside of the momentum interval, and the collision integral loses its gain-minus-loss structure. To prevent this, we additionally set the whole functional I 2↔2 [f ] to zero in such cases to reduce deviations from particle number and energy density conservation. This leads to simplifications of integration boundaries within the collision integral (23), which are further discussed in App. B.
In the numerical algorithm, we first initialize the distribution functionf 0 and its momentum derivativef 0 at the grid pointsp k at initial time. The time step t n → t n+1 follows the explicit Euler method f n+1,k =f n,k − ∆log t αf n,k + βp kf n,k − C[f n ] .

(78)
Since the evaluation of the collision integral at a single momentum point p k neither depends on nor affects the evaluation at other momenta, we parallelize the part of our solver, where the collision integral is computed for each momentum on the grid.
The accuracy of our solution algorithm is mainly limited by the interpolation of the derivative of the distribution functionf . For a typical functional form as in Eq. (45), the relative accuracy for our discretization was up to 10 −2 as compared to the analytical expression of the derivative. Although increasing N p may improve the resolution, the computational costs will grow and we thus found a compromise that still provided sufficiently accurate results.

D. Details on the computation of the scaling function
In Fig. 7 we show the relaxation dynamics of the rescaled distributionf computed by the algorithm introduced above. We start close to the stationary form by choosingf 0 as in Eq. (45) with κ < = 0 and κ > = 3.9. One observes thatf quickly approaches a stationary form, which can be understood as the scaling function f S . Curves at times t ≥ 1.6 are already almost timeindependent. Therefore f S shown in Fig. 6 is computed as the average over these curves, while the error bars are estimated by the standard deviation in this procedure.
The functional form of f S is barely distinguishable from its starting form (45) in Fig. 7, however, small deviations around the scale K S exist. At low momentã p K S and at high momentap K S , the scaling function follows power lawsp −κ< andp −κ> . We have measured the spectral exponents κ i by employing power law fits to the respective regions in the scaling function, κ < = 0 ± 0.01 (sys) κ > = 3.95 ± 0.05 (sys) .
Statistical errors are much smaller than systematic errors, which were estimated to contain possible infrared and ultraviolet cutoff artifacts. 11 To check the stability of these values, we start with slightly larger exponents κ < = 0.5 and κ > = 4.5 for the initial distributionf 0 with the functional form (45). The time evolution of the exponents is shown in Fig. 8, where we use the error estimates of (79). Indeed, one observes that the exponents approach the values (79) of the scaling function.

VIII. CONCLUSION
In this work, we have shown that the large-N kinetic theory at NLO can be applied to highly occupied scalar quantum field theory, which cannot be described by a perturbative kinetic framework. On the other hand, for sufficiently low occupancies or at large momenta it effectively reduces to a perturbative kinetic theory at large N . Hence, the large-N kinetic theory extends perturbative descriptions and constitutes a versatile tool to study the dynamics and the thermalization process of systems out of equilibrium in terms of quasiparticles.
An essential ingredient for the quasiparticle picture, for free movement between collisions and for the proper inclusion of quantum interference effects even at high occupancies is the vertex resummation that appears at NLO in the large-N expansion [6,39,45,50,51]. We analyzed the structure of the effective vertex in detail, which enabled us to show that the mean free path L of quasiparticles at typical momentum modes is larger than their de Broglie wave length. Moreover, the spectral function can be approximated by a quasiparticle form at NLO of the large-N theory since a peak width is suppressed by 1/N . These arguments lead to a welldefined dispersion relation, and thus a consistent quasiparticle picture.
In a second step, we applied the large-N kinetic theory to the highly occupied region of scalar systems at low momenta characterized by a universal self-similar evolution [6,44]. Surpassing former analytical estimates for the scaling exponents α and β of the selfsimilar evolution [6], we solved the effective kinetic equation numerically for the first time. The scaling function obtained, f S (|p|), agrees well with former lattice simulation results, which is a striking confirmation of the applicability of the large-N kinetic theory to highly occupied systems. It reveals a power law behavior ∼ |p| −4 at momenta higher than the typical momentum K S that dominates particle number density, and becomes constant at low momenta.
So far, our study was performed for a non-relativistic scalar field theory under the assumption of a quadratic dispersion relation ω(|p|) ∼ |p| 2 . While this would correspond to a situation where the global condensate is small, this choice is also motivated by recent lattice studies that detected a dynamic exponent z ≈ 2 for the dispersion ω(|p|) ∼ |p| z [57]. Nonetheless, an interesting extension of the current form of the large-N kinetic theory would be to use the full Bogoliubov form [10] for the dispersion. Moreover, our study also captures relativistic scalar field theories with momenta below their effective mass |p| < M because they become effectively non-relativistic and since the studied infrared scaling region is located at momenta below M , which was also observed in lattice simulations [6]. To be able to simulate the complete thermalization process within large-N kinetic theory, including final thermalization, a similar numerical study can be done for relativistic systems.
No explicit assumption on the coupling strength has entered the derivation of the large-N kinetic theory. Therefore, in principle, one could use the large-N kinetic theory also at moderate couplings for finite N .
Indeed, it was shown using 2PI equations to NLO in a 1/N expansion [39,50,51] that the observed selfsimilar regime at low momenta survives to moderate values λ = 1 [60]. However, going to large couplings might become more tricky since we have assumed that the spectral function has a narrow quasiparticle peak that dominates the dynamics, which might not be true for strongly coupled theories. In addition, it was shown within the 2PI framework that for a large coupling λ = 10, inelastic processes may play an important role in the evolution [63]. The large-N kinetic theory in its present form, however, lacks such processes since they would be either off-shell or at next-to-next-to-leading order (NNLO) in a 1/N scheme.
The large-N kinetic theory is an example for a kinetic theory applicable to highly occupied systems, for which conventional kinetic approaches fail. For non-Abelian plasmas, multiple studies indicate strong fields and nontrivial dynamics at low-momentum modes [17][18][19][20][21]29]. An effective description thereof could be an important extension of kinetic approaches [7,[11][12][13]16] to the evolution of weakly coupled non-Abelian plasmas and the thermalization process in ultrarelativistic heavy-ion collisions at high energies.
FIG. 9. The filled region ∆(P ) represents the area of integration in (A8). The figure is taken with adaption from Ref. [6].
We have introduced the momentum difference P ≡ p − q for the effective coupling g 2 eff . In App. A 2 we will show that it only depends on the magnitude P = p 2 + q 2 − 2pq cos(θ p,q ), and hence on the magnitudes of the momenta p = |p| and q = |q| and on the angle θ p,q between them. The three-dimensional momentum integrals in (A1) have been split into radial and solid angle parts according to d 3 q = ∞ 0 dq q 2 dΩ q with dΩ q = 2π 0 dϕ q 1 −1 d cos (θ q ). All angular integrals have been included in dΓ C (p, l, q, r) ≡ dΩ p dΩ l dΩ q dΩ r Except for θ p,q , there is no angular dependence in the residual terms of the collision integral. In the following, the respective integrals are performed analytically.
We start by using the integral representation of the delta function Exploiting dΩ p dΩ q = dΩ q dΩ p,q , where dΩ p,q denotes the angular integration of p around the axis in q direction, and performing integrations over the angles except for θ p,q , we arrive at dΓ C (p, l, q, r) = 2 4 (2π) 5 where we employed 1 −1 dy e ±ixy = 2 sin(x)/x. We have used that the collision integral does not depend on the polar angle θ q of q, even after integrating over θ p,q , such that the integration A change of variables from l and r to u = r + l and v = r−l has been performed in the last step of Eq. (A8) absorbing a factor of 2. Fig. 9 visualizes the region of integration ∆(P ). A second transformation from d cos(θ p,q ) to dP with d cos(θ p,q ) = −(P/pq) dP leads to The next step is the evaluation of the energyconserving delta function. With the quadratic dispersion relation ω p = p 2 /2m, the delta function in the collision integral (A9) reads We use the v-integral to evaluate the delta function according to where R represents every part of the collision integral depending on v. The step function in (A11) can be used to change the integration boundaries of u to ∞ max(P,|p 2 −q 2 |/P ) du. The resulting form of the collision integral is given by Eq. (23).

Retarded self-energy
To also simplify the computation of the effective coupling g 2 eff of (21), we perform the angular integrations of the "one-loop" retarded self-energy Π R in (22) analytically. We first change the integration variable to P − k → k. Introducing y ≡ cos(θ P,k ), where θ P,k is the polar angle between the momenta P and k, one arrives at Π R (t, ω, P ) = lim Note that due to the isotropy of f (t, k), Π R only depends on the absolute value of the momentum P = |P|. Moreover, it vanishes for P = 0 since the integrand becomes identically zero. In addition, one observes that changing the frequency from ω to −ω corresponds to complex conjugation of the whole expression. As a consequence, the "one-loop" retarded self-energy Π R is real for ω = 0. Rewriting (A12) as Π R (t, ω, P ) = lim where we have substituted˜ ≡ m /P k, enables us to use the principal value integral for each fraction in Eq. (A13). This leads to Π R (t, ω, P ) = −mg (2π) 2 P ∞ 0 dk k f (t, k) Γ Π (ω, P, k) (A15) with the kernel Γ Π (ω, P, k) = log P 2 − 2P k 2 − 4m 2 ω 2 For a numerical treatment it is useful to know the singular points of the real part of the integrand Re Γ Π (ω, P, k) and the region where the imaginary part Im Γ Π (ω, P, k) does not vanish. The singularities of the real part are given by where the case P = 0 is excluded since the whole integrand is zero then. The real part can be rewritten as For the imaginary part Im Γ Π (ω, P, k) of the integrand (A16), we will assume ω > 0 since the case of a negative frequency is related to the positive frequency case by complex conjugation as noted above. Then, Im Γ Π (ω, P, k) is zero until the integration variable k has increased sufficiently to fulfill the condition k ≥ |P 2 − 2mω|/2P of the first Heaviside step function but is sill too small to fulfill the condition of the second step function. After exceeding |P 2 + 2mω|/2P , which makes the second step function one as well, the whole expression vanishes again. Altogether, we find − 1 π Im Γ Π (ω, P, k) = 1 if |P 2 −2mω| 2P ≤ k ≤ |P 2 +2mω| 2P 0 else .
(A19) Including also the case of negative frequency, we arrive at the final form of the "one-loop" retarded self-energy in Eq. (24).

Appendix B: Integration boundaries
As explained in Sec. VII C, the collision integral C[f ](t, p) as well as the distribution function f (t, p) are discretized on a grid between the momenta Λ IR and Λ UV in our numerical approach. Outside this domain, we set f (t, p) = 0. Making use of this, the integration boundaries of the collision integral (23) can be further constrained.
Integrals within the real and imaginary part of the "one-loop" retarded self-energy (24)  dk. (B2) The gain-minus-loss part of the integrand of the collision integral yields nonvanishing contributions if To preserve the gain-minus-loss symmetry, we set the whole I 2↔2 to zero if one of the conditions (B4) is not fulfilled, as explained in Sec. VII C. The first constraint is automatically fulfilled since the momentum p is externally set to the correct momentum range. The qintegration is changed to The relation with 2m|ω| ≡ p 2 − q 2 is chosen to be fulfilled for plus and minus sign simultaneously. It can be solved for u to obtain new integration boundaries according to the following procedure: In a first step, we consider the left inequality of (B6), i.e. (u ± 2m|ω|/u)/2 ≥ Λ IR . Using that u > 0, the inequality can be transformed to where the two cases of minus and plus sign on the righthand side are distinguished. For the case of the minus sign, no restriction for u is obtained if Λ 2 IR ≤ 2m|ω|. Otherwise, one obtains In case of the plus sign on the right-hand side of (B7), one finds the constraints The latter condition (B11) excludes u < Λ IR since u has always to be positive. In a second step, we consider the inequality on the right hand side of (B6), which can be transformed to This leads to the four conditions where the condition (B16) is no constraint, since Λ UV − Λ 2 UV + 2m|ω| ≤ 0. All considered cases and the successive constraints are employed at the same time. This allows us to change the integration range for u to