Differential Heavy Quark Distributions and Correlations in Longitudinally Polarized Deep-Inelastic Scattering

We present a first calculation of the heavy flavor contribution to the longitudinally polarized deep-inelastic scattering structure function $g_1^{Q}$, differential in the transverse momentum or the rapidity of the observed heavy quark $Q$ or antiquark $\overline Q$. All results are obtained at next-to-leading order accuracy in QCD within the framework of a newly developed parton-level Monte Carlo generator that also allows one to study observables associated with the produced heavy quark pair such as its invariant mass distribution or its correlation in azimuthal angle. First phenomenological studies are carried out for various heavy quark distributions in a kinematic regime relevant for a future Electron-Ion Collider with a particular emphasis on the expected size of the corresponding double-spin asymmetries and their sensitivity to the still poorly constrained helicity gluon distribution. Theoretical uncertainties associated with the choice of the factorization scale are discussed for selected observables.


I. INTRODUCTION AND MOTIVATION
Heavy flavors, more precisely their contributions to deep-inelastic scattering (DIS) processes or their production in hadron-hadron collisions, have played a rather minor role so far in constraining longitudinally polarized parton distributions (PDFs) [1,2] despite their sensitivity to the elusive gluon helicity density. The prospect of a future high luminosity Electron Ion Collider (EIC) in the U.S. [3], with its vastly extended kinematic reach compared with past and present fixed-target experiments [1,4] will likely change this assessment considerably. For the first time, experiments will be able to probe the helicity structure of nucleons with unprecedented precision at momentum fractions x well below the currently explored range x 0.01 and, at the same time, at energy scales Q above 1 GeV, i.e., in the domain where perturbative QCD should be safely applicable and where contributions from heavy quarks (HQs) are potentially sizable.
Clearly, a proper treatment of HQs in polarized DIS will become a must in any global QCD analysis of helicity PDFs [5][6][7] based on future EIC data [8]. Experiments at the DESY-HERA collider have shown in the past [9] that the charm contribution to unpolarized DIS, i.e., to the structure function F 2 , can reach about a level of 25% at small x. Also, at not too large scales Q, as compared to the HQ mass m, one needs to retain the full dependence on m in all theoretical computations as a massless approximation, i.e., treating charm quarks as massless partons in the proton, is not warranted here. This is primarily also the kinematic region of utmost significance * felix.hekhorn@unimi.it † marco.stratmann@uni-tuebingen.de for investigating polarized DIS at the EIC [3]. The presence of different scales, at least m and Q, significantly complicates all calculations involving HQs in DIS already at the level of next-to-leading order (NLO) accuracy, in particular, the evaluation of the necessary virtual corrections and phase space integrals. This is even more so the case if one considers, in addition, the polarization of the initial-state partons because of the well-known complications due to the presence of γ 5 and the Levi-Civita tensor in the helicity projection operators in n = 4 dimensions [10].
An important step towards completing the already existing suite of NLO HQ production cross sections in longitudinally polarized hadro- [11,12] and photoproduction [13,14] has been made recently in Ref. [15] where the heavy flavor contributions to the inclusive helicity DIS structure function g 1 have been computed more than 35 years after the leading order (LO) expressions were first derived [16]. This paper will build on these NLO results that were obtained with largely analytical methods following closely the corresponding, pioneering calculations in the unpolarized case [17,18].
In the present work we shall develop a new partonlevel Monte Carlo (MC) generator that allows one to study systematically also more exclusive observables in longitudinally polarized DIS, i.e., quantities which are differential in the transverse momentum p T or the rapidity y of the observed heavy (anti)quark Q (Q ) or which are associated with the produced HQ pair such as its invariant mass distribution or its correlation in azimuthal angle. All phase space integrals are computed in four dimensions by means of standard MC techniques [19]. Collinear and soft divergencies at intermediate stages are made explicit by generalized plus distributions within the framework of n-dimensional regularization [20] and cancelled upon adding the renormalized virtual corrections arXiv:2105.13944v1 [hep-ph] 28 May 2021 and appropriate factorization counter terms. Contrary to our results in [15] and thanks to the adopted subtraction method [21] there is no need to actually evaluate integrals in n dimensions.
We benefit greatly from previous applications of the subtraction method to various HQ calculations at NLO accuracy [12,14,[22][23][24], which provide essentially all of the formalism required here. In particular, we follow closely the methods and notations outlined in Refs. [22,24] and shall only very briefly review the main technical aspects in this paper in order to adopt them to the case of helicity-dependent DIS. We will utilize the results of our previous, largely analytical calculation [15] to validate our numerical results and the underlying MC code for fully inclusive and single-differential HQ processes which are amenable to both analytical and numerical methods. We note that our code also provides unpolarized results. They appear in the numerator of experimentally relevant double-spin asymmetries and, hence, are needed in our phenomenological studies. Whenever possible, we use existing unpolarized results in the literature [24][25][26][27] for further validation of our MC code.
Our MC implementation of polarized HQ electroproduction is in principle capable to produce any infrared safe observable, i.e., histograms for exclusive, semiinclusive, or fully inclusive quantities related to any of the outgoing particles or combinations thereof. Experimental cuts can be implemented as well if needed, provided they can be expressed in terms of the available partonic variables. The presented MC program and numerical results are based on the virtual photon-hadron (γ * h) center-ofmass system (c.m.s.) frame, i.e., they make no reference to the lepton in DIS and are given solely in terms of the relevant differential helicity structure function dg 1 . In a next stage, our parton-level MC code can be straightforwardly expanded further, for instance, along the lines of the unpolarized HVQDIS code [28] by interfacing it with appropriate sets of HQ fragmentation functions to model the HQ decay into heavy mesons.
As a first phenomenological application of our MC code, we will explore in this paper some HQ observables of potential interest for the physics program at a future EIC. One of the main goals of the EIC [3] is to collect vital new information about the elusive helicity gluon density ∆g from as many processes as possible and to exploit these data in future global QCD analyses. Since HQ electroproduction is driven by the LO photon-gluon fusion (PGF) process, γ * g → QQ , it is expected to play a major role in this exercise. In Ref. [15] we have studied already the fully inclusive contribution of charm quarks to the DIS structure function g 1 in the kinematic domain x 0.01 and moderate Q 2 relevant for the EIC. Here, we will focus on various differential HQ distributions and correlations, including the experimentally relevant double spin asymmetries, and compute them for the entire suite of helicity PDF uncertainty sets provided by the DSSV collaboration [6].
It will turn out that the single inclusive distribution dg 1 /dp T as well as the invariant mass spectrum of the produced HQ pair, dg 1 /dM , are particularly sensitive to ∆g. HQ correlations which are trivial at LO accuracy, for instance, the back-to-back configuration in azimuthal angle, are also of phenomenological interest as they test higher order QCD corrections most clearly. These types of observables are also particularly sensitive to all-order resummations, i.e., the treatment of multiple soft gluon emission in the vicinity of, say, the back-to-back peak. We should stress that these corrections are not included in our fixed order MC code and, hence, results in these regions of phase space should be taken with a grain of salt; soft gluon resummations in the context of HQ electroproduction have been reported in Refs. [29]. For the phenomenologically most interesting HQ observables we shall also investigate and estimate the remaining theoretical uncertainties at NLO accuracy due to variations of the renormalization and factorization scales. Finally, we wish to remark that for phenomenological applications at a future EIC at small x and for Q not much larger than the HQ mass m, HQ production is most likely best described by retaining the full mass dependence as has been done in our largely analytical calculation [15] and for all numerical results presented here. This also implies that HQs can be only produced extrinsically, for instance, through the PGF mechanism and that the notion of a massless HQ PDF makes no sense 1 . However, in the asymptotic regime, Q m, one might want to consider HQ parton densities which requires to set up some interpolating scheme that matches a theory with n lf light flavors to a theory with n lf +1 massless flavors [31]. Different types of general-mass variable flavor number schemes (GM-VFNS) have been proposed and adopted in global fits of unpolarized PDFs [32,33] but no such prescription has been considered so far for fits of helicity PDFs. Our results, along with the asymptotic limit Q m [34], are necessary ingredients for a polarized GM-VFNS.
The remainder of the paper is organized as follows: in Sec. II we briefly sketch the relevant technical aspects of setting up a parton-level Monte Carlo generator for heavy quark production in longitudinally polarized DIS. The phenomenological results are collected in Sec. III. Here, we first validate our MC code against our previous results based on largely analytical methods. Next, we present detailed results for the transverse momentum and rapidity-differential single-inclusive distributions of the longitudinally polarized DIS charm structure function g their difference in azimuthal angle, their combined transverse momentum, and the cone size variable. The main results of the paper are summarized in Sec. IV.

II. TECHNICAL ASPECTS
In this Section, we will briefly review the main technical aspects pertinent to our goal of setting up a partonlevel MC program for HQ production in polarized DIS. As was already mentioned, we will heavily make use of existing methods adopted for various other calculations of HQ production at NLO accuracy [12,14,[22][23][24]35]. In particular, Refs. [22,24] provide essentially all of the formalism required here, albeit for the unpolarized case.
Photon-gluon fusion, the sole mechanism for HQ electroproduction at LO accuracy [16], receives both virtual and gluon bremsstrahlung O(α s ) corrections, and is supplemented by a genuine NLO light [anti]quarkinduced process where the relevant four-momenta are labeled as q, k 1 , k 2 , p 1 , and p 1 . The photon is virtual, q 2 = −Q 2 , and the heavy quarks are taken to be massive, p 2 1,2 = m 2 , while k 2 1,2 = 0. The spin-dependent matrix elements for all contributing 2 → 2 and 2 → 3 processes have been computed already in our previous paper [15], but rather than performing the needed phase space integrations by largely analytical means we shall now utilize MC methods. This will enable us to study any infrared safe observable in spin-dependent HQ electroproduction up to O(αα 2 s ) even when experimental cuts are imposed. Here, α denotes the electromagnetic and α s the scale-dependent strong coupling. We shall stress already at this point that results obtained with our parton-level MC code for HQ electroproduction will fail whenever fixed-order results are bound to fail, i.e., when they become sensitive to the emission of soft gluons; some examples will be given in Sec. III.
To compute all phase space integrals in four dimensions by means of standard MC techniques [19] it is necessary to cancel any intermediate soft or collinear singularities under the integral sign beforehand. This is achieved by employing a variant of the subtraction method [21,22,24,36] and expressing the matrix elements and the twoand three-particle phase space factors, dPS 2 and dPS 3 , respectively, in terms of the variables x, y, θ 1 , and θ 2 that make the kinematic regions of soft and collinear emissions explicit. x is the reduced invariant mass of the HQ pair, s 5 ≡ s 5 − q 2 = (p 1 + p 2 ) 2 − q 2 scaled by the reduced γ *parton c.m.s. energy squared s ≡ s − q 2 = (q + k 1 ) 2 − q 2 , i.e., ρ * ≡ (4m 2 − q 2 )/s ≤ x ≤ 1, and −1 ≤ y ≤ 1 is the cosine of the angle between q and k 2 in the frame where q+ k 1 = 0. Both θ 1 and θ 2 range between 0 and π and are used to parametrize the spatial orientation of k 1,2 with respect to the plane span by the other three momenta in the c.m.s. frame of the HQ pair. They do not matter for the discussion of singular regions of phase space.
In the matrix element squared for (2) the emission of a soft gluon corresponds to a 1/(1−x) singularity as x → 1 and, likewise, a collinear 1/(1 + y) pole is encountered for y → −1. The latter also appears in the diagrams contributing to process (3) that are proportional to the electric charge of the HQ squared, e 2 Q . The gist of the subtraction method adopted here [22,24] is to tame the singularities by multiplying the squared matrix elements in (2) and (3) by an appropriate Lorentz-invariant factor that vanishes in the soft and/or collinear limit, yielding finite functions that can be integrated numerically over x, y, θ 1 , and θ 2 . The phase space dPS 3 , divided by the same factor, can be expressed in terms of generalized plus distributions within the framework of n = (4 + )-dimensional regularization, and all poles can be cancelled either upon adding the virtual contributions or by applying a factorization counterterm before integration. We will briefly outline the most important steps in what follows. For more details, we refer the reader to Refs. [24,37] Adopting our notation from Ref. [15] for denoting the appropriate projections onto the different HQ contributions (i.e., structure functions) in unpolarized (G, L) and longitudinally polarized (P ) DIS by a subscript k = {G, L, P }, we can define the modified 2 → 3 matrix element squared for PGF (2) schematically as where the the singularities in the limits y → −1 and x → 1 are now regulated.
The specific choice only affects the rate of convergence and stability of the numerical MC integrations. For all our purposes we can adopt the same values as in the corresponding unpolarized HVQDIS code [24,27,28]. dσ (1),s k,g and dσ (1),c− k,g exhibit soft and collinear singularities, respectively, that have been made manifest in n dimensions by introducing theρ and ω prescriptions. To proceed, one notices that the 1/ pole in dσ (1),c− k,g in the limit y → −1 assumes the form dictated by the factorization theorem and, hence, can be absorbed into the definition of the PDFs. To this end, one adds an appropriate "counter cross section" to dσ which is a convolution of the n-dimensional LO gluon-gluon splitting function P where β lf 0 = 11 6 C A − 1 3 n lf with C A = 3 and n lf denoting the number of active light quark flavors. The splitting function P , and a hard (H) part as well as into fourdimensional and O( ) pieces as follows Combining everything, the finite, factorized collinear cross section dσ where dPS x 2 ≡ dPS 2 | s→s 5 dx denotes the relevant twoparticle phase space factor for collinear kinematics with QED is related to the Born matrix element squared for the PGF process, see, e.g., Eq. (13) in Ref. [15]. Equation (15) carries a dependence on the arbitrary factorization scale µ F where the collinear subtraction is performed. We note, that for finite Q 2 there is no collinear divergence associated with the limit y → 1 that requires extra attention in the case of photoproduction.
For k = P the finite soft plus virtual cross section reads The result has no dependence on x and y, nor on θ 2 . The function f P contains all remaining finite contributions, mainly expressed in terms of various logarithms and dilogarithms, but independent ofρ, the scales µ F and µ R , and β lf 0 . Finally, the cross section dσ where we have limited ourselves to k = P as before.
In summary, the PGF process up to NLO accuracy can be expressed as where now all terms are finite and, hence, all phase space integrations can be performed with MC methods in four dimensions. We note that each of the three NLO contributions in (19) individually depends onρ and ω but not their sum. To generate the predictions for HQ structure functions in DIS at the hadronic level, which will be shown and discussed in the next section, one needs to convolute the partonic result in (19) with the appropriate PDF, i.e., in case of k = P with a set of helicitydependent PDFs.
Finally, the genuine NLO contribution from initialstate light quarks (3) is conceptually much simpler to obtain than for the PGF process. To this end, one notices that (3) exhibits only a collinear singularity for y → −1 and, hence, one starts by defining the modified matrix element and appropriately rescaled phase space factor as and dPS 3,q ≡ dPS 3 2s respectively. As before, the divergent factor (1 + y) −1+ that emerges in dPS 3,q can be regularized by introducing a generalized plus distribution. The collinear singularity in dσ is lifted by adding a counter term containing the LO quark-gluon splitting function P (0),ρ k,qg (x) and the dσ (0) k,g (x k 1 ). As a result, the finite light quark induced process at NLO accuracy can be written as and To obtain the contribution of (3) to the hadronic HQ structure functions one needs to perform the convolution of dσ k,q with appropriate combinations of light quark PDFs.

III. PHENOMENOLOGICAL STUDIES
We now turn to a detailed phenomenologically study of various heavy flavor observables in longitudinally polarized DIS that are of potential interest for the physics program at a future EIC [3]. As in our previous study on the inclusive DIS structure function g 1 in Ref. [15], the main focus will be on the sensitivity of the corresponding double-spin asymmetries to the helicity gluon density ∆g in the kinematic range 10 −3 x 5 × 10 −2 that will be experimentally explored for the first time at the EIC.
To this end, we adopt for all our numerical calculations the helicity PDFs of the DSSV group along with their uncertainty estimates [5,6]. In particular, we will highlight the results obtained by utilizing the two most extreme sets that provide the largest excursions from the DSSV best fit gluon PDF in the relevant small x region. Throughout the remainder of the paper, these two sets are labeled as "max ∆g" and "min ∆g", respectively. The DSSV best fit gluon PDF [5,6] is shown in Fig. 1 along with the corresponding uncertainty estimate at 90% confidence level (C.L.) and the two extreme sets. As can be seen, in the small x-region of interest, the sets "max ∆g" and "min ∆g" form the envelope of the uncertainty band. This remains true for all relevant scales µ F used in our study. It should be remarked that for . The solid, dashed, and dot-dashed line represents the best fit, the "max ∆g", and "min ∆g" set, respectively. The shaded area is the DSSV uncertainty estimate for x ∆g at 90% C.L.
those members of the DSSV uncertainty band that exhibit a node at small x, the position of the node moves to smaller values of the momentum fraction with increasing scale µ F .
Since the DSSV sets are only available at NLO accuracy, we have to use them also in all our calculations at LO accuracy and, hence, for the corresponding "Kfactors". To compute all the unpolarized DIS structure functions and HQ distributions that appear in the denominator of the experimentally relevant double-spin asymmetries to be defined below, we adopt the NLO 90% C.L. PDF set of the MSTW group [39]. This set also serves as the reference set in the DSSV global analysis in ensuring compliance with the positivity limit for helicity PDFs. We also use the values of the strong coupling α s (µ R ) at NLO accuracy as determined in the MSTW fit [39]. The choice of the factorization and renormalization scales, µ F and µ R , respectively, depends on the HQ observable under consideration and, hence, will be indicated in each case below. For all fully inclusive HQ distributions we choose [15] We are only going to consider (anti-)charm quark electroproduction throughout this paper.
It should be also mentioned that all kinematic quantities such as transverse momenta, rapidities, and angles of the observed heavy quark and/or antiquark refer to the virtual photon-proton (γ * p) c.m.s. frame, where the positive z-axis coincides with the momentum direction of the hadron. The latter choice differs from the one taken in Ref. [25], where the γ * travels in the z-direction. More details and explicit formulas related to kinematical considerations comprising the choice of momenta, relevant Jacobians and expressions for the various integration limits can be found in Refs. [25,27,37].
As a first important step, we validate our newly developed parton-level MC code against the results already known from our previous calculations based on largely analytical methods [15]. Such a comparison is possible for any single-inclusive HQ observable in DIS that can be expressed in terms of the transverse momentum p T or transverse mass m T = p 2 T + m 2 and the rapidity y of the observed heavy (anti)quark; see, e.g., the appendices in Refs. [25,37] for details concerning the appropriate Jacobians for the required changes of integration variables. As an example, we consider in Fig. 2 the fully inclusive charm contributions F c 1 and g c 1 to the corresponding unpolarized and polarized DIS structure functions F 1 and g 1 , respectively, which have been studied extensively in [15]. The upper and lower panels of Fig. 2 show the relative differences of the numerical results for F c 1 and g c 1 obtained with the largely analytical code of Ref. [15], labeled as "inc", and the MC method described in this paper. Results are presented as a function of x for two different values of Q 2 . We shall note that the necessary multi-dimensional MC integrations are performed with an enhanced adaptive Vegas routine, see Ref. [19]. To obtain numerically stable results, it suffices to use 500k points and 5 iterations at LO accuracy and 4M points and 20 iterations for all computations involving NLO corrections.
As can be seen, the numerical differences for F kinematic range shown in Fig. 2 even close to one permille. Differences for g c 1 are somewhat larger, but usually better than 1 − 2%, except in the vicinity of nodes where g c 1 is numerically very small and changes sign. The comparisons in Fig. 2 reflect the typical numerical accuracy that can be achieved with the chosen number of points and iterations in the Vegas integration. In general, for any single-inclusive HQ observable, see Ref. [37] for more examples, the integrations leading to the unpolarized result are, as expected, noticeably more stable than the corresponding ones for the polarized distributions. The latter may exhibit nodes and, in addition, suffers from significantly larger numerical cancellations between the different contributions to the full NLO result. We wish to stress here, that we have validated our MC code also against all the single-inclusive HQ distributions shown below, with similar numerical differences to the largely analytical code as documented in Fig. 2.
After this prelude, we shall investigate the p T -and ydifferential single-inclusive distributions of the longitudinally polarized DIS charm structure function g c 1 in Figs. 3 and 4, respectively. In both cases we choose a fixed momentum fraction x = 0.001 and set Q 2 = 10 GeV 2 , which are both well within the interesting kinematic range to be explored for the first time at a future EIC. As was mentioned already, the results in this paper have been obtained for a detected heavy anticharm quark. We note that there is in fact a very small asymmetry between observing the heavy quark or the heavy antiquark distributions at NLO accuracy [15,18], which is, however, unlikely to be of any phenomenological relevance at the EIC in case of polarized DIS. For all unpolarized p T -and y-differential distributions that appear in the denominator of the relevant double-spin asymmetries, see Eq. (25) below, we have obtained excellent numerically agreement with the results given in [25] to which we refer the inter- ested reader for a discussion of the most important findings in the unpolarized case. Of course, for all our purposes we have updated the relevant results in Ref. [25] using the MSTW set of PDFs and conventions for the strong coupling, but we refrain from reproducing them here. Figure 3 shows 2x dg c 1 /dp T (left-hand-side) and the corresponding double-spin asymmetry (right-hand-side). The upper and middle row refer to the calculations at LO and NLO accuracy, respectively, and the bottom panels illustrate the K-factor, i.e., the ratio of the NLO and LO results, as a function of p T . We have restricted p T to range where measurements appear to be conceivable at the EIC for the given values of x and Q 2 , and rapidity is integrated over the entire phase space available. In any case, small values of p T correspond to probing ∆g at the smallest momentum fractions kinematically allowed, which is one of the main goals of the EIC spin program. More precisely, for the typical kinematics, x and Q 2 , selected in Fig. 3, the values of momentum fraction in ∆g that are predominantly probed in the convolutional integral with the hard coefficient functions is approximately in the range between 0.015 and 0.05, i.e., about a decade larger than the chosen value of x; see the appendix of Ref. [37] for a collection of all relevant equations entering in this estimate. In addition, the limitation to moderate values of p T has the extra benefit that one does not have to worry about any potentially large logarithms of the type ln(p T /m) and their eventual resummation to all orders in perturbation theory at this stage. This might be an interesting subject to pursue in the future, following the work already available in the unpolarized case [40].
In contrast to the fully inclusive calculation shown in Fig. 2, and following Ref. [25], we now choose for all single-inclusive calculations the scale µ 2 = µ 2 F = µ 2 R = Q 2 + 4(m 2 + p 2 T ) to account for the presence of the additional hard momentum p T . If one were to choose a p T -independent scale, such as µ 2 = Q 2 + 4m 2 , in the calculation, one would observe even larger QCD corrections, in particular, towards higher values of p T , destabilizing the perturbative expansion; a similar observation was made in the unpolarized case, see Ref. [18]. Below, we shall discuss in more detail the impact and relevance of scale variations for dg c 1 /dp T . From the panels in Fig. 3 one can infer that all gluon distributions ∆g that lie within the uncertainty estimate provided by DSSV (shaded bands) produce a node in the p T -distribution both at LO and NLO accuracy somewhere in the range p T ∈ (1.5, 2.0) GeV, naturally accompanied by large QCD corrections in their vicinity. The feature of the common node has no simple kinematic explanation but is deeply rooted in the convolutional integral of the PDFs with the oscillating coefficient functions and their functional properties; see Ref. [15]. In the given range of p T , where the γ * g-subprocess dominates, the two extreme sets min and max ∆g highlighted in Fig. 1 also provide the envelope of the p T -differential structure function. In general, a large positive ∆g leads to the smallest QCD corrections as can be seen best from the bottom panels. However, as was already noted in case of the fully inclusive charm structure function g c 1 in Ref. [15], the NLO corrections are substantial throughout and should be included in any future global QCD analysis including HQ electroproduction data. Since the NLO corrections strongly depend on ∆g and p T one must not resort to any approximations, such as a constant K-factor in future fits.
The experimentally relevant A c 1,p T is larger than in the fully inclusive case, see Ref. [15], and can amount to a few percent. Even an initial measurement with a absolute accuracy of about 0.5% at the EIC would prove to be very valuable in determining the x-shape of ∆g more precisely. This would be true, in particular, if such a measurement could be carried out for at least two different bins of p T , ideally to the left and to the right of the node in A c 1,p T to verify if the sign change predicted by all members of the DSSV uncertainty band for ∆g is indeed realized or not.
Analogously, Fig. 4 presents another single-inclusive distribution of potential interest for the physics program at the EIC: 2x dg c 1 /dy (left-hand-side). As before, the right-hand-side of the plot gives the corresponding double-spin asymmetry A c 1,y , defined analogously to Eq. (25), and the DIS kinematics is again fixed to x = 0.001 and Q 2 = 10 GeV 2 . The transverse momentum of the observed charm antiquark is integrated over the entire phase phase, and the rapidity y is defined in the γ * p c.m.s. frame with negative y pointing in the direction in which the γ * travels. We restrict ourselves to the region of −4 y −2.5 where ∆g is probed at the smallest momentum fractions for any given p T . As for the p T -distribution, different ∆g from the DSSV uncertainty estimate lead to different predictions for 2x dg values of y, the differences in A c 1,y quickly diminish and within the DSSV uncertainty estimate a small negative A c 1,y is expected. Given the anticipated unprecedented, high luminosity of the EIC, it might be even worth studying double-differential distributions in both p T and y for some suitable binning in the future to further optimize the sensitivity of HQ DIS observables to ∆g. We note that the somewhat erratic behavior at NLO observed for y < −4 is due to end-of-phase space effects, i.e., in this region only the emission of soft gluons is kinematically possible. Corresponding large logarithmic terms would need to be resummed to all orders in perturbation theory to stabilize the result. A similar observation was made in the unpolarized case, see Ref. [18].
The combined set of the next two plots, Figs. 5 and 6, gives a flavor of the remaining scale dependence of the single-inclusive distribution dg c 1 /dp T discussed in Fig. 3, which we deem to be the phenomenologically most important one. More specifically, Fig. 5 illustrates the dependence of 2x g c 1 /dp T obtained with the DSSV best fit (solid lines) at LO (top panel) and NLO (middle panel) accuracy under simultaneous variations (shaded bands) As can be seen, the variations for the best fit ∆g are sizable but still significantly smaller than the spread in dg c 1 /dp T introduced by choosing the min or max ∆g uncertainty sets, c.f. Fig. 3. Also the variations in the K-factor, given in the bottom panel of Fig. 5, are rather moderate, perhaps except for the vicinity of the node in the dg c 1 /dp T distribution. At small p T , near the peak of the distribution, the scale variations are significantly smaller then the NLO corrections, as one would hope for. At higher p T the width of the band is, however, rather similar in LO and NLO accuracy. The panels on the left-hand-side are for p T = 0.5 GeV and those on the right-hand-side correspond to setting p T = 5 GeV, i.e., we have picked two representative p Tvalues to the left and to the right of the node in dg c 1 /dp T , c.f. Fig. 3. Most importantly, the extreme min and max ∆g sets still provide the envelope of the DSSV predictions for dg c 1 /dp T in the entire range of µ. For any choice of scale, there is a clear, roughly constant sensitivity to ∆g as discussed above.
In the remainder of this section we turn to observables which involve the detection of both the heavy quark and the heavy antiquark and, hence, require the use of our MC code in their numerical evaluation. Again, we restrict ourselves to charm flavored (anti)quarks. Pre- sumably the most promising quantity is the invariant mass distribution 2x dg c 1 /dM of the HQ pair which is shown in Fig. 7 (left-hand-side) along with the corresponding double-spin asymmetry A c 1,M defined in analogy to Eq. (25). As our default choice of scale we now choose [27] T,Q Q , where p T,Q Q denotes the combined transverse momentum of the produced heavy quark and antiquark pair. In addition to the results for the longitudinally polarized structure function shown in Fig. 7, we have again successfully reproduced the results for the corresponding unpolarized quantities given in Ref. [27] to which we also refer the reader for a discussion of their relevant properties.
As before, we select x = 0.001 and Q 2 = 10 GeV 2 to define the DIS kinematics accessible at the EIC. The top and middle panels of Fig. 7 show the numerical results obtained at LO and NLO accuracy, respectively. The histograms (denoted by different symbols) refer to our usual selection of DSSV sets (best fit, min ∆g, and max ∆g), and the shaded bands correspond to the entire suite of DSSV uncertainty estimates. The lower panels give the respective K-factors for 2x dg c 1 /dM and A c 1,M . In general, we observe a rather similar behavior as for the single-inclusive p T -distributions shown in Fig. 3. All results exhibit a node at around M 7 GeV, which is accompanied by large QCD corrections in its vicinity. The asymmetries peak at around M 4 GeV. Again, initially a set of two measurements at the EIC, one to the left and one to the right of the predicted node, would reveal important new insights into the x-shape of ∆g. As before, all future measurements are required to resolve at least percent-level double-spin asymmetries to be of phenomenological relevance.
Other correlated observables of potential interest are collected in Fig. 8. The top panel shows g c 1 differential in the transverse momentum p T,Q Q of the HQ pair, the middle panel the correlation in its azimuthal angle 2x dg ference in pseudo-rapidity of the charm-anticharm pair, where η = 1 2 ln [(1 + cos(θ))/(1 − cos(θ))] with θ the angle of the HQ relative to the γ * p-axis in the c.m.s. frame. The histograms have been obtained with the MC code for the three sets of DSSV helicity densities (best fit, min ∆g, and max ∆g) and the shaded bands correspond to the range allowed by the DSSV uncertainty estimates. All results have been obtained at NLO accuracy for our standard choice of DIS kinematics, x = 0.001 and Q 2 = 10 GeV 2 , and setting the scale to µ 2 = 4m 2 +Q 2 +p 2 T,Q Q . In this case, we refrain from showing the corresponding double-spin asymmetries as the distributions in Fig. 8 are presumably initially of more limited phenomenological interest at the EIC than the ones presented so far. We note that we have, as usual, successfully reproduced all the respective unpolarized DIS results given in Ref. [24].
All three distributions given in Fig. 8 measure certain features of the DIS HQ structure functions which are present for the first time beyond LO accuracy. The distribution in the transverse momentum of the HQ pair (top panel) reveals the transverse momentum of the additional light-parton jet which recoils against the HQ pair. Since there is no such jet present at LO and in the virtual corrections at NLO accuracy, the p T,Q Q -distribution peaks at small values, and the first bin is dominated by "counter events" in the MC, see also Ref. [24]. Likewise, the distribution in the azimuthal angle ∆ϕ between the produced charm quark and antiquark, measured in the γ * p-c.m.s. frame, peaks at the lowest order back-to-back configuration, i.e., for ∆ϕ = π. The tails are produced by additional radiation of partons present for the first time at NLO accuracy. Similar remarks apply to the cone size variable ∆R shown in the bottom panel, where again the LO distribution is a delta function at ∆R = π.
Finally, we wish to comment on possible future additions to our suite of HQ electroproduction codes. First and foremost, our MC code can be expanded in several directions. Following the unpolarized HVQDIS code [28], the decay of the HQs can be modeled by including hadronization, i.e., additional convolutions with non-perturbative functions such as the charm-to-D-meson fragmentation function. If necessary, one could go even one step further and include also the decay of the heavy meson into experimentally observed electrons and muons; this was done, for instance, in the case of HQ polarized hadroproduction [12]. Another worthwhile addition would be to include also the lepton in the initial state in the MC which would enable one to study all the observables discussed in this paper directly in the laboratory rather than the γ * p c.m.s. frame; again, see the HVQDIS code [28] for an example. Secondly, in close collaboration with experimentalists, one can certainly optimize the HQ DIS observables presented in this paper by scanning for the best ranges in x, Q 2 , p T , etc. Presumably this should wait until the the EIC is further along its way, the machine parameters are confirmed, and experimental collaborations with final detector concepts are formed.

IV. SUMMARY AND OUTLOOK
In this paper we have presented a new parton-level Monte Carlo program at next-to-leading order accuracy in QCD that allows one to study heavy flavor production in longitudinally polarized deep-inelastic scattering in terms of any infrared-safe differential distribution for the structure function g Q 1 . The full heavy flavor mass dependence is retained throughout all calculations which makes this code particularly suited for phenomenological studies at the future Electron-Ion Collider which can explore charm electroproduction at small-to-medium momentum fractions and photon virtualities not much larger than the charm quark mass.
We have validated the Monte Carlo generator against our results for fully inclusive heavy quark production that where obtained previously with largely analytical methods. Within the latter framework it is also possible to compare to single-inclusive distributions differential in the transverse momentum or the rapidity of the observed heavy (anti)quark, albeit with no experimental cuts. In addition, we have verified known unpolarized results for heavy quark electroproduction that are available in the literature.
First phenomenological studies were carried out for various heavy quark distributions in the kinematic regime most relevant for the future Electron-Ion Collider. Particular emphasis was devoted to the expected size of the corresponding double-spin asymmetries and their sensitivity to the still poorly constrained helicity gluon distribution. Apart from single-inclusive distributions, we have also studied observables associated with the produced heavy quark pair including its invariant mass distribution and its correlation in azimuthal angle. Theoretical uncertainties associated with the choice of the factorization scale were studied using the important example of the single-inclusive transverse momentum distribution. Here, it was found that the inclusion of next-to-leading order QCD corrections reduces the scale dependence significantly.
An additional benefit of now having a flexible Monte Carlo generator for heavy quark polarized electroproduc-tion at hand is the possibility to implement and systematically explore relevant experimental cuts that might also help to enhance the sensitivity to the helicity gluon density further. So far, apart from the fully inclusive and transverse momentum differential charm contributions to polarized deep-inelastic scattering, the invariant mass spectrum appears to be the most promising observables in this context.
At present, the code does not include any heavy quark decays but they can be straightforwardly implemented if needed, following already existing frameworks and techniques in other codes. This can include even the decay of the heavy mesons into experimentally observed electrons and muons. The Monte Carlo generator can be also expanded in several other directions most worthwhile might be the inclusion of the lepton beam kinematics which would allow one to study heavy flavor electroproduction not only in the virtual photon-hadron frame but in the laboratory frame as well.
Finally, from a more theoretical perspective, the results presented in this and in our previous paper are also relevant for setting up a general-mass variable flavor number scheme for helicity dependent parton densities which might be more appropriate whenever heavy flavors are studied in some asymptotic region where the relevant scale, for instance, the virtuality of the photon in deep-inelastic scattering, is considerably larger than the heavy quark mass.
We plan to make the codes that were used in this and previous works publicly available.