# Probing the jet transport coefficient of cold nuclear matter in electron-ion
collisions

Peng Ru<EMAIL_ADDRESS>School of Materials and New Energy, South China
Normal University, Shanwei 516699, 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  Zhong-Bo Kang<EMAIL_ADDRESS>Department
of Physics and Astronomy, University of California, Los Angeles, California
90095, USA  Mani L. Bhaumik Institute for Theoretical Physics, University of
California, Los Angeles, California 90095, USA  Center for Frontiers in
Nuclear Science, Stony Brook University, Stony Brook, New York 11794, USA
Enke Wang<EMAIL_ADDRESS>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  Hongxi Xing<EMAIL_ADDRESS>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  Ben-Wei Zhang
<EMAIL_ADDRESS>Key Laboratory of Quark $\&$ 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

###### Abstract

We present a study of the nuclear-medium induced transverse momentum
broadening of particle production in future electron-ion-collision (EIC)
experiments. By considering the multiple scattering between hard partons and
cold nuclear medium within the higher-twist factorization framework in
perturbative QCD, we calculate the transverse momentum broadening of single
hadron production in semi-inclusive measurements, as well as the nuclear
enhancement of the transverse momentum imbalance for di-hadron and heavy-meson
pair productions. In particular, a kinematics dependent non-perturbative jet
transport coefficient $\hat{q}=\hat{q}(x,Q^{2})$ extracted in a global
analysis of the current data, together with its uncertainty determined with a
Hessian method, are input into our calculations and are available for the
community. Significant kinematic and color-state dependences of the nuclear
induced broadening/imbalance are predicted. Our results indicate that the
future EIC measurements are able to provide powerful constraints on the
kinematic dependence of the transport coefficient $\hat{q}$ and thus greatly
facilitate the jet tomography of cold nuclear medium.

##  I Introduction

Exploring at a femtometer scale the properties of the nuclear media in
different matter phases of quantum chromodynamics (QCD), such as cold nucleus,
hadron gas and hot/dense quark-gluon plasma, is one of the main goals of
various high-energy nuclear collisions [1, 2], including lepton-, hadron- and
nucleus-nucleus collisions. Benefited by the factorization in perturbative QCD
[3], the particle(s) produced with a large momentum transfer, such as a parton
jet, can serve as a well-controlled hard probe of the non-perturbative
property of the nuclear medium. The multiple scattering that a hard probe
undergoes in nuclear medium can lead to the transverse momentum broadening and
energy loss of the probe in general [4, 5], which are reflected in the
observed nuclear modifications on the spectra and substructures of jets (or
hadrons) [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. An important medium
property commonly embodied in these effects is the jet transport property [17,
18, 19, 20], quantified as the coefficient $\hat{q}$, which characterizes the
transverse momentum broadening of a (quark) jet per unit propagation length in
the medium and thus measures the strength of the interaction between the probe
and nuclear medium.

Transport coefficient $\hat{q}$ has been an iconic quantity to represent the
medium property seen by jets for a long time, especially in the study of
heavy-ion collisions [20, 21, 22, 23, 24, 25, 26, 27, 28]. Recently, the study
of its dependence on kinematic variables like jet energy and probing scale
became active [27, 28, 29, 32, 33, 34, 31, 36, 35, 37, 38, 30, 39], along with
some related theoretical progress [40, 44, 45, 46, 41, 42, 43]. On this
aspect, the electron-nucleus ($e$A) and proton-nucleus ($p$A) collisions is of
particular importance [32, 33, 34, 36, 35, 37, 44, 45, 46, 47, 48, 49, 50],
since they provide a relative clean environment to delicately study the
kinematic dependence of the transport property of cold nuclear matter and to
test the theoretical framework of the jet-medium interaction, which in turn
can be instructive for the study of nucleus-nucleus (AA) collisions.

In our previous work [32], we performed the first global extraction of the
$\hat{q}$ in cold nuclear matter from the current data in $e$A and $p$A
collisions, mainly on various types of nuclear-induced transverse momentum
broadening. We found that, with a $\hat{q}$ depending on the parton momentum
fraction $x$ and probing scale $Q^{2}$, the theoretical calculations within
the higher-twist expansion formalism can give an overall good agreement with
the world data. The extracted optimal $\hat{q}(x,Q^{2})$ shows significant
enhancements in small- and large-$x$ regions and a mild $Q^{2}$ dependence.
However, the uncertainties of $\hat{q}(x,Q^{2})$ were not yet worked out in
Ref. [32]. Besides, since most of the current data are gathered in the
intermediate $x$ and $Q^{2}$ region, the suggested universality and kinematic
dependence of $\hat{q}(x,Q^{2})$ should be examined in future experiments with
broader $x$ and $Q^{2}$ coverage. These issues motivate our followup study
presented in this paper.

In the first part of this work, we will upgrade our global analysis of
$\hat{q}$ by determining the uncertainties of the extracted $\hat{q}(x,Q^{2})$
with the Hessian matrix method [51, 52], which will generate an uncertainty
set of $\hat{q}(x,Q^{2})$. The Hessian analysis is further upgraded by
including a new data set on $J/\psi$ production published recently [53] to
strengthen the experimental constraints in small $x$ region. Through the
Hessian analysis in this work, we not only learn how the uncertainty of
$\hat{q}(x,Q^{2})$ varies with $x$ and $Q^{2}$, but also give a complete
theoretical prediction with the uncertainty for a related observable, which is
important for the future experimental examination.

The future experiments of electron-ion collisions (EIC) will be an important
place to examine our results [1, 54, 55]. Several future EIC facilities have
been proposed or under construction, e.g., the Electro-Ion Collider in US (US-
EIC) [1], Electron-ion collider in China (EicC) [56] and Jefferson Lab 12 GeV
program (JLab) [57], etc. Since these experiments provide a wide kinematic
($x$ and $Q^{2}$) coverage and high-precision measurements, they are expected
to significantly improve our understanding of the transport property of cold
nuclear matter.

In EIC experiments, the nuclear-induced transverse momentum broadening is
still the type of observable most directly relevant to the transport
coefficient $\hat{q}$. Therefore, in the second part of this work, we study
the transverse momentum broadening of single particle and back-to-back
particle pair productions in future electron-ion collisions. The latter case
is also equivalent to the nuclear enhancement of the transverse momentum
imbalance of the particle pair. Our calculations employ the formalism of
higher-twist expansion, i.e., the generalized factorization in perturbative
QCD, which has been well established [58, 59, 60, 62, 61] and developed [44,
45, 46, 47, 48, 49], and has shown good applicability on describing the
nuclear modification in both cold and hot/dense nuclear matter [47, 65, 20,
32, 63, 64]. In the calculations, we use the $\hat{q}$ together with its
uncertainties extracted from our global analysis as an input. For comparison,
we also provide the predictions with a kinematics independent constant
$\hat{q}$. By preliminarily estimating the experimental uncertainties, we show
the potential of the future EIC measurements to provide powerful constraints
on the $\hat{q}$ in a wide kinematic range.

The rest of this paper is organized as follows. In Sec. II, we briefly review
the framework of our global analysis of $\hat{q}$ and the main results in our
previous analysis (Sec. II.1), and test the extracted $\hat{q}$ with a new
data set from the ALICE experiment at the Large Hadron Collider (Sec. II.2).
We then present our work on Hessian analysis of the uncertainty of $\hat{q}$
(Sec. II.3). The jet energy dependence of $\hat{q}$ is also discussed at the
end of this section. In Sec. III, we present the study of three types of
nuclear-induced transverse momentum broadening/imbalance (Sec. III.1-III.3)
for three EIC facilities: US-EIC, EicC, and JLab, and discuss the advantage of
future EICs in understanding the kinematic dependence of $\hat{q}$ (Sec.
III.4). We give a summary and discussion in Sec. IV.

##  II Global analysis of $\hat{q}$ for cold nuclear matter

In this section we perform an updated global analysis of the jet transport
coefficient $\hat{q}$ for cold nuclear matter (CNM) with the current world
data from electron-nucleus and proton-nucleus collisions. Through the
analysis, a kinematics dependent $\hat{q}$ is extracted and its uncertainties
are estimated with Hessian matrix, which will be used for making predictions
for the EIC observables in the next section.

###  II.1 Framework and previous results of the analysis

First we briefly review the framework of our analysis as well as what have
been done in our previous work [32].

In high-energy $e$A and $p$A collisions, $\hat{q}$ is a key non-perturbative
input in theoretical descriptions of the multiple scattering between the hard
probe and the partons inside the nuclear target. Generally, in absence of the
multiple scattering effects, the cross section of a hard scattering process in
$p$A collisions can be written schematically as follows using leading twist
collinear factorization formalism

$\displaystyle~{}~{}~{}~{}~{}{d\sigma^{S}}=f_{q(g)/p}\otimes f_{q(g)/A}\otimes
d\hat{\sigma}^{\textrm{S}}\otimes D_{h}\,,$ (1)

where the superscript “$S$” denotes the single-scattering process,
$d\hat{\sigma}^{\textrm{S}}$ represents the perturbatively calculable partonic
cross section, $f$ and $D_{h}$ represent the involved parton distribution
functions (PDFs) in initial state and fragmentation functions (FFs) in final
state, respectively, and the subscript ‘${q(g)/A}$’ indicates an incoming
quark (gluon) from the nucleus. The leading twist collinear factorization
formalism underlies the successful global analyses for the PDFs and FFs [66,
67, 68, 52, 69, 70].

In the presence of a large nucleus, the parton multiple scattering effects
become important and can be formulated by generalizing the collinear
factorization in the higher-twist expansion approach [58, 59, 60, 62, 61].
Specifically, let us consider the transverse momentum broadening
$\Delta\langle p_{T}^{2}\rangle$, usually defined as the difference of the
averaged transverse momentum square of the produced particle between $e$A
($p$A) and $ep$ ($pp$) collisions,

$\displaystyle\Delta\langle p_{T}^{2}\rangle_{eA/pA}=\langle
p_{T}^{2}\rangle_{eA/pA}-\langle p_{T}^{2}\rangle_{ep/pp}\,.$ (2)

The leading contribution to the broadening comes from the double scattering
effects enhanced by the nuclear size. For example, in semi-inclusive deep
inelastic scattering (SIDIS), the struck quark that is kicked off by the
virtual photon may experience additional interactions with the partons inside
the nuclear target, resulting in the $\Delta\langle p_{T}^{2}\rangle$ of the
final-state hadrons.

In higher-twist expansion approach, the leading contribution of transverse
momentum broadening can be written generically in the form of a ratio as [71]

$\displaystyle\Delta\langle p_{T}^{2}\rangle\approx\frac{d\langle
p_{T}^{2}\sigma^{D}\rangle}{d\mathcal{PS}}\bigg{/}\frac{d\sigma^{S}}{d\mathcal{PS}}\,,$
(3)

where the denominator ${d\sigma^{S}}\\!/{d\mathcal{PS}}$ is the leading-twist
single scattering cross section in the phase space volume $d\mathcal{PS}$ in
$e$A or $p$A collisions, and the numerator ${d\langle
p_{T}^{2}\sigma^{D}\rangle}/{d\mathcal{PS}}$ is the $p_{T}^{2}$-weighted
double scattering cross section [44, 45, 46, 47, 48]

$\displaystyle\frac{d\langle
p_{T}^{2}\sigma^{D}\rangle}{d\mathcal{PS}}\equiv\int\\!dp_{T}^{2}p_{T}^{2}\frac{d\sigma^{D}}{d\mathcal{PS}dp_{T}^{2}}\,.$
(4)

which can be written as follows [74, 75]

$\displaystyle{d\langle p_{T}^{2}\sigma^{D}\rangle}=f_{q(g)/p}\otimes
T_{ij}\otimes d\hat{\sigma}^{\textrm{D}}\otimes D_{h}\,,$ (5)

where $T_{ij}$ represents the nuclear twsit-4 (T4) parton-parton correlation
functions, which are universal non-perturbative functions encoding the medium
properties characterized by the jet transport coefficient $\hat{q}$. The
$p_{T}^{2}$-weighted cross section, thus the broadening, also depends on the
color representation of the hard probe, i.e., the quark and gluon jets
correspond to the color factors $C_{F}=(N_{c}^{2}-1)/2N_{c}$ and
$C_{A}=N_{c}$, respectively [44, 45, 46, 47, 48].

With the assumption of a loosely bound large nucleus, one can neglect the
momentum and spatial correlations among the nucleons in the nuclear target.
Therefore, the twist-4 matrix element can be effectively factorized in terms
of leading twist PDFs and $\hat{q}$. For example, one can approximate the
twist-4 quark-gluon correlation function $T_{qg}$ as [45]

$\displaystyle
T_{qg}(x,0,0,\mu^{2})\approx\frac{9R_{A}}{8\pi^{2}\alpha_{s}}f_{q/A}(x,\mu^{2})\hat{q}(x,\mu^{2}),$
(6)

where $x$ is the momentum fraction carried by the quark that enters into the
hard interaction, $\mu$ is the factorization scale in perturbative QCD,
$R_{A}$ is the radius of the nucleus with mass number $A$, $f_{q/A}$ is the
parton distribution function of quark $q$ in the nucleus, and
$\hat{q}(x,\mu^{2})$ is the nuclear geometry averaged jet transport
coefficient that we hope to extract. Similarly, for gluon-gluon correlation
function $T_{gg}$, which accesses the process initiated with a gluon in the
nucleus, we assume a same form as Eq. (6) in our study, with $f_{q/A}$
replaced by the gluon distribution $f_{g/A}$. Through our this paper, the
$\hat{q}$ represents the transport coefficient of a quark jet.

With the theoretical framework introduced above, one can perform a global
analysis for $\hat{q}$ similar to what have been done for PDFs and FFs [66,
67, 68, 52, 69, 70]. In particular, since both the twist-4 correlation
functions and the leading-twist PDFs are universal non-perturbative quantities
depend on $x$ and $\mu^{2}$, the $\hat{q}$ from Eq. (6) naturally involves
possible kinematic dependence on momentum fraction and probing scale. Note
that in phenomenological studies [62, 72, 73] $\hat{q}$ was usually assumed to
be a constant value due to the unknown kinematic dependence.

Particle transverse momentum broadening is the type of observable most
directly relevant to $\hat{q}$ for CNM [74, 75]. Our previous global analysis
[32] takes into account the current world data on the transverse momentum
broadening of single hadron production in SIDIS [6] of $e$A collisions, of
Drell-Yan di-lepton production in $p$A collisions [77, 76], and of heavy
quarkonium ($J/\psi$ and $\Upsilon$) production in $p$A collisions [78, 80,
79, 77, 7, 8]. These observables involve multiple scatterings undergone by
quark or gluon jets in the initial or (and) final states of the hard
processes, providing multi-dimensional insight into the transport coefficient
as well as a place to examine the theoretical framework.

Besides the data on transverse momentum broadening, a set of data on the
nuclear modification factor (shadowing effect) of the DIS structure function
[81, 82] is also included in our analysis. In the higher-twist framework, such
a nuclear suppression can be related to the coherent multiple scattering,
which has been calculated by resumming the higher-twist contributions [72],
thus, is also sensitive to the transport coefficient.

In total, there are 215 data points in the analysis [32] from the experiments
at DESY, FNAL, SPS, RHIC and LHC. In particular, different observables or
measurements involve different kinematic regions identified with the momentum
fraction $x$ of the nuclear parton and the probing scale $Q^{2}$, providing
possibility to explore the kinematic dependence of $\hat{q}(x,Q^{2})$.

To address the kinematic dependence of $\hat{q}$ in the global analysis, we
adopt the parametrization form [32]

$\displaystyle\hat{q}(x,\mu^{2})=\hat{q}_{0}\,\alpha_{s}(\mu^{2})\,x^{\alpha}(1-x)^{\beta}\left[\ln(\mu^{2}/\mu_{0}^{2})\right]^{\gamma}\,,$
(7)

which involves four free parameters, $\hat{q}_{0},~{}\alpha,~{}\beta$, and
$\gamma$ to be determined by the experimental data. Such a functional form is
primarily motivated with several physical considerations. For example, in the
small-$x$ region, we expect that $\hat{q}$ depends on the gluon saturation
scale, which exhibits a power-law behavior as $Q_{s}^{2}\propto x^{-1/3}$
[83]. This feature is related to the factor $x^{\alpha}$ in Eq. (7). At large
$x$, the QCD power corrections could be different [84, 85, 86, 87], and may
result in different behavior of $\hat{q}$, which is considered by including
the factor $(1-x)^{\beta}$. Moreover, a logarithmic scale dependence of
$\hat{q}$ is suggested from the radiative corrections [42, 41, 88], thus we
use a factor $[\ln(\mu^{2}/\mu_{0}^{2})]^{\gamma}$ in Eq. (7), where the
exponent $\gamma$ is included to account for potential modification at the
higher-order in perturbative corrections and/or non-perturbative
contributions. Although the QCD scale evolution equation of the twist-4 quark-
gluon correlation function has been derived in a series of previous work [44,
45, 46], it is coupled with the gluon-gluon correlator whose evolution is not
determined. In Eq. (7), the $\alpha_{s}(\mu^{2})$ is introduced to offset the
$\alpha_{s}$ in the denominator of the correlation function in Eq. (6), and
$\mu_{0}=1$ GeV is introduced to make the argument of the logarithm
dimensionless [32]. On the whole, the parametrization form in Eq. (7) has some
generality and is similar to what is usually used in the extraction of other
non-perturbative quantities, such as the parton distribution functions [66].
We note that such a parametrization form and several similar forms of
$\hat{q}$ have been recently applied by other groups, and are shown to work
well in their studies [35, 37].

In the analysis [32], the theoretical calculation of the transverse momentum
broadening is performed at twist-4 level in QCD power expansion and at leading
order (LO) in perturbative expansion with $\alpha_{s}$. A complete NLO
calculation is not yet available. In the calculation, we use the CT14 LO
parton distribution functions with 3 active quark flavors [89], and the DSS
fragmentation functions [67]. The heavy quarkonium production is calculated
with the color evaporation model [47]. Since the transverse momentum
broadening is expressed as a ratio in Eq. (3), it is to some extent
insensitive to the non-perturbative inputs like the PDFs, FFs, and quarkonium
production model. More details of the calculation can be seen in Ref. [32]. In
addition, the possible hadronization of the in-medium jet may weaken the
multiple scattering effects [90], which has not been considered in the current
framework.

The global analysis of $\hat{q}$ starts with finding the optimal $\hat{q}$ by
minimizing the $\chi^{2}$ as a function of the free parameters $\\{a_{j}\\}$
defined as [51, 52]

$\displaystyle\chi^{2}(\\{a_{j}\\})=\sum_{i}\frac{\left[\mathcal{D}_{i}-\mathcal{T}_{i}(\\{a_{j}\\})\right]^{2}}{\sigma_{i}^{2}},$
(8)

where $\mathcal{D}_{i}$ is the value of the $i$-th experimental data point,
$\mathcal{T}_{i}(\\{a_{j}\\})$ is the corresponding theoretical prediction
depending on the values of the free parameters
$\\{a_{j}\\}=\\{\hat{q}_{0},~{}\alpha,~{}\beta,~{}\gamma\\}$ in the $\hat{q}$
parametrization in Eq. (7), and $\sigma_{i}^{2}$ is the statistical and
systematic experimental uncertainties summed in quadrature. The influence of
the possible correlated experimental uncertainties has not been taken into
account in our analysis [91, 52]. The procedure of minimizing the $\chi^{2}$
in the $\\{a_{j}\\}$ space is performed by utilizing the MINUIT package [92].

An optimal $\hat{q}(x,Q^{2})$ within the parametrization form in Eq. (7) is
found at a minimum total $\chi^{2}=260$ ($\chi^{2}$/NDP$=1.21$, where
NDP$=215$ is the number of data points [32]). The optimal $\hat{q}(x,Q^{2})$
exhibits an obvious $x$ dependence, especially in small- and large-$x$
regions, as well as a mild $Q^{2}$ dependence [32] (also seen in Sec. II.3 of
this manuscript). The theoretical results with this $\hat{q}(x,Q^{2})$ show an
overall good agreement with the experimental data, indicating a universal
kinematic dependence of $\hat{q}(x,Q^{2})$ in cold nuclear matter. To further
clarify this kinematic dependence, we also performed the fitting by assuming
$\hat{q}$ is a constant quantity as $\hat{q}=\hat{q}_{0}$, and we found a
minimum total $\chi^{2}=388$ ($\chi^{2}$/NDP$=1.8$), which is apparently
larger than that with the kinematic dependence. Especially, the calculations
with the constant $\hat{q}$ can not give a good description of the data in
small- and large-$x$ regions. For example, the $\chi^{2}$ for the $J/\psi$
broadening at the LHC is 87.3 ($\chi^{2}$/NDP$=7.3$ with NDP$=12$), which is
far from reasonable. In contrast, the result for this part is significantly
improved by using the kinematics dependent $\hat{q}$, reflected by the
$\chi^{2}=4.8$ ($\chi^{2}$/NDP$=0.4$) [32].

The results of the analysis should be examined in future, when more
experimental data with a wider kinematic coverage and a high precision, e.g.,
from EIC, are available. To this end, one needs the uncertainty of the
$\hat{q}(x,Q^{2})$ under the constraints of the current data, which allows to
make a complete theoretical prediction. In the previous work [32], a Lagrange
multiplier method [91, 51] is employed to evaluate the uncertainties of part
of the calculated observable. However, the uncertainties of $\hat{q}(x,Q^{2})$
(varying with $x$ and $Q^{2}$) can not be easily obtained with that method,
which makes the analysis less predictive and thus motivates our reanalysis
with the Hessian matrix method (Sec. II.3). Before presenting that, we show a
test for the extracted $\hat{q}$ with a new data set in the following
subsection.

###  II.2 A test for $\hat{q}(x,Q^{2})$ with new $J/\psi$ data

Figure 1: Transverse momentum broadening $\Delta\langle p_{T}^{2}\rangle$ in
$J/\psi$ production in backward [panel (a)] and forward [panel (b)] regions in
$p$A collisions at the LHC as functions of centrality. Solid curve represents
results with $\hat{q}(x,Q^{2})$ form global analysis [32] and shaded area
corresponds to uncertainty. Dashed line and dotted boundaries represent
results with constant $\hat{q}$ with uncertainty. Circles are new ALICE data
[53] that are not included in previous analysis of $\hat{q}$ [32].

We noted that the ALICE collaboration have published in 2021 a new data set on
the transverse momentum broadening of $J/\psi$ production in $p$-Pb collisions
at $\sqrt{s_{NN}}=8.16$ TeV at the LHC [53], which provides an opportunity to
test the extracted $\hat{q}$ in [32].

In Fig. 1, we show both the theoretical results with the kinematic dependent
$\hat{q}(x,Q^{2})$ and the constant $\hat{q}$ extracted from previous
analysis, confronting with the new data. We find the calculations with the
$\hat{q}(x,Q^{2})$ give a visible rapidity dependence of the broadening,
reasonably description of the data in both backward and forward rapidity
regions. However, the results with the constant $\hat{q}$ obviously
underestimate the broadening in forward (small $x$) region.

A similar rapidity dependence can also be observed in the earlier ALICE
$J/\psi$ data measured at $\sqrt{s_{NN}}=5.02$ TeV (2015) [8], which have been
taken into account in our previous analysis and provided important information
on the $x$ dependence of $\hat{q}$. In this work, we will also include the new
data in the Hessian analysis to provide more constraints on $\hat{q}$.

###  II.3 Uncertainty of $\hat{q}(x,Q^{2})$ from Hessian Matrix

In order to estimate the uncertainties of the kinematics dependent
$\hat{q}(x,Q^{2})$ under the constraints of the current data, we perform an
analysis with a Hessian matrix method. The basic assumption of the Hessian
analysis is that the $\chi^{2}$ can be approximately expressed in a quadratic
form of the free parameters $\\{a_{i}\\}$ in the neighborhood of the minimum
[51, 52] as

$\displaystyle\chi^{2}=\chi_{0}^{2}+\sum_{i,j}H_{ij}y_{i}y_{j}\,,$ (9)

where $\chi_{0}^{2}\equiv\chi^{2}(\\{a_{i}^{0}\\})$ is the global minimum of
the $\chi^{2}$ at the optimal parameter values $\\{a_{i}\\}=\\{a_{i}^{0}\\}$,
$y_{i}=a_{i}-a_{i}^{0}$ is the displacement of $a_{i}$ from its optimal value
$a_{i}^{0}$, and $H_{ij}$ are the elements of Hessian matrix defined as

$\displaystyle H_{ij}=\frac{1}{2}\left(\frac{\partial^{2}\chi^{2}}{\partial
y_{i}\partial y_{j}}\right)_{a_{i}=a_{i}^{0}}.$ (10)

Usually there are interplays among different variables $y_{i}$ (or $a_{i}$) in
the $\chi^{2}$, and the off-diagonal Hessian matrix elements could be non-
zero. This makes the uncertainty estimation, corresponding to a certain
tolerance $\Delta\chi^{2}\equiv\chi^{2}-\chi^{2}_{0}$, not that
straightforward. However, one can disentangle the parameters by defining a new
basis $\\{z_{i}\\}$ of the parameter space, in whose representation the
Hessian matrix is diagonal. With the new set of parameters $\\{z_{i}\\}$, the
$\Delta\chi^{2}$ can be written in a simple form as

$\displaystyle\Delta\chi^{2}=\sum_{i}z_{i}^{2},$ (11)

which means that the contours of the $\chi^{2}$ are spheres in the new basis.
Using this new basis, one can generate the corresponding uncertainty sets of
$\hat{q}$, which can be used to estimate both the uncertainties of $\hat{q}$
and the related theoretical predictions. More details of the Hessian analysis
can be found in Appendix A.

Figure 2: $\Delta\chi^{2}=\chi^{2}-\chi^{2}_{0}$ as a function of each
individual original parameter $a_{i}$ in neighborhood of minimum
$\chi^{2}=\chi^{2}_{0}$. Here
$\\{a_{i}\\}=\\{\hat{q}_{0},\alpha,\beta,\gamma\\}$. Dashed line represents
tolerance $\Delta\chi^{2}=35$.

The Hessian analysis in this work is performed by using the MINUIT package
combined with the ITERATE program [51]. The theoretical framework for
calculating various observables is the same as in our previous analysis
introduced in Sec II.1. After the new ALICE data on $J/\psi$ production (Sec.
II.2) included (now 227 data points in total), a global minimum
$\chi_{0}^{2}=275$ is reached through the analysis.

We first show in Fig. 2 the global $\chi^{2}-\chi_{0}^{2}$ as functions of
each individual original parameters $\\{a_{i}\\}$ in the vicinity of the
minimum. We can see that the experimental data indeed have sensitivities to
all the parameters. The non-zero optimal values of the parameters $\alpha$,
$\beta$ and $\gamma$ suggest the kinematic dependence of $\hat{q}$ on $x$ and
$Q^{2}$, to be further consolidated with the determined uncertainty.

Figure 3: $\Delta\chi^{2}=\chi^{2}-\chi^{2}_{0}$ as functions of new
parameters $\\{z_{i}\\}$ from Hessian analysis. Grey curve shows a quadratic
function $\Delta\chi^{2}=z^{2}$ for reference.
Figure 4: $\hat{q}(x,Q^{2})$ extracted from global analysis, shown as a
function of momentum fraction $x$ of initial-state nuclear parton at
$Q^{2}=$1.2, 10, and 100 GeV2. Green solid curve represents optimal values
($S_{0}$), light-green band represents corresponding uncertainty, and dotted
curves show uncertainty set $\\{S_{\pm k}\\}$ of $\hat{q}(x,Q^{2})$. For
reference, kinematics independent constant $\hat{q}$ and uncertainty extracted
from data are shown as red dashed line and pink band, respectively. Both
uncertainties of $\hat{q}(x,Q^{2})$ and constant $\hat{q}$ correspond to
$90\%$ C.L..  Figure 5: $\hat{q}(x,Q^{2})$ at $Q^{2}=10$ GeV2 and
uncertainties extracted from global analysis with and without ALICE data
(2021) on transverse momentum broadening of $J/\psi$ production. Results of
new analysis are the same as in Fig. 4. Results without $J/\psi$ (21’) data
are shown with a dashed curve with dotted boundaries.

Fig. 3 shows the values of $\chi^{2}-\chi_{0}^{2}$ as functions of the new
parameters $\\{z_{i}\\}$ defined in Eq. (35). They are found to have a very
good agreement with the quadratic form $\Delta\chi^{2}=z_{i}^{2}$, indicating
the good performance of Hessian analysis.

With this Hessian analysis, we obtain the optimal values of the parameters
$\\{a_{i}\\}$ together with their uncertainties corresponding to $90\%$
confidence level (C.L.):

$\displaystyle\hat{q}_{0}=0.0191\\!\pm
0.0061~{}\textrm{GeV}^{2}\\!/\textrm{fm},~{}~{}\alpha=-0.182\pm 0.050,$ (12)
$\displaystyle\beta=-2.85\pm 1.87,~{}~{}~{}\gamma=0.264\pm
0.169\,.~{}~{}~{}~{}~{}~{}$ (13)

We find that with the uncertainties, the global data favor negative $\alpha$
and $\beta$, and a positive $\gamma$ in the parametrization of $\hat{q}$. The
optimal values of the parameters are in good agreement with the results of our
previous analysis [32].

Figure 6: Transverse momentum broadening $\Delta\langle p_{T}^{2}\rangle$, in
Drell-Yan process in $p$A collisions versus nuclear mass number $A$ [panel
(a)], in SIDIS versus Bjorken $x_{B}$ [panel (b)], and in $J/\psi$ production
in backward [panel (c)] and forward [panel (d)] regions in $p$A collisions at
the LHC (5TeV, 2015) as functions of $N_{\textrm{coll}}$. Green solid curve
represents results with optimal $\hat{q}(x,Q^{2})$ and light-green shaded area
corresponds to uncertainty. Red dashed line and dotted band represent results
with constant $\hat{q}$ with uncertainty. Averaged momentum fraction $\langle
x\rangle$ and probing scale $\langle Q^{2}\rangle$ are shown for reference.

In Fig. 4, we show the extracted $\hat{q}(x,Q^{2})$ versus the momentum
fraction $x$ for $Q^{2}=$1.2, 10, and 100 GeV2, including the optimal values
($S_{0}$), uncertainty bands, and uncertainty sets $S_{\pm k}$ of
$\hat{q}(x,Q^{2})$ [defined by Eq. (39)], which can be used to estimate the
uncertainty for theoretical predictions [see Eq. (38)]. For the optimal
values, we can see that the $x$ dependence is noticeable in small- and
large-$x$ regions, and the $Q^{2}$ dependence is relatively mild. The
enhancements of $\hat{q}(x,Q^{2})$ in small- and large-$x$ regions are related
to the negative parameters $\alpha$ and $\beta$, respectively. The $\alpha$
value in Eq. (13) is qualitatively consistent with the growth rate of the
gluon density expected in saturation physics [93], and the negative $\beta$
may indicate an enhancement of the nuclear power correction at large $x$ [84,
85, 86, 87]. Since most of the current data are located in the intermediate
$x$ and $Q^{2}$ regions, the uncertainties of $\hat{q}(x,Q^{2})$ become larger
at small and large values of $x$ or $Q^{2}$, due to the less experimental
constraints. Especially, the uncertainties are dramatically large at $x\gtrsim
0.4$, where no data exist for this kinematic region. For comparison, we also
show in Fig. 4 the $\hat{q}$ extracted by assuming it is a kinematic
independent constant quantity
($\hat{q}=\hat{q}_{0}=0.0150^{+0.0023}_{-0.0025}$ GeV2/fm).

Figure 7: Extracted $\hat{q}$ as a function of jet energy $E_{\textrm{jet}}$.
Solid curve and shaded band represent optimal values and uncertainties,
respectively.

To illustrate the impact of the added new $J/\psi$ data on the extracted
$\hat{q}(x,Q^{2})$, we compare in Fig. 5 the $\hat{q}(x,Q^{2})$ for $Q^{2}=10$
GeV2 extracted with and without the new $J/\psi$ data. A good agreement
between their optimal values is found. What is notable is that, with the
additional constraints from the new data, the uncertainties at small values of
$x$ are reduced to some extent and a more evident $x$ dependence is observed.
More detailed illustrations for the impact of the new $J/\psi$ data on
$\hat{q}(x,Q^{2})$ at different $Q^{2}$ and on the predictions for different
observables can be seen in Appendix B. In addition, for theoretical
predictions, we have compared the uncertainties from the Hessian analysis with
those from our previous analysis via the Lagrange multiplier method, and found
a good agreement between them [94].

###  II.4 Kinematic dependence of $\hat{q}$

The kinematic dependence of $\hat{q}$ is of particular interest since it is
related to the detailed partonic structures of the nuclear matter, similar to
the $x$ and $Q^{2}$ dependence of a parton distribution function. Global
analysis offers an indispensable data-driven understanding of such issues. To
demonstrate how the experimental data determine the kinematic dependence of
the $\hat{q}$ in our analysis, we show in Fig. 6 the results calculated with
the extracted $\hat{q}(x,Q^{2})$ and constant $\hat{q}=\hat{q}_{0}$, for four
representative observables in the analysis: the transverse momentum broadening
in Drell-Yan process (nuclear mass number dependence), in SIDIS (Bjorken $x$
dependence), and in $J/\psi$ production [dependence on the collision-
centrality related $N_{\textrm{coll}}$ in backward and forward rapidity at the
LHC]. For the Drell-Yan [panel (a)] and the backward $J/\psi$ production
[panel (c)], since the involved momentum fraction $x$ of the nuclear parton is
in the intermediate region, the results with the $\hat{q}(x,Q^{2})$ and the
constant $\hat{q}$ are close to each other, while a slightly better agreement
with the Drell-Yan data is given by that with the $\hat{q}(x,Q^{2})$. However,
the theoretical predictions are significantly improved with the kinematics
dependent $\hat{q}(x,Q^{2})$ in the SIDIS [panel (b)] and forward $J/\psi$
production [panel (d)], which correspond to the regions of large and small
$x$, respectively. In particular, the calculation with the constant $\hat{q}$
from the global analysis completely fails to describe the data on the forward
$J/\psi$ production at the LHC. Actually, the enhancements of the
$\hat{q}(x,Q^{2})$ in small and large $x$ regions observed in Fig. 4 stem
largely from the $J/\psi$ production in forward region at the LHC and the
Bjorken-$x$ dependence of SIDIS in the analysis, respectively, which should be
examined through future experiments that involve small and large $x$ regions.

On the other hand, it is also noteworthy that, for the observables that are
taken into account in our analysis, i.e., the transverse momentum broadening
in SIDIS, Drell-Yan process and heavy-quarkonium production, the energy of the
hard probe $E_{\textrm{jet}}$ in the nucleus rest frame can be commonly
expressed with the Lorentz invariant variables $x$ and $Q^{2}$ through the
relationship:

$\displaystyle E_{\textrm{jet}}=\frac{Q^{2}}{2m_{p}x}\,,$ (14)

where the hard scale $Q^{2}$ is the virtuality (or squared invariant mass) of
the virtual photon in SIDIS and Drell-Yan process, and is the squared
invariant mass of the heavy-quark pair that form the quarkonium in the color
evaporation model, $x$ is the momentum fraction of the initial-state nuclear
parton in these processes, and $m_{p}$ is the nucleon mass. With Eq. (14), we
can convert the extracted $\hat{q}(x,Q^{2})$ into the form
$\hat{q}(E_{\textrm{jet}},Q^{2})$, which can be regarded as the jet energy
dependence of $\hat{q}$ in cold nuclear matter. Since the jet energy
dependence is usually discussed in the study of jet quenching in quark-gluon
plasma in relativistic heavy-ion collisions [21, 40, 27], to extract the
$\hat{q}(E_{\textrm{jet}},Q^{2})$ in cold nuclear matter will provide a
reference for the future comparative study. In Fig. 7 we plot the
$\hat{q}(E_{\textrm{jet}},Q^{2})$ with uncertainties. We find that for
$Q^{2}=2-10$ GeV2 the $\hat{q}$ increases with $E_{\textrm{jet}}$ in a wide
range of jet energy corresponding to the small $x$ region. However, for
$Q^{2}=100$ GeV2, the plotted region corresponds to large $x$ values, and the
$\hat{q}$ decreases with $E_{\textrm{jet}}$ with a large uncertainty. Some
similar results in the study of jet quenching can be found in Refs. [21, 40,
27]. Besides, the $Q^{2}$ dependence in $\hat{q}(E_{\textrm{jet}},Q^{2})$ is
more pronounced than that in $\hat{q}(x,Q^{2})$, because the corresponding $x$
will vary with $Q^{2}$ for a certain $E_{\textrm{jet}}$.

Trough the global analysis discussed in this section, a kinematics dependent
transport coefficient $\hat{q}$ in cold nuclear matter and its uncertainties
have been extracted from the current experimental data, which is expected to
motivate the future theoretical and experimental studies to further understand
and constrain the $\hat{q}(x,Q^{2})$, and to consolidate the universality of
the kinematic dependence of $\hat{q}$. Next, we will study related observables
in future electron-ion collisions.

##  III Nuclear induced transverse momentum broadening/imbalance in electron-
ion collisions

The future EIC facilities [1, 56, 57] provide great opportunities to deepen
our understanding of the jet transport property of the cold nuclear medium. In
EIC experiments, the nuclear-medium induced transverse momentum broadening
will continue to be the observable most directly related to the transport
coefficient $\hat{q}$. Since the initial-state projectile is an electron in
EIC, the broadening is induced by the final-state multiple scattering between
the outgoing hard probe and the nucleus.

In this section, we will study the transverse momentum broadening in both
single- and pairwise-particle productions at the EIC. For the latter case, the
nuclear broadening of the particle-pair is equivalent to the nuclear
enhancement of the particle-pair transverse momentum imbalance. Concretely,
using the $\hat{q}$ extracted in the previous section, we will calculate the
transverse momentum broadening of single-hadron production and the enhancement
of transverse momentum imbalance of di-hadron and heavy-meson pair
($D\bar{D}$) productions. These observables will be studied in the kinematic
regions of three proposed EIC facilities: US-EIC, EicC, and JLab (12GeV) [1,
57, 56].

Since we focus on the EIC, it will be useful to give the typical Lorentz-
invariant DIS kinematic variables

$\displaystyle
x_{B}=\frac{Q^{2}}{2p_{N}\cdot\\!q_{\gamma}}\,,~{}~{}y=\frac{q_{\gamma}\cdot\\!p_{N}}{k_{e}\cdot\\!p_{N}}\,,~{}~{}Q^{2}=-q^{2}_{\gamma}\,.$
(15)

Here $x_{B}$ is the Bjorken variable, $y$ is the inelasticity of the
scattering, and $Q^{2}$ is the virtuality of the exchanged photon
$\gamma^{\ast}$. In their definitions, $k_{e}$, $p_{N}$ and $q_{\gamma}$ are
the 4-momenta of the incoming electron, the nucleon and the virtual photon,
respectively. For the final-sate fragmentation process, the hadron momentum
fraction $z_{h}$ is introduced as

$\displaystyle z_{h}=\frac{p_{N}\cdot\\!p_{h}}{p_{N}\cdot\\!q_{\gamma}},$ (16)

where $p_{h}$ is the 4-momentum of the final-state hadron. In addition, the
squared invariant mass of the photon-nucleon ($\gamma^{\ast}$-N) system is
$W^{2}=(q_{\gamma}+p_{N})^{2}$.

Figure 8: Kinematic ($x_{B}$ and $Q^{2}$) regions covered by three EIC
facilities shown as shaded area. Boundaries for US-EIC, EicC, and JLab are
plotted with dashed, dotted-dashed, and dotted curves, respectively. Discrete
circles are sampled momentum fraction $x$ and probing scale $Q^{2}$ from
numerical calculations in current analysis of $\hat{q}$ (Density of the
circles doesn’t represent the density of data points).

In our calculations, the center-of-mass energies ($\sqrt{s}$) of the electron-
nucleon system for the three EIC facilities are taken to be 90 GeV (US-EIC),
10.6 GeV (EicC), and 4.8 GeV (JLab), respectively [1, 57, 56, 95]. The
kinematic range considered in our calculations is: $Q^{2}>1$ GeV2,
$0.01<y<0.95$, $W^{2}>10$ GeV2 ($>4$ GeV2 for JLab), and $0.3<z_{h}<0.8$ [95,
96]. With this kinematic restriction, we plot in Fig. 8 the ranges of Bjorken
$x$ and $Q^{2}$ covered by the three facilities. For comparison, some sampled
values of the $x$ and $Q^{2}$ involved in our global analysis of the current
data are also shown. We can see that the future EIC facilities have the
potential to allow a high-coverage scan on the kinematic dependence of
$\hat{q}$, especially for small- and large-$x$ regions where the current
measurements rarely access.

Besides the wide kinematic coverage, the high precision measurements at future
EIC are expected to provide more powerful constraints on the $\hat{q}$. In
this study, to preliminarily estimate the uncertainty of the measurement in
future EIC, we consider an integrated luminosity as $\mathcal{L}=5$ fb-1 [57,
56, 95], and evaluate the relative statistical uncertainty as
$\delta_{st}=1/\sqrt{\sigma\mathcal{L}}$, where $\sigma$ is the cross section
of the considered process. Due to the lack of the information on the
systematic uncertainty, we simply assume it is on the same order of the
statistical uncertainty, and include an additional factor $\sqrt{2}$ to
estimate the total uncertainty. Now we present the study for three types of
nuclear induced broadening/imbalance as follows.

###  III.1 Single hadron $p_{T}$ broadening in SIDIS

Figure 9: Transverse momentum broadening $\Delta\langle p_{T}^{2}\rangle$ of
single pion production in SIDIS as a function of Bjorken $x_{B}$ at $Q^{2}=$2
and 10 GeV2, for three EIC facilities. Solid curve with shaded area represents
results with $\hat{q}(x,Q^{2})$ with uncertainties. Dashed curve with dotted
band shows results with constant $\hat{q}$ with uncertainties. Circles with
vertical bars represent estimated experimental uncertainties. In calculations,
we have taken $Q^{2}\in[1.5,2.5]$ GeV2 (left panel) and $Q^{2}\in[9,11]$ GeV2
(right panel), and set $x_{B}<0.6$. On top of each panel, we mark the
kinematic ($x_{B}$) ranges covered by three facilities, and overlaps among
them can be seen. The theoretical results for two facilities in their overlap
region is generally similar, since they depend on $x_{B}$ at $Q^{2}$ to a
large extent.  Figure 10: Similar as Fig. 9, but for transverse momentum
broadening $\Delta\langle p_{T}^{2}\rangle$ of single pion production in SIDIS
as a function of $Q^{2}$ within $0.07<x_{B}<0.09$. Results are calculated for
$1<Q^{2}<100$ GeV2.

The transverse momentum broadening of the single hadron production in SIDIS is
a type of observable that has already played an important role in our current
analysis of $\hat{q}$, and will be still important in future EIC. At lowest
order in QCD, the single hadron comes from the fragmentation of the nuclear
struck quark, which can rescatter with the nuclear medium when traversing it.
The leading-twist cross section for the single scattering process can be
written as

$\displaystyle\frac{d\sigma^{S}}{dx_{B}dQ^{2}dz_{h}}\\!=$
$\displaystyle\\!\frac{2\pi\alpha^{2}_{\textrm{em}}}{Q^{4}}\left[1\\!+\\!(1-y)^{2}\right]$
$\displaystyle\times\sum_{q}e_{q}^{2}f_{q/A}(x_{B},\mu^{2})D_{h/q}(z_{h},\mu^{2})\,.$
(17)

The $p_{T}$ broadening in Eq. (3), at given values of $x_{B}$, $Q^{2}$ and
$z_{h}$, can be expressed as [44, 45]

$\displaystyle\Delta\langle
p_{T}^{2}\rangle\\!=\\!\left(\\!\frac{8\pi^{2}\alpha_{s}z_{h}^{2}C_{F}}{N^{2}_{c}-1}\\!\right)\\!\frac{\sum_{q}e_{q}^{2}T_{qg}(x_{B},\\!0,\\!0,\\!\mu^{2})D_{h\\!/\\!q}(z_{h},\mu^{2})}{\sum_{q}e_{q}^{2}f_{q/A}(x_{B},\mu^{2})D_{h\\!/\\!q}(z_{h},\mu^{2})}\,,$
(18)

where the color factor $C_{F}$ corresponds to the transport of a quark jet. In
our calculation, the factorization scale in Eq. (18) is taken to be
$\mu^{2}=Q^{2}$, and $\Delta\langle p_{T}^{2}\rangle$ is obtained by averaging
over a studied kinematic bin. The theoretical inputs, i.e., the PDFs and FFs,
are the same as in our global analysis.

In Fig. 9, we plot the results for the $p_{T}$ broadening of pion production
in SIDIS in electron-lead ($e$-Pb) scattering at three EIC facilities, at
$Q^{2}\\!=\\!2$ GeV2 (left panel) and $10$ GeV2 (right panel), respectively.
The solid curve with shaded band represents the predictions with the kinematic
dependent $\hat{q}(x,Q^{2})$ along with the uncertainties extracted in section
II.3. We find that the $\Delta\langle p_{T}^{2}\rangle$ as a function of
$x_{B}$ clearly reflects the $x$ dependence and the uncertainties of
$\hat{q}(x,Q^{2})$, e.g., the growths of both the broadening and the
uncertainties in small- and large-$x$ regions. For comparison, we also show
the predictions with the kinematics independent constant $\hat{q}$ extracted
from the current data as the dashed curve and dotted uncertainty band. As
expected, the differences between the two theoretical predictions appear
mainly in small- and large-$x$ regions. In the bottom of each panel, we show
the estimated experimental uncertainties as a reference, which are rather
small compared to the theoretical uncertainties. Therefore, we expect the
future EIC experiments should be able to distinguish the two theoretical
predictions and provide powerful constraints on the $\hat{q}$.

Fig. 10 shows the $\Delta\langle p_{T}^{2}\rangle$ as a function of the
probing scale $Q^{2}$ in an intermediate-$x_{B}$ region where the difference
between the extracted $\hat{q}(x,Q^{2})$ and constant $\hat{q}$ is small.
Since the scale dependence of $\hat{q}(x,Q^{2})$ is mild, we only see small
differences between two theoretical predictions in the studied $Q^{2}$
regions. The broadening with the constant $\hat{q}$ slightly decreases with
increasing $Q^{2}$, which is due to the decreasing averaged $z_{h}^{2}$ [a
factor in Eq. (18)] as a result of the scale evolution of the fragmentation
function. Similar as in Fig. 9, the estimated experimental uncertainties are
small for the plotted $Q^{2}$ region.

###  III.2 Nuclear enhancement of di-hadron transverse momentum imbalance

Now we focus on the nuclear-medium enhanced transverse momentum imbalance of
the back-to-back particles production in future EIC. The transverse momentum
imbalance of the particle pair given by
$\vec{p}_{T}=\vec{p}_{1T}+\vec{p}_{2T}$ is equal to the total transverse
momentum of the two particles [49]. Thus, the enhancement of the imbalance is
also the broadening of the total transverse momentum. In this subsection, we
study the imbalance of di-hadron (pion pair) production. The back-to-back
hadron pair comes from the fragmentation of di-jet, which can be produced
through the processes $\gamma^{\ast}q\rightarrow qg$ and
$\gamma^{\ast}g\rightarrow q\bar{q}$ at LO in $\alpha_{s}$.

Figure 11: Similar as Fig. 9, but for nuclear enhancement of transverse
momentum imbalance $\Delta\langle p_{T}^{2}\rangle$ of di-hadron ($\pi\pi$)
production as a function of $x_{B}$. Result for each bin is calculated in
$\gamma^{\ast}$-N system with fixed $Q^{2}$ and $x_{B}$ (central value). Green
square marks averaged momentum fraction $x=x_{B}(1+M_{JJ}^{2}/Q^{2})$ for
corresponding $x_{B}$, with values shown on the right-hand-side vertical axis.
In right panel ($Q^{2}=10$ GeV), dihadron production in JLab is kinematically
prevented/suppressed, thus is not shown.

Following an earlier work [49], the differential cross section of di-hadron
production can be written in the center of mass frame of the virtual-photon-
nucleon ($\gamma^{\ast}$-N) system as

$\displaystyle\frac{d\sigma^{S}}{dy_{1}dy_{2}dp_{1T}dp_{2T}}\\!=\\!\frac{2\pi\alpha_{s}\alpha_{\textrm{em}}}{(W^{2}+Q^{2})^{2}}\\!\sum_{b,c,d}\\!D_{h_{1}/c}(z_{1})D_{h_{2}/d}(z_{2})$
(19) $\displaystyle\times\frac{f_{b/A}(x)}{x}H^{U}_{\gamma^{\ast}b\rightarrow
cd}\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (20)

where we have suppressed the scale $\mu^{2}$-dependence in the PDFs and FFs,
$y_{1(2)}$, $p_{1T(2T)}$ and $z_{1(2)}$ are the rapidity, transverse momentum
and momentum fraction of the hadron $h_{1(2)}$, and
$H^{U}_{\gamma^{\ast}b\rightarrow cd}$ represents the perturbatively
calculable hard function of the partonic subprocess $\gamma^{\ast}b\rightarrow
cd$ [49]. The nuclear enhancement of the transverse momentum imbalance of the
di-hadron, induced by the multiple scattering undergone by the two outgoing
partons $c$ and $d$, is only sensitive to the total color of the two-parton
composite state, which is equal to the color of the initial-state nuclear
parton $b$. Accordingly, the nuclear enhancement at given values of $y_{1(2)}$
and $p_{1T(2T)}$ is expressed as [49]

$\displaystyle\Delta\langle
p_{T}^{2}\rangle=\left(\\!\frac{8\pi^{2}\alpha_{s}}{N^{2}_{c}-1}\\!\right)\frac{\sum_{b,c,d}D_{h_{1}/c}(z_{1},\mu^{2})D_{h_{2}/d}(z_{2},\mu^{2})\frac{1}{x}T_{bg}(x,0,0,\mu^{2})H^{F}_{\gamma^{\ast}b\rightarrow
cd}}{\sum_{b,c,d}D_{h_{1}/c}(z_{1},\mu^{2})D_{h_{2}/d}(z_{2},\mu^{2})\frac{1}{x}f_{b/A}(x,\mu^{2})H^{U}_{\gamma^{\ast}b\rightarrow
cd}}\,,$ (21)

with the hard function $H^{F}_{\gamma^{\ast}b\rightarrow cd}$ written as

$\displaystyle H^{F}_{\gamma^{\ast}b\rightarrow
cd}=\begin{cases}C_{F}H^{U}_{\gamma^{\ast}b\rightarrow
cd}~{}~{}~{}b=\textrm{quark}\\\ C_{A}H^{U}_{\gamma^{\ast}b\rightarrow
cd}~{}~{}~{}b=\textrm{gluon.}\end{cases}$ (22)

Here the color factor $C_{F}$ ($C_{A}$) corresponds to the process initiated
with a nuclear quark (gluon). This is more complicated than the case in SIDIS,
where only the transport of a quark jet is involved at LO in $\alpha_{s}$. In
our calculation, the renormalization and factorization scales are taken to be
the averaged transverse momentum of di-hadron, i.e., $\mu=(p_{1T}+p_{2T})/2$.
We also employ kinematic cuts as $1<p_{1T(2T)}<4$ GeV and $0.1<y_{1(2)}<2.0$.

In Fig. 11, we plot the results for the nuclear enhancement of the di-hadron
(pion pair) imbalance versus Bjorken $x_{B}$, at $Q^{2}\\!=\\!2$ GeV2 (left
panel) and $10$ GeV2 (right panel), respectively. The predictions with both
the kinematics dependent $\hat{q}(x,Q^{2})$ and constant $\hat{q}$, as well as
the estimated experimental uncertainties are shown. Note that the momentum
fraction $x$ carried by the initial-state nuclear parton $b$ can be expressed
as

$\displaystyle
x=\frac{Q^{2}+M_{JJ}^{2}}{2p_{N}\cdot\\!q_{\gamma}}=x_{B}\left(1+\frac{M_{JJ}^{2}}{Q^{2}}\right),$
(23)

where $M_{JJ}$ is the invariant mass of the outgoing two partons (di-jet) $c$
and $d$. Therefore, the $x$ that enters $\hat{q}(x,Q^{2})$ is larger than the
Bjorken variable $x_{B}$. To show this difference, we mark in Fig. 11 the
averaged momentum fraction $\langle x\rangle$ for the corresponding $x_{B}$ as
green square. For a same region of Bjorken $x_{B}$, the di-hadron production
probes the $\hat{q}(x,Q^{2})$ at larger values of $x$ compared to the single-
hadron production.

Figure 12: Similar as Fig. 11, but for nuclear enhancement of transverse
momentum imbalance $\Delta\langle p_{T}^{2}\rangle$ of heavy-meson pair
($D\bar{D}$) production as a function of $x_{B}$.

Since there is a mixture of the processes with color factors $C_{A}$ ad
$C_{F}$, the di-hadron production is expected to suffer stronger multiple
scattering in nuclear medium than the single-hadron production. At the same
time, we can observe the broadening with a constant $\hat{q}$ become weaker
for a greater value of $x_{B}$ (or $x$) due to the decreasing contributions
from the gluon initiated processes. Since our calculation is made in the
$\gamma^{\ast}$-N system, we have simply estimated the experimental
uncertainties with a same luminosity as in $e$-A system to provide a
reference. The measurement of di-hadron imbalance in future EIC is expected to
provide valuable constraints on the behavior of $\hat{q}(x,Q^{2})$ in
large-$x$ region, where the large uncertainties from the current analysis make
the predictions less reliable.

###  III.3 Nuclear enhancement of transverse momentum imbalance of heavy meson
pair

Similar to the di-hadron imbalance, there should be nuclear enhancement of the
transverse momentum imbalance in heavy meson pair production. In this work, we
study that for the $D\bar{D}$ pair from the fragmentation of the produced
back-to-back heavy quark pair $Q\bar{Q}$ ($c\bar{c}$). At lowest order in
$\alpha_{s}$, the $c\bar{c}$ di-jet is produced from the photon-gluon
scattering $\gamma^{\ast}g\rightarrow c\bar{c}$. The differential cross
section for $\gamma^{\ast}+A\rightarrow D(p_{1})+\bar{D}(p_{2})+X$ in the
center of mass frame of the $\gamma^{\ast}$-N system can be expressed in a
similar form of the di-hadron production in Eq. (20) as [49]

$\displaystyle\frac{d\sigma^{S}}{dy_{1}dy_{2}dp_{1T}dp_{2T}}=\frac{2\pi\alpha_{s}\alpha_{\textrm{em}}}{(W^{2}+Q^{2})^{2}}D_{D\\!/\\!Q}(z_{1})D_{\bar{D}\\!/\\!\bar{Q}}(z_{2})$
(24) (25)
$\displaystyle\times\frac{f_{g/A}(x)}{x}H^{U}_{\gamma^{\ast}g\rightarrow
Q\bar{Q}}\,.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (26)

The corresponding nuclear enhancement at given values of $y_{1(2)}$ and
$p_{1T(2T)}$ is

$\displaystyle\Delta\langle
p_{T}^{2}\rangle=\left(\\!\frac{8\pi^{2}\alpha_{s}C_{A}}{N^{2}_{c}-1}\\!\right)\frac{T_{gg}(x,0,0,\mu^{2})}{f_{g/A}(x,\mu^{2})}\,,$
(27)

This simple form is related to the fact that the $c\bar{c}$ di-jet is only
initiated with a nuclear gluon, corresponding to the color factor $C_{A}$. In
our calculation, the fragmentation from $c$ ($\bar{c}$) quark to $D$
($\bar{D}$) meson is described with KKKS08 fragmentation functions [68], and
the renormalization and factorization scales are taken to be the averaged
transverse mass of the $D\bar{D}$ pair, i.e., $\mu=(m_{1T}+m_{2T})/2$. We
employ the same kinematic cuts as in the calculation for di-hadron production.
The momentum fraction $x$ carried by the initial-state nuclear gluon can be
similarly given by Eq. (23). In heavy meson pair production, the averaged
momentum fraction $x$ for a certain value of Bjorken $x_{B}$ can be even
larger than that in di-hadron production, due to the mass of heavy quark.

In Fig. 12, we plot the results for the nuclear enhancement of the imbalance
of $D\bar{D}$ pair production versus the Bjorken $x_{B}$, at $Q^{2}\\!=\\!2$
GeV2 (left panel) and $10$ GeV2 (right panel), respectively. The averaged
$\langle x\rangle$ for each kinematic bin is marked in the plot. It can be
observed that, for a same $x$ value, the $\Delta\langle p_{T}^{2}\rangle$ in
$D\bar{D}$ production is stronger than that in di-hadron production, which is
expected as a result of the greater color factor, to be examined in future EIC
experiments.

The three types of hard probes in future EIC studied in this section are
expected to provide a multi-dimensional understanding of the jet transport
property of cold nuclear matter in a wide kinematic range. In addition,
although the results of $\Delta\langle p_{T}^{2}\rangle$ in this section are
obtained for the electron-lead collisions, the results for other colliding
nuclei can be easily obtained by rescaling with the nuclear radius ($\times
R_{A}/R_{Pb}$), according to Eq. (6).

###  III.4 Advantage of future EIC measurement for understanding kinematic
dependence of $\hat{q}$

Figure 13: Dependence of results on individual parameters, $\hat{q}_{0}$,
$\alpha$, $\beta$, and $\gamma$, represented with circle, diamond, triangle,
and square symbols, respectively. Top row: transverse momentum broadening
$\Delta\langle p_{T}^{2}\rangle$ of single pion production in SIDIS in certain
kinematic regions of three EIC facilities, where $Q^{2}=$10 GeV2 and
$x=x_{1}$, $x_{3}$ and $x_{4}$ for US-EIC, EicC and JLab, respectively. Bottom
row: ratios of $\Delta\langle p_{T}^{2}\rangle$ in different kinematic regions
defined as Eq. (31), where $Q^{2}=$10 GeV2 for left and middle panels, and
$x\in[0.08,0.1]$ for right panel. Blue solid line with grey band represents
theoretical prediction with uncertainty. Abscissa axis shows rescaled
parameter values as $(a_{i}-a_{i}^{0})/\delta a_{i}$, where $\delta a_{i}$ is
uncertainty of $a_{i}$. Yellow circle with vertical bar represents estimated
experimental uncertainty.

From the above study, we can see that the three future EIC facilities jointly
provide a nearly full-kinematics scan of the jet transport property with high
precisions, which will be valuable for determining the kinematic dependence of
$\hat{q}$. To illustrate this, we can note that, in the sense of the
parametrization of $\hat{q}(x,Q^{2})$ in Eq. (7), the kinematic dependence is
determined by the parameters $\alpha$, $\beta$, and $\gamma$, and we have the
ratio for two values of $x$

$\displaystyle
r(x_{1},x_{2})\equiv\frac{\hat{q}(x_{1},Q^{2})}{\hat{q}(x_{2},Q^{2})}=\left(\frac{x_{1}}{x_{2}}\right)^{\alpha}\left(\frac{1-x_{1}}{1-x_{2}}\right)^{\beta}.$
(28)

Assuming that both $\alpha$ and $\beta$ are negative as suggested by our
analysis, we can find that, for small $x$ values ($x_{1,2}\rightarrow 0$ and
$x_{1}<x_{2}$), we have $r(x_{1},x_{2})\rightarrow(x_{1}/x_{2})^{\alpha}$, and
for large $x$ values ($x_{1,2}\rightarrow 1$ and $x_{1}>x_{2}$), we have
$r(x_{1},x_{2})\rightarrow[(1-x_{1})/(1-x_{2})]^{\beta}$. Thus, the behaviors
of $\hat{q}$ in small and large $x$ regions are dominated by $\alpha$ and
$\beta$, respectively. Similarly, the ratio for two values of $Q^{2}$

$\displaystyle
r(Q^{2}_{1},Q^{2}_{2})\equiv\frac{\hat{q}(x,Q^{2}_{1})}{\hat{q}(x,Q^{2}_{2})}=\frac{\alpha_{s}(Q^{2}_{1})}{\alpha_{s}(Q^{2}_{2})}\left[\frac{\ln(Q^{2}_{1}/Q^{2}_{0})}{\ln(Q^{2}_{2}/Q^{2}_{0})}\right]^{\gamma}$
(29)

is only sensitive to the parameter $\gamma$, which dominates the scale
dependence. Accordingly, we can define a secondary observable as a ratio of
$\Delta\langle p_{T}^{2}\rangle$ in EIC, such as

$\displaystyle R_{\textrm{EIC}}(x_{1},x_{2})=\frac{\Delta\langle
p_{T}^{2}\rangle(x_{1},Q^{2})}{\Delta\langle
p_{T}^{2}\rangle(x_{2},Q^{2})}\,,$ (30) $\displaystyle
R_{\textrm{EIC}}(Q^{2}_{1},Q^{2}_{2})=\frac{\Delta\langle
p_{T}^{2}\rangle(x,Q^{2}_{1})}{\Delta\langle
p_{T}^{2}\rangle(x,Q^{2}_{2})}\,,$ (31)

to measure the $x$ and $Q^{2}$ dependence of $\hat{q}(x,Q^{2})$, respectively.
Figure 13 shows, for several observables in SIDIS, the dependence of
theoretical prediction on each parameters $a_{i}$ ($\hat{q}_{0}$, $\alpha$,
$\beta$, and $\gamma$), as a function of the relative parameter displacement
$(a_{i}-a_{i}^{0})/\delta a_{i}$, with $\delta a_{i}$ being the uncertainty of
$a_{i}$. The three panels in the top row show the dependence for
$\Delta\langle p_{T}^{2}\rangle$ at three facilities. We can see a single
measurement of $\Delta\langle p_{T}^{2}\rangle$ is usually sensitive to more
than one parameters. However, as shown in the bottom row, the three jointly
measured $R_{\textrm{EIC}}$ are separately sensitive to the parameters
$\alpha$, $\beta$, and $\gamma$. The future EIC facilities will allow a
precise understanding of the detail of the kinematic dependence of $\hat{q}$,
thanks to the wide kinematic coverage and high-precision measurement.

##  IV Summary and Discussion

To gain foreknowledge on how the future electron-ion-collision experiments can
deepen our understanding of the jet transport property of cold nuclear matter,
in this work, we study the nuclear-medium induced transverse momentum
broadening/imbalance of single-/pairwise-particle production in electron-ion
collisions, for the kinematic regions of three proposed future facilities. Our
theoretical calculations take into account the multiple scattering undergone
by the colored hard probe that traverses the nucleus, within the framework of
the higher-twist expansion, i.e., the generalized factorization in
perturbative QCD.

Particularly, a kinematic dependent jet transport coefficient
$\hat{q}=\hat{q}(x,Q^{2})$, extracted from our global analysis of the current
experimental data, is used in our calculations for EIC. This globally
extracted $\hat{q}(x,Q^{2})$, together with its uncertainty evaluated with a
Hessian matrix method, are available for the community to make theoretical
predictions. Moreover, by adding a new data set on $J/\psi$ production from
the LHC, the Hessian analysis results in reduced uncertainties of
$\hat{q}(x,Q^{2})$ in small-$x$ region. At the same time, we show that this
$\hat{q}(x,Q^{2})$ can be converted into the function of jet energy, i.e.,
$\hat{q}(E_{\textrm{jet}},Q^{2})$, which may be instructive for the study of
the medium modification of jet production in heavy-ion collisions.

The current analysis suggests enhancements of $\hat{q}$ in both small- and
large-$x$ regions, however, the uncertainties in these two regions are still
considerable due to the limited data points therein, which is expected to be
better constrained in future EIC experiments. With the extracted
$\hat{q}(x,Q^{2})$, we study three types of observable in EIC, including the
transverse momentum broadening of single hadron production and the enhancement
of two-particle transverse momentum imbalance in di-hadron and heavy-meson
pair productions. These nuclear induced broadening/imbalance are found to be
sensitive to the color state of the hard probe, and to exhibit a clear
kinematic dependence stemming from $\hat{q}(x,Q^{2})$. Besides, the results
with a constant $\hat{q}$ are also given for comparison. We find that the
future EIC experiments have great potential to provide precise understanding
of $\hat{q}$ in a wide kinematic range and to facilitate the jet tomography of
cold nuclear medium.

Note: The transport coefficient $\hat{q}(x,Q^{2})$ extracted from our global
analysis and the uncertainty set are available for user and can be requested
by email from<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
Acknowledgments.  This research was supported in part by the National Natural
Science Foundation of China (NSFC) under Grants No. 12022512 and No. 12035007,
by Guangdong Major Project of Basic and Applied Basic Research No.
2020B0301030008, by Guangdong Basic and Applied Basic Research Foundation
(Project No. 2022A1515110392, 2022A1515010683), by the National Science
Foundation in US under Grant No. PHY-1945471 (Z.K.), by the China Postdoctoral
Science Foundation under Project No. 2019M652929 (P.R.), and by the MOE Key
Laboratory of Quark and Lepton Physics (CCNU) under Project No. QLPL201802
(P.R.).

##  Appendix A Details of Hessian analysis

We briefly review the Hessian analysis employed in this work. Hessian matrix
analysis is a well-known technique for the uncertainty estimation in a global
analysis [97]. The general idea is to optimize the representation of the
parameter space in the neighborhood of the global minimum $\chi^{2}$, and to
provide an uncertainty set of the parameterized quantity [e.g.
$\hat{q}(x,Q^{2})$ in this work], from which uncertainties of all related
physical quantities can be evaluated.

A Hessian analysis begins with finding in parameter space $\\{a_{i}\\}$ the
coordinate corresponding to the minimum of the global $\chi^{2}$. The basic
assumption of the Hessian analysis is that the $\chi^{2}$ can be approximated
with a quadratic form in the neighborhood of the minimum [51, 52] as

$\displaystyle\chi^{2}=\chi_{0}^{2}+\sum_{i,j}H_{ij}y_{i}y_{j}\,,$ (32)

where $\chi_{0}^{2}\equiv\chi^{2}(\\{a_{i}^{0}\\})$ is the global minimum of
the $\chi^{2}$ at the optimal parameter values $\\{a_{i}\\}=\\{a_{i}^{0}\\}$,
$y_{i}=a_{i}-a_{i}^{0}$ is the displacement of $a_{i}$ from its optimal value
$a_{i}^{0}$, and $H_{ij}$ is the element of Hessian matrix defined as the
second-order partial derivatives of $\chi^{2}$ at the minimum

$\displaystyle H_{ij}=\frac{1}{2}\left(\frac{\partial^{2}\chi^{2}}{\partial
y_{i}\partial y_{j}}\right)_{a_{i}=a_{i}^{0}}$ (33)

Usually there are interplays among different variables $y_{i}$ (or $a_{i}$) in
the $\chi^{2}$, and the off-diagonal Hessian matrix elements could be non-
zero. This makes the uncertainty estimation, corresponding to a certain
tolerance $\Delta\chi^{2}\equiv\chi^{2}-\chi^{2}_{0}$, not that
straightforward.

To disentangle these interplays, one can define a new set of parameters
$\\{z_{i}\\}$, in whose representation the Hessian matrix is diagonal. This
can be achieved by using the complete set of $n$ orthonormal eigenvectors
$V_{i}^{(k)}$ of the symmetric $n\times n$ Hessian matrix ($n$ is the number
of parameters), which satisfy the characteristic equations

$\displaystyle\sum_{j}H_{ij}V_{j}^{(k)}=\lambda_{k}V_{i}^{(k)},$ (34)

with $\lambda_{k}$ being the positive eigenvalues. The new parameters
$\\{z_{i}\\}$ can be expressed as the linear combinations of the original
parameters $\\{y_{i}\\}$ [51, 52]

$\displaystyle z_{i}=\sqrt{\lambda_{i}}\sum_{j}y_{j}V_{j}^{(i)}.$ (35)

With the new set of parameters $\\{z_{i}\\}$, the $\Delta\chi^{2}$ can be
written in a simple form as

$\displaystyle\Delta\chi^{2}=\sum_{i}z_{i}^{2},$ (36)

which means that the contours of the $\chi^{2}$ are spheres in the new basis.
With Eq. (36) and the inverse transformation of Eq. (35), one can define the
uncertainty set $\\{a_{i}^{(\pm k)}\\}$ of the original parameters
$\\{a_{i}\\}$ corresponding to a tolerance $\Delta\chi^{2}$ as [52]

$\displaystyle a_{i}^{(\pm
k)}=a_{i}^{0}\pm\sqrt{\frac{\Delta\chi^{2}}{\lambda_{k}}}V_{i}^{(k)}\,,~{}~{}\textrm{for}~{}~{}k=1,2,\ldots,n.$
(37)

For a quantity $\mathcal{Q}$ as a function of $\\{a_{i}\\}$, whose value
corresponding to $\\{a_{i}^{(\pm k)}\\}$ is $\mathcal{Q}_{\pm
k}\equiv\mathcal{Q}(\\{a_{i}^{(\pm k)}\\})$, its uncertainty can be evaluated
with

$\displaystyle\Delta\mathcal{Q}=\frac{1}{2}\sqrt{\sum^{n}_{k=1}(\mathcal{Q}_{+k}-\mathcal{Q}_{-k})^{2}}.$
(38)

As an example, one can give the uncertainty set of $\hat{q}$

$\displaystyle S_{\pm k}\equiv{\hat{q}}_{\pm k}={\hat{q}}(\\{a_{i}^{(\pm
k)}\\}),$ (39)

and express the uncertainty set for any quantity as a function of $\hat{q}$ as
$\mathcal{Q}_{\pm k}=\mathcal{Q}(S_{\pm k})$.

Figure 14: $\hat{q}(x,Q^{2})$ and uncertainties extracted from global analysis
with and without ALICE data (2021) on transverse momentum broadening of
$J/\psi$ production. Results of new analysis are the same as in Fig. 4.
Results without $J/\psi$ (21’) data are shown with a dashed curve with dotted
boundaries.

In our analysis, the tolerance of the global $\chi^{2}$ is given at a $p\%$
confidence level (C.L.) by simply re-scaling the minimum as [91, 52]

$\displaystyle\Delta\chi^{2}=\chi^{2}_{0}\left(\frac{\xi_{p}}{\xi_{50}}-1\right)\,,$
(40)

where the rescaling parameter $\xi_{p}$ corresponds to the $p$-th percentile
of the $\chi^{2}$ distribution $P(\chi^{2},N)$ and is determined by

$\displaystyle\int_{0}^{\xi_{p}}P(\chi^{2},N)\,d\chi^{2}=p\,\%\,,$ (41)

with $N$ the number of data points and $P(\chi^{2},N)$ defined as

$\displaystyle
P(\chi^{2},N)=\frac{(\chi^{2})^{N/2-1}e^{-\chi^{2}/2}}{2^{N/2}\Gamma(N/2)}.$
(42)

A more thorough and complicated scheme to estimate the tolerance can be found
in Refs. [91, 52]. In our results, we have estimated the uncertainty at $90$%
C.L., corresponding to the tolerance $\Delta\chi^{2}=35$. To get the
uncertainty at any other C.L., one can simply rescale the $\Delta\chi^{2}$ to
get the corresponding uncertainty set $\\{a_{i}\\}$ with Eq. (37).

The Hessian analysis in this work is performed by using the MINUIT package
combined with the ITERATE program [51], which will first search for a global
minimum of the $\chi^{2}$ and then calculate the eigenvectors $V_{i}^{(k)}$
and eigenvalues $\lambda_{k}$ of the Hessian matrix with an iterative method
[51].

Figure 15: Comparison between theoretical results and experimental
measurements for part of the representative observables in the global analysis
of $\hat{q}$, including transverse momentum broadening $\Delta\langle
p_{T}^{2}\rangle$ in SIDIS ($e$A), Drell-Yan process ($p$A) and heavy
quarkonium production ($p$A), and nuclear modification of structure function
$F_{2}(A)/F_{2}(D)$. Theoretical results and uncertainties calculated with
$\hat{q}(x,Q^{2})$ extracted in new analysis with ALICE $J/\psi$ (2021) data
added are plotted with solid curves and shaded areas, while those calculated
with $\hat{q}(x,Q^{2})$ extracted by excluding new $J/\psi$ (2021) data are
plotted with dashed curves and dotted boundaries. Circles represent
experimental data taken from Refs. [6, 77, 76, 78, 80, 79, 7, 8, 53, 81, 82].

##  Appendix B Impact of new ALICE $J/\psi$ data on $\hat{q}$

In comparison with our previous analysis [32], we have added a new data set,
i.e., the transverse momentum broadening of $J/\psi$ in $p$-Pb collisions at
$\sqrt{s_{NN}}=8.16$ TeV at the LHC (ALICE 21’), into the Hessian analysis in
this work. To quantitatively show the impact of the new data, we compare the
results with and without including this data set.

In Fig. 14, we compare the $\hat{q}(x,Q^{2})$ extracted with and without the
new $J/\psi$ data. We also compare theoretical predictions calculated by using
the $\hat{q}(x,Q^{2})$ extracted in the two analyses in Fig. 15. In general,
the new $J/\psi$ data has small impact on the optimal/central values of the
$\hat{q}(x,Q^{2})$ or theoretical predictions, indicating the new analysis is
consistent with our previous work [32].

With the new $J/\psi$ data, the uncertainties of $\hat{q}(x,Q^{2})$ in small
$x$ region are reduced, especially for $Q^{2}=10-100$ GeV2 (squared mass of
$J/\psi$ $\sim 10$ GeV2), as shown in Fig. 14. As a result, we can see in Fig.
15 that, with the $\hat{q}(x,Q^{2})$ extracted in the new analysis, the
uncertainties of the theoretical predictions for $J/\psi$ production in the
forward region at the LHC are obviously reduced (2nd panel in 3rd row and last
panel in 4th row). The theoretical uncertainties for other observables,
related to mid- and large-$x$ regions, have negligible changes.

## References

  * [1] A. Accardi, J. L. Albacete, M. Anselmino, N. Armesto, E. C. Aschenauer, A. Bacchetta, D. Boer, W. K. Brooks, T. Burton and N. B. Chang, et al. Eur. Phys. J. A 52, no.9, 268 (2016),  arXiv:1212.1701 [nucl-ex].
  * [2] Y. Akiba, A. Angerami, H. Caines, A. Frawley, U. Heinz, B. Jacak, J. Jia, T. Lappi, W. Li and A. Majumder, et al. arXiv:1502.02730 [nucl-ex].
  * [3] J. C. Collins, D. E. Soper and G. F. Sterman,  Adv. Ser. Direct. High Energy Phys. 5, 1-91 (1989),  arXiv:hep-ph/0409313 [hep-ph].
  * [4] M. Gyulassy, I. Vitev, X. N. Wang and B. W. Zhang,  arXiv:nucl-th/0302077 [nucl-th].
  * [5] S. Cao and X. N. Wang,  Rept. Prog. Phys. 84, 024301 (2021),  arXiv:2002.04028 [hep-ph].
  * [6] HERMES, A. Airapetian et al.,  Phys. Lett. B684, 114 (2010), arXiv:0906.2478.
  * [7] PHENIX, A. Adare et al.,  Phys. Rev. C87, 034904 (2013), arXiv:1204.0777.
  * [8] ALICE, J. Adam et al.,  JHEP 11, 127 (2015), arXiv:1506.08808.
  * [9] J. Adams et al. [STAR],  Phys. Rev. Lett. 91, 172302 (2003),  arXiv:nucl-ex/0305015 [nucl-ex].
  * [10] S. S. Adler et al. [PHENIX],  Phys. Rev. Lett. 97, 052301 (2006),  arXiv:nucl-ex/0507004 [nucl-ex].
  * [11] G. Aad et al. [ATLAS],  Phys. Rev. Lett. 105, 252303 (2010),  arXiv:1011.6182 [hep-ex].
  * [12] S. Chatrchyan et al. [CMS],  Phys. Rev. C 84, 024906 (2011),  arXiv:1102.1957 [nucl-ex].
  * [13] J. Adam et al. [ALICE],  JHEP 09, 170 (2015),  arXiv:1506.03984 [nucl-ex].
  * [14] S. Chatrchyan et al. [CMS],  Phys. Lett. B 730, 243-263 (2014),  arXiv:1310.0878 [nucl-ex].
  * [15] A. M. Sirunyan et al. [CMS],  Phys. Rev. Lett. 121, no.24, 242301 (2018),  arXiv:1801.04895 [hep-ex].
  * [16] S. Acharya et al. [ALICE],  Phys. Lett. B 802, 135227 (2020),  arXiv:1905.02512 [nucl-ex].
  * [17] R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne, and D. Schiff,  Nucl. Phys. B484, 265 (1997), arXiv:hep-ph/9608322.
  * [18] X.-F. Chen, T. Hirano, E. Wang, X.-N. Wang, and H. Zhang,  Phys. Rev. C84, 034902 (2011), arXiv:1102.5614.
  * [19] A. Majumder and C. Shen,  Phys. Rev. Lett. 109, 202301 (2012), arXiv:1103.0809.
  * [20] K. M. Burke et al. [JET],  Phys. Rev. C 90, 014909 (2014),  arXiv:1312.5003 [nucl-th].
  * [21] F. C. Zhou, G. L. Ma and Y. G. Ma,  Eur. Phys. J. A 56, 70 (2020),  arXiv:1902.00729 [hep-ph].
  * [22] M. Xie, S. Y. Wei, G. Y. Qin and H. Z. Zhang,  Eur. Phys. J. C 79, 589 (2019),  arXiv:1901.04155 [hep-ph].
  * [23] G.-Y. Ma, W. Dai, B.-W. Zhang, and E.-K. Wang,  Eur. Phys. J. C79, 518 (2019), arXiv:1812.02033.
  * [24] L. Chen, G.-Y. Qin, S.-Y. Wei, B.-W. Xiao, and H.-Z. Zhang,  Phys. Lett. B773, 672 (2017), arXiv:1607.01932.
  * [25] C. Andres, N. Armesto, M. Luzum, C. A. Salgado, and P. Zurita,  Eur. Phys. J. C76, 475 (2016), arXiv:1606.04837.
  * [26] A. Kumar, A. Majumder and J. H. Weber,  arXiv:2010.14463 [hep-lat].
  * [27] S. Cao et al. [JETSCAPE],  Phys. Rev. C 104, 024905 (2021),  arXiv:2102.11337 [nucl-th].
  * [28] M. Xie, W. Ke, H. Zhang and X. N. Wang,  arXiv:2206.01340 [hep-ph].
  * [29] E. Bianchi, J. Elledge, A. Kumar, A. Majumder, G. Y. Qin and C. Shen,  arXiv:1702.00481 [nucl-th].
  * [30] S. Shi, J. Liao and M. Gyulassy,  Chin. Phys. C 43, no.4, 044101 (2019),  arXiv:1808.05461 [hep-ph].
  * [31] A. Kumar, A. Majumder and C. Shen,  Phys. Rev. C 101, 034908 (2020),  arXiv:1909.03178 [nucl-th].
  * [32] P. Ru, Z. B. Kang, E. Wang, H. Xing and B. W. Zhang,  Phys. Rev. D 103, L031901 (2021),  arXiv:1907.11808 [hep-ph].
  * [33] P. Ru, Z. B. Kang, E. Wang, H. Xing and B. W. Zhang,  Nucl. Phys. A 1005, 121798 (2021),  arXiv:2004.00027 [nucl-th].
  * [34] M. Arratia, Y. Song, F. Ringer and B. V. Jacak,  Phys. Rev. C 101, 065204 (2020),  arXiv:1912.05931 [nucl-ex].
  * [35] F. Arleo and C. J. Naïm,  JHEP 07, 220 (2020),  arXiv:2004.07188 [hep-ph].
  * [36] Y. Y. Zhang, G. Y. Qin and X. N. Wang,  Phys. Rev. D 100, 074031 (2019),  arXiv:1905.12699 [hep-ph].
  * [37] T. X. Bai and C. G. Duan,  Eur. Phys. J. Plus 136, 1181 (2021),  arXiv:2011.14350 [hep-ph].
  * [38] M. Gyulassy, P. M. Jacobs, J. Liao, S. Shi, X. N. Wang and F. Yuan,  Nucl. Phys. A 1005, 121938 (2021),  arXiv:2012.06151 [hep-ph].
  * [39] Y. Y. Zhang and X. N. Wang,  Phys. Rev. D 105, 034015 (2022),  arXiv:2104.04520 [hep-ph].
  * [40] J. Casalderrey-Solana and X.-N. Wang,  Phys. Rev. C77, 024902 (2008), arXiv:0705.1352.
  * [41] J.-P. Blaizot and Y. Mehtar-Tani,  Nucl. Phys. A929, 202 (2014), arXiv:1403.2323.
  * [42] E. Iancu,  JHEP 10, 095 (2014), arXiv:1403.1996.
  * [43] T. Liou, A. H. Mueller, and B. Wu,  Nucl. Phys. A916, 102 (2013), arXiv:1304.7677.
  * [44] Z.-B. Kang, E. Wang, X.-N. Wang, and H. Xing,  Phys. Rev. Lett. 112, 102001 (2014), arXiv:1310.6759.
  * [45] Z.-B. Kang, E. Wang, X.-N. Wang, and H. Xing,  Phys. Rev. D94, 114024 (2016), arXiv:1409.1315.
  * [46] Z.-B. Kang, J.-W. Qiu, X.-N. Wang, and H. Xing,  Phys. Rev. D94, 074038 (2016), arXiv:1605.07175.
  * [47] Z.-B. Kang and J.-W. Qiu,  Phys. Rev. D77, 114027 (2008), arXiv:0802.2904.
  * [48] Z.-B. Kang and J.-W. Qiu,  Phys. Lett. B721, 277 (2013), arXiv:1212.6541.
  * [49] H. Xing, Z. B. Kang, I. Vitev and E. Wang,  Phys. Rev. D 86, 094010 (2012),  arXiv:1206.1826 [hep-ph].
  * [50] M. Alrashed, D. Anderle, Z. B. Kang, J. Terry and H. Xing,  Phys. Rev. Lett. 129, 242001 (2022),  arXiv:2107.12401 [hep-ph].
  * [51] J. Pumplin, D. R. Stump, and W. K. Tung,  Phys. Rev. D65, 014011 (2001), arXiv:hep-ph/0008191.
  * [52] K. Kovarik, A. Kusina, T. Jezo, D. B. Clark, C. Keppel, F. Lyonnet, J. G. Morfin, F. I. Olness, J. F. Owens and I. Schienbein, et al. Phys. Rev. D 93, 085037 (2016),  arXiv:1509.00792 [hep-ph].
  * [53] S. Acharya et al. [ALICE],  JHEP 02, 002 (2021),  arXiv:2008.04806 [nucl-ex].
  * [54] H. T. Li and I. Vitev,  Phys. Rev. Lett. 126, 252001 (2021),  arXiv:2010.05912 [hep-ph].
  * [55] X. Chen, F. K. Guo, C. D. Roberts and R. Wang,  Few Body Syst. 61, 43 (2020),  arXiv:2008.00102 [hep-ph].
  * [56] D. P. Anderle, V. Bertone, X. Cao, L. Chang, N. Chang, G. Chen, X. Chen, Z. Chen, Z. Cui and L. Dai, et al. Front. Phys. (Beijing) 16, 64701 (2021),  arXiv:2102.09222 [nucl-ex].
  * [57] V. D. Burkert,  Ann. Rev. Nucl. Part. Sci. 68, 405-428 (2018)
  * [58] J.-w. Qiu and G. F. Sterman,  Nucl. Phys. B353, 105 (1991).
  * [59] J.-w. Qiu and G. F. Sterman,  Nucl. Phys. B353, 137 (1991).
  * [60] M. Luo, J.-w. Qiu, and G. F. Sterman,  Phys. Rev. D50, 1951 (1994).
  * [61] M. Luo, J.-w. Qiu, and G. F. Sterman,  Phys. Lett. B279, 377 (1992).
  * [62] M. Luo, J.-w. Qiu, and G. F. Sterman,  Phys. Rev. D49, 4493 (1994).
  * [63] E. Wang and X. N. Wang,  Phys. Rev. Lett. 89, 162301 (2002),  arXiv:hep-ph/0202105 [hep-ph].
  * [64] B. W. Zhang, E. Wang and X. N. Wang,  Phys. Rev. Lett. 93, 072301 (2004),  arXiv:nucl-th/0309040 [nucl-th].
  * [65] Z. B. Kang, I. Vitev and H. Xing,  Phys. Rev. D 85, 054024 (2012),  arXiv:1112.6021 [hep-ph].
  * [66] J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky and W. K. Tung,  JHEP 07, 012 (2002),  arXiv:hep-ph/0201195 [hep-ph].
  * [67] D. de Florian, R. Sassot, and M. Stratmann,  Phys. Rev. D75, 114010 (2007), arXiv:hep-ph/0703242.
  * [68] T. Kneesch, B. A. Kniehl, G. Kramer and I. Schienbein,  Nucl. Phys. B 799, 34-59 (2008),  arXiv:0712.0481 [hep-ph].
  * [69] K. J. Eskola, P. Paakkinen, H. Paukkunen and C. A. Salgado,  Eur. Phys. J. C 77, no.3, 163 (2017),  arXiv:1612.05741 [hep-ph].
  * [70] R. Abdul Khalek, J. J. Ethier, J. Rojo and G. van Weelden,  JHEP 09, 183 (2020),  arXiv:2006.14629 [hep-ph].
  * [71] X.-f. Guo,  Phys. Rev. D58, 114033 (1998), arXiv:hep-ph/9804234.
  * [72] J.-w. Qiu and I. Vitev,  Phys. Rev. Lett. 93, 262301 (2004), arXiv:hep-ph/0309094.
  * [73] J.-w. Qiu and I. Vitev,  Phys. Lett. B632, 507 (2006), arXiv:hep-ph/0405068.
  * [74] J.-w. Qiu and G. F. Sterman,  Int. J. Mod. Phys. E12, 149 (2003), arXiv:hep-ph/0111002.
  * [75] J.-W. Qiu,  Eur. Phys. J. C43, 239 (2005), arXiv:hep-ph/0507268.
  * [76] NA10, P. Bordalo et al.,  Phys. Lett. B193, 373 (1987).
  * [77] P. L. McGaughey, J. M. Moss, and J. C. Peng,  Ann. Rev. Nucl. Part. Sci. 49, 217 (1999), arXiv:hep-ph/9905409.
  * [78] D. M. Alde et al.,  Phys. Rev. Lett. 66, 2285 (1991).
  * [79] M. J. Leitch et al.,  Phys. Rev. D52, 4251 (1995).
  * [80] J.-C. Peng,  AIP Conf. Proc. 494, 503 (1999), arXiv:hep-ph/9912371.
  * [81] E665, M. R. Adams et al.,  Phys. Rev. Lett. 68, 3266 (1992).
  * [82] E665, M. R. Adams et al.,  Z. Phys. C67, 403 (1995), arXiv:hep-ex/9505006.
  * [83] K. J. Golec-Biernat and M. Wusthoff,  Phys. Rev. D59, 014017 (1998), arXiv:hep-ph/9807513.
  * [84] Y. L. Dokshitzer, G. Marchesini, and B. R. Webber,  Nucl. Phys. B469, 93 (1996), arXiv:hep-ph/9512336.
  * [85] S. J. Brodsky,  Dynamical higher twist and high x phenomena: A Window to quark quark correlations in QCD,  in Workshop on Nucleon Structure in High x-Bjorken Region (HiX2000) Philadelphia, Pennsylvania, March 30-April 1, 2000, 2000, arXiv:hep-ph/0006310.
  * [86] L. Gamberg, Z. B. Kang, I. Vitev and H. Xing,  Phys. Lett. B 743, 112-120 (2015),  arXiv:1412.3401 [hep-ph].
  * [87] V. M. Braun, A. Vladimirov and J. H. Zhang,  Phys. Rev. D 99, 014013 (2019),  arXiv:1810.00048 [hep-ph].
  * [88] J. Collins, Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 32, 1-624 (2011).
  * [89] S. Dulat et al.,  Phys. Rev. D93, 033006 (2016), arXiv:1506.07443.
  * [90] W. K. Brooks and J. A. López,  Phys. Lett. B 816, 136171 (2021),  arXiv:2004.07236 [hep-ph].
  * [91] D. Stump et al.,  Phys. Rev. D65, 014012 (2001), arXiv:hep-ph/0101051.
  * [92] F. James and M. Roos,  Comput. Phys. Commun. 10, 343 (1975).
  * [93] E. Iancu and R. Venugopalan,  The Color glass condensate and high-energy scattering in QCD,  in Quark-gluon plasma 4, edited by R. C. Hwa and X.-N. Wang, pp. 249–3363, 2003, arXiv:hep-ph/0303204.
  * [94] P. Ru, Z. B. Kang, E. Wang, H. Xing and B. W. Zhang, PoS HardProbes2020, 152 (2021).
  * [95] E. C. Aschenauer, S. Fazio, M. A. C. Lamont, H. Paukkunen and P. Zurita,  Phys. Rev. D 96, 114005 (2017),  arXiv:1708.05654 [nucl-ex].
  * [96] E. C. Aschenauer, I. Borsa, R. Sassot and C. Van Hulse,  Phys. Rev. D 99, 094004 (2019),  arXiv:1902.10663 [hep-ph].
  * [97] J. Pumplin et al.,  Phys. Rev. D 65, 014013 (2001),  arXiv:hep-ph/0101032 [hep-ph].