††thanks: present address is 6

# Linear Boltzmann transport for jet propagation in the quark-gluon plasma:
Inelastic processes and jet modification

Tan Luo<EMAIL_ADDRESS>Key Laboratory of Quark and Lepton Physics (MOE) and
Institute of Particle Physics, Central China Normal University, Wuhan 430079,
China  Instituto Galego de Física de Altas Enerxías IGFAE, Universidade de
Santiago de Compostela, E-15782 Galicia, Spain  Yayun He<EMAIL_ADDRESS>Key Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle
Physics, Central China Normal University, Wuhan 430079, China  Guangdong
Provincial Key Laboratory of Nuclear Science, Institute of Quantum Matter,
South China Normal University, Guangzhou 510006, China  Guangdong-Hong Kong
Joint Laboratory of Quantum Matter, Southern Nuclear Science Computing Center,
South China Normal University, Guangzhou 510006, China  Shanshan Cao
<EMAIL_ADDRESS>Institute of Frontier and Interdisciplinary Science,
Shandong University, Qingdao, Shandong 266237, China  Xin-Nian Wang
<EMAIL_ADDRESS>Key Laboratory of Quark and Lepton Physics (MOE) and Institute
of Particle Physics, Central China Normal University, Wuhan 430079, China
Nuclear Science Division Mailstop 70R0319, Lawrence Berkeley National
Laboratory, Berkeley, CA 94740

###### Abstract

A Linear Boltzmann Transport (LBT) Monte Carlo model has been developed to
describe jet propagation and interaction with the quark-gluon plasma (QGP) in
relativistic heavy-ion collisions. A complete set of elastic scattering
processes and medium-induced gluon emissions based on the higher-twist
formalism are incorporated for both jet shower and medium recoil partons. It
has been employed to describe experimental data on large transverse momentum
hadron and jet spectra, correlation and jet substructures in high-energy
heavy-ion collisions. We document in detail the structure of the model and
validation of the Monte Carlo implementations of the physics processes in LBT,
in particular, the inelastic process of medium-induced gluon radiation. We
carry out a comprehensive examination of the jet-medium interaction as
implemented in LBT through energy loss and momentum broadening of a single
hard parton, the energy and transverse momentum transfer from leading partons
to medium-induced gluons and jet-induced medium excitation, and medium
modification of reconstructed jets in a static and uniform medium. With
realistic and event-by-event hydrodynamic medium in heavy-ion collisions, we
compute and compare to experimental data on the jet cone-size dependence of
the single inclusive jet suppression at both Relativistic Heavy-Ion Collider
(RHIC) and the Large Hadron Collider (LHC), the dijet asymmetry at LHC and
$\gamma$-jet correlation at RHIC. Effects of medium-induced gluon emissions
and jet-induced medium excitation on jet observables are systematically
examined. Rescatterings of the radiated gluons and recoil partons with the QGP
are found essential to account for the enhancement of soft particle yield
toward the edge of the jet cone.

##  I Introduction

Exploring properties of the color deconfined state of nuclear matter, known as
the quark-gluon plasma (QGP), is the central goal of relativistic heavy-ion
collision experiments at the Relativistic Heavy-Ion Collider (RHIC) and the
CERN Large Hadron Collider (LHC) Shuryak (2017); Adams  _et al._ (2005); Adcox
_et al._ (2005); Muller  _et al._ (2012). Among various probes of the nuclear
matter, jet quenching reveals the fine structure of the QGP at short distance
scales, since the momentum transfer between energetic jet partons and the
medium can be much larger than the thermal scale given by the local
temperature Akiba  _et al._ (2015); Mehtar-Tani  _et al._ (2012); Armesto  _et
al._ (2012); Casalderrey-Solana and Iancu (2011); Kumar  _et al._ (2020). Many
observables have been proposed to study how jets are modified by the QGP as
they traverse the hot and dense medium, from the nuclear modification factor
of high transverse momentum ($p_{\mathrm{T}}$) hadrons and fully reconstructed
jets, their anisotropic flow coefficients, to the inner structures of jets and
jet-related correlations Connors  _et al._ (2018); Wang and Gyulassy (1992);
Majumder and Van Leeuwen (2011); Qin and Wang (2015); Blaizot and Mehtar-Tani
(2015); Cao and Wang (2021); Cao and Qin (2023). Over the past three decades,
studies on jets have transitioned from understanding their interaction
dynamics with the hot subatomic matter to utilizing them to extract transport
properties of the QGP Burke  _et al._ (2014); Cao  _et al._ (2021); Xie  _et
al._ (2022) with precision. Related investigations have also been extended
from ultrarelativistic heavy-ion collisions to smaller collisions systems Park
_et al._ (2017); Zhao  _et al._ (2020); Liu  _et al._ (2022a); Li  _et al._
(2021) and collisions at lower center-of-mass energies Adamczyk  _et al._
(2018); Wu  _et al._ (2022).

As jet partons plough through the QGP medium, they endure both elastic Bjorken
(1982); Braaten and Thoma (1991); Djordjevic (2006); Qin  _et al._ (2008) and
inelastic Baier  _et al._ (1997); Zakharov (1996); Gyulassy  _et al._ (2000);
Wiedemann (2000); Arnold  _et al._ (2002); Wang and Guo (2001) scatterings
with the medium. These interactions not only lead to nuclear modification of
jets, but also energy-momentum deposition to the medium and depletion from the
medium due to back-reaction. The latter two are collectively known as “jet-
induced medium excitation” or “medium response”, and are essential in
understanding how the lost energy from jet partons are redistributed and
thermalize. It has now been widely accepted that jet-induced medium excitation
is responsible for some of the enhanced soft hadron production associated with
jets Chen  _et al._ (2018, 2020), the increased jet energy within a given
circular annulus at large angle with respect to the jet axis Casalderrey-
Solana  _et al._ (2017); Tachibana  _et al._ (2017); Kunnawalkam Elayavalli
and Zapp (2017); Luo  _et al._ (2018); Park  _et al._ (2019), and recovering
the momentum balance between the two hemispheres traversed by the leading and
subleading jets in a dijet event Chatrchyan  _et al._ (2011); Casalderrey-
Solana  _et al._ (2017). Unique signatures of medium response, such as the
diffusion wake in the opposite direction to jet propagation Yang  _et al._
(2021, 2022) and the enhanced baryon-to-meson ratio within and around the jet
cone Chen  _et al._ (2021); Luo  _et al._ (2023); Sirimanna  _et al._ (2022)
have been proposed and wait for confirmation from more precise experimental
measurements.

Measuring medium response is challenging because its signals related to soft
hadrons are usually buried in the huge background from QGP hadronization.
Despite delicate background subtraction techniques developed in heavy-ion
experiments, jet constituents (shower partons including from medium-induced
radiation) and medium response still inevitably overlap each other in reality.
Therefore, seeking signatures of medium response relies on close comparisons
between realistic model calculations and experimental data. Sophisticated
event generators based on different model assumptions have been developed to
simulate jet interactions with the QGP, such as the virtuality-ordered medium-
modified splittings of highly virtual partons Armesto  _et al._ (2009); Zapp
_et al._ (2011); Casalderrey-Solana  _et al._ (2014); Cao and Majumder (2020),
the time-ordered transport descriptions of elastic and inelastic scatterings
of on-shell partons through the QGP Schenke  _et al._ (2009); Wang and Zhu
(2013), and their combination according to the virtuality scale of each parton
Cao  _et al._ (2017); Kumar  _et al._ (2022). Different treatments of medium
response are also implemented in these models, including perturbative
approximation using recoil partons Zapp  _et al._ (2011); He  _et al._ (2022),
hydrodynamic response to energy deposition Tachibana  _et al._ (2022), and a
concurrent simulation of jet and QGP evolution Chen  _et al._ (2018).

The linear Boltzmann transport (LBT) model He  _et al._ (2015); Cao  _et al._
(2018) has been developed to describe jet propagation in QGP medium and to
study the nuclear modification of jets and jet-induced medium excitation. It
incorporates perturbative descriptions of both elastic and inelastic
scatterings between hard and medium partons. The guiding principle of LBT is
the treatment on the equal footing of the offsprings of jet partons and the
thermal partons scattered out of the QGP background which are known as
“recoil” partons. The production of these recoil partons is also accompanied
by particles holes left behind inside the QGP, known as back-reaction or
“negative” partons. Recoil and “negative” partons constitute the jet-induced
medium excitation within the picture of linear response. This LBT model
provides a satisfactory description of the nuclear modification of high
$p_{\mathrm{T}}$ hadrons Xing  _et al._ (2020), single inclusive jets He  _et
al._ (2019, 2022) and $\gamma$-triggered jets Luo  _et al._ (2018), and has
been embedded inside the JETSCAPE framework Putschke  _et al._ (2019) as a
module of jet evolution at low virtuality. As a follow-up work of Ref. He  _et
al._ (2015), we document a detailed description of the model structure and
Monte Carlo implementations in LBT and their validation in this paper. In
particular, we focus on how the medium-induced gluon bremsstrahlung process is
modeled and how it affects the phase space structure of jet and thermal
partons. Jet energy loss and momentum broadening are systematically studied in
a static and uniform QGP medium. As a supplement to our earlier publications,
additional jet observables in heavy-ion collisions with hydrodynamical QGP
background and event-by-event fluctuating initial condition will be presented
and compared to data at RHIC and LHC. These include the nuclear modification
of single inclusive jets, $\gamma$-triggered jets and dijet asymmetry, and
their jet-cone size dependence.

The remainder of this paper is organized as follows. In Sec. II, we first
present the structure and numerical implementations of elastic and inelastic
scatterings between hard/recoil partons and the QGP within the LBT model. In
Sec. III, we discuss in detail the inelastic scattering rate and the medium-
induced gluon spectrum implemented in LBT. We then validate this model setup
by numerically simulating the energy and angular distributions of hard partons
after both single and multiple scatterings inside a static medium in Sec. IV,
including their dependences on the parton energy, medium temperature and
energy loss mechanism. In Sec. V, we analyze the spacetime evolution of jet-
induced medium excitation in the static medium, such as the angular and
momentum distributions of medium partons excited by jets. The energy loss of
fully reconstructed jets, and the longitudinal and transverse momentum
distributions of the jet constituents are then investigated within a static
medium in Sec. VI. We combine the LBT model with a hydrodynamic simulation of
the QGP evolution in Sec. VII, and calculate the nuclear modification of both
single inclusive jets, $\gamma$-triggered jets and dijets in realistic heavy-
ion collisions at RHIC and LHC energies. The release of the LBT code package
is presented in Sec. VIII. A summary is given in Sec. IX.

##  II The linear Boltzmann transport model

The evolution of the phase space distribution of an incoming jet parton
(denoted by “$a$”) inside a thermal medium is described in LBT according to
the Boltzmann equation,

$\displaystyle p\cdot\partial f_{a}(x,p)=E\,\mathcal{C}_{\mathrm{el+inel}},$
(1)

where the four-momentum of $a$ is given by $p^{\mu}=(E,\vec{p})$, and the
$2\rightarrow 2+n$ scattering processes denoted by the collision integral
$\mathcal{C}_{\mathrm{el+inel}}$ can be decomposed into a sequence of
$2\rightarrow 2$ elastic scatterings and their induced gluon emissions.

The elastic processes $ab\rightarrow cd$ can be simulated according to the
scattering rate,

$\displaystyle\Gamma_{\mathrm{el}}^{ab\rightarrow cd}$ $\displaystyle=$
$\displaystyle\frac{\gamma_{b}}{2E}\int\frac{d^{3}p_{b}}{(2\pi)^{3}2E_{b}}\int\frac{d^{3}p_{c}}{(2\pi)^{3}2E_{c}}\int\frac{d^{3}p_{d}}{(2\pi)^{3}2E_{d}}$
(2) $\displaystyle\times$ $\displaystyle f_{b}(\vec{p}_{b})\left[1\pm
f_{c}(\vec{p}_{c})\right]\left[1\pm
f_{d}(\vec{p}_{d})\right]S_{2}(\hat{s},\hat{t},\hat{u})$ $\displaystyle\times$
$\displaystyle(2\pi)^{4}\delta^{(4)}(p+p_{b}-p_{c}-p_{d})|\mathcal{M}_{ab\rightarrow
cd}|^{2},$

where $b$ represents a thermal parton inside the medium, with a spin-color
degeneracy $\gamma_{b}$ (6 for quarks and 16 for gluons), $c$ and $d$ are the
final state particles of $a+b$ scattering. The distribution function of a
given particle $i$ is denoted by $f_{i}(\vec{p}_{i})$. The thermal
distributions for $b$, $c$ and $d$ in Eq. (2) are $f_{i}=1/(e^{p\cdot u/T}\mp
1)$, with “$-$” sign for gluons and “$+$” for quarks. The local temperature
$T$ and the fluid velocity $u^{\mu}$ are provided by the hydrodynamic
simulation of the QGP.

In our current study, the initial-state-averaged and final-state-summed (over
spin and color degeneracies) matrix elements $|\mathcal{M}_{ab\rightarrow
cb}|^{2}$ are taken from perturbative QCD (pQCD) calculations at the leading
order for all possible scattering channels Auvinen  _et al._ (2010). The
corresponding elastic cross sections are $d\sigma_{ab\rightarrow
cd}/d\hat{t}=|{\cal M}_{ab\rightarrow cd}|^{2}/(16\pi\hat{s}^{2})$, with
$s,t,u$ being the Mandelstam variables. Light flavor partons are assumed to be
massless, and bare masses for heavy quarks are $m_{c}=1.27$ GeV for charm and
$m_{b}=4.19$ GeV for beauty quarks. For scatterings between massless partons,
$|\mathcal{M}_{ab\rightarrow cd}|^{2}$ diverges at small angle
($\hat{u},\hat{t}\rightarrow 0$). This is regulated by a double-$\theta$
function

$\displaystyle S_{2}(\hat{s},\hat{t},\hat{u})=\theta(\hat{s}\geq
2\mu_{\mathrm{D}}^{2})\theta(-\hat{s}+\mu_{\mathrm{D}}^{2}\leq\hat{t}\leq-\mu_{\mathrm{D}}^{2}),$
(3)

where

$\mu_{\mathrm{D}}^{2}=(N_{c}+\frac{N_{f}}{2})\frac{g^{2}T^{2}}{3}$ (4)

is the Debye screening mass, with $g^{2}=4\pi\alpha_{\mathrm{s}}$ being the
strong coupling constant, $N_{c}=3$ and $N_{f}=3$ being the numbers of colors
and flavors, respectively.

The total elastic scattering rate for parton $a$ is

$\Gamma^{a}_{\rm el}\equiv\sum_{b,(cd)}\Gamma_{\rm el}^{ab\rightarrow cd},$
(5)

where the summation is over all possible scattering channels $a+b\rightarrow
c+d$.

For inelastic scattering, we relate its rate to the medium-induced gluon
spectrum as

$\Gamma^{a}_{\mathrm{inel}}=\int
dzdk_{\perp}^{2}\frac{1}{1+\delta^{ag}}\frac{dN^{a}_{g}}{dzdk_{\perp}^{2}dt},$
(6)

with the gluon spectrum taken from the higher-twist energy loss calculation
Guo and Wang (2000); Zhang  _et al._ (2004); Majumder (2012),

$\displaystyle\frac{dN_{g}^{a}}{dzdk_{\perp}^{2}dt}=\frac{2\alpha_{\mathrm{s}}C_{A}\hat{q}_{a}(x)P_{a}(z)k_{\perp}^{4}}{\pi\left({k_{\perp}^{2}+z^{2}m^{2}}\right)^{4}}\,{\sin}^{2}\left(\frac{t-t_{i}}{2\tau_{f}}\right),$
(7)

where $z$ and $k_{\perp}$ are the fractional energy and the transverse
momentum of the emitted gluon with respect to its parent jet parton,
$\alpha_{\mathrm{s}}$ is the strong coupling constant, $C_{A}=N_{c}=3$, $m$ is
the quark mass, $t_{i}$ is the production time of the jet parton or the last
time when the parton has inelastic scattering with gluon radiation, and

$\tau_{f}=\frac{2Ez(1-z)}{k_{\perp}^{2}+z^{2}m^{2}}$ (8)

is the formation time of the emitted gluon. The medium induced splitting
functions $P_{a}(z)$ for $a+g({\rm medium})\rightarrow a+g$ processes are
Zhang and Wang (2003); Schafer  _et al._ (2007),

$\displaystyle P_{q}(z)$ $\displaystyle=$
$\displaystyle\frac{(1-z)[1+(1-z)^{2}]}{z},$ (9) $\displaystyle P_{g}(z)$
$\displaystyle=$ $\displaystyle 2\frac{(1-z+z^{2})^{3}}{z(1-z)}.$ (10)

To avoid divergence at $z\rightarrow 0$ in evaluating the inelastic scattering
rate, a lower cut-off on the fractional energy
$z_{\mathrm{min}}=\mu_{\mathrm{D}}/E$ is applied to the emitted gluon. In Eq.
(6), a 1/2 factor is applied when we evaluate the scattering rate from the
gluon spectrum if the parent jet parton is also a gluon. Note that in Eq. (7),
the time (or path length) dependence of the gluon spectrum results from the
interference between gluon radiation amplitudes from different scatterings. It
was originally derived within the higher-twist formalism to take into account
the interference between the radiation amplitude induced by the hard
scattering (that generates the jet parton) and that induced by the secondary
scattering Guo and Wang (2000). It is applied to multiple scatterings in the
QGP medium in LBT and takes into account the inferences between gluon
emissions induced by different scatterings. This will be referred to as the
Landau–Pomeranchuk–Migdal (LPM) interference. Since not every scattering can
induce a gluon emission, the radiation probability according to Eq. (7) for
the current scattering depends on the distance since the last inelastic
scattering. If induced gluon radiation does not happen for the current
scattering, this distance for the propagating parton will be updated in the
calculation of the radiation probability according to Eq. (7) for the next
scattering as illustrated in the flow chart of LBT in Fig. 1.

The jet quenching coefficient $\hat{q}_{a}(x)$ characterizes the transverse
momentum broadening squared of the jet parton $a$ per unit length,
$\hat{q}_{a}=d\langle q^{2}_{\perp}\rangle_{a}/dt$, and can be evaluated using
Eq. (2) with a weight of
$q_{\perp}^{2}=[\,\vec{p}_{c}-(\vec{p}_{c}\cdot\hat{\vec{p}})\hat{\vec{p}}\,]^{2}$
on the right-hand side. Here $\hat{\vec{p}}$ is the unit vector along parton
$a$’s momentum before scattering. If the quantum statistics in the final state
of $2\rightarrow 2$ scattering is neglected, the elastic scattering rate and
jet transport coefficient can be expressed in terms of the elastic scattering
cross sections,

$\displaystyle\Gamma_{\rm el}^{ab\rightarrow
cd}(x)=\rho_{b}(x)\sigma_{ab\rightarrow cd},$ (11)
$\displaystyle\hat{q}_{a}(x)=\sum_{b,(cd)}\rho_{b}(x)\int
d\hat{t}q_{\perp}^{2}\frac{d\sigma_{ab\rightarrow cd}}{d\hat{t}}.$ (12)

Since pQCD $2\rightarrow 2$ processes are dominated by small angle $t$-channel
processes, the elastic scattering rate, jet transport coefficient and the
inelastic rate are all proportional to the quadratic Casimir $C_{2}(a)$ of the
parton. They are a factor $C_{A}/C_{F}=9/4$ larger for a gluon than a quark.
In the LBT model, we use Eq. (2) to evaluate the elastic scattering rate and
the jet transport coefficient $\hat{q}$ which can be approximated well by Eq.
(12) as shown in Ref. He  _et al._ (2015).

Figure 1: (Color online) Flow chart of the LBT Monte Carlo model. The green
shaded boxes represent modules that can be replaced by users with
alternatives.

The above elastic and inelastic scattering rates are the basic building blocks
of the LBT model. They are computed and tabulated as functions of the local
temperature $T$, parton energy $E$ and the propagation time $t-t_{i}$. This is
done as a part of the initialization of the LBT model and provided as an input
data file as a part of the LBT code. The structure of LBT built on top of
these scattering rates are illustrated in Fig. 1. The default modules in green
shaded boxes can be replaced by alternatives supplied by users and the final
hadronization module can also be customized. Models such as the AMPT Lin  _et
al._ (2005) or Trento Moreland  _et al._ (2015) model provide the initial
conditions for both the hydrodynamic evolution of the bulk medium and spatial
probability distribution of hard binary collisions $\langle N_{\rm
binary}\rangle(\vec{x}_{0\perp})$ in the transverse plane. The default
hydrodynamic model used in the past LBT and this study is the CLVisc model and
can be replaced by users with other alternatives. The hydrodynamic model with
the given initial condition provides the spatial and time evolution of the QGP
medium: temperature $T$ and flow velocity $u$ as functions of spatial
coordinates and time in the collision frame. Initial hard processes and jet
showering are simulated with Monte Carlo models such as Pythia Sjostrand  _et
al._ (2008) and Sharpa Gleisberg  _et al._ (2009) with their initial
production points $\vec{x}_{0}$ sampled from the binary collision distribution
$\langle N_{\rm binary}\rangle(\vec{x}_{0\perp})$ in the transverse plane. The
initial longitudinal position $x_{03}$ is set to 0 for hard processes in high
energy heavy-ion collisions.

For each jet shower parton with momentum $\vec{p}_{i}$ and energy $E_{i}$ (the
less energetic daughter from splitting) along the shower tree of the hard jet
parton with initial momentum $\vec{P}$ and energy $E$, one can assign a
formation time,

$\tau_{fi}=\frac{2E}{P_{\rm
T}^{2}}+\sum_{j=1}^{i}\frac{2E_{j}}{\left[\vec{p}_{j}-(\vec{p}_{j}\cdot\hat{\vec{P}})\hat{\vec{P}}\right]^{2}},$
(13)

where the summation is over all prior splittings along the branching tree. The
formation time for the leading parton in the final splitting of the tree is
assumed to be the same as the less energetic daughter. This parton is allowed
to interact with the medium only after this formation time starting at the
position $\vec{x}_{i}=\vec{x}_{0}+\tau_{fi}\vec{p}_{i}/E_{i}$.

The propagation and interaction of jet shower partons in the medium are
simulated in small incremental time steps. For each time step $\Delta t$ in
the collision frame, one first gets the elastic and inelastic scattering rates
from the data table at the local temperature $T$ for parton’s energy
$E_{i}^{\prime}$ in the local comoving frame of the QGP.

For small time step $\Delta t$ in the global computational frame, the average
number of elastic and inelastic scatterings during the time interval are then
given by

$\displaystyle\langle\Delta N^{a}_{\mathrm{el}}\rangle(\vec{x}_{i})$
$\displaystyle=$ $\displaystyle\Gamma^{a}_{\mathrm{el}}(\vec{x}_{i})\Delta
t^{\prime},$ (14) $\displaystyle\langle\Delta
N^{a}_{\mathrm{inel}}\rangle(\vec{x}_{i})$ $\displaystyle=$
$\displaystyle\Gamma^{a}_{\mathrm{inel}}(\vec{x}_{i})\Delta t^{\prime},$ (15)

respectively, where $\Delta t^{\prime}=\Delta t\,p_{i}\cdot u/E_{i}$ is the
time interval in the local comoving frame. The probabilities for elastic and
inelastic scattering in this time interval are

$\displaystyle P^{a}_{\mathrm{el}}$ $\displaystyle=$ $\displaystyle
1-e^{-\langle\Delta N^{a}_{\mathrm{el}}\rangle(\vec{x}_{i})},$ (16)
$\displaystyle P^{a}_{\mathrm{inel}}$ $\displaystyle=$ $\displaystyle
1-e^{-\langle\Delta N^{a}_{\mathrm{inel}}\rangle(\vec{x}_{i})}.$ (17)

These scattering probabilities become $\Gamma_{\mathrm{el}}^{a}\Delta
t^{\prime}\ll 1$ and $\Gamma_{\mathrm{inel}}^{a}\Delta t^{\prime}\ll 1$, if
$\Delta t^{\prime}$ is sufficiently small. To combine elastic and inelastic
processes, the total scattering rate is given by
$\Gamma_{\mathrm{tot}}^{a}(\vec{x}_{i})=\Gamma_{\mathrm{el}}^{a}(\vec{x}_{i})+\Gamma_{\mathrm{inel}}^{a}(\vec{x}_{i})$,
and the total number of scattering, elastic and inelastic, is

$\langle\Delta N^{a}_{\mathrm{tot}}\rangle(\vec{x}_{i})=\Gamma_{\rm
tot}^{a}(\vec{x}_{i})\Delta t^{\prime}.$ (18)

The total scattering probability is,

$P_{\mathrm{tot}}^{a}=1-e^{-\langle\Delta
N^{a}_{\mathrm{tot}}\rangle(\vec{x}_{i})}=P^{a}_{\mathrm{el}}(1-P^{a}_{\mathrm{inel}})+P^{a}_{\mathrm{inel}},$
(19)

which can be understood as a sum of the probability of pure elastic scattering
process without gluon emission (first term) and the probability of inelastic
scattering with at lest one gluon emission (the second term).

Using these probabilities within a given time step $\Delta t$, we first decide
whether a given jet parton scatters with the medium thermal partons. If it
does not scatter, it will propagate in the medium for the next time step
$\vec{x}_{i}\rightarrow\vec{x}_{i}+\Delta t\,\vec{p}_{i}/E_{i}$ in the
collision frame. The scattering rates at the new position will be used to
determine the probabilities for scattering. Such propagation is repeated until
a scattering occurs.

We then decide whether this scattering is elastic or inelastic. If an elastic
scattering occurs based on this probability, the branching ratio
$\langle\Delta N_{\mathrm{el}}^{ab\rightarrow cd}\rangle/\langle\Delta
N_{\mathrm{el}}^{a}\rangle$ is used to select a particular scattering channel,
i.e., the species of $b$, $c$ and $d$. With a given channel, the differential
scattering rate based on Eq. (2) is then applied to determine the momenta of
$b$, $c$ and $d$. For light flavor jets, we assign the parton between $c$ and
$d$ with higher energy as the final state jet parton, while the other as the
“recoil parton” scattered out of the medium background. Meanwhile, $b$
represents the energy-momentum depletion from the medium, denoted as “back-
reaction” or “negative parton” in LBT. Their propagation and interaction with
the medium will be simulated in subsequent time steps according to the above
procedure in the LBT model in order to provide a complete picture of jet-
medium interaction and guarantee the energy-momentum conservation of the
entire system.

If an inelastic scattering happens with a given number of time steps based on
the given probability, we sample the number of emitted gluons using the
Poisson distribution with a mean number of emissions $\langle\Delta
N^{a}_{\mathrm{inel}}\rangle$. We then assign the energy-momentum of each
gluon using its differential spectrum given by Eq. (7). Note that the time
$t_{i}$ in Eq. (7) is the time of the last inelastic scattering. Since the
emitted gluons are induced by elastic scattering within our framework, we
decompose each inelastic scattering ($2\rightarrow 2+n$) into an elastic
scattering ($2\rightarrow 2$) plus a multiple gluon splitting ($1\rightarrow
n$) process. In addition, the final state jet parton directly sampled from Eq.
(2) is on-shell and is not able to split, we need to adjust the energy-momenta
of the $2+n$ final state particles obtained from Eqs. (2) and (7) together to
ensure the energy-momentum conservation of the entire $2\rightarrow 2+n$
process. Details on this adjustment are presented in Appendix A. After each
inelastic scattering the initial time $t_{i}$ in Eq. (7) is reset. Note that
in each time step, we determine the momenta of final jet partons, radiated
gluons, recoil and “negative” partons in the local comoving frame first before
boosting them back into the global collision frame. All these partons will
continue to propagate and interact with the medium according to the above
procedure until the local temperature reaches the QCD phase transition
temperature $T<T_{c}$ according to the equation of state (EoS) used in the
hydrodynamic model. In the current version of LBT model, parton-medium
interaction below $T_{\mathrm{c}}$ is neglected. For the study of jets in
heavy-ion collisions, we currently reconstruct jets with the final state
partons using a jet finding algorithm. Since “negative” partons represent the
depletion of energy and momentum in the medium due to the back-reaction, they
should be subtracted from the jet energy within a given jet cone. In
principle, the final jet shower partons, radiated gluons, recoil partons and
“negative” partons can be hadronized according to a hadronization model as
done within the JETSCAPE framework Putschke  _et al._ (2019). This is not
included in the present work and needs to be implemented in our future effort.

In the following sections till the end of Sec. VI, we will first validate the
LBT model for the propagation of an on-shell single jet parton that enters a
static medium at an initial time of $t_{0}=0$, and develops into a jet that
consists of recoil, “negative” partons and radiated gluons from elastic and
inelastic scatterings. The formation time in Eq. (13) does not apply to the
initial single parton in these brick tests. Jet observables in realistic
heavy-ion collisions will be discussed in Sec. VII, in which initial jet
shower partons are produced from Pythia simulations and start interacting with
the QGP at their formation times

##  III The inelastic process

###  III.1 Inelastic scattering rate

The evolution of energetic partons in heavy-ion collisions is mainly governed
by their scattering rates with the QGP. The elastic scattering rates for
different scattering channels, including their energy and temperature
dependences and the angular distribution of the recoil partons, were presented
in detail in an earlier publication He  _et al._ (2015). In this subsection,
we concentrate on the medium-induced gluon spectra, which determine the rates
for inelastic processes within LBT.

For massless partons (light quarks and gluons), one can complete the
integration over the transverse momentum of the radiated gluon in Eq. (7) and
obtain the differential inelastic rate,

$\frac{d\Gamma^{a}_{\rm
inel}}{dz}=\frac{\alpha_{\mathrm{s}}C_{A}\hat{q}_{a}}{1+\delta^{ag}}\frac{(t-t_{i})}{4E}\frac{P_{a}(z)}{z(1-z)}h(v),$
(20)

where the kinetic limit for $k_{\perp}$ is set as $k_{\perp}\leq\min(z,1-z)E$
and

$\displaystyle h(v)$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{0}^{v}dx\frac{\sin^{2}x}{x^{2}},$ (21)
$\displaystyle v$ $\displaystyle=$
$\displaystyle\frac{(t-t_{i})E}{4z(1-z)}\min(z^{2},(1-z)^{2}),$

which is normalized with the limit $h(\infty)=1$. Since

$\left[\frac{P_{a}(z)}{z(1-z)}h(v)\right]_{z\rightarrow
0}=\frac{2}{\pi}\frac{(t-t_{i})E}{2z},$ (22)

the above differential inelastic rate has a logarithmic infrared divergence.
We therefore impose an infrared cut-off $z\geq\mu_{D}/E$ when we calculate the
integrated inelastic scattering rate $\Gamma^{a}_{\rm inel}$. Shown in Figs. 2
and 3 are $h(v)$ as a function of $z$ and $(t-t_{i})E/4$ and $h(v)$ as a
function of $z$ for four different values of $(t-t_{i})E$, respectively. One
can see that it approaches to its asymptotic value $h(\infty)=1$ when
$(t-t_{i})E$ is very large and the effect of the kinetic restriction on the
transverse momentum $k_{\perp}\leq zE$ is negligible.

Figure 2: (Color online) $h(v)$ defined in Eq. (21) as a function of
fractional momentum $z$ and $(t-t_{i})E$. Figure 3: (Color online) $h(v)$
defined in Eq. (21) as a function of fractional momentum $z$ for four
different values of $(t-t_{i})E$.

Figure 4: (Color online) The effective medium-induced parton splitting
function $P_{a}(z)h(v)/[z(1-z)]$ for a quark (upper panel) and a gluon (lower
panel) for four different values of $(t-t_{i})E=1,4,20,60$ (from bottom to
top) as compared to the asymptotic limit (dashed line) when
$(t-t_{i})E\rightarrow\infty$.

The effective medium-induced parton splitting functions
$P_{a}(z)h(v)/[z(1-z)]$ are shown in Fig. 4 for four different values of
$(t-t_{i})E$, which approaches to the asymptotic form (dashed line) when
$(t-t_{i})E$ is large.

Figure 5: (Color online) The rate of gluon emission from a gluon jet parton at
$t-t_{i}=1$ fm/$c$ as a function of temperature $T$ and initial parton energy
$E$. The strong coupling constant $\alpha_{\rm s}$ is set to 0.3.

Integrating over the momentum fraction, we can obtain the total inelastic rate
or the number of radiated gluons per unit time as shown in Fig. 5 for a gluon
at $t=1$ fm/$c$ as a function of the local temperature $T$ and parton energy
$E$. This rate has only a logarithmic energy dependence at large $E$, which is
contributed by $\hat{q}_{a}$. Since the effective parton splitting functions
have an asymptotic limit $P_{a}(z)/[z(1-z)]$ for large $(t-t_{i})E>4$, one can
obtain the integrated inelastic scattering rate analytically,

$\Gamma^{a}_{\rm
inel}=(t-t_{i})\frac{\alpha_{\mathrm{s}}C_{A}\hat{q}_{a}}{2\mu_{\rm
D}}\left[1+\frac{\mu_{\rm D}}{E}\ln\frac{\mu_{\rm D}}{E}+{\cal
O}(\frac{\mu_{\rm D}}{E})\right].$ (23)

For $(t-t_{i})E<4$, the integrated inelastic rate is instead

$\Gamma^{a}_{\rm
inel}\approx(t-t_{i})^{2}\frac{\alpha_{\mathrm{s}}C_{A}\hat{q}_{a}}{4\pi}\ln\frac{E}{\mu_{\rm
D}}.$ (24)

Using small angle approximation for the parton scattering cross section, the
jet transport coefficient is He  _et al._ (2015)

$\hat{q}_{a}=C_{2}(a)\frac{42\zeta(3)}{\pi}\alpha_{\rm
s}^{2}T^{3}\ln\frac{c_{a}ET}{\mu_{\rm D}^{2}},$ (25)

where $\zeta(3)\approx 1.202$ is the Apéry’s constant, $C_{2}(q)=4/3$,
$C_{2}(g)=3$, $c_{q}=5.8$ and $c_{g}=5.6$. The energy and temperature
dependence of $\Gamma^{a}_{\rm inel}$ at large $(t-t_{i})E$ is therefore
determined mostly by that of $\hat{q}_{a}$ for a fixed value of the strong
coupling constant $\alpha_{\rm s}$.

Figure 6: (Color online) Fractional energy distributions of the radiated
gluons in $q\rightarrow q+g$ (upper) and $g\rightarrow g+g$ (lower) processes
with an initial parton energy $E_{0}=20$ GeV in a static medium at temperature
$T=400$ MeV from semi-analytical calculations (black open circles), direct MC
sampling (red dashed lines) and LBT simulations for $2\rightarrow 3$ (green
solid lines) and $2\rightarrow 2+n$ $(n\geq 1)$ (blue solid lines) .

###  III.2 Medium-induced gluon spectra

For each inelastic scattering, we use Eq. (7) to sample the momentum fraction
$z$ and the transverse momentum $\vec{k}_{\perp}$ of the radiated gluon. Since
in each inelastic process in LBT simulations, the medium-induced gluon
emissions are always accompanied by a $2\rightarrow 2$ elastic scattering
process and the jet parton will also experience transverse momentum transfer
(or broadening) and elastic energy loss. The energy-momentum conservation in
this inelastic process with one or multiple gluon emissions should be
considered together. We adopt a procedure as described in Appendix A to
maintain the global energy-momentum conservation in $2\rightarrow 2+n(g)$
processes.

Shown in Fig. 6 are the normalized fractional energy distributions of gluons,

$\frac{dN}{dz}=\int_{t_{\rm min}}^{t_{\rm max}}dt\frac{d\Gamma^{a}_{\rm
inel}}{dz}\left/\int_{t_{\rm min}}^{t_{\rm max}}dt\Gamma^{a}_{\rm
inel}\right.,$ (26)

emitted from a quark (upper panel) or gluon (lower panel) jet parton with an
initial energy of $E_{0}=20$ GeV in a single inelastic process with different
treatments of global energy-momentum conservation. In order to compare to
results from semi-analytical calculations, we consider a static and uniform
medium at a constant temperature $T=0.4$ GeV and set the strong coupling
constant at $\alpha_{\rm s}=0.3$. In the figure, the black open circles
represent results from numerical integration of Eq. (7) over $k_{\perp}^{2}$
and time between $t_{\rm min}=1$ fm/$c$ and $t_{\rm max}=6$ fm/$c$ with
$t_{i}=0$. The red dashed lines represent results from direct Monte-Carlo
sampling according to Eq. (7), which are in exact agreement with the circles
as expected. The gluon spectra become softer and deviate from the spectrum
given by Eq. (7) when energy-momentum conservation is considered for
$2\rightarrow 3$ process as shown by the green curves. This softening is
caused by the elastic energy loss of the jet parton when considered together
with the induced gluon radiation. When multiple gluon emissions are taken into
account in $2\rightarrow 2+n$ processes, the emitted gluon spectra (the blue
lines) are further shifted towards lower energies because of the successive
energy loss of the jet parton when it emits multiple gluons. Comparing the
upper and the lower panel, one observes that the energy spectrum of emitted
gluons peaks at low energy for a quark jet, but peaks at both small and large
energy for a gluon jet. This can be understood with the $\sim 1/z$ splitting
function for the $q\rightarrow qg$ process, while the splitting is symmetric
$\sim 1/[z(1-z)]$ between the emitted gluon and the jet gluon for
$g\rightarrow gg$. The symmetry around $z=0.5$ can be seen in the semi-
analytical result, but is broken when the realistic $2\rightarrow 2+n$ process
is taken into account.

Figure 7: (Color online) The angular distribution of medium-induced radiated
gluon for a typical parton energy and length $\sqrt{\Delta tE}$ and Debye
screening mass $\mu_{D}/E$.

In a static and uniform medium, one can complete the time integral in Eq. (7)
and obtain the angular distribution of the radiated gluon,

$\displaystyle\frac{dN_{g}^{a}}{d\theta}$ $\displaystyle=$
$\displaystyle\frac{2\alpha_{\mathrm{s}}C_{A}\hat{q}_{a}}{\pi\sqrt{E}}\frac{\Delta
t^{5/2}\cos\theta}{(\sqrt{\Delta
tE}\sin\theta)^{3}}\int_{\mu_{D}/E}^{1}dz\frac{P_{a}(z)}{z^{2}}$ (27)
$\displaystyle\times$ $\displaystyle\left\\{1-\frac{\sin\left(\Delta
tEz\sin^{2}\theta/[2(1-z)]\right)}{\Delta
tEz\sin^{2}\theta/[2(1-z)]}\right\\},$

where the polar angle is defined as $\sin\theta=k_{\perp}/(zE)$. Note that
this polar angle is different from the relative angle between the final two
partons after the splitting. For a gluon jet, the angle is defined for the
less energetic gluon in the final state $\sin\theta=k_{\perp}/(\min[z,1-z]E)$
and the result is twice of the integration over $z\leq 1/2$.

Unlike gluon radiation from a virtual parton in vacuum that has a collinear
divergence $dN_{g}/d\theta\sim 1/\theta$, medium-induced gluons tend to be
emitted with a finite angle $\theta\sim 1/\sqrt{\Delta tE}$ with respect to
the direction of the incoming jet parton. Small angle emissions whose
formation times are longer than the medium length are suppressed due to the
LPM interference as shown in Fig. 7. The angular distribution vanishes as
$dN_{g}\sim\theta$ when $\theta\ll 1/\sqrt{\Delta tE}$. At large angles, gluon
emission becomes incoherent without LPM interference for large values of
$\sqrt{\Delta tE}$, and the angular distribution becomes
$dN_{g}/d\theta\sim\cos\theta/\sin^{3}\theta$, corresponding to the transverse
momentum spectrum $dN_{g}/dk_{\perp}^{2}\sim 1/k_{\perp}^{4}$.

Figure 8: (Color online) Polar angle distributions of the radiated gluons in
$q\rightarrow q+g$ process with an initial parton energy $E_{0}=20$ GeV in a
static medium at temperature $T=400$ MeV from semi-analytical calculations
(black open circles), direct MC sampling (red dashed line), LBT simulations
for $2\rightarrow 3$ (green solid line) and $2\rightarrow 2+n$ $(n\geq 1)$
(blue solid line) processes.

In Fig. 8, we show the normalized polar angle distribution of the emitted
gluon with respect to a parent quark from LBT model. Distributions from a
parent gluon is very similar. Our Monte Carlo sampling faithfully reproduce
the semi-analytic results for medium-induced gluon emissions. When elastic
scatterings are considered together with the gluon radiation in the
$2\rightarrow 3$ process, the global energy-momentum conservation that
suppresses the gluon emission with large momentum fraction (as shown in Fig.
6) also enhances the polar angle distribution at small angle slightly. When
multiple gluon emissions are allowed in LBT simulations (the blue curves), the
average polar angles become larger compared to the case where only single
gluon emission is allowed (the green curves) due to radiative broadening.

##  IV Parton energy loss and transverse momentum broadening

###  IV.1 Parton energy loss

Figure 9: (Color online) The time dependence of the fractional energy loss
$\Delta E/E_{0}$ (upper panel) and energy loss per unit length $dE/dz$ (lower
panel) of a quark (up and down triangle) or a gluon (circle and square) with
initial energy $E_{0}=100$ GeV through a static medium with temperature
$T=300$ MeV, with only elastic (black and blue) and elastic + gluon radiation
processes (red and magenta).

Figure 10: (Color online) Quark energy loss per unit length as a function of
time in a static medium with temperatures $T=200$ MeV (upper panel) and
$T=400$ MeV (lower panel), for different initial energies $E_{0}$. Both
elastic and gluon radiation processes are included.

Figure 11: (Color online) The same as Fig. 10 except for a gluon.

To investigate the energy loss and transverse momentum broadening of a
propagating parton within the LBT model, we again consider a static and
uniform medium with no flow at a constant temperature in this section. We
start with the energy loss of a single parton in Fig. 9 as a function of its
propagation time (or distance) in a static medium with $T=300$ MeV. The parton
is initialized with 100 GeV energy along the longitudinal ($z$) direction. In
the upper panel, we present the accumulated fractional energy loss of the
parton, compared between a quark and a gluon jet parton and different
scattering processes. One can clearly observe a stronger energy loss of a
gluon than a quark due to the larger color factor of the former. Compared to
the energy loss from only elastic scatterings, the parton energy loss
significantly increases when inelastic scattering, or medium-induced gluon
radiation, is included. This indicates the dominating contribution from gluon
radiation to parton energy loss. In the lower panel of Fig. 9, we show the
parton energy loss per unit length $dE/dz$. For elastic scattering, a constant
energy loss rate $dE/dz$ is reached after the traversed length of the parton
surpasses its mean free path, resulting in a linear increase of the
accumulated parton energy loss as shown in the upper panel. In contrast, when
induced gluon radiation process is incorporated, the energy loss rate $dE/dz$
increases with the path length almost linearly for a rather long time until
the mean formation time of single gluon emission is reached. It becomes less
than linear in the distance dependence afterwards. As a result, the
corresponding accumulated parton energy loss shown in the upper panel
increases almost quadratically with respect to the path length at early times
and returns to a slower length dependence at later times. This change of the
energy loss rates with time is caused by the energy dependence of the energy
loss. Since energy loss rate increases with parton energy, it should slow down
as the parton loses more energy. In addition, such an energy dependence of the
energy loss also leads to the change of the ratio between the energy loss of a
gluon and a quark. The energy loss rate of a gluon is larger than that of a
quark by a factor of $C_{A}/C_{F}=9/4$ at early times as seen in Fig. 9. At
later times, the energy loss of these hard partons, which is more for a gluon
than a quark, leads to the reduction of further energy loss and narrows the
difference of energy loss between a quark and a gluon. Note that if one does
not reset $t_{i}$ in Eq. (7) to zero after a radiative gluon is formed and
ignores the variation of the parton energy as it propagates through the
medium, the accumulated energy taken by emitted gluons then keeps increasing
quadratically with respect to the path length and the relative difference
between gluon and quark energy loss remains at $C_{A}/C_{F}$ Cao  _et al._
(2018).

The dependence of parton energy loss on the parton energy and the medium
temperature is further illustrated in Fig. 10 for a quark and Fig. 11 for a
gluon. In each figure, the upper plot is for a static medium at $T=200$ MeV
and the lower plot is for $T=400$ MeV. In each plot, results from different
initial parton energies are compared. Contributions from both elastic
scattering and gluon radiation are included here. Similar to the findings
previously presented in Fig. 9, the increase of parton energy loss per unit
length $dE/dz$ is rapid at early times and then significantly slows down and
even starts to decrease with time as multiple gluon emissions commence and the
hard parton’s energy becomes smaller. This variation of the energy loss rate
is more prominent for lower energy partons and higher medium temperature,
considering the larger fractional energy loss in these scenarios. For the same
reason, it is more prominent for gluons than for quarks.

###  IV.2 Transverse momentum broadening

Figure 12: (Color online) Transverse momentum distribution of a quark with an
initial energy $E_{0}=100$ GeV at different evolution times through a static
medium at temperature $T=300$ MeV. The upper panel includes only elastic
scatterings, while the lower panel includes both elastic scatterings and gluon
emissions.

Figure 13: (Color online) The same as Fig. 12 except for a gluon.

As discussed in Sec. II, the transverse momentum broadening of a hard parton,
or the $\hat{q}$ parameter, is a key quantity that governs the inelastic
scattering rate of a parton. To investigate this broadening of transverse
momentum ($k_{\mathrm{T}}$) of a hard parton traversing a QGP medium, we
present the $k_{\mathrm{T}}$ distributions of the leading parton at different
times through a static medium in Fig. 12 for an initial quark and Fig. 13 for
an initial gluon. In each figure, the upper panel shows the $k_{\rm T}$
distributions due to only elastic scatterings, while the lower panel shows the
distributions due to both elastic and inelastic scatterings. Here, the hard
parton is initialized with an energy $E_{0}=100$ GeV along the $z$ direction,
and the medium temperature is set as $T=300$ MeV. Therefore, starting as a
$\delta$-function at $k_{\mathrm{T}}=0$, the transverse momentum distribution
becomes broader as time evolves. Comparing the upper and lower panels, one
observes a faster broadening of this $k_{\mathrm{T}}$ after the gluon
radiation process is included. The non-exponential $k_{\rm T}$ distributions
due to inelastic scatterings is caused by the restriction on small $k_{\rm T}$
gluon radiations that have long formation times due to the LPM interference.
Comparing Fig. 12 and Fig. 13, one observes a broader $k_{\mathrm{T}}$
distribution for a gluon than for a quark due to the difference of their jet
transport coefficients.

##  V Jet induced medium excitation

After the discussion of energy loss and transverse momentum broadening of a
propagating parton in the last section, we will study the evolution of the jet
shower in the phase space from an initial hard parton propagating in a static
and uniform QGP medium. The partons in the jet shower include contributions
from the medium-modified leading parton, gluons from medium-induced
radiations, recoil partons and “negative” partons. As previously discussed,
the recoil and “negative” partons collectively constitute the jet-induced
medium excitation in the LBT model.

###  V.1 Conic structures in jet-induced medium excitation

Figure 14: (Color online) Energy density distribution of a jet shower
originating from a 100 GeV quark at $x=y=z=0$ and propagating along the $z$
direction in a static medium at $T=300$ MeV. Results are shown at $t=6$ fm/$c$
with contributions from only elastic scatterings with rescatterings (upper
panel), elastic + inelastic scatterings without rescatterings (middle panel),
and elastic + inelastic scatterings with rescatterings (lower panel).

Figure 15: (Color online) The same as Fig. 14 except for a gluon.
Figure 16: (Color online) (Color online) The same as the contour plot in Fig.
14 (c) except for the energy density carried by recoil partons, “negative”
partons, radiated gluons and their sum. Note the color scales in upper and
lower panels are different for clarity.

In Fig. 14, we first present the energy density distribution of a jet shower,
containing radiated gluons from parton-medium interaction, recoil and
“negative” partons but excluding the leading parton, initiated by a hard quark
with an initial energy of $E_{0}=100$ GeV when it propagates through a static
medium at a temperature of $T=300$ MeV. The hard quark starts its propagation
at the origin ($x=y=z=0$) with its momentum along the longitudinal ($z$)
direction. The energy density is calculated at the time of $t=6$ fm/$c$ in the
coordinate space ($z$-$r$ plane) where $r$ denotes the transverse distance
$r=\sqrt{x^{2}+y^{2}}$. The contribution from “negative” partons is subtracted
from the total energy of radiated gluons and recoil partons. The energy
density distributions are averaged over many events. In the figure, the upper
panel is generated with the LBT simulations that only include the elastic
scattering processes, in which rescatterings of recoil partons with the medium
are allowed; the middle panel includes both elastic and inelastic processes,
with rescatterings of radiated gluons and recoil partons switched off; while
the lower panel represents the full LBT simulations that include both elastic
and inelastic scatterings and allow rescatterings of radiated gluons and
recoil partons.

In the panels that include rescatterings, conic structures of the energy
density distribution can be observed along the path of the jet parton. In the
upper panel that only includes elastic scatterings, this energy density
distribution can approximate the Mach-cone structure of the medium response
induced by the jet parton, where the wave front is produced by the propagation
and medium interaction of recoil partons, corresponding to the hydrodynamic
evolution of the energy deposited by the jet parton. The negative net
distribution in the direction opposite to the jet propagation, known as the
“diffusion wake”, is caused by the energy depletion when thermal partons are
scattered out of their original phase space by the jet parton. Comparing the
upper and lower panel, we observe significant enhancement of both the
magnitude and the width of shock wave in the energy density distribution when
inelastic scatterings are introduced, because the medium-induced gluons can
further scatter with the medium and make the jet-induced medium excitation
stronger and broader. Comparing the middle and the lower panels, we observe
the wave front of the energy distribution is close to an arc shape when
rescatterings are not included, a conic structure occurs and the shock wave is
more diffused when rescatterings of the radiated gluons and recoil partons
with the medium are allowed. These rescatterings with the medium slow down the
average speed of parton propagation, which is equivalent to an effective
velocity of sound in a hydrodynamic description of Mach-cone evolution that is
smaller than the velocity of light (for massless partons) when interaction is
considered in the calculation of the equation of state (EoS). Rescatterings
also excite additional partons out of the medium and further broaden the shock
wave induced by the jet parton much like the effect of viscosity in the
hydrodynamic description Similar features can also be seen in Fig. 15 for a
jet shower initiated by a gluon. Comparing Fig. 14 and Fig. 15, we see a
stronger medium excitation induced by a gluon jet than by a quark jet due to
the stronger interactions between a gluon and the medium. Note that the conic
structures obtained from LBT simulations here are not rigorously Mach cones.
Mach cones in a rigorous sense are only present in the hydrodynamic response
to jet energy loss.

To better illustrate the energy density deposited into and depleted from the
medium, we present in Fig. 16 the energy density contour plots from recoil
partons, “negative” partons, radiated gluons and their sum. One observes that
the distribution of “negative” partons overlaps with that of recoil partons
and radiated gluons. They reduce the total energy density carried by recoil
partons and radiated gluons around the forward direction of jet propagation.
They become dominant in the backward region where the net total energy density
becomes negative. This flow of negative net energy density is often referred
to as the diffusion wake trailing the propagating jet parton.

Figure 17: (Color online) Energy distribution of partons inside a jet shower
at different times developed from a 100 GeV quark in a static medium at
$T=300$ MeV, with only elastic scatterings with rescatterings (upper panel),
elastic + inelastic scatterings without rescatterings (middle panel), and
elastic + inelastic scatterings with rescatterings (lower panel).

Figure 18: (Color online) The same as Fig. 17 except for an initial gluon.

###  V.2 Energy distribution

In the LBT simulations, by subtracting the energy-momentum of “negative”
partons from that of radiated gluons and recoil partons, we can rigorously
respect the energy-momentum conservation in each scattering between a hard
parton and a thermal medium parton. Apart from the energy taken by the medium-
induced gluon radiations, the energy lost from the leading parton is
transferred to the jet-induced medium excitation, i.e., the energy of recoil
partons subtracted by that of “negative” partons. To illustrate how the energy
is transferred, we present in Fig. 17 the energy distribution of all partons
within a jet shower as time evolves. Here, the jet shower is initiated by a
100 GeV quark through a static medium with a temperature of 300 MeV. From the
top to the bottom panel, we compare simulations including only elastic
scatterings with rescatterings of recoil partons, both elastic and inelastic
scatterings without rescatterings of recoil partons and radiated gluons, and
both elastic and inelastic scatterings with rescatterings of recoil partons
and radiated gluons in the full LBT model. At early times, the energy
distribution peaks at both the high energy and low energy ends, with the
former contributed by the leading parton close to its initial energy and the
latter contributed by medium response around the temperature scale. As time
evolves, further energy loss of the leading parton suppresses the peak at the
high energy end and broadens the high-energy peak towards the lower energy
region. Meanwhile, successively induced medium excitation enhances the peak at
the low energy end. Comparing the upper and lower panel, we observe the gluon
emission process significantly accelerates the reduction and broadening of the
high energy peak, while enhances the low energy peak from the medium response.
Comparing the middle and lower panel, we observe rescatterings of the radiated
gluons and recoil partons also enhance the effect of medium response at the
low energy end, though have no impact on the leading parton distribution at
the high energy end. In Fig. 18, we present similar energy distribution of
partons inside a jet shower initiated by a hard gluon. Due to the stronger
interaction of a gluon with the medium than a quark, the energy spectrum of
partons inside a gluon jet is broadened faster than inside a quark jet.

Figure 19: (Color online) Projection from azimuthal angle $\varphi$ relative
to the jet parton direction $y$ to the azimuthal angle $\phi$ in the
transverse plane of the collision frame.

###  V.3 Angular distribution

To have a more detailed illustration of how the energy is transferred from the
leading parton to soft components inside a jet, and how a jet gets broadened
as it evolves through the medium, we look at the azimuthal angle distributions
of partons inside a medium-induced jet shower. In heavy-ion collisions, the
most considered case of jet propagation is in the transverse direction with
jet’s rapidity $\eta$ set to 0. To imitate such propagation and the
corresponding azimuthal angle distribution of the produced partons, one can
imagine the longitudinal direction in a static and uniform medium is along
such a transverse direction ($y$) in real heavy-ion collisions. One then can
project a symmetric azimuthal angle ($\varphi$) distribution with respect to
the initial jet parton direction to the angular ($\phi$) distribution in the
transverse plane in the cylindrical frame of the heavy-ion collisions as
illustrated in Fig. 19. This project can be written as

$\frac{dN}{d\phi}=\int d\theta d\varphi\frac{dN}{d\theta
d\varphi}\delta\left(\phi-\arctan[\tan\theta\cos\varphi]\right),$ (28)

where $\theta$ and $\varphi$ are the polar and azimuthal angles, respectively,
relative to the initial jet parton direction $y$. We refer to the angle $\phi$
in the transverse plane of the collisions as the projected azimuthal angle and
the distribution as the projected azimuthal distribution.

Shown in Figs. 20 \- 22 are the projected azimuthal angle distributions of
jet-induced medium partons (radiated gluons and recoil partons with “negative”
partons subtracted) within different ranges of the parton energy $p_{\rm
T}=\sqrt{p_{x}^{2}+p_{y}^{2}}$. The jet parton originates from a 100 GeV quark
and evolves inside a static medium with a constant temperature of $T=300$ MeV.
In each figure, the upper panel shows the angular distributions at $t=2$
fm/$c$ and the lower panel at $t=6$ fm/$c$.

Figure 20: (Color online) Projected azimuthal angle distribution of partons
within different energy intervals inside a jet shower developed from a 100 GeV
gluon in a static medium at $T=300$ MeV. The upper panel shows results at time
$t=2$ fm/$c$, and the lower panel at $t=6$ fm/$c$. Only elastic scatterings
are included, with rescatterings of recoil partons.

Figure 21: (Color online) The same as Fig. 20 except that both elastic and
inelastic scatterings are included without rescatterings of radiated gluons
and recoil partons.

Figure 22: (Color online) The same as Fig. 20 except that both elastic and
inelastic scatterings are included with rescatterings of radiated gluons and
recoil partons.

In Fig. 20, we first study the angular distributions of soft partons produced
by pure elastic scatterings. One can observe that the net distribution of very
low energy partons ($p_{\rm T}$=0 – 1 GeV) is mostly negative, because this
range is dominated by the “negative” partons left behind the jet parton inside
the medium when medium partons are scattered out into recoil partons. In the
momentum space, these “negative” partons center around the direction opposite
to the jet. In contrast, the distribution of recoil partons peaks around the
jet parton direction in the range (1 – 2 GeV) that is slightly higher than the
thermal energy scale considering the energy they gain from scatterings with
the leading parton. In a small angle scattering approximation, the average
energy transfer per scattering is He  _et al._ (2015),

$\langle\delta E_{\rm el}\rangle=\frac{1}{\Gamma^{a}_{\rm
el}}\left\langle\Gamma^{a}_{\rm
el}\frac{q_{\perp}^{2}}{2\omega}\right\rangle\approx\frac{18\pi^{3}}{84\zeta(3)}\alpha_{\rm
s}T\ln\frac{2.6ET}{\mu_{\rm D}^{2}}.$ (29)

This energy scale for recoil partons should increase logarithmically with the
jet parton’s energy. For $\alpha_{\rm s}=0.3$, $E=100$ GeV and and $T=300$
MeV, $\delta E_{\rm el}\approx 2.5$ GeV, consistent with what we see in Fig.
20.

This distribution of recoil partons exhibits a double peak structure in the
projected azimuthal angle at early times, which results from the typical
transverse momentum transfer $\mu_{\rm D}$ between the hard parton and the
medium constituent in the LBT model simulations. As time evolves, more recoil
and “negative” partons are generated. Meanwhile, the two peaks in the
projected azimuthal angle distribution of recoil partons are smeared and in
the end merge due to the successive rescatterings of both the jet parton and
recoil partons with the medium. Although the impact of $\mu_{\mathrm{D}}$ on
the particle distribution here can be easily diminished by other factors of
jet-medium interactions, it worths further investigation considering that
$\mu_{\mathrm{D}}$ is a fundamental quantity characterizing the medium
properties, as recently explored in Ref. Yang  _et al._ (2024).

When the gluon radiation process is included in the simulations, as shown in
Figs. 21 and 22, significantly more soft particles are produced through
medium-induced radiation as the jet evolves through the medium. Clearly, the
energy deposition (positive value of the distribution function) can be
observed in the forward direction ($\phi\approx 0$) even for the softest
energy range (0 – 1 GeV), though energy depletion can still be seen in the
backward direction ($\phi\approx\pi$). Since gluons from medium-induced
radiation are dominantly produced in the collinear direction with respect to
the hard parton, the angular distributions of partons within all energy ranges
peak at $\phi=0$. Since gluons from medium-induced radiation are dominantly
produced in the collinear direction with respect to the hard parton, the
angular distributions of partons within all energy ranges peak near $\phi=0$.
Note that the LPM suppression previously shown in Fig. 7 at very small angle
($\theta<1/\sqrt{\Delta tE}\approx 1/25$ for a 100 GeV parton traveling
through a 1 – 2 fm medium) cannot be observed in these angular distributions
due to contributions from recoil partons, multiple gluon emissions and
rescatterings they experience. Recent studies show, however, a new subjet
observable known as energy-energy correlator can resolve the angular structure
caused by the LPM interference of gluon emissions induced by multiple
scattering in the QGP medium Andres  _et al._ (2023); Yang  _et al._ (2024).
Including rescatterings of radiated gluons and recoil partons (Fig. 22)
generates much broader distributions of the soft partons than not including
these rescatterings (Fig. 21).

##  VI Jet modification in a static medium

In this section, we study properties of fully reconstructed jets originating
from a single hard jet parton in a static and uniform medium. Constituents of
these reconstructed jets include the leading jet parton, medium-induced
gluons, recoil partons, and “negative” partons whose energies are subtracted
from the jet energy within the jet cone. A jet in high-energy collisions is
defined as a cluster of collimated hadrons or partons reconstructed within a
cone,

$\sqrt{(\eta-\eta_{\mathrm{jet}})^{2}+(\phi-\phi_{\mathrm{jet}})^{2}}\leq R,$
(30)

based on a given jet finding algorithm Cacciari  _et al._ (2012) and a
particular background subtraction scheme. In the expression above, $\eta$
($\eta_{\mathrm{jet}}$) and $\phi$ ($\phi_{\mathrm{jet}}$) are the
pseudorapidity and azimuthal angle of the final state parton/hadron (jet),
respectively, and $R$ is a pre-defined jet cone size. In our present study, we
assume the effect of hadronization on the reconstructed jet energy is small
and calculate reconstructed jet energy using the final state partons in LBT
simulations. For other observables such as the jet fragmentation functions one
should include hadronization in the future.

We use a modified version of the FASTJET package Cacciari  _et al._ (2012) to
reconstruct jets with the anti-$k_{\mathrm{T}}$ algorithm, in which $\eta$ and
$\phi$ of “negative” partons are calculated in the same way as for regular
partons, while their energy-momentum is subtracted from that of the regular
partons in each iteration of the jet finding algorithm. All jet partons from
the LBT simulations, including the leading parton, gluons from medium-induced
radiations, recoil partons and “negative” partons, are fed to the FASTJET
package. This corresponds to a perfect subtraction of the unperturbed portion
of the QGP background.

Figure 23: (Color online) Accumulated fractional energy loss $\Delta E/E_{0}$
of a jet with cone size $R=0.4$ as a function of the evolution time including
only elastic scatterings with rescatterings (solid), elastic + inelastic
scatterings without rescatterings (dotted), and elastic + inelastic
scatterings with rescatterings (dot-dashed). The jet is initiated with a 100
GeV quark (upper panel) or gluon (lower panel), propagating through a static
medium at $T=300$ MeV.

###  VI.1 Jet energy loss

During the propagation of a fast parton, it loses energy and experiences
transverse momentum broadening through multiple scatterings with the medium
partons. The lost energy is carried by radiated gluons and jet-induced medium
partons which continue propagating through the medium and interacting with the
medium partons. While some fraction of the lost energy from the leading parton
will still remain inside the jet as long as its carriers are inside the jet
cone, some fraction will be transported outside the jet cone via large-angle
scatterings, large-angle emissions and successive diffusion of the jet-induced
medium partons. The latter results in a reduction of the energy of a
reconstructed jet. Shown in Fig. 23 is the accumulated energy loss of a
reconstructed jet with a cone size of $R=0.4$ as a function of the evolution
time. The jet is initiated by either a quark (upper panel) or a gluon (lower
panel) with an initial energy of 100 GeV and travels through a static medium
at a temperature of 300 MeV. In each panel, one may observe the jet energy
loss is significantly enhanced after the processes of medium-induced gluon
radiation are included and a significant fraction of the radiated gluons are
emitted outside the jet cone as the LPM interference suppresses collinear
gluon radiations. The jet energy loss also becomes larger when rescatterings
of radiated gluons and recoil partons are allowed, since these rescatterings
further transport partons from small to large angles with respect to the jet
axis. For each of the three model settings in each panel, the jet energy loss
appears to increase approximately linearly with time (or distance). This is
different from the feature of the energy loss for a single parton previously
shown in Figs. 10 and 11 due to the different energy loss mechanisms for a
single parton and a reconstructed jet. As expected, the energy loss of a gluon
jet (lower panel) is larger than that of a quark jet (upper panel). Effects of
different components of jet-induced medium excitation on the jet energy loss
were investigated in our earlier work He  _et al._ (2019), where we found
including contributions of recoil partons reduces the jet energy loss, while
including “negative” partons increases the jet energy loss.

Figure 24: (Color online) Fractional longitudinal momentum
($z_{\mathrm{jet}}=p_{\mathrm{L}}/E_{\mathrm{jet}}$) distribution of
constituent partons within a jet with a cone size $R=0.4$ at different
evolution times, with (solid) and without (dashed) rescatterings of radiated
gluons and recoil partons. The jet is initiated with a 100 GeV quark (upper
panel) or gluon (lower panel), propagating through a static medium at $T=300$
MeV with both elastic and inelastic scatterings switched on.

Figure 25: (Color online) The same as Fig. 24 except that the fractional
momentum $z_{0,\mathrm{jet}}$ is defined with respect to the initial parton
energy.

###  VI.2 Jet fragmentation function

Apart from the energy loss of reconstructed jets, one can also study how the
jet energy is redistributed among its constituents – leading jet parton,
radiated gluons, recoil and “negative” partons – in jet substructure
observables. In Fig. 24, we present the fractional longitudinal momentum
($z_{\mathrm{jet}}=p_{\mathrm{L}}/E_{\rm jet}$) distribution within a quark
jet (upper panel) and a gluon jet (lower panel). This is also known as the jet
fragmentation function, albeit in terms of partons instead of hadrons in
experimental measurements. Here, $p_{\mathrm{L}}$ denotes the longitudinal
momentum of a parton inside the jet with respect to the jet direction, and
$E_{\rm jet}$ denotes the energy of the reconstructed jet. Since the jet
develops from a single parton here, this $z_{\mathrm{jet}}$ distribution
starts with a $\delta$-function at $z_{\mathrm{jet}}=1$ before any interaction
occurs in the medium. As the jet evolves through a medium, a new peak near
$z_{\mathrm{jet}}=0$ appears and continues increasing, due to newly emitted
gluons and generated recoil partons at the thermal energy scale. Meanwhile,
the peak at $z_{\mathrm{jet}}=1$ first decreases and then disappears as the
leading parton loses its energy. A similar pattern was reported in an earlier
study on the $\gamma$-triggered jet fragmentation function in realistic heavy-
ion collisions Wang and Zhu (2013). The subtraction of the “negative” partons
from both the parton distribution and the reconstructed jet energy reduces the
soft parton distribution at small $z_{\mathrm{jet}}$. We observe that
rescattering enhances the fragmentation function at very low $z_{\rm jet}$ due
to soft partons from medium response (recoil minus “negative” partons) as well
as at large $z_{\rm jet}$ due to trigger bias. One can also plot the
coincident fragmentation function as shown in Fig. 25 in which the momentum
fraction $z_{0,{\rm jet}}=p_{\rm L}/E_{0}$ is defined relative to the initial
parton energy $E_{0}$ instead of the final jet energy $E_{\rm jet}$. The
enhancement due to rescattering at low $z_{0,{\rm jet}}$ remain approximately
the same. However, there is no enhancement at large $z_{0,{\rm jet}}$. Note
that the medium modification of the fragmentation function in terms of the
momentum fraction relative to the initiating parton energy (or the $\gamma$’s
energy in $\gamma$-triggered jets) is more sensitive to jet-medium
interactions than that in terms of the momentum fraction relative to the
reconstructed jet energy, since in the latter case, the reduced energy of the
reconstructed jet will rescale $z_{\rm jet}$ to a larger value and therefore
enhances the parton distribution at $z_{\mathrm{jet}}\approx 1$. Comparing the
two panels of Fig. 24 or Fig. 25, we observe that the population of the
distribution in the $z_{\mathrm{jet}}<1$ region is faster for a gluon jet than
for a quark jet. In each panel, we also notice that rescatterings of radiated
gluons and recoil partons accelerates the smearing of the distribution
function towards the small $z_{\mathrm{jet}}$ region.

###  VI.3 Jet shape

Figure 26: (Color online) Transverse energy distribution within a jet with a
cone size $R=0.4$ at different times during the jet propagation with (solid
lines) and without (dashed lines) rescatterings of radiated gluons and recoil
partons. The jet is initiated with a 100 GeV quark (upper panel) or gluon
(lower panel), propagating through a static medium at $T=300$ MeV with both
elastic and inelastic scatterings switched on.

In addition to the jet fragmentation function, another jet substructure
observable is the jet shape. It measures the jet energy distribution
transverse to the jet axis, defined as

$\rho(r)=\frac{1}{\Delta
r}\frac{1}{N^{\text{jet}}}\sum_{\text{jets}}\frac{E(r-\Delta r/2,r+\Delta
r/2)}{E(0,R)},$ (31)

where $r=\sqrt{(\eta-\eta_{\mathrm{jet}})^{2}+(\phi-\phi_{\mathrm{jet}})^{2}}$
is the radius to the jet axis in the $\eta$-$\phi$ plane, $E(r_{1},r_{2})$
denotes the energy carried by all partons in the jet shower within the annulus
between $r_{1}$ and $r_{2}$, $R$ is the jet cone size, and $N^{\mathrm{jet}}$
is the total number of jets being analyzed. Shown in Fig. 26 is the shape
functions of a jet at different times during the propagation and evolution in
a static medium. The jet originates from a single quark (upper panel) or a
gluon (lower panel) with initial energy $E_{0}=100$ GeV, and is reconstructed
with a cone size of $R=0.4$ at a given time. Starting as a $\delta$-function
at $r=0$, the jet shape smears towards the large $r$ region as the jet-medium
interactions transport the jet energy into large angles. With both elastic and
inelastic scatterings switched on in the LBT simulations, we observe that the
rescatterings of radiated gluons and recoil partons with the medium accelerate
the broadening of the jet. Comparing the upper and lower panel, we see a
faster broadening of a gluon jet than a quark jet. The medium modification of
the single inclusive and $\gamma$-triggered jet shapes in realistic heavy-ion
collisions were previously explored in Refs. Tachibana  _et al._ (2017) and
Luo  _et al._ (2018), where the jet-induced medium excitation is shown to
significantly enhance the jet shape at large $r$. Considering that the hard
core of the jet is less affected by rescattering than the soft partons at
larger $r$, and the jet shape is defined as a self-normalized function,
introducing rescattering may also enhance the jet shape at very small $r$ when
the intermediate $r$ region is depleted.

##  VII Nuclear modification of jets in heavy-ion collisions

In this section, we use the LBT model to study jet propagation and evolution
in an expanding QGP in realistic heavy-ion collisions following the procedure
as outlined in Sec. II and illustrated by the flow chart in Fig. 1. This is a
supplement to our earlier studies using LBT and we compare to experimental
data that only become available recently. We will especially focus on the
cone-size dependence of the jet modifications.

###  VII.1 Simulating jet quenching event-by-event

As the first step in the setup of the LBT model, we use Pythia 8 Sjostrand
_et al._ (2008) to generate the initial jet showers from nucleon-nucleon
collisions at the given center-of-mass energy within specified regions of
kinematics. To achieve enough statistics, we sometimes generate a given number
of events of initial jet production within different bins of kinematics (e.g.
transverse momentum exchange in hard scatterings) and then average the final
results over these kinematic bins with the corresponding initial jet
production cross sections as weights. The initial spatial distribution of
these jets in heavy-ion collisions is determined according to the hard
scattering locations from the AMPT Lin  _et al._ (2005) model. We also use the
AMPT model to consistently generate the initial energy density distributions
of the QGP by coarse graining the energy-momentum tensor of produced partons
in the initial stage of nuclear collisions Pang  _et al._ (2012). The
subsequent evolution of the QGP is then simulated using the CLVisc
hydrodynamic model Pang  _et al._ (2012, 2015, 2018). With the starting time
of the hydrodynamic evolution set as $\tau_{0}=0.6$ fm/$c$, the specific shear
viscosity set as $\eta/s=0.08$ and hadronic freeze-out temperature $T_{f}=137$
MeV, the hydrodynamic calculation is able to provide a good description of the
soft hadron spectra and their anisotropic flow coefficients at both RHIC and
LHC. Using the local temperature and flow velocity information given by these
hydrodynamic simulations at each time step, we first boost a given hard parton
in the jet shower into the local rest frame of the expanding QGP, in which the
linear Boltzmann equation [Eq. (1)] is solved by simulating the parton
scattering with its surrounding medium. The final state jet partons, together
with its companions (radiated gluons, recoil and “negative” partons) if
scatterings happen during the time interval, are then boosted back to the
global frame, where they stream freely to the locations of the next time step.
This process is iterated for all hard partons until they travel outside the
QGP regime, i.e., to locations where the local temperature drops below the
phase transition temperature $T_{\mathrm{c}}$ (taken as 160 MeV in this work).
We neglect the jet-medium interaction in the hadronic phase where the jet
transport coefficient becomes negligibly small Chen  _et al._ (2010). To
incorporate effects of event-by-event fluctuations, we prepare 200
hydrodynamic events for each collision setup (each centrality for each
colliding system), and simulate a large number of events of jet-medium
interactions through each hydrodynamic profile until the statistical errors of
the extracted observables are sufficiently small. This LBT model has been
successfully applied to study the nuclear modification of single inclusive
hadrons Cao  _et al._ (2018), jets He  _et al._ (2019, 2022) and photon-
triggered jets Wang and Zhu (2013); Luo  _et al._ (2018). As a complementary
study to these earlier work, we present additional results from LBT and
compare to new experimental data in this section and discuss how jet-medium
interactions affect different jet observables in heavy-ion collisions.

In some of the published studies with LBT and in the following subsections, we
have only considered the first term in Eq. (13) for the formation time of the
initial jet shower partons from Pythia simulations. Values of the effective
strong coupling constant $\alpha_{\rm s}$ from fitting to experimental data
are somewhat small, ranging from 0.15 to 0.2, depending on the colliding
energy. If we include the second term in the formation time for the initial
jet shower partons, we can still fit the experimental data by readjusting the
value of $\alpha_{\rm s}$. Since this will increase the formation time for
some of the initial jet shower partons and delay the parton-medium interaction
to later times, the effective value of $\alpha_{\rm s}$ will be larger, from
$\alpha_{\rm s}=0.2$ to 0.3, in order to describe the data. We should
emphasize that this is only an effective and constant coupling constant in LBT
model which should be considered as an adjustable parameter. In the most
recent extension of the LBT model, the quasi-particle LBT (QLBT) model Liu
_et al._ (2022b), one assumes that the QGP consists of thermal quasi-particles
with an effective Debye screening mass that depends on a temperature-dependent
strong coupling $g(T)$ that also controls the interaction between exchanged
gluons and medium partons. The effective strong coupling constant between jet
shower partons and the medium through Bayesian analyses of experimental data
on heavy flavor mesons ranges from 0.2 to 0.3, depending on the temperature
and parton energy. While readjusting $\alpha_{\mathrm{s}}$ can fit the jet
$R_{\mathrm{AA}}$ data for different applications of the parton formation
time, deviations may remain in more detailed substructures of jets that are
more sensitive to the evolution of their soft components. Thus, we utilize the
full Eq. (13) in our current and future calculations.

Figure 27: (Color online) LBT simulations of the nuclear modification factor
$R_{\mathrm{AA}}$ for single inclusive jets in 0-10% Au+Au collisions at
$\sqrt{s_{\mathrm{NN}}}=200$ GeV with jet cone sizes $R=0.2$ (top panel), 0.3
(middle panel) and 0.4 (bottom panel) as compared to STAR data Adam  _et al._
(2020).
Figure 28: (Color online) The ratio of the inclusive jet spectra with jet cone
size $R=0.2$ and 0.4 in 0-10% Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=200$
GeV from the LBT calculations as compared to the STAR data Adam  _et al._
(2020) and the corresponding $p+p$ result from Pythia simulation. The black
dashed line, red solid line and green solid line represent LBT results with
minimal $p_{\mathrm{T}}$ trigger bias as 0, 5 GeV/c and 7 GeV/c for the
leading parton respectively. The purple solid line represents the average of
LBT results between minimal $p_{\mathrm{T}}$ trigger bias of 5 GeV/c and 7
GeV/c for the leading parton, while the purple dashed line represents the
Pythia result for $p+p$ collisions.

###  VII.2 Cone-size dependence of single inclusive jet suppression

The most widely used observable that quantifies the QGP effect on jets is
their nuclear modification factor, as first proposed in Ref. Wang (1998),

$R_{\mathrm{AA}}=\frac{d\sigma^{\rm AA}_{\rm jet}}{\langle N_{\rm bin}\rangle
d\sigma^{pp}_{\rm jet}},$ (32)

where $\langle N_{\rm bin}\rangle$ is the number of binary collisions in the
given centrality bin of A+A collisions. The LBT model was shown to describe
the centrality and $p_{\rm T}$ dependence of $R_{\mathrm{AA}}$ at the LHC
energies and the effect of jet-induced medium response is very important He
_et al._ (2019). Shown in Fig. 27 is the jet $R_{\mathrm{AA}}$ in 0-10%
central Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=200$ GeV as compared to
the recent STAR experimental data Adam  _et al._ (2020). From the top to
bottom panel, we present results for jets with cone sizes of $R=0.2$, 0.3 and
0.4. In the STAR measurement, a high $p_{\mathrm{T}}$ hadron (above 5 GeV)
inside the jet is required to suppress the background. Since the hadronization
process has not been implemented in the LBT calculations presented in this
study, we select jet events with energetic partons and vary the
$p_{\mathrm{T}}$ cut to study its effect on the jet $R_{\mathrm{AA}}$. As
shown in the figure, increasing the $p_{\mathrm{T}}$ cut from 0 (denoted as
“without trigger bias”) to 5 and 7 GeV/$c$ excludes certain jet events at low
$p_{\mathrm{T}}$. Since a large number of jet events at low $p_{\mathrm{T}}$
are generated in A+A collisions after jet-medium interactions, effect of the
trigger bias is more significant in A+A than in $p+p$ collisions. This leads
to a smaller jet $R_{\mathrm{AA}}$ at low $p_{\mathrm{T}}$ when the trigger
$p_{\mathrm{T}}$ cut is increased. However, no visible effect is observed on
the jet $R_{\mathrm{AA}}$ above $p_{\rm T}=20$ GeV/$c$. With a strong coupling
constant set as $\alpha_{\mathrm{s}}=0.2$, the LBT results (shown by the pink
band between 5 and 7 GeV/$c$ cuts for the leading hadron $p_{\mathrm{T}}$) are
consistent with the experimental data. Since the leading hadron usually comes
from the fragmentation of a higher $p_{\mathrm{T}}$ parton, a slightly higher
$p_{\mathrm{T}}$ cut for partons is expected compared to that used for
hadrons.

Shown in Fig. 28 is the ratio between the single inclusive jet spectra with
cone size $R=0.2$ and 0.4 at the RHIC energy. As one uses a larger jet cone,
more particles are included in the jet reconstruction and thus the yield of
jets increases. Our LBT model result with a lower $p_{\mathrm{T}}$ cut on the
triggered parton between 5 and 7 GeV/$c$ agrees with the STAR data. Comparing
$p+p$ and Au+Au collisions, we observe this ratio for the former is slightly
higher than the latter at high $p_{\mathrm{T}}$, while the opposite is seen at
low $p_{\mathrm{T}}$. This cone-size dependence of the jet spectrum in $p+p$
collisions affects the cone-size dependence of the jet $R_{\mathrm{AA}}$ shown
in Fig. 27.

Figure 29: (Color online) LBT results on the nuclear modification factor
$R_{\mathrm{AA}}$ for single inclusive jets in 0-10% Pb+Pb collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02$ TeV, with different jet cone sizes and two
different values of $\alpha_{\mathrm{s}}=$ 0.13 (red) and 0.15 (blue) as
compared to the CMS Sirunyan  _et al._ (2021) and ATLAS Aaboud  _et al._
(2019) data.

CMS experiment Sirunyan  _et al._ (2021) has recently measured $R_{\rm AA}$ of
large $p_{\rm T}$ single inclusive jets at LHC with very large jet cone size
and found a weak jet cone-size dependence of $R_{\rm AA}$ which was not
reproduced by many models, including LBT. The jet $R_{\mathrm{AA}}$ as a
function of $p_{\mathrm{T}}$ from LBT with different jet cone-sizes is shown
in Fig. 29 together with the CMS data for central Pb+Pb collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02$ TeV. The LBT results at LHC with $\alpha_{\rm
s}=0.15$ (blue dashed lines) are tuned to reproduce the ATLAS data Aad  _et
al._ (2015); Aaboud  _et al._ (2019) on single inclusive jet suppression at
both $\sqrt{s_{\mathrm{NN}}}=2.76$ and 5.02 TeV He  _et al._ (2019). This
value of $\alpha_{\mathrm{s}}$ at LHC is smaller than the value we used to fit
the STAR data at RHIC. A possible cause of this colliding energy dependence is
the lower average temperature of the QGP and thus stronger jet-medium
interactions at RHIC than LHC. With this value of $\alpha_{\rm s}$, LBT
however fails to reproduce the CMS data for $R\lesssim 0.4$.

In Fig. 29, we note that the jet $R_{\mathrm{AA}}$ from LBT calculations
increases as the jet cone becomes larger. This is a general feature due to the
jet-induced medium excitation: when more recoil partons scattered into large
angles are included in jet reconstruction as one increases the cone size, the
jet $R_{\mathrm{AA}}$ should increase. On the other hand, such cone-size
dependence becomes much weaker if one only considers the medium-induced gluon
emission process which prefers small angle (shown in Fig. 8). Similar
conclusions have also been drawn by other model calculations, e.g. JEWEL Zapp
_et al._ (2011, 2013); Zapp (2014); Kunnawalkam Elayavalli and Zapp (2017) and
jet-fluid Tachibana  _et al._ (2017); Chang  _et al._ (2020), that take into
account of the effects of medium response induced by jet-medium interaction.

We note that with the jet cone size $R=0.4$, where measurements are available
from both CMS and ATLAS experiment, tension exists between CMS and ATLAS data:
the jet $R_{\mathrm{AA}}$ observed by CMS seems larger than ATLAS. If we
adjust the strong coupling constant to $\alpha_{\mathrm{s}}=0.13$ in LBT to
fit CMS data with $R=0.4$ instead, the agreement between LBT and CMS data
improves for $R<0.4$. But LBT results are still about one sigma blow the CMS
data for small jet cone size $R=0.2$. It is clear that resolving the tension
between ALTAS and CMS data is necessary to draw a more solid conclusion on the
cone size dependence of the nuclear modification on jets.

Figure 30: (Color online) The nuclear modification factor $I_{\mathrm{AA}}$
from LBT simulations for $\gamma$-triggered jets in 0-15% Au+Au collisions at
$\sqrt{s_{\mathrm{NN}}}=200$ GeV with different jet cone sizes as compared to
the STAR data Anderson (2023).

###  VII.3 $\gamma$-jets and dijets

Electro-weak-boson-triggered jets are considered a “golden channel” to
quantify the nuclear modification of jets, since these bosons do not
participate in strong interactions with the QGP and thus serve as a reliable
reference of the initial energy of its associated jets. Suppression and
modification of $\gamma/Z$-triggered jets in Pb+Pb collisions at the LHC
energies can be described well by the LBT model Luo  _et al._ (2018); Zhang
_et al._ (2018). In Fig. 30, we study the nuclear modification factor of jets
triggered by photons within $15<E_{\mathrm{T}}^{\gamma}<20$ GeV
($I_{\mathrm{AA}}$) in central Au+Au collisions at RHIC. Two cone sizes,
$R=0.2$ and 0.5, are used for jet reconstruction. Both the STAR experimental
data Anderson (2023). and our LBT calculations indicate that a larger cone
size captures more energy in the reconstructed jets, thus results in a higher
$I_{\mathrm{AA}}$. This is consistent with our earlier findings for the
suppression factor $R_{\mathrm{AA}}$ for single inclusive jets in LBT
simulations.

Figure 31: (Color online) Dijet asymmetry ($x_{\mathrm{J}}$) distribution in
$p+p$ (blue) and Pb+Pb collisions (red) at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV
from LBT calculations as compared to the ATLAS data Aad  _et al._ (2023). The
transverse momentum of the leading jet is selected within
$158<p_{\mathrm{T,1}}<178$ GeV/$c$ in the upper panel and
$398<p_{\mathrm{T,1}}<562$ GeV/$c$ in the lower panel.

Similar to photon-triggered jets, one may also quantify the jet energy loss by
comparing the transverse momenta of the subleading and leading jet in a dijet
event Qin and Muller (2011). In Fig. 31, we present the distribution of
subleading-to-leading-jet-$p_{\mathrm{T}}$ ratio $(x_{\mathrm{J}}=E_{T}^{\rm
sublead}/E_{T}^{\rm lead})$ of dijets produced at the LHC energy, compared
between $p+p$ and Pb+Pb collisions. One observes that while the
$x_{\mathrm{J}}$ distribution peaks around one in $p+p$ collisions, it shifts
to lower $x_{\mathrm{J}}$ values in Pb+Pb collisions due to the unbalanced
energy loss between subleading and leading jets as they travel through the
QGP. As one increases the $p_{\mathrm{T}}$ range of the leading jet (from the
upper to the lower panel), dijets appear more symmetric in both $p+p$ and
Pb+Pb collisions. Our Pythia simulations for $p+p$ collisions and Pythia+LBT
simulation for Pb+Pb collisions agree well with the corresponding data from
the ATLAS collaboration.

##  VIII Code availability

The source code of the LBT model can be downloaded from
“https://github.com/lbt-jet”, which includes the main modules, the data tables
for the elastic and inelastic scattering rate, a sample data file of the
hydrodynamic profile and the modified version of the FASTJET package that can
take into account of the “negative” partons.

The LBT code is also embedded as a jet energy loss module in the JETSCAPE
framework, which can be downloaded at “https://github.com/JETSCAPE”.

##  IX Summary

We have developed a linear Boltzmann transport (LBT) model for studying jet-
medium interaction in relativistic heavy-ion collisions. This article provides
a detailed description of the model structure, physics processes, Monte Carlo
implementation and numerical validation.

Both elastic and inelastic scatterings between jet partons and the QGP are
implemented in LBT. The former is based on the leading-order pQCD scatterings
between a jet shower parton and a medium parton at the thermal scale, while
the latter is based on the medium-induced gluon emission process within the
higher-twist energy loss formalism. With these scattering processes, we track
the phase space evolution of not only leading partons and their emitted
gluons, but also the thermal partons being scattered out of the QGP background
known as recoil partons and the energy-momentum holes left inside the medium
(“negative” partons). This guarantees the energy-momentum conservation in each
scattering and allows us to study both medium modification of jets and jet-
induced medium excitation. All partons from the LBT simulations are fed into
the FASTJET package for jet finding in the end for reconstructing partonic
jets. The FASTJET routine is modified such that the energy-momentum of
“negative” partons is subtracted from that of regular ones throughout the jet
finding procedure. Within this framework, we present a systematical study on
the evolution of energy and transverse momentum distributions of single hard
partons, jet-induced medium partons, and fully reconstructed jets inside a
static and uniform medium at a fixed temperature.

We show that the accumulated energy loss of a single hard parton increases
linearly with time (or path length) in elastic scatterings, but increases
quadratically at early times after the gluon emission process is included, as
they should according to the formalism implemented in LBT. As the parton
energy decreases along its path, its energy loss rate (per unit length) first
increases with time, then saturates and decreases in our simulations, with its
inflection point becoming earlier in a hotter medium and for a lower energy
parton. These are all the consequences of the energy dependence of the parton
energy loss rate. A broadening in the transverse momentum distributions is
observed for the hard partons, with the broadening being accelerated after the
radiation process is taken into account. By incorporating all partons into the
analysis, including leading partons, gluons from medium-induced radiation and
jet-induced medium excitation (recoil and “negative” partons), a Mach-cone-
like structure can be clearly observed in the energy density distribution of
these partons with an enhanced wave front when rescatterings between recoil
partons and radiated gluons with the medium are allowed. Propagation of
radiated gluons and recoil partons and their further interaction with the
medium convert the energy deposition from jets into the wave front followed by
a diffusion wake in the direction opposite to the jet propagation due to the
energy depletion via “negative” partons. Compared to the pure elastic
scattering processes, the Mach-cone-like structure is significantly enhanced
when medium-induced gluon emission is introduced. By analyzing the parton
spectra within different energy ranges, one can observe the transport of
energy from the hard parton at small angle to soft partons at large angle with
respect to the initial parton. This energy flow results in the energy loss of
fully reconstructed jets, which is found to increase approximately linearly
with time. The energy flow also causes the shift of the jet fragmentation
function towards the smaller momentum fraction region and the enhancement of
the jet shape at large radius. Rescatterings of recoil partons and radiated
gluons with the medium are found essential in the nuclear modification of
jets. They significantly broaden the shock wave of the parton energy density
distribution, and accelerate the energy transport from hard to soft partons
and from small to large angles. Similar features of nuclear modification are
observed for quark jets and gluon jets, except that the modification of the
latter is stronger due to its larger color factor in jet-medium interaction.
This documentation of the physics implementations in the LBT model, in
particular inelastic processes, together with a thorough validation and a
detailed study of the contributions from different model components to jet
energy loss, momentum broadening and its induced medium excitation, will help
to advance our understanding of jet evolution inside a dense QCD medium.

To simulate jet-medium interactions in realistic heavy-ion collisions, the LBT
model is further combined with the Pythia 8 generator for the initial jet
production in nucleon-nucleon collisions, the AMPT model for the spatial
distribution of the hard scattering vertices in nucleus-nucleus collisions,
and the CLVisc hydrodynamic model for the space-time evolution profile of the
QGP with the initial conditions from the same AMPT model calculations.
Complimentary to jet observables presented in our earlier publications, we
calculate the jet cone-size dependence of the jet suppression factor in
different collision systems at different beam energies, the $\gamma$-jet
suppression factor and dijet asymmetry within different kinematic ranges in
the present study, and compare to more recent experimental data at RHIC and
LHC.

###### Acknowledgements.

We thank Wei Chen, Zhong Yang and Weiyao Ke for helpful discussions. This work
is supported in part by National Natural Science Foundation of China (NSFC)
under Grant Nos. 11935007, 11221504, 11861131009, 11890714, 12075098,
12175122, 2021-867 and 12147134, by EU ERDF and H2020 grant 82409, ERC grant
ERC-2018-ADG-835105, Spanish AEI grant PID2020-119632GB-I00 and
CEX2020-001035-M, Xunta de Galicia Research Center accreditation 2019-2022, by
Guangdong Major Project of Basic and Applied Basic Research No.
2020B0301030008, by Guangdong Basic and Applied Basic Research Foundation No.
2021A1515110817, by Science and Technology Program of Guangzhou No.
2019050001, by the Director, Office of Energy Research, Office of High Energy
and Nuclear Physics, Division of Nuclear Physics, of the U.S. Department of
Energy under Contract No. DE-AC02- 05CH11231, by the US National Science
Foundation under Grant No. OAC-2004571 within the X-SCAPE Collaboration.
Computations in this study are performed at the NSC3/CCNU and the National
Energy Research Scientific Computing Center (NERSC), a U.S. Department of
Energy Office of Science User Facility located at Lawrence Berkeley National
Laboratory and operated under Contract No. DE-AC02-05CH11231.

##  Appendix A Energy-momentum conservation of the $2\rightarrow 2+n$ process

In this Appendix, we explain the procedure to guarantee energy-momentum
conservation for a general $2\rightarrow 2+n$ process in LBT.

Figure 32: Illustration of a $2\rightarrow 3$ process.

For a $2\rightarrow 3$ process, as illustrated in Fig. 32, we first sample a
$2\rightarrow 2$ process according to Eq. (2) which by itself respects the
energy-momentum conservation. Then we sample a medium-induced gluon “5” from
parton “3” in Fig. 32 based on Eq.(7). To maintain the energy-momentum
conservation of the entire $2\rightarrow 3$ process, we choose to adjust the
energy-momentum of partons “3” and “4”.

For convenience, we assume parton “1” travels along the $z$-direction with
four-momentum $(E_{1},p_{1x}=0,p_{1y}=0,p_{1z})$, and fix the momentum of
parton “2” from the $2\rightarrow 2$ process as
$(E_{2},p_{2x},p_{2y},p_{2z})$, and the momentum of “5” from the medium-
induced gluon spectrum as $(k_{0},k_{x},k_{y},k_{z})$. The momentum of the
exchanged gluon between “1” and “2”, $(q_{0},q_{x},q_{y},q_{z})$, can also be
calculated from the $2\rightarrow 2$ process. However, since the direction of
this exchanged gluon can be correlated with that of the emitted gluon, we only
keep the magnitude $q_{\perp}=\sqrt{q_{x}^{2}+q_{y}^{2}}$ and re-sample its
direction in the transverse plane with respect to parton “1”. Components
$q_{0}$ and $q_{z}$ will be re-evaluated to recover the energy-momentum
conservation of this $2\rightarrow 3$ process. With this setup, the momenta of
“3” and “4” are then

$\displaystyle p_{3}=\,$
$\displaystyle(E_{1}-q_{0}-k_{0},\,p_{1x}-q_{x}-k_{x},\,p_{1y}-q_{y}-k_{y},$
$\displaystyle p_{1z}-q_{z}-k_{z}),$ (33) $\displaystyle p_{4}=\,$
$\displaystyle(E_{2}+q_{0},\,p_{2x}+q_{x},\,p_{2y}+q_{y},\,p_{2z}+q_{z}),$
(34)

respectively. The two on-shell conditions for “3” and “4”,

$\displaystyle p_{3}^{3}=m^{2},$ (35) $\displaystyle p_{4}^{2}=0,$ (36)

will be solved to determine $q_{0}$ and $q_{z}$. We require that the exchanged
gluon is space-like $(q^{2}<0)$ for the solutions of $q_{0}$ and $q_{z}$. If
there is no solution for Eqs. (35) and (36), or their solutions do not satisfy
the space-like condition, we re-sample the direction of $q_{\perp}$ or the
four-momentum of the emitted gluon “5”. This process is iterated until an
appropriate solution is found or the number of trials exceeds a given large
value. For the latter case, we consider the gluon emission process is not
kinematically allowed by the associated $2\rightarrow 2$ process and reject
its formation. Note that there could be two solutions from Eqs. (35) and (36)
which both satisfy the space-like condition of $q$. In this situation, we
currently select the solution with a smaller $q_{0}$, assuming that the energy
exchange between “1” and “2” is small when they exchange momentum. One may
verify that the same procedure here is also applicable for the scenario where
the momentum exchange ($q$) occurs between the emitted gluon “5” and the
thermal parton “2”.

Figure 33: Illustration of multiple gluon emissions.

For a multiple gluon emission process, we first determine a $2\rightarrow 3$
process as described above and fix the energy-momentum of the recoil parton
“4” in Fig. 32. As illustrated by Fig. 33, if one more gluon “6” is emitted,
we need to adjust the energy-momentum of the previous gluon “5” and the final
state jet parton “3” to respect the energy-momentum conservation of this
$1\rightarrow 3$ splitting process. This is sufficient to guarantee the
energy-momentum conservation of the entire $2\rightarrow 4$ process.

We denote the energy-momentum of “$1^{\prime}$” (state of the jet parton
before emitting the previous gluon “5”) as
($E_{0},\vec{p}_{0\perp}=0,p_{0z}$). This is determined from the $2\rightarrow
3$ process as described earlier and rotated into the coordinate where
“$1^{\prime}$” travels along the $z$-direction for convenience. The energy-
momentum of the new gluon “6” is sampled using Eq. (7) from “3” (before
adjusting its momentum here) as $(k_{6e},\vec{k}_{6\perp},k_{6z})$. The
energy-momentum of the previous gluon is represented by
$(k_{5e},\vec{k}_{5\perp},k_{5z})$. Although it is already given by the
$2\rightarrow 3$ process, we choose to only keep the value of $k_{5e}$ but re-
calculate the $k_{5z}$ component. The transverse component is then re-
evaluated as $k_{5\perp}=\sqrt{k_{5e}^{2}-k_{5z}^{2}}$ given $k_{5z}$ and its
direction is sampled randomly with respect to $\vec{k}_{6\perp}$.

Within this setup, the energy-momentum conservation of the $1\rightarrow 3$
splitting process is achieved using the on-shell condition for the final state
jet parton “3”,

$\displaystyle p_{3}$ $\displaystyle=$
$\displaystyle(E_{0}-k_{5e}-k_{6e},\,-\vec{k}_{5\perp}-\vec{k}_{6\perp},p_{0z}-k_{5z}-k_{6z}),$
$\displaystyle p_{3}^{2}$ $\displaystyle=$ $\displaystyle m^{2}$ (37)

which can be further expanded as

$\displaystyle(E_{0}-k_{5e}-k_{6e})^{2}-(k_{5\perp}^{2}+k_{6\perp}^{2}+2k_{5\perp}k_{6\perp}\cos\theta)$
$\displaystyle-(p_{0z}-k_{5z}-k_{6z})^{2}=m^{2}.$ (38)

Using a randomly sampled $\theta$, we solve the above equation for $k_{5z}$.
Two conditions are examined for the solution: (1) $k_{5e}>|k_{5z}|$ and (2)
the intermediate jet parton between the two emitted gluons
(“$1^{\prime\prime}$” in Fig. 33) is time-like: $p_{1}^{\prime\prime
2}>m^{2}$. If both solutions of Eq. (A) satisfy these conditions, we choose
the solution $k_{5z}$ that is closer to its original value obtained from the
previous $2\rightarrow 3$ process. If there is no appropriate solution, we
iterate the procedure above by re-sampling the angle $\theta$ or the four-
momentum of “6” until the solution is found or the number of trials exceeds a
preset large value.

If more gluons are emitted, we repeat the procedure above to sample the four-
momentum of the last gluon from Eq. (7) and adjust that of the second to last
gluon to ensure the energy-momentum conservation of the entire $2\rightarrow
2+n$ process.

## References

  * Shuryak (2017) E. Shuryak, Rev. Mod. Phys. 89, 035001 (2017), arXiv:1412.8393 [hep-ph] .
  * Adams  _et al._ (2005) J. Adams _et al._ (STAR), Nucl. Phys. A 757, 102 (2005), arXiv:nucl-ex/0501009 .
  * Adcox  _et al._ (2005) K. Adcox _et al._ (PHENIX), Nucl. Phys. A 757, 184 (2005), arXiv:nucl-ex/0410003 .
  * Muller  _et al._ (2012) B. Muller, J. Schukraft, and B. Wyslouch, Ann. Rev. Nucl. Part. Sci. 62, 361 (2012), arXiv:1202.3233 [hep-ex] .
  * Akiba  _et al._ (2015) Y. Akiba _et al._ , (2015), arXiv:1502.02730 [nucl-ex] .
  * Mehtar-Tani  _et al._ (2012) Y. Mehtar-Tani, C. A. Salgado, and K. Tywoniuk, Phys. Lett. B 707, 156 (2012), arXiv:1102.4317 [hep-ph] .
  * Armesto  _et al._ (2012) N. Armesto, H. Ma, Y. Mehtar-Tani, C. A. Salgado, and K. Tywoniuk, JHEP 01, 109 (2012), arXiv:1110.4343 [hep-ph] .
  * Casalderrey-Solana and Iancu (2011) J. Casalderrey-Solana and E. Iancu, JHEP 08, 015 (2011), arXiv:1105.1760 [hep-ph] .
  * Kumar  _et al._ (2020) A. Kumar, A. Majumder, and C. Shen, Phys. Rev. C 101, 034908 (2020), arXiv:1909.03178 [nucl-th] .
  * Connors  _et al._ (2018) M. Connors, C. Nattrass, R. Reed, and S. Salur, Rev. Mod. Phys. 90, 025005 (2018), arXiv:1705.01974 [nucl-ex] .
  * Wang and Gyulassy (1992) X.-N. Wang and M. Gyulassy, Phys. Rev. Lett. 68, 1480 (1992).
  * Majumder and Van Leeuwen (2011) A. Majumder and M. Van Leeuwen, Prog. Part. Nucl. Phys. 66, 41 (2011), arXiv:1002.2206 [hep-ph] .
  * Qin and Wang (2015) G.-Y. Qin and X.-N. Wang, Int. J. Mod. Phys. E24, 1530014 (2015), arXiv:1511.00790 [hep-ph] .
  * Blaizot and Mehtar-Tani (2015) J.-P. Blaizot and Y. Mehtar-Tani, Int. J. Mod. Phys. E 24, 1530012 (2015), arXiv:1503.05958 [hep-ph] .
  * Cao and Wang (2021) S. Cao and X.-N. Wang, Rept. Prog. Phys. 84, 024301 (2021), arXiv:2002.04028 [hep-ph] .
  * Cao and Qin (2023) S. Cao and G.-Y. Qin, Ann. Rev. Nucl. Part. Sci. 73, 205 (2023), arXiv:2211.16821 [nucl-th] .
  * Burke  _et al._ (2014) K. M. Burke _et al._ (JET), Phys. Rev. C 90, 014909 (2014), arXiv:1312.5003 [nucl-th] .
  * Cao  _et al._ (2021) S. Cao _et al._ (JETSCAPE), Phys. Rev. C 104, 024905 (2021), arXiv:2102.11337 [nucl-th] .
  * Xie  _et al._ (2022) M. Xie, W. Ke, H. Zhang, and X.-N. Wang, arXiv:2206.01340 (2022), arXiv:2206.01340 [hep-ph] .
  * Park  _et al._ (2017) C. Park, C. Shen, S. Jeon, and C. Gale, Nucl. Part. Phys. Proc. 289-290, 289 (2017), arXiv:1612.06754 [nucl-th] .
  * Zhao  _et al._ (2020) W. Zhao, C. M. Ko, Y.-X. Liu, G.-Y. Qin, and H. Song, Phys. Rev. Lett. 125, 072301 (2020), arXiv:1911.00826 [nucl-th] .
  * Liu  _et al._ (2022a) Y.-F. Liu, W.-J. Xing, X.-Y. Wu, G.-Y. Qin, S. Cao, and H. Xing, Phys. Rev. C 105, 044904 (2022a), arXiv:2107.01522 [hep-ph] .
  * Li  _et al._ (2021) S.-Q. Li, W.-J. Xing, X.-Y. Wu, S. Cao, and G.-Y. Qin, Eur. Phys. J. C 81, 1035 (2021), arXiv:2108.06648 [hep-ph] .
  * Adamczyk  _et al._ (2018) L. Adamczyk _et al._ (STAR), Phys. Rev. Lett. 121, 032301 (2018), arXiv:1707.01988 [nucl-ex] .
  * Wu  _et al._ (2022) J. Wu, S. Cao, and F. Li, arXiv:2208.14297 (2022), arXiv:2208.14297 [nucl-th] .
  * Bjorken (1982) J. D. Bjorken, FERMILAB-PUB-82-059-THY (1982).
  * Braaten and Thoma (1991) E. Braaten and M. H. Thoma, Phys. Rev. D44, 2625 (1991).
  * Djordjevic (2006) M. Djordjevic, Phys. Rev. C 74, 064907 (2006), arXiv:nucl-th/0603066 .
  * Qin  _et al._ (2008) G.-Y. Qin _et al._ , Phys. Rev. Lett. 100, 072301 (2008), arXiv:0710.0605 [hep-ph] .
  * Baier  _et al._ (1997) R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne, and D. Schiff, Nucl. Phys. B483, 291 (1997), hep-ph/9607355 .
  * Zakharov (1996) B. G. Zakharov, JETP Lett. 63, 952 (1996), arXiv:hep-ph/9607440 .
  * Gyulassy  _et al._ (2000) M. Gyulassy, P. Levai, and I. Vitev, Nucl. Phys. B571, 197 (2000), arXiv:hep-ph/9907461 .
  * Wiedemann (2000) U. A. Wiedemann, Nucl. Phys. B588, 303 (2000), arXiv:hep-ph/0005129 .
  * Arnold  _et al._ (2002) P. Arnold, G. D. Moore, and L. G. Yaffe, JHEP 06, 030 (2002), hep-ph/0204343 .
  * Wang and Guo (2001) X.-N. Wang and X.-F. Guo, Nucl. Phys. A696, 788 (2001), arXiv:hep-ph/0102230 .
  * Chen  _et al._ (2018) W. Chen, S. Cao, T. Luo, L.-G. Pang, and X.-N. Wang, Phys. Lett. B 777, 86 (2018), arXiv:1704.03648 [nucl-th] .
  * Chen  _et al._ (2020) W. Chen, S. Cao, T. Luo, L.-G. Pang, and X.-N. Wang, Phys. Lett. B 810, 135783 (2020), arXiv:2005.09678 [hep-ph] .
  * Casalderrey-Solana  _et al._ (2017) J. Casalderrey-Solana, D. Gulhan, G. Milhano, D. Pablos, and K. Rajagopal, JHEP 03, 135 (2017), arXiv:1609.05842 [hep-ph] .
  * Tachibana  _et al._ (2017) Y. Tachibana, N.-B. Chang, and G.-Y. Qin, Phys. Rev. C95, 044909 (2017), arXiv:1701.07951 [nucl-th] .
  * Kunnawalkam Elayavalli and Zapp (2017) R. Kunnawalkam Elayavalli and K. C. Zapp, JHEP 07, 141 (2017), arXiv:1707.01539 [hep-ph] .
  * Luo  _et al._ (2018) T. Luo, S. Cao, Y. He, and X.-N. Wang, Phys. Lett. B782, 707 (2018), arXiv:1803.06785 [hep-ph] .
  * Park  _et al._ (2019) C. Park, S. Jeon, and C. Gale,  _Proceedings, 27th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions (Quark Matter 2018): Venice, Italy, May 14-19, 2018_ , Nucl. Phys. A982, 643 (2019), arXiv:1807.06550 [nucl-th] .
  * Chatrchyan  _et al._ (2011) S. Chatrchyan _et al._ (CMS), Phys. Rev. C 84, 024906 (2011), arXiv:1102.1957 [nucl-ex] .
  * Yang  _et al._ (2021) Z. Yang, W. Chen, Y. He, W. Ke, L. Pang, and X.-N. Wang, Phys. Rev. Lett. 127, 082301 (2021), arXiv:2101.05422 [hep-ph] .
  * Yang  _et al._ (2022) Z. Yang, T. Luo, W. Chen, L.-G. Pang, and X.-N. Wang, arXiv:2203.03683 (2022), arXiv:2203.03683 [hep-ph] .
  * Chen  _et al._ (2021) W. Chen, S. Cao, T. Luo, L.-G. Pang, and X.-N. Wang, Nucl. Phys. A 1005, 121934 (2021).
  * Luo  _et al._ (2023) A. Luo, Y.-X. Mao, G.-Y. Qin, E.-K. Wang, and H.-Z. Zhang, Phys. Lett. B 837, 137638 (2023), arXiv:2109.14314 [hep-ph] .
  * Sirimanna  _et al._ (2022) C. Sirimanna, I. Soudi, G. Vujanovic, W.-J. Xing, S. Cao, and A. Majumder, arXiv:2211.15553 (2022), arXiv:2211.15553 [hep-ph] .
  * Armesto  _et al._ (2009) N. Armesto, L. Cunqueiro, and C. A. Salgado, Eur. Phys. J. C63, 679 (2009), arXiv:0907.1014 [hep-ph] .
  * Zapp  _et al._ (2011) K. C. Zapp, J. Stachel, and U. A. Wiedemann, JHEP 07, 118 (2011), arXiv:1103.6252 [hep-ph] .
  * Casalderrey-Solana  _et al._ (2014) J. Casalderrey-Solana, D. C. Gulhan, J. G. Milhano, D. Pablos, and K. Rajagopal, JHEP 10, 019 (2014), [Erratum: JHEP09,175(2015)], arXiv:1405.3864 [hep-ph] .
  * Cao and Majumder (2020) S. Cao and A. Majumder, Phys. Rev. C 101, 024903 (2020), arXiv:1712.10055 [nucl-th] .
  * Schenke  _et al._ (2009) B. Schenke, C. Gale, and S. Jeon, Phys. Rev. C80, 054913 (2009), arXiv:0909.2037 [hep-ph] .
  * Wang and Zhu (2013) X.-N. Wang and Y. Zhu, Phys. Rev. Lett. 111, 062301 (2013), arXiv:1302.5874 [hep-ph] .
  * Cao  _et al._ (2017) S. Cao _et al._ (JETSCAPE), Phys. Rev. C 96, 024909 (2017), arXiv:1705.00050 [nucl-th] .
  * Kumar  _et al._ (2022) A. Kumar _et al._ (JETSCAPE), arXiv:2204.01163 (2022), arXiv:2204.01163 [hep-ph] .
  * He  _et al._ (2022) Y. He, W. Chen, T. Luo, S. Cao, L.-G. Pang, and X.-N. Wang, Phys. Rev. C 106, 044904 (2022), arXiv:2201.08408 [hep-ph] .
  * Tachibana  _et al._ (2022) Y. Tachibana, C. Shen, and A. Majumder, Phys. Rev. C 106, L021902 (2022), arXiv:2001.08321 [nucl-th] .
  * He  _et al._ (2015) Y. He, T. Luo, X.-N. Wang, and Y. Zhu, Phys. Rev. C91, 054908 (2015), arXiv:1503.03313 [nucl-th] .
  * Cao  _et al._ (2018) S. Cao, T. Luo, G.-Y. Qin, and X.-N. Wang, Phys. Lett. B777, 255 (2018), arXiv:1703.00822 [nucl-th] .
  * Xing  _et al._ (2020) W.-J. Xing, S. Cao, G.-Y. Qin, and H. Xing, Phys. Lett. B 805, 135424 (2020), arXiv:1906.00413 [hep-ph] .
  * He  _et al._ (2019) Y. He, S. Cao, W. Chen, T. Luo, L.-G. Pang, and X.-N. Wang, Phys. Rev. C99, 054911 (2019), arXiv:1809.02525 [nucl-th] .
  * Putschke  _et al._ (2019) J. H. Putschke _et al._ , (2019), arXiv:1903.07706 [nucl-th] .
  * Auvinen  _et al._ (2010) J. Auvinen, K. J. Eskola, and T. Renk, Phys. Rev. C82, 024906 (2010), arXiv:0912.2265 [hep-ph] .
  * Guo and Wang (2000) X.-F. Guo and X.-N. Wang, Phys. Rev. Lett. 85, 3591 (2000), arXiv:hep-ph/0005044 .
  * Zhang  _et al._ (2004) B.-W. Zhang, E. Wang, and X.-N. Wang, Phys. Rev. Lett. 93, 072301 (2004), arXiv:nucl-th/0309040 [nucl-th] .
  * Majumder (2012) A. Majumder, Phys. Rev. D85, 014023 (2012), arXiv:0912.2987 [nucl-th] .
  * Zhang and Wang (2003) B.-W. Zhang and X.-N. Wang, Nucl. Phys. A720, 429 (2003), arXiv:hep-ph/0301195 .
  * Schafer  _et al._ (2007) A. Schafer, X.-N. Wang, and B.-W. Zhang, Nucl. Phys. A793, 128 (2007), arXiv:0704.0106 [hep-ph] .
  * Lin  _et al._ (2005) Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, Phys. Rev. C72, 064901 (2005), arXiv:nucl-th/0411110 [nucl-th] .
  * Moreland  _et al._ (2015) J. S. Moreland, J. E. Bernhard, and S. A. Bass, Phys. Rev. C 92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Sjostrand  _et al._ (2008) T. Sjostrand, S. Mrenna, and P. Z. Skands, Comput. Phys. Commun. 178, 852 (2008), arXiv:0710.3820 [hep-ph] .
  * Gleisberg  _et al._ (2009) T. Gleisberg, S. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert, and J. Winter, JHEP 02, 007 (2009), arXiv:0811.4622 [hep-ph] .
  * Yang  _et al._ (2024) Z. Yang, Y. He, I. Moult, and X.-N. Wang, Phys. Rev. Lett. 132, 011901 (2024), arXiv:2310.01500 [hep-ph] .
  * Andres  _et al._ (2023) C. Andres, F. Dominguez, R. Kunnawalkam Elayavalli, J. Holguin, C. Marquet, and I. Moult, Phys. Rev. Lett. 130, 262301 (2023), arXiv:2209.11236 [hep-ph] .
  * Cacciari  _et al._ (2012) M. Cacciari, G. P. Salam, and G. Soyez, Eur. Phys. J. C72, 1896 (2012), arXiv:1111.6097 [hep-ph] .
  * Pang  _et al._ (2012) L. Pang, Q. Wang, and X.-N. Wang, Phys. Rev. C86, 024911 (2012), arXiv:1205.5019 [nucl-th] .
  * Pang  _et al._ (2015) L.-G. Pang, Y. Hatta, X.-N. Wang, and B.-W. Xiao, Phys. Rev. D91, 074027 (2015), arXiv:1411.7767 [hep-ph] .
  * Pang  _et al._ (2018) L.-G. Pang, H. Petersen, and X.-N. Wang, Phys. Rev. C 97, 064918 (2018), arXiv:1802.04449 [nucl-th] .
  * Chen  _et al._ (2010) X.-F. Chen, C. Greiner, E. Wang, X.-N. Wang, and Z. Xu, Phys. Rev. C81, 064908 (2010), arXiv:1002.1165 [nucl-th] .
  * Liu  _et al._ (2022b) F.-L. Liu, W.-J. Xing, X.-Y. Wu, G.-Y. Qin, S. Cao, and X.-N. Wang, Eur. Phys. J. C 82, 350 (2022b), arXiv:2107.11713 [hep-ph] .
  * Adam  _et al._ (2020) J. Adam _et al._ (STAR), Phys. Rev. C 102, 054913 (2020), arXiv:2006.00582 [nucl-ex] .
  * Wang (1998) X.-N. Wang, Phys. Rev. C58, 2321 (1998), arXiv:hep-ph/9804357 [hep-ph] .
  * Sirunyan  _et al._ (2021) A. M. Sirunyan _et al._ (CMS), JHEP 05, 284 (2021), arXiv:2102.13080 [hep-ex] .
  * Aaboud  _et al._ (2019) M. Aaboud _et al._ (ATLAS), Phys. Lett. B 790, 108 (2019), arXiv:1805.05635 [nucl-ex] .
  * Aad  _et al._ (2015) G. Aad _et al._ (ATLAS), Phys. Rev. Lett. 114, 072302 (2015), arXiv:1411.2357 [hep-ex] .
  * Zapp  _et al._ (2013) K. C. Zapp, F. Krauss, and U. A. Wiedemann, JHEP 03, 080 (2013), arXiv:1212.1599 [hep-ph] .
  * Zapp (2014) K. C. Zapp, Eur. Phys. J. C74, 2762 (2014), arXiv:1311.0048 [hep-ph] .
  * Chang  _et al._ (2020) N.-B. Chang, Y. Tachibana, and G.-Y. Qin, Phys. Lett. B801, 135181 (2020), arXiv:1906.09562 [nucl-th] .
  * Anderson (2023) D. Anderson (STAR), Acta Phys. Polon. Supp. 16, 55 (2023), arXiv:2212.09202 [nucl-ex] .
  * Zhang  _et al._ (2018) S.-L. Zhang, T. Luo, X.-N. Wang, and B.-W. Zhang, Phys. Rev. C98, 021901 (2018), arXiv:1804.11041 [nucl-th] .
  * Aad  _et al._ (2023) G. Aad _et al._ (ATLAS), Phys. Rev. C 107, 054908 (2023), arXiv:2205.00682 [nucl-ex] .
  * Qin and Muller (2011) G.-Y. Qin and B. Muller, Phys. Rev. Lett. 106, 162302 (2011), [Erratum: Phys. Rev. Lett.108,189904(2012)], arXiv:1012.5280 [hep-ph] .