Unpolarised Transverse Momentum Dependent Distribution and Fragmentation Functions from SIDIS Multiplicities

The unpolarised transverse momentum dependent distribution and fragmentation functions are extracted from HERMES and COMPASS experimental measurements of SIDIS multiplicities for charged hadron production. The data are grouped into independent bins of the kinematical variables, in which the TMD factorisation is expected to hold. A simple factorised functional form of the TMDs is adopted, with a Gaussian dependence on the intrinsic transverse momentum, which turns out to be quite adequate in shape. HERMES data do not need any normalisation correction, while fits of the COMPASS data much improve with a $y$-dependent overall normalisation factor. A comparison of the extracted TMDs with previous EMC and JLab data confirms the adequacy of the simple Gaussian distributions. The possible role of the TMD evolution is briefly considered.


INTRODUCTION
Tranverse momentum dependent parton densities (TMD-PDFs) and fragmentation functions (TMD-FFs), often collectively referred to as TMDs, have recently attracted and keep receiving a huge amount of interest. TMD-PDFs show remarkable spin and angular momentum correlation properties of quarks and gluons and allow a direct connection to the internal 3-D partonic structure of hadrons [1]. TMD-FFs reveal universal features of the hadronization process and couple to TMD-PDFs in several physical observables.
Studies of TMDs are mostly performed in polarised Semi Inclusive Deep Inelastic Scattering (SIDIS) processes, N → h X, in the framework of the TMD factorisation scheme [2][3][4][5], according to which the SIDIS cross section is written as a convolution of TMD-PDFs, TMD-FFs and known elementary interactions. The role of TMDs in such processes has been definitely established by the observation and interpretation of single spin asymmetries, which confirmed the existence of polarised TMDs like the Sivers distribution [6,7] and the Collins fragmentation function [8].
The role of TMDs is already evident in unpolarised cross sections, simply by looking at the transverse momentum, P T , distribution of the final hadron in the γ * − N centre of mass frame, or, at order 1/Q, at the azimuthal dependence of the hadron around the γ * direction. In Ref. [9] a first investigation of SIDIS unpolarised cross sections was performed, mainly based on the EMC Collaboration experimental data [10], gathered from SIDIS experiments at different energies and off different targets. This analysis assumed a simple factorised and Gaussian parameterisation of the TMDs, which was later confirmed by the independent study of Ref. [11], based on data from JLab [12,13] and HERMES [14].
However, very recently, plenty of new data on SIDIS multiplicities have been made available by the HERMES [15] and COMPASS [16] Collaborations. Both of them have performed multivariate analyses of their measurements, providing extremely rich data sets, arranged in independent bins of the kinematical variables, which give a unique opportunity to extract new and more detailed information on the unpolarised TMDs. Also the theoretical investigation of TMDs has much progressed lately [2,3,5,[17][18][19][20] with the study of the QCD evolution of the Sivers and unpolarised TMDs -the so-called TMD evolutionwith the first phenomenological applications [21][22][23][24][25][26][27].
In this paper we analyse the latest HERMES and COMPASS data on unpolarised multiplicities aiming at improving our knowledge of the unpolarised TMDs. We reconsider, with the support of the new data, the first extraction of Ref. [9]; somewhat surprisingly, it turns out that the simple factorised form of the TMDs with the original, flavour independent, Gaussian parameterisation, still works rather well. However, the observed (Gaussian) dependence of the SIDIS cross section on the hadron transverse momentum, P T , is generated by a combination of the (Gaussian) dependences in the quark TMD-PDF and TMD-FF; thus, it is rather difficult to fix separately the parameters of the two Gaussians by studying only unpolarised multiplicities.
Although the HERMES and COMPASS data cover similar Q 2 regions (1 ≤ Q 2 ≤ 10 GeV 2 ), they differ in the experimental set-up, in the statistics, in the binning choices and in the explored x B range; in addition, there seems to be some discrepancy between the two measurements. We then fit the HERMES and the COMPASS multiplicities separately. A simultaneous fit of both sets of data would lead to poor results and is not presented here.
Recently, another study of the unpolarised TMDs has appeared [28], which follows a procedure somehow similar to that of this work, but which considers only the HERMES set of experimental data and does not include any attempt to check for signs of scale evolution.
After a short Section II devoted to the formalism, we present our main results in Section III. In Section IV we briefly discuss the possible role, and look for possible signs, of TMD evolution. In Section V we compare our present results with those of previous analyses [9,11] and check their consistency with other measurements of SIDIS cross sections and P T -distributions [10,12,13,29] which were not included in our fits. Further comments and concluding discussions are presented in Section VI.

II. FORMALISM
The unpolarised + p → h X, SIDIS cross section in the TMD factorisation scheme, at order (k ⊥ /Q) and α 0 s , in the kinematical region where P T k ⊥ Q , reads [30,31]: In the γ * − p c.m. frame the measured transverse momentum, P T , of the final hadron is generated by the transverse momentum of the quark in the target proton, k ⊥ , and of the final hadron with respect to the fragmenting quark, p ⊥ . At order k ⊥ /Q it is simply given by As usual: and the variables x, z and p ⊥ are related to the final observed variables x B , z h and P T and to the integration variable k ⊥ . The exact relations can be found in Ref. [9]; at O(k ⊥ /Q) one simply has The unpolarised TMD distribution and fragmentation functions, f q/p (x, k ⊥ ) and D h/q (z, p ⊥ ), depend on the light-cone momentum fractions x and z and on the magnitudes of the transverse momenta k ⊥ = |k ⊥ | and p ⊥ = |p ⊥ |. We assume these dependences to be factorized and we assume for the k ⊥ and p ⊥ dependences a Gaussian form, with one free parameter which fixes the Gaussian width, The integrated PDFs, f q/p (x) and D h/q (z), can be taken from the available fits of the world data: in this analysis we will use the CTEQ6L set for the PDFs [32] and the DSS set for the fragmentation functions [33]. In general, the widths of the Gaussians could depend on x or z and might be different for different distributions: here, we first assume them to be constant and flavour independent and then perform further tests to check their sensitivity to flavour, x, z and Q 2 dependence. The constant Gaussian parameterisation, supported by a number of experimental evidences [11] as well as by dedicated lattice simulations [34], has the advantage that the intrinsic transverse momentum dependence of the cross section can be integrated out analytically. In fact, inserting Eqs. (5) and (6) into Eq. (1), one obtains where Notice that k 2 ⊥ and p 2 ⊥ will be taken as the free parameters of our fit. According to COMPASS [16] notation the differential hadron multiplicity is defined as: while HERMES [15] definition is where the index n denotes the kind of target. The Deep Inelastic Scattering (DIS) cross section has the usual leading order collinear expression, Inserting Eq. (1), (7) and (11) into Eq. (9) we have a simple explicit expression for the COMPASS and HERMES multiplicities: with P 2 T given in Eq. (8). Notice that, by integrating the above equation over P T , with its magnitude ranging from zero to infinity, one recovers the ratio of the usual leading order cross sections in terms of collinear PDFs and FFs. Its agreement with experimental data has been discussed, for instance, in Refs. [15] and [28].

III. RESULTS
As mentioned in the introduction, the most recent analyses of HERMES and COMPASS Collaborations provide (unintegrated) multivariate experimental data, presented in bins of x B = x, Q 2 , z h = z and P T . The HERMES multiplicities refer to identified hadron productions (π + , π − , K + , K − ) off proton and deuteron targets, and are presented in 6 bins of definite Q 2 and x B values, each for several different values of z h and P T , for a total of 2 660 data points. The selected events cover the kinematical region of Q 2 values between 1 and 10 GeV 2 and 0.023 ≤ x B ≤ 0.6, with a hadronic transverse momentum P T < 2 GeV and a fractional energy z h in the range 0.1 ≤ z h ≤ 0.9.
Instead, the COMPASS multiplicities refer to unidentified charged hadron production (h + and h − ) off a deuteron target ( 6 LiD), and are presented in 23 bins of definite Q 2 and x B values, each for several values of z h and P T , for a total of 18 624 data points. The Q 2 and z h regions covered by COMPASS are comparable to those explored by the HERMES experiment, while they span a region of smaller x B values, 0.0045 ≤ x B ≤ 0.12, and cover a wider P T region (reaching lower P T values). Moreover, the binning choices are very different and COMPASS statistics is much higher than that of HERMES.
For all these reasons, we consider the two data sets separately and perform individual fits.

A. Fit of the HERMES multiplicities
The first step of our analysis consists in using the simple Gaussian parameterisation of Eqs. (5) and (6) and the expression (12), to perform a two parameter fit of the HERMES multiplicities M h n (x, Q 2 , z, P T ). The values of the best fit parameters, the Gaussian widths k 2 ⊥ and p 2 ⊥ , will fix the TMD distribution and fragmentation functions respectively. We do not introduce any overall normalisation constant.
To make sure we work in the region of validity of our simple version of TMD factorisation, Eq. (12), we further restrict the kinematical range explored by the HERMES experiment. In fact, previous studies of the HERMES Collaboration [15] showed that the LO collinear SIDIS cross sections (obtained by integration of Eq. (12) over P T ), agree reasonably well with data only in regions of moderate values of z. The collinear distribution and fragmentation functions which perform best are the CTEQ6L PDF set [32] and the DSS [33] FFs, which we use here. We then consider two possible data selections: z < 0.7 and z < 0.6. Notice that these choices also avoid contaminations from exclusive hadronic production processes and large z re-summation effects [35]. We also fix the same minimum Q 2 as in the CTEQ6L analysis, Q 2 > 1.69 GeV 2 , that amounts to excluding the first two HERMES Q 2 bins. Low P T HERMES data show peculiar deviations from the Gaussian behaviour, which instead are not visible in the COMPASS and JLab [12,13] data: for this reason we prefer not to consider the lowest P T bin in order to explore the regions which exhibit the same kind of behaviour for all experiments. Finally, we apply an additional cut on large P T , requiring P T < 0.9 GeV, as multiplicities at large P T values fall in the domain of the onset of collinear perturbative QCD [36]. In the considered Q 2 range, this implies P T /Q < 0.7. Notice that recent analyses of the same experimental data [26,28] have adopted similar cuts.
Summarising, we limit the analysis of HERMES SIDIS data to the kinematical regions defined by: Moreover, in our fit we do not include the kaon production data points; in fact, the precision and accuracy of the kaon data sample, at present, do not help in constraining the values of the free parameters. When taken into account, the kaon data have little or no impact on the fit and are compatible with the assumption of the same Gaussian width as for pion production. This will be explicitly shown below by computing, using the parameters extracted from pion data, the kaon multiplicities and comparing them with the HERMES results. The above selections reduce the number of fitted HERMES data points to either 576 for z < 0.7, or 497 for z < 0.6. TABLE I: χ 2 values of our best fits, following Eqs. (12) and (8), of the experimental HERMES measurements of the SIDIS multiplicities M h n (x B , Q 2 , z h , PT ) for π + and π − production, off proton and deuteron targets. We show the total χ 2 dof and, separately, the χ 2 point for π + and π − data. CTEQ6 PDFs and DSS FFs are used. Notice that the errors quoted for the parameters are statistical errors only, and correspond to a 5% variation over the total minimum χ 2 . The details of the fits are presented in Table I, where we show the χ 2 per degree of freedom (χ 2 dof ), the χ 2 per number of points (χ 2 point ) for π + and π − production and the resulting values of the two free parameters of the fit, k 2 ⊥ and p 2 ⊥ , with some statistical errors, as explained below. It is worth noticing again that we do not have to use any overall normalisation constant as an extra free parameter; our computations agree well in magnitude with the experimental multiplicities, which are normalised to the collinear DIS cross section.

HERMES
Before drawing hasty conclusions on the numerical values of the parameters, some comments might be helpful.
• Our lowest value of χ 2 dof is obtained by using the kinematical cuts of Eq. (14) with z < 0.6, χ 2 dof = 1.69 for a total of 497 fitted pion data points. The corresponding widths of the Gaussians representing the k ⊥ and p ⊥ dependences of the distribution and fragmentation functions, are: However, if we relax the cut in z to z < 0.7, Eq. (13), then the total χ 2 of the fit becomes larger, χ 2 dof = 2.62, and the value of the extracted k 2 ⊥ Gaussian width significantly decreases while that of p 2 ⊥ increases, as shown in the second row of Table I. This large value of χ 2 reflects the fact that, at large values of z, P 2 T deviates from the assigned linear behaviour in z 2 . Morever, as we already pointed out, the large z region suffers from our lack of knowledge on the collinear fragmentation functions.
• The errors quoted for the free parameters of our fit are obtained from a ∆χ 2 corresponding to a 5% variation over the total minimum χ 2 : following Ref. [37], we relax the usual choice of ∆χ 2 = 1, corresponding to a purely statistical error, in order to include in the quoted errors other, major sources of uncertainty in our fit, which mainly originate from the inaccuracy in the determination of the fragmentation functions. We have checked that, indeed, other choices of collinear PDFs and FFs lead to such uncertainties. Moreover, in reading the errors, one should keep in mind that the parameters are strongly correlated.
The multiplicities obtained from our best fit parameters, with the kinematical cuts of Eq. (14), are compared with the HERMES measurements off a proton target in Figs. 1 and 2 and off a deuteron target in Figs. 3 and 4, separately for positive and negative pions. The shaded uncertainty bands are computed according to Ref. [38].
We have also performed a series of tests to study the effect of kaon data on the extraction. While the optimal parameters do not significantly change when including these data in the fit, the value of χ 2 dof reduces from 1.69 to 1.25, which could naively be interpreted as an improvement in the quality of the fit. However, this is just the result of the large error bands in the kaon subset. In fact, a fit of the kaon data alone would yield χ 2 dof = 0.64, which signals that the errors on these measurements are too large to allow a reliable extraction of the free parameters. This is shown very clearly in Figs. 5-8 where the kaon multiplicities, computed according to Eqs. (12) and (8) with the parameters of Eq. (15) -extracted from the HERMES measurements of pion production only -are compared with the HERMES data.
A careful look at the plots in Figs. 1-4 shows that the description of the HERMES measurements is indeed satisfactory: the Gaussian parameterisation embeds the crucial features of the data, both in shape and size, over a broad kinematical range. The resulting value of χ 2 dof , still a bit sizeable, is somehow expected, given the uncertainties on the collinear fragmentation functions: as stated before, the HERMES analysis [15] showed that the agreement, for the integrated multiplicities, between SIDIS data and collinear LO theoretical computations is not perfect, and that the currently available fragmentation function sets still need further refinements, especially at large z, and for π − production. In fact, including larger values of z in the fit sizeably increases the total χ 2 , as shown in the second line of Table I.
As the HERMES Montecarlo event generator, as well as many phenomenological models, propose a possible dependence of p 2 ⊥ on z, we have also attempted a fit with a z-dependent p 2 ⊥ = N z a (1 − z) b . However, it turns out that this parameterisation cannot be seriously tested by the data selection we have used for our reference fit; in fact, with the cuts of Eqs. (13) and (14), and in particular for the z < 0.6 range, it is quite hard for the best fitting procedure to find a proper convergence. Consequently, one obtains a and b parameters affected by huge statistical errors; this large uncertainties include the zero value and make the resulting parameters hardly significant. Moreover, the total χ 2 dof improves only marginally, down from 1.69 to 1.63.   (12) and (8), with the parameters of Eq. (15), are compared with HERMES measurements for π + SIDIS production off a proton target [15]. The shaded uncertainty bands correspond to a 5% variation of the total χ 2 .  (12) and (8), with the parameters of Eq. (15), are compared with HERMES measurements for π − SIDIS production off a proton target [15]. The shaded uncertainty bands correspond to a 5% variation of the total χ 2 .   (12) and (8), with the parameters of Eq. (15), are compared with HERMES measurements for π + SIDIS production off a deuteron target [15]. The shaded uncertainty bands correspond to a 5% variation of the total χ 2 .  (12) and (8), with the parameters of Eq. (15), are compared with HERMES measurements for π − SIDIS production off a deuteron target [15]. The shaded uncertainty bands correspond to a 5% variation of the total χ 2 .   (12) and (8), with the parameters of Eq. (15), are compared with HERMES measurements for K + SIDIS production off a proton target [15]. Notice that these data are not included in the fit; the shaded uncertainty bands correspond to a 5% variation of the total χ 2 obtained by fitting pion data only.  (12) and (8), with the parameters of Eq. (15), are compared with HERMES measurements for K − SIDIS production off a proton target [15]. Notice that these data are not included in the fit; the shaded uncertainty bands correspond to a 5% variation of the total χ 2 obtained by fitting pion data only.   (12) and (8), with the parameters of Eq. (15), are compared with HERMES measurements for K + SIDIS production off a deuteron target [15]. Notice that these data are not included in the fit; the shaded uncertainty bands correspond to a 5% variation of the total χ 2 obtained by fitting pion data only.  (12) and (8), with the parameters of Eq. (15), are compared with HERMES measurements for K − SIDIS production off a deuteron target [15]. Notice that these data are not included in the fit; the shaded uncertainty bands correspond to a 5% variation of the total χ 2 obtained by fitting pion data only.
To try and obtain more stringent constraints on this unknown z-dependence, we included also all the large-z bins, thus increasing the number of fitted data points from 497 to 657: in this case the minimisation procedure converges easily and the sizes of the errors on a and b become acceptable, excluding the zero values. However, the χ 2 we obtain, χ 2 dof = 2.92, becomes large, reflecting the large uncertainty in the collinear fragmentation functions at large z, and does not allow a reliable interpretation of the physical meaning of the extracted parameters. For this reason, we believe that these measurements are unable to fix the precise z-dependence of the Gaussian width p 2 ⊥ . A study of the possible dependence of k 2 ⊥ on x, similar to that of Ref. [28], shows no significant effect for the HERMES data. Notice that some x dependence of k 2 ⊥ will be further considered in Section IV, in the framework of the TMD evolution.
Similarly, no significant improvement is obtained by allowing flavour dependent k 2 ⊥ Gaussian widths. However, it is interesting to notice that, within a parameterisation in which we allow for two different free parameters for the favoured and disfavoured fragmentation Gaussian widths, HERMES data seem to prefer a configuration in which the favoured p 2 ⊥ fav is smaller than the disfavoured p 2 ⊥ disf by about 15%: we obtain p 2 ⊥ fav = 0.12 ± 0.01 and p 2 ⊥ disf = 0.14 ± 0.01, together with k 2 ⊥ = 0.59 ± 0.08, corresponding to a total χ 2 dof = 1.60. These results are in agreement with those obtained in Ref. [28]. We can therefore conclude that the HERMES multiplicity measurements, in the selected kinematical regions, can satisfactorily be described by modelling the intrinsic transverse momentum dependence with a very simple Gaussian shape with a constant width. Only a very slight indication for a flavour dependence of the fragmentation width is found, with the disfavoured Gaussian slightly wider than the favoured one. We also notice that, presently, the HERMES data do not show any sign of Q 2 evolution. The scale dependence will be discussed in Section IV.

B. Fit of the COMPASS multiplicities
We now consider the COMPASS SIDIS multiplicities of Ref. [16] and best fit them separately. We apply the same cuts as those used for the HERMES set of data, Eq. (13) and (14), with the same motivations, which reduce the number of fitted data points to 6 284 for z < 0.7 and 5 385 for z < 0.6. Our resulting χ 2 dof values are presented in Table II and turn out to be much larger (χ 2 dof = 9.81 and 8.54) than those obtained by fitting the HERMES multiplicities.
Taken at face value the parameters extracted from our best fit (first line of Table II), indicate that COMPASS data lead to values of k 2 ⊥ consistent, within errors, with those obtained in the HERMES fit, Eq. (15), while p 2 ⊥ seems to be slightly larger (considering the tiny error in the determination of this parameter). As in the HERMES fit, by widening the range in z the total χ 2 increases, although in this case the values of the extracted parameters are much more stable.
Our best fit curves (z < 0.6), for the COMPASS multiplicities M h ± D (x B , Q 2 , z h , P T ), are shown in Figs. 9 and 10, where it is clear that the quality of such fits appear to be rather poor.
Attempts of improving the description of the COMPASS data assuming more elaborated modelings of the Gaussian widths, with x and z dependence, do not help. We have also allowed for flavour dependent values of k 2 ⊥ (valence and sea) associated to two (favoured and disfavoured) fragmentation Gaussian widths p 2 ⊥ . However, even by adding a remarkably large number of free parameters, the gain in χ 2 is not satisfactory and the description of the data improves only marginally.
By carefully inspecting Figs. 9 and 10, and the χ 2 contributions of each individual bin, we realised that major improvements of our description of COMPASS data cannot actually be achieved by modifying the Gaussian widths, nor by making it more flexible. Our simple Gaussian model can actually reproduce the shape of the data, even over a large kinematical range; rather, the difficulties of the fit seem to reside in the normalisation of some of the bins, in particular those corresponding to large values of y.
To analyse this effect more closely, we performed the following study: we fitted the same sets of COMPASS data, allowing different normalisation constants N for each individual (x B , Q 2 ) bin (i.e. for each individual panel of Figs. 9 and 10) and we plotted all the obtained values versus x B , Q 2 and y. Although no particular x B or Q 2 trend could be detected, it was immediately evident that the obtained N coefficients showed a linear dependence on the kinematical variable y, as shown in Fig. 11. This suggested that a better description of COMPASS multiplicities could possibly be achieved by assuming a y-dependent normalisation parameter. Consequently, we re-performed a fit of the COMPASS data by adding to the multiplicities of Eq. (12) an overall multiplicative normalisation factor, linearly dependent on y, parameterised as which implies two additional parameters A and B, assumed to be universal and flavour independent. With this parameterisation the quality of our best fit improves very significantly, resulting in a total χ 2 dof of 3.42 (z < 0.6), corresponding to A = 1.06 ± 0.06 and B = −0.43 ± 0.14 and only very slightly different values of the Gaussian widths with respect to those previously obtained, Eq. (16), as shown in the third entry of Table II. A similar improvement is obtained for z < 0.7.
The results of our best fit, for positively and negatively charged hadronic production, are presented in Fig. 12 and 13 respectively. By comparing these figures with the analogous Figs. 9 and 10, in particular the plots on the left sides, corresponding to large y (and low x B ) bins, one can see that the huge gain in χ 2 is due to the fact that only with this second approach we can reproduce the normalisation of the data. Further comments on this issue will be made in the conclusions.
To complete our analysis, we finally performed a fit in which we allowed for two individual Gaussian widths for the favoured and disfavoured fragmentation functions, similarly to what was done for the HERMES measurements, assuming that the final hadrons are mostly pions. However, in this case, these two parameters turn out to be roughly the same. For z < 0.6, we find: k 2 ⊥ = 0.60 ± 0.15, p 2 ⊥ fav = 0.20±0.04, p 2 ⊥ disf = 0.20±0.05, for the Gaussian widths and A = 1.06±0.06, B = −0.43±0.14, for the normalisation constants. These parameters correspond to χ 2 dof = 3.42. In conclusion, we have found that COMPASS data show the need for an overall y-dependent normalisation; having fixed that, then the multiplicities appear to be in agreement with a Gaussian dependence,  12) and (8), with the parameters of Eq. (16), are compared with COMPASS measurements for h + SIDIS production off a deuteron target [16]. The shaded uncertainty bands correspond to a 5% variation of the total χ 2 .
although the resulting value of χ 2 dof remains rather large. Notice that this normalisation issue is not observed in the HERMES multiplicities and its origin, at present, cannot easily be explained and deserves further studies.
Some general comments on COMPASS results, inspired and guided by our grouping of the data in the panels of Figs. 9 and 10 and by the study presented in Fig. 11, could help to understand the origin of the large values of χ 2 dof . Let us consider, for example, the data in the different panels of the same row in Fig. 9. The multiplicity data grouped there have all very similar values of Q 2 and are separated in bins of z; one can notice, going from left to right, that data with very close value of Q 2 and z, still show a sharp x dependence. This can hardly be reproduced by Eq. (12), even considering eventual higher order corrections. Similar considerations apply to Fig. 10.
The large χ 2 which persists even in the case in which we correct with N y , is mainly due to some particular subsets of data, as one can see from Figs. 12 and 13 looking at the rightmost lower panels. These data, if compared with those in the panels to their immediate left (which have very similar values of the binned kinematical variables) show a sudden sharp change, which our smooth Gaussian parameterisation is unable to describe. Such a sharp change corresponds to the first, lowest y point, in Fig. 11.  12) and (8), with the parameters of Eq. (16), are compared with COMPASS measurements for h − SIDIS production off a deuteron target [16]. The shaded uncertainty bands correspond to a 5% variation of the total χ 2 .

IV. SCALE EVOLUTION
In this Section we perform a series of simple tests to check whether the SIDIS multiplicity measurements show any indication of transverse momentum scale evolution. So far, in our model, the Q 2 scale evolution has been taken into account only through the DGLAP evolution of the collinear distribution and fragmentation functions, Eqs. (5) and (6), while no Q 2 dependence is inserted in the transverse momentum distributions.
Multivariate measurements of SIDIS multiplicities offer a unique chance to study the scale evolution of their transverse momentum spectra. A complete description of this evolution should be done in the context of a full TMD evolution scheme, where much progress has recently been achieved [2,3,5,17]. However, we do not use here the complete expression of the evolution equations, which are quite complex and can differ according to the scheme used; rather, we exploit some simple modelling which embeds only those features of the scale evolution which are most relevant in our studies. In particular we expect that, in the kinematical ranges covered by the HERMES and COMPASS experiments, the leading contribution should be given by the non-perturbative, model dependent, part of the TMD evolution equation, which consists essentially in a scale dependence of the PDF and FF Gaussian widths.
For our SIDIS studies we choose a parameterization of k 2 ⊥ and p 2 ⊥ , Eqs. (5) and (6), inspired to that used for Drell-Yan and e + e − processes [39,40] in the context of the Collins-Soper-Sterman resummation scheme [41]: with g 1 , g 2 , g 3 , g 1 , and g 2 free parameters to be extracted by fitting the experimental data. Incidentally, the choice of the argument of ln(10 e x) is just a simple way to assign to the g 3 parameter a suitable normalisation, such that k 2 ⊥ = g 1 + g 3 at x = 0.1 and Q 2 = Q 2 0 . As already noticed, the SIDIS cross section is not sensitive to the individual contributions of k 2 ⊥ and p 2 ⊥ , but only to their linear combination, P 2 T , see Eqs. (7) and (8). Therefore, we take where the g 2 parameter used here replaces the sum of g 2 and g 2 in Eq. (18) and (19). Notice that this is the reason why the value of g 2 extracted by fitting SIDIS data, in general, cannot be directly applied to other processes, like Drell-Yan or e + e − scattering, where knowledge on the individual contribution of either the distribution or the fragmentation TMD Gaussian width is needed. For HERMES data we found no sensitivity to these parameters, confirming the good agreement of their measured multiplicities with the most simple version of our Gaussian model.
The description of COMPASS data is much more involved. First of all, we have to introduce the y-dependent normalisation, as in Eq. (17); then, we study the scale dependence of the P 2 T Gaussian width. We find that this width is not constant and can actually depend on Q 2 and x: by allowing the dependence of Eq. (20), the obtained χ 2 dof decreases from 3.42 to 2.69. However, we realised that a better improvement can be achieved by choosing, for P 2 T , the parametric form P 2 T = a 1 + a 2 ln(10 y) + z 2 [a 1 + a 2 ln(10 y)] , which amounts to replacing the dependence on Q 2 and x by a dependence on their ratio y = Q 2 /(x s).
In this case we obtain χ 2 dof = 2.00. A further decrease in the total χ 2 can be achieved by adding a square root term, g 3f √ y, in the Gaussian width of the fragmentation function: in this case we obtain χ 2 dof = 1.81, close to the value found for the HERMES data, indicating a significant improvement with respect to the results obtained without scale dependence.
Although such an ad-hoc y-dependence can only be considered as a phenomenological attempt to improve the fits, it is consistent with the W 2 dependence of P 2 T already found by the COMPASS Collaboration itself in their data, see Fig. 9 of Ref. [16], and by the EMC Collaboration, see Figs. 6 and 8 of Ref. [10].  17) are compared with the COMPASS measurements for h + SIDIS production off a deuteron target. The shaded uncertainty bands correspond to a 5% variation of the total χ 2 .

V. COMPARISON WITH PREVIOUS EXTRACTIONS AND OTHER EXPERIMENTS
The SIDIS multiplicity data used in our present fits result from the most recent analyses of the HER-MES and COMPASS Collaborations. They represent, so far, the only multivariate analyses available.
Additional measurements are provided by the early EMC results of Ref. [10] or by the more recent SIDIS studies of JLab CLAS [12] and HALL-C [13,29] Collaborations. As we will explain below, these data are not best suited for the extraction of the free parameters of our fit and we have not used them. However, it is worth and interesting to check whether or not the parameters extracted here are consistent with the available EMC and JLab measurements.
• The EMC Collaboration [10] measured P T -distributions in eleven different runs presented in one merged data set, averaging over four different beam energies, three different nuclear targets, without any identification of the final hadrons (not even their charges), and arranging the data in three different bins of z and several ranges of W 2 . In Ref. [9] we exploited these measurements, together with the EMC measurements [42] of the azimuthal dependence of the SIDIS cross section, for a preliminary study of the Gaussian widths of the unpolarised distribution and fragmentation functions. The values found there are slightly different from those we determine in the present fit. Fig. 14 shows the EMC multiplicities [10] as functions of P 2 T , for three bins of z, 0.1 < z < 0.2, 0.2 < z < 0.4 and 0.4 < z < 1.0, and of the invariant mass, W 2 < 90, 90 < W 2 < 150 and 150 < W 2 < 200 (in GeV 2 ). These data are compared with our predictions, computed at two different beam energies, E lab = 120 GeV and E lab = 280, using the values of k 2 ⊥ and p 2 ⊥ extracted from the best fit of the COMPASS multiplicities (see the second entry of Table II). The slope of our P 2 T distributions, which represents the Gaussian width P 2 T , qualitatively agrees with the EMC measurements. The merging of these measurements in one data set, makes any quantitative conclusion quite lax (as one can see by comparing, for instance, the upper and lower panels of Fig. 14 which refer to two very different energies).
• Jlab-HALL-C Collaboration [29] provides unpolarised cross sections, for SIDIS scattering off 1 H (proton) and 2 H (deuteron) targets and for π + and π − final production, as a function of P 2 T , at fixed values of x = 0.32, Q 2 = 3.2 GeV 2 and z = 0.55. At fixed z, these distributions are only sensitive to the total Gaussian width P 2 T = p 2 ⊥ + z 2 k 2 ⊥ , and not to the separate values of k 2 ⊥ and p 2 ⊥ . In Fig. 15 we compare these JLab measurements [29] with the P T distribution calculated in our scheme by using the parameters extracted from a best fit of the HERMES multiplicities, Eq. (15). We use these parameter values as the kinematical region explored at JLab is similar to that spanned in the HERMES experiment. An overall factor N = 1.5, common to all four data sets, has been introduced. Notice that P 2 T , which here represents the only parameter to which these data are sensitive, appears both as the distribution width and as a normalisation factor, see Eq. (12). At a qualitative level, we can conclude that the Gaussian model and the extracted values of k 2 ⊥ and p 2 ⊥ are not in conflict with the JLab cross section measurements of Ref. [29]. However, different normalizations could require different widths. The shaded uncertainty bands are computed  FIG. 14: The EMC multiplicities [10] are plotted as functions of P 2 T for three bins of z, 0.1 < z < 0.2, 0.2 < z < 0.4 and 0.4 < z < 1.0, and of invariant mass, W 2 < 90, 90 < W 2 < 150 and 150 < W 2 < 200. These data are compared to the predictions, computed at two different beam energies, E lab = 120 GeV (upper panel) and E lab = 280 (lower panel), by using the k 2 ⊥ and p 2 ⊥ extracted from the best fit of the COMPASS multiplicities (see the second entry of Table II) .
by propagating the error on the two free parameters, as explained in Appendix A. Notice that the data extend to a region (P 2 T < 0.04 GeV 2 ) excluded in our best fit analysis. • In Fig. 16 we show the CLAS extraction of ( H 1 +H 2 )(P T )/( H 1 +H 2 )(P T = 0), defined in Eqs. (1) of Ref. [12], as a function of P 2 T at the fixed values of x = 0.24 and z = 0.30; in our model, Eq. (12), this ratio corresponds to the actual Gaussian e −P 2 T / P 2 T . From the theoretical point of view, this quantity is more suitable for a comparison with our predictions, as it is unambiguously normalised. However, from the experimental point of view, its determination requires the extrapolation of ( H 1 + H 2 )(P T ) down to P T = 0. The plot shows that the line we obtain by using the Gaussian widths extracted from HERMES multiplicities is qualitatively in agreement with the data taken at Q 2 = 2.37, showing the right slope. Furthermore, CLAS finds quite a strong scale dependence (not shown here), as their slopes change considerably at Q 2 = 2.00 and Q 2 = 1.74, a dependence which we cannot account for in our scheme. Instead, differences among slopes at different values of Q 2 could be obtained by applying the appropriate kinematical cuts to the integration over k ⊥ , rather than integrating over the full k ⊥ range, as suggested, for example, in Ref. [43].

VI. CONCLUSIONS
We have considered a large amount of unpolarised SIDIS multiplicity data, recently made available by the HERMES [15] and COMPASS [16] Collaborations, and their P T dependence. Such a dependence, in selected kinematical regions, is believed to be generated by the intrinsic motion of quarks in a nucleon (information encoded in the TMD-PDFs), and of the hadron in the fragmentation process (encoded in  JLab-Hall C PT -distributions for π + and π − production from SIDIS scattering off a 1 H (proton) and 2 H (deuteron) targets, plotted as a function of P 2 T . Experimental data, corrected for vector meson production, are from Ref. [29]. The lines are the prediction of our model obtained by using the parameters extracted from a best fit of the HERMES multiplicities, see Table I. The shaded uncertainty bands are computed by propagating the error on the two free parameters, as explained in Appendix V.
the TMD-FFs). The separation of the data in different bins of x, Q 2 , z and P T allows to study the P T dependence in appropriate independent kinematical regions, obtaining safe information on the TMDs.
We have parameterised the unknown TMDs in a most simple way, Eqs. (5) and (6). Our results can be summarised as follows.
Our simple Gaussian parameterisation delivers a satisfactory description of the HERMES data points over large ranges of x, z, P T and Q 2 , selected according to Eqs. (13) and (14). These measurements are well described by a TMD Gaussian model with constant and flavour independent widths, k 2 ⊥ and p 2 ⊥ , which we extract as (the only two) free parameters of our fit. There is no need of an overall normalisation. HERMES multiplicities do not show any significant sensitivity to additional free parameters: the fits do not improve by introducing a z-dependence in the Gaussian widths of the TMD-FFs or by allowing a flavour dependence in the Gaussian widths of the TMD-PDFs. We only find a slight improvement in χ 2 by using different (constant) Gaussian widths in the TMD-FFs; the disfavoured fragmentation functions show a preference for a width slightly wider than that of the favoured fragmentation functions.
Fitting COMPASS data turns out to be more difficult. One would expect their high statistics and fine accuracy to provide strong constraints on our model. However, while the Gaussian shape of the P T dependence is qualitatively well reproduced, there are some unresolved issues with their relative overall normalisation. By performing a bin-by-bin analysis, we realised that different normalisation for π + and π − production from SIDIS scattering off proton and deuteron targets, are compared to the predictions obtained from our model by using the parameters extracted from the best fit of the HERMES multiplicities, see Table I. constants are required for different y-bins. In particular, we found that this N y normalisation factor is roughly 1 for very small y, while decreasing linearly with growing y: this implies almost a factor two difference between the largest and the smallest y bin, even within very close values of Q 2 and x, which would be very difficult to accommodate in a QCD driven scheme, even considering scale evolution. This can clearly be seen, for instance, by comparing the panels of each row in Fig. 9 going from left to right. We are presently unable to determine the origin of this effect, which indeed needs further investigation, both on the theoretical and experimental sides.
The COMPASS fit returns a p 2 ⊥ TMD-FF Gaussian width slightly larger than that extracted from the HERMES multiplicities, while it delivers similar k 2 ⊥ values. Comments on the scale evolution sensitivity of COMPASS multiplicitites will be presented in the next item.
Notice that this analysis has been performed on the 2004 run data, when the COMPASS detector was not yet completely set up and no RICH was installed for final hadron separation. Future analyses of more recent COMPASS data with hadron identification might help to clarify the situation.
In our fit, with the phenomenological parameterisation of Eqs. (5) and (6), the only dependence on Q 2 is included in the collinear part of the TMD, i.e. in the collinear PDF or FF factor. The width of the Gaussian, which gives the k ⊥ (p ⊥ ) dependence of the TMDs, does not include any scale dependence. Moreover, the Gaussians are normalised to one (if one integrates over the full k ⊥ range) and consequently our P T integrated cross section corresponds to the usual LO collinear SIDIS cross section.
As an alternative, while retaining this normalisation, we have tried new parameterisations which allow for a Q 2 and/or x-dependence of the Gaussian widths, Eq. (20), or even some y-dependence, Eq. (21). For the HERMES data we did not find any significative x or Q 2 dependence in the trasnverse momentum spectra. For the COMPASS data, instead, some improvement in the quality of the fit can actually be obtained. However, due to the unresolved issues discussed above, we cannot give a clear interpretation of this sensitivity and prefer not to draw, at this stage, any definite conclusion.
Several general considerations should not be forgotten. It is quite possible that the span in Q 2 of the available SIDIS data is not yet large enough to perform a safe analysis of TMD evolution based only on these data. Another related issue is that, always considering the SIDIS data set, the values of P T , while being safely low, are sometimes close to Q and corrections to the TMD factorisation scheme might be still relevant.
Only a combined analysis of SIDIS and Drell-Yan data, that would span in Q 2 from 1 to approximately 80 GeV 2 , would allow a comprehensive study. Such an analysis is beyond the scope of this paper and will be published elsewhere.
• Comparison with other measurements.
We have compared the multiplicities and P T -distributions obtained by using the parameters extracted from our best fit with EMC [10] and JLab [12,29] experimental data. These data are not best suited for fitting and only a qualitative comparison is possible. We found that EMC data are described reasonably well using the parameter values extracted from COMPASS multiplicities while the JLab data seem to be compatible with predictions based on the parameters extracted from HERMES.
Summarising, from this study we find that the Gaussian widths describing the k ⊥ distribution of the unpolarised TMD PDFs and the p ⊥ distribution of the unpolarize TMD FFs span respectively the approximate ranges Indeed, the actual values found here, as in Ref. [28], based on unpolarised multiplicity data, still have large intrinsic uncertainties. Once again, it is important to point out that multiplicities are sensitive to P 2 T only, and consequently they do not strictly allow a separate determination of k 2 ⊥ and p 2 ⊥ . Although the relation P T = z k ⊥ + p ⊥ is dictated by kinematics and is therefore model indipendent, our relation between p 2 ⊥ , k 2 ⊥ and P 2 T , Eq. (8), does depend on the parametrization chosen. While the Gaussian dependence of the TMDs is well in agreement with multiplicity data, a precise determination of the separate values of k 2 ⊥ and p 2 ⊥ would require the analysis of other data, like the azimuthal dependences: in the Cahn effect, for instance, we are sensitive to the ratio k 2 ⊥ / P 2 T [9]. The study of azimuthal dependences, for the new HERMES and COMPASS data, will be performed in a forthcoming paper.