QCD resummation on single hadron transverse momentum distribution with the thrust axis

We derive the transverse momentum dependent (TMD) factorization and resummation formula of the unpolarized transverse momentum distribution ($j_T$) for the single hadron production with the thrust axis in electron-positron collision. Two different kinematic regions are considered, including small transverse momentum limit $j_T \ll Q$, and joint transverse momentum and threshold limit $j_T \ll Q(1-z_h) \ll Q$, where $Q$ and $z_h$ are the hard scattering energy and the observed hadron momentum fraction. Using effective theory methods, we resum logarithms $\ln(Q/j_T)$ and $\ln(1-z_h)$ to all orders. In the end we present the differential cross sections and Gaussian widths calculated for the inclusive charged pion production and find that our results are consistent with the measurements reported by the Belle collaboration.

The electron-positron collider provides a clean environment to study TMD FFs using the inclusive hadron production, since there is no hadronic contamination from the initial states. The standard process to probe TMD FFs in e + e − collisions is the aforementioned back-to-back hadron pair production, e + e − → h 1 h 2 X, which probes the same TMD FFs as those in SIDIS process, e − p → e − hX 1 . Recently, the single-hadron differential cross section for the process, e + e − → hX, is reported by the Belle collaboration [48], where the hadron cross section is studied as a function of the event-shape variable called thrust T , fractional energy z h , and the transverse momentum j T with respect to the thrust axis. The j T distribution shows the usual Gaussian shape and this gives the hope that through such a new measurement, one would gain better understanding of the same TMD FFs. With such an assumption, some phenomenological work has been performed in [49,50].
In this paper we perform a detailed theoretical study for the Belle observable and we develop a TMD factorization formalism for describing such a j T distribution. The plane perpendicular to the thrust axis splits the full phase space into two hemispheres. One only measures the hadron j T in one hemisphere, and the other hemisphere is unmeasured. Such type of measurements are termed as non-global observables [51], which are sensitive to radiation in only a part of phase space. The factorization and resummation formula for non-global observables have a very different structure from the standard global observable [52]. For example, the leading-order evolution equation for Non-Global Logarithms (NGLs) is named as BMS equation [53], which is a non-linear evolution equation. The TMD factorization formalism for the non-global hemisphere event shape has been studied by one of the authors in [54], where they find that the rapidity logarithms evolution does not constitute an essential complicated structure, since it is tied with a universal transverse momentum dependent jet function which also appears in the global observables. Besides, after comparing with the data at the LEP, they also find that the leading non-perturbative effects are related to the Collins-Soper kernel 2 .
We mainly consider the kinematic region with j T Q, where a TMD factorization can be developed which resums ln(Q/j T ). Here Q is the virtuality of the intermediate photon in e + e − → γ * . In this region, Belle collaboration finds that the cross sections can be well described by Gaussians in j T , and that the width of the Gaussians shows an initially rising, then decreasing z h -dependence when z h → 1. Because of this, we further consider the threshold ln(1 − z h ) resummation in the z h → 1 limit. We apply SCET to develop a TMD factorization formalism. Using renormalization group evolution techniques, we resum logarithmic terms to next-to-leading logarithmic (NLL) accuracy, including NGLs. The experimental data are shown as comparison and in good agreement with our theoretical predictions.
The remainder of this paper is organized as follows. In Sec. 2, we present a factorized framework, which only resums the so-called global logarithms. This section would allow us to develop intuition for our framework and understand connection to the standard TMD FFs. In Sec. 3, we present the full factorization formalism, which allows us to resum both global and NGLs. In Sec. 4, numerical results of differential cross sections for pion pro-duction in e + e − → π ± X are presented, as a function of energy fraction z h and transverse momentum j T . We also present the Gaussian width for the j T distribution as computed from our theoretical formalism, and compare them with the Belle experimental data. Finally, conclusions are given in Sec. 5.

TMD formalism: global structure
We consider the process, e + + e − → h (z h , j T ) + X, in e + e − annihilation. The center-ofmass (CM) energy of the e + e − collisions is given by s = Q 2 = (p e + + p e − ) 2 , and the hadron momentum fraction z h = 2p h · q/Q 2 = 2E h /Q is measured. In addition, the hadron's transverse momentum j T is measured with respect to the so-called thrust axisn, which maximizes the event-shape variable thrust T [55]: with the momenta p i of the particles measured in the e + e − CM frame. For convenience, we align the thrust axis to be along +z-direction, and define light-like vectors n µ = (1,n) = (1, 0, 0, 1) andn µ = (1, −n) = (1, 0, 0, −1). We expand any momentum p µ in the light cone frame as p µ = (p + , p − , p T ) with p + = n · p = p 0 − p z and p − =n · p = p 0 + p z . It is important to emphasize that even though we measure the hadron transverse momentum j T with respect to the thrust axis, our cross section is not differential in the thrust variable T . In other words, we consider the cross section for the hadron production, which is differential only in z h and j T : That is to say, the only purpose of the thrust measurement is to provide the thrust axisn and we sum over the entire thrust region 0.5 < T < 1. For the cross section that is further differential in the thrust variable T , see e.g. Refs. [56,57]. We find that such an observable in Eq. (2.2) has a better connection to the standard TMD FFs. The plane perpendicular to the thrust axis divides the full space into two hemispheres: the one on the right (along +z-side) is referred to as the right hemisphere, while the one on the left is the left hemisphere. Note that the observed hadrons are always measured in the right hemisphere, while no measurement is performed for the left hemisphere. Because of this, the kinematics in the left hemisphere are unconstrained, our observable in Eq. (2.2) is a non-global observable [58]. Such observables will involve non-global structures which can not be captured by the traditional exponential formula [52]. In this paper we will apply the jet effective theory [59,60] to derive the factorization and resummation formula. Since the full factorization structure is quite complicated which we save for the next section, in this section, we will for the moment ignore the NGLs that arise from the non-global structure, and write down a factorized formalism to resum the global logarithm and build our intuition. The black lines represent the energetic partons in the unmeasured left hemisphere, while the hadron is measured in the right hemisphere. Vertical dashed line represents a plane that is perpendicular to the thrust axis and that divides the space to left and right hemisphere. The red curves denote soft radiations from the energetic partons with the virtuality of j T . The blue lines in the left panel describe collinear radiations along the thrust axis, while the purple ones in the right panel give collinear-soft (c-soft) radiations .

TMD factorization formalism
We consider the kinematic region where the transverse momentum is small j T Q, and thus is sensitive to TMD physics. Setting the usual power expansion parameter λ = j T /Q, we find that the relevant momentum modes in this region are given by • collinear: p c ∼ Q(λ 2 , 1, λ) • soft: The different modes are illustrated in Fig. 1 (left). The hard modes encode energetic radiations in the left hemisphere: since the hadron is observed in the right hemisphere and has j T Q, any energetic radiation in the right hemisphere will lead to a large transverse momentum for the hadron and thus move the hadron out of the kinematic j T Q region; consequently such radiation is not allowed in the right hemisphere. On the other hand, soft and collinear modes have the same transverse momentum of j T , and thus both contribute to our observable. The difference is that collinear modes encode energetic radiations along the thrust axis, while soft modes describes large angle long wave radiations. Based on the mode analysis, the factorization formalism is given as where D h/i (z h , k T , µ, ν) is the usual TMD FF with k T the transverse momentum of the hadron h with respect to the fragmenting parton i. On the other hand, S i (λ T , µ, ν) is the soft function, with µ and ν renormalization and rapidity scales, respectively. The leading-order (LO) cross section is given by with α em the fine structure constant. Note that the factorization in Eq. (2.3) neglects the power corrections from the ratios j 2 T /Q 2 , which is small in the kinematic j T Q region we consider. Nevertheless, in the region of j T ∼ Q one can include such power corrections from the fixed-order calculations [61]. This is usually referred to as the Y -term in the CSS formalism [28,29]. Figure 2: Three configurations that contribute to the hard function: (a) virtual correction; (b) quark q is on the right hemisphere, while both anti-quarkq and gluon g are on the left hemisphere; (c) gluon g is on the right hemisphere, while both quark q and anti-quarkq are on the left hemisphere. Note that the observed hadron is on the right hemisphere.
It is important to emphasize that the above TMD formalism is already different from the earlier conjectures used in [49,50]. In particular, at leading power, our formalism depends on both quark and gluon TMD FFs, while the previous conjecture contains only quark TMD FFs. To convince that this has to be the case, the easiest way is to look at the Feynman diagram configurations that contribute to our observable at the next-to-leading order (NLO), from which we also derive the hard functions H i with i = q (q), g. At LO, we produce back-to-back quark q and anti-quarkq, each in their corresponding left or right hemisphere, and our hard function is normalized to be H = 1 at this order. At NLO, we receive three contributions as shown in Fig. 2. Here, Fig. 2 (a) is the virtual correction to the LO process e + e − → qq and q (orq) later on fragments into the observed hadron h, and thus this contribution is associated with the quark TMD FFs Fig. 2 (b) and (c) describe the hard scattering with three partons in the final state, where two hard partons are emitted in the left hemisphere and one parton i in the right hemisphere. Here, for three-particle final states the thrust axisn is determined by the direction of the most energetic parton. For (b), it isq and g on the left hemisphere, while q on the right hemisphere, which fragments into the hadron h and thus we have quark TMD FF D h/q . For (c), it is q andq on the left hemisphere, while g on the right hemisphere which fragments into the hadron h and thus we have gluon TMD FF D h/g in Eq. (2.3). We emphasize again that no hard radiations are allowed in the right hemisphere to maintain j T Q.
Direct calculations give us the following expressions for the corresponding bare hard functions at NLO, where we use the notation H i m with index m = 2, 3 at NLO representing the number of final-state partons, while the subscripts (a), (b), (c) correspond to the configurations in Fig. 2. We include the LO result into H i 2 , and we note that H g 3 starts at O(α s ) order, which is free of any divergence. Note that the function H q 2 is the standard dijet hard function, that arise in e.g. back-to-back hadron pair production [17]. If one ignores the non-global structure, the renormalization group (RG) equation for the hard function can be easily obtained from the above expressions. However, the structure for the full RG equations is much more complicated and will be shown in the next section.

TMD formalism in coordinate space
TMD formalism in Eq. (2.3) involves convolution over the transverse momentum k T and λ T . We apply the Fourier transform to go into the coordinate b-space and thus the convolution becomes a simple product. To get started, realizing and we thus can write Eq. (2.3) in the following form where the b-space TMD FF and soft function are defined as Both TMD FFs and soft function suffer from rapidity divergence, which was regularized via the rapidity regulator in [46,62]. As a consequence we have rapidity poles in 1/η and the associated rapidity scale ν, besides the usual poles in 1/ in dimensional regularization and the associated renormalization scale µ. In order to resum relevant logarithms, one can use transitional CSS formalism [28], or effective theory approaches [45,46,62] 3 . The NLO perturbative expressions for TMD FFs are well-known [23], and we list here for completeness: where the splitting functions are given by The NLO soft function S q (b, µ, ν) can also be computed easily. Since at NLO, only soft radiation that is emitted in the right hemisphere contributes to the hadron transverse momentum j T , this will put a constraint for the soft gluon momentum k in the soft function, i.e., k z > 0 or k − > k + .
It might be instructive to point out that the above soft function is exactly half of the standard soft function for the back-to-back hadron pair production in e + e − collisions, as well as those in SIDIS and Drell-Yan processes. This difference is precisely introduced by the constraint k z > 0 for the radiated soft gluon. This situation is similar to the case where one measures the transverse momentum of hadrons inside a jet with a jet radius R, as studied in [23], where soft functions in these two situations are related to each other by a boost along the z-direction.
With the explicit expressions for TMD FFs and soft function at NLO given above, one can easily obtain their corresponding µ and ν evolution equations: Here the relevant anomalous dimensions are given by (2.17) We have the first few coefficients given by , and n f represents the quark flavor number.
It is important to realize that the rapidity divergences between TMD FF D h/q (z h , b, µ, ν) and soft function S q (b, µ, ν) cancel between them. This is to be compared with the standard case, e.g., back-to-back hadron pair production in e + e − collisions, where the rapidity divergences cancel between one TMD FF D h/q (z h , b, µ, ν) and the square-root of the standard soft function S q (b, µ, ν), see e.g. [17,29,62,63]. Following the modern formulation of TMD FFs, we combine them as the so-called properly-defined TMD FFs [29,63] as follows: Using the evolution equations in Eq. (2.13), one can obtain the evolved TMD FFs D TMD to be at the hard scale µ h ∼ Q and thus resum the relevant logarithms ∼ ln(Q 2 /j 2 T ). For example, the standard exercise is to evolve S q from its characteristic scales µ s ∼ µ b and ν s ∼ µ b , and D h/q from its natural scales µ D ∼ µ b and ν D ∼ Q, to the hard scale µ h ∼ Q and a common rapidity scale ν, from which one obtains is the rapidity anomalous dimension [46,62] or Collins-Soper kernel [29,63]. Combine the above evolution equations, we thus obtain where we have On the other hand, the exponent of the evolution factor, i.e. the perturbative Sudakov factor S pert (µ b , µ h ) resums all the global logarithms and is given by Finally when the scale µ b Λ QCD , one can further match the TMD FFs D TMD h/q (z h , b, µ b ) onto the collinear FFs D h/q (z h , µ b ): where C i←q (z, µ b ) = δ iq δ(1 − z) at LO and the higher-order expressions can be found in [17,29,64,65]. On the other hand, when µ b ∼ Λ QCD , one has to introduce non-perturbative contributions, for which we apply the usual b * -prescription to include the TMD evolution in the large b region. Here we have b * defined as where b max is chosen [17] to be 1.5 GeV −1 . At the same time we include non-perturbative function S NP (b, Q 0 , Q), which is given by [23,36] with Q 2 0 = 2.4 GeV 2 , g 2 = 0.84 and g h = 0.042. We choose to work at the next-toleading logarithmic (NLL) level, we thus include two-loop cusp and one-loop normal anomalous dimension, and tree-level matching coefficients. Then plugging in the above results for D TMD h/q (z h , b, µ h ) in Eqs. (2.22) and (2.26), along with the non-perturbative function S NP (b, Q 0 , Q) in Eq. (2.28), into the differential cross section in Eq. (2.7), we obtain the all-order resummation formula where the Bessel function J 0 arises after integrating the angle between b and j T . We have chosen the following scales Such a formalism in Eq. (2.29) resums all the global logarithms in ln(Q 2 /j 2 T ).

TMD formalism at threshold z h → 1
Belle collaboration finds that the hadron cross sections can be well described by Gaussians in j T in the small j T region, and that the width of the Gaussians shows an initially rising for small to intermediate z h , while a decreasing z h -dependence for large z h 1. In the region z h → 1 region, the threshold logarithm of ln(1 − z h ) would become important and thus has to be resummed. In our phenomenological section, we find that the joint threshold and TMD resummation will be able to describe well such a z h -dependence for the Gaussian width. We develop theoretical formalism in this section for this purpose. As we will show below, in the threshold region, an additional mode, so-called collinear-soft (c-soft) mode [57,66,67] is relevant. Such a mode is shown as the purple curves in Fig. 1 (right), and the corresponding momentum scaling is given by Let us start our discussion with the fixed-order result of the perturbative TMD FFs D q/q and D g/q in Eq. (2.10) in the threshold limit. By taking the limit z h → 1, we find at NLO where we keep the overall factor of 1/z 2 h as a convention. Note that in this limit, one can drop the mixing term D g/q in comparison with the more singular terms δ(1 − z h ) and 1 (1−z h ) + in D q/q , i.e. only the flavor diagonal q → q channel contributes. In the threshold limit, we can refactorize the TMD FF D h/q as where S q is a collinear-soft (c-soft) function [57,66,67] that takes into account the soft radiation along the direction of the thrust axis, i.e. the c-soft mode mentioned above. At NLO, it can be computed as follows (2.33) Note that the c-soft function S q (z, b, µ, ν) has the same rapidity anomalous dimension as the TMD FF D h/q (z h , b, µ, ν), which is cancelled after combining soft function S q in Eq. (2.12) and the c-soft function in Eq. (2.33). On the other hand, we also have the collinear FFs at the threshold limit, whose perturbative expression is given by To perform the resummation in the threshold limit, one can perform the Mellin transform or Laplace transformation [68], whose purpose is to convert the above convolution in z-space into a simple product in the corresponding transformed space. Here we choose to perform the Laplace transformation [69], wherez h = 1 − z h . Using the following relation in the threshold limit Note that we have also extended the integration from 0 <z < 1 to 0 <z < ∞ in the threshold limit approximation. The NLO expressions forS q (κ, b, µ, ν) andD h/q (κ, µ) in the Laplace space are given bỹ (2.39) From the above results, one can derive the RG equations for bothS q and D q/q in the Laplace space, where the normal anomalous dimensions γ i expanded as γ i = n=1 γ i n−1 (α s /4π) n with i =S q , f q , and The above RG equations allow us to evolve c-soft functionS q (κ, b, µ, ν) from its natural scale µ S ∼ µ b and ν S ∼ κQ, and the FF D h/q (κ, µ) from initial scale µ F , up to the hard scale µ h and a rapidity scale ν, we obtaiñ Combining the evolution for the global soft function S q (b, µ, ν) in Eq. (2.21), we obtain the following evolution for the properly-defined TMD FFs in the Laplace space, .

(2.46)
Here the perturbative Sudakov factorS pert (µ b , µ h ) is given bỹ where the first integral represents the evolution of c-soft function and the global part of the soft function from µ b to µ h , and the second one is collinear fragmentation function from factorization scale µ F to µ h in the threshold limit. Performing the inverse Laplace transform, we obtain the following expression for TMD FFs in the threshold limit Here we derive the above formula using the first line of Sudakov factor in Eq. (2.47) and setting µ F = µ h , and the parameter η is defined as On the other hand,Ŝ pert (µ b , µ h ) in the momentum space in the threshold limit is given bŷ where the argument in the logarithm is given by (1 − z)Q instead of the usual Q in the threshold limit. Finally using the above result, one can obtain the resummed formalism for the differential cross section at the NLL level where we choose the non-perturbative Sudakov factorŜ NP (b, Q 0 , Q) to have the following formŜ Here motivated by the argument in the perturbative Sudakov function in Eq. (2.50), we replace Q by (1 − z h )Q in the usual non-perturbative function S NP (b, Q 0 , Q) in Eq. (2.28) to obtainŜ NP (b, Q 0 , Q) in the threshold limit.

Factorization and Resummation: full story
The plane perpendicular to the thrust axis divides the full space into two hemispheres. One measures the transverse momentum of hadron h in the right hemisphere, and the left hemisphere is inclusive. As we have emphasized above, hadron transverse momentum with respect to the thrust axis is a non-global observable, since the left hemisphere are unobserved. Such type of observables will involve non-global structures which can not be captured by the traditional exponential formula [52]. In this section, we apply formalism developed in [59,60] for the jet effective theory to derive the factorization and resummation formula. Such a formalism will enable us to resum both global and non-global logarithms. The global logarithmic structure has been discussed in the previous section. Here in this section, we pay more attention to the NGLs [58]. Very recently, a similar structure is also mentioned in [70].
In the standard TMD region where j T Q, following the development in [59,60], we can write the factorization formalism as follows dσ where H, S, and D h/i correspond to hard, soft and TMD FFs, respectively. Besides, different from the formalism in the previous section that resums only global logarithms, the hard and soft functions are now matrices in the color space, so we take color averaging as Tr c [· · · ]/N c after multiplying them and integrating out the solid angles {n} = {n 1 , n 2 , · · · } of the hard partons, where the angular integration is expressed by the symbol ⊗. The index m denotes the number of energetic partons inside the hard function that is defined in [54]. The index m in soft function then represents the number of Wilson lines, and the momentum space soft function is defined as Here X R denotes soft states in the right hemisphere, and one only measures the contributions from soft radiations in the right hemisphere. It is precisely because of such a multi-Wilson line structure that makes the hard and soft function matrices in the color space. After performing Fourier transformation for the observed j T , the factorization formula is given as For the NGLs resummation, we use the same methods in [59] to perform renormalization for the multi-Wilson-line operators. The renormalization constants of the hard and soft function are matrices in the color space, which are given as The expressions for the one-loop global anomalous dimensions have been given in the previous section. After solving the RG equations, we can obtain all-order resummation formula. At the NLL accuracy it has the form as In comparison with the resummed formalism in Eq. (2.29), we have the non-global evolution function U NG , which is given as where U NG is the evolution function for the non-global parts. At the LL accuracy and the large-N c limit, one can calculate it using the parton shower algorithms in [51,71] or the numerical solution of the BMS equations [53]. For the convenience of our numerical calculations in the next section, however, we choose the parametrization given in [51] with a = 0.85 C A , b = 0.86 C A , c = 1.33, and where β 0 = 11 3 C A − 4 3 T F n f , with T F = 1/2. For the differential cross section in the threshold limit, we find that the same nonglobal evolution function U NG arises. We thus write the resummed formalism at the NLL in the threshold limit z h → 1 as

Numerical results
In this section, we will study the differential cross sections and Gaussian widths of the transverse momentum distribution for the single inclusive charged pion production (sum of π + and π − ) in electron-positron annihilation process, e + e − → π ± + X, based on the factorization and NLL resummation formula given in Sec. 3. For parton-to-pion fragmentation functions, we use 2014 DSS analysis [72], where the uncertainties were determined Figure 3: Differential cross sections for the charged pion as a function of j T and different z h bins with ∆z h = 0.05 interval in the region 0.15 < z h < 0.6. In each plot, TMD resummation is applied and data points (black dots and histogram with error bars) from Belle collaboration [48] are shown for comparison. The error bands correspond to 90% C.L. uncertainty obtained from DSS fit [72].
based on the standard iterative Hessian method. Note that Belle data [48] was originally presented in different thrust bins, in 0.5 < T < 0.7, 0.7 < T < 0.8, 0.8 < T < 0.9, 0.9 < T < 0.95 and 0.95 < T < 1.0. Since the theoretical formalism we have developed in this paper is inclusive in the thrust variable, we thus combine the experimental data to obtain the results for the entire region 0.5 < T < 1.0. The data shown in this section are all the ones obtained via such a combination procedure. The errors of the data sets are also combined weighted by corresponding thrust bins. Fig. 3 shows the differential cross sections for pion production in e + e − collision as a function of the pion transverse momentum j T , in different z h bins at the center-of- Figure 4: Transverse momentum j T distribution given by TMD resummation (red band) and the joint TMD and threshold resummation (blue band) for the charged pion production with z h bins 0.65 < z h < 0.7, 0.7 < z h < 0.75, 0.75 < z h < 0.8 and 0.8 < z h < 0.85 from left to right. Results are shown in comparison with Belle data [48] in each pad. Error bands represent 90% C.L. uncertainty. mass energy √ s = 10.58 GeV. The error bands correspond to 90% confidence level (C.L.) uncertainty of parton-to-pion FFs determined in [72]. The hadron transverse momentum with respected to thrust axis are given in 0 < j T < 1.0 GeV for each plot. The energy fraction region 0.1 < z h < 0.65 is divided into eleven sub-regions with ∆z h = 0.05 for each panel. As seen clearly in the figure, for the intermediate z h region (z h 0.5), the evaluations based on TMD resummation in Eq. (3.15) are in good agreement with the data 4 . On the other hand, as z h becomes relatively large (z h 0.5) and thus approaches threshold limit, the agreement becomes worse, which indicates the potential importance of the threshold resummation effect.
In Fig. 4 we compare the differential cross sections obtained by using two resummations schemes: transverse momentum resummation (shown in red curves) and joint transverse momentum and threshold resummation (shown in blue curves). Similarly as before, the error bands correspond to 90% confidence level (C.L.) uncertainty of parton-to-pion FFs. The hadron transverse momentum with respected to thrust axis are given in 0 < j T < 1.0 GeV region. The energy fraction regions are 0.65 < z h < 0.7, 0.7 < z h < 0.75, 0.75 < z h < 0.8 and 0.8 < z h < 0.85 from left to right. In Fig. 4, z h bins are larger than those in Fig. 3 where the threshold logarithms are making some difference, thus compared to TMD resummation, we see that joint resummation has a better performance in these z h bins, especially in the small j T region. As z h gets larger, the consistency between joint resummation results and data gets better with a decreasing Gaussian width. The jointly resummed differential cross section decreases faster, indicating the a smaller Gaussian width value, which is more consistent with experimental data compared to the results with only transverse momentum resummed, where shapes are almost the same for the four z h bins in such a large z h region. To see the change of j T width as a function of z h , we fit the cross section dσ/dz h d 2 j T as a function of j 2 T , dσ and reconstruct the Gaussian width σ 2 j T for both theory and experimental data. We compute Gaussian width as a function of fractional energy z h using both TMD resummation (red curve) and joint resummation (blue curve). In Fig. 5, for small z h region (z h < 0.5), the logarithmically enhanced contribution origins from ln(Q/j T ), thus transverse momentum resummed cross section σ 2 j T fits the data well. As the value of z h is increased, for the TMD factorization theorem in Eq. (3.15), dependence on z h becomes weak, leading to a plateau at the tail region. On the other hand, for the factorization theorem with joint resummation in Eq. (3.19), where transverse momentum and threshold logarithms are jointly resummed, the cross section sharply decreases as z h increases, indicating a better fit for this region. Generally speaking, for kinematic regions distinguished by z h bins, adopting TMD resummation in intermediate z h regions while making use of joint resummation for large z h bins can lead to excellent agreement with measurement for e + e − → πX process, suggesting our factorization and resummation formula results in a reasonable approach for describing single inclusive hadron production at the electron-positron colliders.

Conclusion
Single inclusive hadron production at the e + e − colliders provide a new opportunity to study transverse momentum dependent fragmentation functions (TMD FFs), which are important to understand the 3D structure for the hadrons and the non-perturbative QCD. Belle collaboration has performed the first measurement for this observable, e + e − → h(z h , j T ) + X, where z h is the energy fraction for the hadron, while the hadron's transverse momentum j T is measured with respect to the thrust axis determined by the hadronic event shape. We develop a TMD factorization formalism for such an observable, which resums logarithm of ln(Q/j T ). Realizing the non-global nature of the observable, our factorization formalism involves a multi-Wilson line structure, which allows us to resum both global and non-global logarithms. Besides, as the increasing of the energy fraction z h of the hadron, the threshold soft gluon enhancement effects become more and more important, which require us to perform joint TMD (∼ ln(Q/j T )) and threshold (∼ ln(1 − z h )) resummation. We apply the formalism proposed in [57] based on SCET+ framework [66] to obtain factorization and resummation formula in the joint limit.
In the end we find that TMD resummation formula give a good description for the j T distribution as z h < 0.65. For large z h > 0.65 region, in order to describe the data we need to include threshold resummation effects. Especially, we find that the Gaussian width of the j T distribution given by the TMD formalism freeze to a certain value which is not consistent with the measurement. While after including joint threshold and TMD resummation effects, the theoretical predictions are consistent with data very well.
In the present work we obtained the perturbative resummed cross section at the nextto-leading logarithmic (NLL) accuracy. In the future work, we will include higher order resummation effects using method developed in [73]. Especially, in this case beyond the NLL level, the gluon TMD FF will also contribute to the cross section as shown in (3.3), it will be interesting to study its effects.