SStudy of the η to π Ratio in Heavy-Ion Collisions
Yuanjie Ren ∗ and Axel Drees † Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA Department of Physics and Astronomy, Stony Brook University, Stony Brook, New York 11790, USA (Dated: February 11, 2021)We demonstrate that the p T dependence of the η/π ratio is universal within a few percent forhigh energy p + p , p +A and d +A collisions, over a broad range of collision energies. The η/π ratioincreases with p T up to 4 to 5 GeV/ c where it saturates at a nearly constant value of 0.487 ± p T = 5 GeV/ c the same constant value is also observed in A+A collisions independent ofcollision system, energy, and centrality. At lower p T , where accurate η/π data is absent for A+Acollisions, we estimate possible deviations from the universal behavior, which could arise due to therapid radial hydrodynamic expansion of the A+A collision system. For A+A collisions at RHIC wefind that possible deviations are limited to the p T range from 0.4 to 3 GeV/ c , and remain less than20% for the most central collisions. I. INTRODUCTION
Photons are generally considered ideal probes to studythe quark gluon plasma (QGP) created in heavy ion colli-sions [1], since they have a long mean free path and leavethe collision volume without final state interactions. Ofparticular interest are low momentum or thermal photonswith energies of up to several times the temperature ofthe QGP. The measurement of thermal photons has onlyrecently been possible with the advance of the heavy ionprograms at RHIC [2–4] and LHC [5].One of the experimental key challenges for these mea-surements is to estimate and subtract photons fromhadron decays that constitute the bulk of photons mea-sured in experiments. The two major contributions ofphotons result from π → γ + γ and η → γ + γ decays.Precise knowledge of the parent π and η p T spectrais necessary to estimate the decay photon background.While spectra of pions from heavy ion collisions are wellmeasured at RHIC and LHC, less data exists for η spec-tra, in particular below p T of 2 GeV/ c . Therefor exper-iments need to make assumptions how to model the η spectra below 2 GeV/ c , which leads to sizable system-atic uncertainties. Frequently, experiments have basedthis extrapolation on the hypothesis of transverse mass m T scaling of meson spectra [3–5]. However, it is knownsince the late 1990’s [6] and was recently pointed outagain [7] that m T scaling does not hold below 3 GeV forthe η meson.In this paper we propose a new empirical approach tomodel the η spectrum that is based on the universalityof the η/π ratio across collision systems, beam energies,and centrality selections in heavy ion collisions. With agood understanding of the η/π ratio as function of trans-verse momentum p T and measured π spectra, which arereadily available for many collision systems, one can con-struct a more accurate p T distribution for η mesons.The paper is organized as follows. In the next sec-tion we elaborate more on the failure of m T scaling. In ∗ [email protected] † [email protected] section III we will discuss two empirical fits and a Gaus-sian Process Regression (GPR) to describe the η/π ratiofor p + p and p +A collisions, and document in section IVthe universality of η/π across different collision systems( p + p , p +A, A+A), energies, and collision centrality. InSection V, we estimate possible deviation from the uni-versal trend at low p T due to radial flow in heavy ioncollisions. We provide our result for η/π for RHIC andLHC energies with systematic uncertainties in the finalpart. II. THE FAILURE OF TRANSVERSE MASSSCALING
For measurements of direct photons from heavy ioncollisions, the photons from η and heavier meson decaysare frequently estimated using measured π spectra inconjunction with the m T scaling hypothesis. A typicalimplementation of this method [8] starts with a fit tothe π spectra with a functional form like a modifiedHagedorn function [9]:12 πp T d N d y d p T = A ( M X ) (cid:18) e − ag ( p T ) − bg ( p T ) + g ( p T , M X ) p (cid:19) − n (1)with M X being the meson mass and g ( p T , M X ) = (cid:112) p T + m X − m π . In this implementation the spectraof the η and heavier mass mesons follow the same distri-bution with respect to transverse mass m T ≡ (cid:112) m + p T as the π . The normalisation constant A ( M X ) is the onlyfree parameter, all other parameters are fixed by the fitto the π data. A ( M X ) is fitted to experimental datawhenever such data exists.Fig. 1 compiles available data of η/π for p + p [10–14]and p +A [6, 15] collisions. Also shown on the figure isthe result of m T scaling for two different normalisationconstants A ( M η ) [11, 12, 16] and the expectation from a pythia -6 calculation from [12, 17]. While pythia andthe m T scaling hypothesis agree well, a significant devi-ation from the data is seen at low p T . This was origi-nally discovered at the CERN SPS by CERES/TAPS [6] a r X i v : . [ nu c l - e x ] F e b [GeV/c] T p -
10 1 10 p / h - ALICE p+p 8 TeV (2018)ALICE p+p 2.76 TeV (2017)ALICE p+p 7 TeV (2012)PHENIX p+p 200 GeV (2010)PHENIX p+p 200 GeV (2006)Systematic uncertainties
ALICE p+Pb 5.02 TeV (2018)CERES p+Au 29.1 GeV (1998)CERES p+Be 29.1 GeV (1998))=0.476 ¥ ( p / h scaling, p+p 7 TeV, T m )=0.5 ¥ ( p / h scaling, p+p 200 GeV, T m FIG. 1. The η/π ratio of p+p and p+A collisions. Alsoplotted are the η/π determined by m T -scaling and from a pythia calculation. more than 20 years ago and recently confirmed by AL-ICE at the LHC [10]. Clearly the m T scaling hypothesisis not correct and should not be used to extrapolate me-son spectra to low p T for systems where no data exists. III. DESCRIPTION OF THE η/π RATIO FORP+P AND P+A COLLISIONS
The quantitative agreement of the η/π data shownin Fig. 1 is striking, consider the data covers more than2 orders of magnitude in collision energy. In this sec-tion we will test different methods to obtain an empiri-cal description of η/π . The first two methods (A,B) fita functional shape of the ratio, while the third method(GPR) is a Gaussian Process Regression that does notassume a specific functional shape. All methods yieldsimilar results below 10 GeV/ c , at larger p T the devia-tions are sizable and we will include these deviations inour evaluation of systematic uncertainties. A. Empirical fit A
Method A starts with a ratio of two functions of theform given in Equation 1. The m T -scaling hypothesis isused to reduce the number of parameters: R η/π ( p T ) = R ∞ (cid:16) e − a · g ( p T ) − b · g ( p T ) + g ( p T ) p (cid:17) − n (cid:16) e − ap T − bp T + p T p (cid:17) − n . (2)The advantage of this method is that it preserves arealistic functional form for the p T spectra with an expo-nential decrease at low p T and power law shape at high p T . In principle, this ensures that at high p T the η/π ratio approaches a constant value R ∞ . However, unlikestarting from the π spectrum, the parameters are fittedto the η/π ratio from p+p and p+A collisions shown inFig. 1. We achieve a good fit, though the values of thefit parameters are nonphysical and do not describe theindividual p T spectra. The result is depicted in Fig. 2.The band represents the total uncertainty of the fitfunction from two sources, the uncertainty of fit param-eters, and the systematic uncertainties from data points.The former can be calculated analytically thanks to theexplicit fit function while the latter can be obtained via a“data shuffling approach” which uses a Monte Carlo tech-nique to vary individual data sets within their systematicuncertainties. This approach is discussed in Appendix B.The total uncertainty shown on the figure represents thequadratic sum of statistical and systematic uncertainties. B. Empirical fit B
The second empirical fit function has a very similarform, except that normalization of the exponential andpower law component in the numerator are decoupledby introducing an additional parameter. This is imple-mented such that R ∞ remains the asymptotic value athigh p T . R η/π ( p T ) = A (cid:18) e − a · g ( p T ) − b · g ( p T ) + (cid:0) R ∞ A (cid:1) − n g ( p T ) p (cid:19) − n (cid:16) e − ap T − bp T + p T p (cid:17) − n . (3)The handling of fit and the calculation of the uncer-tainties is identical to Method A. The result is also shownin Fig. 2. In contrast to Method A, which only gradu-ally approaches the asymptotic value at high p T , MethodB reaches the constant at p T of about 5 GeV/ c and ata lower R ∞ = 0 . ± .
024 value, which will be usedas a reference throughout this article. We note that thechange to the constant value is rather abrupt.
C. Gaussian Process Regression (GPR)
Both previous methods have a built-in assumption thatthe η/π has a constant asymptotic value at high p T .However, the data suggest that there might be a maxi-mum around 8 GeV followed by a decrease towards higher p T . In order to avoid any assumptions about the shapewe resort to a machine learning technique called GaussianProcess Regression (GPR), which possesses no physicalknowledge but gives full trust to the data it is given. De-tails about the GPR can be found in [18], and commentsabout the specific implementation we use are summarisedin Appendix A. In general the GPR works best in theregion where many consistent data points are available.Less data points or inconsistent data sets lead to larger [GeV/c] T p -
10 1 10 p / h FIG. 2. Data for the η/π ratio from p + p and p +A collisionscompared to three different methods to describe the data witha universal shape: empirical fit A, empirical fit B, and GPR. [GeV/c] T p p / h Maximal cover of the three p+p methodsp+p 200 GeV, PYTHIA v6.1 )=0.476 ¥ ( p / h scaling, p+p 7 TeV, T m )=0.5 ¥ ( p / h scaling, p+p 200 GeV, T m FIG. 3. Result of combining the three empirical methodsto one universal estimate of η/π as function of p T . Alsoshown for reference are the estimates based on the m T scalinghypothesis and the result of a pythia calculation, both fromFig. 1. uncertainties, and unlike the fitting methods the GPRcan not reliably extrapolate much beyond the range cov-ered by data.The result of the GPR is presented in Fig. 2, withthe band indicating the uncertainties. Over most of the p T range the GPR gives an equally good description ofthe data compared to Methods A and B. As expected, itfollows the data and peaks near 8 GeV/ c . Towards higher p T η/π from the GPR decreases. Whether the drop athigh p T is physical or an artefact of different data setswith different p T ranges not being perfectly consistent inthe range from 3 to 10 GeV/ c will only be resolved withmore precise data. TABLE I. References and systems quoted in this article arecollected in this table. For each A+A system, if differentcentralities have different p T ranges, the one of the minimumbias is presented.System Experiment √ s NN p T range [GeV/c] Ref.p+p CERN WA70 23 GeV 4 − . − . . − . − . . − . −
18 [15]Cu+Au PHENIX 200 GeV 2 . −
19 [21]U+U PHENIX 192 GeV 2 . −
13 [22]d+Au PHENIX 200 GeV 2 . −
11 [12]Au+Au PHENIX 200 GeV 2 . − . . −
17 [23]Pb+Pb ALICE 2.76 TeV 1 . − . Since we do not know the correct functional form of η/π , in particular at high p T , we combine the resultsobtained with the three methods as our best estimate fora universal η/π ratio for p + p and p +A collisions. Thisis achieved by assigning every p T value the minimum ofthe lower uncertainty range of the three methods as thelower bound and the maximum as the upper bound. Theaverage of the lower and upper bound is used as centralvalue. In the following we will use ( η/π ) mcpp to refer tothis combined result, with the superscript mc referring tomaximal coverage of uncertainties. The result is given inFig. 3 and compared to the m T -scaling prediction as wellas the pythia calculation already shown in Fig. 1. Onecan see that all of the theoretical predictions overestimatethe ratio for p T below 3-4 GeV/ c . IV. UNIVERSALITY OF η/π RATIO SYSTEMSAT HIGH p T In the previous section we established that the η/π ratios measured in p + p and p +A collisions are consis-tent with being constant at high p T with a value of R ∞ = 0 . ± .
024 (Section III.B). Here we demonstratethat all available data from p + p , p +A, and A+B colli-sions listed in Table I are consistent with this R ∞ valueindependent of the collision energy, collision system, orcollisions centrality.For this demonstration we adopt the functional formfrom Eq. 3 (empirical Fit B). The parameters are fixedusing the simultaneous fit to the p + p and p +A data tothe following values: a = − . b = 0 . p = 4 . [GeV] NN s10 ) ¥ = T ( p p / h Cu+Au, 200 GeV, M.B., PRC 98, 054903, 2018U+U, 192 GeV, M.B., arXiv:2005.14686, 2020Au+Au, 200 GeV, M.B., PRC 87, 034911, 2013d+Au, 200 GeV, M.B., PRC 75, 024909, 2007Au+Au, 200 GeV, M.B., PRC 75, 024909, 2007p+Pb, 5.02 TeV, EPJ C 78, 624, 2018p+Au, 29.1 GeV, EPJ C 4, 249, 1998p+Be, 29.1 GeV, EPJ C 4, 249, 1998p+p, 23 GeV to 8 TeV) of empirical fit B ¥ R( – )=0.487 ¥ ( p / h FIG. 4. Values of R ∞ = η/π ( p T → ∞ ) as a function of √ s NN for the minimum bias p + p , p +A and A+B data sets.Statistical errors are shown as bars, systematic uncertaintiesas bands. Also shown is a band representing 0 . ± . p + p and p +A data. Note that the A+B data at 200 GeV are offset in √ s NN to avoid overlap of data sets. n = 5 .
07, and the composite parameter R ∞ /A = 2 . R ∞ is determined individually foreach data set using the data shuffling method. For eachdata set we vary the points many times within their sys-tematic uncertainties, as discussed in Appendix B, andcreate an ensemble of R ∞ and σ R ∞ values. The meanof the R ∞ ensemble is used as the measurement of R ∞ for η/π and the standard deviation is quoted as the sys-tematic uncertainty. The mean of the σ R ∞ ensemble isquoted as the statistical uncertainty.Fig. 4 shows the results as a function of the nucleon-nucleon center of mass energy √ s NN for the minimumbias data samples of all collision systems. Also shown onthe figure is the R ∞ value obtained from the combined fitto the p + p and p +A data sets using method B. Withinuncertainties all data sets are consistent with this valueand there is no evidence for a √ s NN dependence of R ∞ .For most publications of η/π from heavy ion colli-sions, the data was also presented for centrality selectedevent classes. In order to include these in the compari-son, we plot R ∞ as a function of the number of producedparticle dN ch /dη (cid:12)(cid:12) η =0 . The dN ch /dη values used are sum-marized in Tab. II.The results are given in Fig. 5. Again all values are con-sistent with a universal value within uncertainties. Thisanalysis strongly suggest that R ∞ does not depend onthe collision systems, √ s NN , or the centrality of the col-lisions and that any apparent differences are likely dueto systematic effects specific to individual data sets. /dy ch dN1 10 ) ¥ = T ( p p / h Cu+Au, 200 GeV, PRC 98, 054903, 2018U+U, 192 GeV, arXiv:2005.14686, 2020Au+Au, 200 GeV, PRC 87, 034911, 2013d+Au, 200 GeV, PRC 75, 024909, 2007Au+Au, 200 GeV, PRC 75, 024909, 2007Pb+Pb, 2.76 TeV, PRC 98, 044901, 2018p+Pb, 5.02 TeV, EPJ C 78, 624, 2018p+Au, 29.1 GeV, EPJ C 4, 249, 1998p+Be, 29.1 GeV, EPJ C 4, 249, 1998p+p, 23 GeV to 8 TeV) of empirical fit B ¥ R( – )=0.487 ¥ ( p / h FIG. 5. Values of R ∞ = η/π ( p T → ∞ ) as a function of dN ch /dη . The presentation is identical to Fig. 4, however,for A+B collisions results from different centrality classes areshown rather than for the minimum bias sample. V. THE EFFECT OF RADIAL FLOW
We have shown that η/π can be described by onecommon function for all p + p and p +A collisions over themeasured p T range from 0.1 to 20 GeV/ c . Furthermore,above p T =5 GeV/ c the same function describes all datafrom heavy ion collisions. Whether this universal func-tion also describes heavy ion data at lower p T can not betested due to the absence of accurate experimental data.However, there are reasons to believe that this universal-ity does not hold at low p T .Evidence for strong collective motion of the bulk ofthe produced particles has been observed in all high en-ergy heavy ion collision. This motion is consistent witha Hubble like hydrodynamic expansion of the collisionvolume, with a linear velocity profile in radial direction.In this velocity profile heavier particles gain more mo-mentum than lighter ones. Radial flow effectively de-pletes the particle yields at low p T and enhances themin an intermediate p T range, which is determined by themass of the particle. For p T much larger than the par-ticle’s mass radial flow becomes negligible. Fig. 6 showsthe effect schematically by comparing π, K , and p spec-tra at RHIC energies. The spectra shown are roughlyto scale and consistent with experimental data from 200GeV Au+Au collisions. They are normalized per particleof the corresponding type at mid rapidity.Since the η meson has about the same mass as thekaon, one would expect that in the momentum rangefrom a few hundred MeV/ c to a few GeV/ c radial flowincreases the yield of η mesons significantly more thanthat of π . This in turn would increase the η/π ratioin heavy ion collisions compared to that observed in p + p and p +A collisions.To quantify the size of the modification due to radial TABLE II. Values for d N ch / d η at mid-rapidity for all colli-sion systems and centrality selections used in this work. Forp+p collisions, the numbers correspond to the inelastic p + p cross section as given in [25]. For all other cases, whenever areference is given, the values are taken directly for the publi-cation. For PHENIX data we use data tabulated in [26]. Thesymbol * in the reference indicates that the value was extrap-olated beyond what was tabulated. All minimum bias values(MB) that are marked by ** were calculated from the central-ity selected data sets for the same system. For all data theuncertainties were calculated assuming that the values quotedin the reference are fully correlated. Reference [6] does notgive an uncertainty on the multiplicity value.System √ s NN Centrality d N ch / d η Ref.p+p √ s – α ( √ s/ GeV) δ [25]p+Au 29.1 GeV – 4.7 [6]p+Be 29.1 GeV – 3.0 [6]p+Pb 5.02 TeV – 16 . ± . −
20% 17 . ± . −
40% 12 . ± . −
60% 8 . ± . −
88% 1 . ± . . ± . −
20% 268 ±
20 [26]20% −
40% 131 ±
10 [26]40% −
60% 54 ± −
93% 12 . ± . ± −
20% 519 ±
26 [26]20% −
60% 156 ±
11 [26]60% −
92% 16 . ± ±
11 **U+U 192 GeV 0% −
20% 636 ±
51 [26]20% −
40% 268 ±
21 [26]40% −
60% 79 ± −
80% 18 . ± ±
35 **Pb+Pb 2.76 TeV 0% −
10% 1448 ±
55 [28]20% −
50% 445 . ±
10 [28] flow we will use a double ratio R flow defined as follows: R flow ≡ (cid:0) ηπ (cid:1) C i (cid:0) ηπ (cid:1) p + p ≈ (cid:16) K ± π ± (cid:17) C i (cid:16) K ± π ± (cid:17) p+p ≡ (cid:16) R K ± AA (cid:17) C i (cid:0) R π ± AA (cid:1) C i , (4)where we take advantage of the fact the momentum boostfrom radial flow is mostly determined by the particlemass and that m K ± ≈ m η . Also charged pions are usedinstead of neutral pions, since π ± and kaons are typicallymeasured simultaneously with the same detector systemsand thus most systematic uncertainties on the measure-ment cancel in the double ratio. The subscript C i refersto a specific collision system, energy and centrality selec- [GeV/c] T p [ a . u .] T d N / dp -9 -8 -7 -6 -5 -4 -3 -2 -1 Hadron spectra at RHICAu+Aup+p p K /10p /100
FIG. 6. Schematic comparison of π, K, p spectra from Au+Auand p+p collisions at √ s NN =200 GeV. All spectra are ap-proximately normalize to their rapidity density at mid rapid-ity. Different particle types are separated by factors of 10 forclarity. [GeV/c] T p ) –p ( AA ) / R – ( K AA = R f l o w R - - Baseline="1.0"-0.0, 0-10% flow
R -0.2, 10-20% flow
RBaseline="1.0"-0.0, 0-10% flow
R -0.2, 10-20% flow
R -0.4, 20-40% flow
R -0.6, 40-60% flow
R -0.8, 60-92% flow
R -0.4, 20-40% flow
R -0.6, 40-60% flow
R -0.8, 60-92% flow R by PHENIX Au+Au 200 GeV – p and – From K
FIG. 7. Double ratio, obtained by (cid:0) R KAA (cid:1) C i / ( R πAA ) C i . tion.Fig. 7 presents R flow for different centrality classes ofAu+Au collisions at 200 GeV. The values were calculatedfrom data published by PHENIX [29]. The data coverthe p T range from 0.5 to 2 GeV/ c and GPR is used toextrapolate somewhat beyond the measured range. Ac-cording to this estimate the η/π ratio is enhanced incentral collisions in a p T region from 0.4 to 3 GeV/ c witha maximum of about 25% near 1 GeV/c. The enhance-ment is reduced for more peripheral collisions and nearlyvanishes for the 60-92% selection.In Fig. 8 we depict the estimate for R flow for Pb+Pb [GeV/c] T p1 10 f l o w R ), 0-10% – p ( AA ) / R – (K AA R , 0-10% ppmc ) p / h / ( PbPb ) p / h ( ALICE Pb+Pb 2.76 TeV
FIG. 8. Double ratio (the flow ratio) for K ± /π ± and η/π inPb+Pb collisions at 2.76 TeV. data at 2.76 TeV calculated from K and π ± data mea-sured by ALICE [24, 30]. Shown are results for a 0-10%centrality selections. Only statistical uncertainties areshown. The flow effect is significantly larger at LHCthan at RHIC: the p T range affected is extended to 5-6 GeV/ c , it reaches its maximum at higher p T around 3GeV/ c , and the maximum has increased to about 50%.All indicates that radial flow effects increase with beamenergy, which is consistent with a higher initial pressureand a longer lifetime of the system at the LHC comparedto RHIC.ALICE also has published η/π for Pb+Pb collisionsat 2.76 TeV [24] down to 1 GeV/ c , which can be usedto verify the validity of the R flow estimate from K/π .For this we have divided Pb+Pb data by the universal( η/π ) mcpp from Fig. 3. The result is also shown in Fig. 8,error bars represent the combine uncertainty of ( η/π ) mcpp and the statistical uncertainty of ( η/π ) P bP b . The ansatzthat ( K ± /π ± ) cent ( K ± /π ± ) pp ≈ ( η/π ) cent ( η/π ) pp is consistent with the data.To construct an η/π ratio for a specific collision sys-tem and centrality selection we modify the universalshape ( η/π ) mcpp determined from p + p and p +A data (seeFig. 3 from section III) with R flow for the selected heavyion sample: (cid:16) ηπ (cid:17) C i = (cid:16) ηπ (cid:17) pp × R flow ≈ (cid:16) ηπ (cid:17) pp × (cid:16) K ± π ± (cid:17) C i (cid:16) K ± π ± (cid:17) pp . (5)Since R flow may be available only in a limited p T re-gion, for example from 0.4 to 2 GeV/ c in Fig. 7, wepropose the following procedure that can be applied to any A+B collisions system if π ± and K data are avail-able for the p T range affected by radial flow. In thefirst step we create pseudo data for η/π by multiply-ing R flow point-by-point with the ( η/π ) mcpp up to p cutT where R flow ( p T ) ≈
1. This range extends to 1.4 or 4.5GeV/c for Au+Au at 200 GeV at 60-92% centrality andPb+Pb at 2.76 TeV, respectively. To ensure that ourflow estimate has the correct asymptotic behavior weadd a second set of pseudo data with constant valuesof η/π = 0 . ± . . × p cutT , which ever islarger. The combine pseudo data are processed througha GPR to obtain a smooth curve. Finally, in order toaccount appropriately for the systematic uncertainties athigh p T we merge the GPR describing the flow effect with( η/π ) mcpp above p cutT . The uncertainty band at low p T isalso taken to be whichever is larger.In Fig. 10 the construction η/π is presented step bystep for three examples: 0-20%, 60-92% Au+Au at 200GeV and 0-10% Pb+Pb at 2.76 TeV, in panels (a) to (c)respectively. The pseudo data generated are representedby points, which are then processed through a GPR re-sulting in the hashed green bands. They are contrastedwith ( η/π ) mcpp , the blue band, and merged with it above p cutT to create the final green envelope representing our η/π estimates. As discussed above the largest flow ef-fect is observed for central Pb+Pb collisions at the LHC(panel (c)). For central Au+Au collison at RHIC (panel(a)) a much smaller effect is observed, and finally periph-eral collisions of the same system are consistent with noflow effect (panel (b)).These best estimates are compared to data in Fig. 9.For the comparison we selected data sets with similarcharged particle densities, so that despite the differencein collision system, centrality or √ s NN matter was cre-ated under similar conditions and evolved the same waywith time. In all three cases our best estimates are con-sistent with the η/π data. VI. SUMMARY DISCUSSION
We find a universal p T dependence of η/π for all p + p and p +A collisions independent of the center of mass en-ergy from √ s NN =23 GeV to 8 TeV. We note that likeoriginally discovered in [6], below 3 GeV/ c the univer-sal ratio is significantly below m T scaling extrapolationsfrom higher p T .That there is no √ s NN dependence is surprising asthe p T spectra of all particles vary strongly with √ s NN and particle production from jet fragmentation becomesincreasingly prevalent at higher energies. None-the-lessthere seems to be no impact on the relative yield at which η and π are produced. This may hint at a largely uni-versal hadronisation process in which hadrons are alwayscreated under the same conditions, even if the underly-ing mechanism is considered different, for example bulkparticle production or jet fragmentation. [GeV/c] T p p / h pseudo dataGPR result of pseudo data in p+p collisions p / h maximal cover of (c) [GeV/c] T p p / h pseudo data for 60-92%GPR result of pseudo data in p+p collisions p / h maximal cover of (b) [GeV/c] T p p / h pseudo data for 0-10%GPR result of pseudo data in p+p collisions p / h maximal cover of (a) FIG. 9. Estimate of the effect of radial flow on η/π for 0-20%, 60-92% Au+Au collisions at 200 GeV and 0-10% Pb+Pbcollisions at 2.76 TeV. Details are discussed in the text. For heavy ion collisions, η/π has the same univer-sal behavior at high p T , independent of collision species,collision energy, or collision centrality. For lower p T wefind evidence for modifications of the relative particleyields due to radial flow. One might speculate that thesame universal harmonization process is at work but thathadrons are produced in a moving reference frame.We have quantified the modification of the η/π ratio due to radial flow using the double ratio R AA ( K ) /R AA ( π ). This assumes that the change of the p T spectra depends entirely on the particle mass, but it [GeV/c] T p p / h T ; maximal cover for high p T GPR(pseudo data) for low pPb+Pb, 2.76 TeV, 0-10%, PRC 98,044901, 2018 (c) [GeV/c] T p p / h T ; maximal cover for high p T GPR(pseudo data) for low pp+Pb, 5.02 TeV, EPJC 78, 624, 2018Au+Au, 200 GeV, 60-93%, PRC 87, 034911, 2013 d+Au, 200 GeV, 0-20%, PRC 87, 034911, 2013U+U, 192 GeV, 60-80%, PRC 102, 064905, 2020Cu+Au, 200 GeV, 60-90%, PRC 98, 054903, 2018 (b) [GeV/c] T p p / h T ; maximal cover for high p T GPR(pseudo data) for low pPb+Pb, 2.76 TeV, 20-50%, PRC 98,044901, 2018Au+Au, 200 GeV, 0-20%, PRC 87, 034911, 2013U+U, 192 GeV, 0-20%, PRC 102, 064905, 2020Cu+Au, 200 GeV, 0-20%, PRC 98, 054903, 2018 (a)
FIG. 10. Estimate of the effect of radial flow on η/π for 0-20%, 60-92% Au+Au collisions at 200 GeV and 0-10% Pb+Pbcollisions at 2.76 TeV. Details are discussed in the text. does not make any assumptions about the similarity of η and kaon spectra themselves. We note that our ap-proach may overestimate the modification due to flow,since kaon production or generally strange quark produc-tion is enhanced in heavy ion collisions. In our estimatethe modification increases with √ s NN . At 200 GeV atRHIC the maximum increase of η/π is estimated to be25% around 1 GeV/ c , in contrast at 2.76 TeV at the LHCthe maximum increase is nearly 50% and occurs at higher p T between 2 and 3 GeV/c.With our original motivation in mind, which was toreduce systematic uncertainties on the measurement ofdirect photons, we proposed a new methodology to cre-ate η/π ratios. This method is more accurate than fre-quently used extrapolations to lower p T based on m T scaling, and does not suffer from the frequent lack ofstatistics for η measurements. Our approach can be ap-plied to all systems for which K/π ± is measured in the p T range affected by radial flow. The method does notrequire actual measurements of η production for a givensystem.We have tested this method for two specific collisionsystems. For Au+Au collisions at 200 GeV the devia-tions due to flow are found to be within ±
15% of theminimum bias values. Even for central collisions the η/π underlying the estimate of photon from hadron de-cays used in direct photon measurements [2–4] is abovewhat we propose. As a consequence direct photon yieldshave been slightly under estimated, though the differ-ences are within quoted systematic uncertainties. Forcentral Pb+Pb collisions at 2.76 TeV the flow modifica-tions are larger and coincidentally bring η/π much closerto the m T scaling assumption used in the measurementof direct photons published by ALICE [5]. ACKNOWLEDGMENTS
We acknowledge the support from the Office of NuclearPhysics in the Office of Science of the Department ofEnergy.
Appendix A: Gaussian Process Regression
In this section we discuss the implementation of theGaussian Process Regression (GPR) used in our analy-sis. Full details about the GPR can be found in [18].We start with a selection of N data points x i , y i , and σ i . In our case this is typically, N values of log 10( p T ), η/π ( p T ) and its variance. We use a Square-Exponential(SE) kernel to describe the correlation between points,which is given by: k SE ( x i , x j ) = σ p exp (cid:18) − ( x i − x j ) l (cid:19) . (A1)Here σ p gives the strength of the correlation betweeny values and l is a length scale that determines the rangein x over which y values are correlated.We introduce the vectors X and Y , which have dimen-sion N and elements x i and y i , i.e. the data. The corre-lations between y values is then defined by a covariancematrix K xx which has the elements( K xx ) ij = k SE ( x i , x j ) + δ ij σ i , (A2)with δ ii = 1 and δ ij = 0 for i (cid:54) = j . The term δ ij σ i addsnoise to the diagonal elements to account for the uncer- tainty on the measured y values. In order to determine σ p and l , we maximize the log likelihood function:log p ( Y | σ p , l ) = − n π − Y T [ K xx ] − Y −
12 log det( K xx ) . (A3)Once the parameters σ p and l are set, we can predict yvalues for any given x value. For this we introduce a vec-tors X ∗ and Y ∗ of dimension R and elements x ∗ i for whichwe want to predict y ∗ i , with R typically much larger thanN. We introduce two more matrices, one of dimension R × N with elements ( K x ∗ x ) ij ≡ k SE ( x ∗ i , x j ), and one ofdimension R × R with elements ( K x ∗ x ∗ ) ij ≡ k SE ( x ∗ i , x ∗ j ).The predicted values Y ∗ and their covariance matrixCov( Y ∗ ) are then calculated as follows: Y ∗ = K x ∗ x [ K xx ] − Y, (A4)Cov( Y ∗ ) = K x ∗ x ∗ − K x ∗ x [ K xx ] − K Tx ∗ x . (A5)The diagonal elements of Cov( Y ) give the variance of Y ∗ due to the statistical uncertainty on the data Y . Werefer to this as vector S ∗ stat . We also consider the fit un-certainty on σ p and l . The variance can be calculated bythe covariance matrix M the fitting procedure providesthrough error propagation of Eq. A4. : S fit =( ∂ l y ∗ ) M ll + 2 ∂ l y ∗ ∂ σ p y ∗ M lσ p + ( ∂ σ p y ∗ ) M σ p σ p , ∀ y ∗ ∈ Y ∗ (A6)with ∂ l and ∂ σ f being the partial derivatives of Y ∗ withrespect to l and σ p .In addition, we incorporate the systematic uncertain-ties using the data shuffling method discussed in Ap-pendix B. We create a large ensemble of different { Y ∗ λ } for the same X ∗ by varying each data set by a Gaussianrandom number (cid:15) ∼ N (0 ,
1) multiplying systematic un-certainties. The pointwise variance of ensambles { Y ∗ λ } ,which we call S ∗ sys , is used as measure of the systematicuncertainty.In all figures that show results from the GPR the centerline represents Y ∗ and the vertical width of the band is (cid:113) S stat + S fit + S ∗ Sys , pointwise.
Appendix B: Data-shuffling method
The data-shuffling method is a Monte Carlo simula-tion approach that allows to estimate the effect of sys-tematic uncertainties on the result of a fit of a func-tion to data. To illustrate how the method works wefirst consider the case of one data set and assume thatthe systematic uncertainties are fully correlated. Herefully correlated means that the correlation matrix is ρ ij = 1 , ∀ i, j . Suppose each data point is described by a4-tuple ( x i , y i , σ stati , σ sysi ). One first defines a Gaussianrandom variable (cid:15) ∼ N (0 , y by a small quantity to y (cid:48) i = y i + σ sysi (cid:15) ac-cordingly. Then in each simulation, one fits with theseshifted data, and gets one fit result. This is repeated L times, which generates L sets of fit parameters. Foreach set of fit parameters one can devide the x valuesinto R bins. Both L and R are usually large numbers.This results in a L -by- R matrix of y λr values. For a fixed r , the mean and standard deviation of { y λr : 1 ≤ λ ≤ L } are calculated. The standard deviation is assigned as systematic uncertainty of the fit for the given r .The method is expanded to multiple data sets by gen-erating independent Gaussian random variables for eachdata set. In principle, more complex correlations of un-certainties for an individual data set can be decoded in ρ ij , however, for the data at hand these correlations arenot known and thus can not be implemented.One can choose as the final y value for a given r ei-ther the mean from data-shuffling, or the fit result of theoriginal data (i.e., the fit result when the Gaussian vari-ables are zero). The difference between them is usuallynegligible. [1] E. V. Shuryak, Quark-Gluon Plasma and Hadronic Pro-duction of Leptons, Photons and Psions, Sov. J. Nucl.Phys. , 408 (1978).[2] A. Adare et al. (PHENIX), Enhanced production of di-rect photons in Au+Au collisions at √ s NN = 200 GeVand implications for the initial temperature, Phys. Rev.Lett. , 132301 (2010), arXiv:0804.4168 [nucl-ex].[3] A. Adare et al. (PHENIX), Centrality dependence oflow-momentum direct-photon production in Au+Au col-lisions at √ s NN = 200 GeV, Phys. Rev. C , 064904(2015), arXiv:1405.3940 [nucl-ex].[4] A. Adare et al. (PHENIX), Beam Energy and Central-ity Dependence of Direct-Photon Emission from Ultra-relativistic Heavy-Ion Collisions, Phys. Rev. Lett. ,022301 (2019), arXiv:1805.04084 [hep-ex].[5] J. Adam et al. (ALICE), Direct photon production in Pb-Pb collisions at √ s NN = 2.76 TeV, Phys. Lett. B ,235 (2016), arXiv:1509.07324 [nucl-ex].[6] G. Agakichiev et al. , Neutral meson production in p Beand p Au collisions at 450-GeV beam energy, Eur. Phys.J. C , 249 (1998).[7] L. Altenk¨amper, F. Bock, C. Loizides, and N. Schmidt,Applicability of transverse mass scaling in hadroniccollisions at energies available at the CERN LargeHadron Collider, Phys. Rev. C , 064907 (2017),arXiv:1710.01933 [hep-ph].[8] A. Adare et al. (PHENIX), Detailed measurement of the e + e − pair continuum in p + p and Au+Au collisions at √ s NN = 200 GeV and implications for direct photon pro-duction, Phys. Rev. C81 , 034911 (2010), arXiv:0912.0244[nucl-ex].[9] R. Hagedorn, Statistical thermodynamics of strong in-teractions at high-energies, Nuovo Cim. Suppl. , 147(1965).[10] S. Acharya et al. (ALICE), π and η meson productionin proton-proton collisions at √ s = 8 TeV, Eur. Phys. J.C , 263 (2018), arXiv:1708.08745 [hep-ex].[11] B. Abelev et al. (ALICE), Neutral pion and η me-son production in proton-proton collisions at √ s = 0 . √ s = 7 TeV, Phys. Lett. B , 162 (2012),arXiv:1205.5724 [hep-ex].[12] S. Adler et al. (PHENIX), High transverse momentum η meson production in p + p , d + Au and Au+Au collisionsat √ s NN = 200-GeV, Phys. Rev. C , 024909 (2007),arXiv:nucl-ex/0611006.[13] S. Acharya et al. (ALICE), Production of π and η mesons up to high transverse momentum in pp colli-sions at 2.76 TeV, Eur. Phys. J. C , 339 (2017), arXiv:1702.00917 [hep-ex].[14] A. Adare et al. (PHENIX), Cross section and double he-licity asymmetry for η mesons and their comparison toneutral pion production in p + p collisions at √ s =200GeV, Phys. Rev. D , 032001 (2011), arXiv:1009.6224[hep-ex].[15] S. Acharya et al. (ALICE), Neutral pion and η mesonproduction in p-Pb collisions at √ s NN = 5 .
02 TeV, Eur.Phys. J. C , 624 (2018), arXiv:1801.07051 [nucl-ex].[16] S. Adler et al. (PHENIX), Mid-rapidity neutral pionproduction in proton proton collisions at √ s = 200-GeV, Phys. Rev. Lett. , 241803 (2003), arXiv:hep-ex/0304038.[17] T. Sjostrand, P. Eden, C. Friberg, L. Lonnblad, G. Miu,S. Mrenna, and E. Norrbin, High-energy physics eventgeneration with PYTHIA 6.1, Comput. Phys. Commun. , 238 (2001), arXiv:hep-ph/0010017.[18] C. Rasmussen and C. Williams, Gaussian Processes forMachine Learning (The MIT Press, 2005).[19] M. Bonesini et al. (WA70 Collaboration), High transversemomentum η production in π − p, π + p, and pp interac-tions at 280 gev/c, Z. Phys. C , 527 (1989).[20] L. Apanasevich et al. (Fermilab E706), Production of π and η mesons at large transverse momenta in pp and pBe Interactions at 530 and 800 GeV/c., Phys. Rev. D , 052001 (2003), arXiv:hep-ex/0204031.[21] C. Aidala et al. (PHENIX), Production of π and η mesons in Cu+Au collisions at √ s NN =200 GeV, Phys.Rev. C , 054903 (2018), arXiv:1805.04389 [hep-ex].[22] U. Acharya et al. (PHENIX), Production of π , η , and K S mesons in U+U collisions at √ s NN = 192 GeV, Phys.Rev. C , 064905 (2020), arXiv:2005.14686 [hep-ex].[23] A. Adare et al. (PHENIX), Neutral pion production withrespect to centrality and reaction plane in Au+Au col-lisions at √ s NN =200 GeV, Phys. Rev. C , 034911(2013), arXiv:1208.2254 [nucl-ex].[24] S. Acharya et al. (ALICE), Neutral pion and η me-son production at mid-rapidity in Pb-Pb collisions at √ s NN = 2.76 TeV, Phys. Rev. C , 044901 (2018),arXiv:1803.05490 [nucl-ex].[25] J. Adam et al. (ALICE), Charged-particle multiplicitiesin proton–proton collisions at √ s = 0 . , 33 (2017), arXiv:1509.07541 [nucl-ex].[26] A. Adare et al. (PHENIX), Transverse energy productionand charged-particle multiplicity at midrapidity in vari-ous systems from √ s NN = 7 . , 024901 (2016), arXiv:1509.06727 [nucl-ex].[27] B. Abelev et al. (ALICE), Pseudorapidity density of charged particles in p + Pb collisions at √ s NN =5 .
02 TeV, Phys. Rev. Lett. , 032301 (2013),arXiv:1210.3615 [nucl-ex].[28] K. Aamodt et al. (ALICE), Centrality dependence of thecharged-particle multiplicity density at mid-rapidity inPb-Pb collisions at √ s NN = 2 .
76 TeV, Phys. Rev. Lett. , 032301 (2011), arXiv:1012.1657 [nucl-ex].[29] A. Adare et al. (PHENIX), Spectra and ratios of identified particles in Au+Au and d +Au collisions at √ s NN = 200 GeV, Phys. Rev. C , 024906 (2013),arXiv:1304.3410 [nucl-ex].[30] J. Adam et al. (ALICE), Centrality dependence of thenuclear modification factor of charged pions, kaons, andprotons in Pb-Pb collisions at √ s NN = 2 .
76 TeV, Phys.Rev. C93