The transverse momentum distribution of hadrons within jets

We study the transverse momentum distribution of hadrons within jets, where the transverse momentum is defined with respect to the standard jet axis. We consider the case where the jet substructure measurement is performed for an inclusive jet sample pp → jet + X. We demonstrate that this observable provides new opportunities to study transverse momentum dependent fragmentation functions (TMDFFs) which are currently poorly constrained from data, especially for gluons. The factorization of the cross section is obtained within Soft Collinear Effective Theory (SCET), and we show that the relevant TMDFFs are the same as for the more traditional processes semi-inclusive deep inelastic scattering (SIDIS) and electron-positron annihilation. Different than in SIDIS, the observable for the in-jet fragmentation does not depend on TMD parton distribution functions which allows for a cleaner and more direct probe of TMDFFs. We present numerical results and compare to available data from the LHC.


Introduction
In recent years, studies of jets and their internal structure have played increasingly important roles in testing the fundamental properties of Quantum Chromodynamics (QCD), and in searching for new physics beyond the Standard Model [1,2]. This is the case in particular in the era of the Large Hadron Collider (LHC), where collimated jets of hadrons are abundantly produced.
In this paper we study the transverse momentum distribution of hadrons h within fully reconstructed jets in pp collisions, pp → (jeth)X, as illustrated in Fig. 1. Specifically we study the ratio F (z h , j ⊥ ; η, p T , R) = dσ pp→(jet h)X dp T dηdz h d 2 j ⊥ dσ pp→jetX dp T dη , (1.1) the center-of-mass (CM) frame in pp collisions. The large light-cone momentum fraction of the jet carried by the hadron h is denoted by z h and j ⊥ is the transverse momentum of the hadron with respect to the standard jet axis. Throughout this paper, bold letters represent two-dimensional transverse momentum vectors, whereas the magnitude of these vectors is referred to as, for example, j ⊥ = |j ⊥ |. This observable has been measured at the LHC in pp collisions for a wide range of jet transverse momenta p T [3]. In addition, it has been measured in both unpolarized pp and transversely polarized p ↑ p collisions at the Relativistic Heavy Ion Collider (RHIC) [4][5][6]. It was proposed in [7] that the latter case can be used to probe azimuthal spin correlations in the fragmentation process, in particular, the so-called Collins function [8].
R h jet j ? Figure 1. The Distribution of hadrons inside a fully reconstructed jet. Here, j ⊥ is the transverse momentum of hadrons with respect to the standard jet axis, and R is the jet radius.
In this work, we develop the theoretical framework to study the above observable F (z h , j ⊥ ; η, p T , R). We consider the case where the jet substructure measurement is performed for an inclusive jet sample pp → jet + X, different than the study in [9] where an exclusive jet sample was studied in the context of heavy quarkonium production. As the experimental measurements [3] were performed for inclusive jet samples, our approach facilitates a direct comparison with the experimental data. In particular, we concentrate on the region of the hadron transverse momentum where j ⊥ p T R. Here, j ⊥ is defined with respect to the standard jet axis, rather than a recoil-free axis, specifically the winner-takeall jet axis as discussed in [10]. While a recoil-free axis can be advantageous for various applications for collider physics, it turns out that there is only a direct relation to the standard transverse momentum dependent fragmentation functions (TMDFFs) when the standard jet axis is used. The standard TMDFFs are also probed in the traditional processes semi-inclusive deep inelastic scattering (SIDIS) and back-to-back hadron pair production in electron-positron annihilation.
Following earlier work on the longitudinal momentum distribution of hadrons inside jets [11][12][13][14][15][16], we can write down the factorized form of the cross section in pp collisions as follows (for more details, see Eq. (3.1) below) dσ pp→(jet h)X dp T dηdz h d 2 j ⊥ = a,b,c f a (x a , µ) ⊗ f b (x b , µ) ⊗ H c ab (x a , x b , η, p T /z, µ) ⊗ G h c (z, z h , ω J R, j ⊥ , µ) . (1.2) Here, f a,b denote the parton distribution functions (PDFs) in the proton with the corresponding momentum fraction x a and x b , respectively. The hard functions H c ab describe the production of an energetic parton c in the hard-scattering event. In addition, the functions G h c (z, z h , ω J R, j ⊥ , µ) are the semi-inclusive TMD fragmenting jet functions (siTMDFJFs) which we describes the production of a jet in the final state with the observed hadron inside. We define this new function in Sec. 2 below. We further demonstrate that the siTMDFJFs can be refactorized in terms of hard matching functions, soft functions, and the transverse momentum dependent fragmentation functions. It is evident that the in-jet fragmentation of hadrons considered in this work provides a very sensitive probe especially for the gluon TMDFF which is so far only poorly constrained by the traditional processes.
The remainder of this paper is organized as follows. In Sec. 2, we provide operator definitions for the siTMDFJFs G h c , which is the essential component in describing the hadron transverse momentum distribution inside jets. We derive a factorization formalism for siTMDFJFs in terms of hard functions, soft functions and the TMDFFs. We calculate all these relevant functions in the factorized expression to next-to-leading order (NLO) and derive their renormalization group equations. We solve the resulting renormalization group (RG) equations in order to resum all the relevant large logarithms. In Sec. 3, we provide a first numerical estimate for the hadron transverse momentum distribution inside jets for LHC kinematics, and we compare with experimental results. We summarize our paper and provide further discussions in Sec. 4.

The semi-inclusive TMD fragmenting jet function
In this section, we introduce the definition of the semi-inclusive TMD fragmenting jet functions G h c (z, z h , ω J R, j ⊥ , µ), which are the relevant new ingredients in order to describe the hadron transverse momentum distribution within jets produced in pp collisions. The siT-MDFJFs describe the fragmentation of a hadron h inside a jet that is initiated by a parton c. We first provide their operator definitions, and we then derive the factorization formalism in terms of hard functions, soft functions and TMDFFs. We derive the relevant RG equations and their solutions. In addition, we work out the relation to standard TMDFFs probed in SIDIS and electron-positron annihilation.

Definition
Following the convention for the definition of the semi-inclusive fragmenting jet function in [16], the siTMDFJFs are defined for quark and gluon jets as follows where N c is the number of the colors for quarks, and j ⊥ is the transverse momentum of the hadron h with respect to the standard jet axis. The large light-cone momentum components of the initial parton c, the jet, and the hadron h are given by ω, ω J , ω h , respectively. We choose to express the results of our calculation in terms of the following ratios of these variables as well as ω J which is related to the transverse momentum of the reconstructed jet. In Eq. (2.1), we have n µ = (1,n) andn µ = (1, −n) where the spatial componentn is chosen along the standard jet axis. In addition, we have n 2 =n 2 = 0 and n ·n = 2. The gauge invariant collinear quark and gluon fields within Soft Collinear Effective Theory (SCET) [17][18][19][20][21] are denoted by χ n and B µ n⊥ as usual, and P is the label momentum operator. The sum over states |X runs over all final state particles except for the observed jet J with the identified hadron h inside.
In [16,22], it was found that the characteristic momentum scale for the jet dynamics with a jet radius R is given by We would like to point out that for the standard jet algorithms in pp collisions, one can simply make the replacement ω J tan(R/2) → p T R, where the jet size R is defined in the (η, φ) plane, see e.g. Ref. [23,24]. Depending on the relative scaling of j ⊥ , µ J and Λ QCD one finds different factorization theorems for the hadron-in-jet cross section. First, we consider the case when j ⊥ is of the same order as the characteristic jet scale µ J , i.e., Λ QCD j ⊥ ∼ p T R. In this case, the siTMDFJF G h c as defined in Eq. (2.1) can be factorized into standard collinear fragmentation functions D h/i (z h , µ) convolved with Wilson coefficients in z h [9,25]. The Wilson coefficients are functions of z, ω J and j ⊥ , and can be calculated perturbatively.
Second, in the region where j ⊥ is much smaller than the characteristic jet scale µ J , i.e., Λ QCD j ⊥ p T R, the perturbative expansion is plagued with large logarithmic corrections of the form ln (p T R/j ⊥ ), and the standard collinear factorization breaks down. Therefore, in the small j ⊥ regime, a new factorization formalism -TMD factorization [26], is required to recover reliable predictions within QCD perturbation theory. This is the kinematic region that we address in this work.

Factorization theorem
We focus on the kinematic region with the relative scalings Λ QCD j ⊥ p T R, referred to as a TMD region. In this region, since the transverse momentum j ⊥ inside the jet is parametrically small, only collinear radiation within the jet with the momentum scaling p c = (p − c , p + c , p c⊥ ) ∼ p T (1, λ 2 , λ) where λ ∼ j ⊥ /p T , and the soft radiation of order j ⊥ are relevant to leading power of the cross section. 1 Harder emissions are only allowed outside the jet and, therefore, do not affect the hadron transverse momentum j ⊥ . Since j ⊥ is defined with respect to the jet axis, any radiation outside the jet will only influence the determination of the jet axis but have no impact on the j ⊥ spectrum. Therefore, we have the following factorized form for G h c derived within SCET where besides the usual renormalization scale µ, the scale ν arises due to the so-called rapidity divergences in the relevant functions to be discussed in detail below. Here, H c→i (z, ω J R, µ) are hard matching functions related to out-of-jet radiation. The soft functions S i (λ ⊥ , µ, νR) take into account soft radiation inside the jet and D h/i (z h , k ⊥ , µ, ν) are the usual TMDFFs, which characterize the collinear degrees of freedom inside the jet. The delta function relates the observed transverse momentum of the hadron j ⊥ to be the total transverse momentum of soft and collinear radiations. Note that λ ⊥ is multiplied by z h inside the delta function to account for the difference between the fragmenting parton and the observed hadron with respect to the jet axis. All the ingredients in the factorization formula will be calculated up to NLO, which determines their RG evolutions. All large logarithms of the form ln R and ln (p T R/j ⊥ ) are resummed by solving the obtained RG equations and by running each component from their natural scales to the hard scale µ ∼ p T at which the cross section is evaluated.

Hard functions
The hard matching functions H c→i (z, ω J R, µ) encode radiation with a virtuality of order O(p T R) outside of the jet. They describe how an energetic parton c from the hard-scattering event produces a jet initiated by parton i with energy ω J and radius R, and can be computed through the matching relation in Eq. (2.4). Up to NLO, they are obtained by the out-ofjet radiation diagrams for inclusive jet (substructure) observables [16]. The same hard matching functions were found in the context of central subjets measured on an inclusive jet sample in [24]. For anti-k T jets, the renormalized expressions are given by where the logarithm L is defined as 6) and the standard splitting functions are also provided here for reference The analogous hard matching coefficients for cone jets can be found in [24]. The RG equations of the functions H i→j take the following form Note that this is a set of four coupled equations with a DGLAP type structure. The anomalous dimensions γ ij (z, ω J R, µ) are given by The second terms of the γ ij (z, ω J R, µ) are the standard DGLAP evolution kernels. Instead, the first term contains a logarithm L and the functions Γ i (ω J R, µ) are given at this order by To summarize, the RG equations encountered here resum double logarithms whereas the DGLAP equations are always associated with the resummation of single logarithms. Evidently, the natural scale for the hard matching coefficients H i→j (z, ω J R, µ) is the same as the jet scale µ J as defined in Eq. (2.3), i.e., (2.14) Thus, by solving the above RG equations and by evolving the hard matching functions from the scale µ J ∼ p T R to the hard-scattering scale µ ∼ p T where the cross section is evaluated, we are resumming large logarithms of jet radius ln R.

TMD fragmentation functions
Now we focus on the TMDFFs D h/i (z h , k ⊥ , µ, ν), which are defined as where ω is the light-cone energy of the initiating quark or gluon. The TMDFFs contain rapidity divergences. We choose to employ the analytic rapidity regulator of [29] which introduces a dependence on the associated rapidity scale ν. Traditionally, TMDs are studied conveniently in Fourier transform space or b-space. Following the standard convention of Ref. [26,30], we define the TMDFFs in b-space as At the same time, through the inverse Fourier transform, we obtain the TMDFFs in momentum space as

Perturbative results
In perturbation theory, the bare TMDFFs suffer from infrared (IR), ultra-violet (UV), and rapidity divergences. To understand the features of these divergences, it is instructive to study the perturbative results for the TMDFFs in the region where k ⊥ Λ QCD . In the following, we consider perturbative splittings i → jk at the parton level, where j refers to the identified parton whose transverse momentum k ⊥ is measured. We denote the corresponding TMDFFs at the parton level by D j/i (z h , k ⊥ , µ, ν). Up to NLO, we find the following results Note that the results for the TMDFFs here are independent of the jet radius parameter R. The allowed collinear radiation inside the jet is so collimated along the jet axis in the kinematical limit that we are considering, j ⊥ /p T R, such that it is insensitive to the jet boundary. We now calculate the Fourier transform of these results according to Eq. (2.16). After expanding around η → 0 and then → 0, we obtain the following results in b-space Here we introduced the scale µ b which is defined as

Renormalization
In this section we perform the renormalization of the TMDFFs and derive the resulting RG equations. We observe that the poles of D q/q and D g/g in the second lines of Eqs. (2.19a) and (2.19c) are UV poles. Therefore, they are subtracted via the usual renormalization procedure. The bare and renormalized TMDFFs are related as For the relevant renormalization constants, we find We thus obtain the associated RG equation and the rapidity renormalization group (RRG) equation where the µand ν-anomalous dimensions are given by (2.27)

Matching onto collinear FFs
After the UV poles are removed via renormalization, the TMDFFs D j/i (z h , b, µ, ν) only contain IR poles which have the expected structure ∼ −1/ P ji (z h ). In the perturbative region 1/b Λ QCD , the TMDFFs can be further matched onto the standard collinear FFs D h/i (z h , µ). With the help of this matching procedure, the remaining IR poles can be subtracted. The matching relation is given by where the matching coefficients are denoted byC j←i . The perturbative results at the parton level for the collinear FFs D j/i (z h , µ) in the MS scheme are given by Together with the expressions for the TMDFFs in Eq. (2.19), we find that the matching coefficients in b-space are given bỹ

Soft functions
The soft function S i (λ ⊥ , µ, νR) is defined as where Y n(n) denotes the soft Wilson line, and P ∈J ⊥ indicates the fact that only the soft radiation inside the jet contributes to the hadron transverse momentum with respect to the jet axis. The soft functions S i (λ ⊥ , µ, νR) also contain rapidity divergences and, thus, the rapidity scale ν arises. The calculation of the soft functions in the perturbative region λ ⊥ Λ QCD is very similar to the standard global soft function that arises in the processes SIDIS, Drell-Yan, and electron-positron annihilation, except that now we restrict the soft radiation to be inside the jet. The final result up to NLO in momentum space is given by where we keep the leading contribution in the limit R 1. The color factors are given by , respectively. After taking the Fourier transform to b-space, we obtain Note that the same result was obtained in [24] in the context of central subjets. Similar to the renormalization of the TMDFFs discussed above, we subtract the UV poles of the soft functions. The renormalized and bare soft functions S i (b, µ, νR) are related by where the multiplicative renormalization constants Z S i are given by The associated RG and RRG equations are given by with the µand ν-anomalous dimensions (2.39)

Solution of the evolution equations and resummation
In this section, we provide the details about how to solve the RG and RRG equations derived above for three functions of the siTMDFJFs, i.e. the hard matching functions, the TMDFFs and the soft functions. The resummation of all large logarithms is obtained by the following two step process. First, we evaluate all fixed order results at their natural scales which eliminates all large logarithms. Second, we evolve all three functions from their natural scales to a common scale µ ∼ p T . Effectively, this procedure resums all large logarithms in the fixed-order results derived above. We are going to find that it is numerically more convenient to evolve the TMDFFs and the soft functions to the jet scale µ ∼ p T R and to combine the result at this scale with the hard matching functions to obtain the siTMDFJFs. Then, we evolve the thus obtained siTMDFJFs from p T R → p T using the RG equations for the combined siTMDFJFs rather than using the RG equations for the three separate functions. We are going to find that the siTMDFJFs satisfy the timelike DGLAP evolution equations like their collinear analogous, the semi-inclusive fragmenting jet functions (siFJFs) as studied in [16]. We show that both approaches for solving the RG equations are equivalent. Besides numerical simplifications, using a combined evolution for the siTMDFJFs also makes the relation to traditional TMDFFs more clear. Before discussing the details of the resummation, we start by introducing the traditional definition of "proper" TMDs that allow for a parton model interpretation of TMD sensitive observables.

Proper TMD definitions
A crucial feature of the results for the TMDFFs D h/i and the soft functions S i is that both have rapidity divergences, but their product D h/i S i is free of rapidity divergences as they exactly cancel. This can be seen clearly from the NLO expressions for D h/i in Eq. (2.19) and S i in Eq. (2.33). The same 1/η poles appear in both expressions but with opposite signs. Following the usual TMD phenomenology [26,31,32], we thus define the "proper" in-jet TMDFFs D R h/i as the product where the superscript R reminds us that it represents the hadron distribution within a jet of the radius R. The cancelation of rapidity divergences for D R h/i can be traced back to the fact that the soft radiation is restricted to be only inside the jet. Note that the TMDFFs D h/i for the in-jet calculation turned out to be the same as for other TMD sensitive observables [26,31,32] and they do not depend on R, as discussed above. However, the soft functions are different in the sense that the soft radiation is restricted to be only inside the jet. Instead, for the "global" soft functions that are relevant for SIDIS and electron-positron annihilation [33][34][35], there is no such phase space constraint. The additional phase space restriction encountered here cuts off half of the rapidity divergences compared to the global soft functions. This leads to the cancelation of the rapidity divergences in Eq. (2.40) for the product D R h/i = D h/i S i . To be more specific, we present results for the global soft functions as well. We usê S i (b, µ, ν) to denote the global soft function in Fourier transform space. Without any phase space constraints on the soft radiation, we obtain the following expression for the global soft function in momentum space [29,36] After taking the Fourier transform to b-space as in Eq. (2.33) and expanding around η, → 0, we findŜ Comparing this result with the S i (b, µ, νR) in Eq. (2.33), we find that the O(α s ) terms differ by an overall factor of 2 and ν ↔ ν tan(R/2). The "proper" standard TMDFFsD h/i as they appear in SIDIS and electron-positron annihilation are then defined aŝ This product is also free of rapidity divergences allowing for a parton model type interpretation of TMD sensitive observables [26,31,32]. It is important to work out the exact relation between the in-jet TMDFFs D R h/i considered in this work and the standard TMDFFsD h/i . We will discuss this relation in more detail after deriving the solution of the RG and RRG equations in the next section.

Hard matching functions
We start with the RG equations for the hard matching functions H i→j , see Eq. (2.11). Note that the anomalous dimensions γ ij (z, ω J R, µ) in Eq. (2.12) contain a purely diagonal piece δ ij δ(1 − z)Γ i (ω J R, µ) and the Altarelli-Parisi splitting functions P ji (z) similar to the timelike DGLAP. We are going to separate these two parts of the anomalous dimensions and the associated evolution. The purely diagonal or non-DGLAP pieces of γ ij (z, ω J R, µ) are going to cancel with the respective terms of the anomalous dimensions of the TMDFFs and the soft functions yielding a standard DGLAP evolution equation for the siTMDFJFs.
To that extend, we start by writing the functions H i→j as where the C i→j (z, ω J R, µ) satisfy evolution equations where the anomalous dimensions are given only by the Altarelli-Parisi splitting functions Note that these evolution equations are similar to DGLAP equations but here we still have four coupled equations. Only the combined siTMDFJFs are going to satisfy the standard timelike DGLAP evolution equations, see Eq. (2.71) below.
The functions E i (ω J R, µ) satisfy multiplicative RG equations where the Γ i (ω J R, µ) are given in Eq. (2.13). The solution for the multiplicative RG equations can be written as The fixed-order results for E i (ω J R, µ) can be obtained from Eq. (2.5) and are given by By choosing µ J = p T R, we obtain E i (ω J R, µ J ) = 1 as the initial condition for the evolution in Eq. (2.46). Using this result in Eq. (2.44) above, we can write the hard matching functions as The functions C i→j (z, ω J R, µ) still need to be evolved from µ ∼ µ J = p T R to µ ∼ p T using the evolution equations in Eq. (2.45) above. Their fixed order expressions are given by

Relation between in-jet and standard TMDFFs
We are now going to derive the solution of the evolution equations for the TMDFFs and the soft functions in order to obtain the proper in-jet TMDFFs D R h/i (z h , b; µ) as defined in Eq. (2.40). For comparison, we also show the results for the standard TMDFFsD h/i (z h , b; µ) as in Eq. (2.43). We start by evolving D h/i (z h , b, µ, ν), S i (b, µ, ν, R) andŜ i (b, µ, ν) using their RG and RRG equations. From the perturbative calculations above, we find that the natural scales for the TMDFFs D h/i , and the two soft functions S i ,Ŝ i are given by To be consistent with the standard Collins-Soper-Sterman (CSS) formalism [37], we evolve from the natural scales of D h/i and S i ,Ŝ i as given in Eq. (2.51), to a common scale µ and ν. In terms of the "proper" TMDFFs in (2.40), the initial conditions for the evolution is given by It might be instructive to point out that the "proper" TMDs chosen as such are equal perturbatively when evaluated at their natural scales, This can be directly verified from the perturbative expressions given above. At this point, it might be instructive to point out that according to Eq. (2.51), the natural rapidity scales for two soft functions S i andŜ i are quite different, ν S /νŜ = 1/ tan(R/2) 1 in the small jet radius limit R 1. Since the "proper" TMDs do not contain rapidity divergences anymore, the ν-dependence will naturally disappear in the end when one evolves to the common rapidity scales. After solving the corresponding RG and RRG equations, we can write the final result in the following form Here the cusp anomalous dimension Γ i cusp and the non-cusp γ i allow for a perturbative evaluation as Γ i cusp = n Γ i n−1 αs π n and likewise for γ i . The first coefficients can be obtained from our calculation and are given by The higher-order expressions can be found for example in Ref. [38]. The "proper" in-jet TMDs D R h/i in Eq. (2.55) may be further expressed in terms of the "proper" standard TMDŝ D h/i in Eq. (2.56) as To obtain the second line we made use of Eq. (2.54). In other words, the evolved "proper" TMDs obtained for the hadron distribution inside jets D R h/i (z h , b; µ) at scale µ is related to the standard TMDsD h/i (z h , b; µ J ) evaluated at scale µ J multiplied by an overall factor. This overall factor is given by an exponential involving an integration over µ from scales µ J to µ. Since scales µ J and µ are both in the perturbative regime, µ J , µ Λ QCD , we find that the relation between the in-jet TMDs D R h/i and the standard TMDsD h/i is purely perturbative.

Solution for the siTMDFJFs
We proceed by combining the above results in order to obtain the evolved siTMDFJFs G h c (z, z h , ω J R, j ⊥ , µ). Starting from Eq. (2.4) and by using the relation Now we can plug in the result of the hard matching coefficients H c→i (z, ω J R, µ) in Eq. (2.49) where we separated and solved the RG equations for the functions E i (ω J R, µ). In addition, we use the results of the evolved in-jet TMDs D R h/i (z h , b; µ) in Eq. (2.59). We find that the siTMDFJFs may eventually be expressed as It is important to note that here we are able to write the result in terms of the standard TMDsD h/i . This is possible since the evolution between the scales µ J and µ of E i (ω J R, µ) cancels with the overall multiplicative factor found in Eq. (2.59) for the in-jet TMDs when written in terms of the standard TMDs. Specifically, we have The result in Eq. (2.62) constitutes the most important part of our work. It explicitly demonstrates that the hadron transverse momentum distribution within jets is related to the standard TMDFFs (as measured in SIDIS and electron-positron annihilation) probed at the jet scale µ J ∼ p T R. Eventually, we can write the result as where we used the inverse Fourier transform as defined in Eq. (2.17) to obtain the TMDFFs in momentum spaceD

Final expression for the siTMDFJFs
In the perturbative region where 1/b Λ QCD , one can further match the TMDFFŝ Using the coefficientsC j←i given in Eq. (2.30) and the perturbative expressions for the soft functions, one obtains It might be instructive to point out that the above matching coefficients are computed in the standard MS scheme, which differs from the simplest minimal subtraction scheme by inserting a factor S for each loop in the counter-terms with S = (4πe −γ E ) . However, in the so-called Collins-11 definition of TMDs, this factor was changed to S JCC = (4π) /Γ(1 − ) [26]. We refer to the latter scheme as MS JCC , in which the π 2 terms are absent in So far, we have discussed the evolution of the siTMDFJFs in the perturbative region, i.e. for 1/b Λ QCD . It is well-known that the evolution of TMDs contains a non-perturbative component in the large-b region. We treat the large-b region by adopting the usual b *prescription [37]. Alternative approaches can be found in [39][40][41][42][43]. One defines b * as where the quantity b max is introduced such that b * → b at small b b max , whereas it approaches the limit b → b max in the large b-region. Using this prescription and the matching coefficients in Eq. (2.67), we can write the evolved TMDFFs in Eq. (2.65) aŝ Here, S i pert (b * , µ J ) is the perturbative Sudakov factor and S i NP (b, µ J ) is the non-perturbative Sudakov factor. We will discuss them in detail in the phenomenological Sec. 3 below.
With all these relevant ingredients available, we may then compute the siTMDFJFs following Eq. (2.64). By using the evolution equations for the functions C i→j in Eq. (2.45) and the expressions for siTMDFJFs in Eq. (2.64), we find that the siTMDFJFs satisfy the standard timelike DGLAP evolution equations This result for the evolution equations of the siTMDFJFs was to be expected. Following our factorization expression for the differential cross section in Eq. (1.2), the product H c ab G h c should be µ-independent order by order analogously to single inclusive hadron production [44,45]. Thus, it is natural that the G h c follow the same DGLAP evolution equations as those for the usual collinear FFs D h/c [16,22].
We would like to summarize again the following aspects of the evolution structure of the siTMDFJFs. The TMD part of the evolution between µ b an µ J is governed by the same TMD evolution equations that have been obtained for the standard TMDs as well. The hard matching functions C c→i follow RG equations where the anomalous dimensions are given by the usual Altarelli-Parisi splitting functions. The evolution of C c→i is carried out between the jet scale µ J and the hard scale µ allowing for the resummation of logarithms in the jet size parameter R. The structure and resummation of single logarithms α n s ln n R becomes more apparent when combining both contributions to obtain the siTMDFJFs. They follow the standard DGLAP structure as it is usually associated with the resummation of single logarithms in the jet size parameter [16,22,46,47]. The obtained structure for the siTMDFJFs provides a convenient method to perform the resummation of all relevant large logarithms. First, we are going to evolve the standard TMDFFsD h/i (z h , j ⊥ ; µ) from the scale µ b to µ J . Second, at the jet scale µ J the TMDFFsD h/i (z h , j ⊥ ; µ J ) will be combined with the remaining hard matching functions C c→i (z, ω J R, µ J ) in Eq. (2.50) to compute the Then, we use the DGLAP evolution equations for the siTMDFJFs in Eq. (2.71) to evolve G h c from the scale µ J to µ and, thus, resum logarithms in the jet size parameter R.
3 Phenomenology for pp → (jeth)X In this section, we present numerical results for the transverse momentum distribution of hadrons inside jets for LHC kinematics. We consider an inclusive jet sample pp → jet + X where a hadron h is identified inside the reconstructed jet. Following [10,14,16,22], the factorization theorem for the process pp → (jeth)X can be written as where f a and f b denote the parton distribution functions (PDFs) in the proton with the corresponding momentum fraction x a and x b , respectively. For all numerical calculations in this work, we choose the CT14 NLO set of PDFs [48]. The hard functions H c ab describe the production of an energetic parton c in the hard-scattering event. They have been calculated analytically up to NLO in [45,49]. The variablesŝ,p T andη denote the partonic CM energy, and the transverse momentum and rapidity of parton c, respectively. They are related to their hadronic analogues aŝ where z is the momentum fraction transferred from parton c to the observed jet. The lower integration limits x min a , x min b and z min can be found for example in [14,16]. Finally, the functions G h c (z c , z h , ω J R, j ⊥ , µ) in Eq. (3.1) are the siTMDFJFs as discussed in Section 2. We would like to stress that the cross section does not depend on TMDPDFs but only on the standard collinear PDFs. Unlike the TMDPDFs, collinear PDFs are very well constrained by data and have been determined in global fits in the literature. Therefore, different than for SIDIS, the hadron in-jet fragmentation considered in this work provides an opportunity to disentangle the effects of TMDPDFs and TMDFFs.
In order to perform numerical calculations, we have to parameterize the non-perturbative Sudakov factors for both quark and gluon TMDFFs. Unfortunately, the quark TMDFFs are not very well constrained so far. The main information for the extraction of quark TMDFFs are obtained from multiplicity distributions of hadrons measured in SIDIS from both the HERMES [50] and COMPASS [51] experiments. These measurements were performed at relative low momentum scales, with photon virtualities Q 2 of several GeV 2 . Thus, there are potential problems when interpreting the data in terms of the usual leading-twist TMD factorization formalism [52]. In addition, since the factorization for SIDIS involves a convolution of TMDPDFs and TMDFFs, the unambiguous extraction of both functions separately is not straightforward. Therefore, current extractions of quark TMDFFs are subject to large uncertainties. Keeping in mind the remaining large uncertainties in our calculation, we are nevertheless going to present numerical estimates for the hadron transverse momentum distribution within jets and compare to LHC measurements. We choose to use the following parametrization of the non-perturbative Sudakov factor following [53,54] with Q 2 0 = 2.4 GeV 2 , b max = 1.5 GeV −1 , g 2 = 0.84, and g h = 0.042. Other parametrizations for the non-perturbative Sudakov factor for TMDFFs have been discussed in [30,55].
Furthermore, we note that the non-perturbative Sudakov factor for the gluon TMDFF is not constrained at all so far. For our numerical calculations, we are going to follow [56][57][58] and adopt a parameterization of the gluon non-perturbative Sudakov factor similar to that for quarks as In comparison to the quark parametrization, the coefficient of the term ∼ ln µ J is enhanced by a color factor C A /C F , whereas the intrinsic part ∼ g h is kept unchanged. In addition, we use the coefficients C j←i in Eq. (2.67) up to the order of α s for the TMDFFs, and keep Γ i 0,1 and γ i 0 in the perturbative Sudakov factor in Eq. (2.70), which is often referred to as the next-to-leading-logarithm prime (NLL ) accuracy. For both C j←i in Eq. (2.67) and C i→j in Eq. (2.50), we use the expressions in the MS JCC scheme as explained in Sec. 2.7, since the TMDFFs that we use for our numerical studies were extracted within this scheme [53,54]. We use the DSS07 parametrization of collinear fragmentation functions for light charged hadrons [59]. Together with the choices for the relevant nonperturbative inputs above for both quark and gluon TMDFFs, we are now going to present first numerical estimates for the transverse momentum distribution of hadrons inside jets and compare to the data provided by the ATLAS collaboration [3]. We choose the following jet kinematics consistent with the available data at a CM energy of √ s = 7 TeV. The jets are reconstructed using the anti-k T algorithm with jet size parameter R = 0.6, and the jet rapidity is integrated over |η| < 1.2. The detailed numerical implementation is very similar to the longitudinal momentum distribution of hadrons inside jets [16], using several numerical techniques developed in the literature [60][61][62]. The RG evolution of the various parts of the cross section is performed as outlined at the end of the last section.
In Fig. 2, we present the comparison of our numerical results and the LHC data for the hadron j ⊥ -distribution inside jets. We make the default scale choices of µ = p T and µ J = p T R. We explore the scale uncertainty by varying µ and µ J independently by a factor of two around their default values and by taking the envelope of these variations. As an example, we choose the jet transverse momentum bins 25 < p T < 40 GeV (left) and 400 < p T < 500 GeV (right). The experimental data are presented for the z h -integrated hadron distribution, i.e. with z h integrated from 0 to 1. This fact hinders a more direct and transparent comparison of our results with the data, since the collinear FFs are only constrained in a finite region z min h < z h < 1 with z min h 0.05 [59,63]. Any z h < z min h is not constrained and can only be obtained by extrapolation. We choose the value for z h in our calculations as the average value z h that are provided in the experimental publication [3], with z h = 0.08 and 0.03 for 25 < p T < 40 GeV and 400 < p T < 500 GeV, respectively. With this caveat in mind, we find that our calculations based on TMDFFs extracted in the literature at low energy scales of several GeV give a reasonable description of the experimental data. The height of the peak is roughly consistent with the data but our results have a broader j ⊥ -distribution than the experimental data. We note that at low jet p T our current numerical estimates agrees somewhat better with the data than in the high p T region.  Figure 3. Breakdown of the hadron j ⊥ -distributions inside jets (blue) into quark initiated (red) and gluon initiated (green) TMDFF channels.
In Fig. 3, we separate the hadron j ⊥ -distribution into quark and the gluon TMDFF components. This separation is valid in the TMD region. We find that at lower jet p T , as shown in the left panel of Fig. 3, the gluon channel dominates over the quark channel due to the overwhelmingly abundant gg initiated events in the pp collisions. The quark TMD fragmenting contribution is suppressed in this region. Therefore, the low jet p T region provides a "golden channel" to extract the gluon TMDFF. At large jet p T (right panel of Fig. 3), where qg initiated events start to dominate, the quark and the gluon TMDFF contributions therefore become comparable to each other. However, due to the difference in the color charges carried by quarks and gluons, the quark TMD fragmenting process peaks at smaller j ⊥ . Away from the peak region, the quark contribution drops more dramatically and exhibits a relatively narrow spectrum compared with the gluon contribution. Therefore, the region away from the peak of the j ⊥ -spectrum will generally be more sensitive to the gluon TMDFF.
To conclude this section, we provide further discussions of our numerical estimates. First, we would like to emphasize that at the moment we only concentrate on the TMD region and we present numerical results without matching onto NLO fixed-order calculations. In other words, we have not considered the effect of the so-called Y -term which can also affect the low j ⊥ -region, as advocated recently in [64,65]. Second, in the TMD region, the non-perturbative parts of the quark TMDFFs have large uncertainties as they have only been constrained from SIDIS data so far, while the gluon TMDFF has not been extracted at all. Third, so far we did not take into account the effect of non-global logarithms (NGLs) [66,67]. They first arise at next-to-next-to leading order due to the hierarchies caused by different constraints in different phase space regions, and affect our results at the current logarithmic accuracy we are considering. The factorization and resummation of NGLs have been studied recently in great detail, see for example Refs. [68][69][70][71][72][73][74]. We expect to obtain significant improvements of our results in the comparison with the experimental data once all these additional factors are taken into account. A dedicated study including all additional effects will be presented in a forthcoming publication.

Conclusions
In this work, we have studied the hadron transverse momentum j ⊥ -distribution within jets, where j ⊥ is defined with respect to the standard jet axis. We set up a factorization formalism that allows for systematic studies of this distribution. As a first step we calculated all the components of the factorization theorem to NLO, and we further resummed all the associated large logarithms ln R and ln(p T R/j ⊥ ). We demonstrated the universality of the TMDFFs that arise for this jet substructure observable and the traditional TMDFFs probed in SIDIS and electron-positron annihilation. We further showed that the hadron distribution within jets produced in pp collisions provides a unique opportunity to study the TMDFFs, especially the gluon TMDFF. For SIDIS and electron-positron annihilation, the gluon TMDFF is usually difficult to access. More specifically, we showed that different than for SIDIS, the j ⊥ spectrum within jets only depends on TMDFFs. There is no dependence on TMDPDFs which allows for a more direct extraction of TMDFFs. Furthermore, we found that at LHC energies we are able to control the sensitivity to different TMDFFs by selecting different values of the jet p T . We observed that the low jet p T region is the ideal region to extract the gluon TMDFF. For large jet p T , the region away from the peak of the j ⊥ -spectrum can also be sensitive to the gluon TMDFF.
In the future, several extensions of this work are possible. For instance, in order to extend our calculations to the region where j ⊥ ∼ p T R, we need to match the resummed result onto fixed order calculations. Such a matching calculation includes the full NLO corrections to this spectrum which may also affect the TMD region. In addition, it will be important to study the numerical impact of NGLs. Besides improvements of the perturbative calculation, a more careful study of the non-perturbative Sudakov evolution will be necessary to determine whether the agreement with the data in the region j ⊥ < 1 GeV can be improved. Also given the relative simple structure of the TMDFFs and the soft functions considered here, a next-to-next-to leading order calculation is possible which will further push forward the accuracy of the theoretical predictions. In this work, we only considered the phenomenology of pp collisions but the formalism developed here is also directly applicable to ep scattering relevant for a future Electron-Ion Collider (EIC). Other phenomenological studies may include for example a global fit of TMDFFs using also data from SIDIS and electron-positron annihilation. Finally, we are also planning to extend our formalism to the polarized case which is crucial to probe the Collins function, Hyperon polarization inside jets and other types of jet substructure observables.