Several Topics on Transverse Momentum-Dependent Fragmentation Functions

The hadronization of a high-energy parton is described by fragmentation functions which are introduced through QCD factorizations. While the hadronization mechanism per se remains uknown, fragmentation functions can still be investigated qualitatively and quantitatively. The qualitative study mainly concentrates on extracting genuine features based on the operator definition in quantum field theory. The quantitative research focuses on describing a variety of experimental data employing the fragmentation function given by the parameterizations or model calculations. With the foundation of the transverse-momentum-dependent factorization, the QCD evolution of leading twist transverse-momentum-dependent fragmentation functions has also been established. In addition, the universality of fragmentation functions has been proven, albeit model-dependently, so that it is possible to perform a global analysis of experimental data in different high-energy reactions. The collective efforts may eventually reveal important information hidden in the shadow of nonperturbative physics. This review covers the following topics: transverse-momentum-dependent factorization and the corresponding QCD evolution, spin-dependent fragmentation functions at leading and higher twists, several experimental measurements and corresponding phenomenological studies, and some model calculations.


Introduction
Quantum chromodynamics (QCD) [1] is known as the fundamental theory of strong interaction in the framework of Yang-Mills gauge field theory [2]. As a key property of QCD, the color confinement prohibits direct detection of quarks and gluons, the fundamental degrees of freedom, with any modern detectors. The emergence of color neutral hadrons from colored quarks and gluons is still an unresolved problem and has received particular interest in recent years [3]. With the progress of QCD into the precision era, unraveling the hadronization mechanism in the high-energy scattering processes has become one of the most active frontiers in nuclear and particle physics.
Due to the nonperturbative nature of QCD, it is still challenging to directly calculate the hadronization process from first principles. Similar to the parton distribution functions (PDFs) [4,5], which were originally defined as the probability density of finding a parton inside the parent hadron, the concept of fragmentation functions (FFs) was introduced by Berman, Bjorken, and Kogut [6] right after the parton model to describe the emergence of a system of the hadron from a high-energy parton isolated in the phase space. An alternative name, the parton decay function, has also frequently been used in early literature.
The modern concept of FFs in QCD was first introduced to describe the inclusive production of a desired hadron in the e + e − annihilation [7,8], which is still the cleanest reaction currently available to investigate the fragmentation process. Within the QCDimproved parton model, the FF has its foundation in the factorization theorem [9,10], in which the differential cross section is approximated as a convolution of short-distance hard scattering and long-distance matrix elements with corrections formally suppressed by inverse powers of a hard scale, e.g., the center-of-mass (c.m.) energy Q = √ s in the e + e − annihilation. The predictive power of this theoretical framework relies on the control of the hard probe, which can be achieved by our ability to calculate the partonic cross section order by order in the perturbation theory, and the universality of the long-distance functions, such as the FFs, to be tested in multiple high-energy scattering processes.
For a single-scale process, e.g., e + e − → hX, where h represents the identified hadron in the final state and X denotes the undetected particles, the process is not sensitive to the confined motion of quarks and gluons in the hadronization process, and one can apply the colinear factorization with the emergence of the detected hadron described by a colinear FF D f →h (z), where the subscript f stands for the parton flavor and z is the longitudinal momentum fraction carried by the hadron h with respect to the fragmenting parton. If two hadrons are identified in a process, e.g., e + e − → h A h B X, where h A and h B are detected hadrons in the final state, the reaction becomes a double-scale problem with one scale Q given by the hard probe and the other scale provided by the transverse momentum imbalance, |p A⊥ + p B⊥ |. When the second scale is much smaller than Q, i.e., the two hadrons are nearly back to back, one needs to use the transverse-momentum-dependent (TMD) factorization. The emergence of each of the hadrons is described by a TMD FF D f →h (z, k ⊥ ), where k ⊥ is the transverse momentum of the fragmenting parton with respect to the observed hadron [8,11]. When the two scales are compatible, the reaction effectively becomes a single-scale process, and one can again use the colinear factorization. The matching between the two regions has been developed. The TMD FFs defined in the e + e − annihilation also play an important role in the study of nucleon three-dimensional structures via the semi-inclusive deep inelastic scattering (SIDIS) process [12]. Instead of identifying two hadrons in a reaction, one can also access TMD FFs in the single-hadron production process by reconstructing the thrust axis, which provides the sensitivity to the transverse momentum of the observed hadron, as proposed in recent years [13][14][15][16].
Taking the parton spin degree of freedom into account, one can define polarized or spin-dependent TMD FFs. They essentially reflect the correlation between parton transverse momentum and its spin during the hadronization process and result in rich phenomena in high-energy scattering processes. For example, the Collins fragmentation function H ⊥ 1 (z, k ⊥ ) [17], naively interpreted as the probability density of a transversely polarized quark fragmenting into an unpolarized hadron, can lead to a single spin asymmetry (SSA) in the SIDIS process with a transversely polarized target [18]. This asymmetry is a key observable for the determination of the quark transversity distribution, the net density of a transversely polarized quark in a transversely polarized nucleon. It also leads to azimuthal asymmetries in e + e − annihilation as measured by Belle, BaBar, and BESIII. The progress of experimental techniques to determine the spin state of produced hyperons, such as Λ and Ω, and vector mesons, such as ρ and K * , offer us the opportunity to extract additional information from FFs. This is far beyond a trivial extension since the spin has been proven to be a powerful quantity to test theories and models, especially in hadron physics. The recent measurement of the spontaneous polarization of Λ from unpolarized e + e − annihilation is such an instance [19]. This observation can be explained by a naively time-reversal odd (T-odd) TMD FF D ⊥ 1T (z, k ⊥ ) and has received interests from various groups [15,[20][21][22][23][24][25][26][27][28][29][30][31][32].
In addition to the leading-twist FFs, which usually have probability interpretations, the high-twist FFs have been found o be much more important than expected in recent years for understanding precise experimental data . Although the colinear factorization at subleading power was demonstrated some time ago, the TMD factorization beyond the leading power is still under exploration, and some approaches have been proposed [67][68][69][70][71][72][73][74][75]. Although high-twist contributions are formally power suppressed, their contributions to the cross section might not be negligible and may have significant effects in certain kinematics or observables. The inclusion of high-twist FFs will also modify the evolution equation and consequently affect the leading-twist FFs. The TMD factorization at subleading power was recently explored with different approaches. Overall, many efforts, both theoretical and experimental, are still required to understand the hadronization process and the upcoming data from future electron-ion colliders.
The remainder of this review is organized as follows. In Section 2 we use e + e − → h A h B X as an example to present the flow of deriving the TMD factorization and the QCD evolution equation of TMD FFs. In Section 3, we present the FFs up to the twist-4 level for spin-0, -1/2, and 1 hadron productions. In Section 4, we summarize the experimental measurements towards understanding the spin-dependent FFs. In Section 5, we briefly lay out some model calculations. A summary is given in Section 6.

Factorization and Evolution
The modern concept of FFs has established on the QCD factorization theorems, which can be derived either from calculating traditional Feynman diagrams in perturbative field theory [10,11,[76][77][78][79][80][81] or in effective theories [82][83][84][85][86]. In the former approach, one first identifies a collection of Feynman diagrams that offers the leading contribution through the Libby-Sterman analysis [87,88]. In this method, the leading contribution is represented by the reduced diagrams.
Taking e + e − → h A h B X process with h A , h B traveling along almost back-to-back directions as an example [11], the leading regions are presented in Figure 1. The cross section is the product of various ingredients, such as the hard part H, the soft part S, and the colinear parts J A , J B . We work in the light-cone coordinate, so that a four momentum p can be written as follows: In the kinematic region where TMD factorization applies, the transverse momentum is considerably small compared with that along the longitudinal direction. Therefore, the momenta of the almost back-to-back hadrons A and B scale as p A ∼ Q(1, λ, √ λ) and p B ∼ Q(λ, 1, √ λ), where Q is the large momentum scale and λ ≪ 1 is a small parameter. The hard part H computes the cross section of interaction among hard partons whose momenta scale as Q(1, 1, 1) in perturbative field theory. The contribution from colinear partons whose momenta are colinear with the final state hadrons A and B are evaluated in the colinear function J A/B . This process results in the gauge invariant bare FFs. The soft part calculates the contribution from soft gluons whose momenta typically take the form of Q(λ, λ, λ). They will be absorbed into the definition of TMD FFs eventually and convert the bare FFs into the renormalized ones.
The interactions between different parts can be eliminated via applying appropriate kinematic approximations and the Ward identity. Finally, the cross section is given by a convolution of those well-separated parts, and we arrive at the factorization theorem of this process.
Depending on the physics of interest, we may derive either colinear factorization or transverse-momentum-dependent (TMD) factorization theorems. For the differential cross section of e + e − → h A h B X as a function of the relative transverse momentum between h A and h B , the TMD factorization theorem applies.
In the single-photon-exchange approximation, the differential cross section of this process can be written as the production of a leptonic tensor and a hadronic tensor. It reads as follows: [43] dσ where, α is the coupling constant, N c = 3 is the color factor, Q is the center-of-mass energy of the colliding leptons, y = (1 + cos θ)/2 with θ the angle between incoming electron and the outgoing hadron h A , z A and z B are light-cone momentum fractions of h A and h B , and P A⊥ is the transverse momentum of h A with respect to the direction of the h B momentum. For the unpolarized lepton beams, the leptonic tensor L µν is given by the following: with l 1 and l 2 being the momenta of colliding leptons. The hadronic tensor W µν contains nonperturbative quantities and is laid out as follows: where q ⊥ = −P A⊥ /z A , and H f (Q, µ) is the hard scattering factor that can be evaluated in the perturbative QCD. Here, D h A 1q (z A , p A⊥ ; µ, ζ A ) is the TMD FF with p A⊥ the transverse momentum of hadron with respect to the fragmenting quark direction, µ is a renormalization scale, and ζ A is a variable to regularize the rapidity divergence. Notice that k i,⊥ is the relative transverse momentum of the fragmenting parton with respect to the hadron momentum. Therefore, we have p i,⊥ by p i,⊥ = −z i k i,⊥ . Please also notice the difference between P A⊥ and p A⊥ . The three-dot symbol stands for various spin-dependent terms which are not explicitly shown.
It is more convenient to perform the TMD evolution in the coordinate space than in the momentum space. Therefore, we need the Fourier transform, to translate the TMD FF into coordinate space one. The hadronic tensor then becomes the following: The TMD FF in the coordinate space is defined as the product of transition matrix elements between the vacuum and the hadronic final states.
Before presenting the final definition for the TMD FF in the coordinate space, we first show the unsubtracted version, which appears in the LO calculation. For the production of hadron A, it reads as follows: where the position vector x = (0, x − , b T ) contains only minus and transverse components, is the rapidity of hadron A, Tr C is a trace in the color space, and Tr D is a trace in the Dirac space. The direction of the Wilson lines in the FF of hadron A is specified by the direction of hadron B which is denoted as n B and vice versa. Notice that the rapidity parameters y A → +∞ and y B → −∞ are introduce, so that n A = (1, −e −2y A , 0 T ) and n B = (−e 2y B , 1, 0 T ) are slightly space-like. Please also notice the difference between y p A and y A . The Wilson line starting from the position x is defined as follows: with a and b being the color indices, and g 0 and A α (0) the bare coupling and the bare gluon field.
Taking the y A → +∞ and y B → −∞ limit and absorbing the soft factors into the unsubtracted TMD FF, we arrive at the final definition of the TMD FF: where y n is an arbitrary rapidity introduced to separate ζ A ≡

Evolution Equations for TMD FFs
To regularize the ultraviolet (UV) and rapidity divergences, the energy scale µ and √ ζ are introduced. As a consequence, the TMD FFs differ at different energy scales. The evolution effects are important for phenomenological studies. The QCD evolution for TMD FFs with respect to ζ is controlled by the Collins-Soper (CS) equation [10,11], which is given as follows: withK(b T ; µ) being the CS evolution kernel. The scale dependence of the evolution kernel is governed by where γ K (µ) is the anomalous dimension. It is given by γ K (µ) = 2C F π α s (µ) with C F = 4/3 being the color factor and α s being the running coupling at the LO accuracy [12].
The µ dependence of the TMD FF is then given by where γ D is another anomalous dimension. At the LO accuracy [12], it is given as follows: . The TMD FF defined by Equation (8) is actually calculable in the colinear factorization approach in the small-b T regime. However, at the large b T region, the discrepancy of these two approaches grows in terms of Λ QCD b T . This region is usually referred to as the nonperturbative regime since a large coordinate corresponds to a small energy scale. The perturbative treatment of the QCD evolution in this region is no longer reliable. To have a consistent formula, the b * -prescription is usually adopted in phenomenology. By we can separate the perturbative part from the nonperturbative part in the QCD evolution. Here, γ E is the Euler constant, and b max is an infrared cutoff which is properly chosen to guarantee that µ b ≫ Λ QCD . Employing the b * prescription, the QCD evolution is always performed in the realm of the perturbative QCD. Therefore, this approach underestimates the contribution from the nonperturbative regime. This part of the contribution can be reintegrated into the final prescription by the introduction of a nonperturbative factor.
Ultimately, we arrive at [12] where the last line is the nonperturbative function that returns the nonperturbative effect that has been deliberately removed from the QCD evolution in the b * -prescription. There is no theoretical approach that can evaluate this nonperturbative function other than the one that extracts it from experimental data [89][90][91][92][93][94][95][96][97][98]. Notice that [99] present a different method to address the nonperturbative physics. Here, is the FF at the initial scale. In the phenomenology, it is usually chosen to coincide with the colinear FF D h 1q (z, µ f ) with the factorization scale specified by µ f = µ b . Similar to the PDF case, QCD evolution tends to broaden the k T distribution width at higher energy scales. Both unpolarized and spin-dependent FFs show such a behavior [100].

TMD Factorization at the Higher Twist
In a semi-inclusive process, normally we can find two energy scales: the typical transverse momentum q ⊥ and the hardest energy scale Q. In the region of Q ≫ q ⊥ ≫ Λ QCD , the TMD factorization framework at the leading twist usually works very well. When q ⊥ ∼ Q ≫ Λ QCD , we should fall back to the colinear factorization. However, in between, there is still a large phase space where q ⊥ is smaller than Q but not much smaller. This is the kinematic region where both the TMD factorization and the colinear factorization can approximately apply. However, the prediction from the TMD factorization deviates from the experimental measurements when q ⊥ /Q becomes not very small, calling for the inclusion of higher twist corrections. The higher twist corrections are also usually referred to as power corrections since they provide contributions in terms of (q ⊥ /Q) n . In addition, twist-3 contributions usually introduce new asymmetries that do not appear at the leading twist level. A comprehensive study on the higher twist contributions is thus vital in phenomenology. Contributions from higher twist TMD PDFs and FFs were studied some time ago [42], but few advances have been made in the systematic derivations of the TMD factorization formula at the higher twist gain. Various theoretical methods have been applied to derive the TMD factorization scheme at the twist-3 level, such as the TMD operator expansion technique [67][68][69], the soft-colinear effective theory approach [70], factorization from functional integral [71][72][73][74], and a very recent work from [75], etc. The TMD factorization at the higher twist level is far from being completed, which requires further theoretical efforts.

Spin-Dependent TMD FFs
In semi-inclusive reactions, the experimental observables are usually different azimuthal asymmetries. In the kinematic region of TMD factorization, they are directly linked to TMD PDFs or FFs. The transverse momenta of partons and hadrons are often entangled with their polarizations. As a consequence, there are abundant polarization-dependent azimuthal asymmetries that can be measured in the experiment. This is particularly true for the transverse polarization. It is thought to provide only subleading power contributions compared to the longitudinal polarization at high energy; however, it often generates leading power contributions when correlated with the transverse momenta. In this section, we summarize the definition of the spin-dependent TMD FFs for hadrons with different spins. The following discussion only applies at the LO level since the TMD factorization for higher twist contributions is still far from being concluded. Therefore, we remove the scale dependence from TMD FFs.

The Intuitive Definition of TMD FFs
FFs represent the momentum distribution of a hadron inside of a hadronic jet produced by the fragmenting high-energy parton. We use D q→h (k; p) to denote the probability density of producing a hadron h with momentum p from a quark with momentum k.
In the high-energy limit, we can safely neglect the quark and hadron mass. Therefore, we have k 2 = p 2 = 0. In the naive parton model picture, the hadrons move colinearly with the parent quark. We thus have p = zk where p is the hadron momentum, k is the quark momentum, and z is the momentum fraction. In this case, the FF is only a scalar function of z. We have where D q→h 1 (z) is simply the unpolarized FF. With the spin degree of freedom being taken into account, the FFs will also depend on additional parameters which characterize the polarization of the final state hadron or the fragmenting quark. For example, for the production of spin-1/2 hadrons, we need to introduce λ q and λ to describe the helicities and introduce ⃗ s Tq and ⃗ S T to describe the transverse polarizations of the quark and the hadron. With more available parameters, we can construct two additional scalar structures, λ q λ h and ⃗ s Tq · ⃗ S T , according to the parity conservation. Therefore, the complete decomposition of the FF is given by where G 1L (z) and H 1T (z) are the longitudinal and transverse spin transfers from the quark to the hadron, respectively. The physical interpretations of these probability densities coincide with those of the leading twist FFs in the colinear factorization approach. In some cases, the transverse momentum of the final state hadron with respect to the quark momentum becomes relevant to the observable of interest. The interplay between the transverse momentum p ⊥ and the polarization parameters induces considerably intriguing phenomena. Again, we use the spin-1/2 hadron production as an example. From the parton model, we obtain the following eight TMD probability densities: Here, we have dropped the q → h superscript for simplicity. These TMD FFs correspond to the eight leading twist TMD FFs defined in the TMD factorization approach. Among them, we notice in particular the famous Collins function H ⊥ 1 [17] and the Siverstype FF D ⊥ 1T [101,102]. They are usually referred to as the naive-T-odd FFs. In neglecting the interaction among the final state hadrons and the gauge link (which will be explained below), the time-reversal invariance demands that these two functions disappear. However, the time-reversal operation converts the "out" state to the "in" state. The interaction among hadrons suggests that one cannot find a simple relation between the "in" and "out" states any longer. Therefore, the time-reversal invariance actually poses no constraints on FFs. This feature can be fully appreciated in the context of parton correlators in the next subsection. Furthermore, we use H to denote FFs accompanied with the transverse polarization of the fragmenting quark⃗ s Tq . They are chiral-odd FFs. The reason for this will also be explained later.

The Definition of TMD FFs from the Parton Correlators
In the language of quantum field theory, the quark FFs are defined via the decomposition of parton correlators, such as the quark-quark correlator and the quark-gluon correlator. Usually, we need to define the gauge-invariant quark-quark correlators in the very beginning. From [8,42,43,50,64], we have the following: (17) where ξ is the coordinate of the quark field, k and p denote the 4-momenta of the fragmenting quark and the produced hadron, respectively; S denotes the hadron spin; and L(ξ; ∞) is the gauge link that ensures the gauge invariance of the definition of the correlator. We use i and j to represent one component of the corresponding spinor. Therefore,Ξ (0) ij (k; p, S) is actually one element in a 4 × 4 matrix which is denoted byΞ (0) (k; p, S).
As for the TMD FFs, we can integrate the above master correlator over the k − component and obtain the following TMD quark-quark correlator: where z = p + /k + is the longitudinal momentum fraction of the hadron, and k ⊥ is the transverse momentum of the fragmenting quark with respect to the hadron momentum. Unlike the discussion in the previous sections, it is more convenient to express the parton correlators as a function of k ⊥ instead of p ⊥ . Nonetheless, since we have the approximation k ⊥ = −p ⊥ /z, these two methods are equivalent.
Although the TMD quark-quark correlator is a nonperturbative object, we can still discuss some general features from the definition. For instance, it possesses hermiticity, parity invariance, and charge-conjugation symmetry. As will be shown below, these properties will constrain the structures of the correlator. However, unlike the case for PDFs, the time-reversal invariance does not mean much for FFs. Furthermore, the quark-quark correlator is a 4 × 4 matrix in the Dirac space. Therefore, it can always be decomposed in terms of 16 Γ-matrices, i.e., The coefficient functions αβ are given by the trace of the corresponding Γ-matrix with the correlator. These coefficient functions can further be decomposed into the products of scalar functions with basic Lorentz covariants according to their Lorentz transformation properties. The basic Lorentz covariants are constructed in terms of the available kinematic variables used in the reaction process. The scalar functions are the corresponding TMD FFs. We will present the detailed decomposition in the following subsections. Notice that the TMD quark-quark correlator given by Equation (18) satisfies the constraints of hermiticity and parity conservation. This will limit the allowed Lorentz structures of the parton correlator.
Higher twist TMD FFs also receive contributions from quark-gluon correlators [42,43,50,64] in addition to the quark-quark correlator mentioned above. For example, the complete decomposition of twist-3 TMD FFs also involves contributions from the following correlator: where D ρ (y) ≡ −i∂ ρ + gA ρ (y) and A ρ (y) denote the gluon field. However, the twist-3 TMD FFs defined via these quark-gluon correlators are not independent from those defined via the quark-quark correlator [64,65]. They are related to each other by a set of equations derived using the QCD equation of motion γ · D(y)ψ(y) = 0. Therefore, we will only show the explicit decomposition of the TMD FFs from the quark-quark correlator in the following subsections.

The Spin Dependence
With the spin degree of freedom being taken into account, the basic Lorentz covariants in the decompositions of the coefficient functions in Equation (19) depend on not only momenta but also parameters describing the hadron polarization. The hadron polarization is defined in the rest frame of the hadron and is described by the spin density matrix.
For spin-1/2 hadrons, the spin density matrix is given by where ⃗ σ is the Pauli matrix, and ⃗ S is the polarization vector in the rest frame of the hadron. The covariant form of the polarization vector reads as follows: Here, M is the hadron mass, λ is the helicity, and S µ T is the transverse polarization vector of the hadron. We have employedn µ to represent the light-cone plus direction and n µ to denote the minus direction. For spin-1/2 hadrons, an additional pseudo-scalar λ and an axial-vector S µ T are at our disposal for constructing the basic Lorentz tensors. For spin-1 hadrons, such as the vector mesons, the polarization is described by a 3 × 3 density matrix, which is usually given as [47] Here, Σ i is the spin operator of spin-1 particle. The rank-2 tensor polarization basis Σ ij is defined by where the second term subtracts the diagonal elements from the product in the first term to give the relation This can be easily seen from the square of the spin-1 operator, i.e., (23), we find that a polarization tensor T is required to fully describe the polarization of a vector meson besides the polarization vector S. The polarization vector S is similar to that of spin-1/2 hadrons. It takes the same covariant form as laid out in Equation (22). The polarization tensor T ij = Tr(ρΣ ij ) has five independent components that consist of a Lorentz scalar S LL , a Lorentz vector S The Lorentz covariant form for the polarization tensor is expressed as [47] T µν = 1 2 where we have used the shorthand notation For spin-3/2 hadrons, such as the decuplet baryons, the polarization is described by a 4 × 4 density matrix which is given by [103,104] Here, Σ i is the spin operator of the spin-3/2 particle, and S i is the corresponding polarization vector. Similar with that for spin-1 case, (Σ ij ) is the polarization tensor basis which has five independent components. It can be constructed from Σ i and is given by Notice that the square of the spin-3/2 operator is given by ∑ i (Σ i ) 2 = 3 2 ( 3 2 + 1)1 = 15 4 1. The rank-2 tensor polarization basis for spin-3/2, Σ ij , is also chosen to be traceless as laid out by Equation (25). Therefore, the second term in Equation (29) is different from that in Equation (24) for spin-1 hadrons. The corresponding polarization tensor T ij also has five independent components which are the same as those for spin-1 hadrons. The rank-3 tensor polarization basis Σ ijk is unique for spin-3/2 hadrons. It has seven independent components which can be constructed as follows: where the symbol {· · · } stands for the sum of all possible permutations. The corresponding rank-3 spin tensor R ijk is defined as follows: Meanwhile, the Lorentz covariant form is given as follows:

Decomposition Result for Spin-Dependent TMD FFs
The results for TMD FFs of spin-1 hadrons defined via quark-quark correlator exist up to twist-4 level in the literature [65]. The leading twist TMD FFs for spin-3/2 hadrons have also been presented in [104]. In this section, we summarize the general decomposition of the quark-quark correlator in terms of TMD FFs for the unpolarized part, polarizationvector-dependent part, rank-2-polarization-tensor-dependent part, and rank-3-polarizationtensor-dependent parts. To describe the production of pseudoscalar mesons, we only need the unpolarized part. To describe the production of baryons, we need to combine the unpolarized and the polarization-vector-dependent parts. The description of the spin-3/2 hadron production requires all four parts. However, it should be noted that different conventions are employed in different works.
The notation system for TMD FFs in this review are laid out here. We use D, G, and H to denote FFs of unpolarized, longitudinally polarized, and transversely polarized quarks, respectively. They are obtained from the decomposition of the γ µ , γ 5 γ µ and γ 5 σ µν terms of the quark-quark correlator. Those FFs defined from the decomposition of the 1 and γ 5 terms are denoted as E. We use the numbers 1 and 3 in the subscripts to denote the leading twist and twist-4 FFs, respectively. Other FFs without numbers in the subscripts are at the twist-3 level. The polarization of the produced hadron will be specified in the subscripts, where L and T represent longitudinal and transverse polarizations, and LL, LT, and TT stand for the rank-2-tensor polarizations. The symbol ⊥ in the superscript implies that the corresponding basic Lorentz structure depends on the transverse momentum k ⊥ .
The decomposition for the unpolarized part is given by the following: Here,k ⊥α ≡ ε ⊥µα k µ ⊥ denotes the transverse vector orthogonal to k ⊥α , with ϵ ⊥µν being defined as ε ⊥µν ≡ ε µναβn α n β . There are eight TMD FFs for the unpolarized part. Among them, the number density D 1 and the Collins function H ⊥ 1 are at the leading twist. They both have twist-4 companions i.e., D 3 and H ⊥ 3 , respectively. The other four are twist-3 FFs. The TMD FFs D ⊥ 1T ,G ⊥ , H ⊥ 1 , H, and H ⊥ 3 are usually referred to as the naive T-odd FFs. The reader may have already discerned that the T-odd FFs are always associated with the Levi-Civita tensor, ε µναβ . It should be noted that T-odd PDFs can only survive thanks to the gauge link. However, for the FFs, the final state interactions between the produced hadrons in the hadronization process can also contribute to the T-oddness. This difference has a more important impact on the polarization-vector-dependent T-odd PDFs and FFs, which are discussed below.
The decomposition for the vector polarized part is given by the following: There are in total 24 polarization-vector-dependent TMD FFs. Of these, 6 contribute at the leading twist, 12 at twist-3, and remaining 6 at twist-4. Among the six leading twist FFs, G 1L is the longitudinal spin transfer, H 1T and H ⊥ 1T are transverse spin transfers, G ⊥ 1T is the longitudinal to transverse spin transfer, H ⊥ 1L is the transverse to longitudinal spin transfer, and D ⊥ 1T induces the transverse polarization of hadrons in the fragmentation of an unpolarized quark. We note in particular that the D ⊥ 1T FF resembles the Sivers function in PDFs [101]. It is responsible for the hadron transverse polarization along the normal direction of the production plane in high-energy collisions. It is also a naive Todd FF. However, as mentioned above, the T-oddness has little meaning in the context of hadronization. The T-odd PDFs arise solely from the gauge link. Therefore, it has been proven theoretically that there is a sign-flip between the Sivers functions in SIDIS and Drell-Yan [105][106][107]. However, the T-oddness of FFs can also be generated from the interaction among final state hadrons. Therefore, there is no such similar relation for the D ⊥ 1T FF between different processes. Besides D ⊥ 1T , there are seven other T-odd FFs, namely, The rest are T-even. All of T-odd FFs are accompanied by the Levi-Civita tensor except for E L and E ′⊥ T . The decomposition for the rank-2-polarization-tensor-dependent part is given as follows: We have used the shorthanded notations such as S kk TT ≡ S αβ TT k ⊥α k ⊥β . There are in total 40 tensor polarization-dependent TMD FFs; of these.10 contribute at the leading twist, 20 contribute at twist-3, and the remaining 10 contribute at twist-4. The 24 TMD FFs defined from the decomposition ofΞ T(0) α and Ξ T(0) ρα are naive T-odd. Among these TMD FFs, we notice in particularly that the S LL dependent TMD FF D 1LL , which is responsible for the spin alignment of the produced vector meson, is decoupled from the quark polarization. This suggests that the vector meson spin alignment can also be observed in the unpolarized high-energy collisions [108][109][110]. Besides, D 1LL also survives the k ⊥ -integral. Therefore, it also appears in the colinear factorization.
The rank-3-polarization-tensor-dependent TMD FFs are unique for spin-3/2 (or higher) hadrons. A complete set of leading twist quark TMD FFs for spin-3/2 hadrons has been given in [104]. There are in total 14 rank-3-polarization-tensor-dependent TMD FFs that can be defined at the leading twist level. We refer interested readers to [104] for a detailed discussion.

TMD FFs of Antiquarks and Gluons
One can define antiquark TMD FFs by replacing the fermion fields in the correlator of quark TMD FFs with the charge-conjugated fields. Therefore, it is easy to find that the traces of the correlator with Dirac matrices I, iγ 5 and γ µ γ 5 will have an opposite sign between quark and antiquark cases, while the traces with γ µ and iσ µ γ 5 are the same [42,43,56]. The definition and parameterization of the antiquark TMD FFs are then full analogous to those of quark TMD FFs.
The gluon FFs are defined through the gluon correlator given by [8,111] Γ µν;ρσ (k; p, S) = ∑ X d 4 ξ (2π) 4 e ik·ξ ⟨0|F ρσ (ξ)|p, S; X⟩⟨p, S; X|U (ξ, 0)F µν (0)|0⟩, (48) where F ρσ (ξ) ≡ F ρσ,a T a is the gluon field field strength tensor, and U (ξ, 0) is the Wilson line in the adjoint representation that renders the correlator gauge invariant. Under the assumption that the fragmenting parton moves in the plus direction, an integration over the k − component is carried out to give the TMD gluon correlator. At the leading twist, we need to consider where i and j are transverse Lorentz indices in the transverse directions. For the spin-1/2 hadron production, there are eight leading twist gluon TMD FFs which are given by the decomposition of the TMD gluon correlator [111]. We have the following: Γ U ,Γ L , andΓ T stand for the unpolarized, the longitudinal, and transverse polarized parts for the hadron production, respectively. Analogously to the quark FFs, we have used D to represent FFs of the unpolarized gluons, G to represent the FFs of the circularly polarized gluons, and H to represent the FFs of the linearly polarized gluons. Higher twist gluon TMD FFs are also discussed in [111], who further detail the parameterizations.

Experiment and Phenomenology
In high-energy experiments, the polarization of final state hadrons is usually measured from the angular distribution of their decay products. It is very challenging to acquire accurate experimental data. In light of a considerably large amount of free parameters, the spin-dependent FFs are not well-constrained experimentally. Compared with the case for unpolarized PDFs or FFs, the quantitative study of spin-dependent FFs is still immature. That said, there are already quite a few phenomenological studies making full use of the available experimental data. In this section, we summarize the available experimental data and the corresponding phenomenological studies.

Λ Hyperons
The polarization of Λ 0 hyperons is usually measured from the angular distribution of the daughter proton in the parity-violating Λ 0 → p + π − decay channel. In the rest frame of Λ 0 , the normalized angular distribution of the daughter proton reads as follows: where α = 0.732 ± 0.014 is the decay parameter of Λ [112], P is the polarization of Λ along a specified direction, and θ * is the angle between the proton momentum and the specified direction to measure the Λ polarization. The LEP experiment is an e + e − collider at the Z 0 -pole. Due to the parity violation in the weak interaction, the produced quark and antiquark are strongly polarized along the longitudinal direction. The longitudinal polarizations of those final state quarks and antiquarks in e + e − annihilation at different collisional energies can be easily computed at the LO level and are explicitly shown in [113]. At the Z 0 -pole, the longitudinal polarization of the final state down-type quarks can reach 0.9. That of the up-type quarks is a bit smaller but is still about 0.6 ∼ 0.7. Based on the SU(6) spin-flavor symmetry, the polarization of Λ 0 is determined by the polarization of the s quark. It is thus proposed in [114] that the final state Λ 0 hyperons are also strongly polarized at LEP, and the measurement of this polarization can probe interesting information on the hadronization mechanism. In the language of QCD factorization, the LEP experiment is the ideal place to study the longitudinal spin transfer G 1L (z), which represents the number density of producing longitudinally polarized Λ 0 hyperons from longitudinally polarized quarks. It is the p T -integrated version of the TMD FF G 1L (z, p ⊥ ).
At the leading order and leading twist, the longitudinal polarization of Λ 0 reads as follows: [108,113] where λ q (y) = ∆ω q (y)/ω q (y) is the helicity of the fragmenting quark with ∆ω q (y) and ω q (y) being defined as follows: Here, y = (1 + cos θ)/2 with θ is the angle between the outgoing Λ and the incoming electron. The coefficient functions are given as A(y) = (1 − y) 2 + y 2 , B(y) = 1 − 2y, A are the coupling constants of the vector current and axisvector current parts of the quark/electron, with Z 0 . M Z being the mass of Z 0 and Γ Z being the width. Notice that λq(y) = −λ q (1 − y). The quark helicity and antiquark helicity have the opposite sign. This is in line with the sign flip in Section 3.5. Therefore, the polarization ofΛ 0 is expected to have the opposite sign with that of Λ 0 .
Since the quark helicity λ q (y) and the production weight ω q (y) are calculable in quantum field theory, the measurement of the longitudinal polarization of final state Λ 0 as a function of z can directly provide information of the longitudinal spin transfer. Such experiments were eventually carried out by ALEPH and OPAL collaborations at LEP in the 1990s [115,116]. As shown in Figure 2, the longitudinal polarization increases monotonically with increasing z, which provides a hint on how to parameterize the longitudinal spin transfer.  [115] and OPAL [116] collaborations at LEP. We have combined the statistical and systematic errors. Neglecting the mass of Λ hyperons in the high-energy limit, the definitions of z in these two experiments are the same as those of the momentum fraction in the light-cone coordinate currently used in the QCD factorizations.
Following the release of these experimental data, many phenomenological studies [108,[117][118][119][120][121][122][123][124] were carried out to understand the longitudinal spin transfer G 1L (z). Among them, the de Florian-Stratmann-Vogelsang (DSV) parameterization [118] offers three scenarios. The first scenario is based on the naive parton model, which assumes that only the s quark contributes to the longitudinal spin transfer at the initial scale. The second scenario assumes that the u and d quarks contribute to negative G 1L (z) at the initial scale. The third scenario assumes that u, d, and s contribute equally. All three can describe the experimental data reasonably well. A more recent Chen-Yang-Zhou-Liang (CYZL) analysis [108] also obtained a good description of the experimental data utilizing the LO formula. The ambiguity again highlights the difficulties in the quantitative study of FFs. It can only be removed through a global analysis of the experimental data in various highenergy reactions. Therefore, many works have also made predictions for the longitudinal polarization of Λ produced in polarized SIDIS [117,[125][126][127] and pp collisions [128][129][130][131].
The inclusive DIS process with the polarized lepton beam has been used to probe the spin structure of the nucleon [132][133][134]. In this process, only the momentum of the final state lepton was measured. Therefore, we can gain information on the nucleon structure but lose those on the hadronization. To restore the access to (spin-dependent) FFs, we have to rely on the semi-inclusive process and measure (the polarization of) one final state hadron (There are two fragmentation regimes in SIDIS, namely the current fragmentation and the target fragmentation. Although the target fragmentation function is also currently a hot topic, it is beyond the scope of this review. We only focus on the study of the current fragmentation function). However, it is not a simple task to do so in the real world. Despite the difficulties, early attempts from the E665 [135] and HERMES [136] collaborations were still successfully performed. Recent measurements from HERMES [137] and COMPASS [138] collaborations have also elevated the quality of experimental data to a level that sheds light on phenomenological studies. These experiments measure the spin transfer coefficient D LL (z) (it is important to not get confused with the spin-alignment-dependent FF D 1LL (z) of vector mesons) which, at the leading order and leading twist approximation, is given by [42] with f 1,q (x B ) being the unpolarized PDF. Due to the presence of the unpolarized PDF of proton/nucleus, the polarized SIDIS experiment favors more contributions from the u and d quarks at large x B than from the e + e − collider. We show the HERMES data set as a function of z (integrating over x B ) and the COMPASS data set as a function of x B (integrating over z) in Figure 3. They are approximately the same at the high-energy limit. However, at the HERMES energy, they are not parallel to each other.
The experimental data from E665 [135] suggested a difference between D LL for the Λ 0 production and that for theΛ 0 production. This was later confirmed by the COMPASS [138] experiment. In [139][140][141], it was shown that such a difference serves as a flavor tag in the study of the G 1L FF. More studies on the flavor dependence of PDFs/FFs have been performed [122,[142][143][144][145][146][147][148][149]. The NOMAD collaboration also carried out similar measurements in the neutrino SIDIS experiment [150,151]. Because of the flavor-changing feature of the charged weak interaction, this experiment opens more opportunities for quantitative research on the flavor dependence of spin-dependent FFs. A sophisticated investigation was presented in [152].
RHIC is the first and, so far, the only polarized proton-proton collider. The helicity of the incident protons can be transferred to that of the partons through the longitudinal spin transfer g 1L (x) of PDFs. Therefore, it also has the capability of probing G 1L (z) of the fragmentation. The first measurement was performed in 2009 [153], while an improved analysis was presented in 2018 [154]. These experiments measure the spin transfer coefficient D LL which is defined as follows: The + symbol in the superscript denotes the helicity of the corresponding proton or Λ hyperon. The updated experimental data from the STAR collaboration [154] at RHIC are shown in Figure 4. This experimental data tend to favor the first and second scenarios in the DSV parameterization [118]. However, it cannot concretely rule out any scenario yet due to the large uncertainties. Moreover, the Xu-Liang-Sichtermann approach [131] based on the SU(6) spin-flavor symmetry can also describe this data well. Moreover, RHIC also measured the transverse spin transfer coefficient D TT , which is sensitive to the convolution of the transversity PDF and the transversity FF [155].  [154]. Data points are taken from [154]. The systematic error and the statistical error have been combined.
The polarizations of partons participating the same hard scattering are strongly correlated. The helicity amplitudes of different partonic processes have been evaluated and summarized in [156]. Thus, [157] proposes the dihadron polarization correlation as a probe to the longitudinal spin transfer G 1L in e + e − annihilations at low energy where the fragmenting quarks are not polarized. Recently, this idea was further investigated and applied to the unpolarized pp collisions in [158]. By measuring the longitudinal polarization correlation of two almost back-to-back hadrons, we also gain access to the longitudinal spin transfer in unpolarized pp collisions. Since this observable avoids the contamination from the longitudinal spin transfer g 1L in PDFs. which is also poorly known, [158] innovated a means to investigating the longitudinal spin transfer G 1L in FFs at RHIC, Tevatron, and the LHC. Furthermore, this work can also be used to constrain the FF of circularly polarized gluons.
Recently, the Belle collaboration measured the transverse polarization of Λ hyperons in e + e − annihilations [19], sparking considerable theoretical interest [15,[20][21][22][23][24][25][26][27][28][29][30][31][32]. In this experiment, one first defines the hadron production plane and then measures the transverse polarization along its normal direction. Since there are two transverse directions, we refer to the polarization along one as P N and the other one as P T . The hadron production plane can be defined in two ways. The first one is defined by the thrust axis and the Λ momentum. In the second, the thrust axis is replaced by the momentum of a reference hadron (in the back-to-back side). Therefore, this experiment is dedicated to probing the D ⊥ 1T (z, p ⊥ ) FF. While the p T -differential experimental data of P N contain sizable uncertainties, the p Tintegrated version is quite precise, as shown in Figure 5. Employing the Trento convention [159] for the definition of D ⊥ 1T (z, p ⊥ ), the p ⊥ integrated transverse polarization is given by [24,25,28,160] whereP ⊥Λ is the unit vector along the direction of P ⊥Λ . The integral in the denominator simply reduces to the product of two colinear FFs. However, to evaluate the numerator, we need to first parameterize the p ⊥ and p h⊥ dependence at the initial scale, which then evolves to the TMD factorization scale through use of the Collins-Soper-Sterman evolution equation. Nonetheless, since the collisional energy at Belle is not very high, a Gaussian ansatz is already a good approximation. More sophisticated approaches incorporating the p ⊥ dependence can be found in [29,31,32].      Figure 5. Transverse polarization of Λ in e + e − annihilation measured by the Belle collaboration [19]. Data points are taken from [19]. Statistical and systematic errors are combined in quadrature.
As shown in Figure 5, the distinct difference between P N measured in the Λ + π + (or K + ) and Λ + π − (or K − ) processes offers an opportunity to explore the flavor dependence. Early attempts, such as the D'Alesio-Murgia-Zaccheddu (DMZ) [24] and Callos-Kang-Terry (CKT) [25] parameterizations, adopted the strategy that valence parton FFs differ from each other and that parton FFs are the same, i.e., D ⊥Λ 1T,u ̸ = D ⊥Λ 1T,d ̸ = D ⊥Λ 1T,s ̸ = D ⊥Λ 1T,sea . However, this approach violates the isospin symmetry, which is one of the most important features of strong interaction. Furthermore, a model calculation [27] based on the strict SU(6) spin-flavor symmetry failed to describe the experimental data. However, it was first shown in [28] that the isospin symmetric Chen-Liang-Pan-Song-Wei (CLPSW) parameterization can still describe the experimental data well as long as the artificial constraint on sea parton FFs is released. This perspective was further investigated in Ref. [31] recently, which concluded that one can obtain good fit to the Belle data with and without implementing the isospin symmetry constraint after taking into account the charm contribution. This confirms that the current Belle dataset does not represent an isospin symmetry violation in the hadronization. Furthermore, [160] proposed to test the isospin symmetry at the future EIC experiment. By comparing the transverse polarizations in ep and eA scatterings at large x, we can ultimately check the difference between D ⊥ 1T,u and D ⊥ 1T,d . The future EIC is a polarized electron-proton/ion collider with unprecedentedly high luminosity. It will open a new window for the quantitative study of spin-dependent FFs. Several works [161][162][163][164] have proposed and made predictions for different observables at the future EIC with polarized proton beams. These observables are sensitive to various combinations of spin-dependent PDFs and FFs. Therefore, the future measurement will reveal information on both hadron structure and hadronization. A recent work [160] also proposed a method to study spin-dependent PDFs/FFs in unpolarized experiments.
The key idea is that the polarizations of the final state quark and initial state parton are correlated. Thanks to the Boer-Mulders function in the PDFs, the initial state quark are transversely polarized although the polarization depends on the azimuthal angle. This transverse polarization can further propagate into final state observables through chiral-odd FFs. By measuring the azimuthal-angle-dependent longitudinal and transverse polarizations of final state Λ, we can probe H ⊥ 1L and H ⊥ 1T even in the unpolarized SIDIS process. Moreover, we can also measure the azimuthal-angle-dependent polarizations in e + e − annihilations to probe combinations of the Collins function and spin-dependent chiral-odd FFs [160]. This idea is akin to those explored in [30,165].

Vector Mesons
Most vector mesons decay through parity-conserving strong interactions. Their polarization vector does not enter the angular distribution of the daughter hadrons. Therefore, it is not possible to measure their polarization vector. In contrast, the tensor polarization does play a role in the angular distribution and therefore can be measured. Among them, spin alignment, which quantifies the deviation from 1/3 of ρ 00 in the spin-density matrix, has received the most attention. Several collaborations [166][167][168][169] at LEP have measured the spin alignment of different vector mesons produced in the e + e − annihilation at the Z 0 -pole. We show the spin alignment of K * 0 and ρ 0 measured by the OPAL [166] and DELPHI [167] collaborations in Figure 6. The off-diagonal matrix elements were also measured in some of the experiments. Thereafter, the NOMAD collaboration measured the vector meson spin alignment for the first time in the neutrino DIS experiment [170]. These measurements offer more information on the hadronization mechanism and have led to several phenomenological studies [171][172][173][174][175][176][177][178][179]. Figure 6 shows that ρ 00 is consistent with 1/3 (i.e., no spin alignment) at the small-z region. However, at large z, a clear spin alignment is observed. This pattern is similar to that for the longitudinal polarization of Λ also measured at LEP [115,116] (shown in Figure 2).
As mentioned above, the quarks produced at LEP and also those at NOMAD are strongly polarized. Therefore, it is tempting to attribute the tensor polarization of final state vector mesons to the longitudinal polarization of the fragmenting quarks. However, a simple tensor structure analysis [47,63,64,175] shows that this is not the case. The spin alignment of the final state mesons is not coupled with the quark polarization. Instead, it is coupled with the quark-polarization-summed cross section. The vector meson spin alignment in e + e − collisions is given by where, ω q is defined to be the same as that for the Λ production in the previous section, and D 1LL (z) is the corresponding FF that is responsible for the vector meson spin alignment. As shown in the above equation, the longitudinal polarization of the fragmenting quark does not play a role here. It was thus first proposed in [63] that the vector meson spin alignment can also be observed in other high-energy collisions with unpolarized quarks fragmenting. Fitting to the experimental data from LEP, other work [108,109] extracted D 1LL (z) and made predictions for the spin alignment of high p T vector mesons in unpolarized pp collisions at RHIC and the LHC [109]. Furthermore, from the same mechanism, there will be a significant spin alignment for vector mesons produced in the unpolarized SIDIS. Measuring vector meson spin alignment at the future EIC will cast new light on the quantitative study of the D 1LL (z) FF. Notice that the spin alignment of low-p T vector mesons in AA collisions has also been measured at RHIC [110,180] and LHC [181] recently. These low-p T hadrons in relativistic heavy-ion collisions are produced through a different hadronization mechanism than those of fragmentation. Their tensor polarization originates from a different source.

Model Calculation
The PDFs and FFs are defined in terms of quark-gluon correlators as laid out in Section 3. Owing to the nonperturbative nature of the hadron state, we cannot directly evaluate them theoretically. Thus far, several proposals for computing quantities that can be related to the PDFs in the lattice QCD approach have been put forward [182][183][184][185][186]. However, it is not possible to study FFs in the lattice QCD yet. In the current stage, the quantitative information is mainly extracted from the experimental data.
However, due to the limited amount of experimental data, the TMD PDFs and FFs are not yet well constrained. As a complementary tool, model calculations have usually been employed to compute different PDFs over the past decades [106,. These investigations offer quantitative insight into the hadron structure and therefore are indispensable for phenomenological study. The same also goes for the FFs. Most of the models can be used to evaluate both PDFs and FFs. We make a nonexclusive brief summary on FF calculations.
There are quite a few models that can be categorized as a spectator model [220,221]. Among them, the quark-diquark model is a simple one which provides the quark-baryondiantiquark vertex,so that the baryon FFs can be easily evaluated. In [222], the colinear baryon FFs were calculated at the leading twist using the quark-diquark model, while in [189,[223][224][225][226], the TMD FFs were further computed at the leading and subleading twists. To compute the meson FFs or the gluon FFs, we need an improved version which offers the vertex among the fragmenting parton, hadron, and spectator. In [227], the Collins function was calculated for pions and kaons in this method. Recently, in [228], approach to calculate leading twist gluon TMD FFs was presented. The chiral invariant model [229] investigated the chiral symmetry and the spontaneous breaking with an effective Lagrangian of quarks, gluons, and goldstone bosons. It can also be classified into the spectator model category. Utilizing this model, several authors [230][231][232][233][234][235] calculated the pion and kaon FFs. In [236], an extended version was also developed to compute the vector meson FFs. Furthermore, several works [237][238][239][240] have evaluated the FFs of different hadrons using a parameterized quark-hadron coupling.
The Nambu-Jona-Lasinio (NJL) model originates from [241,242] who developed an effective theory describing the quark-hadron interaction. It has been employed to evaluate PDFs of different hadrons [191,[243][244][245][246][247][248]. Incorporating with the Feynman-Field model (also known as the quark-jet model) established in [249,250], the NJL=-jet model has been employed to calculate both colinear and TMD FFs of different hadrons [235,[251][252][253][254][255][256][257]. Recent works have also computed FFs of gluon [258] and charm quark [259] with this approach. The Feynman-Field model relates the total FF to the first rank FF. However, it does not specify how to compute the first rank FF. Therefore, in principle, it can be hybridized with another model which provides with the first rank FF to prolong the applicability of the corresponding model.
We make a final remark on the model calculation to conclude this section. All the above-mentioned models compute FFs employing the effective Lagrangian of partons and hadrons of interest. While these calculations offer quantitative insight into the hadronization scheme, we should draw a line between conclusions that are model-dependent and those that are model-independent.

Summary
There are a multitude of topics within the subject of FFs. In this review, we constrain ourselves in a very limited scope that we are familiar with. First, we briefly summarized the derivation of the TMD factorization and the establishment of the QCD evolution equation at the leading twist level. The TMD factorization and the corresponding evolution at the higher twist level are still ongoing topics. Second, we are particularly interested in the spin-related effects. With the spin degree of freedom being taken into account, the interplay between the transverse momentum and the hadron/quark polarization presents a highly intriguing phenomena that can be investigated in experiments. As a result, we need to define more TMD FFs to fully describe the fragmentation process. In quantum field theory, TMD FFs are introduced in the decomposition of parton correlators. We summarized the final results up to the twist-4 level for spin-0, spin-1/2, and spin-1 hadron productions. Finally, although all the TMD FFs have clear definitions in terms of parton fields and hadron states, they are nonperturbative quantities that cannot be directly evaluated from quantum field theory. In contrast to TMD PDFs, FFs cannot be computed even in the lattice QCD approach. The quantitative investigation thus mainly concentrates on the extraction from experimental measurements and model calculations. We summarized several spin-related experiments conducted over the past decades and the corresponding phenomenological studies. In the last section, we also briefly presented several model calculations.
The study of TMD FFs is still a very active field, and many mysteries remain to be explored. The Electron-Ion Collider (EIC) and the Electron-Ion Collider in China (EicC) have been proposed to be built as the new high-energy colliders in the next generation. They will provide new experimental data for the quantitative study of TMD FFs and can significantly boost our understanding of the hadronization mechanism.