Probing gluon helicity with heavy flavor at the EIC

We propose a new measurement of the heavy flavor hadron double spin asymmetry in deep-inelastic scattering at a future Electron-Ion Collider (EIC) to constrain the polarized gluon distribution function inside the proton. Statistical projection on $D^0$ meson double spin asymmetry is calculated with an EIC central detector using an all-silicon tracker and vertexing subsystem. A first impact study was done by interpreting pseudo-data at next-to-leading order in QCD. The sensitivity of the experimental observable in constraining gluon helicity distribution in a wide range of parton momentum fraction $x$ has been investigated considering different beam energy configurations. This measurement complements the inclusive spin-dependent structure function measurement and provides an opportunity to constrain the gluon helicity distribution in the moderate $x$ region.


I. INTRODUCTION
The spin structure of the nucleon has been of fundamental interest in modern hadronic physics ever since the EMC spin puzzle [1]. The current understanding of the structure of the nucleon spin is that it consists of the contributions from quark and gluon helicity distributions, as well as their orbital angular momenta [2,3]. For a longitudinally polarized nucleon, its spin can be decomposed into [3] where ∆Σ(x) and ∆g(x) denote the non-perturbative (longitudinally-)polarized Parton Distribution Functions (PDFs) for the quark singlet and gluons, L q and L g are the orbital angular momenta of quarks and gluons, and x is the momentum fraction carried by the quarks or gluons. The dependence of the PDFs on the factorization scale is left implicit here and for the most part in the rest of the paper. * dpa@m.scnu.edu.cn † felix.hekhorn@unimi.it ‡ yxzhao@impcas.ac.cn After more than 40 years of both experimental and theoretical efforts, the precision of the unpolarized PDFs is reaching higher and higher accuracy [4], but the polarized PDFs (pPDFs) are still not well constrained [5]. In addition, compared to the polarized quark distributions the polarized gluon distribution ∆g is a less known quantity. This is because the gluon only participates via QCD evolution and higher-order corrections to fully inclusive Deep-Inelastic Scattering (DIS) that most of the polarized experiments measure. In the last decade or so, the polarized proton-proton experiments at RHIC [6] have provided a stronger constraint on the polarized gluon distribution due to the fact that ∆g enters the differential cross section at the leading order in hard proton-proton scattering processes. Evidence for a positive ∆g in the region of x > 0.05 has been found. Due to the limited kinematic coverage one still cannot draw a decisive conclusion on the gluon spin contribution to the proton spin [7].
The planned Electron-Ion Colliders (EIC) [8,9] offer unprecedented opportunities for spin physics, especially to constrain the ∆g contribution to the proton spin via the scaling violations of the polarized structure functions [8]. Recently, there have been a series of studies to demonstrate the EIC impact on the pPDFs through various processes including semi-inclusive light hadron production [10,11] and jet production [12]. In this paper we provide a systematic study of charm quark production to constrain the polarized gluon distribution in a wide range of x at the future EIC. Compared to the inclusive DIS measurements, the polarized charm structure function provides direct access to the ∆g from leading order. This will complement the fully inclusive DIS measurements in several important ways, for example, offering a new ingredient on the ∆g determination in addition to the inclusive DIS and providing sensitivity in the moderate x region.
A similar measurement was proposed and performed by the COMPASS collaboration [13]. However, the COMPASS measurement yielded only one data point and was interpreted based on the approximation of photoproduction [14].
Our study is based on two important developments in the recent years. First, the next-to-leading order perturbative QCD formalism for heavy flavor production in polarized DIS has been derived [15]. This will help to achieve high precision from theory side in constraining ∆g from the experiments. Second, an all-silicon tracker conceptual design has been applied and demonstrated in various EIC simulations [16,17]. It also plays an essential role in our analysis of this paper, since it enables high precision measurement of heavy flavor hadrons through their hadronic decay channels.

II. THEORETICAL CALCULATION
In the theory calculation, we focus on electroproduction of inclusive open charm particles. The hadronization effects of the charm quark into the D mesons, electro-weak corrections, intrinsic charm components, and target mass corrections are currently not considered. The cross-sections for the unpolarized and polarized DIS processes are given in terms of three independent structure functions F 1,2 (x, Q 2 ) and g 1 (x, Q 2 ): from which we can define the charm double-spin asymmetry where the superscript c refers to the charm component of the structure functions, dσ ++ and dσ +− are the charm production cross-sections for electron and proton beam spin orientation to be parallel and anti-parallel to each other, respectively, and D(y) is the depolarization factor of the virtual photon depending on the inelasticity y. The target mass as well as the cross-section from longitudinal photon polarization are ignored in the above equations.
In the context of collinear factorization, the structure functions can be computed as a convolution of (p)PDFs (∆)f j and perturbative coefficient functions (∆)c k,j : where z max = Q 2 /(4m 2 + Q 2 ) is the kinematic boundary to create a charm quark pair in the final state with m the charm quark mass. Note that the argument of the PDF is x/z where x is the Bjorken-x and z is the convolution variable. The perturbative next-to-leading order (NLO) calculation of the partonic coefficient functions (∆)c k,j is known in the unpolarized case [18] for quite a while. Their polarized counterparts have become available only recently [15] after the previous leading order (LO) computation [19]. Heavy flavor production can constrain the gluon PDF since at LO the only contribution is photon-gluon-fusion (PGF) (see Fig. 1a) and in the case of unpolarized PDFs this is an established technique [20]. At NLO three different types of contributions have to be considered: real gluon radiation (see Fig. 1b), virtual corrections to PGF (see Fig. 1c), and light quark initiated contributions (see For the factorization and renormalization procedure we use the MS m scheme as described in [14] and choose the respective scales to be µ 2 F = µ 2 R = 4m 2 + Q 2 . We use a fixed flavor scheme with n l = 3 light flavors (uds) together with a pole mass prescription for the heavy charm quark. The actual value for the pole mass m and the prescription for the running coupling α s (µ 2 R ) is provided by the LHAPDF interface [21]. For the calculations of this paper, we use both the NNPDFpol1.1 [22] and DSSV14 [23] polarized PDF sets alongside with their unpolarized counterpart NNPDF23 nlo [24] and MSTW2008nlo [25], respectively.

III. PROJECTIONS FOR THE EXPERIMENTAL OBSERVABLE
Experimentally, the double spin asymmetry in the e + p → e + D 0 + X DIS process can be measured at an EIC as where N ++ and N +− are the luminosity-normalized counts for electron and proton beam spin orientation to be parallel and anti-parallel to each other, respectively, and P e (P p ) is the electron (proton) beam polarization. The beam polarization is assumed to be 80% for the electron beam and 70% for the proton beam at the EIC [8]. Therefore, one has where A c 1 can be calculated as discussed in the Section II. To demonstrate the general size of the double-spin asymmetry A c 1 , we show a representative plot in Fig. 2. An all-silicon tracking detector design [16] at an EIC enables the D 0 reconstruction with a very good signalto-background ratio. Moreover, the large acceptance and high luminosity available at an EIC allows the measurement to be done in a broad kinematic coverage in Bjorken-x and Q 2 .
A simulation study has been performed to obtain uncertainty projections of the experimental observable [16]. The geometry of a silicon tracking system has been implemented in GEANT4 and studied within the full Monte-Carlo framework for detector simulation. The full simulation yields the detector response tables for momentum resolution, single track pointing resolution, tracking efficiency and primary vertex resolution. Afterward, the resolution tables were implemented in a fast smearing simulation framework to allow for the generation of sufficient statistics to carry out detailed studies for physics projections. In our study, three beam energy configurations have been used for electron-proton collisions: 18 GeV×275 GeV, 5 GeV×100 GeV and 5 GeV×41 GeV.
The data were generated by pythiaeRHIC (PYTHIA V6.4) and then fed into the fast smearing framework to accommodate detector response within a 3 T magnetic field. We take advantage of the Kπ two-body-decay to identify the D 0 (D 0 ). Three decay topological distributions, namely, Kπ pair-DCA (distance of closest approach), D 0 Decay-Length rφ in the transverse plane, and the cosθ rφ where θ is the angle between the D 0 pointing direction with respect to the primary vertex and the momentum vector of the Kπ pair, were investigated to obtain the data sample with a good signal-tobackground ratio. In addition to the D 0 -decay topology cuts, the following kinematic cuts in the e + p collision including the squared momentum transfer of the electron Q 2 , the inelasticity y, and the invariant mass of the produced hadronic system W were used in the analysis: Q 2 > 2 GeV 2 , 0.05 < y < 0.8, and W 2 > 4 GeV 2 . The pion/kaon identification was assumed to be feasible up to the momentum limits 10 GeV/c, 6 GeV/c and 50 GeV/c in pseudo-rapidity regions (-3,-1), (-1,1) and (1,3), respectively.
After all the selection requirements had been applied, the data were binned in Bjorken-x. In each bin the reconstructed Kπ invariant mass spectrum was fit with a Gaussian function for signal plus a linear background to extract the number of D 0 signal and background (as shown in Fig. 3). Hence, the uncertainty of A c 1 can be calculated bin by bin, as shown in Fig. 4 for three different beam energy configurations. The weighted center for each data point is according to the Bjorken-x and Q 2 axis, while the size of the error is according to the scale on the right-side vertical axis. The integrated luminosity is corresponding to 100 fb −1 for each beam energy configuration. The electron (proton) beam polarization is assumed to 80% (70%). The uncertainties become larger in the lower beam-energy configuration is due to the decrease of the production cross-section for the charm    The integrated luminosity is 100 fb −1 for each configuration in this plot. The electron (proton) beam polarization is assumed to be 80% (70%). The position of each data point in the plot is defined by the weighted center of Bjorken-x and Q 2 for each particular bin. The uncertainty indicated for each data point should be interpreted using the scale shown on the right-side vertical axis of the plot.
quark. For each beam energy configuration, the uncertainty becomes larger in the higher x region due to the smaller depolarization factor D(y).

IV. PROJECTIONS FOR THE PARTON HELICITY DISTRIBUTIONS
In global data fits, parton distribution functions are expressed in terms of some functional form depending on a number of free parameters whose values are constrained by the experimental data. A practical way of transferring this information to parton distribution functions themselves is to express them as a set of so-called "replicas" generated by means of Monte Carlo sampling of the pa-rameter space. Central values and uncertainties of PDFs become then simple statistical mean values and standard deviations of the full replica set. To assess the impact of the new measurement proposed in this paper on the quark and gluon helicity distributions without performing a full refit, we use the well-established reweighting method of [26]. By exploiting Bayesian inference, the information contained in the new set of data can be incorporated directly into the probability distribution of the initial helicity PDF replicas. More specifically, this is achieved by assigning a weight to each replica measuring its consistency with the new data. The resulting new PDFs containing the information of the new data set preserve the statistical rigor of the original set as long as not too many replicas become suppressed by vanishing weights. A small number of surviving replicas would mean that the impact of the new data is too significant for the reweighting method to work and a full fit is necessary. In our study, two commonly used polarized PDF sets were used separately for the reweighting analysis: NNPDFpol1.1 [22] and DSSV14 [23]. NNPDFpol1.1 is given in the form of 100 replicas while DSSV14 is given with 1000 replicas. Although replicas from both sets are generated using Monte Carlo sampling methods, the two sets profoundly differ in the way the shape of PDF replicas are parametrized. NNPDFpol1.1 uses functional forms provided through the use of neural networks, i.e. with a high number of free parameters. DSSV14 uses a more traditional analytical, although flexible, functional form with much fewer free parameters.
Without a real measurement, one doesn't know the central values of the data points. Therefore, the pseudodata were generated by randomly displacing the theoretical central values using the projected uncertainty. Reweighting is performed using pseudo-data generated for the above-mentioned three different energy configurations and corresponding integrated luminosity of 100 fb −1 . Depending on the energy configuration, after reweighting the surviving replicas are around ∼ 70 and ∼ 850 for NNPDFpol1.1 and DSSV14 respectively, which is a large enough number to justify the use of reweighting in this analysis.
In Fig. 5, we show the impact of the EIC pseudo-data on the uncertainties of the singlet quark helicity distribution ∆Σ(x) and the gluon helicity distribution ∆g(x) at Q 2 = 10 GeV 2 with 100 fb −1 of integrated luminosity. The gray band represents the original NNPDFpol1.1 or DSSV14 absolute uncertainty band before reweighting. Each colored band represents the effect of reweighting using pseudo-data generated with one of the three possible collision energy configurations and all of them combined.
In the bottom area of the plots we also show the ratio between the uncertainties before and after reweighting. For the study with NNPDFpol1.1, one can clearly see the en-ergy dependence of the pseudo-data impact on the gluon helicity distribution: higher collision energy data offer more constraints on the pPDF uncertainty in the lower values of x, and vice-versa. In the quark sector, the two lower center-of-mass energy configurations have less constraining power. While for the study with DSSV14, the impact of 18 GeV×275 GeV and 5 GeV×100 GeV pseudodata sets is similar in both quark and gluon sectors. Interestingly, the impact on ∆g is in the x > 0.1 region for the 5 GeV × 41 GeV configuration, which is a novelty of this measurement. The difference for the resulted impact of EIC pseudo-data by using NNPDFpol1.1 or DSSV14 is mainly due to two reasons. The first one is due to the different data sets included in different global fits. NNPDF-pol1.1 fit contains the world data with x extension to about 4 × 10 −3 including DIS data, open-charm production data from the COMPASS experiment at CERN and high-p T inclusive jet and π 0 , as well as W ± production data from the STAR and PHENIX experiments at RHIC. In addition to the DIS data, the DSSV14 also includes SIDIS data, inclusive jet and identified hadron production measurements from polarized proton-proton collisions at RHIC. The second reason is because of different parameterizations for the shape of quark and gluon helicity distributions. One obvious effect of different parameterizations is the determined uncetainty band beyond the coverage of existing world data, for example in the very low-x region, where the NNPDFpol1.1 shows a significant larger uncertainty compared to the DSSV14 case in both quark and gluon sectors. That is why the high center-ofmass energy configurations show less impact in the relatively low-x region for the DSSV14 study compared to the NNPDFpol1.1 case. Note also that the x value where the impact is large extends beyond the Bjorken-x reach of the data (see Fig. 4), this is because of the shift be- tween Bjorken-x (determined by the virtual photon) and parton-x (x/z in Eqs. (7) and (8)) in the PGF process. The increase of the uncertainty band after reweighting for some values of x in either quark or gluon sector is due to the fact that the reweighting procedure is favoring replicas with an appropriate shape, which is mainly determined by the fixed parameterizations. Similar reason also applies to the situation while looking at the impact by combining all the three pseudo-data sets, compared to the impact of individual pseudo-data set. Resolving such bias is outside the scope of this study and is better left to future fits with real EIC data.
When it comes to helicity PDFs, the first moments, i.e. their integrals over parton momentum fraction (see Eq. (1)) are quantities of great interest. They represent the net quark and gluon contribution to the proton spin [3]. In every practical case, where the x of the data is limited and the integration cannot be performed down to x → 0, the truncated first moments are usually used to represent the contribution to the proton spin down to the x min fraction of proton momentum accessible by experimental data. In Fig. 6 and Fig. 7 we present, for NNPDFpol1.1 and DSSV14 PDF sets respectively, the impact of different EIC pseudo-data sets on truncated first moments for the quark and gluon helicity distributions as a function of the lower integration limit x min . In the same way, we also show the missing contribution to the proton spin which is usually associated to the quark and gluon orbital angular momenta spin contributions [3]. The bottom panel of each plot shows the uncertainty improvement by including a particular set of pseudo-data. Those plots are very instructive as they explicitly show which contribution to the integral is mostly affected by the data. From Fig. 6 we can observe that the   two lower energy configurations are able to target specific regions of the gluon spin contribution to the proton spin. Choosing the 5 GeV × 41 GeV configuration, the precision of our knowledge of the contribution to the proton spin from gluons with momentum fraction down to the intermediate-x region is increased up to a factor of 1.5. If we choose 5 GeV × 100 GeV configuration, this goes up to a factor of 2 around 10 −3 x min 10 −2 . On the other end, the higher energy configuration is able to constrain the quark and gluon spin contributions to the proton spin in the x min region below 10 −3 . Using DSSV14 PDFs set (c.f. Fig. 7) leads to slightly different results. In particular, compared to the NNPDFpol1.1 case, we observe for the 5 GeV × 41 GeV configuration a larger impact in the high-x min and a general smaller impact of the 18 GeV × 275 GeV configuration for the quark sector. This can be well understood by looking at Fig. 5: if we compare the uncertainties of Fig. 5a versus Fig. 5b, the difference of impact of the lowest configuration energy in the quark sector between Fig. 6 and Fig. 7 can be traced back to the difference in the original PDF uncertainties in the high-x region. Similarly, the smaller impact in the small-x min region of the 18 GeV × 275 GeV data is related to the largely different uncertainties of NNPDF-pol1.1 and DSSV14 in the low-x region. As gluon and quark sectors are typically correlated, we also present the first moments of gluon and quark helicity distributions in a two-dimensional plot. The plots in Fig. 8 and Fig. 9 show the first moments of the gluon helicity distribution at Q 2 = 10 GeV 2 as a function of the quark helicity distribution for the three collision energy configurations. For each plot x min is given by the lowest Bjorken-x accessible to the specific energy configuration setting. The black ellipse represents the original 68% C.L. uncertainty boundary before reweighting. The red ellipse represent the 68% C.L. boundary when including pseudo-data with integrated luminosity of 100 fb −1 . Each point in the plot is associated to a specific replica and their colors represent the magnitude of the weights associated after including pseudo-data. The dots in yellow represent the replicas dominating the distributions. They are clustered around the range of values that become relevant once EIC pseudo-data are inserted. The central position of this cluster is not physically relevant at this stage, as it is directly correlated with the distributions' central value shift which, as discussed above, takes meaning only once real experimental data are considered. We can notice that the cluster gets gradually squeezed as the beam energy goes higher at an EIC. Moreover, for the study with DSSV14 PDF sets, the ellipses after including the EIC pseudo-data have different angles comparing to the original distribution, which means that the new observable offer independent ingredients, more specifically gluon-sensitive inputs, into the world data in the DSSV14 global fit.

V. SUMMARY
We have proposed a new measurement on longitudinal double spin asymmetries in the e + p → e + D 0 + X DIS process at an EIC to constrain the gluon helicity distribution ∆g. We would like to emphasize that the classic g 1 measurements at the EIC will play the dominant role to constrain the gluon helicity distribution and its contribution to the proton spin [8]. Our proposal will provide complementary constraints on the gluon helicity distribution. As we show in the impact study, in some kinematics, e.g., moderate x region, heavy flavor production will offer a unique opportunity. Especially, with a lower center-of-mass energy machine, one can improve the precision of gluon helicity distribution in the x > 0.1 region.
Moreover, the theoretical calculation shows that the A c 1 is sizable at the level of 10-20% in the moderate x region, which is an advantage experimentally. To achieve the measurement with good signal significance for the D 0 reconstruction, a state-of-the-art silicon pixel tracking detector system with wide pseudo-rapidity coverage, excellent momentum and spatial resolutions as well as a low mass budget is essential. It has been shown in our study and the work in [16,17] that the proposed all-silicon tracker conceptual design is well suited for the measurement.