Quark Sivers Function at Small $x$: Spin-Dependent Odderon and the Sub-Eikonal Evolution

We apply the formalism developed earlier for studying transverse momentum dependent parton distribution functions (TMDs) at small Bjorken $x$ to construct the small-$x$ asymptotics of the quark Sivers function. First, we explicitly construct the complete fundamental"polarized Wilson line"operator to sub-sub-eikonal order: this object can be used to study a variety of quark TMDs at small-$x$. We then express the quark Sivers function in terms of dipole scattering amplitudes containing various components of the"polarized Wilson line"and show that the dominant (eikonal) term which contributes to the quark Sivers function at small $x$ is the spin-dependent odderon, confirming the recent results of Dong, Zheng and Zhou. Our conclusion is also similar to the case of the gluon Sivers function derived by Boer, Echevarria, Mulders and Zhou (see also the work by Szymanowski and Zhou). We also analyze the sub-eikonal corrections to the quark Sivers function using the constructed"polarized Wilson line"operator. We derive new small-$x$ evolution equations re-summing double-logarithmic powers of $\alpha_s \, \ln^2 (1/x)$ with $\alpha_s$ the strong coupling constant. We solve the corresponding novel evolution equations in the large-$N_c$ limit, obtaining a sub-eikonal correction to the spin-dependent odderon contribution. We conclude that the quark Sivers function at small $x$ receives contributions from two terms and is given by \begin{align} f_{1 \: T}^{\perp \: q} (x, k_T^2) = C_O (x, k_T^2) \, \frac{1}{x} + C_1 (k_T^2) \, \left( \frac{1}{x} \right)^0 + \ldots \end{align} with the function $C_O (x, k_T^2)$ varying slowly with $x$ and the ellipsis denoting the sub-asymptotic and sub-sub-eikonal (order-$x$) corrections.

with the function CO(x, k 2 T ) varying slowly with x and the ellipsis denoting the sub-asymptotic and sub-sub-eikonal (order-x) corrections. The Sivers function [6,7] is crucial for our understanding of the internal structure of the proton (and other hadrons) in terms of its quark and gluon degrees of freedom. It encodes the information about the orbital angular momentum of quarks and gluons in the proton, including spin-orbit coupling effects which lead to an asymmetric distribution of the partonic transverse momentum [8][9][10][11][12][13]. The quark Sivers function is one of the two quark transverse momentum dependent parton distribution functions (TMDs) which is odd under time reversal, the other being the Boer-Mulders function [14]. This leads to the celebrated sign change between the quark Sivers function for semi-inclusive deep inelastic scattering (SIDIS) and Drell-Yan (DY) lepton pair production [15][16][17][18][19] In Sec. III we take the operator definition of the quark Sivers TMD and rewrite it in the small-x/saturation formalism in terms of dipoles containing different parts of the new 'Wilson line' operator. We separately investigate the eikonal and sub-eikonal contributions. In Sec. III A, working at the eikonal level, we show that the leading contribution to the quark Sivers function comes from the eikonal odderon exchange which has a known small-x evolution [50,51,56,58] and asymptotics [54,56]. The odderon intercept is known to be exactly zero at the leading [54,56] and next-to-leading [74] logarithmic approximation in 1/x, and also at any order in the strong coupling α s in the large-N c limit [75]. In addition, the odderon intercept appears to be zero in the strong-coupling calculations based on the anti-de Sitter/Conformal Field Theory (AdS/CFT) correspondence [76][77][78]. We thus conclude, with the perturbative and non-perturbative accuracy of the odderon calculations [54,56,[74][75][76][77][78], that the leading eikonal small-x asymptotics of the quark Sivers function is f ⊥q 1T (x, k 2 T ) ∼ 1/x for both the SIDIS and DY cases. In Sec. III B we explicitly calculate the odderon-generated quark Sivers function in the scalar diquark model of the proton. This contribution dominates at small x as it is proportional to 1/x (see Eq. (66)) and receives no corrections to this power of x from evolution, as described above. It can be compared to the lowest fixed-order (one-loop) contribution to the quark Sivers function in the scalar diquark model of the proton, which, for massless quarks, gives f ⊥q 1T (x, k 2 T ) ∼ x and, therefore, falls off rather rapidly at small x [79]. The sub-eikonal small-x asymptotics of the quark Sivers function is studied in Sec. III C. There we identify the sub-eikonal part of the polarized Wilson line operator (77a) responsible for the Sivers function and construct a new evolution equations (100) (or, at large N c , Eqs. (118) or (121)). These evolution equations re-sum powers of α s ln 2 (1/x): we will refer to this as the double-logarithmic approximation (DLA) (cf., e.g., [1,2,24,31]). Solving these novel evolution equations in the large-N c limit we arrive at the sub-eikonal small-x asymptotic term (129) for the quark Sivers function. Combining this with the spin-dependent odderon we arrive at the eikonal and sub-eikonal small-x asymptotics of the quark Sivers function given in Eq. (1), which we will repeat here for completeness, This is the main physical result of this work. The function C O (x, k 2 T ) is slowly-varying with x and can be exactly determined by the odderon evolution equation (see [65]).
We conclude in Sec. IV by summarizing our results and by describing possible future applications of our technique, such as the determination of the small-x asymptotics for various other TMDs.

II. QUARK S-MATRIX IN THE BACKGROUND FIELD
At small x, the quark and gluon helicity TMD [1] and the quark transversity TMD [2] can be expressed in terms of the 'polarized dipole operators' given by correlators of polarized Wilson lines with the usual eikonal Wilson lines. These polarized dipole operators obey small-x evolution equations, which differ for different TMDs as dictated by the structure of the corresponding part of the polarized Wilson line operator. The polarized fundamental Wilson line up to sub-sub-eikonal order was partially constructed in [2], where, for the transversity calculation, the only contributing portion of the Wilson line operator was the term which coupled the transverse spin of the proton to the transverse spin of a probed quark. Here we want to construct the full polarized Wilson line (quark S-matrix) which applies for scattering of any-polarization quark on any-polarization nucleon, including terms up to (and including) the sub-sub-eikonal order. The calculation will be carried out in the transverse spin basis, with the conversion of our results into the helicity basis being easy to accomplish with the help of Table I below. We consider high-energy scattering of a quark moving in the x − light-cone direction with the large 'minus' momentum p − 2 on a proton moving in the x + light-cone direction with the large 'plus' momentum p + 1 . The lightcone coordinates are defined by such that a space-time 4-vector is written in terms of light-cone coordinates as x µ = (x + , x − , x). The inner product of two 4-vectors is The transverse portion of 4-vectors is denoted by x = (x, y) while its magnitude we is labeled by |x| = x ⊥ . We will make one exception for the transverse momentum k, whose magnitude will be denoted by k T = |k|. The fundamental Wilson line which sums up the eikonal scattering of the quark at transverse position x moving in the x − direction on the gluon field of the proton is [80] with P the path ordering operator, A µ = a A a µ t a the background gluon field, g the strong coupling constant, and the integral running over the light-cone path of the quark. In the case of an infinite path we will abbreviate the Wilson line as We are interested in constructing the generalization of the Wilson line operator which includes all possible sub-eikonal and sub-sub-eikonal interactions of a quark with the target irrespective of the quark and the target polarizations. To this end, we define an S-matrix for the quark-target scattering by where A(p out , p in ) is the scattering amplitude for a quark on a target with p in and p out the incoming and outgoing quark transverse momenta, respectively, while σ and σ are the outgoing and incoming quark polarizations in helicity basis. The amplitude A is normalized such that A = M/(2s) [73], where M is the standard textbook scattering amplitude and s is the center-of-mass energy squared. Extending the notation of [1,2,24], we will denote by V pol x,y;σ ,σ the entire non-eikonal part of the quark scattering S-matrix (7), independent of whether the terms it contains depend on quark or target polarizations or not: x,y;σ ,σ .
From the earlier calculations [1,2,34,35,42] we know that sub-eikonal and sub-sub-eikonal corrections come in as insertions of one or more sub-and/or sub-sub-eikonal operators anywhere along the x − path of the quark. For an insertion of a local operator we have while for an insertion of a bi-local operator (cf. [1,2]) we write Here and below, the two-dimensional integrals denote integration over transverse components of the vector: for instance, d 2 z = dz x dz y . When space allows we will also refer to such integration measure as d 2 z ⊥ . Note that, in general, the operators O pol σ ,σ (z − , z) and O pol σ ,σ (z − 2 , z − 1 ; z 2 , z 1 ) may contain derivatives with respect to z and z 1 , z 2 , respectively, acting on the delta-functions in Eqs. (9) and (10). For the helicity and transversity 'polarized Wilson lines' such derivatives were absent in [1,2] and one also had and V pol for the local and bi-local operators, respectively. For brevity we have suppressed the polarization indices in Eqs. (11) and (12). From the analysis in [1,2,41] we know that the local operator O pol σ ,σ (z − , z) in Eq. (9) can be obtained by calculating the diagrams containing the sub-eikonal and/or sub-sub-eikonal gluon field insertion shown in Fig. 1 below, while the non-local operator O pol σ ,σ (z − 2 , z − 1 ; z 2 , z 1 ) arises due to two quark field insertions with the adjoint Wilson line connecting them, as illustrated diagrammatically in Fig. 2. Note that at the sub-sub-eikonal level the term with two insertions of O pol σ ,σ (z − , z) needs to be included as well. Below we will explicitly calculate the non-eikonal operators O pol σ ,σ (z − , z) and O pol σ ,σ (z − 2 , z − 1 ; z 2 , z 1 ).

A. Gluon insertion operator
Here we calculate the gluon exchange operator O pol G in Feynman gauge (∂ µ A µ = 0). The non-eikonal scattering in Fig. 1 is shown by a gluon coupling to the quark line at the top via a black-circle vertex, where the background field transfers a momentum k to the x − -direction moving quark. Following the approach of [1,2,34], which involves replacing the numerators of quark propagators by the polarization sums, we obtain the (potentially non-eikonal) operator in momentum space The normalizing factor of 1/ 2 p − 2 (p − 2 + k − ) results from the square roots of the denominators of the residues for the two quark propagators adjacent to the non-eikonal vertex integrated over their light-cone 'plus' momentum components, see e.g. Eq. (21) below. The normalization for the operators in [1,2,34] was chosen to be 1/2p − 2 , with the difference between that and the normalization in Eq. (13) not affecting the helicity and transversity operators considered in those references.
FIG. 1: The diagram representing the gluon field contribution to the non-eikonal operator. The black circles on the vertices designate potentially spin-dependent sub-and sub-sub-eikonal scattering.
We take the spinors here for the x − moving quark to be the ± reversed, transversely polarized Brodsky-Lepage (BL) spinors [81] used in [2], with spin quantized along the x-direction. They are related to the helicity-basis ± reversed BL spinors by [21] where χ = ±1 and the helicity-basis ±-reversed BL spinors are [1,2] u σ (p) = 1 with p µ = p 2 +m 2 p − , p − , p and Employing the spinors from Eq. (14) and performing a direct calculation, one can show that the Dirac matrix elements in Eq. (13) arē with S the unit vector in the direction of transverse spin quantization (for the proton polarized along the x-axis we have S =x). The cross-product is defined by Plugging the matrix elements (17) into Eq. (13) we find (for S =x) Since k − ∼ 1/p + 1 is small (with p + 1 the large light-cone momentum component of the parton in the target generating the gluon field), k − p − 2 , we expand some of the terms in the powers of k − /p − 2 , obtaining Note that the "eikonality" of different terms in Eq. (19) is manifest: there is the eikonal A + term, not suppressed by powers of p − 2 , there are the sub-eikonal O(1/p − 2 ) terms, and the sub-sub-eikonal O(1/(p − 2 ) 2 ) terms. There is one caveat about the momentum p − 2 in the sub-eikonal (order-1/p − 2 ) terms in Eq. (19): while in subeikonal calculations one does not need to distinguish the minus momentum component of different quark propagators in Fig. 1, at the sub-sub-eikonal order considered here one has to take into account that the interactions of the quark with the background gluon field result in transfer of a small momentum q − to the quark. If we label P − 2 the momentum of the incoming quark (see Fig. 1), before any interaction, then we can write p − 2 = P − 2 + q − . The difference between 1/p − 2 and 1/P − 2 is of the sub-sub-eikonal order. The difference between 1/(p − 2 ) 2 and 1/(P − 2 ) 2 is of the sub-sub-sub-eikonal order, and is discarded in our calculation. In the coordinate space we replace q − → i∂ − with the derivative acting on everything to the right of (and, hence, earlier than) the insertion of the non-eikonal operator at hand, such that q − is really a total minus momentum transferred to the quark by the target at the light-cone time x − . The order of the operators also matters: we will replace 1/p − 2 → 1/(P − 2 + i∂ − ) with the operator placed to the right of the non-eikonal gluon field, while a similar replacement 1/(p − 2 + k − ) → 1/(P − 2 + i∂ − ) is different from the former by the placement of the operator to the left of non-eikonal field. Strictly-speaking, for consistency we need to expand 1/(P − 2 + i∂ − ) in the powers of i∂ − /P − 2 : we chose to keep the term as is, with the understanding that the overall expression along with this term are correct only up to sub-sub-eikonal accuracy. While the gluon field of the target is often taken to be a function of x − and x only, see e.g. [42], below, in Eqs. (83), we will see that the x + -dependence comes in through the phase of the field, which usually constitutes a higher-order (in eikonality) correction to the leading-order observables generated by the field. Hence, the partial derivative ∂ − with respect to x + is non-trivial, even when applying to a regular light-cone Wilson line (6), if the sub-eikonal x + -dependent phase is included in the otherwise eikonal gluon field A + .
To Fourier-transform the expression (19) into coordinate space we also replace k → −i ∇, k − → i∂ − (both derivatives acting on the gluon field), p 2 → −i ∇, and (p 2 + k) → i ∇, where the operator ∇ acts only on the objects to the right of the operator O G χ ,χ , while ∇ acts only to the left, and neither of them acts directly on A µ , somewhat similar to the notation in [42]. (The ∇ operator is defined by ∇ i = −∂ i = ∂ i .) This yields, with sub-sub-eikonal accuracy, . While the operator in Eq. (20) is also a function of x µ = (x + , x − , x), in the future, by analogy to the light-cone Wilson line (6), we will use , that is, we will just put x + = 0 in all expressions, while implying that the ∂ − -derivatives are applied before putting x + to zero.
The expression (20), while including all the eikonal, sub-eikonal and sub-sub-eikonal corrections to the quark scattering due to an insertion of the background gluon field, is missing the non-eikonal corrections to the free quark propagator. Consider the propagator of the p 2 quark line in Fig. 1, concentrating on the denominator of the propagator. (The numerators in our calculation are written as polarization sums plus the instantaneous terms, just like in the light-cone perturbation theory (LCPT) [81,82].) The integral over p + 2 yields (for p − 2 > 0) where the ellipsis denote the higher-order (sub-sub-eikonal and beyond) corrections, which are obtained by further expanding the exponential in the powers of the phase. Those corrections correspond to multiple insertions of the sub-eikonal operator −i(p 2 2 + m 2 )/(2p − 2 ) in Eq. (21). The sub-eikonal phase corrections were employed two decades ago in the medium jet quenching calculations for heavy ion collisions [83][84][85][86][87].
The above considerations for the sub-sub-eikonal minus momentum transfer apply here as well: therefore, we need to replace 1/p − 2 → 1/(P − 2 + i∂ − ) in the phase operator from Eq. (21). Free quark propagation preserves the quark polarization, and hence the phase operator −i(p 2 [42]). The D i · D i term contains the A i A i product, which diagrammatically arises from the insertion of two non-eikonal gluon fields, with the instantaneous part of the quark propagator inserted between them [42]. Here we have restored this term by simply requiring gauge-covariance of the operator O G χ ,χ , which would lead to gauge-invariance of the quark scattering amplitude. The same philosophy allowed us to replace ∂ − → D − in the denominators, with D µ = ∂ µ − igA µ the covariant derivative. While the terms linear in the field A µ and the terms independent of A µ coming from the first term on the right of Eq. (22) are explicitly present in Eq. (20) already, the remaining terms were again restored by requiring gauge-covariance of the expression. Note that the A i A − A i term appears to require a calculation of a diagram with two consecutive insertions of the instantaneous term in the quark propagator.
In arriving at Eq. (22) we have also completed ∇ × A = −(∂ x A y − ∂ y A x ) → −F 12 , with F µν the full non-Abelian gluon field strength tensor (cf. [1,34]). The F 12 -terms corresponds to the gluon helicity operator found in [1,34], since the δ χ,−χ structure in the transverse spinor basis corresponds to σδ σσ in the helicity spinor basis. Similarly, (22). The reason we can make this replacement is due to the commutator [A − , A i ] being further energy-suppressed, since the transverse components of the gluon field A i are at most sub-eikonal, while A − is sub-sub-eikonal. (This is the case only if the gluon field is emitted by a source, and does not represent the incoming source gluon [41]. In the latter case one needs to perform a calculation with two gluon field insertions and the instantaneous term in the quark propagator we mentioned above, similar to [42]).
The A + -term in Eq. (22) will just yield the usual eikonal Wilson line: we will remove it from the operator O G χ ,χ , labeling the remainder of the operator by O pol G χ ,χ . Finally, replacing all the remaining ∇ i → D i with the sub-sub-eikonal accuracy of our calculation, we arrive at a completely gauge-covariant operator Once again let us point out that the second-to-last term in the curly brackets of Eq. (23) corresponds to transversity [2], while the F 12 -term corresponds to helicity [1,34].
In the helicity basis for spinors the operator (23) can be re-written with the help of the conversion Table I. This yields In the helicity basis Eq. (24) becomes The diagram representing the quark field contribution to the sub-eikonal and sub-sub-eikonal operators.
Now we turn to the two quark exchanges shown in Fig. 2. The part of this operator related to transversity at the sub-sub-eikonal order was calculated in [2]. Here we follow similar steps, but this time retaining all the other terms at the sub-eikonal and sub-sub-eikonal orders. We begin with the left quark exchange in Fig. 2, which, in the A − = 0 gauge, after expansion in k − /p − 2 in some of the terms, gives Here the partial derivatives ∂ z1 and ∂ z2 are acting on U ba z 2 [z − 2 , z − 1 ] δ 2 (z 2 − z 1 ) only, while the partial derivatives with arrows, ∂ z1 and ∂ z2 , are only acting onψ(z − 1 , z 1 ) and ψ(z − 2 , z 2 ), respectively (as indicated by the direction of the arrows). Additionally, α, β = 1, 2, 3, 4 are the Dirac spinor indices.
When comparing Eq. (31) to Eq. (16) in [2] one has to take into account a different definition of the light-cone coordinates (4) used here and that now we do not put p 2 = 0: even then, Eq. (31) differs from Eq. (16) in [2] by several minus signs. In addition, Eq. (16) in [2] was partially projected onto χ = χ , with some transverse polarization-independent terms discarded, as not contributing to transversity, calculating which was the main goal of [2]. Hence it is difficult to fully compare Eq. (31) to Eq. (16) in [2]. However, the transverse spin-dependent parts (∼ χδ χ,χ ) of the two expressions are the same: we, therefore, confirm the polarized Wilson line operator used in the transversity calculation of [2]. The next step is to rewrite the expression (31), derived in the A − = 0 gauge, in a gauge-covariant form, similar to Eq. (23). Following [2], this is easily accomplished by "promoting" the partial derivatives into covariant derivatives, except now the derivatives acting on U ba will become adjoint covariant derivatives. (In the A − = 0 gauge, the A i field of the target is energy-suppressed: therefore, replacing transverse derivatives in Eq. (31) by covariant derivatives constitutes adding sub-sub-sub-eikonal corrections and is allowed by the accuracy of our sub-sub-eikonal approximation.) We thus write acting on spinorsψ and ψ with the color matrix in A i inserted betweenψ α t a and between t b ψ β , respectively, while D ab we recast Eq. (32) as Once again let us stress that the transverse spin-dependent χδ χ,χ term in Eq. (34) is identical to that in Eq. (22) of [2], if one takes into account the difference between the light cone coordinate definitions in that reference and here.
Inserting the operator O pol qq χ ,χ once or twice into the quark S-matrix we obtain up to and including the sub-sub-eikonal order. The corresponding expression in the helicity basis can be constructed with the help of Table I. Since all the derivatives in O pol qq χ ,χ are acting only on the objects within that operator, and not on the functions multiplying it, we can integrate out most of the transverse coordinates in Eq. (35), writing

C. Quark and gluon insertion operator
It is tempting, at this point, to assume that the full sub-sub-eikonal polarized Wilson line is obtained by simply adding equations (23) and (34), and inserting the sum of the two operators once or twice into the Wilson lines, similar to Eqs. (24) and (36), keeping only the terms up to and including the sub-sub-eikonal order.
The diagram representing the quark, anti-quark and sub-eikonal gluon field insertions contributing to the quark S-matrix at the sub-sub-eikonal order.
This is almost correct: the two contributions which would be missed by such a strategy are depicted in Fig. 3 (see also diagram D in Fig. 7 of [2]) and Fig. 4. The first operator, coming from Fig. 3, combines an exchange of the quark and anti-quark, with a sub-eikonal gluon exchange: together this gives a sub-sub-eikonal contribution. The eikonal gluon exchanges between two quark exchanges are already included in Fig. 2 and in the corresponding operator (34) via the adjoint Wilson line U .
The quark exchanges have been found above in Eq. (28). The sub-eikonal gluon exchange with the quark projectile can be obtained from Eq. (23) above. However, we need the expression for the gluon projectile. A calculation along the above lines yields the following sub-eikonal gluon insertion contribution, along with the sub-eikonal phase (cf. [1,39,41]): Here F 12 = a F a 12 T a is the gluon field strength tensor in the adjoint representation, while the adjoint covariant derivatives D ab = ∂ δ ab − gf acb A c and D ab = ∂ δ ab + gf acb A c act on the operators to be added to the left and to the right of the operator (38). Including the sub-eikonal quark exchange contributions we arrive at the following gauge-covariant contribution of the diagrams like that shown in Fig. 3 to the quark S-matrix at the sub-sub-eikonal level: Since the (covariant) derivatives in Eq. (39) are only acting on z 2 , we can integrate out z 1 and z 3 , obtaining The diagram for a double quark exchange following the double anti-quark exchange, with two quark field operators or two anti-quark field operators inserted instead of a quark and anti-quark field operator pair from Fig. 2 inserted twice.
The second operator we need to add, shown in Fig. 4, is similar to a double insertion of the sub-eikonal part of the quark-anti-quark exchange operator in Eq. (34), as given by the second line of Eq. (36), but with the quark and anti-quark fields interchanged in the middle. That is, instead of two quark-anti-quark field operator pairs (ψψ) 2 , in Fig. 4 we have a pair of quark field operators inserted after a pair of anti-quark field operators, ψψψψ, with the appropriate fundamental and adjoint light-cone Wilson lines inserted between the operators. Taking only the sub-eikonal contributions from each pair of (anti-)quark exchanges, a similar calculation to that which yielded Eq. (34) now gives Again, in the last line of Eq. (41) we defined the new operator O pol qqqq χ ,χ , for compactness of the future notation.

E. Full sub-sub-eikonal fundamental polarized Wilson line
We can now assemble all the pieces of the quark S-matrix together to construct the full fundamental polarized Wilson line to the sub-sub-eikonal order. We need to add Eqs. (24), (36), (40), and (41) along with the "interference terms" for the insertions of the O pol G χ ,χ and O pol qq χ ,χ operators. We arrive at with the operators O pol G χ ,χ and O pol qq given by equations (23) and (34), respectively. Note again, that the corresponding expression in the helicity basis can be constructed with the help of Table I.
Equation (42) is the main formal result of this work. It gives the quark S-matrix for the scattering on an arbitrary background quark and gluon fields up to and including the sub-sub-eikonal order. This result can be used to construct small-x asymptotics of all leading-twist quark TMDs, and, probably, of some higher-twist quark distributions as well. Below we will apply (a part of) it to find the small-x asymptotics of the quark Sivers function.

A. Quark Sivers Function at Small x: a General Expression
Having constructed the full polarized Wilson line operator, we can now investigate the small-x evolution of the quark Sivers function. The operator definition for the unpolarized quark TMDs for quarks with longitudinal momentum fraction x and transverse momentum k is (as in [79]) with f q 1 the unintegrated quark distribution, f ⊥ q 1 T the quark Sivers function, M P the proton mass, S P the proton spin, k T = |k|, and U[0, r] the 'staple' gauge link, which we will take here to be the future pointing link used in SIDIS. We work in A − = 0 light cone gauge of the projectile so that the gauge link is just a product of two light-cone Wilson lines, We employ the standard definition of the quark Sivers TMD in Eq. (43). Alternatively, we could calculate the corresponding spin asymmetry A Sivers U T in the SIDIS process at small-x and extract the Sivers TMD by matching that asymmetry onto the TMD factorization expression. In the case of helicity TMDs, the matching on SIDIS was done in [24], leading to the result equivalent to that obtained later in [1] by taking the standard TMD definition and simplifying it at small x. Since the steps in the calculations performed in [24] and [1] are rather general (see also [23]), here we assume that they apply to all TMDs, including the Sivers function, such that the TMD defined in Eq. (43) would appear in the physical processes, such as SIDIS, even at small x.
Using the saturation/CGC formalism, we can rewrite the operator using semi-classical averaging in the proton's wave function as (see [1] and Appendix A of [37]) where we sum over a complete set of partonic states |X , the large angle brackets denote averaging over the proton's wave function [37,[88][89][90][91][92][93], and we have semi-infinite Wilson lines connecting the positions of the quark field operators to infinity. We can take the sum over the intermediate states to be a final state cut, and represent the possible contributions diagrammatically in Fig. 5 as in [1] at the lowest perturbative order (order-α s ). The contributions are organized by whether ζ − and ξ − are positive, negative, or zero, with x − = 0 being the position of the shock wave (assumed to be very thin in the x − direction). For the left-right asymmetric diagrams B, E, and F, their mirror images need to be added: those are not shown in Fig. 5. Each diagram class is represented only by one diagram in Fig. 5. Note that the diagrams in class C also include interactions with the gluons coming from the shock wave and a long-lived anti-quark produced in the final state (see Fig. 6 below): such diagrams are also not shown for brevity in Fig. 5.
At the eikonal order, only the diagrams B, C and D contribute both to the unpolarized quark TMD and to the Sivers function (see e.g. the unpolarized quark TMD f q 1 (x, k 2 T ) calculations in [23,94]). The diagram D does not include any interaction with the target: hence it does not depend on the spin of the target, and cannot contribute to the Sivers function. We are left with the diagrams B and C at the eikonal level, which are shown in Fig. 6 in more detail below.
At the sub-eikonal order, the contributing diagrams for helicity TMD was found in [1] to be the class B diagrams, which is shown in more detail in Fig. 8 below. At that order the diagrams of classes C, D, and F do not contribute to the TMDs dependent on the target proton polarizations, such as helicity, transversity, and the Sivers function at hand, as argued in [1]: for class C the sub-eikonal target interactions with the anti-quark line cancel when the quark (or gluon) exchanges are taken across the final state cut, the diagram D has no spin-dependent contributions since it contains no interaction with the target shockwave, and the diagrams in class F are further energy suppressed since the gluon emission and absorption take place inside the shockwave. The usual leading-order diagram calculation for the Sivers function [15,16,19,79,95,96] comes from diagrams of class A, taking the interference between a re-scattering inside the shockwave and a purely tree level scattering to generate the phase needed to make the transverse spin dependence real. The resulting lowest-order Sivers function is sub-sub-eikonal, , and will not be considered in our calculations below, which are done at the eikonal and sub-eikonal orders. Such a lowest-order diagram can be evolved by adding soft gluon emissions to both class A diagrams and class E diagrams, as depicted in Fig. 5. However, in Appendix A of [1] it was shown that the contributions of graphs A and E to the double logarithmic approximation (DLA) evolution cancel at the sub-eikonal order. (The double logarithmic approximation is defined as the re-summation of the parameter α s ln 2 (1/x) at small x.) We conclude that the sub-sub-eikonal x-dependence of the lowest-order quark Sivers function, f ⊥ q 1 T (x, k 2 T ) ∼ x, is unaffected by the evolution in the DLA, and persists down to small x for the contributions to the Sivers function coming from the diagrams A and E. However, such a sub-sub-eikonal contribution is much smaller than the eikonal and sub-eikonal contributions of interest to us here and will not be considered further.
To summarize, in the eikonal calculations of the quark Sivers function we only need to consider diagrams B and C, while for the sub-eikonal calculations we only need the diagram B. We begin with the eikonal calculations of the diagrams B and C, depicted in Fig. 6 in more detail. 6: Diagrams of classes B and C with kinematics specified. In the diagram B the antiquark propagates from ζ with momentum k 1 , undergoes a transverse spin-dependent interaction with the proton target at the transverse position w, then propagates to ξ with momentum k 2 . In the diagram C the anti-quark interacts with the shock wave again to the right of the cut, before arriving at the point ξ with momentum k 2 .
We begin with the class B diagram in Fig. 6, where the light-cone Wilson line crosses the shockwave from ζ to the final state to the left of the cut, but connects ξ to the final state on the right side of the cut without crossing the shock wave, and is, hence, an identity on the right of the cut. An antiquark with momentum k 1 propagates from ζ through the shockwave, undergoing an interaction at transverse position w, then propagating with momentum k 2 to the final state and connecting on the other side of the cut to ξ. Taking the shockwave to be at light cone time x − = 0, the class B diagrams have ζ − < 0, ξ − > 0 (as in Fig. 6) or ζ − > 0, ξ − < 0 (in the conjugate diagrams of those like B from Fig. 6). For the B-class diagrams the quark field operator at ζ is before the shockwave and the anti-quark operator at ξ is after the shockwave on the other side of the cut, allowing us to write the scattering as a dipole composed of a fully infinite Wilson lines connecting the quark field operators at ξ and at ζ. To evaluate the diagram B we will need the quark propagator through the shock wave [1], where the interaction of the anti-quark with the shock wave is denoted by the Dirac and color matrix V † w ji . Here α, β are the Dirac spinor indices, while i, j are the quark color indices. For simplicity we assume that the anti-quark is massless, m = 0, since the light quark mass is usually negligible in the Sivers function calculations.
Using the propagator (45) in Eq. (44), integrating over ξ, k 2 , ζ − and ξ − , and inserting spinor polarization sums [1], we see that the diagram B and its complex conjugate where c.c. denotes the complex conjugate terms and the intermediate state |X became an anti-quark state |q as only an anti-quark is produced in the final state.
Using the ± reversed transverse Brodsky-Lepage spinors Eq. (15) we evaluate the matrix element for massless anti-quarks asv For the eikonal calculation at hand we only need the Wilson line contribution to the interaction of the anti-quark with the target. We, therefore, writē with the ellipsis denoting the sub-eikonal terms and beyond. Note that the Sivers function couples an unpolarized quark to the transverse polarization of the proton: hence we only need the unpolarized quark interaction with the shock wave. This means we will only need the δ χ1,χ2 term on the right of Eq. (48), even at the sub-eikonal order we consider below. Employing Eqs. (47) and (48) in Eq. (46) and summing over transverse polarizations χ 1 , χ 2 we arrive at where we have explicitly added in the complex conjugate term. Diagram C from Fig. 6 can be calculated similarly, except all the interactions with the anti-quark cancel to the left and to the right of the cut, such that one needs to use the free anti-quark propagator, which can be obtained from Eq. (45) by using V † w ji = γ − δ ij in it. In the end one obtains Combining Eqs. (49) and (50), replacing ξ → w in the latter, we arrive at We have also inserted a time-ordering sign into the last correlator, as justified in [1] (see also [97]).
To extract the Sivers function, which is odd under time reversal, from Eq. (51) we need the part of its right-hand side which changes sign under k → −k. A quick inspection, along with taking k 1 → −k 1 , ζ ↔ w on the right-hand side of the expression, shows that the part of the correlators in Eq. (51) contributing to the Sivers function is an odd function under ζ ↔ w. One can write [56,58] 1 with S ζw the ζ ↔ w-symmetric, C-even Pomeron exchange operator and O ζw the ζ ↔ w-antisymmetric, C-odd odderon. The combination we have in Eq. (51) is The odderon amplitudes being antisymmetric under the ζ ↔ w exchange, O ζw = −O wζ , are the only ones contributing to the Sivers function. We thus write This is our main result for the eikonal-order quark Sivers function at small x.
The odderon is a correlator of eikonal Wilson lines, so the small-x Sivers function we have found is proportional to 1/x, and should actually grow as we decrease x. Moreover, small-x evolution likely leaves this growth unchanged, as the odderon is known to have an intercept equal to zero at the leading [54,56,58] and next-to-leading [74] orders in α s (see also [53] for an earlier solution with a smaller intercept), at all orders in α s in the large N c limit [75,98], and at strong coupling in N = 4 supersymmetric Yang-Mills theory [76][77][78].
Similar conclusion of the odderon dominance in the small-x quark Sivers function was reached in [3] using the expression for the small-x unpolarized quark TMD (unintegrated quark distribution) from [94,99] (see also [23,100]). While qualitatively our result and that in [3] agree on the odderon-driven quark Sivers TMD, a more detailed comparison between the two calculations is left for future work. The importance of C-odd correlators in the Sivers asymmetry was also pointed out in [101,102].

B. Quark Sivers Function at the Eikonal Level: the Spin-Dependent Odderon Contribution
We have shown that the eikonal small-x evolution for the Sivers function in the operator treatment comes from the spin-dependent odderon. But we have not established that the odderon contribution is nonzero. In [5] it was shown that the odderon survives phase space integration if the proton has an asymmetric parton distribution using a diquark model calculation. Here we will show that the odderon is able to generate a nonzero Sivers function in the same diquark model of the proton. The Lagrangian for this model has a spinor field ψ P for the (point-like) proton, the usual spinor fields for the quarks ψ q , and a complex scalar diquark field ϕ which has mass M roughly equal to the proton mass M ≈ M P and the color quantum numbers of an antiquark. The interaction term between these fields is a Yukawa coupling L int = G ϕ * iψi q ψ P +c.c., with effective coupling constant G for the splitting of the proton into a quark and a diquark. We begin by constructing the odderon amplitude O ζw to be used in Eq. (54). At the lowest order it is illustrated in Fig. 7, with the proton target splitting into a quark (solid line) and diquark (dashed line) in order to interact with the odderon (the 3-gluon exchange at the lowest order). The odderon comes from summing over all possible connections of the three gluons in Fig. 7, where the polarized 'target' proton (thicker solid line) splits into a dipole consisting of a quark and diquark, which exchanges three gluons in the symmetric color configuration d abc = 2Tr[t a , {t b , t c }] with the 'projectile' quark-antiquark dipole. We take the 'target' quark to be at transverse position 0 and the diquark at r, as shown in Fig. 7. We want the coupling of the odderon to the transverse spin of the proton, so we need to insert the odderon exchange between the light-cone wave functions of the proton splitting into the quark-diquark pair, and the same wave function, but complex conjugate, describing the merger of the pair back into the proton. The light-cone wave function in the transverse spin basis for a proton at transverse position u and with transverse polarization X, with the quark carrying the fraction γ of the proton p + 1 momentum and the transverse polarization X (as labeled in Fig. 7), is [43] ψ X,X (r, 0, u, γ) = in the limit of massless light quarks (cf., e.g., [79]). The odderon exchange amplitude in the dipole-dipole scattering depicted in Fig. 7 is [21,56] with [21,[58][59][60] One obtains the spin-dependent odderon amplitude O ζw by convoluting a square of the wave function (55) with the three-gluon exchange amplitude (56). Since the amplitude (56) is odd under the r ↔ 0 interchange, only the ∼ r part of the wave function squared contributes: this is exactly the part of the wave function dependent on the proton polarization X [43]. Putting the wave function squared and the amplitude (56) together we have Substituting the odderon amplitude (58) into Eq. (54) we get where the factor of S P comes from the spin quantization axis S multiplied by the proton polarization X. For simplicity we assume that the r-integral is dominated by the perturbatively short distances,m γ r ⊥ 1, such that we can expand the modified Bessel functions to the lowest non-trivial order and obtain In arriving at Eq. (60) we have also integrated over k − 1 .
Next we argue that ln(1/|r|m γ ) is a slowly-varying function compared to powers and exponentials of r ⊥ present in the integrand and approximate this logarithm by one, ln(1/|r|m γ ) ≈ 1. The integral over r then becomes straightforward: performing it and the integral over x we find Next we integrate over l 1 and l 2 while regulating the singularities by the infrared (IR) cutoff Λ. We get To integrate over k 1 analytically we, again, have to neglect a logarithm, by putting ln[(k + k 1 ) 2 /Λ 2 ] ≈ 1. The remaining k 1 -integral is dominated by the singularity at k 1 = −k, which we regulate by the IR cutoff Λ. Keeping only the most singular part in the expansion in 1/Λ we obtain The IR cutoff Λ can, in principle, depend on γ since the natural mass scale for an IR cutoff would bem γ . However, at sufficiently small γ one must replace this cutoff with the QCD confinement scale Λ QCD , since, formally, in the quark-diquark model the transverse distances are bound by r 1/m γ and may get much larger than the confinement scale as γ → 0. For simplicity we take the IR cutoff to be some momentum scale independent of γ, as this will only change an overall numerical factor. Takingm γ = γM P , and integrating over γ we have the Sivers function due to the odderon exchange in the quark-diquark model We have obtained a nonzero, eikonal contribution to the Sivers function corresponding to the same spin-dependent odderon as studied in [5]. Both the eikonality and the dependence on k and M P in Eq. (66) are different from the Born-level results calculated in the same diquark model (see Eq. (A8) in [79]). In particular, the Sivers function in Eq. (66) has a very interesting behavior in the Λ QCD , M P → 0 limit. Since the IR cutoff Λ must be proportional to Λ QCD , we conclude that Λ ∼ Λ QCD ∼ M P . Therefore, M 2 P /Λ 2 ratio will remain constant in the Λ QCD , M P → 0 limit and the Sivers function Eq. (66) will not vanish in the limit of zero proton mass. This may be a feature unique to the quark-diquark model of the proton, but it calls for further analysis.

C. Sivers Function at the Sub-Eikonal Level: a New Evolution
In this Section we will construct a sub-eikonal correction to Eq. (66). Part of the motivation for this is that the effects of the odderon have been historically hard to find in the data. The recent announcement of the odderon discovery [62] was many years in the making and required very careful measurements and extrapolation in energy. It is, therefore, important to understand the background for the odderon contribution (66) in order to discover the latter in the future measurements of the Sivers function. In addition, the existing data on the hadronic single transverse spin asymmetry (STSA) A N measured in p ↑ + p collisions [103][104][105][106] exhibits no evidence for the odderon (with the possible exception of the AnDY Collaboration data [107]): assuming that the x-dependence of the Sivers function is related to that in A N , that is, that A N ∼ x f ⊥ 1 T , we see that the odderon contribution (66) would predict an x-independent A N . However, most of the data [103][104][105][106] show A N which falls off with decreasing x of the transversely polarized proton. This behavior of the data can be roughly approximated as A N ∼ x, which translates into the Sivers function f ⊥ 1 T being independent of x. Since f ⊥ 1 T ∼ x 0 is the sub-eikonal scaling, it is clear that a sub-eikonal study of the Sivers function is in order.

Sub-Eikonal Sivers Function
The analysis of Sec. III A also applies at the sub-eikonal level [1]. At the sub-eikonal level only the diagram B from Fig. 8 contributes [1]. One essential difference we have here as compared to, say, the helicity calculation in [1], is that the sub-eikonal and sub-sub-eikonal terms in the quark S-matrix operator (42) are non-local in the transverse plane. Thus, the calculation in Sec. III A has to be augmented by introducing two different coordinates of the anti-quark to the left and to the right of the shock wave, w and z, as shown in Fig. 8. This is taken into account by replacing along with e ik·(w−ζ) → e ik·(z−ζ) in Eq. (49). We arrive at the Sivers function given by the diagram B and its complex conjugate, where, in the last step, we have extracted the DLA contribution only, and defined with At the sub-eikonal order, the polarized Wilson line in Eq. (68) contains only the ∼ δ χ,χ contribution, since the Sivers function couples an unpolarized quark to the proton spin. 3 We thus write, neglecting the quark mass (m = 0) and employing Eqs. (23) and (34) in Eq. (42), truncated to the sub-eikonal order, keeping the δ χ,χ term only, and neglecting the eikonal contribution already accounted for above, In the A − = 0 gauge we are working in, the transverse components of the gluon field, A, are sub-eikonal. Using this to simplify the gluon part of Eq. (70), while keeping only sub-eikonal terms dependent on the target transverse polarization, we get Further, integrating over z yields We see that while the quark sector operator in Eq. (72) is local in the transverse plane, the gluon sector operator is non-local due to the derivative of the δ-function. Substituting the gluon part of Eq. (72) into the first line of Eq. (68) we arrive at, after multiple integrations by parts, It appears that the non-locality of the gluon operator in Eq. (72) translates into the factors of k 1 and k in Eq. (73). We conclude that the sub-eikonal contribution to the quark Sivers function is The part of the expression on the right-hand side of Eq. (74) that contributes to the Sivers function should change sign after k → −k. Hence, to make (74) an equality, we need to anti-symmetrize its right-hand side under k → −k. Simultaneously changing k 1 → −k 1 and w ↔ ζ we arrive at Further, define Equation (76) becomes Next we define two polarized dipole amplitudes Apart from the transverse positions of the Wilson lines, the amplitudes depend on the longitudinal momentum fraction z, which can be roughly thought of as the smallest of the momentum fractions of the quark and anti-quark lines (see [24,46] for a more precise definition of the argument z). Note that w,ζ (z) = −F [2] ζ,w (z).
The DLA small-x evolution of a sub-eikonal operator only couples it to sub-eikonal operators at the next step, otherwise the evolution would not generate longitudinal logarithms of energy. Hence, the small-x evolution of the operator(s) in Eqs. (75) and (77) (or, equivalently, the evolution for the polarized dipole amplitudes in Eqs. (80)) will only couple to the same operators.

Initial Conditions
For a single quark target, the lowest-order gluon field dependent on the transverse polarization of the quark at b = 0, b − = 0 and an arbitrary b + position in A − = 0 gauge is where X and X are the quark's polarizations before and after the gluon field emission, respectively, and M P is the "quark" mass. This field is indeed sub-eikonal (∼ 1/p + 1 ) at the leading order and may contribute to the operators in Eqs. (77). Note that it does not contribute to the F 12 sub-eikonal term in Eq. (23). Hence, the transverse spin dependence at the sub-eikonal level only contributes to the gluon operator in Eq. (70), and, consequently, the operators in Eqs. (75) and (77).
Note that the field (83) depends on x + and on the rapidly varying b + (with b + even not specified for a quasiclassical target [88][89][90]). However, this dependence only appears at the sub-sub-eikonal order, and, even at that order, it does not affect the field strength tensor F −i needed for the transversity operator in Eq. (23). The x +dependence may appear in any gluon field if we take into account the energy-suppressed phase, as in Eq. (83). However, at the leading orders in eikonality in each channel, in the operators entering Eq. (23), the x + -dependence does not appear.
Another puzzling feature of the field in Eqs. (83) is that the sub-eikonal (∼ 1/p + 1 ) terms in it appear to come in with a factor of i, that is, they are imaginary. This feature, while requiring further interpretation in the future, can be attributed to the conventional wisdom that the transverse spin dependence, Xδ X,X , usually comes in with a factor of i associated with it. This is the well-known i which needs a complex phase to give a real contribution to the Sivers function [15,16,19,110]. Here it comes in through the transverse spin-dependent part of the gluon field.
To construct the initial conditions for our evolution, we will work with the single quark target. Substituting the field from Eq. (83b) along with the polarization-independent eikonal field into Eq. (80a) we immediately see that the two-gluon exchange contribution vanishes, and the first non-trivial polarized dipole amplitudes arise at the three-gluon level. We obtain, after integrating over the impact parameters b between the dipole and the target quark, with c 0 given by Eq. (57) above. Note that the three gluons have to be in the d abc color state, similar to the odderon. One may, therefore, think of Eq. (85) as of the lowest-order contribution to the sub-eikonal spin-dependent odderon. Similar calculation of the initial condition for F [2] w,ζ from Eq. (80b) due to the fields in Eqs. (83b) and (84) of a single quark target readily gives zero in the gluon sector to any order in the gluon exchanges. There is also a quark sector, as one can see from Eq. (77b). However, the quark part of the operator in Eq. (77b) contains only the γ + Dirac matrix and, like other operators in Eqs. (77), is local in the transverse plane: therefore, this operator is similar to that in the unpolarized quark parton distribution function (PDF), and cannot couple to the transverse spin of the target proton. Therefore, it appears that one can discard the quark part of the operator (77b), both in the initial conditions and in evolution.
This can also be seen by the order-by-order evaluation of the quark part of F [2] w,ζ . At the lowest order, the quark operator in Eq. (77b) comes in with a two-quark exchange in the t-channel: this contribution to F [2] w,ζ in Eq. (80b) is zero, since, after averaging over the quark impact parameters, each trace in Eq. (80b) will become transverse coordinate independent and the difference of the two traces will be zero. Adding eikonal gluon exchanges will make each trace in Eq. (80b) a function of w − ζ. However, since the quark operator in Eq. (77b) is local in the transverse plane and comes in with γ + only, it will not generate any dependence on the transverse spin, such that each trace in Eq. (80b) will be a function of |w − ζ|, and their difference will again cancel.
We thus conclude that the initial condition for the second polarized dipole amplitude is zero, Below we will show that the DLA evolution of F [2] w,ζ is closed in the gluon sector, it does not mix with F i w,ζ : therefore, the zero initial conditions (86) imply that F [2] w,ζ = 0 even after evolution.

Small-x Evolution: General Expression
The small-x evolution of the polarized dipoles in Eqs. (80) will ultimately be calculated in the large-N c limit. We will, therefore, assume that the evolution is gluon-driven and neglect the quark field insertion operator in the second line of Eq. (77b). (The same approximation was done for the large-N c limit of the small-x helicity evolution in [24]: see also a discussion of this approximation in [40].) Above we have argued that the initial condition for the entire F [2] w,ζ is zero: the only way it may become non-zero is by mixing with F i w,ζ through evolution. The mixing of the quark part of (77b) in F [2] w,ζ with F i w,ζ , even if non-zero, has to include an (eikonal) interaction of the unpolarized (anti-)quark in the dipole with the target, otherwise the contribution to the dipole amplitude F [2] w,ζ from Eq. (80b) will be zero, since the latter has to be an odd function under w ↔ ζ interchange, and, therefore, has to be a function of both transverse positions w and ζ. An interaction of the unpolarized (anti-)quark in the original parent dipole with the target would make the quark operator evolution sub-leading in N c and can be discarded. Therefore, we will neglect the quark operator from the second line of Eq. (77b) in our evolution analysis below.
To construct the small-x evolution of the operators in Eqs. (77) we follow [80] and rewrite the gluon field as a sum of the background field B µ and the quantum field a µ , and integrate out the quantum fields a µ . To do so, we will need the propagator with the sub-eikonal polarized gluon Wilson line operator (cf. Eqs. (38) and (72)) The propagator is obtained by the technique outlined in [34], employing the δ λ,λ part of the the operator in Eq. (38) to describe the interaction with the shock wave. While, strictly speaking, the gluon field on the right-hand side of Eq. (89) should be B µ , we denote it A µ since at the next step of evolution it will again be separated into the background and quantum fields. Substituting Eq. (89) with s = p + 1 k − into Eq. (88) and integrating over One will also need the standard eikonal propagator a + a + , which is the same as for the unpolarized evolution, We begin with the gluon part of the V [2] pol w from Eq. (77b), which we rewrite as Substituting Eq. (92) into Eq. (80b) and employing the decomposition (87) we obtain a number of contractions. They are diagrammatically represented in Fig. 9, where only representative graphs from most diagram classes are shown and only for the first term in Eq. (92), for brevity.
Since the initial condition for F [2] w,ζ is zero, per Eq. (86), to get a non-zero F [2] w,ζ we need to find evolution steps mixing it with F i w,ζ . To this end we note that in the regime opposite to that in Eq. (90), that is, for The first term in the curly brackets of Eq. (93) has a different sign compared to that in Eq. (90). The second terms are the same in both equations. Since we are interested in evolution mixing F [2] w,ζ with F i w,ζ , below, when discussing the evolution of F [2] w,ζ , we will only talk about the second term in Eqs. (90) and (93), which is the same in both expressions.
To give a longitudinal logarithmic integral dk − /k − , the ⊥ + gluon propagator, like those in Eqs. (90) and (93), has to cross the shock wave [24]: hence, only the diagrams with the gluon crossing the shock wave are shown in Fig. 9 for the ⊥ + propagator. The quantum field a µ can only be involved in contractions (see, e.g., diagram I in Fig. 9), as it is integrated out. The B µ -field outside the shock wave is put to zero. This way the diagrams II, III and V, along with other similar diagrams where we have a factor of the gluon field "outside" the shock wave and of the gluon propagator, are all zero. One can further argue that the diagrams I, IV, and VI cancel between the two terms in Eq. (92). Once again we are talking about the term in the ⊥ + gluon propagator which is the same in Eqs. (90) and (93): the other terms do not cancel, but do not mix F [2] w,ζ with F i w,ζ . This leaves us with the diagrams VII and VIII, where the gluon emission is of the eikonal type, similar to the unpolarized small-x evolution [80,93,[111][112][113][114][115]. These diagrams represent all possible eikonal gluon emissions and absorptions. Eikonal gluons do not have to cross the shock wave to generate longitudinal logarithms: hence, diagram VIII is allowed, along with other virtual diagrams. The eikonal evolution leaves the operator (92) intact in the shock wave, as denoted by the two adjacent boxes inside the shock wave shown in the diagrams VII and VIII: again, no mixing with F i w,ζ is generated. Summarizing the analysis we have just made, we see that the DLA evolution for the dipole amplitude F [2] w,ζ in the large-N c limit couples only to the same amplitudes F [2] w,ζ and does not mix with F i w,ζ . That is, the evolution is of the following type, with K some integral kernel. Since, by Eq. (86), the initial condition (the inhomogeneous term in Eq. (94)) for this evolution is zero, F [2](0) w,ζ = 0, this means that F [2] w,ζ = 0 (95) with the DLA and large-N c accuracy. We thus discard the dipole amplitude F [2] w,ζ and proceed by constructing the evolution for F i w,ζ . The diagrams contributing to the evolution of the polarized dipole amplitude F i w,ζ from (80a) are shown in Fig. 10. Again, we do not show all the eikonal diagrams explicitly. For brevity, we also depict the diagrams contributing only to one of the traces in Eq. (80a). We begin by calculating the sub-eikonal diagrams α, β, γ, and δ in Fig. 10.
By analogy to Eq. (77a), we define an adjoint polarized Wilson line of the same type by With the help of this definition, repeating the steps outlines in more detail in [1,34], we arrive at the following contributions of the diagrams α and β Similarly, diagrams γ and δ give In obtaining Eqs. (97) and (98) we have defined z = k − /p − 2 with the upper cutoff on the z integral given by z, the minus momentum fraction at the previous step of the evolution. The lower limit of the z integral, Λ 2 /s, involves an IR cutoff Λ for the transverse momenta. We have also been using an abbreviated notation for the light-cone Wilson lines V 1 = V x 1 , U 2 = U x 2 , etc. Transverse vectors are defined by x ij = x i − x j with x ij = |x ij |. Most importantly, in deriving Eqs. (97) and (98) we have neglected the first term in the propagator (90) as being the gluon analogue of the (gluon part of the) operator (77b), which we have established to be zero, per Eq. (95).
The eikonal diagrams contribution is well-known. The diagrams , ζ, etc., yield [1,34,80,93,[111][112][113][114][115] Combining equations (97), (98) and (99) we arrive at the evolution equation for the dipole amplitude F i 10 in the operator form, with the inhomogeneous term given by Eq. (85) (for the impact-parameter integrated version of Eq. (100)). Here +(1 ↔ 0) applies to everything on the right, except for the inhomogeneous term. We have also dropped the time-ordering signs for brevity: they are still implied in all correlation functions.
Combining Eqs. (100) and (95) with Eq. (82) we write the sub-eikonal contribution to the Sivers function as

Small-x Evolution in the Large-Nc Limit
Similar to the unpolarized evolution [80,93,[114][115][116][117][118][119][120][121] the evolution equation (100) is not closed: not all the operators on its right-hand side are the same as the dipole operator on its left-hand side (see also (80a)). To obtain a closed evolution equation we will take the large-N c limit. Employing along with the Fierz identity, one can readily show that (cf. [1]) Using Eqs. (102) and (103) in Eq. (100) we obtain where the ellipsis denote the N c -suppressed terms. Taking the large-N c limit yields where all the correlators on the right are functions of z , with this dependence not shown explicitly. We are interested in the solution of this equation outside the saturation region: therefore, we linearize it by replacing tr V 0 V † 2 and tr V 2 V † 1 by N c . This gives the linearized evolution − ln 1 We are interested in the double-logarithmic evolution, for which the transverse integrals in the kernel reduce to logarithms [24]. Whether the transverse integrals are logarithmic in the IR or in the ultraviolet (UV) depends, for most of the terms, on the transverse distance dependence of the correlators on the right-hand side of Eq. (106). Inspired by the inhomogeneous term (85), we will assume that where the function F (x 2 10 , z) may include logarithms of x 2 10 or perturbatively small powers of x 2 10 (e.g. x const √ αs 10 ), along with the z-dependence, but no order-one powers of x 2 10 . Note that Eq. (107) along with Eq. (80a) imply that For the correlators scaling as shown in Eq. (107) we can analyze different terms in the kernel of Eq. (106), extracting the logarithmic contribution.
• The eikonal kernel x 2 10 /(x 2 21 x 2 20 ) is logarithmic in the UV when x 2 → x 1 and in the IR when x 21 ≈ x 20 x 10 . It is not logarithmic when x 2 → x 0 , since the two terms in the square brackets multiplying this kernel in Eq. (106) cancel. We thus approximate, with the DLA accuracy and after integrating over the impact parameters, where the last term on the right is logarithmic because we assume that the correlator multiplying is scales ∼ x 2 21 , per Eq. (107). The IR (upper) limit of the integral in the second term comes from the light-cone (x − -) lifetime ordering condition [24,39], which is essential for the DLA. Similarly, the UV limit of the x 2 21 -integral in the first term on the right of Eq. (109) comes from requiring that the emitted gluon's lifetime is longer than the extent of the shock wave, z x 2 21 > 1/s. One may also think of 1/(z s) as the shortest transverse distance in the problem. There is one important caveat left. Consider the first trace in the last line of Eq. (109). It describes the amplitude in the dipole 10. The subsequent emissions in that dipole will have their lifetimes capped by z x 2 21 , with x 21 being the size of the "neighbor" dipole, not the one we continue the evolution in. We, therefore, define the "neighbor" dipole amplitude Γ(x 2 10 , x 2 21 , z) [24,31,39] by where the lifetime dependence is shown explicitly in the argument. Equation (109) becomes • Next consider the first expression in the second term in the kernel of Eq. (106), that is, This kernel has no UV divergences at either x 2 → x 1 or x 2 → x 0 , since the potentially divergent terms vanish after angular averaging. (Strictly-speaking angular averaging would lead to delta-functions δ 2 (x 21 ) and δ 2 (x 20 ) in the kernel: however, zero-size daughter dipoles generated this way would have zero lifetimes, or, more precisely, the delta-functions would imply that z x 2 21 = 1/s or z x 2 20 = 1/s, not allowing for any further DLA evolution due to impossibility of imposing lifetime ordering like (110) beyond zero lifetime: such contribution may need to be added to the inhomogeneous term, but is not included in the evolution.) Note that employing Eq. (107) one can show that ln(1/Λ) vanishes in Eq. (113): hence the logarithm present in the integrand does not necessarily make the result of the integration logarithmic. Lastly, to determine the IR asymptotics, x 21 ≈ x 20 x 10 , one has to split Eq. (113): in each term we expand in x 10 /x 21 or x 10 /x 20 and average over the angles of x 21 or x 20 , depending on whether the impact-parameter integrated correlator depends on x 2 21 or x 2 20 , respectively. We thus get (after the impact parameter integration) such that Eq. (113) in DLA approximates to F (x 2 20 , z ).
• Finally, let us consider the last term in the kernel of Eq. (106), This term contains a UV divergence at x 2 → x 1 , for the second term in the square brackets. There is no UV divergence at x 2 → x 0 . The contribution coming from the IR region, x 21 ≈ x 20 x 10 , can be evaluated using the above technique. In the end we obtain ≈ −π Using (124) in Eqs. (123) yields confirming the scaling ansatz from Eq. (124). The "neighbor" dipole amplitude is defined only for x 10 > x 21 , that is, ζ > ζ . Repeating the steps from [33], as detailed in Appendix A, one arrives at the solution of Eqs. (125) in the integral form, The leading high-energy asymptotics of F (ζ) is given by the saddle points at ω = ±i √ 3 in the exponent. This is illustrated in Fig. 11, which also shows the entire complex-plane structure of the integrand of Eq. (126a) (cf. [122]). In distorting the original integration contour (the straight vertical line on the right of Fig. 11) to the steepest descent path we pick up the pole at ω = 1. However, the contribution of this pole falls of exponentially with ζ, that is, it scales as ∼ e −2 ζ , and can be safely discarded for ζ 1. We thus see that the integral in Eq. (126a) is indeed dominated by the contribution of the steepest descent contour. (Note that this situation is the exact opposite of the helicity distribution in [33], where the contribution of the right-most pole dominated over the steepest descent path.) Distorting the integration contour to run along the steepest descent path, and integrating over the regions near the saddle points at ω s.p. = ±i √ 3 by expanding , respectively, with some small real parameter ρ (which is then integrated from −∞ to ∞), one arrives at such that F (s 10 , η) = F (ζ) ≈ 3 1/4 8 √ π ζ 3/2 sin 2 √ 3 ζ − π 4 = 3 1/4 8 √ π (η − s 10 ) 3/2 sin 2 We see that the dipole amplitude F (s 10 , η) is an oscillating function of its arguments: such oscillations are interesting and reminiscent of the oscillations found in [40] for quark helicity in the large-N c &N f limit (with N f the number of quark flavors).

Re[ω]
Im[ω] We cross-checked the result (127) by solving Eqs. (123) numerically and found a very good agreement between the analytic and numerical solutions for ζ = η − s 10 3. Furthermore, to test our assumption that Eqs. (118) and Eqs. (121) give the same small-x asymptotics for the Sivers function, we solved Eqs. (118) numerically, obtaining the solution which exhibited the oscillating behavior with the same period and similarly decreasing-with-ζ amplitude of the oscillations as Eq. (127).
Employing the result from Eq. (128) in Eqs. (120) and (101) we see that the z-integral in the latter is dominated by its lower limit, which gives a contribution independent of the center-of-mass energy squared s, and, consequently, of x. Therefore, we conclude that the sub-eikonal small-x asymptotics of the quark Sivers function at large N c and in DLA is given by a constant, The approach to the constant asymptotics of Eq. (129) should be oscillatory with decreasing amplitude of such oscillations, due to the form of the amplitude in Eq. (128). This way, in principle, some residual effects of the oscillations from Eq. (128) may be observable experimentally.
Let us point out that the result (129) is, in a way, similar to the case of the odderon: while, unlike the odderon case, the (DLA) evolution at the sub-eikonal order does significantly affect the dipole amplitude F (x 2 10 , z), the x-dependence of the Sivers function is almost unaffected by the evolution, just like it was for the eikonal odderon contribution. The same conclusion (129) can be obtained by using the initial amplitude (119) in Eqs. (120) and (101).

D. Small-x Asymptotics of the Quark Sivers Function: a Summary
We conclude this Section by summarizing the results of our calculations. The quark Sivers function at small x receives contributions at the eikonal and sub-eikonal order. The eikonal contribution is coming from the spindependent odderon and is given by f ⊥ q 1 T ∼ 1/x with an almost non-perturbative accuracy (see Eqs. (54) and (66)), in agreement with [3,4]. The sub-eikonal contribution to quark Sivers function is calculated in this work for the first time. At large N c and in DLA it is given by Eq. (129). We, therefore, conclude that, at small-x, one can describe the small-x asymptotics of the quark Sivers function as with some functions C O (k 2 T , x) and C 1 (k 2 T ), which can be obtained from the above results. The function C O (k 2 T , x) also depends on x, but in a much slower way than the powers explicitly shown in Eq. (130) (see, e.g., Eq. (26) in [65]). The ellipsis in Eq. (130) denote the order-x sub-sub-eikonal corrections along with the powers of ln(1/x)-suppressed pre-asymptotic corrections to the saddle-point asymptotics shown in Eq. (130).
The situation we have found is qualitatively similar to what was suggested for the h 1 structure function in [123]: a sum of the odderon and DLA contributions. (Note, however, that for a related quantity, the valence quark transversity TMD, only the DLA contribution was found in [2].)

IV. CONCLUSIONS AND OUTLOOK
In this paper we have accomplished several results. In Eq. (42) we have constructed the full sub-sub-eikonal polarized Wilson line/quark S-matrix operator which can also be used to obtain the small-x asymptotics of a number of quark TMDs which have not been studied at small x yet. We employed this operator to study the Sivers function f ⊥ 1 T up to the sub-eikonal order. In the future, one TMD that can be studied with the help of this new operator is the Worm-Gear function g ⊥ 1T coupling the proton's transverse spin to the quark helicity: such coupling is given by the δ χ,−χ terms which can be extracted from Eq. (42). One can also apply the same operator treatment to study the small-x asymptotics of other leading-twist quark TMDs, such as the Boer-Mulders (h ⊥ 1 ) function, by employing the quark polarization-independent ∼ δ χ,χ sub-eikonal corrections from the quark S-matrix operator constructed here. Similarly one can study Pretzelosity (h ⊥ 1T ) which couples transversely polarized quarks (∼ χ δ χ,−χ ) to the transverse spin of the proton, with both spins orthogonal to each other. One can also study the other Worm-Gear function h ⊥ 1L coupling the transversely-polarized quark (∼ χ δ χ,χ ) to the longitudinal spin of the proton. Then, together with the known results for the unpolarized quark TMD at small x [23,94,99,100], quark helicity [1,24,33] and transversity [2] one would have obtained the small-x asymptotics for all the leading-twist quark TMDs. Leading-twist gluon TMDs can also be analyzed in a similar way. This would allow one to make predictions for the data of future small-x experiments studying the proton's spin structure, such as those at the upcoming EIC.
The standard Collins-Soper-Sterman (CSS) [124,125] equations usually applied to TMDs evolve them in Q 2 and not in x. As a consequence, the CSS evolution cannot predict the x-dependence of TMDs, particularly at small x. Constructing small-x evolution of TMDs, which is able to predict the x dependence of TMDs, like what was done in this paper for the quark Sivers function, is thus vital for making predictions for the future TMD measurements at the EIC and for completing our understanding of TMDs in the important small-x region of phase space probed in high energy collisions. In addition, understanding the small-x asymptotics of TMDs would provide a unique angle on the proton momentum and spin structure in the regime dominated by sea quarks and gluons.
To illustrate our method we have constructed the small-x quark Sivers function using the operator formalism introduced in [1] to study the quark helicity TMD and used again in [2] to construct the valence quark transversity TMD. We have reproduced the conclusion of [3] that the spin-dependent odderon dominates in the small-x asymptotics of the quark Sivers function. In addition, we have found the sub-eikonal correction to this odderon-dominated asymptotics. Perhaps naturally, the sub-eikonal contribution to the Sivers function comes from the gauge-covariant operator representing the sub-eikonal phase (21), since existence of a phase is essential for the Sivers function. The two terms, eikonal and sub-eikonal, are summarized in Eq. (130). For the STSA A N observable, if we conjecture that A N ∼ x f ⊥ 1 T as far as the x-dependence is concerned, our prediction (130) implies A N ∼ C O + x C 1 . Since the data [103][104][105][106] appears to be closer to A N ∼ x scaling than to A N ∼ const, we see that the sub-eikonal correction may turn out to be more important for the description of the existing data on STSA. One could speculate that the sub-eikonal correction somehow dominates over the spin-dependent odderon contribution in the experimentallyprobed x-region, possibly just numerically if C 1 C O for some (probably non-perturbative 4 ) reason. Note that the lowest-order perturbative Sivers function scales as f ⊥ 1 T ∼ x at small x [79], leading to A N ∼ x 2 , which seems to disagree with the data [103][104][105][106], though a more detailed analysis of the x-dependence of A N in the data is needed to draw firm conclusions. The sub-eikonal correction we found, once better quantified, may also provide a background for the future spin-dependent odderon searches. These results give an exciting possible new direction for such future experimental studies, particularly in light of the recent announcement by D0 and TOTEM collaborations of odderon detection through the asymmetry between pp and pp collisions [62]. Future experiments such as those to be conducted at the EIC will be able to probe transverse spin asymmetries at small-x and potentially observe both the spin-dependent odderon contribution and the sub-eikonal correction derived in this work.