Transverse single-spin asymmetry in the low-virtuality leptoproduction of open charm as a probe of the gluon Sivers function

We study the low-virtuality inclusive leptoproduction of open charm, $p^\uparrow l\rightarrow D^0+X$ as a probe of the gluon Sivers function. We perform the analysis in a generalised parton model framework. At leading order, this process is sensitive only to the gluon content of the proton. Hence any detection of a transverse single-spin asymmetry in this process would be clear indication of a non-zero gluon Sivers function (GSF). Considering COMPASS and a future Electron-Ion Collider (EIC), we present predictions for asymmetry using fits for the GSF available in literature. Predictions for peak asymmetry values lie in the range of 0.8\% to 13\%. We also present estimates of the upper bound on the asymmetry as obtained with a maximal gluon Sivers function. Further, for the case of the Electron-Ion Collider, we evaluate the asymmetry in the muons decaying from the $D$-meson and find that the asymmetry is well preserved in the kinematics of the muons. Peak values of the muon asymmetry are close to those obtained for the $D$-meson and lie in the range $0.75\%$ to 11\%.


I. INTRODUCTION
Transverse single-spin asymmmetries (SSA) can provide crucial information on the threedimensional structure of hadrons and have hence been a subject of increasing interest in the past two decades. While such asymmetries have been observed since the mid-70s in the hadroproduction of pions, i.e., p ↑ p → π + X [1][2][3], the past few years have provided a large amount of high quality data on SSAs in a wide variety of processes such as p ↑ p → π + X, p ↑ p → K ± + X, p ↑ p → J/ψ + X, lp ↑ → π + X, lp ↑ → K + X etc., (see Refs. [4,5] for reviews on the subject). A theoretical approach based on factorisation in terms of a hardpart and transverse momentum dependent (TMD) parton distribution functions (PDFs) and fragmentation function (FFs) has been formally established for processes which have two scales: a hard, high energy scale such as the virtuality of the photon in the Drell-Yan process and a relatively soft scale of the order of Λ QCD , such as the transverse momentum of the Drell-Yan lepton-pair. Another approach based on factorisation in terms of twist-3 parton correlators has been shown to be valid for the description of SSAs in processes with a single hard scale such as the transverse momentum of a pion in hadronic collisions.
Despite the absence of a formal proof, a lot of work has been done on a TMD description of single hard-scale processes under the assumption of factorisation, in what is generally referred to as the generalised parton model (GPM) approach. This approach has been quite succesful in describing unpolarised cross-sections in the hadroproduction of pions [6]. The leading-order (LO) GPM is able to describe (upto a K-factor) experimental data on pion production in high energy hadron-hadron collisions better than either the LO or the NLO collinear pQCD. It is also able to provide a good description of data on SSA in p ↑ p → π+X at widely different c.o.m energies [4]. One of the important transverse-momentum-dependent distributions which can lead to SSAs is the Sivers distribution [7,8], which encodes the correlation between the azimuthal anisotropy in the distribution of an unpolarised parton and the spin of its parent hadron. This anisotropy in the parton's transverse momentum distribution can lead to an azimuthal anisotropy in the distribution of the inclusive final state, i.e., a SSA. Fits of the u and d quark Sivers functions (QSF) obtained using data on A N (p ↑ p → π + X) at E704 ( √ s = 19.4 GeV), do a good job of describing the main features of the asymmetry observed at STAR ( √ s = 200 GeV) [4]. While the quark Sivers functions have been widely studied over the years, the gluon Sivers function (GSF) still remains poorly understood.
A first indirect estimate of the gluon Sivers function in a GPM framework was obtained in Ref. [9]. The analysis consisted of fitting the GSF to midrapidity data on SSA in pion production at RHIC. In this analysis, the quark contribution to the SSA was calculated using quark Sivers functions (QSF) as extracted from semi-inclusive deep inelastic scattering data.
The GSF fits obtained by this analysis predict asymmetries much smaller than allowed by the positivity bound on the GSF, viz. twice the unpolarised TMD gluon distribution. On the other hand, a recent study of large-p T hadron pair production in COMPASS indicates a substantial negative gluon Sivers asymmetry for both proton and deuteron targets [10].
Since large-p T hadron pairs are produced through the photon-gluon fusion, this process is indeed sensitive to the GSF. However the final state also receives contributions from the QCD Compton process, which is quark initiated. Hence an extraction of the GSF using large-p T hadron production is contaminated by the quark contributions to the SSA and would therefore depend on the extent to which these different processes can be separated in a data sample. With this being the first significant evidence for a non-zero GSF, it is important to study processes such as closed and open charm production which probe the gluon channel cleanly and directly.
In this work, using a GPM approach, we study the low-virtuality leptoproduction (Q 2 ≈ 0) of open-charm as a possible probe of the poorly understood gluon Sivers function (GSF).
At the leading-order (LO) of this process, the production of open-charm happens only via photon-gluon fusion (PGF), making it a direct probe of the gluon content of the proton. The GPM study of open-charm production as a probe of the gluon Sivers function was first proposed in Ref. [11] for the process p ↑ p → D 0 + X. In that study they considered two extreme scenarios for the GSF: zero and maximal. The term 'maximal' here refers to the Sivers function with its positivity bound of twice the unpolarised TMD using the fits of Ref. [9] and found that these fits predict sizeable, measurable asymmetries.
The low-virtuality leptoproduction of J/ψ has also been suggested as a probe of the GSF [13][14][15]. However, leptoproduction of open-charm may have some more advantages over the above mentioned processes. Firstly, unlike the case with p ↑ p → D 0 + X, one need not worry about possible factorisation breaking initial state interactions. p ↑ l → D 0 + X would have the same initial/final state interactions as SIDIS, for which TMD factorisation has been established. A study of SSA in this process might therefore complement studies of SSA in SIDIS and lp ↑ → h + X [16,17] by providing an additional handle on the gluon Sivers function. Secondly, open-charm production is free from dependence on production model, as is the case with closed-charm [18].
We therefore consider the process lp ↑ → D 0 + X in the low-virtuality regime in a GPM framework and see how it could serve as a probe of the GSF at both the COMPASS experiment and a future Electron-Ion Collider (EIC). While present data on open-charm production in COMPASS is limited due to statistics, the proposed Electron-Ion Collider [19] would have a significantly higher luminosity and should be able to provide better data on open/closed charm production. For both experiments, we present estimates for the maximum magnitude of SSA as obtained using the saturated GSF, and also the expected values of SSA obtained using the fits of Ref. [9].
In section II, we give the expressions for the relevant quantities in the GPM framework.
In section III, we give parametrisations for the different TMDs used, and in section IV, we discuss results for both the COMPASS and the EIC kinematics.

II. THE GPM FORMALISM
In this work, we are concerned with the single-spin asymmetry in the low-virtuality leptoproduction of open charm, where dσ ↑(↓) is the invariant differential cross-section for the process p ↑(↓) l → D + X with the spin of the transversely polarised proton being aligned in the ↑(↓) direction with respect to the production plane. Here, ↑ would be the +y direction in a frame where the polarised proton is moving along the +z direction and the meson is produced in the xz plane. Note that this is the convention for collider experiments (such as EIC). In case of fixed target experiments (such as COMPASS), by convention, the polarised proton would be considered to be moving along the −z direction with everything else remaining the same (cf. Fig. 1).
Following the treatment of open-charm hadroproduction [11], we can write the denomi-nator and numerator of Eq. 1 as, and In the above expressions, x g(γ) is the light-cone momentum fraction of the incoming gluon (photon) with the z-axis along the parent proton (lepton) direction, z = p + D /p + c is the lightcone momentum fraction of the D-meson with the z-axis along the fragmenting charm quark direction, k g(γ) is the intrinsic transverse momentum of the gluon (photon) with respect to the parent particle direction, k D is the transverse momentum with which the meson fragments from the charm quark,p is the unit vector along the heavy quark direction, m c is the charm quark mass, andŝ,t andû are the Mandelstam variables for the photon-gluon fusion process γg → cc.
∆ N f g/p ↑ (x, k ⊥ ) and f g/p (x, k ⊥ ) stand for the gluon Sivers function and unpolarised TMD respectively. f γ/l (x, k ⊥ ) is the transverse-momentum-dependent distribution of quasi-real photons in an unpolarised lepton, and D D/c (z, k D ) is the transverse-momentum-dependent fragmentation function. We will discuss the functional forms for all these distributions in Sec. III.
As mentioned earlier, the Sivers function, ∆ N f i/p ↑ (x, k ⊥ ; Q) describes the azimuthal anisotropy in the transverse momentum distribution of an unpolarised parton in transversely polarised hadron, and we have where k ⊥ = k ⊥ (cos φ ⊥ , sin φ ⊥ ). In a generalised parton model (GPM) description of this process, the only possible source of an asymmetry would be a non-zero gluon Sivers function.
Since photon-gluon fusion results in unpolarised final state quarks, there cannot be any contribution from the Collins effect, which would require transversely polarised final state quarks.
The partonic cross-section for photon-gluon fusion into a heavy quark pair is given by [20], where the Mandelstam variables are defined in the usual way, The factor C(x g , x γ , z, k D ) in Eqs. 2 and 3 contains the parton flux and the Jacobian relating the partonic phase-space to the mesonic phase-space. It is give by, The on-shell conditionŝ +t +û = 2m 2 c in Eqs. 2 and 3, gives a quartic equation in z. z can then be fixed by using this equation as shown in Ref. [12].
The delta function δ(k D ·p c ) in Eqs. 2 and 3 ensures that the region of integration for k D is confined to the two-dimensional plane perpendicular to the direction of the charm quark i.e., where k ⊥D represents values of transverse momenta on the allowed plane.
An outline of the treatment of the parton level kinematics and the TMD fragmentation is given in Appendix A.

III. PARAMETRISATION OF THE TMDS
Since we give predictions using the GSF fits of Ref. [9], for consistency we have to use the unpolarised gluon TMD and Sivers function used therein. We use standard factorised Gaussian form for the unpolarised gluon TMD, For the photon distribution f γ/l (x, k ⊥ ), we consider two cases: 1. The Weizsacker-Williams distribution of quasi-real photons with a Gaussian transverse momentum spread [13][14][15], where the Weizsacker-Williams distribution is given by [21][22][23], and the width of the Gaussian is k 2 ⊥γ = 0.1 GeV 2 .
2. The leading order result for the TMD distribution of photons in a lepton from Ref. [24], where m is the mass of the lepton.
The first choice, which we will refer to as Gaussian WW, was used in earlier studies of lowvirtuality leptoproduction by us [13][14][15] (and also in an analysis of low-Q 2 contributions to ep ↑ → h + X [16], but without the Gaussian spread) when first-principles result for the photon TMD distribution was not available. The second choice, which we will refer to as Photon TMD, is the first analytical result available in literature for the transversemomentum-dependent distribution of photons in a lepton [24]. Here we present results using both choices for completeness.
As with the unpolarised densities, we take the transverse-momentum-dependence of the FF to be Gaussian, with k 2 ⊥D = 0.25 GeV 2 . The gluon Sivers function is parametrised as follows [9], with 0 < ρ < 1. N g (x) parametrises the x-dependence of the GSF and is generally written SIDIS1 N g = 0.65 α g = 2.8 β g = 2.8 ρ = 0.687 It must obey |N g (x)| < 1 in order for the Sivers function to satisfy the positivity bound, In this work, for the predictions we consider two options for the gluon Sivers function: 1. the Sivers function with the positivity bound saturated, i.e., N g (x) = 1 and ρ = 2/3.
As mentioned in the introduction, we will refer to the first choice as the 'saturated' Sivers function. It would give an upper bound on the asymmetry for a fixed width k 2 ⊥ . The parameter ρ is set to 2/3 in order to maximize the first k ⊥ -moment of the Sivers function, following Ref. [25]. It must be kept in mind though, that this cannot be treated as giving an absolute upper bound on A N -an increased width k 2 ⊥ , for a fixed value of ρ, naturally would result in an increased asymmetry since the effects of the parton transverse momenta are more pronounced.
The SIDIS1 and SIDIS2 GSFs from Ref. [9] are the first (and so far, only) available extractions of the GSF in a GPM framework. They were obtained by fitting to PHENIX data on A N in inclusive pion production in the midrapidity region at RHIC. The QSFs used in these extractions, also labelled SIDIS1 and SIDIS2 respectively, were fit to data on semi-inclusive deep inelastic scattering. The SIDIS1 QSF set [26] (which was used in the extraction of the SIDIS1 GSF) was fitted to data on pion production in HERMES and positive hadron production in COMPASS with fragmentation functions by Kretzer [27]. It contains only the u and d quark Sivers functions since the data was not sensitive to sea quark contributions. The SIDIS2 QSF set [28] was fitted to flavour segregated data on pion and kaon production from HERMES and COMPASS and hence included sea quark Sivers functions as well. It used fragmentation functions by de Florian, Sassot and Stratmann (DSS) [29].
Both QSF sets give a good description of their respective SIDIS data sets. Furthermore both the GSFs (taken along with their associated QSF sets) describe the data on A N in midrapidity pion production equally well. Despite this the two fits show very different xdependencies, with SIDIS1 being larger in the moderate-x region and SIDIS2 being larger in the low-x region. The values of the parameters of the two GSF fits are given in Table I.

IV. RESULTS
In this section we present results on the unpolarised cross-section and SSA for COMPASS and EIC kinematics. Before going into the results, we should first make a note on the differing kinematic conventions of the two experiments: As COMPASS is a fixed target experiment, by convention the lepton is taken to be along the +z direction. This means that, in the definition of A N in Eq. 1, keeping the conventions for proton spin direction and production plane the same, positive x F and η correspond to the backward hemisphere of the proton, whereas negative x F and η correspond to the forward hemisphere of the proton. Note that this convention differs from that adopted Sec. II where, following the RHIC convention, the proton is taken to be moving along the +z direction. Since EIC, like RHIC, is also a collider experiment, we shall use the same convention for it. In the interest of clarity, the conventions used for the two experiments are illustrated in Fig. 1. Please note that since we are interested in quasi-real photoproduction, we have put a cut, Q 2 < 1 GeV 2 , where Q 2 = −(P l − P l ) 2 is the photon virtuality. This was motivated by the COMPASS antitagging cuts. In regions of large photon transverse momenta, the lepton-photon vertex becomes hard and the photon becomes off-shell. Hence one cannot use hard-parts defined for on-shell initial and final states. The cut on the photon virtuality Q 2 can be implemented by considering its relation to k ⊥γ and x γ , Q 2 = k 2 ⊥γ 1 + xγ 1−xγ . In this work, all results associated with both COMPASS as well as EIC were obtained with the Q 2 < 1 GeV 2 cut. When using the Gaussian WW approximation, this cut does not make a huge difference since the steeply falling k ⊥γ -dependence for k ⊥γ > k 2 ⊥γ prevents large contributions from regions of large virtuality. However, this is not the case with the Photon TMD as it has a much longer tail due to its 1/k 2 ⊥γ dependence. One must also note here that the opposite is true in the low-k ⊥γ region. At low k ⊥γ the Photon TMD falls off very sharply with increasing k ⊥γ whereas the Gaussian WW approximation, by virtue of being a Gaussian, has a flat k T -dependence at very low k T .
The numerical results were obtained using the GRV98LO set for the collinear gluon density and for the collinear part of the FF, the LO parametrisation of the c → D 0 fragmentation function by Kniehl and Kramer [30] was used. The QCD scale was chosen to be

A. COMPASS
The COMPASS experiment is a fixed target experiment involving a 160 GeV muon beam colliding on a proton target with a centre of mass energy with increasing P T in the range 0.5 < P T < 3.0 GeV. The cross-section obtained with the Photon TMD is smaller by roughly 30-40% over the entinre P T range. For a larger value of the width of unpolarised gluon TMD, viz. k 2 ⊥ = 1 GeV 2 instead of 0.25 GeV 2 , the crosssection at fixed P T is not affected much whereas, the cross-section at fixed pseudorapidity spreads out in P T somewhat -becoming smaller by 6% at P T = 0.5 GeV and larger by 40% at P T = 3 GeV -as one would expect. Overall, cross-section estimates for COMPASS are not very sensitive to the unpolarised gluon TMD width. Varying the width of the TMD FF in the range 0 < k 2 ⊥D < 0.25 GeV 2 also does not have any significant effect on the cross-section. Fig. 3 shows estimates for the maximal value of the magnitude of the asymmetry |A max N |, obtained by using the saturated gluon Sivers function viz., N g (x) = 1, ρ = 2/3. The results are presented as a function of x F at fixed P T = 1 GeV (left panel) and as a function of P T at a fixed pseudorapidity η = 1 (right panel). At fixed P T , estimates of |A max N | range from a minimum of about 12% at x F ≈ 0.2 − 0.3 to upto 24% at x F = 0.8. At fixed η, |A max N | shows a general increase with the meson transverse momentum, ranging from around 8% at P T = 0.5 GeV to 24% at P T = 3 GeV. Both the Gaussian WW distribution and the Photon TMD give similar results with the former being slightly smaller at low x F /P T and vice versa. Fig. 4 shows the asymmetries obtained using the SIDIS1 and SIDIS2 fits [9]. As was the  case with the saturated asymmetry, the results obtained with the Gaussian WW approximation and Photon TMD are generally similar, with the former being slightly smaller at low x F /P T and vice versa. Both fits give asymmetry predictions much smaller than allowed by the positivity bound with SIDIS2 giving significantly smaller asymmetries than SIDIS1. This is because the kinematic regions we consider probe the region 0.08 < x g < 0.5, where SIDIS2 is much smaller than SIDIS1, as can be seen from the numbers in Table I. At fixed P T , SIDIS1 gives a peak asymmetry of 4.2% at x F = 0 and SIDIS2 gives a peak asymmetry of 0.8% at x F = 0.8. At fixed η = 1, SIDIS1 gives a peak asymmetry of 7% and SIDIS2 gives a peak asymmetry of almost 1%, both at P T = 3.0 GeV. We have verified that changes in the width of the TMD FF in the range 0 < k 2 ⊥D < 0.25 GeV 2 do not alter the results for either SIDIS1 or SIDIS2 substantially and the general features of the A N predictions stay the same.

B. EIC
The Electron-Ion Collider (EIC) is a proposed experiment with colliding electron and proton/ion beams, with the possibility of both being polarised. It is meant to be capable of attaining high luminosities, with a centre of mass energy of upto 140 GeV in the ep configuration.
In Fig. 5, we show results for the unpolarised invariant cross-section using both the Gaussian WW approximation and the Photon TMD with the same Q 2 < 1 GeV 2 cut as used for COMPASS. We show the cross-section as a function of x F at fixed P T = 1.5 GeV (left panel) and as a function of P T at a fixed pseudorapidity η = 3 (right panel). At fixed P T , in the forward region, the cross-section decreases with increasing x F by more than six orders of magnitude in the range 0 < x F < 0.7. In contrast, in the backward region, the decrease in the cross-section with increasing |x F |, is only around one order of magnitude. This is because, for x F < 0 the gluon density is being probed in the small-x region and both the Weiszacker-Williams distribution and the Photon TMD are being probed in the moderateto-large-x region. In the small-x region the gluon density rises faster with decreasing x than both photon distributions, which behave as 1/x. Further, in the moderate-to-large-x region, the photon distributions fall much less steeply with increasing x than the gluon density.
These two effects combine to give the widely differing behaviour of the cross-section in the backward and forward regions. In the forward region (x F , η > 0), both the Photon TMD and the Gaussian WW approximation give almost identical results, whereas in the backward region the Photon TMD gives slightly smaller results for moderate values of negative x F . This is similar to what was observed at COMPASS. At fixed η, the cross-section decreases with increasing P T by four orders of magnitude in the range 0.5 < P T < 3.0 GeV. With a larger value of the unpolarised TMD width k 2 ⊥ = 1.0 GeV 2 , the cross-section at fixed P T is found to be unaffected in the forward region, but shows a decrease in the backward region of about 40% on average. The increase in the P T -spread of the cross-section is also observed at fixed η, but the effect is very small. The cross-section is also found to be insensitive to changes in the width of the fragmentation function in the range 0 < k 2 ⊥D < 0.25 GeV 2 . In Fig. 6, we show estimates for the maximal value of the magnitude of the asymmetry |A max N | obtained by using the saturated gluon Sivers function, as a function of x F at fixed P T = 1.5 GeV (left panel) and as a function of P T at fixed pseudorapidity η = 3 (right panel).
With the fairly large centre of mass energy of the EIC, we find that the general features of |A max N | are similar to what was obtained for proton-proton collisions at RHIC [11,12]. At fixed η = 3, for the Photon TMD, the asymmetry peaks at 21% at P T = 2 GeV. At fixed P T , large asymmetries are allowed in the forward region, with estimates being upto almost 25% at x F = 0.8. Overall, in the forward region (x F , η > 0) the Photon TMD gives results that are upto 18% larger than the what is obtained with the Gaussian WW approximation. This difference can be understood qualitatively, from the much smaller values of k ⊥γ contributing to production in the case of the Photon TMD and the resultant change in the values of x g , k ⊥g and x γ which contribute for a given P T and x F . As is the case for calculations at RHIC energy and kinematics [11,12], the asymmetry is suppressed in the backward hemisphere of the proton (x F < 0). This is because, in the backward region, the hard-part dσ/dt depends very weakly on the azimuthal angle of the gluon transverse momentum φ ⊥g . This weak dependence, along with the cos φ ⊥g term that is present in the Sivers function (see Eq. 4) leads to a suppression when the azimuthal angle is integrated over. The same has been observed in Ref. [11]. It must be mentioned however, that this feature is energy dependent and the suppression is weaker at lower centre of mass energies. This can be seen from the large values of |A max N | at x F 0.3 for COMPASS shown in Fig. 3. Fig. 7 shows the asymmetries obtained using the SIDIS1 and SIDIS2 fits [9]. As was the case for COMPASS kinematics, both fits give asymmetries much smaller than allowed by the positivity bound, with SIDIS1 giving the larger results of the two. As was found at COMPASS energy, the Photon TMD gives results that are a few percent larger than those obtained using the Gaussian WW approximation. At fixed P T , in the forward region, SIDIS1 gives a peak asymmetry of 13% at x F = 0.4 and SIDIS2 gives a peak asymmetry of 0.8% at x F = 0.3. In the backward region x F < 0, D production gets contributions mainly from from x g < 0.08, where SIDIS2 is larger than SIDIS1. However the overall values of both fits in this region are very small. Combined with the azimuthal suppression, this makes the asymmetries from both fits almost negligible. At fixed pseudorapidity, SIDIS1 gives a peak asymmetry of around 11.5% at P T = 2.5 GeV and SIDIS2 gives a peak asymmetry of around 0.8% at P T = 1.6-2.2 GeV.

C. Single-Spin Asymmetry in open-charm decay muons
So far we have considered the SSA in terms of the D-meson kinematics. It would also be interesting to consider the SSA in terms of the kinematics of the decay muons. A detector such as the proposed ePHENIX [31] would be able to study open heavy flavour production through the leptonic decay channels. With this in mind, we consider the semileptonic decay of the D's in order to obtain the SSA for the decay muons, A µ N , where dσ P l→D+X→µ+X is the Lorentz-invariant inclusive decay-muon cross-section, In keeping with the conventions for the D-meson asymmetry defined in Eq 1, we take the muon to be produced in the xz plane with the proton moving along the +z direction and its spin parallel or antiparallel to the y-axis. Using the narrow width approximation, the expression for the Lorentz-invariant decay-muon cross-section can be written for a general n-body decay channel as follows: where the x i are the n − 1 decay products produced along with the muon. Γ total is the total decay width of the D-meson and BR stands for the branching ratio for the considered n-body decay channel. The above expression consists of the meson production cross-section, the decay matrix element and a phase-space integral over all the decay products except the muon. It makes use of the fact that the decay of a scalar meson can be treated as independent of its production, allowing a factorised form involving the meson cross-section convoluted with the differential decay rate. This can be shown to be true by using the narrow width approximation.
To account for all possible open-charmed meson decays in to muons through all possible channels would be a complex task. To simplify things, we make the following assumptions: First, we consider only the decay of the D 0 to muons. Muons can also be produced through the decay of other charmed mesons states but we do not consider those decays here. For the D 0 we consider the two major semileptonic decay channels, D 0 → K − µ + ν µ which as a branching ratio of 3.33% and D 0 → K * (892) − µ + ν µ which has a branching ratio of 1.92%. Second, in the calculation of the three-body decays, we set the decay matrix elements |M D 0 →K − µ +ν µ | and |M D 0 →K * − µ +ν µ | to 1, and only account for the phase-space kinematics.
In Eq. 19, the momentum p D must be integrated over the entire region of phase-space from which a D-meson can decay to produce a muon of given momentum p µ . The derivation of the closed form expression for decay-muon invariant cross-section, E µ dσ P l→D+X→µ+X d 3 p µ and the integration limits for the momentum of the D-meson, p D is given in Appendix B.
The results for the decay-muon invariant cross-section and SSA, A µ N are presented in Fig. 8. The asymmetry is shown for the case of the saturated GSF and the SIDIS1 and SIDIS2 [9] fits. A µ N is presented as a function of x F µ ≡ 2P Lµ / √ s, with the muon transverse momentum P T µ = 1.5 GeV. It appears that an azimuthal anisotropy in D production would be retained significantly in the decay-muons. The general dependence of the A µ N on x F µ is similar to the dependence of the D-meson SSA on x F . As with the D-meson, the muon SSA is also suppressed in the backward hemisphere. Peak values of the A µ N are close to those obtained for the meson: With the Gaussian WW approximation and SIDIS1 GSF, A µ N has a peak value of 11% at x F µ = 0.3 whereas A N has a peak value of almost 12% at x F = 0.4.
With the SIDIS2 GSF, A µ N has a peak value of 0.75% at x F µ = 0.23 whereas A N has a peak value of 0.8% at x F = 0.3.

V. CONCLUSIONS
In this work, we have presented results for SSA in the low-virtuality leptoproduction of open-charm at both COMPASS and a future Electron-Ion Collider. We find that an asymmetry of upto around 25% is allowed by the saturated gluon Sivers function at both COMPASS and EIC. We also find that, for EIC kinematics, the asymmetry is significantly retained in the distribution of the decay muons. In calculating the asymmetry we used two different forms for the TMD distribution of quasi-real photons in the lepton. The first was the Weizsacker-Williams distribution with a Gaussian transverse-momentum spread (Gaussian WW) and the second was the LO analytical result for the TMD distribution of photons in a lepton (Photon TMD) from Ref. [24]. At COMPASS energy the two forms give similar results, whereas at EIC energy, the Photon TMD gives slightly larger asymmetries in the forward region. The differences in the result for the two distributions can be attributed to the interplay of different k ⊥γ , x γ , x g and and k ⊥g values that get sampled in the two cases for a given value of x F /η and P T .
The two GSF fits of Ref. [9] for which we give predictions, were extracted from data on midrapidity pion production at RHIC. As mentioned earlier, the two differ in the flavour structure of the QSFs used as well as the light quark fragmentation functions used in the extraction. SIDIS1 was obtained using an extraction of the QSFs [26] that included only the u and d flavours and used fragmentation functions by Kretzer [27]. SIDIS2 was obtained using an extraction of the QSFs that also included sea quarks [28] and used more recent fragmentation functions by de Florian, Sassot and Stratmann [29]. While both the GSF fits, taken along with their associated QSF sets, describe the input data on SSA in midrapidity pion production equally well, they have widely differing x-dependencies. This indicates that pion production in the midrapidity region at RHIC is only weakly sensitive to the gluon Sivers function.
In this work, we find that the low-virtuality leptoproduction of open-charm, which probes the gluon content of the proton directly, is able to discriminate well between these two fits. Thus we see that this process offers a good probe of the gluon Sivers function and In this work, we have considered inclusive single-particle leptoproduction in the low virtuality regime. This allows us to handle the process in terms of a TMD distribution of quasi-real photons in a lepton, not unlike a TMD distribution of partons in a hadron. The treatment of parton kinematics here is thus similar to the treatment of transverse-momentum-dependent parton kinematics for inclusive single-particle hadroproduction, which can be found in quite a few places [32,33] including Ref. [6], where heavy meson final states have been considered.
The momenta of the proton, lepton and the D-meson can be written in the p-l centre of mass frame as, where the masses of the proton and lepton have been neglected.
The gluon and the quasi-real photon carry light-cone momentum fractions x g = P + g /P + P , x γ = P − γ /P − l and transverse momenta k g and k γ respectively. Their momenta are given by, where φ ⊥g and φ ⊥γ are the azimuthal angles of gluon and photon transverse momenta.
The heavy quark is produced through photon-gluon fusion gγ → cc and then fragments into the heavy meson. The momentum of the heavy quark is described by z, the light-cone momentum fraction of the heavy meson and k D , the transverse momentum of the meson with respect to direction of heavy quark. In a choice of coordinates where the heavy quark momentum, p c is along the z-axis, the D-meson momentum can be written as where the first term on the right is the component along the heavy quark direction and the second term is the component transverse to it. Here, k D is simply (k Dx , k Dy , 0) = (k D ⊥ , 0).
In the lab coordinates however, k D can have all three components non-zero and is specified as, and the orthogonality condition k D .p c = 0 ensures that k D lies in a plane perpendicular to p c . The light-cone momentum fraction z is given by, This gives us the expression for the energy of the heavy quark, The expression for p c can be obtained from the fact that it is collinear with p D − k D and that the unit vector constructed out of both must therefore be equal, Eqs. 24 and 25 relate the energy and momentum of the observed D-meson with that of the fragmenting heavy quark for given values of k D and z.
The term d 3 k D δ(k D ·p c ) in Eqs. (2) and (3) ensures that the k D integration is only over momenta transverse to the fragmenting parton: where, Limits on k D can be obtained by requiring | cos φ 1 | ≤ 1, Furthermore, the inclusion of intrinsic transverse momenta in the kinematics calls for the following constraints: a) the energy of the incoming parton should not be greater than that of its parent particle, E g(γ) ≤ E p(l) and, b) the energy of the D-meson should not be greater than the energy of the heavy quark E D ≤ E c . The first leads to the following bound on the transverse momenta of the incoming partons, k ⊥g(γ) < √ s min[x g(γ) , x g(γ) (1 − x g(γ) )].
The second constraint, E D ≤ E c is inherently fulfilled by Eq. 24. However, this alone does not ensure that the heavy quark is more energetic than the D-meson in the photon-gluon c.o.m frame. By demanding E c > E D in the γ-g c.o.m frame, we get a lower bound onŝ, s ≥ 2P D .(P γ + P g ).
In our earlier work on open charm production (Ref. [12]) we had not implemented this bound in our calculations. We find that the inclusion of this bound significantly improves the convergence of the integral in the close vincinity of x F = 0.

B. Derivation of Lorentz-invariant decay-muon cross-section
With the decay matrix element set to unity, the three body decay width can be written as follows, where s = (p x 1 + p x 2 ) 2 = (p D − p µ ) 2 is a Lorentz-invariant quantity. Here, x 1 and x 2 are the decay products produced along with the µ. We take x 1 to be K − or K * − and x 2 to bē ν µ . Here we consider the muon to be massless. Since we know the decay factorises, we can write the production and decay as a convolution, where dΓ is the decay-width for an infinitesimal muon momentum region d 3 p µ . Then using Eq. 31 we finally have, The allowed phase-space region for p D can be obtained determined by considering the decay in the rest frame of the D-meson. In it, one can see that the allowed values of muon energy lie in the range 0 < E D com µ < (m 2 D − m 2 K )/2m D , where m K is the mass of the kaon. This constraint on the muon energy in the rest frame of the D, can be translated into the following constraint on the the Lorentz-invariant quantity constructed from the D-meson and muon four-momenta: The expression for s in terms of the momenta involved is, where we have assumed that the muon is massless and is in the xz-plane, P T µ = (P T µ , 0).
Here P T µ and P Lµ are the x and z components of the muon momentum respectively. The lower inequality in Eq. 34 can be cast as follows: We will call the quantity on the right hand side of the above expression, Y . This gives us a constraint on the angle of the D-meson transverse momentum, Naturally, we also require Y ≤ 1, since the lower bound on a cosine term can't be greater than 1. Demanding this gives us upper and lower limits on P L : Demanding Y ≤ 1 also gives us a lower bound on P T :