Calculating TMDs of a Large Nucleus: Quasi-Classical Approximation and Quantum Evolution

We set up a formalism for calculating transverse-momentum-dependent parton distribution functions (TMDs) using the tools of saturation physics. By generalizing the quasi-classical Glauber-Gribov-Mueller/McLerran-Venugopalan approximation to allow for the possibility of spin-orbit coupling, we show how any TMD can be calculated in the saturation framework. This can also be applied to the TMDs of a proton by modeling it as a large"nucleus."To illustrate our technique, we calculate the quark TMDs of an unpolarized nucleus at large-x: the unpolarized quark distribution and the quark Boer-Mulders distribution. We observe that spin-orbit coupling leads to mixing between different TMDs of the nucleus and of the nucleons. We then consider the evolution of TMDs: at large-x, in the double-logarithmic approximation, we obtain the Sudakov form factor. At small-x the evolution of unpolarized-target quark TMDs is governed by BK/JIMWLK evolution, while the small-x evolution of polarized-target quark TMDs appears to be dominated by the QCD Reggeon.


I. INTRODUCTION
Over the past decade quark and gluon transverse momentum-dependent parton distribution functions (TMDs) [1,2] have become an integral component of our understanding of the momentum-space structure of the nucleon. At the same time the main principles for calculating TMDs in the perturbative QCD framework have remained essentially the same over the years: one parameterizes the initial conditions at some initial virtuality Q 2 = Q 2 0 and then applies the Collins-Soper-Sterman (CSS) evolution equation [1] to find the TMDs at all Q 2 . The initial conditions are non-perturbative, and have to be constructed using models of the non-perturbative QCD dynamics (see [3] and references therein). Since very little is known about non-perturbative effects in QCD, often one uses the same form for the parameterization of the initial conditions for several different TMDs, frequently assuming a Gaussian dependence of the TMDs on the parton transverse momentum k T [4]. It is clearly desirable to have a better control of our qualitative and quantitative understanding of TMDs.
To this end one can employ the recent progress in our understanding of small-x physics and parton saturation [5][6][7][8][9][10][11] to put constraints on the TMDs and even calculate them in the high-energy limit. When calculating TMDs one often works either in the s ∼ Q 2 ≫ k 2 T (large-x) or in the s ≫ Q 2 ≫ k 2 T (small-x) regimes, where s is the center-of-mass energy and x = Q 2 /(s + Q 2 ) if one neglects the proton mass. In either case the energy s is large, and the techniques of high-energy QCD should apply. The degrees of freedom in saturation physics are infinite Wilson lines along (almost) light-like paths. The definition of the TMDs involves light-cone Wilson lines as well [12,13], although the integration paths are semi-infinite, forming the so-called "light-cone staple". We can see that there are both similarities and differences between saturation physics and the physics of TMDs. The interface of these two sub-fields of quantum chromodynamics (QCD) has been explored in [14][15][16][17][18][19][20][21][22][23][24]. In the past some success has been achieved in applying saturation physics to study the Sivers function [25,26] both in semi-inclusive deep inelastic scattering (SIDIS) and in the Drell-Yan process (DY). In [27] the Sivers function was constructed by generalizing the quasi-classical Glauber-Gribov-Mueller (GGM) [28] /McLerran-Venugopalan (MV) [29][30][31] approximation of a heavy nucleus with atomic number A ≫ 1. The presence of the atomic number generates a resummation parameter α 2 s A 1/3 [32,33] allowing a systematic resummation of multiple rescatterings, which are essential for the Sivers function. This picture can also be applied to the proton if one models it as a large "nucleus." This largenucleus approximation is known to work well in describing the data from deep inelastic scattering (DIS) experiments on a proton at low-x; it is therefore possible that it would give a reasonable description for proton TMDs as well. The result of [27] was an explicit form of the Sivers function in the s ∼ Q 2 ≫ k 2 T regime, which was different from a simple Gaussian in k T and which can be used as the initial condition for its Collins-Soper-Sterman (CSS) evolution [1]. Another important result of [27] was the realization that the Sivers function can be produced via two different channels: one is the standard "lensing" mechanism of [34][35][36] with additional momentum broadening due to multiple rescatterings in the nucleus, while the other mechanism was due to the orbital angular momentum (OAM) of the nucleus combined with multiple rescatterings. This latter channel has not been reported before [27], and it dominates mainly in the regime where multiple rescatterings are important, or, more precisely, for k T Q s / √ α s , where Q s is the saturation scale of the nucleus. It appears that applications of saturation physics to the calculation of TMDs may lead to qualitatively new channels of generating the relevant observables.
The aim of this work is twofold. First of all we want to generalize the approach of [27] to the calculation of any TMD in the quasi-classical approximation and for s ∼ Q 2 ≫ k 2 T . This is accomplished in Sec. II for the case of an unpolarized nucleus; the generalization to the polarized case is straightforward and is left for future work. The quasi-classical TMD calculation is accomplished using the factorization given in Eq. (18) (and, in more detail in Eq. (39)), which is constructed by analogy with the particle production cross section: the TMD is a convolution of the classical Wigner distribution W describing the nucleons in the nucleus with the TMD of one of the nucleons (φ N ) and the Wilson lines S of the "staple" evaluated in the GGM/MV approximation. The functional form of the nucleon Wigner distribution in the nucleus can be constrained using the symmetries of the nuclear ground state. Considering for simplicity an unpolarized nucleus we arrive at the parameterization of Eq. (55). With the help of this unpolarized nuclear Wigner distribution we construct the unpolarized quark distribution in Eq. (61) and the Boer-Mulders distribution [37] in Eq. (78).
Note that the nuclear effects on the TMDs are not limited to the k T -broadening as is commonly assumed in the literature [38]. Other TMDs for the (transversely or longitudinally) polarized nucleus and for the unpolarized nucleus could be constructed by analogy, but we limit the examples to the two TMDs mentioned above along with the Sivers distribution constructed in [27].
An interesting consequence of the calculation is a mixing between different TMDs of the nucleus and of the nucleons. For instance, as was already observed in [27], the nuclear Sivers The second aim of this work is to understand TMD evolution using the methods of saturation physics. In Sec. III we consider TMD evolution at large and small values of Bjorken-x. The large-x regime, corresponding to s ∼ Q 2 ≫ k 2 T , is considered in Sec. III.1. There we show that the TMD evolution can be thought of as the evolution of the light-cone Wilson lines in the "staple" (if one works in the light-cone gauge of the probe), as shown in Fig. 6. In the large-x regime the Wilson lines combine to form an operator corresponding to the propagation of a (quark or gluon) dipole from the nucleus to positive light-cone infinity (SIDIS) or from negative infinity to the nucleus (DY). (A fundamental semi-infinite dipole is used for quark TMDs, while a semi-infinite adjoint gluon dipole should be used for the gluon TMDs.) In the double-logarithmic approximation resumming powers of α s ln 2 (Q 2 /k 2 T ) we recover the standard Sudakov form factor [39] which one also obtains from the CSS evolution.
In the low-x regime (s ≫ Q 2 ≫ k 2 T ) the diagrams giving the leading contribution to the TMDs are slightly different from the large-x graphs, as illustrated e.g. in Fig. 8. For unpolarized-target quark TMDs in the leading single-logarithmic approximation, which resums powers of α s ln(1/x), we show in Sec. III.2 that TMD evolution is driven by the evolution of a fundamental dipole scattering amplitude on a nucleus, that is by the evolution of the dipole made out of infinite Wilson lines. The evolution of this object is well-known: in the large-N c limit it is given by the Balitsky-Kovchegov evolution equation [40][41][42][43], while beyond that limit it can be found using the Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) evolution equation [44][45][46][47]. For the polarized-target TMDs the situation is more subtle, as the evolution has to carry the polarization information. We argue that polarized quark TMDs at small-x are governed by the nonlinear version [48] of the QCD Reggeon evolution equations [49][50][51][52][53][54][55] (see Fig. 13 and Eq. (94)), which employs the dipole amplitude obtained from the BK/JIMWLK equations.
We have thus constructed quasi-classical expressions for TMDs at large-x in the GGM/MV approximation, and constructed evolution of TMDs at both large-and small-x. The results of this work can be used for constructing a global TMD fit.

II. QUASI-CLASSICAL INITIAL CONDITIONS
In this Section we expand and generalize the approach developed in [27]. While the quasi-classical method detailed below was used in [27] to calculate the Sivers function of a large nucleus, here we will apply it to calculating TMDs of an unpolarized nucleus, such as the unpolarized quark distribution and the Boer-Mulders function. The transverse-momentum-dependent (TMD) quark correlation function in a hadronic state |h(P, S) with momentum P and spin S is defined by where α, β are Dirac indices, the separation vector between the quark fields is r µ = (0 + , r − , r), and the "staple-shaped" gauge link U[0, r] extends to future/past light-cone infinity (±∞ − ) depending on the process under consideration. The abbreviated notation for the integration measure is d 2− r = d 2 r ⊥ dr − . We work in light-front coordinates where g +− = 1 and g +− = 2 are two common choices of the light-front metric. We will use a frame in which the nucleus or the proton have a momentum predominantly in the Projecting the TMD correlator (1) onto various Dirac matrices gives the distribution of quarks with zero, longitudinal, or transverse polarizations, and the independent spinspin and spin-orbit correlations are parameterized in terms of boost-invariant TMD parton distribution functions (the "TMDs"). At leading power in the large momentum P + , the TMD decomposition can be written as (see, e.g. [37,56,57]) where (S, S L ) denote the transverse and longitudinal components of the spin vector, m is the hadron mass, and k × S ≡ k x S y − k y S x with (x, y) the coordinates in the plane transverse to the beam. The eight leading-twist quark TMDs shown here are functions of x and k T = |k|.
At lowest order in α s such that the gauge link U[0, r] can be neglected, the correlator (1) has a simple interpretation as the expectation value of the quark density in the hadronic state |h(P, S) : where Ω ≡ g +− p + d 2− b is a boost-invariant volume factor which normalizes the plane-wave states |h(P, S) . The gauge link U[0, r] is needed to make the definition gauge invariant, and it reflects the distortion of the quark distribution in the presence of initial-or final-state gauge fields [58,59].
By tracing (5) with one of Γ ∈ {γ + , γ + γ 5 , γ 5 γ + γ j ⊥ }, we project out the distributions corresponding to unpolarized quarks (U), longitudinally-polarized quarks (L), and transverselypolarized quarks in directionĵ (labeled T j ), respectively. This can be explicitly verified by evaluating the associated spinor products in (5) using the spinor conventions of [60]: These projections select out the corresponding set of TMDs from the decomposition (3), so that f 1 represents the azimuthally symmetric distribution of unpolarized quarks in an unpolarized hadron, the Sivers function f ⊥ 1T represents the azimuthally asymmetric distribution of unpolarized quarks in a transversely polarized hadron, and so on.

II.1.2. Quasi-Classical Factorization in a Heavy Nucleus
As derived in [27], a heavy nucleus with A ≫ 1 nucleons in the quasi-classical approximation α 2 s A 1/3 ∼ O (1) admits a decomposition of its TMD quark correlator Φ A in terms of the TMD quark correlator φ N of its nucleons, Wilson lines, and the nuclear Wigner distribution. The heavy nucleus in this parametric limit justifies two powerful simplifications of the definition (1) as follows. First, the large number of nucleons justifies a mean-field description of the nucleus |A(P, S) in terms of the light-front wave functions of single-nucleon states |N(p, σ) . Second, because of the large number of nucleons, the initial-/ final-state rescattering described by the gauge link U[0, r] is more likely to occur on the many (A − 1) spectator nucleons, rather than on the same nucleon from which the quark distribution is taken. Schematically, we can write this as The resummation α 2 s A 1/3 ∼ O (1) [33] is a systematic way of calculating the multiple rescattering in the weak-coupling saturation framework, and the accuracy to leading order in A −1/3 means that the sensitivity of the gauge link u[0, r] to the active nucleon is limited to O (α s ).
In this way, the TMD correlator (1) is factorized into a convolution of 3 factors: a wave function piece describing the distribution of nucleons in the nucleus, the TMD correlator of the nucleon itself, and a piece describing the multiple rescattering of the gauge link on the spectator nucleons. This essential picture should apply not only to a heavy nucleus, but to any system -such as a proton at high energies -in which the density of color charges becomes sufficiently large. For semi-inclusive deep inelastic scattering in Bjorken kinematics (see Fig. 1), the quasiclassical factorization formula is [27] Φ A αβ (x, k; P, S) = A wherex ≡ x P + p + is the quark momentum fraction with respect to the active nucleon, d 2+ p = d 2 p ⊥ dp + , and b µ ≡ (0 + , b − , 1 2 (x + y)) is the position of the struck nucleon. The Wigner distribution of nucleons |N(p, σ) inside the nucleus |A(P, S) is is the average momentum of the struck nucleon in the amplitude and complex-conjugate amplitude and σ is its spin. The Wigner distribution is normalized such The semi-infinite dipole scattering amplitude D xy [∞ − , b − ] describes the final-state rescattering on the fraction of nucleons at depths greater than b − and is given in covariant gauge by where is a Wilson line in the fundamental representation and T a is the fundamental generator of SU(N c ). The angle brackets in Eq. (14) denote averaging in the nuclear wave function.
The origin of the nuclear TMD decomposition (11) is in a similar decomposition for the quark production cross section in the SIDIS case [27], illustrated in Fig  The quark production cross section is [27] dσ γ * +A→q+X Note that the light-cone Wilson lines describing the produced quark in the amplitude and in the complex conjugate amplitude (in the above-mentioned gauges) have to be at different transverse positions x ⊥ and y ⊥ since the transverse momentum of the quark k ⊥ is fixed (see e.g. [61]). Together these Wilson lines contribute to the "staple" U in the TMD definition Eq. (11) can be simplified further by replacing the dipole scattering amplitude D xy with its symmetric part S xy ≡ 1 2 (D xy + D yx ) since the anti-symmetric (odderon) part is suppressed by A −1/3 [19]. If an analytic expression is desired, this symmetric part of the dipole scattering amplitude can be evaluated in the quasi-classical GGM/MV multiple scattering approximation to give [28] Here ) is the fraction of nucleons which participate in the final-state rescattering. The logarithm with infrared (IR) regulator Λ, which is often neglected, is only important for recovering the perturbative large-k T (small-|x − y| T ) asymptotics. Alternatively, the Wilson line operators can be evaluated on a configurationby-configuration basis if desired using Monte Carlo methods along the lines of [62,63]. By replacing D xy with S xy , we can rewrite the quasi-classical factorization formula (11) as where we have changed variables to b ≡ x+y 2 and r ≡ x − y and shifted the integration variable k ′ → k ′ +xp.

II.1.3. Lorentz-Covariant Spin Decompositions
The quasi-classical factorization formula (18)  Let us first write a decomposition of the Wigner distribution (12) for nucleons with an arbitrary spin state |S in terms of the light-cone helicity basis |± . For simplicity, we will restrict ourselves to an unpolarized nucleus, but which may have polarized nucleons, which we write in the compact form where δp ≡ p − p ′ ,p ≡ 1 2 (p + p ′ ), and b µ ≡ (0 + , b − , b). As described in [64], an arbitrary spin state can be decomposed in terms of the light-cone helicity basis as |S ≡ cos θ 2 |+ + sin θ 2 e iφ |− ; (20) in the rest frame (R.F.) of the nucleon, the eigen-axis of spin projections for this state is given by the three-dimensional spin vector with λ 2 + S 2 T = 1. The Lorentz-covariant spin vector S µ is obtained by defining (21) as the spatial part in the rest frame and boosting it to a frame in which the nucleon momentum is Note that, by construction, p µ S µ (p) = 0. We have also made use of the fact that the nucleons do not interact with each other in the quasi-classical approximation at hand [66][67][68].
By forming the outer product of (20) and using (22) , we obtain |S S| R.F.  [65,66], obtained by boosting (0, S R.F. ) µ out of the rest frame with a single rotationless boost. This can also be regarded as defining S µ ≡ W µ /m with W µ the Pauli-Lubanski four-vector, and obtaining S R.F. via Eq. (22). Alternatively, one may use light-front boosts, obtaining an expression for S µ in terms of S ′ R.F. , which is the spin in a rest frame rotated with respect to the rest frame used in the second line of (22) (see [65] and the footnote 35 on page 213 of [66]).
Ultimately, for the situations we consider (nonrelativistic nucleon motion and an ultrarelativistic boost along the +z-axis), both definitions of the spin three-vector coincide. We thank Cédric Lorcé for helping us clarify these subtleties.
where we have introduced the Pauli matrices and unit matrix. If we define a Lorentzcovariant set of Pauli matrices analogous to (22) where again p µσ µ (p) = 0 by construction, then we have S µ (p)σ µ (p) = − S R.F. · σ in any frame, and Inserting this into (19) we obtain with where we have also introduced the Hermitean (2 × 2) matrix The projections W unp and −S µŴ µ pol select out the Wigner distributions of nucleons which have zero polarization, or polarization S µ , respectively. And although they were derived here in terms of the light-cone helicity basis |± , they are invariant under a change of basis; any unitary rotation of the spin states transforms the matrix W λλ ′ , but it also rotates the Pauli matrices such that the traces (27) remain invariant. We can also write a decomposition of the matrix W λλ ′ in terms of the complete basis {1, σ} to invert Eqs. (26), (27): From (24) we obtain the trace identity which allows us to recover (27) from (29) by tracing the latter withσ ν or 1 and using the fact thatp µŴ Note that the normalization of the matrix Wigner distribution W λλ ′ (p, b) is such that In a similar way, the quark correlator (1) of the nucleon can be expanded in a basis of spin states, and the polarized and unpolarized parts can be projected using the Pauli matrices: where and the matrix form is Note that S µ andσ µ above depend on the (averaged) momentum of the nucleonp, which In the derivation of (18), a sum over the spins of the intermediate nucleons occurs independently in the amplitude and complex-conjugate amplitude (see Fig. 1). Thus the spins entering the quasi-classical factorization formula are necessarily non-diagonal, which can be compactly expressed in terms of the matrices W λλ ′ and φ λλ ′ : where we have restricted the discussion to an unpolarized nucleus. From (29), (30), and (33), the trace of the two matrices gives which sums over all four independent spin states of the intermediate nucleons. Thus we recover a form similar to (18) in which the spins of the nucleons are diagonal (i.e., the same in the amplitude and complex-conjugate amplitude), as illustrated in Fig. 3. The sum now effectively runs over the four polarizations: U = unpolarized, L = longitudinally-polarized, and T j = transversely polarized in theĵ direction.
This formulation makes the projections defined by (26), (33) independent of the choice of spin basis and expresses the spin dependence in a manifestly Lorentz-covariant form.

II.2. Parameterization of the Wigner Distribution
The quasi-classical factorization formula (39) describes how the properties of the Wigner distribution and the multiple rescattering build up the TMDs of the nucleus from the various TMDs of the nucleons. One powerful feature of (39) is that the Wigner distribution is constructed purely from the light-front wave functions (see, e.g. (28)), without contamination from the scattering dynamics embodied in the gauge link. It is therefore invariant under several symmetries which the TMDs themselves are not, such as P T . We would like to use all the applicable symmetries to parameterize the types of structures which can occur in the Wigner distribution.
The Wigner distribution also does not know about the direction of the collision axis; that information enters operatorially through the gauge link U[0, r] or diagrammatically through the interaction with the virtual photon. We may choose to describe the functional dependence of the Wigner distribution using variables which are appropriate for a collision along the z-axis (e.g. p + , p), but the distribution itself should not treat z as a "special" direction. Naively, one would expect that, when viewed from the rest frame of the nucleus, the Wigner distribution W (p, b) should possess three-dimensional rotation symmetry due to the lack of a preferred collision axis. This is in contrast to the factors φ(x, k) and S(r T , b T ) containing the gauge link, which have at most two-dimensional rotational invariance due to the special role of the longitudinal direction.
There are some subtleties about the application of full rotational symmetry to the Wigner distribution constructed from light-front wave functions; these are addressed most clearly using the "covariant light-front formalism" of [69]. We discuss these issues in Appendix A.
The conclusion is that, if the nucleons move nonrelativistically in the nucleus, then their wave functions -and hence the Wigner distribution in (A19) and (A20) -do possess the desired rotational invariance in the nuclear rest frame. The boost-friendly variables α, b − are given in terms of the rotation-friendly variables p, b through (A21). We can now take advantage of this symmetry to constrain the functional form of the Wigner distribution.

II.2.1. Wigner Distribution of an Unpolarized Nucleus
For an unpolarized nucleus with a non-relativistic distribution of nucleons, the Wigner distribution (A20) possesses manifest rotational invariance in the rest frame and depends on the nucleon spin only linearly. Moreover, because the Wigner distribution is built purely from the nucleon wave functions, it possesses the unbroken discrete symmetries of (ordinary) parity and time reversal. The P and T symmetries of the wave functions translate into the Wigner distribution in the expected way: We can therefore constrain the form of the Wigner distribution using these discrete symmetries and rotational invariance.
For the unpolarized distribution, requiring rotational invariance yields but the quantity p · b is odd under time reversal, so W unp can only depend on it through even powers: Now let us change variables to the boost-invariant quantities using (A21): where α = p + /P + is the longitudinal momentum fraction of the nucleon in the nucleus. In terms of these quantities, the unpolarized distribution depends on where we have separated out the dependence on the longitudinal and transverse variables.
There is one interesting term which appears to mix the longitudinal and transverse variables, arising from the cross term (p b )(p·b) of ( p· b) 2 . In principle such a term is not prohibited by the symmetries of Wigner distribution; however, looking back at the quasi-classical factorization formula (39), we see that this term will not contribute to the cross-section. The reason is that the only other dependence on b arises from the dipole scattering amplitude (17). This factor does not possess the same three-dimensional rotation invariance as the Wigner distribution (A20), because the multiple scattering depends on the depth b − of the struck nucleon. It does, however, retain a residual two-dimensional rotation invariance in the transverse plane, depending only on b − , b T . Thus a term like would integrate out to zero in (39) because of antisymmetry under b → −b. Therefore the terms of the Wigner distribution which survive the d 2 b integral can only depend on the square of this cross-term, , and the dependence on these squared factors are already taken into account. Thus we can replace without missing out on any terms which would survive the d 2 b integral in (39). We use an arrow here rather than an equality to emphasize that we are now dropping contributions to the Wigner distribution which are permitted by symmetry, but would not survive to contribute to the nuclear TMD.
In fact, we can carry this argument one step farther: after the integration is performed, all sensitivity to the direction of the two-vector b drops out. Thus we can in the Wigner distribution now, without loss of generality in the types of terms which can contribute to the quasi-classical factorization formula (39). Doing so further reduces the number of terms which can contribute, since Thus, from all of these simplifications, we can write the relevant part of the unpolarized Wigner distribution as By considering the constraints due to parity, time reversal, and three-dimensional rotation invariance of the Wigner distribution, together with the two-dimensional rotational invariance of the multiple scattering factor (gauge link) to which it couples, we have reduced the unpolarized distribution down to a form which is manifestly even under p → −p, b → −b, These symmetries strongly constrain the interplay between the factors W (p, b), φ(x, k ′ ), S(r, b) entering the quasi-classical factorization formula (39).
Similarly, we can impose rotational invariance on the polarized part of the Wigner distribution: The spin-dependent factors ( S · b) and ( S · p) are odd under parity, while all of the arguments p 2 , b 2 , p· b are P -even; therefore these types of spin dependence cannot enter into the polarized Wigner distribution. The only spin dependence which is permitted by parity is the factor ( S · L) describing spin-orbit coupling, where L = b × p is the orbital angular momentum of the nucleon. Since ( S · L) is even under both parity and time reversal, the T -odd quantity ( p · b) can only occur in even powers, just like in the unpolarized distribution: Now let us again change back to the boost-invariant variables (A21). In addition to the expressions in (43), we must account for the spin-dependent factor to write As before, we now replace b i ⊥ b j ⊥ → 1 2 b 2 T δ ij to keep only the terms which can contribute to the nuclear TMD. In addition to reducing (p · b) 2 T → 1 2 p 2 T b 2 T , this also simplifies the spin dependence for the factors which couple S to b: We see that the spin dependence of the form λ(b × p) has dropped out completely, and the spin dependence of the form (α − 1 A )(S × b) has reduced down to the form (P + b − )(p × S). Thus, without loss of generality, we can identify the structure (P + b − )(p × S) as the only spin dependence which survives the d 2 b integral to contribute to the nuclear TMDs: Again, consideration of the relevant symmetries has reduced the form of the polarized Wigner distribution significantly. The structure that survives contains a prefactor linear in the transverse spin which is odd under b − → −b − and p → −p, times a function which is designation "OAM" reflects the fact that this structure originated from the presence of L · S coupling in the rest frame.
Altogether, (48) and (54) allow us to parameterize the most general structure of the Wigner distribution which is consistent with parity, time reversal, and 2D and 3D rotation symmetries. The structure which can contribute to the quark TMD of an unpolarized nucleus is or, in the notation of (26), Eq. (55) is illustrated in Fig. 4. The parameterization given here was derived for the distribution of nucleons with polarization S N = (S, λ) moving nonrelativistically in an unpolarized nucleus, for which the only nontrivial spin dependence entered as ( L · S N ). It is straightforward to extend this parameterization to describe a polarized nucleus as well. In addition to the equivalent spin-orbit coupling ( L · S A ) which was the source of the "orbital angular momentum channel" to generate the nuclear Sivers function f ⊥A 1T in [27], there are spin-spin couplings like ( S N · S A ) and spin-orbit-spin couplings like ( S N · p)( p · S A ). We will leave this generalization to the full structure of the Wigner distribution of a polarized nucleus for future work.

II.3. Quasi-Classical TMDs of an Unpolarized Nucleus
The TMD quark correlator Φ A unp of an unpolarized nucleus from (36) contains just two leading-twist TMDs: the unpolarized quark distribution f A 1 and the Boer-Mulders distribution h ⊥A 1 : Because f A 1 corresponds to unpolarized quarks and h ⊥A 1 corresponds to transversely-polarized quarks, these two TMDs can be projected from the correlator as in (9) by tracing the Dirac structure with Γ = γ + or Γ = γ 5 γ + γ j ⊥ , respectively: These nuclear TMDs are related to the TMDs of the nucleons through the quasi-classical factorization formula with the Wigner distribution of nucleons parameterized by (56). Since the nucleons can be polarized, the active quarks can come from one of several nucleonic TMDs depending on the polarization, as in (36). Combining the Wigner distribution for a given nucleon polarization with the associated TMDs, we obtain where α = p + /P + andx = x/α. The projection [Γ] of the nuclear TMD quark distribution (58)  The projection Γ = γ + in (58) selects out the unpolarized quark distribution f A 1 in the nucleus and the unpolarized quark distribution f N 1 and quark Sivers function f ⊥N 1T of the nucleon in (60) (see also (9)) . Using this in (59) gives which shows that ( L· S N ) spin-orbit coupling given by W OAM can result in a mixing between the unpolarized quark distribution f A 1 of the nucleus and the quark Sivers function f ⊥N 1T of the nucleons. This is the mirror effect to the OAM channel that was found in [27], in which ( L · S A ) coupling in a polarized nucleus gave rise to a mixing between the Sivers function f ⊥A 1T of the nucleus and the unpolarized f N 1 distribution of the nucleons. The unpolarized quark distribution f 1 is a P T -even function and therefore process independent (see [70] and others); that is, it should give the same result when calculated with a future-pointing gauge link like (17) appropriate for SIDIS as when calculated with a past-pointing gauge link like appropriate for the Drell-Yan process (DY). In the quasi-classical approximation, the difference between the future-pointing gauge link of (17) describing the fraction of nucleons which participate in the scattering. Here T (b T ) is the number density of nucleons per unit transverse area at impact parameter b T . Because f A 1 is P T -even and process-independent, it should give the same result whether evaluated with the future-pointing gauge link of (17) or the past-pointing gauge link of (62). This is clear for the trivial f N 1 channel because both f A 1 and f N 1 are invariant under the P T transformation, and the longitudinal b − integral is also preserved by a simple change of variables. The Sivers function f ⊥N 1T , on the other hand, is P T -odd and changes sign between the two processes, but its associated longitudinal b − integral also changes sign. Thus both channels of (61) are overall even under P T , yielding the process-independent nuclear TMD f A 1 as required. Similar to [27], the ( L · S N ) spinorbit channel here requires at least one additional rescattering to be nonzero; in the limit r 2 T Q 2 s ≪ 1, corresponding to either large transverse momentum or low densities, the multiple scattering factor S(r T , b T ) can be neglected, and the b − integral (65) vanishes. The depth dependence on b − of the spin-orbit factor (P + b − )(p × S), coupled to the angular dependence of the Sivers function, provides a second P T -odd factor (P + b − )(p · k ′ ) which enables a P Todd nucleonic TMD like f ⊥N 1T to contribute to a P T -even nuclear TMD like f A 1 . As we will show, the same ( L · S N ) spin-orbit coupling in the nucleus also generates contributions to the P T -odd nuclear Boer-Mulders function from various T -even nucleonic TMDs.
To better understand our result, let us evaluate Eq. (61) for specific models of the Wigner distribution. First let us consider the case when the internal motion of the nucleons in the nucleus is negligible, that is the nucleons are static with respect to the nuclear center-of-mass (the "standard" GGM/MV model in saturation physics). The corresponding unpolarized Wigner distribution is (cf. Eq. (65) in [27]) with the constant nucleon number density (per (The normalization for W unp is defined in Eq. (32).) Since the nucleons are static, W OAM = 0 and the second term on the right-hand side of Eq. (61) vanishes. Finally we take for the unpolarized quark TMD of the nucleon the lowest-order perturbative expression for a quark target [48,57] which we for simplicity consider in the k T ≫ m N limit. As usual C F = (N 2 c − 1)/(2N c ) is the fundamental Casimir operator of SU(N c ). Using (68), (66) and (17) in Eq. (61) yields The low-x (and lowx = Ax) limit of Eq. (69) is in complete agreement with the previous calculation of a similar quantity in [48] (see Eqs. (22) and (23) there keeping in mind that unintegrated quark distribution in [48] is defined per dk 2 T dx/x unit of phase space, while TMDs are conventionally defined per d 2 k T dx). The quark saturation scale in Eq. (69) is Note that the agreement between our formula (61) for the unpolarized quark distribution with the results obtained from a different method in [48] validates our approach and suggests that other observables, like the unpolarized gluon distribution, when calculated in our technique would agree with the results in the literature [18,32,61,71].
As is well-known, putting the logarithm in the exponent of Eq. (69) to one (as a negligibly slow varying function) casts the equation in a form where the integral over r T can be carried out analytically. This gives This approximation is valid only for k T not much larger than Q s .
To get a feeling for the correction to Eq. (69) resulting from the spin-orbit coupling term in Eq. (61) assume that where ξ is a parameter responsible for the strengths of spin-orbit coupling and we neglect the small change in p + due to the nonrelativistic orbital motion. Just like in [27] we approximate the nucleon Sivers function by the quark target expression, assuming again that k T ≫ m N : Employing Eqs. (71), (72) and (17) in Eq. (61) we get the following correction to the unpolarized quark TMD of a large nucleus: where ρ = 3A 4πR 3 is the nucleon number density in the rest frame and we have made use of the m N ≪ Q s condition and put ln(1/r T Λ) ∼ 1 in the exponents. Comparing Eqs. (73) and (70) we see that the spin-orbit contribution is parametrically different by a factor of ξ/α s .
Without the knowledge of the (most likely) non-perturbative parameter ξ it is impossible to say whether the factor of ξ/α s indicates suppression or enhancement. If we take ξ ∼ α 4 s by analogy to atomic spin-orbit coupling in quantum electrodynamics (QED), then ξ/α s ∼ α 3 s which is parametrically very small. However there is no good reason to justify such a choice in QCD. Further work is needed to fully quantify the role of spin-orbit coupling in the unpolarized quark TMD.

1
Similarly, the projection Γ = γ 5 γ + γ j ⊥ in (58) (9)). Using this in (59) gives It is convenient to separate the diagonal part δ jℓ of the quantity in parentheses, where is the transversity TMD and h ⊥N 1T is the pretzelosity. Then Then we can contract the free index j on both sides with −k n ⊥ ǫ nj T and solve for the nuclear Boer-Mulders function to obtain Thus we see that the same L · S N coupling which allows the quark Sivers function f ⊥N  (64). But there are also P T -even channels -the nucleonic transversity h N 1 and pretzelosity h ⊥N 1T -which mix into the P T -odd nuclear Boer-Mulders function through the role of orbital angular momentum. Just like in the case of the nucleonic Sivers function f ⊥A 1T mixing into the unpolarized quark distribution f A 1 of the nucleus, the necessary reversal of P T symmetry is provided by the depth dependence (P + b − ) in the longitudinal integral (65). And just like before, this mixing requires at least one rescattering to provide this extra P T -odd factor, so it vanishes in the limit of large k T or low charge densities.

III. EVOLUTION
Our goal in this Section is to include the quantum evolution corrections to the quasiclassical TMDs like those calculated above. We will separately consider the cases of x = O(1) (large-x) and x ≪ 1 (small-x).

III.1. Large-x Evolution
The TMDs calculated above and in [27] were found for x ∼ O(1): here we need to find the quantum evolution corrections to a generic TMD given by Eq. (39). Note that we are working in the s ∼ Q 2 ≫ k 2 ⊥ regime, where s is the center-of-mass energy squared per nucleon; hence our x in this section (and above) is order-one, but does not closely approach unity. Let us begin by concentrating on summation of leading logarithms of energy, that is on powers of α s ln(s/k 2 T ). To do this we can use the established formalism of saturation physics [40][41][42][43][44][45][46][47][72][73][74].
We work in the frame where the incoming nucleus has a large P + momentum, P µ = (P + , 1 2 g +− M 2 A P + , 0), while the incoming virtual photon has a large positive q − along with a comparable but negative q + momentum, q µ = (− 1 2 g +− Q 2 q − , q − , 0). (For definitiveness we consider the SIDIS process, though all our discussion and conclusions apply to DY as well.) We will work in the A − = 0 light-cone gauge. Such emissions are easy to sum up. Using crossing symmetry, the correction to a light-cone Wilson line at x ⊥ in the amplitude and another light-cone Wilson line at y ⊥ in the complex conjugate amplitude, as appears in the cross section, can be accounted for by calculating corrections to a pair of Wilson lines in the amplitude: one at x ⊥ another one at y ⊥ . This is illustrated in Fig. 6 and is the basis for expressions like Eq. (39), though in the quasi-classical case one resums multiple rescattering diagrams. The crossing symmetry in Fig. 6 is valid both for multiple rescatterings and for gluon emission [75]. We conclude that all we have to do is evolve a semi-infinite dipole stretching in the x − direction from 0 to +∞. (In the DY case the dipole would stretch from −∞ to 0, but the evolution equation would be the same.) In the leading-logarithm of energy approximation the evolution for this object is simply given by the (half of) virtual corrections to the standard small-x evolution of an infinite dipole [40][41][42][43][44][45][46][47][72][73][74] and reads here. All we know is that this logarithm of energy has a cutoff of the order of k T .
Eq. (79) is easy to solve [72]. Integrating over z ⊥ on the right-hand side yields with ρ some ultraviolet (UV) cutoff having dimensions of distance. Since the shortest transverse distance in the problem is of the order of 1/Q ∼ 1/s (since we assume that s ∼ Q 2 ), we make another approximation: keeping only the leading logarithms of Q 2 we replace ρ → 1/Q. In addition, remembering that we have s ∼ Q 2 we see that In the resulting double-logarithmic approximation (DLA) (in ln Q 2 ) we rewrite Eq. (80) as The solution of Eq. (82) reads with the initial conditions S xy [∞ − , b − ](Q 2 0 ) at Q 2 0 = 1 (x−y) 2 given by Eq. (17). The exponent in Eq. (83) is the well-known Sudakov form-factor [1,39], previously derived in the saturation literature in [76][77][78][79]. For gluon TMDs, where the Wilson lines are adjoint, one would have to replace C F → N c in Eq. (83). While Eq. (83) is written for a fixed coupling constant α s , running coupling corrections to it can be included following [80][81][82].
Note that while we started out this derivation by studying evolution of a semi-infinite dipole with energy, we have been working at large-x where s ∼ Q 2 and evolution with ln(s/k 2 T ) is equivalent to evolution in ln(Q 2 /k 2 T ), as follows from Eq. To summarize, the QCD-evolved large-x quark TMDs of a large unpolarized nucleus in the saturation picture can be obtained from the quark-quark correlator with the semi-infinite dipole scattering amplitude given by Eq. (83) in the DLA and 2 given by Eq. (17). For gluon TMDs the Sudakov formfactor would only be different by the C F → N c substitution.

III.2.1. Evolution of Unpolarized-Target TMDs
Now let us consider TMD evolution in the small-x regime, s ≫ Q 2 ≫ k 2 T . We begin with the unpolarized proton or nucleus TMDs, concentrating on the unpolarized quark TMD f A 1 first. As usual, the dominant contribution at small-x comes from diagrams which are order-α s suppressed compared to the channel shown in Fig. 2 above. Concentrating on the SIDIS process again, the dominant small-x contribution is pictured in Fig. 7. Comparing   Fig. 7 to the lowest-order (no multiple rescatterings) part of Fig. 2, one immediately sees that the former is order-α s suppressed: this is why the channel in Fig. 7 was not considered above. However, at small-x the channel in Fig. 7 is dominant, being enhanced by a power of 1/x compared to the channel in Fig. 2. One may ask why the diagram in Fig. 7 was not included along with the evolution corrections considered in Sec. III.1. Indeed the quark loop in Fig. 7 generates a factor of ln(Q 2 /k 2 T ), such that the suppression of Fig. 7 compared to the lowest-order part of Fig. 2 is diminished to a factor of α s ln(Q 2 /k 2 T ). However it is well-known that the quark loop cannot generate a logarithm of energy; hence the diagram in Fig. 7 is a single-logarithmic correction, which is beyond the DLA accuracy of the evolution in Eq. (83) and can be neglected at large-x.
The splitting of a virtual photon into a qq pair in Fig. 7 happens typically long before the interaction with the target nucleus, denoted by a shaded oval in Fig. 7. Moreover, in the unpolarized SIDIS case all interaction with the target is eikonal and, hence, spinindependent. No interaction with one given nucleon is a special "knockout" interaction, where spin-dependence could be transferred from the target to the probe, as in [27]. Therefore the unpolarized SIDIS cross section at small-x, and the corresponding TMDs, do not have the b − dependent effects akin to those at large-x as shown in Eq. (39) leading to TMD mixing in Eqs. (61) and (78). We also do not need to worry too much about the details of nucleon dynamics encoded in the nuclear Wigner distribution. With this in mind, summing up the diagrams in Fig. 8, where the thick vertical band represents the shock wave and we explicitly show the four contributions where the γ * → qq splitting happens either before or after the photon passes through the shock wave on either side of the cut, we use the light-cone perturbation theory [60] to write the SIDIS cross section as (see [11] for a detailed presentation of a calculation similar to this one, along with the appropriate references) where Ψ γ * →qq is the well-known light-cone wave function for the γ * → qq splitting (normalized as in [11]) given by To extract the unpolarized quark TMD we first impose the Q 2 ≫⊥ 2 condition on Eq. (85), with ⊥ denoting any transverse momentum in the problem (this is part of the s ≫ Q 2 ≫⊥ 2 condition defining our small-x regime). It is well-known that the large-Q 2 asymptotics of Eq. (85) comes from the aligned-jet configurations dominated by either z → 1 or z → 0 regions of phase space [87,88]. Since we are interested in producing a quark with k − ≈ q − we will only take the z ≈ 1 region into account. Note that the transverse photon polarizations lead to a leading power of Q 2 . The large-Q 2 limit of Eq. (85) with z ≈ 1 averaged over transverse photon polarizations is (see e.g. [89] for details of taking the large-Q 2 limit). We have put m f = 0 for simplicity.
To read off the quark TMD f A 1 we compare Eq. (87) to the lowest-order SIDIS process, modeling the target as a single quark. The dominant contribution to this process is illustrated in Fig. 9, resulting in the lowest-order SIDIS cross section where we have also made use of the small-x limit of Eq. (68) for f quark 1 .
Comparing Eq. (88) to Eq. (87) and using s ≈ Q 2 /x we read off the small-x unpolarized quark TMD for a large nucleus to be (see [90] for a similar result) Having the expression (89) for the unpolarized quark TMD we can now determine its evolution at small-x. It is driven by the small-x evolution of the dipole S-matrix. The small-x evolution of S in the leading-logarithmic approximation (LLA) resumming powers of α s ln(1/x) is well known and is given by [7,40,41,[44][45][46][47]91] x,y Y with Y = ln 1/x ≈ ln(s/Q 2 ). In Eq. (90) we have separated S xy (Y ) = Ŝ into the Wilsonline operatorŜ and the averaging in the nuclear wave function evolved up to rapidity Y (see Eq. (14)). Unfortunately Eq. (90) is not a closed integro-differential equation: the object ŜŜ on its right-hand side is different from Ŝ on the left-hand side. The equation closes in the large-N c limit, where it becomes [40][41][42][43] where we again dropped the angle brackets. The initial condition for Eq. (91) is given by (cf. Eq. (17)) Eqs. (89), (90) and (91) give us the expression for and the small-x evolution of the unpolarized quark TMD f A 1 .
x ⊥ y ⊥ The small-x evolution of unpolarized gluon TMD is constructed in a similar manner, with some important differences. Instead of DIS with a photon, which only couples to quarks, consider "DIS" with a scalar current j = −(1/4)F a µν F a µν [28] which couples to gluons. A gluon dipole diagram contributing to the evolution of the unpolarized gluon TMD is shown in Fig. 10. Note that the gluon loop brings in a logarithm of 1/x: hence the diagram is only a small-x evolution correction to the gluon analog of the lowest-order graph in Fig. 9. To include the small-x evolution into gluon TMDs, one has to start with such a lowest-order diagram and first include the Wilson line staple, obtaining the Weizsäcker-Williams (WW) gluon distribution in the classical approximation [61,71,[92][93][94]. Rewriting the semi-infinite adjoint Wilson lines in terms of (derivatives of) the infinite Wilson lines, one can apply the JIMWLK evolution obtaining the resulting evolution equation governing the small-x asymptotics of the WW distribution, as was done in [17,95].
One may worry whether the diagram in Fig. 10 (and other such corrections) needs to be included at large-x, since the gluon loop may give both logarithms of Q 2 and s. The resolution of this question is in the fact that the energy logarithm coming from the gluon loop in Fig. 10 is ln(s/Q 2 ), which is not large in the s ∼ Q 2 regime of Sec. III.1. Hence the gluon ( Fig. 10) and quark (Fig. 8) loop corrections at large-x are suppressed by a power of α s not enhanced by a large logarithm, and thus are outside of the precision of our approximation.
The construction of the small-x evolution of other TMDs of unpolarized proton or nucleus, such as the quark or gluon Boer-Mulders function, can be done similarly to the above, and are also governed by the nonlinear evolution of the correlators of fundamental and/or adjoint Wilson lines.

III.2.2. Evolution of Polarized Target TMDs: an Outline
The evolution of polarized target TMDs at small-x is somewhat different from the unpolarized case. Here we will only give the general outline of this evolution, with a more detailed description left for future work.
Let us concentrate on the quark TMDs. Again the diagrams of the type shown in Fig. 7 need to be considered, except now the quark-gluon vertices are not eikonal, allowing for spindependence to be transferred from the target nucleus to the produced quark. Note that, as a consequence, the contribution of the diagram has a factor of x suppression compared to its eikonal contribution to the unpolarized TMDs. These diagrams are supplemented by the same-order (both in α s and in x) graphs like that shown in Fig. 11. To include the small-x evolution correction into a generic polarized-target quark TMD one has to start by including the effects of GGM/MV multiple rescatterings. The relevant diagrams are shown in Fig. 12. There the red bar represents the shock wave again. Note that we are treating the rescattering on a nucleon carrying the spin information about the target separately, and draw it explicitly in the graphs. (The spin-dependent scattering may also contain gluon exchanges in the t-channel.) Placing such spin-dependent rescattering inside the shock wave rectangle, as shown in Fig. 12, implies that multiple rescatterings may happen both before and after the spin-dependent scattering. The same arguments apply to small-x evolution corrections to the diagram in Fig. 11 for the polarized TMD at hand. In the end one is left with the diagrams where all the evolution and interaction with the target happen entirely either to the left or to the right of the cut, as illustrated in Fig. 13. For simplicity we only show the linear evolution, realized through the quark ladder exchange shown in the left panel of Fig. 13, and the non-eikonal (non-BFKL) gluon ladder exchange shown in the right panel. The nonlinear evolution corrections are represented by interaction with the shock wave. 13. Evolution corrections to the polarized-target SIDIS process at small-x. The gluon ladder in the right panel is of the non-BFKL type [53,54]. The red bands represent the shock wave.
In the left panel of Fig. 13 the anti-quark in the amplitude emits hard gluons, while itself cascading down to lower and lower x until it interacts with the target. After the interaction with the target all the emitted gluons recombine back with the anti-quark. The same applies for the gluon ladder depicted in the right panel of Fig. 13. Combining each diagram with its complex conjugate we write a general polarized-target quark TMD as where λ and Σ are the quark and the target spin projections correspondingly (taken in either longitudinal or transverse basis) and Y ≈ ln 1/x at small-x.
Eq. (93) contains two main ingredients. The function f σσ ′ (x − z, y − z; λ) has to be separately calculated for every polarized TMD (using the γ * → qq light-cone wave function) by analogy to what was done in the unpolarized target case: this is left for future work. The interaction of the x, z (and y, z) dipole with the target, along with its evolution shown explicitly in Fig. 13, is included through the Reggeon exchange amplitude R σ σ ′ x,z (Y ) (R σ σ ′ * y,z (Y )). The evolution of R σ σ ′ x,z (Y ) involves mixing between the quark and gluon ladders, as shown in the right panel of Fig. 13. It is possible that the evolution of the (non-BFKL) gluon ladder depends on the TMD in question. For example, it is conceivable that the evolution of the quark transversity could be different from the quark helicity. While further investigation of this evolution for various polarized TMDs is left for future work, let us give an example of such evolution by considering only the quark sector.
By the quark sector we mean considering only the diagram in the left panel of Fig. 13, without any mixing with the gluon ladders. The quark ladder evolution is built up out of the q → qG splittings with the quark after the splitting being the softer particle. Corresponding small-x evolution kernel is diagonal in quark helicity. Hence the large-N c nonlinear evolution of R σ σ ′ x,z (Y ) for a ladder with quarks in the t-channel is independent of the anti-quark helicities σ, σ ′ on either side of the cut, and in the large-N c limit can be simply read off from the evolution of the flavor non-singlet structure function in [48] (see also [49][50][51][52][53][54][55]): x,z (y).
Here Y i is some initial rapidity for the evolution and r σ σ ′ x,z (Y ) is the initial condition. This initial condition is given by [48]  Eq. (94) is double-logarithmic in energy, that is, it resums powers of α s ln 2 s. Hence it contains R, which resums double-logarithms of energy, and S, resumming single logarithms [48].
Indeed, for consistency, the evolution for R has to be augmented by the single-logarithmic term: this is left for the future work. Eqs. (93) and (94) complete our outline for the general form of the small-x evolution of the polarized quark TMDs.
Note again that the full QCD Reggeon evolution (even at the linear level) also includes gluon ladders, mixing with the quark ladders, as shown in the right panel of Fig. 13. There one of the gluons carries the information about the polarization of the target [53,54]. Hence Eq. (94) is incomplete, and has to be augmented by a mixing term with the gluon ladder along with a separate evolution for the gluon ladder (also mixing with the quark ladder), similar to how it was done in [53,54]. While this is outside the scope of this work, here we note that the x-dependence found in [54] scales as

IV. CONCLUSIONS
In this work, we have presented calculations for the quark TMDs of an unpolarized target in the quasi-classical approximation and with leading-logarithmic quantum evolution. At face value the calculation applies to a heavy nucleus under the resummation α 2 s A 1/3 ∼ O (1), but its core features should also be valid for a high-energy proton with a large parton density generated by quantum evolution. Attempts to model the dense proton as a "nucleus" have been successfully applied in the past to phenomenology in inclusive deep inelastic scattering, so this approach may be a valuable tool in studying the TMDs of the proton as well. Regardless, the formalism presented here makes it possible to perform first-principles calculations of TMDs and spin-orbit structure within a controlled resummation of QCD.
The primary results of our calculation are as follows. The quasi-classical factorization formula (39) expresses the relationship between the quark TMDs of the nucleus, the quark TMDs of the nucleons, the distribution of nucleons in the nucleus, and the multiple rescattering on spectator nucleons. The parameterization (56) expresses the most general form of the nucleon Wigner distribution which is consistent with rotational invariance, parity, and time reversal. Together, these yield the decomposition of the unpolarized quark distribution f 1 (61) and the quark Boer-Mulders distribution h ⊥ 1 (78) in terms of the TMDs of the nucleons. The quasi-classical expressions can be taken as initial conditions for the subsequent quantum evolution. In the large-x regime, the leading evolution (83) is doublelogarithmic, resumming powers of α s ln 2 Q 2 in agreement with the Sudakov form factor and the Collins-Soper-Sterman evolution equation. In the small-x regime, the unpolarized quark distribution is given by (89), and its evolution is governed by BK-JIMWLK equations (90) and (91). The small-x evolution suggested by the polarized TMDs appears to be more complex, involving the familiar BK-JIMWLK evolution as an ingredient for the more intricate Reggeon evolution (94) (augmented by the mixing with the gluon ladders, as shown in Fig. 13).
Clearly further investigation of the perturbative QCD Reggeon evolution is warranted in the future. If the Reggeon evolution leads to polarized TMDs that may grow at small-x faster than 1/x, as Eq. (96) with the power of x taken from [54] appears to suggest, one should try to include saturation effects into this evolution. The inclusion of saturation effects may make the integral of the resulting TMD over small-x convergent, and would allow one to assess the size of the contribution of such an integral. A large contribution coming from the small-x region may help resolve the proton spin puzzle by identifying the phase space region containing the missing spin: this possibility is important and has to be explored in the future.
In our formalism, essentially all of the model dependence is encapsulated into the nonperturbative Wigner distribution of nucleons inside the nucleus. This distribution, however, is highly constrained by symmetry, making it possible to identify distinct channels which couple the TMDs of the nucleus to the TMDs of the nucleons. Our approach is also highly amenable to modeling and phenomenology. First one chooses an ansatz for the Wigner distribution such as (66); this then determines the form of the quasi-classical multiple scattering factor through the integral (63). With these ingredients, one can evaluate the intermediate integrals over p, b, and r which couple the nuclear TMDs to the nucleonic ones. If desired, one can then also choose initial conditions for the nucleonic TMDs, such as the quark target model or scalar diquark model of [57], allowing for explicit analytic or numerical calculation of the nuclear TMDs. These functional forms can then be evolved with the appropriate quantum evolution equations to the kinematics appropriate for comparison with experimental data. In future work, we would like to extend our computation of TMDs in the dense limit to the full set of leading-twist quark and gluon TMDs. This would provide another theoretical benchmark for understanding the range of spin-orbit physics permitted by QCD and would be well-suited to a global fit of available TMD data.
One novel feature of our calculation is that, due to the possibility of spin-orbit coupling in the distribution of nucleons, there can be mixing between different TMDs at the level of the nucleus and the nucleons. In (61) the nucleonic Sivers function mixed into the nuclear quark distribution, and in (78) the nucleonic transversity and pretzelosity mixed into the nuclear Boer-Mulders function. All three of these mixings are due to the same L · S N spin-orbit correlation with strength given by W OAM . The fact that the same underlying correlation is responsible for multiple mixings is a testable prediction of the theory: in principle, if the unpolarized quark distribution f N 1 and Sivers function f ⊥N 1T of the nucleon are known from experiment, then a further measurement of the unpolarized quark distribution f A 1 of the nucleus can allow an extraction of the spin-orbit coupling term W OAM which is present in the nucleus. This would then provide a prediction for the amount of mixing that should occur between the nucleonic transversity h N 1 or pretzelosity h ⊥N 1T and the Boer-Mulders function h ⊥A 1 of the nucleus. Such an extraction would require good coverage of the p T dependence of the nuclear TMDs, and is thus likely to only be accessible at a future electron-ion collider.
The TMD mixing observed here couples the P T -even and P T -odd sectors, mediated by the spin-orbit coupling mechanism and the depth dependence of the multiple scattering. This is a very general feature of TMDs in the dense limit which goes beyond the specific L · S N coupling that can occur in an unpolarized nucleus. A similar mechanism L · S A was already observed for the case of a transversely polarized nucleus, resulting in the mixing of the nucleonic unpolarized quark distribution into the nuclear Sivers function [27]. Similar features should again occur in the general case of a polarized nucleus, with a relatively small number of possible spin-orbit correlations generating a relatively large number of mixings between the nuclear and nucleonic TMDs. Indeed, one can imagine performing the same sort of analysis for the "generalized transverse-momentum-dependent parton distribution functions" (GTMDs) which are analogous operators to (1), taken between off-forward matrix elements of hadronic states (see, for example, [96]). The GTMDs are the "mother functions" from which one can obtain both the TMDs and the generalized parton distributions (GPDs).
If the same type of mixing occurs at the level of the GTMDs, it raises the interesting possibility of a single spin-orbit coupling term being responsible for different mixings in both the TMD and GPD sectors. The methodology presented here, which has at its foundation a genuine resummation of QCD, offers a new approach to the calculation of TMDs and related quantities, with bountiful applications to both theory and phenomenology.
symmetries "light-cone parity" P ⊥ and "light-cone time reversal" T ⊥ which preserve the z-axis and satisfy P ⊥ T ⊥ = P T [98].
To describe the properties of the light-front wave functions under three-dimensional rotations without invoking the interaction Hamiltonian, it is necessary to use the "covariant light-front dynamics" of [69] in which the quantization axis ω µ is kept arbitrary. If the theory is quantized at fixed ω · x, then rotations which transform not only the physical vectors, but also the quantization axis ω µ correspond to "kinematic" transformations which do not invoke the interaction Hamiltonian. Consequently, a light-front wave function for spinless particles is a Lorentz scalar when the quantization axis ω µ is also transformed. In general, the spin indices of a light-front wave function rotate according to their group representation, but they transform more simply in terms of the covariant spin vector S µ (22).
Since the "longitudinal direction" is now defined byn, the longitudinal momentum fraction α of a nucleon with momentum p µ in a nucleus with momentum P µ is and the vector R µ ≡ p µ − αP µ , defined so that ω · R = 0, selects out the momentum transverse ton: The variables α, R 2 are thus scalars under the generalized transformations which include ω µ .
The rotation properties of the wave functions are seen most clearly in the "constituent rest frame" (CRF). In light-front perturbation theory, there is no "rest frame" in the conventional sense, since one component of the momentum (p − ) is not conserved. It is thus impossible to simultaneously set the net spatial momentum P of the nucleus and i p i of the nucleons to zero. Instead, the momenta satisfy a modified conservation law, where τ is the deviation from the "energy shell." In the constituent rest frame, one sets the net spatial momentum of the constituents to zero, i p i = 0, with the parent particle forced to have nonzero momentum, P = τn.
In the constituent rest frame, the variables α, R 2 simplify: where we used (A8) and ω 2 = 0. The wave functions ψ(α, R 2 ) can therefore be expressed as ψ( p 2 , p ·n) which are equivalent invariants but have a simple meaning in the constituent rest frame. The covariant formulation thus makes explicit the dependence on the preferred directionn; this dependence is a feature of the relativistic wave function and in general cannot be avoided. However, there is one limit in which the dependence on the special directionn disappears: the non-relativistic limit p i 2 ≪ m 2 or c → ∞. In this limit, the light-front quantization condition ω · x = ct + ( x ·n) = const. reduces to the equal-time condition ct = const. and the dependence onn drops out. Consequently, in the non-relativistic limit the light-front wave functions should no longer depend onn: ψ( p 2 , p ·n) → ψ( p 2 ), restoring the rotational invariance of the physical vectors in the rest frame. For example, consider the simple wave function for the splitting of massive scalar particles shown in Fig. 14. Taking the tri-scalar vertex as λ, the standard calculation of the wave function gives ψ(k, p) = 1 2g where z = k + /p + is the momentum fraction of the scalar m 1 . In the constituent rest frame and expanding k 2 ≪ m 2 1 , m 2 2 , M 2 in the nonrelativistic limit (i.e., ( kc) 2 ≪ (mc 2 ) 2 for c → ∞), we obtain ψ( k 2 , k ·n) ≈ −λ 3m 2 1 − 4 k 2 3m 2 = ψ( k 2 ) (A12) as desired. Note that in the limit of nonrelativistic constituent motion, the rest frames of the parent particle and the constituents coincide, and we no longer need to distinguish between the "constituent rest frame" and the ordinary rest frame.
Indeed, the non-relativistic limit is precisely the limit which is relevant for the orbital motion of nucleons in a heavy nucleus. If we expand (A3) to lowest order in the nucleon velocities β = v/c = | k|/m N ≪ 1, then we have The Fourier factor exp[−i(p − p ′ ) · b] then becomes and in the rest frame such that Once again projecting out the polarized and unpolarized structures using (26), we obtain W ( p, b, S) which recovers the naive expectation of manifest three-dimensional rotation invariance in the nuclear rest frame.
The key is the change of variables between the boost-invariant quantities likeᾱ, P + b − and the 3-vectors p, b: where M A ≈ Am N . Using the dictionary (A21), we can impose rotational invariance in the rest frame through p, b to constrain the form of the Wigner distribution. Then we can change variables back to the boost-invariant quantitiesᾱ, P + b − which describe the highenergy collision.