Transverse parton momenta in single inclusive hadron production in ${e^ + }{e^ - }$ annihilation processes

We study the transverse momentum distributions of single inclusive hadron production in ${e^ + }{e^ - }$ annihilation processes. Although the only available experimental data are scarce and quite old, we find that the fundamental features of transverse momentum dependent (TMD) evolution, historically addressed in Drell-Yan processes and, more recently, in Semi-inclusive deep inelastic scattering processes, are visible in ${e^ + }{e^ - }$ annihilations as well. Interesting effects related to its non-perturbative regime can be observed. We test two different parameterizations for the $p_\perp$ dependence of the cross section: the usual Gaussian distribution and a power-law model. We find the latter to be more appropriate in describing this particular set of experimental data, over a relatively large range of $p_\perp$ values. We use this model to map some of the features of the data within the framework of TMD evolution, and discuss the caveats of this and other possible interpretations, related to the one-dimensional nature of the available experimental data.


Introduction
Transverse momentum dependent distribution and fragmentation functions (TMDs) are fundamental tools to understand the structure of nucleons in terms of their elementary constituents, quarks and gluons. TMDs are non-perturbative quantities which embed important correlations among partonic and hadronic intrinsic properties, like spin or orbital angular momentum, and their internal transverse motion.
TMD parton distribution functions can be interpreted, at leading twist, as number densities of partons carrying a lightcone momentum fraction x of the parent nucleon momentum P. Unpolarized and polarized TMD distribution functions have extensively been studied in Drell-Yan processes and semiinclusive deep inelastic scattering (SIDIS) in the past; newgeneration, dedicated experiments are currently running (like Drell-Yan at COMPASS or at RHIC) or are being planned (like the Electron-Ion Collider in the US and the AFTER proposal at CERN-LHC).
Of equal importance are the TMD fragmentation functions, which embed fundamental information on the hadronization process, where a hadron h, carrying a light-cone fraction z of the fragmenting parent parton, is produced. TMD fragmentation functions can be measured in single-or double-inclusive hadron production in e + e − annihilation processes or, in a more involved way, in SIDIS, where they necessarily couple to a TMD distri-Email addresses: boglione@to.infn.it (M. Boglione), jogh@jlab.org (J. O. Gonzalez-Hernandez), r.taghavi@stu.yazd.ac.ir (R. Taghavi) bution function. Even with the best SIDIS data presently available, several complications remain to be solved. Extensive recent studies can be found, for example, in Refs. [1][2][3][4][5] In e + e − collisions, at c.m. energies below the Z 0 mass, the electron and positron predominantly annihilate to form a single virtual photon, which can subsequently produce a qq pair. The quark and anti-quark will then convert into hadrons. At sufficiently high energies, these multi-hadronic events are expected to form two back-to-back jets (due to the limited transverse momentum along the original quark direction). Single inclusive distributions in variables relative to the jet direction, which is expected to be the quark direction, will therefore give information about the fragmentation of quarks into hadrons. In particular, the dependence of the inclusive distributions in momentum transverse to the jet axis will provide the golden channel toward the phenomenological extraction of TMD FFs. An illustration of this process is given in Figure 1.
While much effort has been put into measuring and extracting unpolarized and even polarized TMD PDFs (the Sivers function is a well-known example), little or no experimental information on TMD FFs is presently available. BELLE and BaBar Collaborations have recently presented new multidimensional data analyses of the Collins asymmetry in e + e − → h 1 h 2 X processes, which have allowed a first glance to the intrinsic transverse motion of partons inside hadrons through the extraction of the Collins function, a polarized, chirallyodd TMD fragmentation function. Relevant recent literature on this and other e + e − related subjects are presented in Refs. [6][7][8][9]. However, no modern measurements exist of the unpolarized TMD FFs, although a thorough knowledge of this function would be of fundamental importance for any TMD study.
While waiting for up-to-date, high statistics and (possibly) multidimensional results on the p ⊥ distribution of e + e − unpolarized cross sections (the BELLE Collaboration has already presented some of their preliminary Monte Carlo simulations at SPIN 2016 [10]), we concentrate on a set of rather old measurements of single inclusive hadron production in e + e − annihilation processes, e + e − → h X, by the TASSO Collaboration at PE-TRA (DESY) [11,12]: p ⊥ distributions were provided for four different c.m. energies between 14 and 44 GeV, corresponding to charged particle production summed over all charges and all particle species, with no flavor separation. Cross sections are given as functions of p ⊥ , integrated over the energy fraction z h = 2E h / √ s of the detected hadron h. Note that, up to corrections of order p 2 ⊥ /s, z h coincides with the light-cone momentum fraction z. Although no information is offered about possible cuts applied to z h , average values are provided for each c.m. energy set, as summarized in Table 1: they correspond to rather low values, ranging from z h = 0.13 at 14 Gev to z h = 0.08 at 44 GeV. Together with TASSO data, we also consider the analogous MARKII Collaboration measurements [13], collected at the SLAC storage ring PEP, at a fixed c.m. energy of 29 GeV, and PLUTO data on the average transverse momentum square [14], collected at PETRA (DESY).
Crucial to all these data is the correct determination of the jet axis, to which the p ⊥ distributions are most sensitive, beside proper treatment of geometric acceptance effects, trigger bias, kinematics cuts and radiative corrections. These corrections, obtained by comparison to Monte Carlo simulations, are somehow model dependent. Clearly all these issues introduce very large uncertainties which, according to modern standards, were largely underestimated.
Although these data are affected by several limitations (no hadron separation, limited coverage, z h integration, difficulties in reconstructing the jet axis, etc), they represent an extremely interesting example of a direct measurement of intrinsic transverse momenta. In fact, as mentioned above, the jet axis resulting from each e + e − scattering identifies the direction of the fragmenting qq pair (q andq should be back to back in the e + e − c.m. frame if radiative effects are appropriately subtracted), and the detected p ⊥ represents a direct measurement of the transverse momentum of the final hadron with respect to the fragmenting parent parton, see Figure 1.
The purpose of this article is two-folded: first, we will devise and test an appropriate functional form to describe the p ⊥ distributions measured by TASSO and MARKII, achieving as much information as we possibly can on the TMD unpolarized FF; then, a careful interpretation of our results will be provided focusing on the features related to TMD factorization within a TMD evolution scheme, in the non-perturbative regime.
As these measurements offer quite limited information, we will not be able to perform a detailed extraction of the TMD fragmentation functions as done in the past [2,8,[15][16][17]. However, we will observe that interesting signatures of TMD factorization in the non perturbative regime can be detected in these data sets. In particular, we will find indications that a powerlaw p ⊥ behavior, different from the Gaussian parametrization of the TMD FFs we usually used in our previous analyses, might reproduce these data more successfully, especially as p ⊥ grows larger then a few hundred MeV and we enter the region in which TMD evolution effects start to become more visible.

Formalism
Similarly to the collinear case [18], the e + e − → hX cross section can be casted in the following form: Within TMD-factorization, up to power suppressed terms, the hadronic tensor W can be expressed as The first term on the right hand side of Eq. (2), W µν T MD , corresponds to the region of small transverse momenta, while the second, W µν coll , is calculable within collinear factorization and contains corrections that become important at larger values of p ⊥ . In the case at hand, as there is only one observed final hadron, one may write for the TMD term: where Q is the hard scale of the process, z and p ⊥ are the observed hadronic variables, µ is the renormalization scale, while ζ D is a regulator for the light-cone divergences that arise in TMD-factorization. The hard factor H f can be calculated in perturbation theory, while the TMD FF is defined by the relations Notice the conjugate variable to b ⊥ is k ⊥ rather than p ⊥ (k ⊥ is the component of the fragmenting parton momentum transverse to the final hadron, in a reference frame in which the latter is purely longitudinal). Details on the various ingredients of Eq. (5) can be found in Refs. [18,19]. For our purposes it suffices to note that in the first three factors of Eq. (5) all the ingredients are calculable within perturbation theory, except for the collinear fragmentation function d h/ j . Notice that the Wilson coefficients,C j/ f , the integral of the anomalous dimensions, γ D and γ K , and the the Collins-Soper (CS) kernel,K, depend on b ⊥ only through the quantity b * , which is set to remain smaller than some maximum b max . The last factor in Eq. (5) contains only non-perturbative information.
There is some freedom in the definition of Eq. (5), encompassed in the arbitrary quantities µ, D . For our analysis we adopt the usual choices [19] With these choices the CS kernelK vanishes at order α s , so the third factor in Eq. (5) reduces to one. To this same order, the factor containing the anomalous dimensions γ D and γ K can be expressed in a closed analytic form. For our purposes, it will be useful to classify the result of this integral in terms of its dependencies on Q and b * , namely for which we have used the results of the Appendices in Ref. [3]. The quantities N Γ (Q), f Γ (b * , Q 0 ) and λ Γ (b * ) are flavor independent functions that encode the most prominent perturbative effects in the definition of the TMD FF. The last of these three functions, is the most interesting since it correlates b * with Q, which means it has the effect of modifying the shape of the TMD FF under evolution. With these considerations, one may write for the TMD FF: Therefore, except for the overall normalization factor N Γ (Q), the effects of evolution at order α s can be mapped to either the non-perturbative function g K (b ⊥ ), or the perturbative quantity In the region where TMD effects dominate (see Eqs. (1)-(3)), flavor independence of g K and λ Γ in Eq. (8) implies where the symbol F −1 indicates the two-dimensional inverse Fourier transform, from momentum to impact parameter space, and the transformation b ⊥ → z b ⊥ is needed to account for the extra factor of z that appears in the definition of Eq. (3), compared to the TMD term in the hadronic tensor in Eq. (2). In the following section we will use relations (7)-(9) to make an interpretation of our results.

Data fitting and results
A full analysis within a TMD-evolution scheme should include all of the contributions in Eq. (2), as well as a matching prescription to interpolate between regions of small and large p ⊥ , as originally prescribed in Ref. [20]. It is important to stress that all of these ingredients provide crucial constraints that any full analysis should include. Some recent studies related to the complications involved in the matching prescriptions for SIDIS processes are presented in Refs. [4,21]. However, for such type of analyses, multidimensional data sets are most suitable, where one can completely disentangle the effects of different kinematics variables. In the case of the measurements of Ref. [11][12][13][14], the large systematic uncertainties induced by z-integration can only render limited insight on TMD-effects. However, it is still interesting to investigate what information about evolution can be extracted from these measurements. In fact, even for these z-integrated cross sections, one may expect the shape of the p ⊥ distributions to be affected by TMD-evolution effects.
In order to address this question it becomes essential to make an estimate of where the large transverse momentum corrections start becoming important in the TASSO and MARKII data sets. To do so, we start by considering the errors of the TMD approximation, which are of order O(k ⊥ /Q) = O(p ⊥ /zQ). In general, one expects that at p ⊥ ∼ zQ, the cross section should receive contributions from the collinear term in Eq. (2). Using the average values of Table 1, one may estimate that this contributions should be significant at p ⊥ ∼ 2 GeV. We identify this as the matching region. In this article we will attempt to extract information only about the non-perturbative evolution carried by g K in Eq. (8), so we will constrain our analysis to the region of p ⊥ < 1.0 GeV. Note that we also impose a lower cut on p ⊥ , such that p ⊥ > 0.03 GeV; this amounts to excluding the first data point in the TASSO data sets.
In order to relate the shape of the data to possible TMDevolution effects, one needs a model that can reproduce the transverse momentum distribution of the final hadron h. Then, by looking at the b ⊥ -space, one may connect the parameters of the model to some of the information contained in the definition of Eq. (5), as discussed in Ref. [3].
At least two functional forms have been shown to appropriately describe transverse momentum distributions [2,3].
• Gaussian form: it is the most commonly used parametrization for phenomenological studies, it has been shown to reproduce experimental data very successfully in both Drell-Yan and SIDIS processes. It has the advantage of being easy to integrate analytically.
• Power-law: it is very flexible, even with a limited number of free parameters. Not only can it appropriately reproduce the behavior of the cross section at small p ⊥ , but it can also incorporate its tail at larger p ⊥ values.
We will consider both of these functional forms in order to describe the low-p ⊥ region of the data. Our aim is to focus on the kinematics ranges where TMD-effects are dominant. As discussed above, one may estimate that perturbative effects will start to become important roughly around p ⊥ ∼ 2 GeV, but could in fact be non-negligible at even smaller values of transverse momentum, especially as the data we consider are integrated over z. Our working hypothesis is that for p ⊥ < 1 GeV the TMD-term in Eq. (2) is the largest contribution to the cross section.
We model the structure function F 1 in Eq. (1) so that which extends the leading order expression for the collinear cross section. In Eq. (10), the sum runs over all q andq flavors, and we have assumed that the TMD fragmentation function may be written as where D h q (z) is the collinear, unpolarized FF (which we take from Ref. [22]). The function h(p ⊥ ) incorporates all of the p ⊥ dependence of the TMD FF, it is flavor independent and it is normalized so that it integrates to unity.
The TASSO collaboration provides cross sections differential in p ⊥ , normalized to the fully inclusive cross section, which at leading order read so that our fits will involve the expression A note of caution is necessary at this point. The TASSO distributions in p ⊥ at different energies cannot be described by simply using the model of Eq. (13). Instead, one must incorporate a treatment for their normalizations at different values of Q. While this may be in conflict with a possible probabilistic interpretation of the function h(p ⊥ ), it should not affect our conclusions regarding TMD evolution, since for that, we will focus on aspects of the p ⊥ distributions that regard their shape, but not the precise values of their maxima. We note that accounting for all the features of these data may be challenging even within a full TMD analysis, as large systematic errors may translate into out-of-control normalizations when dealing with z-integrated data. The fits in the next subsections are performed on the TASSO p ⊥ -distributions only. We will use the MARK II p ⊥ -dependent normalized cross section and the PLUTO measurements of p 2 ⊥ to cross-check our results. When appropriate, error bands corresponding to a 2σ confidence level are provided, obtained by generating random points in the parameter space for which χ 2 i ∈ [χ 2 0 , χ 2 0 + ∆χ 2 ], where χ 2 0 is the minimal value given by the fit and ∆χ 2 depends on the number of parameters of the model; the relevant cases in our fits involve either 6 or 7 free parameters, which correspond to ∆χ 2 values of 12.85 and 14.34, respectively.
3.1. Gaussian shape at low p ⊥ We start by applying a Gaussian model with a constant width: As mentioned earlier, in order to describe the data we need an appropriate treatment for the normalization. Unexpectedly, we find that using different multiplicative constants, one for each energy, is not enough to obtain a good fit, even in the limited region of p ⊥ < 0.5 GeV (see the first entry in Table 2). Instead, by also introducing a Q-dependent shift for the cross sections, so that one has one can obtain a good description of the data, as seen in the second entry of Table 2 and on the upper panel of Figure 2.  This unconventional prescription to deal with the normalization, while leaving little room for a partonic interpretation for the function h(p ⊥ ), allows us to verify quantitatively that as far as the width of the p ⊥ distribution is concerned, no significant change can be observed with growing Q. One may not conclude, however, that TMD-effects do not appear in these transverse momentum ranges, but rather that the data analyzed do not have the necessary accuracy to show possible width changes in this region. Thus, one must try to extend the description of the data to larger values of p ⊥ . For p ⊥ > 0.5 GeV, a noticeable dependence of the distributions on the c.m energy suggests that a Q-dependent width may be appropriate. However, even with this extension of the model, describing the data past this point turns out to be extremely difficult. To illustrate this, we consider the functional form which allows for some comparison to previous phenomenological studies [23,24] where, within a CSS evolution scheme, a Gaussian behavior of the non-perturbative Sudakov factor in the b ⊥ space was assumed. This time, we use a multiplicative constant for each value of Q (δ = 0 in Eq. (15)), and fix Q 0 = 1.6 GeV. The last entry of Table 2 shows the results of the fit for p ⊥ < 1 GeV, using the extended Gaussian model corresponding to Eq. (16). The obtained minimal χ 2 value points to a rather low quality of the description of the data. The corresponding plot, in the lower panel of Figure 2, shows that there is some tension between a successful description of the peak at low p ⊥ and an equally good description of the tail at larger p ⊥ values. As we will show in the next subsection, the power-law can in fact accommodate for both of these features of the data.
3.2. Power-law shape at low and moderate p ⊥ So far, we have tried the Gaussian ansatz for the p ⊥ distributions. We have found that the data favors a constant width, for values up to p ⊥ = 0.5 GeV. However, this can only be achieved by introducing a Q-dependent shift, which cannot be easily interpreted within a partonic picture. The Gaussian class of models does not seem to be appropriate for larger values of p ⊥ .
In order to describe the data up to p ⊥ = 1 GeV, we test a power-law parametrization, given by where the factor 2(α − 1)M 2 (α−1) is set so that h(p ⊥ ) integrates to unity. The two-dimensional inverse Fourier transform of this function has an exponentially decaying asymptotic behavior, consistent with what would be expected from general arguments within quantum field theory, as discussed in Refs. [18,25]. In impact parameter space the power-law becomes where K 1−α (b ⊥ M) is the modified Bessel function of the second kind. In what follows, we use an independent normalization for each value of Q. As discussed before, we are interested in the shape of the distributions, rather than on their overall normalizations.
In the power-law parametrization of the p ⊥ -differential cross section, the parameter M 2 is related to the position of its peak, p 0 ⊥ , by the relation M 2 = 2 α − 1 p 0 ⊥ 2 . Since the studied distributions reach their maximum at roughly the same value of transverse momentum, p ⊥ ≈ 0.212 GeV, for all values of the c.m. energy, in our main analysis we have imposed the conditions that We have verified that setting M 2 free, does in fact satisfy Eqs. (19) within errors. It is, however, useful to reduce the number of parameters by directly imposing these relations.  First, to test that the power-law can appropriately describe the data, we conducted simple independent fits for each value of the c.m. energy, which renders one value of α for each Q. The most interesting aspect of this preliminary fit is that it shows a clear dependence of the parameter α on Q, despite the use of independent normalizations. The trend of the values for α is displayed in Figure 3, which shows a decrease of the its optimal value with Q. Due to the large uncertainties in the determination of α, there are likely several functional forms that can accommodate the observed behavior. In order to make an educated guess for a suitable Q-dependence in α, we assume that the integration over z does not alter the structure of the relation in Eq. (9), namely for some functiong(b ⊥ ). Thus, one can see that a logarithmic behavior for α may be appropriate by looking at the asymptotic limit of the power-law in b ⊥ -space, Eq. (18). First, for the values α 0 and α that describe the data at two given values Q 0 and Q, if (20) holds, one should have in the large-b ⊥ limit. This can be achieved if for some valuesα and ν, which in turn provides the relation We have implemented Eq. (23) into a fit and confirmed that in fact it reproduces well the TASSO data. Results are shown on the third panel of Table 3, and in the top plot in Figure 4, which includes errors corresponding to a 2σ confidence level.
We compare these results to the MARK II data set in the bottom panel of Figure 4. The argument presented above, leading to a logarithmic Qdependence for the power α, has some caveats. First, it depends on whether Eq. (20) is approximately correct. Furthermore since it considers only the asymptotic large-b ⊥ behavior of Eq. (18), Eq. (23) does not need to hold for values larger than p ⊥ ∼ M. Therefore, even if one can describe the data, any interpretation of the logarithmic behavior of α, in terms of the ingredients that define the TMD FF, Eq. (5), should be taken with great caution. Nonetheless, it is worthwhile to explore this possibility.
Notice that the functiong acquires its behavior from λ Γ and g K . Since the first varies slowly with b ⊥ , and in fact does freeze to a constant value at large enough b ⊥ (see Eq. (7)), one may see Eq. (23) as the manifestation of a logarithmic large-b ⊥ trend for g k (b ⊥ ), analogous to Eq. (22). A behavior consistent with discussions in Refs. [3,25]. However, we stress that one may well reproduce the data by using different assumptions. As a counter example we consider a simple picture where the cross section takes the form and where one accounts for the integration over z by some average value z , leading to where the parameters β 1 , β 2 andM are to be determined by a fit. In this case, it is indeed possible to obtain a good description of the data by using the experimental average values of z of Table 1, since they exhibit a seemingly logarithmic trend. However, it is difficult to make a connection to TMD evolution since the values in Table 1 are in general affected by correlations between Q and z of different origin. In fact, notice that the condition of Eq. (19) implies a logarithmic behavior for our fit parameter M 2 . It is possible that, for instance, the effects of the TMD evolution, encoded in g K , result in changes in the power α that fit the data, while the changes in M 2 are the result of correlations induced by the integration over z. It seems, however, that the lack of information about the z-dependence of the TMD FF in the TASSO and MARK II measurements hinders a more solid conclusion about TMD evolution effects in these data sets.
Finally, it is useful to test the model of Eqs. (17) and (23) for larger values of p ⊥ . In fact, this is necessary to calculate p 2 ⊥ , since it implies integration over the full range of p ⊥ . For this, we keep the first of the constraints (19), but free the parameter p 0 ⊥ . First, we perform a fit including data up to p ⊥ = 2 GeV. As seen in Figure 5 and the last entry of Table 3, the model can successfully accommodate this extended range of p ⊥ . This range is, however, not enough to reproduce to a good accuracy the corresponding TASSO measurements of p 2 ⊥ . Thus, we further extend the range analysis of TASSO data to values up to p ⊥ = 3 GeV and use the resulting minimal parameters to estimate p 2 ⊥ . Figure 6 shows our estimate and the data from TASSO and PLUTO.

Final remarks
TMD FFs embed the essence of hadronization, one of the most important manifestations of QCD in the non-perturbative regime. It is therefore important to gather as much information as possible on these soft quantities, which cannot be computed, but have to be inferred from experiment. Over the last few years, several analyses have been performed to extract the polarized TMD FFs, like the Collins function, using the measurements of the Collins asymmetries in e + e − → h 1 h 2 X processes provided by BELLE and BaBar Collaborations, which delivered multidimensional data (in bins of z 1 , z 2 , p ⊥1 , p ⊥2 ) with impressive statistics and very high precision. Unfortunately, no analogous data have been provided on the p ⊥ -distributions of the unpolarized cross sections, or multiplicities, to allow for  17) and (23), in the range 0.03 GeV < p ⊥ < 1.0 GeV. To avoid overlapping and provide a clear display of the 2σ error bands, we plot the distributions for different energies applying an arbitrary shift. In the bottom panel we compare the results from the fit to TASSO data to the MARK II cross section. Note that a different normalization and its uncertainties have to be determined independently for this data set. We dont display the first bin, centered at p ⊥ = 0.05 GeV.
the extraction of the unpolarized TMD FFs. The absence of these fundamental bricks encumbers the analysis of any other polarized, as well as unpolarized, process. At present, the only available data are some old (and almost forgotten) measurements of e + e − → hX cross sections from TASSO and MARKII Collaborations. Although these are one-dimensional data and are affected by large uncertainties, as discussed in Section 1, they have the unique advantage of delivering measurements at different c.m. energies; therefore, they can provide a valuable starting point not only to learn about the unpolarized TMD FFs, but also to study the physiognomy of their TMD evolution.
In this article we assess the extent to which the effects of TMD evolution can be observed in these data. For this purpose, our main tool is an analysis based on a simple partonic picture in which the cross section is factorized, as in Refs. [2,8,[15][16][17]. While these types of analysis typically use a Gaussian form, we also test a power-law behavior to describe the data, as suggested in Ref. [3]. We extend this class of models to the case in which the parameterization is supplemented by some Q dependence, and use these results to provide an interpretation within a TMD    17) and (23) compared to TASSO p ⊥ -dependent distributions [12], in the range 0.03 GeV < p ⊥ < 2.0 GeV. Error bands are computed using a 2σ-confidence level, as explained in Section 3. evolution framework, see Eq. (5), discussing the caveats related to the limited amount of information provided by these data. We start by modeling the p ⊥ dependence of the cross section by a Gaussian shape and fit the four sets of TASSO cross section data (corresponding to four different c.m. energies) to extract the corresponding free parameters, see Table 2. Our analysis shows that the Gaussian distribution can only describe the data up to p ⊥ ∼ 0.5 GeV, provided the cross sections are adjusted with an ad-hoc, additive term δQ. In this region no  Figure 6: Estimation of the transverse momentum mean value, p 2 ⊥ , obtained by using the parameters extracted by fitting the TASSO p ⊥ distributions up to 3.0 GeV. Empty squares correspond to PLUTO data [14] while filled diamonds correspond to TASSO measurements [12]. The shaded area represents the uncertainty of our calculation and is computed as explained in Section 3.
Q-evolution effects can be observed in the Gaussian width of the p ⊥ -distributions. The difficulties related to the interpretation of these results, however, leads us to consider a different parameterization.
We then focus on a power-law parametrization of the p ⊥ dependence of the cross section, which shows to be more appropriate than the Gaussian model and provides a successful description of the TASSO p ⊥ -distributions over a much larger range of p ⊥ values. As a consistency check, we compare the results of our main fits on TASSO data, reported in the last two panels of Table3, with the MARK II and PLUTO data. We find a reasonable agreement of our model in both cases. Finally, we provide an argument to explain that the Q dependence in the power-law can be consistent with a logarithmic behavior, in the large b ⊥ limit, of the function g K , which encodes the nonperturbative evolution effects in the definition of the TMD FF.
The nature of these data forces us to be cautious with the interpretation of our results. In fact, it seems unlikely that these data by themselves would allow to disentangle the TMD effects from other Q-dependence in the data. This is related to the integration over z, which induces a degree of ambiguity in the possible interpretations. Thus, one should further test any conjecture with multi-dimensional data.
In the foreseeable future, unpolarized single-hadron production at fixed energies by BELLE and BaBar Collaborations, differential in both z and p ⊥ , may indeed provide enough constraints on the z-dependence of the fragmentation process, allowing for the possibility of a full TMD analysis when combined with the TASSO and MARK II data.