Accessing the Gluon Sivers Function at a future Electron-Ion Collider

In this work, we present a systematic study on the feasibility of probing the largely unexplored transverse momentum dependent gluon Sivers function (GSF) in open charm production, high p$_T$ charged di-hadron and di-jet production at a future high energy, high luminosity Electron-Ion Collider (EIC). The Sivers function is a measure for the anisotropy of the parton distributions in momentum space inside a transversely polarized nucleon. It is proposed that it can be studied through single spin asymmetries in the photon-gluon fusion subprocess in electron proton collisions at the EIC. Using a well tuned Monte Carlo model for deep inelastic scattering, we estimate the possible constraints of the GSF from the future EIC data. A comparison of all the accessible measurements illustrates that the di-jet channel is the most promising way to constrain the magnitude of the GSF over a wide kinematic range.

In this work, we present a systematic study on the feasibility of probing the largely unexplored transverse momentum dependent gluon Sivers function (GSF) in open charm production, high pT charged di-hadron and di-jet production at a future high energy, high luminosity Electron-Ion Collider (EIC). The Sivers function is a measure for the anisotropy of the parton distributions in momentum space inside a transversely polarized nucleon. It is proposed that it can be studied through single spin asymmetries in the photon-gluon fusion subprocess in electron proton collisions at the EIC. Using a well tuned Monte Carlo model for deep inelastic scattering, we estimate the possible constraints of the GSF from the future EIC data. A comparison of all the accessible measurements illustrates that the di-jet channel is the most promising way to constrain the magnitude of the GSF over a wide kinematic range.

I. INTRODUCTION
In recent years, an important forefront in hadron physics is to explore the 2+1 dimensional partonic structure of nucleons by including information on the internal parton transverse momentum and coordinate space distributions. These extensions have significantly broadened our understanding to the nucleon structure compared to the 1d picture in the longitudinal momentum space. The transverse momentum structure of nucleons is encoded in the transverse momentum dependent (TMD) [1] parton distribution functions (TMDs), which contain information on both the longitudinal momentum fraction x and the transverse (sometimes called intrinsic) motion k ⊥ of quarks and gluons inside a fast moving nucleon.
When including spin degrees of freedom, TMDs link information on the intrinsic spin of a parton (s q,g ) and their transverse motion (k q,g ⊥ ) to the spin direction of the parent nucleon. At leading twist the most general spin dependent TMD can be denoted by f q,g 1 (x, k q,g ⊥ ; s q,g , S). At leading order, there are eight such combinations, yielding eight independent TMDs [2]. The Sivers function f ⊥ 1T [3], which encapsulates the correlations between a parton's transverse momentum inside the proton and the spin of the proton, has received the widest attention both phenomenologically and experimentally among all TMDs. It * zhengliang@cug.edu.cn † elke@bnl.gov ‡ jhlee@bnl.gov § bxiao@mail.ccnu.edu.cn ¶ zbyin@mail.ccnu.edu.cn was found that the Sivers function is not universal in hard-scattering processes [4], which has its physical origin in the rescattering of a parton in the color field of the remnant of the polarized proton [5]. Proving experimentally the process dependence of the Sivers function is a very important test of the non-Abelian nature of quantum chromo-dynamics (QCD) in TMD factorization.
Experimentally, the quark Sivers function has been measured in semi-inclusive deep inelastic scattering (SIDIS) by the HERMES, COMPASS, and JLab Hall A collaborations [6][7][8]. However, due to the limited statistics precision and the narrow kinematic coverage of the SIDIS data, only the valence quark Sivers function at moderate to high x could be constrained in phenomenological fits [9]. The quark Sivers function has also been studied in polarized proton-proton collisions by the STAR and PHENIX collaborations at the Relativistic Heavy Ion Collider (RHIC) [10,11]. There are first indications both from STAR through the W -boson measurement [12] and COMPASS in Drell-Yan (DY) production [13] for the non-university of the Sivers function [5] if measured in hadron-hadron collisions or SIDIS, but the still challenging statistical precision limits a definite conclusion. Both STAR and COMPASS will soon increase the statistical precision of these measurements by including recent high statistics data. At the future high energy, high luminosity Electron-Ion Collider (EIC) [14], the quark Sivers function can be well constrained over a very wide kinematic range (x, Q 2 , z and p T ) in SIDIS with exquisite precisions. It has been systematically investigated in the one and two hadron final states in Ref. [15] with a modified PYTHIA event generator that includes the quark Sivers effect. is barely known at the current stage [16]. Presently, the major theoretical constraint for the GSF comes from the Burkardt sum rule [17], which requires the total transverse momentum of all partons in a transversely polarized nucleon to vanish. The only direct constrain of the GSF comes from the left-right asymmetry A N data in p ↑ p → π 0 X within the so-called TMD generalized parton model (GPM) framework [18]. This analysis found that the gluon Sivers function is not large [19]. However, the gluon Sivers obtained in the GPM may differ from the gluon Sivers function in the TMD framework [16]. At this moment the only experimental constraint to the gluon Sivers function in the TMD framework comes from the recent SIDIS measurement of high-p T hadron pairs off transversely polarized deuterons and protons at COM-PASS [20]. This analysis found that the gluon Sivers asymmetry is negative at large x B within statistical uncertainties. Interestingly, this finding is in qualitative agreement with results from the calculation based on the light-cone spectator model [21].
Accounting for the different gauge link structures involved in deep-inelastic scattering (DIS) and hadronic collisions, the gluon Sivers function is expected to be process dependent. A test of this non-universality of the gluon Sivers function is of equal importance as for the quarks and is currently not validated. Similar to the sign change of the quark Sivers function, the gluon Sivers function accessed in ep ↑ → e ccX is also predicted to be related to that in p ↑ p → γγX by an overall sign change as shown in [22]. Consequently, the study of the gluon Sivers function at an EIC will provide a unique test of the fundamental non-perturbative QCD effects through complementary information to the proposed gluon Sivers function observables at RHIC and LHC [23,24]. In addition, as pointed out in Ref. [25], there are two different type of gluon TMDs, namely, the Weizsäcker-Williams and the dipole gluon distribution. This is a direct consequence of the different gauge link dependences. By comparing the gluon Sivers functions extracted from DIS and pp collisions, one can test this gauge link dependence since the Weizsäcker-Williams and dipole type T-odd gluon TMDs are expected to behave differently [16,22].
In DIS, the key to study the gluon Sivers function is to tag the Photon-Gluon Fusion (PGF) subprocess. It has been shown in Ref. [25,26] that the gluon transverse momentum distribution can be mapped through quark-antiquark jet correlations for the PGF subprocess γ * g → qq. Ref. [22] suggests that the spin asymmetries measured in heavy quark pair and di-jet production at an EIC can be used to study the Weizsäcker-Williams gluon TMDs including the Sivers function. The open charm production in electron-proton scattering ep ↑ → e ccX is argued to be an ideal probe to tag the PGF process, and can be investigated at a future EIC. A model study has been carried out in [2] and the related experimental considerations for tagging charm quark production through D-mesons in the final state for the PGF subprocess are discussed in [27]. In this paper, we will provide detailed information on EIC projections for open charm production with attainable experimental conditions. Alternative methods tagging the gluon channel through the production of high-p T hadron pairs and di-jets are also studied. The advantages and disadvantages of the different channels will be discussed. Table I shows the definitions of the kinematic variables used in this paper. Virtuality of exchanged photon xB Bjorken x y Energy fraction of virtual photon with respect to incoming electron W Center of mass energy of the γ * p system x Longitudinal momentum fraction of the quark/gluon from the polarized proton involved in the hard interaction z h,q Energy fraction of of a hadron or quark with respect to virtual photon in target rest frame k ⊥ Initial transverse momentum of gluons inside the proton in γ * p center of mass frame k 1⊥ , k 2⊥ Transverse momentum of the two outgoing partons in γ * p center of mass frame p h1⊥ , p h2⊥ Transverse momentum of the trigger/associate particle in γ * p center of mass frame pT Transverse momentum of final state hadron/jet with respect to virtual photon η Pseudorapidity of final state hadron/jet in γ * p center of mass frame PT Transverse momentum scale of final state particle/jet pair with respect to virtual photon kT Vector sum of the transverse momentum for the final state hadron/jet pair in the final state φ kS Sivers angle, the azimuthal angle difference of kT and the proton spin direction p T Lab , p Lab Transverse momentum/momentum of final state hadron in the laboratory frame η Lab Pseudorapidity of the final state hadron/jet in the laboratory frame The remainder of the paper is organized as follows: in Sec. II we discuss the theoretical framework used to build the connection of gluon Sivers function and the size of the single spin asymmetry (SSA). The Monte Carlo setup is described in Sec. III. A detailed description of the results and their projected precision are presented in Sec. IV. We summarize in Sec. V.

II. SINGLE SPIN ASYMMETRY (SSA) ARISING FROM THE GLUON SIVERS EFFECT
The Sivers function describes the distribution of unpolarized partons with flavor a inside a transversely polarized proton with mass M p and can be expressed following the Trento convention in Ref. [28] as: The first term represents the axially symmetric contribution from the unpolarized parton distribution, while the second term ∆ N f a/p ↑ (x, k ⊥ ) generates a distortion away from the center in the number density of unpolarized partons with an intrinsic transverse momentum k ⊥ . The azimuthal dependence of this distortion is given by S · (ˆ P ×ˆ k ⊥ ), where P and S are the polarized proton three momentum and spin polarization vector, respectively. The notation ∆ N f a/p ↑ (x, k ⊥ ) is related to the Sivers function denoted as f ⊥a [29].
The production of high transverse momentum charged hadron pairs or di-jets in DIS through γ * g → qq is dominated by gluons, although it may also have some contribution from the quark channel depending on the process measured. The cross section can be calculated in an effective k t factorization framework at leading order as shown in Ref. [30]. If k 1 and k 2 are the four momenta of the outgoing quarks, one can obtain the dihadron cross section as a generalization of the unpolarized case [31] with the transverse momentum imbalance defined as k ⊥ = | k 1⊥ + k 2⊥ | and the transverse momentum scale as where z q is the momentum fraction of produced quark q with respect to the incoming virtual photon and H γ * g→qq tot gives the combined hard factor that incorpo- (ŝ+Q 2 ) 4 (û t +t u ) of the virtual photon. Eq. 2 can be further simplified using the condition k ⊥ P ⊥ known as the correlation limit [30]. Eq. 2 can thus be expressed as in which 2 f is related to Q 2 as 2 f = z q (1 − z q )Q 2 . Choosing the center of mass frame of the exchanged virtual photon and the proton, in which the proton beam with momentum P is moving in the −z direction, one can obtain an explicit form of the mixed vector product in Eq. 1 as S · (ˆ P ×ˆ k ⊥ ) = sin(φ k − φ S ) with φ k being the azimuthal angle of k ⊥ . A factorized Gaussian parameterization has been adopted for the transverse momentum dependent unpolarized parton distribution function and fragmentation func- There exists a strong correlation between the kinematics of the back-to-back hadron pair and its parent quarks. Therefore, one can use the following variables measurable hadron level P T = | p h1⊥ − p h2⊥ |/2 and k T = | p h1⊥ + p h2⊥ | to access the underlying parton kinematic variables P ⊥ and k ⊥ . A schematic illustration of the encoded kinematic variables is shown in Fig. 1. The GSF can be studied in the single spin asymmetry (SSA) for di-hadron production as follows: where subscript "U" represents the unpolarized electron beam and "T" indicates the transverse polarization of the proton beam. φ kS = φ k T − φ S stands for the angular difference between the total di-hadron transverse momentum k T and the polarized proton spin direction S ⊥ . The amplitude of the SSA is proportional to the corresponding Sivers function divided by the unpolarized parton distributions.

III. MONTE CARLO SIMULATION SETUP
In this section, we will describe the setup for our event generation. We use the PYTHIA-6.4 Monte Carlo (MC) program [32] to simulate the unpolarized cross section as expected at an EIC. The PYTHIA generator has been found to reproduce the charged and open charm particle production in the electron proton collisions at HERA. The comparison of the HERA data [33,34] and the output of the tuned PYTHIA MC for charged particles and D * mesons is shown in Fig. 2 and Fig. 3. Based on this reasonable description of the unpolarized DIS cross section, we will discuss our strategy to obtain the SSA based on weighting the unpolarized results from PYTHIA.
In the simulation, we model the amplitude of the asymmetry as an incoherent superposition of all contributing subprocess on the event-by-event basis. For every event, a weighting factor is obtained according to the kinematics and parton flavor as follows: At the end, the Monte Carlo asymmetry can be understood as the weighted sum of the asymmetry weights from signal (gluon initiated channels) and background (quark initiated channels) processes similar to the strategy used in Ref. [35]: in which N g and N q indicate the number of gluon and quark initiated events in the analyzed event sample. The corresponding event fraction is thus obtained as (Color online) Charged particle transverse momentum distributions for 0 < η < 1.5 defined in the virtual photon-hadron center of mass frame. The HERA data [33] for 5 GeV 2 < Q 2 < 10 GeV 2 , 0.0005 < xB < 0.002 with a beam energy 27.6 GeV × 920 GeV are compared to the tuned PYTHIA results.  [34] for 5 GeV 2 < Q 2 < 100 GeV 2 , 0.02 < y < 0.7 with the beam energy 27.6 GeV × 920 GeV are compared to the tuned PYTHIA results.
In the experiment, it is very hard to reliably separate different subprocesses. Therefore, the fractions of events from different subprocesses are modeled using PYTHIA in this analysis.
The parameterization of the Sivers function is given in a factorized form as in which f a/p (x a , Q 2 ) is the unpolarized parton distribu- describe the x and k ⊥ dependence of the Sivers function. The magnitude of the asymmetry from background contributions is calculated from the quark Sivers function from these recent fits [36] listed as: For the gluon Sivers function we utilize two models as input to our study. The first model is the SIDIS1 set obtained in the fit [18], which follows a similar parameterization form as the quark Sivers function with the parameters given by: The second gluon Sivers model relies on the positivity bound assumption used in [37]: in which the positivity limit is saturated when k ⊥ = 0.8 GeV. We will use 10% and 5% of the positivity bound to study quantitatively the measurability of the gluon Sivers function. We calculate the weight of every event according to the inputs discussed here to obtain the magnitude of the asymmetry in the final state. An example of the first k ⊥ moment of the input Sivers distribu- Fig. 4. Fig. 4(a) shows the quark Sivers functions used to estimate the background contribution, while the gluon Sivers functions are shown in Fig. 4(b). For the current parameterizations the quark Sivers functions are maximum for x >0.1 for the valence quarks and become negligible in the small x regime. The magnitude of the seaquark sivers functions is small over the entire x-range. It is noticeable that the gluon Sivers functions based on the positivity bound assumption and SIDIS1 set have quite a different functional form in x.
We provide in Fig. 5 a comparison of the charged hadron asymmetry measured by the COMPASS experiment [38] with the asymmetry obtained from weighting PYTHIA events according to the method described above with the quark Sivers functions. It is not surprising to see that radiation effects modeled by the parton shower mechanism in PYTHIA are quite weak at the COMPASS energy. The comparison also shows that one can describe both positive and negative charged hadron asymmetries from COMPASS with the event weighting method.
It should be explicitly noted that the parameterizations of the Sivers asymmetry discussed here are not accounting for any effects due to the QCD scale dependence of TMDs. The QCD evolution of the Sivers function can be calculated in the QCD resummation formalism following the Collins-Soper-Sterman method [39,40] by applying the correct Sudakov factor to the spin dependent parton distributions [41][42][43]. The needed precise phenomenological inputs to determine the QCD scale evolution of Sivers asymmetry is not yet available, while they can be obtained from the future RHIC and EIC measurements. We will therefore not address the evolution of the gluon Sivers function in this paper but leave it for future work.
In order to estimate the statistical uncertainty of the SSA in our simulation, we use δA = 1 P 2 N − A 2 N from Ref. [44], where N represents the count of selected pairs in a certain kinematic bin, and P indicates the polarization of the proton beam. In this work, we assume a polarization P = 70% for the EIC beam energy configuration of 20 GeV × 250 GeV with an integrated luminosity L int = 10 fb −1 . Ref. [38]. Radiation effects are estimated by turning on (WithPS) and off (NoPS) parton shower mechanism in PYTHIA shown with dotted and solid lines in this comparison.

IV. RESULTS AND DISCUSSIONS
In this study, the event kinematics has been restricted to 0.01 < y < 0.95 and 1 GeV 2 < Q 2 < 20 GeV 2 with the electron and proton beam energy configuration of 20 GeV × 250 GeV. A detector system especially designed for EIC with a wide acceptance −4.5 < η Lab < 4.5 for measuring charged particles [45] has been assumed, in which case the event kinematics can be well reconstructed from the scattered electron. This selection gives an average event kinematics for the minimum bias events as x B = 0.0012, Q 2 = 2.5 GeV 2 , W = 54.6 GeV. The wide kinematic reach at this high center-of-mass energy ( √ s=141 GeV) makes it possible to study the evolution of the Sivers asymmetry and to access the region dominated by gluons.
As discussed in Sec. II, there is a correlation between the vector sum of the transverse momentum k T for the selected hadron pairs or di-jets with the initial transverse motion of gluons. It then follows naturally to investigate the Sivers asymmetry through the sine-modulation for the angle φ kS = φ k T − φ S , which defines the difference between k T and the spin direction of the proton. To tag the gluon distributions and study the Sivers asymmetries we will study the D meson pair, charged di-hadrons and di-jet production. The detailed experimental cuts for each channel are listed as follows: • D meson pair production D 0 → πK, |η π/K Lab | < 3.5, p π/K T Lab > 0.2 GeV, p D T > 0.7 GeV, z D > 0.1 • charged di-hadron production −4.5 < η Lab < 4.5, p h T > 1.4 GeV, z h > 0.1 • di-jet production π, K, p, γ with p T Lab > 0.25 GeV, |η Lab | < 4.5 used for the jet reconstruction with anti-k T algorithm and a cone radius R = 1 ,the trigger jet has p jet1 T > 4.5 GeV and the associate jet p jet2 T > 4 GeV

A. The Gluon SSA in Open Charm Production
The heavy flavor production has been proven to be very useful for measuring the gluon Sivers asymmetry in proton-proton collisions as shown in Ref. [37,46]. Similar to the case in hadron-hadron reactions, it is well accepted that the heavy flavor production in DIS is a very clean channel to directly probe the gluon distributions. In this section we demonstrate the possibility to measure the gluon Sivers function in open charm production γ * g → cc with D 0 -mesons in the final state. Open charm production has the advantage that quark initiated process are strongly suppressed and one becomes essentially only sensitive to gluon initiated subprocesses. The pseudorapidity distribution of the K-meson from the D 0 decay can be found in Fig. 6. The K momenta are typically a few GeV in the central rapidity region, and extend to 10 GeV at rapidities |η| > 1. The distribution for the π-mesons from the D 0 decay is found to be very similar to the one from Ks. The D 0 meson decay products are required to be in |η π/K Lab | < 3.5 and to have transverse momenta p π/K T Lab > 0.2 GeV to be reconstructed and identified. The p T Lab correlation of the D 0 meson decay products is shown in Fig. 7, most of the π, K products pass the transverse momentum cut. The kinematics of a directly produced D 0 -meson and for one from the decay of a heavier charm-mesons is basically the same, therefore all D 0 -mesons with p T > 0.7 GeV and z h > 0.1 are included in this study. To capture the full charm anti-charm quark pair kinematics, we select events with DD pairs in the final state. Fig. 8 shows that gluon initiated processes account for about 90% of the total selected events over a wide range in x B for two Q 2 bins. For x B > 0.1 quark initiated subprocesses become slightly more important. The sensitivity of the DD pair measurement to the magnitude of the gluon Sivers function is shown in Fig. 9(a). The statistical uncertainty is based on an integrated luminosity of L int = 10 fb −1 . The solid curve represents the parton level asymmetry. The Sivers asymmetry based on the scenario with 10% of the positivity bound of the gluon Sivers function is indicated by the black filled symbols. The SSA for the background quark initiated Sivers effect is consistent with zero and not shown here. The limited statistical precision for the DD final state due to the small branching ratio (3.87% D → Kπ) makes it challenging to precisely determine the gluon Sivers function on the level of 10% of the positiv-ity bound. Therefore, we also investigated the sensitivity to the magnitude of the gluon Sivers function requiring only one D meson. The Sivers angle φ kS is calculated replacing k T with the D meson transverse momentum. Fig. 9(b) depicts the SSA based on the 10% of the gluon Sivers positivity bound assumption, which can be well distinguished from the background SSA due to quark Sivers effects. Comparing Fig. 9(b) and Fig. 9(a), one can observe that the initial parton level asymmetries are similar, but the magnitude of the final state asymmetry for single D mesons is reduced since the transverse momentum of one D meson is not a good proxy for the initial gluon transverse momentum. A similar approach to the open charm production is to select K + K − pairs in the final state, which enhances as underlying process γ * g → ss. We find in our study this measurement is also statistically limited and can only resolve a gluon Sivers signal up to 10% of the positivity bound. Since the global features of this measurement are similar to the di-hadron case, which will be discussed in the next section, we will not provide more information.

B. Gluon SSA through Charged Di-hadron Pairs
In SIDIS production, the Leading Order DIS (LODIS) process γ * q → q is accounting for a large fraction of the charged particle productions. The LODIS process can be largely suppressed by requiring a pair of high p T charged hadrons.
The acceptance for the charged hadrons is required to fit the EIC detector design |η Lab | < 4.5 and the hadron pair candidates must have a transverse momentum p h T > 1.4 GeV and z h > 0.1. To select hadron pairs from back-to-back jets, we ask k T < 0.7P T with P T = | p h1 T − p h2 T |/2. This way, one can eliminate the contribution of two hadrons fragmented from the same parton. The event fractions affected by the gluon and quark initiated processes are shown in Fig. 10. Around 80% of the high p T di-hadron events are generated from gluon initiated processes in the small x B region. The fraction of quark initiated processes grows rapidly as x B approaches 0.1 and with increasing Q 2 . This behavior with Q 2 can be understood that more high p T hadrons are generated through QCD radiation, which has an increased probability with increasing Q 2 . It is especially noted that an understanding of the gluon Sivers function requires to measure its dependence on x B and Q 2 . Figure 11 compares the SSA for charged hadron pairs assuming the magnitude of the gluon Sivers function 5% of its positivity bound (solid circles) and the SIDIS1 set (solid triangles) as well as no gluon Sivers contribution, but a quark Sivers contribution (open circles). The study shows that a gluon Sivers function of a magnitude of 5% of the positivity bound can be measured at an EIC. We also find that the initial parton level asymmetry is significantly diluted, a factor 3 as indicated comparing the black curve and the solid circles. This dilution is larger than the one (factor 2) shown in Fig. 9(a) for DD pairs. The increase can be explained by the stronger smearing due to the fragmentation of light quarks to hadrons of the correlation between the parton and the hadron p T . Figure 11 shows the angular modulation of the Sivers function is only weakly dependent on Q 2 , because of the missing TMD evolution in the current framework, but much more sensitive on x B . The dependence on x B is a natural consequence of the behavior of the Sivers function parameterization with x.

C. Gluon SSA in Di-jet Production
Comparing the hadron level observables with jets, it can be clearly seen that jets provide a more precise reconstruction of the initial gluon kinematics. In the following we study the sensitivity to the gluon Sivers function in di-jet production. The jets are reconstructed from charged hadrons (π, K and protons) measured in the central tracker together with photons accepted in the calorimeter requiring a minimum transverse momentum p T Lab > 0.25 GeV and |η Lab | < 4.5. The jet radius parameter is assumed to be R = 1 in the anti-k T jet reconstruction algorithm. Di-jet events are defined with the trigger jet p jet1 T > 4.5 GeV and the associate jet p jet2 T > 4 GeV. Similar to the di-hadron channel, we use the vector sum of the transverse momentum for the two jets k T = | p jet1 T + p jet2 T | as the proxy to access the underlying gluon dynamics. We present the fractions for quark and gluon initiated processes in Fig. 12. The quark fraction is substantial (close to 30%) for low Q 2 events even for x B ∼ 10 −4 . The fraction of gluon initiated channels is maximized at small x B and drops below the quark fraction if x B is close to 0.01 or 0.1 depending on the Q 2 range. Unlike for the di-hadron case, the gluon event fraction increases with Q 2 . In Fig. 13, it is observed that a gluon Sivers function with the size of 5% of the positivity bound or the SIDIS1 set can be well separated from a SSA based on the quark Sivers effect at large x B for an integrated luminosity of L int = 10 fb −1 . Despite that the initial parton asymmetry for the dijet process is smaller than that for the di-hadron channel, a larger fraction of the initial asymmetry survives in the di-jet channel. The shape of the initial parton level asymmetry is largely preserved in the di-jet asymmetry in all kinematic variables. This is of great advantage to explore dependence of the gluon Sivers function on the hard scattering kinematics. Due to the strong correlation between the momentum of a jet and its mother parton, it is possible to reconstruct the momentum fraction of the parton participating in the hard interaction from the di-jet momentum information: x rec parton = (p jet1 T e −η jet1 + p jet2 T e −η jet2 )/W . The SSA as a function of x rec parton is shown in Fig. 13(c). The initial functional form of the gluon Sivers function on x is well reproduced in the measured SSA as function of x rec parton . The SSA based on a gluon Sivers function with a magnitude of 5% of the positivity bound drops while the one based on the SIDIS1 set increases with x rec parton . The shape of the initial gluon Sivers function is largely the same in the measured SSA. This will allow to distinguish different gluon Sivers models.
With projected high statistics for the di-jet channel, the evolution of the GSF with Q 2 and x B can be studied utilizing a multi-dimensional binning. Figure 14 shows the SSA based on the gluon Sivers function from the SIDIS1 set as function of Q 2 in three x B bins. The differential features of the SIDIS1 gluon Sivers function are: the asymmetry increases in the high x B bins, which is consistent with the behavior in x B . The decrease of the SSA as function of Q 2 is seen especially in the high x B bins. This signature can be utilized to study the evolution of the gluon Sivers function in the di-jet channel.
A more detailed analysis on the main cause of the smearing from the parton level asymmetry to the measured asymmetry in di-jet production is shown in Fig. 15. We present a comparison of the observed Sivers asymmetry in di-jet measurement with different hadronization assumptions to the probed parton level asymmetry. To study the influence of p T in the hadronization as well as the effect due to the decay of particle resonances, we have turned both processes consecutively off in the simulation by setting the respective PYTHIA parameters (PARJ (21) and MSTJ(21)) to zero. The solid red and black curves in Fig. 15 represent the parton level and measured asymmetry, respectively. Comparing the dotted curve (fragmentation p T off) and the dashed curve (particle resonance decay and fragmentation p T off) shows clearly the dominant effect is due to resonance decay in the fragmentation. As is shown in Fig. 15 of the asymmetry varying with x B , we observe the resonance decay becomes slightly more important in the high x B region. The remaining dilution of the parton level asymmetry is caused by QCD radiation, the p T dependence of the hard scattering process and the ability to measure a small transverse momentum imbalance with high p T di-jets. We also perform a study on the impact of using different algorithms and cone sizes in the jet reconstruction procedure. It is found that the effect by changing the jet reconstruction algorithm from anti-k T to k T or the cone size from R = 1 to R = 0.7 is barely visible. We present the result of di-jet asymmetry from cone size R = 0.7 with the dash-dotted line in Fig. 15. The choice of a smaller cone size only leads to a slightly smaller asymmetry. It is implied in this comparison that the di-jet asymmetry size is rather robust in spite of the jet reconstruction methods. An important aspect in comparing the different channels is the coverage in the gluon momentum fraction x g . Figure 16 shows the x g -distributions for the different channels are complementary. Both the heavy flavor DD mesons and the di-hadrons probe lower values of x g ( x g ∼ 0.03) and the di-jet channel is probing the larger x g range ( x g ∼ 0.1).

V. SUMMARY
We have performed a study on the feasibility of measuring the gluon Sivers function in high p T charged dihadrons, heavy flavor mesons and di-jet production in polarized ep ↑ collisions at an EIC. Scanning different assumptions of the magnitude of the gluon Sivers function provides a systematic study on the sensitivity of the different channels.
It is found that although the heavy flavor DD meson production is the cleanest channel to tag gluon initiated processes, it is at the same time also the most statistically challenging process and therefore the sensitivity to small gluon Sivers effects is limited. An alternative method using inclusive D mesons provides sensitivity to a gluon Sivers function with a magnitude of 10% of the positivity bound for a nominal integrated luminosity of L int = 10 fb −1 . But the smearing between parton level and the measured asymmetry is significantly increased. The high p T charged di-hadron channel is statistically more favorable and can resolve a magnitude of the gluon Sivers function down to 5% of the positivity bound. The most precise analyzer for the gluon Sivers effects at an EIC is the di-jet channel, due to its statistical advantage it provides the best sensitivity even for the small Sivers effects and can span the largest Q 2 -range to study TMD evolution effects. Due to its tight correlation between parton and jet kinematics, it has the smallest dilution between the parton level and measured asymmetries.
Overall it is thus the most promising experimental channel to determine and study all features of the gluon Sivers effect at the future EIC.