The $g_T(x)$ contribution to single spin asymmetries in SIDIS

Motivated by a novel origin of transverse single spin asymmetry (SSA) in semi-inclusive Deep Inelastic Scattering (SIDIS) uncovered by some of us, we quantitatively investigate its impact on the theoretical understanding of the mechanism responsible for SSA. This new contribution from the quark-initiated channel first appears in two-loop perturbation theory and involves the $g_T(x)$ distribution. We point out another entirely analogous piece from the gluon-initiated channel proportional to the gluon helicity distribution $\Delta G(x)$. Both contributions are solely expressed in terms of twist-two polarized parton distribution functions and twist-two fragmentation functions in the Wandzura-Wilczek approximation, such that they can be unambiguously evaluated without introducing free parameters. We make predictions for measurements of the asymmetries $A_{UT}$ at the future Electron-Ion Collider (EIC), and find that $A_{UT}$ associated with the $\sin (\phi_h-\phi_S)$, $\sin \phi_S$ and $\sin (2\phi_h-\phi_S)$ harmonics can reach up to 1-2\% even at the top EIC energy.


I. INTRODUCTION
Recently, three of us, together with D. J. Yang, have proposed a novel mechanism for generating transverse singlespin asymmetry (SSA) in semi-inclusive deep inelastic scattering (SIDIS) ep ↑ → e hX [1]. It has been demonstrated that an imaginary phase necessary for SSA can be produced purely within a parton-level cross section starting at two loops. The spin-dependent part of the cross section at high transverse momentum P hT > 1 GeV (measured with respect to the virtual photon direction) can be schematically written as where g T (x) is the twist-three parton distribution function (PDF) associated with a transversely polarized proton, D 1 is the unpolarized twist-two fragmentation function (FF) for the observed hadron h, and H is the hard kernel starting at O(α 2 s ) (see also an earlier related work [2]). The terms omitted in (1) are proportional to the 'genuine twist-three' quark-gluon correlation functions ∼ ψ gF ψ commonly called the Efremov-Teryaev-Qiu-Sterman (ETQS) functions [3,4]. As is well known, the g T distribution can be written as the sum of the Wandzura-Wilczek (WW) part [5] and the genuine twist-three part where ∆q(x) is the twist-two polarized (helicity) quark PDF. It is a consistent truncation of the result in Ref. [1] to keep only the WW part in (2). The new source of SSA can then be entirely expressed in terms of the twist-two PDFs ∆q(x) and the twist-two FFs D 1 (z). This is a remarkable observation in striking contrast to the prevailing view that SSA at high-P hT is explained by the ETQS functions and certain twist-three fragmentation functions (see a review [6] and references therein). Unlike these higher-twist distributions, twist-two distributions are very well constrained by global QCD analyses. Therefore, the mechanism proposed in [1] offers a unique part of SSA that can be unambiguously calculated without introducing free parameters. Moreover, in the transverse-momentumdependent (TMD) PDF framework valid in the low-P hT region (P hT 1 GeV), a new source of SSA proportional to the g 1T (x, k ⊥ ) distribution (the TMD version of g T (x)) has been identified, along with more than a dozen of new contributions involving various twist-three TMDs and FFs and hard kernels up to two loops. Again, this calls into question the prevailing view in the community (see, e.g., [7]) that SSA at low-P hT is entirely attributed to the Sivers and Collins functions.
The purpose of this paper is twofold. First, we extend the analysis of [1] to gluon-initiated channels. There exists a gluonic counterpart of g T (x), the twist-three G 3T (x) distribution [8,9] for a transversely polarized proton. Its WW part is related to the twist-two polarized gluon PDF ∆G(x). In complete analogy to (1), we find the structure arXiv:2109.05440v1 [hep-ph] 12 Sep 2021 which again consists only of twist-two distributions after the WW approximation. We shall identify the two-loop diagrams that go into the hard kernel H g and study their gauge invariance and infrared safety. Equation (3) is a novel gluon-initiated source of SSA in SIDIS to be considered jointly with the previously known mechanism which involves genuine twist-three, three-gluon correlators F gF F [10][11][12]. Second, we perform a detailed numerical analysis of SSA and make predictions for its measurements at the future Electron-Ion Collider (EIC) [13,14]. In doing so, we neglect the 'usual' contributions from the ETQS and twist-three FFs, which have been intensively anallyzed, and focus on the new contributions in order to explore their importance clearly. The results can be viewed as a baseline for future EIC measurements of SSA. Deviations from our predictions, if observed, may be attributed to genuinely twist-three effects.
This paper is organized as follows. In Section II, we describe the SIDIS setup and introduce kinematic variables. In Section III, we first review the result of [1] obtained for the quark-initiated channel, and then propose an analogous, but novel contribution to SSA in the gluon-initiated channel. In Section IV, we perform a detailed analysis of the two-loop diagrams and calculate the hard coefficients in all the partonic channels. In Section V, we implement the obtained formulas numerically and make predictions according to the kinematic coverage of the EIC. We also present results relevant to the COMPASS experiments [15]. Finally, we discuss our findings and conclude in Section VI. Appendices are devoted to a technical proof of the infrared finiteness of the factorization formulas at two-loop level.

II. SIDIS KINEMATICS
In this section we give a brief review of polarized SIDIS e(l)p(P ) → e(l )h(P h )X and introduce involved kinematic variables. We have in mind light hadron production specifically for h = π ± . Heavy-quark production will be studied in a separate work. The spin-dependent part of the differential cross section is given by where is the leptonic tensor, W µν is the hadronic tensor, and ν and µ are the polarization indices of the virtual photon in the amplitude and the complex-conjugate amplitude, respectively. The Bjorken variable is denoted as . We shall work in the so-called hadron frame, where the virtual photon and the proton move in the z direction with The incoming and outgoing leptons have the momenta where φ is the azimuthal angle relative to the z axis, and We adopt the standard variables with the relation In the present frame, the transverse part of q µ t reads q tT = − P hT /z f . The momentum of the final state hadron can then be parametrized as For the transverse spin of the incoming proton we choose In terms of the above variables, the differential cross section (4) takes the following Lorentz invariant form In practice, instead of φ and χ, it is more convenient to define the hadron and spin angles relative to the lepton plane, in accordance with the Trento conventions [16]. The cross section is then a function of φ h and Φ S − χ = φ h − φ S . The dependence on φ h can be factored out by decomposing the hadron tensor W µν using the following set of vectors [17], which form nine independent tensors, V µν k (see [17] for explicit expressions), and their inverses,Ṽ µν k . Among them, the following six symmetric tensors [17] contribute to the decomposition of W µν , With these tensors we can write where have the explicit expressions We are thus led to the representation (see for example [18]) with The Fourier components sin(φ h − φ S ) and sin(φ h + φ S ) are referred to as the Sivers and Collins asymmetries, respectively. While we continue to use this nomenclatures, we emphasize that the new mechanism, which contributes to these asymmetries and will be studied in detail below, has nothing to do with the Sivers and Collins functions, or their collinear twist-three counterparts.

III. NEW CONTRIBUTIONS TO SSA
In this section we first recapitulate the g T (x) contribution to SSA discussed in [1], and apply the so-called Wandzura-Wilczek (WW) approximation to simplify the result. We then derive another new contribution to SSA due to the gluonic counterpart of g T .

A. Quark-initiated channel
In [1], it has been shown that the imaginary phase necessary for SSA in SIDIS can come from the hard kernel in perturbation theory starting at two loops, and all the relevant two-loop diagrams have been identified. However, only the quark (and antiquark) initiated channel was considered there. In this channel, motivated by the structure (1), we factorize the fragmentation function out of the hadronic tensor W µν as where we have taken into account the fact that the observed hadron can also come from the fragmentation of a radiated gluon in the final state. The result of [1] reads (suppressing the label a for simplicity) in which M N is the proton mass, and Tr denotes trace over colors and Dirac indices. Our conventions are 0123 = +1, γ 5 = iγ 0 γ 1 γ 2 γ 3 and αP nS = αβρλ P β n ρ S T λ with the light-like vector n µ satisfying n 2 = 1 and n · P = 1. The g T (x) distribution function is defined as dλ 2π while G F (x 1 , x 2 ) andG F (x 1 , x 2 ) are the ETQS functions (We follow the notation of Ref. [19] where explicit definitions can be found). The 'kinematical' distributionsg(x) and g T (x) are related through the QCD equation of motion The hard matrix elements S (0) µν (xP ) and S (1) µνα (x 1 P, x 2 P ) are computable in perturbation theory. As observed in [1], the first nonzero contribution to S (0) appears at two loops, S (0) ∝ α 2 s , and S (1) is obtained from S (0) by attaching an extra gluon in all possible ways. Two representative diagrams contributing to S (0) are shown in Fig. 1. The crosses denote on-shell lines that lead to an imaginary phase via the Cutkosky rules. We note that Ref. [2] arrived at essentially the same structure as (22), but did not specify the hard kernels S (0,1) .
We shall compute w µν in the WW approximation, namely, by systematically neglecting genuine twist-three distributions everywhere. This is a consistent approximation in the sense that it preserves both QED and QCD gauge invariance. In this approximation, we may write 1 where ∆q(x) and ∆q(x) are the standard twist-two polarized quark and antiquark distributions. Moreover, the first two lines in (22) can be combined into We thus arrive at the structure mentioned in the introduction, The above formulas hold for each quark flavor. In practice, we must sum over flavors weighted by the quark electromagnetic charge. In physical cross sections, we thus apply where g T f is given by (25) for each quark flavor f .
Let us compare (27) with the conventional contribution from the ETQS function [19][20][21] which schematically reads 2 Since H ∼ O(α 2 s ) and H ∼ O(α s ), naively the former is parametrically suppressed by a factor of α s . However, the definition G F ∼ ψ gF ψ explicitly contains the coupling g which actually comes from perturbative diagrams. That is, some suppression associated with the coupling g goes into G F in the convention (28). As for the soft part, both conceptually and practically, we have a far better grasp of twist-two distributions than twist-three distributions: ∆q and D 1 have been well constrained thanks to a wealth of experimental data and global QCD analyses, whereas the ETQS functions are still poorly constrained. These considerations make (27) a new and attractive source of SSA that can be unambiguously calculated without introducing any free parameters. The main goal of this paper is to carry out such a calculation, both analytically and numerically. But before doing so, let us point out that an entirely analogous contribution exists in the gluon-initiated channel.

B. Gluon-initiated channel
The gluonic counterpart of g T (x) for a transversely polarized proton is defined as [8,9,22] dλ 2π Similar to g T (x), it can be written as the sum of the WW part and the genuine twist-three part, where the WW part is related to the polarized (helicity) gluon PDF ∆G(x), and the genuine twist-three part consists of three-gluon correlators ∼ F F F . Their full expressions can be found in [9]. The G 3T distribution appeared in the previous calculation of the double spin asymmetry A LT in proton-proton collisions p → p ↑ → hX [23]. The cross section formula derived in [23] can be straightforwardly adapted to the case of single spin asymmetry in SIDIS ep ↑ → hX. Writing the hadronic tensor as we find (see (17) and (25) of [23]) where g µν T = g µν − P µ n ν − n µ P ν and ω µν = g µν − P µ n ν are the projectors to the transverse space.g(x) is again a kinematical function with its precise definition given in [9] (see also [22] where it is called ∆G (1) T (x)). M F denotes the three-gluon correlators F F F (see (19) of [23]). The hard part also starts at two loops, S (0) ∼ α 2 s , whose diagrams have the same topology as in the quark-initiated channel. A representative diagram is displayed on the left hand side of Fig. 2, and S (1) is obtained by attaching a gluon to this diagram in all possible ways. The diagram on the right, which is an analog of the right diagram in Fig. 1, does not contribute due to Furry's theorem. Note that in the computation of A LT in Ref. [23], the imaginary phase comes from the definition of ∆G(x) for the longitudinally polarized proton. That is, the non-pole part of the hard kernel was calculated. In the present case, the imaginary phase comes from propagator poles in the hard kernel S µν αβ , and this is why two-loop diagrams are needed. In the WW approximation, we may write and neglect the third line of (32). We thus arrive at a new contribution to SSA of the form (3)

IV. COMPUTATION OF THE HARD PART: QUARK-INITIATED CHANNEL
In this and the next sections, we embark on an analysis of the two-loop diagrams for the quark and gluon initiated channels, respectively. The calculation is rather involved, especially because nontrivial cancellations of infrared divergences are in demand. In the end, we shall have infrared safe formulas that can be straightforwardly evaluated numerically.

A. Quark-fragmenting channel
In Fig. 1, either a quark or a gluon in the final state fragments into the observed hadron. For definiteness, we focus on the former process below. The treatment of the latter is basically analogous, and will be included only in the final formulas. The hard factor S (0) µν for the quark-initiated and quark-fragmenting channel explicitly reads [1] where k and p q are the the momenta of incoming and outgoing quarks, respectively, and Equation (34) represents the sum of 12 = 2 × 3 × 2 diagrams, two of which are shown in Fig. 1. One can easily recognize the part of diagrams each piece of Feynman rules corresponds to. In (35) we have performed a color trace, while the Dirac trace is yet to be done. There are three δ-functions, one for the unobserved gluon in the final state δ((k + q − p q ) 2 ), and the other two come from the poles of internal propagators (denoted by the crosses in Fig. 1). Eventually we shall take the collinear limit k → p 1 ≡ xP in these expressions and introduce shorthand notations p 2 ≡ p 1 + q and l 1 ≡ p 2 − p q , the latter being the momentum of the unobserved gluon in Fig. 1. However, this has to be done with some care because the limit does not commute with the k T -derivative acting on the hard kernel in (26). Let us define where S (0) µν (k) can be read off from (34). We first convert the k T -derivatives of the δ-functions to the x-derivatives as [24] S α and then use integration by parts to shuffle the x-derivatives from the δ-functions to the hard factor S (0) µν . From the term ∂δ(l 2 1 )/∂x in (40), we get a term with ∂g T /∂x and a term with From the term ∂δ (p 2 − l 2 ) 2 /∂x in (41), we get a term with ∂δ(l 2 1 )/∂x, a term with ∂g T /∂x and a term with We further convert ∂δ(l 2 1 )/∂x to a term with ∂g T /∂x and a term like (42). The two resulting terms with ∂g T /∂x cancel. In total, we are led to for which l 1 is fixed through momentum conservation as l 1 = p 2 − p q . A general proof on both QED and QCD gauge invariance of the hadronic tensor (22) was given in [1]. It was also realized that the first two terms in (22) contain infrared divergence separately when the momentum l 2 becomes collinear to the incoming quark line (see Fig. 1), but the divergences cancel exactly. Now that we have written the original formula in a significantly different form (44), it is a nontrivial task to check that (44) is gauge invariant and divergence free. In Appendix A, we show that this is indeed the case, but only after summing all the lines of (44). Knowing where divergences are hidden in intermediate expressions greatly helps a numerical analysis.

B. Calculation of the hard coefficients
With (16), (21) and (44), the polarized cross section (12) takes the following form where the common notationsx ≡ x B /x andẑ ≡ z f /z have been introduced. We have included a flavor summation with explicit charges e 2 f , as commented after (27). Equation (45) contains two δ-function constraints with the first one δ q 2 T Q 2 − · · · originating from δ l 2 1 . Solving the conditions l 2 1 = (p 2 − p q ) 2 = 0 and p 2 q = (p 2 − l 1 ) 2 = 0, we find two roots Recalling the definition The second δ-function sets (p 2 − l 2 ) 2 = 0, which, together with the condition l 2 2 = 0, give two roots We are allowed to perform the l 2 integrals as In the last equality, we have switched to the polar coordinate and changed the integration variable from l 2T to ∆ 2 . This facilitates the computation significantly because we do not have to integrate over rational functions involving square roots.
Next, we compute the Dirac traces using Feyncalc [25] and apply the x-and k T -derivatives to the 3rd, 4th and 5th lines of (45), which have to be done carefully. Note that the x-derivative acts outside the l 2 integral in the 3rd line. We can only evaluate δ (p 2 − l 2 ) 2 and the evaluation of δ l 2 1 cannot be performed before taking the x-derivative. In the 4th (5th) line the x (k α T )-derivative is within the l 2 integral and so both p 2 − l 2 and l 1 are put on-shell after the derivatives are taken.
The subsequent integrals over φ 2 and ∆ 2 are the most cumbersome part of the entire calculation. The nontrivial angular dependence comes from the propagator denominators leading to a cos(φ 1 −φ 2 ) term, while in the numerator, after taking Dirac traces, we are left with powers of cos(φ 1 −φ 2 ) as well as linear terms of the forms sin(φ 2 − Φ S ) and cos(φ 2 − Φ S ) arising from P nl2S T and l 2 · S T , respectively. We list the formulas used to carry out such integrals in Appendix C. After the φ 2 integration, the k = 1, 2, 3, 4 terms are proportional to where 12 = − 21 = 1, and the k = 8, 9 terms are proportional to In the individual lines of (45), the integral over the modulus l 2T has a singularity when l 2T → 0 (or when ∆ 2 → 1), and when a 2 = −1. However, as we will demonstrate in Appendix A, the total expression is finite because of the QCD Ward identity. Therefore, we first compute the φ 2 integrals for each line separately, sum up the results from all the lines and perform the l 2T (∆ 2 ) integration afterwards. One notable feature is that the loop integration yields in principle a different expression for each of the four combinations of the roots (a 1 , a 2 ). However, we have found that after the summation over a 2 the results are independent of a 1 . This is an important consistency check as it effectively ensures that, after all, the split (47) is not necessary and we are back to the ordinary z integral over a complete domain allowed by kinematics.
The above discussion is for the quark-initiated and quark-fragmenting channel. We have repeated the whole procedure for the quark-initiated and gluon-fragmenting channel. Adding the two pieces, we finally arrive at the total result where S k = sin(Φ S − χ) for k = 1, 2, 3, 4 and S k = cos(Φ S − χ) for k = 8, 9. Note that we may substitute x∂g T f /∂x ≈ −∆q f (x) in the above expression. The hard coefficients in the quark-fragmenting (qq) channel are given by and in the gluon-fragmenting (qg) channel by Let us briefly comment on the analytic structure of the above results. The hard kernels depend on the virtual quark propagator 1/(p 1 − l 1 ) 2 ≈ 1/(−2p + 1 l − 1 ). We have parametrized the fragmenting parton momentum as p q = p 2 − l 1 = P h /z and l 1 = P h /z in the q → q and q → g channels, respectively. Since P − h = z f q − and p − 2 = q − , we have l − 1 = (1−ẑ)q − in the former case, and l − 1 =ẑq − in the latter case. This is why the hard cross sections for the q → q and q → g channels contain the factor 1/(1 −ẑ) and 1/ẑ, respectively. When z f 1, both factors 1/(1 −ẑ) = z/(z − z f ) and 1/ẑ = z/z f become large as z is varied between z f and 1. When z f → 1, only the former becomes large around the endpoint z z f , so the q → q channel dominates over the q → g channel. This observation will be confirmed in the later numerical analysis. The denominator q T hints that higher-order corrections will introduce the large Sudakov logarithms ln 2 (Q/q T ) at low q T , whose resummation should be implemented in principle. This is however beyond the scope of this work.

V. COMPUTATION OF THE HARD PART: GLUON INITIATED CHANNEL
The gluon initiated channel is somewhat simpler, since the right diagram in Fig. 2 does not contribute due to Furry's theorem as already pointed out. We thus consider only four diagrams: the left diagram in Fig. 2 and its crossing diagrams with the photon and gluon attachments being interchanged. Considering the quark-fragmenting channel for definiteness, we sum the four diagrams and their complex-conjugates in the form Here k represents the initial gluon momentum, p q is the observed quark, with the unobserved antiquark carrying the momentum k + q − p q (equal to l 1 = p 2 − p q in the collinear limit), and l 2 is the loop momentum. The derivative ∂/∂k λ in (32) can be performed along the steps analogous to the case of the quark initiated channel. Defining where S (0)µναβ (k) can be read off from (58), we find It will be useful to write the 2nd line as Similar to (45) for the quark-initiated channel, the individual lines in (63) contain infrared divergences which must be canceled in the sum over all the lines. We will prove this cancellation in Appendix B.
The hard coefficients can be obtained in complete analogy to the quark initiated channel. We have the same sets of roots (a 1 , a 2 ) as before (see (46) and (48)). As we will show in Appendix B, a divergence comes neither from the p 1 − l 2 propagator (corresponding to the choice a 2 = −1 for the root) nor from the q − l 2 propagator (corresponding to the choice a 2 = 1) in (61). Therefore, the l 2T loop integral is finite, which can be performed analytically. We have also confirmed that, similarly to the previous case, we obtain an expression independent of the choice of the roots for l + 1(a1) after the loop integral and after the sum over the l + 2(a2) roots. All in all, the final result for the cross section can be written in a compact way as with x∂G 3T /∂x ≈ −∆G(x)/2 in the present approximation. The hard coefficients are given explicitly by In this computation we have chosen the quark to be observed in the final state (p q → P h /z) while the antiquark goes unobserved. As a nontrivial check we have verified that taking the antiquark as the observed final state (l 1 → P h /z) and the quark as the unobserved one, we recover exactly the same hard coefficients.

VI. NUMERICAL RESULTS
With all the analytical results presented in the previous sections, we are now ready to make predictions for EIC measurements. Specifically, we will numerically compute the asymmetries from the following definition where dσ(φ h , φ S ) is a shorthand for The numerator of (68) is proportional to the O(α 2 s ) polarized cross section we calculated. In terms of the Fourier coefficients (19), we have with the definition and the abbreviations The integration variables x and z are in the ranges As for the unpolarized cross section in the denominator, we use the leading-order O(α s ) formula [17], summarized by Eqs. (54)-(59) in [19]. The angular decomposition can be cast in the following form Since we integrate over the lepton angle (see below), only the first term F 1 is relevant with the explicit expression where G(x) is the unpolarized gluon PDF and the summation over f includes both quarks and antiquarks. The hard factors are given byσ Using F 1 (75) and the relation (20), we obtain from (68) In practice, we show our results as functions of P hT , z f or x B by integrating out all the other variables in the numerator and denominator. Instead of Q 2 , it is convenient to use y = Q 2 /(x B S ep ) so that we have dQ 2 = x B S ep dy and the relations 1 + cosh 2 ψ = 2 1 + (1 − y) 2 y 2 , There are general kinematical constraints on the integration ranges of x B , z f and y. The condition x min < 1 (see (73)) leads to Requiring the upper bound of x B to be positive, we find a condition on y, Similar constraints can be obtained from z min < 1, which are however not very restrictive because P 2 hT S ep . In actual experiments, A U T is integrated over conveniently chosen bins in x B , z f and y, and we will follow suit.
Since we are using the leading-order cross sections for both the numerator and denominator, one may ask a legitimate question about the effect of higher order corrections, in particular when P hT Q and the resummation of the Sudakov logarithms is required. While such a procedure is well established for unpolarized cross sections, that for transversely polarized cross sections is poorly understood. On a general ground, we expect that the impact of resummation largely cancels in the ratio, but this has to be checked, and will be left for a future work. As for the scale µ of the QCD coupling constant α S (and also of PDFs and FFs), we argue that the lower scale µ = P hT is more appropriate than the larger one µ = Q in the typical kinematic region P hT Q we are considering. This is understood simply from the aspect of the Sudakov (k T ) resummation usually done in the Fourier conjugate impact parameter space b T . The running of the coupling tends to pick up a dominant contribution from the large b T region under the inverse Fourier transformation (for which some prescription is needed to avoid the Landau pole [26]). Therefore, the choice of a lower scale µ = P hT ∼ O(1/b T ) fits the above all-order picture better.
The computation is performed with the most recent NNPDF and JAM global fits. For the NNPDF sets we use the helicity PDFs from [27] and FFs from [28]. For the JAM sets we use the helicity PDFs and FFs from a simulatenous fit in [29]. g T f (x) and G 3T (x) are deduced from the helicity PDFs according to the formulas (25) and (30) in the WW approximation. The uncertainties in PDFs (FFs) in NNPDF and JAM fits are quantified by the Monte Carlo replica method to generate a variance according to the normal distribution. In all the plots below the band represents a combination of 1-σ uncertainty due to the replica method and also uncertainty in the scale choice according to 0.5 < µ/P hT < 2.0, added in quadrature. Note that the edge µ = 0.5P hT starts at P hT = 2 GeV.

A. Calculation for COMPASS kinematics
Though our approach is most naturally and legitimately applied to the kinematics for EIC, let us first present the results relevant to the COMPASS experiment [15]. Admittedly, the collinear factorization may not be applicable to the COMPASS kinematics since most of the data points have P hT below 1 GeV. There is, however, one published data point with P hT ≈ 1.5 GeV. We thus only show the P hT distribution for P hT > 1 GeV, integrating out the other variables over the following coverage [15] 0.003 ≤ x B ≤ 0.7 , as well as where With a 160 GeV muon beam colliding on a fixed proton target, the center of mass energy is S µp ≈ 17.4 GeV. The P hT distributions are shown in Fig. 3 for both π + and π − productions. We see that the Sivers asymmetry for π + is smaller than 0.5% in magnitude using the NNPDF fits and about ∼ 0.5% − 1% in magnitude using the JAM fits. In either case, the sign is opposite to the highest P hT COMPASS data point (see the top-right plot in Fig. 9 of [15]). Although the significant experimental uncertainty makes a meaningful comparison difficult, the result does indicate the importance of other sources of SSA, such as the ETQS function. However, P hT 1 GeV is the borderline between the collinear and TMD approaches. Therefore, our analysis implies that not only the Sivers function but also the new higher-twist contributions found in [1] need to be included in the global determination of nonperturbative inputs in this regime. As for the Collins asymmetry, our result is negligibly small. The data show nonvanishing central values at P hT = 1.5 GeV (see the top-right plot in Fig. 6 of [15]), but they are consistent with zero after the large error bars are taken into account.

B. Calculation for EIC kinematics
We now present our results for the EIC kinematics. Figure 4 shows the z f distribution of the π + Sivers asymmetry for S ep = 45 GeV integrated over the window 0.1 < x B < 0.9 and 0.01 < y < 0.95 and P hT > 1 GeV. The upper bound for the integral over P hT is obtained from (81) by placing the remaining kinematic variables at their extremal values in the above kinematic window. We also impose the conditions Q 2 > 1 GeV 2 and W 2 = (P + q) 2 > 25 GeV 2 . In addition to the total asymmetry ('tot'), respective contributions from different channels (qq, qg, gq) are shown. The asymmetry is largest in the forward region z f → 1, at most 1.5% in magnitude, and decreases towards zero as z f decreases. As we discussed at the end of Section IV, the large z f region is dominated by the quark-fragmenting channel, while the gluon-fragmenting channel becomes important at small z f . Since the final state quark and gluon are back-to-back, this explains the sign change for the two channels. A somewhat larger asymmetry is observed from the JAM fit than from the NNPDF fit. This is in fact a general feature seen also for example in Fig. 3, but most directly understood from the z f -distributions in Fig. 4 where the qg channel contribution is dying off more rapidly as z f → 1 for the JAM fits. Consequently, the cancellation between the qq and the qg channels is less effective using the JAM fits. The underlying reason is the smaller g → π + FF in the JAM fit than in the NNPDF fit. Figure 5 gives the P hT distributions of the π + Sivers asymmetry in low z f (0.05 < z f < 0.4, left) and high z f (0.5 < z f < 0.9, right) bins. We can again see the role of the g → π + FF: at low z f the Sivers asymmetry can even become positive (albeit rather small in magnitude) using the NNPDF fits, while in the large z f region the qg channel quickly dies off so the JAM fits predict a larger (negative) Sivers asymmetry, around 0.5% ∼ 1.5% in magnitude.
Note that the gluon-initiated (gq) channel is negligibly small, almost invisible in the plots. A closer look reveals that the contribution to A N from this channel is less than 10 −3 in all the bins we have studied. We have anticipated that the gluon-initiated channel gives a small contribution to light-hadron production. However, the suppression is stronger than expected, and we attempt to explain the reason in the concluding section.
Further predictions for the P hT distributions of the π + Sivers asymmetry across three bins in x B and Q 2 , using the NNPDF and JAM fits, are exhibited in Fig. 6. We find that the Sivers asymmetry can reach up to 2% in magnitude for the JAM fit covering both large x B (0.1 < x B < 0.7) and moderate x B (0.05 < x B < 0.1) bins for the lowest Q 2 bin. Going from moderate to small x B , the Sivers asymmetry drops to a sub-percent level, as seen in the last x B bin with 0.001 < x B < 0.05. This suppression at small x B in fact has the same origin as the smallness of the gluon-initiated channel mentioned above (see the discussion in the concluding section). Figure 7 covers the Sivers asymmetry for two collision energies (S ep ) and two bins in x B . The results show very mild dropping of the Sivers asymmetry as S ep is increased from 45 GeV to the top EIC energy of 140 GeV (see also [30]). The reason is that the energy dependence mainly comes from the y dependence, which roughly cancels out between the numerator and denominator. Compared with an earlier prediction for EIC in the TMD framework at low momentum P hT < 1 GeV (see Fig. 21 of [31]), our result for the Sivers asymmetry is similar or somewhat smaller in magnitude but opposite in sign, although a detailed comparison is not possible because there is no overlap in the plotted P hT ranges. This suggests that there are cancellations between different mechanisms which may lead to a sign change. However, we emphasize again that when P hT 1 GeV, other sources of asymmetries from various twist-three TMDs found in [1] should be added to the contribution from the Sivers function.
Finally, in Fig. 8, we present a full set of moments introduced in (20) and (79) for three different bins in x B and for fixed bins in Q 2 and z f (1 < Q 2 < 10 GeV 2 , 0.5 < z f < 0.9). We find that two additional moments A sin(φ S ) U T and A sin(2φ h −φ S ) U T reach up to 2% in magnitude in the highest x B bin 0.7 < x B < 0.1. In the TMD framework for lowP hT , the sin(φ S ) and sin(2φ h − φ S ) asymmetries are known to be generated by various twist-three TMDs [32]. We have just demonstrated that the g T distribution (or its TMD counterpart g 1T by extension) is also a potentially significant source of these asymmetries. Indeed, our prediction 1-2 % at P hT = 1 GeV is comparable to previous TMD-based calculations [33,34].

VII. DISCUSSIONS AND CONCLUSIONS
In this paper we have performed the complete analytical and numerical evaluation of the novel two-loop contributions to SSA proportional to the twist-two polarized PDFs ∆q(x) and ∆G(x). Our results indicate that, at the EIC, A U T for pions can reach 1-2% for the three harmonics sin(φ h −φ S ), sin(φ h ) and sin(2φ h −φ S ). On the other hand, contributions from the gluon-initiated channel are negligibly small. Since we are dealing with higher-order perturbative diagrams, we have anticipated that the resulting asymmetry would be small. However, the stronger-than-expected suppression we observed, especially in the gluon sector, calls for an explanation. Parametrically, the asymmetry behaves as s M N P hT (x∆q(x) or x∆G(x)) α s (q(x) or G(x)) .
In addition to the obvious factor of α s , A U T is suppressed by the smallness of polarized PDFs as compared to unpolarized ones. In particular, the gluon-initiated channel is expected to be important for x 1, but there, ∆G(x) ∼ xG(x) as a rule of thumb. On top of this, there is a somewhat unexpected extra factor of x in the numerator which comes from the rewriting dxg T (x) = dx x xg T (x) and dxG 3T (x) = dx x xG 3T (x). Of course the same factor exists in the unpolarized cross section in the denominator, which is, however, accompanied by P + , and the product xP + = p + 1 goes into the hard part and gets absorbed. Therefore, our new contribution, especially in the gluon initiated channel, is strongly suppressed like ∼ x 2 at low x, or in more practical terms, as the selected kinematic bin is sensitive to the low x B region. This tendency has been clearly shown in Fig. 8. In the literature, gluon-initiated channels are usually ignored in the calculation of A U T for light-hadrons (see, however, an attempt in pp collisions [35]), partly because it is believed to be small, but also because nothing is known about the strength of the three-gluon correlators F gF F . For the first time, we have presented a reliably calculable piece of the gluon initiated contributions, and found very small values. After all, our main interest in the gluon initiated processes focuses on A U T in heavy systems such as open charm and J/ψ. This will be studied elsewhere.
It is worthwhile to compare (85) with the well-known parametric estimate of SSA where m q ∼ a few MeV is the current quark mass. This formula has been inferred from the argument in [36], and is often quoted in order to emphasize the smallness of SSA in perturbation theory and the necessity to introduce new nonperturbative distributions. The factor of α s is because one needs loop diagrams such as in Fig. 1 to get an imaginary part, and the factor of m q is because one needs a helicity flip. However, this suppression by m q is illusory for the proton initial state. As is clear from the definitions of g T and G 3T in (23) and (29), m q is replaced by the proton mass M N (see also a related argument in [37]). Thus the correct argument in the DIS case would be that, naively A U T ∼ αsM N P hT is large, but the coefficient is suppressed due to the above-mentioned factor x 2 , resulting in SSA of about 1% as we have shown. In SIDIS at P hT > 1 GeV, this should be comparable to other nonperturbative origins of SSA.
Precisely measuring A U T in the sub-percent region is challenging at the EIC. Conversely, if the future data on A U T turn out to be consistently larger than 1%, most likely genuine twist-three effects are at work. But our result must be subtracted when one tries to extract various twist-three distributions. The distinct kinematical features of our contribution, such as the suppression in low z f and low-x B regions, may be useful to isolate this purely perturbative 'background'. At lower P hT < 1 GeV, predictions based on the Sivers function are available [30,31]. However, in the TMD regime P hT < 1 GeV, there are many other sources of the sin(φ h − φ S ) asymmetry which are unrelated to the Sivers function [1], that must be taken into consideration in order to reliably extract the Sivers function.