text:
# Multiple Parton Scattering in Nuclei: Twist-Four Nuclear Matrix Elements and Off-Forward Parton Distributions

## I Introduction

The success of perturbative QCD (pQCD) in describing hard processes in high-energy collisions is mainly attributed to the asymptotic freedom of QCD  at short distances and to factorization theorems. Specifically, the cross sections of processes that involve large momentum transfers can be factorized into a convolution of perturbative hard scattering cross sections and nonperturbative parton distributions and fragmentation functions that contain long distance physics. Even though they are not calculable within pQCD, these parton distributions and fragmentation functions can be rigorously defined in QCD independently of any specific process and measured in many different experiments. Such factorization has been proven up to next-to-leading twist (twist-four)  for hard processes involving both hadrons and nuclei. We will refer to this as the generalized factorization.

The leading twist-four contributions in hard processes in nuclei normally involve multiple scattering with partons from different nucleons. They generally depend on twist-four nuclear parton matrix elements such as

$`{\displaystyle \frac{dy^{}}{2\pi }𝑑y_1^{}𝑑y_2^{}}`$ $`e^{ix_1p^+y^{}+ix_2p^+(y_1^{}y_2^{})}\theta (y_2^{})\theta (y^{}y_1^{})`$ (1)
$`\times `$ $`{\displaystyle \frac{1}{2}}A|\overline{\psi }_q(0)\gamma ^+F_\sigma ^+(y_2^{})F^{+\sigma }(y_1^{})\psi _q(y^{})|A,`$ (2)
which describes the quark-gluon correlation in a nucleus. This matrix element also appears in both lepton-nucleus deeply-inelastic scattering (DIS)  and in Drell-Yan cross section of $`pA`$ collisions . We will work in the infinite momentum frame, where the four-momentum of the virtual photon and the nucleus (atomic number $`A`$) have the form

$`q`$ $`=`$ $`[Q^2/2q^{},q^{},\stackrel{}{0}_{}],`$ (3)
$`p_A`$ $`=`$ $`A[p^+,0,\stackrel{}{0}_{}],`$ (4)
respectively. The Bjorken variable is then $`x_B=Q^2/2p^+q^{}`$. Our convention for four-vectors is $`k^\mu =[k^+,k^{},\stackrel{}{k}_{}]`$, where

$$k^+\frac{k^0+k^3}{\sqrt{2}},k^{}\frac{k^0k^3}{\sqrt{2}}.$$
(5)
Assuming that the two gluon fields in the rescattering process associated with Eq. (2) come from the same nucleon in the nucleus due to color confinement, it has been argued  that the above twist-four nuclear matrix elements are enhanced by a factor of $`A^{1/3}`$ as compared to the leading twist quark distributions in a nucleus,

$$f_q^A(x)=\frac{dy^{}}{2\pi }e^{ixp^+y^{}}\frac{1}{2}A|\overline{\psi }_q(0)\gamma ^+\psi _q(y^{})|A,$$
(6)
for $`A1`$. For processes involving a large transverse momentum scale $`\mathrm{}_T^2>Q^2`$, the ratio of the twist-four contribution and the leading twist one is therefore proportional to $`\alpha _sA^{1/3}/\mathrm{}_T^2`$. For large values of $`A`$, where the above analysis is valid, this quantity can be related to an expansion parameter. In this sense, the above matrix element is the leading higher-twist contribution to hard processes involving multiple parton scattering in nuclei.

In a recent study , Guo and Wang extended the generalized factorization approach to the problem of parton energy loss and modified quark fragmentation in DIS off a nuclear target due to gluon bremsstrahlung induced by secondary quark-gluon scatterings. Because of the Landau-Pomeranchuck-Midgal (LPM) interference effect, gluon bremsstrahlung with small transverse momentum, or large formation time ($`\tau _fQ^2/M\mathrm{}_T^2`$ in the nucleus rest frame, where $`M`$ is the nucleon mass), is suppressed. This limits the available phase space of the transverse momentum to $`\mathrm{}_T^2>Q^2/MR_AQ^2/A^{1/3}`$, ensuring the validity of the leading logarithmic approximation in the study of jet fragmentation for $`\mathrm{}_T^2Q^2`$ in a large nucleus with $`A^{1/3}1`$. The twist-four contribution to the modified fragmentation function in this case is proportional to $`\alpha _sA^{2/3}/Q^2`$, which depends quadratically on the nuclear size. Such a novel quadratic nuclear size dependence has recently been verified by the HERMES experiment .

Similar to other twist-four processes in a nucleus, the nuclear modification to the fragmentation function is also proportional to twist-four nuclear parton matrix elements. The LPM interference effect is explicitly embedded in the combined twist-four parton matrix elements. The quadratic nuclear size dependence of the modification to the fragmentation function is based on a generalized assumption that the twist-four parton matrix elements factorize into twist-two parton distributions in nucleons . The same approximation has been assumed for the momentum-crossed twist-four parton elements of a nucleus. This is a crude assumption at best, and does not specify the condition of validity nor provide any insight into the the relationship between nuclear and nucleonic parton distributions, which should depend on the nucleon wavefunction inside a nucleus.

It is difficult to determine the validity of this approximation within the framework of pQCD since nucleons do not appear explicitly in this theory. However, many hybrid models which employ nonrelativistic quantum mechanics to define a wavefunction for nucleons in a nucleus have been developed in the literature . In these models, the wavefunction allows one to decompose nuclear parton distributions into nucleonic ones. Phenomenologically, these models assert that scattering processes involving nuclei can be understood as weighted averages over scattering processes involving nucleons. Interactions among nucleons are reflected in the wave function. The simplest of these models  arrives at the relation

$$f_{a/A}(x)=A_x^A\frac{d\alpha }{\alpha }\rho (\alpha )f_{a/N}(x/\alpha )$$
(7)
between the leading-twist distribution of $`a`$-type partons in a nucleus of size $`A`$ and the same distribution found in a nucleon. The correspondence is made via a light-cone nucleon density function $`\rho (\alpha )`$, which is the probability of finding a nucleon in a nucleus with longitudinal momentum fraction $`\alpha `$, normalized to 1. Experimentally, Eq.(7) is approximately satisfied for $`x>0.1`$. At smaller values of $`x`$, the phenomenon of shadowing prevents realization of this naive model .

For twist-four matrix elements, one expects similar results. As twist-four objects are associated with partonic correlations, one expects two types of contributions : effects associated with partonic correlations within a single nucleon and those associated with nucleonic correlations within the nucleus. The former effects, which involve nucleonic twist-four distributions, are a simple extension of the convolution model; one simply substitutes new distributions for the twist-two ones. Contributions from nucleonic correlatons imply multiple scattering within the nucleus, and are inherently new effects. In particular, the distribution $`\rho (\alpha )`$ in Eq.(7) will be replaced by a more complicated distribution describing two-nucleon correlations, or the momentum-sharing between nucleons in the nucleus. Since these latter effects involve two nucleons, they will be enhanced by factors of the nuclear radius relative to the nucleonic higher-twist effects. This makes them the dominant contributions in the limit of large $`A`$.

The purpose of this paper is to show that an analysis of these effects in the spirit of the convolution model reveals contributions from a more general class of twist-two matrix elements. Specifically, the decomposition necessarily contains off-forward parton distributions (OFPD’s)  rather than just the simple diagonal matrix elements. At present, deeply-virtual Compton scattering (DVCS) is the only known process which involves the OFPD’s explicitly . We will demonstrate that these elusive objects could in principle also be probed in multiple parton scattering processes in a nucleus.

This paper is organized as follows. In Section II we will give a brief review of multiple parton scattering in nuclei and modified fragmentation functions, focusing on their relation to the twist-four parton matrix elements in nuclei. Section III presents the main result of this paper, in which we derive a relation between a twist-four nuclear matrix element and a convolution of two nucleonic OFPD’s. In Section IV, we discuss the properties of our result in certain limits. We will find that an analytic relationship between the simple factorized expression mentioned above and our result is not obvious, but numerical models show that they produce the same results in the limit of sharply-peaked nuclear wave functions. Section V contains some discussion about and conclusions of our results.

## II Multiple parton scattering

DIS on a nuclear target is the simplest environment in which to study the problem of multiple parton scattering in a nucleus. In this case, a quark is struck by an energetic virtual photon and then scatters again with partons from other nucleons inside the nucleus. The rescattering will induce gluon bremsstrahlung by the propagating quark and cause the leading quark to lose energy. Such radiative energy loss will be manifested in the modification of the quark fragmentation function as compared to the one measured in DIS off a nucleon target, where there is no such rescattering. The gluon bremsstrahlung will interfere destructively with the final-state radiation of the quark-photon scattering. This LPM interference effect will give rise to some novel nuclear effects in the modified quark fragmentation function.

Applying the generalized factorization of twist-four processes to the exclusive process of hadron production in DIS on a nuclear target and performing a collinear expansion with respect to initial parton transverse momentum, one can obtain an effective modified quark fragmentation function with leading higher-twist contributions :

$`\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)`$ (8)
$`+`$ $`\mathrm{\Delta }\gamma _{qgq}(z,x,x_L,\mathrm{}_T^2)D_{gh}(z_h/z)].`$ (9)
Here, $`D_{qh}(z_h,\mu ^2)`$ is the usual renormalized twist-two quark fragmentation function in vacuum that satifies the normal DGLAP  QCD evolution equation. The additional terms are the leading higher-twist contributions from multiple parton scattering and induced gluon bremsstrahlung. These contributions are very similar in form to the normal gluon radiation in vacuum except that the modified splitting functions,

$`\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^2N_cf_q^A(x)}},`$ (10)
$`\mathrm{\Delta }\gamma _{qgq}(z,x,x_L,\mathrm{}_T^2)`$ $`=`$ $`\mathrm{\Delta }\gamma _{qqg}(1z,x,x_L,\mathrm{}_T^2),`$ (11)
depend on the twist-four two-parton correlation function

$`T_{qg}^A(x,x_T,x_L)`$ $``$ $`{\displaystyle \frac{dy^{}}{4\pi }𝑑y_1^{}𝑑y_2^{}\theta (y_2^{})\theta (y^{}y_1^{})(1e^{ix_Lp^+(y_1^{}y^{})})(1e^{ix_Lp^+y_2^{}})}`$ (13)
$`\times e^{i(x+x_L)p^+y^{}+ix_Tp^+(y_1^{}y_2^{})}A|\overline{\psi }_q(0)\gamma ^+F_\sigma ^{a+}(y_2^{})F_a^{+\sigma }(y_1^{})\psi _q(y^{})|A,`$
where $`x_L=\mathrm{}_T^2/2p^+q^{}z(1z)`$ and $`x_Tk_T^2/2p^+q^{}z`$. The virtual corrections supply the $`\delta `$-function contribution to the ‘+’-function, along with the explicit end point contribution

$`\mathrm{\Delta }T_{qg}^A(x,\mathrm{}_T^2)`$ $``$ $`{\displaystyle _0^1}𝑑z{\displaystyle \frac{1+z^2}{(1z)_+}}T_{qg}^A(x,x_L)`$ (14)
$`=`$ $`{\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],`$ (15)
required for conservation of quark flavor.

The twist-four parton matrix elements are in principle not calculable and can only be measured in experiments, just like twist-two parton distributions. However, under certain assumptions, one can use some nuclear model to relate them to twist-two parton distributions in nucleons. Along the way, one obtains the $`A`$-dependence of these nuclear matrix elements.

If we assume that the nuclear wave function can be expressed as a multiple-nucleon state, with each nucleon a color singlet, the two gluon fields must operate on the same nucleon state inside the nucleus. Consider the dominant case where the quark and gluon fields operate on different nucleons inside the nucleus. The integration over $`y_1^{}`$ and $`y_2^{}`$ in Eq. (13) should give the length scale $`r_NR_A`$, where $`r_N`$ is the nucleon radius and $`R_A1.12A^{1/3}\mathrm{fm}`$ is the radius of the nucleus. The twist-four nuclear parton matrix elements should then be approximately proportional to $`A^{4/3}`$. If the quark and gluon come from the same nucleon, the matrix elements will only be proportional to $`A`$, which is subleading.

The twist-four two-parton correlation function that enters the modified quark fragmentation function has not only normal parton matrix elements representing direct terms in the square of the amplitude, but also those representing interference. The former has momenta flowing directly along two parton fields separately while the latter has momenta flowing across two different parton fields. These two different contributions were called ‘diagonal’ and ‘off-diagonal’ matrix elements, respectively, in Ref. . We call them ‘direct’ and ‘crossed’ here to avoid confusion with truly off-diagonal matrix elements, in which the momenta of the external states are different. The relative signs between these two kinds of matrix elements reflect the physics of the LPM interference effect in the processes of induced gluon radiation. As illustrated by the central cut-diagram in Fig. 1, the gluon radiation can either be produced as final state radiation of the photon-quark hard scattering or initial state radiation of the quark-gluon rescattering. In the former case, the energy of the radiated gluon is provided by the initial quark with $`x=x_B+x_L`$. The quark-gluon rescattering can be very soft since the momentum fraction carried by the gluon, $`x_3=x_T=k_T^2/2p^+q^{}z`$, is very small when $`k_T^20`$. This is normally referred to as a hard-soft process. In the latter case, however, the initial gluon must carry a finite momentum fraction $`x_3=x_T+x_L`$ to induce the gluon radiation. Such a process is called double hard scattering. Contributions from central-cut diagrams such as Fig. 1 contain both of these processes as well as their interference. Since initial and final state radiation amplitudes have a phase difference of $`\pi `$ , the final result of the sum of these contributions is the dipole-like form-factor in the radiation spectrum which is now absorbed into the definition of the two-parton correlation $`T_{qg}^A(x,x_T,x_L)`$. As $`x_Lp^+0`$, the effective gluon radiation spectrum vanishes because the interference becomes complete. This is exactly the LPM effect , which is now embedded in the effective two-parton correlation function.

Because of the LPM interference effect in induced bremsstrahlung, the effective two-parton correlation function that enters the modified fragmentation function essentially contains four independent twist-four nuclear parton matrix elements. The two direct ones correspond to gluon radiation associated with photon-quark and quark-gluon scattering. In the corresponding forward scattering processes, shown in Fig. 1, momentum flows separately along the quark and gluon lines. The twist-four parton matrix elements can then have the interpretation of a two-parton joint distribution inside the nucleus. The two crossed matrix elements are related to the inteference between the two different radiation processes. In this case, there is actually a momentum flow in the amount $`x_Lp^+`$ between quarks and gluons in the forward scattering amplitude. Such crossed matrix elements do not have the interpretation of a normal parton distribution. Since we only consider the case where the quark and gluon are from different nucleons inside the nucleus, these contributions should be related to the off-forward (or skewed) parton distrbution functions of a nucleon. In the limit of vanishing skewedness ($`x_L0`$), the crossed matrix elements approach the direct ones. One can see by inspection that the direct matrix elements are real, while the crossed ones are complex.

In the rest of this paper, we study these twist-four nuclear parton matrix elements within the convolution model and relate them to generalized nucleonic parton distributions. The matrix element at issue,

$`K(x_1,x_2,x_L)`$ $`=`$ $`{\displaystyle \frac{dy^{}}{4\pi }𝑑y_1^{}𝑑y_2^{}\theta (y_2^{})\theta (y^{}y_1^{})}`$ (17)
$`\times e^{ix_1p^+y^{}}e^{ix_2p^+(y_1^{}y_2^{})}e^{ix_Lp^+y_1^{}}A|\overline{\psi }(0)\gamma ^+F_\sigma ^{a+}(y_2^{})F_a^{+\sigma }(y_1^{})\psi (y^{})|A,`$
describes the removal of a quark with momentum $`x_1p`$ and a gluon of momentum $`(x_2+x_L)p`$ from our nuclear state $`|A`$, and the subsequent replacement of a gluon with momentum $`x_2p`$ and a quark of momentum $`(x_1+x_L)p`$, as illustrated in Fig. 2. This nuclear parton matrix element is useful in constructing the physical combinations that appear in many nuclear scattering processes , as well as the in-medium evolution of the parton fragmentation functions . In particular, the correlation function (13) appearing in the medium-modified quark fragmentation function (9) can be written as

$`T_{qg}^A(x,x_T,x_L)`$ $`=`$ $`K(x+x_L,x_T,0)K(x,x_T,x_L)`$ (19)
$`K(x+x_L,x_T+x_L,x_L)+K(x,x_T+x_L,0).`$
## III The Convolution Model

In order to relate the matrix element (17) to the nucleonic degrees of freedom of the nucleus, we must define these degrees of freedom quantitatively. Our formalism is based on light-cone perturbation theory , in which one expands the physical state under consideration in terms of the free-particle states of its constituents. In our case, we would like to express our nuclear state in terms of free nucleonic states. Phenomenologically, we know that this basis is not complete and therefore cannot be guaranteed to span our nuclear state space.

Contributions from higher Fock states are required in quantum field theory to generate effective nonlocal interactions from the underlying contact terms of the Lagrangian. However, most of the effect of higher Fock states is already included in the definition of the twist-two nucleonic parton matrix elements. The additional effects of higher Fock states in a nucleus are induced by nucleonic interactions. However, these effects are limited by the effective interaction energy of the system under consideration. Here, we are mainly concerned with nuclear systems whose interaction energy is small compared with the energy per nucleon. This allows us to consider only the lowest Fock state, that of $`A`$ nucleons. We mention here that although we neglect interactions in the form of higher Fock states, nucleon correlations will still occur through the wave function in our Hilbert space. These correlations lead to nontrivial relations between the nuclear matrix element and the nucleonic distributions.

In light of the above approximation, we write

$$|A𝑑\mathrm{\Pi }_A\varphi (\{p_i\})|\{p_i\}\mathrm{\hspace{0.17em}2}p^+(2\pi )^3\delta ^{(3)}\left(p_A\underset{i=1}{\overset{A}{}}p_i\right),$$
(20)
where $`|\{p_i\}`$ represents the state of $`A`$ free nucleons of momenta $`\{p_i\}`$ normalized as

$$\{p_i^{}\}|\{p_i\}=\underset{i}{}2p_i^+(2\pi )^3\delta ^{(3)}(p_ip_i^{}),$$
(21)
and

$$d\mathrm{\Pi }_A\underset{i=1}{\overset{A}{}}\left\{\frac{d^3p_i}{(2\pi )^3}\frac{\theta (p_i^+)}{2p_i^+}\right\}$$
(22)
represents the differential phase-space of $`A`$ nucleons.

The nucleon states are specified by their ‘+’ and ‘$``$’ momentum components; $`d^3p_idp_i^+d^2p_i`$. The ‘–’ component of each momentum is determined by the on-shell condition,

$$p_i^{}=\frac{p_i^2+M^2}{2p_i^+},$$
(23)
and is not conserved (i.e. $`_ip_i^{}p_A^{}`$). Again, $`M`$ is the nucleon mass. The normalization of our nuclear state, $`A^{}|A=2p^+(2\pi )^3\delta ^{(3)}(p_Ap_A^{})`$ implies

$$𝑑\mathrm{\Pi }_A|\varphi (\{p_i\})|^2\mathrm{\hspace{0.17em}2}p^+(2\pi )^3\delta ^{(3)}\left(p_A\underset{i=1}{\overset{A}{}}p_i\right)=1.$$
(24)
If we were to consider higher Fock states as well, the 1 on the right-hand side of this equation would be replaced by the probability of finding our state in this lowest Fock state.

Our wave function $`\varphi `$ contains many distribution functions. In particular, we can define the one-nucleon density

$`\rho _1(k){\displaystyle 𝑑\mathrm{\Pi }_{A1}|\varphi (k,\{p_i\})|^2}`$ (25)
$`\times 2p^+(2\pi )^3\delta \left(Ap^+k^+{\displaystyle \underset{i=1}{\overset{A1}{}}}p_i^+\right)\delta ^{(2)}\left(k_{}+{\displaystyle \underset{i=1}{\overset{A1}{}}}p_i\right),`$ (26)
which represents the probability of finding a nucleon of momentum $`k`$ in our nucleus irrespective of the momenta of the other nucleons. In light of Eq. (24), $`\rho _1`$ satisfies

$$\frac{d^3k}{(2\pi )^3}\frac{\theta (k^+)}{2k^+}\rho _1(k)=1.$$
(27)
Hence the light-cone nucleon distribution function in Eq.(7) can be written as

$$\rho (\alpha )=\frac{1}{4\pi \alpha }\frac{d^2k_{}}{(2\pi )^2}\rho _1(k),$$
(28)
where $`\alpha k^+/p^+`$. Eq.(7) comes directly from the substitution of (20) into the definition of the nuclear quark distribution function $`f_{q/A}(x)`$ in Eq. (6).

The two-nucleon correlator,

$`\rho _2(k_1;k_2;\mathrm{\Delta }){\displaystyle 𝑑\mathrm{\Pi }_{A2}\varphi ^{}(k_1\mathrm{\Delta }/2,k_2+\mathrm{\Delta }/2,\{p_i\})\varphi (k_1+\mathrm{\Delta }/2,k_2\mathrm{\Delta }/2,\{p_i\})}`$ (29)
$`\times \mathrm{\hspace{0.33em}2}p^+(2\pi )^3\delta \left(Ap^+k_1^+k_2^+{\displaystyle \underset{i=1}{\overset{A2}{}}}p_i^+\right)\delta ^{(2)}\left(k_1+k_2+{\displaystyle \underset{i=1}{\overset{A2}{}}}p_i\right),`$ (30)
contains information about the sharing of momentum by nucleons inside the nucleus and appears in the double scattering process we consider. The two-nucleon density, $`\rho _2(k_1;k_2;0)`$, represents the probability of finding two nucleons with the specified momenta within the nucleus, and is normalized as

$$𝑑\mathrm{\Pi }_2\rho _2(k_1;k_2;0)=1.$$
(31)
The use of this function rather than the two-parton correlator, $`\rho _2`$, leads to the expectation of diagonal parton distributions in twist-four nuclear matrix elements.

Due to the $`\theta `$-functions, our matrix element cannot readily be interpreted as a product of twist-two distributions as it stands. These $`\theta `$-functions order the fields along the light-cone axis to make the multiple scattering process physical. Employing the representation

$$\theta (y)=\frac{1}{2\pi i}𝑑z\frac{1}{zx\pm i\epsilon }e^{i(zx)y}$$
(32)
for the $`\theta `$-function, our nonperturbative distribution takes the form

$`K(x_1,x_2,x_L)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{dz_1}{z_1\omega _1i\epsilon }\frac{dz_2}{z_2\omega _2+i\epsilon }}`$ (35)
$`\times {\displaystyle }{\displaystyle \frac{dy^{}}{2\pi }}{\displaystyle \frac{dy_1^{}}{2\pi }}{\displaystyle \frac{dy_2^{}}{2\pi }}e^{ip^+y^{}(x_1+z_1\omega _1)}e^{ip^+y_1^{}(x_2+x_Lz_1+\omega _1)}e^{ip^+y_2^{}(z_2\omega _2x_2)}`$
$`\times A\left|\overline{\psi }(0)\gamma ^+\psi (y^{})F_\sigma ^{a+}(y_2^{})F_a^{+\sigma }(y_1^{})\right|A.`$
Here, $`\omega _1`$ and $`\omega _2`$ are arbitrary real variables that will be chosen later to simplify the final results.

Substitution of our approximate nuclear state (20) into (35) leads to matrix elements of the form

$$\{p_i^{}\}|\overline{\psi }(0)\gamma ^+\psi (y^{})F_\sigma ^{a+}(y_2^{})F_a^{+\sigma }(y_1^{})|\{p_i\}.$$
(36)
Assuming that the color correlation length along the light-cone within our nucleus is not larger than the nucleon size and neglecting the effects of direct multi-nucleonic correlations as higher twist, we can factorize this expression into a product of single-particle Hilbert space amplitudes:

$`Ap_1^{}|\overline{\psi }(0)\gamma ^+\psi (y^{})F_\sigma ^{a+}(y_2^{})F_a^{+\sigma }(y_1^{})|p_1{\displaystyle \underset{i=2}{\overset{A}{}}}2p_i^+(2\pi )^3\delta ^3(p_i^{}p_i)`$ (37)
$`+\left({\displaystyle \genfrac{}{}{0pt}{}{A}{2}}\right)p_1^{}|\overline{\psi }(0)\gamma ^+\psi (y^{})|p_1p_2^{}|F_\sigma ^{a+}(y_2^{})F_a^{+\sigma }(y_1^{})|p_2{\displaystyle \underset{i=3}{\overset{A}{}}}2p_i^+(2\pi )^3\delta ^3(p_i^{}p_i).`$ (38)
Our dismissal of direct multi-nucleonic correlations is one of the main approximations in this paper, and will be discussed in more detail in the conclusions. For now, we note that these higher-twist corrections are suppressed by powers of $`Q^2`$ and as such can be neglected at large scales.

The contributions to (38) can be interpreted in terms of the multiple parton scattering picture in DIS. Matrix elements related to $`\overline{\psi }(0)\gamma ^+\psi (y^{})`$ represent the probability that a quark in a certain nucleon is struck by our probe. The struck quark then propagates through the nucleus, encountering another parton at some point during its journey. If this rescattering occurs while the struck parton is still in its parent nucleon, the effect is represented by a twist-four nucleonic matrix element convoluted with the single-nucleon density. This is essentially the first term in Eq.(38). Considering the coherence length of the scattering to be of the order of the nucleon size, this term should be proportional to $`A`$.

If the rescattering is with a parton in another nucleon, the probability is related to the twist-two parton distributions in each nucleon convoluted with the two-nucleon correlator. We can call this double-factorized rescattering. Since the nuclear radius grows with $`A^{1/3}`$, one expects this double scattering process to be proportional to $`A^{4/3}`$, enhanced by $`A^{1/3}`$ over the double scattering within the same nucleon. The first term in the above decomposition is therefore suppressed relative to the second as $`A`$ increases.

Using the above decomposition, and keeping only the contribution of the double-factorized rescattering, $`K(x_1,x_2,x_L)`$ becomes

$`K(x_1,x_2,x_L)2\left({\displaystyle \genfrac{}{}{0pt}{}{A}{2}}\right){\displaystyle \frac{dz_1}{z_1\omega _1i\epsilon }\frac{dz_2}{z_2\omega _2+i\epsilon }\frac{d^2k_1}{(2\pi )^2}\frac{d^2k_2}{(2\pi )^2}\frac{d^2\mathrm{\Delta }_{}}{(2\pi )^2}}`$ (39)
$`\times {\displaystyle }{\displaystyle \frac{d\alpha _1}{2\pi }}{\displaystyle \frac{d\alpha _2}{2\pi }}{\displaystyle \frac{d\zeta }{2\pi }}{\displaystyle \frac{\theta (\alpha _1+\zeta )}{2(\alpha _1+\zeta )}}{\displaystyle \frac{\theta (\alpha _1\zeta )}{2(\alpha _1\zeta )}}{\displaystyle \frac{\theta (\alpha _2+\zeta )}{2(\alpha _2+\zeta )}}{\displaystyle \frac{\theta (\alpha _2\zeta )}{2(\alpha _2\zeta )}}`$ (40)
$`\times \delta (x_L2\zeta +z_2z_1\omega _2+\omega _1)\rho _2(k_1;k_2;\mathrm{\Delta })`$ (41)
$`\times {\displaystyle }{\displaystyle \frac{dy^{}}{2\pi }}e^{ip^+y^{}(x_1+z_1\omega _1+\zeta )}k_1{\displaystyle \frac{\mathrm{\Delta }}{2}}\left|\overline{\psi }({\displaystyle \frac{y^{}}{2}})\gamma ^+\psi \left({\displaystyle \frac{y^{}}{2}}\right)\right|k_1+{\displaystyle \frac{\mathrm{\Delta }}{2}}`$ (42)
$`\times {\displaystyle }{\displaystyle \frac{dy_d^{}}{2\pi }}{\displaystyle \frac{1}{p^+}}e^{ip^+y_d^{}(x_2+x_Lz_1+\omega _1\zeta )}k_2+{\displaystyle \frac{\mathrm{\Delta }}{2}}\left|F_\sigma ^{a+}({\displaystyle \frac{y_d^{}}{2}})F_a^{+\sigma }\left({\displaystyle \frac{y_d^{}}{2}}\right)\right|k_2{\displaystyle \frac{\mathrm{\Delta }}{2}},`$ (43)
where $`y_d=y_1y_2`$, $`\alpha _ik_i^+/p^+`$ and $`\zeta \mathrm{\Delta }^+/2p^+`$.

If we could take $`\mathrm{\Delta }=0`$, this expression would be reduced to the ordinary parton distributions probed in deeply-inelastic scattering. As it is, we can employ the off-forward parton distributions (OFPD’s)  \[Note that the distributions used here are slightly different than those defined in . In particular, $`G(x,\xi ,t)=xF_g(x,\xi ,t)`$\]. Thus,

$`{\displaystyle \frac{dy^{}}{2\pi }e^{ip^+y^{}x}k_1\frac{\mathrm{\Delta }}{2}\left|\overline{\psi }\left(\frac{y^{}}{2}\right)\gamma ^+\psi \left(\frac{y^{}}{2}\right)\right|k_1+\frac{\mathrm{\Delta }}{2}}`$ $`=`$ $`2Q({\displaystyle \frac{x}{\alpha _1}},\zeta ,{\displaystyle \frac{t_1}{M^2}}),`$ (44)
$`{\displaystyle \frac{dy^{}}{2\pi }e^{ip^+y^{}x}k_2+\frac{\mathrm{\Delta }}{2}\left|F_\sigma ^{a+}\left(\frac{y^{}}{2}\right)F_a^{+\sigma }\left(\frac{y^{}}{2}\right)\right|k_2\frac{\mathrm{\Delta }}{2}}`$ $`=`$ $`2\alpha _2p^+G({\displaystyle \frac{x}{\alpha _2}},\zeta ,{\displaystyle \frac{t_2}{M^2}}),`$ (45)
where

$$t_i=\frac{4\zeta ^2M^2}{\alpha _i^2\zeta ^2}\frac{(2\zeta \stackrel{}{p}_i+\alpha _i\stackrel{}{\mathrm{\Delta }}_{})^2}{\alpha _i^2\zeta ^2}$$
(46)
represents the squared four-momentum transfers for the two matrix elements. Note $`t_1t_2`$ since the ‘-’-components of our four-vectors are not conserved in this version of perturbation theory.

In terms of these OFPD’s, we have

$`K(x_1,x_2,x_L)`$ $``$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \genfrac{}{}{0pt}{}{A}{2}}\right){\displaystyle \frac{d\alpha _1}{2\pi }\frac{d\alpha _2}{2\pi }\frac{d\zeta }{2\pi }\frac{\theta (\alpha _1+\zeta )}{\alpha _1+\zeta }\frac{\theta (\alpha _1\zeta )}{\alpha _1\zeta }\frac{\theta (\alpha _2+\zeta )}{\alpha _2+\zeta }\frac{\theta (\alpha _2\zeta )}{\alpha _2\zeta }}`$ (50)
$`\times {\displaystyle }{\displaystyle \frac{d^2k_1}{(2\pi )^2}}{\displaystyle \frac{d^2k_2}{(2\pi )^2}}{\displaystyle \frac{d^2\mathrm{\Delta }_{}}{(2\pi )^2}}\rho _2(k_1;k_2;\mathrm{\Delta })`$
$`\times {\displaystyle }{\displaystyle \frac{dz_1}{(z_1\omega _1i\epsilon )}}{\displaystyle \frac{dz_2}{(z_2\omega _2+i\epsilon )}}\delta (x_L2\zeta +z_2z_1\omega _2+\omega _1)`$
$`\times \alpha _2Q({\displaystyle \frac{x_1+z_1\omega _1+\zeta }{\alpha _1}},\zeta ,{\displaystyle \frac{t_1}{M^2}})G({\displaystyle \frac{x_2+x_Lz_1+\omega _1\zeta }{\alpha _2}},\zeta ,{\displaystyle \frac{t_2}{M^2}}).`$
This complicated expression can be reduced to a more enlightening form through a few plausible assumptions on the form of our nucleonic correlation $`\rho _2`$ and the OFPD’s. To begin with, we assume that $`\rho _2`$ is peaked around $`\alpha _i=1`$ and $`\zeta =\mathrm{\Delta }_{}=k_i=0`$, with widths that are governed by the nuclear radius, $`R_A`$. This ansatz is dictated by the expectation that the nucleons are confined within the nuclear radius in position space. Specifically, we write

$$\rho _2(\alpha _1,k_1;\alpha _2,k_2;2\zeta ,\mathrm{\Delta }_{})=R_A^4/x_A^2r_2(\frac{\alpha _11}{x_A},R_Ak_1;\frac{\alpha _21}{x_A},R_Ak_2;\frac{2\zeta }{x_A},R_A\mathrm{\Delta }_{}),$$
(51)
with $`x_A=1/(MR_A)`$ and $`r_2`$ approximately independent of $`A`$. The behavior of the normalization $`R_A^4/x_A^2`$ as a function of $`A`$ can be determined from the normalization condition (31) on our two-particle density. Using the $`\delta `$-function to perform the $`\zeta `$ integration and making the changes

$`\nu _i`$ $`=`$ $`(\alpha _i1)/x_A`$ (52)
$`\stackrel{}{v}_i`$ $`=`$ $`R_A\stackrel{}{k}_i`$ (53)
$`\stackrel{}{\delta }_{}`$ $`=`$ $`R_A\stackrel{}{\mathrm{\Delta }}_{}`$ (54)
$`z`$ $`=`$ $`(z_1+z_2)/2x_A`$ (55)
$`u`$ $`=`$ $`(z_1z_2)/2x_A`$ (56)
in variables, we arrive at the expression

$`K(x_1,x_2,x_L)`$ $``$ $`{\displaystyle \frac{1}{4\pi R_A^2}}\left({\displaystyle \genfrac{}{}{0pt}{}{A}{2}}\right){\displaystyle 𝑑u𝑑z\frac{1}{z+u+\xi i\epsilon }\frac{1}{zu\xi +i\epsilon }}`$ (61)
$`\times {\displaystyle }{\displaystyle \frac{d\nu _1}{2\pi }}{\displaystyle \frac{d\nu _2}{2\pi }}{\displaystyle \frac{\theta (1+x_A(\nu _1u))}{1+x_A(\nu _1u)}}{\displaystyle \frac{\theta (1+x_A(\nu _1+u))}{1+x_A(\nu _1+u)}}`$
$`\times {\displaystyle \frac{\theta (1+x_A(\nu _2u))}{1+x_A(\nu _2u)}}{\displaystyle \frac{\theta (1+x_A(\nu _2+u))}{1+x_A(\nu _2+u)}}`$
$`\times {\displaystyle }{\displaystyle \frac{d^2v_1}{(2\pi )^2}}{\displaystyle \frac{d^2v_2}{(2\pi )^2}}{\displaystyle \frac{d^2\delta _{}}{(2\pi )^2}}r_2(\nu _1,\stackrel{}{v}_1;\nu _2,\stackrel{}{v}_2;2u,\stackrel{}{\delta }_{})(1+x_A\nu _2)`$
$`\times Q({\displaystyle \frac{x_1+x_A(\xi +z)}{1+x_A\nu _1}},x_Au,\stackrel{~}{t}_1)G({\displaystyle \frac{x_2+x_A(\xi z)}{1+x_A\nu _2}},x_Au,\stackrel{~}{t}_2),`$
for (50). For simplicity, we have chosen $`\omega _1=\omega _2=x_L/2`$ and defined the parameters $`\xi x_L/(2x_A)`$ and

$$\stackrel{~}{t}_i\frac{x_A^2}{(1+x_A\nu _i)^2x_A^2u^2}\left[4u^2+((1+x_A\nu _i)\stackrel{}{\delta }_{}2x_Au\stackrel{}{v}_i)^2\right].$$
(62)
Working from this form of $`K(x_1,x_2,x_L)`$, it is easy to see the large-$`A`$ enhancement of the multiple-scattering contribution. Since $`r_2`$ depends only weakly on $`A`$, this contribution to $`K`$ scales like $`A^{4/3}`$ as expected. We can simplify our expression further by assuming that $`r_2`$ is sharply peaked in the sense that all of its moments are finite. In this case, we can expand our integrand about the peak in $`r_2`$ in the formal limit $`A\mathrm{}`$ and drop all non-leading terms :

$`K(x_1,x_2,x_L)`$ $``$ $`{\displaystyle \frac{A^2}{8\pi R_A^2}}{\displaystyle 𝑑u𝑑z\frac{1}{z+u+\xi i\epsilon }\frac{1}{zu\xi +i\epsilon }\stackrel{~}{r}_2(u)}`$ (64)
$`\times Q(x_1+x_A(\xi +z),x_Au,0)G(x_2+x_A(\xi z),x_Au,0),`$
where

$$\stackrel{~}{r}_2(u)\frac{d\nu _1}{2\pi }\frac{d\nu _2}{2\pi }\frac{d^2v_1}{(2\pi )^2}\frac{d^2v_2}{(2\pi )^2}\frac{d^2\delta _{}}{(2\pi )^2}r_2(\nu _1,v_1;\nu _2,v_2;2u,\delta _{}).$$
(65)
Expression (64) is the main result of this paper. Its derivation requires only the assumptions that the lowest Fock state dominates the nuclear wave function and that the nucleonic correlator is sharply peaked with a width dictated by the nuclear radius. Strictly speaking, this expression is valid only in the formal limit $`A\mathrm{}`$. For any finite value of $`A`$, one must investigate the size of the derivatives of the OFPD’s in relation to $`A^{1/3}`$. While this investigation must be done within a specific model, the contributions are expected to be small as long as the singular regions are avoided. Our implicit assumption that the OFPD’s are analytic functions of the virtuality of the momentum transfer, $`t_i`$, is supported by studies of these functions . On the other hand, the kernels dictating the evolution of these functions (cf ) imply that the OFPD’s are not analytic functions of their second argument. This is why we have not expanded our integrand about $`x_Au=0`$.

To put this expression into a form suitable for numerical evaluation, we write

$$\frac{1}{z+u+\xi i\epsilon }\times \frac{1}{zu\xi +i\epsilon }=\mathrm{P}\left[\frac{1}{2z}\right]\left[\frac{1}{z+u+\xi i\epsilon }+\frac{1}{zu\xi +i\epsilon }\right]$$
(66)
and explicitly separate the $`u`$-integration into its pole and principal value parts. After a change of variables, this leads to

$`K(x_1,x_2,x_L){\displaystyle \frac{A^2}{8\pi R_A^2}}\{{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{2z}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dw}{w}}[F(z,wz\xi )F(z,w+z+\xi )`$ (67)
$`+F(z,wz+\xi )F(z,w+z\xi )]`$ (68)
$`+i\pi {\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{2z}}[F(z,z+\xi )F(z,z\xi )]\},`$ (69)
where

$`F(z,u)[Q(x_1+x_L/2+x_Az,x_Au,0)G(x_2+x_L/2x_Az,x_Au,0)`$ (70)
$`+Q(x_1+x_L/2x_Az,x_Au,0)G(x_2+x_L/2+x_Az,x_Au,0)]\stackrel{~}{r}_2(u).`$ (71)
In deriving this result, we have used the facts that $`\stackrel{~}{r}_2(u)=\stackrel{~}{r}_2(u)`$ and that the OFPD’s are even functions of their second argument, as can be seen by inspection.

Several features of this expression are worth pointing out. First, we note that the imaginary part of our matrix element is odd in $`\xi `$. This causes the combination in Eq.(19) to be real, as required. One can check that the two-parton correlation $`T_{qg}^A(x,x_T,x_L)`$, as expressed in Eq. (19), depends only on the real part of $`K(x_1,x_2,x_L)`$. In addition, this combination is necessarily positive. The fact that these consistency requirements are satisfied is gratifying, but not unexpected. On the other hand, the relationship between our expression and the naive expectation

$$K(x_1,x_2,0)Q(x_1)G(x_2),$$
(72)
where $`Q(x)Q(x_1,0,0)`$ and $`G(x)G(x,0,0)`$ are the ordinary diagonal distributions, is entirely unclear. Ignoring, for the moment, the non-analytic nature of the OFPD’s, one can imagine expanding each term in the integrand about the peak of the nuclear sharing function $`\stackrel{~}{r}_2`$ and dropping all higher-order terms. This reduces the OFPD’s to ordinary parton distribution functions, bringing us closer to (72). However, we cannot formally reproduce this simple dependence because of a remaining convolution between the distributions. This convolution causes our matrix element to sample the entire parton distribution functions, regardless of the values of $`x_1`$, $`x_2`$, and $`x_L`$. Its presence is a direct consequence of the correlative nature of the matrix element. Mathematically, it comes from the $`\theta `$-functions.

We could attempt to remove the correlation by expanding our result about $`x_A=0`$. However, corrections to the leading term in our expansion diverge, indicating nonanalyticity at $`x_A=0`$. This does not mean that the leading term is not a good approximation to the full solution, only that the corrections cannot be expressed in terms of powers of $`x_A`$. Nevertheless, assuming that the corrections are small for some range of $`x_A`$ near zero, one can compute Eq. (64) via contour integration directly in the limit of $`x_A0`$, but with fixed value of $`\xi =x_L/2x_A`$. The result,

$`K(x_1,x_2,x_L)|_{x_A0}`$ $``$ $`K_0(x_1,x_2,x_L)`$ (73)
$``$ $`{\displaystyle \frac{\pi A^2}{8R_A^2}}Q(x_1+x_L/2)G(x_2+x_L/2)`$ (74)
$`\times `$ $`\left\{\stackrel{~}{r}_2(\xi )+{\displaystyle \frac{i}{\pi }}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{du}{u}}\left[\stackrel{~}{r}_2(u+\xi )\stackrel{~}{r}_2(u\xi )\right]\right\},`$ (75)
is of the form of the naive expectation. In order to determine the validity of this approximation, we must calculate $`K`$ numerically in some model and compare the results to the above approximation.

## IV Numerical Results

In this section, we explore some of the properties of expression (64) in a specific model. While this model is certainly not expected to conform to the quantitative details of the realistic nuclear wave functions, we expect its general features to be echoed in more realistic treatments.

For simplicity, we assume a Gaussian form for $`\stackrel{~}{r}_2`$ :

$$\stackrel{~}{r}_2(u)=\frac{R_A^2\mathrm{\Delta }_{}^2}{4\pi }e^{u^2}.$$
(76)
The normalization is determined via (31) in conjunction with (51) by assuming a Gaussian dependence of $`r_2`$ on $`\delta _{}`$. The constant $`R_A^2\mathrm{\Delta }_{}^2`$ is a measure of the transverse momentum-sharing among nucleons in our nucleus, and is expected to be of order one. Expanding our integrands about the peak in $`\stackrel{~}{r}_2`$ and performing the integral over $`u`$, we arrive at the expression

$`K(x_1,x_2,x_L)`$ $``$ $`{\displaystyle \frac{A^2}{32\pi ^{3/2}R_A^2}}R_A^2\mathrm{\Delta }_{}^2{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{z}}f(z,0)`$ (77)
$`\times `$ $`\left\{\left[D(z+\xi )+D(z\xi )\right]+{\displaystyle \frac{i\sqrt{\pi }}{2}}\left[e^{(z+\xi )^2}e^{(z\xi )^2}\right]\right\}.`$ (78)
Here,

$`f(z,0)Q(x_1+x_L/2+x_Az)G(x_2+x_L/2x_Az)`$ (79)
$`+Q(x_1+x_L/2x_Az)G(x_2+x_L/2+x_Az)`$ (80)
and

$$D(x)_0^x𝑑te^{t^2x^2}$$
(81)
is Dawson’s integral. Since $`D(x)1/2x`$ as $`x\mathrm{}`$, $`K(x_1,x_2,x_L)`$ is obviously non-analytic in $`x_A`$ at $`x_A=0`$, as mentioned above. However, it can be checked explicitly that this expression reduces to (75) when $`x_A0`$.

In deriving Eq. (78), we have assumed $`x_Au1`$ in the expansion around the peak of $`\stackrel{~}{r}_2(u)`$. Such an approximation does not necessarily represent the leading behavior of the integral. The pole contributions occur in the region $`u\xi `$, where $`x_Aux_L/2`$ can in fact be of the same order or larger than the first argument of our OFPD’s,

$$Q(x_1+x_L/2\pm x_Az,x_L/2,0)G(x_2+x_L/2x_Az,x_L/2,0).$$
(82)
In this case, the full generalized distributions are sampled. However, in the relevant region of the $`z`$-integration in Eq. (69) with a Gaussian form of $`\stackrel{~}{r}_2(u)`$, $`|z|<\xi `$. This causes the first variables in the OFPD’s to be bounded by $`x_1`$ and $`x_2`$, respectively. According to model studies of the OFPD’s , the OFPD’s can be approximated by the ordinary parton distributions when the first argument is larger in magnitude than the second. In addition, since the OFPD’s are continuous and $`G(x,\xi ,t)`$ is expected to be positive definite, the variation of $`G(x,\xi ,t)`$ is small in the region $`x<|\xi |1`$. Hence we expect very small deviations from the ordinary gluon momentum distribution even in the region $`x_2<<x_L<<1`$. The singular nature of the quark distribution function $`Q(x)`$ at $`x=0`$ leads to large variations of the associated OFPD in the region $`x_1<x_L`$. This will essentially be the limitation of our approximation here. Outside of our applicability region, the full generalized parton distributions are needed to predict the behavior of the matrix element.

Using the CTEQ parameterization of parton distributions from data , we can calculate our expression numerically and see how well the results follow the factorized form $`K_0`$. Since only the real parts of $`K`$ enters the two-parton correlation $`T_{qg}^A`$, we will concentrate on the real part of our matrix element. The nature of (78) is such that our matrix element samples the parton momentum distribution at all values of $`x`$ rather than just those close to $`x_{1,2}`$. Since the parton distributions are not known for $`x0`$, we assume a simple extrapolation with a constant value for gluons and set the quark distributions to zero beyond the region of parametrization. The errors introduced by such extrapolation are negligible since the contributions to the integral from this region is very small.

Figure 3 shows the ratio $`K(0.3,0.01,x_L)/K_0(0.3,0.01,x_L)`$ as a function of $`x_A/x_L`$ for three different values of $`x_A`$. Since $`x_L`$ measures the momentum sharing between nucleons in our matrix element, it should be smaller than the characteristic momentum fraction, $`x_A`$, in our nucleus. If $`x_L`$ becomes too large, the momentum transfer is suppressed beyond the simple exponential suppression in $`K_0`$. This explains the behavior of the curve for small $`x_A/x_L`$. As $`x_A/x_L`$ increases, the ratio approaches the direct contribution represented by the dashed line in the figure. As can easily be seen from the graph, there is some residual dependence on the nuclear size. However, the saturation ratio changes by less than 15% as one changes $`A`$ from 208 to 32.

To study the residual dependence of our matrix element on the nuclear size, we plot the ratio $`K/K_0`$ at $`x_L=0`$ as a function of $`x_A`$ for three different values of $`x_1`$ in Figure 4. The ratio is seen to drop monotonically as $`x_A`$ is increased, but the close proximity of the curves $`x_1=0.2`$ and $`x_1=0.3`$ indicates only slight dependence on $`x_1`$ when $`x_1`$ is moderate. As $`x_1`$ becomes smaller, the dependence on both $`x_1`$ and $`x_A`$ becomes far more dramatic. In either case, the $`x_A`$ dependence is approximately linear. Since we have already dropped many terms of this order in $`x_A`$, this behavior is not at all surprising. However, it can lead to large corrections to the factorized form for real nuclei (where $`x_A`$ is of order 0.04).

The $`x_1`$ dependence is clearly illustrated in Figure 5. Here, we see that the factorized form $`K_0`$ is indeed a good approximation for moderate $`x_1`$. In the region $`0.15<x_1<0.75`$, the ratio changes by approximately 15% for <sup>208</sup>Pb, while in the more restricted region $`0.2<x_1<0.7`$, the change is less than 8%. For smaller nuclei, the region of ‘moderate $`x_1`$’ is more restricted. The behavior of our curves for small $`x_1`$ is due to the fact that the denominator of our ratio diverges as $`x_10`$, while the numerator samples the distribution function and smears the divergence. For large $`x_1`$, the behavior is quite similar to that of the ordinary convolution model, Eq.(7), which is due to the fact that the denominator vanishes as $`x_11`$ while the numerator remains finite because the partons share momentum. This effect is usually referred to as Fermi motion. The dependence of the ratio on nuclear size is not too large (on the order of 15% for moderate $`x_1`$) in the range $`0.032x_A0.059`$ considered. The tremendous dependence of the placement of these curves on the value of $`x_2`$ can be attributed to the turnover of the gluon momentum distribution at small $`x`$. Since this turnover occurs at small $`x`$, we expect large variation only when $`x_2`$ is extremely small.

To see this explicitly, we plot the ratio $`K/K_0`$ as a function of $`x_2`$ in Figure 6. Here, we can clearly see that the large variations are confined to $`x_2<0.01`$. Changing the scale of the input distributions allows us to explore the behavior of our ratio as a function of $`Q^2`$. In particular, the turnover of the gluon momentum distribution at $`Q^2=5.0`$ GeV<sup>2</sup> occurs below our cut-off at $`x=10^4`$. Hence the increase observed at $`Q^2=2.5`$ GeV<sup>2</sup> is no longer present. Instead, we see a decrease induced by the increasing value of the gluon momentum distribution. The numerator is approximately constant in the region of small $`x_2`$ due to the smearing of the convolution. Since our ratio actually depends on the combination $`x_2+x_L/2`$ rather than $`x_2`$, the small-$`x_2`$ behavior is stabilized by taking $`x_L`$ finite. As displayed in Figure 7, our ratio saturates as $`x_2`$ is reduced when $`x_L`$ is finite. In addition, its behavior is far more moderate.

Our numerical analysis has shown that in a very general region the ‘naive’ expectation, Eq.(72), for our matrix element is quite a good approximation. As long as $`x_A`$ and $`x_L/x_A`$ are small enough, $`x_2+x_L/2`$ is not too small, and $`0.2<x_1<0.7`$, corrections run in the 10-15% level or less. We can get even better approximations by including an $`x_A`$-dependent constant of order $`0.8`$. However, we must be very careful when using this approximation generally. It is easy to see from the above plots that the ratio varies quite quickly as one leaves the region of validity.

## V Phenomenology

In previous studies , one has assumed a Gaussian spatial nuclear distribution $`\mathrm{exp}[(y^{}/2R_A^{})^2]`$, which leads to a phenomenological form for the two-parton correlation in nuclei,

$$T_{qg}^A(x,x_T,x_L)=\stackrel{~}{C}MR_AAQ(x)(1e^{x_L^2/x_A^2}).$$
(83)
Compared with the recent HERMES experimental data  on nuclear modification of quark fragmentation function in DIS, one has extracted  the constant $`\stackrel{~}{C}\alpha _s^20.00065`$ GeV<sup>2</sup>. The strong coupling constant $`\alpha _s`$ should be evaluated at a scale $`Q^22.8`$ GeV<sup>2</sup>. Such a value is also consistent with what is extracted from the nuclear transverse momentum broadening , $`\stackrel{~}{C}\alpha _\mathrm{s}^2=0.00086`$ GeV<sup>2</sup>.

Assuming now that $`K_0(x_1,x_2,x_L)`$ is a good approximation of $`K(x_1,x_2,x_L)`$, we can express the two-parton correlation in nuclei of Eq.(19) as

$`T_{qg}^A(x,x_T,x_L)={\displaystyle \frac{A^2}{32R_A^2}}R_A^2\mathrm{\Delta }_{}^2`$ $`[Q(x+x_L)G(x_T)+Q(x)G(x_T+x_L)`$ (85)
$`2Q(x+x_L/2)G(x_T+x_L)e^{\xi ^2}].`$
Here we assume a Gaussian distribution for the two-nucleon correlation function $`\stackrel{~}{r}_2`$. For moderate values of $`x>>x_L<x_A<<1`$, one can assume $`Q(x+x_L)Q(x)`$, the above expression can be approximated as

$`T_{qg}^A(x,x_T,x_L)`$ $`=`$ $`{\displaystyle \frac{A^2}{32R_A^2}}R_A^2\mathrm{\Delta }_{}^2Q(x)\left[G(x_T)+G(x_T+x_L)\right]`$ (86)
$`\times `$ $`\left[1{\displaystyle \frac{2G(x_T+x_L/2)}{G(x_T)+G(x_T+x_L)}}e^{x_L^2/4x_A^2}\right].`$ (87)
The above derived factorization form is very close to the phenomenological one in Eq.(83), especially for $`x_Lx_T`$. For large $`x_Lx_T`$, when $`G(x_T)G(x_T+x_L)`$, the coefficient in front of the exponential factor will have additional suppression as compared to the phenomenological model. However, for not so large $`Q^2`$ and $`x_L<x_A`$, $`G(x_T)G(x_T+x_L)`$. In this kinematic region, one can then relate the parameters in the phenomenological form to our result within the convolution model

$$\stackrel{~}{C}MR_A\frac{A}{4R_A^2}R_A^2\mathrm{\Delta }_{}^2G(x_T).$$
(88)
Here, we have reduced the nuclear radius by 1/2 in Eq. (87) in order to match the phenomenological form in Eq. (83). With the value of $`\stackrel{~}{C}\alpha _s^20.00065`$ from HERMES experiment and $`G(x_T)3`$ at $`x_T=x_Bk_T^2/Q^20.01`$ ($`x_B=0.124`$, $`Q^2=2.8`$ GeV<sup>2</sup>, $`k_T^20.25`$ GeV<sup>2</sup>), we have

$$R_A^2\mathrm{\Delta }_{}^21.65,$$
(89)
which is on the order of 1 and independent of the $`A`$ as expected. One can consider this as a qualitative agreement between our calculated two-parton correlation in nuclei and the experimental measurements.

## VI Conclusion and Discussion

In this paper, we have studied the generalized twist-four nuclear parton matrix elements $`K(x_1,x_2,x_L)`$, both direct and momentum-crossed ones, in the framework of a Fock hadronic state expansion of the nuclear states. These twist-four nuclear parton matrix elements determine the effect of multiple parton scattering in hard processes involving nuclei, e.g., the nuclear modification of the fragmentation functions. Assuming that the contributions of higher Fock states induced by nucleonic interaction are small, we have shown that the leading contribution to the twist-four nuclear parton matrix elements can be expressed as a convolution of twist-two nucleonic off-forward parton distributions and the two-nucleon correlation function inside a nucleus. In the limit of extremely large nuclei, $`A\mathrm{}`$, or $`x_A=1/MR_A0`$, we have also shown that the twist-four nuclear parton matrix elements can factorized into the product of twist-two nucleonic parton distributions. However, we demonstrated that the matrix elements $`K`$ are not analytic in $`x_A`$ at $`x_A=0`$ (Corrections around $`x_A=0`$ cannot be expanded as powers of $`x_A`$).

To verify the factorization approximation, we have evaluated the twist-four nuclear matrix elements numerically as the convolution of twist-two nucleonic parton distributions and two-nucleon correlation functions inside a nucleus, assumed to have a simple Gaussian form. For $`x_L<x_A`$ and moderate $`x_1`$, we found that the factorization is a good approximation within $``$20% for large nuclei. However, the deviations become very significant for small $`x_1`$ or $`x_2+x_L/2`$, large $`x_A`$, or $`x_L>x_A`$. The corrections at small $`x_2+x_L/2`$ are particularly large when the gluon momentum distribution $`G(x_2+x_L/2)`$ is large. Furthermore, for $`x_{1,2}x_L/2`$, one can no longer express the nuclear matrix elements as the convolution of twist-two nucleonic parton distributions. In this region, off-forward nucleonic parton distributions, which are not known experimentally, are needed for the numerical evaluation. Therefore, one could conceivably use the measured nuclear effects in this kinematic region to constrain the nucleonic OFPD’s.

Another important nuclear effect on the parton matrix elements that we have not considered so far is similar to the nuclear shadowing of the parton distributions or depletion of the effective parton distribution per nucleon. This effect is known experimentally to be large for small $`x_B`$, at least for the quark distributions . One can understand nuclear shadowing as the consequence of coherent initial scattering processes involving multiple nucleons  when one calculates the nuclear parton distributions in terms of nucleonic parton distributions. In the framework of our lowest Fock state expansion, these multiple scattering processes will involve direct multi-nucleonic correlations. Therefore, such multiple scattering effects in principle are higher-twist contributions and should be suppressed at very large $`Q^2`$. In our calculation of the twist-four nuclear parton matrix elements, such higher-twist contributions from multiple scattering processes involving more than two nucleons should also contribute, but they are suppressed at large $`Q^2`$.

###### Acknowledgements.

The work is supported in part by the Director, Office of Science, Office of High Energy and Nuclear Physics, and by the Office of Basic Energy Sciences, Division of Nuclear Sciences of the U.S. Department of Energy under Grant No. DE-AC03-76SF-00098 and in part by NSFC under project No. 19928511