Factorization of $e^+e^- \to H \, X$ cross section, differential in $z_h$, $P_T$ and thrust, in the $2$-jet limit

Factorizing the cross section for single hadron production in $e^+e^-$ annihilations is a highly non trivial task when the transverse momentum of the outgoing hadron with respect to the thrust axis is taken into account. We work in a scheme that allows to factorize the $e^+e^- \to HX$ cross section as a convolution of a calculable hard coefficient and a Transverse Momentum Dependent (TMD) fragmentation function. The result, differential in $z_h$, $P_T$ and thrust, will be given to all orders in perturbation theory and explicitly computed to Next to Leading Order (NLO) and Next to Leading Log (NLL) accuracy. The predictions obtained from our computation, applying the simplest and most natural ansatz to model the non-perturbative part of the TMD, are in exceptional agreement with the experimental measurements of the BELLE Collaboration. The factorization scheme we propose relates the TMD parton densities defined in 1-hadron and 2-hadron processes, restoring the possibility to perform global phenomenological studies of TMD physics including experimental data from semi-inclusive deep inelastic scattering, Drell-Yan processes, $e^+e^- \to H_1 H_2 X$ and $e^+e^- \to HX$ annihilations.


Introduction
Recently, the BELLE collaboration has measured e + e − → H X cross section, as functions of P T , the transverse momentum of the detected hadron with respect to the thrust axis [1] at a c.m. energy of about 10 GeV. This can be considered a break-through measurement in the investigation of the 3D structure of hadrons, as it offers a direct glance to the transverse motion of the fragmenting partons. Unfortunately the factorization mechanism, which was first devised by Collins, Soper and Sterman (CSS) for Drell-Yan processes and later proven to work for Semi-Inclusive Deep Inelastic Scattering (SIDIS) [2][3][4][5][6] and e + e − → H 1 H 2 X Ref. [7], cannot be directly applied to e + e − → H X processes. In fact, having only one hadron detected in the final state makes it impossible to cast the cross section in a form that allows to define the Transverse Momentum Dependent parton densities (TMDs) in the conventional way, i.e. by including part of the soft radiation generated in the process inside the TMD itself. This is a serious impediment which endangers the possibility to exploit the valuable information encoded in the BELLE experimental data in the phenomenological study of the mechanisms which lead to hadronization and, ultimately, to the structure of hadrons in terms of their elementary constituents, quarks and gluons.
In the attempt to overcome this obstacle, in Ref. [8] we introduced a factorization scheme, closely following the CSS methodology of Ref. [7], in which the definition of the TMD is slightly modified in such a way to leave out all its soft content. As a result, the final cross section is written as the convolution of a hard coefficient, which represents the partonic core of the process and is calculable in pQCD at any desired order, with a Transverse Momentum Dependent Fragmentation Function (TMD FF), which carries non-perturbative information and provides a direct probe of the 3D-dynamics of the hadronization mechanism. The advantage of this scheme is that it allows to relate the TMD fragmentation function extracted in e + e − → H X, which is a 1-hadron class process, to the TMD parton densities extracted from SIDIS, Drell-Yan and e + e − → H 1 H 2 X, i.e. in any 2-hadron class process. This is the main point of strength of the proposed scheme: the immense effort already dedicated in the extraction of TMD distribution and fragmentation functions in the past two decades does not go wasted, as the relation to those function is uniquely determined and preserved. This point will be discussed in more detail in Section 4 where we will present a basic phenomenological application of our final result.
These issues are also treated in a in the context of Soft and Collinear Effective Field Theories (SCET) , see for instance Refs. [9][10][11]. Very recently two papers, Refs. [12] and [13], have been dedicated to the study of e + e − → H X cross sections as derived in a SCET framework.
This article is structured as follows. Section 2 shows the computation of the hard part of the process, i.e. the partonic cross section: two separate subsections are dedicated to the cases in which the fragmenting parton is either a gluon or a fermion, respectively. The computation at Next Lowest Order (NLO) is performed explicitly. In Section 3 we revise the structure of the TMD fragmentation function and its general perturbative and non perturbative content, giving the explicit expressions of all relevant quantities to Next to Leading Log (NLL) accuracy. In Section 4 we build the final expression of the e + e − → HX cross section, we discuss its properties and comment on issues related to its resummation. In Section 5 we show how our final formula can be applied to reproduce the experimental cross section as measured by the BELLE Collaboration. Even without a proper fit the agreement with data turns out to be excellent, confirming the validity and strength of the formalism we propose. In Section 6 we draw our conclusions.

Kinematics of the e + e − → H Xprocess and cross section
In the process measured by the BELLE collaboration an electron and a positron, with momenta l 1 and l 2 respectively, annihilate in a virtual boson with momentum q (since the c.m. energy is around 10 GeV, it can only be a photon). The observed hadron, of momentum P , belongs to a jet initiated by a parton produced in the e + e − annihilation. The BELLE cross section is sensitive to the fractional energy z h of the detected hadron H with respect to the total energy Q = q 2 available in the c.m.: The parton frame is defined by fixing the z-axis of the c.m. frame to coincide with the axis of the jet, which is also the direction of the fragmenting parton. Then, the momentum P of the detected hadron can be written as: P = P + , P − p , P p, T with P + = z h q + , (1.2) and the measurement of P p, T probes the transverse motion of the partons inside the detected hadron. However, the determination of the jet axis is not easy. In what follows, we will identify it with the reconstructed thrust axis n, which is the direction that maximizes thrust, T , defined as: where the sum is over all the final state particles i. The variable T describes the topology of the final state and it ranges from 0.5 to 1, where the lower limit corresponds to a spherical distribution of final state particles, while the upper limit realizes a pencil-like event. Therefore, in this paper we will interpret the transverse momentum measured by BELLE as P T ≡ | P p, T | and we will provide a cross section differential in the three variables: z h , P T and T . Following the scheme presented in Ref. [8], the final cross section will be written schematically as: where M H is the mass of the detected hadron and the integration variable z is the light-cone momentum fraction of the detected hadron with respect to its parent (fragmenting) parton. Kinematics constrain z to lie in the range z h ≤ z ≤ 1. We will focus on spinless hadrons, hence only the unpolarized TMD FFs D 1 appears in the cross section. The TMD FFs are defined in the Fourier conjugate space of the transverse momentum of the fragmenting parton with respect to the direction of the detected hadron. A proper definition of this transverse momentum is obtained by boosting (or equivalently by rotating) the parton frame defined in Eq. (1.2) in order to re-define the z-axis so that it coincides to the direction of the detected hadron. According to Ref. [7], we will refer to the resulting frame as the hadron frame, where the momentum k of the fragmenting parton can be written as: The transverse momenta k h, T and P p, T are related by the following expression: Therefore, in Eq. (1.4) the unpolarized TMD FFs have to be intended as: Notice that the Fourier transform acts like an analytic continuation of the TMD in momentum space and consequently it is totally inadequate when used in the range of large values of P T , where TMDs lose their physical meaning.
The cross section in Eq. (1.4) can also be written as the contraction of two tensors: the leptonic tensor L µν and the hadronic tensor W µν H , which encode the initial and final state information, respectively. The leptonic tensor is defined as the lowest order of the electromagnetic vertex e + e − → γ and it is given by: where l 1 and l 2 are the momenta of the electron and the positron, respectively, and θ is the zenith angle of the electron with respect to the thrust axis n: (1, u) ; (1.9) with u · n = cos θ. The final cross section is insensitive to the value of θ, hence we will integrate over the whole solid angle specified by the leptons. On the other hand, the hadronic tensor is defined by: W µ ν H (z h , T, P T ) = 4π 3 X δ (4) (p X + P − q) × × 0| j µ (0)| P, X, out T T P, X, out|j ν (0)| 0 = = 1 4π X d 4 z e iq·z 0| j µ (z/2) | P, X, out T T P, X, out|j ν (−z/2) | 0 , (1.11) where the factor 1/(4π) coincides with the normalization factor used in Ref. [7], while the final states vectors are labeled by "T" as a reminder of the event topology, which has to be selected according to the value of thrust. Finally, the cross section is given by: where φ is the azimuthal angle of the electron with respect to the x-axis of the c.m. frame. Since we are considering unpolarized leptons, the cross section does not depend on this variable. For practical purposes, the hadronic tensor can be written by using the structure functions obtained by projecting onto its relevant Lorentz structures: We can easily compute the projections: (1.14) (1. 15) Such decomposition allows to define the transverse (T) and the longitudinal (L) component of the cross section as: where we used the Born cross section: This remarkably simplifies the integration over θ, in fact we have: (1.19) Therefore it follows straightforwardly: Another simplification occurs when the projection with respect to P µ P ν (see Eq. (1.15)) is zero (or can be neglected), making the longitudinal cross section to vanish. In this case, the two structure functions are not independent anymore and in fact F 2, H = − 2 z h F 1, H . As a consequence, the hadronic tensor can be written as: where we defined the transverse tensor:

Partonic Cross Section
In this section, we will compute the partonic cross section of e + e − → H X in the 2-jet limit, providing its explicit expression to 1-loop (NLO) precision and also its all-order, resummed, formulation. According to the scheme presented in Ref. [8], the partonic cross section is the hard part of the factorized full cross section and encodes its short-distance contributions. Therefore, it is fully computable in perturbative QCD as it represents the partonic analogue of the whole process. In other words, the partonic cross section describes the process e + e − → j X, where j indicates the parton of type j that fragments into the detected hadron H; j can either be a gluon or a fermion (quark/antiquark) of flavor f . However, in the 2-jet limit, the contribution of the fragmenting gluon is strongly suppressed by kinematics (we will present the explicit computation in Sect. 2.1) so that the only relevant contribution is given by the fragmenting fermions. The topology corresponding to a quasi 2-jet configuration is obtained by requiring the thrust T to be almost 1 or, analogously, τ = 1 − T ∼ 0. Since the lowest order is an exact 2-jet configuration (pencillike), the first genuinely non trivial term is generated at 1-loop, and is due to the emission of a real parton crossing the final state cut. Practically speaking, we will introduce a sharp cut-off τ MAX that restricts the range of τ values to 0 ≤ τ ≤ τ MAX . Then, the 2-jet-like topology for the final state will correspond to the limit τ MAX → 0.
The reference frame for the computation is the partonic analogue of the hadron frame, see Eq. (1.5), where the momentum k 1 of the outgoing parton lies along the plus direction: All the formulas valid for the whole cross section (see Section 1) naturally extend to their partonic version. Clearly the partonic cross section cannot merely be obtained by summing over all the relevant Feynman graphs in the small-τ limit. In fact, in the final result, the overlapping between the hard and the collinear momentum regions has to be appropriately subtracted, while renormalization is achieved through an Ultra-Violet (UV) counterterm. According to Ref. [8], the subtracted and renormalized final state tensor of the partonic cross section is given, order by order in perturbation theory, by: here is the dimensional regularization parameter, where the dimension of space-time is D = 4 − 2 . In the following, we will separately compute the contributions of a fragmenting gluon and a fragmenting fermion, in Sections 2.1 and 2.2 respectively. Figure 1. The only 1-loop Feynman graph contributing to W µ ν g , when the gluon is emitted by the quark. The emission from the antiquark line is analogous.

Fragmenting Gluon
In a 1-loop computation, the configuration in which the detected hadron is produced by the fragmentation of a gluon can only occur through the emission of a real gluon from either the quark or the antiquark line. Then, to compute the final state tensor W µ ν g one simply considers the Feynman graph in Fig. 1. If the 2-jet limit is applied to the final state topology, then W µ ν g is suppressed. In fact, after the emission of the gluon, the fermion cannot deviate drastically from its original direction (otherwise it would generate a third jet) and hence it can only proceed almost collinearly to the gluon. In principle, a 2-jet configuration may be achieved also if the fermion turns soft or if it reflects backwards after the emission of the gluon. However, such configurations are suppressed by power counting (see Chapter 5 of Ref. [7]). As a consequence, the only relevant kinematic configuration in the 2-jet limit is given by the emitting fermion being collinear to the fragmenting gluon. However, this is exactly the same configuration that has to be subtracted out in the final result, in order to avoid double counting due to the overlapping with the collinear momentum region. This simple argument is enough to conclude that the subtracted final state tensor associated to the fragmenting gluon has to be suppressed in a quasi 2-jet configuration. Nevertheless, we will perform the explicit computation of W µ ν g anyway, for pedagogical reasons: it is very simple compared to the more important case of the fragmenting fermion (only one Feynman diagram, no rapidity cut-offs, etc ...) and it will serve as a benchmark to introduce the main features associated with the computation of the final state tensor. Section 2.1.1 will present the explicit computation of the partonic cross section, while in Section 2.1.2 we will describe the mechanism of subtraction.

Unsubtracted Final State Tensor
The unsubtracted final state tensor is obtained by applying the hard approximation (T H in Ref. [7]) to the Feynman graphs contributing to the desired order in perturbation theory. This implies considering their massless limit and replacing every soft-collinear divergent quantity with its lowest order approximation. In momentum space, the lowest order partonic TMD is simply a delta function, both in z and in the transverse momentum of the fragmenting parton k T , therefore in the Fourier conjugate space (where we will work) the unsubtracted final state tensor does not depend on b T (see Section 6 of Ref. [8]).
In the following, we will assume that the fragmenting gluon is emitted by the quark, as in Fig. 1. The case of emission by an antiquark is perfectly analogous. The final state hosts three particles: the fragmenting gluon, the quark and the antiquark. The squared amplitude and the phase space will depend on all the possible combination of scalar products among their momenta. Then, by labeling k 2 the momentum of the antiquark and k 3 that of the quark, it is extremely useful to express all the quantities in terms of the following variables: subjected to the constraint i y i = 1, due to the momentum conservation q = k 1 + k 2 + k 3 . We will use dimensional regularization, hence the dimension of spacetime is D = 4 − 2 .
According to standard conventions (Ref. [7]) the polarization vector e(k 1 , λ) of the on-shell fragmenting gluon is defined to have zero plus and minus components 1 . Then, the squared amplitude associated to the Feynman graph in Fig. 1, summed over all the polarization values λ, is given by: where e f is the charge fraction of the (anti)quark generated by the virtual photon. The third line of the previous equation has been obtained by using λ e µ (k 1 , λ)e ν (k 1 , λ) = −g µ ν .
The denominator can easily be written in terms of the variables y i defined in Eq. (2.3): On the other hand, we can compute the projections of the numerator onto its relevant Lorentz structures as in Eqs. (1.14) and (1.15): where the last line of Eq. (2.7) has been obtained by using transverse momentum conservation, k 2, j = −k 3, j . Furthermore: Therefore, the projection of the full squared amplitude onto the metric tensor gives: (2.10) The last line has been obtained by expressing the ratio k 2 3, T /Q 2 in terms of the y i variables and by defining: where H 0, f is the constant factor appearing in front of the lowest order final state tensor, see Eq. (A.5). Notice that the expression in Eq. (2.10) only diverges if y 2 = 0, which corresponds to the configuration in which the quark is collinear to the fragmenting gluon. On the other hand, in y 1 = 0 the function is regular, in agreement with the power counting prediction about the soft and the collinear-to-antiquark configurations. The other projection Figure 2. The phase space available for the three final state massless particles g, q andq, with momenta k 1 , k 3 and k 2 , respectively. The dashed (red) bands represent the quasi 2-jet configurations, where y 1 and/or y 2 are zero. The sub-regions R i correspond to a value of thrust given by τ = y i . gives: (2.13) Notice that this expression is regular in both y 1 = 0 and y 2 = 0, consequently it is suppressed in the 2-jet limit.
The phase space available for the three final state massless particles can be represented by the triangle in Fig. 2. The phase space cannot extend beyond the edge given by y 2 = 1 − y 1 (or y 3 = 0), due to momentum conservation. Since the fragmenting gluon does not cross the final state cut, the phase space integral only involves k 2 and k 3 , and it reduces to a single integration over y 2 after applying the momentum conservation delta. It is given by: where S is a shorthand notation for (4π) /Γ(1 − ), as in Ref. [7]. The (4π) factor in front of the phase space integral simplifies with the normalization chosen for the hadronic tensor, see Eq. (1.11). Since the phase space integration involves only y 2 , it can be very useful to change variables and write all the quantities in terms of the following: Notice that here we are not interested to the whole phase space, but only to the (very narrow) region corresponding to a quasi 2-jet configuration, where y 1 and/or y 2 tend to zero. This region is indicated by a dashed red band in Fig. 2. The integration over these regions is well defined once we appropriately limit the range of thrust, T = 1 − τ , within the interval 0 ≤ τ ≤ τ MAX . At the partonic level, τ can be computed explicitly and its value changes inside the triangle of Fig. 2. In particular, τ always coincides with the minimum among y 1 , y 2 and y 3 , allowing the definition of three sub-regions inside the triangle, R 1 , R 2 and R 3 , according to the value assumed by τ . Then, the 2-jet limit of each region is obtained by considering the limit of vanishing τ MAX and keeping only the leading divergence in τ = 0. Notice that a pure 3-jet configuration is the point of intersection of all three sub-regions, where y 1 = y 2 = 1/3.
The integration over R 1 can be derived from Eq. (2.14) by letting y 2 vary between y 1 and 1 − 2y 1 , and by imposing y 1 < 1/3. Furthermore, in this case, τ = 1 − y 1 ≤ τ MAX . Notice that in the 2-jet limit the sub-region R 1 approaches the configurations where the quark is either soft or collinear to the antiquark. Then, its contribution must be suppressed by some power of τ . In the following we will verify this power counting prediction. Let us consider the R 1 -contribution to the projections of the unsubtracted final state tensor in the 2-jet limit. In terms of the set of variables defined in Eq. (2.15) we have: having defined the integrals: Since τ MAX is the maximum value that τ can reach, beyond which the final state goes into a 3-jet configuration, we surely can consider τ MAX < 1/3 and drop the first theta function in Eq. (2.16). Furthermore, since τ is limited to be small, in the interval [0, τ MAX ], we can take the following approximation, valid when τ MAX → 0: with Re < min(−a + 1, −b + 1). (2.18) which in our case means: Finally, the upper limit on τ translates into a lower limit on z, due to the Dirac delta. Hence we have 1 − τ MAX ≤ z ≤ 1. In the limit of vanishing cut-off, the contribution of the sub-region R 1 has to be proportional to δ(1 − z). In fact, the integration with a test function T (z) gives: therefore: Finally, we have: The above equation shows that the whole 2-jet contribution of the sub-region R 1 is of order O(τ − ). Since the limit τ → 0 has to be taken before expanding in powers of and since Re < 0, then the whole expression of Eq. 2.22 is suppressed in the 2-jet approximation, as expected from power counting. The other projection of the final state tensor gives an analogous result: Notice that the previous expression is regular for any and vanishes if Re < 0 when τ goes to zero. In conclusion, R 1 does not contribute to the unsubtracted final state tensor in the 2-jet limit.
Let's now consider the contribution of R 2 . In the 2-jet limit, this region approaches the configuration where the quark is either soft or collinear to the fragmenting gluon. While the soft configuration has to be suppressed by some power of τ , the contribution of the quark collinear to the gluon has to become larger and larger as τ goes to zero. Starting from Eq. (2.14), if y 1 > 1/3, then 0 < y 2 < (1 − y 1 )/2, otherwise if y 1 < 1/3, then 0 < y 2 < y 1 . In terms of the variables defined in Eq. (2.15) we have: Notice that the approximation in Eq. (2.26) forces z to be very close to 1 as τ MAX goes to zero, resulting in a configuration in which the quark is soft. The factor (1 − z) − suppresses its contribution as Re < 0, confirming the power counting prediction. Finally: where J [1] g/q is the 1-loop jet thrust function associated to the fragmenting gluon, divergent in τ = 0 and defined in Eq. (B.19). The other projection gives a contribution of O(τ − ) and hence it is suppressed in the 2-jet limit as Re < 0.
The final contribution comes from the R 3 region. Fig. 2 shows that this region reaches the 2-jet configuration (red dashed bands) only in the wedges very close to the axes, hence it is suppressed compared to R 1 and R 2 in the 2-jet limit. The projections give: Terms as (1 − τ z) −n− can be expanded in powers of τ , therefore the expression in the previous equation is suppressed in the limit of small τ . The same happens for the other projection. Therefore, we can neglect the whole contribution of this region. Finally, the only non zero contribution to the unsubtracted final state tensor associated to the fragmenting gluon in the 2-jet limit is given by region R 2 in Eq. (2.27): where the partonic structure function F [1] g is: As expected, W µ ν g ( F g ) is proportional to the thrust jet function of the fragmenting gluon, equipped with a cut-off τ MAX that constrains it to a 2-jet-like final state topology. The expression in Eq. (2.29) has been computed to 1-loop accuracy, but the computation could be repeated in a perfectly analogous way at any order in perturbation theory.
The effect of the combination of the sharp cut-off τ MAX and the function J g/q can be obtained equivalently by allowing for a slight modification of the original definition of the thrust jet function of the fragmenting gluon given in Eq. (B.19). Such definition involves a delta function that relates the value of τ to the transverse momentum k T of the fermion (e.g. the quark) that flows collinearly to the fragmenting gluon. If this transverse momentum is small (compared to Q), T is close to 1 ensuring a 2-jet-like final state topology. However, the original definition also involves an integration over k T that stretches the jet function well beyond the region where it is defined 2 . Therefore, if k T is not allowed to become too large, then the function J g/q is strictly limited to describe the 2-jet limit of the final state. This can be achieved by introducing a cut-off that constrains the integration range to the quasi 2-jet configuration by setting an upper limit for the transverse momentum k T . Indeed, the new definition of J g/q will have to coincide to its original definition, Eq. (B.19), at small values of τ . A proper cut-off for k T is the power counting energy scale λ << Q, used to size the collinear momenta. With this choice, we have: Actually, the integration does not extend to the whole spectrum of transverse momenta, since it must satisfy the implicit condition τ ≤ 1, i.e. T > 0. The maximum value of kT is of order Q and, to 1-loop, it is k We can relate this expression with that appearing in Eq. (2.29). In fact: Notice that this choice sets τ MAX = 1 only in a thin slice of the phase space, where z is very close to 1. Furthermore, τ MAX → 0 consistently implies λ 2 /Q 2 → 0 (as showed in Fig. 3) and we can rewrite Eqs. (2.29) and (2.30) as follows: and: The -expansion of the previous expressions is not straightforward. The thrust jet function for the fragmenting gluon as defined in Section B.3 can easily be expanded in powers of by using Eq. (B.8), however in this case the presence of an explicit cut-off makes the computation a bit trickier. In fact, a direct application of Eq. (B.8) gives: The interplay between the plus distribution and the theta function requires some extra care, especially in the limit of λ 2 /Q 2 → 0. Being a distribution, it is best studied by considering its action on a test function T (τ ). In order to make the following expressions more readable, we define r = λ/Q. Then, by definition of plus distribution and by using Eq. (2.32), we have: The last line (term in square brackets) of the previous expression is suppressed in the limit of r → 0. In fact, both integrals are well defined thanks to the plus distribution prescription; this allows to perform an expansion in powers of r 2 , where the first non-zero term appears at order O r 2 . Therefore it can be neglected, leading to the final result: Integration of the above expression over τ leads to the same result obtained by integrating the non-expanded initial function, apart from power suppressed corrections. In other words, by using Eq. (2.37) the operations of integration and the operation of -expansion commute. Then, we can write the thrust jet function of the fragmenting gluon, equipped with the cut-off as in Eq. (2.31), as follows: The previous expression can be further simplified, in fact the combination of the logarithmic terms can be rewritten as: This result follows by considering the combination of logarithms in the first line of the previous expression as a distribution of z and then integrating it with a test function T (z). Since all the integrals are well defined, they can be expanded in powers of λ 2 /Q 2 to give the final result of Eq. (2.39). Then, finally:

Subtraction Mechanism
The unsubtracted final state tensor found in the previous section describes the fragmentation of a gluon resulting from a e + e − scattering at partonic level. However, this information is also encoded in the TMD FF of the gluon, which appears in the final cross section convoluted with the partonic cross section obtained in Eq. (2.29). Therefore, W µ ν g [1] must be appropriately subtracted in order to remove all contributions that overlap with the partonic version of the TMD FF of the gluon, in order to avoid double counting. This subtraction will be performed in the b T -space, where TMDs are defined explicitly in terms of operators, see e.g. Refs. [7,8].
In momentum space, the partonic version of the 1-loop gluon-from-quark TMD FF is given by (see e.g. Ref. [7]): In principle, its Fourier transform involves an integral over the whole spectrum of k T transverse momenta, from zero to infinity, going far beyond the actual region of overlapping.
However, the presence of the cut-off λ ensures that the Fourier transform of the gluon TMD FF, Eq. (2.41), will overlap with the final state tensor only up to λ, matching the same range of k T that is involved in the thrust jet function of the fragmenting gluon equipped with the cut-off (see Eq. (2.31)). The incomplete Fourier transform of Eq. (2.41) is computed with the help of the following expression: Then, we define the Fourier transform of the gluon-from-quark partonic TMD FF, equipped with the cut-off λ, as obtained from the incomplete Fourier transform. This gives: The order-by-order formula for the subtracted, renormalized partonic cross section is given in Eq. (2.2). Notice that the pole in Eq. (2.43) embodies the collinear divergence associated to the TMD FF. In fact, in the case of a fragmenting gluon, the TMD FF is not UV divergent and hence no UV counterterm has to be added to Eq. (2.43) in order to obtain a renormalized quantity. Therefore, the final expression for the partonic cross section has to be subtracted but not renormalized. Its expression follows from the 1-loop version of Eq. (2.2):  Then, ultimately the subtracted partonic cross section for the case of a fragmenting gluon is suppressed in the 2-jet limit, i.e. when λ 2 /Q 2 → 0: Notice that, being λ the IR energy scale that sizes the almost on-shell collinear momenta, our initial prediction about the power counting suppression of the fragmentation of a gluon is totally confirmed by Eq. (2.46). This also supports the intuition about a gluon-initiated jet in a 2-jet-like final state: in such a topology, one jet is initiated by the quark, the other by the antiquark. Finally, in Eq. (2.46) we can relate λ directly to the measured value of thrust τ meas. . In fact, according to the approximation of thrust in the 2-jet limit, the measured thrust τ meas. is well approximated by the sum of the invariant masses of the two jets (see Eq. (B.1)). According to power counting, such invariant masses are of order λ 2 . Therefore, we simply have O(τ meas. ) = O λ 2 /Q 2 . This rather naive argument can be made more specific by exploting the definition of thrust (Eq. (1.3)) and the relation between the transverse momenta of the fragmenting parton and the detected hadron (Eq. (1.6)). In fact, it is not difficult to prove that (see Ref. [13]): Therefore, the choice λ 2 = τ meas. Q 2 is supported by the kinematical bounds of the process in a total natural way. Hence we can interpret the suppression of the partonic cross section in Eq. (2.46) as due to the topology of the final state of the e + e − scattering.

Fragmenting Fermion
In this case, the detected hadron is produced by the fragmentation of a fermion of flavor f , which we will assume to be a quark. The case of a fragmenting antiquark is totally analogous. At 1-loop, one of the two fermionic legs emits a gluon, which can be either virtual or real, as represented by the Feynman graphs in Fig. 4. If the emitted gluon is virtual, then the final state topology is a perfect pencil-like event, as at LO (see Appendix A). Instead, if the gluon is real, a 2-jet-like configuration can be obtained only in the following three situations: • The gluon is soft. Then the final state tensor is dominated by the soft thrust function S as defined in Appendix B.1.
• The gluon is collinear to the fermion which does not fragment (the antiquark in our case). Then we expect the final state tensor to be dominated by the backward thrust function J B as defined in Appendix B.2.
• The gluon is collinear to the fragmenting fermion. In this configuration the largest contribution to the final state tensor is given by the thrust jet function corresponding to the fragmenting quark J q/q , as defined in Appendix B.3. Notice, however, that this coincides to the configuration that has to be subtracted out in the final result, hence its computation has to be performed carefully.
All previous possibilities are allowed by power counting and all of them are expected to be dominant in the 2-jet limit. Note that the case of fragmenting fermion is way more complicated then the case of fragmenting gluon treated in Section. 2.1, not only because it implies more squared matrix elements, but mostly due to the fact that the subtraction mechanism is made considerably more difficult by the presence of a rapidity cut-off in the partonic quark-from-quark TMD FF. Moreover, the subtracted final state tensor will require the addition of a proper UV counterterm that renormalized its UV divergences. We will deal with all those issues in the following sections. In particular, in Section 2.2.1 we will present the explicit computation of the unsubtracted final state tensor W µν f , while in Section 2.2.2 we will perform the subtractions that lead to the final expression for the NLO partonic cross section. The resummed cross section will be presented in Section 2.2.3.

Unsubtracted Final State Tensor
Similarly to the case of the fragmenting gluon, the unsubtracted final state tensor is obtained by applying the hard approximation (T H in Ref. [7]), which sets all the masses to zero and all the soft-collinear divergent quantities to their lowest order. According to the discussion at the beginning of Section 2.1.1, we will work in the Fourier conjugate space.
Virtual Emission. When the emitted gluon is virtual, the final state hosts two particles: the outgoing quark, of flavor f and momentum k 1 , and the antiquark crossing the final state cut, of momentum k 2 . Momentum conservation sets q = k 1 + k 2 . The 1-loop squared amplitude is given by: This expression can be properly simplified by decomposing the Dirac structure in its scalar, vector and tensor parts, by using momentum conservation and the Passarino-Veltman reduction formula [14]. This leads to: which simply asserts that the 1-loop squared matrix element for the virtual emission of a gluon is the lowest order M [0] f , computed in Eq. (A.2), "dressed" with the vertex factor V . As a consequence, the corresponding contribution to the final state tensor will be simply proportional to the lowest order, computed in Eq. (A.5). Hence: The 1-loop vertex factor is given by: where we introduced the integrals I 0 and I 0 defined as: (2.53) Therefore: Real Emission. In this case, the emitted gluon is real and there are in total three particles in the final state: the fragmenting quark, the gluon and the antiquark. The squared amplitudes depend on all possible combinations of scalar products of the final state particle momenta, encoded in the variables y 1 , y 2 and y 3 defined in Eq. (2.3). Momentum conservation ensures that i y i = 1. The squared amplitudes in which the gluon crosses diagonally the final state cut will be labeled by "diag.". It is given by: The projections give: where H 0, f has been defined in Eq. (2.11). Then we have to compute the squared amplitude in which the gluon attaches to both the upper fermionic lines. It will be labeled by "up" and it gives: Its projections are: Notice that in this case there is a non vanishing contribution from the projection with k 1, µ k 1, ν . However, this is not divergent neither when y 1 → 0, nor when y 2 → 0 and hence we expect that it will be suppressed in a 2-jet-like final state configuration. The last squared amplitude involves the gluon attached to both the lower fermion lines and hence it will be labeled as "down". It is given by: f, R , whose projections are: As for the case of the fragmenting gluon, the available phase space for the three final state particles is represented by the dashed (red) vertical and horizontal boundaries of the triangle in Fig. 2. It is given by Eq. (2.14), where the delta function sets the values of thrust according to the sub-regions R 1 , R 2 and R 3 , while the theta function forces τ to lie in the neighborhood of zero [0, τ MAX ]. Furthermore, we will use the same change of variables adopted in Section 2.1.1, defined in Eq. (2.15).
The integration over the sub-region R 1 is totally analogous to Eq. (2.16). The 2-jet limit of this sub-region corresponds to the left, vertical edge of the triangle in Fig. 2, which coincides with having the gluon collinear to the antiquark except in the vertex, where the gluon turns soft. Therefore, it is reasonable to expect that the result of the integration over R 1 will be related to the backward thrust function J B and to the soft function S, defined in Eqs. (B.13) and (B.7), respectively. It is given by: (2.67) As expected, the 2-jet limit of sub-region R 1 is dominated by the backward thrust function J B and by the soft thrust function S. Notice that only half of the soft function appears in the final result, since "half" of the vertex is shared with sub-region R 2 . The other projection (with respect to k 1,µ k 1,ν ) is suppressed in the 2-jet limit, since it is proportional to the integral τ − I 0, 1 ∼ O(τ − ) that vanishes if Re < 0. Let's now consider the contribution of sub-region R 2 . Its contribution to the 2-jet limit of the process corresponds to the lower edge of the triangle in Fig. 2, that coincides with the configuration in which the gluon is collinear to the fragmenting quark, except in the vertex, where the gluon turns soft. Hence, the final result should be related to the thrust jet function of the fragmenting quark J q/q and to the soft thrust function S, defined in Eqs. (B.20) and (B.7), respectively. It is given by: (2.68) The 2-jet limit is obtained by expanding in powers of τ all the non-divergent terms in τ = 0 and by using Eqs. (2.25) and (2.26) to properly approximate the theta functions. This gives: Eq. (2.26) allows us to select out the terms that are divergent as z approaches 1. Furthermore, in this approximation the limit τ MAX → 0 coincides with the limit z → 1, therefore the last line of Eq. (2.69) has to be proportional to δ(1 − z). In fact, integration with a test function T (z) gives: Therefore, the vanishing cut-off limit can be written as: The final result is: where J [1] q/q is the thrust jet function for the fragmenting fermion, defined in Eq. (B.20). Notice that, as expected, the previous expression provides the missing half of the soft function when it is summed to the contribution of the sub-region R 1 , presented in Eq. (2.67). The other projection (with respect to k 1,µ k 1,ν ) is suppressed in the 2-jet limit, in fact its contribution is of order O (τ − ) and vanishes as Re < 0.
Finally, sub-region R 3 does not contribute to the 2-jet limit. Its contribution can be computed similarly to that of sub-region R 2 , but in this case thrust must satisfy the condition α = 1 − τ /z (as in Eq. (2.28)). As a consequence, the final result will be suppressed by powers of order O (τ − ), vanishing as Re < 0.
Therefore, the full contribution of the real emission to the unsubtracted final state tensor is given by: Following the same argument presented in Section 2.1.1, the cut-off τ MAX can be related to a cut-off on the transverse momentum of the fragmenting parton, as in Eq. (2.32). With this choice, we can drop the theta function anytime it appears multiplied by δ(1 − z). Therefore, we will rewrite the result of Eq. (2.73) as: where: Notice that this result is perfectly consistent with the intuitive prediction about the functions that participate to the 2-jet limit of the real emission unsubtracted final state tensor.
The whole unsubtracted final state tensor is the sum of the virtual emission (Eq. (2.50)) and the real emission term (Eq. (2.74))). This gives: where: The -expansion of V [1] is given in Eq. (2.54), while S [1] and J  is not straightforward, due to the non-trivial interplay between the theta function, that involves the λ cut-off , and the distributions in τ and z in J [1] q/q . We have already considered the product of the theta function with the plus distribution (1/τ ) + in Eq. (2.36). In this case, we also have to face the problem of multiplying the theta function by the product of two plus distributions, (1/(1 − z)) + and (1/τ ) + . This particular combination originates from the -expansion of: where the "+" label can be dropped since z never reaches 1 as long as r = λ/Q is different from zero. The integration with a test function T (z) allows us to write the previous result in terms of distributions in z: where the coefficients c n are rational numbers. The second step of the previous equation has been found by exploiting the analyticity properties of the test function T (z). In the first term in the last line of Eq. (2.80) we recognize the action of the plus distribution, (1/(1 − z)) + , while the last term in Eq. (2.80) can be neglected in the small-r limit, as we only consider the dominant (divergent) parts. Finally we have:

Subtraction Mechanism
The unsubtracted final state tensor of Eq. (2.76) describes the fragmentation of a quark at partonic level. In the final cross section, the TMD FF of the quark encodes part of this same information; hence, the overlapping content has to be properly removed from W µν f , similarly to what we did in Section 2.1.2. In this case, however, the procedure is complicated by the fact that the quark-from-quark TMD FF is itself a subtracted quantity, to avoid overlapping with the soft momentum region (see Ref. [7] and [8]). Furthermore, it needs to be renormalized by a proper UV counterterm and such renormalization affects the subtracted final state tensor, which will need an UV counterterm as well. The subtraction will be performed in b T -space, where the TMDs are explicitly defined in terms of operators (see Ref. [8] and Ref. [7]).
In momentum space, the partonic version of the 1-loop quark from quark TMD FF is given by (see Ref. [7]): Here, the term proportional to δ(1 − z) is the subtraction term of the partonic TMD FF that removes the overlapping with the soft momentum region. It involves a rapidity cut-off ζ that acts as a lower bound for the rapidity of the gluon emitted by the fragmenting quark. Its presence ensures that the gluon is actually collinear, with low transverse momentum and large rapidity. Following the argument of Section 2.1.1, the Fourier transform of Eq. (2.83) will be performed only up to k T = λ, in order to match the same collinear momentum region described by the unsubtracted final state tensor and, in particular, by J [1], (λ) q/q , defined in Eq. (2.75). The incomplete Fourier transform of 1/k 2 T was computed in Eq. (2.42). The incomplete Fourier transform of 1/k 2 T log ζ/k 2 T is given by: The incomplete Fourier transform of (2.83) follows straightforwardly. The result gives the definition of the quark-from-quark TMD FF, equipped with the cut-off λ, analogously to Eq. (2.43). We have: The label "0" reminds that the quantity in the previous expression needs to be renormalized with an UV counterterm. This is the same UV counterterm that renormalizes the quarkfrom-quark TMD FF obtained by performing a complete Fourier transform, running over the full spectrum of transverse momenta. At 1-loop it is given by (see Ref. [7]): (2.90) The result of Eq. (2.90) can easily be generalized to all orders. A 2-jet-like topology of the final state implies that the partonic final state tensor can be factorized as depicted in Fig. 5. This means: In the final result we have removed S . This fixes the renormalization scheme as the MS scheme. where the subtracted, renormalized thrust jet function for the fragmenting fermion is related to its unsubtracted analogue by: This relation can be inverted by using an inverse Mellin transform: for a suitable real c. Notice that in Eq. (2.91), all the non-trivial z dependence is encoded in the subtracted, renormalized thrust jet function J (λ) q/q , i.e. in the contribution that describes the radiation collinear to the fragmenting fermion. Finally, the partonic cross section at all order can be written as:

94)
This result closely resembles the structure of the thrust distribution associated to e + e − events, see e.g. Refs. [15][16][17][18]. Clearly, this structure will be transferred to the final cross section (see Eq. (4.2)). In particular, the (finite part 4 of the) functions V , J B and S are well-known objects that have been widely studied in the past, and they have been resummed to obtain the correct small-τ behaviour of the e + e − thrust distribution, see e.g. Refs. [19][20][21]. An analogous resummation of the new function J (λ) q/q , subtracted and renormalized, is presently not known, and its computation goes beyond the purposes of this paper. In order to try to bypass this difficulty, in the next Section we will propose a CSS-like resummation procedure for the whole partonic cross section, by deriving the proper evolution equations for each of the scales appearing in our final result. We will omit the labels "sub" and "R" since, from now on, we will always refer to subtracted and renormalized quantities.

Evolution and Resummation
Since the UV counterterm that renormalizes the subtracted partonic cross section is exactly the inverse of the UV counterterm of the TMD FF for a fragmenting fermion, the subtracted, renormalized partonic cross section obeys the following RG-evolution equation: where γ D is the anomalous dimension of the TMD FF associated to a fragmenting fermion, see Eq. (3.4). At 1-loop, it is given by: Then, the Eq. (2.95) can easily be verified at 1-loop by using the result of Eq. (2.90). The derivative with respect to the rapidity cut-off ζ plays the role of the CS-evolution equation for the TMDs, see Eq. (3.5). It is given by: Notice that, analogously to the CS-evolution for the TMDs, the kernel K does not depend on the rapidity cut-off. At 1-loop, we can compute K from Eq. (2.90): The analogy with the soft kernel K is evident if we consider the RG-evolution of K. In fact, by combining Eqs. (2.95) and (2.97) with Eq. (3.7), we have: which has to be compared with Eq. (3.8). The function γ K is the anomalous dimension of the soft kernel (see Eq. (3.8)). At 1-loop: in agreement with Eq. (2.98). Notice that, while the TMD depends only on two energy scales (µ and ζ), the subtracted, renormalized partonic cross section depends also on λ, the cut-off limiting the transverse momentum range of the fragmenting parton. Therefore, we define the following λ-evolution equation: where G is the evolution kernel. At 1-loop, its value can be computed directly from Eq. (2.90): This kernel is RG-invariant. In fact, as the r.h.s of Eq. (2.95) does not depend on λ, we can easily write: Finally, the CS-evolution for G can be found by combining Eqs. (2.97) and (2.101): The evolution equations presented in Eqs. (2.95), (2.97) and (2.101) can be solved in order to obtain an expression for the partonic cross section that resums the dependence on µ, ζ and λ. This is particularly important for the transverse momentum cut-off. In fact, all the results obtained so far are valid in the limit of λ 2 /Q 2 → 0, that corresponds to the 2-jet limit of the final state topology. On the other hand, a perturbative expansion can only be trusted when the parameter of the expansion is small. Here, even if α S is a small number, the same does not necessarily hold for the product α S log λ 2 /Q 2 , especially in the limit of vanishing cut-off. As a consequence, a resummed formula is necessary not only to obtain an all-order expression for the partonic cross section, but also to provide an adequate description of the 2-jet-like topology. We have: where we have used the RG-invariance of the kernel G. The label "ref." indicates that the energy scales are fixed to their the reference values: µ = Q, ζ = Q 2 and λ = Q. With this choice, all logarithms involving energy scale ratios vanish in the structure function F 1, f of Eq. (2.89). Therefore, the cross section computed at the reference scales can be expanded in powers of α S (Q), which now can be considered a small parameter, and the perturbative expansion is reliable to any order. The expression to NLO is given by the LO result of Eq. (A.7) added to the 1-loop result presented in Eq. (2.90) in which all logs are set to zero. This gives: It is important to notice that all terms containing the non-trivial dependence on thrust τ are multiplied by δ(1 − z), while all terms containing the non-trivial dependence on z are multiplied by δ(τ ). In Eq. (2.105), all the dependence on λ has been confined in the G-term of the exponent in the second line. Now the limit λ/Q → 0 is not an hazard anymore: in fact, as the cut-off decreases, the partonic cross section becomes more and more suppressed. This exponent can easily be computed in perturbation theory, since the strong coupling α S is evaluated at the scale Q and, consequently, it is not involved in the integration. We can use the 1-loop expressions for K [1] and G [1] in Eqs. (2.98) and (2.102) to obtain: Therefore, for any ζ = 0, the partonic cross section goes to zero when the λ cut-off vanishes. As observed in the final part of Section 2.1.2, this cut-off is related to the measured value of thrust. In particular, since τ meas. = λ 2 /Q 2 the suppression is associated to the topology of the final state. This correctly implies that for perfectly pencil-like event (τ meas. = 0), the partonic cross section vanishes. Clearly, this is a very unrealistic case, in which the radius of the jet is zero and the fragmentation process develops along the direction of the fragmenting quark. The probability of finding a τ meas. = 0 event is then exactly zero.
Eq. (2.105) represents the exact expression of the partonic cross section, to all orders. In particular, the exponent in the second line of Eq. (2.105) correctly expresses its allorder resummation with respect to λ (and ζ). Instead, in Eq. (2.107) the exponent is truncated to 1-loop accuracy in the perturbative expansion. This has some effects on the log counting. In fact, a rigorous NLL computation would also include terms of order α 2 S that we explicitly neglected in the NLO computation of the partonic cross section. In particular, the correlated emission of two gluons, each going in one of the two hemispheres defined by the thrust axis, can produce non-global logarithms (NGLs) of λ 2 /Q 2 . A proper treatment of such contributions goes beyond the purpose of this paper, since our partonic cross section is computed up to NLO, i.e. considering only single gluon emissions. The non-global effects of e + e − → H X have been investigated in Refs. [12] and [13]. General, non-global observables have been studied e.g. in Refs. [22][23][24][25]. The inclusion of the contributions beyond the NLO in the resummed partonic cross section generates difficulties associated to the determination of the resummed expression for the subtracted, renormalized function J (λ) q/q (see the discussion in the end of the previous Section). The resummation proposed above allows to bypass such difficulties by providing a result that can be trusted at fixed order in perturbation theory.

TMD FF for the fragmenting quark
The long-distance behavior of the cross section for e + e − → H X presented in Eq. (1.4) is encoded in the unpolarized TMD FF, which provides a 3D description of the fragmentation of a parton of type j into a spinless hadron H. TMD FFs depend on the transverse momentum of the parton with respect to the direction of the produced hadron and reflect the dynamics of hadronizing partons. In the computation of the partonic cross section carried out in Section 2, we have found that only the TMDs associated to fragmenting fermions survive in the sum over parton-types in Eq. (1.4), as the contribution of the fragmenting gluon is suppressed in a 2-jet-like topology, see Eq. (2.46). Therefore, in the following we will define the TMD FF of a fragmenting fermion of flavor f , providing a review of all the explicit expressions relevant to its computation to NLL accuracy.
The TMD FFs are defined in the Fourier conjugate space of the transverse momentum k T of the fragmenting parton. The conjugate variable is commonly denoted as b T . Furthermore, they are equipped with a rapidity cut-off y 1 which constrains the rapidity of the parton described by the TMDs. We will use the following definition for D 1 (see Ref. [7] and Ref. [8]): where y P is the rapidity of the detected hadron H. This definition holds for any fragmenting parton of type j. The function Z j is the UV counterterm that renormalizes the TMD FF.
Its dependence on j originates from the color representation of the fragmenting parton (gluon or fermion). Eq. (3.1) shows how the overlapping between the TMD and the soft momentum region is removed by a subtraction mechanism similar to that used in Section 2 to cancel the double counting between the hard and the collinear momentum region. In the case of a fragmenting fermion of flavor f , at numerator we have the unsubtracted TMD FF defined as: where x = (0, x − , b T ) and P is the momentum of the outgoing hadron H. The (renormalized) quark field is ψ f , while the operators W q are Wilson lines associated to a fermion field along the light-like minus direction w − . The label "NO S.I." indicates that the Wilson line self energies must not be considered, while the label "(0)" refers to the fact that the definition of Eq. (3.2) holds for the bare unsubtracted TMD, defined with bare fields. For this reason, we multiplied the whole expression by Z 2 , which is the wave-function renormalization factor for the quark field. This fixes the UV counterterm Z q as the renormalization factor of the bare TMD, defined with bare fields. Notice that since the Wilson lines in D 1, H/f are defined along a light-like direction, the unsubtracted TMD is affected by unregulated rapidity divergences. The denominator in Eq. 3.1 is a 2-h soft factor 5 and it represents the overlapping contribution with the soft momentum region. Its subtraction removes the double counting and regulates the rapidity divergences. In fact, S 2-h represents all the soft gluons emitted by the fragmenting fermion along the jet direction, i.e. all the gluons with a small transverse momentum and a rapidity no larger than the rapidity cut-off y 1 . The 2-h soft factor is a color singlet. Then S 2-h is defined as the coefficient of the identity matrix in color space: where the trace is over all color indices. Also in this case, the label "(0)" reminds that Eq. (3.3) defines the bare soft factor, in terms of bare fields. Finally, let us point out that the TMD FF D 1 as defined above represents all the collinear radiation of the jet initiated by a fermion of flavor f that fragments into the hadron H. By collinear, we mean that all partons described by the TMD FF have a small transverse momentum k T and a large rapidity y that lies in the range y 1 ≤ y ≤ y P . The TMD FF D 1 defined in Eq. (3.1) obeys the following evolution equations: 5 See the classification in Ref. [8].
CS-evolution, (3.5) where, as usual, we defined: and k + is the plus component of the momentum of the fragmenting parton. The function γ D is the anomalous dimension of the TMD FF, while K is the rapidity-independent kernel of the CS-evolution. They, in turn, solve the following equations: where γ K is the anomalous dimension of the soft kernel K. The solution to Eqs. (3.4) and (3.5) is given by [26,27]: In the previous expression, the reference scales are µ = µ b and ζ = µ 2 b , where µ b is defined as: where the b T prescription allows to separate the perturbative small-b T behavior of the TMD from its non-perturbative large b T content.
We also introduce a minimum value b min , that allows to recover the collinear FFs by integrating over the transverse momentum of the fragmenting parton. Therefore we will adopt the modified b T prescription defined as: (3.12) In this paper, we use standard choices for the minimum and maximum value of b T involved in the previous equations. In particular, we set b max = 1 GeV −1 and b min = C 1 /Q, where It is important to underline that the TMD definition adopted here, Eq. (3.1), differs from the usual definition commonly used in the factorized cross sections of 2-h class processes like SIDIS and e + e − → H 1 H 2 X, with the two hadrons almost back to back (see Refs. [8] and, e.g., Ref. [7]). In such processes, the final cross section presents a soft factor connecting the target and the detected hadrons, in SIDIS, and the two back-to-back hadrons, in e + e − → H 1 H 2 X. This soft factor cannot be fully computed in perturbative QCD, since it encodes non-perturbative information about the soft radiation that flows through the two collinear parts. But it cannot be directly extracted from experimental data either, since it always appears in connection to the two collinear parts. Notice that, despite the similarities, this soft factor is not the same object that we encountered in the derivation of the partonic cross section (see e.g. Eq. (2.94) or Fig. 5). In fact, as shown to 1-loop accuracy, S( ; τ ) is totally predicted by perturbative QCD, see Eq. (B.7). To overcome the difficulties induced by the presence of a non-trivial soft term in the final cross section, the prescription of Ref. [7] would include part of this soft factor in the definition of the TMDs. This has the advantage of freeing the cross section of the explicit presence of the soft factor, but creates a hadron-class dependence in the TMDs, lowering their degree of universality. In Ref. [8] we investigated the relationship between this currently accepted definition, referred to as "square root definition", and the "factorization definition", adopted in this paper and free of any external soft contribution. The two definitions are totally equivalent at small-b T , in the perturbative region, but they show a rather different large-b T , non-perturbative, behavior. It is not by chance that the "square root definition" is obtained by multiplying the TMD defined in Eq. (3.1) by the square root of M S (b T ), the non-perturbative function that models the large-b T behavior of the soft factor: The "square root definition" is optimal for cross sections corresponding to the 2-h class (Drell-Yan, SIDIS, e + e − → H 1 H 2 X) in which two collinear parts appear, each associated to one of the two reference hadrons. However, this definition lowers the degree of universality of TMDs, since it includes some extra soft physics information typical of the 2-h class. In contrast, the "factorization definition" of Eq. (3.1) allows to define a totally universal TMD, which can be applied equally well to processes like e + e − → H X that do not belong to the 2-h class. Notice that Eq. (3.13) is of crucial importance from a phenomenological point of view, as it relates the TMDs obtained from data analyses based on the square root definition (widely used in the last decade) to the TMDs extracted using the factorization definition. This implies that all previous work on the extraction of polarized and unpolarized TMDs could still be easily exploited in global analyses. The universality breaking effects generated in processes belonging to different hadron classes, as addressed in Ref. [8], have also been investigated in Ref. [28]. Here TMD factorization for 3-h class processes, like dijet and heavy-meson pair production in SIDIS, is presented in a SCET framework, where the introduction of a new soft function is required.
In the following, we will focus on the main three ingredients of Eq. (3.9) separately. In particular, in Section 3.1 we will review the NLO expression for the TMD FF at the reference scale, in Section 3.2 we will compute the perturbative Sudakov factor at NLL precision and finally, in Section 3.3 we will consider the non perturbative content of the TMD FF.

Operator Product Expansion at NLO
The first term in Eq. (3.9) is the TMD FF at the reference scales µ = µ b and ζ = µ 2 b , where µ b has been defined in Eq. (3.10). It is related to the short-distance, small-b T -behavior of D 1 and therefore computable in perturbation theory. The application of the factorization procedure to the TMD itself shows that D 1 can be expressed in this region as an Operator Product Expansion (OPE), in which the basis of the operators is given by the collinear FFs d h/k . The Wilson coefficients C k/f of the OPE are fully computable in perturbative QCD and can be expanded in powers of α S (µ b ). Having set the scales to the reference values, no dangerous logarithms will affect the perturbative expansion, that can be considered reliable at any order. In particular, at LO the Wilson coefficients are just delta functions: (3.14) To 1-loop we have the following expressions, see e.g. Refs. [7,29,30]: (3.15) Notice that the plus prescription in Eq. (3.15) can be dropped. Then, to NLO, the OPE is given by: q/q (ρ) + where we used Eqs. (3.14), (3.15) and (3.16).

Perturbative Sudakov Factor at NLL
The second term in Eq. (3.9) is the perturbative part (small-b T ) of the Sudakov factor, which originates from the resummation of the TMD with respect to the scales µ and ζ.
The exponent cannot be expressed as a fixed order expansion, since the logarithms of the energy scales multiply α S at any order. Therefore, the proper way to expand this quantity consists in counting the power of such logs and performing a NLL expansion (this guarantees consistency with the NLO precision of the OPE, see Section 3.1). This operation is easily done by separating out the part that depends on the rapidity cut-off ζ from the rest of the exponent appearing in the perturbative Sudakov factor. This gives: The recipe to obtain the previous quantity at NLL accuracy is the following: • The anomalous dimension γ K of the soft kernel is expanded up to 2-loops. The 1-loop coefficient can be found in Eq. (2.100), while at 2-loops we have, see Refs. [29,30]: where n f is the total number of fermion fields that we consider.
• All the other quantities in the exponent are expanded up to 1-loop. Their expressions are given by, see Refs. [29,30]: where Eq. (3.20) refers to the 1-loop coefficient of the soft kernel computed at the reference scale.
Therefore, the ζ-independent part of the Sudakov at NLL is: where: 23) and the functions g 1 and g 2 are given by: where β 0 and β 1 are the coefficients of the beta functions up to 2 loop: This result is in agreement with the Sudakov factor computed e.g. in Ref. [31]. On the other hand, the term depending on the rapidity cut-off ζ is given by: where x has been defined in Eq. (3.23), while the functions g K 2 and g K 3 are: K + 1 x K [1] . (3.30)

Non-Perturbative content
The last term in Eq. (3.9) encodes the whole non-perturbative content of the TMD fragmentation function, D 1 . Clearly, this cannot be predicted by perturbative QCD and hence it has to be extracted from experimental data, through a phenomenological analysis. It involves two functions. The first is g K , which describes the long-distance behavior of the soft kernel K, defined as: In most phenomenological applications it is assumed to behave quadratically: with a ∼ 0.01 ÷ 0.1 GeV 2 .
The other non-perturbative function is the model M D , which represents the very heart of the TMD, since it identifies D 1 uniquely. Besides b T , it may depend on the flavor of the fragmenting parton, on its light-cone momentum fraction z and also on the detected hadron H (in the following, however, we will assume that it is just a function of b T , dropping any other dependence). Its functional form is arbitrary, as long as it goes to 1 as b T → 0 and goes to zero in the large-b T -limit, b T → ∞. It is commonly modelled as a Gaussian or as a (second kind) Bessel function, which in the Fourier conjugate momentum space becomes a power-law, reminding of a squared propagator for an on-shell fermion. This will be discussed in more detail in Section 5.

Final Cross Section
We are now ready to assemble the expression for the final cross section of e + e − → H X in the 2-jet limit, as sketched in Eq. (1.4) in its general factorized form. It is differential in the fractional energy of the detected hadron z h , defined in Eq. (1.1), in the transverse momentum P T of the detected hadron with respect to the thrust axis, and also in thrust, T , as defined in Eq. (1.3), The factorization procedure allows to separate out the shortdistance contributions, encoded in the partonic cross section (Section 2), from the longdistance behavior, described by the TMD FFs (Section 3). Referring to the very general expression presented in Eq. (1.4), we should keep in mind that: • The contribution of the fragmenting gluon is suppressed by orders of O λ 2 /Q 2 in the 2-jet limit (see Section 2.1.2). Therefore, the sum over parton types, j, runs only over the flavors f of the fragmenting fermions.
• The final cross section is RG-invariant. In fact, the anomalous dimension of the partonic cross section, Eq. (2.95), is exactly equal and opposite to the anomalous dimension of the TMD FF associated to a fragmenting fermion, Eq. (3.4). Therefore, we are allowed to set µ = Q.
• The final cross section depends on two cut-offs, λ and ζ, so far considered as being independent. However, they have been introduced for intimately related reasons: they both ensure that the partonic cross section describes all particles involved in the process, except those that are collinear i.e. with low transverse momentum (k T ≤ λ) and a very large rapidity (y ≥ y 1 ), which are included in the TMD FF. In fact, a neat and consistent separation of the partonic cross section from the TMD fragmentation function guarantees the successful factorization of the cross section. Furthermore, both cut-offs have to be intended in the limit in which they vanish. On one hand, ζ → 0 (i.e. y 1 → +∞, according to Eq. (3.6)) is the limit in which the factorization procedure of Ref. [7] can be applied. On the other hand, λ → 0 corresponds to the 2-jet limit for the topology of the final state of the process, according to the relation λ 2 /Q 2 = τ meas. , as we discussed at the end of Section 2.1.2. The above considerations suggest that λ and ζ are closely related. The subtraction mechanism presented in Section 2.2.2 was built in such a way that the range of transverse momentum of the collinear emission covered by the unsubtracted hard part of Eq. (2.76) would exactly match the k T -range covered by the subtraction term of Eq. (2.87). In fact, in both cases the collinear gluon can have a transverse momentum k T no larger than λ. As far as rapidity is concerned, however, the subtraction mechanism needs to be applied with special care. Let's consider the 1-loop computation of the partonic cross section associated to a fragmenting fermion. On one hand, the subtraction term describes a collinear gluon that, when τ = 0, has a rapidity larger than y 1 , linked to ζ by the definition of Eq. (3.6). On the other hand, the rapidity of the collinear gluon described by the unsubtracted hard part depends on k T and τ , and when τ = 0 it cannot be smaller than y λ , defined exactly as y 1 but with ζ = λ 2 . Therefore, in order to cover the same range of rapidities in both terms, we have to set y 1 = y λ , which means ζ = λ 2 . Any different value of ζ would lead to an incorrect separation of the collinear and the hard contributions in the final cross section. In fact, if ζ > λ 2 , then y 1 < y λ and the hard part would also include the contribution of the collinear particles with a rapidity in the range y 1 ≤ y ≤ y λ , which is already covered by the TMD FFs. Conversely, if ζ < λ 2 , then y 1 > y λ and the contribution of the collinear particles with rapidity in the range y λ ≤ y ≤ y 1 would not be present at all in the final result. Consequently, the only possible choice is ζ = λ 2 .
Such considerations lead us to the following formula: dσ dz h dT meas. dP 2 not exactly coincide with the variable T = 1 − τ that we used in Section 2 in computing the partonic cross section, since the thrust distribution calculated at partonic level does not encode the information about the non-perturbative hadronization process. According to Ref. [15], the observed distribution of "hadronic" thrust, T meas. , is related to its partonic counterpart, T , by a correlation function C(T meas. , T ). Therefore, the partonic cross section written in terms of T meas. will be obtained as: having used the expression for the partonic cross section as obtained in Eq. (2.105). Monte Carlo simulations show that the correlation function is sharply peaked around T meas. ∼ T and that it is not strongly sensitive to the parton distribution [15]. Therefore, a formula suitable for phenomenology (that clearly requires an expression written in terms of functions of T meas. ) can be obtained by approximating the correlation function with a (smooth) delta function and then neglecting all the contributions proportional to C(T meas. , 1), associated to strictly pencil-like events. In other words, the partonic cross section written in terms of τ meas. can be obtained from the result of the fixed order computation, expressed in terms of distributions of τ , by simply neglecting all contributions proportional to δ(τ ), by removing all the "plus" that label the plus-distributions and finally by replacing τ with its hadronic counterpart. At NLO this operation gives: Of course, this is a rather rough approximation which reduces the potential of a formula like Eq. (2.106). In particular, neglecting the τ = 0 terms not only eliminates the contributions of pencil-like events, but it also drastically modifies the behaviour of the cross section in the region where T is very close to 1. Furthermore, by applying this approximation we loose part of the information on the z h -dependence of the cross section, since all z-distributions that are coupled to δ(τ ) are removed. As a consequence, we expect that the dependence on z h will be compromised, especially at large values of thrust. The z h -dependent neglected part, associated with the last line of Eq. (2.106), involves two terms: an integral from z h to 1 and a contribution dominated by log (1 − z h ) 2 . Therefore, the discrepancy due to the brute approximation of Eq. (4.4) should be more important at small and large values of z h , where these two terms have a non negligible impact on the final cross section. We will come back to this point in the next Section.
Finally, we require λ 2 /Q 2 = τ meas. , and we remove the label "meas", which has now become redundant. At NLO, NLL accuracy, the cross section for e + e − → H, X differential in z h , T and P T is then given by: where, according to Section 3: with x defined as in Eq. (3.23). For a convenient and straightforward application of Eq. (4.6) we recall that the Wilson coefficients C [1] q/q and C [1] g/q are presented in Eqs. (3.15) and (3.16), the functions g 1 , g 2 and g K 2 , g K 3 contributing to the Sudakov are computed in Eqs. (3.24), (3.25), (3.29) and (3.30), and the non-perturbative functions g K and M D are defined in Eqs. (3.32) and (5.1). Notice that the cross section in Eq. (4.5) is not resummed in thrust. A proper resummation in T is beyond the purpose of this paper. Such resummation must also include a correct treatment of the dependence on z h , by considering the terms that have been neglected in the approximation used in Eq. (4.4). Clearly, this is strictly connected to the difficulties in finding a fully resummed expression for the subtracted, renormalized function J (λ) q/q (see the discussion at the end of Section 2.2.2) and, ultimately, for the whole second line of Eq. (4.2).
Very recently, the factorization of the e + e − → HX cross section, as measured by the BELLE Collaboration (Ref. [1]), has been investigated in two papers, both based on the SCET formalism. In Ref. [12], the authors propose to integrate out the thrust T and to reproduce the experimental cross section by combining all the measured thrust bins, within the range [0.5−1.0]. A cross section differential in z h , P T and T , is presented in Ref. [13]; it results from matching three different kinematical regions, each associated with a different factorized expression for the final cross section. The phenomenological application of this formula is not shown.

A basic phenomenological application
In this conclusive Section we provide a basic but realistic example of the strong phenomenological potential of Eq. (4.5). We stress that this is not a proper and complete data-analysis of the BELLE measurements, as in what follows there is no fit of the experimental data and no fine tuning of the free parameters. However, the agreement of the predictions obtained from Eq. (4.5) with the experimental data is astonishing, even without a proper fit.
As anticipated in Section 3.3, to model the non-perturbative part of the TMD FF we will exploit a parameterization that, in the b T space, corresponds to a Bessel-K function, normalized in such a way that it is 1 at b T = 0: where K −1+p is the modified Bessel function of the second kind. This model was successfully applied in Ref. [32] and, more recently, in Ref. [8]. It is known as the "power law model", since its Fourier transform is given by: which mimics the power p of a function closely reminiscent of the propagator of an on-shell fermion (see dedicated discussion in Ref. [33]). As mentioned above, we do not perform a proper fit but simply fix the two non perturbative parameters of the k T -distribution to reasonable and perfectly standard values, inspired by pure common sense: m = 1 GeV, to represent a typical hadronic mass, and p = 2 which translates, in k T -space, in a k T -distribution similar to the square of the Feynman propagator, see Eq. (5.2). One additional free parameter, a, is embedded in the non perturbative g K (b T ) function, as defined in Eq. (3.31). Its value depends on the choice of b M AX ; for b M AX = 1 GeV −1 it is expected to be in the range 0.01 − 0.1 GeV, while considerably larger values are found at smaller b M AX . Interesting discussions and reference values of a can be found, for example, in Refs [34][35][36][37]. We set a = 0.05 GeV −1 , as an intermediate value. There is clearly some degree of correlation among the three parameters a, m and p, so that different values of a can easily be compensated by slight variations of m and p. Reassuringly, the final results are not particularly sensitive to the precise value of a, after m and p have been fixed. It is important to point out, here, that we do not need to introduce any overall normalization parameter.
The collinear fragmentation functions used for the implementation of Eq. (4.5) are the NNFF10 NLO reference set [38], for π + production.
To study the z h and P T behaviour of the cross section at fixed values of thrust, T , we consider the data subset that the BELLE Collaboration has chosen to present in Fig.  9 of their article, [1]. The differential cross sections corresponding to pion production are shown as a function of P T for 15 z h -bins, covering the range [0.10 − 0.85] and one thrust bin, 0.85 < T < 0.9 (we fix T = 0.875). The error bars correspond to statistical and systematical errors added in quadrature, as reported by the BELLE Collaboration. For Figure 6. The e + e − → hX cross section computed to NLO and NLL accuracy, corresponding to 15 z h -bins as measured by the BELLE Collaboration. The central values of the z h -bins is indicated in each panel. Thrust is fixed to one bin, corresponding to 0.85 < T < 0.9. Uncertainty bands are not shown as no fit is involved in this computation. The collinear fragmentation functions used for the implementation of Eq. (4.5) are the NNFF10 NLO reference set [38], for π + production. Uncertainties due to the errors on the collinear fragmentation functions are not shown. the BELLE experiment, Q = 10.58 GeV, while the P T bins span a range from 0.06 GeV to 2.5 GeV. We expect our TMD description to hold at low P T 's (more properly where the ratio P T /z h Q is small). In our case P T should certainly be no larger than 1.2 − 1.5 Gev/c, depending on the corresponding value of z h . Fig. 6 shows the comparison between the BELLE experimental measurements and our predictions, represented by the solid (red) line, obtained by applying Eq. (4.5), giving the e + e − → HX cross section, to 1-loop and NLL accuracy, in the 2 jet limit. Our results are presented without uncertainty bands as there is no fitting involved in this calculation. The uncertainty induced by the error on the collinear FFs is not shown in this plot. The agreement of our results with experimental data is extraordinary, and gives strong support to the naturalness with which Eq. (4.5) works. In particular, our plots show how the full range of z h values explored by BELLE measurements can be described to an unprecedented level of quality.
To analyse the T behaviour of the e + e − → HX cross sections we will refer to a different subset of the BELLE experimental data, large enough to span the whole range of values of T , z h and P T where Eq. (4.5) can safely be applied. In particular we will consider three bins in T , corresponding to 0.70 < T < 0.80, 0.80 < T < 0.85 and 0.85 < T < 0.90. In each of these T bins, we will consider three bins in z h , corresponding to 0.35 < z h < 0.40, 0.55 < z h < 0.6 and 0.75 < z h < 0.80.
The results obtained by using Eq. (4.5) are shown in Fig. 7, where they are compared to the BELLE experimental data corresponding to the nine T and z h bins described above (individual plot legends clearly indicate their values in each panel). As in Fig. 6, error bars corresponds to statistical and systematical errors added in quadrature. The solid lines correspond to our prediction for the e + e − → HX cross section to NLO and NLL accuracy, in the 2 jet limit: they are presented without uncertainty bands as there is no fitting involved in this calculation. The uncertainty induced by the error on the collinear unpolarized FF is not shown. As for the previous plots, the agreement with the experimental data is extremely good, especially if one considers that it has been obtained without a proper fit or a fine tuning of the parameters. This confirms the remarkable success of Eq. (4.5), which can reproduce the experimentally measured e + e − → H X cross section in a perfectly natural way, and with a very mild dependence on the specific values of the free parameters.
Outside of the considered range of z h and T values, the quality of the description starts deteriorating slowly. This is more evident in Fig. 8 where our predictions of the e + e − → H X cross section are plotted as functions of thrust, T , in three z h bins and three representative P T bins (one very low, one central and one rather large value of P T , within the TMD range). While the last T -bin is always greatly overestimated, the description of data deteriorates at extreme values of z h . This is somehow expected. In fact, as we discussed in Section 4, a combined resummation of Eq. (4.5) in T and z h would allow us to take into account terms which are presently neglected and restore a good agreement down to the lowest values of z h . Such a resummation is a project of its own, which we might consider undertaking in the future.
The high quality of the description of the BELLE experimental data as a function of thrust represents our most important achievement. For the first time after their publication, these valuable measurements can be fully explained, in their z h , P T and T behaviour, within a sound and well founded scheme of factorization, with minimal bias from non-perturbative free parameters.

Conclusions
In this paper we have constructed the factorized e + e − → HX cross section in the 2-jet limit, providing its explicit expression to NLO and NLL accuracy. Working in a CSS-like framework, factorization has been built, step by step, by means of power counting rules and applying an appropriate subtraction mechanism which allowed us to neatly separate the short distance, low-b T behaviour, from the long-distance behaviour. In this way, the cross section results in the convolution of the partonic cross section, which in our formalism includes all of the hard contributions, with a TMD fragmentation function, defined according to the scheme presented in Ref. [8]. Here TMDs are purely collinear parts, universal and totally process independent.
The 2-jet limit strongly constrains the final state topology to an extreme configuration, in which the available phase space is strongly reduced, as shown in Figs. 2 and 3. In our Figure 7. e + e − → hX cross section computed to NLO and NLL accuracy, corresponding to 9 of the bins measured by the BELLE Collaboration. The value of z h increases by moving down in each column of panels: the upper row of panels corresponds to z h = 0.375, the middle row to z h = 0.575, the lowest row to z h = 0.775. Moving left to right, instead, corresponds to increasing values of thrust. In the leftmost column of panels T = 0.750, in the central column T = 0.825, in the rightmost column T = 0.875. Uncertainty bands are not shown as no fit is involved in this computation. The collinear fragmentation functions used for the implementation of Eq. (4.5) are the NNFF10 NLO reference set [38] for π + production.Uncertainties due to the errors on the collinear fragmentation functions are not shown. treatment, this is realized by introducing in the computation of the partonic cross section a cut-off τ MAX , ultimately related to the measured value of thrust. As a consequence, the TMD FFs appearing in the final cross section acquire a dependence on the thrust, through their rapidity cut-offs.
In Section 5 we have performed a basic phenomenological analysis in which we compare the BELLE experimental data to the results obtained by using of our final formula, Eq. (4.5), which is differential in z h and P T as well as in thrust, T . This is the first phenomenological analysis of the thrust dependence of the e + e − → HX cross section. The agreement of our predictions with experimental data is excellent, even without performing a proper fit. In fact, in our calculation we fix the value of the three free parameters which regulate the non-perturbative behaviour of the TMD FF (one for g k and two for the model M D ) by completely general considerations, based on pure common sense.
The results obtained are very promising, and will represent the basis on which we will construct a more refined, global analysis, simultaneously considering the unpolarized cross sections from SIDIS and e + e − → HX processes. This will be possible thanks to Figure 8. e + e − → hX cross section computed to NLO and NLL accuracy plotted as a function of thrust, T , in three z-bins and three (representative) P T -bins, as indicated in the plot legends.
In the upper row of panels z = 0.775, in the central row z = 0.575 and in the lower row of panels z = 0.375. Uncertainty bands are not shown as no fit is involved in this computation. The collinear fragmentation functions used for the implementation of Eq. (4.5) are the NNFF10 NLO reference set for π + production [38].Uncertainties due to the errors on the collinear fragmentation functions are not shown. the formalism proposed in Ref. [8], that allows to define any TMD parton density in a completely universal way, separating out the soft and process dependent part, and finally restoring the possibility to perform global analyses including data from processes belonging to different hadron-classes. of the frame used to compute the partonic cross section, i.e. the partonic analogue of the hadron frame, where the outgoing parton has zero transverse momentum. All quantities will be computed in dimensional regularization, with D = 4 − 2 being the spacetime dimension.

B.1 Soft Thrust Function
The 1-loop soft thrust function S( , τ ) is obtained as a result of the application of the soft approximation (T S in Ref. [7]) to the Feynman graphs corresponding the 1-loop partonic cross section. The case of the fragmenting gluon is suppressed in this approximation, therefore it can be neglected. On the other hand, when the fragmenting parton is a fermion, the emitted gluon can be either virtual or real. The virtual contribution vanishes in full dimensional regularization, hence the soft thrust function only originates from the real emission graphs. In order to compute the r.h.s. of Eq. (B.1), we will label the momenta of the three final state particles as in Section 2.2.1, where k 1 refers to the fragmenting quark, k 2 to the antiquark and k 3 to the gluon. In this case we have to use their soft approximations: where q + = q − = Q/ √ 2. Notice that the soft approximation leaves k 3 unchanged in its functional form, but the size of all its components are now of order λ 2 /Q, where λ is a very low energy scale according to power counting. Depending on the hemisphere in which the soft gluon has been emitted, we have two different, equally probable, configurations: τ = y Therefore, by defining the usual light-like directions w 1 = (1, 0, 0 T ) and w 2 = (0, 1, 0 T ), the 1-loop soft thrust function is obtained as: Notice that all the divergences of S( , τ ) are regulated by dimensional regularization, and we did not encounter any unregulated rapidity divergence. This is due to the explicit presence of τ , that acts as a regulator. In fact, the expression in Eq. (B.7) is divergent if either or τ vanish. The -expansion of Eq. (B.7) follows from its integrability property.
In fact, S [1] ( , τ ) is integrable with respect to τ , but its trivial expansion in powers of does not. This problem is solved by considering an expansion in terms of τ distributions instead of functions. Then, one can easily prove that:

B.2 Backward Thrust Function
The 1-loop backward thrust function J B ( , τ ) is obtained by applying the collinear approximation to the backward emission (analogue to T B in Ref. [7]) in the 1-loop Feynman graphs. Here "backward" refers to the direction of emission, which goes along n (negative z-axis), opposite to the thrust axis. As for the previous case, the contribution of the fragmenting gluon is suppressed by power counting and virtual graphs are zero in full dimensional regularization. In the following, we will assume that the fragmenting parton is a quark, hence the approximation concerns the configuration where the gluon is collinear to the antiquark. The case of a fragmenting antiquark is perfectly analogous. The backward approximation applied to the momenta (labeled as usual) gives: As in the previous case, the backward approximation leaves k 3 unchanged in its functional form, but resizes its components as l + ∼ λ 2 /Q, l − ∼ Q and l T ∼ λ, where λ is a very low energy scale according to power counting. The emission of the collinear gluon along the direction of the fragmenting quark does not give any contribution, therefore we only have to consider the emission in the antiquark hemisphere: as the total momentum entering into the antiquark hemisphere, the 1-loop backward thrust function is: