Transverse Momentum Distributions of Heavy Hadrons and Polarized Heavy Quarks

We initiate the study of transverse momentum-dependent (TMD) fragmentation functions for heavy quarks, demonstrate their factorization in terms of novel nonperturbative matrix elements in heavy-quark effective theory (HQET), and prove new TMD sum rules that arise from heavy-quark spin symmetry. We discuss the phenomenology of heavy-quark TMD FFs at $B$ factories and find that the Collins effect, in contrast to claims in the literature, is not parametrically suppressed by the heavy-quark mass. We further calculate all TMD parton distribution functions for the production of heavy quarks from polarized gluons within the nucleon and use our results to demonstrate the potential of the future EIC to resolve TMD heavy-quark fragmentation in semi-inclusive DIS, complementing the planned EIC program to use heavy quarks as probes of gluon distributions.


Introduction
Hadronization -the nonperturbative mechanism that confines quarks and gluons produced in high-energy collisions into the experimentally observed color-singlet mesons and baryons -is a key aspect of virtually any process involving Quantum Chromodynamics (QCD), but its fundamental description from first principles remains elusive [1]. In the quest for this fundamental understanding, the fragmentation of bottom and charm quarks to heavy mesons can play a vital role because the mass of the heavy quark imprints as a perturbative scale on the otherwise nonperturbative dynamics of hadronization. The unique properties of heavy quarks as color-charged, but perturbatively accessible objects make them ideally suited as probes of the hadronization cascade, effectively serving as a static color source coupling to the light degrees of freedom. An improved field-theoretic understanding of heavy-quark fragmentation will also benefit the description of heavy flavor in Monte-Carlo generators for the LHC [2], where many key searches and Higgs coupling measurements involve final-state charm or bottom quarks. A rigorous field-theoretic framework in which hadronization can be studied in detail is that of transverse momentum-dependent (TMD) fragmentation functions (FFs), for which all-order factorization theorems have been established [3]. Like collinear fragmentation functions, TMD FFs depend on the longitudinal momentum fraction z H that the hadron retains from its parent quark. In addition, they describe the transverse momentum that the hadron picks up by recoiling against other fragmentation products, including the full quantum correlations with the quark polarization, which provides a three-dimensional picture of the fragmentation cascade. For processes with initial-state hadrons, TMD FFs are complemented by TMD parton distribution functions (PDFs) describing the three-dimensional motion of quarks and gluons inside the nucleon. The TMD dynamics of light quarks and gluons are a well-established field of experimental study [4][5][6][7][8][9][10][11][12][13][14], phenomenological analysis (see e.g. refs. [15][16][17][18]), and progress towards first-principle calculations using lattice field-theory [19][20][21][22][23][24][25]; for a recent comprehensive overview, see ref. [26]. Precision TMD measurements are a key physics target of the future Electron-Ion Collider (EIC) [27].
In this paper we study, for the first time, the TMD FFs of heavy quarks to heavy hadrons. Our theoretical tool to analyze the fragmentation of heavy quarks is (boosted) Heavy-Quark Effective Theory (bHQET) [28][29][30][31][32][33][34][35][36], which has previously been applied to the well-understood collinear (or longitudinal) heavy-quark FFs [37][38][39][40]. We demonstrate that applying bHQET to TMD FFs gives rise to novel, universal matrix elements describing the nonperturbative transverse dynamics of light QCD degrees of freedom in the presence of a heavy quark (i.e., a static color source). While a large part of this work is devoted to developing this new theoretical formalism, we will also consider the phenomenology of heavy-quark TMD FFs in two distinct processes, e + e − collisions and semi-inclusive deep inelastic scattering (SIDIS), which are illustrated in figure 1: (a) Heavy quarks are copiously produced in e + e − collisions. We are interested in the case where the quarks are produced relativistically, as is the case for charm quarks at existing B factories, such that their fragmentation processes are independent. The  Figure 1: (a) Heavy-quark pair production in the back-to-back limit in e + e − collisions, which gives access to heavy-quark TMD FFs. (Note that arrows indicate momentum flow, not fermion flow.) (b) Semi-inclusive deep inelastic scattering with an identified heavy hadron in the final state. This process gives access to individual heavy-quark TMD FFs convolved with a heavy-quark TMD PDF that can be perturbatively computed in terms of the collinear gluon PDF. (c) Heavy-quark pair production in electron-nucleon collisions, which we do not consider in this work. This process can probe the gluon TMD PDF, but is only sensitive to longitudinal fragmentation dynamics.
production rate is largest in the back-to-back region, where the cross section differential in the small hadron transverse momenta factorizes in terms of two unpolarized TMD FFs D 1 H/Q (z H , k T ). Furthermore, the heavy-quark pair is produced in an entangled transverse polarization state in central events. This entanglement imprints on the distribution of final-state hadrons as an azimuthal modulation known as the Collins effect, whose strength schematically is given by where H ⊥ 1 H/i is called the Collins function. The light-quark Collins effect has been measured in detail in pion and kaon samples by the Belle [41][42][43] and BaBar collaborations [44,45], but despite a proposal in ref. [44], a measurement (or search) of the heavy-quark Collins effect has not yet been performed. We will find that the heavy-quark Collins function encodes intricate nonperturbative physics, which motivates a dedicated measurement. 1 (b) A previously overlooked aspect of heavy-quark phenomenology in electron-nucleon collisions at the future EIC is that heavy quarks can be pair-produced in initial-state gluon splittings at small transverse momentum, with one quark e.g. going down the beam pipe and the other undergoing hard scattering and subsequent fragmentation into a heavy hadron, which in turn is reconstructed in a semi-inclusive measurement, as illustrated in figure 1 (b). Crucially, this process can be described by the standard TMD factorization for semi-inclusive deep inelastic scattering (SIDIS) when both the transverse momentum and the mass of the quark are small compared to the hard scattering energy Q. It thus retains the full sensitivity to the initial and final-state transverse momentum distribution (encoded in heavy-quark TMD PDFs and FFs) with respect to the photon direction.
To make phenomenological predictions for heavy-quark SIDIS, we fill a gap in the literature and compute all polarized heavy-quark TMD PDFs by perturbatively matching them onto collinear twist-2 nucleon PDFs, extending the analysis of the unpolarized heavy-quark TMD PDF in refs. [46,47]. 2 Interestingly, while many polarized TMD PDFs are strongly suppressed for heavy quarks, we find a nonzero leading result for the so-called worm-gear L function h ⊥ 1L , encoding the production of transversely polarized (heavy) quarks from linearly polarized light quarks and gluons, which is possible because the quark mass violates chirality. Indeed, the expectation that transverse quark polarization effects are suppressed by quark masses goes back to the early days of QCD [49], and we find that heavy-quark TMDs provide an arena to make this statement precise within twist-2 collinear factorization. For future phenomenology at the EIC, this nonzero conversion rate provides an exciting avenue to observe the heavy-quark Collins function because the factorized SIDIS cross section contains an azimuthal spin asymmetry, Importantly, this asymmetry unlike eq. (1.1) involves only a single Collins function multiplying the perturbatively predicted worm-gear L function, making it possible to extract its sign.
We point out that an extensive heavy-flavor physics program is already being planned for the Electron Ion Collider (EIC), which will leverage heavy-quark pair production as a hard probe of gluon TMDs [50][51][52][53][54][55] and of cold nuclear matter [56,57], as illustrated in figure 1 (c). We stress that this is not the case we consider in this paper: In case (c), the transverse momentum imbalance, production rate, and distribution of the heavy-quark pair are sensitive to the initial-state gluon TMD PDF (or the nuclear collinear gluon PDFs), but on the fragmentation side are at most sensitive to the well-understood longitudinal mass. We contrast these claims with our findings in section 4.1.1. 2 Ref. [47] also computed all secondary quark mass effects on light-quark distributions, including mass effects in the Collins-Soper kernel, which drive the renormalization of many of the objects we introduce here beyond leading-logarithmic order. The much more involved O(α 2 s ) secondary quark mass effects in the gluon TMD PDF were recently calculated in ref. [48]. momentum (z H ) distribution at leading power [53]. In contrast to this, the TMD processes in figure 1 (a) and (b) are directly sensitive to the nonperturbative transverse dynamics of heavy-quark fragmentation. 3 We note that the TMD fragmentation of light quarks to quarkonia has been studied in ref. [60], in that case by matching onto nonrelativistic QCD, and similarly for light-quark TMD dynamics in hard quarkonium production and decay in ref. [61,62].
The remainder of this paper is structured as follows: In Section 2, we analyze heavyquark TMD FFs and identify the new bHQET matrix elements and perturbative matching coefficients that characterize the fragmentation dynamics. In Section 3, we discuss the allorder structure of matching polarized heavy quark TMD PDFs onto collinear PDFs and explicitly compute the O(α s ) matching onto gluon PDFs. In Section 4, we use our results from the previous two sections to outline the prospects for heavy-quark TMD phenomenology at e + e − colliders and the future EIC.

Important note on conventions:
In an attempt to upset everybody equally, we will use QCD fields for writing down TMD correlators in this paper (thus hopefully making it accessible without background knowledge of soft-collinear effective theory), but will consistently make use of lightcone coordinate conventions as more commonly used in the SCET community (and also, for example, in bHQET). Specifically, we decompose fourvectors p µ in terms of lightlike vectors n µ ,n µ with n 2 =n 2 = 0 and n ·n = 2, such that e.g. p 2 = p − p + + p 2 ⊥ . As another part of this convention, collinear momenta near the mass shell will typically have large components p − ≫ p ⊥ ≫ p + . We always take transverse vectors with subscript ⊥ to be Minkowskian, p 2 ⊥ ≡ p ⊥ · p ⊥ < 0, and denote their magnitude by p T = −p 2 ⊥ . Finally, we define the metric and antisymmetric tensor in transverse space as Our convention for the antisymmetric tensor is ϵ 0123 = +1. 3 A very interesting middle ground is occupied by ref. [58], which analyzed the transverse momentum spectrum of heavy headrons inside groomed jets, also using bHQET. The resulting bHQET matrix elements in ref. [58] are predominantly longitudinal: Depending on the precise parametric regime, the dominant contribution to the transverse momentum either comes from perturbative collinear-soft modes stopping the soft-drop grooming algorithm [59], which factorize from the heavy-quark dynamics, or from bHQET modes subject to a primary soft-drop criterion and a secondary measurement on perturbative transverse momenta. (Compared to our results in this work for inclusive heavy-quark TMD FFs sensitive to the full transverse structure of the fragmentation cascade, this effect of the grooming is also responsible for the striking absence of Collins-Soper scaling reported in ref. [58].) Despite these differences in the experimental observable and theoretical structure, we find that it should be possible to establish a powerful connection between our results here and those of ref. [58]. Specifically, if O(Λ 2 QCD ) power corrections from intrinsic transverse hadronic dynamics can be resolved even within groomed jets, as outlined around their eq. (5.9), this would access a second moment of the unpolarized bHQET TMD fragmentation factor we will introduce in section 2.3. If this connection can be made precise, it would suggest that the bHQET fragmentation factors we consider in section 2.3 can appear also in scenarios where the mass has been integrated out together with additional observables (i.e., the soft-drop criterion in this case).
2 TMDs for heavy quark fragmentation into a heavy hadron

Calculational setup and parametric regimes
We consider the fragmentation of a (possibly polarized) heavy quark Q into a hadron H that contains the heavy quark and carries momentum P µ H . For this paper, we assume that the heavy hadron polarization is not experimentally reconstructed. We work in QCD with n f = n ℓ + 1 flavors, where the n ℓ massless quark flavors are denoted by q and the heavy quark Q has a pole mass m ≡ m c , m b ≫ Λ QCD . We decompose P µ H in terms of lightcone momenta as boosted in the frame of the hard scattering and by definition P H,⊥ = 0, coinciding with the "hadron frame" for fragmentation [3].
We are interested in the dependence of the fragmentation process on the total transverse momentum of additional hadronic radiation X into the final state, which is equal to the initial quark transverse momentum k ⊥ by momentum conservation, and Fourier conjugate to the transverse spacetime separation b ⊥ between quark fields. In position space, the TMD quark-quark correlator describing this fragmentation process is defined as where z H is the fraction of the quark's lightcone momentum retained by H, β, β ′ are the open spin indices of the quark fields, Tr denotes a trace over fundamental color indices, and b ≡ (0, b + , b ⊥ ). We have kept a sum over the possible hadron helicities h H , which are not experimentally resolved, implicit in the constrained sum over states, i.e., The Wilson line W (x) is defined as an anti-path ordered exponential of gauge fields extending to positive infinity along the lightcone directionn µ , For simplicity, we have suppressed the rapidity regulator, the soft factor, and transverse gauge links at infinity in eq. (2.2). The Wilson lines only depend on the direction ofn µ and are thus invariant undern µ → e αnµ . Taking P µ H andn µ to define n µ and tracking the α dependence through the definition of ∆ H/Q (z, b ⊥ ) implies that the "good components" of ∆ H/Q [26,63], by which we mean the components of the fermion fields that are picked out by the projector / n/ n/4 acting on ψ Q and that appear in leading-power factorization theorems, transform as ∆ H/Q → e −α ∆ H/Q under this relabeling. 4 In terms of the correlator in eq. (2.2), the bare unpolarized (D 1 H/Q ) and Collins fragmentation function (H ⊥(1) 1 H/Q ) from the introduction are defined in position space as 5 where tr denotes a trace over spin indices. The unpolarized TMD FF encodes the total rate for producing an unpolarized hadron from an unpolarized quark, while the Collins TMD FF describes the strength of the correlation between the quark's transverse polarization and the direction of the hadron transverse momentum. The leading TMD fragmentation functions have been proven to be universal between processes [68], i.e., they are independent of whether the Wilson line points to the future (e + e − → hadrons) or the past (SIDIS). Note that these scalar projections of the TMD fragmentation correlator are invariant under n µ → e αnµ by construction.
Since Λ QCD ≪ m, the nonperturbative dynamics in the fragmentation process are constrained by heavy-quark symmetry in all cases, but differences arise depending on the hierarchy between these two parametric scales and the magnitude k T of the transverse momentum or, equivalently, the inverse of the transverse distance 1/b T ∼ k T . Broadly speaking, we will consider the two cases illustrated in figure 2. In case (a), which we analyze in section 2.3, k T ∼ Λ QCD is generated during the nonperturbative fragmentation process itself, while perturbative emissions at the scale m ≫ k T are suppressed. In this case, the heavy hadron carries almost all the longitudinal momentum provided by the initial heavy quark, while the k T dependence is carried by universal nonperturbative functions describing how the "brown muck" separates from other light hadronic final states. In this regime, "disfavored" fragmentation functions where the valence content of the identified heavy hadron does not match the initial heavy quark, e.g. Q →H, Q → h, or q, g → H, are forbidden by heavy-quark symmetry. To simplify the analysis, we will count i.e., we will assume that the z H measurement does not probe the precise longitudinal momentum distribution near the endpoint, which is the case e.g. if z H is integrated over 4 In SCET this symmetry under the simultaneous relabeling n µ → e −α n µ andn µ → e αnµ is known as type-III reparameterization invariance [64,65] and is a manifest symmetry of the entire correlator because the bad components (/ n/ n/4) ψQ have been integrated out. 5 We use the same symbol for transverse momentum distributions in kT space and their Fourier transforms in bT space throughout this paper, as the meaning will always be clear from the context. Our conventions for Fourier transforms and the spin decomposition of TMD correlators follow ref. [66]. Note the superscript (1) on the bT -space Collins function indicating a bT derivative that arises from integrating a term / k ⊥ in the momentum-space correlator by parts, and that is specifically required due to the conventional normalization to the hadron mass [67]. For reference, the stated definition of the Collins function in position space is equivalent to a term tr γν ] appearing in the spin decomposition of the correlator.  : Parametric regimes for the fragmentation of a heavy quark Q (green) into a single heavy hadron H with a measurement on transverse momentum. In case (a), the hadron picks up transverse momentum relative to its parent quark as the "brown muck" (brown) coalesces around the quark and splits from other soft hadronic radiation into the final state (orange). In case (b), the transverse and longitudinal momentum distributions are dominated by perturbative emissions at the scale of the heavy quark mass (teal). In this regime, the nonperturbative hadronization process is encoded in a normalization factor, while its effect on the shape of the momentum distribution is subleading. a sufficiently large bin or when taking Mellin moments of the z H distribution. While relaxing this modification on the fragmentation side would be straightforward and leads to interesting fully-differential bHQET fragmentation shape functions, this would also require consistently modifying the description of the opposite collinear sector by reinstating the transverse momentum dependence in the formalism of ref. [40] for e + e − collisions, and by crossing the jet function in that reference into the initial state using the methods of ref. [69] for SIDIS, all of which we leave to future work.
In case (b), which we consider in section 2.4, the distributions in transverse and longitudinal momentum are determined by perturbative dynamics at the scale k T ∼ m ≫ Λ QCD , while the dynamics of the nonperturbative bound state only contribute a normalization factor. In the case of the unpolarized TMD FF, this normalization factor admits an interpretation as the total probability for Q to fragment into H, as is well known for inclusive heavy hadron production cross sections [37,38,70]. In this case, the disfavored fragmentation functions for Q →H or q → H, and Q → h or g → H, are perturbatively suppressed by O(α 2 s ) and O(α s ) at the scale µ ∼ m, respectively. We stress that we continue to assume eq. (2.6) also in this regime, so longitudinal momentum distributions remain perturbative.

Review of Boosted Heavy-Quark Effective Theory
The appropriate theory that describes the dynamics at the scale Λ QCD in either case and makes the heavy-quark symmetry manifest is boosted Heavy-Quark Effective Theory (bHQET) [35,36], i.e., the application of HQET [71] to heavy quarks produced in an energetic collision. The effective theory is constructed by integrating out the off-shell fluctuations of the heavy quark field at the scale m; these in particular include its antiquark component with energy gap 2m. The dynamic degrees of freedom are heavy-quark fields h v (x) that are labeled by the timelike direction v µ , which we choose to be the velocity of the heavy hadron, The tree-level matching of the massive QCD quark field onto h v at µ ∼ m reads: The h v (x) are implemented as Dirac spinors satisfying the projection relations For external states, the matching reads |H, h H ; X⟩ = √ m |H v , h H ; X⟩, and we use a nonrelativistic normalization convention for the bHQET states. In addition, the effective theory contains light-quark and gluon degrees of freedom that have isotropic momentum p µ ∼ Λ QCD in the rest frame of the heavy hadron. The tree-level matching for these is trivial; in particular, the Wilson line W (x) takes the same form in the effective theory, but consists of gluon fields that only have support on a restricted set of modes. For reference, the leading HQET Lagrangian is given by [71] where L light is a copy of the QCD Lagrangian with n ℓ massless quark flavors. The spin degrees of freedom of the heavy-quark can be explicitly decoupled from the light dynamics at leading power in 1/m by performing a field redefinition involving static Wilson lines Y v (x) [72,73], In this way, the heavy-quark Lagrangian becomes that of a free theory, in all external operators, acting as a static source of soft gluons. Specifically, the action of h v (x) on a product state in the decoupled theory is given by where s Q = 1 2 and h Q = ± 1 2 are the spin and helicity of the heavy quark, u(v, h Q ) = u(mv, h Q )/ √ m is an HQET spinor, and s ℓ , h ℓ , and f ℓ are the total angular momentum, helicity, and flavor content of the light degrees of freedom inside the hadron. (We will specify a helicity axis in the following section.) Note that the interpolating field for the light state on the right-hand side formally contains a future-pointing Wilson line to form an overall color singlet, which we suppress. Physical hadron states of definite angular momentum s H and helicity h H also have definite s ℓ , which is a good quantum number in the heavy-quark limit. In general, they involve a coherent sum over the helicity eigenstates in eq. (2.12), where we suppressed the common f ℓ and X, and ⟨s Q , h Q ; s ℓ , h ℓ |s H , h H ⟩ is a Clebsch-Gordan coefficient. (We take the coefficient to vanish for h Q + h ℓ ̸ = h H , i.e., one sum is always eliminated in practice by helicity conservation.) For the case of inclusive fragmentation, it has been known for a long time [38,71] that the factorized form of eq. (2.12) together with parity and eq. (2.13) implies relations between the fragmentation probabilities to different hadron states within the same heavy-quark spin symmetry multiplet, i.e., with the same s ℓ = 1 2 , 1, 3 2 , . . . . As an example, at the strict leading order in 1/m, an unpolarized charm quark is exactly three times as likely to fragment into an excited spin-1 vector meson (D * ) than into the corresponding pseudoscalar state (D). The physical reason for this is that the light dynamics do not see the heavy-quark spin at leading power, and thus the same nonperturbative matrix elements with given s ℓ , h ℓ appear in several cases. This analysis is simplified by the fact that for unpolarized or linearly polarized heavy quarks, light amplitudes for different helicities cannot interfere. One key goal of the next section will be to work out the consequences of heavy-quark spin symmetry for transverse momentumdependent fragmentation functions, where transverse quark polarization will let us access this interference for the first time.

Tree-level matching and discrete symmetries
Using the tree-level matching onto bHQET in eq. (2.8) at the leading order in Λ QCD /m ∼ k T /m, the correlator in eq. (2.2) evaluates to where F H is a bHQET bispinor defined by Note that we have evaluated the matrix element at b + = 0, which is justified at leading order in 1/m. This is easiest to see by using momentum conservation on the first matrix element to translate the fields back to b + = 0, where k − is the total residual minus momentum of the final state, and then expanding away k − ∼ v − Λ QCD compared to the leading O(P − H ) term P − H /z H − mv − in the Fourier phase in eq. (2.14). 6 Using The leftover boost factor 1/(n · v) = 1/v − ∼ m db + is a consequence of the nonrelativistic normalization of the external state, and ensures that projections of the above result onto good components are indeed invariant under the rescaling transformationn → e αnµ discussed below eq. (2.4).
To analyze the spin structure of F H (b ⊥ ), it is convenient to define the auxiliary vector which defines a unit z axis oriented along the spatial component ofn in the rest frame. Written out explicitly, F H (v, z, b ⊥ ) depends on the three four vectors v µ , i.e., the label direction in bHQET, corresponding to P µ H in the full theory, the spacelike vector z µ parameterizing the Wilson line directionn µ relative to v µ , and the spatial separation b µ ⊥ of the fields (with direction . As these three are orthogonal, they define a unique fourth unit direction y ρ = ϵ µνρσ v µ x ν z σ with y 2 = −1.
There are three applicable symmetries constraining the form of F H . First, the correlator populates only the particle-particle components of the bispinor, Second, from its definition, the correlator transforms under hermitian conjugation as =h v simplifies for the particle components. On the second identity we translated both matrix elements by −b ⊥ , exploiting that the resulting phase cancels between the states. Third, since parity is conserved in QCD (and thus in bHQET), the correlator satisfies where P denotes the unitary operator implementing parity in the rest frame of the heavy meson. Note that time reversal is broken by the presence of the out states, and thus is not a good symmetry of fragmentation functions [74].
in terms of two real-valued scalar coefficient functions χ 1,H (b T ) and χ ⊥ 1,H (b T ) that can only depend on v 2 = −z 2 = 1 and b 2 ⊥ = b 2 T . By performing the traces in eq. (2.5), we can identify these two functions with the unpolarized and Collins TMD FF, respectively, These results hold at the leading order in the heavy-quark expansion 7 and up to perturbative corrections at the scale µ ∼ m (which we reinstate in section 2.3.3), but capture the exact nonperturbative dependence on k T ∼ Λ QCD within the χ 1,H and χ ⊥ 1,H . For reference, we can also take suitable traces of eq. (2.22) to obtain explicit definitions of χ 1,H and χ ⊥ 1,H in terms of bHQET matrix elements, which we dub heavy-quark TMD fragmentation factors. Note that finding a nonzero result for the Collins fragmentation function during this step crucially relies on the presence of the Wilson line, which distinguished a nontrivial reference direction z µ in the rest frame. From our analysis so far, we can conclude that both leading-power TMD fragmentation functions for unpolarized hadrons are allowed by the discrete symmetries of QCD in the heavy-quark limit at leading order in 1/m. Written as in eq. (2.24), they are also manifestly independent of the flavor and mass of the heavy quark. However, the physical interpretation of χ 1,H (b T ) and χ ⊥ 1,H (b T ) is still fairly unclear at this point, and in fact we have not yet made use of heavy-quark spin symmetry. We address these questions in the next section, where we will derive an intuitive physical picture of the TMD fragmentation factors in terms of the individual constituents of the heavy hadron, and will derive powerful relations within spin symmetry multiplets.

Heavy-quark spin symmetry
We now return to the full correlator F H (b ⊥ ) defined in eq. (2.15) and analyze its heavyquark spin symmetry properties, which are particularly transparent when working with sterile fields. To do so, we first decompose the out states as in eq. (2.13), For definiteness, we take the magnetization axis defining the helicity eigenvalues to be the spatial component z µ of the Wilson line direction in eq. (2.18), which points back to the hard collision. Crucially, different helicities of the quark and the light degrees of freedom in the amplitude (h Q , h ℓ ) and the complex conjugate amplitude (h ′ Q , h ′ ℓ ) can interfere with each other, as only their sum is constrained to be equal to a common h H by helicity conservation. Note that we cannot use a completeness relation for the Clebsch-Gordan coefficients because they are only summed over the hadron helicity h H , but are not summed over the total hadron angular momentum s H because we assume that the experimental measurement can e.g. tell apart D and D * mesons. Acting on these out states with sterile heavy-quark fields as in eq. (2.12) yields On the last line we defined a generalized spin-density matrix ρ H,h Q h ′ Q (b ⊥ ) for the heavyquark helicities. Its entries are determined by the soft dynamics, the experimentally reconstructed values of s H , s ℓ (i.e., H), and angular momentum conservation.
To proceed, it is useful to express the outer product of spinors as where Σ µ Q is the heavy-quark spin operator acting on the nonrelativistic spin Hilbert space Evaluating the traces in eq. (2.24) then yields where y µ with y 2 = −1 is orthogonal to both b µ ⊥ and the Wilson line direction z µ , see the discussion below eq. (2.18). As expected from its relation to the unpolarized TMD FF, the Fourier transform of χ 1,H is simply the total conditional probability to produce H at rest given an initial quark momentum k T transverse to the direction of the static color source provided by the Wilson line. The bHQET analogue χ ⊥ 1,H of the Collins function on the other hand can be interpreted as a conditional density of quark spin Σ Q with respect to a magnetization axis defined by the static color source and the final-state transverse momentum. These physical interpretations roughly correspond to those of the associated relativistic fragmentation functions, which are often written in a form similar to eq. (2.28). Crucially, however, the meaning of the spin space on which the density matrix is defined is different in the heavy-quark case: For light quarks, whose spin after hadronization is an ill-defined concept, it only refers to the initial spin state in which the light quark is prepared. Heavy quarks at leading power in 1/m, by contrast, retain their initial spin state throughout the fragmentation process, and thus their spin density matrix together with the hadron spin measurement instead probes the angular momentum distribution of the final-state light constituents of the heavy hadron. 8 To make this fully explicit, let us introduce a shorthand for the spin-density matrix ρ ℓ of the light degrees of freedom ℓ ≡ {s ℓ , f ℓ }, The fact that the same light spin density matrix ρ ℓ appears for all hadrons within the same spin symmetry multiplet (same s ℓ and f ℓ , but different s H ) leads to relations between their TMD FFs in the heavy-quark limit. While it is interesting to ask how many independent nonperturbative functions the constraints on ρ ℓ from parity and hermiticity in eqs. (2.20) and (2.21) leave in principle, and how many of them are observable when reconstructing the hadron spin in addition, we now push on towards the combinations that are relevant for an unpolarized hadron and that contribute to the two fragmentation factors at hand.

Unpolarized TMD FF:
We begin with the unpolarized quark case and perform the trace in eq. (2.28), which sets h Q = h ′ Q and thus h ℓ = h ′ ℓ , To illustrate this, consider the pseudoscalar case, where and we have written helicities as ± ≡ ± 1 2 for short. We see that the unpolarized TMD FF encodes information about the magnitude of the amplitude for producing a given light 8 A related difference is that relativistic fragmentation functions have to be interpreted as number densities rather than probabilities due to the semi-inclusive measurement acting also on hadronic states at higher multiplicity. By contrast, the unpolarized bHQET TMD fragmentation factor upon Fourier transform, and up to renormalization effects [75], has an interpretation as a probability density because additional pair production of heavy quarks is power suppressed by kT ≪ m. helicity state. Summing over all hadrons H within the same spin symmetry multiplet M ℓ (i.e., all hadrons with identical light spin and flavor state ℓ), we further define where we used the completeness relation of the Clebsch-Gordan coefficients. Note that in the same way, this sum reduces the quark spin density matrix and the correlator to By evaluating the partial sums in eq. (2.30), it is easy to see that in terms of this baseline, the individual unpolarized TMD fragmentation factors are given by where for the purpose of this equation we used H as a shorthand for D (B) when Q = c (b), and similarly for the excited and higher-spin states. These relations are textbook knowledge in the inclusive fragmentation case (b T = 0) [71]. Our analysis shows, for the first time, that they hold without modification and point by point in the distribution when resolving the hadron transverse momentum. We anticipate that a transverse momentumdependent version of the Falk-Peskin parameter [38] appears when resolving individual (linear) hadron polarizations h H of hadrons with s ℓ = 3/2. Generalizations of the above to higher multiplets are trivial. Given that the sum over fragmentation factors within a multiplet encodes the full information on each individual one, it is interesting to ask what form the total fragmentation factor takes. Performing the complete sum over states, we have We see that the total unpolarized heavy-quark TMD FF reduces to a vacuum matrix element of two staple-shaped Wilson-line configurations along the lightlike and timelike direction, respectively. (Recall that we have suppressed transverse gauge links at infinity; the Wilson lines in the interpolating fields for the out states cancel.) This is reminiscent of the heavy quark-antiquark form factor proposed in ref. [20] to extract the TMD soft function on the lattice, but in contrast to that proposal directly relates to a physical observable, i.e., the total TMD cross section for producing heavy hadrons in e + e − collisions. We contend that this makes eq. (2.35) the theoretically simplest TMD observable possible in QCD, since it is entirely given in terms of vacuum matrix elements of Wilson loops. Collins TMD FF: A naive expectation from heavy-quark spin symmetry might be that the Collins FF should be suppressed by 1/m because it encodes a correlation between the initial quark transverse polarization vector and the transverse momentum of hadronic final-state radiation. In the case of light quarks, this correlation arises directly from the nonperturbative dynamics of the QCD Lagrangian, as illustrated in figure 3 (a), but in the heavy-quark case it naively seems to require a suppressed magnetic interaction with the heavy-quark spin. We will now see that this is not the case. As illustrated in figure 3 (b), the angle between the final-state heavy-quark and light transverse polarization vectors (i.e, the relative phase between their helicity states) determines which hadron in the spin symmetry multiplet is produced, even without a dynamical heavy-quark spin interaction taking place. Reconstructing this information experimentally thus induces a correlation between the heavy-quark and the light spin state. Crucially, spin symmetry ensures that the final-state heavy-quark spin state is identical to the one it was prepared in. The light spin state in turn is in general correlated with the transverse momentum k ⊥ of hadronic final-state radiation, since they both arise from the same nonperturbative dynamics of the light degrees of freedom, leading to a Collins effect at the leading order in 1/m. To illustrate this, it is again instructive to consider the case of the pseudoscalar meson, As expected, the Collins FF in the heavy-quark limit contains information about the strength of the interference, and hence the relative nonperturbative phases, of amplitudes for different light helicities.
As a corollary, we conclude that the Collins FF must vanish at leading order in 1/m when summing over all the hadrons in the spin symmetry multiplet, This is immediate to see from the diagonal form of the summed quark spin density matrix or the full correlator in eq. (2.33). Concretely, this means that the Collins FF vanishes altogether for s ℓ = 0 baryons, χ 1,Λ Q = 0. For the next few multiplets and using the same notation as in eq. (2.34), the explicit relations are Discussion: The spin symmetry relations in eqs. (2.34) and (2.38) are the main results of this section. They hold for all values of b T , which means that they also hold point by point in k T upon Fourier transform. Furthermore, they are unaffected by renormalization, as we discuss in the next section. This makes them substantially stronger than the known sum rules for relativistic TMD fragmentation functions. For the light-quark Collins function in particular, the Schäfer-Teryaev sum rule [76] has only been rigorously proven [77] in the bare case, and requires a sum over all possible hadrons, an integral over z h , and a weighted integral over k T . This is in contrast to the novel sum rule in eq. (2.37) that we have derived for the heavy-quark limit, which only requires summing the Collins function over a subset of hadrons and holds at any value of k T and z H . (We will explicitly extend these spin symmetry relations to k T ∼ m in section 2.4.) It thus implies the Schäfer-Teryaev sum rule for heavy quarks as a corollary. We caution that as in the inclusive case, the validity of the spin symmetry relations rests on the assumption that spin symmetry violation is negligible during the entire fragmentation process. This assumptions breaks down e.g. for those H, H * produced from the decays of H 1 , H * 2 whose spin symmetry-violating mass splitting is comparable to their widths [38].

All-order matching and renormalization
Using known results [78,79] for the perturbative matching of SCET onto bHQET, it is in fact straightforward to generalize eq. (2.23) to all orders in perturbation theory, Here the matching coefficient C m = 1+O(α s ) arises from separately matching the collinear ("unsubtracted") and soft contributions to the TMD FFs onto bHQET and QCD with n ℓ light flavors, respectively.
Importantly, the matching is diagonal in spin and forces z H = 1 because real radiation is parametrically forbidden at the scale µ ∼ m due to k T ≪ m, and thus we can separately match the collinear fields in the two matrix elements in eq. (2.2) onto bHQET. Starting at two loops, the matching coefficient features rapidity logarithms of the Collins-Soper scale ζ over the mass as a consequence of the large boost separating the heavy hadron rest frame and the frame where the soft radiation is isotropic. The appearance of the Collins-Soper scale can most transparently be understood by organizing the matching of the collinear sector onto bHQET in terms of gauge-invariant building blocks W † ψ Q with definite large lightcone momentum ω = √ ζ, as commonly done in SCET. The soft matching is nontrivial because starting at O(α 2 s ), vacuum polarization diagrams involving the heavy quark contribute to the expectation values of soft Wilson line operators in the n ℓ +1 theory. The two-loop result for C m was obtained in ref. [78], and our notation relates to theirs as Here the dependence on the rapidity scale ν cancels between the individual matching coefficients on the right-hand side, leaving behind the dependence on ζ/m 2 . The renormalization properties of χ 1,H and χ ⊥ 1,H follow from eq. (2.39) by consistency with the bHQET matching. Making the renormalization explicit requires introducing a rapidity regulator into the bHQET matrix element definitions, e.g. by modifying the lightlike Wilson lines in a standard fashion, and canceling rapidity divergences using the known TMD soft factor in a theory with n ℓ light quarks, which leaves behind an anomalous dependence on the boost factor √ ζ/m =n · v governed by the Collins-Soper kernel. (In the following we find it convenient to use the shorthand ρ =n · v for the third argument of the TMD fragmentation factors, which physically is given by the exponential of the hadron's rapidity in the frame of the hard collision.) This proceeds in close analogy to the standard relativistic case, so we leave an explicit check to a future perturbative calculation. Note that the Wilson lines in eq. (2.15) are in fact still lightlike up to possible off-lightcone regulators (and thus feature standard rapidity divergences) despite the presence of the mass because the opposite collinear sector is boosted close to the lightcone from the point of view of the bHQET rest frame. Conversely, the bHQET dynamics inside either collinear sector are also boosted from the point of view of the central soft modes. Since TMD renormalization is multiplicative in b T and independent of the hadronic state, it acts in the same way on all terms in the spin symmetry relations in eqs. (2.34) and (2.38), which implies that they also hold for the renormalized objects point by point in b T (or k T ).

Relation to bHQET fragmentation probabilities for Λ QCD ≪ k T
An important property of the TMD fragmentation factors we defined above is their limiting behavior as k T ≫ Λ QCD or, equivalently, b T → 0. In this limit, the unpolarized TMD fragmentation factor χ 1,H is related to the total probability χ H for the quark to fragment into H, which has previously been analyzed in HQET [37][38][39], where the matrix-element definition of χ H [40] is equal to χ 1,H (b T = 0) at the bare level, Because χ H is not renormalized [70], we generally expect a perturbative Wilson coefficient T at leading order in the strong coupling, 9 and is analogous to the behavior of the usual TMD PDF as b T → 0 where it approaches its total momentum-space integral given by the collinear PDF, up to renormalization effects and radiative corrections [75]. Here we also assumed without detailed proof that corrections to this relation are quadratic in b T based on the azimuthal symmetry of χ 1,H (b T ). From H χ H = 1, it follows that the total TMD fragmentation factor χ 1 (b T ) defined in eq. (2.35) is purely perturbative in this regime, In contrast to eq. (2.41), the Collins TMD fragmentation factor χ ⊥ 1,H must vanish at least linearly as b T → 0 because there is no leading bHQET matrix element it could match onto in this limit. This is easiest to see by repeating the symmetry analysis of the bHQET correlator F H in section 2.3.1 at b T = 0, which only admits F H (v, z, b ⊥ = 0) = χ H (1 + / v)/2. As we will see by comparing to the limit Λ QCD ≪ k T ≲ m in section 2.5, the expansion indeed starts at the linear order, and the matrix-element definition of the relevant nonperturbative parameter at O(Λ QCD b T ), as well as its tree-level Wilson coefficient, can all be inferred from consistency.

Matching TMD FFs onto bHQET for Λ QCD ≪ m ≲ k T
We next consider case (b) in figure 2. In this regime, the transverse and longitudinal momentum distributions are determined by dynamics at the scale µ ∼ m ∼ k T and are fully perturbative. The nonperturbative dynamics in this case are encoded in bHQET matrix elements that involve additional gluon fields or derivatives and that can be nonlocal along the lightcone, but in contrast to the previous section are local in the transverse direction. Similar to a standard twist expansion, these bHQET matrix elements are categorized by their mass dimension, which determines their scaling as Λ QCD ≪ m, k T , i.e., their mass dimension ∼ Λ n QCD is compensated by powers of b T or 1/m in the Wilson coefficient. This story plays out differently for the unpolarized vs. the Collins TMD FF, which scale as O(1) and O(Λ QCD b T ), respectively, so we will go through the two cases separately in the 9 It is well known that the formal OPE of relativistic fragmentation functions is ambiguous due to an unconstrained choice of boundary condition at lightcone infinity [80,81]. While this fundamental issue remains present here, it is interesting to ask whether the case of bHQET TMD fragmentation factors, which are Wilson loops, can provide additional insight into this issue.
following. We note that the expansion of TMD FFs in terms of bHQET operators differs from a standard twist expansion insofar as the HQET field h v encoding the interactions with the heavy valence quark remains present in all low-energy matrix elements.

Unpolarized TMD FF
For the unpolarized TMD FF, the unique bHQET matrix element that can arise in the infrared at the leading order in 1/m is the total fragmentation probability χ H as defined in eq. (2.42), which follows from symmetry arguments similar to those below eq. (2.43): Importantly, we have again made use of the assumption in eq. (2.6) that we are sufficiently far away from (or have fully integrated over) the endpoint regime z H → 1, as otherwise there would be a nontrivial bHQET shape function on the right-hand side [39,40]. The unique matching coefficient of χ H , which we dub the partonic heavy-quark TMD FF d 1 Q/Q (z, b T , µ, ζ), is a new object that, to our knowledge, appears in our analysis for the first time. 10 It is independent of the precise hadronic final state, carries the exact dependence on b T m ∼ 1, and can be calculated perturbatively by evaluating eqs. (2.2) and (2.5) for partonic final states including at least one heavy quark, i.e., Its rapidity renormalization is governed by the Collins-Soper kernel of a theory with n ℓ massless and one massive quark degree of freedom [47]. We leave a dedicated NLO calculation of d 1 Q/Q (z, b T , µ, ζ) to future work. Since the dependence on the hadronic final state is purely encoded in χ H , which satisfies the same spin symmetry relations as in eq. (2.34), we conclude that the unpolarized heavy-quark TMD FF satisfies for all values of b T (or k T ), including 1/b T ≳ m, up to corrections of O(Λ QCD /m). Eq. (2.44) continues to be valid for k T ≫ m, but features large perturbative logarithms of b T m ≪ 1 in this limit. Their resummation is enabled by further factorizing the physics at those two scales. To do so, we can first match the heavy-quark TMD FF onto twist-2 heavy-quark collinear FFs at the scale µ ∼ m [46],

47)
10 Curiously, the perturbative transverse dynamics of heavy-quark fragmentationb → Bc have previously been evaluated in refs. [82,83]. The complete tree-level result given in the first reference, which starts at O(α 2 s ), can be considered a very specific subset of the NNLO corrections to the TMD FF we define here if we sum over final states. If we tag on the charm instead, their result corresponds to a different perturbative TMD FF d 1 bc/b ∼ α 2 s whose renormalization, by our analysis, is governed by standard (massive) TMD evolution.
where the sum runs over i = q,q, g. This matching takes the same form as the standard matching of light-quark TMD FFs onto twist-2 FFs at µ ∼ Λ QCD , except that the highest IR scale here is given by m. The Wilson coefficients J i/q (z, b T , µ, ζ) encode the perturbative process q → i in a theory with n ℓ + 1 massless flavors at the scale µ ∼ k T , with the quark retaining a fraction z of the parent's lightcone momentum, and are known to N 3 LO [84,85]. In a second step, we perform the well-known [37][38][39][40] matching of the collinear FF of a massive quark onto bHQET to separate Λ QCD ≪ m, This refactorization condition for d 1 Q/Q can serve as a cross check on future perturbative calculations, and in addition enables resumming logarithms of k T /m ≫ 1.

Collins TMD FF
To identify the low-energy bHQET matrix element that the Collins TMD FF matches onto in the limit Λ QCD ≪ m ∼ k T , we use a two-step matching procedure formally valid for the hierarchy Λ QCD ≪ m ≪ k T . (We will later show that the result is correct for either hierarchy.) As for the unpolarized TMD FF above, this lets us make use of well-known results for the matching of light-quark TMD FFs onto collinear FFs, which we can then further match onto bHQET. We start from the diagrammatic small-b T expansion of the bare Collins TMD FF for light quarks, which is valid for Λ QCD ≪ k T and given by [74,88,89] whereĤ h/q is a twist-3 collinear fragmentation matrix element at the scale µ ∼ Λ QCD , 51) 11 In the literature, this relation is more commonly given as a tree-level equality betweenĤ h/q and a weighted kT integral over the bare momentum-space Collins FF. Using eq. (2.65) and integrating by parts, it is easy to see that this reduces to the derivative of the bT -space Collins FF at bT = 0. The O(αs) corrections to eq. (2.50) were evaluated at finite kT > 0 in ref. [89] and involve twist-3 matrix elements that depend on two independent momentum fractions and reduce toĤ h/q in certain limits by use of the equation of motion. We anticipate that matching these more general matrix elements onto bHQET will reduce the number of independent (residual) momenta to one because the heavy-quark momentum is fixed.
where σ µ− = i 2 [γ µ , / n] and we have defined as a shorthand for the insertion of a gluon field strength tensor G µν anywhere along the lightcone, with W (x, y) a straight Wilson-line segment connecting x and y.
We now consider the heavy-quark Collins FF and at first assume the hierarchy Λ QCD ≪ m ≪ k T . For the matching at the scale µ ∼ k T , the mass is an infrared scale, and thus the twist expansion in eq. (2.50) immediately carries over. The collinear matrix element H H/Q takes the same form as eq. (2.51), but is now defined at the scale µ ∼ m. To implement the separation of scales Λ QCD ≪ m, we matchĤ H/Q onto bHQET. At tree level, this amounts to a replacement of the quark fields as in eq. (2.8), and after expanding the momentum-conserving phase results in where χ H,G ∼ Λ QCD is a novel subleading bHQET matrix element defined by Similar to the total fragmentation probability χ H defined in eq. (2.42), χ H,G no longer depends on b ⊥ , but is simply a constant that depends on the identified hadron H. Note that a nonzero value of χ H,G is compatible with all the symmetries of bHQET: Its defining spin correlator X(v, z) (dropping the spin trace) is hermitian by construction, satisfies P + X = XP + = X, and under parity transforms as PX(v, z)P = X(v, −z). Repeating the analysis in section 2.3.1, it is therefore proportional to 1 + / v, which has nonzero trace. In the last step, we combine eqs. (2.50) and (2.53) to arrive at our final result for the tree-level matching of the heavy-quark Collins TMD FF onto bHQET: Because this derivation assumed Λ QCD ≪ m ≪ k T , eq. (2.55) a priori is only valid up to power corrections in mb T . However, since we found a nonzero result at our tree-level working order and power corrections in mb T can only arise from real radiation in the calculation of the Wilson coefficient, eq. (2.55) as written also holds when integrating out both scales simultaneously. We note that additional low-energy matrix elements will in general be generated when performing the matching at higher orders in α s , but leave a dedicated construction of the basis of bHQET operators at this order in Λ QCD to future work. We point out that an observation of the heavy-quark Collins function in this regime would provide interesting insight into novel gluon correlations in the heavy-quark fragmentation process that are encoded in χ H,G . More specifically, χ H,G encodes a correlation between the gluon field polarization and the transverse polarization of the light constituents of the heavy hadron in the final state, which as in section 2.3.2 is indirectly resolved by reconstructing the total hadron spin, e.g. by distinguishing D and D * mesons. Conversely, χ H,G must vanish when summing over all hadrons in the spin symmetry multiplet M ℓ , This result is straightforward to prove along the lines of section 2.3.2 by decoupling the heavy quark fields in eq. (2.54) and exploiting the completeness relation of the Clebsch-Gordan coefficients, which leaves a trace of the form tr σ βα (1 + / v)/2 = 0. Combining these results at large k T ∼ m with those in eq. (2.38) we conclude that the Collins TMD FF satisfies the following relations for all values of b T (or k T ), which we have proven here up to corrections of O(Λ QCD /m) and up to radiative corrections at the scale µ ∼ k T ∼ m for large k T . We conjecture that the additional bHQET matrix elements generated by the matching at higher orders in α s will involve the same Dirac structure as eq. (2.54), i.e., an additional insertion of γ µ ⊥ , and thus will also satisfy eq. (2.56), but leave a detailed all-order analysis in this regime to future work. (We recall that for k T ≪ m eq. (2.57) holds to all orders in the strong coupling, see section 2.3.2.)

Consistency between regimes for Λ QCD ≪ k T ≪ m
Our results in the previous two sections share a common domain of validity when the transverse dynamics are already perturbative, Λ QCD ≪ k T , but still subject to heavyquark symmetry, k T ≪ m. In this section we analyze the consistency relations that arise from this overlap and relate the perturbative bHQET fragmentation factors to the partonic heavy-quark TMD FFs.
We start with the unpolarized case. Comparing eqs. (2.39) and (2.41), which are valid for Λ QCD ≲ k T , to eq. (2.44), valid for k T ≲ m, we find the following all-order refactorization relation for the partonic heavy-quark TMD FF in the limit k T ≪ m, Here we have canceled off the common nonperturbative factors of χ H . To interpret the z dependence, eq. (2.58) says that counting 1 − z ∼ 1, d 1 Q/Q must approach δ(1 − z) up to an overall factor at the distributional level for b T m → ∞, i.e., all Mellin moments of d 1 Q/Q must become equal in this limit. We expect that eq. (2.58) will provide a powerful consistency check of future perturbative calculations of d 1 Q/Q . It also enables the resummation of large perturbative logarithms of k T /m ≪ 1, complementing the factorized result in eq. (2.49) for the opposite limit. For the Collins TMD FF we compare eq. (2.39) to eq. (2.55) and use C m = 1 + O(α s ). Canceling off the z dependence, which is trivial at tree level, this yields which can be interpreted as the leading linear term in a small-b T expansion of χ ⊥ 1,H , as anticipated in section 2.3.4. As for the Collins function at k T ∼ m, we leave a dedicated higher-order matching calculation to future work, which will involve nontrivial Wilson coefficients integrated against at least one additional O(Λ QCD ) bHQET matrix element.

Model functions and numerical results
For our numerical results we will assume a simple Gaussian model for the unpolarized TMD fragmentation factor, where κ H ∼ Λ QCD has units of GeV. Eq. (2.60) is valid at initial scales µ 0 ∼ m ρ 0 ∼ 1/b T of the TMD evolution and satisfies eq. (2.41) up to corrections in α s (µ 0 ). To be specific, we apply a µ * prescription [90,91] (also known as a "local" b * prescription) starting at O(b 4 T ) to ensure that µ 0 stays perturbative without polluting nonperturbative corrections at O(Λ 2 QCD b 2 T ) [75], where b 0 = 2e −γ E ≈ 1.12292 and we take µ min = 1 GeV. We take ζ 0 ≡ mρ 0 to always be equal to its canonical value, ζ 0 = (b 0 /b T ) 2 . We then use leading-logarithmic (LL) perturbative TMD evolution U q (µ 0 , ζ 0 , µ, ζ) to evolve eq. (2.60) to the overall scales µ ∼ √ ζ ∼ Q, with Q the hard scattering energy. 12 This order is sufficient for the exploratory phenomenology we have in mind, and in particular lets us use TMD evolution and β functions in QCD with n f = 5 massless flavors at all scales since the quark decoupling only induces next-to-leading logarithms of b T m. Specifically, we ignore the decoupling relations and NNLL power-like secondary quark mass corrections to the Collins-Soper kernel γ q ζ (b T , µ) that were determined in ref. [47]. We also ignore nonperturbative contributions to the Collins-Soper kernel, since they are orthogonal to the effects we are interested in here. Overall, this results in the following expression for the evolved unpolarized heavy-quark TMD FF, where for definiteness we considered the integral over z cut ≤ z H ≤ 1. To our working order, the right-hand side of eq. (2.62) is independent of z cut as long as 1 − z cut ∼ 1 in order to satisfy eq. (2.6), and also holds for any truncated z H moment of the TMD FF. Note that the single-parameter model in eq. (2.62) is also accurate at large k T ≳ m, cf. eq. (2.44), where it reduces to χ H and thus is correct up to radiative corrections. We assume a similar model for the Collins TMD fragmentation factor, but have to account for the suppression at small b T by modifying the Gaussian, see eq. (2.59), where we find it convenient to express the overall effect strength in terms of λ H⊥ = χ H,G /χ H ∼ Λ QCD , i.e., relative to the total fragmentation probability χ H . The parameter κ H⊥ ∼ Λ QCD controls the relative impact of higher power corrections and is in general distinct from κ H in eq. (2.60). Combining this with NLL n f = 5 TMD evolution as above, we find, for the evolved heavy-quark Collins function in position space, Taking appropriate Bessel integrals [67], we finally transition to momentum space, To evaluate the TMD evolution and the Bessel integrals, we use the numerical implementation of TMD anomalous dimensions, QCD renormalization-group solutions, and doubleexponential oscillatory integration in SCETlib [92].
Our results for the z H -integrated heavy-quark TMD FFs are shown as a function of k T for different values of the model parameters in figure 4. We use α s (m Z ) = 0.118 GeV as the input value for the strong coupling. We note that due to heavy quark flavor symmetry, the charm and bottom-quark TMD FFs are exactly equal at small k T ≪ m. In other words, they only depend on the universal Gaussian parameters κ H (for the unpolarized TMD FF), κ H⊥ (for the Collins TMD FF), and the Collins effect strength λ H⊥ . At large k T ∼ m, the TMD FFs remain independent of the heavy quark mass up to radiative corrections of O(α s ), which we ignore at our LL working order. These plots are thus identical for both flavors we consider. We point out that the Collins function can in general take any sign, as indicated by the yellow band scanning various values of the effect strength λ H⊥ . The effect of varying the size of higher-power corrections (κ H , κ H⊥ ) decreases as k T increases for both TMD FFs, as expected.

Calculational setup
In this section we consider the production of a heavy quark Q with pole mass m = m c , m b ≫ Λ QCD from light partons within a polarized nucleon N . The nucleon has momentum with P − N ≫ P + N = M 2 N /P − N in the rest frame of the hard scattering that the heavy quark participates in. Note that we again take the large component of the hadron (nucleon) momentum to be along the n µ direction to make this section self contained, but the case of ann-collinear incoming hadron (as would be consistent with the n-collinear outgoing hadron we considered in section 2) follows from n µ ↔n µ .
This time, we are interested in the transverse momentum k ⊥ of the heavy quark with respect to the nucleon beam axis, which is again Fourier conjugate to the transverse spacetime separation b ⊥ between quark fields. The bare TMD quark-quark correlator between forward nucleon states that describes this process is where b ≡ (0, b + , b ⊥ ), W was defined in eq. (2.4), x is the lightcone momentum fraction carried by the heavy quark, and we have suppressed the rapidity regulator, the soft factor, and transverse gauge links at infinity for simplicity. For the explicit perturbative calculations in this section, it will also be useful to define the momentum-space version of the above correlator, The spin decomposition of eq. (3.3) in terms of scalar TMD PDFs is well known [66,93], where M N is the nucleon mass, S L is the longitudinal nucleon polarization in the Trento frame [94], and in our convention Φ Q/N (x < 0, k ⊥ ) decomposes in the same way in terms of the antiquark TMD PDFs f 1Q/N (|x|, k T ), etc. We have suppressed terms ("bad components") that do not contribute to leading-power TMD factorization theorems. As we will see in the next section, the terms proportional to the transverse nucleon polarization S ⊥ vanish for heavy quarks to all orders in the strong coupling when matched onto the leading (twist-2) collinear PDFs. We will also find that the twist-2 matching for the Boer-Mulders function h ⊥ 1 vanishes at O(α s ). The remaining TMD PDFs on the first line, for which we will find nonzero results at O(α s ), are the unpolarized TMD PDF f 1 , the helicity TMD PDF g 1L , and the so-called worm-gear L function h ⊥ 1L ; the latter will be of particular significance, and encodes the production of a transversely polarized quark from a linearly polarized nucleon. For reference, the explicit Hankel transforms relating scalar TMDs in b T and k T space read 13

Matching onto twist-2 collinear PDFs
Heavy-quark TMD PDFs are different from their TMD FF counterparts because the heavy quark cannot be part of the initial-state nucleon wave function at the scale µ ∼ Λ QCD at leading power in Λ QCD /m, 14 whereas in the fragmentation case the heavy quark is always part of the final-state heavy hadron until its eventual weak decay. This means that heavy quarks must be pair-produced in initial-state gluon splittings at the scale µ ∼ m instead. In particular, this means there is at least one perturbative emission with transverse momentum ≳ m setting the scale of k T ≳ m, while the region of k T ≪ m can only be populated by several emissions with small net recoil, which is a power-suppressed configuration. In field theory terms, this means that heavy-quark TMD PDFs can be computed by perturbatively matching them onto collinear twist-2 nucleon PDFs in a theory with n ℓ light flavors, which are the only nonperturbative piece of information in this case. The matching onto twist-2 collinear PDFs is well developed for light quark and gluon TMDs, with notable results including all unpolarized quark matching coefficients through O(α 3 s ) [98,99] and results for 13 We continue to distinguish momentum and position-space functions by their argument, see footnote 5.
For the meaning of the superscript (1), see also there. 14 Power corrections of this kind, which are known as "intrinsic charm" and have received substantial recent interest on the collinear PDF side [95,96], would be an interesting subject to explore in the TMD case in the future. Very recently, the TMD PDFs for charm quarks within Λc baryons, which are leading valence contributions and do not have to be produced from gluons, have been evaluated in a lightfront Hamiltonian model in ref. [97]; while these are phenomenologically inaccessible, it would be interesting to analyze these valence dynamics in the heavy-quark limit as we did for TMD FFs in section 2.
polarized TMDs through O(α 2 s ) [100,101], and many of the following steps are standard, see e.g. [3]. Likewise, the O(α s ) matching of the unpolarized heavy-quark TMD PDF onto gluon collinear PDFs has been given in refs. [46,47]. We nevertheless aim for a selfcontained description, giving us the opportunity to point out the ways in which (polarized) heavy-quark TMD PDFs behave differently.
The bare light-quark and gluon twist-2 collinear correlators are defined as where b ≡ (0, b + , 0) in this case and W (b, 0) denotes a straight Wilson line segment. The collinear correlators are conventionally decomposed as [3] in terms of the unpolarized (helicity) quark and gluon PDFs f i/N (g i/N ) and the transversity quark PDF h q/N . The contribution ∝ S ⊥ to the gluon correlator (i.e., the transversity gluon PDF) vanishes identically for spin-0 and spin-1/2 hadrons in the initial state due to helicity conservation [102]. The matching relation between heavy-quark TMD PDFs and twist-2 collinear PDFs holds at the operator level, and constitutes the leading term in the OPE of the former. Taking nucleon matrix elements of the bare operators, the relation for general spin indices reads where p − is the lightcone momentum carried by the light parton extracted from the collinear PDF and the sum runs over the n ℓ light quark flavors. In pure dimensional regularization, the bare matching coefficients are given by the partonic diagrams where z = xP − N /p − is the fraction of p − injected into the hard scattering process and we have indicated the heavy quark lines in red. The gray-shaded circles denote the sum of all possible QCD diagrams with these external legs, including gluon attachments to the Wilson lines that are part of the operators indicated by ⊗. We have included the respective lowest-order diagram for illustration. As is standard, matching relations between individual scalar TMD and collinear PDFs follow by inserting eq. (3.7) into eq. (3.8) and tracing the resulting Dirac bispinors (. . . ) ββ ′ against the relevant Dirac structures.
Flavor conservation in QCD implies that a single fermion line has to connect the external light-quark states in eq. (3.9). It follows that contractions with the quark transversity PDF involve an odd number of Dirac matrices on the light-quark line and vanish to all orders, i.e., flavor conservation and chirality for light quark flavors imply that all terms ∝ S ⊥ vanish at twist-2 level in eq. (3.4). This is distinct from e.g. the light-quark transversity TMD PDF, which receives a tree-level contribution from the transversity collinear PDF of the same flavor. As in the case of light-quark TMD PDFs, Lorentz covariance further implies that only unpolarized (helicity) collinear PDFs can contribute to the unpolarized and Boer-Mulders (helicity and worm-gear L) TMD PDFs, matching the dependence on S L in the spin decomposition. These conclusions are not modified by the inclusion of the soft factor, the rapidity renormalization, and the UV renormalization of the TMD PDFs, all of which are orthogonal to the spin structure. They are likewise unaffected by the renormalization of the collinear PDFs, which acts autonomously on the unpolarized and longitudinally polarized sectors. Passing to renormalized objects, this altogether leaves us with the following four nontrivial matching relations for heavy-quark TMD PDFs onto collinear PDFs, Here the subscripts λ, λ ′ = ∅, ∥, ⊥ on C Q λ /j λ ′ (z, k T , µ, ζ) label the polarization of the heavy quark and the light parton j, the sum runs over gluons and the n ℓ flavors of light quarks and antiquarks, and we have included a factor of k T /M N on the left-hand side as needed to ensure that the matching coefficient is independent of the hadronic state. We have also changed integration variables from p − in eq. (3.10) to z, exploiting the fact that projections of the matching coefficients onto good components can only depend on z by reparameterization invariance. Note that in a crucial difference to the light-quark case, the heavy-quark worm-gear L TMD PDF, which involves an odd number of Dirac matrices on the heavy-quark line in eq. (3.9) is allowed at twist-2 level because the quark mass breaks chirality. The same is true for the Boer-Mulders function. In both cases, the original argument of ref. [103] why the twist-2 matching for these functions vanishes to all orders in the light-quark case critically relied on chirality. 15 Conversely, the respective matching coefficients must vanish linearly as m → 0 to afford the helicity flip, Lastly, note that to all orders it is only the gluon PDF f g (x) and the quark singlet PDF i=q,q f i (x) that contribute to the sum f 1 Q/N + f 1Q/N due to the invariance of eq. (3.9) under the n ℓ light flavor symmetry, and similarly for the two polarized cases. The difference f 1 Q/N − f 1Q/N of heavy quark and antiquark TMD PDFs receives a nonzero contribution s ) due to the relative orientation of the color flow along the fermion lines in eq. (3.9), as in the light-quark case [84,85,98].
Inverting the Hankel transforms in eq. (3.5), we find the b T -space matching relations where the matching coefficients are given by (n = 1 for λ =⊥ and n = 0 otherwise) For the dimensionless b T -space matching coefficients, eq. (3.11) simply reads (3.14) 15 An earlier version of this manuscript incorrectly stated that the twist-2 matching for the heavy-quark Boer-Mulders function should vanish to all orders based on its transformation behavior under time reversal, which however only constrains the leading O(αs) diagram we consider in the next section. We thank Markus Diehl for pointing this out to us.

One-loop evaluation of matching coefficients
At O(α s ), only the gluon diagram in eq. (3.9) is nonzero. Using standard QCD Feynman rules, we find the leading-order result where p = (p − , 0, 0) is the momentum of the external gluon and ℓ is defined as indicated (in the direction of fermion flow). The ℓ + integral is straightforward to evaluate by contours, which amounts to setting an emitted antiquark on shell, ℓ 2 = m 2 . Note that the diagram is finite in four dimensions and without a rapidity regulator because the quark mass cuts off infrared singularities. This is expected, as the UV and rapidity renormalization only become nontrivial at the next order. Dotting eq. (3.15) into the gluon collinear PDF correlator in eq. (3.7) and projecting onto quark spin structures, we find individual momentum-space matching coefficients with leading-order coefficient functions given by (3.17) As a nontrivial check, we have confirmed that using massive SCET Feynman rules [104] results in the same expressions after performing the spin traces and integrating over the loop momentum. Note that the projection of the O(α s ) twist-2 matching diagram onto the Boer-Mulders function remains zero even for finite quark masses. This is expected because the Boer-Mulders function is odd under time reversal [105], i.e., it changes sign depending on whether the Wilson lines in the operator point to the future (SIDIS) or the past (Drell-Yan). The diagram in eq. (3.15) does not yet feature gluon attachments to the Wilson lines that could resolve their direction, and thus its projection onto the Boer-Mulders function has to vanish. Starting at O(α 2 s ), the matching coefficient can in general receive nonzero contributions from the absorptive part of real-virtual diagrams because chirality is broken by the quark mass, and it would be interesting to investigate these contributions further.
Evaluating the inverse Hankel transforms in eq. (3.13), we find the position space matching coefficients which at this order only depend on the dimensionless combination b T m and are given by where K 0 and K 1 are modified Bessel functions of the second kind. These are the main analytic results of this section. The unpolarized matching coefficient C Q/g has been computed long ago [46], and we agree with the b T -space expression given in that reference as well as with the k T -space result in ref. [47]. The results for the polarization-dependent matching coefficients are new.

Consistency with the light-quark limit
For Λ QCD ≪ m ≪ k T , heavy-quark TMD PDFs can be determined using a two-step matching [47]. First, the TMD operators at the scale µ ∼ k T are matched onto collinear PDFs in a theory with n ℓ + 1 massless quark flavors, which results in the standard massless TMD matching coefficients. In a second step, the n ℓ +1-flavor PDFs are matched onto those in a theory with n ℓ flavors at the scale µ ∼ m. At fixed order, this implies the following consistency relation for the unpolarized and linearly polarized massive TMD matching coefficients, where M j λ /k λ denotes the PDF matching function, the sum runs over all light degrees of freedom, and the subscript λ = ∅, ∥ again labels the polarization of the heavy quark and the light partons j and k. Perturbatively expanding the matching functions as these relations simplify for our dimensionless O(α s ) coefficient functions in b T space, q λ /g λ (z, m, µ) , (3.22) where the µ dependence has to cancel within the matching coefficient. For the unpolarized case, this relation has previously been verified in refs. [46,47]. At NLO, the polarized PDF matching function relevant for our case is given by [106] M (1) The massless matching coefficient for the quark helicity TMD PDF onto the collinear gluon helicity PDF was calculated in ref. [100], (3.24) , it is straightforward to see that our result in eq. (3.19) indeed satisfies eq. (3.22).
By contrast, the worm-gear L matching coefficient is suppressed by one power of the mass, see eq. (3.14), and therefore cannot be reproduced by a leading-power PDF matching at the scale µ ∼ m. Interestingly, it contains a logarithm of mb T at subleading power, refs. [107][108][109][110][111], and it would be interesting to understand whether the logarithm in eq. (3.25) might be amenable to similar techniques.

Numerical results for TMD PDFs
For numerics, we evaluate eq. (3.12) at the boundary scales µ 0 ∼ √ ζ 0 ∼ 1/b T given in and below eq. (2.61), perform the TMD evolution back to µ = √ ζ = Q as described around eq. (2.62), and finally take a numerical Fourier transform as in eq. (3.5). E.g., we have for the evolved unpolarized heavy-quark TMD PDF, and similarly for the other cases. For the input collinear gluon PDFs we use the NNPDF31 nnlo as 0118 unpolarized proton PDF set [112] together with the NNPDFpol11 100 set for the polarized case [113]. Our input values for the strong coupling and the quark pole masses were given in section 2.6.
In figure 5, we show our numerical results for the heavy quark TMD PDFs for producing a charm or bottom quark from a longitudinally polarized proton as a function of k T and  x, respectively. The bottom quark TMD PDFs have a wider peak in k T compared to the charm because of its larger mass, as can be understood from the fact that the expressions in eq. (3.19) only depend on mb T up to RG effects. Note also that the worm-gear L function (after including a Jacobian 2πk T ) is quadratic in the small k T region with a coefficient proportional to 1/m 3 , whereas the unpolarized and helicity TMD PDFs are linear in k T . As this approximation is valid to higher k T in the case of the bottom quark than that of the charm, the bottom-quark TMD PDF has a numerically smaller value over a wide range. As x decreases, the unpolarized heavy-quark TMD PDF rises much more rapidly than the polarized ones, as expected from the smaller gluon polarization fraction at smaller x. We point out that the unpolarized TMD PDF changes sign at very high x ≥ 0.6, indicating a need for resumming subleading-power threshold logarithms of 1 − x using e.g. the tools of refs. [114].

Accessing heavy-quark TMDs in e + e − collisions
In e + e − collisions, TMD fragmentation functions may be accessed from double-inclusive measurements with two identified hadrons, e + e − → H a H b X. For instance, the six-fold differential cross section for this process in the TMD limit P a,T , M a,b ≪ Q is given by [115,116] where cos θ and ϕ are the spherical coordinates of hadron H b with respect to the incoming beams in the center-of-mass frame, z a and z b are the lightcone momentum fractions of the two hadrons, and ⃗ P a,T is the transverse momentum of hadron H a . On the right-hand side, α em is the fine-structure constant, Q is the center-of-mass energy of the collision, y = (1 + cos θ)/2, and ϕ 0 is the azimuthal angle of ⃗ P a,T measured relative to the plane spanned by H b and the beams. The hadronic structure functions factorize into TMD FFs, where F ee denotes a weighted sum over flavors and a convolution of two TMD FFs (i.e., a product in b T space) at total partonic transverse momentum q T = P a,T /z a , and the hard function describing the pair production of quarks is given by Here we have kept the contribution from Z boson exchange and Z-photon interference, as relevant for measurements on the Z pole, where P Z (Q 2 ) = Q 2 /(Q 2 − m 2 Z + iΓ Z m Z ) is the reduced Z propagator and e f (v f , a f ) are the electromagnetic charge (vector, axial coupling to the Z) of a fermion f . We may assume that the experimental measurement involves an integral over symmetric ranges in cos θ such that the forward-backward asymmetry and an associated odd Collins effect in eq. (4.1) drop out.
Crucially, the TMD factorization theorems in eqs. (4.1) and (4.2) only assume that the hard scale Q ∼ z a Q ∼ z b Q is large compared to all other scales, i.e., all masses and transverse momenta, and therefore hold for both light-quark and heavy-quark fragmentation at z a,b ∼ 1 without modification. In particular, the heavy quarks are approximately  massless at the scale µ ∼ Q at which they are produced, and their polarization states are thus fully entangled. The hard function in eq. (4.4) could be modified to account for the effect of perturbative spin flips, but this amounts to retaining power corrections in m/Q further suppressed by powers of α s . Importantly, this means that a characteristic cos(2ϕ 0 ) modulation (the Collins effect) is present both for light and for heavy quarks at leading power and at tree level. As is commonly done for light quarks, the Collins effect strength can be accessed by taking suitable ratios of weighted cross sections, which we here take to be integrated over z a and z b as likely relevant for an initial study of the heavy-quark Collins effect. In figure 6 we show the predicted e + e − → DDX or BBX cross sections as a function of hadron transverse momentum P a,T , and the Collins effect strength R cos(2ϕ 0 ) as a function of q T . The universality for charm and bottom quarks follows along the same lines as for figure 4, and holds as long as the center-of-mass energy is sufficient to produce the quark-antiquark pair in a boosted state. This is the case for charm mesons at typical continuum center-of-mass energies at existing B factories, and for both charm and bottom mesons at higher values of Q such as at the Z pole. The Collins effect is smaller at higher center-ofmass energies because χ ⊥ 1,H is linearly suppressed in b T compared to the unpolarized, which means it predominantly contributes at larger values of b T where the Sudakov suppression at higher energies tends to be stronger.
We show the results of varying κ H (κ ⊥ H ) for the unpolarized (Collins) TMD FF, and illustrate the variation of λ H⊥ by the yellow band, exactly as in figure 4. Note that the information about the absolute sign of the Collins function is lost in e + e − collisions, i.e., for two charge-conjugate hadrons we end up with a positive effect strength for any value of λ H⊥ = λH ⊥ since the effect is proportional to the square of the Collins function. One may nevertheless extract the relative factor between e.g. the D and D * Collins function, which heavy-quark spin symmetry predicts to be exactly minus one, see eq. (2.57), by measuring the Collins effect separately for e + e − → DDX and e + e − → D * D X. Explicitly, our prediction from heavy-quark spin symmetry reads We point out that for generic O(Λ QCD ) model parameters, the Collins effect strength reaches the several-percent level for continuum open charm production at existing B factories, in line with our expectation of an effect strength that is comparable to the light quark case, making a future dedicated measurement (or search) appear very feasible.

Comment on claims regarding a mass suppression of the Collins effect
In e + e − collisions, the "intrinsic" heavy-quark Collins effect we analyzed above has been disregarded so far. Note that this effect is in general distinct from the large background contribution of DD weak decays to e.g. a measurement of the Collins effect on a KK sample. This contribution is indeed considered in experimental analyses [41,42,44,45] and subtracted as a background using Monte-Carlo simulations and heavy-quark enriched samples, but cannot be immediately interpreted as a sign of a (nonperturbative) Collins effect since the progenitor DD pair in this case is not constrained to be near the back-toback limit by the measurement, meaning that e.g. perturbative gluon emissions can also induce azimuthal correlations on the DD pair and thus their weak decay products.
Ref. [44] mentions that it would be possible to look for the intrinsic heavy-quark Collins effect with some further improvements to their analysis, but also incorrectly expects that the Collins effect should be parametrically suppressed by the mass of heavy quarks. The argument sketched in that reference (see beginning of their section IV) is that helicity flips should wash out the spin correlation between the heavy quark and the antiquark. This is not the case, as we have argued above: The quarks are approximately massless at the scale µ ∼ Q at which they are produced, and thus are produced with fully entangled spin states, such that there is no suppression by the mass from physics at this scale. Similarly, in our detailed analysis of the Collins FF at the scale µ ≤ k T ≤ m, we find no suppression of the effect by the mass, and the Collins effect in particular is fully allowed by heavy-quark symmetry when accounting for the presence of lightlike Wilson lines. Note that this is not contradictory to the fact that we do find a suppression of the Collins effect by Λ QCD /k T at large k T , since this suppression is exactly commensurate with the twist suppression of the two Collins functions in the light-quark case, which has been mapped out extensively [41][42][43][44][45]. We conclude that the prospects for a measurement of the intrinsic, nonperturbative heavy-quark Collins effect at B factories are even better than anticipated in ref. [44].

Accessing heavy-quark TMDs at the future EIC
TMD fragmentation functions may also be accessed from single-inclusive measurements with one identified hadron in electron-nucleon collisions, e − (ℓ)+N (P ) → e − (ℓ ′ )+H(P H )+ X, where the scattering is mediated by an off-shell photon with momentum q = ℓ − ℓ ′ (and Q 2 ≡ −q 2 > 0). The fully differential cross section for this process in the TMD regime reads [3,66,93,117,118] On the left-hand side, x = Q 2 /(2P · q), y = (P · q)/(P · ℓ), z H = (P · P H )/(P · q), and ⃗ P H,T is the outgoing hadron's transverse momentum relative to ⃗ q in the Breit frame. On the right-hand side,  [93]. The beam polarization information is encoded in the lepton beam helicity λ e and the covariant nucleon spin vector S µ = (0, S T cos ϕ S , S T sin ϕ S , −S L ) as decomposed in the Trento frame. We have dropped terms proportional to S T , which cannot be populated by leading-power heavy-quark TMD PDFs, see section 3. We have also dropped terms proportional to the Boer-Mulders function, whose twist-2 matching in the heavy-quark case is suppressed by at least one additional power of α s . The hadronic structure functions factorize in terms of one TMD PDF and one TMD FF each,  Table 1: Total cross sections in picobarn for producing charm (left two columns) or bottom-quark hadrons (right two columns) in the TMD region at the future 18 × 275 GeV 2 EIC for different cuts on x > x min , Q > Q cut , q T = P H,T /z < q cut T . See the text for further details on the acceptance cuts we consider.
where the convolution in transverse momentum may be written in position space as [67] (4.10) and the hard function for scattering a quark off a virtual photon is As for e + e − collisions, the TMD factorization theorems in eq. (4.9) only assume that the hard scale Q ∼ zQ is large compared to all low scales, and thus hold for both light and heavy hadron production without modification. Again, the heavy quark is approximately massless at the hard scale such that helicity is conserved during the hard scattering. This means that while the production mechanisms for longitudinally or transversely polarized heavy quarks from an incoming nucleon are different from light quarks (and are fully perturbative), the way they imprint on the distribution of final-state hadrons is the same, leaving nonzero spin asymmetries . (4.12) In particular, the sin(2ϕ H ) modulation induced by a nucleon beam polarization flip gives direct access to the heavy-quark Collins function including its sign, which is not accessible in e + e − collisions.
To assess the statistical power of the future EIC to constrain charm and bottom quark TMD dynamics, we first estimate the expected sample size of heavy hadrons in electronproton collisions. To do so, we consider the total cross section for producing a heavy quark in the TMD region summed over beam polarizations, where Θ DIS (x, y) denotes DIS acceptance cuts given by x > x min , 0.01 < y < 0.95 , (4.14) We consider the EIC at beam energies E e = 18 GeV and E N = 275 GeV. Any experimental cuts on z H > z cut and the additional prefactor of z 3 H in eq. (4.13) are irrelevant at our working order because the heavy quark carries all the energy in all regimes, i.e., z H = 1 either at leading power in the heavy-quark expansion or at the leading perturbative order, see the comments below eq. (2.62). For this estimate we set κ H = 0 in the unpolarized heavy-quark TMD FF, since the total integral of the cross section up to q cut T ≫ Λ QCD is independent of it up to corrections of O(Λ 2 QCD /q cut T ) [75], and sum over all heavy hadrons containing the heavy quark, exploiting H χ H = 1. This means that the total rate at which heavy quarks are produced is predicted fully perturbatively, as expected.
Our results for the expected total charm and bottom-quark TMD cross sections are given in table 1 for Q cut = 4 GeV and Q cut = 10 GeV, where higher Q cut allows for mapping out the TMD region to higher q T before encountering power corrections, but at the cost of much lower rates. (We have also adjusted q cut T accordingly in each case.) Scaled to an integrated luminosity of 10 fb −1 , we expect a total charm quark sample of 35 × 10 3 in the TMD region for the loose cut on Q and in the region x > 0.1 where polarization effects are expected to be most pronounced, see figure 5, and where a measurement of the sin(2ϕ H ) asymmetry is the most promising. This suggests that even with this limited integrated luminosity, percent-level asymmetries should be statistically resolvable.
In figure 7 we show the results for the unpolarized SIDIS cross section with a D (B) meson in the final state, and for the two spin asymmetries defined in eq. (4.12). Note that the effect of different κ H in the unpolarized TMD fragmentation function is negligible in the cross section and the A LL asymmetry, which as expected are dominated by perturbative physics. The A LL asymmetry is very sizable at ∼ 30% at the chosen value of x = 0.2. On the other hand, the A U L asymmetry is substantially smaller (1 − 2%) for the generic O(Λ QCD ) parameters we picked here due to the smaller value of both h ⊥ 1L compared to g 1L and H ⊥ 1 compared to D 1 in most of the contributing TMD region, see figures 4 and 5 and the surrounding discussion. The numerically smaller value of h ⊥ 1L for bottom quarks discussed around figure 5 is likewise reflected in the size of the asymmetry for bottom compared to charm quarks. We emphasize that a measurement of A U L , compared to the Collins effect in e + e − collisions, has the unique benefit of accessing the absolute sign  of the heavy-quark Collins function. Resolving this sign should well be possible within the expected statistics at the future EIC. While we leave the study of systematic effects (such as luminosity uncertainties) to future work, we note that the requirements that the established heavy-flavor/gluon distribution program of the EIC places on instrumentation have already been analyzed in depth in ref. [27]. Among these requirements are secondary vertex reconstruction capabilities and the momentum resolution on soft pions from D decays, all of which will also benefit the kind of differential measurements of semi-inclusive heavy-quark fragmentation that we propose here.

Summary of main results and outlook
In this paper, we have studied the transverse momentum-dependent (TMD) dynamics of bottom or charm quarks with mass m ≡ m c , m b ≫ Λ QCD fragmenting into heavy hadrons for the first time. We considered two parametric regimes for the transverse momentum k T , (a) Λ QCD ≲ k T ≪ m, where the hadron transverse momentum k T is determined by nonperturbative soft radiation into the final state, and (b) Λ QCD ≪ m ≲ k T , where k T is set by perturbative emissions. We assumed throughout that the heavy quark is produced at a hard scale Q ≫ m, k T , i.e., it is boosted in the frame of the hard scattering, such that standard TMD factorization applies at the scale Q and only the low-energy TMD matrix elements are modified by the heavy quark dynamics. In both regimes, the dynamics at scales below the heavy quark mass are constrained by heavy-quark symmetry and encoded in novel low-energy matrix elements in boosted Heavy-Quark Effective Theory (bHQET): • We showed that in regime (a), the unpolarized and Collins TMD fragmentation functions (FF) match onto new, universal nonperturbative bHQET matrix elements χ 1,H (k T ) and χ ⊥ 1,H (k T ), which we dubbed TMD fragmentation factors.
• In regime (b), we made use of the twist expansion for light-quark TMD FFs and combined it with the matching collinear FFs onto bHQET to identify the relevant leading or subleading bHQET matrix elements.
• An important new ingredient in this analysis is the unpolarized partonic heavy quark TMD FF d 1 Q/Q , a perturbative Wilson coefficient that appears in our analysis for the first time and that we expect to appear also in other contexts like flavor-tagged energy-energy correlators in the back-to-back limit.
• We find that the Collins TMD FF scales as Λ QCD /k T at k T ≫ Λ QCD -but is not suppressed by the quark mass -and identified the coefficient as a new subleading bHQET matrix element probing gluon correlations within the fragmentation process.
• We used heavy-quark spin symmetry to express the TMD fragmentation factors in terms of the underlying spin density matrix of the light hadron constituents.
• For the unpolarized TMD FF D 1 H/Q (z H , k T , µ, ζ), this allowed us to prove the following relations between renormalized TMD FFs within spin symmetry multiplets: These relations are a powerful generalization of known results for inclusive heavyquark fragmentation, demonstrating for the first time that they hold point by point in transverse momentum.
• We showed that the Collins function arises from correlations between hadronic radiation into the final state and the transverse polarization of the light constituents, which in turn is correlated with the heavy-quark spin by the experimental reconstruction of e.g. D vs. D * mesons. This new picture of spin correlations in the heavy-quark limit allowed us to prove the following novel sum rule for the renormalized heavy-quark Collins function, which holds up to corrections of O(Λ QCD /m) and up to radiative corrections in α s at the scale µ ∼ m, and likely beyond, To extend our analysis to the possible phenomenology at the future Electron-Ion Collider (EIC), we also considered the production of polarized heavy quarks from a polarized nucleon, which is encoded in all-order matching relations between heavy-quark TMD PDFs and twist-2 collinear light-parton PDFs: • We find that terms proportional to the transverse nucleon polarization vanish for heavy quarks at twist-2 to all orders in the strong coupling due to chirality and flavor conservation in the light-quark sector.
• In contrast to the light-quark case, transverse quark polarization states are populated from unpolarized and linearly polarized nucleons because the quark mass breaks chirality.
• We find nontrivial matching coefficients at O(α s ) for the heavy-quark worm-gear L and helicity TMD PDFs onto the gluon helicity collinear PDF, both of which we computed explicitly for the first time. We anticipate that the heavy-quark Boer-Mulders function will receive a contribution from the twist-2 collinear gluon PDF starting at O(α 2 s ), where it becomes allowed by time-reversal invariance.
Combining the standard TMD factorization theorems for e + e − to hadrons and SIDIS with simple numerical models for the new nonperturbative functions we identified, we provided predictions for unpolarized heavy-quark TMD cross sections, the Collins effect strength for heavy quarks at e + e − colliders (and in particular for cc continuum production at current B factories), as well as for the relevant spin asymmetries at the future EIC: • We find that a measurement of the intrinsic heavy-quark Collins effect is well within reach of existing B factories, and is motivated by the rich nonperturbative structure of the heavy-quark Collins function that our analysis revealed.
• The fact that transversely polarized heavy quarks are produced from linearly polarized nucleons at a significant rate, as encoded in the worm-gear L matching coefficient, in addition provides a clean avenue for probing the heavy-quark Collins functions in heavy-quark SIDIS at the future EIC, including its absolute sign.
The theoretical framework we developed in this paper paves the way for many promising future applications: • While we only considered the case of unpolarized heavy hadrons in this work, an immediate next application of our framework are polarized vector mesons or baryons containing heavy quarks. This gives access to a larger set of transverse-momentum dependent polarized fragmentation functions [115,119,120] which in the heavy-quark case resolve the light spin density matrix in even greater detail and obey additional sum rules.
• Another promising prospect is to consider heavy-quark TMD fragmentation within jets, which makes its rich physics accessible in hadron collisions. This extension is in fact straightforward because our results for the heavy-quark TMD FFs hold independent of the factorization theorem they appear in. This makes it possible to insert them into the hadron-in-jet frameworks of refs. [121,122] in a plug-andplay fashion as long as Q ∼ p jet T R ≫ m, k T . Yet another possibility, which could mitigate the effect of nonglobal logarithms that can become nonperturbative in our regime of interest, would be to apply grooming to the jet and study the hadron transverse momentum spectrum with respect to the groomed jet axis [58,123], see also footnote 3. We look forward to the attendant phenomenology, which may in addition serve as a vacuum baseline for TMD interactions of open charm and bottom quarks with the quark-gluon plasma in heavy-ion collisions.
• Other natural extensions are higher-order calculations of the various new partonic matching coefficients we introduced in this paper, which will reduce the perturbative uncertainties on the lowest-order theory predictions we provided here. This will also involve analyzing the renormalon structure and optimizing the choice of quark mass scheme. In addition, one could consider the matching onto subleading bHQET fragmentation matrix elements (for TMD FFs) or onto twist-3 collinear PDFs (for TMD PDFs, extending the work of ref. [124] to the massive case), which would make it possible to interpret phenomenological extractions in terms of higher-point correlation functions. Higher-order resummed predictions for heavy-quark TMD spectra then immediately follow from our factorization results by solving the attendant renormalization group equations, and will serve as powerful, highly differential benchmarks of the heavy-quark physics encoded in present and future parton showers, including their interface with hadronization models, on which our field-theory analysis of the nonperturbative dynamics places rigorous constraints.
In conclusion, our analysis reveals that a wealth of information on the all-order and nonperturbative structure of QCD resides in the transverse momentum dependence of heavyquark fragmentation. An experimental exploration of this new subfield of TMD physics is in immediate reach of existing B factories and will be an exciting addition to the planned heavy-flavor physics program of the future EIC.