Simulation of ATLAS Collaboration’s Measurement Results of Jet P_T Correlations in Pb + Pb Collisions at √S_NN= 2.76 TeV

Simai Jia¹, Zhanzhi Zhang²

¹School of Physics, Peking University, Beijing, China.
²Second High School Attached to Beijing Normal University International Campus, Beijing, China.

DOI: 10.4236/ns.2019.1110031                

Abstract

In this article, a model based on Glauber Monte Carlo is built to simulate the procedure of jet quenching in Quark-Gluon Plasma (QGP). In this model, energy loss of jets in QGP is parametrized by two quantities: path length of jets in QGP, L, and the initial transverse momentum of the jet P_T. The path length of each branch of the jet are labeled L₁ and L₂. As input data, original jet energy data of p + p collisions were obtained from CMS measurement. After being processed by our model, simulated Pb + Pb jet energy data could be given and were compared to the data of ATLAS’s measurement in Ref. [1]. Distributions of (1/N)dN/dx_J  where x_J=P_T2/P_T1, also noted as “frequency”, are presented as a function of P_T1 and collision centrality. As the final result, two different forms of energy loss formula were found, both of which have good adaptation to certain centrality and P_T ranges.

Keywords

Glauber Monte Carlo Model, Quark-Gluon Plasma, Jet Quenching, Energy Loss Simulation

1. Introduction

Relativistic heavy-ion collisions can produce a medium made up of free quarks and gluons, known as the Quark-Gluon Plasma (QGP). Hard scattering processes occurring in these collisions produce high transverse momentum ( $p_T$ ) partons called jets which can be one or multiple. Due to the reason that jets carry color charges, they strongly couple with the QGP, thus lose energy as they propagate through the medium, resulting in the phenomenon of “jet quenching”, which is widely concerned and studied by various scholars [2 - 4].

A number of research groups have already published their measurement results related to energy of jets produced in these heavy-ion collisions at various collision energy levels. Some experimental data were compared to theoretical results [5 , 6] while others remained unexplained by theory. In this article, a model based on Glauber Monte Carlo is built to simulate the measurement results of dijet events $p_T$ correlations in Ref. [1] by the ATLAS collaboration, successful theoretical simulation of which has not been published before. As the final result of this model, the energy loss formula reveals the specific way that energies of jets are modified by QGP. The completion of this model has filled the blank of theoretical prediction of jet energy distributions in this specific energy level of collision. Following ATLAS’s way of illustrating the data, the distribution of jet energy is plotted against a variable defined as $x_J=p_{T2}/p_{T1}$ where $p_{T1}$ and $p_{T2}$ are the transverse momenta of the two jets and $p_{T2}\le p_{T1}$.

In the Glauber Monte Carlo model section of this article, frequently used terms include impact parameter, b defined as the distance between the center of two colliding heavy-ions and number of nucleons that are hit by at least one nucleon in heavy-ion collisions, $N_{part}$. These nucleons are also called participants.

2. Glauber Monte Carlo Model

2.1. Probability Density as a Function of Impact Parameter

In the real circumstance of a collider, two beams of heavy nuclei move in opposite direction. They randomly distribute in a cross-section area (of mm in dimension) far greater than the cross-section area of one nucleon (of fm in dimension).

To determine the probability density of a collision taking part at impact parameter range $b\sim b+\text{d}b$, first switch to a reference frame where one nucleon in one of the two beams is located at the origin. The positive direction of z axis is selected as the same as the velocity of this nucleon. In this reference frame, incident nuclei are evenly distributed in the whole xy plane.

The number of incident nuclei that fall on the thin ring corresponding to impact parameter range $b\sim b+\text{d}b$ is

$\text{d}N=\sigma 2\pi b\text{d}b$

where $\sigma$ is the area density of incident nucleon. Then the probability density function

$f\left(b\right)=\frac{\text{d}N}{\text{d}b}=2\pi \sigma b\propto b$ (1)

2.2. Applying Monte Carlo Glauber Model

The specific procedure of generating a Monte Carlo Glauber Model is well explained in Ref. [7]. Fermi distribution is used in our model, with the expression of nuclear charge density

$\rho \left(r\right)\propto \frac{1+w{\left(r/R\right)}^2}{1+\exp \left(\frac{r-R}{a}\right)}$ (2)

In the computer program, parameter r should be randomly generated according to the probability density function

$f\left(r\right)\propto 4\pi r^2\frac{1+w{\left(r/R\right)}^2}{1+\exp \left(\frac{r-R}{a}\right)}$ (3)

where the specific ration index does not make a difference. For Pb + Pb collision, the Fermi parameters are listed in Table 1.

2.3. Results of Glauber Monte Carlo Model

Implementing the Monte Carlo Glauber, we can plot a 2-D graph showing the distribution of all nucleons of two Pb nuclei on the transverse plane, as shown in Figure 1.

After that, we can get the number of participants—impact parameter and number of collisions— number of participants relationship, as shown in Figure 2.

2.4. Relationship of Impact Parameter to Centrality Ranges

In the procedure of Monte Carlo Glauber model, impact parameter b, the distance between the two centers of Pb nucleon, is an independent variable. In the ATLAS measurement Ref. [1], however, the variable centrality is used to represent the degree of overlapping between two nucleons. To utilize the MC Glauber model for simulating the results of ATLAS’s experiments, quantitative relationship of impact parameter to centrality needs to be determined.

The definition of centrality cannot be written explicitly in formula. It is based on the division of the histogram of $N_{part}$ as shown in Figure 3.

The curve of $N_{part}$ histogram and x-axis form a certain shape which corresponds to a certain area. In this graph, we can find a series of critical values $x_0\mathrm{,}{x}_1\mathrm{,}{x}_2,\cdots ,x_9\mathrm{,}{x}_{10}$ such that

Table 1. Fermi parameters of ²⁰⁸Pb.

Figure 1. An event for Pb + Pb collision at b = 7.0 fm. Participants are indicated as solid circles, while spectators are dotted circles. The two big dotted circles are the outlines of two Pb nuclei.

Figure 2. Histogram of $N_{part}-b$ and $N_{coll}-N_{part}$ for 10k events. The red line indicates average value.

Figure 3. Histogram of $N_{part}$ for 50k events. Different colors indicate different centrality ranges.

$x_0={\left(N_{part}\right)}_{\max }=416$

${\displaystyle {\int }_{x_1}^{x_0}}f\left(x\right)\text{d}x=10\%{\displaystyle {\int }_0^{x_0}}f\left(x\right)\text{d}x$

${\displaystyle {\int }_{x_2}^{x_1}}f\left(x\right)\text{d}x=\mathrm{10\; <mi>\; \%\; </mi>}{\displaystyle {\int }_0^{x_0}}f\left(x\right)\text{d}x$

$\cdots$

${\displaystyle {\int }_{x_9}^{x_{10}}}f\left(x\right)\text{d}x=\mathrm{10\; <mi>\; \%\; </mi>}{\displaystyle {\int }_0^{x_0}}f\left(x\right)\text{d}x$

$x_{10}={\left(N_{part}\right)}_{\min }=2$

Because $N_{part}$ is a discrete variable. The integrals above degenerate into summation. The specific values of these critical points are listed in Table 2.

3. Methodology

3.1. Calculation of Path Length

To simulate the energy loss of jets produced in QGP, the path length of each branch of the jet first needs to be calculated. Note that the QGP forms as a result of a number of nucleons colliding together, the shape of QGP is likely to be irregular and its boundary does not have a clear definition. It is thus more practical to define the path length based on position of participants, the origin of the jet and the equation of the jet line.

However, there still exists various ways to define jet length. The definition we adopt is the number of nucleon that each branch of the jet passes through. Under such definition, the path length is a discrete function. Another possible definition of the path length is the geometric distance from the origin of the jet and the point where the jet leaves the QGP. This definition was finally abandoned because of the following reasons. For one thing, as mentioned above, the boundary of the QGP is not clear. More importantly, it

Table 2. Critical $N_{part}$ values that divide the centrality ranges.

has already been proved that the density in each part of QGP is not homogeneous [8]. The density is higher where there are originally more participants and vice versa. So, the energy loss of the jet per unit length is not constant.

The specific criterion to judge whether a participant contribute to the jet and if it does, which branch of the jet it contribute to can be explained by Figure 4.

In this graph, all participants in the orange shadow area make one-unit contribution to upper branch of the jet and all participants in the blue shadow area make one-unit contribution to lower branch of the jet. The width of the shadow area is D.

In each centrality range, the histogram of path length of each branch of the jet, noted as $L_1$ and $L_2$, is shown in Figure 5.

3.2. Possible Forms of Energy Loss Formula

As the final step of our computer model, the form of the energy loss formula needs to be determined. As a jet propagates through QGP, the amount of energy it loses may be correlated to the path length and the initial energy of the jet. It is worth mentioning that the jets studied in our investigation have momenta of 100 GeV order of magnitude and thus are all ultra-relativistic. Their energy is approximately proportional to their momenta:

$E=pc.$

Thus energy loss is directly proportional to momentum loss.

Before the energy loss formula is finally determined, we can only come up with possible forms of it.

$\Delta p_T=Cf\left(L\right)g\left(p_T\right)R$ (4)

where C is a constant parameter to adjust the overall magnitude of energy loss, L is path length and R is a randomly generated parameter to simulate fluctuations. Each term of the formula may have the following forms:

$f\left(L\right)=L^{\alpha },\alpha =1,2,3,\cdots$ (5)

$g\left(p_T\right)=\{\begin{array}{l}{p}_T\\ M\left[\ln \left(\frac{p_T}{p_0}\right)\right]\\ \cdots \end{array}$ (6)

Figure 4. The graph helping to interpret the criterion to calculate the path length. The red point is the origin of the jet and the purple line indicates the rough boundary of the QGP.

Figure 5. Thermal diagrams of $L_1$ and $L_2$ in centrality range. (a) 0 - 10%; (b) 10% - 20%; (c) 20% - 30%; (d) 30% - 40%; (e) 40% - 60% and (f) 60% - 80% for 100k events.

$R\{\begin{array}{l}=\text{const}\\ \text{follows}\text{uniform}\text{distribution}\\ \text{follows}\text{Gaussian}\text{distribution}\\ \cdots \end{array}$ (7)

where the sign $M\left[f\left(x\right)\right]$ is defined as $M\left[f\left(x\right)\right]=f\left(x\right)\eta \left(x\right)=f\left(x\right)$ if $f\left(x\right)\ge 0$, $M\left[f\left(x\right)\right]=0$ if

$f\left(x\right)<0$. Note that if $g\left(p_T\right)=\ln \left(p_T/p_0\right)$, $g\left(p_T\right)<0$ when $p_T<p_0$, this sign is introduced to prevent

“negative” energy loss. A jet can never gain energy from QGP as it propagates through the medium. As shown in Ref. [9], $f\left(L\right)$ is expected to have the form $f\left(L\right)=L$ and $g\left(p_T\right)$ is expected to have the

form $g\left(p_T\right)=M\left[\ln \left(p_T/p_0\right)\right]$. On the other hand, the form $f\left(L\right)=L^2$ is derived in Ref. [10]. In our

investigation, $f\left(L\right)=L^{\alpha }$ is adopted, leaving the power of path length as a free variable. This reduces the total number of parameters to 3, C, $\alpha$ and $p_0$. The minima of this 3-dimensional parameter space is relatively easy to look for, and we search for it by traversing the parameter space, evenly selecting points first in larger intervals and then in smaller intervals.

3.3. Implementing the Energy Loss Formula to pp Data

As input data to our model, original data of jet events in pp collisions at $\sqrt{s_{NN}}=2.76\text{TeV}$ are obtained from CMS collaboration’s experimental results. Jets produced in pp collisions are considered to be the initial, unquenched state of jets produced in Pb + Pb collisions. These data are first selected based on the same selection rule of ATLAS’s measurement [1]: psudorapidity $\left|\eta \right|<2.1$, and difference in azimuthal angle $\Delta \varphi >7\pi /8$, where $\Delta \varphi =\left|{\varphi }_1-{\varphi }_2\right|$.

A pair of path length $L_1$ and $L_2$ corresponding to a certain centrality range is randomly generated by the Glauber Monte Carlo model for each jet-jet event. Given the energy loss formula, the amount of energy loss and the final energy of each jet can be calculated. Then the jets are selected based on final $p_T$ ranges and classified into different $x_J=p_{T2}/p_{T1}$ where $p_{T2}<p_{T1}$ ranges. The total number of events in each $x_J$ range is counted and normalized by dividing the total number of events which leads to the eventual simulated result.

3.4. Effects of Some Parameters on the Shape of the Curve

Even though the energy loss formula consists of several independent parameters, and they act on energy of unquenched jets as a whole. It is still discovered that different parameters have different modulation effect on the shape of simulated curve.

Firstly, it is found that the greater the variance of the random term R, the less obvious the peak at centrality range 0 - 10% and $p_T$ range 100 - 126 GeV is. Let R follow Gaussian distribution with different variances and $R=1$ as average value, the result of the curve is shown in Figure 6.

Secondly, we found that $g\left(p_T\right)=\ln \left(\frac{p_T}{p_0}\right)$ can usually make the predicted formula fit to Pb + Pb

curves in both different centrality ranges and different $p_{T1}$ ranges, even though $g\left(p_T\right)=p_T$ can

sometimes outperforms $g\left(p_T\right)=\ln \left(\frac{p_T}{p_0}\right)$ in a fixed $p_{T1}$ range.

Thirdly, it is found that for the form $g\left(p_T\right)=\ln \left(\frac{p_T}{p_0}\right)$, the position of the curve is mostly dominated

by the parameter $p_0$ when C is fixed, as shown in Figure 7.

3.5. Adjustment of Parameters

In the procedure of adjusting the energy loss formula, $f\left(L\right)$ is usually decided at first. After that, two parameters C and $p_0$ are adjusted to match the data in Ref. [1]. To describe the degree of resemblance of our predicted Pb + Pb data to that in ATLAS’s measurement, peaks in each plots of our predicted data and that in ATLAS’s measurement were found and compared. Also, we define a new variable “deviation” as:

$D={\displaystyle \underset{i=1}{\overset{10}{\sum }}}\left|y_i^{\text{Measure}}-y_i^{\text{Simulate}}\right|$

Figure 6. Simulated Pb + Pb curve given by energy loss formulas. $\Delta p_T=20{\left(L/40\right)}^{5.35}ln\left(\frac{p_T}{50}\right)R$ under centrality 0 - 10%. (a) $p_R\left(x\right)=\frac{1}{\sqrt{2\pi }{\sigma }^2}{e}^{-\frac{{\left(x-20\right)}^2}{2{\sigma }^2}}$, standard deviation $\sigma =0.1$. (b) $p_R\left(x\right)=\frac{1}{\sqrt{2\pi }{\sigma }^2}{e}^{-\frac{{\left(x-20\right)}^2}{2{\sigma }^2}}$, $\sigma =5.0$. (c) $p_R\left(x\right)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\left(x-20\right)}^2}{2{\sigma }^2}}$, σ = 10.0. from left to right respectively.

Figure 7. Simulated Pb + Pb curve under energy loss formulas of different p₀. $p_0=20,40,50,60,70,90$ GeV respectively. As $p_0$ increases, the peak shifts from small $x_J$ to large $x_J$.

which is the sum of the distances between simulated data points and experimental data points in 10 $x_J$ ranges. Also, the average deviation, $\bar{D}$, is defined as:

$\bar{D}$ = average value of D in six centrality ranges

The simulated data is considered to better resemble ATLAS’s experimental results when average deviation $\bar{D}$ is smaller.

In our investigation, there are basically two ways to select the parameters. One is to minimize average deviation of centrality ranges 0 - 10% and 10% - 20% at $p_T$ range 100 - 126 GeV, recorded as ${\bar{D}}_2$. The conspicuous peak of ATLAS’s measurement is located in this region.

The other is to minimize average deviation of all 9 different combinations of centrality ranges and $p_T$ ranges, recorded as ${\bar{D}}_9$.

The reason why these two ways of parameter selection were developed is that when ${\bar{D}}_9$ reaches its minumum, tendency of the simulated curve of centrality ranges 0 - 10% and 10% - 20% deviates significantly from the measured result, failing to match the most significant feature of ATLAS’s measurement. As shown in the Result section, these two ways of parameter selection will lead to two different forms of energy loss formula with better fitness at different centrality and $p_T$ ranges.

4. Results

By the end of our investigation, the two ways of parameter selection did not give similar forms energy loss formula, leading to two two different forms of energy loss formula with better fitness at different centrality and $p_T$ ranges instead. And the formula that has the best fitness under two criteria, minimum ${\bar{D}}_2$ and ${\bar{D}}_9$ are shown separately.

4.1. Optimization Result Focusing on Small Centrality Ranges

For the first way of parameter selection, the relationship of minimum ${\bar{D}}_2$ to $\alpha$ is shown and the dependence of ${\bar{D}}_2$ on C and $p_0$ at $\alpha =5.35$ is shown in Figure 8.

The energy loss formula with best fitness is:

(a) (b)

Figure 8. (a) The relationship of minimum ${\bar{D}}_2$ to $\alpha$, fitted by cubic spline function and (b) thermal diagram showing dependence of ${\bar{D}}_2$ on C and $p_0$ at $\alpha =5.35$.

$\Delta p_T=205{\left(L/40\right)}^{5.35}M\left[\ln \left(p_T/62\right)\right]$ (8)

corresponding to deviation ${\bar{D}}_2=1.87$ and ${\bar{D}}_9=2.45$. Predicted data is shown in Figure 9 and Figure 10.

Figure 9. The frequency distributions for jet pairs of $100<p_{T1}<126GeV$ for different centrality ranges under energy loss formula $\Delta p_T=205{\left(L/40\right)}^{5.35}M\left[ln\left(p_T/62\right)\right]$.

Figure 10. The frequency distributions for jet pairs of centrality 0 - 10% for different $p_T$ ranges under energy loss formula $\Delta p_T=205{\left(L/40\right)}^{5.35}M\left[ln\left(p_T/62\right)\right]$.

4.2. Optimization Result Considering All Centrality and p_T Ranges

For the second way of parameter selection, the relationship of minimum ${\bar{D}}_9$ to $\alpha$ and the dependence of ${\bar{D}}_9$ on C and $p_0$ at $\alpha =0.8$ is shown in Figure 11.

The energy loss formula with best fitness is:

$\Delta p_T=11{\left(L/40\right)}^{0.8}M\left[\ln \left(p_T/1.3\right)\right]$ (9)

corresponding to deviation ${\bar{D}}_9=1.84$ and ${\bar{D}}_2=3.56$. Predicted data is shown in Figure 12 and Figure 13.

4.3. Discussion of the Results

For the first way of optimization, the simulated curves deviate from experimental data at centrality ranges 20% - 30% and 30% - 40% most significantly and for the second way of optimization, simulated curves fail to resemble the conspicuous peak at centrality range 0 - 10% and when $100<p_{T1}<126\text{GeV}$.

It is not satisfying enough that an energy loss formula that can simulate all centrality and $p_T$ ranges very well has not been found. So possibility remains that 1) the procedure of jet quenching in QGP cannot be described by a simple form of energy loss formula in Equation (4), and 2) the unfolding procedure done by the ATLAS collaboration cannot fully represent physical reality as the peculiar peak at low centrality ranges was not reported in measurements done by other experimental groups [11].

5. Conclusion

In conclusion, our investigation has simulated the procedure of jet quenching in Quark-Gluon Plasma using Glauber Monte Carlo Model, using energy data of jets produced in pp collision obtained from the CMS collaboration. The research has characterised the energy loss of jets in QGP with path length of jets, L, and their initial transverse momentum, $p_T$. Simulated data were compared with ATLAS’s measurement in Ref. [1] to determine best combinations of parameters in our model. Two different forms of energy loss formula were found to adapt for certain centrality ranges and $p_T$ ranges well. Our result

(a) (b)

Figure 11. (a) The relationship of minimum ${\bar{D}}_9$ to $\alpha$, fitted by cubic spline function and (b) thermal diagram showing dependence of ${\bar{D}}_9$ on C and $p_0$ at $\alpha =0.8$.

Figure 12. The frequency distributions for jet pairs of $100<p_{T1}<126GeV$ for different centrality ranges under energy loss formula $\Delta p_T=11{\left(L/40\right)}^{0.8}M\left[ln\left(p_T/1.3\right)\right]$.

Figure 13. The frequency distributions for jet pairs of centrality 0 - 10% for different $p_T$ ranges under energy loss formula $\Delta p_T=11{\left(L/40\right)}^{0.8}M\left[ln\left(p_T/1.3\right)\right]$.

gives the most possible form of interaction between QGP and hadron jets. With our findings, the method of using jets to probe the features of QGP can be more effective as how they interact with each other is determined. We expect our findings can find their application in future studies of the features of Quark- Gluon Plasma.

Acknowledgements

We gratefully acknowledge Professor Gunther Roland for his instruction on this paper, and Yan Lu for his early stage contribution to this investigation.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Aaboud, M., Kupco, A., Davison, P., Webb, S., Sekula, S., Huston, J., et al. (2017) Measurement of Jet Correlations in Pb + Pb and pp Collisions at √SNN = 2.76 TeV with the ATLAS Detector. Physics Letters B, 774, 379-402. https://doi.org/10.1016/j.physletb.2017.09.078
[2]	Aaboud, M., et al. (2018) Observation of Centrality-Dependent Acoplanarity for Muon Pairs Produced via Two-Photon Scattering in Pb + Pb Collisions at √SNN = 5.02 TeV with the ATLAS Detector. Physical Review Letters, 121, Article ID: 212301. https://doi.org/10.1103/PhysRevLett.121.212301
[3]	Brewer, J., Milhano, J.G. and Thaler, J. (2019) Sorting out Quenched Jets. Physical Review Letters, 122, Article ID: 222301. https://doi.org/10.1103/PhysRevLett.122.222301
[4]	Brewer, J., Rajagopal, K., Sadofyev, A. and Van Der Schee, W. (2018) Evolution of the Mean Jet Shape and Dijet Asymmetry Distribution of an Ensemble of Holographic Jets in Strongly Coupled Plasma. Journal of High Energy Physics, 2, 15. https://doi.org/10.1007/JHEP02(2018)015
[5]	Aaboud, M., et al. (2019) Measurement of the Nuclear Modification Factor for Inclusive Jets in Pb + Pb Collisions at √SNN = 5.02 TeV with the ATLAS Detector. Physics Letters B, 790, 108-128. https://doi.org/10.1016/j.physletb.2018.10.076
[6]	Sirunyan, A.M., et al. (2018) Study of Jet Quenching with Isolated-Photon + Jet Correlations in PbPb and pp Collisions at √SNN = 5.02 TeV. Physics Letters B, 785, 14-39. https://doi.org/10.1016/j.physletb.2018.07.061
[7]	Alver, B., Baker, M., Loizides, C. and Steinberg, P. (2008) The PHOBOS Glauber Monte Carlo.
[8]	Alver, B. and Roland, G. (2010) Collision Geometry Fluctuations and Triangular Flow in Heavy-Ion Collisions. Physical Review C, 81, Article ID: 054905. ([Erratum: Physical Review C, 82, Article ID: 039903 (2010)]) https://doi.org/10.1103/PhysRevC.82.039903
[9]	d’Enterria, D. (2010) Jet Quenching. Landolt-Bornstein, 23, 471. https://doi.org/10.1007/978-3-642-01539-7_16
[10]	Djordjevic, M. (2009) Theoretical Formalism of Radiative Jet Energy Loss in a Finite Size Dynamical QCD Medium. Physical Review C, 80, Article ID: 064909. https://doi.org/10.1103/PhysRevC.80.064909
[11]	Chatrchyan, S., et al. (2012) Jet Momentum Dependence of Jet Quenching in PbPb Collisions at √SNN = 2.76 TeV. Physics Letters B, 712, 176-197. https://doi.org/10.1016/j.physletb.2012.04.058