Numerical solution of the nonlinear boson diffusion equation for gluons

TL;DR

The paper investigates rapid gluon thermalization in the early stages of relativistic heavy-ion collisions using the nonlinear boson diffusion equation (NBDE) for occupation numbers $n(\epsilon,t)$. It develops a robust numerical solver in Julia via the method of lines and validates it against analytic constant-coefficient solutions, then extends the model to energy-dependent drift and diffusion coefficients to capture realistic evolution. By introducing an ultraviolet peak in the energy dependence, infrared thermalization slows while ultraviolet thermalization speeds up, removing the condensate-formation overshoot seen with constant coefficients and enabling control over initialization and equilibration times $\tau_{\mathrm{ini}}$ and $\tau_{\mathrm{eq}}$. The results yield a more physically realistic description of time-dependent gluon condensation and thermalization in overoccupied systems, providing a tractable framework for elastic-gluon dynamics in early-time heavy-ion phenomenology and guiding comparisons with kinetic-theory approaches.

Abstract

The nonlinear boson diffusion equation is taken as a basis to account for the fast thermalization of gluons in the initial stages of relativistic heavy-ion collisions. For constant drift and diffusion coefficients with schematic initial conditions, this equation has previously been solved exactly. In order to achieve a more realistic time evolution towards thermalization, energy-dependent transport coefficients are introduced, requiring numerical solutions of the nonlinear equation. Their accuracy is tested against the exact analytical results in the limit of constant coefficients. The consequences for transient gluon-condensate formation through elastic scatterings in overoccupied systems are discussed.

Figures (5)

• Figure 1: Energy-dependent transport coefficients for an equilibrating system of massless gluons. The upper solid curve is the diffusion coefficient $D(\epsilon)$, the lower long-dashed curve is the drift term $v(\epsilon)$, for an initial distribution $\theta(1-p/Q_\text{s})$ with $p=|{p}|\equiv\epsilon$ for massless gluons, the gluon saturation momentum $Q_\text{s}\simeq 1$ GeV, and an equilibrium temperature $T=-D(\epsilon)/v(\epsilon)=443$ MeV that is constant across the whole energy domain (short-dashed line, right scale).
• Figure 2: Comparison of NBDE-solutions for constant (short-dashed at short times, long-dashed at larger times) and energy-dependent transport coefficients (dot-dashed). The energy dependence has the exponential form shown in Fig. \ref{['fig1']} with its peak at $\epsilon_\mathrm{peak} ={1.7}$ GeV and is chosen such that the mean value of the energy-dependent transport coefficients equals the value for constant coefficients: $\langle D(\epsilon) \rangle = D = {0.67}$ GeV$^2$/fm and $\langle v(\epsilon) \rangle = v = {-1.51}$ GeV/fm. In the top frame, the thermalization in the infrared region is seen to be faster for constant transport coefficients. At the saturation momentum $Q_\mathrm{s}$, both solutions thermalize almost equally fast, while for the ultraviolet region, the energy-dependent transport coefficients provide a much faster thermalization than the constant ones. This is clearly visible in the double logarithmic plot, bottom frame.
• Figure 3: Nonequilibrium evolution of a massless gluon system using the solutions of the nonlinear boson diffusion equation with constant transport coefficients (a), (c) are compared to the numerical solutions (b), (d) of the NBDE with energy-dependent transport coefficients. The initial state is $\theta(1-p/Q_\text{s})$ with the gluon saturation momentum $Q_\text{s}\simeq 1$ GeV. The constant transport coefficients are $D\equiv\langle D\rangle=0.67\,\text{GeV}^2/\text{fm}$ and $v\equiv\langle v\rangle=-1.51\,\text{GeV}/\text{fm}$. The equilibrium temperature is $T=-D/v=443\,\text{MeV}$, with the corresponding Bose--Einstein distribution (solid curves). Time-dependent single-particle occupation-number distribution functions are shown at the timesteps shown in (a), with increasing dash lengths, same steps in all four frames. The calculations for constant coefficients in (a), (c) show identical results, but (c) zooms into UV energies; accordingly for the numerical results with energy-dependent coefficients shown in in (b), (d).
• Figure 4: Gluon-condensate fraction in an overoccupied system using constant coefficients (solid curve) compared with the result for two different energy-dependent transport coefficients with $\epsilon_\mathrm{peak} = 1.7$ GeV, an exponential distribution of the energy dependence of $D$ and $v$ as shown in Fig. \ref{['fig1']} (dashed curve) and a rational distribution (dot-dashed curve). The condensate fraction is calculated by imposing particle-number conservation and using Eqs. (\ref{['eq:CondensateFraction']}) and (\ref{['nopart']}). The mean values of the distributions of the energy-dependent coefficients are equal to $D = 0.67$ GeV$^2$/fm and $v = -1.51$ GeV/fm, as used for constant coefficients. The energy dependencies of $D$ and $v$ are adjusted to counteract the fast (slow) thermalization in the infrared (ultraviolet) region for constant transport coefficients and thus prevent an overshoot above the equilibrium limit $0.369958$ (dotted line) of the condensate fraction.
• Figure 5: Effect of different peak locations $\epsilon_\mathrm{peak}$ in the energy distribution of the transport coefficients (Fig \ref{['fig1']}) on the condensate fraction. The further the maximum lies in the ultraviolet region $\epsilon_\mathrm{peak}>Q_\mathrm{s} = 1\GeV$, the smaller the overshoot gets: The thermalization in the infrared region is slowed down and the thermal tail in the ultraviolet develops faster, causing a balanced thermalization across the whole energy domain and preventing an overshoot. The mean values for the transport coefficients are in all cases $\langle D\rangle = 0.67$ GeV$^2$/fm and $\langle v\rangle = -1.51$ GeV/fm.