Effects of TMD evolution and partonic flavor on $e^+e^-$ annihilation into hadrons

We calculate the transverse momentum dependence in the production of two back-to-back hadrons in electron-positron annihilations at the medium/large energy scales of BES-III and BELLE experiments. We use the parameters of the transverse-momentum-dependent (TMD) fragmentation functions that were recently extracted from the semi-inclusive deep-inelastic-scattering multiplicities at low energy from HERMES. TMD evolution is applied according to different approaches and using different parameters for the nonperturbative part of the evolution kernel, thus exploring the sensitivity of our results to these different choices and to the flavor dependence of parton fragmentation functions. We discuss how experimental measurements could discriminate among the various scenarios.


Introduction
Transverse momentum dependent (TMD) parton distribution functions (PDFs) and fragmentation functions (FFs) depend on the longitudinal and transverse components of the momentum of partons with respect to the parent hadron momentum, as well as on their flavor and polarization state. The TMD PDFs and TMD FFs enlarge the amount of nonperturbative information carried by ordinary integrated PDFs and FFs because they open the window on explorations of the multi-dimensional structure of hadrons in momentum space in terms of their QCD elementary constituents. For example, in the last years several data for single-and double-spin asymmetries in semi-inclusive deep-inelastic scattering (SIDIS) have been accumulated and can be interpreted as originating from the effect of specific combinations of (polarized) TMD PDFs and TMD FFs (for a review, see Refs. [1][2][3][4]).
The TMD PDFs and TMD FFs can be defined only by a careful selection of physical observables that are sensitive to processes with two separate scales. In addition one needs to study the appropriate factorization theorems for these observables. For example, the appropriate factorization theorem for SIDIS holds true if the hard photon virtuality is accompanied by transverse momenta of the order of nucleon mass [5,6], which then are observed as a mismatch of collinear momenta. It is necessary that the definition of TMD functions includes all factorizable long-distance contributions to the physical cross section. These nonperturbative contributions, related to collinear gluon radiation, are summed into socalled gauge links that make the TMD functions color gauge invariant objects. Gauge links provide also the necessary phase to generate the above mentioned spin asymmetries [7][8][9]. Because initial-state and final-state gluon interactions are summed into different gauge links, the TMD functions may be process dependent, although parity and time-reversal invariance can simplify this non-universality to a simple proportionality factor [10,11]. To account for scale dependence, the TMD functions obey evolution equations that generalize the standard Renormalization Group Evolution (RGE) to a multi-scale regime in hard processes. TMD evolution equations have been derived for unpolarized TMD PDFs and TMD FFs [12,13], and for polarized ones only in a limited number of cases [14][15][16]. But despite these recent achievements, the phenomenological implementation of these effects is still under active debate [16][17][18][19]. From the experimental point of view, only few data sets are available with enough statistics that allows for a multidimensional analysis and a direct access to transverse momentum distributions [20,21]; in other cases, the studies were limited in the multidimensional coverage and by the restricted variety of targets and final-state hadrons [22][23][24][25][26].
In a preceding paper [27], the dependence of the intrinsic transverse-momentum distribution of both unpolarized TMD PDFs and TMD FFs upon the flavor and the longitudinal momentum of the parton involved was discussed using the recently published data from the Hermes collaboration [20] on multiplicities for pions and kaons produced in SIDIS off proton and deuteron targets. Although the flavor-independent fit of the data was not statistically excluded, a clear indication was found that different quark flavors produce different transverse-momentum distributions of final hadrons, especially when comparing different species of final hadrons. This feature corresponds quite naturally to the well known strong flavor dependence of integrated PDFs [28][29][30][31], and to indications from some models [32][33][34][35][36][37] and lattice calculations of TMD objects [38]. The SIDIS process is useful because it gives simultaneous access to TMD PDFs and TMD FFs. But the factorized cross section always involves a convolution of transverse momenta of the initial and the fragmenting partons: anticorrelation hinders a separate investigation of the two intrinsic distributions. Moreover, the Hermes data were collected at such a limited range in the hard scale that the statistical analysis of Ref. [27] was reasonably performed even without involving modifications due to evolution effects.
In this paper, we consider the semi-inclusive production of two back-to-back hadrons in electron-positron annihilations. In analogy with the SIDIS process, we define the multiplicities in e + e − annihilations as the differential number of back-to-back pairs of hadrons produced per corresponding single-hadron production. Then, we study their transverse momentum distribution at large values of the center-of-mass (cm) energy, starting from an input expression for TMD FFs taken from the analysis of Hermes SIDIS multiplicities at low energy performed in Ref. [27]. In this framework, we can extract clean and uncontaminated details on the transverse-momentum dependence of the unpolarized TMD FF, which is a fundamental ingredient of any spin asymmetry in SIDIS and, therefore, it affects the extraction also of polarized TMD distributions. Moreover, we can make realistic tests on the sensitivity to various implementations of TMD evolution available in the literature, since the hard scales involved in e + e − annihilations are much larger than the average values explored in SIDIS by Hermes , which is assumed as the starting reference scale.
An important difference between PDFs and FFs is the role of the gauge links arising mostly from resummation of gluons with collinear polarizations. T-odd effects for PDFs enter through the operator definitions of the PDFs after inclusion of appropriate gauge links having also transverse pieces. For FFs T-odd effects are contained in the hadronic states and as a consequence there are less universality-breaking effects for FFs [6,[39][40][41].
In Sec. 2, we outline the theoretical tools needed to work out the cross sections for annihilations in two hadrons and define the e + e − multiplicities. In Sec. 3, we introduce the QCD evolution of TMD FFs as the action of an evolution operator on input fragmentation functions, we describe some procedures to separate perturbative from nonperturbative domains of transverse momenta, and we provide some prescriptions to parametrize the nonperturbative contributions to the evolution kernel and the resummation of soft gluon radiation. In Sec. 4, we introduce the flavor decomposition of fragmentation processes. In Sec. 5, we make predictions for the spectrum in transverse momentum of e + e − multiplicities for production of two back-to-back hadrons, focusing on the sensitivity of results to the flavor of the fragmenting parton and to the different prescriptions for describing TMD evolution. Final comments and remarks are summarized in Sec. 6.
2 Multiplicities for e + e − annihilation into two hadrons e e + q P 1 P 1? f P 2 Figure 1. Kinematics for the e + e − annihilation leading to two back-to-back hadrons with momenta P 1 and P 2 .
We consider the process e + e − → h 1 h 2 X depicted in Fig. 1. An electron e − and a positron e + annihilate producing a vector boson with time-like momentum transfer q 2 ≡ Q 2 ≥ 0. A quark and an antiquark are then emitted, each one fragmenting into a residual jet containing a leading hadron that for simplicity we will consider unpolarized: the hadron h 1 with momentum and mass P 1 , M 1 , and the hadron h 2 with momentum and mass P 2 , M 2 . The two hadrons belong to two back-to-back jets, i.e. we have P 1 ·P 2 ≈ Q 2 . In the following, we will limit Q 2 values to a range where the vector boson can be safely identified with a virtual photon. Using the standard notations for the light-cone components of a 4-vector, we define the following invariants where is the electron momentum. The z 1 is the fraction of parton momentum carried by the hadron h 1 , and similarly for z 2 referred to the hadron h 2 . Covariantly, we can define the normalized time-like and space-like directionŝ Correspondingly, we can define the projector into the space orthogonal toẑ andt: The lepton momentum is then given by whereˆ µ ⊥ = µ ⊥ /| ⊥ | and µ ⊥ = g µν ⊥ ν . The g µν ⊥ projects onto the space orthogonal to q and P 2 . The projector onto the space orthogonal to P 1 and P 2 , namely in the hadron cm frame where P 1 and P 2 have no transverse components, is given by where the non-collinearity is defined as In the electron-positron cm frame of Fig. 1, we define the angle θ = arccos( ·ẑ/| |) whereẑ = −P 2 . It is related to the invariant y ≈ (1 + cos θ)/2. In analogy to the Trento conventions [42], we define the azimuthal angle so that P µ 1 = (0, |P 1⊥ | cos φ, |P 1⊥ | sin φ, 0) in this frame, and in any frame obtained from this one by a boost alongẑ. In general, the covariant definition is cos φ = −q T ·ˆ ⊥ /|q T |.
The cross section for the e + e − annihilation into back-to-back pairs of unpolarized hadrons can be written in a factorized formula at low transverse momenta [12,16,43,44]: where q T ≡ |q T | and A(y) = 1 2 − y + y 2 . The H is the hard annihilation part. The D q h 1 (z, b T ; ζ, µ) is the TMD FF in impact parameter space for an unpolarized quark with flavor q fragmenting into an unpolarized hadron h and carrying light-cone momentum fraction z and transverse momentum conjugated to b T [45]. Both H and D q h 1 are separated at the renormalization/factorization scale µ and evolve with it through renormalization group equations. The D q h 1 depends also on the scale ζ (with ζ 1 ζ 2 = Q 4 ) and evolves with it via a process-independent soft factor. The term Y (q 2 T /Q 2 ) ensures the matching with perturbative calculations at large transverse momenta.
In this paper, we will consider a kinematics where q 2 T /Q 2 ) term and corrections from higher twists of order M 2 /Q 2 or higher will be neglected. Moreover, the soft gluon radiation is here resummed into the TMD FF at the Next-to-Leading-Log level (NLL). It implies that the hard annihilation part is consistently calculated at leading order (LO) in α s , namely H(Q 2 , µ) ≈ 1. Equation (2.8) then simplifies to (2.9) In Sec. 5, we present our results for the q T spectrum of hadron pair multiplicities in e + e − annihilation. In strict analogy with the SIDIS definition [20], we construct the e + e − multiplicities as the differential number of back-to-back pairs of hadrons produced per corresponding single-hadron production after the e + e − annihilation. In terms of cross sections, we have where dσ h 1 h 2 is the differential cross section of Eq. (2.9). The dσ h 1 describes the production of a single hadron h 1 from the e + e − annihilation and it is obtained from the previous cross section by summing over all hadrons produced in one emisphere [43]:

TMD evolution of fragmentation functions
In the following, we describe in more detail the dependence of the fragmentation functions D q h 1 of Eq. (2.9) upon the renormalization/factorization scale µ and the scale ζ. Different scenarios are possible according to the choice of the initial starting value for the factorization scale, and of the low-energy model describing the nonperturbative part of the evolution kernel. We first describe the structure of the input D q h 1 at the starting scale.

Input fragmentation functions at the starting scale
We consider the unpolarized TMD FF extracted by fitting the hadron multiplicities in SIDIS data at low energy from Hermes [20]. The assumed functional form displays a transverse-momentum dependent part which is described in impact parameter space by the following fixed-scale flavor-dependent Gaussian ansatz 1 : where P 2 ⊥ a h (z) with a = q,q, is the flavor-and z-dependent Gaussian width at some starting scale Q 2 0 [27,46,47]. The choice of having separate Gaussian functions for different flavors is motivated by the significant differences displayed by the Hermes data between pion and kaon final-state hadrons [20]. The factorized collinear dependent part d a h 1 (z; Q 2 ) is described by using the DSS parametrization of Ref. [48].
Following Refs. [49,50], a possible energy dependence of the Gaussian distribution was taken into account introducing the logarithmic term with g 2 a free parameter. Choosing Q 2 0 = 1 GeV 2 , it was soon realized that the best-fit value for g 2 was compatible with zero. As a matter of fact, the Q 2 range spanned by Hermes is small and the obtained experimental data for multiplicities are not sensitive to evolution effects. For this reason, the fit was performed by using Eq. (3.1) at a scale fixed to the experimental average value, namely Q 2 = Q 2 0 = 2.4 GeV 2 . With this choice, the possible energy dependence of Eq. (3.2) is automatically eliminated.
In summary, the input to our studies on the evolution of D a h 1 with the scales µ and ζ is referred to the expression in Eq. (3.1) to be considered at the starting scale Q 2 0 = 2.4 GeV 2 . However, depending on the choice of the initial value of the factorization scale this identification is not always straightforward, as will be explained in the following sections.

The µ b prescription
As shown in Eq. (2.9), the TMD FFs generally depend on the factorization scale µ and on the scale ζ, that for convenience we name the rapidity scale. The TMD FFs satisfy evolution equations with respect to both of them [12,13]. The evolution with respect to µ is determined by standard RGE equations, whereas the evolution in ζ is determined by a process-independent soft factor [12,13].
The functional form of TMD FFs at small b T can be calculated in perturbative QCD. Conversely, the nonperturbative part at large b T must be constrained by fitting experimental data. At the medium/large energies of the Bes-III and Belle experiments, the perturbative tail of TMD FFs needs to be taken into account. Using the technique of Operator Product Expansion (OPE), it can be represented as a convolution of (perturbatively calculable) Wilson coefficients C with the (nonperturbative) collinear fragmentation functions d 1 [12,13]: The convolution is defined as The dependence of the coefficients upon both factorization and rapidity scales can be represented in a factorized form: where µ b is defined as and γ E is the Euler constant. The K function in Eq. (3.5) 2 arises from the processindependent soft factor that is necessary to proof the factorization theorem leading to the definition of the TMD FFs; it drives the evolution of TMD FFs in the ζ variable. The convolution in Eq. (3.4) is only valid for small b T , namely b T 1/Λ QCD . Moreover, the expression of the C coefficients consists in a power series in α s ln (µ 2 /µ 2 b ) (including also double logarithms of the same argument). The OPE is valid only when the logarithms do not diverge; this is accomplished, e.g., by choosing µ = µ b or q T , so that the series converges. Accordingly, if we choose µ = µ b we can write the TMD FF as The evolution of this fragmentation function from µ b to another value of µ (e.g., µ = Q) is driven by RGE equations. Instead, the evolution from an initial rapidity scale ζ i to ζ is controlled by the K function. The final expression of the TMD FF at the scales µ = Q and ζ is where the anomalous dimension γ F F reads and Γ cusp and γ V are also power series in α s in the MS scheme [16].
The above procedure is valid up to a maximum value of b T , that we name b max , beyond which we do not trust the perturbative calculation. Hence, it is convenient to reconsider the OPE by introducing the new variableb T that freezes at b max when b T becomes large: (3.10) For b T b max , the evolution in ζ is controlled by the function K(b T ; µb), where The nonperturbative part at large b T is defined as what is left over [51]: By adding the intrinsic transverse distribution at the starting scale (see Eq. (3.1)), Eq. (3.8) becomes . (3.13) If we insert ζ i = µ 2 b and ζ = µ 2 = Q 2 , the above equation reduces to (3.14) Hence, the net effect of evolution can be represented as the action of an evolution operator R on the input TMD FF evaluated at the scale µb, which is running withb T . This peculiar feature grants that there is a smooth matching between the perturbative domain at small b T and the nonperturbative domain at large b T . It is interesting to remark that from Eqs. (3.13) and (3.14) we deduce that modelling the nonperturbative part affects the whole b T spectrum, not only the large b T region.
In this paper, we resum the soft gluon radiation up to NLL contributions in ln (µ/µ b ), which corresponds to include terms linear in α s in the perturbative expansion of K and γ V , and quadratic in the expansion of Γ cusp [16]: are the usual Casimir operators for the gluon and fermion representations of the color group SU(N c ) with N c colors, and T F = n f /2 with n f the number of active quark flavors. Consistently, the coefficients C are computed at LO in α s , namely they reduce to δ functions such that Eq. (3.14) simplifies to The definition of g np (b T ) in Eq. (3.12) obviously implies that this function depends on b max , i.e. on the value of the impact parameter that sets the separation between the perturbative and nonperturbative regimes. Indeed, by perturbatively expanding K(b T ; µ b ) at lowest order we have [51] For b T b max , this expression recovers the quadratic parametrization 1 2 g 2 b 2 T adopted in the fits of Refs. [50] and [52], and it suggests that the parameter g 2 is not free but anticorrelated to b max , and proportional to b 2 max through a perturbatively calculable coefficient. The g np function accounts for the radiation of soft gluons emitted from a parton. A small (large) value of b max implies that the QCD perturbative description is valid up to relatively small (large) b T values. Consequently, the amount of soft gluons emission is larger (smaller) and we expect a large (small) value for g 2 . More generally, this anticorrelation is motivated by the fact that both the exact function K(b T ; µ b ) and the TMD FF itself must not depend on the arbitrary choice of b max . So, b max should not be regarded as a free parameter to be fitted to data, but it should be considered as an arbitrary scale that separates perturbative from nonperturbative regimes: changing b max implies a rearrangement of all terms in Eq. (3.13) such that the TMD FF does not change [51].
For the purpose of this work, we will consider anticorrelated pairs of values for {b max , g 2 }, inspired to the values adopted in Refs. [50] and [52]. We will also explore different expressions for each one of theb T and g np functions. Forb T , our first choice is the socalled "b-star" prescription [12,50] The second choice is based on the exponential function that is steeper and it approaches the asymptotic constant b max more quickly. For g np , we choose a linear function of b 2 T similarly to Refs. [49,50,52] (see also Eq. (3.2)): The second choice is suggested by Eq. (3.17): This expression was considered also in Ref. [53], and it reduces to Eq. (3.20) for small b T . In principle, we have four different combinations of prescriptions: However, after some preliminary exploration we realized that some of them were producing redundant results. Therefore, they have been neglected. In summary, the transverse-momentum spectrum of the multiplicities in Eq. (2.10) will be analyzed by varying the anticorrelated pair of parameters {b max , g 2 }, and by considering only the two combinations {b * T , g lin np } and {b † T , g log np }. Finally, we remark that if we choose Q = µb in Eq. (3.16), i.e. if we switch off evolution effects, we should recover the Gaussian model expression of Eq. (3.1) for the TMD FF at the initial scale Q 0 . Formally, this is not the case because in the second line the collinear d 1 is evaluated at µb and the term ( However, the Gaussian model of Eq. (3.1) is deduced by fitting the Hermes SIDIS data, whose kinematics overlaps the domain of very at the scale Q 0 and at very large b T values, or equivalently for very small parton transverse momenta.

The fixed-scale prescription
In Eq. (3.14), we have expressed the evolved TMD FF at a scale Q as the result of an evolution operator R acting on the same TMD FF evaluated at the scale µb running with b T . Alternatively, we can fix the initial scale at the value With this choice, it is not possible to apply the OPE for calculating a perturbative tail to which the TMD FF should match at low b T , as it was done in Eq. (3.3): we need a model input over the whole b T spectrum. In our case, it is now very easy to identify the input TMD FF at the starting scale Q i with the Gaussian parametrization of Eq.
The contribution from the g np term in the input distribution does not appear because of the choice of the starting scale ζ i = Q 2 i = Q 2 0 . The choice µ i = Q i of identifying the starting factorization scale with a fixed scale for the whole b T spectrum has important consequences also on the function K. From Eq. (3.15), we can expand K in powers of ln (µ/µ b ): if µ i = µ b , the series may not converge. One possible workaround is to apply the resummation technique to the K function itself [13]. Here, we will discuss two different prescriptions: computing K from Ref. [12] at a fixed order in α s ; or dressing K by resumming large logarithms of the kind ln (µ/µ b ) [13]. In the first case, K is expanded in powers of α s ; in the second case, the expansion is in α s ln (µ/µ b ). If µ i = µ b , the two expansions are the same.
We will refer to the first choice as the "fixed-scale" prescription. Contrary to the prescription described in the previous section, there is no need to define an arbitrary scale b max to separate perturbative from nonperturbative regimes. However, the function K is evolved from µb to Q i through its anomalous dimension: where g np (b T ) can get either the expression in Eq. (3.20) or in Eq. (3.21). The perturbative contributions are calculated at NLL as in Eq. (3.15), according to which we have The second choice is connected to the results of Ref. [13], because we resum all large logarithms of the kind ln (µ/µ b ) in the perturbative part as where D R is the resummed contribution computed in Ref. [13], and b T,c is the convergence radius of the perturbative expression. Apart from the resummation of logarithms, the main difference with Eq. (3.24) is the presence of the θ functions: nob T prescription is used to connect the perturbative and nonperturbative domains. And the nonperturbative contribution acts differently: while g np in Eq. (3.24) applies to the whole b T spectrum, in Eq.(3.25) it does only for b T > b T,c . Hence, we use the notationḡ np to account for this difference. For example, the K function must be at least continuous at b T = b T,c . We can match this constraint by defining the nonperturbative contribution at b T > b T,c as where g np can be again either the g lin np prescription of Eq. (3.20) or the g log np prescription of Eq. (3.21).
For Q i Q, the perturbative component D R in Eq. (3.25) diverges for b T < b T,c . Hence, its contribution to the evolution of the fragmentation function becomes negligible, being of the kind (Q/Q i ) −D R . Since K is a smooth function in b T , also the contribution of the nonperturbative partḡ np for b T > b T,c becomes numerically negligible [13]. However, this result cannot be generalized to any value of Q. Since we will make explorative calculations also at the Bes-III scale Q = √ 14.6 GeV which cannot be considered to be much larger than the initial scale Q 0 = √ 2.4 GeV of our input TMD FF, we will consider only the "fixed-scale" prescription of Eq. (3.24).

Summary of evolution kernels
In summary, we consider two possible ways of evolving the TMD FF, according to the choice of the initial factorization scale µ i . It is understood that all formulae are computed at the NLL level of accuracy, according to Eq. (3.15).
B The "fixed-scale" prescription: then where γ F F , µb, g np are defined in the same equations as above, while K is given in Eq. (3.24).

Flavor dependence of fragmentation functions
The flavor sum in Eq. (2.9) can be made explicit and further simplified using the symmetry upon charge-conjugation transformations: At the starting scale Q 0 , we distinguish the favored fragmentation where the fragmenting parton is in the valence content of the final hadron h. All the other channels are classified as unfavored fragmentation and are characterized by the fact that the detected hadron is produced by exciting more than one qq pair from the vacuum. If the final hadron is a kaon, we further distinguish a favored fragmentation initiated by an up quark/antiquark from the one initiated by a strange quark/antiquark. We limit the sum to three flavors u, d, s, and the corresponding antiquark partners.

Favored and unfavored fragmentation to different hadron species
For the final hadron pair being (h 1 , h 2 ) = (π + , π − ), the flavor sum in Eq. (2.9) becomes and Using the charge-conjugation symmetry of Eq. (4.1), it is simple to prove that the result for (h 1 , h 2 ) = (π − , π + ) is identical to the above one in Eq. (4.2).
If the final pions have the same charge, Again, because of charge-conjugation symmetry we get the same result for (h 1 , h 2 ) = (π − , π − ). If the final hadron pair is (h 1 , h 2 ) = (K + , K − ), the flavor sum becomes where and Charge-conjugation symmetry grants the same result for (h 1 , h 2 ) = (K − , K + ). where As before, we get the same result for (h 1 , where and charge-conjugation symmetry applied to Eq. (4.11) gives 14) and grants that the same result in Eq. (4.12) holds also for (h 1 , h 2 ) = (π − , K + ).

Flavor dependent Gaussian ansatz
The starting input to our analysis are the TMD FFs extracted by fitting the hadron multiplicities in SIDIS data from Hermes at Q 2 0 = 2.4 GeV 2 [27]. The assumed functional form displays a transverse-momentum dependent part which is described in impact parameter space by the following flavor-dependent Gaussian ansatz: The cross section of Eq. (2.9) (and, in turn, the multiplicity in Eq. (2.10)) is then a sum of Gaussians, and thus no longer a simple Gaussian. The width of the Gaussian depends also on the fractional momentum z, as done in several model calculations or phenomenological extractions [36,[54][55][56][57][58]. The chosen functional form is [27] where β, δ, γ, are fitting parameters and P 2 ⊥ q h ≡ P 2 ⊥ q h (ẑ), withẑ = 0.5.
Isospin and charge-conjugation symmetries suggest four different Gaussian shapes [27]:  Therefore, we have three independent favored fragmentations, and two independent unfavored fragmentations: The remaining favored channeld → π + is then given by The h = π − and h = K − channels can be deduced from the above ones using the chargeconjugation symmetry of Eq. (4.1). By inserting the Gaussian ansatz with the above assumptions in the expressions of Sec. 4.1, we get

Predictions for TMD multiplicities
In this section, we present our results as normalized multiplicities for the hadron pair (h 1 , h 2 ), where M h 1 h 2 (z 1 , z 2 , q 2 T , y) is defined in Eq. (2.10). In such way, we are able to directly compare the genuine trend in q 2 T for each different case. If not explicitly specified, we choose y = 0.2. For selected values of {z 1 , z 2 }, the results are displayed as a function of P 2 1⊥ = z 2 1 q 2 T . Hence, the useful range in P 2 1⊥ depends on z 1 in order to fulfill the condition q 2 T Q 2 . The range obviously depends also on the choice of the hard scale; we consider Q 2 = 100 GeV 2 , as in the Belle experiment, and Q 2 = 14.6 GeV 2 , as in the Bes-III one. For each specific case, the results are displayed as uncertainty bands: they represent the 68% of the envelope of 200 different values for the intrinsic parameters in Eqs. (4.16)-(4.20) for the D 1 (z, b T ; Q 2 0 ) at the starting scale Q 2 0 , obtained by rejecting the largest and lowest 16% of them. The 200 values are obtained by fitting 200 replicas of SIDIS multiplicities measured by the Hermes collaboration [20]. If the 200 values for each parameter were distributed as a Gaussian, the 68% band would correspond to the usual 1σ confidence interval (for more details, see Ref. [27]).
The results are organized as follows. In Sec. 5.1, we show the sensitivity of the normalized multiplicity to different values of the evolution parameters {b max , g 2 } described in Sec. 3.2 for a final hadron pair (h 1 h 2 ) = (π + π − ). In Sec. 5.2, we compare normalized multiplicities for the two different evolution schemes described in Secs. 3.2 and 3.3. In Sec. 5.3, we discuss the capability of discriminating among the various prescriptions illustrated in Sec. 3.2 for the nonperturbative evolution effects. In Sec. 5.4, we concentrate on the sensitivity of the normalized multiplicities upon varying the fractional energy z of final hadrons. In Sec. 5.5, we show how the results get modified when lowering Q 2 from the Belle scale to the Bes-III scale. Finally, in Sec. 5.6 we discuss the sensitivity of ratios of normalized multiplicities for different final states to the flavor structure of the intrinsic transverse-momentum-dependent part of the input TMD FF at the starting scale of evolution.

Sensitivity to nonperturbative evolution parameters
As already remarked in Sec. 3.2, for a specific evolution scheme the nonperturbative part of the TMD evolution depends on the choice of a prescription for describing the transition from perturbative to nonperturbative regimes, which in turn depends on the two parameters b max and g 2 . In this section, we explore the sensitivity of our predictions to different values of the pair {b max , g 2 }. We adopt as limiting cases the choices {b max = 1.5, g 2 = 0.18} and {b max = 0.5, g 2 = 0.68}, that were deduced in Refs. [52] and [50], respectively, by fitting the transverse-momentum distribution of lepton pairs produced in Drell-Yan processes. If not explicitly specified, the first choice is described by uncertainty bands with dot-dashed borders while the second choice is linked to bands with solid borders. As explained in Sec. 3.2, the two parameters are anticorrelated. In the following, we show results also for the interpolating choice {b max = 1, g 2 = 0.43}. The corresponding results are displayed as uncertainty bands with dashed borders.
In Fig. 2, the normalized multiplicity is shown as a function of P 2 1⊥ = z 2 1 q 2 T ≡ (0.5) 2 q 2 T at the Belle scale Q 2 = 100 GeV 2 for the "µ b scale" evolution scheme and with the {b * T , g lin np } prescription for the transition to the nonperturbative regime, as explained in Sec. 3.2. The explored range in P 2 1⊥ is such that for z 1 = 0.5 the maximum q 2 T satisfies the condition q 2 of 7%. We fix it by propagating to the normalized multiplicity the typical experimental error of 3% for single-hadron production data in e + e − annihilations at Q 2 = 100 GeV 2 and z = 0.5, from which the collinear d q 1 (z; Q 2 ) are extracted [48]. This experimental error of 7% seems small enough to discriminate among predictions produced with different choices of {b max , g 2 }.
Two additional light-gray bands are shown, which are partially overlapped (dot-dashed borders) or completely overlapped (dashed borders) to the band with solid borders corresponding to the choice {b max = 0.5, g 2 = 0.68}. These bands reproduce the outcome of calculations performed in the same conditions but for different (arbitrary) choices of the scale µ b . If the band with solid borders corresponds to calculations with the choice of Eq. (3.6) for µ b , then the light-gray band with dot-dashed borders corresponds to the choice µ b /2, and the one with dashed borders to 2µ b . The almost complete overlap of these results shows that for the selected observable, the normalized multiplicity, the theoretical uncertainty in determining the matching scale µ b (that describes the transition from perturbative to nonperturbative regimes) is negligible with respect to the sensitivity to the parameters describing the nonperturbative effects in the evolution.

Sensitivity to evolution schemes
In this section, we explore the sensitivity of our normalized multiplicity to the choice of the evolution scheme. In Sec. 3, we described two different schemes, the "µ b scale" and the "fixed scale". They differ mainly in the fact that in the latter the whole distribution in impact parameter space b T of the TMD FF D q 1 at beginning of evolution is computed at a fixed scale Q 0 , namely there is no impact parameter that describes the transition from low (perturbative) b T to high (nonperturbative) b T . Actually, one would expect that for small values of g 2 and corresponding not too large values of b max (i.e., where the perturbative description of the evolution of the b T distribution is still applicable and gives the predominant contribution) the predictions from the different schemes should tend to a common result, determined mainly by a fully perturbative calculation. However, the complexity of the evolution kernels, described in Secs. 3.2 and 3.3, indicates that this is too a naïve expectation.
in the same conditions and with the same notation as in Fig. 2, but for the "fixed scale" evolution scheme. The additional light-gray bands with dot-dashed and solid borders are the result related to the "µ b scale" evolution scheme for {b max = 1.5, g 2 = 0.18} and {b max = 0.5, g 2 = 0.68}, respectively.
It is evident that for the maximum (minimum) b max (g 2 ) the band with dot-dashed borders in the "fixed scale" scheme is not similar to the light-gray band with dot-dashed borders in the "µ b scale" scheme. Actually, all the results in the "fixed scale" scheme show a much larger distribution in P 2 1⊥ , somewhat pointing to stronger evolution effects of perturbative origin that seem to be absent in the "µ b scale" scheme (where the scale choice minimizes the effect of large logarithms in the perturbative coefficients). It is important to notice that there is a significant overlap between the band with dot-dashed borders in the "fixed scale" scheme and the light-gray band with solid borders in the "µ b scale" scheme.
Apparently, the normalized multiplicity seems not to be enough sensitive to discriminate among different evolution schemes, since two different choices of them can produce similar results with different evolution parameters {b max , g 2 }. However, this result is observed at a specific value of fractional energies of the final hadrons, namely z 1 = z 2 = 0.5.  In Fig. 4, we show the P 2 1⊥ distribution of normalized multiplicities calculated in the same conditions, notation and conventions as in the previous figure, but at z 1 = 0.3 and z 2 = 0.5. The band with dot-dashed borders in the "fixed scale" scheme can now be easily separated from the light-gray band with solid borders in the "µ b scale" scheme if the indicated hypothetical experimental error is around 7%. Therefore, only when combining the study of both the z and P 2 1⊥ dependencies in the normalized multiplicity we may be able to discriminate among different TMD evolution schemes.

Sensitivity to prescriptions for the transition to nonperturbative transverse momenta
We now focus on exploring the possibility of discriminating among different prescriptions that describe the functional dependence in b T of the nonperturbative Sudakov evolution factor (see Eqs. (3.20) and (3.21)) or the transition from the perturbative low−b T domain to the nonperturbative high−b T one (see Eqs. (3.18) and (3.19)). In Fig. 5, the normalized multiplicity of Eq. (5.2) is shown as a function of P 2 1⊥ = z 2 1 q 2 T ≡ (0.5) 2 q 2 T at the Belle scale Q 2 = 100 GeV 2 with the {b † T , g log np } prescription. Again, as in Fig. 3 there are two groups of uncertainty bands. The former one displays the results for the "fixed scale" evolution scheme in the standard notation, i.e. for {b max = 1.5, g 2 = 0.18} (dot-dashed borders), {b max = 1, g 2 = 0.43} (dashed borders), and {b max = 0.5, g 2 = 0.68} (solid borders). The two additional light-gray bands correspond to the results with the "µ b scale" evolution scheme for {b max = 1.5, g 2 = 0.18} (dot-dashed borders) and {b max = 0.5, g 2 = 0.68} (solid borders). So, also for the {b † T , g log np } prescription we find the same ambiguity as for the {b * T , g lin np } one in Fig. 3: the overlap of the light-gray band with solid borders and of the band with dot-dashed borders indicates that two different evolution schemes give similar results with different evolution parameters {b max , g 2 }. Hence, we wonder if this similar trend suggests that it might not be possible to distinguish between the two schemes. Again, the possible way out is to look at the dependence of the results upon the fractional energy of the final hadrons.  In Fig. 6, the normalized multiplicity of Eq. (5.1) is shown as a function of P 2 1⊥ = z 2 1 q 2 T at the Belle scale Q 2 = 100 GeV 2 for the "µ b scale" evolution scheme. Also in this plot, there are two groups of uncertainty bands. A group displays the results for the {b † T , g log np } prescription in the standard notation, i.e. for {b max = 1.5, g 2 = 0.18} (dot-dashed borders), {b max = 1, g 2 = 0.43} (dashed borders), and {b max = 0.5, g 2 = 0.68} (solid borders). The group of two light-gray bands correspond to the results with the {b * T , g lin np } prescription for {b max = 1.5, g 2 = 0.18} (dot-dashed borders) and {b max = 0.5, g 2 = 0.68} (solid borders). If we focus on the left panel where calculations are performed at z 1 = z 2 = 0.5, the two bands with dot-dashed borders are substantially overlapped, thus reinforcing the suspect that it might not be possible to discriminate between the {b * T , g lin np } and {b † T , g log np } prescriptions. But if we now turn to the right panel, where the same calculation is performed at z 1 = 0.3, z 2 = 0.5, we may hope to have a sufficiently small experimental error that discriminates between the two bands with dot-dashed borders. Unfortunately, the plot suggests also that this option seems possible only for the {b max = 1.5, g 2 = 0.18} case. And further explorations show that the same calculation, when performed in the "fixed scale" evolution scheme, produces more confused results. In summary, a combined study of the z and P 2 1⊥ dependencies in the normalized multiplicity might be able to discriminate among different prescriptions for the nonperturbative effects in the evolution only for a selected set of evolution parameters and schemes.

Sensitivity to hadron fractional-energy dependence
In the previous sections, we found that in several occasions only the combined study of the z and P 2 1⊥ dependencies of the normalized multiplicity allows for discerning results obtained from different parametrizations and prescriptions in the description of nonperturbative effects in the TMD evolution. This is not accidental. With the approximations adopted in this work, the main difference between the two considered evolution schemes lies in fact in the z dependence of the collinear fragmentation function d 1 , as it can be deduced by comparing Eqs.  Figure 7. The normalized multiplicity at z 2 = 0.5 as a function of P 2 1⊥ = z 2 1 q 2 T at the Belle scale Q 2 = 100 GeV 2 for the evolution parameters {b max = 1.5, g 2 = 0.18} and with the {b * T , g lin np } prescription for the transition to the nonperturbative regime (see text). Uncertainty band with dot-dashed borders for z 1 = 0.3, with dashed borders for z 1 = 0.5, with solid borders for z 1 = 0.7. The squared box with error bar corresponds to an experimental error of 7%. Left panel for the "µ b scale" evolution scheme, right panel for the "fixed scale" one.
The plots in Fig. 7 seem to confirm this finding. In the left panel, the normalized multiplicity of Eq. (5.1) is shown as a function of P 2 1⊥ = z 2 1 q 2 T at the Belle scale Q 2 = 100 GeV 2 for the "µ b scale" evolution scheme, the {b * T , g lin np } prescription, and the choice {b max = 1.5, g 2 = 0.18}. The bands display results for the values z 1 = 0.3, z 2 = 0.5 (band with dot-dashed borders), z 1 = z 2 = 0.5 (dashed borders), and z 1 = 0.7, z 2 = 0.5 (solid borders). In the right panel, we show the results of the calculations performed in the same conditions but for the "fixed scale" evolution scheme. It is quite evident that the latter scheme produces P 2 1⊥ distributions that are systematically larger for any combination of {z 1 , z 2 }. This finding holds true also for other choices of the evolution parameters {b max , g 2 } and for the {b † T , g log np } prescription.

5.5
Sensitivity to the hard scale: from Belle to Bes-III All previous results have been obtained at the Belle scale of Q 2 = 100 GeV 2 . We may wonder what happens when reducing the "evolution path" to lower scales, like, e.g., the Bes-III scale Q 2 = 14.6 GeV 2 .  In Fig. 8, the normalized multiplicity of Eq. (5.2) is shown as a function of P 2 1⊥ = z 2 1 q 2 T ≡ (0.5) 2 q 2 T in the same conditions and notation as in Fig. 2 but at the Bes-III scale Q 2 = 14.6 GeV 2 . By comparing these results with the ones in Fig. 2, we deduce that the net effect is a systematic enlargement of the uncertainty bands. This finding occurs also for other combinations of evolutions schemes and nonperturbative prescriptions. Hence, we deduce that working at the Bes-III scale is not useful if we want to discriminate among different evolution parameters {b max , g 2 }, or between the {b * T , g lin np } and {b † T , g log np } prescriptions, or between the "fixed scale" and "µ b scale" evolution schemes.
However, we recall that each uncertainty band is the envelope of the 68% of 200 different curves, each one corresponding to a specific replica of the intrinsic parameters entering the Gaussian widths P 2 ⊥ a h (z) of Eq. (3.1) for the b T distribution of the D a 1 at the starting scale in the evolution. Then, we might envisage that the experimental error is sufficiently smaller than the band width such that it is able to discriminate some of the replicas, in order to narrow the uncertainty on the intrinsic parameters. In any case, this goal will be achieved only by performing additional more precise measurements of SIDIS multiplicities for different final hadron species and on different targets.

Sensitivity to partonic flavor
The sensitivity to the nonperturbative intrinsic parameters, that describe the b T distribution of the TMD FF at the initial scale of evolution, is an important issue. The analysis of SIDIS multiplicities at low Q 2 suggests that some of these parameters are different for different flavors [27]. Hence, we expect that also the distribution in transverse momentum space of the evolved TMD FF will depend on the flavor of the fragmenting partons. However, the cross section in Eq. (2.9) mixes all flavors in the sum. Therefore, it is useful to define an observable that is well suited to explore the effect of flavor in the TMD evolution.
In the following, we will show results for the P 2 1⊥ distribution of ratios of normalized multiplicities corresponding to different final states: In Fig. 9, we show the ratio of Eq. (5.3) between the normalized multiplicity for {π + π − } and the one for {K + K − } at z 2 = 0.5 and y = 0.2 as a function of P 2 1⊥ = z 2 1 q 2 T at the Belle scale Q 2 = 100 GeV 2 for the "fixed scale" evolution scheme, for the evolution parameters {b max = 1.5, g 2 = 0.18}, and with the {b * T , g lin np } prescription for the transition to the nonperturbative regime.
If we suppose to switch off the flavor dependence of the intrinsic parameters, the b T distribution of the TMD FF in Eq. (3.28) is controlled by the same Gaussian width P 2 ⊥ (z) for all channels. This feature remains valid when performing the Bessel transform to momentum space, such that the q 2 T distribution of the cross section can be factorized out of the flavor sum. Therefore, if we take the ratio of normalized multiplicities at the same z 1 we expect the latter to be independent of P 2 1⊥ = z 2 1 q 2 T . This is indeed the result displayed in the left panel of Fig. 9. It is a systematic feature of the "fixed scale" evolution scheme: it holds true for other values of z 1 , as shown in the panel, but also for other combinations of nonperturbative evolution parameters and nonperturbative prescriptions.
If we account for the flavor dependence of the Gaussian widths P 2 ⊥ q→h (z), then the b T distribution is different for the {π + π − } final state from the one for {K + K − }, as it can be realized by inspecting Eqs. (4.29)-(4.35). Consequently, the ratio of normalized multiplicities has a specific P 2 1⊥ = z 2 1 q 2 T distribution that, of course, changes with z 1 . This is indeed the content of the right panel in Fig. 9: the uncertainty band of the 68% of 200 replicas of Gaussian widths with dot-dashed borders corresponds to z 1 = 0.3, the band with dashed borders to z 1 = 0.5, the band with solid borders to z 1 = 0.7.
Almost all the ratios are smaller than unity because in our approximations the fragmentation into kaons has two favoured channels while the fragmentation into pions only one (see Eqs. (4.29) and (4.34)), and the P 2 1⊥ distribution of the fragmentation into kaons seems to be larger than the corresponding one for pions (see the analysis of Ref. [27]). In any case, we believe that the inspection of the P 2 1⊥ distribution of ratios of normalized multiplicities for different final hadrons produced in future e + e − annihilation experiments is a useful tool to discriminate among different scenarios in TMD evolution. For example, if future data for this observable will lie well above unity, the "fixed scale" evolution scheme would be ruled out, independently of the flavor dependence of the intrinsic parameters in the TMD FF at the initial scale of evolution.
In Fig. 10, in the two panels of the upper row we show the same ratio of normalized multiplicities in the same conditions and notation as in the previous figure but for the "µ b scale" evolution scheme. The left panel still corresponds to the case when the flavor dependence of the intrinsic parameters is neglected. However, in the "µ b scale" scheme the b T distribution of the TMD FF is influenced also by the collinear part of the fragmentation function: the d q→h  (3.19). Hence, when performing the Bessel transform of D q 1 in the cross section, the resulting q 2 T distribution depends on the flavor of the fragmenting parton even if the intrinsic parameters do not. This "perturbative" flavor dependence, induced by RGE acting on the evolved collinear part of the TMD FF, mixes with the possible flavor dependence of the intrinsic parameters, making it rather difficult to disentangle the two effects. The left panel in the upper row shows the ratio of normalized multiplicities as a function of P 2 1⊥ = z 2 1 q 2 T for three different values of z 1 . As in the previous figure, the band with dot-dashed borders corresponds to z 1 = 0.3, the band with dashed borders to z 1 = 0.5, and the band with solid borders to z 1 = 0.7. Suprisingly, all the ratios are larger than unity. When including also the flavor dependence in the intrinsic parameters, the uncertainty bands become larger because there is a marked sensitivity to all possible replica values of the intrinsic parameters themselves. Again, as in the previous section we can argue that experimental data will have a sufficiently small error to discriminate among the various replicas.
A further constraint can be achieved by considering a different combination of final state hadrons in the ratio of normalized multiplicities in Eq. (5.3). The lower panel in Fig. 10 shows the results for the ratio between a {π + π − } final state and a {π + K − } final state when neglecting the flavor dependence of intrinsic parameters of the TMD FF at the initial 18 Figure 10. Upper panels: same as in previous figure but for the "µ b scale" evolution scheme. Lower panel: the ratio between the normalized multiplicities M π + π − (z 1 , z 2 = 0.5, q 2 T , y = 0.2)/M π + π − (z 1 , z 2 = 0.5, 0, y = 0.2) and M π + K − (z 1 , z 2 = 0.5, q 2 T , y = 0.2)/M π + K − (z 1 , z 2 = 0.5, 0, y = 0.2) as a function of P 2 1⊥ = z 2 1 q 2 T at the Belle scale Q 2 = 100 GeV 2 in the same conditions and with the same notation as in the upper panels, but for flavor independent intrinsic parameters of input TMD FF (see text).
scale. The notation and conventions are the same as in the other panels. All the ratios are now lower than unity. Hence, combining this result with the content of the upper left panel could represent a very selective test of the "µ b scale" evolution scheme. In fact, when neglecting the flavor dependence of intrinsic parameters the P 2 1⊥ distribution of normalized multiplicities for the {π + π − } final state should be larger than the one for {K + K − } at any z 1 , while at the same time it should turn out narrower than the one for {π + K − } at any z 1 . Moreover, if future data for the {π + π − } back-to-back production in e + e − annihilation will display a much narrower P 2 1⊥ distribution than for the {K + K − } production, at least by 20%, this will represent a further selective test for calculations performed in this evolution scheme, as it can be deduced by combining the results in the panels of the upper row.
Finally, we notice that because of charge conjugation symmetry (see Sec. 4.1) we predict that the ratio between normalized multiplicities leading to (π + , K − ) and (π − , K + ) final states should be equal to unity, irrespective of the choice of evolution schemes, nonperturbative evolution parameters and prescriptions. It would be interesting to cross-check this prediction by measuring this ratio as a function of P 2 1⊥ .

Conclusions
In this paper, we consider the semi-inclusive production of two back-to-back hadrons in electron-positron annihilations. We study the transverse momentum distribution of such pairs of hadrons by observing the mismatch between their collinear momenta, and we focus on charge-separated combinations of pions and kaons. We conveniently define the multiplicities in electron-positron annihilations as the differential number of back-to-back pairs of hadrons produced per corresponding single-hadron production, in analogy to the definition of multiplicity in SIDIS process. In particular, we analyze the multiplicities normalized to the point of vanishing transverse momentum in order to extract clean and uncontaminated details on the transverse momentum dependence of the functions describing the fragmentation process (transverse-momentum dependent fragmentation functions -TMD FFs). The normalized multiplicities are advantageous also because they turn out to be almost insensitive to the theoretical uncertainty related to the arbitrary choice of the renormalization scale.
We consider electron-positron annihilations at large values of the center-of-mass (cm) energy, namely in the experimental conditions of the Belle and Bes-III experiments. We study how TMD FFs evolve with the hard scale. The input expression for TMD FFs is taken from a previous analysis of SIDIS multiplicities measured by Hermes at low energy, which is assumed as the starting scale. Since the hard scale in annihilation processes is much larger, we perform realistic tests on the sensitivity to various implementations of TMD evolution available in the literature.
We find that within a specific evolution scheme the transverse momentum distribution of normalized multiplicities at the Belle scale can be very sensitive to the choice of the parameters describing the nonperturbative part of the evolution kernel. A hypothetical 7% error in such data (compatible with the observed experimental error in collinear backto-back emissions in electron-positron annihilations) could discriminate among different choices of parameters that are justified and adopted in the literature.
But we observe also that at the same Belle scale different evolution schemes with different nonperturbative parameters can give overlapping transverse momentum distributions. Our global results indicate that different evolution schemes can be discriminated only by considering the combined dependence of normalized multiplicities on both the transverse momentum and the fractional energy carried by the final hadrons. And this finding holds true (with some limitations) also for the purpose of discriminating among different prescriptions for describing the transition from nonperturbative to perturbative regimes in transverse momentum.
The dependence on the fractional energy of the final hadrons is contained in the collinear part of the TMD FFs. Different evolution schemes produce different evolution effects also in the collinear fragmentation functions, which in turn emphasize the differences in the final transverse momentum distribution of evolved TMD FFs. The dependence on the fractional energy is contained also in the average squared transverse momenta that describe the width of the input distribution of the TMD FFs at the starting scale. Therefore, by studying this dependence it may be possible to reduce the uncertainty on the intrinsic parameters that describe these input distributions.
To this purpose, focusing on the normalized multiplicities at the Bes-III scale looks more promising. In fact, we observe that in stepping down from Belle to Bes-III scale the transverse momentum distributions of normalized multiplicities become much more sensitive to the details of the input distribution at the starting scale. The uncertainty in the determination of the intrinsic parameters needed to fit the Hermes SIDIS multiplicities reflects in a larger spread of normalized multiplicities as functions of transverse momentum. At the Bes-III scale, a hypothetical experimental error of 7% does not discriminate among results coming from different nonperturbative evolution parameters or from different evolution schemes. But within a specific choice of evolution scheme it can discriminate among results that come from different values of the intrinsic parameters.
The Hermes results al low energy show significant differences between SIDIS multiplicities for final-state pions and kaons. Hence, these data were fitted using transverse momentum distributions for the input TMD FFs that contain flavor dependent parameters. Here, we explore also how the final results for normalized multiplicities at Belle and Bes-III scales are sensitive to the details of this flavor dependence at the starting scale. In doing so, we find that the most convenient observable is represented by the ratio of normalized multiplicities for different final hadron species, particularly at the Belle scale.
The most striking evidence is for evolution schemes where the flavor dependence is strictly localized only in the intrinsic parameters of the input TMD FFs at the starting scale. If we switch off such flavor dependence, the transverse momentum distribution of normalized multiplicities is always the same, irrespective of the species of final hadrons. So, if we select for example pions and kaons, the ratio of the corresponding normalized multiplicities is constant and equal to unity. If the flavor dependence of the intrinsic parameters is switched on, then the ratio deviates to values (mostly) lower than unity, in agreement with general expectations that kaons have a larger distribution in transverse momentum.
The situation is more confused for evolution schemes where the flavor dependence is indirectly contained also in the initial conditions of the evolution equations through the (flavor dependent) collinear part of the fragmentation functions. In this case, this effect mixes up with the flavor dependence contained in the intrinsic transverse momentum distribution, and it is difficult to disentangle one from the other. At variance with the previous class of evolution schemes, in this case the ratio of normalized multiplicities for pions with respect to kaons turns out to be (mostly) larger than unity. Fortunately, more selective criteria are offered by considering a variety of species of final hadrons. If we consider ratios of normalized multiplicities for pions with respect to mixed pion-kaon pairs, the results are (mostly) lower than unity. By combining the results for various final states all together, one would hope to constrain the arbitrary ingredients of TMD FFs as much as possible.
We conclude by stressing that all the results and remarks above refer to the unpolarized TMD FFs that describe the fragmentation of an unpolarized parton into an unpolarized hadron. However, this function is an essential ingredient in all the (spin) azimuthal asym-metries extracted in hard processes like electron-positron annihilation, hadronic collision, and SIDIS. Hence, a better control on the transverse momentum dependence of unpolarized TMD FFs implies also a better knowledge of polarized TMD FFs as well as of (un)polarized TMD parton distributions. For this reason, we are looking forward to a multidimensional analysis of data accumulated by the Belle and Bes-III collaborations, possibly including a study of normalized multiplicities for various hadron species as suggested in this work.