QCD evolution of the gluon Sivers function in heavy flavor dijet production at the Electron-Ion Collider

Using Soft-Collinear Effective Theory, we develop the transverse-momentum-dependent factorization formalism for heavy flavor dijet production in polarized-proton-electron collisions. We consider heavy flavor mass corrections in the collinear-soft and jet functions, as well as the associated evolution equations. Using this formalism, we generate a prediction for the gluon Sivers asymmetry for charm and bottom dijet production at the future Electron-Ion Collider. Furthermore, we compare theoretical predictions with and without the inclusion of finite quark masses. We find that the heavy flavor mass effects can give sizable corrections to the predicted asymmetry.


Introduction
In recent years, one of the most important forefronts of hadron physics has been the exploration of the three-dimensional (3D) partonic structure of nucleons in momentum space. Such 3D information is encoded in the so-called transverse-momentum-dependent parton distribution functions (TMD PDFs), which can further inform us about the confined motion of partons in the nucleon, as well as the correlation between their spins, momenta, and the spin of the nucleon [1]. Thanks to semi-inclusive deep inelastic scattering (SIDIS), a great deal of progress has been made in probing and extracting the TMD PDFs of quarkshowever, information regarding those of gluons is still largely unknown experimentally. Exploring and measuring gluon TMD PDFs is one of the primary goals for the future Electron Ion Collider (EIC).
Among the gluon TMD PDFs, the so-called gluon Sivers function is regarded as one of the "golden measurements" at the future EIC [1]. The gluon Sivers function encapsulates the quantum correlation between the gluon's transverse momentum inside the proton and the spin of the proton, thus providing 3D imaging of the gluon's motion. Quite a few processes have been proposed to probe the gluon Sivers function at the EIC, including heavy quark pair production [2], heavy quarkonium production [3][4][5][6][7], and quarkonium-jet production [8], as well as back-to-back dihadron and dijet production [9]. The feasibility of measuring the gluon Sivers function in the above scenarios has been studied in [10], where the authors use the PYTHIA event generator [11] and the reweighting method of [12] to investigate the spin asymmetry. They conclude that dijet production is the most promising channel for probing gluon Sivers effects, where the selection of a sufficiently small-x value suppresses the contribution of the quark channel and the corresponding quark Sivers function. In this paper, we discuss spin asymmetry in the process of heavy flavor (HF) dijet production, where the contribution of the quark Sivers function is further suppressed compared to that of the light flavor dijet case.
An intriguing feature common to both quark and gluon Sivers functions is that they depend non-trivially on the processes in which they are probed. A well-known example of the process-dependence of the quark Sivers function is its sign change between SIDIS and Drell-Yan processes [13][14][15]. Similarly, it has been demonstrated that the gluon Sivers function for the process of back-to-back diphoton production in p + p collisions, p ↑ p → γγX, carries a sign opposite to that of dijet production in e + p collisions, ep ↑ → e jjX: [2]. In [16], it was demonstrated that the gluon Sivers function in any process can be expressed in terms of two "universal" functions with calculable color coefficients for each partonic subprocess. We briefly discuss such a process-dependence for HF dijet production below. For a comprehensive review on gluon TMD PDFs, see [17,18].
So far, studies of the gluon Sivers function at the EIC are mostly performed within the leading-order (LO) parton model, without considering the impact of QCD evolution. The effects of resummation for back-to-back light flavor dijet production in the unpolarized DIS process have been investigated in [19], where the authors apply the p T -weighted recombination scheme [20] in defining the jet axis to avoid the theoretical complexity arising from non-global logarithm (NGL) resummation [21]. A similar idea is used to study single inclusive jet production in the Breit frame at the EIC in [22,23]. Recently, following the same Soft-Collinear Effective Theory (SCET) framework utilized in [24][25][26][27][28], the TMD factorization formula for light flavor dijet production at the EIC has been derived [29], where the azimuthal-angle-dependent soft function, describing the interaction between two final-state jets through the exchange of low-energy gluons, is analytically calculated at one-loop order. For HF jet production in the kinematic region of comparable jet and heavy quark masses, a new effective theory framework is needed. In this work, we provide such a framework and derive the TMD factorization formula.
The remainder of this paper is organized as follows. In Sec. 2, we detail the factorization framework required to carry out resummation in the back-to-back region where the transverse momentum imbalance of the HF dijet is small. In Sec. 3, we present numerical results for charm and bottom dijet production in both unpolarized and transverselypolarized-proton-electron scattering. We summarize our findings and give an outlook for future investigations in Sec. 4.

Factorization and resummation formula
In this section, we start with the kinematics for HF dijet production in e + p collisions. We then provide the TMD factorization formalism with explicit expressions for all the relevant factorized ingredients. Figure 1. HF dijet production in electron-proton collisions, as stated in Eq. (2.1).

Kinematics
As shown in Fig. 1, we consider HF dijet production in the polarized-proton-electron scattering process where S T is the transverse spin of the polarized proton with momentum P and ( ) is the momentum of the incoming (outgoing) electron. At LO, HF dijets are produced via the γ * g → QQ process. The HF quark Q and antiquarkQ initiate the observed HF jets J Q and JQ with momentum p J and p J , respectively. In this paper, we choose to work in the Breit frame so that both the virtual photon (with momentum q = − ) and the beam proton scatter along the z-axis. For convenience, we define the following variables commonly used in DIS, We may further note that Q 2 = x B y S P , where S P = ( + P ) 2 denotes the electronproton center-of-mass energy. In a fashion analogous to SIDIS, we also define the kinematic variable z = P · p J /P · q, which gives the momentum fraction of the photon carried by the jet J Q . At LO, the four-momenta of the incoming and outgoing particles are expressed as where we have introduced two light-like vectors, n µ = (1, 0, 0, 1) andn µ = (1, 0, 0, −1), and define p µ t such that p µ t p tµ = −p 2 T with p T = p T (cos φ J , sin φ J ). We denote transverse momenta relative to the photon-proton beam by the subscript T , while that relative to the jet direction is given the subscript ⊥. Here, we assume p 2 T m 2 Q and take p 2 J = p 2 J = 0. This allows us to derive the factorized cross section in the following section. Lastly, the parton-level Mandelstam variables can be defined aŝ where x is the momentum fraction of the proton carried by the gluon, and is given by

Factorization formula
In the Breit frame, we define the dijet imbalance as q T = p JT + p JT . For this paper, we examine the back-to-back configuration where q T p JT ∼ p JT ≡ p T . Furthermore, we work in the kinematic regime where m Q p T R p T , with R denoting the jet radius. Overall, in the region with the scale hierarchy as q T R q T m Q p T R p T , the factorized expression for the proton-spin-independent cross section is given by Above, y J is the rapidity of the HF jet J Q and is related to the kinematic variable z through the relation z = e y J p T /Q. In the factorization formula Eq. (2.8), S denotes the soft function while f g/N is the unpolarized gluon TMD PDF. Their perturbative oneloop expressions can be found in Sec. 2.4. In the third line of Eq. (2.8), J Q and S c Q are the massive quark jet and collinear-soft functions, which differ from the corresponding functions utilized in light jet production [24][25][26][27][28]. In Secs. 2.5 and 2.6, we present their explicit calculations at next-to-leading order (NLO). The variables k T , λ T , and l T label the transverse momenta associated with the collinear, soft, and collinear-soft modes. Finally, µ and ν are the factorization and rapidity scales, respectively, while ζ is the Collins-Soper parameter [30,31]. In the derivation of the above factorization formula we apply the narrow jet approximation with R 1. However, as shown in [32][33][34][35] this approximation works well even for fat jets with radius R ∼ O(1), and the power corrections of O(R 2n ) with n > 0 can be obtained from the perturbative matching calculation.
Fourier transforming to b-space, the factorized cross section becomes where soft function S and the gluon TMD PDF f g/N both depend on the rapidity scale ν. However, the soft function can be written as where S nn (b, µ, ν) is the soft function for Higgs production in p + p collisions [36,37], and the function S(b, µ) on the right-hand side no longer depends on the rapidity scale ν. Upon making this replacement, the factorized expression for the cross section can be written in terms of the properly-defined TMD gluon distribution [30] by noting that Here, f TMD g/N (x, b, µ, ζ) on the right-hand side is defined as This is the properly-defined gluon TMD PDF probed in Higgs production in p + p collisions [37] and is thus the counterpart of the quark TMD PDF as probed in Drell-Yan lepton pair production. Finally, Eq. (2.9) can be expressed in the following form In the following sections we calculate the one-loop expressions for all the above functions. An important physical requirement is that the factorized cross section must be independent of the scale µ-we verify this factorization-scale-independence in Sec. 2.7. Next, if one considers the scattering of an electron with a transversely-polarized proton with spin S T , Eq. (2.8) can be generalized. In this case, the spin-dependent cross section is given by the sum where dσ U T depends on the gluon Sivers function. The full expressions for the leading twist gluon distributions are given in [38]. Using these results, we can obtain the expression for the polarized cross section by simply replacing the unpolarized gluon TMD PDF in Eq. (2.8) with the gluon Sivers function, namely Here, it is important to note that there exist both f -and d-type gluon Sivers functions, which are associated with different color configurations in the three-gluon correlator, i.e., involving the antisymmetric f abc and symmetric d abc structure constants of SU (3), respectively. For details, see for instance [16,39]. In Eq. (2.15), we have denoted the gluon Sivers function with the superscript f , which is used to indicate that it is f -type. We note that at LO, this process is only sensitive to the f -type function. Further details on this matter are provided in Sec. 2.3. After making this substitution, the factorized cross section then reads where H Sivers denotes the hard function for the polarized process, and this expression can once again be written as a Fourier transform by defining Finally, the factorization formula for the polarized differential cross section becomes Here, we have applied the redefinition Eq. (2.10) to obtain the rapidity-scale-independent gluon Sivers function

Hard function
In the unpolarized process, the LO hard function is determined by the tree-level cross section for dijet production in DIS, which is expressed as [40,41] where α em is the fine structure constant and Q f denotes the fractional charge of the HF quark. On the right-hand side, the first superscript U indicates that the incoming gluon is unpolarized, while the second superscripts {U + L, L, I, T } correspond to the different helicity states of the off-shell photon. Explicitly, the functions H U,i are expressed as One immediately sees that the functions H U,I and H U,T vanish upon integrating out the azimuthal angle φ J of the jet. As such, these contributions do not play a role in our numerical calculations. The expression for the hard anomalous dimension can be obtained from the calculation of the 3-jet process γ * → qqg at lepton collisions [42,43]. The hard anomalous dimension can also be read from the general structures in [44,45] and is given as where γ cusp is the cusp anomalous dimension, while γ q and γ g represent the single logarithmic anomalous dimensions for the quark and gluon, respectively. With this anomalous dimension, one can then perform resummation by solving the following renormalization group (RG) equation for the hard function where, for brevity, we maintain only the scale µ-dependence in the hard function. We note that in order to perform the evolution at next-to-leading logarithmic (NLL) accuracy, the cusp anomalous dimension is needed at two-loop order and the single logarithmic anomalous dimensions are needed at one-loop order. The values of these expressions are where we have organized the perturbative expansion of each anomalous dimension as For the polarized process, we must consider the process-dependence of the corresponding gluon Sivers functions [16]. Such process-dependence can be computed via the attachment of an additional gluon originating from the gauge link in the definition of the Sivers function. This additional gluon is responsible for the soft pole that generates Sivers asymmetry. This method is widely used in computing the process-dependence of the quark Sivers function, see e.g. [46], which gives the same results as shown in [16,47]. In Fig. 2, the soft poles are represented by red lines. We note that for both the polarized and unpolarized cases, the hard functions can be expressed as matrices in color space [28]. For more complicated processes, the relationship between the polarized and unpolarized hard matrices is non-trivial. However, for the γ * q → QQ process, the color space is one-dimensional and, therefore, the polarized hard function can be simply written as (2.26) Here, C 1 and C 2 are the color factors for the polarized hard process associated with the attachment of the additional gluon to the HF quark and anti-quark [3,48,49], respectively. The function h(Q, y, p T , y J , µ) is the kinematic part of the hard function. For the unpolarized case, the hard function can be written as where factor C u is the color factor associated with the unpolarized hard process. We find that at LO for this process, the attachment of this additional gluon originating from the gauge link in the definition of the Sivers function produces color configurations which are proportional to (−if abc ). This analysis indicates that while there are both d-and f -type gluon Sivers functions [16], this process at LO is only sensitive to the f -type gluon Sivers function. In addition to this, while the term (−if abc ) in Fig. 2 appears in the hard function, this term should be absorbed into the definition of the gluon Sivers function as it originates from the Wilson line in the adjoint representation. Similarly, the term δ ac in that figure should be absorbed into the definition of the unpolarized TMD PDF. For this process, we find that (C 1 + C 2 ) = C u . As a result, the polarized and unpolarized hard functions are equal and the hard anomalous dimension is unchanged for the polarized case.

TMD PDFs and the global soft function
In order to regulate the rapidity divergences in the TMD PDFs and the global soft function, we use the η regulator of [36]. The expression for the unsubtracted gluon TMD in this regularization scheme can be obtained from [50], and is given by where we have set ζ = (n · p g ) 2 [31], the scale µ b is defined as Here the anomalous dimensions are given by with γ fg 0 = γ g 0 , which is given in Eq. (2.24). On the other hand, we define the soft function as where the soft integrals I ij are given by (2.37) Here, i and j are either B, J, or J, where B denotes the beam direction while J and J denote the directions of the jet initiated by the HF quark and anti-quark, respectively. The expressions for the φ b -dependent beam-jet soft function integrals are given in [24] as with c bJ = cos(φ b − φ J ) and where the rapidity of jets are Here, the terms marked by "fin" in their superscripts denote the finite contributions of their respective functions. While the divergent pieces are required for the purposes of resummation, the finite pieces are only needed at NLO. Since we perform our analysis at NLL accuracy, these terms are not needed for this study. Recently in [29], the authors derive the following the jet-jet soft function integral which contains no rapidity divergence. This integral can be written as with the one-loop single logarithmic anomalous dimension Thus, the product of f g/N (x, b, µ, ζ/ν 2 ) and S(b, µ, ν) is ν-independent, and we can construct the properly-defined gluon TMD PDF as in Eq. (2.11).
Here we find that the soft function depends not only the magnitude but also the direction of the vector b. As shown in Sec. 2.6 a similar structure also shows up in the collinear-soft function, and the φ b dependence in the anomalous dimensions will cancel out between these two functions. However, after taking into account the evolution between the soft and collinear-soft function, one finds that the φ b integral is divergent in some phase space region. In order to avoid such divergences we apply methods in [24,28] where one first performs an averaging over the φ b angle in both soft and collinear-soft function. We note that this method does not change the RG invariance as shown in Eqs. (2.73) and (2.74). In addition, as discussed in [27] no significant numerical effects between different methods are observed in the NLL resummation calculation.
The φ b -averaged soft function can be constructed from Eqs. (2.35) and (2.36) by replacing the soft integrals with where we have placed a bar over these integrals in order to distinguish them from the φ b -dependent ones and have defined ∆y = y J − y J . These results are the same as the soft function calculated in [51]. Therefore, the anomalous dimensions for the averaged case arē where the one-loop single logarithmic anomalous dimensions is By comparing Eqs. (2.43) and (2.50), one can see the rapidity anomalous dimension is unchanged. Therefore in the φ b -averaged case, we can once again write the factorized expression in terms of the properly defined TMD PDFs.

Massive quark jet function
In this section, we discuss the calculation of the massive quark jet function at NLO. The massive quark jet function has been investigated in detail for various observables. For example, the factorization formula for the massive event shape distribution involves such a jet function, as the jet and heavy quark masses are of similar magnitude [52][53][54][55]. The corresponding jet function has been calculated to two-loop order [56]. Furthermore, the semi-inclusive massive quark jet fragmentation function has been calculated at NLO and applied to inclusive jet production [57,58]. Recently, the one-loop expression for the socalled unmeasured massive quark jet function has been presented in [59]. The global jet anomalous dimension can be obtained from the divergent terms of the unmeasured massive quark jet function. As shown in Fig. 3, the one-loop calculation involves two types of diagrams: J NLO,V Q and J NLO,R Q , where J NLO,V Q contains only single cut propagators and is thereby unconstrained by the jet algorithm. Explicitly, it is written as where the heavy quark mass m Q is the only physical scale involved. Since the real contribution J NLO,R Q is constrained by the jet algorithm, it will depend on the jet scale p T R in addition to m Q . In this work, we define the HF quark four-momentum q µ with q 2 = m 2 Q , which is known as the M-scheme [55]. We note that in the hierarchy of scales we are considering, the constraint of the anti-k T algorithm [60] is independent of the HF quark mass m Q and is in fact identical to that for massless partons [57], namely where q µ = (q + , q − , q ⊥ ) is the four-momentum of the HF quark and ω J is the large component of the jet four-momentum. The jet scale p T R emerges in Eq. (2.53) upon noting ω J = 2 p T cosh y J . In the phase space integral, we expand the integrated momentum q along the jet direction with q + = (m 2 Q + q 2 ⊥ )/q − given by the power counting requirement p T R ∼ m Q . Explicitly, we have where only a single divergence is exhibited, as the heavy quark mass m Q acts as a regulator of the overlapping soft and collinear regions of phase space. After combining the real and virtual contributions, the logarithmic dependence on the quark mass m Q cancels out. The one-loop global jet renormalization constant then reads where, again, we observe that the heavy quark mass m Q only affects the single pole structure. We further note that as m Q → 0, the massive quark jet renormalization constant reduces to that of the massless jet, Z J Q → Z Jq . This gives us the following expression for the global jet anomalous dimension with the one-loop single logarithmic anomalous dimension as where the first term in the brackets is shared by the massless quark jet function and the second term constitutes the finite quark mass correction. Finally, the renormalized HF jet function is given by the following where the function F(p T R, m Q ) can be expressed as This expression for the HF jet function is equivalent to the semi-analytic form presented in [59], and one can see that as m Q → 0, we have F → 0 and, therefore, J Q → J q . Hence, the massive quark jet function behaves as expected in the massless limit. Explicitly, the bare NLO collinear-soft function is given by where the collinear-soft integrals w αβ are defined in b-space as (2.62) Upon performing the k-integration, we obtain the following expressions for w αβ We see that the finite quark mass corrections only enter into the single pole structure of the collinear-soft function. This is analogous to the observation made in Sec. 2.5 in analyzing the massive quark jet function, and can be understood through the same physical reasoning. The finite terms are given by (2.66) Here, we note that wn J v J reduces to the massless wn J n J function [24,27] as m Q → 0, while w v J v J vanishes. Given the expression for S c Q , we can calculate the renormalization constant Z S c Q , which is given by This renormalization constant leads to the following formula for global collinear-soft anomalous dimension where the one-loop single logarithmic anomalous dimension is (2.69) The anomalous dimension for the collinear-soft function associated with the anti-quark is given by For phenomenological purposes, we utilize the φ b -averaged collinear-soft function, which can be obtained through Eq. (2.63) by making use of the following integrals The resulting anomalous dimension for the φ b -averaged collinear-soft function is denoted byΓ cs Q withγ (2.72) Upon integrating over φ b , we find thatΓ csQ (α s ) =Γ cs Q (α s ) and, therefore, the two averaged collinear-soft functionsS c Q andS cQ behave identically under QCD evolution.

Renormalization group consistency
Armed with the anomalous dimensions of each component, we are now positioned to demonstrate the RG consistency of our factorization framework. Inspection of Eqs. (2.56) and (2.68) reveals that all mass corrections cancel exactly in the sum Γ j Q + Γ cs Q , making the RG consistency of our formalism identical to the massless case. A similar observation is made in [59]. This general physical behavior has also been observed in the context of inclusive HF jet production [57], where the authors offer the intuitive argument that as the heavy quark mass m Q constitutes IR information, it thus does not affect the UV behavior of the semi-inclusive jet function. In the present context, we see that the UV evolution behavior of the product of the jet and collinear-soft functions is insensitive to the IR scale introduced by the heavy quark mass. However, in Sec. 2.8, we will see how the heavy quark mass enters non-trivially and crucially into the evaluation of the differential cross section.
Therefore, upon combining Eqs. (2.73) Furthermore, we note that this consistency is preserved under the operation of φ b -averaging

Resummation formula
Utilizing our EFT framework, all-order resummation is achieved through RG evolution. The resulting all-order expression for the HF dijet production cross section is given at NLL 1 by where J 0 is the zeroth order Bessel function of the first kind. In this expression, µ h , µ j , and µ cs are the hard, jet, and collinear-soft scales, respectively. We have also performed the usual operator product expansion (OPE) of the unpolarized gluon TMD PDF f g/N (x, b, µ, ζ) in terms of the collinear gluon PDF f g/N (x, µ) at the initial scales ζ i = µ 2 i = µ 2 b * , and have kept the coefficient at LO to be consistent with NLL accuracy. The matching coefficient at higher-orders can be found in e.g. [37,[62][63][64][65][66][67]. The function S NP parameterizes the contribution from non-perturbative power corrections which are enhanced for q T ∼ Λ QCD . Explicitly, we apply the formula given in [68], which reads We assign the following values to the parameters: g 1 = 0.106 GeV 2 , g 2 = 0.84 and Q 2 0 = 2.4 GeV 2 , and use the standard b * -prescription [69] b * = b in order to regularize the Landau singularity as b → ∞. Moreover, the spin-dependent cross section is expressed as Here, we have expected a similar OPE for the gluon Sivers function f ⊥ 1T,g/N (x, b, µ, ζ) at the initial scales ζ i = µ 2 i = µ 2 b * and simply expressed the corresponding collinear function at LO as f ⊥,f 1T,g/N (x, µ b * ) for simplicity. In principle, the corresponding collinear functions in the OPE expansion would be the twist-3 three-gluon correlation functions defined in [70,71]. To the best of our knowledge, detailed OPE calculations for the corresponding coefficient functions are not available in the literature. An expansion of the gluon Sivers function in terms of the collinear twist-3 quark-gluon-quark correlator, or the so-called Qiu-Sterman function [72,73], in transverse momentum space is performed in [3]. On the other hand, the coefficient functions for the expansion of the quark Sivers function in terms of the threegluon correlation functions are provided in [74,75]. The computation of the coefficient functions for expanding the gluon Sivers function in terms of the three-gluon correlation functions is essential for a full understanding of the QCD evolution of the gluon Sivers function. We leave this to future work.
Our knowledge about gluon Sivers functions, especially in the proper TMD factorization formalism, is rather limited. At the present moment, the only experimental constraint on the gluon Sivers function, in the TMD framework, comes from the SIDIS measurement of back-to-back hadron pairs off transversely-polarized deuterons and protons at COM-PASS [76]. However, as of yet, there has been no theoretical extraction of the gluon Sivers function from such data. On the other hand, an important theoretical constraint on the gluon Sivers function comes from the Burkardt sum rule [77]. For the phenomenological purposes of the next section, we adopt the non-perturbative parameterization utilized by [78,79] 2 . Specifically, for the non-perturbative Sudakov, we take where the g 2 -dependent term is spin-independent and is, therefore, the same term occurring in Eq. (2.76), while the term ∝ g 1 ρ can be connected to the Gaussian width in transverse momentum space [81] for the gluon Sivers function. For the collinear part of the gluon Sivers function, f ⊥,f 1T,g/p (x, µ) in Eq. (2.78), we take At this point, it is important to note that while the mass corrections in sum of the anomalous dimensions for the collinear-soft and massive jet functions cancel, the massdependence of Γ j Q contributes to the differential cross section. By examining Eqs. (2.75) and (2.78), we see that the mass corrections enter into the evolution between the scales µ j and µ cs . We will see in the following section that this can significantly affect both the q T -distributions and spin asymmetries for HF dijet production at the EIC.

Numerical results
In this section, we present numerical results for HF dijet production in unpolarized and transversely-polarized-proton-electron collisions at the future EIC. We set the energies of the electron and proton beam to be 20 GeV and 250 GeV, respectively. These beamenergy values yield a electron-proton center-of-mass energy of √ S P = 141 GeV. For the all-order resummation formulae in Eqs. (2.75) and (2.78), the renormalization scales for each function are chosen to be Here, note that the Landau singularity associated with the collinear-soft scale is also regularized by the b * -prescription. As given in the calculation of the jet function, we consider HF jets constructed using the anti-k T algorithm with radius R = 0.6. The corresponding kinematic cuts for charm and bottom jets in the Breit frame are charm jets : 5 GeV < p T < 10 GeV, |y J | < 4.5 , 2 Note that the gluon Sivers function in [78], and its updated version [80], is constrained to their study of the p ↑ p → πX process. Technically, this is not subject to a TMD factorization framework, but it serves as a starting point for our numerical study, following [10].
In Fig. 5, we display the normalized unpolarized cross section, 1/σ dσ/dq T , as a function of the imbalance q T . In Fig. 6, the Sivers spin asymmetry A sin(φq−φs) U T is presented as a function of q T /p T following [83], for both charm (left panel) and bottom (right panel) jets, respectively. For both plots, the solid curves are the results obtained using the resummation formula, while the dashed curves represent the resummation prediction using the evolution kernel without finite quark mass corrections. For both the unpolarized q T and A sin(φq−φs) U T distributions, we find that the effects of the finite quark masses are modest for charm jets and quite sizable for bottom jets. This can be attributed to the sizes of the charm and bottom masses relative to their associated jet scales p T R. As discussed in Secs. 2.5 and 2.6, we have that J Q → J q and S c Q → S c q as m Q → 0, making them analytic functions of m Q in the neighborhood of zero mass. Since Eqs. (2.75) and (2.78) carry their mass-dependence through the anomalous dimensions for the jet and collinear-soft functions, Eqs. (2.56) and (2.68), one sees that the massive versions of these functions are connected to the massless versions by the ratio m Q / (p T R)-it is in fact this dimensionless parameter that controls the physical size of the mass corrections. With this in mind, one sees that Eq. m Q / (p T R)) from light flavor jet-pairs than it does charm dijets. This relative positioning is then clearly displayed in Figs. 5 and 6. In order to estimate the theoretical uncertainties, in both Figs. 5 and 6 we also show the uncertainties from scale variations, which are given by the red and blue bands. Here we vary the hard and jet scales by a factor of two around their default values as defined in (3.1), and the total uncertainty bands are obtained by the envelope of all the variations. Since the non-perturbative Sudakov factor in Eq. (2.79) is fitted at the canonical scale µ b * , we do not include theory uncertainties from µ b * and µ cs variations. We find that the scale uncertainty is compatible with the finite quark mass corrections in charm dijet process, while its impact on the bottom dijet process is smaller than the mass correction. Therefore in order to identify the finite quark mass effects in the charm dijet process it is essential to reduce the scale uncertainties. Our factorization and resummation formula provides a clear structure to improve the perturbative accuracy, which makes scale uncertainty further reduction possible. We leave the higher-order perturbative calculations in future work.

Conclusion
A major priority of the future EIC is to explore the gluon TMD PDFs. In this paper, we have investigated the use of back-to-back HF dijet production in transversely-polarized target DIS as a means of probing spin-dependent gluon TMD PDFs. We have calculated the expressions for the mass-dependent jet and collinear-soft functions at next-to-leading order. Using these expressions, as well as Soft-Collinear Effective Theory, we resum the large logarithms associated with these expressions at next-to-leading logarithmic accuracy. We then provide a factorization theorem for this process with QCD evolution in the kinematic region where heavy quark mass m Q p T R p T , with p T and R being the transverse momentum and the radius of the jet, respectively. Furthermore, we generate a prediction for the Sivers asymmetry for charm and bottom dijets at the EIC, which can be used to probe the gluon Sivers function. We carefully study the effects of the HF masses by comparing our mass-dependent predicted asymmetry against the asymmetry in the massless limit. We find that, in the kinematic region we consider, the HF masses generate modest corrections to the predicted asymmetry for charm dijet production but sizable corrections for the bottom dijet process. Furthermore, we also consider the theoretical uncertainties from the scale variation. We find that the scale uncertainty can be compatible with the corrections from finite quark mass effects, especially for charm dijets production. In order to identify the mass effects and reduce the scale uncertainties one has to include higherorder corrections in the matching coefficients and the corresponding anomalous dimensions in Eqs. (2.75) and (2.78), and we leave the detailed perturbative calculations in future work.