Skewness and kurtosis of mean transverse momentum fluctuations at the LHC energies

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Skewness and kurtosis of mean transverse momentum fluctuations at the LHC energies © 2024 The Author(s). Published by Elsevier B.V. Funded by SCOAP³ Published version ALICE Collaboration ALICE Collaboration. (2024). Skewness and kurtosis of mean transverse momentum fluctuations at the LHC energies. Physics Letters B, 850, Article 138541. https://doi.org/10.1016/j.physletb.2024.138541 2024 Phys. Lett. B 850 (2024) 138541 Available online 23 February 2024 0370-2693/© 2024 The Author(s). Published by Elsevier B.V. Funded by SCOAP³. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Physics Letters B journal homepage: www.elsevier.com/locate/physletb Letter Skewness and kurtosis of mean transverse momentum fluctuations at the LHC energies .ALICE Collaboration⋆ A R T I C L E I N F O A B S T R A C T Editor: M. Doser Dataset link: https:// www .hepdata .net /record /ins2692420 The first measurements of skewness and kurtosis of mean transverse momentum (⟨𝑝T⟩) fluctuations are reported in Pb–Pb collisions at √𝑠NN =5.02 TeV, Xe–Xe collisions at √𝑠NN =5.44 TeV and pp collisions at √𝑠=5.02 TeV using the ALICE detector. The measurements are carried out as a function of system size ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5, using charged particles with transverse momentum (𝑝T) and pseudorapidity (𝜂), in the range 0.2 <𝑝 T<3.0GeV/cand |𝜂| <0.8, respectively. In Pb–Pb and Xe–Xe collisions, positive skewness is observed in the fluctuations of ⟨𝑝T⟩for all centralities, which is significantly larger than what would be expected in the scenario of independent particle emission. This positive skewness is considered a crucial consequence of the hydrodynamic evolution of the hot and dense nuclear matter created in heavy-ion collisions. Furthermore, similar observations of positive skewness for minimum bias pp collisions are also reported here. Kurtosis of ⟨𝑝T⟩fluctuations is found to be in good agreement with the kurtosis of Gaussian distribution, for most central Pb–Pb collisions. Hydrodynamic model calculations with MUSIC using Monte Carlo Glauber initial conditions are able to explain the measurements of both skewness and kurtosis qualitatively from semicentral to central collisions in Pb–Pb system. Color reconnection mechanism in PYTHIA8 model seems to play a pivotal role in capturing the qualitative behavior of the same measurements in pp collisions. 1. Introduction The properties of the hot and dense nuclear matter created in heavyion collisions at relativistic energies can be studied using event-by-event fluctuations of different quantities like multiplicity, net-charge, mean transverse momentum (⟨𝑝T⟩), etc. [1–4]. The analysis of event-by-event fluctuations of these variables offers valuable means to probe the dynamical fluctuations that originate from the production of a quark– gluon plasma (QGP) [5]phase during heavy-ion collisions. Fluctuations in the thermodynamic quantity of temperature, which are associated with the phase transition in the quantum chromodynamics (QCD) phase diagram, can manifest themselves in the fluctuations of the ⟨𝑝T⟩of the final-state particles [6]. A non-monotonic behavior of ⟨𝑝T⟩fluctuations as a function of centrality or incident energy was suggested as one of the possible signals of the QGP [3]. However, ⟨𝑝T⟩fluctuations are also affected by non-thermodynamic variations in the initial geometry of the collision that include fluctuations in the initial size, shape and orientation of the colliding nuclei, and the fluctuating number of nucleons participating in the collision. The ⟨𝑝T⟩fluctuations measured at RHIC did not show any beam energy dependence, and the non-monotonic behavior with centrality was also not observed [7]. The previous measurements of ⟨𝑝T⟩fluctuations in Pb–Pb collisions at center-of-mass ⋆E-mail address: alice -publications @cern .ch. energy per nucleon pair, √𝑠NN =2.76 TeV in ALICE [8] suggested a connection of the observed fluctuations of the ⟨𝑝T⟩to the fluctuations in the initial state of the collision. Qualitatively, the data obtained from model calculations using the string melting approach in A Multi-Phase Transport (AMPT) model [9], where partons rescatter and recombine through a hadronic coalescence scheme, exhibit agreement with the observed results [8]. The string melting version was introduced in the model [10]to address and account for the effects of flow originating from the entire partonic system in the overlap volume of heavy-ion collisions, as opposed to considering only the contributions from minijet partons in the default version. By incorporating the string melting approach, the model was able to simulate the collective behavior of the partonic system during the early stages of the collision, leading to a more accurate description of the experimental observations of elliptic flow. The ⟨𝑝T⟩fluctuations were also calculated in the color glass condensate (CGC) [11,12] formulation, where they have been related to initial spatial fluctuations of glasma flux tube via their coupling to a collective flow field. The comparison of these calculations with data showed a good agreement in the semicentral and central collisions [13]. The space–time evolution of the QGP phase produced in the relativistic heavy-ion collisions is well described by the relativistic viscous hydrodynamics [14–16]. In Ref. [17], it has been proposed that the https://doi.org/10.1016/j.physletb.2024.138541 Received 22 September 2023; Received in revised form 18 January 2024; Accepted 19 February 2024 Physics Letters B 850 (2024) 138541 2 ALICE Collaboration skewness of ⟨𝑝T⟩fluctuations can serve as an essential probe of the hydrodynamic behavior of the system created in heavy-ion collisions. The ⟨𝑝T⟩of the particles emitted at freeze-out was found to be correlated to the initial energy of the system at the beginning of the hydrodynamic evolution, instead of the energy of the system at freeze-out [18]. The fluctuations of ⟨𝑝T⟩are therefore found to be related to the fluctuations of initial energy density in an effective hydrodynamic description [19]. It was shown that the skewness of ⟨𝑝T⟩fluctuations is driven by the skewness of the initial energy density fluctuations, which implies that the ⟨𝑝T⟩fluctuations arise from the same collective dynamics in the QGP phase that give rise to anisotropic flow. Furthermore, employing initial conditions from the TRENTo model [20]and evolving them with the V-USPHYDRO viscous hydrodynamic [21] simulations predicted positive skewness of ⟨𝑝T⟩fluctuations surpassing expectations in independent particle emission scenarios. Additional information regarding these simulations can be found in Ref. [17]. The fluctuations in ⟨𝑝T⟩of charged particles can be influenced by various physical effects, such as collective behavior of the system formed in the collisions, fluctuations in the number of participating nucleons, or the presence of jets and resonance decays. Examining the higher-order terms of ⟨𝑝T⟩fluctuations will enable us to delve deeper into the intricate mechanisms underlying the observed fluctuations and acquire valuable insights. In this article, the first experimental study of skewness and kurtosis of ⟨𝑝T⟩fluctuations that represent the thirdand fourth-order fluctuations of ⟨𝑝T⟩are reported at the Large Hadron Collider (LHC) energies. The observables used in this analysis are introduced in Sec. 2. Here, the multiparticle 𝑝Tcorrelators used to study the ⟨𝑝T⟩fluctuations are defined. A brief description of the subsystems of ALICE detector relevant to this analysis is given in Sec. 3. The analysis technique and the method of estimating statistical and systematic uncertainties for the measurements are explained in Sec. 4. The skewness and the kurtosis of ⟨𝑝T⟩fluctuations as a function of system size and the interpretation of results using theoretical models are presented in Sec. 5. The major findings of the analysis are summarized in Sec. 6. 2. Observables In this analysis, the fluctuations in event-by-event ⟨𝑝T⟩of charged particles are investigated using multiparticle 𝑝Tcorrelators. The eventby-event ⟨𝑝T⟩is defined as ⟨𝑝T⟩=∑𝑁ch 𝑖=1 𝑝T,𝑖 𝑁ch ,(1) where 𝑝T,𝑖 is the transverse momentum of the ith particle and 𝑁ch is the total number of charged particles in the event. Alternatively, one can employ the standard moment method for event-by-event analysis of ⟨𝑝T⟩fluctuations. The moment method calculates various order moments of the ⟨𝑝T⟩distribution, providing a comprehensive evaluation of the total fluctuation accounting for both the statistical and dynamical (non-statistical) parts. The advantage of employing multiparticle 𝑝Tcorrelators [22,23]is that they yield zero values for events with randomly sampled particles, thereby effectively isolating the non-statistical fluctuations of interest. The expressions for the two-particle 𝑝Tcorrelator (⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 ⟩) and three-particle 𝑝Tcorrelator (⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 Δ𝑝T,𝑘⟩)[17] are denoted by Eqs. (3) and (4), respectively, where 𝑄𝑛is defined as 𝑄𝑛= 𝑁ch ∑ 𝑖=1 𝑝𝑛 T,𝑖 .(2) Constraints on the indices in summations in the definition of correlators ensure that self-correlations are eliminated. The algebraic expressions of the correlators in terms of 𝑄𝑛’s where 𝑛 =1, 2, ... were derived to account for the practical limitations of calculating multiparticle 𝑝Tcorrelators with nested loops, especially when dealing with events characterized by large multiplicities. Since, the values of 𝑄𝑛can be obtained with Eq. (2)employing a single loop, it becomes feasible to directly utilize the later expressions with 𝑄𝑛s in Eqs. (3) and (4)for the calculation of twoand three-particle 𝑝Tcorrelators, respectively. On similar lines, the four-particle 𝑝Tcorrelator (⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 Δ𝑝T,𝑘Δ𝑝T,𝑙⟩) is derived as shown in Eq. (5). These two-, three-, and four-particle correlators are related to the variance, skewness, and kurtosis of the ⟨𝑝T⟩distribution, respectively. ⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 ⟩=⟨∑𝑁ch 𝑖,𝑗, 𝑖≠𝑗 (𝑝T,𝑖 −⟨⟨𝑝T⟩⟩)(𝑝T,𝑗 −⟨⟨𝑝T⟩⟩) 𝑁ch(𝑁ch −1) ⟩ev =⟨𝑄2 1−𝑄2 𝑁ch(𝑁ch −1)⟩ev −⟨𝑄1 𝑁ch ⟩2 ev ,(3) ⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 Δ𝑝T,𝑘⟩ =⟨∑𝑁ch 𝑖,𝑗,𝑘, 𝑖≠𝑗≠𝑘 (𝑝T,𝑖 −⟨⟨𝑝T⟩⟩)(𝑝T,𝑗 −⟨⟨𝑝T⟩⟩)(𝑝T,𝑘 −⟨⟨𝑝T⟩⟩) 𝑁ch(𝑁ch −1)(𝑁ch −2) ⟩ev =⟨𝑄3 1−3𝑄2𝑄1+2𝑄3 𝑁ch(𝑁ch −1)(𝑁ch −2)⟩ev −3⟨𝑄2 1−𝑄2 𝑁ch(𝑁ch −1)⟩ev ⟨𝑄1 𝑁ch ⟩ev +2⟨𝑄1 𝑁ch ⟩3 ev ,(4) ⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 Δ𝑝T,𝑘Δ𝑝T,𝑙⟩ =⟨∑𝑁ch 𝑖,𝑗,𝑘,𝑙, 𝑖≠𝑗≠𝑘≠𝑙 (𝑝T,𝑖 −⟨⟨𝑝T⟩⟩)(𝑝T,𝑗 −⟨⟨𝑝T⟩⟩)(𝑝T,𝑘 −⟨⟨𝑝T⟩⟩)(𝑝T,𝑙 −⟨⟨𝑝T⟩⟩) 𝑁ch(𝑁ch −1)(𝑁ch −2)(𝑁ch −3) ⟩ev =⟨𝑄4 1−6𝑄4+8𝑄1𝑄3−6𝑄2 1𝑄2+3𝑄2 2 𝑁ch(𝑁ch −1)(𝑁ch −2)(𝑁ch −3) ⟩ev −4⟨𝑄3 1−3𝑄2𝑄1+2𝑄3 𝑁ch(𝑁ch −1)(𝑁ch −2)⟩ev ⟨𝑄1 𝑁ch ⟩ev +6⟨𝑄2 1−𝑄2 𝑁ch(𝑁ch −1)⟩ev ⟨𝑄1 𝑁ch ⟩2 ev −3⟨𝑄1 𝑁ch ⟩4 ev .(5) In the above equations, angular brackets ⟨…⟩ev denote an average over all events and ⟨⟨𝑝T⟩⟩ =⟨𝑄1∕𝑁ch⟩ev is the event-by-event ⟨𝑝T⟩ averaged over all events. The standardized skewness (𝛾⟨𝑝T⟩) and the intensive skewness (Γ⟨𝑝T⟩) are two measures of the skewness of the ⟨𝑝T⟩ distribution, constructed using the two-particle and three-particle 𝑝T correlators [17], given by 𝛾⟨𝑝T⟩=⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 Δ𝑝T,𝑘⟩ ⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 ⟩3∕2 ,(6) and Γ⟨𝑝T⟩=⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 Δ𝑝T,𝑘⟩⟨⟨𝑝T⟩⟩ ⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 ⟩2.(7) The kurtosis 𝜅⟨𝑝T⟩of ⟨𝑝T⟩distribution is defined as 𝜅⟨𝑝T⟩=⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 Δ𝑝T,𝑘Δ𝑝T,𝑙⟩ ⟨Δ𝑝T,𝑖Δ𝑝T,𝑗 ⟩2.(8) The article focuses on analyzing higher order fluctuations of ⟨𝑝T⟩, particularly the third and fourth order. Three key observables, namely the standardized skewness, intensive skewness, and kurtosis of the ⟨𝑝T⟩ distribution (defined in Eqs. (6), (7), and (8), respectively), are investigated as a function of the system size in Pb–Pb collisions at √𝑠NN = 5.02 TeV, Xe–Xe collisions at √𝑠NN = 5.44 TeV, and pp collisions at √𝑠=5.02 TeV. Physics Letters B 850 (2024) 138541 3 ALICE Collaboration 3. Experimental setup The measurements reported in the article are obtained using the data recorded by the ALICE detector at the LHC. A detailed description of the ALICE detector and its performance can be found in Refs. [24,25]. The primary sub-detectors relevant to this analysis are the Time Projection Chamber (TPC) [26], the Inner Tracking System (ITS) [25], and the V0 detector [27]. The TPC and ITS are used for tracking and reconstructing the primary vertex, while the V0 detector is used for triggering and the default centrality estimation. The V0 detector consists of two scintillator arrays, V0A and V0C, located on both sides of the interaction point, covering the pseudorapidity (𝜂) intervals 2.8 <𝜂<5.1and −3.7 <𝜂<−1.7, respectively [27]. The data analyzed here are obtained from Pb–Pb collisions at √𝑠NN =5.02 TeV, Xe–Xe collisions at √𝑠NN = 5.44 TeV, and pp collisions at √𝑠=5.02 TeV. These datasets were recorded in 2015 (pp and Pb–Pb) and 2017 (Xe–Xe) during LHC Run 2. A minimum bias (MB) trigger condition is applied to select collision events that requires at least one hit in both the V0A and the V0C detectors. Events that pass the MB trigger criteria and have a reconstructed primary vertex position within 10 cm along the beam axis to the nominal interaction point are selected in the analysis. Using information from multiple detectors as described in Ref. [25], events with more than one reconstructed primary interaction vertex (pile-up events) are rejected. In total, 84 million Pb–Pb collisions, 1.2 million Xe–Xe collisions, and 95 million pp collisions pass the above mentioned criteria and are used for the analysis. The charged-particle tracks reconstructed using the ITS and the TPC in the ALICE central barrel are selected in the kinematic range 0.2 <𝑝 T<3.0GeV∕𝑐and |𝜂| <0.8, where 𝑝Tis the track momentum in the plane transverse to the beam axis. The accepted tracks exhibit approximately uniform azimuthal acceptance in this region. The tracks that have at least 70 out of a maximum possible 159 reconstructed space points in the TPC, and at least one hit in the two innermost layers of the ITS (ITS has in total six layers) are selected. This selection assures a resolution better than 300 μm [25]on the distance-of-closest-approach (DCA) to the primary vertex in the plane perpendicular (DCA𝑥𝑦) and parallel (DCA𝑧) to the beam axis for the selected tracks. In order to suppress the contamination from secondary particles, the DCA of the tracks to the primary vertex must be within 1 cm in the longitudinal direction and 0.1 cm in the transverse plane. Moreover, the 𝜒2per degree of freedom of the track fit to the space points in the TPC and the ITS are required to be less than 4 (2.5 for Pb–Pb) and 36, respectively. 4. Data analysis The mean transverse momentum ⟨𝑝T⟩distributions obtained for three different centrality classes in Pb–Pb collisions at √𝑠NN =5.02 TeV are shown in Fig. 1. The distributions are not corrected for detector inefficiencies. The centrality classes are formed by splitting the events based on the measured amplitude distribution in the V0A and V0C detectors as described in Refs. [28,29]. The mean of the ⟨𝑝T⟩distributions is found to increase whereas the width of the distribution decreases from peripheral to central collisions. The larger width of the distribution for peripheral (lower multiplicity) collisions indicates larger fluctuations compared to the central (higher multiplicity) collisions. Similar behavior is observed in Xe–Xe collisions at √𝑠NN =5.44 TeV and pp collisions at √𝑠=5.02 TeV. Before starting the analysis with the ⟨𝑝T⟩correlators presented in Sec. 2, it is important to demonstrate that these fluctuations are not a trivial consequence of fluctuating 𝑁ch in a given centrality class. The definition of ⟨𝑝T⟩in Eq. (1) clearly shows that the varying multiplicity in different events of the same centrality class can affect the ⟨𝑝T⟩distributions and consequently its fluctuations. To examine whether these fluctuations persist after removing the stochastic effects of the multiplicity, a check is performed. For each centrality class, the minimum Fig. 1. Efficiency uncorrected mean transverse momentum distributions for different centrality classes in Pb–Pb collisions at √sNN =5.02 TeV. number of tracks per event in the kinematic acceptance of the analysis, 𝑁min ch is determined. In each event of the centrality class, the 𝑁min ch number of tracks are then selected randomly to calculate the genuine ⟨𝑝T⟩distribution free from biases of multiplicity fluctuations. The left panel of Fig. 2shows the ⟨𝑝T⟩distribution for the 0–5% central events in Pb–Pb collisions at √𝑠NN =5.02 TeV. The distribution in red markers shows the original event-by-event ⟨𝑝T⟩distribution in which 𝑁ch fluctuates whereas the distribution shown with black markers is the ⟨𝑝T⟩ distribution in which 𝑁ch is fixed to 𝑁min ch in each event of the centrality class. A Gaussian function is used to fit the two distributions and the ratio of data to fit is shown in the bottom panel of Fig. 2. The ratio shows that the data points are above the fit in the right-side of the distribution, and below the fit in the left-side, indicating that the ⟨𝑝T⟩distributions are positively skewed. The ⟨𝑝T⟩distributions exhibit substantial deviations from a Gaussian distribution even under the condition of a fixed 𝑁ch, thereby reflecting the analogous traits observed in the original distribution with fluctuating 𝑁ch. The middle and right panels of Fig. 2 visually exemplify the consistent manifestation of these characteristics within the distributions of ⟨𝑝T⟩for both semicentral and peripheral collisions. It is concluded that the positive skewness in the event-by-event ⟨𝑝T⟩distributions is not a trivial consequence of event-by-event statistical fluctuations of 𝑁ch, and therefore new and non-trivial information can be extracted from higher-order moments of ⟨𝑝T⟩fluctuations. The ⟨𝑝T⟩correlators described in Eqs. (3), (4), and (5)are calculated using charged-particle tracks in different centrality (multiplicity) classes. For Pb–Pb collisions, the observables are analyzed in 18 equal-width centrality classes ranging between zero to ninety percent, i.e., 0–5%, 5–10%,..., 85–90%. For Xe–Xe collisions, the total number of events are divided into four centrality classes, 0–20%, 20–40%, 40–60%, and 60–80%. The centrality classes are taken wider for Xe– Xe collisions because of the limited data sample. In pp collisions, the events are classified into ten multiplicity classes, 0–1%, 1–5%, 5–10%, 10–15%, 15–20%, 20–30%, 30–40%, 40–50%, 50–70%, and 70–100%. The method of defining these multiplicity classes can be found in Ref. [30]. The ⟨𝑝T⟩correlators are calculated in the unit multiplicity bins of a given centrality class, and then merged using centrality bin width correction formula [31]. In this way, both the effects of multiplicity fluctuations, and the necessity of using non-trivial multiplicity weight to determine all-event averages from single-event averages, are eliminated. The statistical uncertainties are evaluated using the bootstrap method. By repeatedly resampling the dataset and analyzing each resample, the bootstrap method provides a robust and computationally feasible approach for estimating the uncertainties associated with the observed results. It allows for a comprehensive assessment of the vari- Physics Letters B 850 (2024) 138541 4 ALICE Collaboration Fig. 2. Comparison of event-by-event mean transverse momentum ⟨𝑝T⟩distribution from actual analysis and modified analysis in which 𝑁ch in each event is fixed to 𝑁min ch (𝑁min ch is the minimum number of charged particles in a given centrality class) and 𝑁min ch number of particles are selected in each event randomly for the 0–5% (left), 30–40% (middle), and 50–60% (right) centrality classes. ability in the data, taking into account the inherent fluctuations and dependencies present in the measurements. The systematic uncertainties on the observables are estimated separately for each collision system by varying event and track selection criteria. The uncertainties related to the event selection include the variation of the accepted vertex position along the beam direction and pileup cut where applicable. The uncertainties due to track selection include the variation of the selection criteria on DCA𝑥𝑦, DCA𝑧, the number of reconstructed space points in the TPC, and the quality of the track fit from their nominal values. Systematic uncertainties obtained from the contribution of each of these sources are considered as uncorrelated and the total systematic uncertainty on the observables is obtained by adding them in quadrature. Table 1shows the summary of the contributions to the total systematic uncertainty on standardized skewness, intensive skewness, and kurtosis in Pb–Pb, Xe–Xe, and pp collisions. The standardized skewness, intensive skewness, and kurtosis reported in this draft are considered robust and independent of detection efficiencies [17]. The detector inefficiencies cancel to leading order within these ratios, despite their potential presence in individual 𝑝T correlators. Consequently, efficiency correction is not carried out for these ratios. The robustness of the measured ratios is further affirmed by performing a Monte Carlo (MC) closure test. The MC closure test uses simulations based on different event generators for producing particles that corresponds to generated (true) results and GEANT3 [32] for the transport of particles through the geometry of ALICE detectors. For Pb–Pb and Xe–Xe collisions, the MC events are generated using HIJING (Heavy-Ion Jet Interaction Generator) [33]whereas for pp collisions, the events are produced by PYTHIA8 [34,35]event generator with Monash2013 tune. The experimental conditions prevailing during the data taking are accounted for in the generated events using GEANT3 by reproducing the actual configuration of different detectors during the runs. Results obtained from the generated events are compared with their corresponding reconstructed ones (without applying efficiency corrections), and they show a good agreement within uncertainties. The ratio of the reconstructed to the generated results for the observables is fitted with zeroth-order polynomial functions to quantify the amount of closure obtained. The percentages of agreement for standardized skewness, intensive skewness, and kurtosis are 97.3%, 99.6% and 99.1%, respectively, in pp collisions at √𝑠=5.02 TeV. Similar study is also performed for the other two collision systems, Pb–Pb collisions at √𝑠NN =5.02 and Xe–Xe collisions at √𝑠NN =5.44 TeV. The percentages of closure for the observables (standardized skewness, intensive skewness, kurtosis) in Pb–Pb (Xe–Xe) collisions are found to be 97.7% (97.4%), 96.4% (97.7%), and 99.7% (98.7%). These differences are added to the systematic uncertainties and reported in Table 1. To study how the skewness and kurtosis of ⟨𝑝T⟩distributions vary with the size of the collision system and to have a comparison between heavy-ion collisions (Pb–Pb and Xe–Xe) and small systems like pp collisions, expressing the centrality (multiplicity for pp collisions) classes as average charged particle multiplicity density (⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5) is more suitable. The conversion of centrality to ⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5is performed using the measured values of ⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5for certain centrality classes from Ref. [36]for Pb–Pb collisions. It should be noted that the centrality classes used in this analysis are narrower and different from that in the reference. A linear relation is observed between average number of charged particles (⟨𝑁acc⟩) in the kinematic range (|𝜂| <0.8 and 0.2 <𝑝 T<0.3GeV/c) of this analysis and ⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5over the full centrality range. The extracted fit parameters from the linear relation between ⟨𝑁acc⟩and ⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5were used to assign a value of ⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5for any value of ⟨𝑁acc⟩. Hence, by calculating the ⟨𝑁acc⟩for the centrality classes of this analysis, the corresponding value of ⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5is obtained. This method was previously used in the study of second-order fluctuations of ⟨𝑝T⟩in Ref. [8]. Similar approach is followed for Xe–Xe collisions using ⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5values of centrality classes from Ref. [37]. For pp collisions, the measured values of ⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5are taken from Ref. [30]since the multiplicity classes are the same. 5. Results 5.1. Standardized skewness The standardized skewness of charged particles produced in Pb–Pb collisions at √𝑠NN =5.02 TeV and Xe–Xe collisions at √𝑠NN =5.44 TeV as a function of system size denoted by the cubic root of average charged particle multiplicity density, ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5is shown in Fig. 3. The femtoscopic radii associated with the radius of the fireball at freeze-out, scale linearly with ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5and therefore it is used as a proxy for the system size [38]. The uncertainties on ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5are obtained by propagating the uncertainty on ⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5. It is observed that 𝛾⟨𝑝T⟩decreases with increasing ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5for both Pb–Pb and Xe–Xe collisions. This is generally expected from fluctuations in a dilution scenario caused by superposition of partially independent particle-emitting sources [7,8]. If the fluctuations in ⟨𝑝T⟩have been purely statistical and generated by finite multiplicity (𝑁), then 𝛾⟨𝑝T⟩would have a dependence of 1∕√𝑁. Aslight Physics Letters B 850 (2024) 138541 5 ALICE Collaboration Table 1 Contributions to the total systematic uncertainty on standardized skewness 𝛾⟨𝑝T⟩, intensive skewness Γ⟨𝑝T⟩, and kurtosis 𝜅⟨𝑝T⟩in Pb–Pb collisions at √𝑠NN =5.02 TeV, Xe–Xe collisions at √𝑠NN = 5.44 TeV, and pp collisions at √𝑠=5.02 TeV. Ranges are given where the uncertainties depend on centrality or multiplicity. Observables Sources of systematic uncertainty Pb–Pb Xe–Xe pp Vertex 𝑧-position 0.7–2.7% 1.5–9.1% 0.3–0.8% Pileup 0.5–6.4% – – 𝜒2 TPC∕𝑛TPCclusters 0.5–2.6% 3.4-6.5% 0.1–0.4% 𝜒2 ITS∕𝑛ITSclusters 0.5–2.2% 3.6–11.3% <0.1% 𝛾⟨𝑝T⟩𝑛TPCcrossedrows 0.6–3.8% 2.6–6.8% 0.9–1.9% DCA𝑥𝑦 0.6–4.5% 2.5–7.3% 0.5–1.1% DCA𝑧0.4–1.8% 2.8–5.9% <0.1% MC closure 2.3% 2.6% 2.7% Total 2.9–8.7% 9.1–17.2% 2.9–3.6% Vertex 𝑧-position 0.6–2.7% 1.7–9.1% 0.3–0.8% Pileup 0.5–6.3% – – 𝜒2 TPC∕𝑛TPCclusters 0.4–2.7% 3.6–6.9% 0.1–0.3% 𝜒2 ITS∕𝑛ITSclusters 0.4–2.2% 3.7–11.4% <0.1% Γ⟨𝑝T⟩𝑛TPCcrossedrows 0.8–3.8% 3.0–7.0% 0.3–1.7% DCA𝑥𝑦 0.7–4.7% 0.4–4.2% 0.3–1% DCA𝑧0.4–1.9% 2.7–6.0% <0.1% MC closure 3.6% 2.3% 0.4% Total 4–9% 9.5–17.6% 0.8–2.1% Vertex 𝑧-position 0.1–3.6% 0.9–4.5% 0.1–1.4% Pileup 0.2–4.5% – – 𝜒2 TPC∕𝑛TPCclusters 0.1–1.8% 0.9–4.6% 0.1–1.1% 𝜒2 ITS∕𝑛ITSclusters 0.1–0.9% 1.2–3.9% <0.1% 𝜅⟨𝑝T⟩𝑛TPCcrossedrows 0.1–1.3% 3.0–7.0% 0.3–1.7% DCA𝑥𝑦 0.1–2.3% 0.3–0.9% 0.5–3.3% DCA𝑧0.1–2.9% 0.7–1.9% <0.1% MC closure 0.3% 1.3% 0.9% Total 0.5–7.3% 2.7–11.7% % 1.8–4.7% Fig. 3. Standardized skewness 𝛾⟨𝑝T⟩shown as a function of ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5 in Pb–Pb collisions at √𝑠NN =5.02 TeV, Xe–Xe collisions at √𝑠NN =5.44 TeV, and pp collisions at √𝑠=5.02 TeV. The predictions from different event generators and hydrodynamic model calculations from V-USPHYDRO [17]and MC-Glauber+MUSIC [40,41]are denoted by colored bands. The statistical (systematic) uncertainties are represented by vertical bars (boxes). enhancement of 𝛾⟨𝑝T⟩in the most central (⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5>11) Pb– Pb collisions is also observed within uncertainty. One of the probable reasons could be the reduction of two-particle correlations ⟨Δ𝑝𝑖Δ𝑝𝑗⟩ towards the central collisions reported previously in Ref. [8]. Hydrodynamic model calculations from Ref. [17]that uses TRENTo initial conditions evolved by hydrodynamic code V-USPHYDRO referred as VUSPHYDRO model in the figures are compared with the data for both Pb– Pb and Xe–Xe collisions. A small specific shear viscosity, 𝜂∕𝑠=0.047 is used in the V-USPHYDRO model. As shown in Fig. 3, the V-USPHYDRO model calculations capture the general system size dependence of the measurement, but fail to describe it quantitatively. In the model, 𝛾⟨𝑝T⟩ is larger in Pb–Pb collisions compared to Xe–Xe collisions for almost all multiplicities, however data do not show such system dependence. Model calculations using HIJING [33], which incorporates various phenomena like multiple minijet production with initial and final state radiation or nuclear effects such as parton shadowing and jet quenching, are shown in the figure. HIJING considers nucleus–nucleus collisions as superposition of independent binary collisions of wounded nucleons. The latest version of the model, HIJING∕𝐵𝐵 v2.0 with shadowing and jet quenching effects [39]is used. The 𝛾⟨𝑝T⟩in HIJING model exhibits a strong dependence with ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5and reproduces the measurement in the semiperipheral to semicentral region. Calculations performed using the MC-Glauber+MUSIC model, which adopts the Glauber Monte Carlo [40] approach to generate initial conditions for the collisions, followed by the MUSIC hydrodynamic model [41]with 𝜂∕𝑠= 0.1 are presented. Notably, both the HIJING and MC-Glauber+MUSIC model utilize the Monte Carlo Glauber initial conditions, with the primary difference lying in the subsequent evolution. While HIJING lacks the implementation of collective phenomena, MUSIC incorporates such effects. Interestingly, the results obtained for the 𝛾⟨𝑝T⟩in Pb–Pb collisions remain same for both models from semiperipheral to semicentral collisions. This prompts the question of whether the 𝛾⟨𝑝T⟩is primarily sensitive to the details of the initial conditions. The results in Pb–Pb and Xe–Xe are also compared with the measurements carried out in pp collisions at √𝑠=5.02 TeV. The 𝛾⟨𝑝T⟩in pp collisions shows a steeper decrease with increasing ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5 compared to heavy-ion collisions. The observed 𝛾⟨𝑝T⟩is also found smaller in pp collisions than in Pb–Pb collisions in the overlapping Physics Letters B 850 (2024) 138541 6 ALICE Collaboration Fig. 4. Intensive skewness Γ⟨𝑝T⟩shown as a function of ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5in Pb– Pb collisions at √𝑠NN =5.02 TeV, Xe–Xe collisions at √𝑠NN =5.44 TeV, and pp collisions at √𝑠=5.02 TeV. The colored dashed lines represent the independent baseline for each system. The predictions from different event generators and hydrodynamic model calculations from V-USPHYDRO [17]and MC-Glauber+MUSIC [40,41]are denoted by colored bands. The statistical (systematic) uncertainties are represented by vertical bars (boxes). ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5region. The measurements in pp collisions are compared with PYTHIA8 [35]model, which is a pQCD based MC event generator. The PYTHIA8 model is found to describe numerous experimental results in pp collisions at LHC energies quite successfully [42–44]. This model can involve color reconnection (CR) mechanism [45]in pp collisions, which could contribute to explaining certain collective effects that have been observed in these collisions. Comparison of the measurements in pp collisions with PYTHIA8 model for both cases, with and without enabling CR mechanism (referred as PYTHIA8 CR ON and PYTHIA8 CR OFF), is presented. PYTHIA8 CR OFF model calculations utterly fail to describe the measurements both qualitatively and quantitatively. On the other hand, PYTHIA CR ON calculations reproduce the dependence of 𝛾⟨𝑝T⟩with ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5although with a larger slope. 5.2. Intensive skewness Fig. 4shows the intensive skewness studied with respect to ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5in Pb–Pb collisions at √𝑠NN =5.02 TeV and Xe–Xe collisions at √𝑠NN =5.44 TeV. The colored dashed lines indicate the independent baselines Γindependent evaluated for each system separately. The formula [17]used to calculate the independent baseline for intensive skewness is Γindependent =⟨(𝑝T−⟨𝑝T⟩)3⟩⟨𝑝T⟩ ⟨(𝑝T−⟨𝑝T⟩)2⟩2,(9) where ⟨𝑝T⟩is the average transverse momentum in given centrality class. The 𝑝Tspectra for different centrality classes are the input distributions for the calculation of Γindependent as a function of centrality. The second and third central moments of the 𝑝Tdistributions, ⟨(𝑝T−⟨𝑝T⟩)2⟩ and ⟨(𝑝T−⟨𝑝T⟩)3⟩are evaluated using the 𝑝Tspectrum for each centrality class and the ratio of these moments gives the Γindependent for that centrality class. For Pb–Pb and Xe–Xe collisions, the 𝑝Tspectra are taken from Refs. [46] and [47], respectively. Consecutively, the centralities are converted into ⟨d𝑁ch∕d𝜂⟩|𝜂|<0.5and ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5. A positive intensive skewness, larger than the independent baseline, is observed in both Pb–Pb and Xe–Xe collisions as a function of system size. This observation is in accordance with the predictions from the V-USPHYDRO hydrodynamic model calculations. The V-USPHYDRO model results describe the data qualitatively but do not show quantitative agreement. The HIJING model calculations though seem to capture the decreasing trend of the measurements in the semiperipheral to semicentral region (3.8 ≤⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5≤8.9) but are unable to explain the rise of Γ⟨𝑝T⟩in the measurements from the semicentral to central region (9.0 ≤⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5≤12.4). The observed increase in both 𝛾⟨𝑝T⟩ and Γ⟨𝑝T⟩for ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5>11 in Pb–Pb collsions might be indicative of the onset of thermalization [48,49]. This phenomenon is distinct from hydrodynamization, which could occur in smaller systems and at shorter timescales. As a direct consequence of thermalization, a significant increase of the skewness of ⟨𝑝T⟩fluctuations is predicted for most central collisions. If the system thermalizes, ⟨𝑝T⟩increases with density which is a characteristic phenomenon in relativistic fluid dynamics. At fixed multiplicity, larger density corresponds to smaller volume and correspondingly larger impact parameter which suggests that the rise of skewness results from a reduction in impact parameter fluctuations within the most central collisions. The MC-Glauber+MUSIC model calculations exhibit similar behavior observed in the measurements of Γ⟨𝑝T⟩in Pb–Pb collisions and provide better agreement quantitatively than the V-USPHYDRO model results. The increase of Γ⟨𝑝T⟩for central collisions in the measurements, seems to be depicted by both of the hydrodynamic model calculations with different initial conditions and therefore demands theoretical inputs to understand this rise. The measurement performed in pp collisions at √𝑠=5.02 TeV is also shown along with its independent baseline. For pp collisions, the baseline is calculated using the 𝑝Tspectra from Ref. [50]. The hydrodynamic prediction of positive intensive skewness above of its baseline is pronounced in the measurements in pp collisions as well. The distinct non-monotonic behavior observed in Pb–Pb collisions as a function of ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5is noticeably attenuated in pp collisions. The Γ⟨𝑝T⟩exhibits a subtle, monotonic decrease as ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5increases in pp collisions. In the PYTHIA8 CR OFF model, the results for Γ⟨𝑝T⟩exhibit an increase as ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5increases, which is in direct contrast to the observations from the PYTHIA8 CR ON model. While the PYTHIA8 CR ON model appears to be in closer agreement with the experimental measurements compared to CR OFF, neither model adequately captures the qualitative behavior observed in the measurements. 5.3. Kurtosis The kurtosis evaluated using Eq. (5)as a function of ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5 in Pb–Pb collisions at √𝑠NN =5.02 TeV and Xe–Xe collisions at √𝑠NN =5.44 TeV is shown in Fig. 5. The 𝜅⟨𝑝T⟩for a Gaussian distribution, which serves as a baseline for independent particle production, is shown with dotted line. The 𝜅⟨𝑝T⟩is found to decrease with increasing system size in Pb–Pb and Xe–Xe collisions and approaches the Gaussian baseline towards the most central Pb–Pb collisions. The measurements are compared with HIJING model calculations for both Pb–Pb and Xe–Xe collisions. While the HIJING model successfully captures the decreasing trend observed in the measurements with respect to ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5, it exhibits a pronounced and rapid decrease in 𝜅⟨𝑝T⟩compared to the data points. This is related to trivial system size dependence of 𝜅⟨𝑝T⟩in HIJING that should follow 1∕𝑁dependence. Throughout the entire range of ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5, the HIJING model overestimates the measurements. The MC-Glauber+MUSIC calculations reproduce the measurements in the central to midperipheral region (5.4 ≤⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5≤12.4). The measurements in pp collisions at √𝑠=5.02 TeV are also shown in the Fig. 5. In pp collisions, the 𝜅⟨𝑝T⟩decreases with ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5and remains consistently above the Gaussian baseline even for the highest value of ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5. Calculations of 𝜅⟨𝑝T⟩in PYTHIA8 show similar behavior with ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5as in data when CR mechanism is turned ON. PYTHIA8 CR OFF model completely fails to describe the measured 𝜅⟨𝑝T⟩. Physics Letters B 850 (2024) 138541 7 ALICE Collaboration Fig. 5. Kurtosis 𝜅⟨𝑝T⟩shown as a function of ⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5in Pb–Pb collisions at √𝑠NN =5.02 TeV, Xe–Xe collisions at √𝑠NN = 5.44 TeV and pp collisions at √𝑠=5.02 TeV. The predictions from different event generators and hydrodynamic calculations from MC-Glauber+MUSIC model [40,41]are represented by colored bands. The dashed line indicates the Gaussian baseline. The statistical (systematic) uncertainties are represented by vertical bars (boxes). 6. Summary In summary, first results on higher order fluctuations of mean transverse momentum (⟨𝑝T⟩) of charged particles, skewness, and kurtosis in Pb–Pb, Xe–Xe, and pp collisions at LHC energies are presented. Two measures of skewness, standardized skewness and intensive skewness, are studied using threeand two-particle 𝑝Tcorrelators as a function of system size (⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5). The standardized skewness is found to decrease with increasing system size in all three collision systems. Positive intensive skewness, larger than the baseline and consistent with the predictions from the hydrodynamics studies in Ref. [17], is observed in both Pb–Pb and Xe–Xe collisions. The anticipated positive intensive skewness for A–A collisions is likewise confirmed in pp collisions at similar collision energy. HIJING model calculations, which do not incorporate hydrodynamic evolution, are able to exhibit the qualitative behavior of intensive skewness in the region of 3.8 ≤⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5≤ 8.9. The striking rise of intensive skewness in the most central region is unexplained by HIJING. This effect is, however captured in the hydrodynamic calculations of both V-USPHYDRO and MC-Glauber+MUSIC models. The MC-Glauber+MUSIC model shows good agreement with the measurements of the standardized skewness and intensive skewness in the region of 4.8 ≤⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5≤12.4. In contrast, the V-USPHYDRO model, which utilizes TRENTo initial conditions, overestimates both measures of skewness. The discrepancy between the results obtained from the two hydrodynamic models suggests that these measurements can provide valuable insights into the initial stages of the collision. Also, the standardized skewness in HIJING and MCGlauber+MUSIC models exhibits similar behavior, suggesting that its sensitivity may primarily be attributed to the details of the initial conditions rather than the subsequent evolution. In pp collisions, the PYTHIA8 model with the color reconnection mechanism enabled is able to qualitatively reproduce the measurements. The kurtosis measured using the twoand four-particle 𝑝Tcorrelators shows a decrease with the system size. This decreasing behavior is reproduced by the HIJING model calculations in Pb–Pb and Xe–Xe collisions, which, however, do not show a quantitative agreement with the measurements. Hydrodynamic calculations from MC-Glauber+MUSIC model explain the measurement well in the region of 4.8 ≤⟨d𝑁ch∕d𝜂⟩1∕3 |𝜂|<0.5≤12.4. The value of kurtosis from the MC-Glauber+MUSIC model has a similar value to that from a Gaussian distribution. The measurements of kurtosis in Pb–Pb collisions also agreeing with the Gaussian baseline for highest values of the system size, hint about the production of a locally thermalized system in most central Pb–Pb collisions. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability This manuscript has associated data in a HEPData repository at: https://www .hepdata .net /record /ins2692420. Acknowledgements The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration gratefully acknowledges the resources and support provided by all Grid centres and the Worldwide LHC Computing Grid (WLCG) collaboration. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the ALICE detector: A. I. Alikhanyan National Science Laboratory (Yerevan Physics Institute) Foundation (ANSL), State Committee of Science and World Federation of Scientists (WFS), Armenia; Austrian Academy of Sciences, Austrian Science Fund (FWF): [M 2467-N36] and Nationalstiftung für Forschung, Technologie und Entwicklung, Austria; Ministry of Communications and High Technologies, National Nuclear Research Center, Azerbaijan; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Financiadora de Estudos e Projetos (Finep), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Universidade Federal do Rio Grande do Sul (UFRGS), Brazil; Bulgarian Ministry of Education and Science, within the National Roadmap for Research Infrastructures 20202027 (object CERN), Bulgaria; Ministry of Education of China (MOEC), Ministry of Science & Technology of China (MSTC) and National Natural Science Foundation of China (NSFC), China; Ministry of Science and Education and Croatian Science Foundation, Croatia; Centro de Aplicaciones Tecnológicas y Desarrollo Nuclear (CEADEN), Cubaenergía, Cuba; Ministry of Education, Youth and Sports of the Czech Republic, Czech Republic; The Danish Council for Independent Research | Natural Sciences, the Villum Fonden and Danish National Research Foundation (DNRF), Denmark; Helsinki Institute of Physics (HIP), Finland; Commissariat à l’Énergie Atomique (CEA) and Institut National de Physique Nucléaire et de Physique des Particules (IN2P3) and Centre National de la Recherche Scientifique (CNRS), France; Bundesministerium für Bildung und Forschung (BMBF) and GSI Helmholtzzentrum für Schwerionenforschung GmbH, Germany; General Secretariat for Research and Technology, Ministry of Education, Research and Religions, Greece; National Research, Development and Innovation Office, Hungary; Department of Atomic Energy, Government of India (DAE), Department of Science and Technology, Government of India (DST), University Grants Commission, Government of India (UGC) and Council of Scientific and Industrial Research (CSIR), India; National Research and Innovation Agency - BRIN, Indonesia; Istituto Nazionale di Fisica Nucleare (INFN), Italy; Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) and Japan Society for the Promotion of Science (JSPS) KAKENHI, Japan; Consejo Nacional de Ciencia (CONACYT) y Tecnología, through Fondo de Cooperación Internacional en Ciencia y Tecnología (FONCICYT) and Dirección General de Asuntos del Personal Académico (DGAPA), Mexico; Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO), Netherlands; The Research Council of Norway, Norway; Commission on Science and Technology for Sustainable Development in the South (COMSATS), Pakistan; Pontificia Universidad Católica del Perú, Peru; Ministry of Education and Physics Letters B 850 (2024) 138541 8 ALICE Collaboration Science, National Science Centre and WUT ID-UB, Poland; Korea Institute of Science and Technology Information and National Research Foundation of Korea (NRF), Republic of Korea; Ministry of Education and Scientific Research, Institute of Atomic Physics, Ministry of Research and Innovation and Institute of Atomic Physics and Universitatea Nationala de Stiinta si Tehnologie Politehnica Bucuresti, Romania; Ministry of Education, Science, Research and Sport of the Slovak Republic, Slovakia; National Research Foundation of South Africa, South Africa; Swedish Research Council (VR) and Knut & Alice Wallenberg Foundation (KAW), Sweden; European Organization for Nuclear Research, Switzerland; Suranaree University of Technology (SUT), National Science and Technology Development Agency (NSTDA) and National Science, Research and Innovation Fund (NSRF via PMU-B B05F650021), Thailand; Turkish Energy, Nuclear and Mineral Research Agency (TENMAK), Turkey; National Academy of Sciences of Ukraine, Ukraine; Science and Technology Facilities Council (STFC), United Kingdom; National Science Foundation of the United States of America (NSF) and United States Department of Energy, Office of Nuclear Physics (DOE NP), United States of America. In addition, individual groups or members have received support from: European Research Council, Strong 2020 - Horizon 2020 (grant nos. 950692, 824093), European Union; Academy of Finland (Center of Excellence in Quark Matter) (grant nos. 346327, 346328), Finland. References [1] E.V. Shuryak, Event per event analysis of heavy ion collisions and thermodynamical fluctuations, Phys. Lett. B 423 (1998) 9–14, arXiv :hep -ph /9704456. [2] M.A. Stephanov, K. Rajagopal, E.V. Shuryak, Event-by-event fluctuations in heavy ion collisions and the QCD critical point, Phys. Rev. D 60 (1999) 114028, arXiv : hep -ph /9903292. [3] H. Heiselberg, Event-by-event physics in relativistic heavy ion collisions, Phys. Rep. 351 (2001) 161–194, arXiv :nucl -th /0003046. [4] ALICE Collaboration, The ALICE experiment – a journey through QCD, CERN-EP2022-227, arXiv :2211 .04384 [nucl -ex], 11.2022. [5] E.V. Shuryak, Quark-gluon plasma and hadronic production of leptons, photons and psions, Phys. Lett. B 78 (1978) 150. [6] L. Stodolsky, Temperature fluctuations in multiparticle production, Phys. Rev. Lett. 75 (1995) 1044–1045. [7] STAR Collaboration, J. Adams, et al., Incident energy dependence of pt correlations at RHIC, Phys. Rev. C 72 (2005) 044902, arXiv :nucl -ex /0504031. [8] ALICE Collaboration, B. Abelev, et al., Event-by-event mean 𝑝𝐓fluctuations in pp and Pb-Pb collisions at the LHC, Eur. Phys. J. C 74 (2014) 3077, arXiv :1407 .5530 [nucl -ex]. [9] Z.-W. Lin, C.M. Ko, B.-A. Li, B. Zhang, S. Pal, A multi-phase transport model for relativistic heavy ion collisions, Phys. Rev. C 72 (2005) 064901, arXiv :nucl -th /0411110. [10] C. Zhang, L. Zheng, F. Liu, S. Shi, Z.-W. Lin, Update of a multiphase transport model with modern parton distribution functions and nuclear shadowing, Phys. Rev. C 99 (2019) 064906, arXiv :1903 .03292 [nucl -th]. [11] F. Gelis, E. Iancu, J. Jalilian-Marian, R. Venugopalan, The color glass condensate, Annu. Rev. Nucl. Part. Sci. 60 (2010) 463–489, arXiv :1002 .0333 [hep -ph]. [12] H. Weigert, Evolution at small x(bj): the color glass condensate, Prog. Part. Nucl. Phys. 55 (2005) 461–565, arXiv :hep -ph /0501087. [13] S. Gavin, G. Moschelli, Fluctuation probes of early-time correlations in nuclear collisions, Phys. Rev. C 85 (2012) 014905, arXiv :1107 .3317 [nucl -th]. [14] S. Bhadury, M. Kurian, V. Chandra, A. Jaiswal, Second order relativistic viscous hydrodynamics within an effective description of hot QCD medium, J. Phys. G 48 (2021) 105104, arXiv :2010 .01537 [hep -ph]. [15] U. Heinz, C. Shen, H. Song, The viscosity of quark-gluon plasma at RHIC and the LHC, AIP Conf. Proc. 1441 (2012) 766–770, arXiv :1108 .5323 [nucl -th]. [16] P. Romatschke, New developments in relativistic viscous hydrodynamics, Int. J. Mod. Phys. E 19 (2010) 1–53, arXiv :0902 .3663 [hep -ph]. [17] G. Giacalone, F.G. Gardim, J. Noronha-Hostler, J.-Y. Ollitrault, Skewness of mean transverse momentum fluctuations in heavy-ion collisions, Phys. Rev. C 103 (2021) 024910, arXiv :2004 .09799 [nucl -th]. [18] F.G. Gardim, G. Giacalone, M. Luzum, J.-Y. Ollitrault, Effects of initial state fluctuations on the mean transverse momentum, Nucl. Phys. A 1005 (2021) 121999, arXiv :2002 .07008 [nucl -th]. [19] F.G. Gardim, G. Giacalone, M. Luzum, J.-Y. Ollitrault, Thermodynamics of hot strong-interaction matter from ultrarelativistic nuclear collisions, Nat. Phys. 16 (2020) 615–619, arXiv :1908 .09728 [nucl -th]. [20] J.S. Moreland, J.E. Bernhard, S.A. Bass, Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions, Phys. Rev. C 92 (2015) 011901, arXiv :1412 .4708 [nucl -th]. [21] J. Noronha-Hostler, G.S. Denicol, J. Noronha, R.P.G. Andrade, F. Grassi, Bulk viscosity effects in event-by-event relativistic hydrodynamics, Phys. Rev. C 88 (2013) 044916, arXiv :1305 .1981 [nucl -th]. [22] STAR Collaboration, S.A. Voloshin, Multiplicity and mean transverse momentum fluctuations in Au+Au collisions at RHIC, AIP Conf. Proc. 610 (2002) 591–596, arXiv :nucl -ex /0109006. [23] S.A. Voloshin, Mean p(t) fluctuations from two particle and four particle correlations, arXiv :nucl -th /0206052. [24] ALICE Collaboration, K. Aamodt, et al., The ALICE experiment at the CERN LHC, J. Instrum. 3 (2008) S08002. [25] ALICE Collaboration, B. Abelev, et al., Performance of the ALICE experiment at the CERN LHC, Int. J. Mod. Phys. A 29 (2014) 1430044, arXiv :1402 .4476 [nucl -ex]. [26] J. Alme, et al., The ALICE TPC, a large 3-dimensional tracking device with fast readout for ultra-high multiplicity events, Nucl. Instrum. Methods A 622 (2010) 316–367, arXiv :1001 .1950 [physics .ins -det]. [27] ALICE Collaboration, E. Abbas, et al., Performance of the ALICE VZERO system, J. Instrum. 8 (2013) P10016, arXiv :1306 .3130 [nucl -ex]. [28] ALICE Collaboration ALICE, Centrality determination in heavy ion collisions, ALICEPUBLIC-2018-011, http://cds .cern .ch /record /2636623, 2018. [29] ALICE Collaboration, Centrality determination using the Glauber model in XeXe collisions at √𝑠NN =5.44 TeV, ALICE-PUBLIC-2018-003, https://cds .cern .ch / record /2315401, 2018. [30] ALICE Collaboration, S. Acharya, et al., Pseudorapidity distributions of charged particles as a function of midand forward rapidity multiplicities in pp collisions at √𝑠=5.02, 7and 13 TeV, Eur. Phys. J. C 81 (2021) 630, arXiv :2009 .09434 [nucl -ex]. [31] X. Luo, J. Xu, B. Mohanty, N. Xu, Volume fluctuation and auto-correlation effects in the moment analysis of net-proton multiplicity distributions in heavy-ion collisions, J. Phys. G 40 (2013) 105104, arXiv :1302 .2332 [nucl -ex]. [32] R. Brun, F. Bruyant, F. Carminati, S. Giani, M. Maire, A. McPherson, G. Patrick, L. Urban, GEANT detector description and simulation tool, CERN-W5013, http:// cds .cern .ch /record /1082634, 10.1994. [33] X.-N. Wang, M. Gyulassy, HIJING: a Monte Carlo model for multiple jet production in p p, p A and A A collisions, Phys. Rev. D 44 (1991) 3501–3516. [34] T. Sjöstrand, S. Mrenna, P.Z. Skands, PYTHIA 6.4 physics and manual, J. High Energy Phys. 05 (2006) 026, arXiv :hep -ph /0603175. [35] T. Sjöstrand, S. Ask, J.R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C.O. Rasmussen, P.Z. Skands, An introduction to PYTHIA 8.2, Comput. Phys. Commun. 191 (2015) 159–177, arXiv :1410 .3012 [hep -ph]. [36] ALICE Collaboration, J. Adam, et al., Centrality dependence of the charged-particle multiplicity density at midrapidity in Pb-Pb collisions at √𝑠NN =5.02 TeV, Phys. Rev. Lett. 116 (2016) 222302, arXiv :1512 .06104 [nucl -ex]. [37] ALICE Collaboration, S. Acharya, et al., Centrality and pseudorapidity dependence of the charged-particle multiplicity density in Xe–Xe collisions at √𝑠NN =5.44 TeV, Phys. Lett. B 790 (2019) 35–48, arXiv :1805 .04432 [nucl -ex]. [38] ALICE Collaboration, K. Aamodt, et al., Two-pion Bose-Einstein correlations in central Pb-Pb collisions at √𝑠𝑁𝑁 =2.76 TeV, Phys. Lett. B 696 (2011) 328–337, arXiv :1012 .4035 [nucl -ex]. [39] V. Topor Pop, M. Gyulassy, J. Barrette, C. Gale, X.N. Wang, N. Xu, Baryon junction loops and the baryon-meson anomaly at high energies, Phys. Rev. C 70 (2004) 064906, arXiv :nucl -th /0407095. [40] B. Alver, M. Baker, C. Loizides, P. Steinberg, The PHOBOS Glauber Monte Carlo, arXiv :0805 .4411 [nucl -ex]. [41] B. Schenke, S. Jeon, C. Gale, (3+1) D hydrodynamic simulation of relativistic heavyion collisions, Phys. Rev. C 82 (2010) 014903, arXiv :1004 .1408 [hep -ph]. [42] ALICE Collaboration, J. Adam, et al., Charged-particle multiplicities in proton– proton collisions at √𝑠=0.9to 8 TeV, Eur. Phys. J. C 77 (2017) 33, arXiv : 1509 .07541 [nucl -ex]. [43] ALICE Collaboration, S. Acharya, et al., Charged-particle multiplicity distributions over a wide pseudorapidity range in proton-proton collisions at √𝑠=0.9, 7, and 8 TeV, Eur. Phys. J. C 77 (2017) 852, arXiv :1708 .01435 [hep -ex]. [44] ALICE Collaboration, B. Abelev, et al., Multiplicity dependence of the average transverse momentum in pp, p-Pb, and Pb-Pb collisions at the LHC, Phys. Lett. B 727 (2013) 371–380, arXiv :1307 .1094 [nucl -ex]. [45] T. Sjostrand, M. van Zijl, A multiple interaction model for the event structure in hadron collisions, Phys. Rev. D 36 (1987) 2019. [46] ALICE Collaboration, S. Acharya, et al., Transverse momentum spectra and nuclear modification factors of charged particles in pp, p-Pb and Pb-Pb collisions at the LHC, J. High Energy Phys. 11 (2018) 013, arXiv :1802 .09145 [nucl -ex]. [47] ALICE Collaboration, S. Acharya, et al., Transverse momentum spectra and nuclear modification factors of charged particles in Xe-Xe collisions at √𝑠NN =5.44 TeV, Phys. Lett. B 788 (2019) 166–179, arXiv :1805 .04399 [nucl -ex]. [48] R. Samanta, S. Bhatta, J. Jia, M. Luzum, J.-Y. Ollitrault, Thermalization at the femtoscale seen in high-energy Pb+Pb collisions, arXiv :2303 .15323 [nucl -th]. [49] R. Samanta, J.a.P. Picchetti, M. Luzum, J.-Y. Ollitrault, Non-Gaussian transverse momentum fluctuations from impact parameter fluctuations, arXiv :2306 .09294 [nucl -th]. [50] ALICE Collaboration, S. Acharya, et al., Charged-particle production as a function of multiplicity and transverse spherocity in pp collisions at √𝑠=5.02 and 13 TeV, Eur. Phys. J. C 79 (2019) 857, arXiv :1905 .07208 [nucl -ex].