Bose-Einstein Condensation from the QCD Boltzmann Equation

Abstract

We present a novel numerical scheme to solve the QCD Boltzmann equation in the soft scattering approximation, for the quenched limit of QCD. Using this we can readily investigate the evolution of spatially homogeneous systems of gluons distributed isotropically in momentum space. We numerically confirm that for so-called “overpopulated” initial conditions, a (transient) Bose-Einstein condensate could emerge in a finite time. Going beyond existing results, we analyze the formation dynamics of this condensate. The scheme is extended to systems with cylindrically symmetric momentum distributions, in order to investigate the effects of anisotropy. In particular, we compare the rates at which isotropization and equilibration occur. We also compare our results from the soft scattering scheme to the relaxation time approximation.

1. Introduction

The study of quark-gluon plasma (QGP), the phase of strongly interacting matter formed in relativistic nuclear collisions and consisting of quasi-free quarks and gluons, is of increasing relevance in modern physics [1]. It represents a testing ground for the Standard Model, as well as for finite temperature field theory and possible grand unification theories. It is also of cosmological significance, as the early universe was dominated by this phase of matter.

Experiments at the Super Proton Synchrotron (SPS), Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) allow us to probe the energy scales at which the QGP is produced. Inferring its properties and phenomenological behaviour is a central goal of the heavy ion programs at these facilities. The theoretical tools that have been developed to describe it are manifold, as the various stages of a heavy ion collision represent very different physical regimes that demand similarly diverse mathematical formalisms to describe (see Figure 1).

Prior to the collision, the nuclei are accelerated to near-light speed, with a Lorentz factor on the order of 100. They are therefore subject to strong Lorentz contraction along the beam axis. At these energies, the lifetime of gluons emitted from the valence quarks or other gluons is long enough to allow additional emissions of soft gluons from themselves. This process keeps increasing the number density of gluons until saturation occurs as recombination of gluons becomes non-negligible, forming the state of matter called the Color Glass Condensate (CGC) [2,3,4]. This regime of large gluon number can be approximated by classical dynamics.

In the first stage of a collision, a large number of gluons are liberated from the CGC. These gluons form a dense, off-equilibrium state called the glasma. Extensive hydrodynamic analyses of HIC indicate that as the medium expands, rapid thermalization occurs (characteristic time on the order of 1 fm) and a QGP in local equilibrium forms [5,6,7,8]. This rapid thermalization is indicative of strong interactions. As the medium continues to expand and decrease in temperature, it eventually drops below the deconfinement temperature $(T_c\approx 170\mathrm{MeV})$ and hadronization occurs.

Despite longstanding efforts and various approaches to describe the dynamics of heavy ion collisions (see e.g., [9,10]), the rapid equilibration of the QGP remains to be thoroughly understood. Another question that has received a lot of recent attention is the possible formation of a gluon condensate in heavy ion collisions [11,12]. We will address these two questions by adapting and numerically solving the QCD Boltzmann equation assuming the dominance of soft gluon exchange in binary collisions. In this framework we can describe the evolution of the QGP from the early pre-equilibrium stages through thermalization towards freeze-out.

2. The Boltzmann Transport Equation

The fundamental equation of kinetic theory is the Boltzmann transport equation. It is a non-linear integro-differential equation describing the evolution of a distribution function of particles, for our purposes a dilute gas of gluons “in a box”. (Quarks are omitted for conceptual simplicity and also motivated for systems which are gluon-dominated). For a spatially homogeneous system under the assumption that $2\to 2$ processes dominate, it can be written as

$$\begin{array}{c}\hfill {\partial }_{t}f=\frac{1}{2}\int \frac{d^3{p}_2}{{\left(2\pi \right)}^{3}2E_2}\frac{d^3{p}_3}{{\left(2\pi \right)}^{3}2E_3}\frac{d^3{p}_4}{{\left(2\pi \right)}^{3}2E_4}\frac{|{\mathcal{M}}_{12\to 34}{|}^2}{2E_1}{\left(2\pi \right)}^4\delta (p_1+p_2-p_3-p_4)(f_3{f}_4{\overline{f}}_1{\overline{f}}_2-f_1{f}_2{\overline{f}}_3{\overline{f}}_4).\end{array}$$

(1)

Here $f_i$ is the distribution function of particle i with 4-momentum $p_i=(E_i,{\mathit{p}}_{\mathit{i}})$. As shorthand, we write $\overline{f_i}\equiv 1+f_i$.

The transition amplitude $\mathcal{M}$ of binary gluon scattering reads at tree level

$$|{\mathcal{M}}_{12\to 34}{|}^2=72g^4\left[3-\frac{tu}{s^2}-\frac{su}{t^2}-\frac{st}{u^2}\right],$$

(2)

where s, t and u are the familiar Mandelstam variables and g is related to the QCD coupling constant $\alpha$ by $g^2=4\pi \alpha$.

For small scattering angles, $|t|\ll s$ and expression (2) simplifies to

$$|{\mathcal{M}}_{12\to 34}{|}^2\approx 144g^4\frac{s^2}{t^2},$$

(3)

which is to be regulated, e.g., by making the substitution

$$\frac{1}{t^2}\to \frac{1}{{\left(t-{\mu }^2\right)}^2},$$

(4)

where $\mu$ is the screening mass.

While this equation is a challenge to solve, Boltzmann’s H-Theorem guarantees that regardless of the initial condition, the equilibrium distribution function will be a Bose-Jüttner function [13],

$$f_{eq}(x,p)={\left[e^{\frac{p^{\alpha }{u}_{\alpha }\left(x\right)-\mu \left(x\right)}{T\left(x\right)}}-1\right]}^{-1}.$$

(5)

Here T, u and $\mu$ parameterize the temperature, collective flow velocity and chemical potential, respectively.

There is one caveat; there exist “overpopulated” initial distributions (see Figure 2) which contain more gluons than can be “accommodated” in a Bose-Jüttner distribution while maintaining particle number and energy conservation. It has been argued [11] that under the assumption of approximate gluon number conservation, a transient state close to equilibrium may form with a Bose-Einstein condensate.

3. The Fokker-Planck and Relaxation Time Approximations

Under the assumption that small-angle scattering dominates, the RHS of Equation (1) can be approximated as the divergence of a current in momentum space [11],

$$D_{t}f=-{\mathbf{\nabla }}_p·\mathcal{J}\left(\mathit{p}\right),$$

(6)

which is a Fokker-Planck type equation. Here the components of the current $\mathcal{J}$ read

$${\mathcal{J}}_i\left(\mathit{p}\right)=\frac{9}{4\pi }{g}^4\mathcal{L}{\int }_{\mathit{k}}{\mathcal{V}}_{ij}(\mathit{p},\mathit{k})\left\{f_p\overline{f_p}{\nabla }_k^j{f}_k-f_k{\overline{f}}_k{\nabla }_p^j{f}_p\right\},$$

(7)

where

$${\mathcal{V}}^{ij}=\left(1-\mathit{v}·\mathit{w}\right){\delta }^{ij}+\left(v^i{w}^j+v^j{w}^i\right),$$

(8)

and we define $\mathit{p}\equiv {\mathit{p}}_1$, $\mathit{k}\equiv {\mathit{p}}_2$ and denote the corresponding unit vectors by $\mathit{w}\equiv \mathit{p}/p$ and $\mathit{v}\equiv \mathit{k}/k$.

In Equation (7), $\mathcal{L}$ is the so-called Coulomb logarithm emerging for screened interactions with vector boson exchange, $\mathcal{L}={\int }_{q_{min}}^{q_{max}}\frac{dq}{q}=\ln \frac{q_{max}}{q_{min}}$ where $q_{max}$ and $q_{min}$ are cutoffs of the order of the equilibrium temperature T and the screening mass $\mu$ introduced in expression (4), respectively [11]. We take $\mathcal{L}$ to be a constant of order 1 in our analysis.

It is convenient to rescale the time variable in Equation (6) as $\tau =\frac{9}{4\pi }{g}^4\mathcal{L}t$ to eliminate the constant factor in Equation (7). The integral in (7) can then be performed (see Appendix A), yielding

$$\mathcal{J}\left(\mathit{p}\right)=I_a{\mathbf{\nabla }}_{p}f+I_{b}f\overline{f}\widehat{\mathit{p}}+({\mathbf{\nabla }}_{p}f·\widehat{\mathit{p}})\mathcal{I}+({\mathbf{\nabla }}_{p}f\times \widehat{\mathit{p}})\times \mathcal{I},$$

(9)

where $I_a=\int f\overline{f}$, $I_b=\int \frac{2f}{p}$ and $\mathcal{I}\equiv ({\mathcal{I}}_x,{\mathcal{I}}_y,{\mathcal{I}}_z)=\int \mathit{w}f\overline{f}$ are functionals of the distribution function f.

It will be interesting to compare the Fokker-Planck approximation of the Boltzmann equation to the well-known Relaxation Time Approximation (RTA),

$${\partial }_{t}f=\frac{p^{\mu }{u}_{\mu }}{p_0}\frac{f_{\infty }-f}{\tau }.$$

(10)

The RTA is easily solvable (and convergences to the same equilibrium distribution); however it lacks QCD-specific features and, as we will see, yields qualitatively different behavior to the Fokker-Planck approximation, which we argue is more physically motivated.

4. The Method of ’B-Lines’

We have developed a flux-conservative numerical scheme that allows us to readily solve the Boltzmann equation in the Fokker-Planck approximation, (6) + (9), which we call the method of ‘B-lines’. The name is given in analogy to splines, with the ’B’ referring to an efficient parameterization of the distribution function in terms of piecewise Bose functions. We have implemented it for distribution functions spherically and cylindrically symmetric in momentum space; here we will discuss the scheme for the simpler isotropic case.

For spherically symmetric distribution functions, we discretize the momentum grid into bins of width $\Delta$ and construct a piecewise Bose interpolation of f,

$$f^{\left(i\right)}\left(p\right)=\frac{1}{e^{a_{i}p+b_i}-1}.$$

(11)

The domain of $f^{\left(i\right)}$ is $p\in [i,(i+1\left)\right]\Delta$ for $0\le i<M-1$ and $[\Delta (M-1),\infty )$ for $i=M-1$.

A couple of points are in order. Firstly, it should be noted that this approach is equivalent to a linear interpolation of an expedient transformation of the distribution function, $g\equiv ln(\overline{f}/f)$, i.e., $g^{\left(i\right)}\left(p\right)=a_{i}p+b_i$. One of the reasons that we choose to make this transformation is that a piecewise linear interpolation directly in terms of f would not allow us to describe the formation of the Bose-Einstein condensate. Secondly, for equilibrium distribution functions approaching equilibrium our interpolation scheme becomes exact, which is a nice property. An equilibrium distribution in g-space of course is simply a straight line.

Physically, the $a_i$ correspond to local (in momentum space) inverse temperature parameters, and the $b_i$ correspond to the chemical potential. We determine them by sampling the distribution function at the gridpoints,

$$\begin{array}{cc}\hfill f_0& \equiv f\left(\delta \right),\hfill \\ \hfill f_i& \equiv f(\Delta i),0<i<M.\hfill \end{array}$$

(12)

Here $\delta$ is small relative to $\Delta$ but non-zero in order to avoid singularities at the origin.

Having established the details of the initial interpolation, we now consider the process by which we evolve the distribution function in time. We separate our $M+1$ cells on the p-axis at the momenta

$$\begin{array}{cc}\hfill p_0& =\delta ,\hfill \\ \hfill p_i& =\left(i+\frac{1}{2}\right)\Delta ,0<i<M,\hfill \end{array}$$

(13)

such that cell 0 is $[0,\delta ]$, cell M is $[(M-\frac{1}{2})\Delta ,\infty )$ and intermediate cells with index $0<i<M$ are $[p_{i-1},p_i]$. These cell boundaries are “staggered” with respect to the grid used for the interpolation, with the distribution function in each cell being interpolated by two B-lines. This is because as we will see shortly, the first derivative of our interpolation of the distribution function is required to be continuous at our cell boundaries, which is not in general the case at B-line boundaries.

From the B-lines we can easily calculate the particle number (per volume) in each cell,

$$\begin{array}{cc}\hfill n_0& =\frac{4\pi }{{\left(2\pi \right)}^3}{\int }_0^{\delta }\mathrm{d}p{p}^2{f}^{\left(0\right)}\left(p\right),\hfill \\ \hfill n_i& =\frac{4\pi }{{\left(2\pi \right)}^3}{\int }_{p_{i-1}}^{\Delta i}\mathrm{d}p{p}^2{f}^{(i-1)}\left(p\right)+\frac{4\pi }{{\left(2\pi \right)}^3}{\int }_{\Delta i}^{p_i}\mathrm{d}p{p}^2{f}^{\left(i\right)}\left(p\right),0<i<M,\hfill \\ \hfill n_M& =\frac{4\pi }{{\left(2\pi \right)}^3}{\int }_{p_{M-1}}^{\infty }\mathrm{d}p{p}^2{f}^{(M-1)}\left(p\right).\hfill \end{array}$$

(14)

These integrals are combinations of polylogarithm functions depending on the B-line parameters $a_i,b_i$. For convenience we have set the gluon degeneracy to 1; $d_g=16$ can be reinstated as needed.

Now, recall that the Fokker-Planck equation (6) can be written as a continuity equation (conserving both particle number and energy). Thus the rate of change of particle number in cell i,

$${\partial }_{\tau }{n}_i=\frac{4\pi }{{\left(2\pi \right)}^3}{\int }_{p_i}^{p_{i+1}}\mathrm{d}p{p}^2{\partial }_{\tau }f(p,\tau )=\frac{4\pi }{{\left(2\pi \right)}^3}{p}^2\mathcal{J}\left(p\right){|}_{p_i}^{p_{i+1}},$$

(15)

is given by ${\varphi }_{i+1}-{\varphi }_i$, the net radial flux into cell i, where

$${\varphi }_i\equiv \frac{4\pi }{{\left(2\pi \right)}^3}{p}^2\mathcal{J}\left(p_i\right)=\frac{4\pi }{{\left(2\pi \right)}^3}{p}^2(I_a{\partial }_{p}f+I_{b}f(1+f)){|}_{p=p_i}.$$

(16)

For the zeroth cell, the flux ${\varphi }_0$ at $p=0$ is zero by definition. Similarly, the last cell’s rightmost boundary is at infinity, with zero flux through it.

We thus arrive at the following non-linear system of ODEs,

$$\begin{array}{c}{\dot{n}}_0=\varphi \left(\delta \right),\hfill \\ {\dot{n}}_i={\varphi }_{i+1}-{\varphi }_i,0<i<M,\hfill \\ {\dot{n}}_M=-{\varphi }_M.\hfill \end{array}$$

(17)

Note that the integrals $I_a$ and $I_b$ ($\mathcal{I}$ vanishes in the spherically symmetric case), which determine the flux (16) depend non-linearly on all of the $f_i$ and must be updated at each time step. We can readily solve these ODEs using the forward Euler method. Having updated the particle number in each cell, it is straightforward to find the evolved set of B-line parameters.

One advantage of the modified scheme is that it allows us to compute the number of particles in the condensate. Analysis shows that the flux into the zeroth cell becomes non-zero after finite time for overpopulated systems approaching the equilibrium Bose-Einstein distribution. Note that this flux remains well-defined as we take the limit $\delta \to 0$. For large t, the flux becomes proportional to $I_b-I_a/T$, which vanishes for an equilibrium distribution. The particle number in cell 0 has a “regular” contribution, which vanishes for $\delta \to 0$, as well as the condensate contribution.

5. Results

Figure 3, Figure 4, Figure 5 and Figure 6 show the evolution in the special case of spherically symmetric, CGC-inspired initial conditions of the form

$$f\left(p\right)=f_0\frac{1}{e^{\left(p-Q\right)/C}+1},$$

(18)

where $f_0$ and C are constants and Q sets the overall momentum scale. For these figures, we have chosen $f_0=0.225$ and $C=0.05Q$ (and $Q\to 1$), which is a moderately overpopulated initial condition where some $8\%$ of the particles asymptotically condense. We point out the qualitative difference between the Fokker-Planck and Relaxation Time Approximations; in the latter the condensate begins to form immediately, whereas the Fokker-Planck scheme exhibits a characteristic “lag”. The onset time $t_c$ is fairly short, assuming $\mathcal{L}\approx 1$, $\alpha \approx 0.4$ and $Q\approx 2$ GeV we find $t_c\approx 0.2$ fm/c.

Generalizing from spherically symmetric to cylindrically symmetric initial conditions, we are able to explore the effects on anisotropy on the evolution of the distribution function. It is important to differentiate between isotropic distribution functions just boosted out of their rest frame and distribution functions that are “genuinely” anisotropic, i.e., even in their rest frame.

“Genuinely” anisotropic distributions are often parameterized in the form [14,15]

$$f_{\mathrm{iso}}\left(\sqrt{{\omega }^2+\xi p_z^2}\right),$$

(19)

where $\xi >-1$ specifies the anisotropy. Similarly, we consider as a generalization of (18) initial conditions of the form

$$f(\omega ,p_z)=f_0\frac{\sqrt{1+\xi }}{e^{\left(\sqrt{{\omega }^2+\xi p_z^2}+b{p}_z-Q\right)/C}+1},$$

(20)

where $\omega =|p|$, b is a boost parameter and the numerator is a normalization that is convenient with regard to the particle number density.

We extract the equilibration time by studying the entropy, evaluated towards final equilibrium. In particular we would like to compare it to the time taken for the initially anisotropic distribution function to isotropize. To this end, as a measure of the anisotropy of a distribution function, we define the “anisotropy parameter”

$$\alpha =\frac{T_{LRF}^{22}}{T_{LRF}^{33}},$$

(21)

where $T_{LRF}^{\mu \nu }$ is the energy-momentum tensor in the local rest frame. In the rest frame, for a cylindrically symmetric distribution function with some anisotropy, $T^{11}=T^{22}=P_{\perp }$ is the transverse pressure of the fluid, while $T^{33}=P_z$ is the longitudinal pressure. For an isotropic distribution they are equal; thus the ratio $\alpha$ must also approach 1 as the system isotropizes.

In Figure 7 we plot the evolution of the normalized entropy and anisotropy parameters associated with an initial condition of form (20), with parameters $f_0=0.1,C=0.05,\xi =b=0.2$ and $Q=1$. Figure 8 shows a log plot of this evolution.

From a fit, the gradients of the lines in Figure 7 are identical within the uncertainties, which corroborates that the rates of isotropization and equilibration are strongly correlated.

6. Discussion

In summary, we have developed an efficient numerical scheme to solve the relativistic Boltzmann equation for gluons in the small-scattering approximation under the assumption of spherically/cylindrically symmetric initial conditions and spatial homogeneity. Among our findings, we have reproduced results from [11] regarding the formation of a transient Bose-Einstein condensate state for overpopulated, spherically symmetric initial conditions. We have investigated the rate at which an anisotropic distribution function becomes isotropic and compared it to the rate of thermalization. Further, we have compared these results to the relaxation-time approximation to the Boltzmann equation.

Possible directions for future work might include investigating the timescale of hydrodynamization (i.e., the time at which hydrodynamics becomes applicable). Following [16], it would be interesting to explore the relation between Bose-Einstein condensation and Kolmogorov turbulence in the relativistic case. Another follow-up would be to study the non-equilibrium attractor described by [17] for the relaxation-time approximation, and see if a similar phenomenon can be observed in the Fokker-Planck approximation.

Scope for further extension of our scheme exists, and such an extension is planned. In particular, it is desirable to extend the scheme to remove the assumption of spatial homogeneity and describe systems without symmetry assumptions in which the above scheme would essentially represent a single spatial cell. A challenge is the fact that the computational complexity scales geometrically with each additional degree of freedom - the so-called “curse of dimensionality”. (Boltzmann equation solvers as well as hydro-codes typically rely on assumptions of symmetry, and for good reason).