Large-N kinetic theory for highly occupied systems

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Large-N kinetic theory for highly occupied systems © Authors. Published by the American Physical Society. Funded by SCOAP3. Accepted version (Final draft) Walz, R.; Boguslavski, Kirill; Berges, J. Walz, R., Boguslavski, K., & Berges, J. (2018). Large-N kinetic theory for highly occupied systems. Physical Review D, 97(11), Article 116011. https://doi.org/10.1103/PhysRevD.97.116011 2018 Large-Nkinetic theory for highly occupied systems R. Walz,1,* K. Boguslavski,2,†and J. Berges1,‡ 1Institut für Theoretische Physik, Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany 2Department of Physics, University of Jyväskylä, P.O. Box 35, 40014 University of Jyväskylä, Finland (Received 5 February 2018; published 12 June 2018) 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 Nof the underlying scalar quantum field theory. Extending previous studies, we demonstrate that the large-Nkinetic 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-Nkinetic 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. DOI: 10.1103/PhysRevD.97.116011 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 twith position xand momentum p[1]. Accordingly, the derivation of kinetic theory from the underlying quantum field theory involves a series of crucial assumptions [2–6]. An important condition is that the de Broglie wavelength ∼1=jpjof 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 1 [2]. Likewise, perturbative descriptions exist far from equilibrium [8–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–13]. Much less is known about effective kinetic descriptions for general far-from-equilibrium situations. Pressing applications concern systems in which the occupancies of relevant modes are nonperturbatively 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 Refs. [14,15] and [11–13,16] for current perturbative kinetic descriptions). In this situation, typical gauge boson occupancies can become nonperturbatively large at low momenta below the Debye mass scale. These modes may influence the evolution of important quantities like the longitudinal pressure PLof the expanding plasma, evidence of which was found from real-time lattice simulations [17–19]. Recent studies have found universal scaling *[email protected] †[email protected] ‡[email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. PHYSICAL REVIEW D 97, 116011 (2018) 2470-0010=2018=97(11)=116011(22) 116011-1 Published by the American Physical Society 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 have been extensively discussed in the literature [22–29], as have methods for accessing spectral information at the Debye scale and below [30–33]. Remarkably, longitudinally expanding non-Abelian plasmas and self-interacting scalar field theories are found to share important universal aspects of their far-fromequilibrium evolution [17,34]. Similar self-similar 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–37], and nonrelativistic systems such as ultracold quantum gases or other condensed matter systems after a strong quench [38–41]. Characteristic infrared properties of these highly occupied systems turn out to be quantitatively the same for both relativistic and nonrelativistic 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 mean interactions are repulsive [47]. Perturbative kinetic approaches [48–50] break down at such large typical occupation numbers f≳1=λand are not able to reproduce key features of this low-momentum dynamics [6]. To describe the evolution also for nonperturbatively large occupation numbers, we consider an effective kinetic description for scalar systems that is not based on a weak-coupling or diluteness expansion. Developed in Refs. [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 nonperturbative expansions in the number of species available. For scalar field theories with Nspecies and quartic self-interactions, this results from a large-Nexpansion to next-to-leading order (NLO) based on a two-particle irreducible (2PI) resummation of self-energy diagrams [37,39,43,51,52], which translates to a vertex resummation in the kinetic framework [6,45]. So far, the large-Nkinetic theory at NLO has been successfully applied to analytically compute the selfsimilarity exponents near nonthermal fixed points at low momenta, agreeing well 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-Nkinetic 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 ∼jpj−4tail of the distribution in the regime dominated by number conservation. We analyze in detail the range of validity and quasiparticle picture of the large-Nor “vertex-resummed” 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 nonrelativistic low-energy limit relevant, e.g., also for the description of ultracold Bose gases. Section III summarizes main aspects of perturbative kinetic theory, before we present the large-Nkinetic 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. Vand VI. We present a numerical solution of the large-Nkinetic theory for the description of a nonthermal fixed point in Sec. VII. After concluding in Sec. VIII,we end with two Appendixes on calculational details of the collision integrals (Appendix A) and on integration boundaries (Appendix B). II. RELATIVISTIC AND NONRELATIVISTIC 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 xin three dimensions and quartic selfinteractions. For N¼4, the relativistic model describes the Higgs sector of the Standard Model of particle physics [53]. 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 multicomponent field theories [54]. In a nonrelativistic setting, the Heisenberg magnet for N¼3is a prominent example, and the case N¼2can be used to describe the two real components of a complex Bose field in systems of ultracold atomic gases dominated by s-wave scattering [55]. The considered relativistic quantum theory is described, on a classical level, by the action S½φ¼Zt;x1 2∂μφa∂μφa−m2 2φaφa−λ 4!NðφaφaÞ2ð1Þ with the notation Rt;x≡RdtRd3xand the (renormalized) mass mand coupling parameter λ. Here, summation over repeated Lorentz indices μ¼0;…;3and 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¼kB¼ℏ¼1. R. WALZ, K. BOGUSLAVSKI, and J. BERGES PHYS. REV. D 97, 116011 (2018) 116011-2 For processes with characteristic momenta below the mass scale m, one may expect an effectively nonrelativistic 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 nonrelativistic field theories. One may have in mind the phenomenologically important case of an N¼2-component nonrelativistic field theory, which can equivalently be described in terms of a complex field ϕðt; xÞ. Following standard procedures [56], the effective low-energy description may then be characterized by the nonrelativistic action Snr½ϕ;ϕ¼Zt;xϕi∂tþ∇2 2mϕ−g 2ðϕϕÞ2:ð2Þ Here we also introduced the effective nonrelativistic coupling g, which is no longer dimensionless and may be related to the relativistic parameters as [6] g∼λ m2:ð3Þ For dilute Bose systems, gcan be related to the s-wave scattering length, given by a¼mg=ð4πÞ[56]. For the nonrelativistic quantum theory in Eq. (2), the expectation value of the particle density n¼hϕϕiis conserved, and we will consider spatially homogeneous systems. The density nand scattering length acan be used to define a characteristic “coherence length”of which the inverse is the momentum scale Q¼ffiffiffiffiffiffiffiffiffiffiffiffiffi 16πan p∼ffiffiffiffiffiffiffiffiffi mgn p:ð4Þ We also define the “diluteness parameter” ζ¼ffiffiffiffiffiffiffiffi na3 p∼Qmg; ð5Þ which provides a dimensionless expansion parameter for the nonrelativistic system, similar to the dimensionless coupling λfor the relativistic system. Specifically, with (3), we obtain ζ∼ðQ=mÞλ. We emphasize that for the relativistic theory the particle number is not conserved in general. However, for the highly occupied system considered in Sec. VII, an approximately conserved particle number is dynamically generated at a nonthermal fixed point such that the nonrelativistic theory and the relativistic one can be in the same universality class of infrared scaling phenomena [6]. 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 λ≪1or 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 tand 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 anticommutator expectation value hfφa;φbgi≡hφaφbþφbφaidetermines (in the absence of external forces) the change of the distribution function according to [45] Z∞ 0 dω 2πω∂ ∂thfφa;φbgiðt; ω;pÞ≡∂fðt; pÞ ∂tδab:ð6Þ Here, the frequency ωand spatial momentum parise 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 hfφa;φbgi∼δ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 ∂fðt; pÞ ∂t¼C½fðt; pÞ:ð7Þ 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) lead to the well-known Boltzmann equation with the collision integral for elastic 2↔2scatterings [56], Crel½fðt;pÞ ¼Zl;q;r λ2ðNþ2Þ 6N2I2↔2½fðt;p;l;q;rÞ ×ð2πÞ4δð3Þðpþl−q−rÞδðωrel pþωrel l−ωrel q−ωrel rÞ 2ωrel p2ωrel l2ωrel q2ωrel r ;ð8Þ with the notation Rq≡Rd3q=ð2πÞ3and the relativistic dispersion ωrel p¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m2þp2 q:ð9Þ The functional I2↔2½fcontains the distribution functions fp≡fðt; pÞdescribing the changes by loss or gain through scattering: LARGE-NKINETIC THEORY FOR HIGHLY OCCUPIED …PHYS. REV. D 97, 116011 (2018) 116011-3 I2↔2½fðt; p;l;q;rÞ ¼ðfpþ1Þðflþ1Þfqfr−fpflðfqþ1Þðfrþ1Þ ≈ f≫1ðfpþflÞfqfr−fpflðfqþfrÞ:ð10Þ 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 leads to a quasiparticle form for the spectral function given by the expectation value of the commutator of two fields [8,45]: h½φa;φbiðω;pÞ¼ Oðλ0Þ2πsgnðωÞδðω2−ðωrel pÞ2Þδab:ð11Þ From higher orders in the coupling, the spectral function would receive corrections leading to a mass shift and nonzero 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, we do not consider them here. Since there are only elastic collisions contributing to this order, the particle number is artificially conserved. Inelastic processes can also be taken into account by going to higher order in the coupling [44]. Starting from a general out-ofequilibrium state, such inelastic processes are relevant to describe the approach to thermal equilibrium at late times, since otherwise a thermal distribution with chemical potential for the particle number would appear even in the absence of a conserved number. Essentially, neglecting inelastic processes limits the time until which the approximation can be applied [57]. For the purposes of this section, going beyond the given order is not necessary.1 Similarly, for the nonrelativistic dispersion ωp¼ jpj2=2mat lowest perturbative order, the collision term for the kinetic equation of the theory with action (2) becomes [10] Cnr½fðt; pÞ¼Zl;q;r 2g2I2↔2½fðt; p;l;q;rÞð2πÞ4 ×δð3Þðpþl−q−rÞ ×δðωpþωl−ωq−ωrÞ:ð12Þ The quadratic dispersion relation at this order can also be viewed as arising from the low-momentum limit of the above relativistic collision integral (8) [9,50]. We note that taking into account subleading corrections for nonrelativistic 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], which we do not consider here. The perturbative power counting for λ≪1leading to (8),or(12) for ζ≪1, with elastic 2↔2scatterings as in (10) assumes that the relevant occupancies fpfor typical momenta are not too high. More precisely, only for fp≪1=λin the relativistic, or fp≪1=ζin the nonrelativistic case, the higher-order corrections are parametrically small. Before we discuss this issue in more detail below in Sec. VI, we will introduce in the following an alternative kinetic description based on a large-Nexpansion of the underlying quantum field theory. IV. LARGE-NKINETIC THEORY The standard Boltzmann equation described in the last section is based on a weak-coupling expansion, which restricts the range of validity of the kinetic theory to perturbative problems. However, for the N-component field theory, an alternative kinetic description with an extended range of validity may be derived based on a nonperturbative expansion in N. For details about its derivation from the underlying quantum field theory, we refer to Refs. [6,56]. Here, we give the relevant expressions that are used below to solve the large-Nkinetic 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,52]. In particular, the spectral function at leading order (LO) in a large-Nexpansion reads h½φa;φbiðω;pÞ¼ LOlargeN2πsgnðωÞδðω2−ðωrel pÞ2Þδab;ð13Þ where, in contrast to the lowest-order perturbative Eq. (9), the dispersion for the relativistic theory now contains an effective mass term M2: ωrel p¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M2þp2 q:ð14Þ At LO, the effective mass term is given by the gap equation [6,52] M2¼m2þλ 6Zp fðt0;pÞ ωrel pð15Þ evaluated at some given time t0. 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-Nexpansion there is a well-defined effective kinetic description in terms of scatterings between quasiparticles. Similar to the previous section, we start by considering the 1For instance, for a relativistic scalar field theory describing interacting pions in the context of heavy-ion collisions, conditions for a conserved particle number density and the time interval in which this approximation can be trusted have been discussed in Refs. [58–60]. R. WALZ, K. BOGUSLAVSKI, and J. BERGES PHYS. REV. D 97, 116011 (2018) 116011-4 relativistic theory and extend the discussion to the nonrelativistic case in the end. To discuss subleading corrections in the 1=N expansion, it is convenient to employ the auxiliary field formulation of the same model [52]. For this purpose, we rewrite the original action (1) by introducing an auxiliary field χðxÞas S½φ;χ¼−Zt;x1 2φað□þm2Þφa− 3N 2λχ2þ1 2χφaφa: ð16Þ 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 χðxÞ¼ λ 6NφaðxÞφaðxÞ:ð17Þ While the auxiliary field is not a dynamical degree of freedom, it can be used to conveniently express scattering corrections in terms of the expectation value Dðx−yÞ≡hχðxÞχðyÞi−hχðxÞihχðyÞi:ð18Þ Since χrepresents a two-point function according to (17), the function Dðx; yÞencodes a four-point function or vertex. Specifically, at NLO in the 1=N expansion, scatterings are mediated by (18) [52]. This is indicated in Fig. 1,in which dashed lines represent the two-point function (18) in Fourier space. The modified vertex at NLO is shown for scatterings in the s,t, and uchannels, respectively. Accordingly, the effective kinetic equation at NLO is given by the same kinetic equation (7), however with the different collision term [6] Crel NLO½fðt; pÞ¼Zl;q;r λ2 effðt; p;l;q;rÞ 6NI2↔2½fðt; p;l;q;rÞ ×ð2πÞ4δð3Þðpþl−q−rÞ ×δðωrel pþωrel l−ωrel q−ωrel rÞ 2ωrel p2ωrel l2ωrel q2ωrel r :ð19Þ 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 timeand momentum-dependent effective coupling function λ2 effðt;p;l;q;rÞ≡λ2 31 j1þΠrel Rðt;ωrel pþωrel l;pþlÞj2 þ1 j1þΠrel Rðt;ωrel p−ωrel q;p−qÞj2 þ1 j1þΠrel Rðt;ωrel p−ωrel r;p−rÞj2ð20Þ incorporates the vertex corrections for the different scattering channels according to Fig. 1. The appearance of the renormalized one-loop retarded self-energy Πrel Rðt;ω;pÞ¼ λ 12 Zq fðt;p−qÞ ωrel qωrel p−q1 ωrel qþωrel p−q−ω−iϵ þ1 ωrel q−ωrel p−q−ω−iϵþ1 ωrel q−ωrel p−qþωþiϵ þ1 ωrel qþωrel p−qþωþiϵð21Þ in the denominator of Eq. (20) is the result of a geometric series summation of an infinite number of scattering processes at NLO in the large-Nexpansion [51].We emphasize that Πrel R, and thus also λ2 eff, is time dependent since it depends on the evolving distribution function. From (20), one observes that for jΠrel Rj≪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 (19) essentially describes standard perturbative 2↔2scatterings, however at large N.2In contrast, for high characteristic occupancies with f∼1=λ, the collision term (19) can be strongly modified if Πrel Rstarts to become of order 1. We will discuss the corresponding behavior of the effective coupling in more detail in the following sections, for which we will introduce the nonrelativistic effective kinetic equation relevant at low momenta below. Similar to the lowest-order perturbative kinetic equation of Sec. III, the large-Nkinetic theory at NLO only involves elastic scattering processes. The lack of inelastic processes implies conservation of the particle number density n¼Rpfðt; pÞ¼const, which follows from the kinetic equation (7) and RpCrel NLO½fðt; pÞ¼0. Taking into account inelastic processes is possible by going beyond NLO, which is, however, beyond the scope of the present study ppp lll qqq rr r p+l p-q p-r FIG. 1. Scattering processes at next-to-leading order in the large-Nexpansion, which are mediated by an effective interaction. While the solid lines represent particles with given 4momenta, the dashed line represents the function (18) in Fourier space describing s-, t-, and u-channel exchange. 2The prefactor ðNþ2Þ=ð6N2Þin Eq. (8) becomes 1=ð6NÞat NLO in the large-Nexpansion. LARGE-NKINETIC THEORY FOR HIGHLY OCCUPIED …PHYS. REV. D 97, 116011 (2018) 116011-5 that aims to provide a kinetic description of numberconserving dynamics near nonthermal fixed points [6,44,47]. Moreover, here, the condensate formation time diverges with volume ∼V1=αwith positive scaling exponent αas shown in lattice simulations [6]. Hence, no emergence of a condensate is expected within a finite time for the infinite volume considered, which is consistent with the self-similar evolution that we will observe in numerical calculations of the large-Nkinetic theory in Sec. VII. Therefore, we will consider the collision integral (19) without including a condensate in the following. Further discussions on the distinction between the perturbative and nonperturbative regimes can be found in Ref. [61]. To simplify the following discussion and to make the connection to ultracold atoms, we will restrict ourselves to momenta below the (effective) mass scale jpj≪Mand describe the dynamics in terms of a nonrelativistic quantum field theory. Note that also relativistic theories with m¼0but M>0may be described by the nonrelativistic limit for small momenta. To be consistent with the nonrelativistic theory defined by (2), we will use the symbol mfor the mass. Following along the lines of Sec. III, we consider first the case N¼2to illustrate the effective kinetic equation for a nonrelativistic complex scalar field, i.e., with two real field components. The case of general Nthen proceeds accordingly [62]. For the quadratic dispersion relation, one obtains [6] Cnr NLO½fðt;pÞ ¼Zl;q;r g2 eff½fðt;ωp−ωq;p−qÞI2↔2½fðt;p;l;q;rÞ ×ð2πÞ4δð3Þðpþl−q−rÞδðωpþωl−ωq−ωrÞ:ð22Þ The effective coupling in the collision integral reads g2 eff½fðt; ω;PÞ¼ g2 j1þΠRðt; ω;PÞj2;ð23Þ with the one-loop retarded self-energy ΠRðt;ω;PÞ¼lim ϵ→0þgZk fðt;P−kÞ ×1 ωk−ωP−k−ω−iϵþ1 ωk−ωP−kþωþiϵ ð24Þ and the momentum difference P¼p−q. As for the relativistic theory, the collision integral (22) reduces to its perturbative expression3for small jΠRj≪1. For later use, it is helpful to further evaluate the expressions for Cnr NLO and ΠR. Using magnitudes of momenta p¼jpj, and similarly for q,k, and P, we find for an isotropic system (details are given in Appendix A) C½fðt; pÞ¼ m 32π3pZ∞ 0 dqq Zpþq jp−qj dP ×g2 eff½fðt; ωp−ωq;PÞ ×Z∞ max ðP;jp2−q2j PÞ duu−ðp2−q2Þ2 u3 ×I2↔2½ft; p; 1 2u−p2−q2 u; q; 1 2uþp2−q2 u ð25Þ and ΠRðt; ω;PÞ¼ mg ð2πÞ2PZ∞ 0 dkkfðt; kÞ × log ðkþP 2Þ2−m2ω2 P2 ðk−P 2Þ2−m2ω2 P2 þiπZjP2þ2mωj 2P jP2−2mωj 2P dkkfðt; kÞ;ð26Þ where we have dropped the labels of Cnr NLO to shorten the notation. Accordingly, the effective kinetic equation (7) depends on time tand the magnitude of the momentum p. V. BEHAVIOR OF THE LARGE-NRESUMMED EFFECTIVE VERTEX To discuss the extended range of validity of large-N kinetic theory, we first consider the effective coupling g2 eff appearing in the collision integral (25), which is a function of the difference in energies of the inand outgoing particles ω¼ωp−ωq¼ðp2−q2Þ=2mand of the magnitude of the momentum change P¼jp−qjin 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≫qand thus 2mω≈p2≈P2. Similarly, q≫pleads to the same effective coupling because g2 eff½fðt; ω;PÞis symmetric in ω. To ease the discussion, we will therefore use ω≥0in the following. In the second case, the magnitude of the momentum of the particle undergoing the collision stays at the same order p∼qwhile it changes its direction. The nearly collinear regime is discussed separately and constitutes the third scenario. The effective coupling in (23) can be calculated from the retarded self-energy ΠRðt; ω;PÞgiven in (26). For the 3The factor of 2 difference between the large-NNLO expression (22) and the perturbative collision integral in (12) is due to an omitted term that is of order NNLO. Note that a similar modification is found for the relativistic theory, as commented on in footnote 2. R. WALZ, K. BOGUSLAVSKI, and J. BERGES PHYS. REV. D 97, 116011 (2018) 116011-6 distribution function fðt; kÞthat enters the integrals in ΠR, we assume that for l¼0, 1, 2 integrals of the form Zdkklfðt; kÞ∼Klþ1fðt; KÞð27Þ 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 in which k2fðt; kÞis maximal, k2fðt; kÞjk¼K¼max ðk2fðt; kÞÞ;ð28Þ such that it provides the dominant contributions to the particle number density n∝Zd3k ð2πÞ3fðt; kÞ∼K3fðt; KÞ:ð29Þ For the integrals in (27) to converge between the maximal limits of 0 and ∞, the distribution function fðt; kÞshould fall off faster than k−3at large momenta k≳Kand decrease more slowly than k−1at low momenta k≲K, or not decrease there at all. A. Dispersive regime, P2≈2mω We start with the regime p≫q. Then, one has P≈p and ω≈ωp¼p2=2m≈P2=2m. The second relation states that the energy difference of inand outgoing momenta ω follows the nonrelativistic dispersion relation with momentum (difference) P, and we refer to this regime as dispersive. The one-loop retarded self-energy (26) in this limit reads ΠRt; P2 2m;P¼mg ð2πÞ2PZ∞ 0 dkkfðt;kÞlog kþP k−P þiπZP 0 dkkfðt;kÞ:ð30Þ Its real part involves an integration over all momenta, and we can use (27) in the limiting cases of P≳Kand P≲K. In the first case, the integrand is dominated at low momenta k, and we can approximate logðkþPÞ−log jk−Pj≈ 2k=P þOððk=PÞ3Þ. Similarly, the second case leads to logðkþPÞ−log jk−Pj≈2P=k þOððP=kÞ3Þ. Hence, the limiting expressions are ReΠRt; P2 2m;P ∼ P≳KmgKfðt; KÞK2 P2ð31Þ ReΠRt; P2 2m;P ∼ P≲KmgKfðt; KÞ:ð32Þ For the imaginary part, large and small ingoing momenta P≳Kand P≲Klead to the expressions ImΠRt; P2 2m;P ∼ P≳KmgKfðt; KÞK Pð33Þ ImΠRt; P2 2m;P ∼ P≲KmgPfðt; PÞ:ð34Þ In (34), we used that kfðt; kÞshould be a growing function at low momenta to be consistent with (27). To get the corresponding limiting expressions for the effective coupling, we first assume that for typical soft momenta Kthe occupation number fðt; KÞis sufficiently large such that for the considered momenta Pthe 1 in the denominator of g2 eff in (23) can be neglected, and the effective coupling reads g2 eff ≈g2ððReΠRÞ2þðImΠRÞ2Þ−1. Since Pfðt; PÞis limited by Kfðt; KÞ, the real part (32) dominates at low momenta P≲K. On the other hand, the imaginary part (33) decreases more slowly than the real part at high momenta P≳Kand is thus larger. With this, the effective coupling is parametrically g2 eff½ft; P2 2m;P ∼ P≳K1 ðmKfðt; KÞÞ2 P2 K2ð35Þ g2 eff½ft; P2 2m;P ∼ P≲K1 ðmKfðt; KÞÞ2:ð36Þ Accordingly, it is constant below Kand follows the power law P2beyond K. We note that at even larger momenta it becomes constant ≃g2when the þ1in the denominator of its definition becomes important. As noted above, the regime q≫pleads to the same expressions (35),(36), with P≈q. The frequency argument gets a minus sign −q2=2m, which does not change the values of g2 eff because of its symmetry. B. Momentum-dominated regime, P2≳2mω In the momentum-dominated regime, we consider the situation in which inand outgoing momenta are of the same order p∼q. For further simplifications, we also assume that the frequency difference is small, ω¼ωp−ωq≲P2=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. VC. In the considered regime, the imaginary and real parts of the one-loop retarded self-energy (26) become ReΠRðt;ω;PÞ≈mg ð2πÞ2PZ∞ 0 dkkfðt;kÞ2log 2kþP 2k−Pð37Þ ImΠRðt; ω;PÞ≈mg 4π P 2ft; P 22mω P2:ð38Þ LARGE-NKINETIC THEORY FOR HIGHLY OCCUPIED …PHYS. REV. D 97, 116011 (2018) 116011-7 Expanding the logarithm of the real part and approximating the integral as in (27), one arrives at similar expressions as in the dispersive case ReΠRðt; ω;PÞ∼ P≳2KmgKfðt; KÞð2KÞ2 P2ð39Þ ReΠRðt; ω;PÞ∼ P≲2KmgKfðt; KÞ:ð40Þ A closer look on the logarithm of the original expression in (26) reveals that if 2mω≲2KP is satisfied the estimate for lower momenta (40) is even valid in the collinear regime where 2mωexceeds P2. Hence, the full range of validity for (40) is 2K≳P≳2mω=2K, which may only hold if ð2KÞ2≳2mω. Otherwise, for ð2KÞ2≲2mω, no region with the value (40) exists. To compute the effective coupling, we again assume large occupation numbers and thus neglect the 1 in the denominator in (23), which yields g2 eff ≈g2ðð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ω≲P2and, for small momenta P=2≲K, from ðP=2Þfðt; P=2Þ≲Kfðt; KÞ, while for larger momenta P=2≳K, it results from ðP=2Þ3fðt; P=2Þ≲ K3fðt; KÞ, which are both requirements for fðt; kÞand were formulated below Eq. (27). Hence, the effective coupling parametrically follows g2 eff½fðt; ω;PÞ∼ P≳2K1 ðmKfðt; KÞÞ2 P4 ð2KÞ4ð41Þ g2 eff½fðt; ω;PÞ∼ P≲2K1 ðmKfðt; KÞÞ2:ð42Þ The main difference from the dispersive regime is the steep power law P4at large momenta. Interestingly, the transition between smalland large-momentum expressions proceeds at the slightly larger scale 2K. C. Collinear regime, P2≲2mω The remaining case is when inand outgoing momenta are nearly collinear, cos θpq ∼1, i.e., the case in which P2≲2mω. The logarithm appearing in ReΠRin (26) can be written as log  1−P2ð2kþPÞ2 ð2mωÞ2 −log  1−P2ð2k−PÞ2 ð2mωÞ2 ≈−4P4 ð2mωÞ2 2k P;ð43Þ where we have expanded it in the second line. Assuming that the integral is dominated at momenta k∼Kas 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 ΠRas ReΠRðt; ω;PÞ∼−mgKfðt; KÞð2KÞ2 ðPinvÞ2ð44Þ ImΠRðt; ω;PÞ≈mg 4π Pinv 2ft; Pinv 2;ð45Þ where we have introduced the inverse momentum Pinv ¼2mω=P. Recall that for the real part the remaining situation of 2mω≲2KP for low momenta P≲2Kin the collinear regime has been discussed in Sec. VB, in which it led to the expression (40). From this, we can compute the effective coupling in the collinear regime. Interestingly, the real part is negative, and 1þReΠRin the denominator of the effective coupling (23) 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 (23) can again be neglected, and the real and imaginary parts of ΠRcan be compared in order to estimate g2 eff ≈g2ððReΠRÞ2þðImΠRÞ2Þ−1. At first, we consider Pinv ≳2K. This condition translates to 2mω≳ 2KP and corresponds to the real part as given by (44). Using ðPinvÞ3fðt; PinvÞ≲K3fðt; KÞfor large momenta Pinv, one finds ImΠR≲jReΠRj. Similarly, the case Pinv ≲ 2Ktranslates to 2mω≲2KP, and the real part is then given by (40). With Pinvfðt; PinvÞ≲Kfðt; KÞfor low momenta, one again finds ImΠR≲jReΠRj. Thus, as in the momentum-dominated regime, the real part dominates the effective coupling. At low momenta P≲2K, it leads to (42) for 2mω≲2KP, while in all other cases, in the collinear regime, one has g2 eff½fðt; ω;PÞ∼ 1 ðmKfðt; KÞÞ2ð2mωÞ4 ð2KPÞ4:ð46Þ 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 g2 eff numerically as defined in (23). For the distribution function, we use fðt; pÞ≃ 1 g A ðp=BÞκ<þðp=BÞκ>;ð47Þ which has been suggested in Ref. [6] to approximate the distribution function at low momenta during the self-similar regime. The parameter Bis related to the momentum scale Kdefined in Eq. (28) via K¼Bðð2−κ<Þ=ðκ>−2ÞÞ1=ðκ>−κ<Þ, such that both are of R. WALZ, K. BOGUSLAVSKI, and J. BERGES PHYS. REV. D 97, 116011 (2018) 116011-8 momentum KSand the amplitude fSðKSÞobtained here within large-Nkinetic theory to the overline notation, as KS¼t−β ref KS≈0.0577 ð76Þ and fSðKSÞ¼tα reffSðKSÞ, and compare them to the corresponding quantities in Ref. [6]. Since in our prescription of the large-Nkinetic theory for very high occupation numbers (see footnote 9) 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∼fSðKSÞK3 S¼fSðKSÞKS 3. Hence, it is not suitable for a comparison with the lattice results. On the other hand, being independent of n, the scale KS(or KS) is fixed by the mass parameter m(or by both mand the reference time tref) and enables a quantitative comparison between large-Nand lattice simulation results. In Fig. 3 of Ref. [6], the transition scale KSis located within the momentum range 0.05 ≲KS≲0.08. This is consistent with our result from the large-Nkinetic theory (76). We note that the small deviations between large-N kinetic and lattice results may have different reasons. First of all, we use the large-Nkinetic 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.10 And finally, regarding the discussion below Eq. (67), the observed power law on the lattice could be a superposition of power laws with different origins. Therefore, it was of great importance to pin down the power-law exponent 4 in Eq. (74) that can be associated with the large-Nkinetic theory11 to be able to distinguish it from other possible contributions. We have seen that the large-Nkinetic 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. 6and 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. To solve the fixed-point equation (73), we start again with the full kinetic equation in (7) with the collision integral C½fas given by (25).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] in which 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 fðt; pÞ≡tα˜ fðt; tβpÞ≡tα˜ fðt; ˜ pÞ;ð77Þ with the values for the scaling exponents from (66). Comparing (77) to the self-similar evolution (65), one finds that the scaling function is the stationary limit of this rescaled distribution ˜ fðt; ˜ pÞ→fSð˜ pÞ. Plugging (77) into the kinetic equation (7), one arrives at the rescaled kinetic equation ∂ ˜ fðt; ˜ pÞ ∂logt¼−α ˜ fðt; ˜ pÞþβ˜ p∂ ˜ fðt; ˜ pÞ ∂˜ p−C½ ˜ fðt; ˜ pÞ;ð78Þ since the explicit time factors cancel. Note that Eq. (78) is equivalent to the original kinetic equation (7) but reduces to the fixed-point form (73) when ˜ fbecomes time independent. Hence, a stationary solution of this equation corresponds to the scaling function fSð˜ 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 fSð˜ pÞin time. Moreover, the overall amplitude of ˜ fdrops out of the kinetic equation (78) because we neglect the 1 in the denominator of the effective coupling (23) (see also footnote 9). This corresponds to the high-occupancy limit ˜ fðt; ˜ KÞ→∞with ˜ fðt; ˜ pÞ= ˜ fðt; ˜ KÞkept fixed for each momentum. Here, ˜ Kis the typical (rescaled) infrared momentum scale in which particle number density nis dominated with respect to the distribution ˜ f. Similarly, the coupling constant galso drops out of the kinetic equation. To numerically solve (78), we discretize time logarithmically with constant Δlog t¼log tkþ1−log tk¼const between successive times tkþ1and tk. The time can then be calculated as tk¼ekΔlog t. Moreover, we choose a logarithmically spaced momentum grid for the distribution function ˜ fðt; ˜ pÞ, its derivative, and the collision integral, in order to resolve the distribution function at very low momenta. The grid involves Npmomenta between ΛIR and ΛUV such that the ratio between successive momenta is constant ˜ pkþ1=˜ pk¼const. We employed Δlog t¼0.1, ΛIR ¼0.017,ΛUV ¼17–52, and Np¼100 for the plots of this section. For the computation of the collision integral and for the interpolation of the distribution function, we use methods12 from the GNU Scientific Library [71]. Interpolation is required since the integration methods need continuous 10Similar finite-time artifacts have been observed for nonAbelian gauge theory, in which a kinetic description compared well with lattice results but showed slight differences at late times beyond the time of lattice simulations when being closer to the nonthermal fixed point [26]. 11See also Ref. [61] for an alternative approach to this problem and an outcome consistent with our results. 12From Ref. [71], we frequently use the integration method GSL _ INTEGRATION _ CQUAD , which is particularly suitable for singular integrands. Such occur, for instance, in the real part of the one-loop retarded self-energy ΠR. The employed interpolation type is GSL _ INTERP _ AKIMA . LARGE-NKINETIC THEORY FOR HIGHLY OCCUPIED …PHYS. REV. D 97, 116011 (2018) 116011-15 functions for the integrand, and it is performed based on the sampling points ð˜ pk; ˜ fn;k ≡ ˜ fðtn;˜ pkÞÞ at time tn.For momenta outside the momentum grid ½ΛIR;ΛUV,weset the distribution function to zero. Therefore, some terms in the functional I2↔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 I2↔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 (25), which are further discussed in Appendix B. In the numerical algorithm, we first initialize the distribution function ˜ f0and its momentum derivative ˜ f0 0at the grid points ˜ pkat initial time. The time step tn→tnþ1 follows the explicit Euler method ˜ fnþ1;k ¼ ˜ fn;k −Δlog t½α ˜ fn;k þβ˜ pk ˜ f0 n;k −C½ ˜ fn:ð79Þ Since the evaluation of the collision integral at a single momentum point pkneither depends on nor affects the evaluation at other momenta, we parallelize the part of our solver in which 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 function ˜ f0. For a typical functional form as in Eq. (47), the relative accuracy for our discretization was up to 10−2as compared to the analytical expression of the derivative. Although increasing Npmay 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 distribution ˜ fcomputed by the algorithm introduced above. We start close to the stationary form by choosing ˜ f0as in Eq. (47) with κ<¼0and κ>¼3.9.One observes that ˜ fquickly approaches a stationary form, which can be understood as the scaling function fS. Curves at times t≥1.6are already almost time independent. Therefore, fSshown in Fig. 6is computed as the average over these curves, while the error bars are estimated by the standard deviation in this procedure. The functional form of fSis barely distinguishable from its starting form (47) in Fig. 7; however, small deviations around the scale KSexist. At low momenta ˜ p≲KSand at high momenta ˜ p≳KS, the scaling function follows power laws ˜ p−κ<and ˜ p−κ>. We have measured the spectral exponents κiby employing power-law fits to the respective regions in the scaling function, κ<¼00.01ðsysÞ κ>¼3.95 0.05ðsysÞ:ð80Þ Statistical errors are much smaller than systematic errors, which were estimated to contain possible infrared and ultraviolet cutoff artifacts.13 To check the stability of these values, we start with slightly larger exponents κ<¼0.5and κ>¼4.5for the initial distribution ˜ f0with the functional form (47). The time evolution of the exponents is shown in Fig. 8, where we use the error estimates of (80). Indeed, one observes that the exponents approach the values (80) of the scaling function. VIII. CONCLUSION In this work, we have shown that the large-Nkinetic 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-Nkinetic theory extends perturbative descriptions and constitutes a versatile tool to study the dynamics 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-Nexpansion [6,39,45,51,52]. We analyzed the structure of the effective vertex in detail, which enabled us to show that the mean free path Lof quasiparticles at typical momentum modes is larger than their de Broglie wavelength. Moreover, the spectral function can be FIG. 8. The exponents of approximate power laws ˜ p−κiat low (κ<) and high (κ>) momenta of the rescaled distribution ˜ fas functions of time. The distribution has been initialized as in (47) with A¼27000,B¼1.135,κ<¼0.5, and κ>¼4.5. 13We note that at low and high momenta close to the cutoffs the scaling function starts to show deviations from power-law or constant behavior, which is the main source of error in the powerlaw fits. Increasing the momentum grid to lower and larger momenta reduces the effects of these numerical artifacts but comes at the price of increased numerical costs. R. WALZ, K. BOGUSLAVSKI, and J. BERGES PHYS. REV. D 97, 116011 (2018) 116011-16 approximated by a quasiparticle form at NLO of the largeNtheory since a peak width is suppressed by 1=N. These arguments lead to a well-defined dispersion relation and thus a consistent quasiparticle picture. In a second step, we applied the large-Nkinetic 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 self-similar evolution [6], we solved the effective kinetic equation numerically for the first time. The scaling function obtained, fSðjpjÞ, 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 ∼jpj−4at momenta higher than the typical momentum KSthat dominates particle number density and becomes constant at low momenta. No explicit assumption on the coupling strength has entered the derivation of the large-Nkinetic theory. Therefore, in principle, one could use the large-Nkinetic theory also at moderate couplings for finite N. Indeed, it was shown using 2PI equations to NLO in a 1=N expansion [39,51,52] that the observed self-similar regime at low momenta survives to moderate values λ¼1[68]. However, it was argued within the 2PI framework that for a large coupling λ¼10 inelastic processes may play an important role at all times in the evolution [72].Thelarge-N kinetic theory in its present form at NLO, however, lacks such processes. To be able to simulate a complete thermalization process within large-Nkinetic theory, starting far from equilibrium and the evolution toward a nonthermal fixed point including the subsequent final thermalization dynamics, one would have to go beyond this order to capture inelastic processes. While going to NNLO is challenging, the description may be partially simplified at late times relevant for the final approach to thermal equilibrium since the typical occupancies become smaller such that standard perturbative approximations become available again at least for small enough couplings. The large-Nkinetic 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–21,29].An effective description thereof could be an important extension of kinetic approaches [7,11–13,16] to the evolution of weakly coupled non-Abelian plasmas and the thermalization process in ultrarelativistic heavy-ion collisions at high energies. ACKNOWLEDGMENTS We thank J. P. Blaizot, I. Chantesana, T. Gasenzer, A. Kurkela, T. Lappi, A. Piñeiro Orioli, S. Schlichting, and R. Venugopalan for useful discussions and collaborations on related work. K. B. gratefully acknowledges support by the European Research Council under Grant No. ERC-2015-COG-681707. This work is part of and supported by the DFG Collaborative Research Centre “SFB 1225 (ISOQUANT).” APPENDIX A: TOWARD SOLVING THE LARGE-NKINETIC EQUATION In this Appendix, we solve some of the integrals appearing in the collision integral of the nonrelativistic large-Nkinetic theory (22) analytically for d¼3spatial dimensions. 1. Collision integral We start by averaging the Boltzmann equation (7) over the solid angle of p. Because of the isotropy of the distribution function, the left-hand side does not change, while the righthand side becomes C½fðt; pÞ¼ZdΩp 4πC½fðt; pÞ ¼π ð2πÞ10 Z∞ 0 dll2dqq2drr2ZdΓCðp; l; q; rÞ ×δðωpþωl−ωq−ωrÞ ×g2 eff ½fðt; ωp−ωq;PÞ ×I2↔2½fðt; p; l; q; rÞ:ðA1Þ We have introduced the momentum difference P≡p−q for the effective coupling g2 eff . In Appendix A2,we will show that it only depends on the magnitude P¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2þq2−2pq cosðθp;qÞ qand hence on the magnitudes of the momenta p¼jpjand q¼jqjand on the angle θp;qbetween them. The three-dimensional momentum integrals in (A1) have been split into radial and solid angle parts according to Rd3q¼R∞ 0dqq2RdΩqwith RdΩq¼R2π 0dφqR1 −1dcosðθqÞ. All angular integrals have been included in ZdΓCðp; l; q; rÞ≡ZdΩpdΩldΩqdΩrð2πÞ3 ×δð3Þðpþl−q−rÞ:ðA2Þ 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 ð2πÞ3δð3Þðpþl−q−rÞ¼Zd3xeiðpþl−q−rÞx:ðA3Þ LARGE-NKINETIC THEORY FOR HIGHLY OCCUPIED …PHYS. REV. D 97, 116011 (2018) 116011-17 Exploiting RdΩpdΩq¼RdΩqdΩp;q, where RdΩp;q denotes the angular integration of paround the axis in q direction, and performing integrations over the angles except for θp;q, we arrive at ZdΓCðp;l;q;rÞ¼24ð2πÞ5Z1 −1 dcosðθp;qÞZ∞ 0 dxx2 ×sinðPxÞ Px sinðlxÞ lx sinðrxÞ rx ;ðA4Þ where we employed R1 −1dyeixy ¼2sinðxÞ=x.Wehave used that the collision integral does not depend on the polar angle θqof q, even after integrating over θp;q, such that the integration R1 −1d cosðθqÞ¼2can be performed in the end, which leads to Eq. (A4). With sinðaÞsinðbÞsinðcÞ¼1 4ð−sinða−b−cÞþsinðaþb−cÞ þsinða−bþcÞ−sinðaþbþcÞÞ ðA5Þ and Z∞ 0 dxsinðaxÞ x¼π 2sgnðaÞ;ðA6Þ we obtain ZdΓCp; l; q; rÞ¼ð2πÞ6 lr Z1 −1 d cosðθp;qÞ P ×ðsgnðPþl−rÞþsgnðP−lþrÞ −sgnðP−l−rÞ−sgnðPþlþrÞÞ: ðA7Þ Taking the integrals over land rof the collision integral (A1) into account, the sign functions in expression (A7) can be conveniently evaluated via Z∞ 0 dlZ∞ 0 drðsgnðPþl−rÞþsgnðP−lþrÞ −sgnðP−l−rÞ−sgnðPþlþrÞÞ ¼2Z∞ P drZPþr 0 dlþZP 0 drZPþr P−r dl−Z∞ 0 dlZ∞ Pþl dr ≡2ZΔðPÞ dldr¼Z∞ P duZP −P dv: ðA8Þ A change of variables from land rto u¼rþland v¼ r−lhas been performed in the last step of Eq. (A8), absorbing a factor of 2. Figure 9visualizes the region of integration ΔðPÞ. A second transformation from d cosðθp;qÞto dPwith d cosðθp;qÞ¼−ðP=pqÞdPleads to C½fðt; pÞ¼ 1 64π3pZ∞ 0 dqZpþq jp−qj dPZ∞ P duZP −P dv ×qðu2−v2Þδðωpþωðu−vÞ=2−ωq−ωðuþvÞ=2Þ ×g2 eff½fðt; ωp−ωq;PÞ ×I2↔2½ft; p; u−v 2;q;uþv 2:ðA9Þ The next step is the evaluation of the energy-conserving delta function. With the quadratic dispersion relation ωp¼p2=2m, the delta function in the collision integral (A9) reads δ1 2mp2þ1 4ðu−vÞ2−q2− 1 4ðuþvÞ2 ¼2mδðp2−q2−uvÞ:ðA10Þ We use the vintegral to evaluate the delta function according to 2mZP −P dvRðvÞδðp2−q2−uvÞ ¼2m uRp2−q2 uΘP−jp2−q2j u;ðA11Þ where Rrepresents every part of the collision integral depending on v. The step function in (A11) can be used to change the integration boundaries of uto R∞ max ðP;jp2−q2j=PÞdu. The resulting form of the collision integral is given by Eq. (25). 2. Retarded self-energy To also simplify the computation of the effective coupling g2 eff of (23), we perform the angular integrations of the one-loop retarded self-energy ΠRin (24) analytically. FIG. 9. The filled region ΔðPÞrepresents the area of integration in (A8). The figure is taken with adaption from Ref. [6]. R. WALZ, K. BOGUSLAVSKI, and J. BERGES PHYS. REV. D 97, 116011 (2018) 116011-18 We first change the integration variable to P−k↦k. Introducing y≡cosðθP;kÞ, where θP;kis the polar angle between the momenta Pand k, one arrives at ΠRðt; ω;PÞ¼lim ϵ→0þ g ð2πÞ2Z∞ 0 dkk2fðt; kÞZ1 −1 dy ×1 ðP2−2PkyÞ=2m−ω−iϵ þ1 ðP2−2PkyÞ=2mþωþiϵ:ðA12Þ Note that, due to the isotropy of fðt; kÞ,ΠRonly depends on the absolute value of the momentum P¼jPj. Moreover, it vanishes for P¼0since 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 ΠRis real for ω¼0. Rewriting (A12) as ΠRðt; ω;PÞ¼lim ˜ ϵ→0þ −mg ð2πÞ2PZ∞ 0 dkkfðt; kÞZ1 −1 dy ×1 y−P2−2mω 2Pk þi˜ ϵþ1 y−P2þ2mω 2Pk −i˜ ϵ; ðA13Þ where we have substituted ˜ ϵ≡mϵ=Pk, enables us to use the principal value integral lim ϵ→0þZ1 −1 1 xþDiϵdx¼PV Z1 −1 1 xþDdx ∓iπZ1 −1 δðxþDÞdx ¼log 1þD −1þD ∓iπΘð1−jDjÞ ðA14Þ for each fraction in Eq. (A13). This leads to ΠRðt; ω;PÞ¼ −mg ð2πÞ2PZ∞ 0 dkkfðt; kÞΓΠðω;P;kÞðA15Þ with the kernel ΓΠðω;P;kÞ¼log ðP2−2PkÞ2−4m2ω2 ðP2þ2PkÞ2−4m2ω2 −iπΘ1−jP2−2mωj 2Pk  −Θ1−jP2þ2mωj 2Pk :ðA16Þ 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 ksing ¼P 2mjωj P;ðA17Þ where the case P¼0is excluded since the whole integrand is zero then. The real part can be rewritten as ReΓΠðω;P;kÞ¼log ðk−P 2Þ2−m2ω2 P2 ðkþP 2Þ2−m2ω2 P2 :ðA18Þ For the imaginary part ImΓΠðω;P;kÞof the integrand (A16), we will assume ω>0since 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 khas increased sufficiently to fulfill the condition k≥jP2−2mωj=2Pof the first Heaviside step function but is still too small to fulfill the condition of the second step function. After exceeding jP2þ2mωj=2P, which makes the second step function one as well, the whole expression vanishes again. Altogether, we find − 1 πImΓΠðω;P;kÞ¼1if jP2−2mωj 2P≤k≤jP2þ2mωj 2P 0else : ðA19Þ Including also the case of negative frequency, we arrive at the final form of the one-loop retarded self-energy in Eq. (26). 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 (25) can be further constrained. Integrals within the real and imaginary parts of the oneloop retarded self-energy (26) are simplified to Re ΠR∶Z∞ 0 dk→ZΛUV ΛIR dk; ðB1Þ Im ΠR∶ZjP2þ2mωj 2P jP2−2mωj 2P dk→Zmin ðjP2þ2mωj 2P;ΛUVÞ max ðjP2−2mωj 2P;ΛIRÞ dk: ðB2Þ The gain-minus-loss part LARGE-NKINETIC THEORY FOR HIGHLY OCCUPIED …PHYS. REV. D 97, 116011 (2018) 116011-19 I2↔2½ft; p; 1 2u−p2−q2 u;q;1 2uþp2−q2 u ðB3Þ of the integrand of the collision integral yields nonvanishing contributions if ΛIR ≤p≤ΛUV; ΛIR ≤q≤ΛUV; ΛIR ≤ 1 2ðup2−q2 uÞ≤ΛUV:ðB4Þ To preserve the gain-minus-loss symmetry, we set the whole I2↔2to 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 pis externally set to the correct momentum range. The qintegration is changed to Z∞ 0 dq→ZΛUV ΛIR dq: ðB5Þ The relation ΛIR ≤ 1 2u2mjωj u≤ΛUV ðB6Þ with 2mjωj≡jp2−q2jis chosen to be fulfilled for plus and minus signs simultaneously. It can be solved for uto obtain new integration boundaries according to the following procedure. In a first step, we consider the left inequality of (B6), i.e., ðu2mjωj=uÞ=2≥ΛIR. Using that u>0, the inequality can be transformed to ðu−ΛIRÞ2≥Λ2 IR ∓2mjωj;ðB7Þ where the two cases of minus and plus signs on the righthand side are distinguished. For the case of the minus sign, no restriction for uis obtained if Λ2 IR ≤2mjωj. Otherwise, one obtains u≤ΛIR −ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 IR −2mjωj qif u<ΛIR;ðB8Þ u≥ΛIR þffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 IR −2mjωj qif u>ΛIR:ðB9Þ In case of the plus sign on the right-hand side of (B7), one finds the constraints u≥ΛIR þffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 IR þ2mjωj qif u>ΛIR;ðB10Þ u≤ΛIR −ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 IR þ2mjωj qif u<ΛIR:ðB11Þ The latter condition (B11) excludes u<ΛIR since ualways has to be positive. In a second step, we consider the inequality on the righthand side of (B6), which can be transformed to ju−ΛUVj≤ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV ∓2mjωj q:ðB12Þ This leads to the four conditions u≥ΛUV −ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV −2mjωj qif u<ΛUV;ðB13Þ u≤ΛUV þffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV −2mjωj qif u>ΛUV;ðB14Þ u≤ΛUV þffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV þ2mjωj qif u>ΛUV;ðB15Þ u≥ΛUV −ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV þ2mjωj qif u<ΛUV;ðB16Þ where the condition (B16) is no constraint, since ΛUV −ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV þ2mjωj p≤0. All considered cases and the successive constraints are employed at the same time. This allows us to change the integration range for uto Z∞ maxðP;2mjωj PÞ du →ZΛUVþffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV−2mjωj p maxðP;2mjωj P;ΛIRþffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 IRþ2mjωj p;ΛUV−ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV−2mjωj pÞ du: ðB17Þ To have a nontrivial integration range for u, its integration boundaries should additionally satisfy ΛUV þffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV −2mjωj q≥max P; 2mjωj P;ðB18Þ while the other two arguments of the maximum function in the lower integration boundary of uare always smaller than the upper integration boundary. Equation (B18) imposes constraints on the Pintegration. Thus, we change the boundaries of the Pintegration according to Zpþq jp−qj dP→Zmin ðpþq;ΛUVþffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV−jp2−q2j pÞ max ðjp−qj;jp2−q2j=ðΛUVþffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Λ2 UV−jp2−q2j pÞÞ dP: ðB19Þ Furthermore, uand Pintegrations are performed only if the lower integration boundaries are smaller than the upper ones. Otherwise, these integrals are set to zero. R. WALZ, K. BOGUSLAVSKI, and J. BERGES PHYS. REV. D 97, 116011 (2018) 116011-20 [1] L. P. Pitaevskii and E. M. Lifshitz, Physical Kinetics (Pergamon, New York, 1981). [2] L. Kadanoff and G. Baym, Quantum Statistical Mechanics (Benjamin, New York, 1962). [3] A. H. Mueller and D. T. Son, Phys. Lett. B 582, 279 (2004). [4] J. Berges and S. Borsanyi, Phys. Rev. D 74, 045022 (2006). [5] P. B. Arnold, Int. J. Mod. Phys. E 16, 2555 (2007). [6] A. P. Orioli, K. Boguslavski, and J. Berges, Phys. Rev. D 92, 025041 (2015). [7] J.-P. Blaizot and E. Iancu, Phys. Rep. 359, 355 (2002). [8] V. E. Zakharov, V. S. Lvov, and G. Falkovich, Kolmogorov Spectra of Turbulence I: Wave Turbulence (SpringerVerlag, Berlin, 1992). [9] R. Micha and I. I. Tkachev, Phys. Rev. D 70, 043538 (2004). [10] S. Nazarenko, Wave Turbulence (Springer-Verlag, Berlin, 2011). [11] P. B. Arnold, G. D. Moore, and L. G. Yaffe, J. High Energy Phys. 01 (2003) 030. [12] A. Kurkela and E. Lu, Phys. Rev. Lett. 113, 182301 (2014). [13] A. Kurkela and Y. Zhu, Phys. Rev. Lett. 115, 182301 (2015). [14] A. Kurkela, Nucl. Phys. A956, 136 (2016). [15] K. Fukushima, Rep. Prog. Phys. 80, 022301 (2017). [16] R. Baier, A. H. Mueller, D. Schiff, and D. T. Son, Phys. Lett. B502, 51 (2001). [17] J. Berges, K. Boguslavski, S. Schlichting, and R. Venugopalan, Phys. Rev. D 92, 096006 (2015). [18] J. Berges, K. Boguslavski, S. Schlichting, and R. Venugopalan, Phys. Rev. D 89, 074011 (2014). [19] J. Berges, K. Boguslavski, S. Schlichting, and R. Venugopalan, Phys. Rev. D 89, 114007 (2014). [20] J. Berges, M. Mace, and S. Schlichting, Phys. Rev. Lett. 118, 192005 (2017). [21] M. Mace, S. Schlichting, and R. Venugopalan, Phys. Rev. D 93, 074036 (2016). [22] J.-P. Blaizot, F. Gelis, J.-F. Liao, L. McLerran, and R. Venugopalan, Nucl. Phys. A873, 68 (2012). [23] A. Kurkela and G. D. Moore, Phys. Rev. D 86, 056008 (2012). [24] J.-P. Blaizot, J. Liao, and L. McLerran, Nucl. Phys. A920, 58 (2013). [25] Z. Xu, K. Zhou, P. Zhuang, and C. Greiner, Phys. Rev. Lett. 114, 182301 (2015). [26] M. C. A. York, A. Kurkela, E. Lu, and G. D. Moore, Phys. Rev. D 89, 074036 (2014). [27] J.-P. Blaizot, J. Liao, and Y. Mehtar-Tani, Nucl. Phys. A961, 37 (2017). [28] K. Zhou, Z. Xu, P. Zhuang, and C. Greiner, Phys. Rev. D 96, 014020 (2017). [29] N. Tanji and R. Venugopalan, Phys. Rev. D 95, 094009 (2017). [30] T. Lappi and J. Peuron, Phys. Rev. D 95, 014025 (2017). [31] A. Kurkela, T. Lappi, and J. Peuron, Eur. Phys. J. C 76, 688 (2016). [32] J. M. Pawlowski and A. Rothkopf, Phys. Lett. B 778, 221 (2018). [33] J. M. Pawlowski and A. Rothkopf, EPJ Web Conf. 175, 07001 (2018). [34] J. Berges, K. Boguslavski, S. Schlichting, and R. Venugopalan, Phys. Rev. Lett. 114, 061601 (2015). [35] J. H. Traschen and R. H. Brandenberger, Phys. Rev. D 42, 2491 (1990). [36] L. Kofman, A. D. Linde, and A. A. Starobinsky, Phys. Rev. Lett. 73, 3195 (1994). [37] J. Berges, A. Rothkopf, and J. Schmidt, Phys. Rev. Lett. 101, 041603 (2008). [38] N. G. Berloff and B. V. Svistunov, Phys. Rev. A 66, 013603 (2002). [39] C. Scheppach, J. Berges, and T. Gasenzer, Phys. Rev. A 81, 033611 (2010). [40] B. Nowak, J. Schole, D. Sexty, and T. Gasenzer, Phys. Rev. A85, 043627 (2012). [41] B. Nowak, S. Erne, M. Karl, J. Schole, D. Sexty, and T. Gasenzer, in Proceedings of the International School on Strongly Interacting Quantum Systems Out of Equilibrium (Oxford University, New York, 2016). arXiv:1302.1448 [42] J. Schmiedmayer and J. Berges, Science 341, 1188 (2013). [43] J. Berges and G. Hoffmeister, Nucl. Phys. B813, 383 (2009). [44] G. D. Moore, Phys. Rev. D 93, 065043 (2016). [45] J. Berges and D. Sexty, Phys. Rev. D 83, 085004 (2011). [46] M. Karl and T. Gasenzer, New J. Phys. 19, 093014 (2017). [47] J. Berges, K. Boguslavski, A. Chatrchyan, and J. Jaeckel, Phys. Rev. D 96, 076020 (2017). [48] Y. B. Zel’dovich and E. V. Levich, Soviet J. Exp. Theor. Phys. 28, 1287 (1969). [49] D. V. Semikoz and I. I. Tkachev, Phys. Rev. Lett. 74, 3093 (1995). [50] D. V. Semikoz and I. I. Tkachev, Phys. Rev. D 55, 489 (1997). [51] J. Berges, Nucl. Phys. A699, 847 (2002). [52] G. Aarts, D. Ahrensmeier, R. Baier, J. Berges, and J. Serreau, Phys. Rev. D 66, 045008 (2002). [53] J. F. Donoghue, Dynamics of the Standard Model (Cambridge University Press, Cambridge, England, 1994). [54] L. Kofman, Lect. Notes Phys. 738, 55 (2008). [55] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation (Clarendon, Oxford, 2003). [56] J. Berges, arXiv:1503.02907. [57] A. Arrizabalaga, J. Smit, and A. Tranberg, Phys. Rev. D 72, 025014 (2005). [58] D. N. Voskresensky, Zh. Eksp. Teor. Fiz. 105, 1473 (1994) [J. Exp. Theor. Phys. 78, 793 (1994)]. [59] D. N. Voskresensky, Yad. Fiz. 59, 2090 (1996) [Phys. At. Nucl. 59, 2015 (1996)]. [60] E. E. Kolomeitsev and D. N. Voskresensky, Nucl. Phys. A973, 89 (2018). [61] I. Chantesana, A. P. Orioli, and T. Gasenzer, arXiv: 1801.09490. [62] T. Gasenzer, J. Berges, M. G. Schmidt, and M. Seco, Phys. Rev. A 72, 063604 (2005). [63] F. Halzen and A. Martin, Quarks & Leptons: An Introductory Course in Modern Particle Physics (Wiley, New York, 1984). [64] A. Schachner, A. P. Orioli, and J. Berges, Phys. Rev. A 95, 053605 (2017). [65] T. Gasenzer, B. Nowak, and D. Sexty, Phys. Lett. B 710, 500 (2012). LARGE-NKINETIC THEORY FOR HIGHLY OCCUPIED …PHYS. REV. D 97, 116011 (2018) 116011-21 [66] S. Mathey, T. Gasenzer, and J. M. Pawlowski, Phys. Rev. A 92, 023635 (2015). [67] M. Karl, B. Nowak, and T. Gasenzer, Phys. Rev. A 88, 063615 (2013). [68] J. Berges and B. Wallisch, Phys. Rev. D 95, 036016 (2017). [69] J. Berges and D. Sexty, Phys. Rev. Lett. 108, 161601 (2012). [70] B. Nowak and T. Gasenzer, New J. Phys. 16, 093052 (2014). [71] M. Galassi et al.,GNU Scientific Library Reference Manual, 3rd ed. (2009), ISBN: 0954612078. [72] S. Tsutsui, J.-P. Blaizot, and Y. Hatta, Phys. Rev. D 96, 036004 (2017). R. WALZ, K. BOGUSLAVSKI, and J. BERGES PHYS. REV. D 97, 116011 (2018) 116011-22