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

Exploring the properties of nuclear matter at the femtometer scale across different phases—including cold nuclei, hadron gas, and hot/dense quark-gluon plasma—is a central goal of high-energy nuclear collision experiments, such as lepton–, hadron–, and nucleus–nucleus collisions [1, 2]. In these collisions, perturbatively calculable observables provide unique and controlled probes of the non-perturbative structure of the nuclear medium across a wide range of phenomena [3, 4]. Among these observables, those arising from the multiple scattering of partons propagating through the nuclear medium are of particular interest [5−11], because they encode the medium’s transport properties, typically characterized by the jet transport coefficient $ \hat q $ [12−16]. This parameter quantifies the transverse momentum broadening of a (quark) jet per unit path length in the medium and reflects the gluon density of the medium as ''seen'' by the jet. Thus, the observation that $ \hat q $ for the hot dense medium created in relativistic heavy-ion collisions is orders of magnitude larger than that for cold nuclear matter serves as a signature of the quark-gluon plasma [17−19].

In relativistic nucleus-nucleus (AA) collisions, the transport coefficient $ \hat q $ plays a central role in jet tomography of the quark–gluon plasma. In particular, its dependence on the medium temperature and on kinematic variables such as jet energy and the probing scale [20−31] offers valuable insights into the structure of the medium, a topic that has recently attracted considerable attention and motivated related theoretical developments [32−44]. Regarding the kinematic dependence of $ \hat q $, electron-nucleus (eA) and proton-nucleus (pA) collisions are of particular importance [22−26, 45]. They provide a cleaner environment for detailed studies of the transport properties of cold nuclear matter and for precision tests of theoretical descriptions of jet–medium interactions, which can in turn inform and guide analogous studies in AA collisions.

In previous work [22], we performed the first global extraction of $ \hat q $ in nuclei, based on existing data from eA and pA collisions and within the high-twist (HT) expansion framework, i.e., a generalized factorization in perturbative QCD [46, 47]. The data revealed a non-trivial kinematic dependence of the transport property, and theoretical calculations incorporating $ \hat q $ as a function of both the parton momentum fraction x and the probing scale $ Q^2 $ achieved overall good agreement with world data [22]. However, because most current data are concentrated in the intermediate x and $ Q^2 $ region, the extracted $ \hat q $ remains weakly constrained over a wide kinematic range. The considerable uncertainties in these regions motivate future measurements to reduce them and enable rigorous tests of the universality and kinematic dependence of $ \hat q $.

Upcoming electron–ion collision (EIC) experiments—including those at the U.S. Electron-Ion Collider (US-EIC) [1], the Electron-ion collider in China (EicC) [48, 49], and the Jefferson Lab 12 GeV program (JLab) [50]—will offer an excellent opportunity to address these issues, owing to their broad kinematic coverage and high-precision measurements. In these experiments, nuclear-induced transverse momentum broadening (TMB) in particle production serves as an observable directly linked to the jet transport coefficient $ \hat q $. Measurements of such observables can provide an intuitive picture of the kinematic dependence of $ \hat q(x,Q^2) $ and are therefore of significant interest.

In this work, we study the transverse momentum broadening of single-particle and back-to-back pair production in future electron–ion collisions. Our calculations are based on the high-twist expansion [51−53], which has demonstrated broad applicability in describing nuclear modifications in both cold and hot/dense nuclear matter [13, 22, 54−56]. A central input to these calculations is the jet transport coefficient $ \hat q $, considered both with and without kinematic dependence, as extracted from global analyses performed within the same HT framework. In this work, we further improve the global analysis by incorporating a Hessian method to systematically propagate uncertainties. By comparing the theoretical uncertainties arising from $ \hat q $ with the estimated experimental uncertainties, we demonstrate the potential of future EIC measurements to provide strong constraints on $ \hat q $ over a wide kinematic range.

The remainder of this paper is organized as follows. In Sec. II, we present the main results of our global analysis of $ \hat q $, enhanced by Hessian-based uncertainty quantification. The implications of the results for a possible jet-energy dependence of $ \hat q $ are also discussed. In Sec. III, we investigate three types of nuclear-induced transverse momentum broadening (Secs. III.A–C) across three proposed EIC facilities: the US-EIC, EicC, and JLab. The potential of future measurements to reveal the kinematic dependence of $ \hat q $ is discussed. Finally, a summary is provided in Sec. IV.

In this section, we briefly introduce the framework of our global analysis of $ \hat q $ based on existing data from eA and pA collisions and present the key improvements and main results of the analysis.

In this section, we briefly introduce the framework of our global analysis of $ \hat q $ based on existing data from eA and pA collisions and present the key improvements and main results of the analysis.

Our global extraction of $ \hat q $ is performed by analyzing world data from eA and pA collisions, primarily for the TMB of particle production, within the higher-twist expansion framework of perturbative QCD [22]. In this formalism, the effects of parton multiple scattering in the nuclear medium (target) are incorporated through twist-4 corrections [36, 51]. Taking TMB as an example, it is defined as the difference in the average transverse momentum squared, $ \langle p_T^2\rangle $, of produced particles between eA (pA) and $ ep $ ($ pp $) collisions:

Our calculation of TMB includes the leading-twist single-scattering cross section and the twist-4 double-scattering contribution enhanced by the large nuclear radius. Using pA collisions as an example, at twist-4 and leading order (LO) in $ \alpha_s $, the TMB can be schematically expressed as [36, 51]

where the superscripts ''S'' and ''D'' denote single- and double-scattering processes, respectively. Here, $ \mathrm{d}\sigma $ and $ \mathrm{d}\hat{\sigma} $ represent the differential cross sections at the observable level and for the perturbatively calculable partonic process, respectively. The terms f and $ D_h $ denote the relevant parton distribution functions (PDFs) and fragmentation functions (FFs), with the subscripts $ {q(g)/p} $ and $ {q(g)/A} $ indicating an incoming quark (or gluon) from the proton and nucleus, respectively. In the $ p_T^2 $-weighted double-scattering cross section ${\mathrm{ d}}\langle p_T^2\sigma^D\rangle $, $ T_{ij} $ is a twist-4 parton–parton correlation function that characterizes the relevant non-perturbative structure of the nuclear target in these processes.

At leading order in perturbation theory, the relevant jet transport coefficient in nuclei can be related to the twist-4 parton–parton correlation functions $ T_{ij} $. For example, the quark–gluon correlation function $ T_{qg} $ is usually approximated as [22, 36]

where x is the momentum fraction carried by the parton from the nucleus, μ is the factorization scale, $ R_A $ is the nuclear radius for a nucleus with mass number A, $ f_{q/A} $ is the parton distribution function of a quark q in the nucleus, and $ \hat q(x,\mu^2) $ is the effective nuclear-geometry-averaged jet transport coefficient that we aim to extract. Similarly, for the gluon–gluon correlation function $ T_{gg} $, which is relevant to processes initiated by a gluon in the nucleus, we adopt the same form as in Eq. (3), with $ f_{q/A} $ replaced by the gluon distribution $ f_{g/A} $. Throughout this paper, $ \hat q $ denotes the effective transport coefficient for a quark jet.

With this formalism, a global analysis of $ \hat q(x,Q^2) $ can be performed. It is important to note that the nonperturbative quantities directly constrained in such an analysis are the twist-4 correlation functions $ T_{ij} $. The quantity $ \hat q $ extracted in this manner corresponds, at leading order, to the transverse momentum broadening of a quark jet per unit propagation length, as shown in Eq. (2). Moreover, the relation between $ \hat q $ and $ T_{ij} $ may not hold exactly when higher-order contributions are included [37−43], although recent next-to-leading-order studies suggest that these corrections may be small [57].

On the other hand, a close connection exists between the twist-4 approach and formalisms that resum an arbitrary number of independent rescatterings, such as the Baier–Dokshitzer–Mueller–Peigne–Schiff (BDMPS) fra-mework [12] (see Ref. [12] for a general discussion). The correspondence between the HT formalism and color glass condensate effective theory has also been demonstrated recently [56, 58]. Nevertheless, it is important to emphasize that the interpretation of the variables x and $ Q^2 $ entering $ \hat q $ may differ across theoretical treatments. For example, in the BDMPS approach, x typically denotes the momentum fraction of the rescattered gluon in the nucleus, and $ Q^2 $ represents the scale of the rescattering process [12, 24]. These variables are therefore defined in a way that is more directly related to the intrinsic transport properties of the medium. In contrast, the $ \hat q(x,Q^2) $ factorized from the $ T_{ij} $ via Eq. (3) directly inherits the variables of the twist-4 correlation functions. In this context, x can be understood as the total momentum fraction carried by the two nuclear partons participating in double scattering. This fraction is fixed by the scattering kinematics (e.g., the Bjorken variable) and is related to the momentum fraction of the medium gluon by integration over the variables involved in the twist-4 matrix element $ T(x,x_1,x_2) $ [36, 51]. The scale $ Q^2 $ in $ T_{ij} $ serves as the factorization scale in QCD factorization and is typically taken to be the hard-scattering scale in practical calculations. Although a recent study [24] based on the BDMPS formalism has shown behavior of $ \hat q(x,Q^2) $ similar to our results, it is important to keep this difference in definition in mind when comparing the variable dependence between the two approaches. In addition, because x and $ Q^2 $ in our formalism are kinematically linked to the jet energy, the extracted $ \hat q(x,Q^2) $ naturally encodes the dependence of $ \hat q $ on the jet energy, as will be further discussed in Sec. II.C.

In the global analysis, to properly account for the functional form and possible kinematic dependence of $ \hat q $, we adopt the parametrization introduced in Ref. [22]:

which involves four free parameters—$ \hat q_0, \alpha, \beta $, and γ—that are constrained by experimental data. This parametrization is primarily motivated by a combination of physical considerations and numerical practicality [22], and predictions obtained using this form have been shown to agree well with data within the present theoretical framework. Similar parameterizations of $ \hat q $ have recently been adopted by other groups and have shown good performance in their studies [24, 27].

The global analysis incorporates 227 data points, including world data on the TMB of single-hadron production in semi-inclusive deep inelastic scattering (SIDIS) from eA collisions, Drell–Yan dilepton production in pA collisions, and heavy quarkonium ($ J/\psi $ and Υ) production in pA collisions. In addition, data on the nuclear shadowing effect in the deep inelastic scattering structure function are included [22].

In the analysis [22], the CT14 leading-order parton distribution functions with three active quark flavors [59] and the DSS fragmentation functions [60] are employed in the calculations. Heavy quarkonium production is described using the color evaporation model [51]. Since the TMB is expressed as a ratio in Eq. (2), it is relatively insensitive to non-perturbative inputs such as PDFs, FFs, and the quarkonium production model. Further details of the calculation can be found in Ref. [22]. It should also be noted that possible hadronization effects of the in-medium jet, which may reduce the impact of multiple scattering [61], are not taken into account in the present framework.

Our global extraction of $ \hat q $ is performed by analyzing world data from eA and pA collisions, primarily for the TMB of particle production, within the higher-twist expansion framework of perturbative QCD [22]. In this formalism, the effects of parton multiple scattering in the nuclear medium (target) are incorporated through twist-4 corrections [36, 51]. Taking TMB as an example, it is defined as the difference in the average transverse momentum squared, $ \langle p_T^2\rangle $, of produced particles between eA (pA) and $ ep $ ($ pp $) collisions:

Our calculation of TMB includes the leading-twist single-scattering cross section and the twist-4 double-scattering contribution enhanced by the large nuclear radius. Using pA collisions as an example, at twist-4 and leading order (LO) in $ \alpha_s $, the TMB can be schematically expressed as [36, 51]

where the superscripts ''S'' and ''D'' denote single- and double-scattering processes, respectively. Here, $ \mathrm{d}\sigma $ and $ \mathrm{d}\hat{\sigma} $ represent the differential cross sections at the observable level and for the perturbatively calculable partonic process, respectively. The terms f and $ D_h $ denote the relevant parton distribution functions (PDFs) and fragmentation functions (FFs), with the subscripts $ {q(g)/p} $ and $ {q(g)/A} $ indicating an incoming quark (or gluon) from the proton and nucleus, respectively. In the $ p_T^2 $-weighted double-scattering cross section ${\mathrm{ d}}\langle p_T^2\sigma^D\rangle $, $ T_{ij} $ is a twist-4 parton–parton correlation function that characterizes the relevant non-perturbative structure of the nuclear target in these processes.

At leading order in perturbation theory, the relevant jet transport coefficient in nuclei can be related to the twist-4 parton–parton correlation functions $ T_{ij} $. For example, the quark–gluon correlation function $ T_{qg} $ is usually approximated as [22, 36]

where x is the momentum fraction carried by the parton from the nucleus, μ is the factorization scale, $ R_A $ is the nuclear radius for a nucleus with mass number A, $ f_{q/A} $ is the parton distribution function of a quark q in the nucleus, and $ \hat q(x,\mu^2) $ is the effective nuclear-geometry-averaged jet transport coefficient that we aim to extract. Similarly, for the gluon–gluon correlation function $ T_{gg} $, which is relevant to processes initiated by a gluon in the nucleus, we adopt the same form as in Eq. (3), with $ f_{q/A} $ replaced by the gluon distribution $ f_{g/A} $. Throughout this paper, $ \hat q $ denotes the effective transport coefficient for a quark jet.

With this formalism, a global analysis of $ \hat q(x,Q^2) $ can be performed. It is important to note that the nonperturbative quantities directly constrained in such an analysis are the twist-4 correlation functions $ T_{ij} $. The quantity $ \hat q $ extracted in this manner corresponds, at leading order, to the transverse momentum broadening of a quark jet per unit propagation length, as shown in Eq. (2). Moreover, the relation between $ \hat q $ and $ T_{ij} $ may not hold exactly when higher-order contributions are included [37−43], although recent next-to-leading-order studies suggest that these corrections may be small [57].

On the other hand, a close connection exists between the twist-4 approach and formalisms that resum an arbitrary number of independent rescatterings, such as the Baier–Dokshitzer–Mueller–Peigne–Schiff (BDMPS) fra-mework [12] (see Ref. [12] for a general discussion). The correspondence between the HT formalism and color glass condensate effective theory has also been demonstrated recently [56, 58]. Nevertheless, it is important to emphasize that the interpretation of the variables x and $ Q^2 $ entering $ \hat q $ may differ across theoretical treatments. For example, in the BDMPS approach, x typically denotes the momentum fraction of the rescattered gluon in the nucleus, and $ Q^2 $ represents the scale of the rescattering process [12, 24]. These variables are therefore defined in a way that is more directly related to the intrinsic transport properties of the medium. In contrast, the $ \hat q(x,Q^2) $ factorized from the $ T_{ij} $ via Eq. (3) directly inherits the variables of the twist-4 correlation functions. In this context, x can be understood as the total momentum fraction carried by the two nuclear partons participating in double scattering. This fraction is fixed by the scattering kinematics (e.g., the Bjorken variable) and is related to the momentum fraction of the medium gluon by integration over the variables involved in the twist-4 matrix element $ T(x,x_1,x_2) $ [36, 51]. The scale $ Q^2 $ in $ T_{ij} $ serves as the factorization scale in QCD factorization and is typically taken to be the hard-scattering scale in practical calculations. Although a recent study [24] based on the BDMPS formalism has shown behavior of $ \hat q(x,Q^2) $ similar to our results, it is important to keep this difference in definition in mind when comparing the variable dependence between the two approaches. In addition, because x and $ Q^2 $ in our formalism are kinematically linked to the jet energy, the extracted $ \hat q(x,Q^2) $ naturally encodes the dependence of $ \hat q $ on the jet energy, as will be further discussed in Sec. II.C.

In the global analysis, to properly account for the functional form and possible kinematic dependence of $ \hat q $, we adopt the parametrization introduced in Ref. [22]:

which involves four free parameters—$ \hat q_0, \alpha, \beta $, and γ—that are constrained by experimental data. This parametrization is primarily motivated by a combination of physical considerations and numerical practicality [22], and predictions obtained using this form have been shown to agree well with data within the present theoretical framework. Similar parameterizations of $ \hat q $ have recently been adopted by other groups and have shown good performance in their studies [24, 27].

The global analysis incorporates 227 data points, including world data on the TMB of single-hadron production in semi-inclusive deep inelastic scattering (SIDIS) from eA collisions, Drell–Yan dilepton production in pA collisions, and heavy quarkonium ($ J/\psi $ and Υ) production in pA collisions. In addition, data on the nuclear shadowing effect in the deep inelastic scattering structure function are included [22].

In the analysis [22], the CT14 leading-order parton distribution functions with three active quark flavors [59] and the DSS fragmentation functions [60] are employed in the calculations. Heavy quarkonium production is described using the color evaporation model [51]. Since the TMB is expressed as a ratio in Eq. (2), it is relatively insensitive to non-perturbative inputs such as PDFs, FFs, and the quarkonium production model. Further details of the calculation can be found in Ref. [22]. It should also be noted that possible hadronization effects of the in-medium jet, which may reduce the impact of multiple scattering [61], are not taken into account in the present framework.

By performing a global analysis using the MINUIT package [62], an optimal form of $ \hat q(x,Q^2) $ is obtained, with a minimum $ \chi^2/\mathit{NDP}=1.21 $, where NDP denotes the number of data points (NDP = 227). In previous work [22], the Lagrange multiplier method was employed to estimate the uncertainties of selected calculated observables. However, this approach does not readily provide uncertainties for $ \hat q(x,Q^2) $ as a function of x and $ Q^2 $, thereby limiting the predictive power of the analysis and motivating the present improvement using the Hessian matrix method [63–64].

The Hessian analysis assumes that the $ \chi^2 $ function can be approximated by a quadratic form in terms of the free parameters $ \{a_i\} $ near the minimum:

where $ \chi_0^2\equiv\chi^2(\{a_i^0\}) $ is the global minimum at the optimal parameter values $ \{a_i\}=\{a_i^0\} $, $ y_i=a_i-a_i^0 $ denotes the deviation of $ a_i $ from its optimal value $ a_i^0 $, and $ H_{ij} $ are the elements of the Hessian matrix, defined as

In general, correlations among the parameters lead to nonzero off-diagonal elements in the Hessian matrix, complicating the direct estimation of uncertainties for a given tolerance $ \Delta\chi^2\equiv\chi^2-\chi^2_0 $. This issue can be addressed by diagonalizing the Hessian matrix through a linear transformation to a new basis $ \{z_i\} $, in which the $ \chi^2 $ increment takes the simple form

In this diagonal representation, the $ \chi^2 $ contours become spherical, enabling the straightforward generation of uncertainty sets for $ \hat q $ and related theoretical predictions. Further technical details of the Hessian procedure can be found in Refs. [63–64]. Briefly, by applying the inverse transformation from the new parameters $ \{z_i\} $ back to the original parameters $ \{y_i\} $, one obtains the uncertainty sets $ \{a_i^{(\pm k)}\} $ (with $ k=1,\ldots,n $, where n is the number of parameters) for the original parameters $ \{a_i\} $, corresponding to a given tolerance $ \Delta\chi^2 $. Furthermore, for any observable $ {\cal{Q}} $ that depends on $ \{a_i\} $, its uncertainty can be estimated using

where $ {\cal{Q}}_{\pm k}\equiv{\cal{Q}}(\{a_i^{(\pm k)}\}) $ denotes the uncertainty sets corresponding to the parameter sets $ a_i^{(\pm k)} $. For example, the uncertainty sets for $ \hat q $ are defined as

and the uncertainty sets for any quantity that depends on $ \hat q $ can be expressed as $ {\cal{Q}}_{\pm k}={\cal{Q}}(S_{\pm k}) $.

In the present work, the Hessian analysis is performed using the MINUIT package in combination with the ITERATE program [63]. Figure 1 shows the variation of $ \chi^2-\chi_0^2 $ with respect to each original parameter $ \{a_i\} $ near the minimum. The data are sensitive to all parameters, and the nonzero optimal values of α, β, and γ suggest a non-negligible dependence of $ \hat q $ on x and $ Q^2 $. This kinematic dependence is further examined in light of the uncertainties determined in this analysis. Figure 2 shows $ \chi^2-\chi_0^2 $ as a function of the transformed parameters $ \{z_i\} $. The excellent agreement with the quadratic form $ \Delta\chi^2=z_i^2 $ confirms the validity and reliability of the Hessian approach.

Figure 1.  (color online) $ \Delta\chi^2=\chi^2-\chi^2_0 $ as a function of each original parameter $ a_i $ in the neighborhood of the minimum at $ \chi^2=\chi^2_0 $. Here, $ \{a_i\}=\{\hat q_0, \alpha, \beta, \gamma\} $. The dashed line represents the tolerance $ \Delta\chi^2=35 $.

Figure 2.  (color online) Values of $ \Delta\chi^2=\chi^2-\chi^2_0 $ are shown as functions of the new parameters $ \{z_i\} $ obtained from the Hessian analysis. The gray curve shows the reference quadratic function $ \Delta\chi^2=z^2 $.

Using this Hessian analysis, we obtain the optimal values of the parameters $ \{a_i\} $, together with their uncertainties at the 90% confidence level (C.L.):

In Fig. 3, we show the extracted $ \hat q(x, Q^2) $ as a function of the momentum fraction x for fixed values of $ Q^2 = 1.2, 10, $ and $ 100\; {\rm{GeV}}^2 $. The central curves correspond to the optimal fit ($ S_0 $), and the shaded bands represent the uncertainty envelopes obtained from the eight Hessian eigenvector sets $ S_{\pm k} $. The optimal results exhibit a noticeable x dependence in both the small- and large-x regions, whereas the $ Q^2 $ dependence is relatively mild. Because most of the current data lie in the intermediate-x and intermediate-$ Q^2 $ regions, the uncertainties of $ \hat q $ increase at small and large values of x or $ Q^2 $. This increase is especially pronounced for $ x \gtrsim 0.4 $, where no data currently exist. For comparison, Fig. 3 also includes the result of a fit in which $ \hat q $ is assumed to be a kinematics-independent constant, yielding $ \hat q = \hat q_0 = 0.0150^{+0.0023}_{-0.0025} {\rm{GeV}}^2/{\rm{fm}} $. This constant fit gives a significantly larger minimum $ \chi^2/{\rm{NDP}} = 2.55 $, indicating a clear preference for the kinematics-dependent parametrization.

Figure 3.  (color online) $ \hat q(x,Q^2) $ extracted from the global analysis as a function of the momentum fraction x of the initial-state nuclear parton at $ Q^2= $1.2, 10, and 100 GeV$ ^2 $. The green solid curve represents the optimal values ($ S_0 $), the light-green band represents the corresponding uncertainty, and the dotted curves show the uncertainty set $ \{S_{\pm k}\} $ of $ \hat q(x,Q^2) $. For reference, the kinematics-independent constant $ \hat q $ and the uncertainty extracted from the data are shown as the red dashed line and pink band, respectively.

Since most of the data points are concentrated in the intermediate-x and intermediate-$ Q^2 $ regions, the uncertainty of the constant $ \hat q $, determined with the same tolerance $ \Delta\chi^2=35 $, is primarily driven by the data in those regions. It is therefore significantly smaller than the uncertainty of $ \hat q(x,Q^2) $ over most of the kinematic range.

As an additional test of the Hessian uncertainty evaluation, we compare the results obtained here with those from a fit that excludes the recent LHC data on the TMB of $ J/\psi $ production (denoted $ J/\psi $(21') [65]), as shown in Fig. 4. The best-fit values from the two fits are in good agreement. Notably, including the additional constraints from the $ J/\psi $(21') data significantly reduces the uncertainties at small x and reveals a more pronounced x dependence. This improved constraint arises mainly from data in the forward-rapidity region. This dataset was not included in our previous analysis, in which the small-x behavior of $ \hat q $ was primarily driven by the $ J/\psi $ TMB data published in 2015 [8].

Figure 4.  (color online) $ \hat q(x,Q^2) $ at $ Q^2=10 $ GeV$ ^2 $ and its uncertainties, extracted from a global analysis with and without the ALICE data (2021) on TMB of $ J/\psi $ production. The results of the new analysis are identical to those in Fig. 3. Results without the $ J/\psi $ (21') data are shown as a dashed curve with dotted boundaries.

By performing a global analysis using the MINUIT package [62], an optimal form of $ \hat q(x,Q^2) $ is obtained, with a minimum $ \chi^2/\mathit{NDP}=1.21 $, where NDP denotes the number of data points (NDP = 227). In previous work [22], the Lagrange multiplier method was employed to estimate the uncertainties of selected calculated observables. However, this approach does not readily provide uncertainties for $ \hat q(x,Q^2) $ as a function of x and $ Q^2 $, thereby limiting the predictive power of the analysis and motivating the present improvement using the Hessian matrix method [63–64].

The Hessian analysis assumes that the $ \chi^2 $ function can be approximated by a quadratic form in terms of the free parameters $ \{a_i\} $ near the minimum:

where $ \chi_0^2\equiv\chi^2(\{a_i^0\}) $ is the global minimum at the optimal parameter values $ \{a_i\}=\{a_i^0\} $, $ y_i=a_i-a_i^0 $ denotes the deviation of $ a_i $ from its optimal value $ a_i^0 $, and $ H_{ij} $ are the elements of the Hessian matrix, defined as

In general, correlations among the parameters lead to nonzero off-diagonal elements in the Hessian matrix, complicating the direct estimation of uncertainties for a given tolerance $ \Delta\chi^2\equiv\chi^2-\chi^2_0 $. This issue can be addressed by diagonalizing the Hessian matrix through a linear transformation to a new basis $ \{z_i\} $, in which the $ \chi^2 $ increment takes the simple form

In this diagonal representation, the $ \chi^2 $ contours become spherical, enabling the straightforward generation of uncertainty sets for $ \hat q $ and related theoretical predictions. Further technical details of the Hessian procedure can be found in Refs. [63–64]. Briefly, by applying the inverse transformation from the new parameters $ \{z_i\} $ back to the original parameters $ \{y_i\} $, one obtains the uncertainty sets $ \{a_i^{(\pm k)}\} $ (with $ k=1,\ldots,n $, where n is the number of parameters) for the original parameters $ \{a_i\} $, corresponding to a given tolerance $ \Delta\chi^2 $. Furthermore, for any observable $ {\cal{Q}} $ that depends on $ \{a_i\} $, its uncertainty can be estimated using

where $ {\cal{Q}}_{\pm k}\equiv{\cal{Q}}(\{a_i^{(\pm k)}\}) $ denotes the uncertainty sets corresponding to the parameter sets $ a_i^{(\pm k)} $. For example, the uncertainty sets for $ \hat q $ are defined as

and the uncertainty sets for any quantity that depends on $ \hat q $ can be expressed as $ {\cal{Q}}_{\pm k}={\cal{Q}}(S_{\pm k}) $.

In the present work, the Hessian analysis is performed using the MINUIT package in combination with the ITERATE program [63]. Figure 1 shows the variation of $ \chi^2-\chi_0^2 $ with respect to each original parameter $ \{a_i\} $ near the minimum. The data are sensitive to all parameters, and the nonzero optimal values of α, β, and γ suggest a non-negligible dependence of $ \hat q $ on x and $ Q^2 $. This kinematic dependence is further examined in light of the uncertainties determined in this analysis. Figure 2 shows $ \chi^2-\chi_0^2 $ as a function of the transformed parameters $ \{z_i\} $. The excellent agreement with the quadratic form $ \Delta\chi^2=z_i^2 $ confirms the validity and reliability of the Hessian approach.

Figure 1.  (color online) $ \Delta\chi^2=\chi^2-\chi^2_0 $ as a function of each original parameter $ a_i $ in the neighborhood of the minimum at $ \chi^2=\chi^2_0 $. Here, $ \{a_i\}=\{\hat q_0, \alpha, \beta, \gamma\} $. The dashed line represents the tolerance $ \Delta\chi^2=35 $.

Figure 2.  (color online) Values of $ \Delta\chi^2=\chi^2-\chi^2_0 $ are shown as functions of the new parameters $ \{z_i\} $ obtained from the Hessian analysis. The gray curve shows the reference quadratic function $ \Delta\chi^2=z^2 $.

Using this Hessian analysis, we obtain the optimal values of the parameters $ \{a_i\} $, together with their uncertainties at the 90% confidence level (C.L.):

In Fig. 3, we show the extracted $ \hat q(x, Q^2) $ as a function of the momentum fraction x for fixed values of $ Q^2 = 1.2, 10, $ and $ 100\; {\rm{GeV}}^2 $. The central curves correspond to the optimal fit ($ S_0 $), and the shaded bands represent the uncertainty envelopes obtained from the eight Hessian eigenvector sets $ S_{\pm k} $. The optimal results exhibit a noticeable x dependence in both the small- and large-x regions, whereas the $ Q^2 $ dependence is relatively mild. Because most of the current data lie in the intermediate-x and intermediate-$ Q^2 $ regions, the uncertainties of $ \hat q $ increase at small and large values of x or $ Q^2 $. This increase is especially pronounced for $ x \gtrsim 0.4 $, where no data currently exist. For comparison, Fig. 3 also includes the result of a fit in which $ \hat q $ is assumed to be a kinematics-independent constant, yielding $ \hat q = \hat q_0 = 0.0150^{+0.0023}_{-0.0025} {\rm{GeV}}^2/{\rm{fm}} $. This constant fit gives a significantly larger minimum $ \chi^2/{\rm{NDP}} = 2.55 $, indicating a clear preference for the kinematics-dependent parametrization.

Figure 3.  (color online) $ \hat q(x,Q^2) $ extracted from the global analysis as a function of the momentum fraction x of the initial-state nuclear parton at $ Q^2= $1.2, 10, and 100 GeV$ ^2 $. The green solid curve represents the optimal values ($ S_0 $), the light-green band represents the corresponding uncertainty, and the dotted curves show the uncertainty set $ \{S_{\pm k}\} $ of $ \hat q(x,Q^2) $. For reference, the kinematics-independent constant $ \hat q $ and the uncertainty extracted from the data are shown as the red dashed line and pink band, respectively.

Since most of the data points are concentrated in the intermediate-x and intermediate-$ Q^2 $ regions, the uncertainty of the constant $ \hat q $, determined with the same tolerance $ \Delta\chi^2=35 $, is primarily driven by the data in those regions. It is therefore significantly smaller than the uncertainty of $ \hat q(x,Q^2) $ over most of the kinematic range.

As an additional test of the Hessian uncertainty evaluation, we compare the results obtained here with those from a fit that excludes the recent LHC data on the TMB of $ J/\psi $ production (denoted $ J/\psi $(21') [65]), as shown in Fig. 4. The best-fit values from the two fits are in good agreement. Notably, including the additional constraints from the $ J/\psi $(21') data significantly reduces the uncertainties at small x and reveals a more pronounced x dependence. This improved constraint arises mainly from data in the forward-rapidity region. This dataset was not included in our previous analysis, in which the small-x behavior of $ \hat q $ was primarily driven by the $ J/\psi $ TMB data published in 2015 [8].

Figure 4.  (color online) $ \hat q(x,Q^2) $ at $ Q^2=10 $ GeV$ ^2 $ and its uncertainties, extracted from a global analysis with and without the ALICE data (2021) on TMB of $ J/\psi $ production. The results of the new analysis are identical to those in Fig. 3. Results without the $ J/\psi $ (21') data are shown as a dashed curve with dotted boundaries.

In jet tomography of cold nuclear matter, the kinematic dependence of transport properties is of particular interest because it provides insight into the detailed partonic structure of nuclear matter.

To illustrate how the experimental data constrain the nontrivial x and $ Q^2 $ dependence of $ \hat q $ in our analysis, we compare in Fig. 5 the results obtained using the extracted $ \hat q(x, Q^2) $ with those obtained using a constant $ \hat q = \hat q_0 $ for four representative observables included in the fit: the TMB in the Drell–Yan process, in SIDIS, and in $ J/\psi $ production at backward and forward rapidity at the LHC.

Figure 5.  (color online) Results for TMB in the Drell–Yan process in pA collisions as a function of the nuclear mass number A [panel (a)], in SIDIS as a function of Bjorken $ x_B $ [panel (b)], and in $ J/\psi $ production in the backward [panel (c)] and forward [panel (d)] regions in pA collisions at the LHC (2015) as functions of $ N_{ \rm{coll}} $. The green solid curve represents the results obtained with the optimal $ \hat q(x, Q^2) $, and the light-green shaded area indicates the corresponding uncertainty. The red dashed line and dotted band represent the results obtained with constant $ \hat q $ and its uncertainty. The average momentum fraction $ \langle x\rangle $ and probing scale $ \langle Q^2\rangle $ are shown for reference. Experimental data are taken from Refs. [6, 8, 66].

For the Drell–Yan process [panel (a)] and backward $ J/\psi $ production [panel (c)], the nuclear parton momentum fraction x lies in the intermediate region, and the results obtained with $ \hat q(x, Q^2) $ and with the constant $ \hat q $ are comparable. In contrast, the kinematic-dependent $ \hat q(x, Q^2) $ provides significantly improved theoretical descriptions for SIDIS [panel (b)] and forward $ J/\psi $ production [panel (d)], which probe the large- and small-x regions, respectively. In particular, the calculation using the constant $ \hat q $ from the global analysis fails to describe the forward $ J/\psi $ production data at the LHC. Notably, the enhanced values of $ \hat q(x, Q^2) $ in the small- and large-x regions observed in Fig. 3 originate primarily from forward $ J/\psi $ production at the LHC and the Bjorken-x dependence of SIDIS, respectively.

These enhancements of $ \hat q(x, Q^2) $ are governed by the parameters α and β in the adopted parametrization. The extracted value of α in Eq. (10) is qualitatively consistent with the growth rate of the gluon density anticipated in saturation physics [67], whereas the negative β may indicate an enhancement of nuclear power corrections at large x [68−70]. It is worth noting that, within the present framework, several effects are not yet incorporated, such as higher-order and higher-twist contributions, the dynamics of coherent multiple scattering, and potential violations of collinear factorization. These effects could play a role in certain kinematic regions and may influence the extracted x and $ Q^2 $ dependence of $ \hat q $ or the correlation functions $ T_{ij} $.

For the observables included in our analysis, the energy of the propagating parton $ E_{ \rm{jet}} $ in the rest frame of the nucleus can be estimated from the Lorentz-invariant variables x and $ Q^2 $ via the relation:

where $ Q^2 $ denotes the virtuality (or squared invariant mass) of the virtual photon in SIDIS or the Drell–Yan process, or the squared invariant mass of the heavy-quark pair that forms the quarkonium; x is the corresponding Bjorken variable; and $ m_N $ is the nucleon mass. Using this relation, we convert the extracted $ \hat q(x, Q^2) $ into $ \hat q(E_{ \rm{jet}}, Q^2) $, as shown in Fig. 6. For $ Q^2 = 2 $–$ 10 \,\, {\rm{GeV}}^2 $, $ \hat q $ increases with $ E_{ \rm{jet}} $ over a broad energy range, corresponding to the small-x region. Similar trends have been observed in jet-quenching studies, e.g., in Refs. [14, 28, 32]. In contrast, for $ Q^2 = 100 \,\, {\rm{GeV}}^2 $, the plotted region corresponds to large x, and $ \hat q $ exhibits a decreasing trend with $ E_{ \rm{jet}} $, albeit with large uncertainties. Moreover, the $ Q^2 $ dependence of $ \hat q(E_{ \rm{jet}}, Q^2) $ is more pronounced than that of $ \hat q(x, Q^2) $, because at fixed $ E_{ \rm{jet}} $, the corresponding x necessarily varies with $ Q^2 $. Because the energy dependence of $ \hat q $ is often discussed in the context of jet quenching in the quark–gluon plasma produced in relativistic heavy-ion collisions [14, 28, 30, 32], the present determination of $ \hat q(E_{ \rm{jet}}, Q^2) $ for cold nuclear matter provides a valuable reference for future comparative studies.

Figure 6.  (color online) The extracted $ \hat q $ as a function of the jet energy $ E_{ \rm{jet}} $. The solid curve and shaded band represent the optimal values and associated uncertainties, respectively.

In jet tomography of cold nuclear matter, the kinematic dependence of transport properties is of particular interest because it provides insight into the detailed partonic structure of nuclear matter.

To illustrate how the experimental data constrain the nontrivial x and $ Q^2 $ dependence of $ \hat q $ in our analysis, we compare in Fig. 5 the results obtained using the extracted $ \hat q(x, Q^2) $ with those obtained using a constant $ \hat q = \hat q_0 $ for four representative observables included in the fit: the TMB in the Drell–Yan process, in SIDIS, and in $ J/\psi $ production at backward and forward rapidity at the LHC.

Figure 5.  (color online) Results for TMB in the Drell–Yan process in pA collisions as a function of the nuclear mass number A [panel (a)], in SIDIS as a function of Bjorken $ x_B $ [panel (b)], and in $ J/\psi $ production in the backward [panel (c)] and forward [panel (d)] regions in pA collisions at the LHC (2015) as functions of $ N_{ \rm{coll}} $. The green solid curve represents the results obtained with the optimal $ \hat q(x, Q^2) $, and the light-green shaded area indicates the corresponding uncertainty. The red dashed line and dotted band represent the results obtained with constant $ \hat q $ and its uncertainty. The average momentum fraction $ \langle x\rangle $ and probing scale $ \langle Q^2\rangle $ are shown for reference. Experimental data are taken from Refs. [6, 8, 66].

For the Drell–Yan process [panel (a)] and backward $ J/\psi $ production [panel (c)], the nuclear parton momentum fraction x lies in the intermediate region, and the results obtained with $ \hat q(x, Q^2) $ and with the constant $ \hat q $ are comparable. In contrast, the kinematic-dependent $ \hat q(x, Q^2) $ provides significantly improved theoretical descriptions for SIDIS [panel (b)] and forward $ J/\psi $ production [panel (d)], which probe the large- and small-x regions, respectively. In particular, the calculation using the constant $ \hat q $ from the global analysis fails to describe the forward $ J/\psi $ production data at the LHC. Notably, the enhanced values of $ \hat q(x, Q^2) $ in the small- and large-x regions observed in Fig. 3 originate primarily from forward $ J/\psi $ production at the LHC and the Bjorken-x dependence of SIDIS, respectively.

These enhancements of $ \hat q(x, Q^2) $ are governed by the parameters α and β in the adopted parametrization. The extracted value of α in Eq. (10) is qualitatively consistent with the growth rate of the gluon density anticipated in saturation physics [67], whereas the negative β may indicate an enhancement of nuclear power corrections at large x [68−70]. It is worth noting that, within the present framework, several effects are not yet incorporated, such as higher-order and higher-twist contributions, the dynamics of coherent multiple scattering, and potential violations of collinear factorization. These effects could play a role in certain kinematic regions and may influence the extracted x and $ Q^2 $ dependence of $ \hat q $ or the correlation functions $ T_{ij} $.

For the observables included in our analysis, the energy of the propagating parton $ E_{ \rm{jet}} $ in the rest frame of the nucleus can be estimated from the Lorentz-invariant variables x and $ Q^2 $ via the relation:

where $ Q^2 $ denotes the virtuality (or squared invariant mass) of the virtual photon in SIDIS or the Drell–Yan process, or the squared invariant mass of the heavy-quark pair that forms the quarkonium; x is the corresponding Bjorken variable; and $ m_N $ is the nucleon mass. Using this relation, we convert the extracted $ \hat q(x, Q^2) $ into $ \hat q(E_{ \rm{jet}}, Q^2) $, as shown in Fig. 6. For $ Q^2 = 2 $–$ 10 \,\, {\rm{GeV}}^2 $, $ \hat q $ increases with $ E_{ \rm{jet}} $ over a broad energy range, corresponding to the small-x region. Similar trends have been observed in jet-quenching studies, e.g., in Refs. [14, 28, 32]. In contrast, for $ Q^2 = 100 \,\, {\rm{GeV}}^2 $, the plotted region corresponds to large x, and $ \hat q $ exhibits a decreasing trend with $ E_{ \rm{jet}} $, albeit with large uncertainties. Moreover, the $ Q^2 $ dependence of $ \hat q(E_{ \rm{jet}}, Q^2) $ is more pronounced than that of $ \hat q(x, Q^2) $, because at fixed $ E_{ \rm{jet}} $, the corresponding x necessarily varies with $ Q^2 $. Because the energy dependence of $ \hat q $ is often discussed in the context of jet quenching in the quark–gluon plasma produced in relativistic heavy-ion collisions [14, 28, 30, 32], the present determination of $ \hat q(E_{ \rm{jet}}, Q^2) $ for cold nuclear matter provides a valuable reference for future comparative studies.

Figure 6.  (color online) The extracted $ \hat q $ as a function of the jet energy $ E_{ \rm{jet}} $. The solid curve and shaded band represent the optimal values and associated uncertainties, respectively.

The broad kinematic coverage and high-precision measurements expected from upcoming electron–ion collision experiments present a unique opportunity to deepen our understanding of the structure of cold nuclear matter. In particular, measurements of nuclear-induced TMB provide an excellent means of performing a kinematic scan of the jet transport coefficient $ \hat q $.

In this section, we study TMB in both single-particle and particle-pair production at EICs. For the latter, the TMB of a particle pair is equivalent to the nuclear enhancement of its transverse momentum imbalance. Since the initial-state projectile is an electron, the observed TMB arises from final-state multiple scattering between the outgoing parton and the nucleus.

Specifically, using the $ \hat q $ extracted in the previous section, we calculate the TMB for single-hadron production within the HT framework [51−53], as well as the enhancement of transverse momentum imbalance for di-hadron and heavy-meson pair ($ D\bar D $) production. These observables are investigated in the kinematic regions relevant to three proposed EIC facilities: the US-EIC, EicC, and the Jefferson Lab 12 GeV program [1, 49, 50].

Since we focus on electron–ion collisions, it is useful to introduce the standard DIS kinematic variables:

Here, $ x_B $ is the Bjorken variable, y denotes the inelasticity of the scattering, and $ Q^2 $ is the virtuality of the exchanged photon $ \gamma^* $. In these definitions, $ k_e $, $ p_N $, and $ q_\gamma $ are the four-momenta of the incoming electron, the target nucleon, and the virtual photon, respectively. For the final-state fragmentation process, the momentum fraction $ z_h $ carried by the produced hadron is defined as

where $ p_h $ is the four-momentum of the final-state hadron. Additionally, the squared invariant mass of the $ \gamma^* $–nucleon system is given by $ W^2 = (q_\gamma + p_N)^2 $.

In our calculations, the center-of-mass energies $ \sqrt{s} $ of the electron–nucleon system for the three EIC facilities are set to 90 GeV (US-EIC), 10.6 GeV (EicC), and 4.8 GeV (JLab), respectively [1, 49, 50, 71]. The kinematic region considered in our analysis is defined by the following cuts: $ Q^2 \gt 1 \,\, {\rm{GeV}}^2 $, $ 0.01 \lt y \lt 0.95 $, $ W^2 \gt 10\,\, {\rm{GeV}}^2 $ (with $ W^2 \gt 4\,\, {\rm{GeV}}^2 $ for JLab), and $ 0.3 \lt z_h \lt 0.8 $ [71, 72]. Under these kinematic constraints, Fig. 7 shows the ranges of Bjorken x and $ Q^2 $ accessible at each facility. For comparison, we also overlay representative values of x and $ Q^2 $ from the current data included in our global analysis. Future EIC experiments will enable a comprehensive scan of the kinematic dependence of $ \hat q $, particularly in the small- and large-x regions, where existing measurements are scarce.

Figure 7.  (color online) Kinematic ($ x_B $ and $ Q^2 $) regions covered by three EIC facilities are shown as shaded areas. The boundaries for US-EIC, EicC, and JLab are plotted as dashed, dash-dotted, and dotted curves, respectively. The discrete circles represent the sampled momentum fraction x and probing scale $ Q^2 $ used in the numerical calculations for the current analysis of $ \hat q $ (the density of the circles does not represent the density of data points).

In addition to broad kinematic coverage, the high precision anticipated at future EIC facilities will provide significantly tighter constraints on $ \hat q $. To estimate the measurement uncertainties in this study, we assume an integrated luminosity of $ {\cal{L}} = 5 \, \, {\rm{fb}}^{-1} $ [49, 50, 71] and evaluate the relative statistical uncertainty as $ \delta_{{\rm{st}}} = 1/\sqrt{\sigma {\cal{L}}} $, where σ denotes the cross section of the process under consideration. In the absence of detailed information on systematic uncertainties, we assume that they are comparable in magnitude to the statistical uncertainty and therefore include an overall factor of $ \sqrt{2} $ to estimate the total uncertainty.

The broad kinematic coverage and high-precision measurements expected from upcoming electron–ion collision experiments present a unique opportunity to deepen our understanding of the structure of cold nuclear matter. In particular, measurements of nuclear-induced TMB provide an excellent means of performing a kinematic scan of the jet transport coefficient $ \hat q $.

In this section, we study TMB in both single-particle and particle-pair production at EICs. For the latter, the TMB of a particle pair is equivalent to the nuclear enhancement of its transverse momentum imbalance. Since the initial-state projectile is an electron, the observed TMB arises from final-state multiple scattering between the outgoing parton and the nucleus.

Specifically, using the $ \hat q $ extracted in the previous section, we calculate the TMB for single-hadron production within the HT framework [51−53], as well as the enhancement of transverse momentum imbalance for di-hadron and heavy-meson pair ($ D\bar D $) production. These observables are investigated in the kinematic regions relevant to three proposed EIC facilities: the US-EIC, EicC, and the Jefferson Lab 12 GeV program [1, 49, 50].

Since we focus on electron–ion collisions, it is useful to introduce the standard DIS kinematic variables:

Here, $ x_B $ is the Bjorken variable, y denotes the inelasticity of the scattering, and $ Q^2 $ is the virtuality of the exchanged photon $ \gamma^* $. In these definitions, $ k_e $, $ p_N $, and $ q_\gamma $ are the four-momenta of the incoming electron, the target nucleon, and the virtual photon, respectively. For the final-state fragmentation process, the momentum fraction $ z_h $ carried by the produced hadron is defined as

where $ p_h $ is the four-momentum of the final-state hadron. Additionally, the squared invariant mass of the $ \gamma^* $–nucleon system is given by $ W^2 = (q_\gamma + p_N)^2 $.

In our calculations, the center-of-mass energies $ \sqrt{s} $ of the electron–nucleon system for the three EIC facilities are set to 90 GeV (US-EIC), 10.6 GeV (EicC), and 4.8 GeV (JLab), respectively [1, 49, 50, 71]. The kinematic region considered in our analysis is defined by the following cuts: $ Q^2 \gt 1 \,\, {\rm{GeV}}^2 $, $ 0.01 \lt y \lt 0.95 $, $ W^2 \gt 10\,\, {\rm{GeV}}^2 $ (with $ W^2 \gt 4\,\, {\rm{GeV}}^2 $ for JLab), and $ 0.3 \lt z_h \lt 0.8 $ [71, 72]. Under these kinematic constraints, Fig. 7 shows the ranges of Bjorken x and $ Q^2 $ accessible at each facility. For comparison, we also overlay representative values of x and $ Q^2 $ from the current data included in our global analysis. Future EIC experiments will enable a comprehensive scan of the kinematic dependence of $ \hat q $, particularly in the small- and large-x regions, where existing measurements are scarce.

Figure 7.  (color online) Kinematic ($ x_B $ and $ Q^2 $) regions covered by three EIC facilities are shown as shaded areas. The boundaries for US-EIC, EicC, and JLab are plotted as dashed, dash-dotted, and dotted curves, respectively. The discrete circles represent the sampled momentum fraction x and probing scale $ Q^2 $ used in the numerical calculations for the current analysis of $ \hat q $ (the density of the circles does not represent the density of data points).

In addition to broad kinematic coverage, the high precision anticipated at future EIC facilities will provide significantly tighter constraints on $ \hat q $. To estimate the measurement uncertainties in this study, we assume an integrated luminosity of $ {\cal{L}} = 5 \, \, {\rm{fb}}^{-1} $ [49, 50, 71] and evaluate the relative statistical uncertainty as $ \delta_{{\rm{st}}} = 1/\sqrt{\sigma {\cal{L}}} $, where σ denotes the cross section of the process under consideration. In the absence of detailed information on systematic uncertainties, we assume that they are comparable in magnitude to the statistical uncertainty and therefore include an overall factor of $ \sqrt{2} $ to estimate the total uncertainty.

The TMB of single-hadron production in semi-inclusive deep inelastic scattering (SIDIS) has played an important role in current global analyses of $ \hat q $ and will remain a key observable at future EIC facilities. At leading order in perturbative QCD, the detected hadron originates from the fragmentation of a struck quark, which may undergo multiple scattering as it propagates through the nuclear medium. The leading-twist cross section for the single-scattering process is given by

where the sum runs over all quark flavors, $ e_q $ denotes the quark electric charge, $ f_{q/A} $ is the nuclear parton distribution function, and $ D_{h/q} $ represents the fragmentation function into hadron h.

The TMB, as defined in Eq. (1), can be expressed at fixed $ x_B $, $ Q^2 $, and $ z_h $ as [36, 53]

where the color factor $ C_F $ reflects the quark-jet character of the propagating parton. In our calculation, the factorization scale is set to $ \mu^2 = Q^2 $, and $ \Delta \langle p_T^2\rangle $ is averaged over the kinematic bin considered. The same non-perturbative inputs (PDFs and FFs) as those in the global analysis are used throughout.

In Fig. 8, we present the TMB for pion production in SIDIS in electron–lead (e-Pb) collisions at three proposed EIC facilities, for fixed $ Q^2 = 2 \,\, {\rm{GeV}}^2 $ (left panel) and $ 10 \,\,{\rm{GeV}}^2 $ (right panel). The solid curves and associated shaded bands correspond to predictions obtained using the kinematic-dependent $ \hat q(x, Q^2) $ and its uncertainties. The results show that $ \Delta \langle p_T^2 \rangle $ as a function of $ x_B $ closely reflects both the x dependence and the uncertainty structure of $ \hat q(x, Q^2) $. In particular, both the broadening and its uncertainties increase substantially in the small- and large-x regions. For comparison, predictions based on a kinematic-independent constant $ \hat q $ (dashed curves), together with its uncertainty bands (dotted curves), are also shown. As expected, the main differences between the two sets of predictions arise in the small- and large-x regimes. The estimated experimental uncertainties, shown in the lower part of each panel, are substantially smaller than the theoretical uncertainties. This suggests that future EIC measurements will be able to clearly discriminate between the two scenarios and provide strong constraints on $ \hat q $.

Figure 8.  (color online) 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 GeV$ ^2 $ for three EIC facilities. The solid curve with the shaded area represents results obtained using $ \hat q(x,Q^2) $ with uncertainties. The dashed curve with the dotted band shows results obtained using a constant $ \hat q $ with uncertainties. Circles with vertical bars represent estimated experimental uncertainties. In the calculations, we take $ Q^2\in[1.5,2.5] $ GeV$ ^2 $ (left panel) and $ Q^2\in[9,11] $ GeV$ ^2 $ (right panel), and set $ x_B<0.6 $. At the top of each panel, we indicate the kinematic ($ x_B $) ranges covered by the three facilities, showing the overlaps among them. The theoretical results for two facilities in their overlap region are generally similar because they depend primarily on $ x_B $ at fixed $ Q^2 $.

Figure 9 shows $ \Delta \langle p_T^2 \rangle $ as a function of $ Q^2 $ in an intermediate-$ x_B $ region, where the difference between the kinematic-dependent and constant-$ \hat q $ predictions is small. Owing to the mild scale dependence of $ \hat q(x, Q^2) $, only modest differences between the two predictions are observed over the $ Q^2 $ range considered. The slight decrease in the broadening with increasing $ Q^2 $ in the constant-$ \hat q $ case arises from the decrease in the average value of $ z_h^2 $ [cf. Eq. (15)], driven by the scale evolution of the fragmentation functions. As in Fig. 8, the projected experimental uncertainties remain small across the plotted $ Q^2 $ region.

Figure 9.  (color online) Similar to Fig. 8, but for the transverse momentum broadening $ \Delta \langle p_T^2\rangle $ of single-pion production in SIDIS as a function of $ Q^2 $ in the range $ 0.07<x_B<0.09 $. The results are calculated for $ 1<Q^2<100 $ GeV$ ^2 $.

The TMB of single-hadron production in semi-inclusive deep inelastic scattering (SIDIS) has played an important role in current global analyses of $ \hat q $ and will remain a key observable at future EIC facilities. At leading order in perturbative QCD, the detected hadron originates from the fragmentation of a struck quark, which may undergo multiple scattering as it propagates through the nuclear medium. The leading-twist cross section for the single-scattering process is given by

where the sum runs over all quark flavors, $ e_q $ denotes the quark electric charge, $ f_{q/A} $ is the nuclear parton distribution function, and $ D_{h/q} $ represents the fragmentation function into hadron h.

The TMB, as defined in Eq. (1), can be expressed at fixed $ x_B $, $ Q^2 $, and $ z_h $ as [36, 53]

where the color factor $ C_F $ reflects the quark-jet character of the propagating parton. In our calculation, the factorization scale is set to $ \mu^2 = Q^2 $, and $ \Delta \langle p_T^2\rangle $ is averaged over the kinematic bin considered. The same non-perturbative inputs (PDFs and FFs) as those in the global analysis are used throughout.

In Fig. 8, we present the TMB for pion production in SIDIS in electron–lead (e-Pb) collisions at three proposed EIC facilities, for fixed $ Q^2 = 2 \,\, {\rm{GeV}}^2 $ (left panel) and $ 10 \,\,{\rm{GeV}}^2 $ (right panel). The solid curves and associated shaded bands correspond to predictions obtained using the kinematic-dependent $ \hat q(x, Q^2) $ and its uncertainties. The results show that $ \Delta \langle p_T^2 \rangle $ as a function of $ x_B $ closely reflects both the x dependence and the uncertainty structure of $ \hat q(x, Q^2) $. In particular, both the broadening and its uncertainties increase substantially in the small- and large-x regions. For comparison, predictions based on a kinematic-independent constant $ \hat q $ (dashed curves), together with its uncertainty bands (dotted curves), are also shown. As expected, the main differences between the two sets of predictions arise in the small- and large-x regimes. The estimated experimental uncertainties, shown in the lower part of each panel, are substantially smaller than the theoretical uncertainties. This suggests that future EIC measurements will be able to clearly discriminate between the two scenarios and provide strong constraints on $ \hat q $.

Figure 8.  (color online) 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 GeV$ ^2 $ for three EIC facilities. The solid curve with the shaded area represents results obtained using $ \hat q(x,Q^2) $ with uncertainties. The dashed curve with the dotted band shows results obtained using a constant $ \hat q $ with uncertainties. Circles with vertical bars represent estimated experimental uncertainties. In the calculations, we take $ Q^2\in[1.5,2.5] $ GeV$ ^2 $ (left panel) and $ Q^2\in[9,11] $ GeV$ ^2 $ (right panel), and set $ x_B<0.6 $. At the top of each panel, we indicate the kinematic ($ x_B $) ranges covered by the three facilities, showing the overlaps among them. The theoretical results for two facilities in their overlap region are generally similar because they depend primarily on $ x_B $ at fixed $ Q^2 $.

Figure 9 shows $ \Delta \langle p_T^2 \rangle $ as a function of $ Q^2 $ in an intermediate-$ x_B $ region, where the difference between the kinematic-dependent and constant-$ \hat q $ predictions is small. Owing to the mild scale dependence of $ \hat q(x, Q^2) $, only modest differences between the two predictions are observed over the $ Q^2 $ range considered. The slight decrease in the broadening with increasing $ Q^2 $ in the constant-$ \hat q $ case arises from the decrease in the average value of $ z_h^2 $ [cf. Eq. (15)], driven by the scale evolution of the fragmentation functions. As in Fig. 8, the projected experimental uncertainties remain small across the plotted $ Q^2 $ region.

Figure 9.  (color online) Similar to Fig. 8, but for the transverse momentum broadening $ \Delta \langle p_T^2\rangle $ of single-pion production in SIDIS as a function of $ Q^2 $ in the range $ 0.07<x_B<0.09 $. The results are calculated for $ 1<Q^2<100 $ GeV$ ^2 $.

We now examine the nuclear-medium-induced enhancement of the transverse momentum imbalance in back-to-back particle production at future EIC facilities. The transverse momentum imbalance of a particle pair, defined as $ \boldsymbol{p}_T = \boldsymbol{p}_{1T} + \boldsymbol{p}_{2T} $, corresponds to the total transverse momentum of the two particles [52]. Thus, the enhancement of this imbalance is equivalent to the broadening of the total transverse momentum.

In this subsection, we focus on the imbalance in di-hadron, specifically pion-pair, production. The back-to-back hadron pair originates from the fragmentation of a di-jet system, which, at leading order in $ \alpha_s $, can be produced through the partonic processes $ \gamma^* q \rightarrow q g $ and $ \gamma^* g \rightarrow q \bar{q} $.

Following previous work [52], the differential cross section for di-hadron production can be expressed in the center-of-mass frame of the virtual-photon–nucleon ($ \gamma^* $-N) system as

In this expression, where the $ \mu^2 $ dependence of the PDFs and FFs has been suppressed, $ y_{1(2)} $, $ p_{1T(2T)} $, and $ z_{1(2)} $ denote the rapidity, transverse momentum, and momentum fraction of the hadron $ h_{1(2)} $, respectively. The function $ H^U_{\gamma^{\ast}b\rightarrow cd} $ denotes the perturbatively calculable hard part of the partonic subprocess $ \gamma^{\ast}b\rightarrow cd $ [52]. The nuclear enhancement of the di-hadron transverse momentum imbalance, driven by multiple scattering of the two outgoing partons c and d, is sensitive only to the total color charge of the two-parton composite state, which is equivalent to that of the initial-state nuclear parton b. Accordingly, for given $ y_{1(2)} $ and $ p_{1T(2T)} $, this enhancement is given by [52]

where the hard function $ H^F_{\gamma^{\ast}b\rightarrow cd} $ is written as

This expression involves the color factors $ C_F $ and $ C_A $, which correspond to processes initiated by a nuclear quark and a nuclear gluon, respectively. This introduces additional complexity compared with the SIDIS case, where only quark-jet transport contributes at leading order in $ \alpha_s $. In our calculation, the renormalization and factorization scales are set to the average transverse momentum of the di-hadron system, i.e., $ \mu = (p_{1T} + p_{2T})/2 $. Kinematic cuts of $ 1 \lt p_{1T(2T)} \lt 4 $ GeV and $ 0.1 \lt y_{1(2)} \lt 2.0 $ are also applied.

It is worth noting that, in the back-to-back di-hadron region, the potential role of Sudakov resummation should be considered. For instance, if one considers the differential transverse momentum imbalance distribution without integrating over the imbalance, the back-to-back region is characterized by two distinct scales: the momentum imbalance ($ p_T $) and the individual hadron transverse momentum ($ p_{1(2)T} $), with $ p_T\ll p_{1(2)T} $. In this regime, Transverse Momentum Dependent (TMD) factorization is required, which necessarily entails the implementation of Sudakov resummation [73]. In the present analysis, however, we integrate over the imbalance scale; therefore, the only remaining scale in our observable is the transverse momentum of the produced hadrons. Consequently, the appropriate theoretical framework for our calculation is collinear factorization (specifically, higher-twist factorization), within which Sudakov resummation is not present.

In Fig. 10, we present the nuclear enhancement of the di-hadron (pion pair) transverse momentum imbalance as a function of Bjorken $ x_B $, for fixed $ Q^2 = 2 \,\, {\rm{GeV}}^2 $ (left panel) and $ 10 \,\, {\rm{GeV}}^2 $ (right panel). Predictions obtained with both the kinematics-dependent $ \hat q(x, Q^2) $ and the constant $ \hat q $ are shown, along with estimated experimental uncertainties.

Figure 10.  (color online) Similar to Fig. 8, but for the nuclear enhancement of the transverse momentum imbalance, $ \Delta \langle p_T^2\rangle $, in di-hadron ($ \pi\pi $) production as a function of $ x_B $. The result for each bin is calculated in the $ \gamma^{\ast} $-N system with fixed $ Q^2 $ and $ x_B $ (central value). The green squares indicate the average momentum fraction $ x=x_B(1+M_{JJ}^2/Q^2) $ for the corresponding $ x_B $, with values shown on the right-hand vertical axis. In the right panel ($ Q^2=10 $ GeV), di-hadron production at JLab is kinematically forbidden or suppressed and is therefore not shown.

It is also worth noting that, for this process, the momentum fraction x of the nuclear parton entering $ \hat q(x, Q^2) $ is related to $ x_B $ via

where $ M_{JJ} $ denotes the invariant mass of the outgoing di-jet system. Consequently, the relevant x values are larger than $ x_B $. To illustrate this shift, we indicate in Fig. 10 the average $ \langle x \rangle $ corresponding to each $ x_B $ point using green square symbols. Thus, compared with single-hadron production, di-hadron production probes $ \hat q(x, Q^2) $ at larger momentum fractions for a given $ x_B $ range. This shift also affects the propagation of uncertainties: because the uncertainties of $ \hat q(x,Q^2) $ are associated with the variable x, they are mapped onto the observable at smaller $ x_B $. Consequently, the uncertainties in the small-$ x_B $ region do not increase significantly, as they are predominantly determined by those of $ \hat q(x,Q^2) $ in the intermediate-x region.

From Fig. 10, we observe that the broadening predicted with a constant $ \hat q $ decreases as $ x_B $ (or x) increases. This trend arises from the combination of the color factors $ C_A $ and $ C_F $ associated with the contributing partonic channels, together with the decreasing contribution of gluon-initiated processes at larger x. Since our calculation is performed in the $ \gamma^* $-nucleon center-of-mass frame, the experimental uncertainties are simply estimated using the same integrated luminosity as that for the e-A system, which serves as an illustrative reference.

Overall, future measurements of di-hadron imbalance at the EIC are expected to provide valuable constraints on the behavior of $ \hat q(x, Q^2) $ in the large-x region, where current uncertainties remain substantial and theoretical predictions are therefore less reliable.

We now examine the nuclear-medium-induced enhancement of the transverse momentum imbalance in back-to-back particle production at future EIC facilities. The transverse momentum imbalance of a particle pair, defined as $ \boldsymbol{p}_T = \boldsymbol{p}_{1T} + \boldsymbol{p}_{2T} $, corresponds to the total transverse momentum of the two particles [52]. Thus, the enhancement of this imbalance is equivalent to the broadening of the total transverse momentum.

In this subsection, we focus on the imbalance in di-hadron, specifically pion-pair, production. The back-to-back hadron pair originates from the fragmentation of a di-jet system, which, at leading order in $ \alpha_s $, can be produced through the partonic processes $ \gamma^* q \rightarrow q g $ and $ \gamma^* g \rightarrow q \bar{q} $.

Following previous work [52], the differential cross section for di-hadron production can be expressed in the center-of-mass frame of the virtual-photon–nucleon ($ \gamma^* $-N) system as

In this expression, where the $ \mu^2 $ dependence of the PDFs and FFs has been suppressed, $ y_{1(2)} $, $ p_{1T(2T)} $, and $ z_{1(2)} $ denote the rapidity, transverse momentum, and momentum fraction of the hadron $ h_{1(2)} $, respectively. The function $ H^U_{\gamma^{\ast}b\rightarrow cd} $ denotes the perturbatively calculable hard part of the partonic subprocess $ \gamma^{\ast}b\rightarrow cd $ [52]. The nuclear enhancement of the di-hadron transverse momentum imbalance, driven by multiple scattering of the two outgoing partons c and d, is sensitive only to the total color charge of the two-parton composite state, which is equivalent to that of the initial-state nuclear parton b. Accordingly, for given $ y_{1(2)} $ and $ p_{1T(2T)} $, this enhancement is given by [52]

where the hard function $ H^F_{\gamma^{\ast}b\rightarrow cd} $ is written as

This expression involves the color factors $ C_F $ and $ C_A $, which correspond to processes initiated by a nuclear quark and a nuclear gluon, respectively. This introduces additional complexity compared with the SIDIS case, where only quark-jet transport contributes at leading order in $ \alpha_s $. In our calculation, the renormalization and factorization scales are set to the average transverse momentum of the di-hadron system, i.e., $ \mu = (p_{1T} + p_{2T})/2 $. Kinematic cuts of $ 1 \lt p_{1T(2T)} \lt 4 $ GeV and $ 0.1 \lt y_{1(2)} \lt 2.0 $ are also applied.

It is worth noting that, in the back-to-back di-hadron region, the potential role of Sudakov resummation should be considered. For instance, if one considers the differential transverse momentum imbalance distribution without integrating over the imbalance, the back-to-back region is characterized by two distinct scales: the momentum imbalance ($ p_T $) and the individual hadron transverse momentum ($ p_{1(2)T} $), with $ p_T\ll p_{1(2)T} $. In this regime, Transverse Momentum Dependent (TMD) factorization is required, which necessarily entails the implementation of Sudakov resummation [73]. In the present analysis, however, we integrate over the imbalance scale; therefore, the only remaining scale in our observable is the transverse momentum of the produced hadrons. Consequently, the appropriate theoretical framework for our calculation is collinear factorization (specifically, higher-twist factorization), within which Sudakov resummation is not present.

In Fig. 10, we present the nuclear enhancement of the di-hadron (pion pair) transverse momentum imbalance as a function of Bjorken $ x_B $, for fixed $ Q^2 = 2 \,\, {\rm{GeV}}^2 $ (left panel) and $ 10 \,\, {\rm{GeV}}^2 $ (right panel). Predictions obtained with both the kinematics-dependent $ \hat q(x, Q^2) $ and the constant $ \hat q $ are shown, along with estimated experimental uncertainties.

Figure 10.  (color online) Similar to Fig. 8, but for the nuclear enhancement of the transverse momentum imbalance, $ \Delta \langle p_T^2\rangle $, in di-hadron ($ \pi\pi $) production as a function of $ x_B $. The result for each bin is calculated in the $ \gamma^{\ast} $-N system with fixed $ Q^2 $ and $ x_B $ (central value). The green squares indicate the average momentum fraction $ x=x_B(1+M_{JJ}^2/Q^2) $ for the corresponding $ x_B $, with values shown on the right-hand vertical axis. In the right panel ($ Q^2=10 $ GeV), di-hadron production at JLab is kinematically forbidden or suppressed and is therefore not shown.

It is also worth noting that, for this process, the momentum fraction x of the nuclear parton entering $ \hat q(x, Q^2) $ is related to $ x_B $ via

where $ M_{JJ} $ denotes the invariant mass of the outgoing di-jet system. Consequently, the relevant x values are larger than $ x_B $. To illustrate this shift, we indicate in Fig. 10 the average $ \langle x \rangle $ corresponding to each $ x_B $ point using green square symbols. Thus, compared with single-hadron production, di-hadron production probes $ \hat q(x, Q^2) $ at larger momentum fractions for a given $ x_B $ range. This shift also affects the propagation of uncertainties: because the uncertainties of $ \hat q(x,Q^2) $ are associated with the variable x, they are mapped onto the observable at smaller $ x_B $. Consequently, the uncertainties in the small-$ x_B $ region do not increase significantly, as they are predominantly determined by those of $ \hat q(x,Q^2) $ in the intermediate-x region.

From Fig. 10, we observe that the broadening predicted with a constant $ \hat q $ decreases as $ x_B $ (or x) increases. This trend arises from the combination of the color factors $ C_A $ and $ C_F $ associated with the contributing partonic channels, together with the decreasing contribution of gluon-initiated processes at larger x. Since our calculation is performed in the $ \gamma^* $-nucleon center-of-mass frame, the experimental uncertainties are simply estimated using the same integrated luminosity as that for the e-A system, which serves as an illustrative reference.

Overall, future measurements of di-hadron imbalance at the EIC are expected to provide valuable constraints on the behavior of $ \hat q(x, Q^2) $ in the large-x region, where current uncertainties remain substantial and theoretical predictions are therefore less reliable.

As in di-hadron production, nuclear enhancement of the transverse momentum imbalance is also expected in heavy-meson pair production. In this work, we investigate this effect for $ D\bar{D} $ pairs arising from the fragmentation of a back-to-back heavy-quark pair $ Q\bar{Q} $, specifically $ c\bar{c} $. At leading order in $ \alpha_s $, the $ c\bar{c} $ dijet is produced through the photon–gluon fusion process $ \gamma^* g \rightarrow c\bar{c} $. The differential cross section for $ \gamma^* + A \rightarrow D(p_1) + \bar{D}(p_2) + X $ in the center-of-mass frame of the $ \gamma^* $-nucleon system can be expressed in a form analogous to that given in Eq. (3) for di-hadron production [52]:

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

This simple expression reflects the fact that the $ c\bar{c} $ dijet is initiated solely by a nuclear gluon, corresponding to the color factor $ C_A $. In our calculation, the fragmentation of c ($ \bar{c} $) quarks into D ($ \bar{D} $) mesons is described by the KKKS08 fragmentation functions [74]. The renormalization and factorization scales are set to the average transverse mass of the $ D\bar{D} $ pair, i.e., $ \mu = (m_{1T} + m_{2T})/2 $. We adopt the same kinematic cuts as in the di-hadron production calculation. The nuclear parton momentum fraction x is again given by Eq. (19). In heavy meson pair production, owing to the heavy quark mass, the average momentum fraction x for a given $ x_B $ is typically larger than that in di-hadron production.

In Fig. 11, we present the nuclear enhancement of the transverse momentum imbalance for $ D\bar{D} $ pair production as a function of Bjorken $ x_B $, at fixed $ Q^2 = 2\,\,{\rm{GeV}}^2 $ (left panel) and $ 10\,\, {\rm{GeV}}^2 $ (right panel). The average value $ \langle x \rangle $ for each kinematic bin is also indicated in the figure. For a given momentum fraction x, the broadening $ \Delta \langle p_T^2 \rangle $ in $ D\bar{D} $ production is larger than that in di-hadron production, as expected from the larger color factor associated with gluon-initiated processes. This prediction can be tested in future EIC measurements.

Figure 11.  (color online) Similar to Fig. 10, but for the nuclear enhancement of the transverse momentum imbalance $ \Delta \langle p_T^2\rangle $ in heavy-meson pair ($D\bar{D} $) production as a function of $x_B $.

The three types of transverse momentum broadening at a future EIC studied in this section are expected to provide a multidimensional probe of the jet transport properties of cold nuclear matter over a broad kinematic range. Although the results presented here are obtained for electron–lead collisions, predictions for other nuclear targets can be readily obtained by rescaling with the nuclear radius, i.e., multiplying by $ R_A / R_{{\rm{Pb}}} $, in accordance with Eq. (3).

As in di-hadron production, nuclear enhancement of the transverse momentum imbalance is also expected in heavy-meson pair production. In this work, we investigate this effect for $ D\bar{D} $ pairs arising from the fragmentation of a back-to-back heavy-quark pair $ Q\bar{Q} $, specifically $ c\bar{c} $. At leading order in $ \alpha_s $, the $ c\bar{c} $ dijet is produced through the photon–gluon fusion process $ \gamma^* g \rightarrow c\bar{c} $. The differential cross section for $ \gamma^* + A \rightarrow D(p_1) + \bar{D}(p_2) + X $ in the center-of-mass frame of the $ \gamma^* $-nucleon system can be expressed in a form analogous to that given in Eq. (3) for di-hadron production [52]:

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

This simple expression reflects the fact that the $ c\bar{c} $ dijet is initiated solely by a nuclear gluon, corresponding to the color factor $ C_A $. In our calculation, the fragmentation of c ($ \bar{c} $) quarks into D ($ \bar{D} $) mesons is described by the KKKS08 fragmentation functions [74]. The renormalization and factorization scales are set to the average transverse mass of the $ D\bar{D} $ pair, i.e., $ \mu = (m_{1T} + m_{2T})/2 $. We adopt the same kinematic cuts as in the di-hadron production calculation. The nuclear parton momentum fraction x is again given by Eq. (19). In heavy meson pair production, owing to the heavy quark mass, the average momentum fraction x for a given $ x_B $ is typically larger than that in di-hadron production.

In Fig. 11, we present the nuclear enhancement of the transverse momentum imbalance for $ D\bar{D} $ pair production as a function of Bjorken $ x_B $, at fixed $ Q^2 = 2\,\,{\rm{GeV}}^2 $ (left panel) and $ 10\,\, {\rm{GeV}}^2 $ (right panel). The average value $ \langle x \rangle $ for each kinematic bin is also indicated in the figure. For a given momentum fraction x, the broadening $ \Delta \langle p_T^2 \rangle $ in $ D\bar{D} $ production is larger than that in di-hadron production, as expected from the larger color factor associated with gluon-initiated processes. This prediction can be tested in future EIC measurements.

Figure 11.  (color online) Similar to Fig. 10, but for the nuclear enhancement of the transverse momentum imbalance $ \Delta \langle p_T^2\rangle $ in heavy-meson pair ($D\bar{D} $) production as a function of $x_B $.

The three types of transverse momentum broadening at a future EIC studied in this section are expected to provide a multidimensional probe of the jet transport properties of cold nuclear matter over a broad kinematic range. Although the results presented here are obtained for electron–lead collisions, predictions for other nuclear targets can be readily obtained by rescaling with the nuclear radius, i.e., multiplying by $ R_A / R_{{\rm{Pb}}} $, in accordance with Eq. (3).

To gain insight into how future electron–ion collision experiments can advance our understanding of the jet transport properties of cold nuclear matter, we investigate nuclear-medium-induced transverse momentum broadening (TMB) in single- and pair-particle production within the kinematic regions accessible at three proposed future facilities. Our theoretical calculations account for multiple scattering of partons propagating through the nucleus and are performed within the framework of the higher-twist expansion.

A key ingredient in our study is a jet transport coefficient $ \hat q $ that depends on both the parton momentum fraction x and the probing scale $ Q^2 $. This quantity is extracted from a global analysis of existing experimental data, and its uncertainties are quantified using the Hessian matrix method, enabling robust predictions for future measurements. The extracted $ \hat q(x, Q^2) $ exhibits nontrivial kinematic dependence within the obtained uncertainty bands. Moreover, it is converted into a function of jet energy, which may provide valuable guidance for studies of jet quenching in heavy-ion collisions.

Using this $ \hat q(x, Q^2) $, we investigate three TMB-related observables at the EIC: TMB in single-hadron production and the nuclear enhancement of transverse momentum imbalance in di-hadron and heavy-meson pair ($ D\bar{D} $) production. These nuclear-induced broadening effects are found to be sensitive to the color state of the propagating parton and to exhibit a clear kinematic dependence inherited from $ \hat q(x, Q^2) $. For comparison, results obtained with a constant $ \hat q $ are also presented.

Our findings indicate that future TMB measurements at the EIC have significant potential to provide precise constraints on jet transport properties over a wide kinematic range and to advance tomographic imaging of the non-perturbative structure of cold nuclear matter.

To gain insight into how future electron–ion collision experiments can advance our understanding of the jet transport properties of cold nuclear matter, we investigate nuclear-medium-induced transverse momentum broadening (TMB) in single- and pair-particle production within the kinematic regions accessible at three proposed future facilities. Our theoretical calculations account for multiple scattering of partons propagating through the nucleus and are performed within the framework of the higher-twist expansion.

A key ingredient in our study is a jet transport coefficient $ \hat q $ that depends on both the parton momentum fraction x and the probing scale $ Q^2 $. This quantity is extracted from a global analysis of existing experimental data, and its uncertainties are quantified using the Hessian matrix method, enabling robust predictions for future measurements. The extracted $ \hat q(x, Q^2) $ exhibits nontrivial kinematic dependence within the obtained uncertainty bands. Moreover, it is converted into a function of jet energy, which may provide valuable guidance for studies of jet quenching in heavy-ion collisions.

Using this $ \hat q(x, Q^2) $, we investigate three TMB-related observables at the EIC: TMB in single-hadron production and the nuclear enhancement of transverse momentum imbalance in di-hadron and heavy-meson pair ($ D\bar{D} $) production. These nuclear-induced broadening effects are found to be sensitive to the color state of the propagating parton and to exhibit a clear kinematic dependence inherited from $ \hat q(x, Q^2) $. For comparison, results obtained with a constant $ \hat q $ are also presented.

Our findings indicate that future TMB measurements at the EIC have significant potential to provide precise constraints on jet transport properties over a wide kinematic range and to advance tomographic imaging of the non-perturbative structure of cold nuclear matter.