TMD Fragmentation Functions at N$^3$LO

We compute the unpolarized quark and gluon transverse-momentum dependent fragmentation functions (TMDFFs) at next-to-next-to-next-to-leading order (N$^3$LO) in perturbative QCD. The calculation is based on a relation between the TMDFF and the limit of the semi-inclusive deep inelastic scattering cross section where all final-state radiation becomes collinear to the detected hadron. The required cross section is obtained by analytically continuing our recent computation of the Drell-Yan and Higgs boson production cross section at N$^3$LO expanded around the limit of all final-state radiation becoming collinear to one of the initial states. Our results agree with a recent independent calculation by Luo et al.


Introduction
Highly energetic scattering processes allow us to test our understanding of fundamental interactions with incredible precision.It is the asymptotic freedom of strong interactions of QCD that allows us to contrast our first principle understanding of the interactions with experimental data.The interacting elementary particles of QCD -quarks and gluons -are however concealed in our observations as they form hadronic bound states as the strong interactions confine at long distances.The gateway that bridges the world of partonic interactions, where observables are calculable in perturbative QCD, to the observation in real-life detectors is provided by factorisation theorems.Factorisation theorems split the long range, confining part of a scattering process from the short range collision process of quarks and gluons.The long range part of this factorisation theorems is typically expressed in terms of parton distribution functions (PDFs) and fragmentation functions (FFs).PDFs and FFs are independent of the particularities of the scattering process and are universal, such that they can be measured and used in many different experiments and observables.
Longitudinal FFs are the simplest example of FFs, as they only describe the probability of a quark or a gluon to convert to a hadron that carries a given momentum fraction of the fragmenting parton [1][2][3][4].This notion is expanded by transverse-momentum dependent FFs (TMDFFs) [5][6][7][8][9][10][11][12], which encode the probability of a hadron to arise from a fragmenting parton with a certain fraction of the partons longitudinal momentum and a small transverse momentum relative to the parton.
TMDFFs are intrinsically nonperturbative objects, as they relate the dynamics of partons and hadrons, and as such have been extracted from various experiments [23,32,[46][47][48][49][50][51][52].However, for transverse momenta q T that are much larger than the confinement scale Λ QCD , an operator product expansion in Λ QCD /q T allows one to express each TMDFF in terms of a standard longitudinal FF and a q T -dependent matching kernel.The matching kernels are calculable order by order in perturbation theory and are currently known at next-to-next-to-leading order (NNLO) [53][54][55] in perturbative QCD.In the regime of perturbative q T , they can be used for example in extractions of longitudinal FFs from differential measurements of suitable observables, see for example refs.[56][57][58][59].
In this article we present the calculation of the matching kernels for all unpolarized quark and gluon TMDFFs at N 3 LO.TMD parton distribution functions (TMDPDFs), the initial-state counterparts of TMDFFs, are already known at this order [60,61].Our calculation relies on a recently developed method to expand hadron collider cross sections around the limit where final state QCD radiation is collinear to an incoming parton [62].We demonstrate explicitly how partonic cross sections for the production of electro-weak gauge bosons can be related to DIS cross sections via analytic continuation.We apply this analytic continuation to the collinear limit of the partonic cross section of gluon-fusion Higgs boson production and Drell-Yan production to obtain their DIS counter part.The collinear limit of these production cross sections was recently computed by us for the calculation of the N 3 LO TMDPDFs [60] and N -jettiness beam function [63].We then establish an analytic relation among the TMDFF matching kernels and the newly obtained collinear limit of the DIS cross sections.The combination of this collinear limit with the TMD soft function yields the scheme-independent TMDFF at N 3 LO.With this we are finally able to extract the desired perturbative TMDFF matching kernels.
The paper is structured as follow.In section 2, we setup the kinematics for SIDIS retaining full information on the momentum of the final state hadron.In section 3, we show how to use crossing symmetry and analytic continuation to obtain results for fully differential partonic cross sections in SIDIS from analogous cross sections in proton-proton collision.In section 4, we study the behaviour of the partonic cross section when taking the radiation to be collinear either to the struck proton or to the final state hadron.In section 5, we make use of the framework developed in the previous sections to extract the TMDFFs at N 3 LO by imposing a transverse-momentum measurement to the leading collinear expansion of the cross section.We conclude in section 6.

Setup
In this section we introduce our notation for the description of semi-inclusive deep inelastic scattering (SIDIS), reviewing both the scattering process and providing the definitions of all required kinematic variables and the associated final-state phase space.Finally, we define the transverse momentum observables of interest in this article.

Semi Inclusive Deep Inelastic Scattering
We study cross sections for the production of a hadron H in DIS alongside additional radiation, which we indicate as a multiparticle state X.In particular, we focus on the hadronic part of the DIS cross section that is initiated by the scattering of a proton with momentum P 1 and an electro-weak boson h with the space-like momentum q, P (P 1 ) + h(q) → H(−P 2 ) + X(−k) . (2.1) Here, we take all momenta to be incoming.This process is schematically depicted in figure 1 for the example of a virtual photon as the electro-weak gauge boson.In this article, we will consider DIS with either a virtual photon or a Higgs boson as the electro-weak gauge boson.
We are interested in SIDIS, where we measure an observable O that depends on the final-state hadronic momenta.For perturbative O Λ QCD , the cross section differential q p 2 p 1 k H(P 2 ) P (P 1 ) Figure 1.Schematic picture of the DIS process in eq.(2.1), producing a final-state hadron H in the scattering of an electroweak boson, here a photon, off the incoming proton.
in O can be factorized as Here, the overall normalization σ0 is the Born cross section and the sum runs over parton flavors i, j.In eq.(2.2), ηi+h→j+X is a perturbatively calculable partonic coefficient function encoding the underlying partonic process i + h → j + X, which is convolved with the nonperturbative parton distribution function (PDF) f i and fragmentation function (FF) d H/j .The PDF f i (x) encodes the probability to extract the parton with flavor i and momentum fraction x from the proton, while the FF d H/j (y) describes the fragmentation of a parton of flavor j into a hadron of type H which carries the momentum fraction y of the parent parton.We define the hadronic invariants In analogy, we introduce the partonic variables z and ξ, (2.4) The convolution integrals abbreviated by ⊗ x B and ⊗ x F in eq.(2.2) can now be written explicitly as where the partonic coefficient function is given by (2.6) Here we introduced the normalization factor N i related to the helicity and color average of the incoming particle, which for an incoming quark or gluon takes the value In eq.(2.6), the sum runs over the number m of additional partons in the final state besides the parton of flavor j that fragments into the hadron H, and Φ 1+m is the associated m + 1parton phase space.The δ functions implement the measurements of ξ and O, and the squared matrix element |M i+h→j+m | 2 corresponds to the partonic process of producing the m + 1 partons in the collision of a parton of flavor i with the hard probe h.

Kinematics and Final State Phase Space
We are interested in observables differential in the four momentum P 2 of the final state hadron H, while we are inclusive over all additional final-state radiation.A convenient set of variables to describe the kinematics of the corresponding partonic process is given by Here, k is the sum of all m final state momenta of the particles produced in addition to the parton with momentum p 2 , (2.9) The differential m + 1-parton phase space is given by It can be parameterized using the variables in eq.(2.8) as where (2.12) The kinematic variables are defined in the following domains, We can now express the desired partonic coefficient function defined in eq.(2.6) in terms of the partonic coefficient function differential in the above variables, The second line is the central object in this work, from which all desired observables can be easily projected out.It can be expanded as Here, we have expanded ηi+h→j+X in the strong coupling constant α s /π, and denote the coefficients as η ij for brevity.In the second line we have split off the terms η V ij which arise purely from Born contributions and virtual corrections.The remaining functions η ( ,m,n) ij are separately holomorphic in the vicinity of w 1 = 0 and w 2 = 0.
The benefit of using the variables defined in eq.(2.8) is that together with q 2 they fully specify the momentum p 2 and thus are sufficient to express in ηi+h→j+X differential in p 2 .For example, the Lorentz-invariant momentum fractions defined in eq.(2.4) are given by (2.16)

Transverse Momenta
In SIDS, two particular definitions of transverse momentum play a key role.These two different definitions of transverse momentum are most naturally measured in two different inertial frames.We define the infinite momentum frame (also referred to as Breit frame) and the hadron frame as follows: Infinite Momentum Frame Hadron Frame q = (0, 0, Q) q = (q 0 , q T , q z ) . (2.17) Here, E 1 and E 2 represent the energies of the initial and final state hadrons, respectively.The explicit vectors in the above table are Euclidean vectors.The momentum component | q T | of the momentum q is orthogonal to the plane spanned by the momenta P 1 and P 2 of the hadrons and is most naturally measured in the hadron frame.The momentum component | P 2T | of the momentum P 2 is orthogonal to the plane spanned by the momenta q and P 1 and is most naturally measured in the infinite momentum frame.We express both transverse momenta in terms of Lorentz invariant quantities by Here, is the invariant mass of the dihadron system.Inserting the parametrisation in terms of w 1 , w 2 and x as defined in eq.(2.8), the two transverse momenta of interest can be expressed in a Lorentz-invariant fashion as Below, we will be mostly interested in the limit that | P 2T | 2 , or equivalently | q T | 2 , becomes small.We will approach this limit by considering the limit w 2 → 0, for which one obtains the simple relation lim (2.21)

Crossing from Production to DIS Cross Sections
In the previous section, we introduced the SIDIS process P (P 1 ) + h(q) → H(−P 2 ) + X(−k) for the scattering off an electroweak boson h off the proton P , thereby producing a detected final-state hadron H in association with additional hadronic radiation.The associated cross section is related by a factorization theorem to the partonic process i(p 1 ) + h(q) → j(−p 2 ) + X(−k), which we describe by the partonic coefficient function η ij where we are fully differential in p 1 and p 2 , but integrate over k.
We now want to relate, i.e. cross, this partonic configuration to the one where both partons are in the initial state and produce an outgoing electroweak boson h, which we hence refer to as "production".Concretely, we study the crossing relation where, as always, we choose all momenta as ingoing.
Recently, we have studied this production process refs.[62,64,65].In particular, in ref. [62] we showed that the corresponding partonic coefficient function is given by where the differential phase space for h + n partons is given by Here, all variables are identical to the ones introduced in section 2.2 for SIDIS.In particular, note that the squared matrix elements are identical in the DIS and production case up to . . .< l a t e x i t s h a 1 _ b a s e 6 4 = " o c g i A j p D j k x J r / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " o c g i A j p D j k x J r / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " o c g i A j p D j k x J r / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " o c g i A j p D j k x J r / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " o c g i A j p D j k x J r / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " o c g i A j p D j k x J r x j Z g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " x j Z g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " x j Z g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " r l N 3 7 8 5 r j e u i j j I c w h G c g A s X 0 I B b a E I L C K T w D K / w Z j 1 Z L 9 a 7 9 T E f L V n F z g H 8 g f X 5 A 5 B x k w I = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 s 4 l 5 V H 7 6 U A r X n P 4 H H z K p 9 r l N 3 7 8 5 r j e u i j j I c w h G c g A s X 0 I B b a E I L C K T w D K / w Z j 1 Z L 9 a 7 9 T E f L V n F z g H 8 g f X 5 A 5 B x k w I = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 s 4 l 5 V H 7 6 U A r X n P 4 H H z K p 9 r l N 3 7 8 5 r j e u i j j I c w h G c g A s X 0 I B b a E I L C K T w D K / w Z j 1 Z L 9 a 7 9 T E f L V n F z g H 8 g f X 5 A 5 B x k w I = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 s 4 l 5 V H 7 6 U A r X n P 4 H H z K p 9 the crossing of momenta p 2 and q.Furthermore, in both cases, production and DIS, the final state radiation is integrated over the phase space dΦ m .The dependence of the cross sections on the momenta p 2 and q is fully retained.
In order to relate the partonic coefficient function of DIS to production, or vice versa, we need to understand the analytic structure of the partonic coefficient function.Crossing p 2 and q changes the sign of the numerical value of the invariants s and w 1 , and consequently it is important to understand the analytic branch structure of the partonic coefficient functions at s = 0 and w 1 = 0. Since s is the only variable with explicit mass dimension in our choice of independent variables, it immediately follows that the partonic coefficient function at O(α n s ) depends on s only through the multiplicative factor s −n .The analytic dependence on w 1 was already hinted at in eq.(2.15), but needs to be investigated in more detail.
The partonic coefficient function comprises of amplitudes interfering with complex conjugate amplitudes, integrated over the m-parton phase space.This can be further split into interference of l-loop Feynman diagrams with conjugate k-loop Feynman diagrams, as illustrated in figure 2. Similar to the decomposition of the partonic coefficient function in eq.(2.15), the analytic structure of the depicted interference diagram can be decomposed as Here, the functions f (i 1 ,i 2 ,j 1 ,j 2 ) (w 1 , w 2 , x) do not contain any branch cuts at s = 0, w 1 = 0 or w 2 = 0.When performing a computation of analytic partonic coefficient functions, it is easy and often useful to keep track of the individual functions f (i 1 ,i 2 ,j 1 ,j 2 ) (w 1 , w 2 , x).
The second line in eq.(3.4) differs between DIS and production kinematics due to the different signs of s and w 1 .Explicit phases occur in a given loop amplitude depending on the kinematic configuration of the external momenta.The phases are easily determined by equipping the Lorentz-invariant scalar products s, sw 1 and sw 2 with a definite Feynman prescription, Crossing from DIS to production kinematics then requires us to analytically continue the second line of eq.(3.4).As an example, we consider the case i The same analytic structure as outlined above for the interference of two Feynman diagrams naturally holds for the entire partonic coefficient function as well, dη Once the universal functions η are identified, it is easy to perform the analytic continuation between DIS and production kinematics.The above was observed and explicitly verified for the computation of the ingredients of Higgs and DY production up to N 3 LO in refs.[65][66][67][68][69] and holds in particular for the interference of amplitudes for massless QCD corrections for the processes under consideration.We note that it is of course also possible to relate DIS or production kinematics to partonic cross sections where only the electroweak gauge boson is in the initial state and all partons are in the final state, for example e + e − annihilation.
In addition to the analytic continuation from DIS to production kinematics there are some other, trivial differences in the partonic coefficient functions.First, the overall normalisation factor N i and N production ij differ, which can be trivially accounted for.Second, the phase space measure dΦ 1+m and dΦ production h + m differ by factors depending on the kinematic variables.However, this difference is accounted for by a simple multiplicative factor that does not require any additional analytic continuation.With this we have identified all differences between DIS and production kinematics in bare, partonic coefficient functions and can relate one to the other as long as the required analytic information is retained in the computation of one of them.

Collinear limit of partonic coefficient functions
In this section we briefly review the method introduced in ref. [62] to expand cross sections in the kinematic limit where all final-state radiation becomes collinear to the parton with momentum p 1 or p 2 .In order to illustrate this, it is instructive to decompose the momentum k into its components along these directions, Here, the k ⊥ component is chosen orthogonal to p 1 and p 2 .In order to illustrate the collinear limit with respect to either massless parton we introduce an auxiliary rescaling parameter λ and indicate the collinear limit by The respective limit is then achieved by taking λ → 0. The variables w 1 , w 2 and x defined in eq.(2.8) were chosen such that the action of either collinear rescaling transformation in eq. ( 4.2) on the partonic coefficient function simply amounts to a rescaling of w 1,2 .Specifically, in the p 1 -collinear limit only w 1 is rescaled, while in the p 2 -collinear limit only w 2 is rescaled, while the other variables are not affected, An expansion of our partonic coefficient function in the p 1,2 -collinear limit is thus equivalent to an expansion in w 1,2 .More details on how such an expansion can be performed for multiloop partonic coefficient functions can be found in ref. [62].
A key difference between the p 1 -and p 2 -collinear limit is that the former corresponds to a collinear initial-state singularity, which were already discussed in refs.[60,62,63], while the latter corresponds to collinear final-state singularity.Here, we only only briefly look at the impact of the p 1 -collinear limit on the more familiar variables given in eqs.(2.16) and (2.18), Note, that the p 1 -collinear limit of the phase space is identical for DIS and production kinematics up to the domain of the variables, lim Furthermore, in the strict p 1 -collinear limit, which is defined by only retaining momentum modes in loop integrals where the loop momentum itself is collinear to p 1 [62], none of the partonic coefficient functions require any analytic continuation when crossing between DIS and production kinematics.Thus, up to overall normalization factors the strict p 1 -collinear limit agrees between production and DIS kinematics.Of course, this is an immediate consequence of the universality of collinear dynamics of QCD and the factorization of collinear initial-state singularities.
The limit of all final-state radiation becoming collinear to the momentum p 2 corresponds to collinear final-state singularities, which were not discussed in ref. [62] and are the main focus of this article.In this limit, the familiar variables in eqs.(2.16) and (2.18) become p 2 −collinear : Note, that ξ in the p 2 -collinear limit behaves reciprocal to z in the p 1 -collinear limit, which is a consequence of their definition in eq.(2.4).In contrast to the p 1 -collinear limit, in the p 2 -collinear limit the phase space for DIS and production kinematics differ slightly by lim Furthermore, in order to cross from production to DIS kinematics it is necessary to analytically continue parts of the partonic coefficient function, as outlined in section 3.

Calculation of the TMD Fragmentation Functions
In this section we calculate the TMDFFs at N 3 LO from a perturbative calculation of the SIDIS process P (P 1 ) + h(q) → H(−P 2 ) + X(−k).We will briefly review the required factorization for SIDIS in the limit of small transverse momentum in section 5.1, before showing in section 5.2 how it relates to the kinematic limit where the final-state momenta P 2 and k are collinear to each other.In section 5.3, we discuss our results for the TMDFFs.

SIDIS factorization at small transverse momentum
We consider the unpolarized SIDIS process in eq.(2.1) in a frame where the incoming proton P and outgoing hadron H are aligned along the lightcone vectors in eq.(5.1), i.e.
In this frame, the momentum q µ of the electroweak boson h is given by In particular, it has a nonvanishing transverse momentum q T .Note, that the above coordinates correspond to the hadron frame introduced in sec.2.3.The factorization of the SIDIS cross section in the limit of small transverse momentum, q T Q, was first derived in [11] and elaborated on in refs.[85,86].We follow the notation established in the treatment of TMD factorization within Soft-Collinear Effective Theory (SCET) [87][88][89][90] in the formalism of the rapidity renormalization group equation [91,92].For Drell-Yan like processes, the factorized cross section is given by see appendix A for more details.In eq.(5.5), σ0 is the same Born cross section as before, the sum runs over all parton flavors i, j contributing to the Born process i + h → j, and the hard function H ij encodes virtual corrections to the Born process.As is common, the factorization in eq.(5.5) is expressed in Fourier space, with b T Fourier conjugate to q T .The TMD beam and fragmentation functions Bi (x B , b T ) and DH/j (x F , b T ) encode the effect of radiation collinear to the incoming proton and outgoing hadron, respectively, and are defined below.They depend on b T and the momentum fractions x B,F as defined in eq. ( 2.3).The soft function Sq (b T ) encodes the transverse recoil due to soft radiation, and is independent of the quark flavors i and j.Eq. (5.5) depends not only on the common renormalization scale µ, which we take as usual as the MS scale, but also on the scale ν that arises from the regularization of so-called rapidity divergences [5,86,[91][92][93][94][95][96][97][98], for which we employ the exponential regulator of ref. [92].The momentum fractions ω a,b in eq.(5.5) are defined as the lightcone components (5.6) They are closely related to the Collins-Soper scale ξ a,b ∝ ω 2 a,b [5,6].For gluon-induced processes, the factorized cross section reads dσ dx F d 2 q T = σ0 x 2 F 2H ρσρ σ (q 2 , µ) (5.7) The only difference to eq. (5.7) is the Lorentz structure of Bρσ g and Dρσ H/g , which arises due to the helicity structure of the gluon field, One can decompose the gluon TMDFF as where we suppressed the scales for brevity.The decomposition of Bρσ g has the same structure as eq.(5.8).We will only consider Higgs production, where due to the scalar nature of the Higgs boson and thus we only require the combination where we suppressed all arguments for brevity.Since this structure is very similar to the combination in eq.(5.5), in the following we will always use the form in eq.(5.5), with the implicit understanding that the B g D g term has to be added for Higgs production.Furthermore, since B g = O(α s ) and D g = O(α s ), their NNLO results are sufficient to describe Higgs production at N 3 LO.They have already been calculated in ref. [53], and we will not consider them in our calculation of the N 3 LO TMDFFs.
Before proceeding, we remark that the precise form of the beam, fragmentation and soft functions in eqs.(5.5) and (5.7) depends on the chosen rapidity regulator.In our work, we will use the exponential regulator of ref. [92], and the ensuing rapidity renormalization scale is denoted as ν, but many other rapidity regularization schemes are known in the literature [86,91,[93][94][95]98].The scheme ambiguity can be eliminated by combining the beam and fragmentation functions with the soft function, which in our case reads (5.11) These combinations are manifestly ν independent, reflecting the independence of the rapidity regulator.As is common, we have introduced the so-called Collins-Soper scale ζ a,b = ω 2 a,b , a remnant of the rapidity regularization.Note that while there is an established notation distinguishing TMD beam functions Bi and TMDPDFs fi , so far no such notation exists for the TMDFF.To make clear which function we refer to, we will label the TMDFF including the soft function by an explicit superscript "TMD".
Similar combinations as in eq.(5.11) can be constructed in all viable rapidity regulators (often, these functions are combined at the bare level prior to renormalization, with UV renormalization applied to the product).To calculate the TMDFF itself, one has to calculate collinear and soft matrix elements separately, and hence it is natural to separate the fragmentation functions from the soft function.Thus, we will provide results both for the scheme-dependent DH/j before and the scheme-independent DTMD H/j after combination with the soft function.
The TMDFFs in eqs.(5.5) and (5.7) are well-defined QCD hadronic matrix elements [99].Using SCET notation, the bare fragmentation functions are defined as Here, we make explicit that we regulate UV divergences by working in d = 4−2 dimensions and regulate rapidity divergences using the exponential regulator of ref. [92].In eq.(5.12), the sum is over all additional hadronic final states X, the trace is over color and spin, and P is the momentum of the hadron H.The fields χ n and B µ n⊥ are collinear quark and gluon fields in SCET, with the pair of fields in each equation separated by b µ = (0, b − , b ⊥ ).The matrix elements in eq.(5.12) are defined in the hadron frame as specified in eq. ( 5.3), i.e. the outgoing hadron H defines the lightcone direction nµ , and b T is transverse to it.
For perturbative b T Λ −1 QCD , the TMDFF can be matched perturbatively onto the collinear FF.For the renormalized TMDFF, this relation reads [85,86] where the matching coefficients Cjj are perturbatively calculable.In the second line in eq. ( 5.13) we have replaced z → x F /z, which will be more convenient for our extraction of Cjj .Corrections to eq. ( 5.13) are suppressed as O(b T Λ QCD ).With this we may rewrite the TMDFF of eq.(5.13) in terms of the Mellin convolution (5.15) The TMDFFs in eq.(5.12) are defined in a coordinate system where P µ 2 = P + 2 nµ /2 defines the lightcone direction and has vanishing transverse momentum, and hence b ⊥ is Fourier-conjugate to the transverse momentum of the parton that initiates the fragmentation process.Alternatively we may consider the transverse momentum of the final state hadron P 2T which is naturally defined in the infinite momentum frame, see sec.2.3, and following ref.[54] we denotes this definition of the TMDFF FH/j .Since the two transverse momenta in these two frames are related as P T 2 = −x F q T , see eq. (2.21), the two TMDFFs are related by (5.16) The first relation is an immediate consequence of eq.(2.21).The second equation immediately follows upon Fourier transform in d − 2 dimensions.In ref. [54], the matching relation for the FH/j was written as (5.18)

TMD fragmentation functions from the collinear limit
The TMDFF can be obtained from the collinear limit of SIDIS following the same strategy applied in refs.[60,62] to calculate the TMDPDF from the collinear limit of proton-proton scattering.We start from the cross section differential in the transverse momentum q T , which in the limit of small q T Q is given by the factorization theorem in eq.(5.5), The key insight is that the hard, beam, fragmentation and soft function in eq. ( 5 Here we used that the hard and soft functions are normalized to unity at tree level, while the TMD beam function reduces to the PDF itself.Note that eq.(5.20) is to be understood at the bare level, as only combining it with all other limits will cancel all appearing infrared divergences, and thus q T and b T are treated in d − 2 dimensions.We also used that both photon and Higgs exchange are flavor diagonal to fix j = ī.We want to relate eq.(5.20) to the SIDIS cross section defined in collinear factorization.Combining eqs.(2.5) and (2.14), we obtain where the expressions for ξ and q 2 T are given by eqs.(2.16) and (2.20).In the limit that all final state radiation becomes collinear to P 2 , i.e. w 2 → 0, all required variables becomes p 2 −collinear : Note that in this limit, the partonic coefficient function scales as δ(1 − z), see eq. (2.11), and thus renders the convolution in z trivial.Furthermore, we fix w 1 = −(1 − ξ)/ξ, and obtain lim strict p 2 −coll.
Comparing eqs.(5.20) and (5.23), we can immediately read off the relation between the perturbative matching kernel and take the Fourier transform with respect to q T , DH The perturbative matching kernel as defined in eqs.(5.13) -(5.15) is then given by Knaive The superscript "naive" in eq. ( 5.25) indicates that this is not yet the final result for the (bare) matching coefficient.First, we note that we still have to regulate rapidity divergences that arise as w 2 → 0, or equivalently ξ → 1.In our approach, this regulator must only act on the total momentum k.The only known regulator in the literature that fulfills this constraint is the exponential regulator [92], which amounts to inserting a factor exp[2τ e −γ E k 0 ] into the integral.In our parameterization, this regulating factor reads where in the last step we neglect the w 1 term that is not required to regulate the w 2 → 0 limit and use the momentum fraction ω b of eq.(5.6).Since eq. ( 5.26) vanishes exponentially as ξ → 1 and x → 1, it regulates all rapidity divergences in the p 2 -collinear sector.We identify the rapidity regularisation scale as as τ has inverse mass dimensions.Secondly, the TMDFF is defined as the purely collinear limit of the cross section, but the above matrix element still contains overlap with the soft factor.Its subtraction is referred to as zero-bin subtraction [100].In the case of the exponential regulator, this is that the first four terms in the threshold expansion of the collinear limit of the limit of ηij used here matches the collinear expansion of the threshold expansion of refs.[64,111].
In figure 3 we illustrate our results by showing the three-loop matching kernel K ij (z) in all quark channels (left) and gluon channels (right).The different channels have been rescaled as indicated in the figure to account for their different magnitudes.
For completeness, in appendix B we present the ξ → 0 limit of the kernels, both for the quark and the gluon TMD fragmentation functions.These results are interesting for the study of the high energy behavior of TMDFFs, similar to studies of the small-x behavior of TMDPDFs in refs.[53,60,[112][113][114][115].Note that the timelike TMDFF shows a doublelogarithmic series in ln ξ, such that the N 3 LO coefficient contains up to α 3 s ln 5 ξ, in contrast to the single-logarithmic series observed for the spacelike TMDPDF, where one encounters at most α 3 s ln 2 ξ at this order.

Conclusions
We have computed the perturbative matching kernel relating transverse-momentum dependent fragmentation functions (TMDFFs) with longitudinal fragmentation functions at N 3 LO in QCD, obtaining analytic results for all partonic channels contributing to the quark and unpolarized gluon TMDFF.These results for this matching kernel, defined in eq.(5.15), are provided as ancillary files together with the arXiv submission of this article.Our calculation is based on a simple extension of a framework recently developed by us, that allows to expand differential hadronic cross sections efficiently in the collinear limit [62].This method was developed in detail in ref. [62] for the collinear expansion of differential hadron collider production cross sections.We have demonstrated explicitly how they are related to DIS cross sections via analytic continuation.By analytically continuing our recent computation of the collinear limit of the gluon fusion Higgs boson and DY production cross section to DIS kinematics, we have obtained the TMDFFs in similar fashion as the N -jettiness beam functions and TMDPDFs calculated in refs.[60,62,63].Our new results demonstrate once more the potency of this method obtaining universal ingredients arising in the infrared and collinear limits of QCD to an unprecedented level of precision in perturbation theory.
An important check on our calculation lies in the cancellation of all infrared and ultraviolet poles against suitable counterterms.Since these counterterms can be fully predicted using known anomalous dimension, this provides a highly nontrivial check.In particular, it involves the cancellation of infrared divergences against the QCD mass factorisation counterterm comprised of time-like splitting functions.Thus, as a by product, our calculation confirms the recent results for the NNLO timelike splitting function ref. [74][75][76]80], in particular the correct result in the qg channel first obtained in ref. [80].
There are several phenomenological applications of our results.Firstly, the TMDFFs obtained in this paper constitute the last missing ingredient to describe the singular structure of the transverse momentum distribution of QCD radiation in color-singlet decays at N 3 LO.They also enable the resummation of transverse momentum distributions at N 3 LL accuracy, both in e + e − annihilation and Higgs decay to quarks or gluons as well as in SIDIS.In particular they allow for the calculation of the jet functions for the Energy-Energy Correlator (EEC) and the Transverse EEC jet functions in the back-to-back limit [41,43] at N 3 LO.For the case of the EEC this allows to push the resummation accuracy to N 3 LL which constitutes the most accurate resummation carried for an event shape to date.We carry out this calculation in ref. [42].
Our method to expand cross sections around the collinear limit in the final state can be used to calculate higher order terms in the collinear expansion.Such higher order terms would allow one to study the structure of factorization beyond leading power for IRC safe observables in e + e − annihilation and Higgs decay [116][117][118][119][120][121][122][123][124][125][126][127][128][129][130][131] as well as the appearance of subleading power rapidity divergences [98,132,133].Furthermore, they would provide data to validate the resummation of power suppressed logarithms [130,134].It would also be interesting to explore the application of the methods developed here and in ref. [62] to TMDFFs involving a jet measurement [58,135,136].
Note: While this article was under completion, an independent calculation was made available on the arXiv in ref. [137] based on the method proposed in ref. [80].The authors of ref. [137] provided an important cross check on intermediary results for genuine two loop contributions in the K gq channel that allowed us to track an error in a routine related to the analytic continuation of the partonic coefficient functions.The initial discrepancy was a non-logarithmically enhanced finite and rational term proportional to (C A − C F )ζ 2 ζ 3 in the K gq and K qg channel.After this was resolved, we find perfect agreement among all analytic results.Topical Collaboration.M.E. is also supported by the Alexander von Humboldt Foundation through a Feodor Lynen Research Fellowship.For the quark channels, the high-energy limit z → 0 of the kernels K(3) qi (z) contributing to the quark TMD fragmentation function is given by The expressions for the high energy limit z → 0 up to O(z 40 ), as well as that for the threshold limit z → 1 up to O((1 − z) 40 ), can be found for all channels in electronic form in the ancillary files of this work.
H N U c 6 L 8 + 5 8 L F r X n G L m B P 7 A + f w B y 2 + P Q g = = < / l a t e x i t > . . .< l a t e x i t s h a 1 _ b a s e 6 4 = " o c g i A j p D j k x J r 5 k P W C 1 B B + y V G n k = " > A A A B 7 X i c b V B N S 8 N A E J 3 4 W e t X 1 a O X Y B E 8 l U Q E P R a 9 e K x g P 6 A N Z b P Z t G s 3 u 2 F 3 U i i h / 8 G L B 0 W 8 + n + 8 + W / c t j l o 6 4 O B x 3 s z z M w L U 8 E N e t 6 3 s 7 a + s b m 1 X d o p 7 + 7 t H x x W j o 5

2 < 1 <
r 5 e a 8 N r P u E x S g 5 I t F 4 W p I C Y m 8 7 / J k C t k R k w t o U x x e y t h Y 6 o o M z a d k g 3 B W 3 1 5 n b R r V c + t e v d X l c Z N H k c R z u A c L s G D O j T g D p r Q A g Y j e I Z X e H O E 8 + K 8 O x / L 1 o K T z 5z C H z i f P w E E j Z k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D j S 1 o O A r 1 V R / V a t z f C T q K / y m J g E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K U I 9 F L x 4 r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p I R n U B u W K W 3 U X I O v E y 0 k F c j Q H 5 a / + M G Z p h N I w Q b X u e W 5 i / I w q w 5 n A W a m f a k w o m 9 A R 9 i y V N E L t Z 4 t T Z + T C K k M S x s q W N G S h / p 7 I a K T 1 N A p s Z 0 T N W K 9 6 c / E / r 5 e a 8 N r P u E x S g 5 I t F 4 W p I C Y m 8 7 / J k C t k R k w t o U x x e y t h Y 6 o o M z a d k g 3 B W 3 1 5 n b R r V c + t e v d X l c Z N H k c R z u A c L s G D O j T g D p r Q A g Y j e I Z X e H O E 8 + K 8 O x / L 1 o K T z 5 z C H z i f P w E E j Z k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D j S 1 o O A r 1 V R / V a t z f C T q K / y m J g E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K U I 9 F L x 4 r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p I R n U B u W K W 3 U X I O v E y 0 k F c j Q H 5 a / + M G Z p h N I w Q b X u e W 5 i / I w q w 5 n A W a m f a k w o m 9 A R 9 i y V N E L t Z 4 t T Z + T C K k M S x s q W N G S h / p 7 I a K T 1 N A p s Z 0 T N W K 9 6 c / E / r 5 e a 8 N r P u E x S g 5 I t F 4 W p I C Y m 8 7 / J k C t k R k w t o U x x e y t h Y 6 o o M z a d k g 3 B W 3 1 5 n b R r V c + t e v d X l c Z N H k c R z u A c L s G D O j T g D p r Q A g Y j e I Z X e H O E 8 + K 8 O x / L 1 o K T z 5 z C H z i f P w E E j Z k = < / l a t e x i t > p l a t e x i t s h a 1 _ b a s e 6 4 = " D j S 1 o O A r 1 V R / V a t z f C T q K / y m J g E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K U I 9 F L x 4 r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p I R n U B u W K W 3 U X I O v E y 0 k F c j Q H 5 a / + M G Z p h N I w Q b X u e W 5 i / I w q w 5 n A W a m f a k w o m 9 A R 9 i y V N E L t Z 4 t T Z + T C K k M S x s q W N G S h / p 7 I a K T 1 N A p s Z 0 T N W K 9 6 c / E / r 5 e a 8 N r P u E x S g 5 I t F 4 W p I C Y m 8 7 / J k C t k R k w t o U x x e y t h Y 6 o o M z a d k g 3 B W 3 1 5 n b R r V c + t e v d X l c Z N H k c R z u A c L s G D O j T g D p r Q A g Y j e I Z X e H O E 8 + K 8 O x / L 1 o K T z 5 z C H z i f P w E E j Z k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D j S 1 o O A r 1 V R / V a t z f C T q K / y m J g E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K U I 9 F L x 4r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W yx W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p I R n U B u W K W 3 U X I O v E y 0 k F c j Q H 5 a / + M G Z p h N I w Q b X u e W 5 i / I w q w 5 n A W a m f a k w o m 9 A R 9 i y V N E L t Z 4 t T Z + T C K k M S x s q W N G S h / p 7 I a K T 1 N A p s Z 0 T N W K 9 6 c / E / r 5 e a 8 N r P u E x S g 5 I t F 4 W p I C Y m 8 7 / J k C t k R k w t o U x x e y t h Y 6 o o M z a d k g 3 B W 3 1 5 n b R r V c + t e v d X l c Z N H k c R z u A c L s G D O j T g D p r Q A g Y j e I Z X e H O E 8 + K 8 O x / L 1 o K T z 5 z C H z i f P w E E j Z k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D j S 1 o O A r 1 V R / V a t z f C T q K / y m J g E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K U I 9 F L x 4 r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p I R n U B u W K W 3 U X I O v E y 0 k F c j Q H 5 a / + M G Z p h N I w Q b Xu e W 5 i / I w q w 5 n A W a m f a k w o m 9 A R 9 i y V N E L t Z 4 t T Z + T C K k M S x s q W N G S h / p 7 I a K T 1 N A p s Z 0 T N W K 9 6 c / E / r 5 e a 8 N r P u E x S g 5 I t F 4 W p I C Y m 8 7 / J k C t k R k w t o U x x e y t h Y 6 o o M z a d k g 3 B W 3 1 5 n b R r V c + t e v d X l c Z N H k c R z u A c L s G D O j T g D p r Q A g Y j e I Z X e H O E 8 + K 8 O x / L 1 o K T z 5 z C H z i f P w E E j Z k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D j S 1 o O A r 1 V R / V a t z f C T q K / y m J g E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K U I 9 F L x 4 r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p I R n U B u W K W 3 U X I O v E y 0 k F c j Q H 5 a / + M G Z p h N I w Q b X u e W 5 i / I w q w 5 n A W a m f a k w o m 9 A R 9 i y V N E L t Z 4 t T Z + T C K k M S x s q W N G S h / p 7 I a K T 1 N A p s Z 0 T N W K 9 6 c / E / r 5 e a 8 N r P u E x S g 5 I t F 4 W p I C Y m 8 7 / J k C t k R k w t o U x x e y t h Y 6 o o M z a d k g 3 B W 3 1 5 n b R r V c + t e v d X l c Z N H k c R z u A c L s G D O j T g D p r Q A g Y j e I Z X e H O E 8 + K 8 O x / L 1 o K T z 5 z C H z i f P w E E j Z k = < / l a t e x i t > p l a t e x i t s h a 1 _ b a s e 6 4 = " P G 2 3 2 O

4 L 8 6 7 8 7 F
o L T n F z D H 8 g f P 5 A / 9 x j Z g = < / l a t e x i t > A l < l a t e x i t s h a 1 _ b a s e 6 4 = " e V s y b ot A M E x / O q + U K x e f d C 7 1 Y 0 4 = " > A A A B 9 H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u q G 5 c V 7 A P a o d x J 0 z Y 0 k x m T T K E M / Q 4 3 L h Rx 6 8 e 4 8 2 / M t L P Q 1 g O B w z n 3 c k 9 O E A u u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / n h u b P w U l e F U s F m p m 2 g W I x 3 j k H U s l R g y 7 a f z 0D N y Z p U + G U T K P m n I X P 2 9 k W K o 9 T Q M 7 G Q W U i 9 7 m f i f 1 0 n M 4 N p P u Y w T w y R d H B o k g p i I Z A 2 Q P l e M G j G 1 B K n i N i u h I 1 R I j e 2 p Z E v w l r + 8 S p o X V c + t e g + X l d p t X k c R T u A U z s G D K6 j B P d S h A R S e 4 B l e 4 c 2 Z O C / O u / O x G C 0 4 + c 4 x / I H z + Q P 3 O J I 2 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " e V s y b o t A M E x / O q + U K x e f d C 7 1 Y 0 4 = " > A A A B 9 H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u q G 5 c V 7 A P a o d x J 0 z Y 0 k x m T T K E M / Q 4 3 L h R x 6 8 e 4 8 2 / M t L P Q 1 g O B w z n 3 c k 9 O E A u u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / n h u b P w U l e F U s F m p m 2 g W I x 3 j k H U s l R g y 7 a f z 0D N y Z p U + G U T K P m n I X P 2 9 k W K o 9 T Q M 7 G Q W U i 9 7 m f i f 1 0 n M 4 N p P u Y w T w y R d H B o k g p i I Z A 2 Q P l e M G j G 1 B K n i N i u h I 1 R I j e 2 p Z E v w l r + 8 S p o X V c + t e g + X l d p t X k c R T u A U z s G D K6 j B P d S h A R S e 4 B l e 4 c 2 Z O C / O u / O x G C 0 4 + c 4 x / I H z + Q P 3 O J I 2 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " e V s y b o t A M E x / O q + U K x e f d C 7 1 Y 0 4 = " > A A A B 9 H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u q G 5 c V 7 A P a o d x J 0 z Y 0 k x m T T K E M / Q 4 3 L h R x 6 8 e 4 8 2 / M t L P Q 1 g O B w z n 3 c k 9 O E A u u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / n h u b P w U l e F U s F m p m 2 g W I x 3 j k H U s l R g y 7 a f z 0D N y Z p U + G U T K P m n I X P 2 9 k W K o 9 T Q M 7 G Q W U i 9 7 m f i f 1 0 n M 4 N p P u Y w T w y R d H B o k g p i I Z A 2 Q P l e M G j G 1 B K n i N i u h I 1 R I j e 2 p Z E v w l r + 8 S p o X V c + t e g + X l d p t X k c R T u A U z s G D K6 j B P d S h A R S e 4 B l e 4 c 2 Z O C / O u / O x G C 0 4 + c 4 x / I H z + Q P 3 O J I 2 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " e V s y b o t A M E x / O q + U K x e f d C 7 1 Y 0 4 = " > A A A B 9 H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u q G 5 c V 7 A P a o d x J 0 z Y 0 k x m T T K E M / Q 4 3 L h R x 6 8 e 4 8 2 / M t L P Q 1 g O B w z n 3 c k 9 O E A u u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / n h u b P w U l e F U s F m p m 2 g W I x 3 j k H U s l R g y 7 a f z 0 D N y Z p U + G U T K P m n I X P 2 9 k W K o 9 T Q M 7 G Q W U i 9 7 m f i f 1 0 n M 4 N p P u Y w T w y R d H B o k g p i I Z A 2 Q P l e M G j G 1 B K n i N i u h I 1 R I j e 2 p Z E v w l r + 8 S p o X V c + t e g + X l d p t X k c R T u A U z s G D K 6 j B P d S h A R S e 4 B l e 4 c 2 Z O C / O u / O x G C 0 4 + c 4 x / I H z + Q P 3 O J I 2 < / l a t e x i t > A ⇤ k < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 s 4 l 5 V H 7 6 U A r X n P 4 H H z K p 9 T i

Figure 2 .
Figure 2. Schematic picture of the interference of a l-loop Feynman diagram with a complexconjugate k-loop Feynman diagram.
) and thus our kernels K are identical to their kernels with rescaled arguments, Kjj ξ, b T , µ

Figure 3 .
Figure 3.The N 3 LO TMD fragmentation function boundary term K(3) ij (z) as a function of z.The matching kernels in all channels entering the quark fragmentation (left) and the gluon fragmentation (right) are displayed.For illustration purposes the kernels have been rescaled as indicated.
199) encode different dynamics.The hard function H arises from hard virtual corrections to the Born process, while B, D and S are constructed such that they only arise from p 1collinear, p 2 -collinear and soft momenta in loop integral and real emissions, respectively.It follows that by calculating the strict p 2 -collinear limit, defined such that both loop and real momenta are expanded in the p 2 -collinear limit, only the fragmentation function contributes to eq. (5.19), lim strict p 2 −coll.