Heavy and light flavor jet quenching in different collision systems at the LHC energies

Jet quenching and anisotropic flows are the two most important signatures of the quark-gluon plasma (QGP) created in high-energy nucleus-nucleus collisions performed at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC). At low transverse momentum ($p_{T}$) region, the produced final-state hadrons exhibit large anisotropies in the azimuthal angle distribution. Such anisotropy can be successfully explained by relativistic hydrodynamics simulation Romatschke:2017ejr ; Rischke:1995ir ; Heinz:2013th ; Gale:2013da ; Huovinen:2013wma , which suggests that the produced hot and dense nuclear matter is a strongly-interacting QGP. The strong interaction among QGP constituents can convert initial-state geometric anisotropy into final-state anisotropy of hadron momentum distributions Alver:2010gr ; Petersen:2010cw ; Qin:2010pf ; Staig:2010pn ; Teaney:2010vd ; Schenke:2010rr ; Ma:2010dv ; Qiu:2011iv . At high $p_{T}$ region, the final-state hadron yields are strongly suppressed compared to the binary-scaled nucleon-nucleon collisions. Such phenomenon can be well explained by the energy loss of hard jet partons during their propagation through the soft QGP medium, and is usually called jet quenching Wang:1991xy ; Qin:2015srf ; Blaizot:2015lma ; Majumder:2010qh ; Gyulassy:2003mc ; Cao:2020wlm . The interaction between jets and medium may also lead to the modification of jet-related correlations Qin:2009bk ; Chen:2016vem ; Chen:2016cof ; Chen:2017zte ; Luo:2018pto ; Zhang:2018urd ; Kang:2018wrs , the suppression of full jet yields Qin:2010mn ; Young:2011qx ; Dai:2012am ; Wang:2013cia ; Blaizot:2013hx ; Mehtar-Tani:2014yea ; Cao:2017qpx ; Kang:2017frl ; He:2018xjv ; Chang:2016gjp ; Casalderrey-Solana:2016jvj and the change of jet structure and substructure Tachibana:2017syd ; KunnawalkamElayavalli:2017hxo ; Brewer:2017fqy ; Chien:2016led ; Milhano:2017nzm ; Chang:2019sae . Currently, one important task of heavy-ion physics is to quantitatively extract various transport coefficients of the QGP and jet transport parameters using low and high $p_{T}$ data Bernhard:2019bmu ; Everett:2020xug ; Burke:2013yra ; Cao:2021keo ; Cao:2018ews ; Rapp:2018qla , and ultimately to probe the structure and transport properties of QGP at various energy and length scales.

In recent years, experiments have observed large anisotropic flows in proton-lead collisions at the LHC and deuteron-gold collisions at RHIC Abelev:2012ola ; Aad:2012gla ; Chatrchyan:2013nka ; Adare:2013piz . The successfulness of relativistic hydrodynamics in explaining the collectivity in such small systems implies the possible formation of mini QGP Bozek:2013uha ; Bzdak:2013zma ; Qin:2013bha ; Werner:2013ipa ; Bozek:2013ska ; Nagle:2013lja ; Schenke:2014zha ; Zhao:2019ehg ; Zhao:2017rgg . However, no jet quenching has been confirmed in these small systems so far, e.g., $R_{pA}\approx 1$. On the other hand, the color glass condensate framework shows that the initial-state effect, i.e., the interaction between partons originated from the nucleon projectile and dense gluons inside the target nucleus before the onset of hydrodynamics evolution, can also generate significant amount of collectivity Dusling:2012iga ; Dusling:2017dqg ; Mace:2018vwq ; Davy:2018hsl ; Zhang:2019dth . Nowadays, how to disentangle the contributions from initial and final state effects is still in hot debate.

In order to understand the above intriguing results, the system size scan has been proposed to bridge the gap between large and small collision systems Citron:2018lsq ; Huang:2019tgz . By scanning different sizes of collision systems, we may build a smooth transition from large to small collision systems, and help to disentangle the initial and final state contributions. Along this line, the system size scan predictions for Pb+Pb, Xe+Xe, Ar+Ar, O+O and other collision systems have been performed from relativistic hydrodynamics Sievert:2019zjr ; Lim:2018huo . As for hard probes, there have been studies on the nuclear modification factor $R_{AA}$ for Xe+Xe collisions Zigic:2018ovr ; Shi:2019nyp and O+O collisions Huss:2020dwe ; Huss:2020whe and the system scan study for heavy quarks and heavy flavor mesons Katz:2019qwv . However, a systematic study of jet quenching in different collision systems for both light and heavy flavor partons and hadrons is still lacking. This is the main purpose of our work.

In this work, we perform a systematic study on both heavy and light flavor jet quenching in different collision systems at the LHC energies. The productions of jet partons and hadrons are calculated within a next-to-leading-order perturbative QCD framework, the evolution of heavy and light flavor jet partons inside the QGP is simulated via a linear Boltzmann transport model, and the space-time profile of the QGP fireball is obtained via a (3+1)-dimensional viscous hydrodynamics simulation. Using our state-of-the-art jet quenching model Xing:2019xae , we obtain a nice description of $R_{\rm AA}$ for charged hadrons and $D$ mesons in central and mid-central Pb+Pb and Xe+Xe collisions measured by CMS collaboration. Using the same setup, we further compute the $p_{T}$ and centrality dependences of $R_{AA}$ for charged hadrons, $D$ and $B$ mesons in Pb+Pb, Xe+Xe, Ar+Ar and O+O collisions at the LHC energies. Our study shows a clear system size dependence for both light and heavy flavor hadron $R_{AA}$ across different collision systems. This suggests that $R_{pA}\approx 1$ in proton-nucleus collisions is mainly due to the small size of the collision system. Interestingly, we observe sizable jet quenching effect for both heavy and light flavor hadrons in central O+O collisions at the LHC energies, which can be tested in future experiments.

The paper is organized as follows. In Section II, we briefly present the next-to-leading-order framework for computing the momentum spectra of hard jet partons and high $p_{T}$ hadron. In Section III, we provide the linear Boltzmann transport model for describing the evolution of jet partons interacting with the soft QGP medium. In Section IV, we discuss the setups of the initial conditions and hydrodynamics simulation for different collision systems. In Section V, we provide the numerical results on $R_{AA}$ as a function of $p_{T}$ and centrality for charged hadrons, $D$ mesons, $B$ mesons and heavy flavor decayed electrons in Pb+Pb, Xe+Xe, Ar+Ar and O+O collisions at the LHC energies. Section VI contains our summary.

In this work, we use the next-to-leading-order (NLO) perturbative QCD framework developed in Refs. Jager:2002xm ; Aversa:1988vb to calculate the production of both light and heavy flavor partons and hadrons at high transverse momentum ($p_{T}$) region. In this framework, the differential cross section for high $p_{T}$ hadrons can be written as follows:

	$\displaystyle E_{h}\frac{d\sigma_{pp\to h}}{d^{3}p_{h}}$	$\displaystyle=\frac{1}{\pi S}\sum_{abc}\int\frac{dz_{c}}{z_{c}^{2}}\int\frac{dv}{v(1-v)}\int\frac{dw}{w}$
		$\displaystyle\times f_{a}(x_{a},\mu)f_{b}(x_{b},\mu)\frac{d\sigma_{ab\to c}}{dvdw}D_{h/c}(z_{c},\mu).$		(1)

In the above equation, the sum $\sum_{abc}$ is over all possible partonic species. The integral ranges are set as: $z_{c}\in[1-V+VW,1]$, $v\in[VW/z_{c},1-(1-V)/z_{c}]$, $w\in[VW/(vz_{c}),1]$, where the hadronic variables $S$, $T$, $U$, $V$, $W$ are defined as:

	$\displaystyle V=1+\frac{T}{S},\,W=-\frac{U}{S+T},$
	$\displaystyle S=(p_{A}+p_{B})^{2},\,T=(p_{A}-p_{h})^{2},\,U=(p_{B}-p_{h})^{2}.$		(2)

The corresponding partonic variables are defined as:

	$\displaystyle v=1+\frac{t}{s},\,w=-\frac{u}{s+t},$
	$\displaystyle s=(p_{a}+p_{b})^{2},\,t=(p_{a}-p_{c})^{2},\,u=(p_{b}-p_{c})^{2}.$		(3)

$f_{a}(x_{a})$ and $f_{b}(x_{b})$ are parton distribution functions (PDFs), and $D_{h/c}$ is fragmentation function (FF). Neglecting all masses, one has the following relations:

	$\displaystyle x_{a}=\frac{VW}{vwz_{c}},\,x_{b}=\frac{1-V}{(1-v)z_{c}},$
	$\displaystyle s=x_{a}x_{b}S,\,t=\frac{x_{a}}{z_{c}}T,\,u=\frac{x_{b}}{z_{c}}U.$		(4)

At one-loop level, the partonic cross section takes the following form:

$\displaystyle\frac{d\sigma_{ab\to c}}{dvdw}=\frac{d\sigma_{ab\to c}^{(0)}}{dv}\delta(1-w)+\frac{\alpha_{s}(\mu)}{\pi}\frac{d\sigma_{ab\to c}^{(1)}}{dvdw}.$

(5)

Since the above NLO perturbative QCD framework neglects all masses, it is expected to be applicable in high $p_{T}$ regime, $p_{T}\gg m_{c},m_{b}$. In this work, we focus on heavy and light flavor hadron productions at high $p_{T}$ ($p_{T}>$ 6-8 GeV), where the fragmentation mechanism is also the dominant contribution for hadron production. In our calculation, the parton distribution functions are taken from CTEQ parameterizations Pumplin:2002vw . The cold nuclear matter (CNM) effects on PDFs in nucleus are taken from EPS09 parameterizations Eskola:2009uj . The fragmentation functions are taken from Refs. Kretzer:2000yf ; Kneesch:2007ey ; Kniehl:2008zza for charged hadrons, $D$ mesons, and $B$ mesons, respectively.

In Fig. 1, we show the transverse momentum spectra for charged hadrons, $D$ mesons and $B$ mesons in p+p collisions. The upper two plots are charged hadron spectra in p+p collisions at $\sqrt{s}=5.02$ TeV and $\sqrt{s}=5.44$ TeV, and the lower two plots are $D$ and $B$ meson spectra in p+p collisions at $\sqrt{s}=5.02$ TeV. In the NLO perturbative QCD calculation, the scales in the above equations are all taken as the transverse momenta of hadrons $\mu=p_{T}$. One can see that the NLO perturbative QCD framework can provide a nice description of both heavy and light flavor hadron production at high $p_{T}$ region in nucleon-nucleon collisions, as measured by CMS collaboration Khachatryan:2016odn ; Sirunyan:2017xss ; Sirunyan:2018eqi ; Sirunyan:2017oug . This serves as the baseline to study the hot QGP medium modification effect in high-energy nucleus-nucleus collisions.

We use linear Boltzmann transport (LBT) model to simulate the elastic and inelastic interactions between hard partons and the soft QGP Wang:2013cia ; He:2015pra ; Cao:2016gvr ; Cao:2017hhk ; Xing:2019xae . In the LBT model, one numerically uses the Monte-Carlo method to solve the following Boltzmann equation for the propagating hard parton “1”:

$\displaystyle p_{1}\cdot\partial f_{1}(x_{1},p_{1})=E_{1}(C_{\mathrm{el}}+C_{\mathrm{inel}}),$

(6)

where $C_{\mathrm{el}}$ and $C_{\mathrm{inel}}$ represent the collision integrals for elastic and inelastic interactions between the hard parton and the medium constituents.

To simulate the elastic process ($1+2\rightarrow 3+4$), the key quantity is to calculate the elastic scattering rate:

	$\displaystyle\Gamma_{\rm el}(\vec{p}_{1})$	$\displaystyle=\frac{\gamma_{2}}{2E_{1}}\int\frac{d^{3}p_{2}}{(2\pi)^{3}2E_{2}}\int\frac{d^{3}p_{3}}{(2\pi)^{3}2E_{3}}\int\frac{d^{3}p_{4}}{(2\pi)^{3}2E_{4}}$
		$\displaystyle\times f_{2}(\vec{p}_{2})[1\pm f_{3}(\vec{p}_{3})][1\pm f_{4}(\vec{p}_{4})]S_{2}(s,t,u)$
		$\displaystyle\times(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})|M_{12\rightarrow 34}|^{2}.$		(7)

In the above equation, $\gamma_{2}$ is the degeneracy factor of parton “2”. The condition $S_{2}(s,t,u)=\theta(s\geq 2\mu_{D}^{2})\theta(t\leq-\mu_{D}^{2})\theta(u\leq-\mu_{D}^{2})$ is imposed to avoid possible divergence at small angle $u,t\rightarrow 0$ Auvinen:2009qm ; He:2015pra , with the Debye screening mass taken as $\mu_{D}^{2}=6\pi\alpha_{s}T_{m}^{2}$, with $T_{m}$ the temperature of the QGP medium. The $\pm$ sign represents the Pauli block and Bose enhancement effects for the thermal partons. The factor $(1-f_{3})$ is neglected in the calculation due to the much lower density of hard partons in the QGP. The matrix element $|M_{12\rightarrow 34}|^{2}$ is taken from the leading order perturbative QCD calculation Combridge:1978kx . Since we only consider the leading-order $2\to 2$ scattering processes, the distribution for the exchanged transverse momentum ($q_{\perp}$) between hard partons and QGP typically has a hard tail. To account for the possible contributions from multiple soft scatterings, we impose an effective cutoff $q_{\perp}<10T_{m}$ for exchanged transverse momentum. Such cutoff will reduce the parton energy loss, which can be compensated by shifting effective strong coupling $\alpha_{\rm s}$ to a larger value. In this work, we use $\alpha_{\mathrm{s}}=0.25$ for the interaction vertex connecting to the thermal partons inside the medium. For the vertices connecting to jet partons, we set $\alpha_{s}=4\pi/\left[9\ln(2E_{1}T_{m})/\Lambda^{2})\right]$ with $\Lambda=0.2~{}\text{GeV}$.

To simplify the above elastic scattering rate, we may choose the jet parton momentum $\vec{p}_{1}$ along the $+z$ direction and the medium parton $\vec{p}_{2}$ in the $xz$ plane. Then the elastic scattering rate can be obtained as follows:

	$\displaystyle\Gamma_{\rm el}(\vec{p}_{1})$	$\displaystyle=\frac{\gamma_{2}}{16E_{1}(2\pi)^{4}}\int dE_{2}d\theta_{2}d\theta_{4}d\phi_{4}$
		$\displaystyle\times f_{2}(E_{2},T_{m})(1\pm f_{4}(E_{4},T_{m}))S_{2}(s,t,u)\big{|}{M_{\mathrm{12\rightarrow 34}}}\big{|}^{2}$
		$\displaystyle\times\frac{p_{2}p_{4}\sin\theta_{2}\sin\theta_{4}}{E_{1}+E_{2}-p_{1}\cos\theta_{4}-p_{2}\cos\theta_{24}},$		(8)

where $\theta_{2}$ is the polar angle of the momentum $\vec{p}_{2}$, $(\theta_{4},\phi_{4})$ are the polar and azimuthal angles of the momentum $\vec{p}_{4}$, and $\theta_{24}$ is the angle between $\vec{p}_{2}$ and $\vec{p}_{4}$:

$\displaystyle\cos\theta_{24}=\sin\theta_{2}\sin\theta_{4}\cos\phi_{4}+\cos\theta_{2}\cos\theta_{4}.$

(9)

From energy-momentum conservation, one can find,

$\displaystyle E_{4}=\frac{p_{1}\cdot p_{2}+m_{4}^{2}}{E_{1}+E_{2}-p_{1}\cos\theta_{4}-p_{2}\cos\theta_{24}}.$

(10)

Using the above formula, one can calculate the elastic scattering rate $\Gamma_{\rm el}^{(i)}$ for each channel and the total scattering rate $\Gamma_{\rm el}=\sum_{i}\Gamma_{\rm el}^{(i)}$. If the number of elastic scatterings in a given time step $\Delta t$ follows the Poisson distribution, then the probability for having elastic scatterings in $\Delta t$ is given by:

$\displaystyle P_{\rm el}=1-e^{-\Gamma_{\rm el}\Delta t}.$

(11)

For inelastic radiative process, we use the higher-twist formalism for the medium-induced gluon radiation spectrum Wang:2001ifa ; Zhang:2003wk ,

	$\displaystyle\frac{dN_{g}}{dxdk_{\perp}^{2}dt}=$	$\displaystyle\frac{2\alpha_{s}C_{A}P(x)k_{\perp}^{4}}{\pi(k_{\perp}^{2}+x^{2}M^{2})^{4}}\hat{q}(E_{1},T_{m})$
		$\displaystyle\times\sin^{2}\bigg{(}\frac{t-t_{i}}{2\tau_{f}}\bigg{)},$		(12)

where $x$ and $k_{\perp}$ are the fraction energy and transverse momentum of the emitted gluon, $M$ is the mass of the parent parton (here we use $M_{c}=1.27$ GeV/c and $M_{b}=4.19$ GeV/c), $\alpha_{s}$ is the strong coupling for the splitting vertex, $P(x)$ is the vacuum splitting function, $\tau_{f}=2Ex(1-x)/(k_{\perp}^{2}+x^{2}M^{2})$ is the formation time for gluon radiation, $(t-t_{i})$ denotes accumulative time for a gluon radiation event, i.e., $t_{i}$ is reset to zero after each radiation, $\hat{q}(E_{1},T_{m})$ is the jet transport parameter characterizing the transverse momentum broadening rate from elastic collisions and can be calculated as follows:

	$\displaystyle\hat{q}$	$\displaystyle(E_{1},T_{m})=\frac{\gamma_{2}}{16E_{1}(2\pi)^{4}}\int dE_{2}d\theta_{2}d\theta_{4}d\phi_{4}(\vec{p}_{2\perp}-\vec{p}_{4\perp})^{2}$
		$\displaystyle\times f_{2}(E_{2},T_{m})(1\pm f_{4}(E_{4},T_{m}))S_{2}(s,t,u)\big{|}{M_{\mathrm{12\rightarrow 34}}}\big{|}^{2}$
		$\displaystyle\times\frac{p_{2}p_{4}\sin\theta_{2}\sin\theta_{4}}{E_{1}+E_{2}-p_{1}\cos\theta_{4}-p_{2}\cos\theta_{24}}.$		(13)

The above formula naturally includes the dependences of $\hat{q}$ on jet energy and medium temperature across different colliding systems investigated in this work. Using the above gluon radiation formula, one may calculate the average gluon radiation rate in a give time step $\Delta t$ as follows:

$\displaystyle\Gamma_{\rm inel}=\int dxdk_{\perp}^{2}\frac{dN_{g}}{dxdk_{\perp}^{2}dt},$

(14)

In this work, we impose a lower cut-off $x_{\rm min}=\mu_{D}/E$ for the emitted gluon energy to avoid the possible divergence as $x\to 0$. If the number of radiated gluons follows a Poisson distribution, then the probability for having inelastic scatterings in a given time step $\Delta t$ is given by:

$\displaystyle P_{\mathrm{inel}}=1-e^{-\Gamma_{\rm inel}\Delta t}.$

(15)

To include both elastic and inelastic interactions, the total scattering rate is given by:

$\displaystyle\Gamma_{\rm tot}=\Gamma_{\rm el}+\Gamma_{\rm inel}.$

(16)

The total scattering probability is then obtained as:

$\displaystyle P_{\rm tot}=1-e^{-\Gamma_{\rm tot}\Delta t}=P_{\rm el}+P_{\rm inel}-P_{\rm el}P_{\rm inel}.$

(17)

One can split the total interaction probability into two parts: $P_{\rm el}(1-P_{\rm inel})$ for pure elastic scatterings and $P_{\rm inel}$ for inelastic scatterings.

In this work, we use the (3+1)-dimensional CLVisc Pang:2018zzo hydrodynamics model to simulate the space-time evolution of QGP created in different collision systems: Pb+Pb collisions at 5.02A TeV, Xe+Xe collisions at 5.44A TeV, Ar+Ar collisions at 5.85A TeV and O+O collisions at 6.5A TeV. The initial conditions for hydrodynamics simulation are obtained as follows. The entropy distribution $S(\tau_{0},x,y,\eta_{s})$ at the initial proper time $\tau_{0}$ is factorized into the entropy density $s(x,y)$ in the transverse plane and the longitudinal envelop function $H(\eta_{s})$,

$$S(\tau_{0},x,y,\eta_{s})=Ks(x,y)H(\eta_{s}),$$

(18)

where $K$ is an overall normalization factor which is adjusted to describe the charged hadron spectra in the most central collisions Adam:2016ddh ; Acharya:2018hhy . For Ar-Ar and O-O collisions at the LHC energies, where the experimental data on multiplicity are unavailable, we use the empirical power-law formula: ${2}/{\left<N_{\text{part}}\right>}\left<dN_{\text{ch}}/d\eta\right>|_{|\eta|<0.5}=0.7455(s_{\rm NN}/{\rm GeV}^{2})^{0.1538}$ Acharya:2018hhy , which is fitted to the multiplicity data of RHIC Au+Au and LHC Pb+Pb collisions at various colliding energies. The entropy density $s(x,y)$ in the transverse plane is generated from the TRENTo modelMoreland:2014oya , in which the positions of nucleons in nucleus are sampled according to the following Woods-Saxon distribution:

$$\rho(r,\theta)=\rho_{0}\frac{1+w\frac{r^{2}}{R^{2}(\theta)}}{1+\exp\left(\frac{r-R(\theta)}{d}\right)},$$

(19)

where $\rho_{0}$ denotes the nuclear density at the center, $d$ is a surface thickness parameter, the radius parameter

$\displaystyle R(\theta)=R_{0}[1+\beta_{2}Y_{20}(\theta)+\beta_{4}Y_{40}(\theta)],$

(20)

with $Y_{nl}(\theta)$ being the spherical harmonic functions, and $(w,\beta_{2},\beta_{4})$ controlling the deviations from spherical shape. Table.(1) shows the parameters used the Woods-Saxon distribution for different collision systems.

Using the Glauber model simulation, one may calculate the nuclear thickness function $T_{A}(x,y)$ and $T_{B}(x,y)$ by assigning each participant nucleon a Gaussian smearing function (with width 0.5 fm). Then the local entropy density $s(x,y)$ is constructed from the thickness function according to the following formula,

$$s(x,y)=\left(\frac{T_{A}^{p}+T_{B}^{p}}{2}\right)^{\frac{1}{p}}.$$

(21)

In this work, we choose the parameter $p=0$ which mimics the IP-Plasma or EKRT initial condition models. For the longitudinal direction, we use the following parameterized form for the envelop function $H(\eta_{s})$ Pang:2018zzo ,

$$H(\eta_{s})=\exp\left[-\frac{(|\eta_{s}|-\eta_{0})^{2}}{2\sigma^{2}_{\eta_{s}}}\theta(|\eta_{s}|-\eta_{0})\right].$$

(22)

The values of $\eta_{0}$, $\sigma_{\eta_{s}}$ and the overall normalization $K$ and the initial time $\tau_{0}$ for different collision systems are listed in Table 2.

Using the above setup, we may obtain the initial conditions for hydrodynamics simulation. In this work, for each centrality we average 5000 Trento events to get the smooth initial entropy distribution as the input for hydrodynamics simulation. Then we use CLVisc hydrodynamics model Pang:2018zzo to simulate the space-time evolution of the QGP. In the hydrodynamics model, the following equation of motion for the energy-momentum tensor $T^{\mu\nu}$ and the shear stress tensor $\pi^{\mu\nu}$ are solved in Milne coordinate using the Kurganov-Tadmor (KT) algorithm,

	$\displaystyle\partial_{\mu}T^{\mu\nu}=0,$		(23)
	$\displaystyle\pi^{\mu\nu}=\eta_{v}\sigma^{\mu\nu}-\tau_{\pi}\left[\Delta^{\mu\nu}_{\alpha\beta}u^{\lambda}\partial_{\lambda}\pi^{\alpha\beta}+\frac{4}{3}\pi^{\mu\nu}\theta\right],$		(24)

where $\sigma^{\mu\nu}$ is the symmetric shear tensor, $\theta$ is the expansion rate. In this study, the shear-viscosity-to-entropy-density-ratio is chosen to be $\eta_{v}/s=0.16$ and the relaxation coefficient is $\tau_{\pi}=\frac{3\eta_{v}}{sT_{m}}$, and the equation of state is taken from the s95-pce parameterization. After hydrodynamic evolution, the QGP medium will be converted to hadrons according to the Cooper-Frye formula where the switching temperature is taken as $T_{\rm sw}=137$ MeV. Note that the soft hadrons produced from the freezeout hypersurface of the fluid have negligible contribution to the nuclear modification of hard hadrons with $p_{T}>6-8$ GeV explored in this paper Zhao:2021vmu .

In Fig. 3, we show the nuclear modification factor $R_{AA}$ as a function of transverse momentum $p_{T}$ for charged hadrons in Pb+Pb collisions at 5.02A TeV, Xe+Xe collisions at 5.44A TeV, Ar+Ar collisions at 5.85A TeV and O+O collisions at 6.5A TeV at the LHC. Four panels show four different centralities: 0-10%, 10-30%, 30-50% and 50-80%. The bands show the uncertainties due to the cold nuclear matter effect: solid/dashed/dash-dotted/dash-dot-dotted curves represent our calculation results with CNM while dotted curves are corresponding results without CNM. Note that as there is no nuclear PDF for ⁴⁰Ar from EPS09 parameterizations, we take the nuclear PDF of ⁴⁰Ca as an approximation in the calculation. We have checked that for Pb+Pb collisions, the use of the impact parameter dependent nPDF – EPS09s Helenius:2012wd makes the CNM uncertainty band of $R_{AA}$ narrower in peripheral collisions, but has small effect in central collisions. One possible reason is that central collisions have much larger binary collision number than peripheral collisions, thus have much larger contribution when averaging over the centrality (impact parameter) for the CNM effect.

From the plots, one can see that our jet quenching model can provide a very good description of the transverse momentum and centrality dependences of $R_{AA}$ for charged hadrons in Pb+Pb and Xe+Xe collisions as measured by CMS Collaboration. Using the same setups, we predict charged hadron $R_{AA}$ for smaller collision systems, such as Ar+Ar and O+O collisions, as a function of $p_{T}$ and centrality. One can clearly see the system size dependence for $R_{AA}$, i.e., $R_{AA}$ is larger for smaller systems, due to smaller jet energy loss effect. The system size dependence can also been seen from the centrality dependence: $R_{AA}$ becomes closer to unity when moving from central to peripheral collisions.

It is very interesting to see that the jet quenching effect on $R_{AA}$ is still very significant in central O+O collisions: $R_{AA}$ at $p_{T}\sim$ 8-10 GeV can reach about 0.5. Such quenching effect should be large enough to be observed in the future O+O collisions at the LHC. The confirmation of our result in the future can serve as the bridge to search for the QGP signal in even smaller collision systems like proton-nucleus collisions. It is worth noting that the experimental bias in event selection may lead to sizable deviation from unity for $R_{AA}$ in peripheral collisions Acharya:2018njl ; Loizides:2017sqq .

In order to see the system size dependence more clearly, we plot in Fig. 3 the nuclear modification factor $R_{AA}$ as a function of $N_{\rm part}$ for four different collision systems (as shown by four different curves). The left figure shows the result for $p_{T}=(6.7,7.7)$ GeV and the right figure for $p_{T}=(26,30)$ GeV. Our model can describe the experimental data for Pb+Pb collisions and Xe+Xe collisions quite well. One now sees that $R_{AA}$ is a decreasing function of $N_{\rm part}$, which clearly demonstrates the system size effect on the jet energy loss and nuclear modification factor $R_{AA}$. As one goes from larger to smaller collision systems, $R_{AA}$ increases and tends to reach unity. One also observes that for the same values of $N_{\rm part}$, there exist some differences in $R_{AA}$ from four different collision systems due to different collision energies. We have checked that the effect of the different initial parton spectra is negligible (as the collision energies are not very different). This means that different initial geometries and hydrodynamic evolution profiles play more important roles for such difference.

In Fig. 5, we show $R_{AA}$ of $D$ mesons as a function of $p_{T}$ for Pb+Pb, Xe+Xe, Ar+Ar and O+O collisions at the LHC energies. The setup is the same as Fig. 3. Four different centralities (0-10%, 10-30%, 30-50% and 50-80%) are shown. One can see that our model can can also describe $D$ meson $R_{AA}$ as a function of $p_{T}$ and centrality for Pb+Pb collisions as measured by ALICE and CMS Collaborations. This means that the flavor, energy and system size dependences of jet quenching can be well understood in our jet quenching model. Using the same setups, the predictions for smaller collision systems (Ar+Ar and O+O collisions) are shown. We can see that the jet quenching effect on $D$ meson $R_{AA}$ is also quite significant in central O+O collisions.

In Fig. 5, we show $R_{AA}$ of $D$ mesons as a function of $N_{\rm part}$ in four collision systems for two $p_{T}$ bins: $p_{T}=(6.7,7.7)$ GeV in the left panel and $p_{T}=(26,30)$ GeV on the right. Similar to charged hadron $R_{AA}$, $D$ meson $R_{AA}$ is also a decreasing function of $N_{\rm part}$, which clearly shows that jet quenching effect becomes weaker for smaller system sizes.

Figure 7 shows $B$ meson $R_{AA}$ as a function of $p_{T}$ for four different collision centralities (0-10%, 10-30%, 30-50% and 50-80%) in Pb+Pb, Xe+Xe, Ar+Ar and O+O collisions at the LHC energies. One can see that due to larger masses of $b$ quarks, the values of $B$ meson $R_{AA}$ are larger than $D$ meson and charged hadron $R_{AA}$. We can see that even for $B$ mesons, jet quenching effect on $R_{AA}$ is still sizable in central O+O collisions. Figure 7 shows $B$ meson $R_{AA}$ as a function of $N_{\rm part}$ in four collision systems for $p_{T}=(6.7,7.7)$ GeV (Left) and $p_{T}=(26,30)$ GeV (Right). Similar system size ($N_{\rm part}$) dependence can be seen for $B$ meson $R_{AA}$. Our results on the flavor, mass, energy and system size dependences of jet quenching can be tested by future measurements.

Recently, experiments have also measured the nuclear modification factor $R_{AA}$ of heavy flavor decayed electrons. Using the same setups, we further calculate $R_{AA}$ for electrons decayed from $D$ and $B$ mesons for central 0-10%, 10-30%, 30-50% and 50-80% Pb+Pb, Xe+Xe, Ar+Ar and O+O collisions at the LHC energies, as shown in Fig. 8. One can see that our model can also describe the heavy flavor decayed electron $R_{AA}$ in Pb+Pb collisions at 5.02A TeV measured by ALICE Collaboration. In addition, we also observe similar system size dependence for heavy flavor decayed electron $R_{AA}$.

Using our state-of-the-art jet quenching model, we have performed a systematic study on heavy and light flavor jet quenching in different collision systems at the LHC energies. The momentum spectra of light and heavy flavor hadrons at high $p_{T}$ are calculated within a next-to-leading-order perturbative QCD framework, by taking into account both quark and gluon contributions to high $p_{T}$ light and heavy flavor hadron productions. The evolution and energy loss of light and heavy partons in the QGP is simulated via the linear Boltzmann transport (LBT) model. The space-time evolution of the QGP fireball is simulated using a (3+1)-dimensional viscous hydrodynamics. Using our jet quenching model, we have calculated the nuclear modification factor $R_{AA}$ for charged hadrons, $D$ mesons, $B$ mesons and heavy-flavor decayed electrons as a function of $p_{T}$ and centrality ($N_{\rm part}$) for central 0-10%, 10-30%, 30-50% and 50-80% Pb+Pb, Xe+Xe, Ar+Ar and O+O collisions at the LHC energies. Our model can give a good description of charged hadron, $D$ meson and heavy flavor decayed electron $R_{\rm AA}$’s in central and mid-central Pb+Pb (and Xe+Xe) collisions measured by CMS and ALICE Collaborations. We have also provided the detailed predictions for the transverse momentum and centrality dependences of $R_{AA}$’s for charged hadrons, $D$ meson, $B$ mesons, and heavy flavor decayed electrons in Pb+Pb, Xe+Xe, Ar+Ar and O+O collisions at the LHC energies. Our results show that both light and heavy flavor $R_{AA}$’s have a strong system size dependence. As the system size becomes smaller, jet quenching effect on $R_{AA}$ diminishes: $R_{AA}$ tends to approach unity as the system size shrinks. Our result indicates that $R_{pA}\approx 1$ in proton-nucleus collisions is mainly due to the small size of the nuclear medium produced in the collisions. Very interestingly, our calculation also show significant jet quenching effect on $R_{AA}$’s for charged hadrons and heavy flavor hadrons and their decayed electrons in central O+O collisions at the LHC energies, which can be tested by future experiments. In summary, our study provides an important baseline for studying the flavor, mass, energy and system size dependences of jet quenching. Further studies along this line should help to identify the unique signatures of mini QGP droplet and search for the disappearance of QGP by scanning the system size and collision energy of relativistic nuclear collisions.

This work is supported in part by Natural Science Foundation of China (NSFC) under Grants No. 11775095, No. 11890710, No. 11890711 and No. 11935007. H. X. is supported by the Guangdong Major Project of Basic and Applied Basic Research No. 2020B0301030008, Key Project of Science and Technology of Guangzhou (Grant No. 2019050001), the National Natural Science Foundation of China under Grant No. 12022512, No. 12035007.