Event patterns extracted from transverse momentum and rapidity spectra of Z bosons and quarkonium states produced in pp and Pb-Pb collisions at LHC

Transverse momentum ($p_T$) and rapidity ($y$) spectra of $Z$ bosons and quarkonium states (some charmonium $c\bar c$ mesons such as $J/\psi$ and $\psi(2S)$, and some bottomonium $b\bar b$ mesons such as $\Upsilon(1S)$, $\Upsilon(2S)$, and $\Upsilon(3S)$) produced in proton-proton ($pp$) and lead-lead (Pb-Pb) collisions at the large hadron collider (LHC) are uniformly described by a hybrid model of two-component Erlang distribution for $p_T$ spectrum and two-component Gaussian distribution for $y$ spectrum. The former distribution results from a multisource thermal model, and the latter one results from the revised Landau hydrodynamic model. The modelling results are in agreement with the experimental data measured in $pp$ collisions at center-of-mass energies $\sqrt{s}=2.76$ and 7 TeV, and in Pb-Pb collisions at center-of-mass energy per nucleon pair $\sqrt{s_{NN}}=2.76$ TeV. Based on the parameter values extracted from $p_T$ and $y$ spectra, the event patterns (particle scatter plots) in two-dimensional $p_T$-$y$ space and in three-dimensional velocity space are obtained.


Introduction
High energy nucleus-nucleus collisions at the relativistic heavy ion collider (RHIC) [1][2][3][4] and large hadron collider (LHC) [5][6][7][8] provide excellent environment and condition of high temperature and density [9], where a new state of matter, namely the quark-gluon plasma (QGP) [10][11][12], is expected to form and to live for a very short time. It is regretful that the QGP cannot be directly measured in experiments due to its very short lifetime. Instead, to understand the formation and properties of QGP, the distribution laws of finalstate particles are studied. Because of the complexities and difficulties in experiments, the observables are limited. To obtain more information from limited observables, we need more modelling and theoretical analyses.
Generally, an interacting system of nucleus-nucleus collisions undergoes a few stages which include, but are not limited to, the scattering, color glass condensate, thermaliza-tion, hadronization, chemical equilibrium and freeze-out, kinetic equilibrium and freezeout. The distribution laws of final-state particle transverse momenta and rapidities reflect the situation of interacting system at the stage of kinetic freeze-out, while the feed-down corrected yields and the ratios of those yields of different particles reflect the situation at the stage of chemical freeze-out. Based on the descriptions of transverse momentum (p T ) and rapidity (y) spectra, one can extract some information of transverse excitation and longitudinal expansion of interacting system. Thus, other information such as the event pattern (particle scatter plot) in two-dimensional p T -y space and three-dimensional velocity space can be extracted from parameters fitted to p T and y spectra [13,14].
In peripheral nucleus-nucleus collisions, less nucleons take part in the interactions. The most peripheral nucleus-nucleus collisions contain only two nucleons which are from the two collision nuclei respectively. Proton-proton (pp) collisions are similar to the most peripheral nucleus-nucleus collisions in the case of neglecting the spectator (cold nuclear) effect. As an input quantity and a basic collision process, pp collisions can be used to give comparisons with nucleus-nucleus collisions. We are interested in both pp collisions and lead-lead (Pb-Pb) collisions at the LHC.
To understand the stage of kinetic freeze-out in high energy collisions, we can analyze the p T and y spectra to obtain the probability density functions f p T (p T ) for p T and f y (y) for y. Nevertheless, these probability density functions cannot directly give us a whole and perceptual picture of the interacting system at the stage of kinetic freezeout. In fact, a whole and perceptual picture can help us understand the interacting mechanisms in detail. Fortunately, we can use the Monte Carlo method to extract some discrete values of p T and y based on f p T (p T ) and f y (y). Other quantities such as energy, momentum components, velocity, and velocity components can be obtained according to some definitions and assumptions. Because f p T (p T ) for p T and f y (y) for y are based on descriptions of experimental spectra, the extracted discrete values are independent of models.
In this paper, based on a hybrid model of two-component Erlang distribution for p T spectrum (which results from a multisource thermal model [15][16][17]) and two-component Gaussian distribution for y spectrum (which results from the Landau hydrodynamic model and its revisions [18][19][20][21][22][23][24][25][26]), we analyze together p T and y spectra of Z bosons and quarkonium states (some charmonium cc mesons such as J/ψ and ψ(2S), and some bottomonium bb mesons such as Υ(1S), Υ(2S), and Υ(3S)) produced in pp collisions at center-of-mass energies √ s = 2.76 and 7 TeV, and in Pb-Pb collisions at center-of-mass energy per nucleon pair √ s N N = 2.76 TeV. The modelling results are in agreement with the experimental data measured at the LHC. Based on the parameters extracted from p T and y spectra, the event patterns (particle scatter plots) at kinetic freeze-out in two-dimensional p T -y space and in three-dimensional velocity space are obtained. The structure of the present work is as followings. The model and method are shortly described in section 2. Results and discussion are given in section 3. In section 4, we summarize our main observations and conclusions.

The model and method
Firstly, we need modelling descriptions of p T and y spectra. In the framework of multisource thermal model [15][16][17], we can obtain an Erlang distribution or a two-, threeor multi-component Erlang distribution to fit p T spectrum. According to the model, many (m) emission sources which stay at the same excitation state are assumed to form in high energy collisions. Each (the i-th) source is assumed to contribute to transverse momentum p T i by an exponential function where p T i denotes the average value of p T i , which results in Eq. (1) to be a probability distribution and ∞ 0 f p T i (p T i )dp T i = 1. The m sources which contribute to p T resulting in an Erlang distribution [17] which is the folding of m exponential distributions and has the average transverse momentum p T = m p T i . In the case of considering the two-component Erlang distribution, we have which has the average transverse momentum p T = k 1 m 1 p T i 1 + (1 − k 1 )m 2 p T i 2 , where k 1 and 1 − k 1 denote the relative contributions of the first and second components which contribute in the low-and high-p T regions respectively, and the subscripts 1 and 2 denote the quantities related to the first and second components respectively. Eqs. (2) and (3) are probability distributions which are normalized to 1. When we compare them with experimental data, normalization constant (N p T ) which is used to fit the data is needed. On y spectrum, we choose the Landau hydrodynamic model and its revisions [18][19][20][21][22][23][24][25][26] which are called the revised Landau hydrodynamic model in the present work. In the model, the interacting system is described by the hydrodynamics. The y spectrum can be described by a Gaussian function [25,26] where σ y denotes the rapidity distribution width and y C denotes the mid-rapidity (peak position). In symmetric collisions, y C = 0 is in the center-of-mass reference frame.
In the case of considering the two-component Gaussian function for y spectrum, we have where k B (1 − k B ), y B (y F ), and σ yB (σ yF ) denote respectively the relative contribution, peak position, and distribution width of the first (second) component which distributes in the backward (forward) rapidity region. In symmetric collisions such as pp and Pb-Pb collisions which are considered in the present work, we have k B = 1−k B = 0.5, y B = −y F , and σ yB = σ yF . As probability distributions, Eqs. (4) and (5) are normalized to 1. When we compare them with experimental data, a normalization constant (N y ) which is used to fit the data is needed. Secondly, we need discrete values of p T and y. The related calculation is performed by a Monte Carlo method. Let r i , r 1i , r 2i , and R 1,2,···,7 denote random numbers in [0,1]. Eqs. (2)-(5) results in for the first component in the low-p T region, or for the second component in the high-p T region, and y = σ yB −2 ln R 3 cos(2πR 4 ) + y B for the first component in the backward rapidity region, or y = σ yF −2 ln R 5 cos(2πR 6 ) + y F for the second component in the forward rapidity region, respectively. It should be clarified that the random numbers used above are independent in [0,1]. Through the conversions Eqs.  (11) are accustomed expressions which result from the Gaussian distribution. If we use y = σ y √ −2 ln R 1 instead of Eq. (9), we obtain an accustomed expression which result from the Rayleigh distribution f (y) = (y/σ 2 y ) exp(−y 2 /2σ 2 y ) which is different from the Gaussian function. The energy E is given by where m 0 denotes the rest mass of the considered particle. The x-, y-, and z-components of momentum p are given by p x = p T cos ϕ, p y = p T sin ϕ, p z = p 2 T + m 2 0 sinh y (13) Figure 1 presents (a)(b) transverse momentum spectra, d 2 σ/(dydp T ), and (c)(d) rapidity spectra, dσ/dy, of Z bosons produced in pp collisions at √ s = 2.76 TeV for the (a)(c) dimuon (µµ) and (b)(d) dielectron (ee) decay channels in |y| < 2.0 and |y| < 1.44 respectively, where σ on the vertical axis denotes the cross-section, and the integral luminosity L int = 5.4 pb −1 . The closed squares represent the experimental data of the CMS Collaboration [27], and the error bars are only the statistical uncertainties. The curves are our results calculated by using the two-component Erlang distribution for p T spectrum and the two-component Gaussian distribution for y spectrum, which are the results of the multisource thermal model [15][16][17] and the revised Landau hydrodynamic model [18][19][20][21][22][23][24][25][26], respectively. The values of free parameters [m 1 , p T i 1 , k 1 , m 2 , p T i 2 , y F (= −y B ) and σ yF (= σ yB )], normalization constants (N p T and N y ), and χ 2 /dof are listed in Tables 1 and  2, where the normalization constant N p T (or N y ) is used to give comparison between the normalized curve with experimental p T (or y) spectrum. One can see that the results calculated by using the hybrid model are in agreement with the experimental data of Z bosons produced in pp collisions at √ s = 2.76 TeV measured by the CMS Collaboration.
In some cases, the values of χ 2 are very large due to very small experimental errors.   Collaboration. It should be noted that the ATLAS collaboration has also published similar experimental results [28,29]. We can do similar fits for the spectra of p T and y, which are not presented in the present work by design due to the similarity.

Figure
Type  Table 3. In the calculation, the contributions of the first and second components in the two-component Gaussian distribution for y spectrum are not distinguished by design. Figure 4 is another example of event pattern in three-dimensional velocity (β TeV and (c)(d) Pb-Pb collisions at √ s N N = 2.76 TeV. The total number of particles for each panel is 1000. The red and blue globules represent the contributions of the first and second components in the two-component Erlang distribution for p T spectrum respectively. Most particles appear in the region of low-|β x | and low-|β y |. The values of root-mean-squares ( β 2 x for β x , β 2 y for β y , and β 2 z for β z ) and the maximum |β x |, |β y |, and |β z | (|β x | max , |β y | max , and |β z | max ) are listed in Table 4. One can see that The event pattern in velocity (or coordinate) space looks like a rough cylinder along oz axis (the beam direction) when it emit Z bosons, and the maximum velocity surface (the surface consisted of the maximum velocities in different directions) is a fat cylinder which has the length being 1.6-2.2 times of diameter. As for Figure 3, the contributions of the first and second components in the two-component Gaussian distribution for y spectrum are not distinguished by design.
In Figure 5,  Tables 3 and 4. One can see that p 2 T ( y 2 ) for charmonium cc mesons is less (greater) than that for Z bosons. Once again, β 2 x ≈ β 2 y ≪ β 2 z , and |β x | max ≈ |β y | max < |β z | max . At the same time, β 2 x and |β x | max ( β 2 y and |β y | max , β 2 z ) for charmonium cc mesons are larger than those for Z bosons, and the two types of particles have nearly the same |β z | max . The event pattern in velocity (or coordinate) space looks like a rough cylinder along the beam direction when it emit charmonium cc mesons, and the maximum velocity surface is a fat cylinder which has the length being 1.2-1.4 times of diameter.  Table 3. Values of the root-mean-squares p 2 T for p T and y 2 for y corresponding to the particle scatter plots in  pp, ψ(2S) 3.40 ± 0.09 2.67 ± 0.05 Figure 9(a) pp, Υ(1S) 5.53 ± 0.11 2.26 ± 0.05 Figure 9(b) pp, Υ(2S) 6.03 ± 0.11 2.13 ± 0.04 Figure 9(c) pp, Υ(3S) 6.49 ± 0.11 2.23 ± 0.04 Figure 12(a) Pb-Pb, Υ(1S) 6.26 ± 0.14 1.88 ± 0.04 Figure 12 In the upper panel of Figure 8, the spectra of (a) p T for Υ(1S), Υ(2S), and Υ(3S) and (b) y in backward (forward) region for Υ(1S), Υ(2S), and Υ(3S) produced in pp collisions at √ s = 7 TeV are presented, where σ iS and B iS on the vertical axis denote the crosssection and branching fractions of Υ(iS) (i = 1, 2, and 3), respectively. The closed circles, squares, and triangles represent the experimental data of Υ(1S), Υ(2S), and Υ(3S), respectively, measured by the LHCb Collaboration [35], and the error bars are the total uncertainties which combine the systematic and statistical uncertainties to be equal to their root-sum-square. The curves for p T and y spectra are our results calculated by using the two-component Erlang distribution and the two-component Gaussian distribution respectively. The values of free parameters, normalization constants, and χ 2 /dof are listed in Tables 1 and 2. One can see that the results calculated by using the hybrid model are in agreement with the experimental data of quarkonium states (bottomonium bb mesons Υ(1S), Υ(2S), and Υ(3S)) measured by the LHCb Collaboration. Table 4. Values of the root-mean-squares β 2 x for βx, β 2 y for βy, and β 2 z for βz, as well as the maximum |βx|, |βy|, and |βz| (|βx|max, |βy|max, and |βz|max) corresponding to the particle scatter plots in Figures 4, 7, 10, and 13.

Figure
Type   The closed symbols represent the experimental data of the ALICE [30][31][32] and LHCb Collaboration [33,34] in different conditions shown in the panels, and the error bars are only the statistical uncertainties. The curves are our results calculated by using the hybrid model.
symbols represent different data measured by the ALICE [30], LHCb [35], and CMS collaborations [36,37] Tables 3 and 4. One can see that p 2 T ( y 2 ) for bottomonium bb mesons is close to that for charmonium cc mesons, and less (greater) than that for Z bosons. Once more, β 2 x ≈ β 2 y ≪ β 2 z , and |β x | max ≈ |β y | max < |β z | max . At the same time, β 2 x and |β x | max ( β 2 y and |β y | max , β 2 z ) for bottomonium bb mesons are close to those for charmonium cc mesons, and larger than those for Z bosons. The three types of particles have nearly the same |β z | max . The event pattern in velocity (or coordinate) space looks like a rough cylinder along the beam direction when it emit  Figures 12 and 13, where the event patterns in p T − y and β x − β y − β z spaces are given respectively. The red squares and globules represent the contribution of    1, 2, and 3), respectively. The closed symbols represent the experimental data of the LHCb [35], ALICE [30], and CMS Collaborations [36,37]. The error bars in (a)(b) are the total uncertainties which combine the systematic and statistical uncertainties to be equal to their root-sum-square, and the error bars in (c)(d) are only the systematic uncertainties. The curves are our results calculated by using the hybrid model.  denotes the yield, and T AA denotes the nuclear overlap function which depends on the collision centrality. The closed symbols represent the (preliminary) experimental data of the ALICE [38] and CMS Collaborations [39,40] in different conditions shown in the panels, and the error bars are only the systematic uncertainties. The curves are our results calculated by using the hybrid model. We can give qualitatively a prediction on the prospects to actually get J/ψ rapidity distributions in Figure 11(b) for the given centrality classes listed in Figure 11(a). From Table 1 we see that p T i 1 and p T i 2 for Figure 11(a) increase slightly with decrease of the centrality. As a non-leading particle, the rapidity distribution of J/ψ does not depend obviously on the centrality. These two characteristics render that in the derived event pattern, the particles will scatter slightly in higher p T region in peripheral collisions. From central to peripheral collisions, the increase in p T is slight due to the slight increases in p T i 1 and p T i 2 .
The most important set of parameters in Table 1 are the mean contributions, p T i 1 and p T i 2 , of each source for p T spectrum in the first and second components. Then, the mean contribution of each source in the first+second components is p T i = p T /[k 1 m 1 + (1 − k 1 )m 2 ]. The mean transverse momenta contributed by the first and second components are m 1 p T i 1 and m 2 p T i 2 respectively. The most important set of parameters in Table 2 are the peak position y F (y B ) and distribution width σ yF (σ yB ) for the forward (backward) region in the two-component Gaussian function for y spectrum. To see the dependences of these parameters on m 0 , Figure 14 shows the relations between (a) p T i ( p T i 1 , p T i 2 ) One can see an obvious increase in p T i , p T i 1 , p T i 2 , p T , m 1 p T i 1 , and m 2 p T i 2 , and a decrease in y F and σ yF , when m 0 changes from ∼3 GeV to ∼90 GeV. This mass scale imposes a difference on the observed data from the presence of a quantum chromodynamics (QCD) hard scale. From Table 1 one can see that in some cases k 1 is less than 0.5, and p T i 1 is greater than p T i 2 . These results does not mean a reversal of the first and second components, or a reversal of the low-and high-p T regions. In fact, the first component that contributes in the low-p T region is determined by both m 1 and p T i 1 , and the second component that contributes in the high-p T region is determined by both m 2 and p T i 2 . The mean transverse momenta, m 1 p T i 1 and m 2 p T i 2 , contributed by the first and second components, respectively, play a determinate action. Generally, we restrict m 1 p T i 1 < m 2 p T i 2 , which ensures the first component contributing in the low-p T region and the second component contributing in the high-p T region. The size of k 1 is not related to the relative sizes of m 1 p T i 1 and m 2 p T i 2 .
We would like to point out that Z bosons and quarkonium states are chosen by us in the present work due to their productions being through hard scattering and at initial state, which are different from those of light flavor particles which are mostly produced through soft excitation and at intermediate state. Although we describe uniformly the presence of a variable mass scale from 3 to 90 GeV, we cannot provide further information in theory. In fact, the method used in the present work is independent of models, which is only based on the description of probability distribution of experimental data. If we do not use formulas, but discrete probabilities based on data, similar event patters can be obtained. In addition, our choices on data sets are not sole. In fact, other data sets which contained both p T and y spectra can be used to structure the event patterns.
It should be noted that in the calculation of χ 2 in Figures 1, 2, 5 ,8 and 11, we have to determine the corresponding relation between a theoretical point and an experimental data in a given bin of a variable. According to ref. [41], there are more than one methods to determine the corresponding relation, which include considerations in variable and frequency respectively. In the present work, we have simply used the method of histograms, i.e. presenting the theoretical point and experimental data in the form of histograms when we calculate χ 2 . In most cases, if the bin is narrow enough, to use the closest theoretical point and experimental data is an approximate alternative method.
It should also be noted that Figures 3,4,6,7,9,10,12,13 are supposed to only give pictorial representations of the particle scatter plots at the stage of kinetic freeze-out of the interacting system. These pictorial representations are based on the experimental data themselves. They are independent of models due to Eqs. (3) and (5) being only parameterizations of the experimental data, and Eqs. (7), (8), (10), and (11) being only stochastic extractions based on Eqs. (3) and (5). Any other event generators that describe simultaneously the spectra of p T and y can give similar pictorial representations. As we know, concrete and ready-made pictorial representations of the particle scatter plots at the stage of kinetic freeze-out given by other event generators are not available at present.
It seems that possible effects from kinematic cuts are not considered for the extracted discrete values, then for the event patterns. As we know, all the experimental data are subject to trigger and event selection restrictions, which have subsequent effects on the kinematic region in which the cross sections are defined. In fact, because our analyses are based on the experimental data themselves, the extracted parameters imply these restrictions. This renders that there are considerations on the impacts of these restrictions on the subsequent event shapes. Because the differences in experimental restrictions for pp and Pb-Pb collisions as well as for spectra of Z bosons and quarkonium states are small, the impacts on the values of parameters, discrete values, and event patterns are small in some details. Generally, these small differences in restrictions do not effect largely the global shapes and contour profiles of events.
(particle scatter plots) for different particles in three-dimensional velocity space or other available spaces based on the transverse momentum and rapidity spectra of considered particles. Because the present analyses are based on the experimental data themselves, the extracted parameters imply experimental restrictions. This renders that the discrete values and event patterns extracted by us are accurate to the best of our ability. It is expected that different particles correspond to different event patterns in different shapes with different sizes. To give a better comparison, the same experimental restrictions for different particles are needed in the future.