Multiplicity dependence of charged-particle production in pp, p–Pb, Xe–Xe and Pb–Pb collisions at the LHC

Multiplicity ( N ch ) distributions and transverse momentum ( p T ) spectra of inclusive primary charged particles in the kinematic range of | η | < 0 . 8 and 0 . 15 GeV / c < p T < 10 GeV / c are reported for pp, p–Pb, Xe–Xe and Pb–Pb collisions at centre-of-mass energies per nucleon pair ranging from √ s NN = 2 . 76 TeV up to 13 TeV. A sequential two-dimensional unfolding procedure is used to extract the correlation between the transverse momentum of primary charged particles and the charged-particle multiplicity of the corresponding collision. This correlation sharply characterises important features of the final state of a collision and, therefore, can be used as a stringent test of theoretical models. The multiplicity distributions as well as the mean and standard deviation derived from the p T spectra are compared to state-of-the-art model predictions. Providing these fundamental observables of bulk particle production consistently across a wide range of collision energies and system sizes can serve as an important input for tuning Monte Carlo event generators.

Multi-particle correlations measured in all collision systems, contain an imprint of the initial state of the colliding partners, characterised via their quark and gluon degrees of freedom, allowing the extraction of fundamental properties of the quark-gluon plasma [3,4,18].Hydrodynamic-like (final-state) collective flow is increasingly part of the modelling of small collision systems [19,20] in an interplay with initialstate phenomena (see reviews in Refs.[17,21]).Such collision systems are usually modelled in the framework of the colour glass condensate (CGC) [22], where multi-particle production proceeds via the decay of colour flux tubes that stretch between two colliding hadrons in the longitudinal direction and are coherent in the transverse direction over a range that is inversely proportional to a saturation scale Q s , which appears due to the non-linearity of parton evolution at high energies.
In the PYTHIA8 event generator [23], which describes a broad range of observables in pp collisions, the initial state is determined by parton distribution functions extracted from measurements.With increasing collision energy, the role of multi-parton interactions, i.e., when two or more distinct (hard) parton interactions occur within a pp collision, becomes more and more important [24].The respective colour strings may cut each other or reconnect, which leads to a redistribution of energy from particle production to transverse momentum, and, therefore, a smaller number of particles but with higher average p T .A phenomenon recently implemented in PYTHIA is the interaction between transversely-extended colour strings, exerting pressure on each other [25].This produces effects similar to those resulting from final-state collective dynamics, akin to that of a long-lived quark-gluon medium.Recently, PYTHIA8 was extended with the Angantyr model for heavy-ion collisions [26], which uses a Glauber-based initialstate modelling with Gribov colour fluctuations to determine the number of pp sub-collisions.In current PYTHIA8 implementations with the Angantyr model, there is no colour reconnection between individual pp sub-collisions but between the multiple partonic interactions of the pp sub-collisions.
A broad range of observables is also described successfully in the EPOS family Monte Carlo (MC) event generators [27].The initial state is realised in EPOS through a parton-based Gribov-Regge theory [27] which is a multiple scattering framework, recently augmented with the treatment of saturation effects [28].While in EPOS3 [29], a full hydrodynamic evolution is included for the final state, in EPOS LHC [30] a parameterised hydrodynamic component of a small volume with high density of thermalised matter is used together with a dilute region, i.e. a core plus corona implementation.In both PYTHIA8 and EPOS LHC models, the respective parameters are tuned using the Run 1 data at the LHC [30,31].
The mean transverse momentum, ⟨p T ⟩, of the charged-particle p T spectrum and its correlation with the charged-particle multiplicity N ch carry essential information on the underlying particle production mechanism.This has been studied by many experiments at hadron colliders in pp(p) covering centreof-mass energies from √ s = 31 GeV up to  as well as in Xe-Xe [42] and Pb-Pb [8] collisions at √ s NN = 5.44 TeV and 5.02 TeV, respectively.All experiments observed an increase of ⟨p T ⟩ with N ch in the central rapidity region, explained in the PYTHIA approach as a consequence of nontrivial colour reconnections (see discussion in Ref. [43]).While in the CGC approach ⟨p T ⟩ is a universal function of the ratio of the charged-particle multiplicity to the transverse area of the collision [44], in the EPOS LHC model it is determined by the collective expansion [30].For all collision systems, the ⟨p T ⟩-N ch correlation is a basic observable for tuning or calibrating the theoretical models [19], a simple observable which allows extracting fundamental properties of a deconfined quark-gluon medium [4].
As bulk production at the LHC is driven by a complex interplay of soft and hard QCD processes, the endeavour to find a consistent model description for particle production in all collision systems has not been concluded yet.At the LHC, a large amount of data was recorded in Run 1 and Run 2, covering different collision systems from pp to heavy-ion collisions at various centre-of-mass energies, which allows a detailed study of particle production over a wide range of charged-particle multiplicity.This Letter presents a measurement of the charged-particle multiplicity distributions and the corresponding transverse momentum spectra as a function of N ch in pp, p-Pb, Xe-Xe and Pb-Pb collisions at centreof-mass energies per nucleon pair ranging from √ s NN = 2.76 TeV up to 13 TeV.From these spectra, the mean ⟨p T ⟩ and standard deviation σ (p T ) = ⟨(p T − ⟨p T ⟩) 2 ⟩ within 0.15 GeV/c < p T < 10 GeV/c are extracted.The comprehensive set of measurements presented in this Letter can serve as a high-precision input for tuning phenomenological models towards the goal of understanding particle production in the non-perturbative regime of QCD and its transition to the perturbative regime.This Letter is structured as follows.Section 2 briefly describes the experimental setup.In Section 3, the data used for this analysis and a detailed description of the analysis procedure are presented.Results are discussed in Section 4, and a summary is given in Section 5.

Experimental setup
The measurements reported in this Letter were obtained with ALICE at the Large Hadron Collider.In the following, the detectors relevant for this work are discussed briefly.A more detailed description of the ALICE apparatus in its configurations of LHC Run 1 and Run 2 can be found in Refs.[45,46].The present analysis is based on tracking information from the Inner Tracking System (ITS) [47] and the Time Projection Chamber (TPC) [48].Both detectors are located within a large solenoidal magnet which provides a nominal field strength of 0.5 T for all of the data taking periods analysed in this work, except for the Xe-Xe data taking, where the magnetic field was reduced to 0.2 T in order to extend the kinematic acceptance of the detector to lower transverse momentum.The ITS is comprised of six cylindrical layers of silicon detectors with radii between 3.9 cm and 43 cm.Its two innermost layers are equipped with silicon pixel detectors (SPD), the two intermediate layers consist of silicon drift detectors (SDD), and the two outermost layers are made of double-sided silicon strip detectors (SSD).The large cylindrical TPC has an active radial range from about 85 cm to 250 cm and an overall length along the beam direction of about 500 cm.It covers the full azimuth in the pseudorapidity range |η| < 0.9 and provides track reconstruction with up to 159 space points along the trajectory of a charged particle as well as particle identification via the measurement of specific energy loss dE/dx.The V0 detector, which consists of two scintillator arrays covering the pseudorapidity ranges of 2.8 < η < 5.1 (V0A) and −3.7 < η < −1.7 (V0C), is used for triggering on hadronic collisions and multiplicity measurements at forward rapidity [8,49].Contamination from electromagnetic interactions in Pb-Pb and Xe-Xe collisions was strongly suppressed using signals from two neutron-Zero-Degree Calorimeters (ZDC), positioned on both sides of the interaction point at 114.0 m distance, see [46] for details.

Analysis procedure
This analysis aims to obtain the correlation between primary charged-particle p T spectra and their corresponding event multiplicities N ch , both defined consistently in the kinematic range |η| < 0.8 and 0.15 GeV/c < p T < 10 GeV/c.This kinematic selection ensures optimal momentum resolution and homogeneous tracking efficiency over the entire range.In addition, the multiplicity distributions for N ch > 0 events are reported.For each collision, the number of reconstructed tracks selected for the analysis provides an experimental measure (N meas ch ) for its multiplicity.Due to detector acceptance and reconstruction efficiency, this measured track multiplicity N meas ch is different from the actual number of primary charged particles (N ch ) produced in the kinematic region under study.Secondary particles from weak decays or from interactions with the detector material as well as tracks that are smeared into the acceptance (i.e. from |η| ≥ 0.8 and p T ≤ 0.15 GeV/c, p T ≥ 10 GeV/c) remaining in the sample of selected tracks further contribute to the measured p T spectra and consequently to the measured track multiplicity.The event-by-event fluctuations of reconstruction efficiency and contamination effects lead to a rather broad correlation between N ch and N meas ch that is shown in Fig. 1 for an example data set.As a result, the transverse momentum spectrum associated to a given N meas ch value carries contributions on particle production from a range of different N ch values.In addition, the finite resolution of the transverse momentum measurement itself results in a (small) smearing of the measured transverse momentum p meas T .The correct correlation between the collision-characteristic N ch value and its corresponding p T distribution can be extracted by unfolding the measured quantities.

Data sets and event selection
The data analysed in this work were collected between 2010 and 2018 during the LHC Run 1 and Run 2 data-taking periods.They comprise pp, p-Pb, Xe-Xe, and Pb-Pb collisions at centre-of-mass energies per nucleon pair ranging from √ s NN = 2.76 up to 13 TeV.Hadronic collisions were selected with two different minimum-bias (MB) interaction triggers.In Run 1, a single hit in either of V0A or V0C detectors or in the SPD was required (denoted as V0OR) in coincidence with a crossing of two particle bunches in the LHC at the nominal ALICE interaction point.For the Run 2 data taking period, a signal in both V0A and V0C was necessary to fulfil the MB trigger condition (denoted V0AND).Therefore, the former is also sensitive to single diffractive pp events while the latter almost exclusively selects non-single diffractive interactions.The offline event selection (for details see Ref. [8,49]) is optimised to reject beam-induced background and pileup collisions.Events with no reconstructed vertex and those with poor vertex resolution are rejected.Both the trigger efficiency and the vertex requirements affect mostly low multiplicity events.To ensure full pseudorapidity coverage of the tracking detectors (in particular by the inner ITS layers) and therefore avoid a possible asymmetry in the kinematic selection of the tracks, all collisions are required to have their reconstructed primary vertex located within |V meas z | < 10 cm along the beam direction with respect to the nominal interaction point.Table 1 shows an overview of the data sets and their corresponding number of events selected for this analysis.In order to have comparable results regardless of the trigger setup and event selections, Track selection A primary charged particle [50] is defined as a charged particle with a mean proper lifetime τ larger than 1 cm/c, which is either produced directly in the interaction or from decays of particles with τ smaller than 1 cm/c, excluding particles produced in interactions with the detector material.Charged-particle trajectories are reconstructed using both the ITS and TPC detectors.In order to select only tracks with good p T resolution for analysis, a minimal length in the active detector volume as well as a good agreement of the final track parameterisation with its comprising space points are required.The contamination of the track sample with weakly decaying particles, secondary particles from interactions of primary particles with the detector material and from pileup events is significantly reduced by selecting only tracks originating from a location close to the primary interaction vertex.In previous ALICE publications [8,42] those selection criteria were studied in great detail and optimised for best track quality and minimal contamination from secondary particles.

Particle-composition correction
The measured data are complemented by Monte Carlo simulations implementing a realistic GEANT3 [51] model of the ALICE detector and mimicking the experimental conditions present during the data taking.From these simulations, information about efficiency, secondary contamination, and smearing of N ch as well as p T is obtained.However, it was found in previous measurements [52,53] that the current state-of-the-art MC event generators do not perfectly reproduce the relative particle abundances, in particular for hyperons.Since the detector response, as well as the effect of the track selection, vary for the different particle species (e.g.due to different lifetimes), a purely MC-based correction for efficiency and feed-down contamination of inclusive charged particles would depend on the accuracy of the relative hadron abundances produced by the respective underlying event generator.To take this effect into account, the particle abundances from the event generator are re-weighted by means of a data-driven approach that was already employed in other ALICE analyses [8,42].This particle-composition correction utilises several ALICE measurements of identified (π, K, p, Λ) particle p T spectra as a function of multiplicity (in coarse intervals) for a range of collision systems [9,[53][54][55][56] as input and offers N ch -and p T -dependent correction factors for particle abundances.These data-driven adjustments for the generator bias result in a more accurate description of the detector performance and are applied prior to the unfolding corrections described in the following.

Event-level unfolding
The measurement of the raw charged-particle multiplicity distribution is influenced by several effects.In the experiment, some collisions that occurred within |V z | < 10 cm with respect to the nominal interaction vertex are not detected by the minimum-bias trigger or discarded by the subsequent event selection.Due to the experimental vertex-position resolution, a valid collision event may also be reconstructed outside of |V meas z | < 10 cm and therefore rejected in the analysis.On the other hand, the measured and selected sample of events may contain collisions without any primary charged particle produced within the kinematic range of interest (i.e.events with N ch = 0) or collisions with a true vertex position located outside |V z | < 10 cm.In addition, as discussed before, the value of the measured track multiplicity N meas ch itself is affected by tracking efficiency and track selection as well as contamination with secondaries and particles smeared into the kinematic acceptance, resulting in a broad correlation between the actual number of primary charged particles N ch and the measured track multiplicity N meas ch .Using the particle-composition corrected MC simulation, the measured track multiplicity distribution can be corrected for the efficiency, contamination, and smearing effects by means of an iterative unfolding procedure proposed by D'Agostini [57] and implemented in the RooUnfold [58] software package.
Starting with an initial assumption (prior) for the desired multiplicity distribution, which in this case is taken from the MC simulation, unfolding weights (posterior probabilities) are calculated by combining the prior with the detector response and the measured track multiplicity distribution according to Bayes' theorem.By again applying these posterior probabilities to the measured track multiplicity distribution, an updated and more accurate guess for the prior is calculated.This procedure is repeated at least three times and, in order to avoid overfitting, immediately stopped once the χ 2 /ndf between the resulting multiplicity distributions of two consecutive iterations drops below unity.In this context, the number of degrees of freedom refers to the number of data points in the respective distribution.The procedure is found to be very stable and the resulting unfolded spectrum after convergence is not sensitive to the choice of a prior.
Track-level unfolding While the one-dimensional multiplicity distributions are straightforward to unfold with the iterative D'Agostini method described above, extracting the correlation between p T spectra and the corresponding multiplicity poses a greater challenge.In principle, this two-dimensional (2D) deconvolution could be done in the same manner, i.e., by unfolding the distribution of (N meas ch , p meas T )-pairs and thereby extracting the corrected (N ch , p T )-distribution of primary charged particles.However, for the highly-granular measurement carried out here, this would imply a huge number of possible combinations and therefore, in practice, require an unreasonably large Monte Carlo event sample to sufficiently populate the smearing matrix that represents the detector response.Hence, in this analysis, a new approach was developed aiming to effectively achieve the 2D unfolding in a simpler way, by splitting it into multiple one-dimensional unfolding problems.Starting from the raw yield of charged-particle tracks as a function of measured track transverse momentum p meas T and measured track multiplicity N meas ch , this technique yields the fully corrected transverse momentum spectra of primary charged particles as a function of their corresponding primary charged particle multiplicity N ch .These fully corrected N ch -dependent p T spectra are then normalised to the number of N ch > 0 events obtained from the unfolded multiplicity distributions.In addition, the spectra are divided by the widths ∆p T and ∆N ch of the respective intervals chosen for analysis.In the present work, for pp, p-Pb and AA collisions with N ch ≤ 100, multiplicity intervals of ∆N ch = 1 are used, while for the rest of the range in AA collisions intervals of ∆N ch = 9 are chosen.This choice is driven by optimising the resolution of the unfolding procedure versus computing time.As the bulk of particles are produced at low transverse momentum, the p T intervals are set to have decreasing granularity from low to high p T .
In the experiment, event-level as well as track-level effects influence whether a charged particle with transverse momentum p T from an event of multiplicity N ch is detected with the measured properties p meas T and N meas ch .For the reasons described above, entire collision events, and in consequence all of their corresponding particles, can either be lost if the event is rejected or incorrectly selected and as a result contribute to the background contained in the measured track sample.Further, for an event which is selected for analysis, the p meas T spectrum as well as the measured track multiplicity N meas ch are affected by tracking efficiency, transverse momentum resolution, and contamination from secondaries or particles smeared into the kinematic acceptance of the measurement.While the event-level effects change the p Tintegrated detector response and are intrinsically multiplicity dependent, the track-level detector response mainly depends on the transverse momentum of the respective particle.However, there is still a small p T dependence of contamination and efficiency related to the event selection as the trigger may bias toward specific types of events (e.g., by selectively rejecting diffractive collisions which may have a transverse momentum shape different from that of non-diffractive events with the same multiplicity).On the other hand, also the tracking efficiency and contamination are (slightly) multiplicity dependent, as for instance the particle composition of the event changes with multiplicity, which is relevant in particular for AA collisions.
The basic idea behind the novel procedure employed in this analysis is that the multiplicity dimension (which is mostly affected by event-level effects) and the p T dimension (which is dominated by track-level effects) can be treated in two separate, sequential unfolding stages.In a first step, the N meas ch dependent raw track yield in each p meas T interval is unfolded separately using the event-level efficiencies and contamination as well as the (p T integrated) multiplicity smearing matrix of primary charged particles, which contains the probability for a primary charged particle from a collision with true multiplicity N ch to be found in an event with measured track multiplicity N meas ch .As a result of this deconvolution, the measured tracks are reassigned to the true multiplicities N ch and corrected for track losses related to the event selection and contamination from background events.In a second step, the p meas T dependent track yield is unfolded individually in each N ch interval using the corresponding p T dependent tracking efficiencies and p meas T dependent contamination, as well as the (multiplicity-integrated) transverse momentum smearing matrix of primary charged particles.Note that since the measured track distributions in each of the individually unfolded p meas T or N ch intervals are different, unique posterior probabilities (i.e., unfolding weights) are obtained in each of these cases.For all of the p meas T intervals, the p meas T integrated N ch distribution of measured primary charged particles taken from the MC simulation is used as the initial prior for the unfolding, while for all the N ch intervals the N ch integrated p T distribution of measured primary charged particles is used.
To validate the self-consistency of this sequential unfolding approach, it is applied to a MC sample which includes the transport of particles through the detector and the results are then compared with its underlying generator-level expectation.In Fig. 2 this closure test is shown for the mean and standard deviation of the unfolded N ch dependent transverse momentum spectra simulated for pp, p-Pb and Pb-Pb collisions at √ s NN = 5.02 TeV using the PYTHIA8, EPOS LHC and HIJING [59] event generators, respectively.Note that these are the particle-composition corrected simulations.The ratios in the bottom panels show a very good agreement between the unfolded and generated distributions over the whole range of multiplicities.The non-closure is mostly well below 1% and the remaining relative differences are used as an estimate for the systematic uncertainty of the method.The closure test was validated by cross-checking the largest data set, i.e. pp collisions at √ s = 13 TeV, with EPOS LHC as an alternative MC generator.
In the top panel of Fig. 3, the evolution of spectral shapes with multiplicity obtained with this unfolding procedure is shown for pp collisions at √ s = 5.02 TeV.This double-differential measurement not only allows characteristic properties of the spectra, e.g.⟨p T ⟩, and σ (p T ), to be extracted but also allows a direct comparison of the spectral shape produced in collisions with different multiplicities.The bottom left panel of Fig. 3 shows the ratio of multiplicity-dependent self normalised p T distributions P prim (p T |N ch ) to the multiplicity-integrated self normalised p T distribution P prim (p T ) of primary charged particles.A hardening of the spectra with multiplicity is apparent, which continuously increases with N ch , a trend observed earlier in coarser multiplicity intervals [41] and with different event selection methods [60].In the bottom right panel of Fig. 3 the self normalised multiplicity distribution of primary charged particles  for fixed transverse momentum P prim (N ch |p T ) divided by the p T -integrated N ch distribution P prim (N ch ) is shown as a function of N ch .As expected, it is evident that high p T particles are produced mostly in high-multiplicity events.

Systematic uncertainties
The accuracy of the corrections applied in this work depends on how well the measured track properties are reproduced by the MC simulations.The systematic uncertainties related to a possible disagreement were estimated by varying the track-selection criteria in reasonable ranges.A detailed list of those track quality criteria and their variations can be found in previous ALICE publications [8,42].In order to estimate the systematic uncertainty related to the particle-composition correction, the yields of identified particles are varied within their respective systematic uncertainties and the extrapolation of those spectra to p T = 0.15 GeV/c is performed with different functions.In addition, an uncertainty for the accuracy of the unfolding procedure is assigned that is quantified by the level of agreement between the unfolded simulated measurement and the expected distributions produced by the generator (MC closure test).For each variation, the fully corrected results are calculated and their deviation to the nominal result is considered an uncertainty.Therefore, any effect of the variations on the performance of the analysis procedure is included in the respective systematic uncertainty.To obtain the overall systematic uncertainties, all individual contributions are assumed to be fully uncorrelated and are added in quadrature.
The systematic uncertainties of multiplicity distributions are around 2-5% at low N ch , but increase towards higher multiplicities up to 10-20%, depending on the collision system and energy.While for pp and p-Pb collisions the track quality requirements are the most relevant contributions, in AA collisions the systematic uncertainty related to the particle-composition correction is the most prominent one.The systematic uncertainty of the mean transverse momentum of the unfolded spectra is largely dominated by the contributions from track-selection variations, yet the total systematic uncertainties in most of the reported N ch range are around 1%.At the lowest and highest multiplicities the contribution from the closure tests increases to up to 2%.The contribution from the particle-composition correction is around 0.5%.The systematic uncertainties of the standard deviation of the spectra σ (p T ) are dominated by the track quality requirements and are below 2% on average.At very low and very high N ch , the systematic uncertainties go up to 3-5% due to the MC closure test contribution.

Results and discussion
Multiplicity distributions as well as the mean and standard deviation of charged-particle p T spectra as a function of N ch are presented in comparison with model predictions.Table 2 summarises the mean and standard deviation of both the multiplicity distributions and N ch -integrated p T spectra for pp, p-Pb, Xe-Xe, and Pb-Pb collisions at the various centre-of-mass energies per nucleon pair.
The left panel of Fig. 4 shows the probability density of charged-particle multiplicity N ch for pp (top), p-Pb (middle), and AA (bottom) collisions at different energies per nucleon pair.For all collision systems, the distributions reach a maximum around N ch ≈ 2 and then fall steeply off over several orders of magnitude.In the pp and p-Pb systems, the slope of the decay with N ch decreases with increasing collision energy.This can be attributed to the larger momentum transfer in the initial hard scattering which results in larger multiplicities.The right panel of Fig. 4 presents the data after scaling the probability density and the charged-particle multiplicity with the average number of charged particles ⟨N ch ⟩ according to the Koba-Nielsen-Olesen (KNO) [61] scaling.Figure 5 shows the corresponding ratios of the KNO-scaled multiplicity distributions at various centre-of-mass energies per nucleon pair relative to √ s = 13 TeV, √ s NN = 8.16 TeV and √ s NN = 5.02 TeV for pp, p-Pb and Pb-Pb collisions, respectively.KNO scaling apparently holds in pp and AA collisions within 20%, and in p-Pb collisions within 10%.
In Fig. 6 the mean and standard deviation of the p T spectra are compared for pp, p-Pb and Pb-Pb collisions at the same centre-of-mass energy per nucleon pair of √ s NN = 5.02 TeV.All three collision systems have similar values at N ch = 1 and then coincide until Pb-Pb deviates at N ch ≈ 5 and p-Pb deviates at N ch ≈ 25 from the trend observed in pp.This observation is consistent with an earlier comparison of the ⟨p T ⟩-N ch correlation for the three systems at different centre-of-mass energies [62].Figure 7 shows the mean (left) and standard deviation (right) of the transverse momentum spectra as a function of the charged-particle multiplicity N ch for pp (top), p-Pb (middle), and AA (bottom) collisions at different centre-of-mass energies per nucleon pair.For all collision systems, a clear ordering of ⟨p T ⟩ as well as σ (p T ) with collision energy is observed, which can be attributed to the larger momentum transfers involved at higher √ s NN .For pp collisions at all centre-of-mass energies, the average transverse momentum increases monotonically with an almost linear trend up to N ch ≈ 16 and beyond that continues with an again almost linear dependence on N ch but reduced slope.In p-Pb collisions, a similar multiplicity dependence is observed up to N ch ≈ 25.At higher multiplicities, the increase in ⟨p T ⟩ is slower than in pp collisions.In both pp and p-Pb, σ (p T ) follows a similar trend as ⟨p T ⟩.On the other hand, for AA collisions one observes an increase of ⟨p T ⟩ with multiplicity up to about one third of the measured range, followed by a constant trend for the rest of the N ch range.The σ (p T ) increases for N ch ≲ 100 to a maximum and decreases afterwards.This is unique to large collision systems and is presumably a consequence of flow and jet quenching [8].The high N ch resolution of this measurement makes it possible to spot differences between the spectral evolution with multiplicity in Xe-Xe and Pb-Pb collisions at √ s NN = 5.44 TeV and √ s NN = 5.02 TeV, respectively.The observed difference in the trends might be a result of the slightly deformed Xe nuclei [63].In Fig. 8 both the mean (left) and standard deviation (right) of the p T spectra as a function of N ch are summarised for all data sets (top panels) and then shown as a function of relative multiplicity N ch /⟨N ch ⟩ (middle panels) as well as divided by their respective multiplicity-integrated values (bottom panels).In the latter scaling, the overall energy-dependent increase of average kinematic energy and number of produced particles are accounted for.As a result, the values for each collision system align almost perfectly for ⟨p T ⟩/⟨p T ⟩ incl and σ (p T )/σ (p T ) incl .The left panel of Fig. 9 shows the relative standard deviation of the spectra σ (p T )/⟨p T ⟩ as a function of N ch (left) and as a function of N ch /⟨N ch ⟩ (right).For pp collisions, this relative width of the p T spectra increases with multiplicity.The same trend is also observed for the larger collision systems.However, after around N ch ≈ 20 both for p-Pb and AA collisions, the standard deviation rises at the same rate as the mean, resulting in a flattening of the σ (p T )/⟨p T ⟩ ratio.After this plateau, the spectra in AA collisions become narrower relative to their mean values.The right panel of Fig. 9 shows the relative standard deviation of the spectra σ (p T )/⟨p T ⟩ as a function of the relative multiplicity N ch /⟨N ch ⟩.The plateau observed in p-Pb collisions starts at the relative multiplicity N ch /⟨N ch ⟩ ≈ 1, while the decrease observed in AA collision already begins at lower relative multiplicities of around N ch /⟨N ch ⟩ ≈ 0.2.
Figures 10 and 11 compare measured results for pp and p-Pb collisions with predictions from PYTHIA8 (solid lines) and EPOS LHC (dashed lines).Here, the PYTHIA8.306event generator is used with the Monash-2013 tune [31] for pp collisions and with the Angantyr model [26] for p-Pb collisions.The top left panel represents the ratio of models over measurements for the multiplicity distributions, the top right panel for the respective KNO-scaled multiplicity distributions and the bottom left and right panel show the ratios for ⟨p T ⟩ and σ (p T ), respectively.
In pp collisions, the overall shapes of the multiplicity distribution and KNO-scaled distribution shown in the upper panels of Fig. 10 are better described by EPOS LHC, while PYTHIA8 falls sharply off above N ch /⟨N ch ⟩ ≈ 4. Both models agree with the experimental distributions within 25% with larger deviations at highest multiplicities.For ⟨p T ⟩ and σ (p T ) shown in the bottom panels of Fig. 10, PYTHIA8 underpredicts the experimental data on ⟨p T ⟩ at the lowest values of N ch by up to 4%.The N ch dependent ⟨p T ⟩ values produced by PYTHIA8 increase faster than the measurements with an almost linear dependence up to N ch ≈ 20, after which the ratio shows a flat multiplicity dependence with an offset from unity varying from 0.5% at √ s = 5.02 TeV up to 4% at the highest centre-of-mass energy.EPOS LHC is further off at low multiplicities by up to 5% and increases slower than the measurements, underestimating them by up to 6% around N ch ≈ 9.At higher multiplicities, the increase is faster with a linearly rising ratio up to N ch ≈ 20 − 30, reaching a plateau which describes the measurements within ±2%.The experimental data on σ (p T ) is reproduced by both models within 10% at charged-particle multiplicities N ch > 10, with larger deviations at the lowest multiplicities.Results from model calculations in comparison with measurements from p-Pb collisions are shown in Fig. 11.PYTHIA8/Angantyr predicts the charged-particle multiplicity distribution within 30% (Fig. 11, top left) over the whole multiplicity range.EPOS LHC agrees within 20% for N ch < 70 but fails to describe the measurement at higher multiplicities.The KNO-scaled multiplicity distributions shown in the top right panel of Fig. 11 are described by both models within 20% up to a relative multiplicity of 2.5.Beyond that, both models exhibit increasing deviations from the measurement.PYTHIA8/Angantyr underpredicts ⟨p T ⟩ by about 5% at low multiplicities (Fig. 11, bottom left), N ch < 20, with the deviation increasing as a function of multiplicity, reaching about 25% at N ch = 110.This might result from the missing colour reconnection between the sub-collisions in the model.It is expected that high string density effects, as the recently-introduced shoving mechanism [64], will lead to an increase of ⟨p T ⟩ as a function of the multiplicity.EPOS LHC reproduces the ⟨p T ⟩ and σ (p T ) measurement within 10%.describes the pp measurements very well over the entire multiplicity range, even at the highest multiplicities.However, it significantly underpredicts the p-Pb measurements above N ch > 10, where the deviation increases with multiplicity.In Pb-Pb collisions, the ⟨p T ⟩ is systematically underpredicted by PYTHIA8/Angantyr by around 8% on average.Again, this points to the missing treatment of high string density effects which are not included in the PYTHIA8/Angantyr model, yet [26].The EPOS3 model overpredicts the ⟨p T ⟩ in all systems up to N ch ≈ 10 − 20, and underpredicts p-Pb and Pb-Pb measurements at higher multiplicities, less than PYTHIA8/Angantyr, but also cannot reproduce the ⟨p T ⟩ evolution with N ch .The hydrodynamical model calculations do not describe the measurements well, except for Pb-Pb collisions at the highest multiplicities of N ch > 1000.
In the CGC approach, the average transverse momentum is a universal function of the ratio of chargedparticle multiplicity and transverse area of a collision [44].The transverse area S T (N ch ) is derived from the interaction radius R( 3 dN g /dy) as a function of gluon multiplicity dN g /dy.This interaction radius was calculated within the CGC framework for pp collisions at a reference energy of √ s = 7 TeV = W 0,pp and for p-Pb collisions at √ s NN = 5.02 TeV = W 0,p-Pb [65].Parameterisations of these interaction radii proposed in Ref. [22] are used to calculate the interaction area.Following the arguments in Ref. [44], the ⟨p T ⟩ vs. N ch measurements presented in this Letter are scaled for each collision energy W = √ s NN for the respective collision system (pp, p-Pb) with a factor of (W /W 0 ) λ /(λ +2) .Here, the exponent λ characterises the saturation scale and was determined in Ref. [44] to be λ = 0.22 as the best fit to the transverse momentum distributions measured with ALICE.In order to approximate the gluon density corresponding to a measured final state multiplicity, a proportionality factor γ, defined by the equation dN g /dy = γN ch , is needed.Here, the naive value γ = 3/2 1  ∆η motivated by the ratio of the number of charged particles to all particles, as done in Ref. [22], is used.A weak dependence of the results on γ was observed.Figure 13 shows the ⟨p T ⟩ as a function of (W /W 0 ) λ /(λ +2) N ch /S T for pp and p-Pb collisions at various collision energies (left), and the ratio of those curves to the 13 TeV result (right), which has the highest reach in this scaling observable.The disagreement from the scaling is significant given the measurement's uncertainties, but still within about 10% over the entire range.At comparable values of the scaling variable, the ratio shows a distinct energy ordering, and all ratios exhibit a noticeable peak.In a more recent study [66], the determination of the exponent λ was revisited and it was found that the differential cross sections are better described when using λ ≈ 0.32 instead of λ = 0.22.However, with this updated value for the characteristic exponent λ , the geometrical scaling of the data presented in this Letter agrees only within ±15%.In general, the scaling results are found to be very sensitive to the value of λ and the best agreement is actually found for λ = 0, which effectively removes the energy scaling term (W /W 0 ) λ /(λ +2) proposed in Ref. [44].An approximate geometrical scaling of ⟨p T ⟩ was also observed in AA collisions as discussed in Ref. [67].

Summary and conclusions
A comprehensive study of inclusive charged-particle production at the LHC is presented, spanning a wide range of collision energies in pp, p-Pb, Xe-Xe and Pb-Pb collisions.Multiplicity distributions are compared across centre-of-mass energies for all collision systems, and, in addition, shown in the KNOscaling form.The KNO scaling is observed to hold within about 20%, not only for pp collisions, but also for p-Pb and AA collisions.Transverse momentum spectra are measured as a function of chargedparticle multiplicity in narrow N ch intervals.The mean and standard deviation of these p T spectra as a function of N ch are compared with PYTHIA8 (Angantyr), EPOS LHC, EPOS3, and hydrodynamical model predictions.For pp collisions, the spectral shape evolution with multiplicity is described fairly well by PYTHIA8 and EPOS LHC for all centre-of-mass energies, while EPOS3 and the hydrodynamical model fail to predict this observable.In general, for p-Pb and AA collisions there is a large tension with the data for all considered models, except for EPOS LHC which offers the best model prediction for p-Pb collisions.The geometric scaling of ⟨p T ⟩ proposed within the colour glass condensate framework is found to hold in first order, with deviations at the level of 10%.
Since the study of charged-particle production as a function of multiplicity plays a key role in understanding the properties of strongly-interacting matter created in collision systems of different sizes and energy densities, in the future, this rich high-precision set of multidimensional measurements can help to improve the theoretical modelling of the complex interplay of hard and soft QCD processes that govern particle production at LHC energies.

Figure 2 :
Figure 2: The N ch dependence of the mean (left panel) and standard deviation (right panel) of the p T distributions for (particle-composition corrected) Monte Carlo events in pp, p-Pb and Pb-Pb collisions at √ s NN = 5.02 TeV.Results propagated through a full GEANT model of ALICE and corrected with the sequential 2D unfolding (closed circles) procedure described in the text are compared with the generator-level (open squares) distributions and their ratios are shown in the bottom panels.

Figure 3 :
Figure3: Top panel: the correlation of primary charged particle p T spectra with multiplicity per N ch > 0 event for pp collisions at √ s = 5.02 TeV.Bottom panels: the corresponding relative change of p T (left) and N ch (right) distributions with respect to the inclusive ones.In the left panel, each of the curves represents a single N ch value, ranging from N ch = 1 (blue) to N ch = 55 (red).In the right panel, the colours represent the p T intervals used in this analysis from the lowest in blue to the highest one in red.

Figure 4 :
Figure 4: Probability density of charged-particle multiplicity N ch (left) and the corresponding KNO-scaled distributions (right) for pp (top), p-Pb (middle), and AA (bottom) collisions at different centre-of-mass energies per nucleon pair.Statistical and systematic uncertainties are shown as bars and semi-transparent bands, respectively.

Figure 5 :
Figure 5: Ratios of the KNO-scaled multiplicity distributions at various centre-of-mass energies per nucleon pair relative to √ s = 13 TeV for pp collisions (top panel) and relative to √ s NN = 8.16 TeV and √ s NN = 5.02 TeV for p-Pb and Pb-Pb collisions, respectively (left and right bottom panels).Statistical and systematic uncertainties are shown as bars and semi-transparent bands, respectively.

Figure 6 :
Figure 6: Mean (left) and standard deviation (right) of the charged-particle transverse momentum spectra as a function of the charged-particle multiplicity for pp, p-Pb, and Pb-Pb collisions at a centre-of-mass energy per nucleon pair of √ s NN = 5.02 TeV.Statistical and systematic uncertainties are shown as bars and semi-transparent bands, respectively.

Figure 7 :
Figure 7: Mean (left) and standard deviation (right) of the charged-particle transverse momentum spectra as a function of the charged-particle multiplicity for pp (top), p-Pb (middle), and AA (bottom) collisions at different centre-of-mass energies per nucleon pair.Statistical and systematic uncertainties are shown as bars and semitransparent bands, respectively.

Figure 8 :
Figure 8: Mean (left) and standard deviation (right) of the charged-particle transverse momentum spectra as a function of the charged-particle multiplicity (top) and relative multiplicity N ch /⟨N ch ⟩ (middle, bottom) for pp, p-Pb, Xe-Xe and Pb-Pb collisions at various centre-of-mass energies per nucleon pair.The bottom panels show both quantities relative to their multiplicity-inclusive value.Statistical and systematic uncertainties are shown as bars and semi-transparent bands, respectively.

Figure 9 :
Figure 9: Relative standard deviation of the charged-particle transverse momentum spectra as a function of the absolute (left) and relative (right) charged-particle multiplicity for pp, p-Pb, Xe-Xe and Pb-Pb collisions at various centre-of-mass energies per nucleon pair.Statistical and systematic uncertainties are shown as bars and semitransparent bands, respectively.

Figure 10 :
Figure 10: Ratio of model predictions to data for pp collisions at various energies.The upper panels show it for the multiplicity distributions (left) and their KNO-scaling form (right), the bottom panels represent ⟨p T ⟩ (left) and σ (p T ) (right).The semi-transparent bands indicate the relative systematic uncertainties of the data.

Figure 12
Figure 12 shows a comparison of the measured ⟨p T ⟩ as a function of N ch for pp, p-Pb and Pb-Pb collisions at the same centre-of-mass energy per nucleon pair, √ s NN = 5.02 TeV, with results from three different model calculations: PYTHIA8 (left panel; Angantyr for p-Pb and Pb-Pb), EPOS3 (middle panel), and hydrodynamics with CGC initial conditions [20] (right panel).As shown before, PYTHIA8

Figure 11 :
Figure 11: Ratio of model predictions to data for p-Pb collisions at various energies.The upper panels show it for the multiplicity distributions (left) and their KNO-scaling form (right), the bottom panels represent ⟨p T ⟩ (left) and σ (p T ) (right).The semi-transparent bands indicate the relative systematic uncertainties of the data.

Figure 12 :
Figure 12: Comparisons of ⟨p T ⟩ as a function of N ch at √ s NN = 5.02 TeV for pp, p-Pb, and Pb-Pb collisions to three different model predictions.Statistical and systematic uncertainties are shown as bars and semi-transparent bands, respectively.

Figure 13 :
Figure 13: Average transverse momentum ⟨p T ⟩ as a function of the scaling variable (W /W 0 ) λ /(λ +2) N ch /S T [44] for pp and p-Pb collisions at various energies (left) and the ratio of all data sets to that in pp collisions at 13 TeV (right).The reference energy W 0 corresponds to √ s = 7 TeV for pp and √ s NN = 5.02 TeV for p-Pb collisions.Statistical and systematic uncertainties are shown as bars and semi-transparent bands, respectively.

Table 1 :
Overview of the analysed data sets.The definitions of the two MB triggers are explained in the text.

Table 2 :
Global characteristics of the analysed data sets with corresponding systematic uncertainties.Both the N ch and the p T spectra are consistently defined in the kinematic range |η| < 0.8 and 0.15 GeV/c < p T < 10 GeV/c and only events with N ch > 0 are considered.