Universality-breaking effects in $e^+e^-$ hadronic production processes

Abstract

Recent BELLE measurements provide the cross section for single hadron production in $e^+e^-$ annihilations, differential in thrust and in the hadron transverse momentum with respect to the thrust axis. Universality breaking effects due to process-dependent soft factors make it very difficult to relate this cross sections to those corresponding to hadron-pair production in $e^+e^-$ annihilations, where transverse momentum dependent (TMD) factorization can be applied. The correspondence between these two cross sections is examined in the framework of the Collins-Soper-Sterman factorization, in the collinear as well as in the TMD approach. We propose a scheme that allows to relate the TMD parton densities defined in 1-hadron and in 2-hadron processes, neatly separating, within the soft and collinear parts, the non-perturbative terms from the contributions that can be calculated perturbatively. The regularization of rapidity divergences introduces cut-offs, the arbitrariness of which will be properly reabsorbed by means of a mechanism closely reminiscent of a gauge transformation. In this way, we restore the possibility to perform global phenomenological studies of TMD physics, simultaneously analyzing data belonging to different hadron classes.

1 Introduction

QCD describes hadronic matter through the dynamics of its elementary consituents, quarks and gluons. However, confinement prevents the direct observation of partonic degrees of freedom, which are shaded by the hadronization mechanism.

Recently, the BELLE Collaboration at KEK has measured the $e^+e^- \rightarrow H X$ cross section at a c.m. energy of $Q^2\sim 112$ $\hbox {GeV}^2$ as a function of $P_T$, the transverse momentum of the observed hadron h relative to the thrust axis [1]. The data are binned in $P_T$ and selected in thrust in such a way to ensure that $T \sim 1$, which corresponds to a two jet configuration. This is one of the measurements which go closer to being a direct observation of a partonic variable, the transverse momentum of the hadron with respect to its parent fragmenting parton.

These data have indeed triggered a great interest of the high energy physics community, especially among the experts in the phenomenological study of TMD phenomena and factorization. However, there are difficulties in the analysis of these data as, due to the nature of this process, a TMD factorization as that formulated in Ref. [2] cannot be directly applied. In this case, in fact, collinear factorization would rather be the correct approach.

In this paper we will follow very closely the formulation proposed by J. Collins in Ref. [2] for $e^+e^- \rightarrow H_A \, H_B \, X$ processes and we will give a different definition of TMDs which, by extending their degree of universality, becomes suitable to be applied also to the $e^+e^- \rightarrow H \, X$ process. We will move along the lines suggested, for instance, in Ref. [3].

In this new definition, the soft factor of the process, which is responsible for potential universality breaking effects, is not absorbed in the TMD, to prevent it from influencing its genuinely universal nature. Instead, it appears explicitly in the cross section which acquires a new term, that we will call soft model, $M_S$. After being modelled using a suitable parameterization, it can be extracted from experimental data. While the TMDs are truly universal and can be extracted from any process, $M_S$ is universal only among a restricted number of processes. In other words, $M_S$ is universal only within his hadron class. Later on in the paper we will define what we mean exactly by “class”; for the moment being we anticipate that, for instance, Drell-Yan, Semi Inclusive Deep Inelastic Scattering (SIDIS) and $e^+e^- \rightarrow H_A \, H_B\, X$ processes belong to the same hadron class, while DIS and $e^+e^- \rightarrow H\, X$ belong to a different class.

The advantage of this formulation is that a well defined expression relates TMDs extracted using different definitions. Consequently, all results obtained in past phenomenological analyses can easily be reformulated according to this new framework, with no loss of information.

We stress that the factorization procedure itself will not be altered from its original form. Rather, we introduce a new methodology to implement the phenomenological application of that same scheme, changing the focus on the fundamental ingredients of the phenomenological models.

The paper is organized as follows. In Sect. 1.1 we will outline the basics of TMD and collinear factorization. Section 2 will be devoted to the study of the soft factor and its factorization properties, while in Sect. 3 we will examine the collinear parts of hadronic processes. Here we will define the TMDs and show how a particular transformation, which we will call “rapidity dilation”, allows to consider them invariant with respect to the choice of the rapidity cut-offs introduced by the regularization of the rapidity divergences. In Sect. 4 we will briefly outline a new way of classifying hadronic processes in terms of their “hadron class”. Section 5 will be dedicated to the study of 2-hadron processes and their cross sections, respectively. Finally, in Sect. 6 we will apply this formalism to $e^+e^- \rightarrow H\, X$, giving a simple example of how this scheme can be applied to a phenomenological analysis. Appendices A, B and C will be dedicated to the definition of Wilson lines, to the small $b_T$ behaviour of the soft factor and of the TMDs, and to the kinematics of a $e^+e^- \rightarrow H\, X$ process, respectively.

Throughout the paper we will adopt a pedagogical approach, as we intend to provide a review which could be useful to young beginners as well as to experienced researchers.

1.1 Collinear and TMD factorization

Modern studies of high energy QCD processes are based on factorization, a procedure that allows to separate the cross section of a hadronic process involving a hard energy scale Q into a part which is fully computable in perturbation theory and a non-perturbative contribution, with an error suppressed by powers of m/Q, where m is a typical low energy mass scale. In general, the perturbative part is process dependent but it can be computed, at any given order, for any given process. The non-perturbative part, cannot be computed: it should rather be inferred from experimental data. However, when defined in an appropriate way, it is universal, in the sense that it can be extracted from one process and then used in any other. If factorization applies and universality is preserved, then the theory can be predictive. Nowadays, several different schemes are available to implement factorization. In the following, we will adopt the modern version of the Collins-Soper-Sterman (CSS) scheme [4, 5], often referred to as CSS2, presented in Ref. [2].

When factorization applies, then the cross section of the process will appear as a convolution of contributions which can be classified in terms of the following three categories:

• 1. Hard part It corresponds to the elementary subprocess and it provides the signature of the process, as it identifies the partonic scattering uniquely. It is fully computable in perturbation theory in terms of Feynman diagrams, up to the desired accuracy.

• 2. Collinear parts These contributions are associated to the initial and final state hadrons of the process and contain the collinear divergences related to the massless particles emitted along the hadron direction. Each of them corresponds to a bunch of particles strongly boosted along this direction, which move almost collinearly, very fast. Due to their characteristic divergences, collinear parts cannot be fully computed in perturbation theory: their non-perturbative content has to be extracted from experimental data. Among all the particles in the collinear group, two of them deserve special attention: the reference hadron and the reference parton. If the collinear group refers to the initial state of the process, the reference hadron coincides with the initial hadron and the reference parton is the parton confined inside it that is struck in the hard scattering; if the collinear group refers to the final state, the reference hadron is the detected hadron and the reference parton is the fragmenting parton, i.e. the particle that initiates the hadronization process.

• 3. Soft part It embeds the contribution due to the soft gluon radiation that connects the collinear parts and that flows through the detector. It contains soft divergences and carries non-perturbative information, therefore it cannot be computed in perturbation theory. It cannot be directly extracted from data, either, as the energy of the soft radiation is so low that detectors are not sensitive to it. Since the collinear parts interact among each other only through soft gluons, their contribution can affect the cross section in a non-trivial way. Moreover, the soft part is always associated with the collinear terms and there is no way to extract them separately. This is sometimes referred to as the soft factor problem.

Fig. 1

a Pictorial representation of (the hadronic part of) a DIS process. The struk quark is associated to the collinear part relative to the target hadron, while the radiated gluon is the hard real emission. b Pictorial representation of (the hadronic part of) $e^+e^-\rightarrow H \, X$. The quark line corresponds to the fragmenting quark associated to the collinear part representing the final hadron, while the radiated gluon is the hard real emission

In several cases the contribution of the soft part is trivial. In particular, any time in addition to the collinear partons there are real emissions with hard transverse momentum, the soft factor fully factorizes but its value reduces to unity. In these cases, in fact, the soft gluons are kinematically overpowered and do not correlate the collinear parts anymore: in this way each collinear cluster of partons is totally independent from any other. Technically speaking, in such a situation the soft factor involves an integration over all the components of the total soft momentum so that the soft information is washed out in the integral. Whether there could be a hard real emission or not is determined by kinematics. Hence, it is the hard factor that discriminates among different cases.

Kinematical configurations where hard real emissions are present, see for instance Fig. 1, represent instances in which collinear factorization holds. In these cases it is possible to relate each collinear part with a Parton Distribution Function (PDF) or a Fragmentation Function (FF), depending on whether the associated reference hadron is in the initial or in the final state, respectively. As an example of a collinearly factorized process, one could consider the case of an $e^+e^-$ scattering where two spinless hadrons $H_A$ and $H_B$ are produced in the final state, in a configuration far from being back-to-back in the center of mass frame (which, in this case, corresponds to the lab frame). The resulting cross section is given by (see Eq. (12.84) in Ref. [2]):

$$\begin{aligned}&\frac{d \sigma }{\left( \frac{d^3 {{\varvec{p}}}_A}{E_A}\right) \, \left( \frac{d^3 {{\varvec{p}}}_B}{E_B}\right) } = \sum _{j_A,\,j_B} \int \, \frac{d {\widehat{z}}_A}{{\widehat{z}}_A^2} \, d_{H_A / j_A} ({\widehat{z}}_A) \,\nonumber \\&\quad \times \int \, \frac{d {\widehat{z}}_B}{{\widehat{z}}_B^2} \, d_{H_B / j_B} ({\widehat{z}}_B) \, \frac{d {\widehat{\sigma }}}{(\frac{d^3 {{\varvec{k}}}_A}{\epsilon _A}) \, \left( \frac{d^3{{\varvec{k}}}_B}{\epsilon _B}\right) }\,,\nonumber \\ \end{aligned}$$

(1)

where $d {\widehat{\sigma }}$ is the partonic cross section, i.e. the hard part, while $d_{H_i / j_i} ({\widehat{z}}_i)$, for $i = A, \,B$, are the usual FFs associated to the outgoing hadrons, with momenta ${{\varvec{p}}}_A$ and ${{\varvec{p}}}_B$, and to the fragmenting partons of flavor $j_A$ and $j_B$, corresponding to the two collinear contributions to the cross section of the process.

Configurations in which kinematics forbid hard real emissions, instead, are extremely complex, but still very interesting. Here, the soft factor does not reduce to unity, and soft gluons have a non-trivial impact on the cross section as they correlate the collinear parts. This correlation originates from momentum conservation laws in the transverse direction. In fact, with no hard real emissions and consequently no large transverse momentum entering into the game, the low transverse momentum components of soft and collinear particles cannot be neglected anymore: the information regarding the (total) soft transverse momentum survives and the soft factor results in an integration over the plus and minus (but not over the transverse) components. In these cases it is not possible to associate a PDF or a FF to the collinear contributions: parton densities are now related to different and more general objects, known as Transverse Momentum Dependent (TMD) parton functions, either TMD PDFs or TMD FFs depending on whether they refer to an initial or a final state hadron. In this cases collinear factorization breaks and a different, more involved, factorization scheme has to be applied, commonly referred to as TMD factorization. As an example of a TMD factorized process, we can once again consider the production of two spinless hadrons from an $e^+e^-$ scattering where, this time, the two hadrons are almost back-to-back in the $e^+e^-$ center of mass frame. In this case, there are no hard real emissions and the hadronic part of the cross section is given by (see Eq. (13.31) in Ref. [2]):

$$\begin{aligned}&W^{\mu \,\nu }(Q,\,p_A,\,p_B) \nonumber \\&\quad = \frac{8 \pi ^3 z_A z_B}{Q^2} \, \sum _f H^{\mu \,\nu }_{f,\,{\overline{f}}}(Q) \,\nonumber \\&\qquad \times \,\int d^2{{\varvec{k}}}_{A, \,h \,T} \,d^2{{\varvec{k}}}_{B, \,h \,T} \, {\mathbb {S}}({{\varvec{q}}}_{h \,T} - {{\varvec{k}}}_{A, \,h \,T} - {{\varvec{k}}}_{B, \,h \,T}) \, \nonumber \\&\qquad \times D_{H_A / f} ({{\varvec{k}}}_{A, \,h \,T}) \, D_{H_B / {\overline{f}}} ({{\varvec{k}}}_{A, \,h \,T}), \end{aligned}$$

(2)

where $H^{\mu \,\nu }_{f,\,{\overline{f}}}(Q)$ is the hard part, ${\mathbb {S}}$ represents the soft factor and the functions $D_{H_i / f}$, for $i = A, \,B$, are the TMD FFs associated to the outgoing hadrons and to the fragmenting partons of flavor f and ${\bar{f}}$.

Once again, it is kinematics that determines which fatorization scheme has to be used: if the two hadrons are back-to-back then TMD factorization, Eq. (2), must be applied, otherwise Collinear factorization, Eq. (1), will be appropriate.

2 Soft factor

In this section we focus on those kinematical configurations in which there is no hard radiation, where TMD factorization has to be applied and the soft factor plays a non trivial role.

From the point of view of the soft gluons, each collinear group is simply a bunch of particles strongly boosted in a certain direction. The boost is so strong that the soft gluons are only sensitive to the color charge and to the direction of the collinear particles. As a consequence, the propagation of collinear particles is well approximated by a Wilson line, see Appendix A, in the direction of the corresponding collinear group, usually represented by double lines, see Eq. (3). In the massless limit, the versor which identifies this direction is light-like. However, a light-like Wilson line brings unregulated rapidity divergences. In order to cancel them, it is common to introduce a rapidity cut-off $y_i$ which tilts the corresponding Wilson line away from its original light-like direction. Obviously, the final result for the cross section should not depend on these rapidity cut-offs, which then have to be removed in the final stage of the computation. As explained in Ref. [2], the self-interactions of these Wilson lines should not be included into the definition of the soft factor. If ${{\varvec{k}}}_{S,\,T}$ is the total transverse momentum of the real soft radiation flowing through the detector, then the soft factor of a generic process is defined as [2]:

(3)

where D is the dimension of space time ($D = 4 - 2 \epsilon$ in dimensional regularization), $\mu$ is the renormalization scale and $\{y_i\}$ are Lorentz invariant combinations of the rapidity cut-offs ($i = 1, \dots N$ where N is the number of collinear parts in the process). The dependence on the parton-types $j_1 \dots j_N$ of the partons associated to the collinear parts is only on their Wilson line approximation, which changes according to their color representation (fermions or gluons). The label “NO S.I.” reminds us not to consider the Wilson lines self energies Ref. [2]. This implies that $N = 1$ is excluded, since it would correspond only to a Wilson line self energy-like contribution. Finally, the factor $Z_S$ is a UV-renormalization factor that cancels, order by order, the UV divergences generated when the integration region stretches outside of the soft region. The role of $Z_S$ will become clear later on, when the soft factor will be defined in the Fourier conjugate space, Eq. (4).

It is important to stress that, with this definition, the soft factor is sensitive to the number $N \ge 2$ of collinear groups, each one associated to a reference hadron h. Therefore it is not totally blind to the rest of the process, but carries some residual information about the overall process. For this reason, in what follows we will always add a label “N-h” to the soft factor ${\mathbb {S}}$ in order to take into account this dependence.

It is usually more convenient to define the soft factor in the Fourier conjugate ${{\varvec{b}}}_T$ space of ${{\varvec{k}}}_{S,\,T}$, where the quantities involved in the cross section can be identified through an operator definition. In the following, the Fourier transformed quantities will be labeled by a tilde. In particular, the Fourier transform of the soft factor, $\widetilde{{\mathbb {S}}}_{\mathrm{{N}-h}}$, is a matrix in color space, given by the vacuum expectation value of a product of Wilson lines:

$$\begin{aligned}&\widetilde{{\mathbb {S}}}_{\mathrm{{N}-h}}({{\varvec{b}}}_{T};\, \mu ,\,y_i,\,j_k)\nonumber \\&\quad = \int d^{D-2} {{\varvec{k}}}_{S,\,T} \, e^{i \, {{\varvec{k}}}_{S,\,T} \cdot {{\varvec{b}}}_{T}} \,{\mathbb {S}}_{\mathrm{{N}-h}}({{\varvec{k}}}_{S,\,T},\, \mu ,\,y_i,\,j_k) \nonumber \\&\quad = Z_S(\mu ,\,\,y_i,\,j_k) \langle 0 | \prod _{i = 1}^N \, W_{j_i}(\infty ,\,-{{{\varvec{b}}}_T}/{2}; \, n_i(y_i)\,)^\dagger \,\nonumber \\&\qquad \times \, \prod _{k = 1}^N \, W_{j_k}(\infty , \,{{{\varvec{b}}}_T}/{2};\, n_j(y_j)\,) | 0 \rangle \, \vert _{\mathrm{NO S.I.}}. \end{aligned}$$

(4)

The Wilson line $W_{j_i}( \infty ,\,{{{\varvec{b}}}_T}/{2};\, n_i)$ goes from ${{{\varvec{b}}}_T}/{2}$ towards infinity in the direction of $n_i$, which is not light-like thanks to the rapidity cut-off $y_i$, and has the color representation given by the parton type $j_i$.

We can obtain more information about the soft factor by studying its structure in detail. Since all the collinear information is replaced by spinless eikonal propagators, ${{\varvec{k}}}_{S,\,T}$ is the only vector appearing in the soft factor. Therefore, ${\mathbb {S}}$ is always rotational invariant and depends only on the modulus $|{{\varvec{k}}}_{S,\,T}| = k_{S,\,T}$. This reflects on the Fourier conjugate space, where the dependence on ${{\varvec{b}}}_T$ is only through its modulus $|{{\varvec{b}}}_T| = b_T$. Moreover the natural leading momentum region of ${\mathbb {S}}$ is where all the momenta are soft, with components of size $\lambda _S = {\lambda ^2}/{Q}$, where $\lambda<< Q$ is a very low energy scale. When the soft factor is Fourier transformed, the total transverse soft momentum ${{\varvec{k}}}_{S,\,T}$ is integrated out and its dependence is replaced by ${{\varvec{b}}}_T$. At fixed $b_T$ we can roughly access all momenta with $k_{S,\,T} \le \frac{1}{b_T}$, hence this operation can be regarded as a sort of analytic continuation of the function ${\mathbb {S}}_{\mathrm{{2}-h}}(k_{S,\,T})$ outside of its natural momentum region, since when $b_T$ is small $k_{S,\,T}$ can be very large. This generates UV divergences which will have to be canceled order by order by the UV counterterm $Z_S$.

Fig. 2

Leading regions for the soft factor, ${\mathbb {S}}_{\mathrm{{N}-h}}$, at small $b_T$

The application of the factorization procedure to the soft factor itself gives us the possibility to express $\widetilde{{\mathbb {S}}}_{\mathrm{{N}-h}}$ in terms of perturbative and non-perturbative parts (see Ref. [2]). Leading regions involve hard, collinear and soft subgraphs, as represented in Fig. 2. The hard factor is associated to the external Wilson line vertices and contains hard subgraphs with highly virtual loops. There is a collinear subgraph corresponding to each Wilson line and all of them are connected by the soft subgraph. Furthermore, if the entering transverse momentum $k_s$ is large enough, there can be more hard subgraphs ${\mathbb {C}}_{\alpha }$ with production of final-state jets of high transverse momentum (i.e. hard gluon emissions which cross the cut). In each hard jet there is a fully inclusive sum/integral over final states, hence the sum-over-cuts argument presented in Ref. [2] allows us to consider them as being far off-shell and part of the hard factor. In this case collinear factorization holds and the soft factor is unity. Furthermore, there is no convolution between the hard part and the collinear factors ${\mathbb {C}}_i$, since the cut eikonal propagators that exit from the hard subgraphs do not carry momentum. As a consequence, all the collinear parts are integrated 2-h soft factors and are unity as well. Therefore, the only remaining effective region is the hard factor with all the extra hard jets. It has the same structure of ${\mathbb {S}}_{\mathrm{{N}-h}}$ but now it is fully computable in perturbation theory. In particular, it is a standard result that general soft functions exponentiate and that the exponent can be computed by using web technology, see for example Refs. [6,7,8,9]. Hence at small $b_T$ the soft factor can be written schematically as:

$$\begin{aligned}&{\mathbb {S}}_{\mathrm{{N}-h}} (b_T;\, \mu ,\,y_i,\,j_k) \overset{\hbox {low }b_T}{\sim }\nonumber \\&\quad \int d^{D-2} {{\varvec{k}}}_{S,\,T} \, e^{i \, {{\varvec{k}}}_{S,\,T} \cdot {{\varvec{b}}}_{T}} \, \nonumber \\&\quad \hbox {exp}\left[ \sum _{{\mathcal {W}}} {\mathcal {W}}(k_{S,\,T};\, \mu ,\,y_i,\,j_k) \right] , \end{aligned}$$

(5)

where the sum is extended to all the (multiparton) webs and the sums over the diagrams in each web, corresponding to a certain color mixing matrix, has not been shown for simplicity. In order to separate the small and large $b_T$ behavior of $\widetilde{{\mathbb {S}}}_{\mathrm{{N}-h}}$, we can modify its functional dependence on ${{\varvec{b}}}_T$ by introducing a function ${{\varvec{b}}}_T^\star ({{\varvec{b}}}_T)$ such that it coincides with ${{\varvec{b}}}_T$ at small $b_T$, while at large $b_T$ it is no larger than a certain $b_{\mathrm{{max}}}$. A possible choice, according to Ref. [2, 4, 5], is given by:

$$\begin{aligned} {{\varvec{b}}}_T^\star \left( b_T\right) = \dfrac{{{\varvec{b}}}_T}{\sqrt{1 + {b_T^2}/{b_{\mathrm{{max}}}^2}}} \end{aligned}$$

(6)

Then, by dividing and multiplying $\widetilde{{\mathbb {S}}}_{\mathrm{{N}-h}}$ in Eq. (5) by its small $b_T$ behavior, we easily obtain a factorized expression which holds valid at any value of $b_T$:

$$\begin{aligned} \widetilde{{\mathbb {S}}}_{\mathrm{{N}-h}}(b_T; \, \mu ,\,y_i,\,j_k)&= \widetilde{{\mathbb {S}}}_{\mathrm{{N}-h}}(b_T^\star ; \, \mu ,\,y_i,\,j_k) \, \nonumber \\&\quad \times \, \frac{\widetilde{{\mathbb {S}}}_{\mathrm{{N}-h}}(b_T; \, \mu ,\,y_i,\,j_k)}{\widetilde{{\mathbb {S}}}_{\mathrm{{N}-h}}(b_T^\star ; \, \mu ,\,y_i,\,j_k)}\nonumber \\&= \, \int d^{D-2} {{\varvec{k}}}_{S,\,T} \, e^{i \, {{\varvec{k}}}_{S,\,T} \cdot {{\varvec{b}}}_T^\star } \, \nonumber \\&\quad \times \hbox {exp} \left[ \sum _{{\mathcal {W}}} {\mathcal {W}}(k_{S,\,T};\, \mu ,\,y_i,\,j_k) \right] \nonumber \\&\quad \times M_S (b_T; \, \mu ,\,y_i,\,j_k), \end{aligned}$$

(7)

where $M_S (b_T; \, \mu ,\, \{y_i\}_{i = 1 \dots N}, \,\{j_i\}_{i = 1 \dots N})$ is the fully non-perturbative function that models the N-h soft factor at large $b_T$, while the whole perturbative content is gathered in the webs.

In the t’Hooft limit,¹ the soft factor is strongly simplified. Regarding the perturbative part, the only surviving diagrams are planar and the exponentiation becomes trivial. Furthermore, we can also make some guess on the non-perturbative part which is, in principle, a fully arbitrary function, since there is no way to extract it independently from experiments. In this limit the non-perturbative contribution of $\widetilde{{\mathbb {S}}}_{\mathrm{{N}-h}}$ only regards the incoherent emission of free glueballs, of every possible kind.² The function that models this kind of emission is a Poisson distribution, similarly to what happens for photons in QED.

2.1 2-h soft factor

In the 2-h class, there are two directions for the collinear parts which can be identified to the plus and the minus direction in the c.m. frame. The Wilson lines are tilted with respect to these light-like directions by introducing two rapidity cut-offs $y_1$ and $y_2$. The original plus and minus directions are restored if the cut-offs are removed, i.e. by taking the limits $y_1 \rightarrow +\infty$ and $y_2 \rightarrow -\infty$. In total, there are four Wilson lines, two on each side of the final state cut. The only relevant case for applications involves Wilson lines that replace fermionic collinear partons. Hence, in the following we will drop the dependence on the parton types for simplicity. Furthermore, the 2-h soft factor is color singlet, proportional to the identity matrix in color space, i.e. $( \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} )^i_{\;j} \propto \delta ^i_{\;j}$. Then, $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ is defined as the coefficient in front of the delta function. By using the definition in Eq. (4) we have:

$$\begin{aligned}&\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}({{\varvec{b}}}_{T};\, \mu , \,y_1 - y_2) \nonumber \\&\quad = Z_S(\mu ,\,y_1 - y_2) \nonumber \\&\quad \quad \times \, \frac{\hbox {Tr}_C}{N_C} \, \langle 0 | W(-{{{\varvec{b}}}_T}/{2}, \,\infty ; \, n_1(y_1)\,)^\dagger \, \nonumber \\&\quad \quad \times \, W({{{\varvec{b}}}_T}/{2}, \,\infty ; \, n_1(y_1)\,) \, \nonumber \\&\quad \quad \times \, W({{{\varvec{b}}}_T}/{2}, \,\infty ; \, n_2(y_2)\,)^\dagger \, \nonumber \\&\quad \quad \times \, W(-{{{\varvec{b}}}_T}/{2}, \,\infty ; \, n_2(y_2)\,) | 0 \rangle \, \vert _{\mathrm{NO S.I.}}, \end{aligned}$$

(8)

where $N_C$ is the number of colors available for quarks and antiquarks (3 in QCD). The Eq. (8) describes a loop for the full path outlined by the Wilson lines. It starts (e.g.) from $-{{{\varvec{b}}}_T}/{2}$ and goes to ${{{\varvec{b}}}_T}/{2}$, passing through $\infty$, along the almost plus direction $n_1$. Then it comes back, again passing through $\infty$, along the almost minus direction $n_2$. Notice that the only Lorentz invariant combination for a function depending on two rapidities (e.g. $y_1$ and $y_2$) is their difference (e.g. $y_1 - y_2$). It is possible to write the evolution equation for ${\mathbb {S}}_{\mathrm{{2}-h}}$ in the $b_T$-space with respect to both rapidity cut-offs, $y_1$ and $y_2$, using a single rapidity-independent kernel ${\widetilde{K}}({{\varvec{b}}}_T;\,\mu )$ defined as [2]:

$$\begin{aligned}&\lim _{y_2 \rightarrow -\infty } \frac{\partial \log {\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}({{\varvec{b}}}_{T}; \, \mu , \, y_1 - y_2)}}{\partial y_1} = \frac{1}{2}\, {\widetilde{K}}({{\varvec{b}}}_T;\,\mu ) \end{aligned}$$

(9)

$$\begin{aligned}&\lim _{y_1 \rightarrow +\infty } \frac{\partial \log {\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}({{\varvec{b}}}_{T}; \, \mu , \, y_1 - y_2)}}{\partial y_2} = - \frac{1}{2}\, {\widetilde{K}}({{\varvec{b}}}_T;\,\mu ) \,, \end{aligned}$$

(10)

It has an anomalous dimension $\gamma _K$:

$$\begin{aligned} \frac{d {\widetilde{K}}({{\varvec{b}}}_T;\,\mu )}{d \log {\mu }} = - \gamma _K(\alpha _s(\mu )), \end{aligned}$$

(11)

where $\gamma _K$ depends on $\mu$ through the strong coupling $\alpha _S$ and is independent of $b_T$. Then, ${\widetilde{K}}$ can be written as:

$$\begin{aligned} {\widetilde{K}}({{\varvec{b}}}_T;\,\mu ) = {\widetilde{K}}({{\varvec{b}}}_T;\,\mu _0) - \int _{\mu _0}^\mu \, \frac{d \mu '}{\mu '} \, \gamma _K(\alpha _s(\mu ')) . \end{aligned}$$

(12)

For large values of $(y_1 - y_2)$, the solution to the evolution equations for $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ is given by:

$$\begin{aligned}&\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}({{\varvec{b}}}_T; \, \mu , \, y_1 - y_2) \nonumber \\&\quad = \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}({{\varvec{b}}}_{T}; \, \mu _0, \, 0) \, \nonumber \\&\qquad \times \hbox {exp} \Big \{ \frac{y_1 - y_2}{2} \, {\widetilde{K}}({{\varvec{b}}}_T;\,\mu )\Big \}\nonumber \\&\qquad + {\mathcal {O}}\Big ( e^{-(y_1 - y_2)} \Big ) \nonumber \\&\quad = \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}({{\varvec{b}}}_{T}; \, \mu _0, \, 0) \, \nonumber \\&\qquad \times \hbox {exp} \Big \{ - \frac{y_1 - y_2}{2} \nonumber \\&\qquad \times \Big [ \int _{\mu _0}^{\mu }\, \frac{d \mu '}{\mu '}\, \gamma _K(\mu ) \; - {\widetilde{K}}({{\varvec{b}}}_T;\,\mu _0)\Big ]\Big \} \, \nonumber \\&\qquad +{\mathcal {O}}\Big ( e^{-(y_1 - y_2)} \Big ), \end{aligned}$$

(13)

where the reference values of the RG scale and of the rapidities are chosen to be $\mu _0$ and $y_{1,\,0} = y_{2,\,0}$, respectively. In the solution of the evolution equation, two functions appear: the fixed scale soft factor $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}(b_T; \, \mu _0, \, 0)$ and the soft kernel ${\widetilde{K}}(b_T;\,\mu )$. Both of them can be separated in terms of their perturbative and non-perturbative contents by using the $b^\star$ prescription, similarly to what was done in Eq. (7):

$$\begin{aligned}&\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}(b_T; \, \mu _0, \, 0) = \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}(b_{T}^\star ; \, \mu _0, \, 0) \, M_S^{(0)}(b_T) \,; \end{aligned}$$

(14)

$$\begin{aligned}&{\widetilde{K}}(b_T;\,\mu ) = {\widetilde{K}}(b_T^\star ;\,\mu ) - g_K(b_T) \,. \end{aligned}$$

(15)

Finally, consistency between Eqs. (7) and (13) requires that:

$$\begin{aligned}&\lim \limits _{\begin{array}{c} y_1 \rightarrow +\infty \\ y_2 \rightarrow -\infty \end{array}} \, \int d^{D-2} {{\varvec{k}}}_{S,\,T} \, e^{i \, {{\varvec{k}}}_{S,\,T} \cdot {{\varvec{b}}}_T^\star } \, \nonumber \\&\qquad \times \hbox {exp} \left[ \sum _{{\mathcal {W}}} {\mathcal {W}}(k_{S,\,T};\, \mu ,\,y_1-y_2) \right] \nonumber \\&\quad = \frac{y_1 - y_2}{2}{\widetilde{K}}(b_T^\star ;\,\mu ) = \frac{y_1 - y_2}{2} \nonumber \\&\qquad \times \left[ {\widetilde{K}}(b_T^\star ;\,\mu _0) - \int _{\mu _0}^{\mu }\, \frac{d \mu '}{\mu '}\, \gamma _K(\mu ') \right] ; \end{aligned}$$

(16)

$$\begin{aligned}&\lim \limits _{\begin{array}{c} y_1 \rightarrow +\infty \\ y_2 \rightarrow -\infty \end{array}} \, M_S(b_T; \, \mu ,\,y_1 - y_2) = M_S^{(0)}(b_T) \, e^{-\frac{y_1 - y_2}{2} \, g_K(b_T)} \,;\nonumber \\&\quad \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}(b_T^\star ; \, \mu _0, \, 0) = 1 \,. \end{aligned}$$

(17)

Notice that the non-perturbative function $M_S(b_T; \, \mu ,\,y_1 - y_2)$ loses its dependence on $\mu$ in the large rapidity limit, as $g_K$ does not depend on the RG scale. Since we are only interested in the asymptotic behaviour of $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$, we will drop the label (0) from $M_S(b_T)$ and we will refer to it as the soft model, i.e. the non-perturbative part which will have to be parametrized and treated phenomenologically, possibly taking inspiration from the properties of the soft factor in the t’Hooft limit. The two non-perturbative functions $M_S$ and $g_K$ should not contribute at small $b_T$ by definition, hence we require that $g_K(b_T)\rightarrow 0$ and $M_S(b_T) \rightarrow 1$ when $b_T \rightarrow 0$. Furthermore, since the Fourier transform of $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ has to be well behaved, the contribution of $g_K$ and $M_S$ should be suppressed at large $b_T$. Notice that the factor in front of $g_K$, being proportional to the difference of the rapidity cut-offs, is always large and negative in the large rapidity cut-off limit. In conclusion, the 2-h soft factor in $b_T$ space can be written as:

$$\begin{aligned}&\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}(b_T; \, \mu ,\,y_1 - y_2) \nonumber \\&\quad = e^{\frac{y_1 - y_2}{2}{\widetilde{K}}(b_T^\star ;\,\mu )} \, M_S(b_T) \, \nonumber \\&\qquad e^{-\frac{y_1 - y_2}{2} \, g_K(b_T)} + {\mathcal {O}}\Big ( e^{-(y_1 - y_2)} \Big ) \,. \end{aligned}$$

(18)

This result shows that the soft factor itself can be factorized in a purely perturbative part, process dependent but calculable within pQCD, and a part which is genuinely non perturbative and, inevitably, will have to be committed to a phenomenological model, in this case embedded in the functions $M_S(b_T)$ and $g_K(b_T)$.

Although the definition of Eq. (8) implies that $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} = 1$ at $b_T = 0$, a direct fixed order perturbative computation of ${\widetilde{K}}$ does not reproduce the correct behavior in this region. In this regard, since the soft factor is unity at $b_T = 0$, then ${\widetilde{K}}$ goes to zero at small $b_T$, but an explicit calculation gives instead a larger and larger value as $b_T$ decreases, forcing $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ to vanish in $b_T = 0$. This kind of problems arise because the integrated soft factor can be defined through perturbative QCD only as a bare quantity. A solution can be found by applying some regularization procedure, for instance one can modify the $b^\star$ prescription of Eq. (6) allowing for the introduction of a new parameter $b_{\mathrm{MIN}} \ne 0$ that provides a minimum value for $b_T$ (see Appendix B).

3 Collinear Parts and TMDs

Let’s now consider a generic collinear part. If kinematics forbid hard real emissions, the information about the transverse momentum ${{\varvec{k}}}_{T}$ of the reference parton survives. All the collinear particles are boosted very strongly in the collinear group direction, that we can identify with the plus direction without loss of generality. To them, everything outside of the collinear group is moving very fast in the opposite direction, so fast that the only surviving information is the color charge and the direction. In other words, as seen from the collinear factor, the rest of the process is well approximated by a light-like Wilson line flowing in the direction opposite to that of the collinear group.

Assuming for simplicity that the reference parton is a quark, if ${{\varvec{k}}}_T$ is the total transverse momentum of the collinear group, then the collinear factor (along the plus direction) is defined as in Ref. [2]:

(19)

where the color average ${\hbox {Tr}_{\mathrm{C}}}/{N_C}$ is due to the fact that collinear factors are color singlets and, analogously to the 2-h soft factor, they are defined as the coefficient in front of the delta in color space. The variable $\xi$ is the light-cone fraction of the momentum k of the reference parton, quark of flavor j, with respect to the momentum P of the reference hadron H, $\mu$ is the renormalization scale at which ${\mathbb {C}}$ is evaluated and $y_P$ is the (very large) rapidity of the reference hadron. The definition of $\xi$ in the initial and final state is given by:

$$\begin{aligned} \xi = {\left\{ \begin{array}{ll} x = \frac{k^+}{P^+} &{}\hbox { initial state hadron;} \\ z = \frac{P^+}{k^+} &{}\hbox { final state hadron,} \end{array}\right. } \end{aligned}$$

(20)

To the Wilson line, instead, we can associate a very large and negative rapidity. Similarly to the definition in Eq. (3), $Z_C$ is the UV-counterterm of ${\mathbb {C}}$, while the label “NO S.I.” reminds us not to consider the Wilson lines self interactions. It is important to stress that the collinear factor ${\mathbb {C}}$ is totally blind to the rest of the process, it only depends on its intrinsic variables. As for the soft factor, the operator definition is simpler in the Fourier conjugate space. Here, we have:

$$\begin{aligned}&{\widetilde{{\mathbb {C}}}}_{j,\,H}(\xi ,\, {{\varvec{b}}}_{T}; \, \mu , \, y_P ,\,-\infty ) = \int d^{D-2} {{\varvec{k}}}_T \, e^{i \, {{\varvec{k}}}_T \cdot {{\varvec{b}}}_T} \,{\mathbb {C}}_{j,\,H}(\xi ,\, {{\varvec{k}}}_{T}; \, \mu , \, y_P,\,-\infty ) \nonumber \\&\quad = Z_C(\mu , y_P , -\infty ) \, \frac{\hbox {Tr}_{\mathrm{C}}}{N_C} \int \frac{d x^-}{2 \pi } e^{i k^+ x^-} {\left\{ \begin{array}{ll} \langle P \,(H) | {\overline{\psi }}(-{x}/{2}) W_j(-{x}/{2}, {x}/{2};\, w_{-}) \psi ({x}/{2})| P (H) \rangle \vert _{\hbox {NO S.I.}} &{}\hbox { initial state,} \\ \begin{aligned} \dfrac{1}{\xi } &{}\times \, \sum \limits _{X} \, \langle P\,(H), X;~\hbox {out} | {\overline{\psi }}(-{x}/{2}) W_j(-{x}/{2}, \infty ;\, n_1(y_1)\,)^\dagger | 0 \rangle \\ &{}\quad \times \, \langle 0 | W_j({x}/{2}, \infty ;\, w_{-}\,) \, \psi ({x}/{2})| P (H), X;~\hbox {out} \rangle \vert _{\hbox {NO S.I.}} \end{aligned} \, &{}\hbox { final state,} \end{array}\right. } \nonumber \\ \end{aligned}$$

(21)

where $x = (0, x^-, {{\varvec{b}}}_T)$ and $w_{-}$ is the light-like minus direction of the Wilson line.

The presence of a light-like Wilson line allows for particles with a low, or even a very large negative rapidity, to be considered as part of the collinear group. This contradiction reflects in the computation by inducing the presence of unregulated rapidity divergences. This problem can be solved by subtracting out these unphysical contributions from the collinear factor by using the subtraction method described in Ref. [2]. Since all the non truly collinear contributions are due to the overlapping with the soft region, they can be rearranged in one global term which turns out to be a 2-hadron soft factor. Hence we can use the definitions given in Eqs. (3) and (4) to subtract them out and obtain:

$$\begin{aligned}&{\widetilde{{\mathbb {C}}}}_{j,\,H}^{\mathrm{sub}} (\xi ,\, {{\varvec{b}}}_{T}; \, \mu , \, y_P - y_1) \nonumber \\&\quad = Z_C^{\mathrm{sub}}(\mu ,\, y_P - y_1) \,Z_2 \left( \alpha _S(\mu )\right) \nonumber \\&\qquad \times \lim _{y_{u_2} \rightarrow -\infty } \dfrac{{\widetilde{{\mathbb {C}}}}_{j,\,H}^{(0)} (\xi , \,{{\varvec{b}}}_{T}; \, \mu , \, y_P - y_{u_2})}{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}^{(0)} (b_{T}; \, \mu , \,y_1 - y_{u_2})} \,, \end{aligned}$$

(22)

where $y_1$ is the rapidity cut-off carried by ${\mathbb {S}}_{\mathrm{{2}-h}}$ that, similarly to the case of the soft factor, should be removed in the final result for the cross section. After subtraction, the particles in ${\mathbb {C}}$ can only have a rapidity y such that $y_1< y < y_P \sim +\infty$. Hence if $y_1$ is chosen to be sufficiently large, only strongly boosted particles in the plus direction contribute to ${\mathbb {C}}$, according to the naive physical intuition. The subtracted collinear factor has its own UV-counterterm $Z_C^{\mathrm{sub}}$, for this reason in the previous definition the quantities inside the limit are bare, in the sense that they have to be considered without their UV-renormalization factors. Since the unsubtracted collinear part is defined with renormalized quark fields, see Eq. (21), then if $Z_C^{\mathrm{sub}}$ is the ratio of the renormalized collinear part to the unrenormalized collinear part, we have to multiply explicitly by the wave-function renormalization factor of the quark field, $Z_2$.

Having given the general definition of ${\mathbb {C}}$, TMDs can be obtained straightforwardly. In fact, as ${\mathbb {C}}$ is an operator acting onto the space of Dirac spinors, it belongs to the Clifford algebra built from the Dirac matrices $\left\{ \gamma ^\mu \right\}$. Therefore, we simply expand ${\mathbb {C}}$ on the basis of this algebra. Neglecting all the dependences on partonic and hadronic variables, we have:

$$\begin{aligned} {\mathbb {C}}^{\mathrm{sub}}= & {} {\mathcal {S}} \, {\mathbb {I}} + {\mathcal {V}}^\mu \, \gamma _\mu + {\mathcal {A}}^\mu \, \gamma ^5 \, \gamma _\mu \nonumber \\&+ i \, {\mathcal {P}} \, \gamma ^5 + i {\mathcal {T}}^{\mu \, \nu } \, \sigma _{\mu \, \nu } \gamma ^5. \end{aligned}$$

(23)

Then, the TMDs are related to the coefficients ${\mathcal {S}},{\mathcal {V}}, \dots , {\mathcal {T}}^{\mu \, \nu }$ of the Clifford Algebra expansion and the definition in Eq. (22) naturally extends to TMDs. Such coefficients can be further expanded in terms of all the Lorentz tensors contributing to the leading twist approximation (see e.g. Ref. [10]). This allows to isolate all the dependence on the vector part of ${{\varvec{b}}}_T$ in the coefficients of such expansion, leaving a set of scalar functions depending only on the modulus $b_T$. These scalar functions are the TMDs. For example, the coefficient of $\gamma ^+$ defines the unpolarized TMDs and the Sivers function:

$$\begin{aligned} {\mathcal {V}}^+= & {} {1}/{4}\, \hbox {Tr}_{\mathrm{Dirac}} \, \bigl [ \gamma ^+ {\mathbb {C}}^{\mathrm{sub}} \bigr ] \nonumber \\= & {} {\left\{ \begin{array}{ll} f_1 - \frac{1}{M} |{{\varvec{S}}}_T \times {{\varvec{k}}}_T| f_{1 T}^\perp \, &{}\hbox { initial state} \,,\\ D_1 - \frac{1}{M} |{{\varvec{S}}}_T \times {{\varvec{k}}}_{h,T}| D_{1 T}^\perp \, &{}\hbox { final state}\,. \end{array}\right. } \end{aligned}$$

(24)

Formally, if C is a generic TMD function referring to a collinear factor in the plus direction, then its definition equipped with subtractions is inherithed directly from Eq. (22) and it is given by:

$$\begin{aligned}&{\widetilde{C}}_{j,\,H}^{\mathrm{sub}} (\xi ,\, b_T; \, \mu , \, y_P - y_1) \nonumber \\&\quad = \left( Z_{\mathrm{TMD}}\right) _j(\mu ,\, y_P - y_1) \, Z_2 \left( \alpha _S(\mu )\right) \nonumber \\&\qquad \times \lim _{y_{u_2} \rightarrow -\infty } \dfrac{\frac{1}{4} \hbox {Tr}_{\mathrm{Dirac}} \left[ \varGamma \, {\widetilde{{\mathbb {C}}}}_{j,\,H}^{(0)} (\xi , \,{{\varvec{b}}}_{T}; \, \mu , \, y_P - y_{u_2})\right] ^{\begin{array}{c} \hbox {leading}\\ \hbox {twist coeff.} \end{array}} }{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}^{(0)} (b_T; \, \mu , \,y_1 - y_{u_2})}\, \nonumber \\&\quad = \frac{1}{4} \hbox {Tr}_{\mathrm{Dirac}} \left[ \varGamma \, {\widetilde{{\mathbb {C}}}}_{j,\,H}^{\mathrm{sub}} (\xi , \,{{\varvec{b}}}_{T}; \, \mu , \, y_P - y_{u_2})\right] ^{\begin{array}{c} \mathrm{leading}\\ \hbox {twist coeff.} \end{array}}, \end{aligned}$$

(25)

where $\varGamma$ is the proper Dirac matrix combination to extract the desidered TMD and $Z_{\mathrm{TMD}}$ is its own UV counterterm. The label “leading twist coeff.” means that the TMDs are obtained, after the projection onto the Clifford Algebra, as the coefficients of the expansion at leading twist. The operator definition of TMD as given in Eq. (25), which follows directly from the TMD factorization prescription, will be referred to as the factorization definition. Notice that within this definition, the TMD is a purely collinear object, as all soft sub-divergences have been subtracted out.

3.1 Evolution equations for TMDs

In the factorization definition³ of TMDs, Eq. (25), a 2-h soft factor appears as a consequence of the subtraction mechanism. Therefore, we can use the results of Sect. 2.1 to write the evolution equation (Collins-Soper evolution) for ${\widetilde{C}}$ with respect to the rapidity cut-off $y_1$. On the other hand, the evolution with respect to the scale $\mu$ (i.e. the Renormalization Group evolution) is ruled by the anomalous dimension $\gamma _C$. The equations are given by:

$$\begin{aligned}&\dfrac{\partial \log {{\widetilde{C}}_{j,\,H}(\xi ,\, b_{T}; \, \mu , \,\zeta )}}{\partial \log {\sqrt{\zeta }}} = \frac{1}{2} {\widetilde{K}}(b_T;\,\mu ) \,, \end{aligned}$$

(26)

$$\begin{aligned}&\dfrac{\partial \log {{\widetilde{C}}_{j,\,H}(\xi ,\, b_{T}; \, \mu , \, \zeta )}}{\partial \log {\mu }} = \gamma _C \left( \alpha _S(\mu ),\, \frac{\zeta }{\mu ^2} \right) \,, \end{aligned}$$

(27)

which, for later convenience, have been re-written in terms of a new variable, $\zeta$, defined as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} \zeta = \left( M x\right) ^2 \, e^{2 (y_P - y_1)} &{}\hbox { initial state hadron;} \\ \zeta = \left( \frac{M}{z}\right) ^2 \, e^{2 (y_P - y_1)} &{}\hbox { final state hadron,} \end{array}\right. } \end{aligned}$$

(28)

where M is the mass of the reference hadron, while x and z are the light-cone fractions of the momentum of the reference parton with respect to the hadron. Thanks to the definitions in Eq. (20), in both initial and final states we can write $\zeta \sim Q^2 e^{-2 y_1}$. In addition to the previous evolution equations, we also have the RG evolution of ${\widetilde{K}}$, Eq. (11), and the CS evolution of $\gamma _C$, given by:

$$\begin{aligned} \frac{\partial \gamma _C \left( \alpha _S(\mu ),\, {\zeta }/{\mu ^2} \right) }{\partial \log {\sqrt{\zeta }}} = -\frac{1}{2}\gamma _K(\alpha _S(\mu )) , \end{aligned}$$

(29)

which gives:

$$\begin{aligned} \gamma _C \left( \alpha _S(\mu ),\, {\zeta }/{\mu ^2} \right)= & {} \gamma _C \left( \alpha _S(\mu ),\,1 \right) \nonumber \\&- \frac{1}{4} \gamma _K (\alpha _S(\mu )) \log {\frac{\zeta }{\mu ^2}}. \end{aligned}$$

(30)

With the help of Eqs. (11), (29) and (30), we can rewrite the solution to Eqs. (26) and (27) as [2]:

$$\begin{aligned}&{\widetilde{C}}_{j,\,H}(\xi ,\, b_{T}; \, \mu , \,\zeta ) \nonumber \\&\quad = {\widetilde{C}}_{j,\,H}(\xi ,\, b_{T}^\star ; \, \mu _0, \,\zeta _0)\nonumber \\&\quad \quad \times \, \hbox {exp} \Big \{ \frac{1}{4} \, {\widetilde{K}}(b_T^\star ;\,\mu _0) \, \log {\frac{\zeta }{\zeta _0}} + \int _{\mu _0}^{\mu } \frac{d \mu '}{\mu '} \, \nonumber \\&\quad \quad \times \left[ \gamma _C(\alpha _S(\mu '),\,1) - \frac{1}{4} \, \gamma _K(\alpha _S(\mu ')) \, \log {\frac{\zeta }{\mu '^2}}\right] \Big \} \nonumber \\&\quad \quad \times \left( M_C\right) _{j,\,H}(\xi ,\,b_T) \hbox { exp} \Big \{ -\frac{1}{4} \, g_K(b_T) \, \log {\frac{\zeta }{{\overline{\zeta }}_0}} \Big \} \end{aligned}$$

(31)

where the standard choices for the reference values of the scales are⁴:

$$\begin{aligned}&\mu _0 = \mu _b = \frac{2 e^{-\gamma _E}}{b_T^\star } \,; \end{aligned}$$

(32)

$$\begin{aligned}&\zeta _0 = \mu _b^2 \,; \end{aligned}$$

(33)

$$\begin{aligned}&{\left\{ \begin{array}{ll} {\overline{\zeta }}_0 = \left( M x \right) ^2 &{}\hbox { initial state;} \\ {\overline{\zeta }}_0 = \left( \frac{M}{z} \right) ^2 &{}\hbox { final state.} \end{array}\right. } \end{aligned}$$

(34)

In the solution of the evolution equation the $b_T^\star$ prescription, Eq. (6), has been used in order to separate the perturbative from the non-perturbative content, in complete analogy to what was done for the soft factor in Sect. 2.1. In particular, in Eq. (31), the non-perturbative behavior of the TMD is described by two functions. The first is $g_K$, the same function that appears in Eq. (18) in the asymptotic behavior of $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$. The second is the TMD model function $\left( M_C\right) _{j,\,H}(\xi ,\,b_T)$, that embeds the genuine non-perturbative behavior of the TMD: it depends on the flavor of the reference parton and on the reference hadron associated to the collinear part. By definition, the model should not influence the TMD at small $b_T$. Furthermore, since the Fourier transform of the TMD has to be well behaved, the model should be sufficiently suppressed at large $b_T$.⁵ These properties restrict the behaviour of the non-perturbative function $M_C$ at small and large $b_T$ as follows

$$\begin{aligned}&\lim _{b_T \rightarrow 0} M_C(\xi ,\,b_T) = 1;&\lim _{b_T \rightarrow \infty } M_C(\xi ,\,b_T) = 0 . \end{aligned}$$

(35)

The factorization procedure can be applied either to the full collinear factor or to the TMDs themselves, in order to study their behavior at small $b_T$, outside of their natural collinear momentum region. This is given by a convolution of a finite (calculable in perturbative QCD) hard coefficient ${\mathcal {C}}$ with the TMD integrated over ${{\varvec{k}}}_T$. The proof can be found in Chapter 13 of Ref. [2]. Hence, TMDs at small $b_T$ can be written as Operator Product Expansions (OPE):

$$\begin{aligned}&{\widetilde{C}}_{j,\,H}(\xi ,\, b_{T}; \, \mu , \,\zeta ) \overset{\hbox {low }b_T}{\sim } \;\; \widetilde{{\mathcal {C}}}_j^{\,k}(b_{T}; \, \mu , \,\zeta ) \otimes c_{k,\,H} (\mu ) ,\nonumber \\&\quad = {\left\{ \begin{array}{ll} \left( \widetilde{{\mathcal {C}}}_j^{\;k}(b_{T}; \, \mu , \,\zeta ) \otimes f_{k/H}(\mu ) \right) (x) &{}\hbox { initial state;} \\ z^{-2+2\epsilon } \left( d_{H/k} (\mu ) \otimes \widetilde{{\mathcal {C}}}_{\;j}^{k}(b_{T}; \, \mu , \,\zeta ) \right) (z) &{}\hbox { final state.} \end{array}\right. } \end{aligned}$$

(36)

where $\widetilde{{\mathcal {C}}}_j^{\,k}$ are the Wilson Coefficients of the OPE, which are matrices in the flavor space. A sum over k is implicit. In the second line of Eq. (36) we distinguish the Wilson Coefficients of the initial state from those corresponding to the final state according to the position of their upper and lower flavor indices. The convolution $\otimes$ of two generic functions f and g is defined as

$$\begin{aligned} \left( f \otimes g\right) (x) = \int _x^1 \, \frac{d \rho }{\rho } f({x}/{\rho }) g(\rho ) , \end{aligned}$$

(37)

where we recall that the Wilson Coefficients of the final state have a normalization factor $\rho ^{2-2\epsilon }$ when the convolution is made explicit, see Ref. [11]. The integrated TMDs are indicated by lowercase letters. In the following, $c_{k,\,H}$ will be a generic integrated TMD, while f will label integrated TMD PDFs and d will refer to integrated TMD FFs. Thanks to the OPE, the solution of the evolution equations, Eq. (31), can be rewritten as

$$\begin{aligned}&{\widetilde{C}}_{j,\,H}(\xi ,\, b_{T}; \, \mu , \,\zeta ) \nonumber \\&\quad = \left( \widetilde{{\mathcal {C}}}_j^{\;k}(b_{T}^\star ; \, \mu _0, \,\zeta _0) \otimes c_{k,\,H} (\mu _0) \right) (\xi ) \,\nonumber \\&\quad \quad \times \hbox { exp} \Big \{ \frac{1}{4} \, {\widetilde{K}}(b_T^\star ;\,\mu _0) \, \log {\frac{\zeta }{\zeta _0}} + \int _{\mu _0}^{\mu } \frac{d \mu '}{\mu '} \, \nonumber \\&\quad \quad \times \left[ \gamma _C(\alpha _S(\mu '),\,1) - \frac{1}{4} \, \gamma _K(\alpha _S(\mu ')) \, \log {\frac{\zeta }{\mu '^2}}\right] \Big \} \nonumber \\&\quad \quad \times \left( M_C\right) _{j,\,H}(\xi ,\,b_T) \hbox { exp} \Big \{ -\frac{1}{4} \, g_K(b_T) \, \log {\frac{\zeta }{{\overline{\zeta }}_0}} \Big \} . \end{aligned}$$

(38)

The definition of integrated TMDs coincides with the Fourier transformed TMDs in $b_T = 0$. Perturbative QCD fails to give the right result in $b_T = 0$ because of the new UV divergences introduced by the integral over the whole range of $k_T$. In fact, as we explain in Appendix B.2, ${\widetilde{C}}$ goes to zero as $b_T \rightarrow 0$ (see Eq. (134)) and the usual collinear PDFs and FFs are not recovered. This problem is completely analogous to that encountered in Sect. 2.1 and it can be solved in a similar way, by defining a regularization procedure for the definition of the integrated TMDs (see Appendix B).

3.2 Rapidity dilations

In the operator definition of the TMD, Eq. (25), we introduced a rapidity cut-off $y_1$, required in the subtraction mechanism of the overlapping between soft and collinear momentum regions; $y_1$ acts as a lower bound for the rapidity of the particles described by the TMD, which are supposed to be collinear, hence very fast moving along the reference direction (plus direction) of the jet. Therefore, despite it is a full-fledged arbitrary cut-off, its value has to be chosen “large enough” to preserve the physical meaning of the TMDs. In fact, in a physical observable $y_1$ should be set in the limit in which the factorization procedure holds, i.e. $y_1 \rightarrow \infty$. However, TMDs are not physical observables and they depend on $\mu$ as well as $y_1$, hence they can be considered at a fixed, and finite, value of $y_1$.

The transformation rule for a shift in the rapidity cut-off can be easily obtained from the solution of the evolution equations in Eq. (31). If we shift $y_1$ to ${\widehat{y}}_1 = y_1 - \theta$, where $\theta$ is some real number, then (neglecting the dependence on all variables except $\zeta$) the full TMD transforms as:

$$\begin{aligned}&{\widetilde{C}}(\zeta ) \mapsto {\widetilde{C}}({\widehat{\zeta }}) ={\widetilde{C}}(\zeta ) \, \hbox {exp} \left[ \frac{1}{2} \theta \, {\widetilde{K}}\right] , \end{aligned}$$

(39)

Fig. 3

Pictorial representation of a TMD Fragmentation Function, in which the separation between perturbative and non-perturbative regime is explicitly shown, corresponding to two different values of the rapidity cut-off. In a the rapidity cut-off of the TMD FF is set to a generic value $y_1$. In b the rapidity cut off has been shifted to ${\widehat{y}}_1 > y_1$. The two TMDs represent different physical configurations, as the range spanned by the rapidities of the particles belonging to those TMDs are different. This transformation alters the fragmentation mechanism; in fact at extremely large values of the rapidity cut off, one can reach a quasi 1-dimensional configuration

Therefore, the full effect of this transformation is a dilation factor which depends on the soft kernel ${\widetilde{K}}(b_T,\,\mu )$ and the shift parameter $\theta$. Notice that the transformed TMD describes a different physical configuration, as the rapidities of the collinear particles have been shrinked to a narrower range. In particular, as ${\widehat{y}}_1$ approaches the factorization limit of infinite rapidity, the particles belonging to that TMD become more and more tightly aligned along the reference direction of the collinear group. As a consequence, in the limit of infinite rapidity cut-off, the collinear particle motion is basically 1-dimensional and the 3D picture of the hadron structure is altered (see Fig. 3). Since the non-perturbative information about the 3D structure of the hadrons is encoded into the model $M_C$, we can define a transformation that makes the TMD invariant with respect to the shift of the rapidity cut-off by acting simultaneously on the model. The resulting transformed TMD will describe the same physical configuration of the initial TMD, because the alteration due to the tightened range of rapidity will be totally reabsorbed by the transformed model, that will compensate for the dilation factor $\hbox {exp} \left[ \frac{1}{2} \theta \, {\widetilde{K}}\right]$ in Eq. (39). Then, for any $\theta < 0$ such transformation is defined as:

$$\begin{aligned}&y_1 \mapsto {\mathcal {D}}_{\theta } (y_1) = y_1 - \theta , \end{aligned}$$

(40)

$$\begin{aligned}&M_C(b_T) \mapsto {\mathcal {D}}_{\theta } \left( M_C(b_T)\right) \nonumber \\&\quad = M_C(b_T)\, \hbox {exp} \left[ -\frac{1}{2} \theta \, {\widetilde{K}} (\mu ,\,b_T)\right] , \end{aligned}$$

(41)

where only the dependence on $b_T$ has been shown explicitly in $M_C$. Due to the dilation factor in front of the model, we will refer to the previous transformation ${\mathcal {D}}_{\theta }$ as a rapidity dilation (RD), that makes TMDs invariant with respect to the choice of the rapidity cut-off:

$$\begin{aligned}&{\widetilde{C}}(\zeta ,\,M_C) \mapsto {\mathcal {D}}_{\theta } \left( {\widetilde{C}}(\zeta ,\,M_C) \right) \nonumber \\&\quad = {\widetilde{C}}\left( \zeta e^{2\theta },\,{\mathcal {D}}_{\theta } M_C\right) = {\widetilde{C}}(\zeta ,\,M_C) . \end{aligned}$$

(42)

The transformed model, Eq. (41), acquires the same properties of Eq. (35). In fact, since ${\widetilde{K}}$ goes to zero at small $b_T$, then the dilation factor is 1 for $b_T \sim 0$. Furthermore, since ${\widetilde{K}}$ is basically negative, at large $b_T$ the dilation factor give an additional suppression beside those due to the properties of $g_K$ and $M_C$. Rapidity dilation make TMDs invariant under the choice of the rapidity cut-off $y_1$, which nevertheless has to be considered an arbitrary and large parameter. This is in fact one of the necessary hypothesis at the basis of any factorization formula: all the particles described by a collinear part (and ultimately by a TMD) must have a large and positive rapidity, according to the reference direction of the collinear group. Hence, a rapidity dilation performed with a very large and positive $\theta$ would contradict the initial hypothesis on the validity of factorization itself. The correct way to interpret this transformation is to apply it only after the factorization formula has been derived, in the limit of $y_1 \rightarrow \infty$. Roughly speaking, the model associated with a certain choice describes how collinear particles with rapidity in the range⁶$y_1 \le y < \infty$ behave in the non-perturbative regime. Then, rapidity dilations simply balance the perturbative and non-perturbative information encoded in the TMDs according to the choice of rapidity cut-off, in order to keep their combination invariant. Furthermore, rapidity dilations offer a new point of view which helps in the comparison between different physical configurations. Let’s consider, for example, those depicted in Fig. 3, where panel (a) and panel and (b) represent TMDs described by $C(\zeta ,\,M)$ and $C({\mathcal {D}}_\theta \zeta ,\,M)$, respectively, with ${\mathcal {D}}_\theta \zeta < \zeta$. It is interesting to point out that this interpretation is totally equivalent to considering the two TMDs evaluated within the same range of rapidity but associated to two different non-perturbative models, i.e. interpreting the TMD depicted in panel (a) as described by $C({\mathcal {D}}_\theta \zeta ,\,{\mathcal {D}}_\theta M)$ and the TMD in panel (b) as described by $C({\mathcal {D}}_\theta \zeta ,\,M)$. Notice that rapidity dilations define a group under the multiplication laws:

$$\begin{aligned} {\mathcal {D}}_{\theta _2} \circ {\mathcal {D}}_{\theta _1} = {\mathcal {D}}_{\theta _1+ \theta _2} , \quad {\mathcal {D}}_{\theta } \circ {\mathcal {D}}_{-\theta } = \mathrm {id} \,. \end{aligned}$$

(43)

Rapidity dilations are closely reminiscent of a 1-parameter gauge transformations for the “fields” $y_1$ and $M_C$, that make the TMD invariant. Here the TMD plays the role of the “Lagrangian”. In this sense, rapidity dilations might be considered a symmetry for the TMDs. There is an interesting analogy between the action of rapidity dilations on the rapidity cut-off, $y_1$, and the action of the Renormalization Group (RG) on the energy scale $\mu$. Here, an arbitrary $\mu$ allows to regularize the UV divergences, but it introduces some arbitrariness in the theory, as $\mu$ can be set to any value. Some quantities are independent of the choice of this scale, like cross sections, where the RG-variation of the fields is exactly compensated by the RG-variation of the external on-shell particles (LSZ mechanism). Similarly, rapidity dilations (RD) allow to control the arbitrariness in the choice of the rapidity cut-off, $y_1$, that regularizes the rapidity divergences. In particular, TMDs defined along the plus direction are RD-invariant, as the transformation of the model $M_C$, Eq. (41), exactly compensates for the rapidity shift, Eq. (40). However, there are quantities that are not RD-invariant. As we will show in Sect. 3.2.1, the TMD defined along the minus direction ${\widetilde{C}}_-$ or the 2-h soft factor $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ are examples of such quantities. On the other hand, combinations as ${\widetilde{C}}_+ \, {\widetilde{C}}_- \, \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ are RD-invariant, because the dilation factor in ${\widetilde{C}}_-$ exactly compensates the variation in $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$.

When a TMD appears in a cross section, its non-perturbative content, i.e. the last line of Eq. (38), has to be extracted from experimental data and the result will depend on the choice of the rapidity cut-off, as represented in Fig. 3. Rapidity dilations ensure that the physical information encoded in the TMDs stays unaltered if the rapidity cut-off is moved toward the limit of infinite rapidity. If ${\widetilde{C}}^{\mathrm{NP}}$ denotes the full non-perturbative content of the TMD, then its transformation rule under rapidity dilation is given by:

$$\begin{aligned}&{\mathcal {D}}_{\theta }\left( {\widetilde{C}}^{\mathrm{NP}}_{j,\,H} \right) (\zeta ,\,M,\,g_K) \nonumber \\&\quad = {\widetilde{C}}^{\mathrm{NP}}_{j,\,H} ({\mathcal {D}}_{\theta }\zeta ,\,{\mathcal {D}}_{\theta }M,\,g_K) \nonumber \\&\quad = \left( {\mathcal {D}}_{\theta }M\right) _{j,\,H}(b_T) \hbox { exp} \Big \{ -\frac{1}{4} \, g_K(b_T) \, \log {\frac{{\mathcal {D}}_{\theta }\zeta }{{\overline{\zeta }}_0}} \Big \}\nonumber \\&\quad = \left( M\right) _{j,\,H}(b_T) \hbox { exp} \Big \{ -\frac{1}{4} \, g_K(b_T) \, \log {\frac{\zeta }{{\overline{\zeta }}_0}} \Big \} \,\nonumber \\&\quad \quad \times \hbox { exp} \Big \{ -\frac{1}{2} \,\theta \, {\widetilde{K}}(b^\star _T,\,\mu ) \Big \} \nonumber \\&\quad = {\widetilde{C}}^{\mathrm{NP}}_{j,\,H}(\zeta ,\,M,\,g_K) \, \hbox { exp} \Big \{ -\frac{1}{2} \,\theta \, {\widetilde{K}}(b^\star _T,\,\mu ) \Big \}, \end{aligned}$$

(44)

where the second step is given by Eqs. (15), (40) and (41). In this case, the dilation factor is fully computable in perturbative QCD. Notice that the function $g_K$ is not affected by the rapidity dilation, as it should. In fact $g_K$ is also involved in the definition of the soft factor $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ (Eq. (18)), which must not depend on the extraction of the TMD.

Due to rapidity dilations, the choice of the model depends on the choice of the rapidity cut-off. Therefore, in general two independent extractions of TMDs, that use different values of $\zeta$, will feature different models. However, rapidity dilations allow to relate these independent extractions of TMDs. The main difficulty here is that theory is devised in the $b_T$-space, while measurements are performed in the transverse momentum space.

As a practical example, let’s suppose we want to compare the TMD extractions of two independent research Groups, A and B, that have analyzed the same sample of data. They will provide two TMD functions in transverse momentum space $C^{(A)}$ and $C^{(B)}$. Since they obtained their result fitting the same experimental data, the two functions have to be compatible within the overlapping of the respective uncertainty bands, built by considering all source of errors (collinear PDFs/FFs uncertainties, experimental errors, fitting uncertainties, etc ...). However a meaningful comparison can only be made for small values of transverse momentum, because TMDs turn non-physical at large $k_T$ (see B.2). Both results can be written as the Fourier transform of their $b_T$ counterparts. Schematically:

$$\begin{aligned}&C^{(A)} (k_T,\,\zeta ,\,M^{(A)}) \nonumber \\&\quad = \int \, \frac{d^2 {{\varvec{b}}}_T}{(2\pi )^2} \, e^{i {{\varvec{k}}}_T \cdot {{\varvec{b}}}_T} \, {\widetilde{C}}^{\mathrm{P}}_{\zeta } (b^\star _T) \, {\widetilde{C}}^{\mathrm{NP}}_{\zeta ,\,M^{(A)}} (b_T); \end{aligned}$$

(45)

$$\begin{aligned}&C^{(B)} (k_T,\,\zeta ',\,M^{(B)}) \nonumber \\&\quad = \int \, \frac{d^2 {{\varvec{b}}}_T}{(2\pi )^2} \, e^{i {{\varvec{k}}}_T \cdot {{\varvec{b}}}_T} \, {\widetilde{C}}^{\mathrm{P}}_{\zeta '} (b^\star _T) \, {\widetilde{C}}^{\mathrm{NP}}_{\zeta ',\,M^{(B)}} (b_T). \end{aligned}$$

(46)

where only the dependence on the rapidity cut-off and on the model are shown explicitly and ${\widetilde{C}}^{\mathrm{P}}$ denotes the perturbative content of the TMD. In principle, also the choice of $g_K$ may be different for the two extractions. However, Group A and B have to agree also in the estimate of the ${\mathbb {S}}_{\mathrm{{2}-h}}$, and this gives further constraints. Hence, even if $g_K^{(A)}$ and $g_K^{(B)}$ have different functional forms, they should share more or less the same shape. For simplicity, in the following we will set $g_K^{(A)} \sim g_K^{(B)}$.

The two rapidity cut-off are different but, supposing $\zeta ' < \zeta$, a certain real number $\theta < 0$ must exists such that $\zeta ' = \zeta e^{2\theta }$. Hence, Group A can perform a rapidity dilation in order to comply with the choice of Group B. By using Eq. (44), they can write:

$$\begin{aligned}&C^{(A)} (k_T,\,\zeta ,\,M^{(A)}) \nonumber \\&\quad = \int \frac{d^2 {{\varvec{b}}}_T}{(2\pi )^2} \, e^{i {{\varvec{k}}}_T \cdot {{\varvec{b}}}_T} \, {\widetilde{C}}^{\mathrm{P}}_{\zeta '} (b^\star _T) \, {\widetilde{C}}^{\mathrm{NP}}_{\zeta ',\,{\mathcal {D}}_\theta \, M^{(A)}} (b_T) \nonumber \\&\quad = \int \frac{d^2 {{\varvec{b}}}_T}{(2\pi )^2} \, e^{i {{\varvec{k}}}_T \cdot {{\varvec{b}}}_T} \, {\widetilde{C}}^{\mathrm{P}}_{\zeta '} (b^\star _T) \, {\widetilde{C}}^{\mathrm{NP}}_{\zeta ,\,M^{(A)}} (b_T)\, \nonumber \\&\qquad \times \hbox {exp} \Big \{ -\frac{1}{2} \,\theta \, {\widetilde{K}}(b^\star _T) \Big \}. \end{aligned}$$

(47)

In this way, the two estimates are written with the same perturbative part in $b_T$-space. If TMDs were valid throughout the whole spectrum of $k_T$s, then the comparison between Eqs. (46) and  (47) would lead to ${\mathcal {D}}_\theta \, M^{(A)} \sim M^{(B)}$. However, the Fourier transforms in Eqs. (45) and (46) have to be compatible at small $k_T$, but there is no constraint for larger values. Furthermore, we known that the perturbative content is constant for $b_T$ larger than a certain $b_{\mathrm{SAT}} \ge b_{\mathrm{MAX}}$. At the same time, the non-perturbative content should be of order 1 at small/moderate $b_T$, i.e. up to $b_{\mathrm{SAT}}$, in order not to interfere too drastically on the perturbative information. Therefore, we can roughly split the Fourier transform in two parts as:

$$\begin{aligned} C (k_T)&\sim \int ^{b_{\mathrm{SAT}}} \frac{d^2 {{\varvec{b}}}_T}{(2\pi )^2} \, e^{i {{\varvec{k}}}_T \cdot {{\varvec{b}}}_T} \,\nonumber \\&\quad \times \left( {\widetilde{C}}^{\mathrm{P}}_{\zeta } (b^\star _T) - {\widetilde{C}}^{\mathrm{P}}_{\zeta } (b_{\mathrm{SAT}}) \right) \, {\widetilde{C}}^{\mathrm{NP}}_{\zeta ,\,M} (b_T)\nonumber \\&\quad + {\widetilde{C}}^{\mathrm{P}}_{\zeta } (b_{\mathrm{SAT}}) \, \int \frac{d^2 {{\varvec{b}}}_T}{(2\pi )^2} \, e^{i {{\varvec{k}}}_T \cdot {{\varvec{b}}}_T} \, {\widetilde{C}}^{\mathrm{NP}}_{\zeta ,\,M} (b_T). \end{aligned}$$

(48)

The first part is integrated only up to the saturation value $b_{\mathrm{SAT}}$. It is clearly dominated by perturbative information, since ${\widetilde{C}}^{\mathrm{NP}}$ is not drastically different from 1 in that range, for any choice of $\zeta$ and M. As a consequence, this part is almost the same for both Eqs. (47) and  (46). On the other hand, the second part is simply proportional to the Fourier Transform of the non perturbative content of the TMD. Different choices of $\zeta$ and M can give integrands of the same order up to $b_{\mathrm{SAT}}$ but they can differ on how rapidly they go to zero as $b_T$ goes to infinity: at large $b_T$ they could be very small and at the same time differ for many orders of magnitude. This difference is not evident at small values of $k_T$, since the area under the curve in $b_T$ space, after the saturation value, can be neglected in any case. However, the differences may be consistent at large $k_T$. This is not a problem, since TMDs lose their physical meaning in this region. Hence, Group A can compare its result with that of Group B by Fourier transforming its non-perturbative part after it has been rapidity-dilated. Then, the following relation should hold within uncertainties:

$$\begin{aligned} R^{\mathrm{NP}}(k_T)&= \frac{\displaystyle {\mathcal {FT}} \Big [ {\widetilde{C}}^{\mathrm{NP}}_{\zeta ',\,M^{(B)}} (b_T)\Big ]}{\displaystyle {\mathcal {FT}} \Big [ {\widetilde{C}}^{\mathrm{NP}}_{\zeta ,\,M^{(A)}} (b_T)\, \hbox { exp} \Big \{ -\frac{1}{2} \,\theta \, {\widetilde{K}}(b^\star _T) \Big \} \Big ]}\nonumber \\&\sim 1 \quad \hbox { at small } k_T. \end{aligned}$$

(49)

3.2.1 Rapidity dilation and z-axis reflection

The behaviour under z-axis reflection, which simply exchanges the plus and minus directions, is particularly important for widely studied processes, like SIDIS, Drell-Yan and $e^+e^-\rightarrow H_A\,H_B \,X$, where two TMDs associated to opposite directions are multiplied together. If $R_z$ is the Lorentz transformation that reverses the z-axis, then the rapidity of the reference hadron swaps its sign under the action of $R_z$. On the other hand, the rapidity cut-off is not the rapidity of any real particle. It is just an ad hoc number and hence it is trivially invariant under the action of $R_z$. However, the particles belonging to the collinear group associated to the TMD in the minus direction should have a very large negative rapidity according to the limit $y_1 \rightarrow +\infty$. Therefore, a proper rapidity cut-off would be $y_2 = - y_1$, as if $y_1$ had changed its sign. Summarizing:

$$\begin{aligned} {\left\{ \begin{array}{ll} y_P &{}\mapsto R_z\left( y_P \right) = - y_P; \\ y_1 &{}\mapsto R_z\left( y_1 \right) = y_1 {\mathop {=}\limits ^{def}} - y_2 . \end{array}\right. } \end{aligned}$$

(50)

As a consequence, the variable $\zeta$ for a TMD in the minus direction is obtained by simply replacing $\zeta _{+} \propto \hbox {exp} \left( y_P - y_1\right)$ with $\zeta _{-} \propto \hbox {exp} \left( y_2 - y_P\right)$ and the full TMD transforms as:

$$\begin{aligned} {\widetilde{C}}_{+} (\zeta _{+}) \mapsto R_z\left( {\widetilde{C}}_{+} (\zeta _{+}) \right) = {\widetilde{C}}_{-} (\zeta _{-}) , \end{aligned}$$

(51)

where only the dependence on the rapidity cut-off has been made explicit.

There is a non trivial interplay between z-axis reflection and rapidity dilations, since the two transformations do not commute. In fact, if the rapidity cut-off $y_1$ of $C_+$ is shifted, then the rapidity cut-off $y_2$ of $C_-$ is shifted as well, but with the sign reversed. This can easily be seen by a direct computation, with the help of Eqs. (40) and (50):

$$\begin{aligned} {\mathcal {D}}_{\theta } \left( y_2\right) = {\mathcal {D}}_{\theta } \left( -y_1\right) = -y_1 + \theta = y_2 + \theta . \end{aligned}$$

(52)

Therefore, according to Eq. (41), the model of $C_-$ transforms as:

$$\begin{aligned} {\mathcal {D}}_{\theta } \left( M_{C_-}(b_T)\right) = M_{C_-}(b_T)\, \hbox {exp} \left[ \frac{1}{2} \theta \, {\widetilde{K}}\right] . \end{aligned}$$

(53)

However, in the z-reversed TMD, $C_-$, the rapidity cut-off appears with the opposite sign with respect to $C_+$. Hence, there is no more compensation between the rapidity shift and the transformed model, and $C_-$ is not invariant under rapidity dilations. This can be summarized by saying that the two transformations do not commute:

$$\begin{aligned} {\left\{ \begin{array}{ll} &{}R_z\left( {\mathcal {D}}_{\theta } \left( {\widetilde{C}}_{+} (\zeta _{+}) \right) \right) = R_z\left( {\widetilde{C}}_{+} (\zeta _{+}) \right) = {\widetilde{C}}_{-} (\zeta _{-}) ;\\ &{}{\mathcal {D}}_{\theta }\left( R_z \left( {\widetilde{C}}_{+} (\zeta _{+}) \right) \right) = {\mathcal {D}}_{\theta }\left( {\widetilde{C}}_{-} (\zeta _{-}) \right) = {\widetilde{C}}_{-} (\zeta _{-})\, \hbox {exp} \left[ \theta \, {\widetilde{K}}\right] . \end{array}\right. } \end{aligned}$$

(54)

Finally, lets consider the behavior of the (asymptotic) 2-h soft factor, defined in Eq. (18), under rapidity dilations. Since the soft model is not affected by the transformation and $y_1 \rightarrow y_1 - \theta$, while $y_2 \rightarrow y_2 + \theta$, it transforms as:

$$\begin{aligned}&{\mathcal {D}}_{\theta } \left( \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}(b_T; \, \mu ,\,y_1 - y_2) \right) = \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}(b_T; \, \mu ,\,y_1 - y_2) \, \nonumber \\&\quad \hbox {exp} \left[ -\theta \, {\widetilde{K}}\right] . \end{aligned}$$

(55)

Therefore, by exploiting Eqs. (42), (54), and (55), the combination ${\widetilde{C}}_+ \, {\widetilde{C}}_- \, \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ cross section) is invariant under rapidity dilations.

4 Universality and process classification

Process-independent quantities play the most important role in factorized cross sections. They are universal, which means that once they have been estimated they can be used in any cross section, regardless of the specific process. This is particularly useful for those quantities that carry non-perturbative information. Since they cannot be computed analytically, they have to be extracted from experimental data. However, if they are universal, any process that allows for their presence in the cross section can be exploited, and we can prefer those with a richer amount of data. A lack of universality would undermine the predictive power of QCD itself. In fact, if the non-perturbative quantities had to be extracted again for each individual process, the phenomenological analysis of a hadronic cross sections would be reduced to a mere fit of experimental data.

In general, a factorized cross section is a convolution of three different objects: the hard part, the collinear factors and the soft factor (see Sect. 1).

The hard part is completely process-dependent. However, it can be computed in perturbative QCD and its lack of universality does not affect the predictive power of the theory.

Collinear parts and the TMDs, as defined in Sect. 3 by the factorization definition, Eq. (25), depend only on their internal variables and hence are completely blind to the kinematics of the process. Therefore, they can be really considered universal quantities.

On the other hand, the soft factor, defined in Sect. 2, is not completely process-independent. In fact, it depends on the number N of the collinear factors involved in the factorized cross section, each related to its reference parton of type j and to its reference hadron H. Therefore, they are not insensitive to the kinematics of the process in which they appear, because they depend both on the number of the Wilson lines replacing the collinear parts and also on their color representation, which is fixed by the parton type j and differs from quark and gluons. However, at fixed N and for reference partons of the same kind, soft factors are actually the same object, modulo crossing symmetry. As an example, Drell-Yan scattering with two quark-initiated collinear factors in the initial state, $e^+e^-\rightarrow H_A \, H_B\, X$, with two quark-initiated collinear factors in the final state, and also SIDIS, with one quark-initiated collinear factor in the initial and one in the final state, share the same soft factor ${\mathbb {S}}_{\mathrm{{2}-h}}$ modulo the crossing symmetry that relates the three processes. Notice that in this case there are only two collinear factors and charge conservation allows only two quarks as reference partons.

Since processes with a different number N of collinear factors have a different soft factor in their factorized cross section, it is possible to classify them according to this number. This coincides with the number of reference hadrons participating to the hadronic process. The classes derived with this criterion will be called hadron classes. Formally, a process belongs to the N-h class if it globally involves N collinear parts, which can appear in the initial and/or in the final state, in all possible combinations and for all the allowed kind of reference partons. Therefore, ${\mathbb {S}}_{\mathrm{{N}-h}}$ can be considered universal only within the N-h class, modulo crossing symmetry and the possible color representations of its Wilson lines. This is a weaker kind of universality, that holds only for a limited number of processes. For instance, processes involving one collinear group belong to the ${1}$-hadron class: deep inelastic scattering (DIS), corresponding to one reference hadron in the initial state; $e^+e^-\rightarrow H\,X$ correspoding to one reference hadron in the final state. Processes involving two collinear groups belong to the ${2}$-hadron class: here we have Drell–Yan like scattering, $e^+e^-\rightarrow H_A \, H_B\, X$, and SIDIS.

The classification above has nothing to do with the nature of the factorization adopted (collinear or TMD): it depends only on the specific kinematics of the processes, that can be different case by case. However, it is possible to identify common properties within each hadron-class that allows to determine, a priori, which factorization scheme should be used. Consider for example the ${1}$-hadron class case. In both DIS and $e^+e^-\rightarrow H\,X$ there is at least one hard real emission, since there is always a fermion leg crossing the final state cut (see Fig. 4). The collinear factor associated to the real emission is totally crossed by the final state cut and hence it does not have any reference hadron. Therefore, it can be considered far off-shell and part of the hard factor. All the information about soft transverse momentum is washed away and the collinear factorization scheme has to be applied. The factorized cross section for ${1}$-hadron class processes is then written as a convolution of the collinear part associated to the reference hadron with an hard factor that, once considered together with the hard real emissions, can be interpreted as a partonic cross section, i.e. the partonic counterpart of the process.

Fig. 4

Pictorial representation of DIS (a) and of $e^+e^-\rightarrow H \, X$ (b). In both cases, there is at least one hard real emission, which produces a collinear factor completely crossed by the final state cut. Consequently it can be reabsorbed in the hard factor of the cross section

In the ${2}$-hadron class, instead, the choice of factorization scheme is non-trivial and depends on the specific kinematics of the process. It is dictated by the size of one parameter, namely the ratio between the modulus of the weak boson transverse momentum $q_T$ and the typical energy scale of the process Q (see Ref. [2]). When ${q_T}/{Q} \ll 1$, TMD factorization has to be applied, while if ${q_T}/{Q} \gg 1$, collinear factorization will be appropriate.⁷ The cross sections predicted in these two kinematical ranges, computed within two different approximations, do not automatically match; in fact, the intermediate region, where $q_T\sim Q$, is usually called “matching region”. Several studies have been devoted to the implementation of different algorithms to map these kinematics regions and to match the collinear cross sections to the TMD cross section (a problem known as “matching”), see for example Refs. [12,13,14,15,16].

According to the previous considerations, one can build a hierarchy based on universality. The lowest level is occupied by quantities, like the hard part, that are completely process dependent but usually fully computable in perturbation theory. At the top of the hierarchy we find quantities, like the collinear factors, that are absolutely process independent: they carry non-perturbative information but their universality properties guarantee that they can be extracted from one particular process and then used in any other. In the middle there are quantities which are only universal within their own N-hadron class, like the soft factors. As they carry non-perturbative information, they cannot be computed perturbatively. They too have to be extracted from experimental data, but they can only be used, class by class, for the processes involving the same number of collinear groups and for the same kind of reference partons.

In this sense, it is very important to provide a working scheme where objects with different degrees of universality are neatly separated, in such a way to maximize their perturbative content and their universal parts, while reducing the class-dependent factors to the minimum. For the latter, special experimental efforts will be required in order to gather a large number of high quality data corresponding to several different processes, which will then be analyzed simultaneously in a completely consistent framework. The latest analyses of the BELLE Collaboration and the current plans towards the realization of a new Electron Ion Collider (EIC) are indeed moving towards this direction [1, 17,18,19].

The classification introduced above has to be intended as a criterion to classify processes on the basis of their factorized cross section properties, and of their corresponding soft factor. Therefore, the number N that labels the classes is not the number of all hadrons involved in the process, in general much greater than the number of collinear factors. The difference is more evident when we consider the final state of a scattering process. In general, experimentalists detect a huge number of hadrons, grouped in jets. The number of jets does not correspond to the number of collinear factors, which instead is the number of reference hadrons, i.e. the number of jets in which the hadron is detected in order to study the jet’s fragmentation properties. The actual topology of the event (e.g. the number of jets) is described by event-shape variables, like thrust. An example will be presented in Sect. 6.

5 The 2-hadron class

We will now focus on the ${2}$-hadron class of processes. As mentioned above, in this class the choice of factorization scheme depends on a single parameter, the ratio ${q_T}/{Q}$ (see Sect. 1.1). The ${2}$-hadron class plays a crucial role, as its soft factor ${\mathbb {S}}_{\mathrm{{2}-h}}$ is exactly the same object that appears at denominator in the subtracted collinear factor ${\mathbb {C}}$, Eq. (22) and, consequently, in the general definition of the TMD, Eqs. (23) and (24).

5.1 2-h class cross section

In Sect. 2 we have provided a useful formalism to decompose the 2-h class soft factor and the collinear part ${\mathbb {C}}$ in a fully perturbative (computable) part and a strictly non-perturbative term, which can be modeled through the functions $g_K(b_T)$, $M_S(b_T)$ and $M_C(b_T)$, as shown in Eqs. (18) and (31). At this stage we have achieved all the necessary tools to be able to write an explicit expression for the ${2}$-hadron class cross section. Its generic structure is analogous to that given in Eq. (2) for $e^+e^-\rightarrow H_A\, H_B \, X$, with $H_A$ and $H_B$ in an almost back-to-back configuration:

$$\begin{aligned} d \sigma _{\mathrm{2-h}}&\sim H \, \times \, {\mathcal {FT}} \Big [ \, {\widetilde{C}}_+ \, \times \, {\widetilde{C}}_- \, \times \, \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} \, \Big ]\nonumber \\&\sim H \, \times \, {\mathcal {FT}} \Big [ \, \dfrac{{\widetilde{C}}_+^{\mathrm{unsub}}}{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}} \, \times \, \dfrac{{\widetilde{C}}_-^{\mathrm{unsub}}}{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}} \, \times \, \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} \, \Big ]\,, \end{aligned}$$

(56)

where $C_+$, $C_-$ refers to TMDs defined along the plus and the minus direction, respectively. The soft factor $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ appearing in Eq. (56) is the same object that appears as subtraction factor in the factorization definition of the TMDs. Reorganizing the three $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ factors and reabsorbing them in the TMD, leads to a different definition of TMDs (see e.g. [2, 20]):

$$\begin{aligned}&{\widetilde{C}}_+^{\; \mathrm{sqrt}}(\xi _+,\, {{\varvec{b}}}_{T}; \, \mu , \, y_{P_1} - y_n) \nonumber \\&\quad = \lim \limits _{\begin{array}{c} y_{u_1} \rightarrow +\infty \\ y_{u_2} \rightarrow -\infty \end{array}} \, {\widetilde{C}}_+^{\mathrm{unsub}} (\xi _+, \,{{\varvec{b}}}_{T}; \, \mu , \, y_{P_1} - y_{u_2}) \,\nonumber \\&\qquad \times \sqrt{\dfrac{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} ({{\varvec{b}}}_{T}; \, \mu , \,y_{u_1} - y_n)}{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} ({{\varvec{b}}}_{T}; \, \mu , \,y_{u_1} - y_{u_2}) \, \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} ({{\varvec{b}}}_{T}; \, \mu , \,y_n - y_{u_2})}} \end{aligned}$$

(57)

$$\begin{aligned}&{\widetilde{C}}_-^{\; \mathrm{sqrt}}(\xi _-,\, {{\varvec{b}}}_{T}; \, \mu , \, y_n-y_{P_2}) \nonumber \\&\quad = \lim \limits _{\begin{array}{c} y_{u_1} \rightarrow +\infty \\ y_{u_2} \rightarrow -\infty \end{array}} \, {\widetilde{C}}_-^{\mathrm{unsub}} (\xi _-, \,{{\varvec{b}}}_{T}; \, \mu , \, y_{u_1} - y_{P_2}) \,\nonumber \\&\qquad \times \sqrt{\dfrac{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} ({{\varvec{b}}}_{T}; \, \mu , \,y_n - y_{u_2})}{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} ({{\varvec{b}}}_{T}; \, \mu , \,y_{u_1} - y_{u_2}) \, \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} ({{\varvec{b}}}_{T}; \, \mu , \,y_{u_1} - y_n)}} . \end{aligned}$$

(58)

This definition of TMDs is often referred to as the square root definition.

There are many advantages to it. First of all, a single rapidity cut-off $y_n$ is sufficient to regularize all rapidity divergences, the perturbative computations are much easier and the evolution equation are unified and symmetrized, see Ref. [2]. Moreover, as mentioned above, the square root definition allows to solve the soft factor problem in the ${2}$-hadron class. In fact, according to this definition, the cross section assumes a “Parton-Model”-like structure, where all soft gluons are reabsorbed in the TMD definition, very convenient for phenomenological applications:

$$\begin{aligned} d \sigma _{\mathrm{2-h}}&\sim H \, \times \, {\mathcal {FT}} \Big [ \, {\widetilde{C}}_+ \, \times \, {\widetilde{C}}_- \, \times \, \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} \, \Big ]\nonumber \\&\sim H \, \times \, {\mathcal {FT}} \Big [ \, {\widetilde{C}}_+^{\; \mathrm{sqrt}} \times {\widetilde{C}}_-^{\; \mathrm{sqrt}} \, \Big ]\,, \end{aligned}$$

(59)

As an example, the unpolarized cross section for $e^+e^-\rightarrow H_A\, H_B \, X$ for almost back-to-back spinless hadrons, Eq. (2), becomes:

$$\begin{aligned}&W^{\mu \,\nu }(Q,\,p_A,\,p_B) \nonumber \\&\quad = \frac{8 \pi ^3 z_A z_B}{Q^2} \sum _f H^{\mu \,\nu }_{f,{\overline{f}}}(Q) \nonumber \\&\qquad \times \int d^2 {{\varvec{b}}}_T \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}({{\varvec{b}}}_T) {\widetilde{D}}_{1,H_A / f} (z_A, {{\varvec{b}}}_T) {\widetilde{D}}_{1 H_B / {\overline{f}}}(z_B, {{\varvec{b}}}_T)\nonumber \\&\quad = \frac{8 \pi ^3 z_A z_B}{Q^2} \sum _f H^{\mu \,\nu }_{f,\,{\overline{f}}}(Q) \nonumber \\&\qquad \times \int d^2 {{\varvec{b}}}_T {\widetilde{D}}_{1,H_A / f}^{\; \mathrm{sqrt}} (z_A, {{\varvec{b}}}_T) {\widetilde{D}}_{1, H_B / {\overline{f}}}^{\; \mathrm{sqrt}} (z_B, {{\varvec{b}}}_T) \,. \end{aligned}$$

(60)

Despite its numerous advantages, the square root definition lowers the degree of universality of the TMD, as it relates it to the 2-h soft factor which, by definition, is only universal within its corresponding 2-h class. In other words, the square root definition is optimal for the ${2}$-hadron class, as it beautifully simplifies the 2-h cross section making it suitable for phenomenological applications; its drawback, however, is that it ceases to be valid outside the ${2}$-hadron class. On the other hand, abandoning the square root definition of the TMDs in favor of the factorization definition, Eq. (23), will force us to face the soft factor problem and take a new (and potentially very hard) challenge: reformulating the way we do phenomenology, in terms of newly defined fundamental objects, where the soft factors are modeled explicitly rather than absorbed in the definition of the TMDs.

We will attempt such a strategy, adopting the factorization definition of the TMD, Eq. (25), and relying on the results of Sects. 2 and 3 for the decomposition of the collinear and soft factors in terms of their perturbative and non-perturbative parts.

Using the solution of the evolution equations for the TMDs, Eq. (38), and the soft factor, Eq. (18), it is possible to write the ${2}$-hadron class cross section in terms of perturbative and non-perturbative functions. Apart from the hard factor and a Fourier transform, the relevant structure is given by:

$$\begin{aligned}&{\widetilde{C}}_+(\xi _+,\, {{\varvec{b}}}_{T}; \, \mu , \, \zeta _1) \; {\widetilde{C}}_-(\xi _-,\, {{\varvec{b}}}_{T}; \, \mu , \,\zeta _2) \; \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}\nonumber \\&\qquad \times (b_T; \, \mu ,\,y_1 - y_2) \nonumber \\&\quad = {\widetilde{C}}_+(\xi _+,\, {{\varvec{b}}}_{T}^\star ; \, \mu _0, \,\mu _0^2) \; {\widetilde{C}}_-(\xi _-,\, {{\varvec{b}}}_{T}^\star ; \, \mu _0, \,\mu _0^2) \nonumber \\&\quad \times \hbox { exp} \Big \{ \frac{1}{4} \, {\widetilde{K}}(b_T^\star ;\,\mu _0) \, \log {\frac{\zeta _1 \zeta _2}{\mu _0^4}}\nonumber \\&\qquad + \int _{\mu _0}^{\mu } \frac{d \mu '}{\mu '} \, \left[ 2 \gamma _C(1) - \frac{1}{4} \, \gamma _K(\mu ') \, \log {\frac{\zeta _1 \zeta _2}{\mu '^4}} \right] \Big \} \nonumber \\&\quad \times M_{C_+}(\xi _+,\,b_T) \, M_{C_-}(\xi _-,\,b_T)\nonumber \\&\qquad \hbox {exp} \Big \{ -\frac{1}{4} \, g_K(b_T) \, \log {\frac{\zeta _1 \zeta _2}{\overline{\zeta _1}_0 \overline{\zeta _2}_0}} \Big \} \nonumber \\&\quad \times \hbox { exp} \Big \{ \frac{y_1 - y_2}{2} \left[ {\widetilde{K}}(b_T^\star ;\,\mu _0) - \int _{\mu _0}^{\mu }\, \frac{d \mu '}{\mu '}\, \gamma _K(\mu ) \right] \Big \} \nonumber \\&\quad \times M_S(b_T) \hbox { exp} \Big \{ -\frac{y_1 - y_2}{2} g_K(b_T) \Big \} \,, \end{aligned}$$

(61)

where the reference values of the scales can be set to standard choices, Eqs. (32), (33), (34) and the errors due to the evolution equations are neglected, since they are suppressed by ${\mathcal {O}}\Big ( e^{-(y_1 - y_2)} \Big )$. From Sect. 3.2.1, the product of the two rapidity cut-off gives $\zeta _1 \, \zeta _2 \sim Q^4 e^{-2(y_1-y_2)}$, hence the second and the third lines in Eq. (61) generate contributions that exactly cancel the fourth line and the exponential of the fifth line, respectively. Therefore, we simply have:

$$\begin{aligned}&{\widetilde{C}}_+(\xi _+,\, {{\varvec{b}}}_{T}; \, \mu , \, \zeta _1) \; {\widetilde{C}}_-(\xi _-,\, {{\varvec{b}}}_{T}; \, \mu , \,\zeta _2) \; \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}\nonumber \\&\qquad \times (b_T; \, \mu ,\,y_1 - y_2)\nonumber \\&\quad = {\widetilde{C}}_+(\xi _+,\, {{\varvec{b}}}_{T}^\star ; \, \mu _0, \,\mu _0^2) \;{\widetilde{C}}_-(\xi _-,\, {{\varvec{b}}}_{T}^\star ; \, \mu _0, \,\mu _0^2) \nonumber \\&\qquad \times \hbox { exp} \Big \{ {\widetilde{K}}(b_T^\star ;\,\mu _0) \, \log {\frac{Q}{\mu _0}}\nonumber \\&\qquad + \int _{\mu _0}^{\mu } \frac{d \mu '}{\mu '} \, \left[ 2 \gamma _C(1) - \gamma _K(\mu ') \, \log {\frac{Q}{\mu '}} \right] \Big \} \nonumber \\&\qquad \times M_{C_+}(\xi _+,\,b_T) \, M_{C_-}(\xi _-,\,b_T) \, M_S(b_T) \nonumber \\&\qquad \times \hbox { exp}\Big \{ -g_K(b_T) \, \log {\frac{Q}{\sqrt{\overline{\zeta _1}_0 \overline{\zeta _2}_0}}} \Big \} \,. \end{aligned}$$

(62)

As expected, in the previous equation there is no residual dependence on the rapidity cut-offs $y_1$ and $y_2$, hence we can simply set $\zeta _{1,\,2} = Q^2$. Needless to say, this result is compatible with rapidity dilations since, as shown in Sect. 3.2.1, the combination ${\widetilde{C}}_+ {\widetilde{C}}_- \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ is invariant on the choice of the rapidity cut-off.

The same steps that lead to Eq. (62) could be repeated using the square root definition of the TMDs, Eqs. (57) and (58). The final result will be very similar to Eq. (62), with the exception that the soft model, $M_S$, would not be there. One may therefore be induced to think that adopting the factorization definition instead of the square root definition results in the replacement of the whole 2-h soft factor with its long-distance behavior, given by the sole $M_S$. However, this is not the case. In fact, the 2-h soft factor strongly correlates the two TMDs and combines with them in the structure of the factorization theorem of Eq. (59). In particular, the exponential resulting from the solution of the evolution equations (see Eq. (18)), which encodes the whole dependence on the scale and rapidity cut-offs, combines with its analogous counterparts contained in the two TMDs in such a way that the final cross section does not show any explicit dependence on $y_1$ and $y_2$. The neat effect of applying the factorization definition of the TMD instead of the square root definition can be investigated by direct comparison, as we will show in Sect. 5.2. For the moment being, we may expect that the two TMD definitions will differ in their long-distance behavior. In fact, one might imagine that, in Eq. (62), $M_S$ could be reabsorbed by assigning one of its square roots to $M_{C_+}$ and the other to $M_{C_-}$, giving back the same result that could have been obtained by using the square root definition. This naive intuition will be confirmed in the next Section.

5.2 Factorization definition vs. square root definition

We can now compare the factorization definition with the square root definition of the TMDs. Reference [2] shows that the unsubtracted TMDs ${\widetilde{C}}_i^{\mathrm{unsub}}$ ($i = 1,\,2$), are the same in the two definitions. Hence we can compute their ratio (here we pick the plus direction):

$$\begin{aligned}&\dfrac{{\widetilde{C}}_+^{\; \mathrm{sqrt}}(\xi _+,\, {{\varvec{b}}}_{T}; \, \mu , \, y_{P_1} - y_n)}{{\widetilde{C}}_+ (\xi _+,\, {{\varvec{b}}}_{T}; \, \mu , \, y_P -y_1)}\nonumber \\&\quad = \lim \limits _{\begin{array}{c} y_{u_1} \rightarrow +\infty \\ y_{u_2} \rightarrow -\infty \end{array}} \sqrt{\dfrac{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} (b_{T}; \, \mu , \,y_{u_1} - y_n)}{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} (b_{T}; \, \mu , \,y_{u_1} - y_{u_2}) \, \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} (b_{T}; \, \mu , \,y_n - y_{u_2})}} \,\nonumber \\&\qquad \times \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}} (b_{T}; \, \mu , \,y_1 - y_{u_2}) \nonumber \\&\quad = \sqrt{\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}(b_T; \, \mu _0, \, 0) \, \hbox {exp} \Big ( (y_1 - y_n) \, {\widetilde{K}}(b_T;\,\mu ) \Big )} \nonumber \\&\quad = \sqrt{M_S(b_T)} \times e^{\frac{(y_1 - y_n)}{2} \, {\widetilde{K}}(b_T^\star ;\,\mu )} \, e^{-\frac{(y_1 - y_n)}{2} \, g_K(b_T)} \,, \end{aligned}$$

(63)

where in the second line we used the solution to the evolution equations for the 2-h soft factor, Eq. (13), while in the last step we used Eq. (18) in order to separate the perturbative from the non-perturbative content. Obviously, a perfectly analogous result holds for the TMD relative to the opposite direction, ${\widetilde{C}}_-$. Notice that Eq. (63) is an exact result, since all the corrections to the asymptotic part of the solution to the evolution equations vanish due the limits $y_{u_1} \rightarrow +\infty$ and $y_{u_2} \rightarrow -\infty$, for any finite value of $y_n$, $y_1$.

According to Sect. 3.2 (see e.g. Fig. 3), the physical configuration described by a TMD depends on the cut-off that constrains the rapidity of the particles that move collinearly to the reference direction. Therefore, the comparison between the two different TMD definitions has to be done at the same value of the rapidity cut-off, which implies setting $y_1 = y_n$, to ensure that we compare TMDs that describe the same physics. As a consequence, in Eq. (63) the dependence on the soft kernel ${\widetilde{K}}$ disappears, leaving only a square root of the soft model $M_S(b_T)$. Therefore we have:

$$\begin{aligned}&{\widetilde{C}}^{\; \mathrm{sqrt}}(\xi ,\, {{\varvec{b}}}_{T}; \, \mu , \, y_P - y_n) \nonumber \\&\quad = \sqrt{M_S(b_T)} \times {\widetilde{C}}(\xi ,\, {{\varvec{b}}}_{T}; \, \mu , \, y_P - y_n) \,, \end{aligned}$$

(64)

which clearly holds for both ${\widetilde{C}}_+$ and ${\widetilde{C}}_-$. This is a very important result, as it shows that the choice of TMD definition (square root or factorization definition) only affects the non-perturbative content of the TMDs, while having no impact on the perturbative part. Consequently, $C^{\; \mathrm{sqrt}}$ will differ from C mainly in the small $k_T$ region.

According to Eq. (64), the square root definition is obtained from Eq. (31) by multiplying the TMD defined through the factorization definition by a square root of the soft model. In other words, the contribution of the soft physics just acts on $M_C(\xi ,\,b_T)$:

$$\begin{aligned} M_C^{\; \mathrm{sqrt}} (\xi ,\,b_T)= M_C(\xi ,\,b_T) \times \sqrt{M_S(b_T)}\,. \end{aligned}$$

(65)

To conclude, we can compare the effect of using either one of two different TMD definitions in the cross section. Had we used the square root definition, its net effect in Eq. (62) would have been the replacement:

$$\begin{aligned}&M_{C_+}(\xi _+,\,b_T) \, M_{C_-}(\xi _-,\,b_T) \, M_S(b_T)\nonumber \\&\quad \rightarrow M_{C_+}^{\; \mathrm{sqrt}} (\xi _+,\,b_T) \, M_{C_-}^{\; \mathrm{sqrt}} (\xi _-,\,b_T) \,. \end{aligned}$$

(66)

Clearly the square root definition offers an ideal framework to perform the phenomenological study of the 2-h class of processes: it solves the soft factor problem by reabsorbing it in the TMD definition and allows to extract the model functions $M_{C_{1,\,2}}^{\;{\mathrm{sqrt}}}$ from experimental data. However, this operation makes it impossible to disentangle the non-perturbative soft effects due to $M_S$ which, instead, remains explicit when using the factorization definition for the TMD.

Equation (66) is particularly important from the phenomenological point of view, as it relates the TMDs obtained from data analyses based on the square root definition (which has been very widely used in the last 10 years) to the TMDs extracted using the factorization definition. In this regard, the methodology proposed in this paper allows to profit of the past experience and to benefit of all the results obtained in previous analyses, while extending the scheme to all those processes which could not be considered before, because they belong to a different hadron class. A rather straightforward example of this strategy will be the combined analysis of the BELLE measurements of the polarization of $\varLambda$ hyperons [17] in $e^+e^-\rightarrow \varLambda \pi (K) X$ processes (2-h class), already studied in Ref. [21] within a generalized-parton model approach, and in $e^+e^-\rightarrow \varLambda X$, i.e. in a 1h-class process. This will be presented in a forthcoming paper.

6 Factorization of $e^+e^-\rightarrow H\,X$

In this section we will focus on the $e^+e^-\rightarrow H\,X$ process, which belongs to the ${1}$-hadron class according to the classification of Sect. 4. Here we have only one true collinear part, associated to the reference hadron H, which can be quark- or gluon-initiated; beside, in any case, there is always at least one hard real emission that gives a collinear contribution crossing the final state cut, hence included in the hard factor, which can then be interpreted as a partonic cross section. The soft factor of the process is unity according to the collinear factorization scheme, see Ref. [2].

The thrust T will be included in the derivation of the final result. It is an event-shape variable that describes the topology of the final state, i.e. the number of observed jets. It can take values from 0.5 to 1.0: the lower limit corresponds to a spherical distribution of particles in the final state, while the upper limit indicates an exact two-jet configuration (pencil-like events). Among all jets, only one is related to the collinear part, while the others have to be included in the partonic cross section. Therefore, the value of thrust will determine which Feynman graphs have to be considered in the calculation of the hard part, that will acquire a non-trivial dependence on T.

Similar cross sections have been studied in the framework of Soft Collinear Effective Theories (SCET), within a TMD-like factorization scheme. The relations connecting the SCET to the CSS definition of TMDs have been investigated in Refs. [22, 23], where perfect equivalence has been found with the square root definition, see Eq. 58. In particular, the cross section for $e^+e^-\rightarrow H \hbox {(jet)} \, X$ has been considered in Refs. [24,25,26]. In these papers the dependence on thrust T is not considered; instead the radius R of the jet is introduced as a reference for the transverse momentum of the detected hadron. For the case of $e^+e^-$ event shape angularities, but with no dependence on transverse momentum, one could refer for example to Ref. [27].

In the (thrust dependent) cross section of $e^+e^-\rightarrow H\,X$, the leptonic tensor $L_{\mu \,\nu }$, corresponding to the initial state contribution, is Lorentz contracted with the hadronic tensor $W^{\mu \,\nu }_H$, associated to the final state (see for example Ref. [2]). The cross section is then written as:

$$\begin{aligned} \frac{d \sigma }{{d^3 {\varvec{P}}}/{2 E_{P}} \, dT} = \frac{2 \alpha ^2}{Q^6} L_{\mu \, \nu } W^{\mu \, \nu }_H (T). \end{aligned}$$

(67)

Since the coupling of QED is much smaller than $\alpha _S$, the leptonic tensor can be well approximated by its lowest order:

$$\begin{aligned} L^{\mu \, \nu } = l_1^\mu l_2^\nu + l_2^\mu l_1^\nu - g^{\mu \, \nu } l_1 \cdot l_2, \end{aligned}$$

(68)

where $l_1$ and $l_2$ are the momenta of the incoming electron and positron, and the electron mass is neglected.

The hadronic tensor $W^{\mu \, \nu }_H$ depends on the momentum P of the outgoing hadron and on the momentum q of the boson connecting the initial with the final state. Furthermore, it depends on thrust, T. Its definition is:

$$\begin{aligned}&W^{\mu \, \nu }_H(P,\,q,\,T)\nonumber \\&\quad = 4 \pi ^3 \, \sum _X \delta ^{(4)} \left( p_X +P - q \right) \, \nonumber \\&\quad \quad \times \langle 0 | \, j^\mu (0) | \, P,\,X,\,\hbox {out} \, \rangle _T \, {}_T \langle P,\,X,\,\hbox {out} | j^\nu (0) | \, 0 \rangle \nonumber \\&\quad = \frac{1}{4\pi } \, \sum _X \, \int d^4 z\, e^{i q \cdot z} \langle 0 | \, j^\mu \left( z/2\right) | \, P,\,X,\,\hbox {out} \, \rangle _T \, {}_T \nonumber \\&\quad \quad \times \langle P,\,X,\,\hbox {out} | j^\nu \left( -z/2\right) | \, 0 \rangle , \end{aligned}$$

(69)

where $j^\mu$ are the electromagnetic currents for the hadronic fields and the final states have been labeled by “T” in order to recall that their topology is fixed by the value of thrust. The factor $1/(4\pi )$ in the last line coincides with the normalization choice of Ref. [2]. The hadronic tensor can be decomposed in terms of structure functions:

$$\begin{aligned} W^{\mu \, \nu }_H&= \left( -g^{\mu \, \nu } + \frac{q^\mu q^\nu }{q^2} \right) F_{1,\,H} \, \nonumber \\&\quad + \frac{\left( P^\mu - q^\mu \frac{P \cdot q}{q^2} \right) \left( P^\nu - q^\nu \frac{P \cdot q}{q^2} \right) }{P \cdot q} \, F_{2,\,H}. \end{aligned}$$

(70)

Thanks to this decomposition and by using the definition of the fractional energy $z = 2 \, {P \cdot q}/{Q^2}$, see Eq. (148), we can easily compute the projections:

$$\begin{aligned} -g_{\mu \,\nu }W^{\mu \, \nu }_H&= 3 F_{1,\,H} - \left( \frac{2}{z} \, \frac{M^2}{Q^2} - \frac{z}{2} \right) \nonumber \\ F_{2,\,H}&= 3 F_{1,\,H} + \frac{z}{2} \, F_{2,\,H} + {\mathcal {O}} \left( \frac{M^2}{Q^2} \right) ; \end{aligned}$$

(71)

$$\begin{aligned} \frac{P_\mu P_\nu }{Q^2} W^{\mu \, \nu }_H&= \left( -\frac{M^2}{Q^2} + \left( \frac{z}{2} \right) ^2 \right) \nonumber \\&\quad \times F_{1,\,H} + \left( \frac{2}{z} \frac{P^4}{Q^4} - z \frac{P^2}{Q^2} + \left( \frac{z}{2} \right) ^3 \right) \, F_{2,\,H} \nonumber \\&= \left( \frac{z}{2} \right) ^2 \, \left[ F_{1,\,H} + \frac{z}{2} \, F_{2,\,H} \right] + {\mathcal {O}} \left( \frac{M^2}{Q^2} \right) . \end{aligned}$$

(72)

Fig. 5

a Leading momentum regions for the hadronic tensor $W^{\mu \, \nu }_H$ as they appear in the first step of factorization. b Actual representation of the leading momentum regions of the hadronic tensor $W^{\mu \, \nu }_H$. All the collinear factors corresponding to real emissions have been included into the hard part. The soft factor of the process is equal to one, as expected for a ${1}$-hadron class process (see Sect. 4)

6.1 Factorization of the hadronic tensor

The factorization procedure allows to factorize the hadronic tensor $W^{\mu \, \nu }_H$ into hard, collinear and soft parts, as shown in Fig. 5a. According to Ref. [2] and by using dimensional regularization, we have:

$$\begin{aligned}&W^{\mu \, \nu }_H(P,\,q,\,T) \nonumber \\&\quad = \sum _{N \ge 2} \, \sum _{j_1} \int \frac{d^D k_{1}}{(2 \pi )^D} \, \nonumber \\&\quad \quad \times \sum _{j_2} \int \frac{d^D k_{2}}{(2 \pi )^D} \prod _{\alpha = 3}^N \, \int \frac{d^D k_{\alpha }}{(2 \pi )^D} {\mathbb {C}}_{\alpha }(k_\alpha )_{j_\alpha } \nonumber \\&\quad \quad \times \hbox {Tr}_D \Big \{ {\mathcal {P}}_1 {\mathbb {C}}_{1}(k_1,\,P)_{j_1,\,H} \overline{{\mathcal {P}}}_1 {\mathbb {H}}_{j_1,\dots j_N}^\mu ({\widehat{k}}_1,\dots {\widehat{k}}_N,\,T) \, \nonumber \\&\quad \quad \times {\mathcal {P}}_2 {\mathbb {C}}_{2}(k_2)_{j_2} \overline{{\mathcal {P}}}_2 \, ({\mathbb {H}}^\dagger )_{j_1\, \dots j_N}^\nu ({\widehat{k}}_1,\dots {\widehat{k}}_N,\,T) \Big \} \nonumber \\&\quad \quad \times \int \frac{d^D k_{S}}{(2 \pi )^D} {\mathbb {S}}_{\mathrm{{N}-h}}{}_{j_1\, \dots j_N}(k_S) \, \delta ^{(n)}\nonumber \\&\quad \quad \times \left( q - {\widehat{k}}_1 - {\widehat{k}}_2 - \sum _{\alpha } {\widehat{k}}_\alpha \right) . \end{aligned}$$

(73)

In Eq. (73) the collinear parts are represented by ${\mathbb {C}}_j$: they depend only on the entering total momentum $k_j$ and on the type j of the corresponding parton (either a gluon or a quark/antiquark of flavor j/$\bar{j}$), and they are averaged over the color of the initiating parton. Among them, ${\mathbb {C}}_1$ and ${\mathbb {C}}_2$ are associated to the fermionic legs of the quark and the antiquark, hence they appear associated to the fermionic projectors, ${\mathcal {P}}_j$ and $\overline{{\mathcal {P}}}_j$, which connect them to the hard parts and make the jet partons on-shell. Since the hard part and the collinear parts are computed in the same frame (the h-frame, as defined in Appendix C) the expressions for these projectors are simply

$$\begin{aligned} {\mathcal {P}} = \frac{\gamma ^- \, \gamma ^+}{2} \quad \text {and}\quad \overline{{\mathcal {P}}} = \frac{\gamma ^+ \, \gamma ^-}{2}. \end{aligned}$$

(74)

Furthermore, by charge conjugation, $j_2 = {\overline{j}}_1$. The projectors defined above will be fundamental in extracting the leading twist FFs of the quark and the anti-quark in the cross section. All the other collinear parts, ${\mathbb {C}}_\alpha$, are generated by gluons. In this case, the role of the fermionic projectors of Eq. (74) is played by a gluon density matrix $\rho _{j'\,j}$ that encodes the information about the gluon polarization. In the following, we will consider the case of a fragmenting quark, corresponding to the collinear part ${\mathbb {C}}_1$ as depicted in Fig. 5.

In Eq. (73) the hard parts are represented by ${\mathbb {H}}$ and its hermitian conjugate, ${\mathbb {H}}^\dagger$: they encode the kinematics of the process. Momentum conservation is ensured by the appropriate delta function. However, in the hard contributions, the parton momenta are approximated, in that only their leading components are considered, as stressed by the “$^\wedge$” hats on them. In practice, the momentum ${\widehat{k}}_\alpha$ is $k_\alpha$ projected onto the (unknown) direction of its corresponding collinear part:

$$\begin{aligned} {\widehat{k}}_\alpha = w_{\alpha } \, \frac{k_\alpha \cdot {\widetilde{w}}_\alpha }{w_{\alpha } \cdot {\widetilde{w}}_\alpha }, \end{aligned}$$

(75)

where $w_{\alpha }$ and ${\widetilde{w}}_\alpha$ are the light-like vectors corresponding to the plus and the minus directions, respectively, in the reference frame of ${\mathbb {C}}_\alpha$. The approximated momentum of the fragmenting quark is simply:

$$\begin{aligned} {\widehat{k}}_1 = \left( {\widehat{k}}_{1,\,h}^+,\,0,\,{{\varvec{0}}}_T \right) _h, \end{aligned}$$

(76)

where ${\widehat{k}}_{1,\,h}^+ = k_{1,\,h}^+$, since the reference frame of the fragmenting parton corresponds (by definition) with the hadron frame. Furthermore, kinematics impose constraints on the possible values that $k_{1,\,h}^+$ can assume, since $P_h^+< k_{1,\,h}^+ < {P_h^+}/{z}$ (see Eqs. (144), (146) and (147)).

Finally, the soft contribution is a N-h soft factor, where N is the total number of partons exiting the hard scattering. It depends on the collinear parton type only through their color representation. Notice that the total soft momentum $k_S$ cannot be involved in the kinematics of the process, since it is washed out by the real hard emission (at least one, ${\mathbb {C}}_2$). In fact, none of the $k_S$ components appear in the conservation delta. As a consequence, ${\mathbb {S}}_{\mathrm{{N}-h}}$ is integrated over all the components of $k_S$ and its contribution becomes trivial. As expected for a process belonging to the ${1}$-hadron class, the soft factor is unity and can be omitted in the leading region representation. This is a common feature in a collinear factorization scheme, see Ref. [2], but it requires a more careful treatment in this case. In fact, the cross section we are dealing with is also sensitive to the value of thrust, which provides a precise indication about the topology of the final state. In other words, we know exactly the total number N of partons that generate the collinear blobs ${\mathbb {C}}_i$ in Fig. 5a. For instance, if $T \sim 1$, then the topology of the final state is 2-jet like, and only ${\mathbb {C}}_1$ and ${\mathbb {C}}_2$ are allowed to be represented as leading regions. Despite its strong influence on the whole process, the explicit information about T is solely encoded into the hard parts. Therefore, the soft factor ${\mathbb {S}}_{\mathrm{{N}-h}}$ is totally blind to the precise value of T (see the definition of Eq. (4)) hence, ultimately, it can be integrated out. It is important to stress, however, that this does not imply that the soft radiation does not contribute to the final cross section. Rather, it means that the contribution of the soft gluons has to be considered in association to the information on thrust. In fact, it turns out that a soft function S(T), explicitly depending on T and totally predicted by perturbative QCD, will be one of the ingredients appearing in the final cross section. An explicit computation of such soft function is presented in Ref. [28], for the case of a 2-jet like topology of the final state.

A similar argument applies for all the collinear parts except ${\mathbb {C}}_1$, which is the only one associated to the final, detected hadron H that makes the sum over the final states incomplete. In all other cases, the blobs ${\mathbb {C}}_{i \ne 1}$ are totally crossed by the final state cut (see Fig. 5a). According to a collinear factorization scheme, such contributions are suppressed in the collinear region they are supposed to describe, and can actually be considered hard contributions. Again, this does not imply that the radiation collinear to other direction than that of the detected hadron does not give any contribution to the final cross section. Rather, it means that such contributions cannot be disentangled from the information about thrust. Indeed, they will be represented in the final cross section by jet functions $J_i(T)$, explicitly depending on T and totally predicted by perturbative QCD. An explicit example is presented in Ref. [28], for the case of a 2-jet like topology of the final state.

Following the previous argument, all the soft and collinear contributions (except ${\mathbb {C}}_1$) can be recast, together with the pure hard vertices ${\mathbb {H}}$ and ${\mathbb {H}}^\dagger$, into one single factor, as depicted in Fig. 5b. The final result is then given by:

$$\begin{aligned}&W^{\mu \, \nu }_H(P,q,T) \nonumber \\&\quad = \sum _{j_1} \int d {\widehat{k}}_{1,\,h}^+ \int \frac{d k_{1,\,h}^- \, d^{D-2}{{\varvec{k}}}_{1,\,T,\,h}}{(2 \pi )^D} \, \hbox {Tr}_D \nonumber \\&\qquad \times \Big \{ {\mathcal {P}}_1 {\mathbb {C}}_{1}(k_1,\,P)_{j_1,\,H} \overline{{\mathcal {P}}}_1 {\mathcal {H}}_{j_1}^{\mu \, \nu }(Q,\,{\widehat{k}}_{1,\,h}^+,\,T) \Big \}. \end{aligned}$$

(77)

In the above equation, all the contributions that can be totally predicted by perturbative QCD have been collected in the hard coefficient ${\mathcal {H}}^{\mu \, \nu }$. Notice that, here, “hard” indicates that ${\mathcal {H}}^{\mu \, \nu }$ is completely computable by using perturbative techniques; hence it can be considered a “hard contribution” according to the naive classification introduced in Sect. 1.1. In fact, despite the label “hard”, it encodes all the thrust dependent soft and collinear terms that have been recast into the large blob of Fig. 5b. Notice that, while the collinear part ${\mathbb {C}}_1$ depends on all the components of $k_1$, the hard contribution depends only on its leading component, $k_{1,\,h}^+$. Then, ${\mathbb {C}}_1$ and ${\mathcal {H}}^{\mu \, \nu }$ are not completely disentangled, because a convolution over $k_{1,\,h}^+$ will survive. In the following we will drop the index “1” related to the fragmenting parton, which has become redundant. Applying the fermionic projectors and parity conservation, the only surviving contribution in the case of $e^+e^-\rightarrow H\,X$ is given by the coefficient of $\gamma ^-$ in the expansion of Eq. (23):

$$\begin{aligned} {\mathcal {P}} {\mathbb {C}}(k,\,P)_{j,\,H} \overline{{\mathcal {P}}} = \gamma ^- \; \frac{\hbox {Tr}_D}{4} \Big \{ \gamma ^+ {\mathbb {C}}(k,\,P)_{j,\,H} \Big \}. \end{aligned}$$

(78)

The Dirac trace of $\gamma ^+ {\mathbb {C}}(k,\,P)_{j,\,H}$ defines two TMD FFs (as in Eq. (24)):

$$\begin{aligned}&\frac{1}{{\widehat{z}}} \, \int \frac{d k_{h}^-}{(2 \pi )^D} \frac{\hbox {Tr}_D}{4} \Big \{ \gamma ^+ {\mathbb {C}}(k,\,P)_{j,\,H}\Big \} \nonumber \\&\quad = D_{1,\,H/j}({\widehat{z}},\,\vert -{\widehat{z}} \, {{\varvec{k}}}_{h,\,T} \vert ) \nonumber \\&\qquad - \frac{1}{M} |{{\varvec{S}}}_T \times {{\varvec{k}}}_{h,\,T}| D_{1T,\,H/j}^{\perp }({\widehat{z}},\,\vert -{\widehat{z}} \, {{\varvec{k}}}_{h,\,T} \vert ), \end{aligned}$$

(79)

where M and ${{\varvec{S}}}_T$ are the mass and the transverse spin of the detected hadron, while ${\widehat{z}} = {P_h^+} / {k_h^+}$. The function $D_{1,\,H/j ^\perp }$ is the unpolarized TMD FF, while $D_{1T,\,H/j}$ is the Sivers-like TMD FF. For simplicity, in the following we will collectively indicate with $D_{j,\,H}({\widehat{z}},\,-{\widehat{z}} \, {{\varvec{k}}}_{h,\,T})$ the sum of the two contributions in the r.h.s. of Eq. (79). Therefore:

$$\begin{aligned} W^{\mu \, \nu }_H(P,\,q,\,T)&= \sum _{j} \int d {\widehat{k}}_h^+ \; {\widehat{z}} \nonumber \\&\quad \times \int d^{D-2}{{\varvec{k}}}_{h,\,T} \; D_{j,\,H}({\widehat{z}},\,-{\widehat{z}} \, {{\varvec{k}}}_{h,\,T})\, \nonumber \\&\quad \times \hbox {Tr}_D\Big \{ \gamma ^- {\mathcal {H}}^{\mu \, \nu }(Q,\,{\widehat{k}}_h^+,\,T)_j\Big \} \nonumber \\&= \sum _{j} \int _z^1 \, d {\widehat{z}} \; {\widehat{W}}_j^{\mu \, \nu }({z}/{\widehat{z}},\,Q,\,T) \nonumber \\&\quad \times \int d^{D-2}{{\varvec{k}}}_{h,\,T} \; D_{j,\,H}({\widehat{z}},\,-{\widehat{z}} \, {{\varvec{k}}}_{h,\,T}) , \end{aligned}$$

(80)

where the kinematics constraints over ${\widehat{z}}$ have been taken into account. The role of the hard factor in the previous equation is played by the function ${\widehat{W}}_j^{\mu \, \nu }$, which is the partonic analogue of the full hadronic tensor $W^{\mu \, \nu }_H$. It is defined as:

$$\begin{aligned} {\widehat{W}}_j^{\mu \, \nu }({\widehat{k}},\,q,\,T) = \hbox {Tr}_D \Big \{ {\widehat{k}}_h^+ \gamma ^- {\mathcal {H}}_j^{\mu \, \nu }(Q,\,{\widehat{k}}_h^+,\,T)\Big \}, \end{aligned}$$

(81)

Notice that since the approximated parton momentum has only a plus component, we can write . Therefore, ${\widehat{W}}_j$ is the algebraic expression corresponding to the pictorial representation given in Fig. 6a. Its actual definition has to be equipped with the subtraction of the double counting due to the overlapping with the collinear momentum region (see Sect. 6.3). As already stressed, the label “hard” associated to ${\widehat{W}}_j$ only refers to its perturbative nature, since it is totally predicted by perturbative QCD. It encodes all the contributions related to soft and collinear radiation, properly expressed by soft and jet functions that explicitly depend on the thrust T. In fact, the factorization procedure can be further applied to ${\widehat{W}}_j$, in order to expose all the soft/collinear terms which are implicitly written in Eq. (81). This is shown in Fig. 6b. In the following, we will not use the factorized form for the partonic final state tensor, since we are not interested in any particular final state topology. An explicit example, valid in the 2-jet case, is presented in Ref. [28].

Fig. 6

a Pictorial representation of ${\widehat{W}}_j^{\mu \, \nu }$. b Factorization of ${\widehat{W}}_j^{\mu \, \nu }$, in which soft and collinear contributions are explicitly represented, according to the topology of the final state associated to the value of the thrust T. All the quantities are totally predicted by perturbative QCD

In Eq. (80), the dependence on the parton transverse momentum is only in the collinear part and, in principle, the integrand of $W^{\mu \, \nu }_H$ could be defined as the hadronic tensor differential in ${{\varvec{k}}}_{h,\,T}$. However, although the parton transverse momentum is not a physical observable, kinematics relates ${{\varvec{k}}}_{h,\,T}$ with the transverse momentum ${{{\varvec{P}}}}_{p,\,T}$ of the outgoing hadron in the parton frame, i.e. measured with respect to its final state jet axis, that we identify with the thrust axis (see Appendix C). This can be measured (as it has been done by the BELLE Collaboration, Ref. [1]) and the definition of the hadronic tensor differential in $P_{p,\,T}$ is obtained by the change of variables ${{\varvec{k}}}_{h,\,T} = -\frac{1}{{\widehat{z}}} \; {{{\varvec{P}}}}_{p,\,T}\left[ 1 + {\mathcal {O}}( \frac{P_{p,\,T}^2}{Q^2} )\right]$ (see Eq. (154)). Therefore:

$$\begin{aligned}&\frac{d W_H^{\mu \, \nu }(z,\,Q,\,T)}{d^2 {{{\varvec{P}}}}_{p,\,T}} \nonumber \\&\quad = \sum _{j} \int \frac{d {\widehat{z}}}{{\widehat{z}}^{2}} \, {\widehat{W}}_j^{\mu \, \nu }({z}/{\widehat{z}},\,Q,\,T) \, D_{j,\,H}({\widehat{z}},\, {{{\varvec{P}}}}_{p,\,T}) \nonumber \\&\qquad \times \left[ 1 + {\mathcal {O}}\left( \frac{P_{p,\,T}^2}{Q^2} \right) \right] . \end{aligned}$$

(82)

The differential of ${{{\varvec{P}}}}_{p,\,T}$ carries information about two variables: the modulus $P_{p,\,T}$ and the azimuthal angle $\beta$ in the $x\,y$-plane of the parton frame. While the first can be measured, the angle $\beta$ cannot be determined experimentally. In fact, an angular dependence in the TMD contribution $D_{j,\,H}$ can originate from the Sivers-like contribution $|{{\varvec{S}}}_T \times {{{\varvec{P}}}}_{p,\,T}|$ (see Eq. (79)). However, as explained in Ref. [17], the transverse spin of the hadron is orthogonal to its transverse momentum with respect to the axis of the jet, identified with the thrust axis. Hence $|{{\varvec{S}}}_T \times {{{\varvec{P}}}}_{p,\,T}| = \pm S_T \, P_{p,\,T}$ for any choice of the x-axis in the parton frame. Therefore, the integration over $\beta$ is trivial and results just in a $2\pi$ factor on the r.h.s of Eq. (82):

$$\begin{aligned}&\frac{d W_H^{\mu \, \nu }(z,\,Q,\,T)}{d P^2_{p,\,T}} \nonumber \\&\quad = \pi \, \sum _{j} \int \frac{d {\widehat{z}}}{{\widehat{z}}^{2}} \, {\widehat{W}}_j^{\mu \, \nu }({z}/{\widehat{z}},\,Q,\,T) \, D_{j,\,H}({\widehat{z}},\, P_{p,\,T}) \nonumber \\&\qquad \times \left[ 1 + {\mathcal {O}}\left( \frac{P_{p,\,T}^2}{Q^2}\right) \right] , \end{aligned}$$

(83)

where:

$$\begin{aligned} D_{j,\,H}({\widehat{z}},\, P_{p,\,T})&= D_{1,\,H/j}({\widehat{z}},\,P_{p,\,T}) \nonumber \\&\quad {\mp } \frac{{\widehat{z}}}{M} S_T \, P_{p,\,T} \, D_{1T,\,H/j}^{\perp }({\widehat{z}},\,P_{p,\,T}). \end{aligned}$$

(84)

6.2 Factorized cross section

The full cross section is obtained by contracting the hadronic tensor of Eq. (83), and its partonic counterpart, Eq. (81), with the leptonic tensor, as in Eq. (67):

$$\begin{aligned}&\frac{d \sigma }{({d^3 {\varvec{P}}}/{2 E_{P}}) \, dP^2_{p,\,T} \,dT} \nonumber \\&\quad = \pi \sum _{j} \int _z^1 \frac{d {\widehat{z}}}{{\widehat{z}}^{2}} \, \frac{d {\widehat{\sigma }}_j}{d^3 {\varvec{\widehat{k}}}/{2 E_{\widehat{k}}}\,dT} \, D_{j,\,H}({\widehat{z}},\, P_{p,\,T}) \,\nonumber \\&\qquad \times \left[ 1 + {\mathcal {O}}\left( \frac{P_{p,\,T}^2}{Q^2} \right) \right] , \end{aligned}$$

(85)

where the dependence on thrust has been made explicit. Let’s focus on the r.h.s. of the previous equation. The only non-zero component of the approximated parton momentum ${\widehat{k}}$ is in the plus direction, as defined in Eq. (76). Therefore, its Lorentz invariant phase space measure can only be written as:

$$\begin{aligned} d^3 {\varvec{\widehat{k}}}/{2 E_{\widehat{k}}}&= \frac{1}{2} d |{\widehat{{{\varvec{k}}}}}| \, |{\widehat{{{\varvec{k}}}}}|\, d\cos {\theta } \, d \phi \nonumber \\&= \frac{Q^2 }{8} \,\frac{z}{{\widehat{z}}} \, d \left( \frac{z}{{\widehat{z}}} \right) \, d\cos {\theta } \, d \phi , \end{aligned}$$

(86)

and carries information about the polar angle $\theta$ and the azimuthal angle $\phi$ with respect to the beam axis (LAB frame, see Appendix C). On the l.h.s the same variables have to appear explicitly. Hence, the Lorentz invariant phase space of the detected hadron is written in the LAB frame as well:

$$\begin{aligned}&{d^3 {\varvec{P}}}/{2 E_{P}} = \frac{Q^2 }{8} \,z \, d z \, d\cos {\theta } \, d \phi \, \left[ 1 + {\mathcal {O}} \left( \frac{M^2}{Q^2} \right) \right] . \end{aligned}$$

(87)

Finally, the cross section is given by:

$$\begin{aligned}&\frac{d \sigma }{dz \,d\cos {\theta } \,d \phi \, dP^2_{p,\,T} \,dT} \nonumber \\&\quad = \pi \sum _{j} \int _z^1 \frac{d {\widehat{z}}}{{\widehat{z}}} \, \frac{d {\widehat{\sigma }}_j}{d({z}/{{\widehat{z}}})\,d\cos {\theta } \,d \phi \,dT} \, D_{j,\,H}({\widehat{z}},\, P_{p,\,T}) \, \nonumber \\&\qquad \times \left[ 1 + {\mathcal {O}}\left( \frac{P_{p,\,T}^2}{Q^2} ,\; \frac{M^2}{Q^2}\right) \right] . \end{aligned}$$

(88)

There are five independent observables:

• 1. The fractional energy $z = {2 |{{{\varvec{P}}}}|}/{Q}$.

• 2. The polar angle $\theta$ of the outgoing hadron with respect to the electron.

• 3. The azimuthal angle $\phi$ of the outgoing hadron with respect to the x-axis in the LAB frame. This is significant only if such axis can be defined unambiguously, as in the case of polarized leptons. Otherwise, we can simply drop $d\phi$ on both sides of Eq. (88) as a result of integration, which is our case.

• 4. The thrust T, defined in Eq. (156).

• 5. The (modulus of the) transverse momentum of the outgoing hadron $P_{p,\,T}$ with respect to its final state jet axis, that we identify with the thrust axis.

Common scenarios are those in which experiments provide two or three of the variables listed above:

• z and $\theta$ are measured, but the thrust axis is not reconstructed, hence $P_{p,\,T}$ is unknown. In this case, in addition to the integration over T, the previous cross section has to be integrated over all possible values of the transverse momentum $P_{p,\,T}$, restoring the integrated TMDs as in Eq. (80):

  $$\begin{aligned} \frac{d \sigma }{dz \,d\cos {\theta }}&= \sum _{j} \int _z^1 \frac{d {\widehat{z}}}{{\widehat{z}}} \, \frac{d {\widehat{\sigma }}_j}{d {z}/{{\widehat{z}}} \, d\cos {\theta }} \, d_{j/H}({\widehat{z}})\,\nonumber \\&\quad \times \left[ 1 + {\mathcal {O}}\left( \frac{M^2}{Q^2}\right) \right] , \end{aligned}$$

  (89)

  where we used the results of Appendix B.2. This result coincides with the cross section presented in Chapter 12 of Ref. [2]. The convolution over ${\widehat{z}}$ is between renormalized quantities, as we did for the OPE of TMDs at small $b_T$ in Eq. (140).

  The dependence on $\theta$, both in the partonic and in the full cross section, can expressed in terms of longitudinal (L) and transverse (T) contributions:

  $$\begin{aligned} \frac{d \sigma }{d x \, d \cos {\theta }} = \frac{3}{8} (1 + \cos ^2{\theta }) \frac{d \sigma _T}{d x} + \frac{3}{4} \sin ^2{\theta } \frac{d \sigma _L}{d x}\,, \end{aligned}$$

  (90)

  where x can be z in the full cross section, or ${z}/{{\widehat{z}}}$ in its partonic counterpart. The structure functions are related to the transverse and the longitudinal component of the cross section as follows:

  $$\begin{aligned}&\frac{d \sigma _T}{d x} = \frac{4 \pi \alpha ^2}{3 Q^2} x F_1(x, Q^2), \end{aligned}$$

  (91)
  $$\begin{aligned}&\frac{d \sigma _L}{d x} = \frac{\pi \alpha ^2}{3 Q^2} \left[ 2 x F_1(x, Q^2) + x^2 F_2(x, Q^2)\right] . \end{aligned}$$

  (92)

• z and $P_{p,\,T}$ are measured, but the polar angle $\theta$ of the outgoing hadron with respect to the beam axis is integrated over. Indeed, the measurement of the transverse momentum $P_{p,\,T}$ has to be done with respect to the jet axis, which for our purposes coincides with the thrust axis. Therefore, if the cross section is differential in $P_{p,\,T}$, it also has to be differential in T (or in an analogous variable that allows to determine the axis of the jet).⁸

  The integration of the partonic cross section with respect to $\theta$ is straightforward and follows from Eqs. (90), (91) and (92):

  $$\begin{aligned}&\int _{-1}^1 \, d \cos {\theta } \, \frac{d {\widehat{\sigma }}_j}{d x \, d \cos {\theta } \, dT} = \frac{d \sigma _T}{d z} + \frac{d \sigma _L}{d z} \nonumber \\&\quad = \frac{4 \pi \alpha ^2}{3 Q^2} \, x \, \left( \frac{3}{2} \, F_{1,\,j}(x,\,Q^2,\,T) + \frac{x}{4} \, F_{2,\,j}(x,\,Q^2,\,T)\right) . \end{aligned}$$

  (93)

  For simplicity, in the following we will collectively indicate with ${d\sigma }/{d x}$ the sum of the two contributions on the r.h.s. of Eqs. (93). Therefore:

  $$\begin{aligned} \frac{d \sigma }{dz\, dP^2_{p,\,T} \,dT}&= \pi \sum _{j} \int _z^1 \frac{d {\widehat{z}}}{{\widehat{z}}} \, \frac{d {\widehat{\sigma }}_j}{d({z}/{{\widehat{z}}})\,dT} \, D_{j,\,H}({\widehat{z}},\, P_{p,\,T}) \,\nonumber \\&\quad \times \left[ 1 + {\mathcal {O}}\left( \frac{P_{p,\,T}^2}{Q^2} ,\; \frac{M^2}{Q^2}\right) \right] . \end{aligned}$$

  (94)

  Since TMDs are defined in the Fourier conjugate space, see Eq. (25), it is more convenient to write the cross section using their $b_T$-space counterparts:

  $$\begin{aligned}&\int d^{D-2} {{{\varvec{P}}}}_{p,\,T} \; e^{i \, {{{\varvec{P}}}}_{p,\,T} \cdot {{\varvec{b}}}_T} D_{j,\,H}(z, {{{\varvec{P}}}}_{p,\,T}) \, \nonumber \\&\qquad \times \left[ 1 + {\mathcal {O}}\left( \frac{P_{p,\,T}^2}{Q^2}\right) \right] \nonumber \\&\quad = z^{D-2} \int d^{D-2} {{\varvec{k}}}_{T,\,h} \; e^{i \, {{\varvec{k}}}_{T,\,h} \cdot (- z {{\varvec{b}}}_T)} D_{j,\,H}(z, -z {{\varvec{k}}}_{T,\,h}) \nonumber \\&\quad = z^{D-2} {\widetilde{D}}_{j,\,H}(z, \, - z \, {{\varvec{b}}}_T), \end{aligned}$$

  (95)

  where $D_{j,\,H}$ is actually only a function of the modulus of ${{{\varvec{P}}}}_{p,\,T}$. Hence:

  $$\begin{aligned}&D_{j,\,H}(z, P_{p,\,T}) \, \left[ 1 + {\mathcal {O}}\left( \frac{P_{p,\,T}^2}{Q^2}\right) \right] \nonumber \\&\quad = \int \frac{d^{2} {{\varvec{b}}}_T}{(2 \pi )^{2}} \, e^{i \, \frac{{{{\varvec{P}}}}_{p,\,T}}{z} \cdot {{\varvec{b}}}_T} \, {\widetilde{D}}_{j,\,H}(z, \,b_T). \end{aligned}$$

  (96)

  Notice that all definitions in Eqs. (25) and (38) hold for the Fourier transformed TMD FFs ${\widetilde{D}}_{j,\,H}$. Finally, the cross section in its final form is given by:

  $$\begin{aligned} \frac{d \sigma }{dz\, dP^2_{p,\,T} \,dT}&= \pi \sum _{j} \int _z^1 \frac{d {\widehat{z}}}{{\widehat{z}}} \, \frac{d {\widehat{\sigma }}_j}{d({z}/{{\widehat{z}}})\,dT} \,\nonumber \\&\quad \times \int \frac{d^{2} {{\varvec{b}}}_T}{(2 \pi )^{2}} \, e^{i \, \frac{{{{\varvec{P}}}}_{p,\,T}}{{\widehat{z}}} \cdot {{\varvec{b}}}_T} \, {\widetilde{D}}_{j,\,H}({\widehat{z}}, \,b_T) \, \nonumber \\&\quad \times \left[ 1 +{\mathcal {O}}\left( \frac{M^2}{Q^2}\right) \right] . \end{aligned}$$

  (97)

  In this final cross section, all the non-perturbative information about the hadronization is contained in the TMD FFs. On the other hand, the partonic cross section is totally predicted by perturbative QCD, despite it encodes the contributions associated to both the soft radiation and the radiation collinear to the directions of the observed jets. In fact, the factorization of the partonic final state tensor depicted in Fig. 6b automatically extends to the full partonic cross section. In this regard, shortly after the publication of this work, the transverse momentum spectrum of single-hadron production in $e^+e^-$ has been investigated in two papers, Refs. [29, 30], both based on a SCET treatment. In this framework, all soft/collinear contributions appear explicitly in the final cross section, in contrast to Eq. (97) in which they are implicitly contained in the partonic cross section. The underlying structure of the cross section formulae presented in these articles are consistent with our results: an explicit comparison up to NLO may be obtained by using the results of Ref. [28].

  As for the cross section in Eq. (89) the convolution in Eq. (97) is between renormalized quantities, as we will discuss in the next Section. Furthermore, in contrast to Eq. (94), the cross section written in terms of the Fourier transform is a function defined for any value of $P_{p,\,T}$ and in fact the errors are only sized as ${M^2}/{Q^2}$. However, the physical meaning is lost for large values of the transverse momentum of the outgoing hadron, as the TMDs themselves become non physical in the large $P_{p,\,T}$ region (see the discussion at the end of Appendix B.2). Hence, the cross section of Eq. (97) can only be trusted where $P_{p,\,T} \ll Q$ or, more precisely, where $P_{p,\,T} \ll P^+ = z \, {Q}/{\sqrt{2}}$, which is the actual condition that allows to consider the outgoing hadron as a collinear particle, according to the power counting rules.

6.3 Subtraction mechanism

As it is clear from Eq. (97), the $e^+e^-\rightarrow HX$ cross section, differential in z, in thrust T and in the transverse momentum of the detected hadron with respect to the thrust axis, $P_{p,\,T}$, offers a direct probe of the transverse motion of partons. Recently the BELLE Collaboration has provided high statistics experimental data corresponding to such cross section [1]. Although the final result of Eq. (97) is simply the Fourier Transform of the convolution of a TMD FF and a thrust-dependent hard factor, i.e. the partonic cross section integrated over $\theta$, the phenomenological application of Eq. (97) requires special care.

The final cross section is RG invariant if the anomalous dimension of the hard factor is exactly equal and opposite to that of the TMD FF, order by order in perturbation theory. This argument applies to renormalized quantities, i.e. functions provided of the proper UV counterterm. Furthermore, the hard factor in Eq. (97) has to be properly subtracted to avoid double counting due to the overlapping with the collinear momentum region. Therefore, the hard factor of the final cross section is defined in two steps: first it is equipped with subtractions, then it is renormalized.

The unsubtracted analogue of the hard factor in Eq. (97) is the partonic version of the full cross section. Being a partonic quantity, it is completely unaware of the outgoing hadron. It describes the process at partonic level, which means $e^+e^-\rightarrow f \, X$, where f is a parton of type f that replaces the detected hadron. The most convenient frame where to compute ${\widehat{\sigma }}^{\mathrm{unsub}}$ is the analogue of the hadron frame, where the momentum of f lies along the plus direction. The expression of its final state tensor ${\widehat{W}}_f^{\mu \,\nu }{}^{\mathrm{, unsub}}$ is obtained from the integrand of Eq. (80). For a given value of T, the phase space available for real emissions is restricted, because only the final state topology associated with that particular value of thrust can be reached. A simple way to force the phase space to describe only the region of interest is introducing sharp cut-offs that shrink the available range of values of T. For instance, if we are interested in the quasi 2-jet limit, we can force T to remain in the neighborhood of 1 by defining a minimal value of T, $T_{\mathrm{MIN}} \le 1$. In practice, the unsubtracted final state tensor is obtained by computing the contribution of all the Feynman graphs needed when T lies in the range defined by the topology cut-offs, in the massless limit and by setting all the soft/collinear divergent quantities to their lowest order (in the language of Ref. [2], this corresponds to the application of the hard approximator $T_H$, modified to include the introduction of the cut-offs for thrust). Since the lowest order for the collinear part is just a product of delta function that sets the momentum of the fragmenting parton to be equal to that of the outgoing parton, and the lowest order for the soft factor is unity, in $k_T$-space we simply have:

$$\begin{aligned}&\frac{d {\widehat{W}}_f^{\mu \,\nu }{}^{\mathrm{, unsub}}(\epsilon ;\,z,\,T;\,T_c)}{d^{2-2\epsilon } {{\varvec{k}}}_T} \nonumber \\&\quad =\sum _j \, \int _z^1 \frac{d {\widehat{z}}}{{\widehat{z}}} \, {\widehat{W}}_j^{\mu \,\nu }{}^{\mathrm{, unsub}}(\epsilon ;\,{z}/{{\widehat{z}}},\,T;\,T_c) \, \delta _{j\,f} \, \delta (1-{\widehat{z}}) \, \delta ^{2-2\epsilon }({{\varvec{k}}}_T), \end{aligned}$$

(98)

where $T_c$ stands for a generic topology cut-off for thrust. As a consequence, the Fourier transform of the previous expression does not depend on $b_T$. In Eq. (98) we explicitly showed the dependence on $\epsilon$, which is the regulator used in dimensional regularization (where the spacetime dimension is set to $D = 4 - 2\epsilon$). In fact, the unsubtracted final state tensor is collinear divergent and presents poles in $\epsilon$. The standard subtraction procedure removes the overlapping with the momentum region described by the collinear part, that represents the boundary of the phase space corresponding to the emissions along the direction of the outgoing parton f, and hence it also removes the collinear divergences of ${\widehat{W}}_f^{\mu \,\nu }{}^{\mathrm{, unsub}}$. The subtraction term is obtained by considering the collinear approximation of the unsubtracted final state tensor and it coincides with the partonic version of the (bare) TMD FFs, that will be denoted by $D^{(0)}_{f/j}$ (in the language of Ref. [2], the subtraction term is obtained by applying $T_A T_H$, i.e. applying both the approximator collinear to the outgoing particle and the hard approximator). However, due to the presence of topology cut-offs, represented by $T_c$ in Eq. (98), we have to slightly modify the subtraction mechanism. In particular, it is not the whole partonic TMD has to be subtracted out, but only the part that actually overlaps. This coincides with the contribution given by the transverse momenta that lie in the power counting momentum region, i.e. where $k_T$ is at most of order $\lambda$, with $\lambda$ being some IR energy scale. Order by order we can relate $\lambda$ and $T_c$ by a precise kinematic relation. As a consequence, the Fourier transform of the partonic TMD cannot cover the whole range of $k_T$, but it stops when $k_T$ reaches $\lambda$. The resulting quantity does not depend on $b_T$, nevertheless it shows an explicit dependence on the transverse momentum cut-off and, ultimately, on $T_c$. Summarizing, the (bare) Fourier transformed subtraction term is defined by:

$$\begin{aligned}&{\widetilde{D}}^{(0),\,(\lambda )}_{f/j} (\epsilon ;\,z,\,\lambda (T_c),\,\zeta ) \nonumber \\&\quad = \int d^{2-2\epsilon } {{\varvec{k}}}_T \, e^{-i \, {{\varvec{k}}}_T \cdot {{\varvec{b}}}_T} \, D^{(0)}_{f/j} (\epsilon ;\,z,\,k_T,\,\zeta )\, \theta \left( \lambda (T_c) - k_T\right) \end{aligned}$$

(99)

This quantity is both collinear and UV divergent. Since the poles of ${\widetilde{D}}^{(0),\,(\lambda )}_{f/j}$ are the same that would be obtained by a complete Fourier transform, the UV divergence is renormalized by using the same UV counterterm that heals the UV divergence in the usual TMDs (i.e. those defined without a cut-off on transverse momentum). In this regard, as a consequence of the mechanism of subtractions, the subtracted final state tensor acquires the same UV divergences of ${\widetilde{D}}^{(0),\,(\lambda )}_{f/j}$, but with opposite signs. Then we can easily renormalize it by using the inverse of the TMD UV counterterm. In $b_T$-space we have the following factorization formula:

$$\begin{aligned}&\underbrace{{\widehat{W}}_f^{\mu \,\nu }{}^{\mathrm{, unsub}} (\epsilon ; \,z, \,T, \,T_c)}_\text {coll. divergent} \nonumber \\&\quad = \sum _{j} \int _z^1 \, \frac{d {\widehat{z}}}{{\widehat{z}}} \; \underbrace{{\widehat{W}}_j^{\mu \, \nu }{}^{\mathrm{, sub,} (0)} (\epsilon ;\,{z}/{\widehat{z}},\,T,\,\lambda (T_c),\, \zeta )}_\text {UV divergent} \, \nonumber \\&\qquad \times \underbrace{{\widehat{z}} \; {\widetilde{D}}^{(0),\,(\lambda )}_{j,\,f} (\epsilon ; \, {\widehat{z}},\,\lambda (T_c),\,\zeta )}_\text {UV divergent and coll. divergent} \nonumber \\&\quad = \sum _{j} \int _z^1 \, \frac{d {\widehat{z}}}{{\widehat{z}}} \;\nonumber \\&\qquad \times \left\{ {\widehat{W}}_k^{\mu \, \nu }{}^{\mathrm{, sub,} (0)} (\epsilon ;\,{z}/{\widehat{z}},\,T,\,\lambda (T_c),\, \zeta ) Z_{\mathrm{TMD}}^{-1}{}^k_{j} (\epsilon ;\,\mu ,\,\zeta ) \right\} \nonumber \\&\qquad \times \left\{ {\widehat{z}} \, Z_{\mathrm{TMD}}{}_{j}^{l} (\epsilon ;\,\mu ,\,\zeta )\, {\widetilde{D}}^{(0),\,(\lambda )}_{l,\,f} (\epsilon ; \, {\widehat{z}},\,\lambda (T_c),\,\zeta ) \right\} \nonumber \\&\quad = \sum _{j} \int _z^1 \, \frac{d {\widehat{z}}}{{\widehat{z}}} \; {\widehat{W}}_j^{\mu \, \nu }{}^{\mathrm{, sub}} ({z}/{\widehat{z}},\,T, \,\mu ,\,\lambda (T_c), \, \zeta ) \, \nonumber \\&\qquad \times \underbrace{{\widehat{z}} \; {\widetilde{D}}^{(\lambda )}_{j,\,f} (\epsilon ; \, {\widehat{z}},\,\mu ,\,\lambda (T_c),\,\zeta )}_{\mathrm{coll. divergent}} , \end{aligned}$$

(100)

where we simply used the associative property of convolutions and a sum over repeated upper-lower flavor indices is implicit. Notice that at this stage the renormalized, subtracted final state tensor has acquired a dependence on both the topology cut-off $T_c$ and on the rapidity cut-off $\zeta$. Order by order in perturbation theory, the functions ${\widehat{W}}_j^{\mu \, \nu }{}^{\mathrm{, sub}}$ are determined recursively by using:

$$\begin{aligned}&{\widehat{W}}_j^{\mu \nu ,\,{[n]}}{}^{\mathrm{, sub}} (z,\,T,\,\mu ,\,\lambda (T_c),\, \zeta )\nonumber \\&\quad = {\widehat{W}}_f^{\mu \nu ,\,{[n]}} {}^{\mathrm{, unsub}} (\epsilon ; \,z, \,T, \,T_c)+\nonumber \\&\qquad - \sum _j \, \sum _{m = 1}^n \; \int _z^1 \frac{d{\widehat{z}}}{{\widehat{z}}}\; {\widehat{W}}_{j,\;R}^{\mu \nu ,\,{[n-m]}} {}^{\mathrm{, sub}} ({z}/{{\widehat{z}}},\,T,\,\mu ,\,\lambda (T_c),\, \zeta ) \, \nonumber \\&\qquad \times \left[ {\widehat{z}} \, {\widetilde{D}}_{j,\,f}^{[m],\,(\lambda )} (\epsilon ; \, {\widehat{z}},\,\mu ,\,\lambda (T_c),\,\zeta ) \right] , \end{aligned}$$

(101)

In the previous result, we used the fact that the lowest order of the partonic TMDs equipped with $\lambda$ is just a delta function, ${\widetilde{D}}_{f/j}^{[0],\,(\lambda )}\left( {\widehat{z}} \, \right) = \delta _{f \; j} \, \delta (1-{\widehat{z}})$, as in the case in which there is no cut-off on transverse momenta. From now on, when the labels “sub” and “R” are not explicitly indicated, ${\widehat{W}}_j$ will be implicitly considered both subtracted and renormalized. Once the expression of ${\widehat{W}}_j$ is known, the full subtracted, renormalized cross section is computed straightforwardly through the partonic structure functions ${\widehat{F}}_{1,\,j}$ and ${\widehat{F}}_{2,\,j}$, obtained as in Eqs. (71) and (72).

The partonic cross section resulting from the previous subtraction procedure obeys the following RG evolution equation:

$$\begin{aligned}&\frac{\partial }{\partial \log {\mu }} \log {\left( \frac{d {\widehat{\sigma }}_j \left( \mu ,\,\lambda (T_c),\,\zeta \right) }{dz \, dT}\right) } \nonumber \\&\quad = -\gamma _{D,\,j} \left( \alpha _S(\mu ),\,{\zeta }/{\mu ^2} \right) , \end{aligned}$$

(102)

where $\gamma _{D,\,j}$ is the anomalous dimension of the TMD FF of flavor j. The RG invariance of the full cross section follows straightforwardly, as the anomalous dimensions of the partonic cross section and of the full TMD FF appearing in Eq. (97) are equal and opposite. The derivative of the partonic cross section with respect to the rapidity cut-off $\zeta$ plays the same role of the CS evolution for the TMD (see Eq. (26)). Then, in analogy with the soft kernel ${\widetilde{K}}$, we define the rapidity-independent kernel that determines the CS-evolution for the partonic cross section as:

$$\begin{aligned}&\frac{\partial }{\partial \log {\sqrt{\zeta }}} \log {\left( \frac{d {\widehat{\sigma }}_j \left( \mu ,\,\lambda (T_c),\,\zeta \right) }{dz \, dT}\right) } \nonumber \\&\quad = \frac{1}{2} {\widehat{K}} \left( \alpha _S(\mu ),\, {\mu ^2}/{\lambda (T_c)^2}\right) . \end{aligned}$$

(103)

The kernel ${\widehat{K}}$ has an additive anomalous dimension:

$$\begin{aligned} \frac{\partial }{\partial \log {\mu }} \, {\widehat{K}}\left( \alpha _S(\mu ),\,{\mu ^2}/{\lambda (T_c)^2}\right) = \gamma _K\left( \alpha _S(\mu )\right) . \end{aligned}$$

(104)

This anomalous dimension is equal and opposite to $\gamma _K$, which is associated to the soft kernel ${\widetilde{K}}$ (see Eq. (11)). Finally, the partonic cross section shows an explicit dependence on the scale $\lambda (T_c)$ used to constrain the transverse momentum of the fragmenting parton. The corresponding evolution has no analogue in the TMDs. It is given by:

$$\begin{aligned}&\frac{\partial }{\partial \log {\lambda }}\, \log {\left( \frac{d {\widehat{\sigma }}_j \left( \mu ,\,\lambda (T_c),\,\zeta \right) }{dz \, dT}\right) }\nonumber \\&\quad = G\left( \alpha _S(\mu ),\, {\mu ^2}/{Q^2},\,{\zeta }/{\mu ^2},\, {\mu ^2}/{\lambda (T_c)^2}\right) . \end{aligned}$$

(105)

The $\lambda$-evolution kernel G is RG invariant. Furthermore, it obeys the following CS-evolution:

$$\begin{aligned}&\frac{\partial }{\partial \log {\sqrt{\zeta }}} G\left( \alpha _S(\mu ),\, {\mu ^2}/{Q^2},\,{\zeta }/{\mu ^2},\, {\mu ^2}/{\lambda (T_c)^2}\right) \nonumber \\&\quad = \frac{1}{2} \, \frac{\partial }{\partial \log {\lambda }} {\widehat{K}}\left( \alpha _S(\mu ),\,{\mu ^2}/{\lambda (T_c)^2}\right) . \end{aligned}$$

(106)

Finally, the solution to the evolution equations Eqs. (102), (103) and (105) gives:

$$\begin{aligned}&\frac{d {\widehat{\sigma }}_j (\mu ,\,\lambda (T_c),\,\zeta )}{dz \, dT} = \frac{d {\widehat{\sigma }}_j}{dz \, dT}\Big |_{\begin{array}{c} \text {ref.} \end{array}}\, \nonumber \\&\quad \times \hbox {exp} \Bigg \{ \int _\mu ^Q\,\frac{d\mu '}{\mu '}\, \gamma _D\left( \alpha _S(\mu '),\,{\zeta }/{(\mu ')^2}\right) \Bigg \} \nonumber \\&\quad \times \hbox {exp} \Bigg \{ \frac{1}{4}\, {\widehat{K}}\left( \alpha _S(Q),\,1\right) \, \log \frac{\zeta }{Q^2} \nonumber \\&\quad -\int _{\lambda (T_c)}^Q \,\frac{d\lambda '}{\lambda '}\, G\left( \alpha _S(Q),\, 1,\,{\zeta }/{Q^2},\, {Q^2}/{(\lambda ')^2}\right) \Bigg \}, \end{aligned}$$

(107)

where we have also used the RG-invariance of the kernel G. The label “ref.” stands for the energy scales at the reference values $\mu = Q$, $\zeta = Q^2$ and $\lambda (T_c) = Q$.

So far, we have considered $\lambda$ and $\zeta$ as independent. As a consequence, the evolution of the partonic cross section can be written in perfect analogy to the TMD evolution. Furthermore, this approach makes the RG-invariance of the final cross section explicit, see Eq. (102). However, the correct separation between hard and collinear momentum regions, represented by the partonic cross section and the TMD FFs respectively, can only be obtained by setting $\zeta = \lambda (T_c)^2$. This is due to the presence of an upper boundary for the transverse momentum of the fragmenting parton, which automatically reflects onto a lower limit for the rapidity of the produced particles. Therefore, the only cut-off left in the final cross section is the topology cut-off $T_c$. It enters in the final formula of Eq. (97) through the partonic cross section and through the rapidity cut-off of the TMD FFs. Its value has to be chosen according to the topology of the final state and, ultimately, it depends on the kinematics of the process.

Fig. 7

Amplitude squared for the LO partonic tensor, in the limit $T \rightarrow 1$

6.4 A simple example of the rapidity dilation mechanism

A crude example of the thrust-dependent cross section of $e^+e^-\rightarrow H \, X$ process can be obtained from a basic QCD approximation in the 2-jet limit, i.e. $T \rightarrow 1$. At lowest order, the subtraction mechanism is trivial and the subtracted, renormalized hard coefficient are easily computed from Eq. (101):

$$\begin{aligned} \left( {\widehat{W}}_f^{\mu \, \nu }\right) ^{[0]} (z,\,T) = \left( {\widehat{W}}_f^{\mu \,\nu }{}^{\mathrm{, unsub}} \right) ^{[0]} (z,\,T). \end{aligned}$$

(108)

In the 2-jet limit, the only Feynman diagram contributing to the l.h.s of the previous equation is given by Fig. 7. It is an exact 2-jet configuration, hence $T = 1$. As a consequence, the phase space integration is trivial and does not require the introduction of the topology cut-off $T_c$. The actual computation is easier for the projections (see Eqs. (71) and (72)):

$$\begin{aligned}&-g_{\mu \,\nu } \left( {\widehat{W}}_f^{\mu \,\nu }\right) ^{[0]} (z,\,T)\nonumber \\&\quad =(1-\epsilon ) \, \delta _f^{\;q} e_q^2 \; 2 N_C \, \delta (1-z) \, \delta (1-T) ; \end{aligned}$$

(109)

$$\begin{aligned}&\quad \frac{k_\mu k_\nu }{Q^2} \left( {\widehat{W}}_f^{\mu \,\nu }\right) ^{[0]} (z,\,T)= 0 . \end{aligned}$$

(110)

Notice that the gluon contribution is always suppressed in a 2-jet configuration. Then we can compute the lowest order subtracted, renormalized structure functions:

$$\begin{aligned} {\widehat{F}}^{[0]}_{1,\,f}(z,\,T)&= \delta _f^{\;q} e_q^2 \;N_C \, \delta (1-z) \, \delta (1-T) ; \end{aligned}$$

(111)

$$\begin{aligned} {\widehat{F}}^{[0]}_{2,\,f}(z,\,T)&= -\frac{2}{z} \; {\widehat{F}}^{[0]}_{1,\,f}(z,\,T) . \end{aligned}$$

(112)

Finally, by using Eqs. (91), (92) and (93), the LO subtracted and renormalized partonic cross section appearing in the final result of Eq. (97) is given by:

$$\begin{aligned} \frac{d {\widehat{\sigma }}_f^{[0]}}{dz\,dT}&= a_T \, \frac{4 \pi \alpha ^2}{3 Q^2} \, z \, {\widehat{F}}^{[0]}_{1,\,f}(z,\,T) \nonumber \\&= a_T \, \frac{4 \pi \alpha ^2}{3 Q^2} \, N_C \, \delta _f^{\;q} e_q^2 \; \delta (1-z) \, \delta (1-T), \end{aligned}$$

(113)

where the factor $a_T$ accounts for the limited acceptance in the polar angle $\theta$. In the following, the detected hadron will be considered spinless for simplicity. Hence, the Sivers-like contribution disappears and in the cross section will remain only the unpolarized TMD FF $D_1$. Its crudest estimate is the Leading Log (LL) approximation, given by:

$$\begin{aligned}&{\widetilde{D}}_{1\;j/H}^{LL}(z,\, b_{T}; \, Q, \,\zeta ) \nonumber \\&\quad = \frac{1}{z^2} \, d_j (z,\,\mu _b) \nonumber \\&\quad \quad \times \hbox { exp} \left\{ L_b \, g^{\mathrm{LL}}_1 \left( a_S (Q) L_b \right) \right. \nonumber \\&\quad \quad \left. + g^{\mathrm{LL}}_2 \left( a_S (Q) L_b , \, \log {\left( {\zeta }/{Q^2}\right) } \right) \right\} \nonumber \\&\quad \quad \times \left( M_{D_1}\right) _{j,\,H}(z,\,b_T) \nonumber \\&\quad \quad \times \hbox { exp}\left\{ -\frac{1}{4} \, g_K(b_T) \, \log {\left( z^2 \, \frac{\zeta }{M_H^2} \right) } \right\} , \end{aligned}$$

(114)

where $L_b = \log {({Q}/{\mu _b})}$ and the functions $g^{\mathrm{LL}}_1$, $g^{\mathrm{LL}}_2$ are given in Eq. (137). Notice that, since there is no need for a topological cut-off, the rapidity cut-off $\zeta$ is unconstrained in the LO, LL approximation. This is a consequence of the low degree of information encoded in the perturbative part of the TMD. Basically, all the constraints on the rapidity of the collinear particles are contained in the non-perturbative part of the TMD FFs. Therefore, any modification of $\zeta$ has to be traced back to a modification of the non-perturbative model $M_{D_1}$ that describes the fragmentation mechanism. The rapidity dilation transformation discussed in Sect. 3.2 allows to choose the rapidity cut-off consistently with the choice of the model. Therefore, the LO, LL cross section is written in terms of a generic $\zeta$. It is given by:

$$\begin{aligned}&\frac{ d \sigma _{2{\mathrm{-jet}}} ^{[0],\;LL}}{dz\, dP^2_T \,dT} \nonumber \\&\quad = \pi \sum _{j} \int _z^1 \frac{d {\widehat{z}}}{{\widehat{z}}} \, \frac{d {\widehat{\sigma }}_j^{[0]}}{d({z}/{{\widehat{z}}})\,dT} \, \nonumber \\&\qquad \times \int \frac{d^2 {{\varvec{b}}}_T}{(2 \pi )^2} \, e^{i \, \frac{{{{\varvec{P}}}}_{p,\,T}}{{\widehat{z}}} \cdot {{\varvec{b}}}_T} \, {\widetilde{D}}^{LL}_{j,\,H}({\widehat{z}}, \,b_T,\,Q,\,\zeta ) \,\nonumber \\&\qquad \times \left[ 1 + {\mathcal {O}}\left( \frac{M^2}{Q^2}\right) \right] \nonumber \\&\quad = a_T \, \frac{4 \pi ^2 \alpha ^2}{3 Q^2} \, N_C \, \delta (1-T) \, \sum _q \, e_q^2 \, \nonumber \\&\qquad \times \int \frac{d^2 {{\varvec{b}}}_T}{(2 \pi )^2} \, e^{i \, \frac{{{{\varvec{P}}}}_{p,\,T}}{z} \cdot {{\varvec{b}}}_T} \, \frac{1}{z^2} \, d_q (z,\,\mu _b) \nonumber \\&\qquad \times \hbox { exp} \left\{ L_b \, g^{\mathrm{LL}}_1 \left( a_S (Q) L_b \right) \right. \nonumber \\&\qquad \left. + g^{\mathrm{LL}}_2 \left( a_S (Q) L_b , \, \log {\left( {\zeta }/{Q^2}\right) } \right) \right\} \nonumber \\&\qquad \times \left( M_{D_1}\right) _{q,\,H}(z,\,b_T) \nonumber \\&\qquad \times \hbox { exp} \left\{ -\frac{1}{4} \, g_K(b_T) \, \log {\left( z^2 \, \frac{\zeta }{M_H^2} \right) } \right\} \, \nonumber \\&\qquad \times \left[ 1 +{\mathcal {O}}\left( \frac{M^2}{Q^2}\right) \right] \end{aligned}$$

(115)

Notice that this formula represents the simplest, non trivial approximation beyond the parton model picture. It holds valid to LO in the perturbative expansion and at $T=1$, hence it only has illustrative purposes. A reliable phenomenological analysis should not rely on Eq. (115), but rather on the full NLO expression, with the appropriate accuracy in the order of logarithms, which will soon be presented in Ref. [28]. In the following, we will give a prototypical application of this LO cross section formula to a small sub-sample of the BELLE data [1], which should only serve as an example of the rapidity dilation mechanism discussed in Sect. 3. The simplicity of its usage and the small number of free parameters involved in the fitting procedure are indeed the points of strength of Eq. (115).

Fig. 8

$e^+e^- \rightarrow H X$ LO cross section computed according to Eq. (115) using two different parameterizations for its non-perturbative part, and at two different values of the rapidity cut-off ($\zeta =Q_0^2$ solid black line, $\zeta =Q_0^2/4$ dashed green line). See text for more details

For our example, we will consider only the subset of the BELLE $e^+e^- \rightarrow HX$ cross sections, corresponding to $0.55<z<0.6,0.85<T<0.90$, in 20 $P_T$ bins ranging from 0.06 to 2.5 GeV. For the BELLE experiment $Q=10.58$ GeV. This data sub-sample is shown in Fig. 8. Statistical and systematical errors are added in quadrature. The analysis will be performed using the NNFF10 fragmentation function set at LO [31], and fixing the values of $b_{\mathrm{MIN}}$ and $b_{\mathrm{MAX}}$ as follows: $b_{\mathrm{MIN}}=C1/Q \sim 0.1$ $\hbox {GeV}^{-1}$ and $b_{\mathrm{MAX}}=1$ $\hbox {GeV}^{-1}$.

Let’s now suppose that, somewhere around the globe, Group A performs a phenomenological analysis of the above BELLE data subset using a power-law parameterization of the model in $P_T$-space which, in the $b_T$ space, corresponds to a Bessel-K function, normalized in such a way that it is 1 at $b_T=0$:

$$\begin{aligned} M_A (b_T,\,m,\,p) = \frac{2^{2-p}}{\varGamma (-1+p)} \, (b_T \, m)^{-1+p} \, K_{-1+p}(b_T \,m)\nonumber \\ \end{aligned}$$

(116)

where $K_{-1+p}$ is the modified Bessel function of the second kind. This model was successfully used in Ref. [32] to fit the $e^+e^-\rightarrow HX$ cross sections measured by the TASSO and MARKII Collaborations [33, 34]. Group A knows that the TMD cross section will become unphysical as $P_T$ grows larger, as it is only valid in the TMD region where $P_T<< P^+$ (here $P^+ = z\,Q/\sqrt{2}\sim 4.3$ GeV). Therefore they fix $P_{T,{\mathrm{MAX}}} = 1.8$ GeV. After this point the cross section will rapidly fall to zero and become negative. Having set their rapidity cut-off at $\zeta =Q^2$, Group A best fit returns $m_A=0.35$ and $p_A=3.00$ for their two free parameters. The resulting cross section is shown in Fig. 8 (red, solid line).

Fig. 9

The effect of a rapidity shift from $\zeta =Q^2$ to $\zeta =Q^2/4$ on the cross section extracted by Group A. As expected the cross section is not invariant under this transformation

On the other side of the planet Group B, totally unaware of the work of Group A, performs a fit on the same data sample, but they choose a Gaussian parameterization for the model of their cross section (clearly the perturbative part has the same functional form in both cases)

$$\begin{aligned} M_B(b_T,\,m, \,p) = e^ {-m \, {b_T}^{2} }\,. \end{aligned}$$

(117)

Here there is only one free parameter, m, as the power p has been fixed to 2 to obtain a Gaussian form. They set their rapidity cut-off to $\zeta =Q^2/4$ and decide to be conservative on their TMD-regime requirement, so they fix $P_{T,{\mathrm{MAX}}} = 1.3$ GeV. Their fit has only one free parameter, $m_B$, which the $\chi ^2$ minimization procedure sets to 0.12. The corresponding cross section is shown in Fig. 8, by the green dashed line.

Notice that, in principle, there is at least one more free parameter in both analyses, which is used to model the $g_K$ function, see Eq. (115). As explained in Sect. 3, however, $g_K$ does not depend on the rapidity cut-off, nor on the flavour j of the fragmenting quark. Therefore, it does not play any active role in a rapidity dilation and is not relevant in this example. We will therefore suppose it to be the same for Group A and B and parameterize it as $g_K = a \, b_T^2$ with $a=0.11$, fixed “a priori”.

Fig. 10

Left panel: non perturbative contribution to the LO cross section, corresponding to the same choice of rapidity cut-off. The solid black line represents the extraction of Group B, while the dashed green line is obtained from the extraction of group A by applying a rapidity dilation, i.e. through a transformation that brings $\zeta =Q_0^2$ to $\zeta =Q_0^2/4$ and compensates this variation by changing the value of the free parameters of the model $M_A$. See text for more details

As it is clearly shown in Fig. 8, the results obtained by Group A and B are consistent, within errors, as they fit the same data sample. Similarly, also the TMD fragmentation functions extracted by the two groups will be consistent at small $P_T$, where they carry a truly physical information about the transverse motion of the hadronizing parton. In $b_T$-space, the two TMDs are very similar at small $b_T$ but they may differ in their large $b_T$ behaviour, because of the different choices of models, $M_A(b_T)$ and $M_B(b_T)$.

It is at this point that Eq. (49) becomes crucial: in fact, it allows the two Groups to relate their independent extractions through a rapidity dilation. The two extractions will not correspond to a one-to-one relation in $b_T$-space, nevertheless rapidity dilations preserve the physical meaning of TMDs. First of all, Group A performs a rapidity shift on their extraction: as expected the cross section is not invariant for a variation of the rapidity cut-off. This is illustrated in Fig. 9. However, by applying a full rapidity dilation, i.e. transforming their TMD according to Eq. (42), Group A can match their results to those obtained by Group B, in the range of small $P_T$ where the TMD approximation holds valid and where information from the experimental data is able to constrain the model. In fact, according to Eq. (49), here we have:

$$\begin{aligned}&{\mathcal {FT}} \Big [ \sigma ^{\mathrm{NP}}_{\zeta ',\,M_B} (b_T) \Big ] \sim {\mathcal {FT}} \Big [ \sigma ^{\mathrm{NP}}_{\zeta ,\,M_A} (b_T)\, \exp \Big \{ -\frac{1}{2} \,\theta \, {\widetilde{K}}(b^\star _T) \Big \} \Big ] \nonumber \\&\quad \hbox { at small } P_T. \end{aligned}$$

(118)

Here $\theta =\log 2$.

This is shown in Fig. 10, where the black solid line represents the results of Group B for the non-perturbative contribution to the full cross section (left hand side of Eq. (118)), while the green line corresponds to the results of Group A for the analogous quantity after the application of a rapidity dilation (right hand side of Eq. (118)). Notice that $\sigma ^{NP}_A$ is related to $\sigma ^{NP}_B$ by a factor which is purely perturbative and therefore calculable and totally model independent, see Eq. (118).

Figure 11 shows the ratio of these two curves as a function of $P_T$, $R^{\mathrm{NP}}$. In an ideal world, where all extraction converged to the same model, $R^{\mathrm{NP}}$ would be 1 at all values of $P_T$ (dashed red line). However, in a realistic case $R^{\mathrm{NP}}$ is very close to 1 only at small $P_T$, as it should, and it starts deteriorating as $P_T$ grows larger. It is not by chance that it stays close to 1 up to $P_T \sim 1.3$, which corresponds to the value of $P_{T,{\mathrm{MAX}}}$ set by group B. After that point, the cross section starts to become unphysical and the ratio itself becomes meaningless. In Fig. 11 a thin gray vertical line marks $P_{T,{\mathrm{MAX}}} = 1.3$ GeV. Notice that the invariance under rapidity dilation is considerably powerful: it allows to preserve the physical part of the cross section, embodied by the TMD function at small $P_T$, even in a realistic situation in which a very limited range of $P_T$ is constrained by experimental information, while compensating for the variation of the rapidity cut-off in the perturbative part by a transformation of the non-perturbative model.

Fig. 11

The ratio between the non pertubative contributions to the cross sections calculated according to the extraction of Group B and the rapidity dilated extraction of Group A

7 Conclusions

In this paper we have extended the TMD factorization mechanism to processes belonging to different hadron classes. This is potentially a very powerful tool, as it allows us to exploit the same definition of TMD parton densities in different processes, which up to now could not be used in a simultaneous data analysis. With this extended definition of TMD, in particular, we have been able to apply the TMD formalism to the process of one hadron production from $e^+e^-$ scattering, belonging to the 1-h hadron class. Within this scheme, the TMD FFs extracted from a phenomenological analysis of the $P_T$ dependent $e^+e^-\rightarrow H X$ cross sections, can be related to the analogous TMD FFs as extracted in a 2-h class process, like SIDIS or $e^+e^-\rightarrow H_A \, H_BX$.

Clearly the extension of the factorization scheme comes to a price, a price that in this case turns out to be rather large and two-folded. First of all the soft factor, which is responsible for a (partial) breaking of universality, cannot be included in the definition of the TMD, as it is elegantly done in the standard TMD factorization through the “square-root” TMD definition. Freed by its soft contribution, TMD becomes truly universal and can be used in any class of processes. The soft factor, however, assumes a fundamental role as it becomes a pivotal ingredient of the factorized cross section, where the non-perturbative effects of soft physics are encoded in the soft model $M_S$. It will have to be extracted within its corresponding hadron class and should only be used within that class. The process $e^+e^-\rightarrow H X$ is a slightly exceptional case, as the soft factor here becomes unity, as shown in Sect. 6.

Having recovered a solid and truly universal definition of TMD, we can factorize cross sections as that of $e^+e^-\rightarrow H X$, where there is only one single TMD embodying the long-distance contributions. The all-order expression of this cross section has been obtained in Sect. 6 following a factorization scheme derived from the CSS factorization procedure. As a rude, first estimate of the final cross section we presented the result obtained to leading order (LO) and leading log (LL) accuracy. Here we had to face an additional problem: the arbitrariness in the choice of the rapidity cut-off reflects in the LO result, undermining its predictive power. To make the TMD independent of the choice of the rapidity cut-off, they have to be made invariant under a specific transformation, which we call “rapidity dilation”.

Such transformations regulate how the perturbative and non-perturbative contributions are balanced within the TMD itself. In fact, in physical observables the cut-off has to be taken very large ($y_1 \rightarrow \infty$) but in the TMDs alone (which are not physical observables) there is total arbitrariness in choosing its particular value. Rapidity dilations control this arbitrariness by acting both on the rapidity cut-off and on the model $M_C$. The larger $y_1$ the more $M_C$ is suppressed, and the TMD is, basically, only perturbative. Less extreme values of $y_1$, instead, will correspond to a more dominant non-perturbative contribution.

Separating perturbative and non-perturbative contributions is a highly non-trivial problem, which affects any phenomenological analysis. For example, ambiguities originate when we have to fix the value of $b_{MAX}$, which marks the critical value of the impact parameter at which non-perturbative contributions start becoming non negligible. TMDs are well defined within the approximation in which the partonic $k^+$ is very large while $k_T$ is small (i.e. collinear according to power counting), they should therefore correspond to partons with a very large rapidity and very small transverse momentum with respect to the jet axis. The rapidity cut-off $y_1$, formally, will have to be taken to infinity but, in practice, the specific size of $y_1$ will determine how far we stretch the perturbative content of the TMD and where the non-perturbative contribution will become dominant.

To clarify the practical relevance of rapidity dilation invariance, in the last Section of this paper we have presented a simple example to show how rapidity dilations can offer a tangible help in relating phenomenological analyses performed using different non-perturbative model assumptions and different values of the rapidity cut-off. This provides a solid basis for the interpretation of the results of independent TMD extractions.

Finally, we want to stress that the scheme we are proposing does not require a new start in the phenomenological analysis of all classes of hadronic processes. In fact, we can relate the TMDs obtained from data analyses based on the square root definition to the TMDs extracted using the factorization definition. This allows us to benefit all previous phenomenological analyses and extend them to 1-h class processes. This is indeed the strategy we are planning to pursue in the near future.

Appendices

Appendix A: Wilson Lines

A Wilson line (or a gauge link) is a path-ordered exponential operator defined by:

$$\begin{aligned} W_\gamma = P \left\{ \exp \left[ -i g_0 \, \int _0^1 ds \, {\dot{\gamma }}^\mu (s) A^a_{(0)\,\mu }(\gamma (s)) t_a \right] \right\} , \end{aligned}$$

(119)

where $\gamma$ is a generic path and P denotes the path ordering (i.e. when the exponential is expanded the fields corresponding to higher values of s are to be placed to the left). The coupling constant and the gluon field are bare quantities, as indicated by the label “0”. In the previous formula, $t_a$ are the generating matrices of the gauge group, in the appropriate representation. The Wilson lines guarantee that PDFs and FFs (in both collinear and TMD cases) are gauge invariant, by linking the quark to the anti-quark fields in the definition of the collinear factor (see Eq. (22) ). The Wilson line represents the (all order) propagation of a particle strongly boosted in some direction n. If this direction is a straight line the Wilson line depends only on the endpoints of the path and can be written in a compact way as:

$$\begin{aligned} W_n\left( x_2,\,x_1,\,n\right) = P \left\{ \exp \left[ -i g_0 \int _{x_1}^{x_2} d\lambda n^\mu A^a_{(0)\,\mu }(\lambda n) t_a \right] \right\} , \end{aligned}$$

(120)

If the strongly boosted particle is a quark, the associated Feynman rules are:

(121)

(122)

More details can be found for instance in Chapter 7 of Ref. [2].

Appendix B: Perturbative QCD and small $b_T$ region

Soft factors and collinear parts are well defined functions only over a rather small region in the transverse momentum space, according to power counting rules. The Fourier transform to the impact parameter space can be regarded as a kind of analytic continuation, because at fixed $b_T$ we can roughly access all transverse momenta with $k_T \le \frac{1}{b_T}$, even trespassing the original momentum region. In particular, the small $b_T$ region is associated with large transverse momenta, where perturbative QCD can be applied and a power expansion in $\alpha _s$ allows us to perform explicit calculations. This can be proved by a direct application of the factorization procedure to the small $b_T$ approximation of the Fourier transformed function. For the soft factor this can be found in Sect. 2, while for collinear parts we refer to Chapter 13 of Ref. [2] and to Ref. [35].

Despite the undeniable advantage provided by the possibility to perform explicit calculations in the small $b_T$ region, perturbative QCD is not enough to reproduce integrated quantities, which correspond to the Fourier transformed functions evaluated in $b_T = 0$. These can be recovered from the operator definitions, that obviously give a non-perturbative, all-order point of view. Therefore in $b_T = 0$, Eqs. (4) and (8) simply confirm that the integrated soft factor is the identity matrix, while Eq. (21) reproduces the integrated PDFs and FFs. The failure of perturbative QCD in $b_T = 0$ is due to the fact that the integral over ${{\varvec{k}}}_T$ is intrinsically ill defined, since it extends well beyond the physical momentum region where the TMDs and the soft factor are defined. As a consequence, new UV divergences arise and the counterterms in Eqs. (4) and (21) are not sufficient to cancel them. Therefore, the perturbative approach lead to definition of integrated functions as bare quantities and they need a renormalization in order to acquire physical meaning and reproduce the correct results. In the following, such renormalization procedure will be investigated for both the 2-h soft factor and the TMDs.

Fig. 12

Feynman graphs contributing to the small $b_T$ behavior of the 1 loop soft factor $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$. a Virtual diagrams (zero in dimensional regularization). b Real diagrams

1.1 Appendix B.1: Small $b_T$ behaviour of 2-h Soft Factor

The Feynman graphs in Fig. 12 show that in the small $b_T$ region the (renormalized) 2-h soft factor is given by:

$$\begin{aligned} \widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}(b_T,\,\mu ,\,y_1 - y_2)&= 1 - \frac{\alpha _S(\mu )}{4\pi } \, 8 C_F \,(y_1 - y_2) \nonumber \\&\quad \times \log {\frac{\mu \, b_T}{C_1}} + {\mathcal {O}} \left( \alpha _S^2, \; e^{-(y_1-y_2)} \right) , \end{aligned}$$

(123)

where $C_1 = 2 e^{-\gamma _E}$. The perturbative expansion of the previous equation should be valid at small $b_T$; however in this region $\log (\mu b_T/C_1)$ becomes large and sufficiently near to $b_T = 0$ it completely oversizes $\alpha _S$ so that the expansion becomes meaningless. Resummation in principle solves this problem. The soft kernel can be directly obtained from Eq. (123) by using the definition of Eq. (9) or of Eq. (10):

$$\begin{aligned} {\widetilde{K}}(b_T,\,\mu ) = - \frac{\alpha _S(\mu )}{4\pi } \, 16 C_F \, \log {\frac{\mu \, b_T}{C_1}} + {\mathcal {O}} \left( \alpha _S^2 \right) . \end{aligned}$$

(124)

This expressions implies that ${\widetilde{K}}(b_T,\,\mu )$ is large and positive as $b_T$ decreases. Therefore, the resummed soft factor of Eq. (18) vanishes in $b_T = 0$. An improvement can be reached by using a leading log estimate of ${\widetilde{K}}$ by using its evolution equation solution, Eq. (12). Actually, it is inappropriate to count the logs of a quantity, like the soft kernel, which is already the result of a resummation procedure. Despite this, we can apply the same recipe and set all terms to order $\alpha _S^0$ except $\gamma _K$, which has to be taken to 1 loop. This gives:

$$\begin{aligned} {\widetilde{K}}^{LL}(b_T,\,\mu ) = \frac{\gamma _K^{[1]}}{2 \beta _0} \, \log {\left( 1 - \frac{\alpha _S(\mu )}{4\pi } \log {\frac{\mu \, b_T}{C_1}} \right) }, \end{aligned}$$

(125)

that coincides with Eq. (124) in the limit $\alpha _S \rightarrow 0$. With this estimate, the divergence of ${\widetilde{K}}$ is much less severe but it is still there. calculations. An easy way to solve the problem and ensure that the perturbative QCD computation agrees with the operator definition prediction is to introduce a cut-off that prevents the soft transverse momentum to reach the UV region when it is integrated out. This can be implemented in $b_T$-space by introducing a new parameter $b_{\mathrm{MIN}} \ne 0$ that provides a minimum value for $b_T$. A modification of the $b^\star$ prescription, Eq. (6), is a simple way to insert this cut-off directly in the definition of the soft factor. For example, we can use the modified $b^\star$ prescription of Ref. [14]:

$$\begin{aligned} {{\varvec{b}}}_T^\star \left( b_c (b_T)\right) = {{\varvec{b}}}_T^\star \left( \sqrt{b_T^2+b^2_{\mathrm{MIN}}}\right) . \end{aligned}$$

(126)

Fig. 13

Leading Log (LL) Fourier Transformed 2-h Soft Factor in the small $b_T$ region. The introduction of $b_{\mathrm{MIN}}$ (blue solid line) allows to recover a full agreement with the operator definition of Eq. (8). If the regularization is not introduced (red dashed line), $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ vanishes in $b_T = 0$. Here $b_{\mathrm{MAX}} = 1$ $\hbox {GeV}^{-1}$

Then, the integrated soft factor is given by the unintegrated $\widetilde{{\mathbb {S}}}_{\mathrm{{2}-h}}$ evaluated in $b_T^\star \left( b_c (0)\right) = b_{\mathrm{MIN}}$. If $\mu$ can be considered a large energy scale (e.g. if it can be set equal to the hard energy scale Q of the process) then we can set $b_{\mathrm{MIN}} \propto {1}/{\mu }$. Consequently, all logs in Eqs. (123) and (124) are heavily suppressed and the soft kernel is zero at small $b_T$, while the soft factor is unity, see Fig. 13. Despite this kind of regularization has been devised for the 2-h soft factor, it applies equally well to the general soft factor ${\mathbb {S}}_{\mathrm{{N}-h}}$, where N can be any integer.

1.2 Appendix B.2: Small $b_T$ behaviour of TMDs

Formally, the integrated TMD is the Fourier transformed TMD computed at $b_T = 0$. In order to recover this result from Eq. (38) by applying perturbative QCD, the Fourier transformed TMD has to be renormalized, otherwise it would vanish in $b_T = 0$. This result can be proved by following the procedure described in Ref. [14]. First of all, thanks to the properties of the model $M_C$, Eq. (35), and of $g_K$, we can neglect all the non-perturbative content in Eq. (38) at small $b_T$. Furthermore, in this region we can approximate $b_T^\star$ with $b_T$. Then, it is a standard result that the $\alpha _S$ expansion of Wilson coefficients can be written as:

$$\begin{aligned}&\widetilde{{\mathcal {C}}}_f^{\,j}(\rho ,\, b_T; \, \mu , \,\zeta ) = \sum _{n=0}^{\infty } \left( \frac{\alpha _S(\mu )}{4\pi }\right) ^n \,\nonumber \\&\quad \times \sum _{k=0}^{2 n} \, \sum _{l=0}^{[{k}/{2}]} \widetilde{{\mathcal {C}}}_f^{\,j[n,\,k-l,\,l]}(\rho )\nonumber \\&\quad \times \left( \log {\frac{\mu }{\mu _b}} \right) ^{k-l} \, \left( \log {\frac{\zeta }{\mu _b^2}} \right) ^l, \end{aligned}$$

(127)

where [k/2] denotes the integer part of k. If the scales are fixed according to the standard choices of Eqs. (32) and (33), all the logs disappear and the only $b_T$ dependence in the Wilson coefficients is given by $\alpha _S (\mu _b)$. Since $\mu _b \propto {1}/{b_T}$, when $b_T \rightarrow 0$ the energy scale becomes very large and $\alpha _S$ can be considered a small parameter. For example, at 1 loop:

$$\begin{aligned} \frac{\alpha _S(\mu _b)}{4\pi } \overset{\hbox {low }b_T}{\sim } \frac{1}{2 \beta _0 \log {{\mu _b}/{\varLambda _{\mathrm{QCD}}}}}. \end{aligned}$$

(128)

Then, the Wilson coefficients evaluated at the scales $\mu _b$, $\zeta _b$ are well approximated at small $b_T$ by their lowest order term, which is simply a delta function:

$$\begin{aligned} \widetilde{{\mathcal {C}}}_j^{\;k}(\rho ,\, b_T; \, \mu _b, \,\zeta _b) = \delta _j^{\;k} \delta (1-\rho )+ {\mathcal {O}} \left( \frac{1}{\log {\frac{\mu _b}{\varLambda _{\mathrm{QCD}}}}} \right) \end{aligned}$$

(129)

On the other hand, ${\widetilde{K}}$, which is at exponent, allows for a different number of logarithms in front of each power of $\alpha _S$:

$$\begin{aligned} {\widetilde{K}}(b_T; \, \mu ) = \sum _{n=1} \left( \frac{\alpha _S(\mu )}{4\pi }\right) ^n \, \sum _{l=0}^{n} {\widetilde{K}}^{[n, l]} \, \left( \log {\frac{\mu }{\mu _b}} \right) ^{l}. \end{aligned}$$

(130)

As in the previous case, all the explicit dependence on $b_T$ vanishes if $\mu = \mu _b$ and ${\widetilde{K}}$ at small $b_T$ is well approximated by its lowest order term. However, since in this case the series starts from ${\mathcal {O}} \left( \alpha _S(\mu _b) \right)$, we can simply neglect this contribution. Finally, the anomalous dimensions have a simple expansion in $\alpha _S$:

$$\begin{aligned}&\gamma _{C}(\alpha _S(\mu ), 1) = \sum _{n=1} \left( \frac{\alpha _S(\mu )}{4\pi }\right) ^n \, \gamma _{C}^{[n]} , \end{aligned}$$

(131)

$$\begin{aligned}&\gamma _{K}(\alpha _S(\mu )) = \sum _{n=1} \left( \frac{\alpha _S(\mu )}{4\pi }\right) ^n \, \gamma _{K}^{[n]} . \end{aligned}$$

(132)

The easiest way to study the behavior of their contribution in Eq. (38) at small $b_T$ is to consider its derivative with respect to $\log {b_T}$. Since ${\partial }/{\partial \log {b_T}} = - {\partial }/{\partial \log {\mu _b}}$, we can compute with the help of Eq. (128):

$$\begin{aligned}&\frac{\partial }{\partial \log {b_T}} \int _{\mu _b}^{\mu } \frac{d \mu '}{\mu '} \, \sum _{n=1} \left( \frac{\alpha _S(\mu ')}{4\pi }\right) ^n\nonumber \\&\qquad \times \left[ \gamma _C^{[n]} - \frac{1}{4} \, \gamma _K^{[n])} \, \log {\frac{\zeta }{\mu '^2}}\right] \nonumber \\&\quad = \sum _{n=1} \left( \frac{\alpha _S(\mu _b)}{4\pi }\right) ^n \, \left[ \gamma _C^{[n])} - \frac{1}{4} \, \gamma _K^{[n]} \, \log {\frac{\zeta }{\mu _b^2}}\right] \nonumber \\&\quad = \frac{1}{4} \frac{\gamma _K^{[1]}}{\beta _0} + {\mathcal {O}} \left( \frac{1}{\log {b_T \, \varLambda _{\mathrm{QCD}}}} \right) . \end{aligned}$$

(133)

This behavior affects the whole TMD at small $b_T$, giving:

$$\begin{aligned} {\widetilde{C}}(\xi ,\, b_{T}; \, \mu , \,\zeta ) \overset{\hbox {low }b_T}{\sim } \left( b_T\right) ^{\frac{1}{4} \frac{\gamma _K^{[1]}}{\beta _0}} \times \hbox {log corrections } . \end{aligned}$$

(134)

From Eq. (124) $\gamma _K^{[1]} = 16 \, C_F$, then the TMDs goes to zero when $b_T \rightarrow 0$ with a power-law behavior.

This is also confirmed by a direct computation of the leading log (LL) estimate of Eq. (38). In this approximations, all the quantities are taken at order $\alpha _S^0$, except $\gamma _K$ which instead is computed at 1 loop. The result is:

$$\begin{aligned}&{\widetilde{C}}_{f,\,H}(\xi ,\, b_{T}; \, \mu , \,\zeta ) \overset{\hbox {low }b_T}{\sim } \left( \delta _f^{j} \otimes c_j (\mu _b) \right) (\xi ) \nonumber \\&\quad \times \hbox { exp} \left\{ L_b \, g^{LL}_1 \left( a_S (\mu ) L_b \right) + g^{LL}_2 \left( a_S (\mu ) L_b \right) \right\} \end{aligned}$$

(135)

where $L_b = \log {{\mu }/{\mu _b}}$, $a_S = {\alpha _S}/{4\pi }$ and:

$$\begin{aligned} g^{LL}_1 (x)&= \frac{\gamma _K^{[1]}}{4 \beta _0} + \frac{\gamma _K^{[1]}}{8 \, x \, \beta _0^2} \, \log {(1 - 2 \beta _0 x)}, \end{aligned}$$

(136)

$$\begin{aligned} g^{LL}_2 (x)&= \frac{1}{8 \beta _0} \, \gamma _K^{[1]} \, \log {\frac{\zeta }{\mu ^2}} \, \log {(1 - 2 \beta _0 x)}. \end{aligned}$$

(137)

Notice that the function $g^{LL}_2$ contributes to the LL estimate even if it typically appears at NLL. This is due to the presence of three scales instead of two. In fact, if $\sqrt{\zeta }$ equals either $\mu$ or $\mu _b$, only $g^{LL}_1$ contributes to LL. sufficient to cancel them.

As a consequence of the previous arguments, integrated TMDs are bare quantities when approached perturbatively. Formally:

$$\begin{aligned}&\int d^{D-2} {{\varvec{k}}}_T \, C_{j,\,H}(\xi ,\, k_{T}; \, \mu , \,\zeta ) = c_{j,\,H}^{(0)}(\xi ,\,\mu ) \nonumber \\&\quad = {\left\{ \begin{array}{ll} f_{j/H}^{(0)}(x,\,\mu ) &{}\hbox { initial state;} \\ z^{-2+2\epsilon } d_{H/j}^{(0)}(z,\,\mu ) &{}\hbox { final state.} \end{array}\right. } \end{aligned}$$

(138)

The bare integrated TMDs in the equation above acquire their dependence on $\mu$ through the renormalized fields used to compute them. Real bare quantities are defined through bare fields and are obtained by multiplying by $Z_2$ as in Eq. (25). Notice that integration makes the soft-collinear subtractions trivial, because the ${\mathbb {S}}_{\mathrm{{2}-h}}$ appearing in the factorization definition is unity when integrated over all soft transverse momentum. The required UV counterterm depends on the plus component of the momentum of the reference parton, i.e. on the collinear momentum fraction $\xi$. Hence, the renormalized quantities are not simple products of the bare quantities with the UV counterterm, like in Eq. (25), but rather convolutions

$$\begin{aligned} c_{j,\,H}(\xi ,\,\mu ) = \left( \left( Z_{\mathrm{int}}\right) _{j}^{\,k} (\alpha _S(\mu )) \otimes c_{k,\,H}^{(0)} \right) (\xi ), \end{aligned}$$

(139)

where now $c_{k,\,H}^{(0)}$ denotes a bare quantity computed with bare fields. With this definition, we can interpret the renormalized integrated TMDs as the usual PDFs and FFs used in collinear factorized cross sections as in Eq. (1).

The factorization procedure applied to the TMD at small $b_T$ does not give Eq. (36) directly. Instead, it expresses the final result as a convolution between a collinear part, represented by the unrenormalized integrated TMDs, and a hard factor ${\mathbb {H}}$ which has to be properly subtracted in order to cancel the double counting due to the overlapping between the hard and the collinear momentum region. This subtraction mechanism is completely anologous to that used in the definition of the subtracted collinear part in Eq. (22). Roughly speaking, the UV part of the bare integrated TMDs is (minus) $Z_{\mathrm{int}}$, then the subtracted hard part acquires the divergence induced by the counterterm. Despite this, we can still define a finite hard part by interpreting ${\mathbb {H}}^{\mathrm{sub}}$ as a bare quantity as well, with its renormalized finite counterpart represented by the Wilson Coefficients in the OPE. As a consequence, the required counterterm will be exactly $Z^{-1}_{\mathrm{int}}$. Then, a straightforward application of the convolution property shows that:

$$\begin{aligned}&{\widetilde{C}}_{j,\,H}(b_{T}; \, \mu , \,\zeta ) \overset{\hbox {low }b_T}{\sim } \;\;\; \left( {\mathbb {H}}^{\mathrm{sub}}\right) _j^{\,k} (b_{T}; \, \mu , \,\zeta ) \otimes c_{k,\,H}^{(0)}\nonumber \\&\quad = \left[ \left( {\mathbb {H}}^{\mathrm{sub}}\right) _j^{\,k} (b_{T}; \, \mu , \,\zeta ) \otimes \left( Z^{-1}_{\mathrm{int}}\right) _{k}^{\,l} (\alpha _S(\mu )) \right] \nonumber \\&\qquad \otimes \left[ \left( Z_{\mathrm{int}}\right) _{l}^{\,m} (\alpha _S(\mu )) \otimes c_{m,\,H}^{(0)} \right] \nonumber \\&\quad = \widetilde{{\mathcal {C}}}_j^{\,k} (b_{T}; \, \mu , \,\zeta ) \otimes c_{k,\,H} (\mu ) . \end{aligned}$$

(140)

Therefore, the functions $c_{j,\,H}$ appearing in the OPE are the renormalized integrated TMDs. Notice that the same procedure is used in the cross sections where the usual PDFs and FFs appear.

Different renormalizations of the integral over ${{\varvec{k}}}_T$ are allowed. A common procedure, for instance, is to introduce a cut-off as we did for the 2-h soft factor in B.1, by introducing a new parameter $b_{\mathrm{MIN}} \ne 0$ that provides a minimum value for $b_T$, for istance as in Eq. (126). Then, the integrated TMD is given by the unintegrated TMD evaluated in $b_T^\star \left( b_c (0)\right) = b_{\mathrm{MIN}}$:

$$\begin{aligned}&\int d^{D-2} {{\varvec{k}}}_T \, C_{f,\,H}(\xi ,\, k_{T}; \, \mu , \,\zeta ) = {\widetilde{C}}_{f,\,H}(\xi ,\, b_{\mathrm{MIN}}; \, \mu , \,\zeta )\nonumber \\&\quad \sim \left( \widetilde{{\mathcal {C}}}_j^{\,k} (b_{\mathrm{MIN}}; \, \mu , \,\zeta ) \otimes c_{k,\,H} (\mu ) \right) (\xi ), \end{aligned}$$

(141)

where in the last step we used the OPE expansion valid at small $b_T$. In general, this result does not coincide with $c_{f,\,H}(\xi ,\,\mu )$, but it will do if the Wilson Coefficients can be well approximated by their lowest order. If $\mu$ can be considered a large energy scale (e.g. if it can be set equal to the hard energy scale Q of the process) then we can set $b_{\mathrm{MIN}} \propto {1}/{\mu }$. Then all the logs inside the Wilson Coefficients are heavily suppressed and the lowest order approximation is reliable. Therefore, if $\mu$ is large enough, the cut-off approach gives the same result of the renormalization through the UV counterterm $Z_{\mathrm{int}}$. Thanks to $b_{\mathrm{MIN}}$, the subtraction mechanism implemented in the factorization procedure applied to the TMD at small $b_T$ is now applied to the collinear parts instead of the hard factor. Therefore, we do not have to worry about subtracting the hard part. However, the final result coincides with that of Eq. (36) because, trivially, ${\mathbb {H}}^{\mathrm{sub}} \otimes {\mathbb {C}}^{\mathrm{unsub}} = {\mathbb {H}}^{\mathrm{unsub}} \otimes {\mathbb {C}}^{\mathrm{sub}}$.

The integration over ${{\varvec{k}}}_T$ of the TMD, actually gives the area under the curve designed by the TMD in $k_T$-space. Even with the introduction of an explicit $b_{\mathrm{MIN}}$, the value of such integral is very small. Since in momentum space, at small $k_T$, the TMD is positive (e.g. Gaussian behavior), the small value of the integrand implies that the TMD has to change sign at a certain $k_T$. This is equivalent to say that the TMD loses its physical meaning when $k_T$ becomes too large. In fact, the power counting imposes $k_T \sim \lambda$, where $\lambda$ is some small IR energy scale.

Fig. 14

Momentum flow that determines the kinematics boundaries on ${\widehat{z}}$

Fig. 15

The LAB frame and the h-frame are both c.m. frames, but differ by a spatial rotation

Appendix C: Kinematics

As stressed in Sect. 1.1, kinematics play a crucial role in factorization, as it determines whether we need to apply a TMD or a collinear factorization scheme. The study of kinematics is strictly connected to the choice of the frame. In the case of $e^+e^-\rightarrow H\,X$, three four-vectors underlay the kinematical configuration:

• The momentum k of the fragmenting parton.

• The momentum P of the outgoing detected hadron H of mass M, $P^2 = M^2$.

• The momentum q of the highly virtual time-like photon that makes the partonic state. Its squared momentum gives the square of the center of mass energy $Q>> M$, $q^2 = Q^2$.

Clearly, the choice of the frame is completely arbitrary since the cross section will be Lorentz invariant. Three main frames are useful in deriving the final form of the factorized cross section: in this appendix we will provide a short description of all of them.

• 1. Hadron frame, labeled by h. This is the frame where the outgoing hadron H has no transverse components and it moves very fast along the (positive) $z_h$-direction:

  $$\begin{aligned} {{{\varvec{P}}}}_{T,\,h} = {{\varvec{0}}}_T. \end{aligned}$$

  (142)

  Furthermore, since H is strongly boosted in the plus direction its plus component is very large, of order $\sim Q$. As a consequence, its minus component has to be very small in order to satisfy the on-shell condition $P^2 = 2 P_h^+ \,P_h^- = M^2$. Therefore, in this frame, the full four-momentum P can be written as:

  $$\begin{aligned} P = \left( P_h^+,\,\dfrac{M^2}{2 P_h^+}, {{\varvec{0}}}_T \right) _h \sim Q \left( 1,\,\dfrac{M^2}{Q^2}, \,0 \right) . \end{aligned}$$

  (143)

  The fragmenting parton belongs by definition to the same collinear group of the outgoing hadron, hence it is almost collinear to it: it has a very large plus component, a low transverse momentum and an even lower minus component. It is almost on-shell, with a very low virtuality. Power counting (see Chapter 5 in Ref. [2]) allows us to quantify the sizes of these quantities by introducing a small infrared scale $\lambda<< Q$. Then $k^2 = \lambda ^2$, which means $k_h^+ \sim Q$, $k_h^- \sim {\lambda ^2}/{Q}$ and $k_{h,\,T} \sim \lambda$. Neglecting all the suppressed components, k and P become exactly collinear, i.e. $k \propto P$. This can be made explicit by setting:

  $$\begin{aligned} k_h^+ = \dfrac{1}{{\widehat{z}}} P_h^+, \end{aligned}$$

  (144)

  Therefore $P \sim {\widehat{z}} k$, and

  $$\begin{aligned} k = \left( \frac{ P_h^+}{{\widehat{z}}},\,k_h^-, {{\varvec{k}}}_{T,\,h} \right) _h \sim Q \left( 1,\,\dfrac{\lambda ^2}{Q^2}, \, \frac{\lambda }{Q} \right) . \end{aligned}$$

  (145)

  Since power counting rules are defined in the hadron frame, this is the most appropriate frame where to implement factorization. We can interpret ${\widehat{z}}$ as the collinear momentum fraction that the outgoing hadron takes off the fragmenting parton. Clearly ${\widehat{z}}$ has kinematics boundaries, due to the requirement that all the particles crossing the final state cut are physical, i.e. they have positive energy. With the help of Fig. 14 and by applying the power counting rules, we obtain the following constraints:

  • Positive energy for the final state of the jet

    $$\begin{aligned} (k - P)_h^0 \sim k_h^+ - P_h^+ = P_h^+ \left( \frac{1}{{\widehat{z}}} - 1 \right) \ge 0, \end{aligned}$$

    (146)

    which gives ${\widehat{z}} \le 1$.

  • Positive energy in the hard part of the process (given that $q_h^- >0$)

    $$\begin{aligned} (q - k)_h^0\ge & {} 0 \rightarrow q_h^+ - k_h^+ \nonumber \\= & {} \frac{Q}{\sqrt{2}} \left( 1 - \frac{\sqrt{2} P_h^+}{Q} \, \frac{1}{{\widehat{z}}}\right) \ge 0. \end{aligned}$$

    (147)

    The fractional energy z is defined as

    $$\begin{aligned} z = 2 \, \frac{P \cdot q}{Q^2} = 2 \frac{E_{\mathrm{CM}}}{Q} \sim \frac{\sqrt{2} P_h^+}{Q},\nonumber \\ \end{aligned}$$

    (148)

    where $E_{\mathrm{CM}}$ is the energy of the detected hadron in the center of mass frame. Then Eq. (147) gives the kinematics boundary: ${\widehat{z}} \ge z$, with $z \le 1$.

  The scaling of the components of the four-momentum q is obtained from the momentum conservation relation:

  $$\begin{aligned} q = k + \sum _{\alpha } k_\alpha , \end{aligned}$$

  (149)

  where $k_\alpha$ is the momentum of a generic real emission. As explained in Sect. 1.1, since the process $e^+e^-\rightarrow H X$ belongs to the ${1}$-hadron class, there is always at least one real emission (in this case the anti-quark leg that does not fragment) with a hard momentum, i.e. with all components very large, at least of order Q. As a consequence, the only component of k that survives in Eq. (149) is $k_h^+$, while all the others are strongly suppressed by the large momenta $k_\alpha$.

• 2. c.m. frame, labeled by $\gamma$. In this frame the spatial momentum of q is zero

  $$\begin{aligned} {{\varvec{q}}}_{\gamma } = {{\varvec{0}}}\,, \end{aligned}$$

  (150)

  which means

  $$\begin{aligned} q = \left( Q, \, {{\varvec{0}}} \right) _\gamma = \left( \frac{Q}{\sqrt{2}},\,\frac{Q}{\sqrt{2}}, {{\varvec{0}}}_{T} \right) _\gamma . \end{aligned}$$

  (151)

  Since rotations send null spatial vectors into null spatial vectors, the condition in Eq. (150) is defined modulo a rotation in space. Therefore, if we set the z-axis of this frame to be the direction of the outgoing hadron, we can identify the hadron frame with the c.m. frame and apply power counting and the whole factorization procedure directly in this frame. This is a big advantage, since usually the calculation of the hard part of the cross section is much easier in the c.m. frame but in general it does not coincide with the hadron frame, which on the other hand makes simpler the application of the factorization procedure.⁹ Then we can write the components of q in the h-frame as in Eq. (151). From Eq. (149) it follows that the total transverse momentum of the real emissions exactly cancels the contribution of ${{\varvec{k}}}_{T,\,h}$, hence $|\sum _\alpha {{\varvec{k}}}_{\alpha ,\,T,\,h}| \sim \lambda$.

  Notice that the LAB frame, in which the z-axis coincide with the beam axis, is a valid c.m. frame but it is not the hadron frame, as they differ by a spatial rotation, as shown in Fig. 15.

  The lepton pair is back-to-back in both the frames, but the direction of their spatial momenta is different.

• 3. Parton frame, labeled by p. As explained in Ref. [2], in order to properly define a fragmentation function we need a frame in which the fragmenting parton has zero transverse momentum. This is the parton frame, defined by requiring

  $$\begin{aligned} {{\varvec{k}}}_{T,\,p} = {{\varvec{0}}}_T. \end{aligned}$$

  (152)

  In principle we have two Lorentz transformations available that we can use to reach the parton frame from the hadron frame: a rotation of the (small) angle between the fragmenting parton and the outgoing hadron and a (light) transverse boost in the ${{\varvec{k}}}_{T,\,h}$ direction. By defining ${{\varvec{k}}} = {{{\varvec{k}}}_{T,\,h}}/{k_h^+}$, the angle of the rotation is $\alpha = -\sqrt{2} k$, while the parameter of the boost is ${\varvec{\beta }} =\sqrt{2} {{\varvec{k}}}$. The two choices give the same result:

  $$\begin{aligned} k&= \left( k_h^+, \,k_h^{-} - \dfrac{k_{h,\,T}^2}{2 k_h^+}, {{\varvec{0}}}_T \right) _p \nonumber \\&\quad + {\mathcal {O}} \left( \dfrac{\lambda ^2}{Q^2} \right) \left( 1,\,\dfrac{\lambda ^2}{Q^2}, \,1 \right) ; \end{aligned}$$

  (153)
  $$\begin{aligned} P&= \left( {\widehat{z}} \, k_h^+, \, \dfrac{M^2 + {\widehat{z}}^2 \,k_{h,\,T}^2}{2 {\widehat{z}} \, k_h^+}, \,-{\widehat{z}} \, {{\varvec{k}}}_{T,\,h} \right) _p \nonumber \\&\quad + {\mathcal {O}} \left( \dfrac{M^2,\,\lambda ^2}{Q^2} \right) \left( 1,\,\dfrac{M^2,\,\lambda ^2}{Q^2}, \,1 \right) . \end{aligned}$$

  (154)

  Notice that the plus components remain the same in the two frames (apart from power suppressed corrections). In this frame we can identify the $z_p$-axis as the axis of the experimental jet of hadrons in which H is detected. In fact, all the (almost) collinear particles in the jet have been generated by the same fragmenting parton and hence the sum of their spatial momenta has to be equal to ${{\varvec{k}}}_p = |{{\varvec{k}}}| \, {\widehat{z}}_p$, that lies on the (positive) z direction in this frame. Therefore, measuring $P_{p,\,T}$ gives the transverse momentum of the outgoing hadron with respect the jet axis. By definition, this axis coincides with the partonic thrust axis ${\widehat{n}}_p$, which is the direction that maximizes the partonic thrust $T_p$ defined as

  $$\begin{aligned} T_p = \dfrac{\sum _i |{{\varvec{k}}}_{h,\,i} \cdot {\widehat{n}}_p|}{\sum _i |{{\varvec{k}}}_{h,\,i}|}, \end{aligned}$$

  (155)

  where the sum runs over all the partons produced in the hard scattering, and ${{\varvec{k}}}_{h,\,i}$ is the spatial momentum in the c.m. frame of the ith outgoing hard parton. For example, in the case of two (back-to-back) partons $T_p = 1$ and ${\widehat{n}}_p$ is the axis of the parton pair, while for three partons $T_p = \hbox {max} \{x_1, \,x_2,\,x_3 \} \ge {2}/{3}$, with $x_i = 2 {|{{\varvec{k}}}_{h,\,i}|}/{Q}$, and ${\widehat{n}}_p$ is the direction of the i-th parton. Since $P_{p,\,T}$ is strictly connected to $k_{h,\,T}$, as shown in Eq. (154), its measurement offers powerful information on the partonic variables. However, the experimental measurement is on the transverse momentum of the outgoing hadron with respect to the hadron thrust axis ${\widehat{n}}_h$, which is the direction that maximizes the hadronic thrust $T_h$ defined as

  $$\begin{aligned} T_h = \dfrac{\sum _i |{{{\varvec{P}}}}_{\mathrm{CM},\,i} \cdot {\widehat{n}}_h|}{\sum _i |{{{\varvec{P}}}}_{\mathrm{CM},\,i}|}, \end{aligned}$$

  (156)

  where now the sum runs over all the detected particles in the center of mass frame (e.g. the LAB frame). Its value is close to its partonic counterpart, but they are not the same. As shown in Ref. [36], the observed distribution of hadronic thrust is related to the distribution with respect to the partonic thrust (which can be computed in perturbation theory) by a correlation function $C(T_h,\,T_p)$ that is sharply peaked around $T_h \sim T_p$. Therefore, roughly speaking, we can set $C(T_h,\,T_p) \sim \delta (T_h - T_p)$ and the direction which maximizes the hadronic thrust is approximately the same axis that maximizes the partonic thrust, i.e. ${\widehat{n}}_p \sim {\widehat{n}}_h$. The estimate of how much they differ can be made more quantitative in the simple case of a 2-jet configuration. In fact, in this case we have $T_p = 1$ and $T_h \sim 1 - {(M_1^2 + M_2^2)}/{Q^2}$ (see Ref. [36]), where $M_{1,\,2}$ is the invariant mass of the hadronic jets, hence $T_p - T_h \sim {\mathcal {O}}({M^2}/{Q^2})$. In this paper, we consider $P_{p,\,T}$ as a valid estimate of the transverse momentum of the outgoing hadron with respect to the hadronic thrust axis.