format: text

text:
# Multiple Parton Scattering in Nuclei: Parton Energy Loss

## I Introduction

Hard processes in high-energy strong interactions are always localized in space-time because of the large momentum-energy transfer. The asymptotic behavior of QCD allows one to compute these cross sections perturbatively. Together with the factorization theorem and the experimental information of parton distributions and fragmentation functions, hard processes in hadronic collisions have been well understood . One can then in turn use them as probes of nuclear matter as well as hot quark-gluon plasma which is expected to be formed in high-energy nuclear collisions. In particular, the propagation of an energetic parton and its induced energy loss has been proposed as a probe of the properties of dense matter formed in high-energy nuclear collisions . Based on a model of multiple scattering and induced radiation in QCD proposed by Gyulassy and Wang (GW), recent theoretical studies  show that a fast parton will lose a significant amount of energy via induced radiation when it propagates through a hot partonic matter. The most interesting feature of the result is the quadratical distance dependence of the total energy loss because of the non-Abelian nature of QCD radiation and the Landau-Pomercanchuk-Migdal (LPM)  interference. Such a quadratical distance dependence is also a consequence of the GW static color-screened potential model for multiple scattering where colors are screened but not confined. In this paper we will study parton multiple scattering inside a nucleus where colors are confined to the size of a nucleon. In this case parton propagation will certainly be different from that in a partonic matter and one should expect that the parton energy loss to be related to the nucleon size or confinement scale. Parton energy loss in $`eA`$ DIS has been studied before within various models of intranuclear scattering and the modification of hadronization . In this paper we will study the problem in the framework of multiple parton scattering in perturbative QCD (pQCD).

Unlike the situation in QED, the energy loss of a parton cannot be directly measured because partons are not the final experimentally observed particles. The total energy of a jet as traditionally defined by a cluster of hadrons in the phase space will not change much due to medium induced radiation because a jet so defined contains particles both from the leading parton and from the radiated gluons. This is particularly the case if multiple scattering and induced radiation do not dramatically  change the energy profile of the jet in phase-space. It is also virtually impossible to determine the jet energy event by event because of the large background and its fluctuation in heavy-ion collisions. One then has to resort to particle distributions within a jet and study the effect of parton energy loss by measuring the modification of the particle distribution due to multiple scattering and induced radiation . One such distribution is the fragmentation function of the produced parton, $`D_{ah}(z,\mu ^2)`$, where $`z`$ is the fractional energy of the parton $`a`$ carried by the produced particles $`h`$. Unlike the situation that has been considered in most of the recent theoretical studies  of parton energy loss of an on-shell parton during its propagation through medium, partons produced in hard processes are normally far off-shell characterized by the momentum scale $`\mu ^2`$. Final-state radiation of these off-shell partons in free space, such as in $`e^+e^{}`$ annihilation, leads to the $`\mu `$-dependence of the fragmentation functions as given by Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP)  QCD evolution equations,

$`{\displaystyle \frac{D_{qh}(z_h,\mu ^2)}{\mathrm{ln}\mu ^2}}`$ $`=`$ $`{\displaystyle \frac{\alpha _s(\mu ^2)}{2\pi }}{\displaystyle _{z_h}^1}{\displaystyle \frac{dz}{z}}\left[\gamma _{qqg}(z)D_{qh}(z_h/z,\mu ^2)+\gamma _{qgq}(z)D_{gh}(z_h/z,\mu ^2)\right],`$ (1)
$`{\displaystyle \frac{D_{gh}(z_h,\mu ^2)}{\mathrm{ln}\mu ^2}}`$ $`=`$ $`{\displaystyle \frac{\alpha _s(\mu ^2)}{2\pi }}{\displaystyle _{z_h}^1}{\displaystyle \frac{dz}{z}}\left[{\displaystyle \underset{q=1}{\overset{2n_f}{}}}\gamma _{gq\overline{q}}(z)D_{qh}(z_h/z,\mu ^2)+\gamma _{ggg}(z)D_{gh}(z_h/z,\mu ^2)\right]`$ (2)
where $`\gamma _{abc}(y)`$ are the splitting functions of the corresponding radiative processes . When a parton is produced in a medium, it will suffer multiple scattering and induced radiation that will then lead to medium modification of the DGLAP evolution of the parton fragmentation functions. In this paper we will derive such modified evolution equations for parton fragmentation functions in the simplest case of deeply inelastic $`eA`$ scattering (DIS). Multiple scattering and induced radiation suffered by the leading quark in this case will give rise to an additional term in the DGLAP evolution equations. As a consequence, the modified fragmentation functions become softer. This can be directly translated into the energy loss of the leading quark.

The study here is very similar to that by Luo, Qiu and Sterman (LQS)  on nuclear dependence of jet cross section in deeply inelastic scattering. The difference is that they considered large transverse momentum jet production whereas we will concentrate on soft gluon emission which is responsible for the DGLAP evolution equations of fragmentation functions. Depending on the fractional momentum carried by the second parton in the case of double parton scattering, one can categorize the processes into soft-hard, hard-hard and their corresponding interferences. In their study of large transverse momentum jet production, LQS only considered soft-hard and hard-hard processes, because the interference terms become negligible as we will show. Since we have to consider soft gluon emission, the interference terms are very important and they will cancel contributions from hard-hard and hard-soft processes in the limit of collinear emission (or zero transverse momentum). In this paper, we will include all these processes and treat them in the same manner. Utilizing the generalized factorization of higher-twist parton distributions, we find that each process in the double scattering probes different twist-four parton correlation of the nuclear medium. As a result of the sum of all contributions, the additional term in the modified DGLAP evolution equation is proportional to the combined twist-four matrix element,

$`{\displaystyle \frac{1}{4\pi }}{\displaystyle }`$ $`dy^{}dy_1^{}dy_2^{}\theta (y_2^{})\theta (y^{}y_1^{})(1e^{ix_Lp^+(y_1^{}y^{})})(1e^{ix_Lp^+y_2^{}})`$ (3)
$`e^{i(x_B+x_L)p^+y^{}}A|\overline{\psi }_q(0)\gamma ^+F^{\alpha +}(y_2^{})F_\alpha ^+(y_1^{})\psi _q(y^{})|A,`$ (4)
where

$$x_L=\mu ^2/2p^+q^{}z(1z),$$
(5)
with $`\mu `$ the typical transverse momentum or factorization scale and $`z`$ the fractional momentum of the emitted gluon. In this paper, we denote the four-momentum of the virtual photon and the target nucleon in DIS as

$`q`$ $`=`$ $`[Q^2/2q^{},q^{},\stackrel{}{0}_{}],`$ (6)
$`p`$ $`=`$ $`[p^+,0,\stackrel{}{0}_{}],`$ (7)
respectively. The Bjorken variable is then $`x_B=Q^2/2p^+q^{}`$.

In the above matrix element, one can identify $`1/x_Lp^+=2q^{}z(1z)/\mu ^2`$ as the formation time of the emitted gluons. For large formation time as compared to the nuclear size, the above matrix element vanishes, demonstrating a typical LPM interference effect. This is because the emitted gluon (with long formation time) and the leading quark are still a coherent system when they propagate through the nucleus. Additional scattering will not induce more gluon radiation, thus limiting the energy loss of the leading quark.

The reduction of phase space available for gluon radiation due to the LPM interference effect is critical for applying the LQS formalism to the problem in this paper. In the original LQS approach , the generalized factorization for processes with a large final transverse momentum $`\mathrm{}_T^2Q^2`$ allows one to consider the leading contribution in $`1/Q^2`$ which is enhanced by the nuclear size $`R_AA^{1/3}`$. For large $`Q^2`$ and $`A`$, higher-twist contribution from double parton rescattering that is proportional to $`\alpha _sR_A/Q^2`$ will then be the leading nuclear correction. One can neglect contributions from more than two parton rescattering. In deriving the modified fragmentation functions, we however have to take the leading logarithmic approximation in the limit $`\mathrm{}_T^2Q^2`$, where $`\mathrm{}_T`$ is the transverse momentum of the radiated gluon. Since the LPM interference suppresses gluon radiation whose formation time ($`\tau _fQ^2/\mathrm{}_T^2p^+`$) is larger than the nuclear size $`MR_A/p^+`$ in our chosen frame, $`\mathrm{}_T^2`$ should then have a minimum value of $`\mathrm{}_T^2Q^2/MR_AQ^2/A^{1/3}`$. Here $`M`$ is the nucleon mass. Therefore, the logarithmic approximation is still valid for large nuclei ($`MR_A1`$). In the meantime, the leading higher-twist contribution proportional to $`\alpha _sR_A/\mathrm{}_T^2\alpha _sR_A^2/Q^2`$ can still be small for large $`Q^2`$ so that one can neglect processes with more than two parton rescattering. The parameter for twist expansion of the fragmentation processes inside a nucleus is thus $`\alpha _sA^{2/3}/Q^2`$ as compared to $`\alpha _sA^{1/3}/Q^2`$ for processes with large final transverse momentum as studied by LQS . This is why the nuclear modification to the fragmentation function as derived in this paper depends quadratically on the nuclear size $`R_A`$.

Since our approach allows us to relate the energy loss of a fast parton to the parton correlations in the medium, the final results will be sensitive to the properties of the medium through which a produced parton propagates. In particular, the medium-dependence of the above matrix elements in a deconfined quark-gluon plasma will be very different from that in an ordinary nuclear medium. Therefore, the parton energy loss and its dependence on the medium size in hadronic matter are expected to differ from that in a quark-gluon plasma. This simply reflects the different parton correlations in different types of media.

The results of the present study were already reported in Ref. . We will provide a detailed derivation and discussion in this paper. The rest of this paper is organized as follows: In the next section, we give a brief overview of the framework of our study including generalized factorization of twist-four parton distributions in hard processes. In Section III we will discuss in detail the calculation of different contributions to soft gluon emission involving double parton scattering. To simplify the calculation, we will use the techniques of helicity amplitude with soft gluon approximation ($`z<<1`$). In Section IV, we will consider virtual corrections from unitary constraints. We will then define, in Section V, the effective parton fragmentation functions in deeply inelastic $`eA`$ collisions and derive the modified DGLAP evolution equations. We will also calculate the average energy loss suffered by the leading quark. Finally, in Section VI, we will discuss our results and their implications in other hard processes.

## II General Formalism

Consider the following semi-inclusive process in the deeply inelastic lepton-nucleus scattering,

$$e(L_1)+A(p)e(L_2)+h(\mathrm{}_h)+X,$$
(8)
where $`L_1`$ and $`L_2`$ are the four-momenta of the incoming and the outgoing leptons respectively, $`p`$ is the momentum per nucleon for the nucleus with the atomic number $`A`$, and $`\mathrm{}_h`$ is the observed hadron momentum. The momentum of the virtual photon ($`\gamma ^{}`$) is $`q=L_2L_1`$.

The differential cross section of semi-inclusive processes in DIS with an observed final state hadron $`\mathrm{}_h`$ can be expressed as

$$E_{L_2}E_\mathrm{}_h\frac{d\sigma _{\mathrm{DIS}}^h}{d^3L_2d^3\mathrm{}_h}=\frac{\alpha _{\mathrm{EM}}^2}{2\pi s}\frac{1}{Q^4}L_{\mu \nu }E_\mathrm{}_h\frac{dW^{\mu \nu }}{d^3\mathrm{}_h},$$
(9)
where $`s=(p+L_1)^2`$ is the total invariant mass of the lepton-nucleon system. The leptonic tensor $`L_{\mu \nu }`$ is given by

$$L_{\mu \nu }=\frac{1}{2}\mathrm{Tr}(\gamma L_1\gamma _\mu \gamma L_2\gamma _\nu ).$$
(10)
The semi-inclusive hadronic tensor $`E_\mathrm{}_hdW^{\mu \nu }/d^3\mathrm{}_h`$ is defined as,

$$E_\mathrm{}_h\frac{dW^{\mu \nu }}{d^3\mathrm{}_h}=\frac{1}{2}\underset{X}{}A|J^\mu (0)|X,hX,h|J^\nu (0)|A2\pi \delta ^4(q+pp_X\mathrm{}_h),$$
(11)
where $`_X`$ runs over all possible intermediate states and $`J_\mu `$ is the hadronic electromagnetic (EM) current, $`J_\mu =e_q\overline{\psi }_q\gamma _\mu \psi _q`$.

In the parton model with collinear factorization approximation, one can in general factorize the semi-inclusive cross section into a parton fragmentation function and the partonic cross section. By factorizing the parton fields out of the EM current as

$$J^\mu (0)\underset{c,\alpha }{}\stackrel{~}{J}_{c,\alpha }^\mu (0)\widehat{\varphi }_c^\alpha (0)\underset{c,\alpha }{}\widehat{\varphi }_c^\alpha (0)\stackrel{~}{J}_{c,\alpha }^\mu (0),$$
(12)
one can express the semi-inclusive hadronic tensor as

$$\frac{dW^{\mu \nu }}{dz_h}=\underset{\alpha ,\beta ,c}{}z_h^2\widehat{d}_{ch}^{\alpha \beta }(z_h,\mathrm{}_h)\underset{X}{}A|\stackrel{~}{J}_{c,\alpha }^\mu (0)|XX|\stackrel{~}{J}_{c,\beta }^\nu (0)|A2\pi \delta (\mathrm{}_c^2),$$
(13)
as illustrated by Fig. 1, where $`\alpha ,\beta `$ are indices (spinor for quarks and Lorentz for gluons) of the partonic field of species $`c`$ with momentum $`\mathrm{}_c`$, and $`\widehat{d}_{ch}^{\alpha \beta }(z_h,\mathrm{}_h)`$, related to the fragmentation function, is defined as

$$\widehat{d}_{ch}^{\alpha \beta }(z_h,\mathrm{}_h)\underset{S}{}\frac{d^4\mathrm{}_c}{(2\pi )^4}d^4y0|\widehat{\varphi }_c^\alpha (0)|h,Sh,S|\widehat{\varphi }_c^\beta (y)|0e^{i\mathrm{}_cy}\delta (z_h\frac{\mathrm{}_h^{}}{\mathrm{}_c^{}}).$$
(14)
After such factorization, one only needs to calculate the cross sections of partonic hard processes. The fragmentation functions defined in terms of parton matrix elements in the above equation are nonperturbative and can only be obtained via experimental measurements. In addition, one can also factorize out the parton distributions inside nuclei (or nucleons). Such factorization has been proven to all orders and in the leading twist .

In the infinite-momentum frame as defined by Eq. (6), the dominant momentum component of the leading parton in a single jet event is the minus component. Therefore, we define the momentum fraction carried by a produced hadron as $`z_h=\mathrm{}_h^{}/\mathrm{}_c^{}`$.

### A Leading Twist Contributions

The leading-twist contribution to DIS to the lowest order ($`𝒪(\alpha _s^0)`$) comes from a single hard $`\gamma ^{}+q`$ scattering as illustrated in Fig. 2. The semi-inclusive hadronic tensor $`dW_{\mu \nu }/dz_h`$ can be factorized as

$$\frac{dW_{\mu \nu }^{S(0)}}{dz_h}=\underset{q}{}𝑑xf_q^A(x)H_{\mu \nu }^{(0)}(x,p,q)D_{qh}(z_h),$$
(15)
where $`f_q^A(x)`$ is the quark distribution as defined by

$`f_q^A(x)`$ $``$ $`{\displaystyle \frac{d^4k}{(2\pi )^4}d^4ye^{iky}\delta (x\frac{k^+}{p^+})A|\overline{\psi }_q(0)\frac{\gamma ^+}{2p^+}\psi _q(y)|A}`$ (16)
$`=`$ $`{\displaystyle \frac{dy^{}}{2\pi }e^{ixp^+y^{}}\frac{1}{2}A|\overline{\psi }_q(0)\gamma ^+\psi _q(y^{})|A}.`$ (17)
Here quark fragmentation function $`D_{qh}(z_h)`$ is defined as

$`D_{qh}(z_h)`$ $``$ $`{\displaystyle \frac{z_h^3}{2}}d_{qh}(z_h,\mathrm{}_h){\displaystyle \frac{z_h^3}{4\mathrm{}_h^{}}}\mathrm{Tr}[\gamma ^{}\widehat{d}_{qh}(z_h,\mathrm{}_h)],`$ (18)
$`\widehat{d}_{qh}^{\alpha \beta }(z_h,\mathrm{}_h)`$ $``$ $`{\displaystyle \underset{S}{}}{\displaystyle \frac{d^4\mathrm{}_q}{(2\pi )^4}d^4ye^{i\mathrm{}_qy}\delta (z_h\frac{\mathrm{}_h^{}}{\mathrm{}_q^{}})0|\psi _q^\beta (0)|h,Sh,S|\overline{\psi }_q^\alpha (y)|0}`$ (19)
$`=`$ $`{\displaystyle \underset{S}{}}{\displaystyle \frac{\mathrm{}_h^{}}{z_h^2}}{\displaystyle \frac{dy^+}{2\pi }e^{i\mathrm{}_h^{}y^+/z_h}0|\psi _q^\beta (0)|h,Sh,S|\overline{\psi }_q^\alpha (y^+)|0}.`$ (20)
The hard part of $`\gamma ^{}+q`$ partonic scattering is

$`H_{\mu \nu }^{(0)}`$ $`=`$ $`e_q^2{\displaystyle \frac{1}{2}}\mathrm{Tr}(\gamma p\gamma _\mu \gamma (q+xp)\gamma _\nu )(2\pi )\delta [(q+xp)^2]`$ (21)
$`=`$ $`4\pi e_q^2\left[x_Be_{\mu \nu }^L{\displaystyle \frac{1}{2}}e_{\mu \nu }^T\right]\delta (xx_B),`$ (22)
where momentum conservation gives $`\mathrm{}_q=xp+q`$ and by definition $`\mathrm{}_h^{}=z_h\mathrm{}_q^{}=z_hq^{}`$. The transverse and longitudinal tensors are defined as

$`e_{\mu \nu }^T`$ $`=`$ $`g_{\mu \nu }{\displaystyle \frac{q_\mu q_\nu }{q^2}},`$ (23)
$`e_{\mu \nu }^L`$ $`=`$ $`{\displaystyle \frac{1}{pq}}[p_\mu {\displaystyle \frac{pq}{q^2}}q_\mu ][p_\nu {\displaystyle \frac{pq}{q^2}}q_\nu ].`$ (24)
Note that a sum and an average over color indices are implied in the quark distribution and fragmentation function, respectively. In Eq. (18), we used the collinear approximation in the unpolarized fragmentation functions,

$`\widehat{d}_{qh}(z_h,\mathrm{}_h)`$ $``$ $`{\displaystyle \frac{1}{2}}d_{qh}(z_h,\mathrm{}_h)\gamma \mathrm{}_h+\mathrm{},`$ (25)
$``$ $`{\displaystyle \frac{z_h}{2}}d_{qh}(z_h,\mathrm{}_h)\gamma \mathrm{}_q+\mathrm{},`$ (26)
while other higher-twist terms are neglected. We will use such a collinear approximation throughout this paper even in the case of double scattering. We will also neglect other higher-twist contributions to the fragmentation function that are independent of the nuclear size.

At the next-to-leading order $`𝒪(\alpha _s)`$, the dominant (in leading log approximation) real radiative contribution to the fragmentation process in an axial gauge ($`A^{}=0`$) comes from the final state radiation. Since we are only interested in the collinear behavior of the radiative corrections in order to study the evolution equation of the fragmentation functions, we will keep only the leading (divergent) contribution when the gluon’s transverse momentum vanishes(with respect to $`q+xp`$), i.e., $`\mathrm{}_T0`$. Using the same collinear approximation, one finds the leading contribution from the radiative correction (for $`\mathrm{}_T^2`$ up to a factorization scale $`\mu ^2`$) to the quark fragmentation process,

$$\frac{dW_{\mu \nu }^{S(1)q}}{dz_h}=\underset{q}{}𝑑xf_q^A(x)_{z_h}^1\frac{dz}{z}D_{qh}(z_h/z)H_{\mu \nu }^{(1)q}(x,p,q,z),$$
(27)
where $`z_h=\mathrm{}_h^{}/q^{}`$ and $`z=\mathrm{}_q^{}/q^{}`$ are the fractional momentum carried by hadron and the final quark, respectively. The factorized form in Eq. (27) is quite general for all processes involving gluon radiation in the final state. It is then enough to calculate the partonic hard part $`H_{\mu \nu }^{(1)q}(x,p,q,z)`$ as illustrated in Fig. 3,

$$H_{\mu \nu }^{(1)q}(x,p,q,z)=H_{\mu \nu }^{(0)}(x,p,q)_0^{\mu ^2}\frac{d\mathrm{}_T^2}{\mathrm{}_T^2}\frac{\alpha _s}{2\pi }C_F\frac{1+z^2}{1z},$$
(28)
where the color factor is $`C_F=(N_c^21)/2N_c=4/3`$. Since we neglect the non-leading-log terms, the tensor structure of $`H_{\mu \nu }^{(1)q}`$ remains the same as $`H_{\mu \nu }^{(0)}`$ in Eq. (22). Similarly, if the final hadron comes from the gluon fragmentation, the contribution is

$$\frac{dW_{\mu \nu }^{S(1)g}}{dz_h}=\underset{q}{}𝑑xf_q^A(x)H_{\mu \nu }^{(0)}(x,p,q)_0^{\mu ^2}\frac{d\mathrm{}_T^2}{\mathrm{}_T^2}\frac{\alpha _s}{2\pi }_{z_h}^1\frac{dz}{z}C_F\frac{1+(1z)^2}{z}D_{gh}(z_h/z).$$
(29)
Here $`z=\mathrm{}_g^{}/q^{}`$ and $`D_{gh}(z_h)`$ is the gluon fragmentation function defined as

$`D_{gh}(z_h)`$ $``$ $`{\displaystyle \frac{z_h^2}{2}}\epsilon _{\mu \nu }(\mathrm{}_g)d_{gh}^{\mu \nu }(z_h,\mathrm{}_h)`$ (30)
$`=`$ $`{\displaystyle \frac{z_h^2}{2\mathrm{}_h^{}}}{\displaystyle \underset{S}{}}{\displaystyle \frac{dy^+}{2\pi }e^{i\mathrm{}_h^{}y^+/z_h}0|F^\mu (0)|S,hS,h|F_\mu ^{}(y^+)|0},`$ (31)
where $`d_{gh}^{\mu \nu }(z_h,\mathrm{}_h)`$ is given by Eq. (14) for gluon fields and

$$\epsilon _{\mu \nu }(\mathrm{}_g)=\underset{\lambda =1,2}{}\epsilon _\mu (\mathrm{}_g,\lambda )\epsilon _\nu (\mathrm{}_g,\lambda ),$$
(32)
with $`\epsilon _\mu (\mathrm{}_g,\lambda )`$ being the polarization vector of a gluon in an axial gauge.

There are both infrared and collinear divergences in Eq. (28) and Eq. (29). The infrared divergences come from the phase space where the gluon’s fractional momentum goes to zero. These divergences will be canceled by the virtual corrections in the quark self-energy diagrams as we will discuss later. The collinear divergences when $`\mathrm{}_T0`$ will be absorbed into the renormalized fragmentation functions and are responsible for the evolution of the quark fragmentation function with respect to the factorization scale $`\mu `$. Summing up all the contributions from Eqs. (15), (27) and (29), one has

$`{\displaystyle \frac{dW_{\mu \nu }^S}{dz_h}}`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle 𝑑xf_q^A(x)H_{\mu \nu }^{(0)}(x,p,q)D_{qh}(z_h,\mu ^2)},`$ (33)
$`D_{qh}(z_h,\mu ^2)`$ $``$ $`D_{qh}(z_h)+{\displaystyle _0^{\mu ^2}}{\displaystyle \frac{d\mathrm{}_T^2}{\mathrm{}_T^2}}{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle _{z_h}^1}{\displaystyle \frac{dz}{z}}[C_F{\displaystyle \frac{1+z^2}{1z}}D_{qh}(z_h/z)`$ (34)
$`+`$ $`C_F{\displaystyle \frac{1+(1z)^2}{z}}D_{gh}(z_h/z)]+\mathrm{virtual}\mathrm{corrections},`$ (35)
where the renormalized quark fragmentation function $`D_{qh}(z_h,\mu ^2)`$ satisfies the DGLAP evolution equation in Eq. (1) after inclusion of the virtual corrections. The virtual corrections have to be included to ensure unitarity and infrared safety of the final result. We will discuss the virtual contributions later in this paper.

### B Higher Twist Contributions

In a nuclear medium, the outgoing quark in DIS may experience additional scattering with other partons from the nucleus. The additional scattering may induce additional gluon radiation and cause the leading quark to lose energy. Such induced gluon radiation will effectively give rise to additional terms in the evolution equation leading to modification of the fragmentation functions in a medium. These additional terms from multiple scattering are always non-leading-power contributions. There are many next-leading-twist contributions to the semi-inclusive hadronic tensor in DIS. One type of contribution involves non-leading-twist matrix elements of parton distributions in the nucleus which in general are related to two-parton correlations. Since two-parton correlations in a nucleus can involve partons from different nucleons in the nucleus, they are proportional to the size of the nucleus and thus are enhanced by a nuclear factor $`A^{1/3}`$ as compared to two-parton correlations in a nucleon. In this paper we will consider only those non-leading-twist contributions that are enhanced by the nuclear factor. We will neglect those contributions that are not enhanced by the nuclear medium. They are in general not related to multiple parton scattering in the nucleus.

In the LQS generalized factorization scheme for processes in nuclei with large final transverse momentum ($`\mathrm{}_T^2Q^2`$), contributions from double parton scattering are leading in $`\alpha _sA^{1/3}/Q^2`$. One can similarly apply the generalized factorization to the fragmentation processes in nuclei thanks to the LPM interference effect in gluon radiation induced by multiple parton scattering. As we have argued in the Introduction and will demonstrate later in this paper, the formation time of gluon radiation due to the LPM interference requires the radiated gluon to have a minimum transverse momentum $`\mathrm{}_T^2Q^2/MR_AQ^2/A^{1/3}`$. The nuclear corrections to the fragmentation function due to double parton scattering will then be in the order of $`\alpha _sA^{1/3}/\mathrm{}_T^2\alpha _sA^{2/3}/Q^2`$. For large values of $`A`$ and $`Q^2`$, these corrections are leading and yet the requirement $`\mathrm{}_T^2Q^2`$ for the logarithmic approximation in deriving the modified DGLAP evolution is still valid. Therefore, for large $`Q^2`$ and $`A`$, we can neglect processes with more than two parton rescattering and consider only double scattering.

We will work in the LQS framework  in which the twist-four parton distributions are factorizable. We will simply apply such factorization of twist-four parton distributions to the study of semi-inclusive processes in DIS. In general, the twist-four contributions can be expressed as the convolution of partonic hard parts and four-parton matrix elements . At the lowest order (processes without gluon radiation) in this framework, rescattering with collinear gluons gives the eikonal contribution to the gauge-invariant leading-twist and lowest-order result in Eq. (15), assuming collinear factorization of the quark fragmentation function. As we will show later, rescattering with another quark in the leading order will contribute to the semi-inclusive cross section but not to the renormalization equation because there is no collinear divergency. Therefore, we first concentrate on processes involving rescattering with gluons. For next-to-leading order processes (with gluon radiation) involving a secondary scattering with another gluon from the nucleus, shown in Fig. 4 as an example, the double scattering contributions to $`dW_{\mu \nu }/dz_h`$ can be expressed as

$`{\displaystyle \frac{dW_{\mu \nu }^D}{dz_h}}`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle _{z_h}^1}{\displaystyle \frac{dz}{z}}D_{qh}(z_h/z){\displaystyle \frac{dy^{}}{2\pi }𝑑y_1^{}𝑑y_2^{}\frac{d^2y_T}{(2\pi )^2}d^2k_T\overline{H}_{\mu \nu }^D(y^{},y_1^{},y_2^{},k_T,p,q,z)};`$ (36)
$`\times `$ $`e^{i\stackrel{}{k}_T\stackrel{}{y}_T}{\displaystyle \frac{1}{2}}A|\overline{\psi }_q(0)\gamma ^+A^+(y_2^{},0_T)A^+(y_1^{},y_T)\psi _q(y^{})|A.`$ (37)
Here $`\overline{H}_{\mu \nu }^D(y^{},y_1^{},y_2^{},k_T,p,q,z)`$ is the Fourier transform of the partonic hard part $`\stackrel{~}{H}_{\mu \nu }(x,x_1,x_2,k_T,p,q,z)`$ in momentum space,

$`\overline{H}_{\mu \nu }^D(y^{},y_1^{},y_2^{},k_T,p,q,z)`$ $`=`$ $`{\displaystyle 𝑑x\frac{dx_1}{2\pi }\frac{dx_2}{2\pi }e^{ix_1p^+y^{}+ix_2p^+y_1^{}+i(xx_1x_2)p^+y_2^{}}}`$ (38)
$`\times \stackrel{~}{H}_{\mu \nu }^D(x,x_1,x_2,k_T,p,q,z),`$ (39)
where $`k_T`$ is the relative transverse momentum carried by the second parton in the double scattering. Values of the momentum fractions $`x,x_1,`$ and $`x_2`$ are fixed by $`\delta `$-functions and poles in the partonic hard part. They normally depend on $`k_T`$.

In order to pick up the next-leading-twist contribution, we expand the partonic hard part around $`k_T=0`$,

$`\overline{H}_{\mu \nu }^D(y^{},y_1^{},y_2^{},k_T,p,q,z)`$ $`=`$ $`\overline{H}_{\mu \nu }^D(y^{},y_1^{},y_2^{},k_T=0,p,q,z)`$ (40)
$`+`$ $`{\displaystyle \frac{\overline{H}_{\mu \nu }^D}{k_T^\alpha }}|_{k_T=0}k_T^\alpha +{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2\overline{H}_{\mu \nu }^D}{k_T^\alpha k_T^\beta }}|_{k_T=0}k_T^\alpha k_T^\beta +\mathrm{}.`$ (41)
This is known as collinear expansion , On the right-hand-side of Eq. (41), the first term gives the eikonal contribution to the leading-twist results. It does not correspond to the physical double scattering, but simply makes the matrix element in a single scattering gauge invariant. The second term for unpolarized initial and final states vanishes after being integrated over $`k_T`$. The third term will give a finite contribution to the double scattering process. Substituting Eq. (41) into Eq. (37), and integrating over $`d^2k_T`$ and $`d^2y_T`$, we obtain

$`{\displaystyle \frac{dW_{\mu \nu }^D}{dz_h}}`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle _{z_h}^1}{\displaystyle \frac{dz}{z}}D_{qh}(z_h/z){\displaystyle \frac{dy^{}}{2\pi }𝑑y_1^{}𝑑y_2^{}\frac{1}{2}A|\overline{\psi }_q(0)\gamma ^+F_\sigma ^+(y_2^{})F^{+\sigma }(y_1^{})\psi _q(y^{})|A}`$ (42)
$`\times `$ $`\left({\displaystyle \frac{1}{2}}g^{\alpha \beta }\right)\left[{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{k_T^\alpha k_T^\beta }}\overline{H}_{\mu \nu }^D(y^{},y_1^{},y_2^{},k_T,p,q,z)\right]_{k_T=0},`$ (43)
where $`k_T^\alpha A^+k_T^\beta A^+`$ are converted into field strength $`F^{\alpha +}F^{\beta +}`$ by partial integrations.

## III Quark-gluon double scattering

In this section we calculate the hard part of parton rescattering with gluons. We consider processes with a quark and a gluon in the initial state. In these processes, there is first a hard photon-quark scattering. The produced parton from the first hard scattering then has a second scattering with another initial gluon from the nucleus. We refer to such processes as quark-gluon double scattering processes. We defer the calculation of virtual corrections to the latter part of this paper. There are a total of 23 cut diagrams corresponding to 3 different double scattering processes, their interferences, and the interferences between single scattering and 7 different triple scattering processes (see Appendix). We use the process shown in Fig. 4 with three different cuts as an example to outline the steps that are essential to evaluate the hard part of double scattering.

Since we are interested in the leading log behavior of the fragmentation functions, we will only retain the leading contributions in the limit of $`\mathrm{}_T0`$ ($`\mathrm{}_T`$ being gluons transverse momentum) and neglect contributions that are finite when $`\mathrm{}_T=0`$. This corresponds to the leading log approximation for $`\mathrm{}_T^2Q^2`$.

### A Hard vs Soft Rescattering

Following the collinear approximation involving leading twist-4 parton distributions, as outlined in Ref. , the hard partonic part in Fig. 4 with a central cut is

$`\overline{H}_{C\mu \nu }^D(y^{},y_1^{},y_2^{},k_T,p,q,z)`$ $`=`$ $`{\displaystyle 𝑑x\frac{dx_1}{2\pi }\frac{dx_2}{2\pi }e^{ix_1p^+y^{}+ix_2p^+y_1^{}+i(xx_1x_2)p^+y_2^{}}\frac{d^4\mathrm{}}{(2\pi )^4}}`$ (44)
$`\times `$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}\left[p\gamma \gamma _\mu p^\sigma p^\rho \widehat{H}_{\sigma \rho }\gamma _\nu \right]2\pi \delta _+(\mathrm{}^2)\delta (1z{\displaystyle \frac{\mathrm{}^{}}{q^{}}}).`$ (45)
Here,

$`\widehat{H}_{\sigma \rho }`$ $`=`$ $`{\displaystyle \frac{C_F}{2N_c}}g^4{\displaystyle \frac{\gamma (q+x_1p)}{(q+x_1p)^2iϵ}}\gamma _\alpha {\displaystyle \frac{\gamma (q+x_1p\mathrm{})}{(q+x_1p\mathrm{})^2iϵ}}\gamma _\sigma \gamma \mathrm{}_q\gamma _\rho `$ (46)
$`\times `$ $`\epsilon ^{\alpha \beta }(\mathrm{}){\displaystyle \frac{\gamma (q+xp\mathrm{})}{(q+xp\mathrm{})^2+iϵ}}\gamma _\beta {\displaystyle \frac{\gamma (q+xp)}{(q+xp)^2+iϵ}}\mathrm{\hspace{0.17em}\hspace{0.17em}2}\pi \delta _+(\mathrm{}_q^2),`$ (47)
and

$$\epsilon ^{\alpha \beta }(\mathrm{})=g^{\alpha \beta }+\frac{n^\alpha \mathrm{}^\beta +n^\beta \mathrm{}^\alpha }{n\mathrm{}}n^2\frac{\mathrm{}^\alpha \mathrm{}^\beta }{(n\mathrm{})^2}$$
(48)
is the polarization tensor of a gluon propagator in an axial gauge ($`nA=0`$) with $`n=[1,0^{},\stackrel{}{0}_{}]`$ and $`\mathrm{}`$, $`\mathrm{}_q=q+(x_1+x_2)p+k_T\mathrm{}`$ are the 4-momenta carried by the gluon and the final quark, respectively. $`z=\mathrm{}_q^{}/q^{}`$ is the fraction of longitudinal momentum (the large minus component) carried by the final quark.

To simplify the calculation of the trace part of $`\widehat{H}_{\sigma \rho }`$ and extract the leading contribution in the limit $`\mathrm{}_T0`$ and $`k_T0`$, we again use the collinear approximation,

$$p^\sigma \widehat{H}_{\sigma \rho }p^\rho \gamma \mathrm{}_q\frac{1}{4\mathrm{}_q^{}}\mathrm{Tr}\left[\gamma ^{}p^\sigma \widehat{H}_{\sigma \rho }p^\rho \right].$$
(49)
The two longitudinal components in the $`\mathrm{}`$-integration are fixed by the two $`\delta `$-functions in Eq. (45) and the remainder becomes $`d\mathrm{}_T^2/(4\pi )^2`$. To carry out the integration over $`x,x_1`$ and $`x_2`$, we rewrite the $`\delta `$-function in Eq. (47) as

$$\delta _+(\mathrm{}_q^2)=\frac{1}{2p^+q^{}z}\delta (x_1+x_2x_Lx_Dx_B),$$
(50)
where $`x_B=Q^2/2p^+q^{}`$ is the Bjorken variable and

$`x_L`$ $`=`$ $`{\displaystyle \frac{\mathrm{}_T^2}{2p^+q^{}z(1z)}},x_D={\displaystyle \frac{k_T^22\stackrel{}{k}_T\stackrel{}{\mathrm{}}_T}{2p^+q^{}z}}.`$ (51)
With the help of the on-shell condition $`\mathrm{}_q^2=0`$ and $`\mathrm{}^2=0`$ in the $`\delta `$-functions, one can also simplify the variables of the propagators in Eq. (47):

$`(q+xp)^2`$ $`=`$ $`2p^+q^{}(xx_B),(q+xp\mathrm{})^2=2p^+q^{}z(xx_Lx_B),`$ (52)
$`(q+x_1p)^2`$ $`=`$ $`2p^+q^{}(x_1x_B),(q+x_1p\mathrm{})^2=2p^+q^{}z(x_1x_Lx_B).`$ (53)
The integration over $`x`$, $`x_1`$ and $`x_2`$ can now be cast in the form

$`I_C(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,p,q,z)`$ $`=`$ $`{\displaystyle 𝑑x\frac{dx_1}{2\pi }\frac{dx_2}{2\pi }\frac{e^{ix_1p^+y^{}+ix_2p^+y_1^{}+i(xx_1x_2)p^+y_2^{}}}{(x_1x_Biϵ)(x_1x_Lx_Biϵ)}}`$ (54)
$``$ $`{\displaystyle \frac{\delta (x_1+x_2x_Lx_Dx_B)}{(xx_Lx_B+iϵ)(xx_B+iϵ)}}.`$ (55)
Two of the above integrations can be carried out by contour integration. There are four possible poles in the denominators and two of them are used to obtain the residues of the contour integration. There are four possible choices for the pair of poles representing subprocesses with different kinematics. When we choose the poles

$$x=x_B+x_L,x_1=x_B+x_L,$$
(56)
which is the momentum fraction carried by the initial quark, the momentum fraction carried by the initial gluon in the rescattering is $`x_2=x_D`$, which vanishes when $`k_T0`$ according to the definition in Eq. (51). In this case, the final gluon is produced via the final state radiation of the hard photon-quark scattering. After the final state radiation, the quark becomes on-shell and encounters another scattering with a very soft gluon. This subprocess corresponds to rescattering with a soft gluon after a hard collision and is normally referred to as a hard-soft process. On the other hand, if we choose the poles,

$$x=x_B,x_1=x_B,$$
(57)
the momentum fraction carried by the initial gluon in the rescattering is $`x_2=x_L+x_D`$, which is finite ($`x_L`$) even when $`k_T0`$. We call this type of rescattering hard and refer the corresponding double scattering as double-hard processes. According to such a choice of a pair of poles, the quark becomes on-shell immediately after the hard photon-quark scattering. The gluon radiation is actually induced by the hard secondary quark-gluon scattering. It is thus produced by the initial state radiation of the secondary scattering. The other two combinations of the poles,

$`x`$ $`=`$ $`x_B+x_L,x_1=x_B;`$ (58)
$`x_1`$ $`=`$ $`x_B+x_L,x=x_B,`$ (59)
represent the interferences between hard-soft and double-hard processes. Using the relations imposed by the $`\delta `$-function, $`x_2=x_B+x_L+x_Dx_1`$, and the momentum conservation, $`x_3=x_1+x_2x`$, one can verify that there is a momentum transfer $`x_L`$ between the initial quark and gluon fields. If we consider the quark and gluon come from different nucleons inside the nucleus, this means that there is a momentum transfer between different nucleons in these interferences processes. Such an observation will become important when we study later the parton matrix elements involved.

One can easily find out the residues of the contour integrations with the above four possible choices of poles. After making a change of variable $`x_1+x_2x_Lx_Dx`$, one finds

$`I_C(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,p,q,z)`$ $`=`$ $`{\displaystyle 𝑑x\frac{\delta (xx_B)}{x_L^2}\overline{I}_C(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,x,p,q,z)};`$ (60)
$`\overline{I}_C(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,x,p,q,z)`$ $`=`$ $`e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\theta (y_2^{})\theta (y^{}y_1^{})`$ (61)
$`\times `$ $`(1e^{ix_Lp^+y_2^{}})(1e^{ix_Lp^+(y^{}y_1^{})}).`$ (62)
In the above contributions, the interferences between hard-soft and double-hard processes have the negative sign relative to the hard-soft and double-hard contributions. Previous studies  have considered only hard-soft and double-hard processes, but neglected the interferences. For large $`\mathrm{}_T`$ or $`x_L`$, this is justified because the interference terms will vanish due to the oscillation of the phase factors in the integration of $`y_1^{}`$ and $`y_2^{}`$ over the size of the whole nucleus. Note that in the double-hard term the phase factor $`\mathrm{exp}[ix_Lp^+(y_1^{}y_2^{})]`$ does not vanish after integration over $`y_1^{}`$ and $`y_2^{}`$, because color confinement requires that $`|y_1^{}y_2^{}|`$ remain the size of a nucleon. We will discuss other consequences of color confinement later in this paper. At the $`\mathrm{}_T0`$ limit, one can no longer neglect the interference terms because they exactly cancel the contributions from the hard-soft and double-hard processes. This was pointed out in Ref. . However, our final results, for the first time explicitly demonstrate such cancellation. In the latter part of this section, we will discuss the physical implications of such cancellation and its connection with the LPM interference effect. It will help us to reorganize our final results.

Because of the collinear (or leading log) approximation, the tensor structure of the double scattering contributions is generally the same as in the leading order single scattering,

$$\overline{H}_{\mu \nu }^D(y^{},y_1^{},y_2^{},k_T,p,q,z)=𝑑xH_{\mu \nu }^{(0)}(x,p,q)\overline{H}^D(y^{},y_1^{},y_2^{},k_T,x,p,q,z).$$
(63)
After completing the trace in Eq. (49) and using the results of the contour integration in Eq. (62), we have the contributions from Fig. 4 with a central cut,

$`\overline{H}_C^D(y^{},y_1^{},y_2^{},k_T,x,p,q,z)`$ $`=`$ $`{\displaystyle \frac{d\mathrm{}_T^2}{\mathrm{}_T^2}\frac{\alpha _s}{2\pi }C_F\frac{1+z^2}{1z}}`$ (64)
$`\times `$ $`{\displaystyle \frac{2\pi \alpha _s}{N_c}}\overline{I}_C(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,x,p,q,z),`$ (65)
where $`\overline{I}_C`$ is given in Eq. (62) and $`H_{\mu \nu }^{(0)}`$ in Eq. (22).

In addition to the central-cut diagram, we also need to consider asymmetrical-cut diagrams in Fig. 4 that represent interferences between single and triple scattering. The trace part is exactly the same as in the central-cut diagram. The differences are that one of the quark propagators is replaced by the cut propagator in the central-cut diagram. In addition, each of the asymmetrical-cut diagrams has only two possible pairs of poles, whereas there are four choices in the central-cut diagram. After similar contour integration, we have for the left-cut diagram

$`I_L(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,p,q,z)`$ $`=`$ $`{\displaystyle 𝑑x\frac{dx_1}{2\pi }\frac{dx_2}{2\pi }\frac{e^{ix_1p^+y^{}+ix_2p^+y_1^{}+i(xx_1x_2)p^+y_2^{}}}{(x_1x_Biϵ)(x_1x_Lx_Biϵ)}}`$ (66)
$``$ $`{\displaystyle \frac{\delta (xx_Bx_L)}{(x_1+x_2x_Bx_Lx_Diϵ)(xx_B+iϵ)}}`$ (67)
$`=`$ $`{\displaystyle 𝑑x\frac{\delta (xx_B)}{x_L^2}\overline{I}_L(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,x,p,q,z)};`$ (68)
$`\overline{I}_L(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,x,p,q,z)`$ $`=`$ $`e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\theta (y_1^{}y_2^{})\theta (y^{}y_1^{})`$ (69)
$`\times `$ $`(1e^{ix_Lp^+(y^{}y_1^{})}).`$ (70)
In the second step of the above equation we have made a change of variable $`xx_Lx`$. On the left-hand side of the left-cut diagram, the kinematics of the final state gluon radiation in the single hard photon-quark scattering requires the quark line to be off-shell, i.e., $`x=x_B`$. Similarly, on the right-hand side, only one of the quark propagators connecting the gluon radiation vertex can be on shell. This leads to two possible pair of poles in the propagators,

$`x_1`$ $`=`$ $`x_B+x_L,x_2=x_D;`$ (71)
$`x_1`$ $`=`$ $`x_B,x_2=x_L+x_D,`$ (72)
giving rise to the two contributions in Eq. (70). They correspond to soft and hard first quark-gluon rescattering, respectively, in the triple scattering process. In both cases, $`x_3=x_D`$. So the second quark-gluon rescattering in the triple scattering process on the right-hand side of the left-cut diagram is always soft.

Similarly, one can also obtain for the right-cut diagram of Fig. 4,

$`\overline{I}_R(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,x,p,q,z)`$ $`=`$ $`e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\theta (y_2^{})\theta (y_2^{}y_1^{})`$ (73)
$`\times `$ $`(1e^{ix_Lp^+y_2^{}}).`$ (74)
Other coefficients of the asymmetrical-cut diagrams are exactly the same as in the central-cut diagram given in Eq. (65).

Since QCD is a non-Abelian gauge theory, rescattering can also happen between radiated and target gluons. Shown in Fig. 5 is an example that involves gluon-gluon rescattering. The contribution to the hard part of double scattering from this central cut diagram is

$`\overline{H}_{\mathrm{Fig}.\text{5}}^D(y^{},y_1^{},y_2^{},k_T,x,p,q,z)`$ $`=`$ $`{\displaystyle \frac{d\mathrm{}_T^2}{(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}\frac{\alpha _s}{2\pi }C_A\frac{1+z^2}{1z}}`$ (75)
$`\times `$ $`{\displaystyle \frac{2\pi \alpha _s}{N_c}}\overline{I}_{\mathrm{Fig}.\text{5}}(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,x,p,q,z),`$ (76)
$`\overline{I}_{\mathrm{Fig}.\text{5}}(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,x,p,q,z)`$ $`=`$ $`e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\theta (y_2^{})\theta (y^{}y_1^{})`$ (77)
$`\times `$ $`[e^{ix_Dp^+y_2^{}/(1z)}e^{ix_Lp^+y_2^{}}]`$ (78)
$`\times `$ $`[e^{ix_Dp^+(y^{}y_1^{})/(1z)}e^{ix_Lp^+(y^{}y_1^{})}],`$ (79)
which has a structure very similar to the contribution in Eq. (62) and (65) from the central-cut diagram in Fig. 4. The four different terms correspond similarly to the hard-soft, double-hard scattering and their interferences. In the hard-soft process, the emitted gluon encounters another scattering with a very soft target gluon after it is produced in the final state radiation of the hard photon-quark scattering. Note that the initial quark has momentum fraction $`x=x_B+x_L+x_D/(1z)`$ while the soft gluon carries momentum fraction $`x_2=x_Dx_D/(1z)`$. In the double-hard process, the gluon is induced by the secondary hard scattering, similar to the quark-gluon rescattering, but in this case via a three-gluon vertex. The $`k_T`$ dependence of the cross section is typical of the gluon-gluon interaction.

As another example, we illustrate the interference between quark-gluon and gluon-gluon rescattering processes in Fig. 6. The contribution from the central-cut diagram is,

$`\overline{H}_{\mathrm{Fig}.\text{6}}^D(y^{},y_1^{},y_2^{},k_T,x,p,q,z)`$ $`=`$ $`{\displaystyle 𝑑\mathrm{}_T^2\frac{\stackrel{}{\mathrm{}_T}(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}\frac{\alpha _s}{2\pi }\frac{C_A}{2}\frac{1+z^2}{1z}}`$ (80)
$`\times `$ $`{\displaystyle \frac{2\pi \alpha _s}{N_c}}\overline{I}_{\mathrm{Fig}.\text{6}}(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,x,p,q,z),`$ (81)
$`\overline{I}_{\mathrm{Fig}.\text{6}}(y^{},y_1^{},y_2^{},\mathrm{}_T,k_T,x,p,q,z)`$ $`=`$ $`e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\theta (y_2^{})\theta (y^{}y_1^{})`$ (82)
$`\times `$ $`[e^{ix_Dp^+y_2^{}/(1z)}e^{ix_Lp^+y_2^{}}][1e^{ix_Lp^+(y^{}y_1^{})}].`$ (83)
Notice that all interference terms in Eqs. (70-74) and (83) have the opposite sign as compared to symmetrical diagrams. The calculation of the other diagrams with all possible cuts is tedious but similarly straightforward. The splitting functions in all quark-gluon double scattering and their interferences have the same form in $`z`$ as the $`qqg`$ splitting function in the single scattering case. The residues of the contour integrations in different processes are different and the resultant different phase factors are essential to give rise to the interesting physical consequences as we discuss below.

### B LPM Interference

Before we list the results of all processes, it is helpful to discuss first the underlying physics in gluon radiation induced by quark-gluon double scattering. As we have seen in the calculation of the contributions from diagrams in Figs. 4-6, the so-called hard-soft processes correspond to the case where the gluon radiation is induced by the hard scattering between the virtual photon and an initial quark with momentum fraction $`x_B`$. The quark is knocked off-shell by the virtual photon and becomes on-shell again after radiating a gluon. Afterwards the quark or the radiated gluon will have a secondary scattering with another soft gluon from the nucleus. We denote such a hard-soft process by the diagram in Fig. 7, where the off-shell quark is marked by a filled circle. In the double hard processes, on the other hand, the quark is on-shell after the first hard scattering with the virtual photon. The gluon radiation is then induced by the scattering of the quark with another gluon that carries finite momentum fraction $`x_L+x_D`$. We denote such hard rescatterings by the diagram in Fig. 8 where an off-shell parton is again marked by a filled circle.

The conventional categorization of soft and hard rescattering by LQS is according to the hardness of the virtual partons or the momentum fraction carried by the second initial parton. The hardness of the off-shell parton is characterized by the virtuality that is proportional to $`\mathrm{}_T^2`$. In the limit $`\mathrm{}_T0`$, these two processes become identical. That is why they cancel each other completely as we see from the final results \[Eqs. (62), (70), (74), (79) and (83)\]. Only gluon radiation with finite transverse momentum has non-vanishing contribution to the radiation matrix elements. In the collinear limit $`\mathrm{}_T0`$, the original terms ”soft” and ”hard” used by LQS in the large $`\mathrm{}_T`$ processes becomes ambiguous. But we will keep using them in this paper for historical reasons and for consistency with the large $`\mathrm{}_T`$ limit. Physically, one should interpret the soft rescattering as the processes in which the gluon radiation is induced by the original first hard scattering ($`\gamma `$-quark), while the hard rescattering as the processes in which the gluon radiation is induced by the rescattering, as illustrated in Figs. 7 and 8. The classification of these two different processes according to the pole structure of radiation matrix elements is most convenient in discussing the interference effects.

We define the formation time of the gluon radiation as

$$\tau _f\frac{2q^{}z(1z)}{\mathrm{}_T^2}.$$
(84)
For very collinear ($`\mathrm{}_T0`$) gluon radiation the formation time can become much larger than the relative distance between the two scattering. Then the two radiation processes should have destructive interference, leading to the LPM interference effect. One can see that our results \[Eqs. (62), (70), (74), (79) and (83)\] have a clear manifest of such an LPM interference effect due to the cancellation by the interferences between hard-soft and double hard scattering. Since the gluon’s formation time should be smaller than the nuclear size $`MR_A/p^+`$ in the given frame, the gluon should then have a minimum transverse momentum $`\mathrm{}_T^2Q^2/MR_A`$, as we will demonstrate later in this paper.

We now reorganize the contributions of different processes according to the above classification of hard-soft and double hard rescattering. We will list our complete calculation of quark-gluon double scattering in the following. To help understand the reorganization of different processes we provide in the Appendix all the radiation amplitudes induced by single, double and triple scattering. There we use the approach of helicity amplitudes for high-energy parton scattering and the results are the same as our complete calculation in the limit of soft radiation ($`z_g=1z0`$).

All the central-cut diagrams are shown in Fig. 9. The contributions from double-hard (Fig. 9a) and hard-soft scattering (Fig. 9b) are,

$`\overline{H}_{CH}^D`$ $`=`$ $`{\displaystyle 𝑑\mathrm{}_T^2\frac{\alpha _s}{2\pi }\frac{1+z^2}{1z}e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\frac{2\pi \alpha _s}{N_c}\theta (y_2^{})\theta (y^{}y_1^{})}`$ (85)
$`\times `$ $`C_A{\displaystyle \frac{k_T^2}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}e^{ix_Lp^+(y^{}y_1^{}+y_2^{})},`$ (86)
$`\overline{H}_{CS}^D`$ $`=`$ $`{\displaystyle 𝑑\mathrm{}_T^2\frac{\alpha _s}{2\pi }\frac{1+z^2}{1z}e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\frac{2\pi \alpha _s}{N_c}\theta (y_2^{})\theta (y^{}y_1^{})}`$ (87)
$`\times `$ $`\{C_F{\displaystyle \frac{1}{\mathrm{}_T^2}}+C_A{\displaystyle \frac{1}{(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}e^{ix_Dp^+(y^{}y_1^{}+y_2^{})/(1z)}`$ (88)
$``$ $`{\displaystyle \frac{C_A}{2}}{\displaystyle \frac{\stackrel{}{\mathrm{}_T}(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}[e^{ix_Dp^+y_2^{}/(1z)}+e^{ix_Dp^+(y^{}y_1^{})/(1z)}]\},`$ (89)
respectively. As we have discussed earlier, the gluon radiation in the double-hard scattering is induced by the rescattering of an on-shell quark with another hard gluon. Therefore the contribution has almost the same form as obtained by Gunion and Bertsch (GB) for soft gluon bremsstrahlung induced by a single elastic scattering of on-shell quarks, except for the phase factors and the splitting function which gives the GB result in the soft radiation limit ($`z_g=1z0`$). In the processes we consider here, the gluon radiation can also be induced by the first hard scattering that produces the quark jet. In these so-called hard-soft processes, the final quark or gluon then has a second scattering with another soft gluon from the nucleus. In Eq. (89), the first term comes from quark rescattering, the second term from gluon rescattering and the last two terms come from the interferences between the first two processes.

There are also two interferences between double hard and hard-soft scattering as shown in Figs. 9c and 9d. Their contributions are,

$`\overline{H}_{CI_2}^D`$ $`=`$ $`{\displaystyle 𝑑\mathrm{}_T^2\frac{\alpha _s}{2\pi }\frac{1+z^2}{1z}e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\frac{2\pi \alpha _s}{N_c}\theta (y_2^{})\theta (y^{}y_1^{})}`$ (90)
$`\times `$ $`()\left\{{\displaystyle \frac{C_A}{2}}{\displaystyle \frac{\stackrel{}{k_T}(\stackrel{}{k_T}\stackrel{}{\mathrm{}_T})}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}+C_A{\displaystyle \frac{\stackrel{}{k_T}\stackrel{}{\mathrm{}_T}}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}e^{ix_Dp^+(y^{}y_1^{})/(1z)}\right\}e^{ix_Lp^+y_2^{}}.`$ (91)
$`\overline{H}_{CI_1}^D`$ $`=`$ $`{\displaystyle 𝑑\mathrm{}_T^2\frac{\alpha _s}{2\pi }\frac{1+z^2}{1z}e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\frac{2\pi \alpha _s}{N_c}\theta (y_2^{})\theta (y^{}y_1^{})}`$ (92)
$`\times ()`$ $`\left\{{\displaystyle \frac{C_A}{2}}{\displaystyle \frac{\stackrel{}{k_T}(\stackrel{}{k_T}\stackrel{}{\mathrm{}_T})}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}+C_A{\displaystyle \frac{\stackrel{}{k_T}\stackrel{}{\mathrm{}_T}}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}e^{ix_Dp^+y_2^{}/(1z)}\right\}e^{ix_Lp^+(y^{}y_1^{})},`$ (93)
It is interesting to note that contributions from the double hard processes or gluon radiation induced by secondary hard scattering and the interferences with hard-soft processes all vanish in the collinear limit of the secondary scattering with a gluon ($`k_T0`$). What remains is the radiation spectrum induced by the first hard scattering (photon-quark). As we will show later, combining with interferences between single and triple scattering in the same collinear limit, it gives the eikonal contribution to the next-to-leading order correction of the single scattering \[Eq. (28)\] corresponding to the first term in the collinear expansion in Eq. (41). For finite $`k_T`$ the above radiation spectra exhibit interesting interference patterns which were discussed in detail by Gyulassy, Lévai and Vitev  and Wiedemann  in the framework of the GW static color-screened potential model. In the present paper, we want to study their contribution to the higher-twist corrections to the fragmentation functions.

To complete our calculation we also have to consider all the interferences between single and triple scattering. These correspond to the asymmetrical-cut diagrams. The right-cut and left-cut dialgrams in Fig. 4 are just two examples. After similar reorganization of all different processes, we list in Figs. 10 all the possible right-cut interference diagrams. The contributions from the first three right-cut diagrams in Fig. 10, which involve double hard scattering, are

$`\overline{H}_{RH}^D`$ $`=`$ $`{\displaystyle 𝑑\mathrm{}_T^2\frac{\alpha _s}{2\pi }\frac{1+z^2}{1z}e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\frac{2\pi \alpha _s}{N_c}\theta (y_2^{})\theta (y_2^{}y_1^{})}`$ (94)
$`\times `$ $`(){\displaystyle \frac{\stackrel{}{k_T}(\stackrel{}{k_T}\stackrel{}{\mathrm{}_T})}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}[{\displaystyle \frac{C_A}{2}}e^{i(x_D^0x_Dx_L)p^+(y_1^{}y_2^{})}{\displaystyle \frac{C_A}{2}}`$ (95)
$``$ $`C_Ae^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})}]e^{ix_Lp^+y_2^{}},`$ (96)
where $`x_D^0=k_T^2/2p^+q^{}`$. We have made the variable change $`k_Tk_T`$ in some of the contributions in order to obtain a more compact form of the final result. The contributions from the two hard-soft processes in the right-cut diagrams in Fig. 10 are

$`\overline{H}_{RS}^D`$ $`=`$ $`{\displaystyle 𝑑\mathrm{}_T^2\frac{\alpha _s}{2\pi }\frac{1+z^2}{1z}e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\frac{2\pi \alpha _s}{N_c}\theta (y_2^{})\theta (y_2^{}y_1^{})}`$ (97)
$`\times `$ $`()\{{\displaystyle \frac{1}{\mathrm{}_T^2}}[C_F+C_Ae^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})}]`$ (98)
$``$ $`{\displaystyle \frac{\stackrel{}{\mathrm{}_T}(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}{\displaystyle \frac{C_A}{2}}e^{ix_Dp^+y_2^{}/(1z)}[1+e^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})}]\}.`$ (99)
Similarly, contributions from the left-cut diagrams, which are just the complex conjugates of the right-cut diagrams in Fig.10, are,

$`\overline{H}_{LH}^D`$ $`=`$ $`{\displaystyle 𝑑\mathrm{}_T^2\frac{\alpha _s}{2\pi }\frac{1+z^2}{1z}e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\frac{2\pi \alpha _s}{N_c}\theta (y^{}y_1^{})\theta (y_1^{}y_2^{})}`$ (100)
$`\times `$ $`(){\displaystyle \frac{\stackrel{}{k_T}(\stackrel{}{k_T}\stackrel{}{\mathrm{}_T})}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}[{\displaystyle \frac{C_A}{2}}e^{i(x_D^0x_Dx_L)p^+(y_1^{}y_2^{})}{\displaystyle \frac{C_A}{2}}`$ (101)
$``$ $`C_Ae^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})}]e^{ix_Lp^+(y^{}y_1^{})},`$ (102)
$`\overline{H}_{LS}^D`$ $`=`$ $`{\displaystyle 𝑑\mathrm{}_T^2\frac{\alpha _s}{2\pi }\frac{1+z^2}{1z}e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\frac{2\pi \alpha _s}{N_c}\theta (y^{}y_1^{})\theta (y_1^{}y_2^{})}`$ (103)
$`\times `$ $`()\{{\displaystyle \frac{1}{\mathrm{}_T^2}}[C_F+C_Ae^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})}]`$ (104)
$``$ $`{\displaystyle \frac{\stackrel{}{\mathrm{}_T}(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})}{\mathrm{}_T^2(\stackrel{}{\mathrm{}_T}\stackrel{}{k_T})^2}}{\displaystyle \frac{C_A}{2}}e^{ix_Dp^+(y^{}y_1^{})/(1z)}[1+e^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})}]\}.`$ (105)
### C Collinear Expansion

Following the procedure of extracting the next-leading-twist contributions to the semi-inclusive DIS cross section, as outlined in Sec. II B, we now expand the hard part in $`k_T`$ according to Eq. (41). To simplify our notation, we factor out all the $`\theta `$-functions and $`k_T`$-independent part of $`\overline{H}_{C,R,L}^D`$(the sum of all the contributions from central-cut, right-cut or left-cut diagrams) and define

$`\overline{H}^D`$ $`=`$ $`{\displaystyle 𝑑\mathrm{}_T^2\frac{\alpha _s}{2\pi }\frac{1+z^2}{1z}e^{i(x+x_L)p^+y^{}+ix_Dp^+(y_1^{}y_2^{})}\frac{2\pi \alpha _s}{N_c}}`$ (106)
$`\times `$ $`\left[H_C^D\theta (y_2^{})\theta (y^{}y_1^{})H_R^D\theta (y_2^{})\theta (y_2^{}y_1^{})H_L^D\theta (y^{}y_1^{})\theta (y_1^{}y_2^{})\right].`$ (107)
From Eqs. (86)-(105), we have

$$\overline{H}_C^D(k_T=0)=\overline{H}_R^D(k_T=0)=\overline{H}_L^D(k_T=0)=\frac{C_F}{\mathrm{}_T^2}.$$
(108)
They all come from hard-soft processes. As we have pointed out before, double-hard scattering processes and their interferences all vanish when $`k_T=0`$. Note that the combination of $`\theta `$-functions,

$`{\displaystyle 𝑑y_1^{}𝑑y_2^{}\left[\theta (y_2^{})\theta (y_2^{}y_1^{})+\theta (y^{}y_1^{})\theta (y_1^{}y_2^{})\theta (y_2^{})\theta (y^{}y_1^{})\right]}`$ (109)
$`=`$ $`{\displaystyle _0^y^{}}𝑑y_1^{}{\displaystyle _0^{y_1^{}}}𝑑y_2^{},`$ (110)
is an ordered integral limited by the value of $`y^{}`$. We will refer to any terms that are proportional to the above combination of $`\theta `$-functions as contact contributions. Combining with Eqs. (37) and (63), the first term in the $`k_T`$ expansion in Eq. (107), $`\overline{H}^D(k_T=0)`$, gives the eikonal contribution to the next-leading-order correction of single scattering in Eq. (28). Such eikonal contributions involving non-physical gauge fields do not correspond to any physical double scattering. They only make the final results gauge invariant and can be gauged away. However, we do see the importance of including all the interferences (right-cut and left-cut diagrams) between single and triple scattering. This is essentially a special case of the generalized proof of factorization of leading-twist contributions in DIS  where contributions from any number of soft gluon rescatterings can be eikonalized.

The dominant contributions to the double quark-gluon scattering come from the quadratic term in the $`k_T`$ expansion of $`\overline{H}^D`$ in Eq. (43). We note that these high-twist terms generally involve two additional spatial integrations with respect to $`y_1^{}`$ and $`y_2^{}`$. If no restriction is imposed, these integrations will only be limited by the nuclear size and thus will produce nuclear enhancement as compared to DIS with a nucleon target. Because of the rapid oscillation of $`e^{ixp^+y^{}}`$, any term proportional to Eq. (110) that limits the integration over $`y_1^{}`$ and $`y_2^{}`$ to the value of $`y^{}`$ will not have nuclear enhancement. Neglecting any term that is proportional to Eq. (110) (contact term) and keeping the leading term when $`\mathrm{}_T0`$, we have

$`_{k_T}^2H_C^D|_{k_T=0}`$ $`=`$ $`{\displaystyle \frac{4C_A}{\mathrm{}_T^4}}(1e^{ix_Lp^+y_2^{}})(1e^{ix_Lp^+(y^{}y_1^{})})+𝒪(x_B/Q^2\mathrm{}_T^2),`$ (111)
$`_{k_T}^2H_L^D|_{k_T=0}`$ $`=`$ $`0+𝒪(x_B/Q^2\mathrm{}_T^2),`$ (112)
$`_{k_T}^2H_R^D|_{k_T=0}`$ $`=`$ $`0+𝒪(x_B/Q^2\mathrm{}_T^2).`$ (113)
The four terms in $`H_C^D`$ correspond to hard-soft, double hard scattering and their interferences. Substituting the above in Eq. (63) and (43), we have the leading contribution to the semi-inclusive tensor from double quark-gluon scattering with quark fragmentation

$`{\displaystyle \frac{W_{\mu \nu }^{D,q}}{dz_h}}`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle 𝑑xH_{\mu \nu }^{(0)}(xp,q)_{z_h}^1\frac{dz}{z}D_{qh}(z_h/z)\frac{\alpha _s}{2\pi }C_A\frac{1+z^2}{1z}}`$ (114)
$`\times `$ $`{\displaystyle \frac{d\mathrm{}_T^2}{\mathrm{}_T^4}\frac{2\pi \alpha _s}{N_c}T_{qg}^A(x,x_L)}+(\mathrm{virtual}\mathrm{correction}),`$ (115)
where

$`T_{qg}^A(x,x_L)`$ $`=`$ $`{\displaystyle \frac{dy^{}}{2\pi }𝑑y_1^{}𝑑y_2^{}e^{i(x+x_L)p^+y^{}}(1e^{ix_Lp^+y_2^{}})(1e^{ix_Lp^+(y^{}y_1^{})})}`$ (116)
$`{\displaystyle \frac{1}{2}}`$ $`A|\overline{\psi }_q(0)\gamma ^+F_\sigma ^+(y_2^{})F^{+\sigma }(y_1^{})\psi _q(y^{})|A\theta (y_2^{})\theta (y^{}y_1^{})`$ (117)
is the quark-gluon correlation function which essentially contains four independent twist-four parton matrix elements in a nucleus. Since $`x_L`$ as defined in Eq. (51) depends on $`\mathrm{}_T`$ and $`z`$, the matrix element $`T_{qg}^A(x,x_L)`$ has an implicit dependence on $`z`$ and $`\mathrm{}_T`$. The contribution from gluon fragmentation is similarly

$`{\displaystyle \frac{W_{\mu \nu }^{D,g}}{dz_h}}`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle 𝑑xH_{\mu \nu }^{(0)}(x,p,q)_{z_h}^1\frac{dz}{z}D_{gh}(z_h/z)\frac{\alpha _s}{2\pi }C_A\frac{1+(1z)^2}{z}}`$ (118)
$`\times `$ $`{\displaystyle \frac{d\mathrm{}_T^2}{\mathrm{}_T^4}\frac{2\pi \alpha _s}{N_c}T_{qg}^A(x,x_L)}+(\mathrm{virtual}\mathrm{correction}).`$ (119)
The above contribution from double quark-gluon scattering is very similar to the next-to-leading order single scattering in Eq. (28). Even the splitting function has the same form except for the color factor. The contribution is, however, proportional to the twist-four parton matrix elements. Since there is not much restriction on the spatial integral, such twist-four matrix elements will give the nuclear enhancement we are looking for. We will return later for more discussions on the $`A`$ dependence of the double scattering contribution.

In the process of a collinear expansion, we have kept $`\mathrm{}_T`$ finite and took the limit $`k_T0`$. As a consequence, the gluon field in one of the twist-four parton matrix elements in Eq. (117) carries zero momentum, corresponding to the soft-hard process. Later in this paper we will discuss how to estimate these twist-four matrix elements. Assuming a factorized form, they will be proportional to the product of quark and gluon distributions in a nucleon. However, the gluon distribution $`xf_g(x)`$ at $`x=0`$ is not defined in QCD. This is because we dropped higher order terms of the collinear expansion in Eq. (113). These terms are considered higher-twist contributions in the LQS framework. As a remedy to the problem, we simply resum a subset of the higher-twist terms in the collinear expansion and restore the phase factors such as $`\mathrm{exp}(ix_Tp^+y^{})`$ which is related to the intrinsic transverse momentum of the initial partons, $`x_Tk_T^2/2p^+q^{}z`$. In this case, soft gluon fields in the parton matrix elements will carry a fractional momentum $`x_T`$. Similarly, one can also restore the collinear structure of the scattering amplitude by replacing $`1/\mathrm{}_T^4`$ in Eq. (113) with $`1/\mathrm{}_T^2(\mathrm{}_T^2+k_T^2)`$.

### D Virtual Corrections

So far we have not considered virtual corrections which will ensure the final result to be infrared safe. The calculation of the virtual correction to the single scattering in Fig. 11 is very similar to the real correction in Fig. 3. One can easily find,

$$\frac{dW_{\mu \nu }^{S(v)}}{dz_h}=\underset{q}{}𝑑xf_q^A(x)H_{\mu \nu }^{(0)}(x,p,q)D_{qh}(z_h)_0^{\mu ^2}\frac{d\mathrm{}_T^2}{\mathrm{}_T^2}\frac{\alpha _s}{2\pi }_0^1𝑑zC_F\frac{1+z^2}{1z}.$$
(120)
One can see that the integral over $`z`$ is infrared divergent at $`z=1`$ or $`z_g=0`$. However, such divergency cancels exactly the infrared divergency in Eq. (28) of the radiative contribution to the single scattering process. The sum is then infrared safe. Using the definition of the $`+`$functions ,

$$_0^1𝑑z\frac{F(z)}{(1z)_+}_0^1𝑑z\frac{F(z)F(1)}{1z}$$
(121)
with $`F(z)`$ being any function which is sufficiently smooth at $`z=1`$, one can rewrite  the sum of radiative and virtual correction \[Eq. (35)\] in the following form

$`D_{qh}(z_h,\mu ^2)`$ $`=`$ $`D_{qh}(z_h)+{\displaystyle _0^{\mu ^2}}{\displaystyle \frac{d\mathrm{}_T^2}{\mathrm{}_T^2}}{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle _{z_h}^1}{\displaystyle \frac{dz}{z}}[\gamma _{qqg}(z)D_{qh}(z_h/z)`$ (122)
$`+`$ $`\gamma _{qgq}(z)D_{gh}(z_h/z)].`$ (123)
The splitting functions are defined as

$`\gamma _{qqg}(z)`$ $`=`$ $`C_F\left[{\displaystyle \frac{1+z^2}{(1z)_+}}+{\displaystyle \frac{3}{2}}\delta (1z)\right],`$ (124)
$`\gamma _{qgq}(z)`$ $`=`$ $`\gamma _{qqg}(1z).`$ (125)
Therefore, with the definition of the $`+`$function, the $`\delta `$-function terms in the splitting functions take into account the self-energy virtual correction which cancels the infrared divergences from the radiative processes. The final renormalized quark fragmentation $`D_{qh}(z_h,\mu ^2)`$ satisfies the DGLAP evolution equations in Eq. (1).

When cast into the DGLAP evolution equation in Eq. (1), the splitting function from the real correction can be interpreted as the probability for the quark to radiate a gluon with momentum fraction $`1z`$. Then one must also take into account the probability of no gluon radiation in the evolution to ensure unitarity. Such unitarity requirement gives rise to the same virtual correction as calculated from the diagram in Fig. 11. In the double scattering case, we will use the same unitarity requirement to obtain virtual corrections. The virtual contribution to the quark fragmentation in double scattering processes is, for example,

$$\frac{W_{\mu \nu }^{D(v),q}}{dz_h}=\underset{q}{}𝑑xH_{\mu \nu }^{(0)}(xp,q)D_{qh}(z_h)\frac{\alpha _s}{2\pi }C_A_0^1𝑑z\frac{1+z^2}{1z}\frac{d\mathrm{}_T^2}{\mathrm{}_T^4}\frac{2\pi \alpha _s}{N_c}T_{qg}^A(x,x_L).$$
(126)
One can single out the infrared divergent part by rewriting the integral

$`{\displaystyle _0^1}𝑑z{\displaystyle \frac{1+z^2}{1z}}T_{qg}^A(x,x_L)`$ $`=`$ $`T_{qg}^A(x,x_L)|_{z=1}{\displaystyle _0^1}𝑑z{\displaystyle \frac{2}{1z}}\mathrm{\Delta }T_{qg}^A(x,\mathrm{}_T^2),`$ (127)
$`\mathrm{\Delta }T_{qg}^A(x,\mathrm{}_T^2)`$ $``$ $`{\displaystyle _0^1}𝑑z{\displaystyle \frac{1}{1z}}\left[2T_{qg}^A(x,x_L)|_{z=1}(1+z^2)T_{qg}^A(x,x_L)\right].`$ (128)
The second term is finite since $`T_{qg}^A(x,x_L)`$ is a smooth function of $`z`$. The first term can be combined with the radiative contribution in Eq. (115) to cancel the infrared divergency. With the help of the $`+`$function, the final result can be expressed as

$`{\displaystyle \frac{W_{\mu \nu }^{D,q}}{dz_h}}`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle 𝑑xH_{\mu \nu }^{(0)}(xp,q)\frac{2\pi \alpha _s}{N_c}\frac{d\mathrm{}_T^2}{\mathrm{}_T^4}_{z_h}^1\frac{dz}{z}D_{qh}(z_h/z)}`$ (129)
$`\times `$ $`{\displaystyle \frac{\alpha _s}{2\pi }}C_A\left[{\displaystyle \frac{1+z^2}{(1z)_+}}T_{qg}^A(x,x_L)+\delta (z1)\mathrm{\Delta }T_{qg}^A(x,\mathrm{}_T^2)\right].`$ (130)
Here the implicit $`z`$-dependence of $`T_{qg}^A(x,x_L)`$ plays an important role in the final result. The above integrand will be proportional to the splitting function for single scattering in Eq. (124) if one ignores the $`z`$ dependence of $`T_{qg}^A(x,x_L)`$. Similarly, the final result for contributions from gluon fragmentation is

$`{\displaystyle \frac{W_{\mu \nu }^{D,g}}{dz_h}}`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle 𝑑xH_{\mu \nu }^{(0)}(xp,q)\frac{2\pi \alpha _s}{N_c}\frac{d\mathrm{}_T^2}{\mathrm{}_T^4}_{z_h}^1\frac{dz}{z}D_{gh}(z_h/z)}`$ (131)
$`\times `$ $`{\displaystyle \frac{\alpha _s}{2\pi }}C_A\left[{\displaystyle \frac{1+(1z)^2}{z_+}}T_{qg}^A(x,x_L)+\delta (z)\mathrm{\Delta }T_{qg}^A(x,\mathrm{}_T^2)\right],`$ (132)
where we have used the fact that $`x_L`$ in Eq. (51) is invariant under the transform $`z1z`$ and so is $`T_{qg}^A(x,x_L)`$.

## IV Quark-quark Double Scattering

We have so far only considered quark-gluon double scattering in a nucleus. After the first photon-quark hard scattering, the leading quark can also rescatter with another quark from the nucleus. Such quark-quark double scattering processes also contribute to the semi-inclusive DIS at twist-four. Unlike quark-gluon rescattering, quark-quark double scattering even contributes at the lowest order without gluon radiation. At the lowest order there is only one kind of quark-quark double scattering diagram, as shown in Fig. 12 and its crossing variations. One can easily calculate their contributions and obtain

$`{\displaystyle \frac{dW_{qq\mu \nu }^{D(0)}}{dz_h}}`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle 𝑑x\frac{dy^{}}{2\pi }𝑑y_1^{}𝑑y_2^{}e^{ixp^+y^{}}A|\overline{\psi }_q(0)\frac{\gamma ^+}{2}\psi _q(y^{})\overline{\psi }_q(y_1^{})\frac{\gamma ^+}{2}\psi _q(y_2^{})|A}`$ (133)
$`\times `$ $`{\displaystyle \frac{2\pi \alpha _s}{N_c}}8C_F{\displaystyle \frac{x_B}{Q^2}}H_{\mu \nu }^{(0)}(x,p,q)\{D_{gh}(z_h)\theta (y_2^{})\theta (y^{}y_1^{})`$ (134)
$``$ $`D_{qh}(z_h)[\theta (y^{}y_1^{})\theta (y_1^{}y_2^{})+\theta (y_2^{})\theta (y_2^{}y_1^{})]\},`$ (135)
where to combine all the crossing diagram we use the fact that fields on the light-cone commute with each other . In the above contributions, the term proportional to the gluon fragmentation function is from the central-cut diagrams while the terms containing quark fragmentation function are from left and right-cut diagrams (or interferences). Using Eq. (110) and neglecting the contact term, we have

$$\frac{dW_{qq\mu \nu }^{D(0)}}{dz_h}=\underset{q}{}𝑑xT_{qq}^{A(I)}(x)\frac{2\pi \alpha _s}{N_c}8C_F\frac{x_B}{Q^2}H_{\mu \nu }^{(0)}(x,p,q)\left[D_{gh}(z_h)D_{qh}(z_h)\right],$$
(136)
where

$$T_{qq}^{A(I)}(x)=\frac{dy^{}}{2\pi }𝑑y_1^{}𝑑y_2^{}e^{ixp^+y^{}}A|\overline{\psi }_q(0)\frac{\gamma ^+}{2}\psi _q(y^{})\overline{\psi }_q(y_1^{})\frac{\gamma ^+}{2}\psi _q(y_2^{})|A\theta (y_2^{})\theta (y^{}y_1^{})$$
(137)
is a four-quark matrix element in a nucleus. This twist-four contribution is explicitly suppressed by $`1/Q^2`$ as compared to the single scattering case. Such leading-order and high-twist contributions from quark-quark double scattering essentially mix the quark and gluon fragmentation functions. In this case, there is no induced radiation and thus no energy loss. However, these contributions will change the final differential cross section or semi-inclusive spectrum, since quark and gluon fragmentation functions are different. One thus can consider the case as modification of fragmentation functions without energy loss. Furthermore, these processes do not change the integrated total cross section. One can verify this by using the momentum sum rule $`_h𝑑z_hz_hD_{q,gh}(z_h)=1`$. This is consistent with a general theorem that final state interaction will not change the total cross section.

We also note that contributions from quark-quark double scattering are free of any divergences, especially collinear divergency. Therefore, they will not contribute to the QCD evolution of the effective parton fragmentation functions in the nuclear medium.

To the next order we have to consider radiative corrections to the processes in Fig. 12. In addition we also have to consider the diagram in Fig. 13. Again, calculations of these diagrams are tedious but straightforward. Take the diagram in Fig. 13 for example: The pole structure is exactly the same as the central-cut diagram of quark-gluon double scattering in Fig. 4. The resultant contribution is

$`{\displaystyle \frac{dW_{Fig.\text{13}\mu \nu }^{D(1)}}{dz_h}}`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle _{z_h}^1}{\displaystyle \frac{dz}{z}}D_{qh}(z_h/z){\displaystyle \frac{d\mathrm{}_T^2}{\mathrm{}_T^2}\frac{\alpha _s}{2\pi }C_F\frac{1+z^2}{(1z)^2}\frac{x_B}{Q^2}}`$ (138)
$`\times `$ $`{\displaystyle \frac{2\pi \alpha _s}{N_C}}{\displaystyle 𝑑xT_{qq}^{A(II)}(x,x_L)H_{\mu \nu }^{(0)}(x,p,q)},`$ (139)
where

$`T_{qq}^{A(II)}(x,x_L)`$ $`=`$ $`{\displaystyle \underset{q_i}{}}{\displaystyle \frac{dy^{}}{2\pi }𝑑y_1^{}𝑑y_2^{}e^{i(x+x_L)p^+y^{}}(1e^{ix_Lp^+y_2^{}})(1e^{ix_Lp^+(y^{}y_1^{})})}`$ (140)
$`\times `$ $`A|\overline{\psi }_q(0){\displaystyle \frac{\gamma ^+}{2}}\psi _q(y^{})\overline{\psi }_{q_i}(y_1^{}){\displaystyle \frac{\gamma ^+}{2}}\psi _{q_i}(y_2^{})|A\theta (y_2^{})\theta (y^{}y_1^{}).`$ (141)
The four terms with different phase factors in the above equation correspond to hard-soft, double hard and the interferences. Notice that the two diagonal terms are defined exactly the same as $`T_{qq}^{A(I)}(x)`$. This structure is exactly the same as in the quark-gluon double scattering. However, this contribution is only proportional to $`1/\mathrm{}_T^2`$ versus $`1/\mathrm{}_T^4`$ in the quark-gluon double scattering. It is thus suppressed by $`1/Q^2`$. This remains to be the case for all quark-quark double scattering processes. Since we are interested only in the collinear behavior of the gluon radiation processes in order to study the QCD evolution, we will neglect in this paper all radiative contributions from quark-quark double scattering. Similarly, contributions proportional to $`1/\mathrm{}_T^2`$ from gluon-gluon double scattering can also be neglected to this approximation, e.g., in the collinear expansion in Eq. (113).

## V Modified Quark Fragmentation Function

For a complete result of the semi-inclusive cross section of DIS off a nucleus, one should also include the higher-twist contribution to the quark distributions and their QCD evolution as studied by Mueller and Qiu . Here we will simply replace the leading-twist quark distributions $`f_q^A(x)`$ by $`\stackrel{~}{f}_q^A(x,\mu _I^2)`$ which contains nuclear modification to the quark distributions and their QCD evolution and $`\mu _I^2`$ is the factorization scale for the quark distributions in a nucleus. We should note that even the leading-twist quark distributions $`f_q^A(x)`$ in a nucleus are different from $`A`$ free nucleons \[$`f_q^N(x)`$\].

### A Modified Evolution Equations

Including the higher-twist contributions to the quark distributions and summing up all the leading contributions from single and double scattering processes in Eqs. (123), (130) and (132), we have

$$\frac{dW_{\mu \nu }}{dz_h}=\underset{q}{}𝑑x\stackrel{~}{f}_q^A(x,\mu _I^2)H_{\mu \nu }^{(0)}(x,p,q)\stackrel{~}{D}_{qh}(z_h,\mu ^2)$$
(142)
as the total semi-inclusive tensor in DIS off a nucleus up to twist-four corrections. We define the modified effective quark fragmentation function as

$`\stackrel{~}{D}_{qh}(z_h,\mu ^2)`$ $``$ $`D_{qh}(z_h,\mu ^2)+{\displaystyle _0^{\mu ^2}}{\displaystyle \frac{d\mathrm{}_T^2}{\mathrm{}_T^2}}{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle _{z_h}^1}{\displaystyle \frac{dz}{z}}[\mathrm{\Delta }\gamma _{qqg}(z,x,x_L,\mathrm{}_T^2)D_{qh}(z_h/z)`$ (143)
$`+`$ $`\mathrm{\Delta }\gamma _{qgq}(z,x,x_L,\mathrm{}_T^2)D_{gh}(z_h/z)],`$ (144)
where $`D_{qh}(z_h,\mu ^2)`$ is given in Eq. (123) for leading-twist contributions and

$`\mathrm{\Delta }\gamma _{qqg}(z,x,x_L,\mathrm{}_T^2)`$ $`=`$ $`\left[{\displaystyle \frac{1+z^2}{(1z)_+}}T_{qg}^A(x,x_L)+\delta (1z)\mathrm{\Delta }T_{qg}^A(x,\mathrm{}_T^2)\right]{\displaystyle \frac{C_A2\pi \alpha _s}{(\mathrm{}_T^2+k_T^2)N_c\stackrel{~}{f}_q^A(x,\mu _I^2)}}`$ (145)
$`\mathrm{\Delta }\gamma _{qgq}(z,x,x_L,\mathrm{}_T^2)`$ $`=`$ $`\mathrm{\Delta }\gamma _{qqg}(1z,x,x_L,\mathrm{}_T^2).`$ (146)
The twist-four matrix element $`T_{qg}^A(x,x_L)`$ is given in Eq. (117) and $`\mathrm{\Delta }T_{qg}^A(x,\mathrm{}_T^2)`$ is given in Eq. (128). Since the two-parton correlation function $`T_{qg}^A`$ can involve partons from two different nucleons inside the nucleus, it should have a nuclear enhancement that is proportional to the nuclear size $`R_A`$, or $`T_{qg}^A/\stackrel{~}{f}_q^AA^{1/3}`$, as we will show in the next subsection. Because of the LPM interference structure contained in $`T_{qg}^A(x,x_L)=0`$ ($`x_L\mathrm{}_T^20`$), only those radiated gluons induced by multiple parton scattering that have finite transverse momentum, $`\mathrm{}_T^2Q^2/A^{1/3}`$ contribute to the fragmentation function. Therefore, the nuclear correction to the fragmentation function in Eq. 144 is proportional to $`\alpha _sA^{2/3}/Q^2`$ from double parton scattering. For large $`A`$ and $`Q^2`$, this is the leading nuclear correction. Since $`T_{qg}(x,x_L)/\stackrel{~}{f}_q^A(x)`$ is proprotional to the gluon density per unit transverse area in the nucleus $`\rho _g`$ and $`\alpha _s/Q^2`$ can be considered as the parton rescattering cross section $`\sigma _g`$, the leading nuclear correction to the fragmentation function here is then proportional to $`\sigma _g\rho _gR_A`$. This is exactly the same as the expansion parameter in the opacity expansion approach in Refs. .

We should emphasize here that the factorized form of the semi-inclusive tensor in Eq.(142) only serves to define the effective quark fragmentation function $`\stackrel{~}{D}_{qh}(z_h,\mu ^2)`$. Such a modified fragmentation function for final state particle production has an explicit dependence on the initial parton distribution through the high-twist double scattering processes. Therefore, the factorization for semi-inclusive processes in DIS is broken explicitly at twist-four correction. This is a natural consequence of the non-vanishing parton energy loss at twist-four when the leading quark suffers multiple scattering through the nuclear medium.

Taking the derivative with respect to the collinear factorization scale $`\mu ^2`$, we obtain the modified DGLAP evolution equation in leading order of $`\alpha _s`$ for the modified quark fragmentation function,

$`{\displaystyle \frac{\stackrel{~}{D}_{qh}(z_h,\mu ^2)}{\mathrm{ln}\mu ^2}}`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle _{z_h}^1}{\displaystyle \frac{dz}{z}}[\stackrel{~}{\gamma }_{qqg}(z,x,x_L,\mu ^2)\stackrel{~}{D}_{qh}(z_h/z,\mu ^2)`$ (147)
$`+`$ $`\stackrel{~}{\gamma }_{qgq}(z,x,x_L,\mu ^2)D_{gh}(z_h/z,\mu ^2)].`$ (148)
The modified splitting functions are defined as

$`\stackrel{~}{\gamma }_{qqg}(z,x,x_L,\mu ^2)`$ $`=`$ $`\gamma _{qqg}(z)+\mathrm{\Delta }\gamma _{qqg}(z,x,x_L,\mu ^2)`$ (149)
$`\stackrel{~}{\gamma }_{ggq}(z,x,x_L,\mu ^2)`$ $`=`$ $`\stackrel{~}{\gamma }_{qqg}(1z,x,x_L,\mu ^2),`$ (150)
where $`\gamma _{qqg}(z)`$ given in Eq. (124) is the splitting function in single scattering processes. If one only considers single-jet events in DIS, the gluon fragmentation function does not contribute in the leading order. Before we carry out similar calculations for the higher-twist correction to the DGLAP evolution of the gluon fragmentation functions, we may assume that the gluon fragmentation function which enters into the above equation follows the normal DGLAP evolution,

$$\frac{D_{gh}(z_h,\mu ^2)}{\mathrm{ln}\mu ^2}=\frac{\alpha _s}{2\pi }_{z_h}^1\frac{dz}{z}\left[\underset{q=1}{\overset{2n_f}{}}\gamma _{gq\overline{q}}(z)\stackrel{~}{D}_{qh}(z_h/z,\mu ^2)+\gamma _{ggg}(z)D_{gh}(z_h/z,\mu ^2)\right],$$
(151)
where the normal splitting functions are

$`\gamma _{gq\overline{q}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[z^2+(1z)^2],`$ (152)
$`\gamma _{ggg}`$ $`=`$ $`2C_A\left[{\displaystyle \frac{z}{(1z)_+}}+{\displaystyle \frac{1z}{z}}+z(1z){\displaystyle \frac{1}{12}}(11{\displaystyle \frac{2}{3}}n_f)\delta (z1)\right],`$ (153)
and $`n_f`$ is the number of quark flavors. All the fragmentation functions obey the momentum sum rule,

$$_0^1𝑑z\underset{h}{}zD_{ah}(z,\mu ^2)=1.$$
(154)
One can check that the modified quark fragmentation function $`\stackrel{~}{D}_{qh}(z,\mu ^2)`$ in Eq. (144) still satisfies the momentum sum rule. Such a momentum sum rule might seem count-intuitive since there is momentum transfer of $`x_Lp^+`$ to the fragmentating quark from the second partons in the nucleus in the double-hard processes. However, one should note that the initial quark from the nucleus carries momentum $`x_Bp^+`$ in this case. In single scattering and hard-soft processes, however, the initial quark carries momentum $`(x_B+x_L)p^+`$. So the total momentum transfer from the nucleus to the quark is the same in both single and double scattering. The momentum sum rule for the modified fragmentation functions should still be valid.

Compared to Eq. (1), the renormalization equation for the modified quark fragmentation function in Eq. (148) is similar to the original DGLAP evolution equation for the fragmentation functions in vacuum. However, the modified splitting functions $`\stackrel{~}{\gamma }`$ \[Eq. (149)\] have an extra term $`\mathrm{\Delta }\gamma `$ from induced gluon radiation. Similarly to the ordinary parton cascade in vacuum, the induced radiation will soften the modified quark fragmentation function. As we have argued in the Introduction, such softening will be the only experimentally measurable consequence of quark energy loss in a medium. Another important difference between the modified and the original DGLAP evolution equations is that the induced splitting function $`\mathrm{\Delta }\gamma `$ depends on the twist-four parton matrix elements $`T_{qg}^A(x,x_L)`$ of the nucleus. Because of the interferences between hard-soft and double-hard processes, $`T_{qg}^A(x,x_L)`$ explicitly manifests the LPM effect which modifies the gluon radiation spectra such that the induced splitting functions are also modified from their form in vacuum. The induced splitting functions also explicitly depend on the factorization scale $`\mu `$.

### B Twist-four Parton Matrix Elements

As we have observed, modified DGLAP evolution equations for parton fragmentation functions depend on both the parton distribution $`\stackrel{~}{f}_q^A(x,\mu _I^2)`$ and two-parton correlation functions $`T_{qg}^A(x,x_L)`$. Naturally, they are essential to any numerical solution to the modified DGLAP evolution equations for the fragmentation functions. One also has to know them, especially their dependence on the nuclear size, in order to have any numerical estimate of the effective quark energy loss.

As defined in Eq. (117), the parton correlation function $`T_{qg}^A(x,x_L)`$ essentially contains four independent twist-four parton matrix elements in a nucleus. They are in principle not calculable and can only be measured independently in experiments just like parton distributions. However, under certain assumptions, one might be able to relate the parton correlation functions to single parton distributions in the nuclei and in the meantime obtain the $`A`$-dependence of the parton correlations. There are four parton field operators and three spatial (in the longitudinal direction) integrations in $`T_{qg}^A(x,x_L)`$. We assume that the nuclear wave function can be expressed as a multiple-nucleon state and each nucleon is in a color singlet state. In this case, there should not be long range color correlation between two nucleons due to color confinement. In other words, the interaction between nucleons inside a nucleus can only be mediated via color-singlet objects. Consider first the case where the quark and gluon fields operate on different nucleons inside the nucleus. This means that the double scattering happens with two partons from two different nucleons. Because of the color confinement, $`y_1`$ and $`y_2`$ in Eq. (117) are restricted such that $`|y_1y_2|r_N`$, where $`r_N`$ is the radius of a nucleon. The integrations over $`y_1`$ and $`y_2`$ should, therefore, give the length scale of $`r_NR_A`$, where $`R_A=1.12A^{1/3}`$ is the radius of the nucleus. In this case the correlation functions in Eq. (117) should be approximately proportional to $`A^{4/3}`$. If both the quark and gluon fields operate on the same nucleon state inside the nucleus, $`y_1`$ and $`y_2`$ are both restricted to be within the size of the nucleon. Then the correlation functions should only be roughly proportional to $`A`$. This corresponds to double scattering with two partons in the same nucleon. As we have discussed before, we will neglect these contributions because they yield no nuclear enhancement.

Among the four independent matrix elements in $`T_{qg}^A(x,x_L)`$, there are two diagonal terms, one proportional to a phase factor $`\mathrm{exp}[i(x+x_L)p^+y^{}+ix_Tp^+(y_1^{}y_2^{})]`$ and another to $`\mathrm{exp}[ixp^+y^{}+i(x_L+x_T)p^+(y_1^{}y_2^{})]`$. The first term corresponds to hard-soft double scattering processes in which the gluon carrying momentum fraction $`x_T`$ is soft. The second term comes from double-hard scattering in which the gluon carries large momentum fraction $`x_L+x_T`$. In this paper, we adopt the approximation used in Ref. . Under such an approximation, one simply relates the two-parton correlations to the product of two single parton distributions,

$`{\displaystyle \frac{dy^{}}{2\pi }𝑑y_1^{}𝑑y_2^{}}`$ $`e^{ix_1p^+y^{}+ix_2p^+(y_1^{}y_2^{})}`$ (155)
$`\times `$ $`{\displaystyle \frac{1}{2}}A|\overline{\psi }_q(0)\gamma ^+F_\sigma ^+(y_2^{})F^{+\sigma }(y_1^{})\psi _q(y^{})|A\theta (y_2^{})\theta (y^{}y_1^{})`$ (156)
$`=`$ $`{\displaystyle \frac{C}{x_A}}f_q^A(x_1)x_2f_g^N(x_2),`$ (157)
where $`x_A=1/MR_A`$, $`f_q^A(x)`$ is the quark distribution inside a nucleus as defined in Eq. (17), and $`f_g^N(x)`$ is the gluon distribution inside a nucleon as defined by

$$f_g^N(x)\frac{1}{xp^+}\frac{dy^{}}{2\pi }e^{ixp^+y^{}}N|F_\sigma ^+(0)F^{+\sigma }(y^{})|N.$$
(158)
$`C`$ is assumed to be a constant, reflecting the strength of two-parton correlation inside a nucleus. As we have argued, $`C`$ should be proportional to the nucleon radius $`r_N`$, characterizing the confinement scale. We also assume that the nuclear distribution takes a Gaussian form $`\rho (r)\mathrm{exp}(r^2/2R_A^2)`$. In terms of the light-cone coordinates, we have then $`\rho (y^{})=\rho _0\mathrm{exp}(y_{}^{}{}_{}{}^{2}/2R_{A}^{}{}_{}{}^{2})`$ where $`R_A^{}=\sqrt{2}R_AM/p^+`$ and $`M`$ is the nucleon mass. The distribution is normalized such that

$$_0^{\mathrm{}}𝑑y^{}\rho _0e^{y_{}^{}{}_{}{}^{2}/2R_{A}^{}{}_{}{}^{2}}=1.$$
(159)
With such a distribution, the third integration in Eq. (157) gives $`𝑑y_2^{}R_A^{}1/x_Ap^+`$.

In addition to the two diagonal terms, there are also two off-diagonal matrix elements in $`T_{qg}^A(x,x_L)`$ which come from the interferences between hard-soft and double-hard scattering processes. These two new matrix elements introduced by the interferences have never been studied before in DIS. Similar matrix elements have been discussed before in the context of $`pA`$ collisions . For finite value of $`x_L`$, this involves off-diagonal parton matrix elements and they should be related to the so-called skewed parton distributions of nucleons . For a rough estimate here, we can generalize the approximation on parton correlations in Ref.  to the off-diagonal matrix elements. As we have pointed out in the calculation of the cut-diagrams in Fig. 4, in the interferences between double-hard and hard-soft processes, there is actually momentum flow of $`x_Lp^+`$ between the two nucleons where the initial quark and gluon come from. Without strong long range two-nucleon correlation inside a nucleus, the amount of momentum flow $`x_Lp^+`$ should then be restricted by uncertainty principle. Note that the phase factors in the off-diagonal matrix elements have an additional factor of $`\mathrm{exp}(\pm ix_Lp^+y_2^{})`$ relative to their diagonal counterpart. Assuming the same Gaussian nuclear distribution in light-cone coordinates and using

$$_0^{\mathrm{}}𝑑y^{}\rho _0e^{y_{}^{}{}_{}{}^{2}/2R_{A}^{}{}_{}{}^{2}\pm ix_Lp^+y^{}}=e^{x_L^2/x_A^2},$$
(160)
we should have a similar approximation for the off-diagonal matrix elements:

$`{\displaystyle \frac{dy^{}}{2\pi }𝑑y_1^{}𝑑y_2^{}}`$ $`e^{ix_1p^+y^{}+ix_2p^+(y_1^{}y_2^{})\pm ix_Lp^+y_2^{}}`$ (161)
$`\times `$ $`{\displaystyle \frac{1}{2}}A|\overline{\psi }_q(0)\gamma ^+F_\sigma ^+(y_2^{})F^{+\sigma }(y_1^{})\psi _q(y^{})|A\theta (y_2^{})\theta (y^{}y_1^{})`$ (162)
$`=`$ $`{\displaystyle \frac{C}{x_A}}f_q^A(x_1)x_2f_g^N(x_2)e^{x_L^2/x_A^2}.`$ (163)
Relative to the diagonal matrix elements, the off-diagonal ones are suppressed by a factor $`\mathrm{exp}(x_L^2/x_A^2)`$. Combining all the four terms together, we have

$$T_{qg}^A(x,x_L)=\frac{C}{x_A}[f_q^A(x+x_L)x_Tf_g^N(x_T)+f_q^A(x)(x_L+x_T)f_g^N(x_L+x_T)](1e^{x_L^2/x_A^2}).$$
(164)
Notice that $`\tau _f=1/x_Lp^+`$ is the gluon’s formation time. Thus, $`x_L/x_A=L_A^{}/\tau _f`$ with $`L_A^{}=R_AM/p^+`$ being the nuclear size in the our chosen frame. It is then clear from the above that the interferences between double-hard and hard-soft scattering cancel exactly the double scattering contributions for collinear gluon radiation whose formation time is much larger than the nuclear size. The LPM effect now explicitly manifests itself via the effective matrix elements in the gluon radiation spectrum that is induced by double parton scattering.

The diagonal parton matrix elements that arise from hard-soft and double hard processes are exactly the same as what appear in processes with large final transverse momentum such as di-jet production studied in Ref. . These diagonal matrix elements have an enhancement linear in nuclear size $`R_A`$ relative to the single parton distribution in nuclei as in Eq. (157). The off-diagonal parton matrix elements, on the other hand, come from the interference between hard-soft and double hard processes. Because such parton matrix elements involve transferring momentum fraction $`x_L`$ between two nucleons inside a nucleus, they should be suppressed when $`x_L`$ is larger than what the uncertainty principle allows. For a Gaussian nuclear distribution, these off-diagonal matrix elements are exponentially suppressed for large $`x_L/x_A`$ as in Eq. (163). This is why one can neglect the interferences between soft and hard rescattering in the LQS  study of processes with large final transverse momentum, $`\mathrm{}_T>>Q^2/A^{1/3}`$. However, in our current study of fragmentation processes, the off-diagonal matrix elements in the interferences are very important. The interference terms cancel exactly the normal contributions that contain the diagonal matrix elements as shown in Eq. (164) for $`\mathrm{}_T=0`$. For small but finite $`\mathrm{}_T`$, the cancellation is not complete. The effective parton matrix element (diagonal minus the off-diagonal ones) is then proportional to $`x_L^2/x_A^3`$ which depends cubically on the nuclear size $`R_A`$. It, however, also vanishes quadratically with $`x_L\mathrm{}_T^2`$. Therefore, the interferences effectively cut off the transverse momentum of the radiated gluon at $`\mathrm{}_T^2Q^2/A^{1/3}`$. In terms of the LPM interference, this means that the formation time of induced gluon radiation cannot be larger than the nuclear size $`R_A`$. Such a consequence of the LPM interference effect is critical for us to apply the LQS generalized factorization to the problem of nuclear modification of fragmentation function in this study.

### C Parton Energy Loss

So far we have derived the modification to the quark fragmentation functions and their DGLAP evolution equations. Given the twist-four parton matrix elements, one can then solve the modified DGLAP evolution equations for the modified fragmentation functions. Doing so, one can effectively resum all the leading-log corrections at twist-four. We leave the numerical solutions to future studies. In this paper, we instead will formulate the effective energy loss of leading quarks.

In principle the modification of the fragmentation functions would be the only experimental effect of induced gluon radiation via multiple scattering. One can never directly measure the energy loss of the leading quark. The net effect of the energy loss is the suppression of leading particles on one hand and the enhancement of soft particles on the other, leading to the modification of the fragmentation functions. One can then experimentally characterize the parton energy loss via the momentum transfer from large to small momentum regions of the fragmentation functions. We will discuss this in detail when we numerically calculate the modified fragmentation functions.

Upon a close examination of Eq. (144), we see that the first term is the renormalized fragmentation function in vacuum. The rest are particle production induced by the rescattering of the quark through the nuclear medium. In particular, the last term is particle production from the fragmentation of the gluon which is induced by the secondary scattering. Such particle production is at the expense of the energy loss of the leading quark. We can thus quantify the quark energy loss by the momentum fraction carried by the induced gluon,

$`\mathrm{\Delta }z_g(x_B,\mu ^2)`$ $`=`$ $`{\displaystyle _0^{\mu ^2}}{\displaystyle \frac{d\mathrm{}_T^2}{\mathrm{}_T^2}}{\displaystyle _0^1}𝑑z{\displaystyle \frac{\alpha _s}{2\pi }}z\mathrm{\Delta }\gamma _{qgq}(z,x_B,x_L,\mathrm{}_T^2)`$ (165)
$`=`$ $`{\displaystyle _0^{\mu ^2}}{\displaystyle \frac{d\mathrm{}_T^2}{\mathrm{}_T^2}}{\displaystyle _0^1}𝑑z[1+(1z)^2]{\displaystyle \frac{T_{qg}^A(x_B,x_L)}{\stackrel{~}{f}_q^A(x_B,\mu _I^2)}}{\displaystyle \frac{C_A\alpha _s^2}{N_c(\mathrm{}_T^2+k_T^2)}}.`$ (166)
Using the approximation of the parton correlation functions in Eq. (164) and changing the integration variable, we have

$`\mathrm{\Delta }z_g(x_B,\mu ^2)`$ $`=`$ $`{\displaystyle \frac{C_A\alpha _s^2}{N_c\stackrel{~}{f}_q^A(x_B,\mu _I^2)}}{\displaystyle \frac{x_B}{x_AQ^2}}C{\displaystyle _0^1}𝑑z{\displaystyle \frac{1+(1z)^2}{z(1z)}}{\displaystyle _0^{x_\mu }}{\displaystyle \frac{dx_L}{x_L(x_L+x_T)}}(1e^{x_L^2/x_A^2})`$ (167)
$`\times `$ $`[f_q^A(x_B+x_L)x_Tf_g^N(x_T)+f_q^A(x_B)(x_L+x_T)f_g^N(x_L+x_T)],`$ (168)
where $`x_\mu =\mu ^2/2p^+q^{}z(1z)=x_B/z(1z)`$ if we choose the factorization scale as $`\mu ^2=Q^2`$.

One can numerically evaluate the above in the future with assumed parton distributions in nuclei. It is nevertheless interesting to discuss the nuclear dependence of the parton energy. The exponential factor in the above equation comes from the LPM interferences. It regularizes the integration over $`x_L`$ which is otherwise divergent. Combining with the QCD radiation spectrum, it also limits the integration over $`x_L`$ to $`x_L<x_A`$ and gives $`𝑑x_L/x_L^21/x_A`$. The fractional energy loss by the quark is, therefore, proportional to

$$\mathrm{\Delta }z_g\frac{C_A\alpha _s^2}{N_c}\frac{x_B}{x_A^2Q^2}.$$
(169)
Since $`x_A=1/MR_A`$, the energy loss thus depends quadratically on the nuclear size. The extra size dependence comes from the combination of the QCD radiation spectrum and the modification of the available phase space in $`\mathrm{}_T`$ or $`x_L`$ due to the LPM interferences. One also notes that, though the fractional energy loss is suppressed by $`1/Q^2`$ for a fixed value of $`x_B`$, the total energy loss, $`\mathrm{\Delta }E=q^{}\mathrm{\Delta }z_g`$, is not.

## VI Summary and Discussions

Working in the framework of generalized factorization of higher-twist parton distributions, we have studied multiple parton scattering in deeply inelastic $`eA`$ collisions, in particular the induced gluon radiation and the resultant modification to the final hadron spectra. We have defined modified quark fragmentation functions to take into account the effect of multiple parton scattering. We have presented a detailed derivation of the modified quark fragmentation functions and their QCD evolution equations to the next-leading-twist. The modification is shown to depend on twist-four parton matrix elements, both diagonal and off-diagonal, in nuclei.

Depending on whether the gluon radiation is induced by the secondary scattering, one can categorize the multiple parton scattering as soft or hard, according to the fractional momentum carried by the secondary parton involved. We have considered both soft and hard scattering and their interferences. We have shown that these are exactly the so-called LPM interference effect. We have demonstrated that LPM interferences modify the available phase space in the emitted gluon’s momentum. Coupled with the gluon spectrum in QCD, this leads to the quadratical dependence of the modification of fragmentation functions or the effective parton energy loss on the nuclear size $`R_A`$.

The LPM interference effect in this study is critical for applying the LQS generalized factorization to the fragmentation processes. Because of the LPM interference effect, the gluon radiation with small transverse momentum is suppressed. Equivalently, the formation time of the gluon radiation has to be smaller than the nuclear size $`R_A`$. This requires that the radiated gluon have a minimum transverse momentum $`\mathrm{}_T^2Q^2/A^{1/3}`$. For such a value of transverse momentum, we can still take the leading logarithmic approximation $`\mathrm{}_T^2Q^2`$ in the study of fragmentation function in a large nucleus ($`A1`$). Therefore, the applicability of the our study here is limited to experimental situations in which $`Q^2`$ must be large enough such that $`Q^2Q^2/A^{1/3}\mathrm{\Lambda }_{\mathrm{QCD}}^2`$. In the meantime, the resultant nuclear corrections proportional to $`\alpha _sA^{1/3}/\mathrm{}_T^2\alpha _sA^{2/3}/Q^2`$ from double parton scattering are still leading for large values of $`Q^2`$. This is why the nuclear corrections to the fragmentation function and the induced parton energy loss depend quadratically on the nuclear size $`R_A`$. In the processes with large final transverse momentum $`\mathrm{}_T^2Q^2`$ as studied by LQS , the leading correction is instead proportional to $`\alpha _sA^{1/3}/Q^2`$. This is because one can neglect the interference contribution if $`\mathrm{}_T^2Q^2/A^{1/3}`$.

We have also considered double-quark scattering processes. Though their contributions to the QCD evolution equations can be neglected as compared to quark-gluon scattering, they do have a leading-order contribution which mixes quark and gluon fragmentation functions. Since they involve quark-antiquark correlations in nuclei, the modification to quark and antiquark fragmentation functions will be different. This might give different modification to the spectra for negative and positive hadrons as observed in experiments

There is currently little information on the twist-four parton matrix elements in nuclei, especially the off-diagonal ones. We have only outlined a very crude estimate, assuming factorization of the two parton correlations in nuclei. This enables us to estimate the nuclear dependence of the final results. Future work on this is necessary in order to have any quantitative study of the problem. With that information, one should be able to numerically solve the QCD evolution equation for the modified fragmentation functions.

The method we developed here to study the modification of parton fragmentation functions cannot yet be applied directly to other processes, such as high-energy $`pA`$ and $`AA`$ collisions. One can, however, find important implications from the study in this paper. The $`A^{4/3}`$ nuclear size dependence of the twist-four parton matrix elements relies on the color confinement within the nucleon radius $`r_N`$. So the strength of parton correlation as represented by the constant $`C`$ in Eq. (157) is proportional to $`r_N`$. If we replace the nuclear target in DIS by a droplet of quark-gluon plasma, then the integration over $`y_1^{}y_2^{}`$ in Eq. (157) is no longer restricted to $`r_N`$ but rather to the effective parton correlation length $`1/\mu _D`$ with $`\mu _D`$ being the Debye screening mass in the quark-gluon plasma. So the modification of the parton fragmentation function and the effective parton energy loss will be sensitive to $`\mu _D`$. If there is a dramatic change in the parton correlation during the QCD phase transition, one should then expect a similar behavior in parton energy loss.

## Acknowledgements

We thank J. Qiu and G. Sterman for helpful discussions. This work was supported by the Director, Office of Energy Research, Office of High Energy and Nuclear Physics, Divisions of Nuclear Physics, of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098 and DE-FG-02-96ER40989 and in part by NSFC under project 19928511.

## Appendix

In this Appendix we calculate the hard part of gluon radiation processes associated with multiple quark-gluon scattering in DIS off a nucleus. We use the technique of helicity amplitude for high-energy parton scattering, where one can neglect the transverse recoil induced by the scattering. The final results will agree with the complete calculation in the limit of soft radiation. Such an exercise can help us to cross-check the results of our complete calculation. More importantly, it can help us to understand the physical processes more clearly on the level of amplitudes and to reorganize the final results in this paper according to the nature of different processes.

### A Single Scattering

Assuming the dominant component of a fast quark is $`\mathrm{}_q[0,\mathrm{}_q^{},\stackrel{}{0}_{}]`$. One has,

$$\overline{u}(\mathrm{}_q,\lambda )\gamma ^\mu u(\mathrm{}_q^{},\lambda ^{})2\sqrt{\mathrm{}_q^{}\mathrm{}_q^{}}\delta _{\lambda \lambda ^{}}\underset{¯}{n}^\mu ,$$
(170)
where $`\underset{¯}{n}=[0,1,\stackrel{}{0}_{}]`$, $`\lambda `$ and $`\lambda ^{}`$ are quark helicities. According to the collinear factorization of the parton matrix elements, each initial quark line contributes $`u(p)`$ while each initial gluon line (with Lorentz index $`\sigma `$) contributes $`p_\sigma `$ to the amplitude. Each initial parton line with momentum $`p_i=[x_ip^+,0,\stackrel{}{k}_T]`$ also has a momentum integral $`dx_i/2\pi `$ with a phase factor $`\mathrm{exp}(ix_ip^+y_i^{}+i\stackrel{}{k}_T\stackrel{}{y}_{Ti})`$. We summarize these special Feynman rules for initial parton lines as:

$`\mathrm{initial}\mathrm{quark}`$ $``$ $`u(p){\displaystyle \frac{dx_i}{2\pi }e^{ix_ip^+y_i^{}+i\stackrel{}{k}_T\stackrel{}{y}_{Ti}}},`$ (171)
$`\mathrm{initial}\mathrm{gluon}`$ $``$ $`p_\sigma {\displaystyle \frac{dx_i}{2\pi }e^{ix_ip^+y_i^{}+i\stackrel{}{k}_T\stackrel{}{y}_{Ti}}}.`$ (172)
All the internal and final external parton lines follow the normal Feynman rules.

We will work with an axial gauge $`A^{}=0`$. The final gluon with momentum $`\mathrm{}`$ is assumed to carry $`1z`$ momentum fraction of the struck quark and has polarization $`\epsilon (\mathrm{})`$,

$`\mathrm{}`$ $`=`$ $`[{\displaystyle \frac{\mathrm{}_T^2}{2(1z)q^{}}},(1z)q^{},\stackrel{}{\mathrm{}_T}],`$ (173)
$`\epsilon (\mathrm{})`$ $`=`$ $`[{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{(1z)q^{}}},0,\stackrel{}{ϵ}_T].`$ (174)
According to the above Feynman rules, the amplitude for gluon radiation in the single scattering case of DIS as shown in Fig. 14(a) is

$`M_\mu ^S(y)`$ $`=`$ $`{\displaystyle \frac{dx}{2\pi }\overline{u}(xp+q)\gamma _\mu u(p)\overline{M}^S},`$ (175)
$`\overline{M}^S(y)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}T_ce^{i(x_B+x_L)p^+y^{}},`$ (176)
where we have used the on-shell condition for the final quark

$$(xp+q\mathrm{})^2=2p^+q^{}(xx_Bx_L)=0,$$
(177)
and assumed soft radiation limit $`1z0`$. $`T_c`$ is the color matrix in the adjunct representation with the color index $`c`$ for the radiated gluon.

To calculate the hard part of the parton scattering, one can simply take the square of the amplitude $`M_\nu ^{S(1)}(0)M_{\mu }^{S(1)}{}_{}{}^{}(y)`$, average over the initial quark’s spin, sum over the final gluon’s polarization and integrate over the final partons’ momenta. One should also average over initial and sum over the final partons’ color indices. Note that one of integrations over $`x`$ is carried out by the overall momentum conservation $`2\pi \delta (x_i)`$. In the other integral $`dx/2\pi `$, we absorb the factor $`1/2\pi `$ into the definition of the parton matrix elements associated with the processes. Following this convention, we have the hard part of gluon radiation induced by single scattering,

$`\overline{H}_{\mu \nu }^{(1)}`$ $`=`$ $`{\displaystyle 𝑑x\frac{1}{2}\mathrm{Tr}[p\gamma \gamma _\mu (xp+q)\gamma \gamma _\nu ]2\pi \frac{\delta (xx_B)}{2p^+q^{}}}`$ (178)
$`\times `$ $`{\displaystyle \frac{d^2\mathrm{}_T}{2(2\pi )^3}\frac{dz}{1z}\overline{M}^S(0)\overline{M}_{}^{S}{}_{}{}^{}(y)}`$ (179)
$`=`$ $`{\displaystyle 𝑑xH_{\mu \nu }^{(0)}(x,p,q)\frac{\alpha _s}{2\pi }\frac{d\mathrm{}_T^2}{\mathrm{}_T^2}𝑑zC_F\frac{2}{1z}e^{i(x+x_L)p^+y^{}}}.`$ (180)
One can see that the above is the same as the hard part in Eq. (28) in the soft radiation limit $`z1`$.

### B Double Scattering

For double scattering processes as shown in Fig. 14(b)-(d), there is one loop-integration which is carried out by the contour integration around a pole in one of the two propagators. We also assume that the initial gluon carries transverse momentum $`k_T`$. The on-shell requirement of the final quark is

$$[q+(x+x_1)p+k_T\mathrm{}]^2=2zp^+q^{}(x+x_1x_Bx_Lx_D)=0,$$
(181)
where $`x_L`$ and $`x_D`$ are defined in Eq. (51). The kinematics in the process of Fig. 14(b) only allows one choice of pole. The contour integration around this pole gives

$$\frac{dx}{2\pi }\frac{e^{ixp^+(y^{}y_1^{})i(x_B+x_L+x_D)p^+y_1^{}}}{(xp+q)^2+iϵ}=\frac{i}{2p^+q^{}}e^{ix_Bp^+y^{}i(x_L+x_D)p^+y_1^{}}\theta (y^{}y_1^{}).$$
(182)
The momentum fraction carried by the initial gluon is then $`x_1=x_L+x_D`$. Since this momentum fraction is finite when $`k_T=0`$, the secondary scattering is considered hard. The final gluon in this process is from the final state radiation of the secondary scattering. This is what we refer to as double hard scattering. Using the helicity selection rule in Eq. (170), one can find the amplitude for the process in Fig. 14(b),

$`\overline{M}^{D(1)}(y,y_1)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}T_cT_{a_1}e^{i(x_B+x_L)p^+y^{}ix_Dp^+y_1^{}}`$ (183)
$`\times `$ $`e^{ix_Lp^+(y^{}y_1^{})}ig\theta (y^{}y_1^{}),`$ (184)
where $`a_1`$ is the color index of the initial gluon. Note that there should also be an integral over $`k_T`$. We have absorbed this together with the transverse phase factor $`\mathrm{exp}(i\stackrel{}{k_T}\stackrel{}{y}_{1T})`$ into the definition of parton matrix elements.

In the other two diagrams of double scattering in Fig. 14, there are two possible choices of poles. One choice of poles corresponds to soft rescattering of the final quark \[Fig. 14(c)\] or of the final gluon \[Fig. 14(d)\] where the initial gluon’s momentum fraction $`x_1=x_D`$ \[for Fig. 14(c)\] or $`x_1=zx_D/(1z)`$ \[for Fig. 14(d)\] goes to zero for $`k_T=0`$. The final gluon is induced by the hard photon-quark scattering. This is what we call hard-soft double scattering. The second choice however gives the initial gluon finite momentum $`x_1=x_L+x_D`$ and the gluon radiation is from the initial state radiation of the hard secondary scattering in what we call double-hard scattering. One can find the amplitudes of these two diagrams as

$`\overline{M}^{D(2)}(y,y_1)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}T_{a_1}T_ce^{i(x_B+x_L)p^+y^{}ix_Dp^+y_1^{}}`$ (185)
$`\times `$ $`[1e^{ix_Lp^+(y^{}y_1^{})}]ig\theta (y^{}y_1^{}),`$ (186)
$`\overline{M}^{D(3)}(y,y_1)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}[T_{a_1},T_c]e^{i(x_B+x_L)p^+y^{}ix_Dp^+y_1^{}}`$ (187)
$`\times `$ $`[e^{ix_Dp^+(y^{}y_1^{})/(1z)}e^{ix_Lp^+(y^{}y_1^{})}]ig\theta (y^{}y_1^{}).`$ (188)
Since the gluon radiation induced by the hard rescattering in these two diagrams is initial state radiation, the amplitude has the opposite sign with respect to the hard-soft processes where the gluon is produced by final state radiation of the first hard scattering.

### C Triple Scattering

For gluon radiation associated with triple scattering, there are all together 7 different diagrams as shown in Fig. 15. we list the amplitudes here:

$`\overline{M}^{T(1)}(y,y_1,y_2)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}T_cT_{a_2}T_{a_1}e^{i(x_B+x_L)p^+y^{}ix_D^0p^+(y_1^{}y_2^{})}`$ (189)
$`\times `$ $`e^{ix_Lp^+(y^{}y_2^{})}(g^2)\theta (y_1^{}y_2^{})\theta (y^{}y_1^{}),`$ (190)
$`\overline{M}^{T(2)}(y,y_1,y_2)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}T_{a_2}T_cT_{a_1}e^{i(x_B+x_L)p^+y^{}ix_Dp^+(y_1^{}y_2^{})}`$ (191)
$`\times `$ $`[e^{ix_Lp^+(y^{}y_1^{})}e^{ix_Lp^+(y^{}y_2^{})i(x_D^0x_D)p^+(y_1^{}y_2^{})}]`$ (192)
$`\times `$ $`(g^2)\theta (y_1^{}y_2^{})\theta (y^{}y_1^{}),`$ (193)
$`\overline{M}^{T(3)}(y,y_1,y_2)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}T_{a_2}T_{a_1}T_ce^{i(x_B+x_L)p^+y^{}ix_Dp^+(y_1^{}y_2^{})}`$ (194)
$`\times `$ $`[1e^{ix_Lp^+(y^{}y_1^{})}](g^2)\theta (y_1^{}y_2^{})\theta (y^{}y_1^{}),`$ (195)
$`\overline{M}^{T(4)}(y,y_1,y_2)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}[T_{a_2},T_c]T_{a_1}e^{i(x_B+x_L)p^+y^{}ix_Dp^+(y_1^{}y_2^{})}`$ (196)
$`\times `$ $`[e^{i(x_D^0x_D)p^+(y_1^{}y_2^{})+ix_Lp^+(y^{}y_2^{})}`$ (197)
$``$ $`e^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})+ix_Lp^+(y^{}y_1^{})}](g^2)\theta (y_1^{}y_2^{})\theta (y^{}y_1^{}),`$ (198)
$`\overline{M}^{T(5)}(y,y_1,y_2)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}T_{a_2}[T_{a_1},T_c]e^{i(x_B+x_L)p^+y^{}ix_Dp^+(y_1^{}y_2^{})}`$ (199)
$`\times `$ $`[e^{ix_Lp^+(y^{}y_1^{})}e^{ix_Dp^+(y^{}y_1^{})/(1z)}](g^2)\theta (y_1^{}y_2^{})\theta (y^{}y_1^{}),`$ (200)
$`\overline{M}^{T(6)}(y,y_1,y_2)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}T_{a_1}[T_{a_2},T_c]e^{i(x_B+x_L)p^+y^{}ix_Dp^+(y_1^{}y_2^{})}`$ (201)
$`\times `$ $`\left[e^{ix_Lp^+(y^{}y_1^{})}e^{ix_Dp^+(y^{}y_1^{})/(1z)}\right]e^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})}`$ (202)
$`\times `$ $`(g^2)\theta (y_1^{}y_2^{})\theta (y^{}y_1^{}),`$ (203)
$`\overline{M}^{T(7)}(y,y_1,y_2)`$ $`=`$ $`2g{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\stackrel{}{\mathrm{}}_T^2}}[T_{a_1},[T_{a_2},T_c]]e^{i(x_B+x_L)p^+y^{}ix_Dp^+(y_1^{}y_2^{})}`$ (204)
$`\times `$ $`[1e^{ix_Lp^+(y^{}y_1^{})}]e^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})}`$ (205)
$`\times `$ $`(g^2)\theta (y_1^{}y_2^{})\theta (y^{}y_1^{}).`$ (206)
In $`\overline{M}^{T(4)}`$, $`\overline{M}^{T(6)}`$ and $`\overline{M}^{T(7)}`$, where the radiated gluon is coupled to the second initial gluon through a triple-gluon vertex, we have made variable change $`k_Tk_T`$.

To obtain these amplitudes, two loop-integrations are carried out by contour integrations around two poles. For each diagram, there are sometimes two possible choices of the two poles, and consequently two contributions to the amplitude. They represent soft or hard rescattering processes. However, one of the quark-gluon scatterings must be soft. In Fig. 15(a), for example, the gluon induced by the third (hard) scattering after the quark has a soft second scattering. Fig 15(b) has two contributions. One corresponds to gluon radiation induced by the third scattering, the other by the second scattering. One of the two quark-gluon scatterings must be soft. In Fig 15(c), one contribution comes from gluon radiation induced by the second scattering followed by soft quark-gluon scattering. The other contribution corresponds to gluon radiation from the first photon-quark hard scattering and double soft scattering afterwards. Such analyses of each diagram help us to reorganize the amplitude according to the physical processes.

### D Soft and Hard Rescattering

We can reorganize all the radiation amplitudes of double and triple scattering according to our classification of soft and hard rescattering. The total amplitude for the hard-soft double scattering in Figs. 14(b)-(d), as represented by Fig. 7, is

$`\overline{M}_S^D(y,y_1)`$ $`=`$ $`e^{i(x_B+x_L)p^+y^{}ix_Dp^+y_1^{}}ig\theta (y^{}y_1^{})`$ (207)
$`\times `$ $`2g\left[{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}T_{a_1}T_c{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}[T_{a_1},T_c]e^{ix_Dp^+(y^{}y_1^{})/(1z)}\right].`$ (208)
The amplitude for double hard scattering, as represented by Fig. 8, is,

$`\overline{M}_H^D(y,y_1)`$ $`=`$ $`e^{i(x_B+x_L)p^+y^{}ix_Dp^+y_1^{}}ig\theta (y^{}y_1^{})`$ (209)
$`\times `$ $`2g[T_c,T_{a_1}]\left[{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}\right]e^{ix_Lp^+(y^{}y_1^{})}.`$ (210)
In triple scattering processes, the gluon radiation can be induced by the initial photon-quark hard scattering which is then followed by two soft scatterings. We denote the gluon radiation together with the first soft quark-gluon or gluon-gluon scattering by the effective diagram in Fig. 7. Similarly, the second soft scattering can be either quark-gluon or gluon-gluon as shown in Fig. 10(d) and (e). The amplitudes for these two diagrams with double soft scattering are,

$`\overline{M}_{S(1)}^T(y,y_1,y_2)`$ $`=`$ $`e^{i(x_B+x_L)p^+y^{}ix_Dp^+(y_1^{}y_2^{})}(g^2)\theta (y_1^{}y_2^{})\theta (y^{}y_1^{})`$ (211)
$`\times `$ $`2g\left[{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}T_{a_2}T_{a_1}T_c{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}T_{a_2}[T_{a_1},T_c]e^{ix_Dp^+(y^{}y_1^{})/(1z)}\right],`$ (212)
$`\overline{M}_{S(2)}^T(y,y_1,y_2)`$ $`=`$ $`e^{i(x_B+x_L)p^+y^{}ix_Dp^+(y_1^{}y_2^{})}(g^2)\theta (y_1^{}y_2^{})\theta (y^{}y_1^{})`$ (213)
$`\times `$ $`2g[{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}T_{a_1}[T_{a_2},T_c]e^{ix_Dp^+(y^{}y_1^{})/(1z)}`$ (214)
$`+`$ $`{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}[T_{a_1},[T_{a_2},T_c]]]e^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})},`$ (215)
where $`\overline{M}_{S(1)}^T`$ is the combination of the first term in $`\overline{M}^{T(3)}`$ and the second term in $`\overline{M}^{T(5)}`$ while $`\overline{M}_{S(2)}^T`$ is the combination of the second term of $`\overline{M}^{T(6)}`$ and the first term of $`\overline{M}^{T(7)}`$.

In hard rescattering, the gluon radiation is induced either by the second or the third scattering. The amplitude for induced radiation by the third scattering as shown in Fig. 10(c) is

$`\overline{M}_{H(1)}^T(y,y_1,y_2)`$ $`=`$ $`e^{i(x_B+x_L)p^+y^{}ix_D^0p^+(y_1^{}y_2^{})}()g^2\theta (y_1^{}y_2^{})\theta (y^{}y_1^{})`$ (216)
$`\times `$ $`2g[T_c,T_{a_2}]T_{a_1}\left[{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}\right]e^{ix_Lp^+(y^{}y_2^{})}.`$ (217)
It is the combination of $`\overline{M}^{T(1)}`$, the second term in $`\overline{M}^{T(2)}`$, and the first term in $`\overline{M}^{T(4)}`$. If the gluon is induced by the second scattering, the third soft scattering can either be quark-gluon or gluon-gluon as shown in Figs. 10(a) and (b). Their amplitudes are

$`\overline{M}_{H(2)}^T(y,y_1,y_2)`$ $`=`$ $`e^{i(x_B+x_L)p^+y^{}ix_Dp^+(y_1^{}y_2^{})}()g^2\theta (y_1^{}y_2^{})\theta (y^{}y_1^{})`$ (218)
$`\times `$ $`2gT_{a_2}[T_c,T_{a_1}]\left[{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}\right]e^{ix_Lp^+(y^{}y_1^{})},`$ (219)
$`\overline{M}_{H(3)}^T(y,y_1,y_2)`$ $`=`$ $`e^{i(x_B+x_L)p^+y^{}ix_Dp^+(y_1^{}y_2^{})}()g^2\theta (y_1^{}y_2^{})\theta (y^{}y_1^{})`$ (220)
$`\times `$ $`2g[T_{a_1},[T_c,T_{a_2}]]\left[{\displaystyle \frac{\stackrel{}{ϵ}_T\stackrel{}{\mathrm{}}_T}{\mathrm{}_T^2}}{\displaystyle \frac{\stackrel{}{ϵ}_T(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)}{(\stackrel{}{\mathrm{}}_T\stackrel{}{k}_T)^2}}\right]e^{ix_Lp^+(y^{}y_1^{})}`$ (221)
$`\times `$ $`e^{i(1z/(1z))x_Dp^+(y_1^{}y_2^{})},`$ (222)
respectively. Again, $`\overline{M}_{H(2)}^T`$ is the combination of the first term in $`\overline{M}^{T(2)}`$, the second term in $`\overline{M}^{T(3)}`$ and the first term of $`\overline{M}^{T(5)}`$. $`\overline{M}_{H(3)}^T`$ is the sum of the second term in $`\overline{M}^{T(4)}`$, the first term of $`\overline{M}^{T(6)}`$ and the second term of $`\overline{M}^{T(7)}`$.

With the amplitudes of the re-classified diagrams, we can easily calculate the partonic hard part of double scattering,

$$\overline{H}^D=\frac{d^2\mathrm{}_T}{2(2\pi )^3}\frac{dz}{1z}[\overline{M}^D(0,y_2)\overline{M}_{}^{D}{}_{}{}^{}(y,y_1)+\overline{M}^T(0,y_2,y_1)\overline{M}_{}^{S}{}_{}{}^{}(y)+\overline{M}^S(0)\overline{M}_{}^{T}{}_{}{}^{}(y,y_1,y_2)],$$
(223)
which includes both the double scattering and the interferences of triple and single scattering. The final results are the same as listed in Sec. III B with the splitting functions replaced by their form in the limit of $`z1`$.