Nuclear gluons at RHIC in a multi-observable approach

We explore the possibility of measuring nuclear gluon distributions at the Relativistic Heavy-Ion Collider (RHIC) with $\sqrt{s}=200 \, {\rm GeV}$ proton-nucleus collisions. In addition to measurements at central rapidity, we consider also observables at forward rapidity, consistent with proposed upgrades to the experimental capabilities of STAR and sPHENIX. The processes we consider consist of Drell-Yan dilepton, dijet, and direct photon-jet production. The Drell-Yan process is found to be an efficient probe of gluons at small momentum fractions. In order to fully utilize the potential of Drell-Yan measurements we demonstrate how the overall normalization uncertainty present in the experimental data can be fixed using other experimental observables. An asset of the RHIC collider is its flexibility to run with different ion beams, and we outline how this ability could be taken advantage of to measure the $A$ dependence of gluon distributions for which the current constraints are scarce.


I. INTRODUCTION
Good control over the partonic structure of protons and heavier nuclei has become an indispensable ingredient in modern particle, heavy-ion, and astro-particle physics. For processes involving a large momentum transfer, Q Λ QCD ∼ 200 MeV, the nucleon's relevant degrees of freedom can be described by parton distribution functions (PDFs). Despite the progress in theoretical first-principle methods [1], the PDFs are still most reliably determined through a statistical analysis of a global set of experimental data. Along with the precise data from the Large Hadron Collider (LHC), the list of data types that are included in state-of-the-art PDF analyses has grown, now ranging from traditional inclusive deeply-inelastic scattering to jet, top-quark and heavy gauge-boson production [2,3]. At this moment, global fits of proton PDFs do not use any data from the Relativistic Heavy-Ion Collider (RHIC), and nuclear-PDF fits [4,5] use only inclusive pion data from RHIC [6,7]. The advantage of the lower center-of-mass (c.m.) energies of RHIC, √ s = 200 and 500 GeV, is that the underlying event is not as large as it is at the LHC, and thus e.g. jets can be better resolved at lower transverse momenta (p T ) [8,9]. These jet measurements are compatible with next-to-leading-order (NLO) perturbative QCD calculations [9][10][11] down to p T ∼ 10 GeV (which is the minimum p T of the measurements), so nothing really forbids using them in PDF analysis. Similarly, low-mass Drell-Yan (DY) events can be better resolved from the decays * ilkka.m.helenius@jyu.fi † lajoie@iastate.edu ‡ jdosbo@umich.edu § petja.paakkinen@jyu.fi ¶ hannu.paukkunen@jyu.fi of heavy flavour. In p+p collisions at the higher c.m. energy of √ s = 500 GeV, measurements of W ± bosons also become feasible [12]. These measurements provide complementary constraints on the fixed-target measurements [13] of the u/d ratio.
The current status of the global determination of nuclear PDFs has been recently reviewed e.g. in Refs. [4,5], and the field is rapidly developing. This is mainly driven by the p+Pb measurements at the LHC, but is also motivated by theoretical advances in upgrading global analyses to the next-to-NLO (NNLO) QCD level. Currently, the most recent global NLO fits are EPPS16 [14], nCTEQ15 [15] and DSSZ12 [16]. Out of these, EPPS16 has the widest data coverage and is currently the only one to use LHC measurements. In the current NNLOlevel fits [17,18], the experimental input is restricted to fixed-target data only. In this paper, we report our studies on the future prospects for constraining nuclear PDFs at RHIC, particularly with measurements at central and forward rapidities where forward acceptance corresponds to that proposed for the STAR [19] and sPHENIX experiments [20]. Projections of forward direct photon and Drell-Yan measurements at RHIC on nuclear PDFs have been separately considered e.g. in Ref. [19]. Here, we aim for a more systematic approach by combining multiple observables into a simultaneous analysis and more carefully assessing the experimental normalization uncertainty. We will base our study mainly on the EPPS16 analysis.

II. EXPERIMENTAL DATA PROJECTIONS
Several processes are expected to have an impact on nuclear PDFs at RHIC c.m. energies. Here, we will focus on such double differential measurements which, to leading order, probe PDFs at fixed momentum fractions.

arXiv:1904.09921v1 [hep-ph] 22 Apr 2019
In particular, we construct pseudodata projections for the Drell-Yan di-lepton, dijet, and direct photon-jet processes, differential in the invariant mass M and rapidity y of the produced pair. While the Drell-Yan process (on fixed target) has been used as a constraint for nuclear PDFs already in the pioneering EKS98 analysis [21], the use of dijets [22,23] has been realized only in the recent EPPS16 fit [14]. Currently, there are no available photonjet data to be included in the global analyses though the potential of the process has been discussed [24]. In principle, these processes individually constrain different quark-gluon combinations of PDFs, and can also be used together to limit the effect of normalization uncertainties as will be described later.
To generate our projections, we first impose fiducial acceptance requirements on a barrel detector with forward instrumentation. Uncertainty projections were generated for a barrel covering the full azimuth of 0 < φ < 2π and pseudorapidity acceptance of |η| < 1, where the detector is assumed to have full tracking as well as electromagnetic and hadronic calorimetry such that jets can be robustly measured. In conjunction with the barrel central rapidity detector, a forward spectrometer, incorporating tracking and electromagnetic and hadronic calorimetry with pseudorapidity acceptance of 1.4 < η < 4 and full azimuthal coverage, is also considered.
Projections were determined by taking the cross sections as predicted in the Pythia 6 event generator [25,26] and multiplying them by total integrated luminosity projections at RHIC. Assumed luminosities were 197 pb −1 for p+p collisions and 0.33 pb −1 for p+Au collisions, corresponding to the anticipated sPHENIX run plan for the second and third years of operation in the early 2020's. Estimates of experimental efficiencies are also applied, as described below for each process. The total expected yields were converted to per event yields by dividing by the total p+p cross section times the expected luminosity. Thus, the ratio of the p+Au to p+p yields is always unity and the statistical uncertainties of the ratio are indicative of the actual statistical uncertainties on a measurement of R pA , where R pA is defined as Since the detector has both central and forward instrumentation, there are a number of rapidity regions that each observable can probe. Ideally measurements should be made in each region, as different values of x will be probed in both the proton and the nucleus when measuring at central and/or forward rapidities. Thus, we consider several combinations of observables, for which a summary table is shown in Tab. I.
Drell-Yan data are generated in both the central barrel and the forward arm independently from one another. In each case, the Drell-Yan dilepton pair is measured in the invariant mass range of 4.5 < M + − < 9 GeV/c 2 . Experimental efficiencies for the reconstruction of Drell-Yan pairs were estimated from a full Geant4 [27] simulation of the sPHENIX detector (including forward instrumentation). Tight cuts on the simulated data were used to reduce backgrounds so that the simulated measurement is dominated by Drell-Yan pairs and backgrounds from decays, conversions, etc. were minimal. Overall pair efficiencies vary between 10-15% over the invariant mass range considered for both central and forward measurements. The dijet data are considered in both regions to probe a variety of x values. Jets were determined from final-state Pythia particles with the anti-k T algorithm with R = 0.4 [28]. Two dijet pseudodata samples are constructed, one in which both jets are measured in the central barrel and another where one jet is measured at central rapidity and the other is measured at forward rapidity. Since the p T reach of jets becomes smaller at forward pseudorapidities, the leading jet at central rapidity is required to have p T > 12 GeV/c and the subleading jet at forward rapidity is required to have p T > 8 GeV/c. Experimental efficiencies for the reconstruction of dijet pairs were estimated in a similar fashion as the Drell-Yan data. Jets were reconstructed with both hadronic and electromagnetic calorimeter deposits and efficiencies were found to be approximately 80%, where these efficiencies become smaller towards the edge of the detector acceptance when some fraction of the jet cone lies outside of the detector. At forward pseudorapidity, the efficiencies were generally smaller, varying between 40-70% depending on the pseudorapidity of the jet. The direct photonjet channel is expected to have high impact on the gluon nuclear PDF at small x when the process is measured in the forward direction. However, the photon-jet channel is difficult to measure at forward pseudorapidities due to large backgrounds from π 0 → γγ decays. Thus, we only generate photon-jet pseudodata where both are measured in the central barrel, where previous direct photon measurements at RHIC have been made and where future RHIC experimental upgrades are expected to be able to measure this process. We note that if new forward instrumentation at RHIC were available to separate direct photons from backgrounds at forward pseudorapidities this would add a powerful additional observable to constrain the nuclear gluon PDF at low x, complementary to Drell-Yan [19]. The photon-jet cross sections were generated with p γ T > 10 GeV/c and p jet T > 8 GeV/c. Photon-jet reconstruction efficiencies were evaluated similarly to the dijet and Drell-Yan data, where the efficiency was found to be approximately 70% integrated across the central rapidity of the barrel detector.

A. Generation of pseudodata
From the Pythia simulations for p+p and p+Au collisions we keep the relative statistical error, but construct the pseudodata points for the expected nuclear modifica- where δ uncorr. signifies the total uncorrelated data uncertainty and r is a Gaussian random variable. To obtain δ uncorr. , we add in quadrature the statistical uncertainty in the anticipated yield in p+p and p+Au collisions. A 4% normalization uncertainty is assumed to account for the model dependence in determining N coll used in determining the R pA ratio. The overall scale of this uncertainty is unimportant, however, assuming it is common to all measurements, as we will detail later. In addition, for dijet (photon-jet) measurement, another 5% (4%) uncorrelated bin-to-bin systematic uncertainty is added, corresponding to the residual experimental systematic error that does not cancel in the ratio. For the Drell-Yan case the statistical uncertainty dominates and no additional systematic uncertainty is added. A systematic uncertainty of the order of 5% is clearly smaller than what one can expect to be present in measurements for the absolute cross sections. However, if the p+p and p+Au runs are made soon one after the other (so that the detector configuration and calibration remains unaltered), much of the systematic uncertainty can be expected to cancel. We note that recent dijet measurements by the CMS collaboration [29] quote a systematic uncertainty even less than 5%.
The central values for R pA in Eq. (2) were obtained by NLO-level calculations using the CT14NLO [30] freeproton PDFs and EPPS16 [14] nuclear modifications. For dijets we used meks (v1.0) [31], with the anti-k T algorithm, taking a jet cone R = 0.4, and fixing the QCD scales to the average of the two highest-p T jets. The leading jet was required to have p T > 12 GeV/c and the subleading jet p T > 8 GeV/c. These unequal cuts are necessary to avoid sensitivity to soft-gluon resummation. For photon-jet we used jetphox (v1.3.1) [32,33] where the jet was defined by a k T algorithm with R = 0.4 and the QCD scales were fixed to p T of the photon. No isolation criteria for the photons were imposed. The NLO Drell-Yan cross sections are standard, and were calculated with a private code based on Ref. [34], fixing the QCD scales to the invariant mass of the dilepton pair. Besides the photon-jet process for which there are no data in the EPPS16 analysis, the scale choices are the same as made in the EPPS16 fit.

A. The Hessian reweighting technique in a nutshell
We estimate the impact of the projected data on the EPPS16 nuclear PDFs by the PDF reweighting (also called PDF profiling) method [23,[35][36][37][38]. In this method, one studies the function where the first term describes the behaviour of the original global χ 2 in the EPPS16 analysis, and the second term is the contribution of the new data to the overall χ 2 budget. The central fit of EPPS16 corresponds to z = 0, and the error sets S ± i are defined in the z space by and they are known to increase the original χ 2 function by T = 52 units. The values for δ ± i are given in the EPPS16 paper [14], from which a i and b i coefficients in Eq. (3) can be solved. The contribution from the "new" data is defined as where D i and E i denote the ith data point and it's error. The overall normalization uncertainty is marked by E norm. . The theoretical prediction T i ( z) we write as where the coefficient for β ik and γ ik can be obtained by computing the predictions for T i with all the PDF error sets. As a result, the total χ 2 in Eq. (3) becomes an analytic function of z which we numerically minimize with respect to z and f N . We note that to avoid D'Agostini bias [39], the normalization factor in Eq. (5) f N multiplies the theoretical prediction T i , and not the data value D i . Since the pseudodata are based on EPPS16, the new minimum is, by construction, always very close to z = 0. After finding the parameters z min that correspond to the minimum of Eq. (3), we expand where H is the second-derivative (Hessian) matrix. By diagonalizing the matrix H this becomes where v = P ( z − z min ), where P is the orthogonal matrix that diagonalizes H, i.e. P T HP = 1. The new error sets are then defined as in Eq. (4) assuming that the original tolerance is not altered, i.e. that each new error setŜ ± i still correspond to ∆χ 2 (Ŝ ± i ) = 52.

B. Correlating the overall normalization
The normalization uncertainty in Eq. (5) we discuss here is that of the luminosity determination of the minimum-bias data sample. At the LHC, the p+A luminosities are determined by Van der Meer scans [40]. Alternatively, the measured per-event yields dN pA /N events are converted to cross sections dσ pA by where the average number of binary nucleon-nucleon collisions N coll is estimated from a Glauber-type model [41]. This leads to a model-dependent normalization uncertainty which is difficult to determine. Furthermore, the inelastic proton-nucleon cross section σ inelastic pn appearing in Eq. (9) is very sensitive to the physics at low scales Q 2 ∼ Λ 2 QCD , and it is presumably lower than the values measured in proton-proton collisions due to shadowing/saturation effects. The overall normalization is thus problematic in this approach. However, the normalization issue can be overcome by simultaneously measuring several observables from the same minimum-bias data sample. The reason is that there is only one single normalization uncertainty and in Eq. (5) the index i runs through all data points, not just those belonging to a one single observable. Including data that probe PDFs in a relatively better constrained region thus serves to "calibrate" the overall normalization.
To demonstrate how this works we have performed a PDF-profiling analysis first using only the forward Drell-Yan data, and then supplementing these data with the central-barrel dijet data. When only the Drell-Yan data are used, the constraints appear very weak. This is shown in Fig. 1 where the original EPPS16 uncertainties on the predictions are barely affected by the reweighting. The reason for the inability of these data to provide constraints is that the nuclear modification is predicted to be rather flat at small x and the variations in PDFs around the original central value can be compensated by suitably tuning the normalization f N . The flatness of the predicted Drell-Yan nuclear modification originates, to some extent, from the fit functions used in the EPPS16 analysis, but also from the scale evolution of PDFs which tends to flatten out the nuclear modifications in sea-quarks. Here, we took the normalization uncertainty to be 4%, but if a larger number would have been used (e.g. 10%) even fewer constraints would have been obtained.
The situation changes when the central-barrel dijet projections are also included. These data probe the nuclear PDFs at much higher x than the Drell-Yan data and carry significant sensitivity also to the rather-well constrained sum of valence quarks u A valence +d A valence . The nuclear modifications for the dijets are expected to exhibit some excess (antishadowing) around y dijet ∼ 1 which turns into a suppression (EMC effect) for y dijet ∼ −1. Such a pattern can not be mimicked by the overall normalization and leaves thus less room for f N variation. Since the normalization is now common for the dijet and Drell-Yan data, the Drell-Yan data have much larger impact. This is shown in Fig. 2, which should be compared   to Fig. 1. While the uncertainties for the central-barrel dijet data are only slightly reduced from their original EPPS16 values, the inclusion of these data is crucial in fixing the overall normalization. We have also observed that our results do not significantly depend on the exact value we pick for the normalization uncertainty.

C. Simultaneous analysis of Drell-Yan, dijet and photon-jet pseudodata
Following the observation made in the previous subsection, our strategy is to simultaneously analyze several observables that share the common normalization uncertainty. To separate the effect of forward-arm measurements, we first present the results using only the centralbarrel data, and then include the data simulated with the forward-arm acceptance. In Fig. 3 we summarize the Drell-Yan, dijet and photon-jet pseudodata within the central-barrel acceptance −1 < η < 1. The light-blue bands ("EPPS16") show the original EPPS16 predictions, and the darker bands ("EPPS16+CB") the error bands obtained after the reweighting analysis. We observe a modest improvement in the uncertainty bands for dijet and γ-jet cases. The precision of the Drell-Yan measurements is not expected to be high enough to set constraints as the sea quarks at 10 −2 < x < 10 −1 are already rather well constrained by the fixed-target DIS data. In Fig. 4, in turn, we show the combined pseudodata within the full centralbarrel and forward-instrumentation acceptance [42]. The light blue bands are again the original EPPS16 predictions, and the green bands ("EPPS16+CB+FI") are the uncertainties obtained in the reweighting exercise. The reduction in the PDF uncertainties is now more significant than in the central-barrel-only case shown in Fig. 3.
The impact of both "EPPS16+CB" and "EPPS16+CB+FI" analyses on EPPS16 is shown in Fig. 5 where we plot the average sea-quark modification for Au, , together with the gluon nuclear modification, Here, f p/A i (x, Q 2 ) denotes the parton density in a bound proton and f p i (x, Q 2 ) is the free-proton PDF. We omit here the valence quarks as we found no effects there. The improvement we find in R A u+d+s is rather weak in both cases. In the central-barrel analysis, there is a modest improvement in the gluons across all values of x, though the small-x improvement is merely a consequence of better constrained anti-shadowing regime which is transmitted to small x via momentum conservation and correlations in the EPPS16 fit function. The improvement for R A g (x, Q 2 ) in the full analysis is clearly larger. Thanks to the wider acceptance, the full pseudodata sample is able to provide better direct constraints also at lower x. In particular, the gluon distribution gets significantly better constrained, the level of improvement being of the order of 50% in places. We have found that here the most decisive factor is the forward-arm Drell-Yan data sample. Formally, the Drell-Yan production occurs via qq annihilation, but at small x the behavior of sea quarks is strongly driven by the gluons. Thus, the Drell-Yan data indirectly -but still rather strongly -constrains the gluons. This is in line with Ref. [43] (Sect. 10.4) where a similar phenomenon was observed in the case of W ± production. We note that the dijet and photon-jet pseudodata probe the mid-and high-x part of the nuclear PDFs. The uncertainties for these two observables are dominated by the assumed 5% uncorrelated bin-tobin systematic error and the obtained improvement in nuclear PDFs are dictated by this assumption. If systematic uncertainties like those achieved in p+Pb collisions at the LHC [29] could be reached, the impact would be clearly larger. In addition, the systematic uncertainty of the LHC measurements is almost always of correlated nature, but such correlation is difficult to estimate in advance. All in all, assuming a 5% uncorrelated systematic uncertainty appears thus a reasonable test scenario which should not overstate the constraining power.

IV. CONSTRAINING THE A DEPENDENCE OF NUCLEAR PDFS WITH LIGHTER IONS
The mass-number (A) dependence of the current nuclear PDFs is not well known -direct constraints exist only for large-x valence quarks and intermediate-x seaquarks. On one hand, e.g. in the EPPS16 analysis, the guideline has been that the nuclear effect should be larger for larger nuclei at the parametrization scale Q = m charm which then tends to lead to physically sound A systematics also at larger Q. On the other hand, in the recent nuclear-PDF analysis by the NNPDF collaboration [18] there is less direct control over the A dependence and thus the nuclear effects from one nucleus to another can fluctuate significantly. Due to the p+Pb and Pb+Pb collisions program at the LHC, the near-future improvements on nuclear PDFs are bound to be driven by the Pb nucleus. For example, the dijet [29], D-meson [44] and W ± [45] measurements efficiently constrain [23,46] the gluons in the Pb nucleus, perhaps providing even stronger constraints for Pb than what we have found in the present study for Au. The LHCb fixed-target mode facilitates measurements on lighter noble-gas targets [47], but only the very-high x regime of nuclear PDFs is accessible. However, e.g. in astrophysical applications the relevant nuclei (e.g. oxygen and nitrogen) are much lighter and thus collider measurements involving lighter nuclei would be very much useful [48]. In addition, the study of the onset of jet quenching and saturation phenomena with lighter ions will require nuclear PDFs for nuclei other than Au or Pb. Interests towards light-ion beams at the LHC have been expressed [43] but since the main focus of LHC is still in p+p collisions, it is not clear whether and when this would materialize. Here, the flexibility of RHIC to run with different ions is a clear asset. Indeed, at least p, d, Al, Cu, Ru, Zr and U ions have already been used in physics runs which demonstrates that a proper "A scan" is, in principle, possible. The same multi-ion option would also be available if RHIC is eventually turned into an Electron-Ion Collider [49], where the possibilities to constrain nuclear PDFs are undisputed [50,51]. To highlight the present uncertainties for light ions, Fig. 6 shows the nuclear effects for A = 40 (Ca, Ar) from the EPPS16 and nCTEQ15 [15] global fits of nuclear PDFs. While the uncertainty bands overlap, the shape at intermediate and large x are quite different; while the nuclear effects in nCTEQ15 monotonically rise towards high x, the EPPS16 error band more closely resembles the typical pattern of shadowing, antishadowing and EMC-effect. Fig. 7 shows how this different behaviour would be reflected in dijet production. In the backward direction (y dijet < 0) one is sensitive to the large-x part of nuclear PDFs and the nCTEQ15 prediction tends to be above the EPPS16 one, consistently with Fig. 6. The difference in Fig. 7 is not as marked as in Fig. 6 as towards large x the valence quarks also play an increasingly important role. Towards y dijet 0 the probed x gets lower and, in line with Fig. 6, the nCTEQ15 prediction tends to be at the lower limit of EPPS16. Assuming a similar ∼50% reduction in the gluon PDF uncertainties as found for Au in Fig. 5, it appears reasonable that the measurements would be able to resolve between nCTEQ15 and EPPS16. In an approach like that of the NNPDF collaboration [18], where more freedom for the A dependence is allowed than in nCTEQ15 or EPPS16, the benefit would be even more pronounced. An additional interesting possibility we would like to point out would be to study p+A i collisions of two isobaric nuclei A 1 and A 2 (e.g. A 1 = 96 44 Ru vs. A 2 = 96 40 Zr collisions) with constant A but differences in proton and neutron content. Precision measurements of e.g. (p+Ru)/(p+Zr) ratios for hard processes (like those discussed in this paper) would allow a study of the assumptions made in the present global fits of nuclear PDFs. Indeed, it is currently assumed that the nuclear effects depend only on the mass number A, and not on the mutual balance of neutrons and protons. In addition, the isospin symmetry (i.e. u proton/A = d neutron/A ) is assumed to be exact. Thus, (p+Ru)/(p+Zr) ratios, or other similar constant-A combinations, would test the assumptions made in global analyses at a deeper level and also test other theoretical approaches, e.g. the importance of short-range nucleon-nucleon correlations [52], or the lack of them [53]. In principle, in an optimal situation the neutron-to-proton mixture in the two nuclei should be as different as possible, with (at least nearly) constant A. Such a measurement makes optimal use of the flexibility of the RHIC facility.

V. SUMMARY
Using the Au nucleus as a test case, we have examined the prospects for constraining nuclear gluon PDFs at RHIC with new measurements that assume detector acceptances similar to those proposed for STAR and sPHENIX with forward upgrades. We have found that the Drell-Yan process at low invariant mass is able to significantly constrain the low x gluon distribution with up to 50% reduction in the current EPPS16 uncertainty. The constraints at higher x depend considerably on the assumed systematic uncertainty, which is expected to dominate over the statistical uncertainty for dijet and photon-jet processes. Assuming an order of 5% bin-tobin independent systematic uncertainty leads to modest constraints at mid-and high-x region. Even so, we find the inclusion of additional observables along with the Drell-Yan data of utmost importance to overcome the overall normalization uncertainty in the R pA ratio. Without supplementing the Drell-Yan pseudodata with other observables (here either dijets, photon-jet, or both), we find that the power of the measurement of Drell-Yan to constrain the small-x behavior of the gluon is lost. In addition, we note that Drell-Yan constrains nuclear gluons indirectly, and it is possible that even stronger constraints could be obtained if measurements of forward direct photons could be added to this suite of observables.
While the focus of our analysis was on the Au nucleus, similar constraints can be expected to be obtained for any other nucleus. In this respect we briefly discussed the A dependence of nuclear PDFs and highlight the significant opportunity for improvements that could be attained with a proper A-scan -measuring the same observables with several nuclear beams -for which the RHIC collider provides a unique opportunity.