Progress in the Glauber model at collider energies

We review the theoretical and experimental progress in the Glauber model of multiple nucleon and/or parton scatterings, after the last 10--15 years of operation with proton and nuclear beams at the CERN Large Hadron Collider (LHC) and with various light and heavy colliding ions at the BNL Relativistic Heavy Ion Collider (RHIC). The main developments and the state-of-the-art of the field are summarized. These encompass measurements of the inclusive inelastic proton and nuclear cross sections, advances in the description of the proton and nuclear density profiles and their fluctuations, inclusion of subnucleonic degrees of freedom, experimental procedures and issues related to the determination of the collision centrality, validation of the binary scaling prescription for hard scattering cross sections, and constraints on transport properties of quark-gluon matter from varying initial-state conditions in relativistic hydrodynamics calculations. These advances confirm the validity and usefulness of the Glauber formalism for quantitative studies of QCD matter produced in high-energy collisions of systems, from protons to uranium nuclei, of vastly different size.


INTRODUCTION
The central goal of high-energy nucleus-nucleus (AA) collisions is to study the collective thermodynamic and transport properties of quarks and gluons, the elementary degrees of freedom of the theory of strong interaction (Quantum Chromodynamics, QCD) [1]. Colliding heavy-ions at center-of-mass (c.m.) energies above tens of GeV is the only experimental way known to produce a large multibody system of deconfined partons, the quark-gluon plasma (QGP), predicted by lattice QCD calculations for energy densities exceeding a critical value of c ≈ 0.5 GeV/fm 3 [2,3]. Such collisions provide the only available means to study the thermodynamics of a non-Abelian quantum gauge field theory in the laboratory.
Larger transverse spatial overlaps of the two incoming ions lead to a larger volume of the created system thereby resulting in "mesoscopic" conditions conductive to partial (local) thermalization and subsequent collective behavior of the produced partons. Hence, the interpretation of data from high-energy heavy-ion collisions relies on a detailed knowledge of the initial QGP matter distribution, resulting from the overlap of the two nuclei colliding at a given impact parameter (b), complemented with a theoretical modeling of its subsequent spacetime evolution. Deriving from the data the transverse size of the created QGP, as well as any event-by-event irregularities arising from density and/or color fluctuations, is crucial to interpret the experimental observations, because the total and local depositions of energy in the collision chiefly influence the initial conditions and the evolution of the produced strongly-interacting system.
The generic method to study the properties of QCD matter relies on measuring the distributions of a variety of observables in AA collisions, and compare them to the same measurements in protonproton (pp) and/or proton-nucleus (pA) collisions, where a QGP is not expected to be produced. Common observables include particle production yields (light hadrons, heavy-quarks, quarkonia, jets, photons, etc.) as functions of their transverse momentum (p T ), pseudorapidity (η), azimuthal angle (φ) [4]. Carrying out such quantitative comparisons among colliding systems of different sizes requires appropriate normalization of the measured distributions by using e.g. the number of participating nucleons (N part ), the number of independent binary nucleon-nucleon (NN) collisions (N coll ), the medium transverse area (A ⊥ ), or the system eccentricity ε n (given by the n th moments of its azimuthal spatial distribution), as described below.
Computing the aforementioned quantities typically relies on describing multiple scatterings of the constituents (nucleons, themselves characterized by parton densities) inside the colliding nuclei, using the so-called Glauber formalism, named after Roy Glauber's pioneering work on the use of quantum mechanical scattering theory to calculate cross sections in pA and AA collisions [5,6]. The Glauber formalism is based on the geometric or eikonal approximation, which assumes projectile nucleons to travel along straight lines (i.e. with negligible momentum exchanges compared to their longitudinal momenta) and to undergo multiple independent subcollisions with nucleons in the target. In the original form of the Glauber approach [7][8][9], the total hadronic cross section of a collision of two nuclei (with A and B number of nucleons, respectively) is given by a (2A + 2B + 1)-dimensional integral where b is the collision impact parameter and s denotes a position in the transverse plane. The interaction probability σ(s) is normalized to give the nucleon-nucleon inelastic cross section σ NN = d 2 s σ(s). The nuclear thickness function T A (b) = dz ρ A (b, z) describes the transverse nucleon density by integrating the nuclear density ρ along the longitudinal direction (z).
In the so-called optical limit of the Glauber model derived in [9], local density fluctuations and correlations are ignored so that each nucleon in the projectile interacts with the incoming target as a flux tube described with a smooth density. The total cross section then reduces to where is known as the nuclear overlap function, normalized as d 2 b T AB (b) = A B by integrating over all impact parameters. Expressions (1) and (2) give identical results for large enough nuclei and/or for sufficiently small values of σ NN . In contrast to optical calculations, Monte Carlo Glauber (MCG) calculations [10][11][12][13][14][15][16][17][18][19] evaluate the phase space of Eq. (1) stochastically by distributing the A, B nucleons of nucleus A and B, respectively, in coordinate space according to their corresponding nuclear densities, separated by an impact parameter b sampled from * dσ/db ∝ b. The nuclear transverse profiles are usually approximated with parametrizations of charge density distributions extracted from low-energy electron-nucleus scattering experiments [20,21] that, for large spherical nuclei, are usually described by 2-parameter Fermi (2pF) (also called Woods-Saxon) distributions, ρ(r) = ρ 0 /(1 + exp( r−R a )), with half-density radius R and diffusivity a. Following the eikonal approximation, the collision is treated as a sequence of independent binary nucleon-nucleon collisions, i.e. the nucleons travel on straight-line trajectories and their interaction probability does not depend on the number of collisions they have suffered before. In its simplest form, an interaction takes place between two nucleons if the distance d between their centers satisfies d < σ NN /π. Alternatives to this so-called "black disk" approximation, using for example a Gaussianlike distribution for the nucleon overlap or more complex forms, are also used as discussed below.
One of the most typical phenomenological applications of the Glauber model is to provide the initial conditions for the number density (or entropy or energy densities, related to the former via the QCD equation of state, EoS) of the medium formed in nuclear collisions as input for hydrodynamic calculations of its subsequent space-time evolution. The initial entropy profile in the transverse plane at midrapidity (η = 0) is typically assumed to be proportional to a linear combination of the number density of particles produced in soft (with yields assumed to scale as the N part pairs) and hard (scaling as N coll ) scatterings [22] s 0 ( r ⊥ ) ≡ dS with the relative weight † often taken as α = 0.2 at the LHC. Two typical snapshots in the (x, y) plane of the medium formed in PbPb and pPb collisions at the LHC are shown in Fig. 1. Right after the collision, the spectator nucleons (in grey in the figure) continue undisturbed with the original longitudinal momenta inside the beam line, whereas the hot and dense QGP formed at midrapidity at the LHC expands longitudinally (transversely) at about (0.6 times) the speed of light. Among the quantities derived from Glauber models, the initial azimuthal anisotropy of the QGP system is an important one as it directly propagates to various flow components in the final state through collective hydrodynamic pressure gradients, and is thereby particularly sensitive to the transport and thermodynamic properties of QCD matter [23]. The harmonics eccentricities ε n of the produced QGP (with n = 2 being its ellipticity, n = 3 triangularity, etc.) are theoretically defined as ε n ≡ r n e inφ w(r, φ)rdrdφ r n w(r, φ)rdrdφ , * Hereafter, for simplicity, the impact-parameter variable b (b NN ) is used to denote the distance between nuclei (or nucleons). † The historical wounded nucleon Model (WNM) [8] for low-energy heavy-ion collisions assumes an entropy deposit only for each nucleon that engages in one or more inelastic collisions, and would correspond to α = 0. Figure 1: Events displays of initial transverse density profiles in the (x, y) plane of PbPb collisions at √ s NN = 5.5 TeV with a triangular-like shape (5-10% centrality, left), and pPb at √ s NN = 8.8 TeV (0-1% centrality, right) generated with TGlauberMC [18]. The colored bidimensional surfaces indicate the local density (fm −2 ), following the weighting given by Eq. (4) with α = 0.2. The grey circles indicate spectator nucleons that do not participate in the collision process.
where the weights w(r, φ) are often taken as the initial density of the medium derived from the Glauber model itself via a linear combination of the underlying N part and N coll distributions, as in Eq. (4). Assuming linear hydrodynamics response, the final harmonic flows determined from the azimuthal distributions of the produced hadrons via dN dφ = N 2π (1 + ∑ n (2v n cos(n(φ − ψ n ))) where ψ n is the flow plane angle for harmonic flow v n , are proportional to the initial eccentricity: v n = κ ε n , with κ ≈ 0.2 depending on the EoS, and deviations being sensitive to the ratio of the medium shear viscosity over entropy density η/s [24][25][26][27][28][29]. Whereas flow studies historically focused on the elliptic component (v 2 ), analyses of the higher Fourier harmonics have blossomed in the last years because of the large flow signals observed in the RHIC and LHC data [30][31][32][33][34]. A typical PbPb event at the LHC, leading to a medium with a triangular-like profile, is displayed in Fig. 1 (left).
Another phenomenologically relevant quantity often derived from MCG models is the average QGP path length L(b) traversed by a given perturbative probe, such as an energetic parton, produced in the collision [35][36][37]. The L(b) dependence of the energy loss suffered by a parton going through a QGP provides valuable information on the jet quenching mechanism and on the medium properties [38].
All Glauber quantities mentioned above (T AA , N coll , N part , A ⊥ , ε n , L) depend on the impact parameter b, which is not directly measurable by the experiments, but is monotonically correlated with the overall multiplicity of produced particles in the collision: A smaller impact parameter, i.e. a more central collision, will on average lead to a higher particle multiplicity. The reaction centrality is usually expressed in percentiles of the total inelastic hadronic cross section, with 0% meaning "most central", i.e. fully headon collisions at b = 0 fm, and 100% meaning "most peripheral", i.e. grazing collisions beyond which there is no QCD interaction ‡ (corresponding to b R A + R B , where R A and R B are the nuclear radii). Experimental measurements performed in intervals of multiplicity can be mapped to centrality ranges using Glauber-based simulations, and from there extended to all other relevant quantities.
This review aims at summarizing the main theoretical and experimental progress in the Glauber formalism over the 10-15 years passed since a previous review summarized the status of the field at the BNL RHIC energies [40]. Key differences at the LHC (with TeV colliding energies per nucleon, up to 50 ‡ Significant electromagnetic (e.m.) interactions of the ions can also happen in "ultraperipheral interactions" [39] at impact parameters larger than the sum of their nuclear radii. times larger than those studied at the previous RHIC energy frontier) compared to previous accelerators, are driven by the fact that the nucleon substructure becomes more important, and that even the particle multiplicities produced in small collision systems (pp, pPb) can reach values as large as those measured in CuCu collisions at 200 GeV [41]. An example of the high-density medium formed in pPb collisions at the LHC is schematically shown in Fig. 1 (right). Many of the Glauber model developments have gone in parallel to experimental and theoretical studies of pp, pA, and AA collisions at the CERN LHC, as well as to prospective studies for future facilities such as the Future Circular Collider (FCC) [42]. The document is organized as follows. Section 2 covers measurements and calculations of the proton and nuclear inelastic cross sections. Section 3 discusses improvements in the description of the proton and nucleus transverse profiles. Section 4 reviews phenomenological applications of the Glauber model. Section 5 examines the experimental methods used to determine the collision centrality and discusses their inherent biases. The review ends with a summary of the main conclusions in Section 6. Appendix A provides some basic technical details of MCG simulations.

INELASTIC PROTON AND NUCLEAR CROSS SECTIONS
A key ingredient of all Glauber calculations, via Eqs. (1) and (2), is the inclusive inelastic nucleonnucleon cross section, σ NN , evaluated at the same c.m. energy √ s NN as the pA or AA collision under consideration. The value of σ NN receives contributions from both (semi)hard parton-parton scatterings (aka. "minijets"), computable above a given p T ≈ 2 GeV cutoff by perturbative QCD (pQCD) approaches, as well as from softer "peripheral" scatterings of diffractive nature, with a scale not very far from Λ QCD ≈ 0.2 GeV. Due to the latter nonperturbative contributions, that cannot be computed from first-principle QCD calculations today, one resorts to phenomenological fits of the experimental data to predict the evolution of σ NN as a function of √ s NN . At high c.m. energies, above a few tens of GeV, any potential difference due to the valence-quark structure is increasingly irrelevant, as the bulk of the pQCD cross section proceeds through gluon-gluon scatterings, and all experimental measurements of pp and pp (as well as nn and np) collisions can be combined to extract σ NN .
The collision-energy dependence of the inelastic cross section σ NN is shown in Fig. 2 from all available measurements today, starting off with fixed-target studies in the range √ s NN ≈ 10-30 GeV performed in the 1970s-90s [43], combined with data at pp (UA5 [44] at √ s = 200 and 900 GeV, E710 [45,46] and CDF [47,48]  , or (ii) from measurements of inelastic particle production data in the central detectors collected with MB triggers. The latter measurements are less precise than the former, as they require an extrapolation, dominated by diffractive contributions, to forward regions of phase space not covered by the detectors. Many Glauber MC codes, such as GLISSANDO [11,14,19], use the COMPETE collision-energy parametrization of the NN cross section [65], as provided in the PDG review [43]. However, such a relatively complex multiparameter expression is only required when aiming at a reproduction of the NN cross section data over the full collision energy range experimentally measured, down to √ s ≈ 1 GeV. For the regime of energies relevant at colliders, the √ s dependence of the experimental σ NN results can be fit to a simpler parametrization, as done in [18], resulting in a = 28.84 ± 0.52, b = 0.0458 ± 0.0167 and n = 2.374 ± 0.123 with goodness-of-fit over number of degrees of freedom χ 2 /N dof = 0.7 (and s given in GeV 2 units). The first row of Table 1 lists  , as well as with the typical 42 ± 3 mb used so far in the RHIC literature [40]. At the top FCC energy of √ s NN = 100 TeV, the expected cross section of σ NN = 107.5 ± 6.5 mb is also in agreement with the value 105.1 ± 2.0 mb derived from the average of various model predictions [66].  [67], and σ pPb = 2.06 ± 0.08 b [68] as well as σ pPb = 2.10 ± 0.07 b [69], respectively. At 5.02 TeV, σ PbPb = 7.62 ± 0.15 b for σ NN = 67.3 ± 1.2 mb. The measured ALICE and CMS inelastic pPb and PbPb cross sections provide therefore an inclusive validation of the Glauber model at TeV energies. This fact further justifies the common application of the Glauber approach to derive pp cross section from cosmic-ray proton-air cross section measurements [64].

DEVELOPMENTS IN THE INITIAL STATE
In the past fifteen years, there have been numerous developments related to the description of the nucleon and nuclei radial profiles in Glauber models including e.g. the incorporation of various sources of event-by-event fluctuations (particularly relevant for collisions involving light-ions and/or protons), subnucleonic degrees of freedom, neutron skin effects, and deformed light-and heavy-ion distributions. The most important advances on these topics are summarized in this section.

Proton Transverse Profile
The transverse profile of the proton (or, generically, of the nucleon) is of key importance in many aspects of the Glauber formalism. First, in nuclear collisions, it determines the NN interaction probability, and its realistic description is particularly important to properly describe proton-nucleus results (where intrinsic fluctuations in the proton shape are more relevant than in AA collisions). Second, it is a prime ingredient of MC event generators for pp collisions such as PYTHIA 8 [70] or HERWIG++ [71] in order to reproduce, through the underlying multiparton interactions (MPIs), the properties of inclusive hadron production both in minimum-bias (MB) collisions -as a sum of the particle production activity from all, central to peripheral, pp collisions-and in the so-called "underlying event" (UE) accompanying hard scatterings at the LHC [72]. Third, the proton transverse distribution is also a basic element in calculations of double-and triple-parton scatterings cross sections in proton and nuclear collisions, where the effective N-parton scattering (NPS) cross section bears a simple geometric interpretation in terms of powers of the inverse of the integral of the pp overlap function over all impact parameters, (with n = 2, 3, .., N for double-, triple-, N-scatterings) [73]. Improving the treatment of pp collisions within a Glauber-like approach based on their impact parameter and underlying parton-parton scatterings has attracted a growing interest in the last years to interpret highmultiplicity pp results where collective (QGP-like) phenomena have been observed in the data [74][75][76], a possibility anticipated in [77].
The most simplistic approach used in Glauber models is to consider a fixed proton shape at all colliding energies, disregarding any varying distributions of its parton contents (valence and sea quarks, gluons) and their correlations. The simplest form is a hard-sphere parametrization with uniform density, ρ(r) = 1 4/3 πR 3 Θ (r − R) with radius consistent with electron-proton scattering fits giving a root-meansquare radius, R rms ≈ 0.8 fm [78]. A profile more consistent with the proton charge form factor is given by an exponential, ρ(r) = 1 8πR 3 e −r/R , reproducing to a large extent the spatial distribution of its valence quarks, with R = R rms / √ 12 = 0.234 fm. A single Gaussian ansatz, although not very realistic, makes subsequent calculations especially transparent and hence was used in some analytical approaches. A double Gaussian ansatz, ρ(r) ∝ 1−β a 3 1 e −r 2 /a 2 1 + β a 3 2 e −r 2 /a 2 2 , which corresponds to a distribution with a small core region of radius a 2 containing a fraction β of the total matter embedded in a larger region of radius a 1 , is a common choice, available e.g. in the PYTHIA 8 MC generator [79].
In MC event generators for pp collisions, parametrizations of the overlap function as a function of impact parameter, rather than the individual radial density of each colliding proton, are usually used. In PYTHIA 8, the pp overlap is parametrized as normalized to one, d 2 b T pp (b) = 1, where r p is a characteristic "radius" of the proton, Γ is the gamma function, and the exponent m varies between a more Gaussian-like (m ≈ 2) to a more peaked exponentiallike (m ≈ 1) distribution. The popular PYTHIA 8 "Monash" parameter settings (tune) uses m = 1.85 [80]. In the HIJING [10] and HERWIG++ [71] generators, the pp overlap is approximated by the Fourier transform of the proton e.m. form factor as where µ ∝ 1/r p is a free parameter, and K 3 is the modified Bessel function of the third kind. In HIJING, µ = 3.9/ σ NN /2π is used to describe the dependence of µ on the effective size of the proton. The importance of properly taking into account the proton transverse profile in pp collisions at LHC energies is illustrated in Fig. 3 (left) where the high-p T hadron yield (normalized to its impact-parameterintegrated value) as a function of centrality measured by the ALICE experiment at √ s = 5.02 TeV, is shown. The multiplicity classes are obtained by ordering events according to the charged-particle response in the ALICE VZERO scintillators (2.8 < η < 5.1 and −1.7 < η < −3.7) for events with at least one charged particle produced at midrapidity (|η| < 1). Central (peripheral) pp collisions feature about ten times larger (smaller) yields than MB collisions, consistent with the expectation from the MPI picture [82]. The centrality dependence of the yields is well reproduced by HG-PYTHIA [83] and (slightly less well) by PYTHIA 8 (Monash tune), where the former basically corresponds to PYTHIA with the proton overlap profile given by Eq. (8).
The above discussion highlights the importance in MCG simulations of the choice of the NN collision profile, whereby two nucleons collide if their impact parameter is less than a given distance parameter D ≈ σ NN /π. The simplest collision profiles used are the hard sphere (HS), , with a fitted to reproduce the σ NN inelastic cross section (e.g. a = 0.92 at RHIC energies) [14], although more involved probabilistic ways to model the NN interaction have been known for a long time [84]. At LHC energies, a modification of the collision profile proposed by Refs. [14,85] uses a parametrization based on the Euler Γ(z) and incomplete Γ(α, z) gamma functions where ω interpolates between the HS (ω → 0) and Gaussian (ω → 1) cases. The choice (G, ω) = (1, 0.4) reproduces the measured σ NN ≈ 73 mb and σ pp,el ≈ 25 mb LHC results.  NN gives σ NN , as commonly done for profile functions. One can see that the PYTHIA 8 Monash and HIJING profiles correspond approximately to the ω = 0.6 choice in Eq. (9).
An extension of the nuclear Glauber approach to pp collisions at the parton-level was described in Ref. [77] where, by analogy to the nuclear case, the thickness function of a proton with N g partons is written as T p (x, ) = N g dz ρ(x, ; z), normalized to d 2 b T p (b) = N g . The overlap function of a pp collision at b can be obtained as a convolution over the corresponding thickness functions of each proton normalised to d 2 b T pp (b) = N 2 g . From the partonic cross section σ gg , one can then define the number of binary parton-parton collisions in a pp collision at impact parameter b: N coll,gg (b) = σ gg T pp (b). From this, the probability of an inelastic parton-parton interaction at impact parameter b can be defined as , where the value of σ gg is obtained by requiring that the proton-proton inelastic cross section obtained from the pp overlap function, d 2 b 1 − e −σ gg T pp (b) = σ pp , matches the σ inel pp ≈ 40-80 mb values measured at RHIC and LHC energies ( Table 1). Values of σ gg ≈ 6 mb are used that are consistent with a simplistic perturbative gluon-gluon cross-section of σ gg = K ·(9/2) π α 2 s /p 2 T for α s ≈ 0.5 at a p T -cutoff of the order 1 GeV, where K ≈ 2 is a factor accounting for higher-order pQCD corrections. The final particle multiplicity density in a pp collision follows the same impact-parameter-dependence as that of the number of binary parton-parton collisions, N coll,gg (b), and the average multiplicity in a pp collision, integrated over all impact-parameters, , with the absolute normalisation dN 0 /dη chosen so as to reproduce the MB pp multiplicity of dN MB /dη ≈ 10 measured at midrapidity at LHC energies [86].
Other extensions of the nuclear MCG simulations to account for subnucleonic degrees of freedom, such as e.g. three constituent valence quarks, exist [16,17,19,87,88]. The number of partonic constituents and the way to distribute them in pA and AA collisions -between the two extremes of bound to individual nucleons (according e.g. to a radially exponential nucleon form factor) or freely-distributed over the nucleus (following a global 2pF profile)-must be chosen so as to reproduce basic experimental quantities (inelastic cross sections σ NN , σ pA , σ AA , overall particle multiplicities), and have a different impact on the event ellipticity and triangularity. Accounting for subnucleonic degrees of freedom is particularly important to generate realistic initial conditions for small QGP systems at the LHC, as discussed in Sec. 4.3.

Fluctuations and Correlations
Traditionally, Glauber studies have just dealt with average nucleon and nuclear transverse densities, but in the last decade more and more studies have incorporated event-by-event shape fluctuations generated by various underlying mechanisms. Such developments have been largely motivated by an increasing number of experimental observations that pointed to the need of enhanced eccentricities fluctuations in order to explain e.g. the large azimuthal harmonics -direct [89], elliptic [90], and triangular [27,[30][31][32][33][34] flows-observed in the data.
There are two obvious sources of quantum fluctuations in MCG simulations: fluctuations in the positions of nucleons within the nucleus, and fluctuations at the subnucleonic level. The former are in principle properly accounted for by conventional MCG simulations [91], and we thus focus mostly on the latter. Nucleons are composite quantum-mechanical systems with varying spatial and momentum configurations of their internal quark and gluon constituents, and the overall transverse area occupied by their color fields changes event by event, a phenomenon often referred to as color fluctuations (CFs) [92][93][94]. Such fluctuations can lead to potentially large changes of the effective collision-by-collision nucleon transverse size that are often not accounted for in MCG codes, notably when the NN interaction is approximated by a simple "black disk" approach (see above). Whereas such CFs tend to average out in central AA collisions, they are of importance for pp and pA collisions. The CFs have been evaluated the- oretically in terms of the cross section for inelastic diffractive processes in pN scattering, often referred to as "Glauber-Gribov" approach [95], generalized to the nuclear case [92], and have been phenomenologically encoded into an event-by-event variation of σ NN given by a probability distribution of the form where σ 0 denotes the mean σ NN value, and Ω its width. The normalization C is computed from the provided input (mean σ NN and Ω) requiring σPdσ/ Pdσ = σ 0 , with the dispersion given by the ratio of inelastic diffraction over elastic cross section at t = 0 (zero momentum exchange). Variations of the nucleon interaction strength naturally lead to increased fluctuations in the number of participant nucleons and binary nucleon collisions, e.g. longer tails in the N part and N coll distributions, compared to the pure eikonal picture. Example distributions for pPb collisions at the LHC with varying width values Ω = 0, 0.5, 1, are shown in Fig. 4 (left). Microscopic refinements in the MCG treatment of fluctuations from diffractive scatterings in pA collisions have been discussed in Ref. [96], which proposed a log-normal parametrization instead of Eq. (10). In addition, CFs explain the existence of partonic configurations of the nucleon with smaller-than-average interaction strength that lead to observable differences in the centrality dependence of single jet production in dAu and pA collisions observed in the data at RHIC and the LHC [97,98]. Such events are characterized by a configuration in which a large fraction of the proton momentum is carried by a single parton, e.g. x 0.1, that is more spatially compact than the average, and their interaction strength is thereby reduced with respect to the eikonal limit [99]. The standard MCG models underestimate the associated number of "peripheral" events with low hadronic activity, while they overestimate the "central" ones with a large hadronic activity, although inclusive jet production rates remain unmodified. Most Glauber-like calculations are based upon the pure eikonal approximation and disregard higherorder terms of the expansion Eq. (1) that account for short-range two-, and three-nucleon correlations (many-nucleon correlations are further suppressed). The impact of realistically correlated NN configurations (centrally correlated nucleon configurations, two-body full correlations, three-body chains) on the medium eccentricity has been studied in Ref. [84], where it was found that their combined effects cancel out and bring the results close to the no-correlations case. The impact of NN correlations and their interplay with diffractive effects on N coll have been estimated in Ref. [100]. This latter work found that such correlations slightly decrease (respectively, increase) by a few percent the N coll values (resp., σ pA and σ AA ) compared to the pure Glauber results, but that such effects are cancelled when taking into account diffractive (aka. "Gribov inelastic shadowing") corrections that act in the opposite direction. The number of nucleons that are diffractively excited in the multiple collisions but revert back to their ground state before the scattering process is completed, both increases the nuclear transparency (i.e. reduce the nuclear cross section) and reduces the N coll results back to the values obtained with the conventional Glauber codes. Another possible source of correlations in MCG simulations is due to the recentering procedure [18] by which the MC setup of the initial nuclear profile without nucleons overlaps is done in such a way so that the c.m. of each nucleus is fixed at a given location in each event. Those correlations are found to be small, in particular for large nuclei [91]. Correlations at the parton level due to the interference of same-color gluons from different nucleons, evaluated with the DIPSY generator [101] based on a BFKL resummation of small-x dipoles [102], were found to result only in few percent effects for pA collisions with heavy nuclei.

Neutron Skin and Isospin Effects
The transverse profile of nucleons inside a nucleus is commonly described by a single 2pF distribution. An ansatz, based on electromagnetically probing the charge distribution (protons) of nuclei, that is however not supported by measurements with strongly-interacting probes that prefer instead two nonidentical distributions for protons and neutrons [103], in particular at the surface of heavy stable neutron-rich nuclei, such as 208 Pb with a neutron excess of N/Z ≈ 1.5 [104,105]. These differences appear because protons around the center of the nucleus feel their common electromagnetic repulsion from all directions resulting in an electrostatic equilibrium at a constant charge density, but the outermost protons at r 6 fm, where the nucleon density begins to drop, need additional "skin" or "halo" neutrons in the periphery to counteract the outward Coulomb repulsion and maintain a sufficient nuclear surface tension.
The nominal heavy-ion species at the LHC is 208 Pb, the heaviest stable doubly magic nucleus and one of the most intensively studied isotopes. While the average charge radius of 208 Pb is known to within ±0.001 fm [21,106], past estimates placed the uncertainty in the neutron radius at about ±0.2 fm [107]. Neutron point density parameters of R n = 6.67±0.03 (stat.) fm and a n = 0.55±0.01 (stat.) +0.02 −0.03 (syst.) fm have been measured by the Crystal Ball collaboration via coherent pion photoproduction [105], while the CERN LEAR experiment reports comparable values of R n = 6.684 ± 0.02 (stat.) fm and a n = 0.571 fm derived from antiproton-nucleus interactions coupled with radiochemistry techniques [104]. These data favor a peripheral neutron distribution in the form of a neutron "skin" rather than a neutron "halo", i.e. the neutron distribution is slightly broader than the proton one because of its larger diffusivity (a n − a p ≈ 0.1 fm), but has the same half-radius as the proton distribution (R p ≈ R n ≈ 6.7 fm). The combined point density distribution for proton and neutrons has been implemented in recent MCG simulations [18] via the weighted sum of the individual 2pF distributions. For peripheral PbPb collisions, this results in a maximum ∼4% increase in N coll and approximately half this percentage for pPb collisions, largely driven by the increase of the central radius in the D2pF compared to the 2pF parametrization. If one focuses instead on the transverse distribution of the underlying d quarks, the dominant flavor in neutrons, larger effects are expected for precise phenomenological studies of isospin-dependent gauge boson (γ, W ± , and Z) cross sections in nuclear compared to proton collisions [108][109][110][111]. Figure 4 (right) shows the average fraction of pp, pn, and nn scatterings for PbPb and pPb collisions versus centrality illustrating the increasing relevance of the neutron density for peripheral collisions.

Light and Heavy Deformed Nuclei
For a number of years, the RHIC machine has been providing collisions with a variety ions beyond the nominal gold nucleus, ranging from the lightest species such as deuteron and helium-3 [112] to the heav-iest ones such as uranium [113], as a means to study the system size dependence of various QGP-related signals. As aforementioned, for spherical nuclei, the probability density distribution in MCG models is sampled from the underlying 2pF or D2pF radial probability functions, and taken to be uniform in azimuthal and polar angles, but light species such as the deuteron ( 2 H and b = 1.177 fm −1 , and r = 2r denoting the distance between the proton and neutron, is often employed [114][115][116]. For 3 H and 3 He nuclei, configurations are computed from Green's function MC calculations using the AV18/UIX Hamiltonian, which correctly sample the position of the three nucleons, including their correlations [117]. Similarly, results of wavefunction-based calculations are available for helium-4, carbon, and oxygen [118]. For slightly deformed nuclei, such as e.g. sulfur, the Fermi distribution is modified with an extra parameter w and a Gaussian term, ρ(r) = ρ 0 . Details on all relevant parametrizations can be found in App. A and in Refs. [15,118,119]. The description of the transverse profile of heavy deformed nuclei starts off with the 2pF expressions modified with an expansion of R in spherical harmonics, ρ(r) , with 16 √ π (35 cos 4 (θ) − 30 cos 2 (θ) + 3), and deformation parameters β 2 (quadrupole) and β 4 (hexadecapole) [21]. The higher harmonic eccentricities (ε n ) of the initial QGP produced in collisions of deformed nuclei, such as U, are particularly sensitive to the parametrization of their profiles. A proper description of collisions of heavy deformed nuclei requires also accounting for their relative, tip-on-tip and side-on-side, orientations. Tip-on-tip collisions produce a smaller elliptic flow but larger particle multiplicities (entropy densities), whereas on the contrary, side-on-side collisions generate a larger elliptic flow but a smaller multiplicity. The scaling of v 2,3 flows with multiplicity in ultracentral collisions (0-1% centrality percentile) in small and deformed systems produced in UU, dAu, 9 BeAu, 9 Be 9 Be, 3 He 3 He, and 3 HeAu collisions at RHIC energies, including or not fluctuations from subnucleonic degrees of freedom, has been theoretically studied in [120,121]. This work indicates that such collisions can help discriminate between different initial entropy densities of the QGP medium formed at RHIC and LHC [122,123]. Implications for the extraction of QGP transport properties such as its viscosity, are further developed in Sec. 4.3.

PHENOMENOLOGICAL APPLICATIONS
The Glauber model has many important phenomenological uses in nuclear collisions at colliders. We consider three typical cases here. First, in the definition of the baseline scalings for comparing hardscattering cross sections in pp, pA, and AA collisions. Second, as an underlying framework for MC event generators used in high-energy heavy-ion and cosmic-ray physics. Third, to provide realistic initialstate conditions of the created QGP for subsequent spacetime evolution in hydrodynamics codes. We succinctly review below the basic ideas and latest progress in these three areas.

Binary Scaling for Hard Scatterings
One of the most extended uses of the Glauber model is to properly normalize the fractional cross sections or yields for the production of a given particle in hard-scattering processes (i.e. partonic processes characterized by mass and/or p T scales above a few GeV) in AA and pA collisions, to be able to compare them to those expected in the simpler pp collisions where no QGP formation is, in principle, expected. For pQCD observables that do not suffer any final-state effects, the assumption of binary scaling allows the extraction of modifications of the nuclear parton distribution functions (PDFs) compared to the free proton ones.
It is informative to recall the basic scaling rules for perturbative scatterings in nuclear collisions [124,125]. For a given hard process A+B → h +X, from the generic Eq. (2) for the inclusive AA cross section, one obtains the following relationship between pp and nuclear collisions where the second expression for the inclusive hard cross section is obtained integrating the former over the impact parameter. The associated minimum-bias invariant yield per nuclear collision, N hard AB = σ hard AB /σ AB , for a given hard process in an AB collision compared to that of a pp collision is N hard AB MB = A·B σ AB · σ hard NN , where σ AB is the inclusive inelastic AB cross section given by Eq. (2). The average nuclear overlap function at impact parameter b for minimum-bias collisions is The corresponding expressions for a given impact parameter b can be obtained by multiplying each nucleon in nucleus A with the density along the z direction in nucleus B, integrated over nucleons in nucleus A, i.e.
Similarly, one obtains a useful expression for the probability of an inelastic NN collision or, equivalently, for the number of binary inelastic collisions, N coll , in a nucleus-nucleus collision at impact parameter b: From this last expression, one can see that the nuclear overlap function, can be thought of as the hard-scattering integrated luminosity (i.e. the number of hard collisions per unit of cross section) per AB collision at a given impact parameter.
The expressions above allow writing the standard binary (or point-like) collision scaling formula that relates the hard-scattering yields in nuclear and proton collisions as N AB (b) = N coll (b) · N pp . The nuclear modification factor for hard-scattering processes is thereby defined as the ratio of AA over scaled pp cross sections and/or yields (here, differential in p T and η) as for minimum-bias collisions, and dependent on b as In the absence of any final-and/or initial-state effects, one expects R AB = 1 for any hard-scattering process. The equality of Eq. (15) or (16) to unity, modulo few percent nuclear PDF effects (see below), for colorless hard probes that do not suffer final-state interactions in the produced QGP was confirmed previously in heavy-ion collisions at SPS and RHIC as well as, in the last years, at LHC energies, and constitutes a validation of the basic assumptions of the Glauber model itself. Prominent examples include the production yields of photons [126][127][128][129], and W and Z bosons [130][131][132][133][134][135][136][137] in pPb and PbPb collisions. Figure 5 shows the R AA values measured for isolated photons [129] and Z bosons [135] as functions of p γ T and rapidity |y Z |, respectively, in PbPb collisions at the LHC. Both ratios are around unity with small variations due to nuclear PDF modifications related either to the increased number of d-quarks in the Pb nucleus compared to protons (isospin effects) and/or to few percent (anti)shadowing effects   Figure 5: Nuclear modification factors R AA measured for isolated photons as a function of p γ T (left) [129] and Z bosons as a function of rapidity (right) [135] in PbPb collisions at the LHC, compared to pQCD NLO calculations with nuclear PDFs.
at the large virtualities (Q 2 ≈ p γ T , m Z ) probed in the underlying partonic scatterings. Under the key assumption of binary scaling, the precision of the PbPb data (a few percent experimental uncertainties in the case of electroweak gauge bosons) allows the derivation of the EPPS16 [138], nCTEQ15 [139], and nNNPDF2.0 [140] nuclear PDFs at next-to-leading order (NLO), or more recently also next-to-nextto-leading order (NNLO) [141] accuracy, through global fits combining nuclear deep-inelastic scattering (DIS) and LHC electroweak boson data. Precision electroweak boson measurements in PbPb collisions can also be used to derive a data-driven normalization, in principle independent of the Glauber model, for cross section measurements in PbPb collisions, as explored in [142].

Heavy-Ion Monte Carlo Event Generators
All existing generic event generators of ultrarelativistic pA and AA collisions -such as HIJING 1.0 [10] and 2.0 [143], EPOS-LHC [144], AMPT [145], QGSJET-II [146], DPMJET-III [147], as well as the more recent PYTHIA 8-based [148] ANGANTAYR [149] code-internally rely on a Glauber picture to model the early stage of the collision through the proper computation of the number of inelastic subcollisions for any reaction centrality. The main differences among models arise from their treatment of the underlying (semi)hard scatterings: minijets in the case of the HIJING, ANGANTAYR, and AMPT codes mostly used in collider physics; and Regge-Gribov "cut pomerons" (aka. parton ladders, giving rise to one or two strings spanned between two colliding nucleons, or between a nucleon and another pomeron) in the case of the EPOS-LHC, QGSJET-II, and DPMJET-III codes that are mostly used in cosmic-ray physics [150]. The final hadronization of partons or strings is carried out via (variations of) the Lund fragmentation model [151] in all generators.
The HIJING generator relies on the eikonal approach to determine the number of inelastic subcollisions of two types: soft NN collisions treated as in the FRITIOF approach [152], and (multiple) hard parton-parton collisions treated perturbatively as in PYTHIA. The transverse momentum cutoff p T,0 ≈ 2 GeV that separates hard from soft scatterings increases slowly with collision energy (logarithmically, similar to the inelastic cross section Eq. (6) evolution), so that the total number of minijets per unit transverse area satisfies p 2 T,0 /π > T AA (b) σ hard /(πR 2 A ), where σ hard is the pQCD cross section for 2 → 2 parton scatterings, and T AA (b) is the overlap function of the AA collision. The probability for an inelastic NN collision is given by dσ NN The AMPT code uses directly the Glauber initial conditions generated by HIJING as input for its parton cascade evolution. The most recent heavy-ion event generator is ANGANTAYR, which follows a Glauber approach similar to that of HIJING but further takes into account Glauber-Gribov corrections, by dividing up each inelastic sub-collision as either single-, double-diffractive, or absorptive (i.e. nondiffractive), and with CF effects implemented through a model with fluctuating nucleon radii resulting in a fluctuating NN cross section inspired by the approach of Ref. [93].
At variance with other generators, EPOS-LHC keeps track of how many times a given nucleon interacts with nucleons from the other nucleus, and separates them event-by-event into the "core" (nucleons that collided more than once) and "corona" (nucleons that interacted exactly once). Event-by-event, a fraction of the string segments that do not overlap (corona) fragment into hadrons normally, following the Lund string model, whereas the other clusters with large density of strings, are used to create a (QGPlike) core that can flow and hadronize collectively. Such a two-component core-corona medium leads to a consistent reproduction of the particle-multiplicity dependence of a number of observables ( p T , ratio of different hadron yields, etc.) in pp, pA, and AA collisions [144]. The subsequent collective expansion of the medium, defined by the initial core particle density, is taken care of by relativistic hydrodynamics equations.
Apart from the hadronic MC event generators mentioned above, generators of ultraperipheral (photoninduced) nuclear collisions, such as STARLIGHT [153] and SUPERCHIC 3 [154] also employ a Glauber approach to determine the non-hadronic overlap probability needed to compute the cross sections of purely exclusive final states. In STARLIGHT, the probability of having no hadronic interactions is given by P no.had (b) = e −σ NN T A (b) , where T A (b) is the nuclear thickness function. For large nuclei, applying this probability function is roughly equivalent to imposing a b > R A condition on their photon flux. This probability is also used when the photon is emitted by the proton, leading to an effective γ flux in pA collisions considerably smaller than if it was calculated directly from the e.m. proton form factor.

Initial Conditions for Hydrodynamic QGP Evolution Calculations
The strongly-interacting medium created in AA collisions at the LHC is a dynamical system that expands, cools down, and transforms into a hadron gas at times around τ = 10-15 fm/c. Under such conditions, the extraction of QGP properties can only be accomplished by comparing the experimental measurements of the hadronic final-state to theoretical predictions that include a model of the full spacetime evolution of the heavy-ion collision process [155]. The current state-of-the-art for the QGP evolution is given by 3-D viscous relativistic fluid dynamics calculations [24,27,29,156,157] -where the plasma thermodynamic properties evolve according to the lattice QCD EoS [2,3] with nonideal corrections encoded in the medium shear viscosity ηmatched to a hadron transport cascade (often, the UrQMD code [158]) once the energy density drops below c . The largest source of uncertainty in the extraction of medium properties from data-theory comparisons lies in the description of the initial state of the QGP [26], a topic that is discussed next.
Hydrodynamic models start their evolution at τ 0 from a given input entropy (or energy) density in the transverse plane s 0 ( r ⊥ ). In principle, such initial conditions (ICs) should be consistently derived by solving the non-equilibrium evolution of the matter created in the first parton-parton interactions, but achieving thermalization of the interacting fields in ultrashort O(1 fm) time-scales remains a difficult theoretical problem in heavy-ion physics [159]. Therefore, one assumes that the produced matter has (pre)equilibrated, and the ICs are usually given either by (i) the number density of produced gluons after the primary collisions, via s 0 ( r ⊥ ) ∝ dN g /(τ 0 d 2 r ⊥ dη), derived analytically e.g. in models such as MC-KLN [160,161] and IP-Glasma [162,163] based on the Color-Glass Condensate (CGC) effective theory for parton saturation in heavy-ion collisions [164], or by pQCD NLO calculations with ad hoc parton saturation such as EKRT [165,166]; (ii) Glauber MC profiles such as those shown in Fig. 1 obtained with the ansatz given by Eq. (4) with α adjusted to match the observed multiplicity distributions, or similarly generated with one of the MC event simulations (HIJING, EPOS, etc.) discussed in Sec. 4.2 .
The corresponding local deposition of the entropy density from the underlying (parton or nucleon) collisions, and thereby the transverse area and azimuthal anisotropies of the produced medium, are then model dependent. The T R ENTo approach [167] has parametrized all different approaches with a generic function of the participant target and projectile thickness functions T part A,B , as with a continuous parameter p that effectively interpolates among different entropy deposition schemes. , is equivalent to the generalized mean ansatz with p = 1. The IP-Glasma approach -which combines CGC effects of the incoming gluon distributions, where the gluon thickness function T g (b) is a Gaussian function of the impact parameter (from the center of the probed nucleon) with a width constrained by HERA DIS data [168], with an event-by-event classical Yang-Mills evolution of the produced glasma gluon fields-deposits density following T part A · T part B . The default entropy deposition parameter p = 0 of T R ENTo, derived from a global fit to multiple experimental data, produces similar initial eccentricities to IP-Glasma [167]. In the KLN model, the gluon multiplicity N g can be determined perturbatively in the k T -factorization CGC approach from the parton saturation momenta of each nucleus Q 2 s,A ∝ T A , leading to s ∝ T min 2 + log(T max /T min ) . This would correspond to a parameter p ≈ −0.67 in the generalized ansatz (17). The EKRT approach combines collinear factorized pQCD minijet production with a simple model of gluon saturation, and predicts an energy density given by 0 ∝ K sat π Q 3 sat with saturation momentum Q sat (K sat ; T A , T B ), corresponding to an exponent p ≈ 0.
Smaller, more negative, values of p pull the generalized mean towards a minimum function, and hence correspond to models with more extreme gluon saturation effects, leading to the following schematic hierarchy of more saturated ICs: WNM < IP-Glasma, T R ENTo, EKRT < MC-KLN. Of course, such a hierarchy only accounts for average density effects, and the different approaches feature also key physics differences that lead to more (or less) QGP shape fluctuations, and thereby larger (smaller) eccentricities ε n that impact significantly the extractions of e.g. the medium viscosity-over-entropy η/s ratio from comparisons of azimuthal flows v n to the hydrodynamic predictions. In particular, the IP-Glasma model also generates the full energy-momentum tensor of the medium, with momentum anisotropies with a length scale (of the order of Q −1 s (x) = 0.1-0.2 fm) smaller than those present in other calculations (0.4-1 fm), resulting in a finer structure of the initial entropy density compared to the MC-KLN and MC-Glauber models [162].
At the LHC, elliptic v 2 and triangular v 3 flows have been studied e.g. by the ALICE experiment in XeXe and PbPb collisions at √ s NN = 5.44 and 5.02 TeV, respectively. The ratios of v 2,3 /ε 2,3 as a function of particle transverse density (given by (1/A ⊥ )dN ch /dη), where ε n and A ⊥ are derived from the ICs of various models described above, are shown in Fig. 6. The hydrodynamic expectation is that v n /ε n Figure 6: Comparisons of the ratio v n /ε n (for n = 2, 3) as a function of particle transverse density, (1/A ⊥ )dN ch /dη, in XeXe and PbPb collisions at the LHC, where v n and dN ch /dη have been measured by the ALICE experiment, and the ε n and A ⊥ parameters are derived using the ICs of various models described in the text. Figured adapted from Ref. [169].
increases monotonically with the transverse density across different collision energies and systems, and a violation of such a scaling may indicate an incorrect modeling of the initial transverse area A ⊥ and/or the azimuthal anisotropies ε n . The results of Fig. 6 indicate that the standard MCG using nucleons and the MC-KLN model (first and second panels) fail to reproduce the expected scalings for v 2 (red symbols), whereas the T R ENTo model with p = 0, equivalent to IP-Glasma, as well as the MCG with constituent partons (third and fourth panels) feature better scaling behaviors accross flow coefficients and systems (although a drop at the largest densities is observed). These results illustrate the type of constraints imposed by the data on IC medium models, which suggest, in this case, the need of a higher number of subnucleonic sources in order to achieve a steady increase of v 2,3 /ε 2,3 for more central collisions.
Implications for the extraction of QGP transport properties, such as its viscosity constrained from the experimental ratio of triangular over elliptic flows (v 3 /v 2 ), have also been studied e.g. in [119,120] for UU collisions at RHIC energies. The work [119] shows that a model overestimation of the ε 3 /ε 2 ratio will imply a larger amount of the viscous damping needed in the subsequent theoretical hydrodynamic evolution to match the experimental UU data. A critical comparison of ICs derived from IP-Glasma and MCG models for light systems produced in pAu and dAu collisions at RHIC can be found in Ref. [170].

EXPERIMENTAL DEVELOPMENTS
On the experimental front, the LHC data has provided a wealth of new results that have helped to improve the extraction of relevant quantities from the Glauber approach. We review here two experimental aspects of importance for the determination of the reaction centrality in pA and AA collisions. The centrality determination -a proxy for the (arguably) most important parameter of Glauber models: the collision impact parameter-is found to be subject to stronger biases at the LHC than at lower c.m. energies.

Collision Centrality Estimates
As mentioned earlier, neither the impact parameter nor any derived Glauber quantity can be directly measured experimentally. Instead, average quantities are obtained within Glauber approaches for classes of events whose inclusive particle multiplicities and/or energy distributions can be reproduced by the corresponding calculation over a given b range. Since on average the impact parameter is monotonically related to the overall particle multiplicity, one typically measures multiplicity (or energy) distributions over a suitably large phase space. The mapping to calculated quantities then proceeds in intervals of centrality or centrality classes, which are obtained by binning the distribution in fractions of its total integral. Centrality is then typically defined as the percentile of the per-event multiplicity distribution dN/dM above M i relative to where M ap (< M i ) is the multiplicity value for which the fraction of total cross section was determined at the c ap point. The anchor point (AP) sets the absolute scale of the centrality. Clearly, one would like to achieve M ap close to zero, resulting in c ap close to 100%. However, because of trigger inefficiency and the increasing background contamination from ultraperipheral photonuclear collisions, experiments typically can only set the anchor point between 80 to 90% of the total hadronic cross section. An example is given in Fig. 7, which shows the sum of the amplitudes in the ALICE VZERO scintillators (at 2.8 < η < 5.1 and −1.7 < η < −3.7), called V0M centrality estimator, representing the uncorrected charged particle multiplicity distribution. The vertical lines indicate the typical centrality binning obtained from slicing the distribution in fractions of the total integral starting from 90% § , where smaller fractions refer to more central collisions.
To determine the anchor point, two approaches are typically used. The first one involves simulation of hadronic and e.m. processes, including a detailed description of the detector response, and hence gives direct access to the fraction of hadronic events below M AP . The second one involves MCG modeling together with a simple description for particle production at detector level, to simulate the uncorrected multiplicity distribution. The calculated distribution will describe the data down to the most peripheral events up to some point where background contamination and trigger inefficiency start to matter. The point where data and simulation start to separate can be used to set the AP.
To model multiplicity in heavy-ion collisions, one exploits that the majority of the initial-state NN collisions can essentially be treated like MB pp collisions, with a small perturbation from rarer hard interactions. The charged particle multiplicity in minimum-bias pp collisions at high energy can be described by a negative binomial distribution (NBD) [86], given by where µ is the mean and k is related to the width of the multiplicity distribution. Hence, the multiplicity for heavy-ion events can be approximated as a superposition of many NBDs, quickly approaching the Gaussian limit. A typical approach is to assume that the number of particle emitting sources can be described by the two-component approach, f · N part + (1 − f ) · N coll as in Eq. (4). A minimization procedure is then applied to the distribution of hadronic activity to determine the µ, k, and f parameters, with the values listed in Fig. 7 obtained with a χ 2 /N dof close to unity. The AP can be determined with ∼1% uncertainty (absolute). The centrality resolution is at the level of as few % in peripheral collisions and better than 1% in most central collisions.
Due to the finite kinematic acceptance, trigger inefficiency, and detector resolution, as well as the possible biases of the event selection, the details of the centrality determination differ between experiments, and even among collision systems within a given experiment. A short introduction with references to the approaches at RHIC is given in [40], while for more details on similar approaches at the LHC, see [171,172] for ALICE, [173] for CMS, and [174] for ATLAS. Alternative centrality estimators based on the transverse energy measured at forward rapidities (approx. 3 < |η| < 5 for CMS and ATLAS), as well as far-forward neutral energy in Zero Degree Calorimeters (ZDCs) along the beam line [171,175], have been also employed. The same methods can be applied to determine the centrality in pA or dAu collisions [176][177][178]. However, unlike for collisions of larger nuclei, the centrality determination is often subject to larger biases due to fluctuations in the categorization of events (see next section).

Collision Centrality Biases
As discussed in the previous section, medium effects on the production of perturbative QGP probes are in general quantified by the nuclear modification factor Eq. (16), defined as the ratio of the per-event yield measured in AA collisions over the same yield expected from an incoherent superposition of N coll binary pp collisions. However, event centrality classification involves selection of event samples for which the properties of the underlying binary NN collisions may deviate from those of unbiased pp collisions [83]. In this case, R AA (and R pA ) can deviate from unity even in the absence of nuclear effects. There are two main sources of selection biases. First, the number of hard processes is suppressed for increasingly peripheral AA collisions because of a simple geometrical bias: the probability for collisions increases proportional to b while the nuclear density decreases, leading to an increased probability for more-peripheral-than-average NN collisions. Second, the centrality selection, which relies mostly on measurements dominated by soft particle production -the relative weight of soft-to-hard contributions to the total energy distribution is given by the f ≈ 0.8 parameter in Fig. 7-biases the average multiplicity of individual NN collisions, and hence can affect the normalization of yields of collisions dominated by hard processes. As shown in Fig. 3 (left), hard scatterings are more probable in central NN collisions with large partonic overlap thereby leading to many simultaneous MPIs with lots of UE activity, and since hard processes are dominated by the production of jets that fragment (or heavy resonances that decay) into a large number of final hadrons, a peripheral AA event with a hard scattering often has a hadronic activity much larger than that typical of its centrality class. Since the AA (and pA) centrality determination is based on ordering the measured multiplicity or summed energy in the event, peripheral nuclear events with a hard scattering can thereby be wrongly assigned to a more central class.
which effectively leads to a reduction of N coll in peripheral collisions compared to T NN = δ(b NN ) because of their increased probability for less-central-than-average NN collisions. Standard MCG calculations, which by construction include the geometrical bias, however, do not use information about individual NN collisions. Even though the NN collisions are still modeled as occurring incoherently, the number of hard processes for a given centrality selections is not taken proportional to N coll , but to where N hard NN is the average number of hard scatterings in a NN collision for a given centrality selection and N hard NN is its unbiased average value. The mean number of hard scatterings per collision depends on b NN , and can be written as N hard (b NN ) = σ hard T NN (b NN ), where σ hard is the pQCD cross-section for 2 → 2 parton scatterings. Since the yield of hard and soft processes are correlated via their common b NN , the NN collisions can be biased towards lower or higher than average impact parameters, when ordering the measured multiplicity or transverse energy necessary for the centrality determination. This leads to a selection bias on N hard NN in addition to the inherent geometrical bias. Due to the strong dependence of σ hard on √ s, the selection bias is more relevant at LHC than at RHIC (and negligible at SPS) collision energies.
The relevance of the selection bias induced by the correlation between soft and hard particles is demonstrated in Fig. 8 that shows the nuclear modification factors for charged particle production, integrated above a large enough p T value (8 and 10 GeV), in pPb and PbPb collisions at √ s NN = 5.02 TeV compared to calculations using HG-PYTHIA [83]. The calculation, which uses the HIJING model to determine the distribution of N hard (b NN ) in a nuclear collision and PYTHIA 6.28 (Perugia 2011 tune) to generate the corresponding NN events, purposely does not include nuclear modification effects, unlike most of the models discussed in Sec. 4.2. As for the data, where the V0M estimator was used, the calculation determines the centrality using charged particles in the acceptance of the ALICE VZERO detectors. The calculation describes well the pPb and very peripheral PbPb data indicating that the strong apparent suppression from unity in this region results, in both cases, from the event selection. For more central PbPb collisions, where parton energy loss leads to the known, large suppression of particle production in PbPb compared to pp collisions, the R AA is not affected by such biases.
The hard-soft event selection bias is particularly important when fluctuations of the centrality estimator caused by b NN are of similar size as the dynamic range of N coll , as is the case in pA collisions, and strongly enhanced by trivial autocorrelations if the phase space for the measurement and event categorization are close-by or overlap [176][177][178]. This can be already deduced by computing the ratio between the average multiplicity of the centrality estimator and the average multiplicity per average ancestor of the Glauber fit as demonstrated in Fig. 8 of [177]. In contrast, centrality measurements based on zerodegree energy should not introduce any selection bias, while the geometric bias could still play a role. In the so-called hybrid method, described in [177], the pPb centrality selection is based on ZDC neutral energy in the Pb-going directions (slow neutrons), and N coll is determined from the measured charged particle multiplicity M according to N coll = N coll · M/ M , where N coll and M are, respectively, the centrality-averaged number of collisions and multiplicity. In case soft and hard particle yields are affected in the same way, the selection bias would cancel out in the nuclear modification factor.

SUMMARY
An outstanding topic in the physics of the strong interaction is the understanding of the thermodynamic and transport properties of hot and dense quark-gluon matter accessible to experimental study via highenergy collisions of nuclei. To correctly identify and interpret signals of collective partonic behavior in AA collisions, it is necessary to have a realistic extrapolation of the baseline hadron production properties of pp and pA collisions where, in principle, no quark-gluon plasma (QGP) is expected to be formed. In this work, we have reviewed the main developments and the state-of-the-art of the Glauber approach to describe multiple scatterings in proton and nuclear collisions, after 10 years of operation with pp, pPb, and PbPb collisions at the CERN LHC, as well as with deformed light-and heavy-ions at BNL RHIC. The Glauber model allows for arguably the simplest, most economical, and yet successful, understanding of collisions of extended hadronic objects based on an impact-parameter (b) superposition of independent elementary scatterings, each of which produces particles and thereby defines the local and global density of the precursor QGP. Key derived quantities in Glauber models include the nuclear overlap function T AA (b), number of participant nucleons N part (b), number of binary collisions N coll (b), transverse area A ⊥ (b), eccentricities ε n (b), average path-length L(b), of the strongly-interaction medium produced at different collision centralities, which are fundamental for the extraction of the QGP properties from the data.
The new LHC measurements, performed at 50 times larger c.m. energies than at previous nuclear collisions, and the latest precision RHIC data for a variety of colliding systems, have required to revisit and improve various ingredients of the Monte Carlo Glauber (MCG) simulations. Our review has first provided a new fit of the world measurements of the inclusive inelastic nucleon-nucleon cross sections (σ NN ), a key ingredient of MCG models. The inelastic hadronic pA and AA cross sections measured at the LHC are well reproduced by the corresponding MCG results derived using σ NN , a fact that confirms the overall validity of the Glauber model at the highest c.m. energies ever studied. Second, improved descriptions of the proton and nuclear density profiles, including subnucleonic degrees of freedom and neutron skin effects, and any associated sources of fluctuations and correlations, have been examined. Third, we have reviewed the main applications of the Glauber model for collider studies. The binary scaling prescription to quantitatively compare hard scattering cross sections in pp, pA, and AA collisions, has been validated by measurements of electroweak probes at the LHC whose yields are unaffected by final-state interactions in the QGP. The use of the Glauber formalism in MC event generators for heavyion physics, as well as to provide initial entropy-density profiles as input for relativistic hydrodynamics calculations, have also been discussed. The importance of a realistic description of the medium eccentricities to extract key transport properties, such as the QGP shear viscosity, from comparisons of elliptic and triangular flows measurements at RHIC and LHC to viscous hydrodynamics predictions has been highlighted. Last, the experimental procedures used and the inherent biases introduced by them, in the determination of the collision centrality from the data, which rely on the application of the MCG model, have been briefly discussed.
As an illustrative summary of our review, Fig. 9 shows the impact of different Glauber model ingredients and experimental biases, presented as a ratio as a function of centrality of the different elements over the default standard MCG calculation, for PbPb (left) and pPb (right) collisions at √ s NN = 5.5 and 8.8 TeV, respectively. First, the red long-dashed curves indicate the magnitude of the experimental shift introduced in measured nuclear modification factors when not properly accounting for event selection (multiplicity-and process-dependent) biases introduced by the centrality determination. These biases are significant in the most peripheral centrality classes and need to be carefully modeled and/or corrected for, specially when aiming at precision measurements at large impact parameters. The other curves indicate the ratio of the values of N coll obtained with modified ingredients with respect to standard MCG simulations. Inclusion of Glauber-Gribov fluctuations (via Eq. (10) with Ω = 1) , modified NN collision profiles (via Eq. (9) with ω = 0.4), or neutron skin effects, lead to few percent modifications of the N coll values in different centrality ranges, in principle within the assigned Glauber model systematic uncertainties [18]. An analytical calculation of N coll in the optical Glauber limit leads, however, to significant underestimations of the number of collisions for peripheral PbPb and pPb collisions.
The results of Fig. 9 emphasize large quantitative corrections needed in the Glauber model for the most peripheral AA and pA collisions. The availability of very large data samples of electroweak bosons at the LHC, with cross section measurements with few percent experimental uncertainties similar to or smaller than those of the Glauber model, opens up the possibility of using them to define an alternative experimental proxy for the nuclear overlap function. The ratio N V /(σ V NN N evt ), where σ V NN is the vector boson (V = ∑W ± , Z) production cross section in NN collisions (that can be estimated from pp measurements) and N V /N evt the per-event AA yields, has been suggested [180] as a data-driven T AA (b) proxy that would eliminate the need for Glauber modeling, and reduce corrections for centrality and event-selection effects, while cancelling uncertainties in the determination of quantities such as R AA . Such a proposal would require a higher level of theoretical accuracy (next-to-next-to-leading order) in the nuclear parton distribution functions, and in their centrality dependence, in order to fully exploit the high-precision σ V pp measurements.
All in all, the results summarized in this review show that, despite its simplicity (or, arguably, thanks to it), the Glauber model has stood for over 50 years as an indispensable and useful baseline approach to be able to quantitatively compare collisions of systems of varying size, from protons to uranium nuclei. Its continued exploitation to extract from the data key thermodynamics and transport properties of strongly interacting matter at the highest densities and temperature accessible in the laboratory, remains unchallenged for the years to come.

A Details on MCG simulations and parameters
This appendix provides a concise summary of technical details and parameters of MCG simulations. The MCG calculations for AA collisions are typically performed in two steps: (i) the nucleon positions in each colliding nucleus are determined for a given event, and (ii) the two nuclei are "collided" assuming that the nucleons travel in a straight-line along the beam axis (eikonal approximation). Conventional implementations can be found in the codes of Refs. [11,12,[14][15][16][17][18][19]. Among them, Refs. [16,17] provide a simple approach to include subnucleonic degrees of freedom.

Setup of nuclei
The position of each nucleon in the nucleus is determined according to a probability density function, where one commonly requires a minimum internucleon separation (d min = 0.4 fm) between the centers of the nucleons to mimic hard-core repulsion of nucleons inside the nucleus ¶ . The probability distribution for spherical nuclei is taken to be uniform in azimuthal and polar angles, and the radial probability function is taken from nuclear charge densities extracted in low-energy electron scattering experiments [20,21]. The nuclear charge density is usually parametrized by a Fermi distribution with three parameters where R is the nuclear radius, a is the skin depth, and w is a parameter to describe deviations from a spherical shape. Values used for various nuclei are listed in Table 2. The overall normalization given by the nucleon density (ρ 0 ) plays no role in the MCG simulation.   [20,21] for Eqs. (23) and (25). The suffix 'HN' indicates rescaled values to account for finite nucleon profiles as derived in Ref. [182]. All configurations can be found in the TGlauNucleus::Lookup function of the TGlauberMC code [183].
3. The proton taken from the Hulthén form, Eq. (24), with the neutron placed opposite to it.
For 3 H and 3 He, the configurations have been computed (and stored in a database) from Green's function MC calculations using the AV18 + UIX model interactions, which correctly sample the position of the three nucleons, including correlations, as in Ref. [117]. These, and the results of wavefunctionbased calculations for He-4, carbon, and oxygen derived in [118] are publicly available [183]. There are also rescaled values § § for the radius and skin depth of Cu ("CuHN"), Au ("AuHN") and Pb ("PbHN") to take into account the finite nucleon profile as derived in [182] and [119]. || Parameters R = 5.4 ± 0.1 fm and a = 0.61 +0.07 −0.09 fm for 132 Xe [181] are used, with the radius scaled down by 0.99 = (129/132) 1/3 and a reduced by 0.02 fm, to symmetrize the uncertainty and to approximate the smaller 129 Xe neutron skin. ** Parameters R = 5.32 fm and a = 0.57 fm for Sb(121.7) [20,21] are used, with the radius of Sb scaled up by 1.019 = (129/122) 1/3 using the average of A ≈ 122 between the two isotopes of natural Sb (A = 121 with 57.2%, and A = 123 with 42.8% abundances). The resulting parameters are consistent with those obtained from 132 Xe. † † These values are usually also used for 208 Pb, since the Bessel-Fourier coefficients for the two nuclei are similar [20,21]: the values R = 6.624 fm and a = 0.549 fm for 208 Pb ("Pb*") from [20] are consistent within uncertainties with the 207 Pb values. ‡ ‡ Separate (point-like) profiles for protons and neutrons, as described in [16]. § § Rescaled parameters that include the recentering effect, in addition, are given in [18].  Table 3: Nuclear charge density parameters for deformed nuclei from [20,21]. The suffix '2' for Si2, Cu2, Au2, and U2 refers to the use of Eq. (26). All configurations can be found in the TGlauNucleus::Lookup function of the TGlauberMC code [183].
For sulfur, and a few other nuclei (Si, Ca, Ni,...), a 3-parameter Gaussian-like form is typically used .

Collision of nuclei
After setting up the two nuclei, the impact parameter of the collision is chosen from dN/db ∝ b up to some large maximum b max with b max > R A + R B 20 fm. The centers of the nuclei are shifted to (−b/2, 0, 0) and (b/2, 0, 0), respectively. The reaction plane, defined by the impact parameter and the beam direction, is given by the xand z-axes, while the transverse plane is given by the xand y-axes. The longitudinal coordinate is irrelevant in the calculation as the nucleons move along a straight trajectory along the beam axis.
The "ball diameter" defined as D = σ NN /π defines the interaction strength of two nucleons: Two nucleons from different nuclei are usually assumed to collide if their relative transverse distance is less than the ball diameter, i.e. d < D. Variants to this approach use nucleon overlap profles that differ from the hard-sphere (or "black-disc") approach, as discussed in Sec. 3, see Eq. (9). If no nucleon-nucleon collision is registered for any pair of nucleons, then no nucleus-nucleus collision occurred. Counters for determination of the total (geometric) cross section are updated accordingly. ¶ ¶ Same parameters as for 129 Xe in Table 2 (App.) with β 2 and β 4 from Ref. [184]. *** Same parameters as for 129 Xe in Table 2 (App.) with β 2 = 0.18 ± 0.02 from interpolation between measured deformation parameters for the even-A Xe isotopes [185], and β 4 set to zero. † † † Implementation as in [186].