IP-Glasma Framework

• IP-Glasma is a QCD-based framework that models the initial state in ultra-relativistic heavy-ion collisions using CGC theory and dynamic Yang–Mills evolution.
• It integrates IP-Sat constraints with event-by-event fluctuations to generate realistic energy density profiles and geometric eccentricities for hydrodynamic simulations.
• The framework accurately predicts observables such as gluon multiplicity and flow harmonics, enhancing quantitative QGP phenomenology.

The IP-Glasma framework is a first-principles initial-state model for ultra-relativistic heavy-ion collisions, combining the Color Glass Condensate (CGC) effective theory of high-density QCD matter with an impact parameter–dependent saturation (IP-Sat) model for parton distributions, and event-by-event classical Yang–Mills (CYM) dynamics of gluonic fields. It provides a comprehensive description of the early-time, pre-equilibrium evolution of strongly interacting QCD matter, generating lumpy, fluctuating distributions of energy, momentum, and spatial geometry that serve as event-by-event fluctuating initial conditions for viscous relativistic hydrodynamics. The framework has systematically evolved from a 2+1D (boost-invariant) realization to a full 3+1D (rapidity-dependent) formulation, incorporating state-of-the-art JIMWLK small-$x$ QCD evolution, and underpins much of the current quantitative phenomenology of QGP collectivity, flow harmonics, and early-time dynamics.

1. Theoretical Foundations: CGC, IP-Sat, and Wilson Lines

The IP-Glasma framework is rooted in the CGC effective theory, valid at small Bjorken-$x$, where high gluon occupancy allows the treatment of gluon fields as classical color fields generated by stochastic large-$x$ color charges $\rho^a(x^\mp, \mathbf{x}_\perp)$. The statistical distribution $W_Y[\rho]$ of these charges is evolved in rapidity $Y$ via the JIMWLK equation, implementing quantum small-$x$ evolution. The local color-charge squared fluctuation $\mu^2(\mathbf{x}_\perp)$ is tied to the local saturation scale $Q_s$.

The IP-Sat model provides constraints on $Q_s(x,b)$ from HERA DIS data, parameterizing the dipole–nucleus forward amplitude as

$$N(x,r,b) = 1 - \exp\left[ -\frac{\pi^2}{2N_c} r^2 \alpha_s(\mu^2) xg(x,\mu^2) T_p(b) \right]$$

with $T_p(b)$ the nucleon transverse profile and $xg(x,\mu^2)$ the gluon density.

Wilson lines,

$$V_{A,B}(\mathbf{x}_\perp) = \mathcal{P}\exp\left[ -ig \int dx^{\mp} \frac{\rho_{A,B}(x^{\mp},\mathbf{x}_\perp)}{\nabla_\perp^2 + m^2}\right],$$

encode the eikonal color rotation of high-energy partons traversing the nucleus. The IR regulator $m$ ensures color neutrality.

2. Event-by-Event Fluctuations and Initial Condition Sampling

Each event's geometry is constructed by sampling nucleon positions from a Woods–Saxon distribution, followed by sampling local color charges $\rho^a(\mathbf{x}_\perp)$ as independent Gaussian random variables with variance determined by the local $Q_s^2$ extracted from IP-Sat. Sub-nucleonic fluctuations (“hot spots”) and rapidity correlations are also implemented where required.

The sampled color charge densities are exponentiated to define Wilson lines for each nucleus, typically discretized into $N_y$ longitudinal sheets (with $N_y$ ranging from 10–50). The Gaussian two-point correlator,

$$\langle \rho^a_{k}(\mathbf{x}_\perp) \rho^b_{l}(\mathbf{y}_\perp) \rangle = \delta^{ab} \delta_{kl} \delta^{(2)}(\mathbf{x}_\perp - \mathbf{y}_\perp) \frac{g^2\mu_{A(B)}^2(\mathbf{x}_\perp)}{N_y},$$

determines the quantum fluctuations at the sub-nucleonic scale.

For 3D (rapidity-dependent) implementations, Wilson lines are further evolved in rapidity using the JIMWLK equation, with stochastic Langevin updates incorporating Gaussian noise and the appropriate kernel regulated by $m$.

3. Classical Yang–Mills Evolution and Glasma Field Dynamics

Following initialization, the classical (sourceless) Yang–Mills equations,

$[D_\mu, F^{\mu\nu}] = 0,$

are solved numerically on a transverse (2+1D) or full three-dimensional (3+1D, in $\tau$, $x$, $y$, $\eta$) lattice. Typical parameters are $N_x=N_y=512$–$1024$, lattice spacing $a=0.02$–$0.10$ fm, and infrared regulator $m=0.2$–$0.5$ GeV.

In 2+1D (boost-invariant) simulations, the initial gauge fields $A^i = A^i_A + A^i_B$ and $A^\eta = (ig/2)[A^i_A, A^i_B]$ are evolved with vanishing $A_\eta$ gradients. In 3+1D realizations, genuine rapidity dependence is implemented by constructing pure-gauge pre-collision fields in all three spatial directions. The Gauss law constraint $[D_i, E^i] + [D_\eta, E^\eta] = 0$ is solved numerically for the initial transverse electric field, typically via a covariant Poisson equation for a scalar potential $\phi$.

Time stepping is accomplished via leapfrog integration, updating electric and magnetic fields and the link variables using the lattice Hamiltonian equations in temporal or Fock–Schwinger gauge.

4. Observables: Multiplicities, Geometry, and Event-by-Event Fluctuations

Key event-by-event observables computed in the IP-Glasma framework include:

• Gluon multiplicity: The gluon spectrum is extracted from equal-time field correlators in Coulomb gauge after CYM evolution,

$$\frac{dN_g}{dy} = \frac{2}{(2\pi)^2N_c}\int d^2k_\perp \left[ \frac{g^2}{\tau} \mathrm{tr}(E_i(k)E_i(-k)) + \tau\, \mathrm{tr}(\pi(k)\pi(-k))\right]/\tilde k_\perp.$$

The event-by-event multiplicity distribution is well described by a negative binomial distribution, with parameters directly determined by the fluctuating initial Glasma fields. For instance, $k/\bar n = \mathcal{O}(1)$ in central collisions, decreasing in more peripheral events.

• Energy density: The lattice energy density at time $\tau$ reads

$$\varepsilon(\mathbf{x}_\perp, \tau) = T^{\tau\tau} = \frac{2}{g^2 a^4}(N_c - \mathrm{Re}\,\mathrm{tr}\,U_{12}) + \frac{1}{a^4}\,\mathrm{tr}\, \pi^2,$$

which exhibits spatial granularity on scales $\sim1/Q_s$.

• Eccentricities: Energy-weighted eccentricities up to $n=6$,

$$\varepsilon_n e^{in\Psi_n} = \frac{\int d^2 x_\perp\, r^n e^{in\phi}\, \varepsilon(x_\perp)}{\int d^2 x_\perp\, r^n \varepsilon(x_\perp)},$$

quantify the initial geometry and drive final-state harmonic flow, with fluctuation-dominated odd moments $\varepsilon_3, \varepsilon_5$ systematically larger than in MC-KLN or wounded-nucleon models.

• Longitudinal structure (3D): 3+1D implementations generate rapidity-dependent profiles, enabling computation of rapidity correlations, two-point $C(\eta_1,\eta_2)$, and harmonic decorrelation ratios $r_n(\eta_a,\eta_b)$.

5. Coupling to Hydrodynamics and Hadronic Cascade

The energy-momentum tensor $T^{\mu\nu}$ of the Glasma fields at a switching time $\tau_0$ (typically $0.4$–$0.6$ fm/c) is diagonalized (“Landau matching”) to extract the local energy density $\varepsilon$ and flow velocity $u^\mu$. These serve as initial conditions for event-by-event relativistic viscous hydrodynamics, e.g., the MUSIC code, which is run with constant or temperature-dependent shear viscosity $\eta/s$ (e.g., $0.095$) and bulk viscosity $\zeta/s(T)$ tuned to single-particle spectra and flow harmonics.

Hydrodynamic evolution proceeds until the fluid reaches a switching temperature (e.g., $T_{\rm sw}=145$ MeV), after which a hadronization prescription (e.g., Cooper–Frye) generates particle distributions. The subsequent hadronic rescattering and resonance decays are simulated using a microscopic hadronic transport model such as UrQMD.

6. Key Phenomenological Insights

The IP-Glasma framework, when interfaced with hydrodynamics and hadronic afterburners, correctly predicts:

• The centrality and system dependence of charged-hadron multiplicities at both RHIC and LHC energies, with the inclusion of JIMWLK evolution necessary for reproducing the $\sqrt{s}$-dependence from RHIC to LHC (Mäntysaari et al., 28 Aug 2025).
• Realistic distributions of flow harmonics $v_n$ up to $n=5$ and their fluctuations, due to large initial-state eccentricities. The ratio $\varepsilon_3/\varepsilon_2 \sim 0.7$–$0.8$ for mid-central Au+Au, close to experiment (Schenke et al., 2012).
• A more compact initial transverse energy density ($R_{\rm rms}$) than wounded-nucleon (MC-Glauber) models, leading to stronger radial flow and requiring larger bulk viscosity $\zeta/s(T)$ to match final $\langle p_T\rangle$ (Schenke et al., 2018).
• Nontrivial initial flow fields in both transverse and longitudinal directions, with pre-equilibrium longitudinal flow $u^\eta$ of magnitude $\sim0.5\langle u^\perp\rangle$, not captured in purely boost-invariant (2+1D) treatments (McDonald et al., 2017).
• Genuine rapidity-dependent correlations and flow-vector decorrelation observables $r_n(\eta_a,\eta_b)$ in $3+1$D, controlled directly by JIMWLK evolution, agreeing with CMS factorization-measurements (McDonald et al., 2020, Schenke et al., 2017).

7. Extensions and Ongoing Developments

Recent developments and extensions include:

• Efficient numerical methods for $3+1$D CYM evolution in Milne coordinates with realistic treatment of color currents, boundary conditions, and Gauss-law enforcement (Matsuda et al., 2023).
• Fully event-by-event hydrodynamic simulations initialized by $3+1$D IP-Glasma, facilitating systematic exploration of longitudinal observables and the impact of initial-state longitudinal fluctuations (McDonald et al., 2023).
• Incorporation of heavy-quark transport (e.g., charm) initialized and propagated through the Glasma and QGP using the MARTINI event generator, probing early-stage charm–medium interactions (Gale et al., 23 Sep 2024).

Compared to more parametric or geometric models (e.g., TRENTo, MC-Glauber, EKRT), the IP-Glasma framework is distinguished by its direct implementation of small-$x$ QCD evolution, event-by-event geometry fluctuations, and a principled mapping to hydro and hadronic observables. The predictive power and flexibility of the framework have made it a standard tool for quantitative QGP phenomenology and for interpreting the initial conditions in relativistic heavy-ion collisions.