Attractive vs. repulsive interactions in the Bose-Einstein condensation dynamics of relativistic field theories

We study the impact of attractive self-interactions on the nonequilibrium dynamics of relativistic quantum fields with large occupancies at low momenta. Our primary focus is on Bose-Einstein condensation and nonthermal fixed points in such systems. As a model system we consider O(N)-symmetric scalar field theories. We use classical-statistical real-time simulations, as well as a systematic 1/N expansion of the quantum (2PI) effective action to next-to-leading order. When the mean self-interactions are repulsive, condensation occurs as a consequence of a universal inverse particle cascade to the zero-momentum mode with self-similar scaling behavior. For attractive mean self-interactions the inverse cascade is absent and the particle annihilation rate is enhanced compared to the repulsive case, which counteracts the formation of coherent field configurations. For N>= 2, the presence of a nonvanishing conserved charge can suppress number changing processes and lead to the formation of stable localized charge clumps, i.e. Q-balls.


I. INTRODUCTION
Nonequilibrium quantum fields with large occupation numbers at low momenta are frequently encountered in the context of cosmology. Important examples include the decay of coherent oscillations of the inflaton during the reheating stage after inflation [1] or the production of dark matter axions from the misalignment mechanism and from the decay of axionic strings and domain walls [2][3][4][5][6]. Therefore, understanding the dynamics of such systems is important. Particularly interesting is the possibility of Bose-Einstein condensation, which is the macroscopic occupation of the quantum state with the lowest energy, in such systems. In the case of axions, the resulting collective quantum behavior may have important observational consequences in cosmology (see, e.g., Refs. [7][8][9]) and would leave distinct imprints in direct detection experiments [10][11][12]. While many arguments have focused on the formation rate of a condensate [13][14][15][16][17][18][19], a qualitatively even more important point is the attractive nature of the relevant interactions [20,21] that tends to favor localized structures instead of a spatially constant condensate [21]. Inspired by this, it is one of our main aims to study the impact of attractive interactions and to delineate the differences to the repulsive case.
In recent years, there has been a significant advance in the theoretical understanding of the dynamics of isolated highly occupied quantum fields. Many characteristic properties of the dynamics in such extreme conditions turn out to be insensitive to details of the underlying model and initial conditions. This allows one to classify theories with different microscopic descriptions into universality classes [22,23]. This notion of universality is based on the existence of nonthermal fixed points [24,25] that are nonequilibrium attractor solutions with self-similar scaling behavior of correlation functions and of the particle momentum distribution. Self-similarity in this case is associated to the transport of some conserved quantity in momentum space [26,27]. Self-similar scaling regions can represent transport of energy toward high momenta or particle number transport toward the zero mode. While the first case drives the thermalization process by pushing the typical hard scale to larger momenta [27], the second one may lead to the formation of a Bose-Einstein condensate out of equilibrium [28]. Both cascades have been observed in different regions of the same momentum distribution function in scalar models with positive quartic self-coupling [23,29].
In general, for scalar N-component fields φ a , a ¼ 1;…;N, with self-interactions of the form ∼λ 2n ðφ a φ a Þ n , positive or negative signs of the couplings determine whether they are repulsive or attractive. While repulsive interactions have the tendency to dilute concentrations of energy density over position space, attractive interactions encourage such local concentrations. It is thus not surprising that condensation dynamics is strongly affected by the presence of attractive interactions, since the lowest-energy configurations in this case are localized "clumps." Therefore, how characteristic features of the dynamics depend on the type of selfinteractions is an interesting question.
In the present work, we investigate the influence of attractive self-interactions on the far-from-equilibrium dynamics of scalar fields and, in particular, their impact on Bose-Einstein condensation. For our model, we consider OðNÞ-symmetric relativistic field theory 1 in 3 þ 1 dimensions, exhibiting an attractive quartic self-interaction (λ < 0) and stabilized by a repulsive sextic term so that the potential has an OðNÞ-symmetric global minimum. We are interested in the nonperturbative regime of large occupancies f ∼ 1=jλj for typical momenta, in the weak coupling regime with jλj ≪ 1. Our analysis is based on real-time classical-statistical lattice simulations, using the possibility of mapping the dynamics of highly occupied quantum fields onto a classical-statistical field theory evolution 2 [33,34]. In addition to this, we develop a vertex-resummed kinetic theory, based on the 1=N expansion of the quantum two-particle-irreducible (2PI) effective action in the presence of quartic and sextic interactions to next-to-leading order [35]. This extends well-established kinetic descriptions [26,27] to the considered highly occupied regime, in analogy to the repulsive φ 4 theory [23,36]. We derive the corresponding effective kinetic equation, discuss its properties, and compare it with our observations from classical-statistical simulations.
The outline of this work is as follows. In Sec. II, we describe our model and the considered initial conditions. In Sec. III, we present the results from the classical-statistical simulations for the single-component model. In Sec. IV, we extend that discussion to the case of multicomponent fields and derive the vertex-resummed kinetic description. We conclude in Sec. V. Throughout this paper, we work in Minkowski spacetime with metric η μν ¼ ð1; −1; −1; −1Þ and set ℏ ¼ k B ¼ c ¼ 1.

II. MODEL
In this paper, we study OðNÞ-symmetric scalar field theories in flat, 3 þ 1-dimensional spacetime. We consider potentials that contain an attractive self-interaction but have a global minimum at φ ¼ 0. The simplest form of a potential capturing these properties is given by φ a ðxÞφ a ðxÞ þ λ 4!N ðφ a ðxÞφ a ðxÞÞ 2 þ λ 6 6!N 2 ðφ a ðxÞφ a ðxÞÞ 3 : ð2:1Þ Here, summation over repeated indices is implied, and the factors of 1=N are chosen such that the classical action UðφÞ scales proportionally to N. The parameters m 2 and λ 6 are positive. For completeness, we consider both cases λ > 0 and λ < 0 and refer to those models throughout this paper as repulsive and attractive, respectively. We assume weak couplings, in particular jλj ≪ 1.
After the rescaling prescriptions x ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi jλj m 2 φ a φ a N q and uðxÞ ¼ jλj Nm 4 UðφÞ, the potential can be written as þ sgnðλÞ x 4 4! þ g 2 x 6 6! : ð2:2Þ The dimensionless parameter g 2 , defined as characterizes the strength of the sextic coupling. In the case λ < 0, the potential falls below quadratic near its minimum but stays positive everywhere if g 2 is sufficiently large (g 2 > 5=8). Similarly, it becomes convex for g 2 > 3=2. This is demonstrated in Fig. 1, where the rescaled potential is shown for several values of g 2 , together with the free theory parabola uðxÞ ¼ x 2 =2, which may be viewed as a dividing line between repulsion and attraction. We consider spatially homogeneous and isotropic quantum systems and employ Gaussian initial conditions, which can be formulated in terms of the macroscopic field 3 ϕ a ðtÞ ¼ hφ a ðt; xÞi and the single-particle distribution function fðt; jpjÞ at the initial time t ¼ 0 [37]. To observe the buildup of correlations and condensation, we start with a  vanishing macroscopic field and large initial occupation numbers f ≳ 1=jλj of the following form: These initial conditions correspond to a "box" in momentum space up to the momentum scale p 0 , with the amplitude A 0 ≳ 1. Another type of initial conditions that is also often employed in the literature corresponds to a large macroscopic field ϕð0Þ ∼ 1= ffiffiffiffiffi jλj p in the absence of initial occupation numbers with fð0; jpjÞ ¼ 0 [27]. However, in a transient regime of parametric resonance, the evolution of the field ϕ triggers growth of unstable modes, and the system eventually becomes highly occupied, recovering the fluctuation initial conditions [27,29]. We will thus only employ box-type initial conditions (2.4).
Despite starting with Gaussian initial conditions, selfinteractions render the dynamics of the system non-Gaussian, and higher-order nonfactorizable correlation functions build up for t > 0. Although the notion of particles is not uniquely defined in interacting relativistic field theory, the following definition of fðt; jpjÞ turns out to provide a useful quasiparticle interpretation for scalar systems [24,37]. With the help of the anticommutator expectation value hfφ a ðt; xÞ;φ a ðt 0 ; x 0 Þgi; ð2:5Þ We start our discussion with the simplest case of a singlecomponent field theory. Since we consider systems with large occupation numbers at characteristic momenta, perturbative approaches based on a coupling expansion are inapplicable. On the other hand, the large occupation numbers of typical modes that dominate the energy density strongly exceed unity, such that genuine quantum effects become subdominant and the dynamics of the system becomes essentially classical. More precisely, in this regime, quantum field theory can be mapped onto a classicalstatistical field theory [33,34], which can be simulated numerically. Otherwise, as soon as typical occupation numbers become of order unity, f ≲ 1, the mapping becomes inaccurate [32,[38][39][40], and the system leaves the classical regime, thermalizing eventually [41]. Thus, we restrict ourselves to the case of sufficiently large occupation numbers and study the scalar systems with classicalstatistical simulations in this section.
In such lattice simulations, one samples over the initial conditions and evolves each realization according to the classical field equation of motion. Correlation functions are obtained by taking ensemble averages over the classical trajectories, where quantum observables are translated to their classical counterparts. For instance, the anticommutator expectation value h 1 2 f:; :gi in (2.5) is replaced by the averaged product of fields in the classical-statistical framework.
Numerical simulations are performed on a threedimensional (3D) cubic lattice with periodic boundary conditions with up to 1024 3 points and lattice spacings in the range of 0.0625-0.5 p −1 0 , with p 0 introduced by the initial conditions (2.4). 6 The classical equations of motion are solved using the standard leapfrog algorithm, and the details are presented in Appendix A. The distribution function fðt; pÞ and the dispersion relation ωðt; pÞ are computed by use of Eqs. (2.6) and (2.7), respectively. The latter is also used to calculate the effective mass MðtÞ by fitting ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi p 2 þ M 2 ðtÞ p . Another important quantity is the total particle number density 7 If not stated otherwise, we employ g 2 ¼ 2 and m ¼ 0.7 p 0 . More generally, all dimensionful quantities shown in the 4 Because of the system's homogeneity, the Fourier transformation is performed with respect to the relative coordinates: Fðt; t 0 ; pÞ ¼ R d 3 Δx expð−ipΔxÞFðt; t 0 ; ΔxÞ. 5 A closely related alternative definition of the distribution function, also employed in the literature [23,28], is ðfðt; jpjÞþ 1 2 þ ð2πÞ 3 δ ð3Þ ðpÞn cond ðtÞÞ=ωðt; pÞ ¼ Fðt; t; pÞ. It turns out to yield the same distribution function, if a relativistic form for ω is employed with a suitably chosen effective mass. 6 In the central and lower plots in the left panel of Fig. 2, we combined two data sets with different lattice spacings (0.0625 p −1 0 with 0.25 p −1 0 and 0.1 p −1 0 with 0.25 p −1 0 , respectively) to show a larger momentum range of the distribution function. In the overlapping momentum region, the two simulations agree to very good accuracy. 7 On the lattice, integrals are translated to sums over momenta according to plots of this section are made dimensionless by rescaling with appropriate powers of p 0 . In addition, we plot rescaled functions F ↦ jλjF, f ↦ jλjf, and n ↦ jλjn, since for fixed m 2 , A 0 and g 2 , the classical-statistical dynamics of these combinations does not depend on the absolute value of λ. However, the dynamics crucially depends on its sign. Both cases of λ > 0 and λ < 0 are considered separately.
For a better understanding of the role of the negative coupling, we first briefly discuss the dynamics for positive λ, which corresponds to repulsive self-interactions. The same model without the sextic coupling has been extensively studied previously [23,27,28]. For the range of initial occupancies we consider, we find that the presence of the sextic coupling does not introduce any significant difference for the dynamics.
For the numerical simulations, we employ box initial conditions with A 0 in the range 20-500. In the beginning of the time evolution, the momentum distribution transforms from the initial box into a smoother function, exhibiting a dual cascade. The direct UV cascade populates the initially empty high-momentum modes. At the same time, the inverse cascade at low momenta constantly transports particles to the zero-momentum mode, leading to farfrom-equilibrium Bose-Einstein condensation [28]. After some transient time, the dynamics of each of the cascades becomes self-similar, fðt; jpjÞ ¼ t α f S ðt β jpjÞ; ð3:2Þ with separate sets of real scaling exponents α and β and fixed-point distributions f S . This self-similar dynamics reflects the system being in the vicinity of a nonthermal fixed point [24,25].
In the upper left panel of Fig. 2, we show the rescaled distribution function t −α fðt; jpjÞ as a function of the rescaled momentum t β jpj at different times for the values α ¼ 3=2 and β ¼ 1=2, while the inset gives the same distributions without rescaling. Since the rescaled curves lie on top of each other, the system follows a self-similar evolution with the given values for the exponents α and β.
The fixed point distribution f S exhibits an approximate power-law behavior ∼jpj −κ , with κ ≈ 4.5. The same values for these three exponents have been obtained also for the φ 4 model [23].
The occupancy of the zero-momentum mode grows during the self-similar evolution approximately as a power-law ∼t α . This can be directly understood from (3.2), by setting p ¼ 0. For a finite volume, this growth continues until Fðt; t; p ¼ 0Þ becomes proportional to the volume V, which signals the emergence of the condensate according to Ref. [28] (see also footnote 5) with the replacement ð2πÞ 3 δð0Þ → V. This is demonstrated in the upper right panel of Fig. 2, where the quantity Fðt; t; p ¼ 0Þ=V is shown as a function of time for several volumes. All of this has also been observed in the φ 4 model [23]. There, it was also shown that this dynamics implies that the time required for creating a condensate scales with volume.
The inverse cascade represents the transport of the conserved particle number to the zero mode. Particle number conservation follows from the relation α ¼ 3β, which was discussed above: Moreover, the condensate stays approximately constant after its creation, as can be seen in the upper right panel of Fig. 2. Hence, also, the total particle number density nðtÞ is approximately conserved, since its dominant contribution comes from the infrared modes and the condensate. Therefore, at this stage of the evolution, elastic scatterings are the dominant processes for the dynamics in repulsive scalar systems and lead to the emergence of a long-lived condensate. 8 At high momenta, another self-similar scaling region emerges after some time. There, we have extracted scaling exponents β 0 ≈ −1=4.2 and α 0 ≈ 4β 0 , with error estimates of the order of 10%. This ratio of the exponents corresponds to a transport of conserved energy toward high momenta in the ultrarelativistic regime with ω p ≈ jpj. We have also verified that these exponents are unaffected by the presence of the sextic coupling. 9 The scaling exponents of both cascades thus seem to be independent of the microscopic details of the model. In Sec. IV, we explain the reason for this from the point of view of the 2PI 1=N expansion to next-to-leading order.

B. Attractive (λ < 0)
We now consider the attractive model with λ < 0. In addition, we take the positive sextic coupling characterized by g 2 sufficiently large such that the potential has its global minimum at the origin. Results are presented for g 2 ¼ 2. 8 The decay of the condensate and the infrared modes has been studied in Ref. [42]. There, it was shown that with an increasing mass parameter m such decays become more suppressed. For the values of the mass parameter and the initial occupancy considered here, we confirm that the total particle number density is practically conserved. 9 The scaling exponents of the self-similar energy cascade have been studied in Ref. [27], and for k ↔ l scattering processes, they are given by β 0 ¼ −1=ð2ðk þ lÞ − 1Þ. Our extracted values slightly exceed β 0 ¼ −1=5, which would correspond to 2 ↔ 1 þ soft scatterings, where the soft particle corresponds to the involvement of the condensate [27]. The strength of the fluctuations, characterized by the parameter A 0 , determines which of the two self-interactions is more important for the dynamics. Therefore, we split our discussion into two parts. 10 We first consider the case of relatively weak initial fluctuations, such that the negative, attractive quartic coupling is more important and mean interactions are attractive. In the second regime, initial fluctuations are strong, and the positive repulsive sextic coupling dominates, turning the mean interactions repulsive. In systems with high occupation numbers, we will even encounter a transition between these regimes during the time evolution.

Attractive mean interactions
The inverse particle cascade is absent in this case. This is demonstrated for A 0 ¼ 60 in the central left panel of Fig. 2, which contains several snapshots of the distribution function. While a transient growth of the low-momentum occupancies is observed at early times, it does not lead to an enhancement of infrared modes as strong as in the repulsive models, ∼jpj −4.5 , that were discussed above. For comparison, a power law ∼jpj −2 is additionally shown in the plot.
Because of the absence of an inverse particle cascade, the zero-momentum mode does not develop a condensate part. This is demonstrated in the central right panel of Fig. 2, where the time evolution of the correlation function Fðt; t; p ¼ 0Þ divided by the volume V is shown for several volumes. As was mentioned in the previous section, for a finite volume, the presence of a condensate would lead to a contribution in Fðt; t; p ¼ 0Þ, which scales proportional to V. Such a contribution is clearly absent in the attractive case, which becomes more stringent when compared to the repulsive system depicted in the upper right panel of Fig. 2. In other words, no long-range order is being established. At hard momenta, a self-similar cascade develops after some transient time. The extracted values of the scaling exponents β 0 ≈ −1=5 and α 0 ≈ 4β 0 are slightly smaller than those observed in the repulsive model. 11 As in the repulsive theory, the ratio α 0 =β 0 ¼ 4 is characteristic of a transport of conserved energy density toward high momenta. Hence, the dynamics at high momenta can be also interpreted as a direct energy cascade.
As compared to repulsive systems, number-changing processes appear to be more efficient in the attractive model. The time evolution of particle number density nðtÞ is shown in the upper left panel of Fig. 3, where the lines with A 0 ¼ 12 and A 0 ¼ 60 correspond to the considered regime with mean attractive interactions. After some transient time, nðtÞ decreases approximately as a power law nðtÞ ∼ t β 0 . This exponent can be understood by plugging the self-similar ansatz into Eq. (3.3) in the absence of a highly occupied infrared region. This loss of particles is quite different from the repulsive model discussed above. The repulsive model features an enhanced infrared region, and its total particle number density is approximately conserved, despite the presence of a direct energy cascade with α 0 ¼ 4β 0 at hard momenta.
The absence of long-range order in our considered system can be viewed as a manifestation of the tendency of attractive self-interactions to fragment the field into localized configurations [21]. Importantly, number-changing processes prevent the formation of visible spatial structures. 12

Repulsive mean interactions
If the initial occupancy A 0 is sufficiently increased, with other parameters fixed, the repulsive sextic term becomes more important initially. In this case, an early stage is observed, during which the distribution function exhibits a dual cascade. This is shown in the lower left panel of Fig. 2, where we show several snapshots of the distribution function at different times. At early times, a self-similar inverse particle cascade is observed with the same scaling exponents α ≈ 3=2 and β ≈ 1=2 and an approximate powerlaw behavior ∼jpj −4.5 as in the repulsive λ > 0 case. However, at later times, the dynamics starts to deviate from a self-similar evolution, and eventually the infrared enhancement decays.
A similar evolution is observed for the zero-momentum mode. At early times, condensation occurs as in the repulsive model (especially for sufficiently small volumes, where the condensation time is small enough), as is visible in the lower right panel of Fig. 2. However, at later times, the condensate collapses.
The decay of the high infrared occupancy and the collapse of the condensate are caused by inelastic processes. Indeed, as can be seen in the upper left panel of Fig. 3 10 Qualitatively, the two regimes can be easily understood from Fig. 1, where the blue potential curve lies below the free theory curve for small field values and above it for large field values. We will quantify this for the case of an N component theory in Sec. IV, Eq. (4.13). 11 In the absence of a condensate, one would expect elastic 2 ↔ 2 scatterings to be the dominant processes for the direct cascade at high momenta, which would lead to β 0 ¼ −1=7 [27]. As in the repulsive model, we observe some discrepancy from the expected value of β 0 . Nevertheless, the hierarchy of the two β 0 exponents, that jβ 0 j for the repulsive model exceeds its value in the attractive system, qualitatively coincides with expectations. 12 Although in single simulations relatively long-lived configurations with energy density concentration may be observed, they do not affect the distribution function in any significant way. We have verified this by removing them by explicitly replacing the field value in these regions by a small fixed field value and comparing the distribution function afterward. We note that the importance of such configurations may depend on the chosen parameter region; see, for example, Ref. [43].
for A 0 ¼ 220 and A 0 ¼ 500, the total particle number density decreases with time. This interpretation is further strengthened by the observation of accompanying peaks of the distribution function at momenta jpj ≳ MðtÞ shown in the lower left panel of Fig. 2, where MðtÞ is the effective mass. The locations of these peaks correspond to 2n → 2 inelastic scatterings off the "condensate" or off modes with momenta considerably below MðtÞ. Such decay peaks have also been observed in quartic models [32,42]. For instance, the sharp peak of the t ¼ 8000 curve is located at momentum jpj ≈ ffiffi ffi 3 p M ≈ 1, 13 which corresponds to an inelastic scattering of four particles with jpj ≪ M into two particles with energy ϵ ≈ 2M. The peaks move toward smaller momenta with time. This is a consequence of the particle loss, which leads to the decrease of the effective mass (see the lower right panel of Fig. 3). Moreover, the peaks become more pronounced as the height of the initial box is increased, since in this case the condensate and low-momentum modes contain more particles that can undergo inelastic scatterings.
Importantly, the particle loss generates a transition from repulsive to attractive mean interactions. An important quantity that characterizes the strength of the fluctuations, and visualizes the transition, is the integrated correlator  13 The value of the effective mass at t ¼ 8000 is M ≈ 0.6, which can be deduced from the lower right panel of Fig. 3, using the mass parameter m ¼ 0.7. This correlation function characterizes the strength of the fluctuations and determines shifts to the effective mass and the coupling. As can be seen from its definition (3.4), it is an integral over all momentum modes and is additionally dominated by the large infrared modes, which can be understood from footnote 5. Its decrease is a consequence of particle loss. The time evolution of this quantity is shown in the upper right panel of Fig. 3 for different initial occupancies. As can be seen, for large values of A 0 (220 and 500 in the plot), there is a stage when the correlator decreases rapidly. This stage represents the transition between the two regimes and corresponds to the rapid decrease of the total particle number density, the collapse of the condensate, and the effective mass MðtÞ falling below its bare value m, as seen in the residual panels of the figure. Around this time window, the observed decay peaks in the particle spectrum that were discussed above become more pronounced.
In addition, we show the zero crossing of a modified quartic coupling, as a gray dashed line in the upper right panel of Fig. 3. This time-dependent quartic coupling results from a systematic 1=N expansion to next-to-leading order and will be derived in Sec. IV B. There, we will show that to this order the coupling parameter λ should essentially be replaced byλðtÞ in scattering processes of scalar theories with a large N.
Interestingly, this modified coupling changes its sign during the transition, from positive to negative, as can be seen in the figure. This is a remarkable observation given the fact that the modified coupling is derived in a 1=N expansion while we consider a single-component model here. Similar observations will be made in Sec. IVA for scalar theories with N ≥ 2 field components. In summary, although the dynamics of the infrared fixed point is insensitive to the precise form of the self-interactions, it does depend on whether mean self-interactions are repulsive or attractive. While repulsive self-interactions have the effect of moving the system closer to the infrared nonthermal fixed point, which also leads to condensation, the presence of attractive self-interactions has the opposite effect. Such interactions are accompanied by enhanced annihilation processes, which prevent the formation of large structures. In contrast, an energy cascade to higher momenta occurs in both cases and is, in this sense, a more generic phenomenon.

IV. DYNAMICS FOR MULTICOMPONENT FIELDS (N > 1)
In this section, we study the dynamics of our model in the case in which the number of components is larger than 1. Our analysis includes two different nonperturbative methods. One of them is classical-statistical simulations, similar to those that were employed in the previous section. Many aspects of the dynamics turn out to be universal and do not indicate strong dependence on N. In addition to this, we develop a vertex-resummed kinetic theory for our model in Sec. IV B, based on a 1=N expansion of the 2PI quantum effective action up to the next-to-leading order (NLO). The 1=N expansion relies on a large number of field components and is based on a classification of the contributions to the effective action according to their scaling with N. This provides a controlled expansion parameter, which is not restricted to weakly coupled or low occupied systems. The 2PI 1=N expanded theory to NLO and the vertexresummed kinetic theory have been successfully applied to repulsive quartic OðNÞ-symmetric scalar field theories and were able to describe observed phenomena [35,37,44]. Therefore, these approximations also seem to be a suitable approach to the case of attractive interactions.

A. Classical-statistical simulations
We extend here the analysis of Sec. III to the N > 1 case. More specifically, we consider two-, four-, and eight-component fields and use lattices with up to 512 3 points and lattice spacings in the range of 0.1-0.5 p −1 0 . We checked that the shown results are insensitive to the lattice spacing and grid size. The details of the numerical calculations can be found in Appendix A. We employ g 2 ¼ 2 and m ¼ 0.7 p 0 , with p 0 introduced by the initial conditions (2.4), and rescale all dimensionful quantities in the figures with appropriate powers of p 0 . The repulsive and attractive models are discussed separately.

Repulsive (λ > 0)
As in the quartic theory [23], the basic features of the time evolution of the momentum distribution function fðt; jpjÞ reveal weak dependence on the number of field components. Similar to the N ¼ 1 case, two cascades emerge after some transient time: an inverse particle cascade at low momenta and a direct energy cascade at high momenta. The inverse cascade becomes self-similar, with the same scaling exponents α ≈ 3=2 and β ≈ 1=2 and the same scaling function f S ðpÞ that were observed in the single-component model, and leads to the formation of a condensate. The dynamics at low momenta is dominated by elastic scatterings, so the total particle number (3.1) is approximately conserved.
There are some differences in the spatial structure of the condensate for different N, which, however, do not affect the above-mentioned universal self-similar dynamics. They can be studied directly by looking at the classical field φ a ðt; xÞ in position space of a single simulation, before taking the ensemble average. This has been done in Ref. [42] for the repulsive φ 4 theory. In the singlecomponent model, the condensate in one simulation represents a coherently "back-and-forth" oscillating classical field. For N ≥ 2, more general orbital forms in field component space for coherent oscillations are possible. Circular orbits turn out to be energetically preferred [42], and as a result, in multicomponent models, the coherent oscillations are circular.
This may be further understood by the presence of continuous OðNÞ symmetries for N ≥ 2 that impose additional restrictions on the field [42]. Importantly, there are NðN − 1Þ=2 conserved Noether charges of the form associated with global rotations in the corresponding planes in field space. For our considered initial conditions, each of these charges vanishes (see Appendix A). However, the charges may differ from zero locally. As a result, for N ¼ 2, the condensate splits into different domains, each of which has a positive or negative charge [see Fig. 4(a)] depending on the direction of the circular oscillations, such that the total charge vanishes. These domains are separated from each other by thin walls, also visible in the figure, where the field oscillates linearly and carries no charge. While numberchanging processes within each domain are suppressed because of charge conservation, the presence of domain walls leads to a small leakage of particle number. 14 For N > 2, there is a continuous group of possible oscillations, and domain walls are absent there. Therefore, annihilation processes are even stronger suppressed. We refer to Ref. [42] for a detailed discussion of topological defects for different N.

Attractive mean interactions with λ < 0
When the quartic coupling is negative, two qualitatively different regimes of the dynamics can be observed, depending on which of the two self-interactions is more important.
For relatively weak fluctuations, mean interactions are attractive, while for strong fluctuations, they are repulsive. The dynamics in both regimes proceeds along the same lines as for N ¼ 1.
In the case of mean attraction, the inverse cascade, which was observed in the repulsive model, is absent, and no condensate is formed. The low-momentum modes grow for some transient time and decay slowly afterward. Also, the total particle number density soon decreases approximately as a power-law ∼t β 0 with β 0 ≈ −1=5, being the same as discussed in Sec. III for N ¼ 1. Hence, number-changing processes are important also for N > 1 with mean attractive interactions.
Although in single simulations small clumps, which typically contain positive and negative charge, can develop during the evolution, they quickly decay via particle annihilation and thus are not stable. A typical late-time configuration is shown in Fig. 4(c) for A 0 ¼ 60 at t ¼ 5000. This is a notable difference from the nonrelativistic case, which features a conserved charge (particle number) and accordingly leads to the formation of persistent clumps (see, e.g., Ref. [45]). 15 As already mentioned, our considered initial conditions generate no Oð2Þ charge. However, with a slight modification of the initial conditions, which on the level of a single simulation consists in correlating the initial oscillation phases of each of the field components, a nonvanishing charge Q can be generated. This is explained in Appendix A, where it is also shown that the value of the charge can be controlled by the height A 0 of the initial distribution. The charge density of a typical late-time configuration with these modified initial conditions is shown in Fig. 4(b) for the case of A 0 ¼ 60 at time t ¼ 5000. As can be seen, most of the system's charge is concentrated in an almost spherical object, FIG. 4. Snapshots of the charge density ρðxÞ ¼ φ 1 ðxÞ _ φ 2 ðxÞ − φ 2 ðxÞ _ φ 1 ðxÞ at t ¼ 5000 for two-component field theory with A 0 ¼ 60 and (a) λ > 0 without charge asymmetry, (b) λ < 0 with charge asymmetry and (c) λ < 0 without charge asymmetry. 14 Similar to the single-component case, a nonzero bare mass further suppresses particle number-changing processes [42]. 15 For the case of axions, long-range gravitational interactions are supposed to be the dominant attractive interaction. In this case, it is not a priori clear if the number-changing processes are similarly effective or if there exists an approximate particle number conservation. which corresponds to a so-called Q ball [46]. 16 This is in stark contrast to our usual choice of initial conditions (2.4) without an initial charge and more similar to the nonrelativistic model with a conserved particle number.
Indeed, it is a well-known feature of potentials, which contain an attractive self-interaction and thus grow slower as compared to quadratic in some range of field values, that the presence of a Noether charge associated with an unbroken global symmetry leads to the existence of compact lowest-energy field configurations, i.e., Q balls [46]. They are spherically symmetric in coordinate space, and at each spatial point, the classical field oscillates along circular orbits in field space. Q balls are nontopological solitons since their stability against decay is guaranteed by the conserved charge they carry. For the simplest case of N ¼ 2, a Q ball has the form (in polar coordinates with an appropriate choice of the origin) where ϕðrÞ determines the radial profile. Q balls with a large charge can be described using the thin-wall approximation [46], and their profile function has a form close to a step function, ϕðrÞ ¼ ϕ Q ΘðR − rÞ, with a constant amplitude ϕ Q and with R denoting the radius. In contrast, Q balls with small charge have thick walls, and there is evidence that they become unstable for a sufficiently small charge [47].

Repulsive mean interactions with λ < 0
For sufficiently high initial occupancies, mean interactions become initially repulsive. The dynamics in this regime is very similar to the λ > 0 case. In particular, an inverse particle cascade with the same values of the scaling exponents α ≈ 3=2 and β ≈ 1=2 is again observed and leads to the formation of a condensate. We have seen in Sec. III B that in the single-component model particle number was not conserved and the particle loss gradually increased the role of the (attractive) quartic self-interaction. This eventually generated a transition to the mean attractive regime, accompanied by a fragmentation of the condensate. We observe a similar behavior for N ¼ 2. However, the annihilation rate is suppressed as compared to N ¼ 1, and as a result, the lifetime of the condensate is longer. The structure of the condensate plays an important role for this suppression [42]. As discussed in Sec. IVA 1, for N ¼ 2, the classical field locally carries maximal charge, which means that particle annihilation from the same domain of the condensate would violate charge conservation. The domain walls, that are necessarily present in this case, are therefore expected to be the dominant source of particle evaporation. For larger values of N, such as N ¼ 4 or N ¼ 8, number-changing processes in the repulsive regime are even stronger suppressed, which makes the condensate long lived.
To demonstrate the transition, we plot the integrated correlation function Fðt; t; Δx ¼ 0Þ (3.4) as a function of time in Fig. 5 for several values of the initial occupancy A 0 . The two panels correspond to N ¼ 2 (left) and N ¼ 8 (right). As in the N ¼ 1 case (see Fig. 3), the region where Fðt; t; Δx ¼ 0Þ decreases rapidly corresponds to the collapse of the condensate and approximately separates the that is derived in Sec. IV B and shown in the right y axis of the right panel. 16 After some time several smaller Q-balls add up to a single large Q-ball, that is visible in the figure. mean repulsive from the mean attractive regime. This is also demonstrated by a black dashed line that shows the zero crossing of the modified couplingλðtÞ of Eq. (3.5). The modified coupling indeed changes its sign during the fast transition, similarly to what was observed for the N ¼ 1 case. Since the derivation ofλðtÞ is based on a 1=N expansion, as will be shown in the following section, the modified coupling is expected to provide a good description for theories with many field components. That its zero crossing mostly falls into the transition region also for theories with few field components is an unanticipated but important observation that supports the interpretation in terms of repulsive and attractive mean interactions.
In addition, one observes another interesting phenomenon in Fig. 5. For N ¼ 2, the transition occurs for any of the considered amplitudes of the initial occupancy within the simulation time, occurring later for higher amplitudes. Compared to N ¼ 2 for N ¼ 8, it generally occurs at later times. Remarkably, there seems to be a rather dramatic increase between A 0 ¼ 190 (blue) and A 0 ¼ 220 (green), such that the transition does not occur within the simulation time. The system remains in the mean repulsive regime, and the modified coupling stays positive,λðtÞ > 0. This difference between N ¼ 2 and N > 2 theories is consistent with the argument about enhanced suppression of annihilation processes for N > 2 as compared to N ¼ 2 theories in the end of Sec. IVA 1.
Finally, if the initial conditions are modified to generate a nonvanishing charge density, as explained in Sec. IVA 2 and Appendix A, a condensate is formed for sufficiently high initial occupancies. It consists of a single charge domain, and number-changing processes are forbidden. As a result, the condensate is stable, and no transition to the mean attractive regime occurs in this case. This is consistent with the picture that for very high charge densities the condensate corresponds to a circular orbit in field space with an amplitude such that the field is always in the repulsive regime according to Fig. 1.

B. Vertex-resummed kinetic theory from 2PI 1=N expansion
To understand the above observations from classicalstatistical simulations, we derive in the following an effective kinetic description that is based on a 1=N expansion to NLO. As it turns out, the description is similar to the vertex-resummed kinetic theory that has been used to describe scaling phenomena in repulsive quartic scalar theories [23,36]. The main difference is the emergence of a modified couplingλðtÞ, which incorporates corrections from the presence of the sextic coupling. As mentioned in previous sections, the modified coupling varies with time and may even change its sign, which is then accompanied by a transition from repulsive to attractive mean interactions.

2PI 1=N expansion to NLO
We start with the full quantum field theory the properties of which are encoded in the 2PI effective action Γ½ϕ; G [48]. It is a free energy functional, parametrized in terms of the macroscopic field ϕ a ðxÞ and the time-ordered connected two-point function G ab ðx; yÞ, also referred to as the (full) propagator and given by The time variables are defined on the Schwinger-Keldysh contour C [49], which starts at t ¼ 0, runs forward along the real-time axis, and then returns back to t ¼ 0.
Here, T C denotes the time-ordering operator along the time contour C. The correlation functions ϕ and G satisfy the quantum equations of motion that are given by stationarity conditions of the effective action, and sgn C ðx 0 − y 0 Þ is AE1 depending on whether x 0 is after or before y 0 along the closed time contour. F is the statistical propagator, and ρ denotes the spectral function. Both are real-valued functions with different symmetry properties: F ab ðx; yÞ ¼ F ba ðy; xÞ, ρ ab ðx; yÞ ¼ −ρ ba ðy; xÞ. The contributions to the 2PI effective action can be written as [48] Γ½ϕ; Here, the traces include spacetime integration 17 as well as summation over field indices, e.g., Tr C fG −1 0 ðϕÞGg ¼ R xy;C G −1 0ab ðx; y; ϕÞG ba ðy; xÞ. The function G −1 0 is the inverse classical propagator, defined as 17 The subscript C corresponds to an integration of the time variable along the Schwinger-Keldysh contour C. iG −1 0 ab ðx; y; ϕÞ ¼ The 2PI functional Γ 2 ½ϕ; G includes all diagrams, which do not become disconnected after removing two inner (full) propagator lines. This functional is the sum of all such 2PI diagrams with full propagators G and vertices given by the nonlinear and nonquadratic (with respect to φ) parts of S C ½ϕ þ φ, which we denote by S int ½φ; ϕ. Setting ϕ ¼ 0, the inverse classical propagator and the interaction term for our model (2.1) are given by

ð4:11Þ
In a practical computation, to solve the equations of motion (4.4), the effective action needs to be truncated. This can be done in a nonperturbative way using a 1=N expansion, by classifying the contributions to Γ 2 ½ϕ ¼ 0; G based on their scaling with the number of field components N. The diagrams up to NLO have been derived in Refs. [51,52]. There are only two vacuum graphs at leading order [51], shown in Fig. 6(a). These correspond to diagrams with only one vertex in which propagators G aa ðx; xÞ connect legs with the same indices forming closed bubbles. Because of the implicit summation over the index a, each bubble, and consequently each of the two shown diagrams, scales proportional to N. Therefore, we obtain Both terms in (4.12) contribute only as a coordinatedependent shift to the effective mass, and in particular, at leading order (LO), no scattering processes are included (see Appendix B). The NLO contribution consists of an infinite series of diagrams [52], shown in Fig. 6(b). The first row contains all diagrams without sextic vertices. It is the complete NLO contribution in the φ 4 theory and has been resummed in Ref. [35]. In the two lower rows, diagrams include at least one sextic vertex. The latter, whenever encountered, necessarily contains one closed bubble G aa ðx; xÞ to compensate its additional factor of 1=N, while the rest of its legs form the same propagator chains as in the first row. Note that the two diagrams appearing at LO are also encountered at NLO, however, with a different index structure, as explicitly shown in Fig. 6.
To combine the diagrams of the first row with those containing sextic vertices, it is convenient to introduce a coordinate-dependent "modified" coupling, 18 Fðx; xÞ; ð4:13Þ which plays the role of the quartic coupling at NLO in the presence of the sextic self-interaction. With this, the resummation of the NLO diagrams proceeds in the same way as for the φ 4 theory and yields [52] where we have definedBðx;y;GÞ¼δðx−yÞþi~λ ðxÞ 6N G 2 ðx;yÞ, with G 2 ðx; yÞ ¼ G ab ðx; yÞG ab ðx; yÞ, and the logarithm sums the infinite series Leading-order (a) and next-to-leading-order (b) diagrams 18 In Eq. (4.13), we have used the fact that ρ aa ðx; xÞ ¼ 0, and we have defined Fðx; xÞ ¼ F aa ðx; xÞ=N.
The arguments of this section can be generalized to any higher-order self-interactions of even power. For a general interaction ð2nÞ!N n−1 ðφ a ðxÞφ a ðxÞÞ n , the NLO modified quartic coupling is given bỹ 6ðn − 1Þλ 2n ð2n − 1Þ! ðFðx; xÞÞ n−2 : ð4:16Þ Stated differently, at NLO, the higher power selfinteractions contribute as a coordinate-dependent shift to the quartic coupling, given by (4.13) and (4.16), while topologically the diagrams are the same as for the quartic scalar theory. In general, the modified coupling can be positive or negative, depending on the signs of the coupling parameters and on the magnitude of Fðx; xÞ, which is connected to the system's occupation numbers. We have already shown the time evolution of the modified coupling in previous sections in Figs. 3 and 5 for the sextic theory with attractive quartic interaction, and we have observed qualitatively different dynamics depending on the sign of the modified coupling. To get closer to an interpretation in terms of quasiparticles, we derive the corresponding kinetic theory in the following.

Effective kinetic equation
In a kinetic framework, the time evolution of the distribution function is given by the Boltzmann equation. It has the form ∂fðt; pÞ ∂t ¼ C½fðt; pÞ; ð4:17Þ where the collision integral C½f represents a sum over all possible scatterings with in-or outgoing momentum p. For instance, in a coupling expansion, the leading-order perturbative expression for the 2 ↔ 2 collision integral in φ 4 theory is given by [37] C 2↔2 ½fðt; pÞ with the shortcut notation f p ¼ fðt; pÞ, the Bose enhancement factors f p þ 1, and, writing R q ¼ R d 3 q=ð2πÞ 3 , with the integration measure Z dΩ 2↔2 ½fðt; p; q; r; lÞ Such a perturbative approach is not expected to provide an accurate description of the time evolution of the highly occupied infrared modes, which dominate the dynamics in our model. Because of jλjf ≳ 1, infinitely many 2 ↔ 2 scattering diagrams become as important as the perturbative expression (4.19). Remarkably, kinetic theory can still be extended to the highly occupied regime, with the help of a systematic resummation of scattering processes by use of the 1=N expansion to NLO that was discussed above. This has been done for the φ 4 theory [36,37,53] by deriving the corresponding truncated 2PI equations of motion and recasting them into a form similar to the Boltzmann equation, under additional assumptions such as a comparably smooth time evolution, an effective memory loss, and on-shell quasiparticles. In this section, we generalize that derivation to our case. Since at NLO our model can be viewed as a φ 4 theory with a coordinate-dependent (more specifically timedependent) quartic coupling, the form of the Boltzmann equation is similar to that of the φ 4 theory [37], with the replacement λ ↦λ. The details of the derivation are presented in Appendix B. Here, we state the final expression for the integration measure of the Boltzmann equation: 19   FIG. 7. On the left-hand side, the self-consistent equation for the s-channel part of the effective coupling is depicted. Filled circles represent the NLO effective four-vertex (4.21), while dashed circles correspond to the modified couplingλðxÞ from (4.13), the diagrammatic representation of which is shown on the right-hand side. 19 Note that it includes only particle number conserving 2 ↔ 2 scatterings. The first on-shell number-changing processes, such as 4 → 2 scatterings, appear only at NNLO in the 1=N expansion. Off-shell inelastic processes, such as 3 → 1 processes, exist in the 2PI equations of motion already at NLO [36], where they are found to be important for final thermalization [54]. However, they are not taken into account in the Boltzmann equation. Z dΩ 2↔2 ½fðt; p; q; r; lÞ ¼ Z lqrλ 2 eff ðt; p; q; r; lÞ 6N ð2πÞ 4 δðp þ l − q − rÞ δðω p þ ω l − ω q − ω r Þ 2ω p 2ω l 2ω q 2ω r : ð4:20Þ A time-and momentum-dependent effective four-vertex has been introduced, λ 2 eff ðt; p; q; r; lÞ ¼λ 2 ðtÞ which is a consequence of the infinite series of 2PI diagrams at NLO (4.15). The three terms correspond to different scattering channels, where the s scattering channel (p þ l term) is illustrated in Fig. 7 as a self-consistent diagrammatic equation. Expanding the left equation in the figure, the effective vertex appears to be a geometric series of "chain" diagrams with double propagator lines that connect modified couplingsλðtÞ of Eq. (4.13). Therefore, the expression for the effective vertex in (4.21) can be regarded as an analytic continuation of the geometric series. The "one-loop" retarded self-energy Π R in (4.21) corresponds to two propagator lines and one modified coupling in this picture, i.e., a "chain link," and can be written as [23] Π R ðt;ω;pÞ ¼ lim ð4:22Þ Having found the corresponding Boltzmann equation, we now discuss some of its properties. In the case in which the quartic coupling parameter and the modified coupling are of the same order jλj ≈ jλj, then for sufficiently small f ≪ jλj −1 or very large occupancies f ≫ jλj −1 , such that jΠ R j ≪ 1 or jΠ R j ≫ 1, we have from (4.21) respectively. In both cases, only the square of the coupling enters the equation, and there is no dependence on its sign. In the low-occupancy case, the perturbative Boltzmann equation (4.19) for the φ 4 theory with the prefactor for large N and the replacement λ ↦λ is recovered. In contrast, in the other case of large occupancies, we haveλ 2 eff ðt; p; q; r; lÞ ≪ λ 2 ðtÞ, and therefore the resummation generates an effectively weak coupling. Since jΠ R j contains a factor ofλ, the modified coupling drops out of the entire Boltzmann equation, and it becomes identical to the resummed kinetic equation in φ 4 theory. In Refs. [23,37], self-similar solutions of the form (3.2) have been studied for this regime. Under additional assumptions of particle number conservation and a nonrelativistic dispersion relation ω p ≈ M þ p 2 =2M for typical momenta, the values of the exponents α ¼ 3=2, β ¼ 1=2 have been found. This solution describes the inverse particle cascade that we have observed also for our model in the repulsive and mean repulsive cases in Sec. III with the same scaling exponents.
However, if jΠ R j ∼ 1, the sign of the modified couplingλ becomes important due to the denominators of (4.21). In addition to this, the effective coupling has a pole at Π R → −1. If such a pole is encountered, the approximation based on the 1=N expansion may break down.
To compare with our classical-statistical simulations for negative quartic coupling, we have shown the time evolution of the correspondingly computed modified couplings (4.13) in Figs. 3 and 5, for the 1-, 2-, and 8-component models. As can be seen, the regimes of mean repulsion and mean attraction are consistent with the modified coupling being positive or negative, respectively. WhenλðtÞ ≳ jλj, the system enters the above-mentioned self-similar regime after some time, which is rather well described by the vertex-resummed kinetic description. This does not happen numerically whenλðtÞ ∼ −jλj, where number-changing processes become important. These processes get also enhanced when jλðtÞj ≪ jλj, and the kinetic formulation of the 1=N expansion up to NLO turns out to be insufficient to describe the essential aspects of the dynamics in these regimes.

V. CONCLUSION
We have studied the impact of attractive self-interactions on the nonequilibrium dynamics of relativistic scalar systems with large occupation numbers at low momenta. Focusing on OðNÞ-symmetric field theories in 3 þ 1 dimensions with both quartic and sextic self-interactions, we compared the dynamics of the repulsive (λ > 0) and attractive (λ < 0) models.
In the repulsive theory with quartic and sextic interactions, we recover the same nonthermal attractor that has been observed already for the quartic theory. This observation is in line with the expectation that sextic couplings in 3 þ 1 dimensions represent irrelevant couplings in the renormalization group sense. The initially high population produces in this case a dual cascade in momentum space with self-similar dynamics at low and high momenta. The low-momentum region involves a steep power law and leads to the formation of a Bose-Einstein condensate. Its scaling properties are insensitive to the number of field components.
This evolution is also observed if the classical quartic coupling is negative λ < 0, but the fluctuations are so large that the repulsive sextic interaction dominates and mean interactions are repulsive. In contrast, for smaller fluctuations in the mean attractive regime, the inverse cascade is absent. The particle annihilation rate is enhanced as compared to the repulsive case, and no large spatial structures form, unless some nonvanishing conserved charge is present, which would forbid number-changing processes and lead to the formation of Q balls. In the absence of a charge, we have observed that enhanced particle loss can generate a transition between mean repulsion and mean attraction. This is accompanied by the collapse of the condensate and by a rapid decay of low-momentum modes. We also observed that with an increasing number of field components the system stays longer in the mean repulsive regime.
From the 1=N expansion to NLO, we have seen that the contribution from the sextic interaction, as well as from arbitrary interactions of the form ðφ a φ a Þ n with n ≥ 3, can be viewed as a time-dependent shift of a modified quartic couplingλðtÞ. This explains why the scaling solutions observed in the (mean) repulsive case are unaffected by the presence of higher-order self-interactions. Moreover, the observed transition between mean repulsion and mean attraction occurs approximately when this modified quartic coupling changes its sign.
An important application of our studies concerns the physics of dark matter axions and axionlike particles, which exhibit attractive self-interactions, as well as gravitational attractive interaction. Based on our analysis, we expect inelastic collisions to play an important role in the dynamics of such systems. Studying the impact of the phenomenologically more relevant gravitational interaction is left for future work.
Higher-order correlation functions in the case of a Gaussian state can be expressed in terms of one-and two-points correlators. To avoid quenches at the initial time, we have initialized the modes with an effective mass, obtained by iteratively solving the gap equation To generate a nonvanishing charge in the system, one can modify the conditions that the random numbers c a andc a from (A3) have to satisfy. Especially, if one additionally imposes hc a 0 ;pcb 0 ;p 0 i ¼ Vδ p;−p 0 ðA7Þ for some a 0 > b 0 , then the charge density hQ a 0 b 0 i=V ¼ R p ðf p þ 1 2 Þ will be generated. Because of the large occupation numbers in the employed initial conditions (2.4), we omit the additional vacuum 1=2 in the above formulas. While this does not change our results as long as typical occupation numbers are large f ≫ 1, it allows a continuum extrapolation while avoiding the need to renormalize observables, which is in general not possible in classical-statistical theories [39].

APPENDIX B: TRANSPORT EQUATIONS FROM THE 2PI 1=N EXPANSION
In this Appendix, we describe how the effective kinetic equation (4.20) is derived starting from quantum field theory using the 2PI effective action. Although we consider the potential of (2.1), the discussion can be easily generalized to arbitrary OðNÞ-symmetric polynomial potentials.
In the 2PI formalism, the equations of motion for oneand two-point functions are given by the stationarity conditions (4.4). In the symmetric regime, the macroscopic field is zero, and we are left with the equation for the propagator. The latter can be rewritten, taking into account the decomposition of (4.8), as [37,44] where Σ ab ðx; y; GÞ ¼ 2i is called proper self-energy and is the sum of all oneparticle-irreducible diagrams with two external lines. After decomposing the full propagator G ab into spectral and statistical functions according to (4.5), and making a similar decomposition for the self-energy, with an additional separation of "local" and "nonlocal" parts, ab ðxÞδðx − yÞ þΣ ab ðx; yÞ; ðB3Þ the self-consistent equation (B1) is replaced by an equivalent pair of coupled evolution equations for the statistical propagator and the spectral function [37], where and (4.10) has been used. An effective mass term M 2 ab has been introduced, which in the symmetric regime is given by M 2 ab ðxÞ¼m 2 δ ab þΣ ð0Þ ab ðx;GÞ. These equations are called Kadanoff-Baym equations or 2PI equations of motion [55]. If the expressions for the self-energies are known, one can prescribe the statistical propagator and its derivatives at the initial time (the equaltime spectral function is fixed by the commutation relations) and follow the time evolution of the two-point functions.
Transport equations are obtained from the 2PI equations of motion under several assumptions [37,53]. One of them concerns the effective memory loss of the system about the initial conditions, which allows one to take the limit t 0 → −∞. The next step is switching to relative and center coordinates Furthermore, a gradient expansion to the lowest order is employed, which means that only the lowest-order terms in the number of derivatives with respect to the center coordinates X μ and powers of the relative coordinates s μ are kept. Such a description is expected to be suitable for a comparably smooth time evolution.
Moreover, by virtue of OðNÞ rotations, the propagators are considered to be diagonalized, so that F ab ðx:yÞ ¼ Fðx; yÞδ ab and ρ ab ðx; yÞ ¼ ρðx; yÞδ ab . The two-point functions are Fourier transformed with respect to the relative coordinates where the factor i is included in the definition ofρ to make it real valued. Analogous transformations are done for the self-energies. Taking into account all above-mentioned approximations, the 2PI equations of motion can be recast into the following form [37]: 2p μ ∂FðX; pÞ ∂X μ ¼Σ ρ ðX; pÞFðX; pÞ − Σ F ðX; pÞρðX; pÞ; ðB9Þ Note that, due to the gradient expansion to the lowest order, the effective mass term is not present in these equations. For spatially homogenous systems, the dependence on X reduces to pure time dependence. Therefore, Eq. (B10) implies thatρðpÞ does not depend on time. In other words, in this approximation, the particle spectrum does not change with time. From (B9), the time evolution of the particle momentum distribution can be obtained. The latter can be defined according to [53] fðt; pÞ ¼ where Fðt; pÞ ¼ ðfðt; pÞ þ 1 2 ÞρðpÞ. 20 From (B9), the time evolution of the distribution function defined in that way is given by ∂fðt; pÞ ∂t ¼ Z ∞ 0 dp 0 2π ½Σ ρ ðt; pÞFðt; pÞ − Σ F ðt; pÞρðpÞ: So far, the form of the potential has not been specified. Now, we truncate the self-energies according to the truncation scheme of Sec. IV B 1, which is based on a 1=N expansion to NLO. Taking into account (B2), as well as the expressions (4.12) and (4.14), the LO contribution to the self-energy is given by  20 This definition yields the same on-shell distribution function as our previous definition (2.6) used for lattice simulations if the spectral functionρðpÞ follows a δ-like free-field form.
In the second equation,B −1 sums the "geometric series": As can be seen, the LO contributes to Σ ð0Þ ab according to (B3) and therefore adds only to a coordinate-dependent mass shift to the free evolution, similar to the case of φ 4 theory [37]. To separate the local and the nonlocal parts in Σ NLO , it is useful to split B −1 ðx; yÞ ¼ δðx − yÞ − iIðx; yÞ and to rewrite (B14) as The first term with the ½:: brackets contributes to the effective mass, while the second term represents the nonlocal part of the self-energy, and its diagrammatic representation is shown in Fig. 8. Importantly, the expression of the nonlocal partΣ at NLO precisely coincides with the corresponding expression for the quartic theory [37], with the only difference being that λ is replaced byλðxÞ. By introducing the chain term The last step is plugging the expressions (B19)-(B24) for the self-energies into the equation (B12). We note that in the employed derivative expansion we havẽ λðxÞ ¼λ X þ s 2 ≈λðXÞ þ s μ ∂λðXÞ ∂X μ ; ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f λðxÞλðyÞ q ≈λðXÞ: Therefore, the only difference compared to the φ 4 theory [36,37] is the substitution λ →λðtÞ, where t is the central time coordinate. We state here the final result: C NLO ½fðt; pÞ ¼ Z dΩ 2↔2 ½fðt; p; l; q; rÞ½ðf p þ 1Þðf l þ 1Þf q f r − f p f l ðf q þ 1Þðf r þ 1Þ þ Z dΩ 1↔3 ðaÞ ½fðt; p; l; q; rÞ½ðf p þ 1Þðf l þ 1Þðf q þ 1Þf r − f p f l f q ðf r þ 1Þ þ Z dΩ 1↔3 ðbÞ ½fðt; p; l; q; rÞ½ðf p þ 1Þf l f q f r − f p ðf l þ 1Þðf q þ 1Þðf r þ 1Þ þ Z dΩ 0↔4 ½fðt; p; l; q; rÞ½ðf p þ 1Þðf l þ 1Þðf q þ 1Þðf r þ 1Þ − f p f l f q f r : ðB26Þ In the above expression f p ¼ fðt; pÞ and the integration kernels are given by where Under an additional quasiparticle assumption, using the δ-like free-field form of the spectral functioñ the terms corresponding to off-shell 1 ↔ 3 and 0 ↔ 4 processes vanish, and the term representing elastic 2 ↔ 2 scatterings simplifies to (4.20).