Thermalization of overpopulated systems in the 2PI formalism

TL;DR

This paper investigates how overpopulated scalar systems thermalize using the 2PI formalism, focusing on soft mode growth and transient Bose-Einstein condensate formation, revealing that particle number changing processes prevent BEC formation.

Contribution

It applies the 2PI formalism to study thermalization in overpopulated systems, highlighting the role of particle number changing processes in preventing BEC formation.

Findings

Rapid approach to condensation threshold

Particle number changing processes hinder BEC formation

Soft modes grow significantly during thermalization

Abstract

Based on the two-particle irreducible (2PI) formalism to next-to-leading order in the $1/N$ expansion, we study the thermalization of overpopulated systems in scalar $O(N)$ theories with moderate coupling. We focus in particular on the growth of soft modes, and examine whether this can lead to the formation of a transient Bose-Einstein condensate (BEC) when the initial occupancy is high enough. For the value of the coupling constant used in our simulations, we find that while the system rapidly approaches the condensation threshold, the formation of a BEC is eventually hindered by particle number changing processes.

Equations

$$S[\phi ]=\int d^{4}x\left[\frac{1}{2}{\partial }_{\mu }{\phi }_a{\partial }^{\mu }{\phi }_a-\frac{1}{2}{m}^2{\phi }_a{\phi }_a-\frac{\lambda }{4!N}({\phi }_a{\phi }_a{)}^2\right],$$

$$F_{ab}(x,y)=⟨\{{\phi }_a(x),{\phi }_b(y)\}⟩,$$

$${\rho }_{ab}(x,y)=i⟨[{\phi }_a(x),{\phi }_b(y)]⟩,$$

$$F(x,y)=\int \frac{d^{3}p}{(2\pi {)}^3}F(t,t^{\prime },p)e^{i\vec{p}\cdot (\vec{x}-\vec{y})}.\rho (x,y)=\int \frac{d^{3}p}{(2\pi {)}^3}\rho (t,t^{\prime },p)e^{i\vec{p}\cdot (\vec{x}-\vec{y})},$$

$$\left({\partial }_t^2+p^2+M^2(t)\right)F(t,t^{\prime },p)=$$

$$+{\int }_{t_0}^{t^{\prime }}d{t}^{\prime \prime }{\Sigma }_F(t,t^{\prime \prime },p)\rho (t^{\prime \prime },t^{\prime },p),$$

$$\left({\partial }_t^2+p^2+M^2(t)\right)\rho (t,t^{\prime },p)=$$

$$M^2(t)$$

$${\Sigma }_F(t,t^{\prime },p)$$

$${\Sigma }_{\rho }(t,t^{\prime },p)$$

$$I_F(t,t^{\prime },k)$$

$$-{\int }_{t_0}^{t}d{t}^{\prime \prime }{I}_{\rho }(t,t^{\prime \prime },k)\left(F(t^{\prime \prime },t^{\prime },|\vec{k}-\vec{q}|)F(t^{\prime \prime },t^{\prime },q)-\frac{1}{4}\rho (t^{\prime \prime },t^{\prime },|\vec{k}-\vec{q}|)\rho (t^{\prime \prime },t^{\prime },q)\right)$$

$$\displaystyle\qquad+2\int^{t^{\prime}}_{t_{0}}dt^{\prime\prime}I_{F}(t,t^{\prime\prime},k)\rho(t^{\prime\prime},t^{\prime},|\vec{k}-\vec{q}|)F(t^{\prime\prime},t^{\prime},q)\Big{]},$$

$$I_{\rho }(t,t^{\prime },k)$$

$$\displaystyle\qquad-\int^{t}_{t^{\prime}}dt^{\prime\prime}I_{\rho}(t,t^{\prime\prime},k)\rho(t^{\prime\prime},t^{\prime},|\vec{k}-\vec{q}|)F(t^{\prime\prime},t^{\prime},q)\Big{]}.$$

$$F_{\text{eq}}(t,t^{\prime },p)$$

$${\rho }_{\text{eq}}(t,t^{\prime },p)$$

$$f(t,p)$$

$${\omega }_p(t)$$

$$\displaystyle F(t,t^{\prime},p)\big{|}_{t=t^{\prime}=0}$$

$$\displaystyle\partial_{t}F(t,t^{\prime},p)\big{|}_{t=t^{\prime}=0}$$

$$\displaystyle\partial_{t}\partial_{t^{\prime}}F(t,t^{\prime},p)\big{|}_{t=t^{\prime}=0}$$

$$f_{\text{in}}(p)$$

$$n_{\text{in}}=\frac{A{Q}^3}{6{\pi }^2},{\epsilon }_{\text{in}}=\frac{A{Q}^4}{8{\pi }^2},$$

$$n_{\text{eq}}=\frac{\zeta (3)}{{\pi }^2}{T}^3,{\epsilon }_{\text{eq}}=\frac{{\pi }^2}{30}{T}^4.$$

$$n(t)$$

$$\epsilon (t)$$

$$\displaystyle+\frac{\lambda}{12N}\int_{t^{0}}^{t}\!dt^{\prime}\!\int\!d^{3}z\Big{[}-2I_{F}(t,t^{\prime},|\vec{x}-\vec{z}|)F\rho(t,t^{\prime},|\vec{x}-\vec{z}|)$$

$$\displaystyle\hskip 113.81102pt-I_{\rho}(t,t^{\prime},|\vec{x}-\vec{z}|)\left(F^{2}-\frac{\rho^{2}}{4}\right)(t,t^{\prime},|\vec{x}-\vec{z}|)\Big{]}.$$

$$\int \frac{d^{3}k}{(2\pi {)}^3}A(|\vec{p}-\vec{k}|)B(k)=\frac{1}{32{\pi }^{2}p}{\int }_p^{\infty }du{\int }_{-p}^{p}dv(u^2-v^2)A\left(\frac{u-v}{2}\right)B\left(\frac{u+v}{2}\right).$$

$$M_{\text{ren}}^2(t)=m^2+\frac{\lambda (N+2)}{6N}\int \frac{d^{3}p}{(2\pi {)}^3}(F(t,t,p)-F(0,0,p)).$$

$$f(t,p)\approx \frac{T(t)}{{\omega }_p-\mu (t)},$$

Full text

††preprint: YITP-17-49

I Introduction

Understanding the thermalization process of relativistic heavy-ion collisions remains an outstanding issue in the physics of quark-gluon plasmas (see a recent review Gelis:2016upa and references therein). It is commonly assumed that the initial state of the matter formed after a collision is dominated by gluons whose typical momentum is the saturation scale $Q_{s}$, and whose occupation number $f$ is inversely proportional to the strong coupling constant $\alpha_{s}$. In high energy collisions, such as those at the Large Hadron Collider (LHC), the occupation number may become large, $f\sim 1/\alpha_{s}(Q_{s})\gg 1$, leading to an overpopulated initial state. By this we mean that, given their total energy, the number of the initial gluons, if conserved during the evolution, is too large to be accommodated in the final state in a thermal distribution .

In general, when overpopulated systems are driven towards thermal equilibrium by elastic processes which conserve particle number, the particles in the final state that cannot be accommodated in the Bose-Einstein distribution populate a Bose-Einstein condensate (BEC). The formation of such a condensate has been much studied in various contexts, mainly within kinetic theory (see Semikoz:1994zp ; Semikoz:1995rd ; Lacaze:2001qf ; Blaizot:2013lga ; Epelbaum:2014mfa ; Xu:2014ega ; Meistrenko:2015mda , and references therein for some representative works). Most of these works only consider elastic $2\to 2$ scattering in the Boltzmann equation. Inelastic scattering may qualitatively change the picture Blaizot:2011xf ; Kurkela:2011ti , allowing only for the formation of a transient BEC Blaizot:2011xf . A recent analysis of the specific effect of $2\to 3$ processes, dominated in QCD by collinear splittings, suggests a complete hindrance of BEC formation Blaizot:2016iir .

The validity of kinetic theory to discuss this kind of phenomenon may be questioned. It is known in particular Berges:2004yj , from the way it can be derived from the field equations of motion, that the Boltzmann equation is only an approximate description of the full nonequilibrium evolution, that is best suited for the dilute regime $f(p)\ll 1/\lambda$ where $\lambda\ll 1$ is the coupling constant. Strictly speaking then, the soft momentum regime $f(p\approx 0)\gg 1/\lambda$, which is important for the discussion of the BEC formation, is outside the region of applicability of the Boltzmann equation. While the range of applicability of kinetic equations may be wider than what is suggested by a specific derivation from microscopic physics, it is interesting to explore the dynamics of BEC formation using more complete formalisms. This has been done for instance within classical statistical field theory in Berges:2012us (see also Berges:2012ev ; Berges:2014bba ), where the formation of a transient BEC has been observed, as well as its subsequent decay due to the inelastic processes which are naturally included in this formalism. However, the classical statistical field theory is sensitive to the ultraviolet cutoff, unless the coupling constant is chosen extremely small, $\lambda\ll 1$, corresponding to very large occupations $f\gg 1$ Epelbaum:2014mfa . If one has in mind getting insight into the situation in QCD, where the coupling $\alpha_{s}\sim 0.3-0.4$ is not particularly weak, we should explore other regimes.

In this paper, we investigate the problem of the BEC formation in the two-particle irreducible (2PI) formalism in the scalar $O(N)$ theory to next-to-leading order in the $1/N$ expansion Berges:2001fi ; Aarts:2001yn ; Aarts:2002dj . This framework allows us to explore the entire momentum regime, from the deep infrared where $f(p)\gg 1/\lambda$ all the way to the ultraviolet regime where $f$ exhibits the exponential Boltzmann tail. Moreover, the calculations are not limited to weak coupling, and allow for the study of arbitrary occupations. Unfortunately, in practice, these calculations involve heavy numerics. As a result simulations have been sparse in the literature, often done in 1+1 or 1+2 dimensions Berges:2001fi ; Aarts:2001yn ; Nishiyama:2010wc ; Hatta:2012gq , and in 1+3 dimensions on small lattices Arrizabalaga:2004iw ; Tranberg:2008ae . Very recently, a large volume simulation in 1+3 dimensions has appeared Berges:2016nru , but the focus there was on the so-called non-thermal fixed points and the self-similar behaviours which emerge far-from-equilibrium at small coupling. Instead, we shall work at moderate coupling, so that the system eventually approaches equilibrium in a finite (computer) time, and we focus on whether and how a BEC is formed as we dial the initial occupancy. We shall indeed find that the system approaches condensation for large enough initial occupancies. However, our choice of a relatively large coupling enhances the effect of inelastic processes which eventually hinder the condensation.

II Time evolution equation

The 2PI formalism for non equilibrium dynamics that we shall use in this paper is well developed and extensively discussed in the literature (for a pedagogical introduction see e.g. Berges:2004yj ). In this section, we just summarize the main equations that we solve. We consider a scalar field theory with $O(N)$ symmetry in a non-expanding geometry. The action is given by

[TABLE]

where $a=1,\dots,N$. In the 2PI formalism, the basic building blocks of the non-equilibrium evolution equations are the statistical function

[TABLE]

and the spectral function

[TABLE]

where $\{,\}$ and $[,]$ denote respectively an anticommutator and a commutator. We assume unbroken $O(N)$ symmetry so that the one-point function $\langle\varphi_{a}(x)\rangle$ vanishes at all times, and the 2-point functions take the forms $F_{ab}=\delta_{ab}\,F$, $\rho_{ab}=\delta_{ab}\,\rho$. We also assume that the system is spatially homogeneous and isotropic, and perform a Fourier transformation to momentum space

[TABLE]

where $F(t,t^{\prime},p)$ depends only on the modulus $p=|\vec{p}|$ of the momentum, and similarly for $\rho(t,t^{\prime},p)$. The evolution equations then read

[TABLE]

where $M$ is the effective mass (the local part of the self energy)

[TABLE]

The right hand sides of Eqs. (5) and (6) contain memory effect from the initial time $t_{0}=0$ to the current time $t$ and $t^{\prime}$.

The self energies $\Sigma_{F}$ and $\Sigma_{\rho}$ are given by the sum of 2PI diagrams starting from three loops. We shall however go beyond the skeleton loop expansion, and use here the large-$N$ approximation, and resum the ‘ring diagrams’ which appear at the next-to-leading order in the $1/N$ expansion Aarts:2002dj . In this approximation, the self energies are evaluated as

[TABLE]

where

[TABLE]

The diagrams resummed in this expansion scheme are shown in Fig. 1.

In thermal equilibrium, the statistical and spectral functions take the form:

[TABLE]

where $\omega_{p}$ is the particle energy and $f_{\text{BE}}(\omega_{p})=1/(e^{\beta\omega_{p}}-1)$ is the Bose-Einstein distribution. Since the particle number is not conserved, the chemical potential $\mu$ vanishes in equilibrium. The memory integrals in Eqs. (5) and (6) vanish for these distributions in the long time limit $t,t^{\prime}\to\infty$.

When the system is not in equilibrium, but the quasi-particle picture is valid (which occurs on short time scales), one may define the quasi-particle distribution and the quasi-particle energy from the following equations Berges:2004yj

[TABLE]

where $m_{\text{qp}}$ is the quasi-particle mass.

Motivated by the form of the equilibrium 2-point function (12), we choose the following initial conditions for Eqs. (5) and (6)

[TABLE]

where

[TABLE]

The initial conditions for the spectral function are automatically fixed by the equal-time commutation relation as $\rho(t,t^{\prime})|_{t=t^{\prime}=0}=0$, $\partial_{t}\rho(t,t^{\prime})|_{t=t^{\prime}=0}=1$ and $\partial_{t}\partial_{t^{\prime}}\rho(t,t^{\prime})|_{t=t^{\prime}=0}=0$. The momentum $Q$ is the analog of the ‘saturation momentum’ in the context of heavy-ion collisions. The parameter $A$ characterizes the initial occupancy. We shall be interested in how the nature of the thermalization process changes as we dial $A$ to larger and larger values. A general expectation is that for large enough values of $A$, $A>A_{c}$, a Bose-Einstein condensate is formed. A crude estimate of the critical value $A_{\text{crit}}$ is as follows Blaizot:2011xf . Assuming $m=0$ for simplicity, the initial number density and the energy density are

[TABLE]

per degree of freedom. If the system eventually thermalizes at the temperature $T$ and vanishing chemical potential, we have, for a noninteracting theory,

[TABLE]

Condensation occurs if the dimensionless ratio $n_{\text{in}}/\varepsilon_{\text{in}}^{3/4}$ exceeds the equilibrium value $n_{\text{eq}}/\varepsilon_{\text{eq}}^{3/4}$. This leads to the condition $A\geq A_{\text{c}}\approx 0.15$. It is straightforward to generalize this argument to the case of non interacting massive particles (see also Meistrenko:2015mda ). In Fig. 2, we plot the function $A_{\text{c}}(m)$ which we obtained numerically. Note that this estimate is only an approximation in the case of the fully interacting theory: it ignores indeed the effects of the interactions, aside from that of giving a mass to the excitations.

It is interesting to study whether and how the expectation of BEC formation is borne out in our field theoretic system where the coupling is not too small and the particle number is not conserved. To quantify the latter effect, we follow the evolution in time of the number density $n(t)$, and the time derivative of its integral within a sphere of radius $p$, $J(t,p)$:

[TABLE]

The potential delta-function contribution at $p=0$ is kept in $n(t)$ when we discretize the momentum integral to a sum over different $p$-bins.111We do not introduce a specific field to describe the potential condensate, as done e.g in Ref. Berges:2012us , Eq. (4). If the particle number is conserved, and if no particle accumulates in the state $p=0$, $J$ can be interpreted (to within a factor $4\pi p^{2}$) as the flux of particles in momentum space through the sphere of radius $p$ (counted positively if the particles are towards $p=0$). In the present situation, because of the existence of inelastic processes that do not conserve particle number, $J$ also includes, in addition to the flux just mentioned, contributions from particle creation and annihilation inside the sphere of radius $p$.

We also control the time evolution of the total energy density, given by

[TABLE]

In contrast to $n(t)$, $\varepsilon(t)$ should be strictly conserved, and this serves as a check of our numerical simulations.

III Numerical results

In this section, we present the numerical solutions of Eqs. (5) and (6). We fix $N=4$ and $\lambda=10$ throughout. We choose a rather large value of the coupling constant, $\lambda=10$, in order to achieve thermalization on a relatively short (computing) time Hatta:2012gq . The initial occupancy is chosen uncorrelated to the value of the coupling constant. We consider three different values of $A$: $A=0.1$, $A=1$ and $A=5$. The results depend on two dimensionful quantities: the bare mass $m$ and the parameter $Q$ in the initial distribution function. As we shall explain shortly, we use two different values of the ratio of these parameters, $m/Q=1/5$, or $m/Q=1/2$, depending on the value of $A$.

The calculations are done on a regular lattice in momentum space. We call $a$ the lattice spacing and $\Lambda$ the maximum value of the momentum on the grid. We choose the lattice spacing as $a_{p}=0.2\,m$, and the time step as $a_{t}=0.01/m$. In order to ensure energy conservation, the momentum cutoff $\Lambda$ must be at least several times larger than $Q$, and this requirement becomes more and more severe as $A$ is increased (the final, equilibrium, distribution extends to larger momenta as $A$ grows). We use $\Lambda=2\,Q=10\,m$ (50 grid points) for (a) $A=0.1$ and (b) $A=1$, and $\Lambda=7.5\,Q=15\,m$ (75 grid points) for (c) $A=5$. Note that, according to Fig. 2, the case $A=0.1$ is ‘underpopulated’ and the cases $A=1$ and $A=5$ are ‘overpopulated’. The $A=5$ case is numerically more demanding; in order to make the ratio $\Lambda/Q$ large, we increased the number of grid points and chose a larger ratio $m/Q=1/2$.

For the calculation of the memory integrals in Eqs. (5) and (6), we exploit the fact that the self energies $\Sigma_{F}(t,t^{\prime\prime})$ and $\Sigma_{\rho}(t,t^{\prime\prime})$ decay rapidly as $t-t^{\prime\prime}$ increases, and we introduce a cutoff $t_{c}$ on the $t^{\prime\prime}$ integration: $\int^{t}_{t_{0}}dt^{\prime\prime}\dots\longrightarrow\int^{t}_{t-t_{c}}dt^{\prime\prime}\dots$ with $t_{c}=50/Q$. We have checked that the results are unchanged even if we use the smaller value $t_{c}=25/Q$. The convolution integrals in the collision term is somewhat simplified by using the formula

[TABLE]

Finally, we perform a mass renormalization in order to eliminate the most severe ultraviolet divergences. These appear in the $p$-integral in Eq. (7), which is quadratically divergent due to the vacuum contribution to $F$. At each time step, we perform a simple subtraction and use the renormalized mass

[TABLE]

The same mass renormalization is applied to the total energy (24). Thanks to this and our choice of the cutoff $\Lambda$, the total energy is conserved to better than 1% accuracy in all the results to be presented below. There exist also logarithmic divergences that could be eliminated by a coupling constant renormalization. However these are milder and in practice do not affect the results Arrizabalaga:2004iw .

Finally, we observe that large numbers are generated by a combination of the coupling constant and various factors of $\pi$ emerging from angular integrations. As a very rough estimate, after rescaling all dimensionful quantities by appropriate factors of $Q$, one extracts, in the right hand side of Eqs. (5) and (6) an overall factor $\lambda^{2}Q^{2}/(72N\pi^{4})$ or $\lambda^{2}Q^{2}/(1152N\pi^{6})$ depending on how one estimates the angular integrals222In Eqs. (5) and (6), we estimate the magnitude of $\Sigma$ as follow. In Eq. (8) we isolate the factor $\lambda/3N$, and in Eq. (10) the factor $\lambda/6$, while $F(t,t^{\prime},p)$ is taken to be of order $Q^{-1}$. Collecting these factors together with factors $1/(2\pi^{2})$, or $1/(8\pi^{3})$, coming from the angular integrations, one obtains the estimates mentioned in the text.. We could absorb these factors in a “natural” time scale $\tau$, with $t\approx 16\,\tau$ or $t\approx 210\,\tau$ (in units $Q^{-1}$). These numbers are in fact underestimates, as we shall see, but they indicate that we should expect the dynamics to develop over time units that are typically in the range $100\,Q^{-1}$ or even larger.

III.1 $A=0.1$: underpopulated case

We start our discussion with the underpopulated case, i.e., with $A=0.1$. The results are summarized in Fig. 3. The upper left figure shows the time evolution of the distribution $f(t,p)$, and the right figure shows the temperature $T(t)$ and the chemical potential $\mu(t)$ in units of $Q$ extracted by fitting $f(t,p)$ to the Bose-Einstein distribution in the soft, low-$p$ ($p\lesssim 0.5Q$), region. We find that such a fit is possible already at early times. The smooth evolution of the distribution function, as well as the behaviors of the temperature and chemical potential as a function of time, clearly show the trend towards thermalization. However, it takes a very long time for the system to fully thermalize. This can be seen by noticing that the effective chemical potential $\mu$ is always negative and eventually approaches, very slowly, the equilibrium value $\mu=0$ from below. Also, as shown in the lower figure, the zero mode $p=0$ occupancy grows monotonously and has not reached it equilibrium value at the end of the simulation. The growth of soft modes is mostly due to elastic scattering. In fact the number density (23) is approximately conserved: it only changes by less than 2% during the time span of the simulation.

III.2 $A=1$: overpopulated case

Next we turn to the ‘overpopulated’ case $A=1$, with the main results displayed in Fig. 4. One sees both quantitative and qualitative differences as compared to the $A=0.1$ case. The lower left figure shows that the zero mode occupancy $f(p=0)$ rises sharply at early times, reaches a maximum and then decreases slowly. The particle number is approximately conserved during the growth of the soft modes, and it starts to decrease appreciably at a later time as shown in the lower right figure. By the end of the simulation, the system has lost about 20% of particles. Another visualization of the growth of soft modes is given in the right figure of Fig. 5 where we plot the particle flux. Before $f(0)$ reaches a maximum, the flux in the soft region is positive, meaning that particles are flowing into the soft momentum region. After the peak, the flux turns negative everywhere.

As already pointed out, the growth of soft modes is primarily due to elastic processes, which conserve particle number. Eventually inelastic processes, presumably $2$ to $4$ processes333We have no way at this point to isolate precisely which inelastic processes dominate. In leading order both $2\to 4$ and $1\to 3$ (and their reverse) can contribute. However the $1\to 3$ process can only occur off shell, which is made possible in a heat bath by the thermal width acquired by the quasiparticles. In contrast $2\to 4$ processes can occur on shell, and for this reason we believe these to be the dominant ones., take over and start to eliminate the particles in excess in order to reach the equilibrium distribution. The existence of this “delay” in the action of the inelastic processes is also observed in calculations based on kinetic equations Blaizot:2016iir . A further perspective on this feature is provided by the comparison of the full $1/N$ NLO calculation with a 3-loop 2PI calculation where we keep only the first line of (10) and (11). The corresponding occupancy of the zero momentum mode is shown in Fig. 4, lower left panel. One sees that the inelastic collisions suppress the growth of the low momentum modes sooner than in the full NLO calculation. This we interpret as a result of the screening of the interaction in the full calculation Berges:2008wm , making inelastic processes less effective.

In the upper right figure, we determined the set $(T,\mu)$ from the soft ($p\lesssim 0.5\,Q$) and hard ($p\gtrsim 1.0\,Q$) momentum regions separately. The two determinations converge at late times, indicating that the soft and hard modes approach thermalization. However, it takes an extremely long time (much longer than the simulation time) for the system to reach the true equilibrium state where $\mu_{\text{soft}}=\mu_{\text{hard}}=0$. What is interesting is that, unlike what happens in the $A=0.1$ case, where the chemical potential remains negative, here $\mu_{\text{soft}}$ quickly turns positive and approaches very close to (but never exceeds) the quasiparticle mass $m_{\text{qp}}$ defined by Eq. (15). This is more clearly shown in the left plot of Fig. 5 (magnified version of the upper-right figure in Fig. 4) where $\mu_{\text{soft}}$ and $m_{\text{qp}}$ differ only by less than 2% when they are the closest, which roughly occurs at the time when $f(p=0)$ peaks out. This is reminiscent of the dynamical BEC formation in number-conserving theories Semikoz:1994zp ; Blaizot:2013lga ; Epelbaum:2014mfa ; Meistrenko:2015mda , and in fact, the particle number is approximately conserved up to this point. The spectrum in the infrared is well approximated by the Bose-Einstein, or Rayleigh-Jeans distribution

[TABLE]

but does not exhibit any obvious scaling behavior Semikoz:1994zp ; Berges:2015ixa . Moreover, the subsequent behavior is very different from what it would be if the threshold for BEC was crossed. Since $\mu$ never becomes equal to the quasiparticle mass, $f(p=0)$ does not diverge. Instead, it starts to decrease, and at around the same time, the total particle number $n(t)$ also starts to decrease as a result of inelastic processes, as already discussed. Such a ‘decay of a condensate’ has also been observed in Berges:2012us where $f(p=0)$ exhibits a plateau behavior before starting to decrease. The reason we do not see such a plateau here is presumably due to the strength of the coupling which makes the recombination process more efficient.

III.3 $A=5$: strongly overpopulated case

The results for $A=5$, shown in Fig. 6, do not differ much, qualitatively, from those of the previous case $A=1$. Most comments made for this case could be repeated here. The zero mode occupancy $f(p=0)$ rises more sharply, and faster, and it peaks out earlier than in the $A=1$ case. At around the peak, $\mu_{\text{soft}}$ makes its closest approach to $m_{\text{qp}}$. After that, both $f(p=0)$ and the total particle number start to decrease, the latter dropping by about 40% by the end of the simulation. Given that the value $A=5$ is safely labeled ‘overpopulated’, these results corroborate our discussion above, that the formation of a BEC is hindered by inelastic processes. It is interesting to note that the competition between elastic and inelastic processes forces the system to spend a fair amount of time in the vicinity of the onset of condensation, a situation somewhat reminiscent of that of non thermal fixed points.

IV Discussions and conclusions

Due to technical reasons, it is difficult to go beyond $A=5$ in our present numerical setup. Since our lattice size is limited, already when $A=5$ we had to employ a $Q$ value smaller than in the $A=1$ case in order to make the ratio $\Lambda/Q$ large and ensure cutoff-independence. The case $A=10$, $\lambda=10$ has been studied in the same formalism in Berges:2016nru , and the evolution of $f(t,p)$ reported there is qualitatively similar to the $A=1$ and $A=5$ cases above.444Ref. Berges:2016nru used massless fields $m=0$ and observed that $f(t,p)$ approaches the Bose-Einstein distribution with $\mu=0$. In our simulations, $m$ is kept finite, and this allows us to study in detail the transient regime where $\mu$ approaches $m$ before relaxing to the true equilibrium value $\mu=0$. Actually, in Ref. Berges:2016nru the authors considered the region $1000\geq A\geq 10$ with fixed $A\lambda=100$. When $A\gg 10$, which corresponds to the weak coupling regime, the spectrum can no longer be fitted by the Bose-Einstein or Rayleigh-Jeans distribution. This is a far-from-equilibrium regime where nonthermal fixed points appear and $f(t,p)$ exhibits self-similar scaling behaviors. Our study is complementary to this, being more focused on the regime where one can discuss the formation of a BEC on top of the Bose-Einstein distribution.

In the previous studies of the BEC formation based on the Boltzmann equation, special care has been taken in order to treat the singular behavior near the condensation onset. For instance, Ref. Semikoz:1994zp used a much larger lattice with logarithmic spacings in $p$ to allow soft particles to decay into even softer ones, instead of accumulating them in a single $p=0$ mode as in our simulation. Besides, usually one splits the Boltzmann equation into two equations, one for the zero mode and the other for nonzero modes. It is not clear to us whether such more elaborate treatments can qualitatively change the overall picture that emerges from our calculations. Our result suggests, at least for the range of parameters that we have considered, that the formation of a BEC is hindered by the particle number changing processes. Amusingly though, condensation is approached, and the competition between elastic and inelastic processes forces the system to spend a relatively long time in the vicinity of the condensation onset, an effect somewhat reminiscent of that of non-thermal fixed points observed for other ranges of parameters.

It would be interesting to extend the present analysis to longitudinally expanding systems as in heavy-ion collisions. The BEC formation in the presence of longitudinal expansion has been studied in Epelbaum:2015vxa ; Berges:2015ixa in different approaches. The NLO-2PI simulation in the longitudinally expanding geometry was performed in Hatta:2012gq but in 1+2 dimensions (see also Tranberg:2008ae ). It should be straightforward to generalize the work of Hatta:2012gq to 1+3 dimensions.

Finally, applications of the 2PI formalism to nonequilibrium phenomena in QCD have been so far quite limited, and there are only a few attempts in the literature which however take into account only low-order 2PI diagrams Nishiyama:2010mn ; Hatta:2011ky . Unfortunately, a systematic all-order resummation scheme of 2PI diagrams (like the $1/N$ expansion in the scalar theory) is lacking in gauge theory.

Acknowledgments

We thank Akira Ohnishi for many discussions. This work was initiated when J.-P. B. was a visiting professor at Yukawa Institute for Theoretical Physics, Kyoto University. S.T. is supported by the Grant-in-Aid for JSPS fellows (No.26-3462). J.-P. B thanks François Gelis and J$\ddot{\rm u}$rgen Berges for useful discussions.