Saturation Scale Fluctuations

• Saturation Scale Fluctuations are defined as local, event-by-event variations of the semihard momentum scale Qₛ that marks the onset of nonlinear gluon dynamics.
• They are modeled using Gaussian or log-normal distributions, reflecting both explicit stochastic fluctuations and emergent geometric variations from proton substructure.
• These fluctuations significantly affect particle multiplicity tails, rapidity correlations, and diffraction observables, influencing predictions in high-energy QCD phenomenology.

Saturation scale fluctuations are event-by-event and local-in-transverse-coordinate variations of the semihard momentum scale $Q_s(x,b)$ that controls the onset of nonlinear gluon dynamics at small $x$. In the Color Glass Condensate (CGC), $Q_s$ is the dynamically generated scale in the high-density regime; its magnitude, spatial profile, and fluctuations determine the strength and coherence of the gluon fields that drive particle production, diffraction, and rapidity correlations. Across current phenomenology, these fluctuations are represented either as explicit stochastic variations of $\ln Q_s^2$, often with a Gaussian or log-normal law, or as emergent consequences of fluctuating proton geometry, hotspot structure, color-charge noise, and stochastic small-$x$ evolution (Albacete et al., 2014, Bzdak et al., 2015).

1. Definition and theoretical setting

Several equivalent criteria define the saturation scale. In the dipole picture, it is the transverse scale at which the dipole scattering amplitude becomes order one, $N(r \approx 1/Q_s,x,b)\sim O(1)$, while in occupation-number language saturation corresponds to $f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$, implying strong classical fields. A standard phenomenological form is

$$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$

with $\lambda \approx 0.2\text{--}0.3$ favored by data, and with nuclear enhancement $Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$ (Albacete et al., 2014).

In the McLerran–Venugopalan framework, the valence color sources are sampled from a Gaussian weight,

$x$0

so the local color-charge density $x$1 fluctuates and directly induces fluctuations of $x$2. In practice, MV-like initial conditions feed into BK/JIMWLK evolution and into models such as IP-Glasma, where $x$3 is sampled event by event (Albacete et al., 2014).

Beyond mean-field evolution, BK/JIMWLK determines the average small-$x$4 dynamics, while stochastic generalizations introduce dispersion of the saturation front. A frequently used coarse-grained variable is

$x$5

with fluctuations modeled by

$x$6

This implies that $x$7 is log-normally distributed, with $x$8 controlling the strength of intrinsic saturation-scale fluctuations (Bzdak et al., 2015, McLerran et al., 2015).

A distinct but related viewpoint treats fluctuations not as an externally imposed $x$9 prior, but as emergent from geometry. In impact-parameter dependent dipole models, $Q_s$0 is proportional to the local thickness profile $Q_s$1. In the proton-hotspot picture, event-by-event fluctuations of hotspot number, width, and position generate local patches of larger or smaller $Q_s$2, so saturation-scale fluctuations become a direct consequence of fluctuating geometry rather than an independent stochastic variable (Mäntysaari et al., 2016, Kumar et al., 2024).

2. Statistical representations of $Q_s$3 fluctuations

A widely used implementation introduces fluctuations directly in $Q_s$4. In high-multiplicity $Q_s$5 studies with IP-Glasma, the event-by-event and point-by-point distribution is taken as

$Q_s$6

with the variance parameterized as

$Q_s$7

Fits to multiplicity distributions from $Q_s$8 to $Q_s$9 TeV require $\ln Q_s^2$0, with weak energy dependence, and imply rare events where $\ln Q_s^2$1 exceeds the mean by more than a factor of five (McLerran et al., 2015).

In $\ln Q_s^2$2 multiplicity phenomenology, the same Gaussian law in $\ln Q_s^2$3 is used, but the fitted width is much larger. For $\ln Q_s^2$4Pb at $\ln Q_s^2$5 TeV, a single rapidity-independent width $\ln Q_s^2$6 is reported to describe ATLAS and ALICE pseudorapidity distributions across centrality bins, after a centrality-dependent normalization $\ln Q_s^2$7 is fixed at midrapidity (McLerran et al., 2015). The related geometrical-scaling analysis reaches the same phenomenological conclusion, stating that a sizable width $\ln Q_s^2$8 is needed to describe ALICE $\ln Q_s^2$9Pb multiplicity data (Praszalowicz, 2016).

Other analyses instead embed the fluctuations in proton substructure. In exclusive $x$0 production, the proton is modeled by $x$1 hot spots whose centers fluctuate with width $x$2, and each hotspot carries a Gaussian profile with width $x$3. The local profile $x$4, and therefore the local saturation scale, fluctuates because the hotspot positions and widths vary event by event (Mäntysaari et al., 2016, Salazar et al., 2021). In multiplicity-dependent $x$5 production, this approach is supplemented by two independent log-normal rescalings,

$x$6

with baseline values

$x$7

so the data favor substantial size fluctuations and only modest $x$8-normalization fluctuations (Salazar et al., 2021).

A further formulation characterizes fluctuations through the gluon distribution itself rather than through dipoles. Defining

$x$9

the effective potential for $N(r \approx 1/Q_s,x,b)\sim O(1)$0 in Gaussian small-$N(r \approx 1/Q_s,x,b)\sim O(1)$1 ensembles becomes

$N(r \approx 1/Q_s,x,b)\sim O(1)$2

which can be rewritten in terms of the multiplicative field $N(r \approx 1/Q_s,x,b)\sim O(1)$3 as a Liouville-type potential,

$N(r \approx 1/Q_s,x,b)\sim O(1)$4

In the MV model, fluctuations pile up near $N(r \approx 1/Q_s,x,b)\sim O(1)$5, whereas after JIMWLK evolution they become approximately scale invariant above the absorptive boundary set by $N(r \approx 1/Q_s,x,b)\sim O(1)$6 (Dumitru et al., 2017).

3. Rapidity asymmetry and multi-particle correlations

In symmetric $N(r \approx 1/Q_s,x,b)\sim O(1)$7 collisions, intrinsic fluctuations of the incoming proton saturation scales generate an event-wise forward–backward tilt of the rapidity distribution. A simple parameterization takes

$N(r \approx 1/Q_s,x,b)\sim O(1)$8

with $N(r \approx 1/Q_s,x,b)\sim O(1)$9 independently drawn from the Gaussian law above. The single-inclusive multiplicity is then modeled by the CGC/$f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$0-factorization-inspired expression

$f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$1

For $f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$2, the event distribution is generically tilted in rapidity; averaging restores symmetry, but the odd variance remains nonzero (Bzdak et al., 2015).

The two-particle rapidity correlation function is defined as

$f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$3

and experiments expand

$f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$4

with $f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$5. The coefficient $f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$6 carries the rapidity-odd asymmetry: in symmetric systems $f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$7, while $f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$8 measures the variance of the event-by-event forward–backward tilt (Bzdak et al., 2015).

Near midrapidity, the key analytic relation between the measured odd-parity variance and the fluctuation width is

$f(k_\perp \lesssim Q_s)\sim 1/\alpha_s$9

with the small-$$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$0 behavior $$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$1. Comparison with ATLAS long-range rapidity correlations in minimum-bias $$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$2 at $$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$3 TeV yields $$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$4, consistent with values inferred from multiplicity-distribution analyses (Bzdak et al., 2015).

The same mechanism extends to genuine $$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$5-particle rapidity correlations. Saturation-scale fluctuations induce a leading structure

$$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$6

where $$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$7 denotes the genuine cumulant. Explicitly,

$$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$8

$$Q_s^2(x,b)=Q_0^2(b)\left(\frac{x_0}{x}\right)^\lambda,$$9

Within the color-domain model, the effective variance decreases as

$\lambda \approx 0.2\text{--}0.3$0

leading to $\lambda \approx 0.2\text{--}0.3$1 and to the prediction that the four- and eight-particle cumulants change sign at intermediate multiplicity (Bzdak et al., 2016).

4. Diffraction, geometry, and the separation of fluctuation classes

Exclusive vector-meson production provides a complementary probe because coherent diffraction measures the average gluon profile, while incoherent diffraction measures its event-by-event variance. In the Good–Walker formalism,

$\lambda \approx 0.2\text{--}0.3$2

so the incoherent cross section vanishes for a nonfluctuating proton (Mäntysaari et al., 2016).

Within IPsat, the dipole cross section is

$\lambda \approx 0.2\text{--}0.3$3

with Gaussian proton profile

$\lambda \approx 0.2\text{--}0.3$4

and $\lambda \approx 0.2\text{--}0.3$5. Event-by-event geometric fluctuations are introduced by replacing $\lambda \approx 0.2\text{--}0.3$6 with a sum over constituent-hotspot profiles. Saturation-scale fluctuations are then added as multiplicative log-normal noise,

$\lambda \approx 0.2\text{--}0.3$7

with $\lambda \approx 0.2\text{--}0.3$8, either at the quark level or on a transverse grid of cell size $\lambda \approx 0.2\text{--}0.3$9 fm (Mäntysaari et al., 2016).

The quantitative conclusion from HERA $Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$0 photoproduction is that geometric proton shape fluctuations dominate the incoherent spectrum for $Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$1, while saturation-scale fluctuations are subleading. With $Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$2, the incoherent cross section increases by $Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$3 for $Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$4 relative to the same geometry without $Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$5 fluctuations, but $Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$6 fluctuations alone remain orders of magnitude below the data. Color-charge fluctuations contribute mainly near $Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$7 and are numerically modest in the measured kinematic window (Mäntysaari et al., 2016).

A more recent hotspot-evolution model reaches a similar conclusion from a different direction. There, the proton thickness evolves with the resolution scale $Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$8 through hotspot splitting, and the local saturation scale is defined by

$Q_{s,A}^2\sim A^{1/3}Q_{s,N}^2$9

Although hotspot evolution generates event-by-event fluctuations in hotspot number, width, normalization, and therefore in $x$00, the measured HERA incoherent $x$01 spectrum for $x$02 is described equally well by saturated bSat/IPSat and linearized bNonSat once realistic geometric fluctuations are included. This indicates that, in that channel and kinematic range, the data are dominated by geometry rather than by saturation-scale-specific nonlinear effects (Kumar et al., 2024).

In proton–proton diffraction, the same distinction appears in Good–Walker language. The diffractive excitation cross section is

$x$03

so saturation suppresses observable fluctuations when $x$04. In central $x$05 collisions at high energy, the distribution $x$06 becomes sharply peaked near the black-disk limit, strongly reducing $x$07; diffraction is then dominated by a ring at moderate $x$08, where $x$09 (Flensburg et al., 2010, Flensburg, 2011). A plausible implication is that the effective event-by-event impact of $x$10 fluctuations is clipped by unitarity in central $x$11, even if the underlying Born-level fluctuations remain broad.

5. Multiplicity, geometrical scaling, and heavy-flavor observables

Multiplicity distributions in small systems are particularly sensitive to the high-$x$12 tail. In IP-Glasma, standard MV color-charge fluctuations and impact-parameter fluctuations generate a baseline multiplicity distribution that is too narrow; adding intrinsic event-by-event $x$13 fluctuations with $x$14 broadens the distribution and reproduces the very-high-multiplicity tail from RHIC to LHC energies. The inferred fluctuations are $x$15, and the analysis predicts $x$16 for $x$17 TeV $x$18 with $x$19 mb (McLerran et al., 2015).

In the geometrical-scaling framework, the inclusive cross section is written as

$x$20

with $x$21. Event-by-event fluctuations of $x$22 are again modeled by a Gaussian in $x$23. If $x$24 in a given window, averaging over fluctuations yields

$x$25

so fluctuations largely renormalize the normalization while preserving the scaling form. The phenomenological fit to $x$26Pb multiplicity data requires $x$27, whereas geometrical scaling in $x$28 inelastic cross sections is consistent with $x$29 up to $x$30 GeV/$x$31 (Praszalowicz, 2016).

Multiplicity-dependent $x$32 production in $x$33 and $x$34Pb accesses the interplay between geometric and normalization fluctuations. The rcBK-evolved dipole amplitude is parameterized as

$x$35

with the local nuclear or proton scale implemented through

$x$36

The observable

$x$37

shows a characteristic forward–backward asymmetry in $x$38Pb. Backward rapidity is faster than linear because midrapidity charged hadrons are more strongly affected by saturation than $x$39, while forward rapidity is sublinear or close to linear because the nuclear $x$40 suppresses $x$41 in high-activity events. Artificially weakening saturation removes this asymmetry, and the data prefer substantial hotspot-size fluctuations $x$42 together with modest normalization fluctuations $x$43 (Salazar et al., 2021).

6. Broader dynamical consequences and open directions

In heavy-ion initial-state modeling, saturation-scale fluctuations regulate not only multiplicity tails but also the granularity of the initial energy density. In MC-EKRT, the effective saturation mechanism is implemented by rejecting candidate dijets that overlap within the transverse exclusion radius

$x$44

so the saturation strength fluctuates because minijet multiplicities are Poisson sampled and because hotspot substructure amplifies local nuclear thickness. The resulting local saturation-scale fluctuations are necessary for describing the large-multiplicity tail in Pb+Pb multiplicity distributions, and hotspots enhance the ratio $x$45 by suppressing $x$46 more strongly than $x$47 (Hirvonen et al., 2024).

A distinct recent development concerns event-by-event transverse-momentum correlations in geometrical scaling. There, fluctuations of the saturation momentum $x$48 induce a single-mode structure in the semi-inclusive spectrum,

$x$49

with

$x$50

The two-particle transverse-momentum correlation then factorizes as

$x$51

and the integrated fluctuation measure satisfies

$x$52

This suggests a possible connection between $x$53-correlation observables and fluctuations of the emission region, potentially testable with HBT analyses (Osada, 1 Jun 2026).

At the formal level, the role of rare fluctuations remains important in the saturation regime. For the dipole $x$54-matrix, rare dilute configurations reduce the exponent of the saturation-region solution at LO, become comparatively less important with running-coupling BK, and re-emerge at full NLO because gluon-loop corrections partially compensate the quark-loop suppression. In the full NLO result, the dominant $x$55 term in $x$56 is reduced by a factor $x$57 when rare fluctuations are included, showing that fluctuation effects remain structurally relevant even deep in the black regime (Xiang et al., 2018).

Several limitations recur across the literature. The choice of fluctuation law is model dependent; impact-parameter dependence and proton-size fluctuations remain only partially constrained; and the relative importance of geometric fluctuations, color-charge fluctuations, and explicit $x$58-normalization noise depends strongly on the observable. Diffraction at moderate $x$59 primarily constrains geometry, while long-range rapidity asymmetries and multiplicity tails are more directly sensitive to $x$60 (Mäntysaari et al., 2016, Bzdak et al., 2015). This suggests that no single observable isolates all components of saturation-scale fluctuations.

A consistent picture nevertheless emerges. Intrinsic and geometric fluctuations of the saturation scale are required to account for the high-multiplicity tails of $x$61 events, the centrality dependence of $x$62 rapidity distributions, and the odd rapidity moments measured in two-particle correlations. At the same time, incoherent diffraction shows that local strength fluctuations of $x$63 are not interchangeable with geometric proton-shape fluctuations. A plausible implication is that future constraints will rely on combining observables that project different fluctuation sectors: rapidity-odd cumulants, multiplicity-dependent heavy flavor, incoherent diffraction near and away from $x$64, and system- and energy-dependent scaling tests (McLerran et al., 2015, Bzdak et al., 2016, Salazar et al., 2021).