The spectrum of dileptons produced in heavy ion collisions at high energies has been proposed in the past to provide information about the dynamical evolution of quark-gluon plasmas (see, e.g., the review by Ruuskannen in ref.[@dilepton]). Unlike hadronic probes, that spectrum is sensitive to the earliest moments in the evolution, when the energy densities are an order of magnitude above the QCD confinement scale (1 GeV/fm$^3$). However, the small production cross sections for pairs with invariant mass above 1 GeV together with the large combinatorial background from decaying hadrons necessitate elaborate procedures to uncover the signal from the noise. The PHENIX detector[@CDR], now under construction at the Relativistic Heavy ion Collider (RHIC), is designed to measure $ee,e\mu$, and $\mu\mu$ pairs to carry out this task. One of the important sources of background in the few GeV mass range arises from semileptonic decay of charmed hadrons ($D \bar D$). As shown recently in ref.[@vogt], the expected thermal signals in that mass range may only reach 10% of the number of pairs from open charm decay. Therefore, special kinematical cuts and precise $e\mu$ measurements will be needed to uncover the thermal signals[@CDR].

There have been several attempts[@pcm; @MW] to take advantage of the large open charm background as a probe of the evolving gluon density. Most mid-rapidity charm pairs are produced via gluon fusion[@ccbar]. Therefore, the dileptons from charm decay carry information about the distribution of primordial gluons before hadronization. In fact, the inside-outside cascade nature of such reactions greatly suppresses[@LM; @LMW] all sources of charm production except the initial perturbative QCD source. Therefore, the open charm background is dominated by the [*initial*]{} gluon fusion ($gg\rightarrow
c\bar{c}$) rate.

The initial $c\bar c$ rate depends sensitively on the nuclear gluon structure function[@vogt2], $g_A(x,Q^2)$. The quantity of fundamental interest[@mullerq; @eskolaq] is the gluon shadowing function $$R_{g/A}(x,Q^2)=g_A(x,Q^2)/Ag_N(x,Q^2).
\label{EQ:shad}$$ The point of this Letter, is to show that the $A$ dependence of continuum dilepton pairs in the few GeV mass region provides a novel probe of that unknown gluon structure. We also show that the required measurements of $ee$,$e\mu$, and $\mu\mu$ pair yields in $p+A\rightarrow \ell^+ \ell^- X$ at $\sqrt s=200$ AGeV are not only experimentally feasible at RHIC but also that the open charm signal can be easily extracted via the proposed PHENIX detector.

A key advantage of the continuum dilepton pairs from open charm decay over those from $J/\psi$ decay is that it is possible to test the applicability of the underlying QCD dynamics at a given fixed $\sqrt s$ by checking for a particular scaling property discussed below. In quarkonium production the mass is fixed and the required scaling can be checked only by varying the beam energy. In $p+A\rightarrow J/\psi$ production the required scaling was unfortunately found to be violated in the energy range $20< \sqrt{s} < 40$ AGeV[@moss] thus precluding a determination of $R_A$. Possible explanations for the breakdown of QCD scaling for $J/\psi$ production and its anomalous negative $x_F$ behavior[@vogt2; @moss] include interaction of next-to-leading order quarkonium Fock state with nuclear matter[@ks], nuclear and comover $J/\psi$ dissociation[@huf], and parton energy loss mechanisms[@gavin]. It remains an open question whether the required scaling will set in at RHIC energies $60 < \sqrt s < 200$ AGeV.

Shadowing of the quark nuclear structure functions, $q_A(x,Q^2)$, is well established from deep inelastic $\ell A\rightarrow \ell X$ reactions[@emc]. For heavy nuclei, the quark structure functions are shadowed by a factor, $R_{q/A}(x\ll 0.1,Q^2)\approx 1.1 - 0.1 \;A^{1/3}$, nearly independent of $Q^2$. Perturbative QCD analysis[@mullerq; @eskolaq] predicts that the gluon structure will also be shadowed at $x<0.01$ due to gluon recombination processes. Therefore, a systematic measurement of gluon shadowing is of fundamental interest in understanding the parton structure of nuclei. The gluon nuclear structure is also of central importance in the field of nuclear collision since it controls the rate of mini-jet production that determines the total entropy produced at RHIC and higher energies[@hijingprl]. We focus here on $p+A$ collision rather than $A+A$ to minimize the combinatorial $\pi,K$ decay backgrounds and other final state interaction effects.

To demonstrate the sensitivity of open charm dileptons to gluon shadowing $R_{g/A}(x,Q^2)$, we calculate the pair spectrum from open charm decay in $p+Au$ reactions at 200 GeV/A using for illustration two different gluon shadowing functions. The first assumes that gluon shadowing is identical to the measured quark shadowing, i.e. $R_{g/A}(x,Q^2)=R_{q/A}(x,Q^2=4\;{\rm
GeV}^2)$, and is thus essentially $Q^2$ independent. This is the default assumption incorporated in the HIJING Monte Carlo model[@Hijing]. The second is taken from Eskola[@Eskola], where $R_{g/A}(x,Q^2)$ is computed using the modified GLAP evolution[@mullerq] starting from an assumed $R_{g/A}$ consistent with the measured quark shadowing function at $Q^2=4$ GeV$^2$. This second shadowing function depends on scale and differs for $q$ and $g$. The two shadowing functions are shown in Fig. 1.

The initial charm pair distribution from $B+A$ nuclear collisions at impact parameter $\vec{b}$ is computed as in [@MW]: $$\begin{aligned}
&&\frac {dN^{BA}(\vec{b})}{dp_\perp^2 dy_1 dy_2}
=K \int d^2r_b d^2r_a \delta ^{(2)} (\vec b-\vec r_b -\vec r_a) \nonumber \\
&   &\sum _{b,a} x_b \Gamma_{b/B}(x_b, Q^2, \vec r_b) x_a \Gamma_{a/A}(x_a,
Q^2, \vec r_a) \frac {d \hat \sigma_{ab}}{d \hat t} 
\label{EQ:dncc} \end{aligned}$$ where $\Gamma_{a/A}(x, Q^2,\vec r)=T_A(\vec r) f_{a/N}(x,Q^2) R_{a/A}(x,Q^2,
\vec r)$ is the nuclear parton density function in terms of the known nucleon parton structure functions $f_{a/N}(x,Q^2)$, the nuclear thickness function $T_A(\vec r) = \int dz n_A(\sqrt {z^2+\vec r ^2})$, and the unknown impact parameter dependent shadowing function $R_{a/A}(x,Q^2, \vec r)$. The conventional kinematic variables $x_{a,b}$ are the incoming parton light cone momentum fractions, $\hat t=-(p_b-p_1)^2$, and the final $c,\bar{c}$ have rapidities $y_1,y_2$ and transverse momenta $\vec{p}_{1\perp}=-\vec{p}_{2\perp}
\equiv \vec{p}_\perp$. Next-to-leading order corrections to the lowest order parton cross sections ${d \hat \sigma_{ab}}/{d \hat t}$ are approximately taken into account by a constant $K$ factor[@SV; @vogt3]. We use the recent MRSA[@MRSA] parton structure functions. From a fit to low energy open charm production data in $pp$ collisions, we fix $m_c=1.4$ GeV, $K=3$, $Q^2 =
\hat s/2$ for $g g \rightarrow c\bar c $ and $\hat s$ for $q\bar q \rightarrow
c \bar c$. This choice of parameters leads to a charm pair cross section of $340 \mu b$ for $pp$ collisions at RHIC. For impact parameter averaged collisions, the integral over the transverse vector in eq.(\[EQ:dncc\]) leads to a factor: $BA R_{b/B}(x_b,Q^2) R_{a/A}(x_a,Q^2)/\sigma ^{BA}_{in}$, where $\sigma ^{BA}_{in}$ is the inelastic $B+A$ cross section and $R_{a/A}(x_a,Q^2)$ is the impact parameter averaged shadowing function.

In order to compute the $D\bar D$ pair distribution, a hadronization scheme for $c \rightarrow D$ must be adopted. In $e^+e^-\rightarrow c\bar{c}\rightarrow
D\bar{D}X$ hadronization can be modelled via string fragmentation or fitted for example via the Peterson fragmentation function (see [@vogt4]). The final $D$ carries typically only a fraction $\sim 0.7$ of the original $c$ momentum. However, in $pp$ collisions charm hadronization is complicated by the high density of partons produced during beam jet fragmentation. In this system recombination or coalescence of the heavy quark with comoving partons provides another mechanism which in fact seems to be dominant at least at present energies. The inclusive $pp \rightarrow DX$ data is best reproduced[@vogt; @e769] with the delta function fragmentation $D(z)=\delta
(1-z)$ in all observed $x_f$ and moderate $p_T$ regions. This observation can be understood in terms of a coalescence model if the coalescence radius is $P_c \sim 400$ MeV. We assume that hard fragmentation continues to be the dominant mechanism at RHIC energies, and hence no additional A dependent effects due to hadronization arise. However, this assumption must be tested experimentally. Our main dynamical assumption therefore is that the $pp\rightarrow D$ transverse momentum distributions can be accurately reproduced from the QCD level rates assuming hard fragmentation as at present energies[@e769]. At RHIC this can be checked either via single inclusive leptons[@CDR] or directly via $K\pi$[@star].

With this assumption, the impact parameter averaged $D\bar D$ pair distribution in $p+A$ is given by $$\begin{aligned}
&&\frac {dN^{pA}}{dp_\perp^2 dy_3 dy_4}
= \frac {KA}{\sigma ^{pA}_{in}} \frac {E_3 E_4}{E_1 E_2} \nonumber \\
&&\sum _{b,a} x_b f_{b/N}(x_b, Q^2)  x_a
f_{a/N}(x_a, Q^2) R_{a/A}(x_a,Q^2)\frac {d \hat \sigma_{ab}}{d \hat t} \end{aligned}$$ From this distribution of charmed mesons, we computed the dilepton spectrum via Monte Carlo using JETSET7.3[@jetset] to decay the $D\bar{D}$.

While the lepton pair mass $M$ and rapidity $y$ fluctuate around a mean value for any fixed $M_{c\bar{c}}$ and $y_{c\bar{c}}$, on the average $\langle y
\rangle \approx \langle y_{c\bar{c}} \rangle $, and $\langle M_{c\bar{c}}
\rangle$ can be well approximated by a linear relation $ \langle M_{c\bar{c}}
\rangle \approx \beta M+ M_0$ over the lepton pair mass range $1 \le M \le
4$GeV. For the delta function fragmentation case, $\beta\approx 1.5$ and $M_0
\approx 3.0$ GeV are determined by the D-meson decay kinematics.

Since gluon fusion dominates, there is an approximate scaling in terms of the light cone fraction $x_A$ that one of the gluon carries from nucleus $A$: $$\begin{aligned}
\ln x_A&=& -y_{c\bar{c}}+ \ln( M_{c\bar{c}}/\sqrt s) \nonumber \\
&\approx& -y+ \ln [(\beta M + M_0) /\sqrt s] \label{EQ:shift}\end{aligned}$$ The ratio of the dilepton $dN/dMdy$ spectra in $pA$ to scaled $pp$ for different pair masses is thus expected to scale approximately as $$\begin{aligned}
&&R^{pA}_{e\mu} ( M,y=-\ln x_A + \ln [ (\beta M+M_0)/
\sqrt s ] \;) \nonumber \\ 
&&\equiv \frac {1} {\nu} \frac {dN^{pA}_{e\mu}} {dN^{pp}_{e\mu}}
\approx R_{g/A}(x_A,Q^2\!\approx\! (\beta M\!+\!M_0)^2\!/2\;)
\label{EQ:scaling}\end{aligned}$$ where $\nu \equiv A {\sigma ^{pp}_{in}}/{\sigma ^{pA}_{in}} \sim A^ {1/3}$. In order to test gluon dominance and the accuracy of the above approximate scaling we compare in Fig. 1 the gluon shadowing function to the above dilepton ratio from the Monte Carlo calculation. The solid curves are the input gluon shadowings as a function of the Bjorken $x$. The other six curves are ratios of dilepton $dN/dMdy$ spectra of shadowed $p+Au$ over those from unshadowed $p+Au$ as given by eq.(\[EQ:scaling\]) as a function of the scaling variable $x_A$ for three different dilepton masses.

Overall, the scaled ratio of lepton pair spectra approximates the shadowing function remarkably well. We conclude that this ratio can therefore serve to map out gluon shadowing in nuclei. Note that in the case of Eskola shadowing, even the $Q^2$ dependence of the shadowing function is visible through the rise of the ratio curves with $M$ in the small $x$ region.

Finally, it is important to estimate the dilepton background to see if the proposed signal is experimentally feasible. Therefore we also calculated the signal and background dileptons for the PHENIX detector[@CDR] taking into account the detector geometry and specific kinematical cuts. The electron background is mainly due to $\pi^0$ Dalitz and photon conversions. The background muon arises from random decays of charged pions and kaons. The electrons from $\pi^0$ Dalitz decay can be suppressed by small angle cut on dielectrons[@Zhang], and the background muons are mainly suppressed by reducing the free decay volume. For the detector geometry, the electron arms cover pseudo-rapidity range $-0.35 < \eta _e < 0.35$, and azimuthal angle range $\pm (22.5^\circ ,112.5^\circ )$. The muon arm covers pseudo-rapidity range $1.15 < \eta _\mu < 2.44$ and almost full azimuthal angle. For the kinematical cut, we take $E_e > 1$ GeV, $E_\mu > 2$ GeV, and require that the relative azimuthal angle of the lepton pair $\phi _{\ell^+\ell^-} > 90^\circ $ in order to improve the signal-to-background ratio.

For this background study, we use central $p+Au$ collision with HIJING shadowing. We calculate the $ee$, $e\mu$, and $\mu\mu$ signal and backgrounds coming into the detector, where the background electrons and muons are generated from HIJING Monte Carlo calculation. We find (see Fig. 2) that the signal-to-background ratio for $ee$ is very large, thus the dielectron signal from open charm decay is the easiest to extract. That ratio falls to about 2.5 for $e\mu$, and about 1/4 for $\mu\mu$. The pair rapidity distributions of the signal and backgrounds for central $p+Au$ collision are shown in Fig. 2. Due to the detector geometry, $ee$, $e\mu$, and $\mu\mu$ spectra cover pair rapidity regions centered at about 0, 1, and 2 (One can also study $Au+p$ collision and measure pair rapidity regions centered at $-1$ and $-2$). Like-sign subtraction[@NA38] should significantly reduce the noise, especially in the $\mu\mu$ channel. We conclude that the proposed measurement is feasible.

Further details of the calculations will be reported elsewhere[@long]. For example, the Cronin effect and energy loss are estimated to lead to distortions up to $10\%$ at RHIC energies. From studies of the nuclear dependence of transverse momentum of $J/\psi$ production, the typical increment, $\delta
p_\perp^2(A)$, due to multiple collisions is limited to $\sim 0.34$ GeV$^2$ even for the heaviest nuclei[@gavinmg]. This momentum spread is shared between the $c$ and $\bar{c}$, and distributed approximately as $g(\delta
p_\perp)=e^{-\delta p_{\perp}^2/\Delta ^2}/\pi \Delta ^2$, $\Delta^2 \sim 0.17$ GeV$^2$. For energy loss, we assume that the charm quark loses an energy $\delta E$ in the lab frame, and the charm quark $p_\perp$ is reduced by $(1-\epsilon)$, where $\epsilon=\delta E/[m_\perp \sinh (y+y_0)]$. Combining these two effects, we may write for the charm quark that the final $\vec
p_\perp=(\vec {p_{\perp}^{ini}} + \vec \delta p_\perp) (1-\epsilon)$. Therefore the D-meson spectrum $F(p_\perp) \equiv d^2 N/d m_{\perp}^2$ becomes $F^{\prime}(p_\perp)=\int F({p_{\perp}^{ini}}) g(\delta p_\perp) d\vec \delta
p_\perp$. Taking $F(p_\perp) \propto e^{-\alpha p_\perp}$, $\alpha \simeq 1.3$ GeV$^{-1}$ from HIJING, and expanding the convolution to lowest order, the relative change of the D-meson spectrum is given by $$\begin{aligned}
F^{\prime}(p_\perp)/F(p_\perp) \simeq 1-\alpha p_\perp \epsilon +\alpha^2
\Delta^2 /4 \nonumber\end{aligned}$$ This is also the relative change of the $c\bar c$ pair spectrum when $M_{c\bar
c}=2 m_\perp$ in case of $y_1=y_2$, $p_{\perp 1}=p_{\perp 2}$. For $m_\perp
\sim 3$GeV, $\delta E=10$GeV, $\cosh y_0 \simeq 100$, the relative change in the pair spectrum is estimated to be $F^{\prime}(p_\perp)/F(p_\perp) \simeq
1\;-0.1( e^{-y}- 1)$.

In summary, we calculated the lepton pair spectra from open charm decay in two different shadowing scenarios. By scaling the ratios for different mass ranges according to eq.(\[EQ:scaling\]), we showed that the dilepton rapidity dependence of those ratios on $x_A$ reproduces well the underlying gluon shadowing function defined in eq.(\[EQ:shad\]). Finally we showed that the measurements required to extract the gluon shadowing are experimentally feasible at RHIC. The above analysis depends on the validity of our basic assumption regarding the hard fragmentation of charm quarks, and this can be checked explicitly via the single inclusive measurements of $D$ production.

We conclude by emphasizing the importance of determining gluon shadowing in $p+A$ to fix theoretically the initial conditions in $A+A$. In $A+A$ the open charm decay is regarded as an annoying large background that must be subtracted to uncover the thermal signal. In $p+A$ that background becomes the signal needed to determine the incident gluon flux in $A+A$. The continuum charm dileptons in $p+A$ at RHIC are likely to provide a unique source of information on the low $x_A$ nuclear gluon structure at least until HERA is capable of accelerating heavy nuclei.

Acknowledgements: We thank M. Asakawa, R. Bellweid, S. Gavin, P.E. Karchin, P. McGaughey, S. Nagamiya, M. Tannenbaum, J. Thomas, R. Vogt, X.-N. Wang, G. Young and W. Zajc for useful discussions. This work was supported by the Director, Office of Energy Research, Division of Nuclear Physics of the Office of High Energy and Nuclear Physics of the U.S. Department of Energy under Contract No. DE-FG02-93ER40764.

[20]{} P.V. Ruuskannen, [*Quark-Gluon Plasma*]{}, ed. R.W. Hwa, (World Scientific, Singapore, 1990), Adv. Ser. on Directions in High En. Phys., Vol. 6, p519.

S. Nagamiya [*et al.*]{}, in PHENIX Conceptual Design Report (1993) (unpublished); on internet see [*http://rsgi01.rhic.bnl.gov/ phenix/phenix\_home.html*]{}. Recently a second muon arm is added.

R. Vogt, B. V. Jacak, P. L. McGaughey and P. V. Ruuskanen, Phys. Rev. D [**49**]{}, 3345 (1994).

K. Geiger, Phys. Rev. D [**48**]{}, 4129 (1993).

B. Müller and X. N. Wang, Phys. Rev. Lett. [**68**]{}, 2437(1992).

S. J. Brodsky, P. Hoyer, A. H. Mueller and W. K. Tang, Nucl. Phys. B [**369**]{}, 519 (1992).

Z. Lin and M. Gyulassy, Phys. Rev. C [**51**]{}, 2177 (1995); ERRATUM-ibid. C [**52**]{}, 440 (1995).

P. Lévai, B. Müller and X.-N. Wang, Phys. Rev. C [**51**]{}, 3326 (1995).

S. Liuti, R. Vogt, Phys. Rev. C [**51**]{}, 2244 (1995); S. Gupta and H. Satz, Z. Phys. C [**55**]{}, 391 (1992). R.V. Gavai, R. M. Godbole, Report No. TIFR/TH/93-57 (unpublished).

A. Mueller and J. Qiu, Nucl. Phys. B [**268**]{}, 427 (1986); J. Collins and J. Kwiecinski, Nucl. Phys. B [**335**]{}, 89 (1990).

K. J. Eskola, J. Qiu, and X.-N. Wang, Phys. Rev. Lett. [**72**]{}, 36 (1994).

D. M. Alde [*et al.*]{}, Phys. Rev. Lett. [**66**]{}, 133 (1990). J. Badier [*et al.*]{}, Z. Phys. C [**20**]{}, 101 (1983).

D. Kharzeev and H. Satz, Report No. CERN-TH/95-214.

C. Gerschel and J. Hüfner, Phys. Lett. B [**207**]{}, 253 (1988).

S. Gavin and J. Milana, Phys. Rev. Lett. [**68**]{}, 1834 (1992).

EM Collaboration, M. Arneodo [*et al.*]{}, Nucl. Phys. B [**333**]{}, 1 (1990); NM Collaboration, P. Amaudruz [*et al.*]{}, Z. Phys. C [**51**]{}, 387 (1991); E665 Collaboration, M. R. Adams [*et al.*]{}, Phys. Rev. Lett. [**68**]{}, 3266 (1992).

X.-N. Wang, M. Gyulassy, Phys. Rev. Lett. [**68**]{}, 1480 (1992).

X.-N. Wang and M. Gyulassy, Phys. Rev. D [**44**]{}, 3501 (1991).

K. J. Eskola, Nucl. Phys. B [**400**]{}, 240 (1993).

I. Sarcevic and P. Valerio, Phys. Lett. B [**338**]{}, 426 (1994).

R. Vogt, Report No. LBL-37105 (1995).

A. D. Martin, R. G. Roberts and W. J. Stirling, Univ. of Durham, Report No. DTP/94/78 (1994).

R. Vogt, S. J. Brodsky and P. Hoyer, Nucl. Phys. B [**383**]{}, 643 (1992); R. Vogt and S. J. Brodsky, Nucl. Phys. B [**438**]{}, 261 (1995).

S. Frixione [*et al.*]{}, Nucl. Phys. B [**431**]{}, 453 (1994); P. E. Karchin (E769 Collaboration), Report No. FERMILAB-CONF-95-053-E.

STAR Collaboration, Nucl. Phys. A [**566**]{}, 277c (1994).

T. Sjöstrand, Report No. CERN-TH.6488/92 (1992), on intenet see [*http://heplibw3.slac.stanford.edu/FIND/FREEHEP/NAME/JETSET/FULL.*]{}

P. H. Zhang [*et al.*]{}, in The TALES/SPARHC EXPERIMENT at RHIC, p23.

NA38 Collaboration, Phys. Lett. B [**220**]{}, 471 (1989).

Z. Lin and M. Gyulassy, in preparation.

S. Gavin and M. Gyulassy, Phys. Lett. B [**214**]{}, 241 (1988).
