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

The first measurements of skewness and kurtosis of mean transverse momentum ($\langle p_\mathrm{T}\rangle$) fluctuations are reported in Pb$-$Pb collisions at $\sqrt{s_\mathrm{NN}}$ = 5.02 TeV, Xe$-$Xe collisions at $\sqrt{s_\mathrm{NN}}$ $=$ 5.44 TeV and pp collisions at $\sqrt{s} = 5.02$ TeV using the ALICE detector. The measurements are carried out as a function of system size $\langle \mathrm{d}N_\mathrm{ch}/\mathrm{d}\eta\rangle_{|\eta|<0.5}^{1/3}$, using charged particles with transverse momentum ($p_\mathrm{T}$) and pseudorapidity ($\eta$), in the range $0.2<p_\mathrm{T}<3.0$ GeV/$c$ and $|\eta|<0.8$, respectively. In Pb$-$Pb and Xe$-$Xe collisions, positive skewness is observed in the fluctuations of $\langle p_\mathrm{T}\rangle$ 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 $\langle p_\mathrm{T}\rangle$ 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.


Introduction
The properties of the hot and dense nuclear matter created in heavy-ion collisions at relativistic energies can be studied using event-by-event fluctuations of different quantities like multiplicity, net-charge, mean transverse momentum (⟨p T ⟩), etc. [1][2][3][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 quarkgluon 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 ⟨p T ⟩ of the final-state particles [6].A non-monotonic behavior of ⟨p T ⟩ fluctuations as a function of centrality or incident energy was suggested as one of the possible signals of the QGP [3].However, ⟨p T ⟩ fluctuations are also affected by nonthermodynamic 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 ⟨p 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 ⟨p T ⟩ fluctuations in Pb-Pb collisions at centre-of-mass energy per nucleon pair, √ s NN = 2.76 TeV in ALICE [8] suggested a connection of the observed fluctuations of the ⟨p 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 ⟨p 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][15][16].In Ref. [17], it has been proposed that the skewness of ⟨p T ⟩ fluctuations can serve as an essential probe of the hydrodynamic behavior of the system created in heavy-ion collisions.The ⟨p 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 ⟨p 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 ⟨p T ⟩ fluctuations are driven by the skewness of the initial energy density fluctuations, which implies that the ⟨p T ⟩ fluctuations arise from the same collective dynamics in the QGP phase that give rise to anisotropic flow.Furthermore, employing initial conditions from the T R ENTo model [20] and evolving them with the V-USPHYDRO viscous hydrodynamic [21] simulations predicted positive skewness of ⟨p T ⟩ fluctuations surpassing expectations in independent particle emission scenarios.Additional information regarding these simulations can be found in Ref. [17].
The fluctuations in ⟨p 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 ⟨p 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 ⟨p T ⟩ fluctuations that represent the third-and fourth-order fluctuations of ⟨p T ⟩ are reported at the Large Hadron Collider (LHC) energies.The observables used in this analysis are introduced in Sec. 2. Here, the multiparticle p T correlators used to study the ⟨p 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 ⟨p 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.

Observables
In this analysis, the fluctuations in event-by-event ⟨p T ⟩ of charged particles are investigated using multiparticle p T correlators.The event-by-event ⟨p T ⟩ is defined as where p T,i is the transverse momentum of the ith particle and N 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 ⟨p T ⟩ fluctuations.The moment method calculates various order moments of the ⟨p 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 p T correlators [22,23] is that they yield zero values for events with randomly sampled particles, thereby effectively isolating the nonstatistical fluctuations of interest.The expressions for the two-particle p T correlator (⟨∆p T,i ∆p T, j ⟩) and three-particle p T correlator (⟨∆p T,i ∆p T, j ∆p T,k ⟩) [17] are denoted by Eqs. 3 and 4, respectively, where Q n is defined as 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 Q n 's where n = 1, 2, ... were derived to account for the practical limitations of calculating multiparticle p T correlators with nested loops, especially when dealing with events characterized by large multiplicities.Since, the values of Q n can be obtained with Eq. 2 employing a single loop, it becomes feasible to directly utilize the later expressions with Q n s in Eqs. 3 and 4 for the calculation of two-and three-particle p T correlators, respectively.On similar lines, the four-particle p T correlator (⟨∆p T,i ∆p T, j ∆p T,k ∆p T,l ⟩) is derived as shown in Eq. 5.These two-, three-, and four-particle correlators are related to the variance, skewness, and kurtosis of the ⟨p T ⟩ distribution, respectively.
In the above equations, angular brackets ⟨. ..⟩ ev denote an average over all events and ⟨⟨p T ⟩⟩ = ⟨Q 1 /N ch ⟩ ev is the event-by-event ⟨p T ⟩ averaged over all events.The standardized skewness (γ ⟨p T ⟩ ) and the intensive skewness (Γ ⟨p T ⟩ ) are two measures of the skewness of the ⟨p T ⟩ distribution, constructed using the twoparticle and three-particle p T correlators [17], given by and The kurtosis κ ⟨p T ⟩ of ⟨p T ⟩ distribution is defined as The article focuses on analyzing higher order fluctuations of ⟨p T ⟩, particularly the third and fourth order.
Three key observables, namely the standardized skewness, intensive skewness, and kurtosis of the ⟨p T ⟩ distribution (defined in Eqs.6, 7, and 8, respectively), are investigated as a function of the system size in Pb-Pb collisions at √ s NN = 5.02 TeV, Xe-Xe collisions at √ s NN = 5.44 TeV, and pp collisions at √ s = 5.02 TeV.

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.1 and −3.7 < η < −1.7, respectively [27].The data analyzed here are obtained from Pb-Pb collisions at √ s NN = 5.02 TeV, Xe-Xe collisions at √ s NN = 5.44 TeV, and pp collisions at √ s = 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 0.5 0.6 0.7 0.8 0.9 1) c (GeV/ 〉 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 < p T < 3.0 GeV/c and |η| < 0.8, where p T is 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-ofclosest-approach (DCA) to the primary vertex in the plane perpendicular (DCA xy ) and parallel (DCA z ) 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 χ 2 per 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.

Data analysis
The mean transverse momentum ⟨p T ⟩ distributions obtained for three different centrality classes in Pb-Pb collisions at √ s 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 ⟨p T ⟩ distributions are 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 √ s NN = 5.44 TeV and pp collisions at √ s = 5.02 TeV.
Before starting the analysis with the ⟨p T ⟩ correlators presented in Sec. 2, it is important to demonstrate that these fluctuations are not a trivial consequence of fluctuating N ch in a given centrality class.The definition of ⟨p T ⟩ in Eq. 1 clearly shows that the varying multiplicity in different events of the same centrality class can affect the ⟨p T ⟩ distributions and consequently its fluctuations.To examine whether  Comparison of event-by-event mean transverse momentum ⟨p T ⟩ distribution from actual analysis and modified analysis in which N ch in each event is fixed to N min ch (N min ch is the minimum number of charged particles in a given centrality class) and N 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.
these fluctuations persist after removing the stochastic effects of the multiplicity, a check is performed.For each centrality class, the minimum number of tracks per event in the kinematic acceptance of the analysis, N min ch is determined.In each event of the centrality class, the N min ch number of tracks are then selected randomly to calculate the genuine ⟨p T ⟩ distribution free from biases of multiplicity fluctuations.The left panel of Fig. 2 shows the ⟨p T ⟩ distribution for the 0-5% central events in Pb-Pb collisions at √ s NN = 5.02 TeV.The distribution in red markers shows the original event-by-event ⟨p T ⟩ distribution in which N ch fluctuates whereas the distribution shown with black markers is the ⟨p T ⟩ distribution in which N ch is fixed to N 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 ⟨p T ⟩ distributions are positively skewed.The ⟨p T ⟩ distributions exhibit substantial deviations from a Gaussian distribution even under the condition of a fixed N ch , thereby reflecting the analogous traits observed in the original distribution with fluctuating N ch .The middle and right panels of Fig. 2 visually exemplify the consistent manifestation of these characteristics within the distributions of ⟨p T ⟩ for both semicentral and peripheral collisions.It is concluded that the positive skewness in the event-by-event ⟨p T ⟩ distributions is not a trivial consequence of event-by-event statistical fluctuations of N ch , and therefore new and non-trivial information can be extracted from higher-order moments of ⟨p T ⟩ fluctuations.
The ⟨p 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 variability 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 xy , DCA z , 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 1 shows 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 p 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 are 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 √ s = 5.02 TeV.Similar study is also performed for the other two collision systems, Pb-Pb collisions at √ s NN = 5.02 and Xe-Xe collisions at √ s 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 ⟨p 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 (⟨dN ch /dη⟩ |η|<0.5 ) is more suitable.The conversion of centrality to ⟨dN ch /dη⟩ |η|<0.5 is performed using the measured values of ⟨dN ch /dη⟩ |η|<0.5 for 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 (⟨N acc ⟩) in the kinematic range (|η| < 0.8 and 0.2 < p T < 0.3 GeV/c) of this analysis and ⟨dN ch /dη⟩ |η|<0.5 over the full centrality range.The extracted fit parameters from the linear relation between ⟨N acc ⟩ and ⟨dN ch /dη⟩ |η|<0.5 were used to assign a value of ⟨dN ch /dη⟩ |η|<0.5 for any value of ⟨N acc ⟩.Hence, by calculating the ⟨N acc ⟩ for the centrality classes of this analysis, the corresponding value of ⟨dN ch /dη⟩ |η|<0.5 is obtained.This method was previously used in the study of second-order fluctuations of ⟨p T ⟩ in Ref. [8].Similar approach is followed for Xe-Xe collisions using ⟨dN ch /dη⟩ |η|<0.5 values of centrality classes from Ref. [37].For pp collisions, the measured values of ⟨dN ch /dη⟩ |η|<0.5 are taken from Ref. [30] since the multiplicity classes are the same.

Standardized skewness
The standardized skewness of charged particles produced in Pb-Pb collisions at √ s NN = 5.02 TeV and Xe-Xe collisions at √ s NN = 5.44 TeV as a function of system size denoted by the cubic root of average charged particle multiplicity density, ⟨dN ch /dη⟩ 1/3 |η|<0.5 is shown in Fig. 3.The femtoscopic radii associated with the radius of the fireball at freeze-out, scale linearly with ⟨dN ch /dη⟩ 1/3 |η|<0.5 and therefore it is used as a proxy for the system size [38].The uncertainties on ⟨dN ch /dη⟩ |η|<0.5 > 11) Pb-Pb collisions is also observed within uncertainty.One of the probable reasons could be the reduction of two-particle correlations ⟨∆p i ∆p j ⟩ towards the central collisions reported previously in Ref. [8].Hydrodynamic model calculations from Ref. [17] that uses T R ENTo initial conditions evolved by hydrodynamic code V-USPHYDRO referred as V-USPHYDRO model in the figures are compared with the data for both Pb-Pb and Xe-Xe collisions.A small specific shear viscosity, η/s = 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, γ ⟨p 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/BB v2.0 with shadowing and jet quenching effects [39] is used.The γ ⟨p T ⟩ in HIJING model exhibits a strong dependence with ⟨dN ch /dη⟩ 1/3 |η|<0.5 and 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 η/s = 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 γ ⟨p T ⟩ in Pb-Pb collisions remain same for both models from semiperipheral to semicentral collisions.This prompts the question of whether the γ ⟨p 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 √ s = 5.02 TeV.The γ ⟨p T ⟩ in pp collisions shows a steeper decrease with increasing ⟨dN ch /dη⟩ 1/3 |η|<0.5 compared to heavy-ion collisions.The observed γ ⟨p T ⟩ is also found smaller in pp collisions than in Pb-Pb collisions in the overlapping ⟨dN ch /dη⟩ 1/3 |η|<0.5 region.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][43][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), are presented.PYTHIA8 CR OFF   although with a larger slope.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

Intensive skewness
where ⟨p T ⟩ is the average transverse momentum in given centrality class.The p T spectra 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 p T distributions, ⟨(p T − ⟨p T ⟩) 2 ⟩ and ⟨(p T − ⟨p T ⟩) 3 ⟩ are evaluated using the p T spectrum 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 p T spectra are taken from Refs.[46] and [47], respectively.Consecutively, the centralities are converted into ⟨dN ch /dη⟩ |η|<0.5 and ⟨dN 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 ≤ ⟨dN ch /dη⟩ 1/3 |η|<0.5 ≤ 8.9) but are unable to explain the rise of Γ ⟨p T ⟩ in the measurements  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).
from the semicentral to central region (9.0 ≤ ⟨dN ch /dη⟩ 1/3 |η|<0.5 ≤ 12.4).The observed increase in both γ ⟨p T ⟩ and Γ ⟨p T ⟩ for ⟨dN 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 ⟨p T ⟩ fluctuations is predicted for most central collisions.If the system thermalizes, ⟨p 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 Γ ⟨p T ⟩ in Pb-Pb collisions and provide better agreement quantitatively than the V-USPHYDRO model results.The increase of Γ ⟨p T ⟩ for central collisions in the measurements, seems to be depicted by both of the hydrodynamic model calculations with different initial conditions and therfore demands theoretical inputs to understand this rise.
The measurement performed in pp collisions at √ s = 5.02 TeV is also shown along with its independent baseline.For pp collisions, the baseline is calculated using the p T spectra 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 ⟨dN ch /dη⟩

Kurtosis
The kurtosis evaluated using Eq. 5 as a function of ⟨dN ch /dη⟩  5.The κ ⟨p T ⟩ for a Gaussian distribution, which serves as a baseline for independent particle production, is shown with dotted line.The κ ⟨p 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 capture the decreasing trend observed in the measurements with respect to ⟨dN ch /dη⟩ 1/3 |η|<0.5 , it exhibits a pronounced and rapid decrease in κ ⟨p T ⟩ compared to the data points.This is related to trivial system size dependence of κ ⟨p T ⟩ in HIJING that should follow 1/N dependence.Throughout the entire range of ⟨dN ch /dη⟩

Summary
In summary, first results on higher order fluctuations of mean transverse momentum (⟨p 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 three-and two-particle p T correlators as a function of system size (⟨dN 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 ≤ ⟨dN 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 ≤ ⟨dN ch /dη⟩ 1/3 |η|<0.5 ≤ 12.4.In contrast, the V-USPHYDRO model, which utilizes T R ENTo 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 MC-Glauber+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 two-and four-particle p T correlators 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 ≤ ⟨dN 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.

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 In

Figure 1 :
Figure 1: Efficiency uncorrected mean transverse momentum distributions for different centrality classes in Pb-Pb collisions at √ s NN = 5.02 TeV.

Figure 2 :
Figure2: Comparison of event-by-event mean transverse momentum ⟨p T ⟩ distribution from actual analysis and modified analysis in which N ch in each event is fixed to N min ch (N min ch is the minimum number of charged particles in a given centrality class) and N 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.
model calculations utterly fail to describe the measurements both qualitatively and quantitatively.On the other hand, PYTHIA CR ON calculations reproduce the dependence of γ ⟨p T ⟩ with ⟨dN ch /dη⟩ 1/3 |η|<0.5

Figure 4
Figure 4 shows the intensive skewness studied with respect to ⟨dN ch /dη⟩ 1/3 |η|<0.5 in Pb-Pb collisions at √ s NN = 5.02 TeV and Xe-Xe collisions at √ s 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

1 / 3 |η|<0. 5 3 |η|<0. 5
is noticeably attenuated in pp collisions.The Γ ⟨p T ⟩ exhibits a subtle, monotonic decrease as ⟨dN ch /dη⟩ 1/increases in pp collisions.In the PYTHIA8 CR OFF model, the results for Γ ⟨p T ⟩ exhibit an increase as ⟨dN ch /dη⟩ 1/3 |η|<0.5 increases, 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.
stitute) 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 2020-2027 (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'Energie Atomique (CEA) and

Table 1 :
Contributions to the total systematic uncertainty on standardized skewness γ ⟨p T ⟩ , intensive skewness Γ ⟨p T ⟩ , and kurtosis κ ⟨p T ⟩ in Pb-Pb collisions at √ s NN = 5.02 TeV, Xe-Xe collisions at √ s NN = 5.44 TeV, and pp collisions at √ s = 5.02 TeV.Ranges are given where the uncertainties depend on centrality or multiplicity.