Transverse-momentum-dependent Multiplicities of Charged Hadrons in Muon-Deuteron Deep Inelastic Scattering

A semi-inclusive measurement of charged hadron multiplicities in deep inelastic muon scattering off an isoscalar target was performed using data collected by the COMPASS Collaboration at CERN. The following kinematic domain is covered by the data: photon virtuality $Q^{2}>1$ (GeV/$c$)$^2$, invariant mass of the hadronic system $W>5$ GeV/$c^2$, Bjorken scaling variable in the range $0.003<x<0.4$, fraction of the virtual photon energy carried by the hadron in the range $0.2<z<0.8$, square of the hadron transverse momentum with respect to the virtual photon direction in the range 0.02 (GeV/$c)^2<P_{\rm{hT}}^{2}<3$ (GeV/$c$)$^2$. The multiplicities are presented as a function of $P_{\rm{hT}}^{2}$ in three-dimensional bins of $x$, $Q^2$, $z$ and compared to previous semi-inclusive measurements. We explore the small-$P_{\rm{hT}}^{2}$ region, i.e. $P_{\rm{hT}}^{2}<1$ (GeV/$c$)$^2$, where hadron transverse momenta are expected to arise from non-perturbative effects, and also the domain of larger $P_{\rm{hT}}^{2}$, where contributions from higher-order perturbative QCD are expected to dominate. The multiplicities are fitted using a single-exponential function at small $P_{\rm{hT}}^{2}$ to study the dependence of the average transverse momentum $\langle P_{\rm{hT}}^{2}\rangle$ on $x$, $Q^2$ and $z$. The power-law behaviour of the multiplicities at large $P_{\rm{hT}}^{2}$ is investigated using various functional forms. The fits describe the data reasonably well over the full measured range.


Introduction
A complete understanding of the three-dimensional parton structure of a fast moving nucleon requires the knowledge of the intrinsic motion of quarks in the plane transverse to the direction of motion, both in momentum and coordinate space. While the spatial distributions of quarks in the transverse plane are described by generalised parton distributions (GPDs), the momentum distributions of quarks in the transverse plane are described by transverse-momentum-dependent (TMD) parton distribution functions (PDFs). A precise knowledge of TMD-PDFs is found [1] to be crucial in the explanation of many singlespin effects observed in hard scattering reactions [2,3,4], in addition to the important role they play in spin-independent processes. In a similar way, transverse-momentum-dependent fragmentation functions (TMD-FFs) are crucial for the description of hard scattering reactions involving hadron production. Both PDFs and FFs are non-perturbative quantities that are assumed to be process-independent. The simplest examples are the spin-averaged TMD-PDF f q 1 (x, k T ) and the spin-averaged TMD-FF D h q (z, p h⊥ ), where x is the Bjorken scaling variable, k T is the quark intrinsic transverse momentum, z is the fractional energy of the final-state hadron, and p h⊥ is the transverse momentum of the final-state hadron relative to the direction of the fragmenting quark. After integration over k T and p h⊥ , the TMD-PDFs and TMD-FFs reduce to the standard spin-averaged collinear PDFs and FFs, where collinear means along the direction of the virtual photon. While the knowledge on the collinear PDFs and FFs is quite advanced, very little is presently known about the dependence of TMD-PDFs and TMD-FFs on k T and p h⊥ , as only sparse experimental data are available to date.
One of the most powerful tools to assess TMD-PDFs and TMD-FFs is the semi-inclusive measurement of deep inelastic scattering (SIDIS), N → hX, where one hadron is detected in coincidence with the scattered lepton in the final state. According to the QCD factorisation theorem [5,6] the deep inelastic scattering (DIS) process is considered to proceed via two independent sub-processes, i.e. the elementary QED process q → q is followed by the hadronisation of the struck quark. The outgoing hadrons provide information about the original transverse motion of the quark in the nucleon via their transverse momentum vector P hT . The latter is defined with respect to the virtual photon direction. The SIDIS cross section can be written as a convolution of a 'hard' scattering cross section, which is calculable in perturbative QCD (pQCD), with the non-perturbative TMD-PDFs and TMD-FFs. It depends on five kinematic variables. Two variables describe inclusive DIS, i.e. the negative square of the four-momentum transfer Q 2 = −q 2 and the Bjorken scaling variable x = −q 2 /(2P · q), where q and P denote the fourmomenta of the virtual photon and the nucleon, respectively. Three more variables describe the final-state hadrons, i.e. the fraction of the virtual photon energy that is carried by a hadron, z = (P · P h )/(P · q), the magnitude P hT of the transverse momentum of a hadron and its azimuthal angle φ in the system of virtual photon and nucleon. Here, P h denotes the four-momentum of the hadron. In the present analysis, the dependence on φ is disregarded. When integrating it over, the differential cross section for spinindependent SIDIS reads as follows in twist two "TMD factorisation scheme" [7,8]: with Here, y is the lepton energy fraction that is carried by the virtual photon and s is the centre-of-mass energy, which are related to x and Q 2 through Q 2 = xys. The hadron transverse momentum is related to k T and p h⊥ by P hT = zk T + p h⊥ [7]. An important consequence of the factorisation theorem is that the fragmentation function is independent of x, and the parton distribution function is independent of z, while both depend on Q 2 .
In addition to azimuthal asymmetries in spin-independent SIDIS [9], the most relevant experimental observable to investigate spin-averaged TMD-PDFs and TMD-FFs is the differential hadron multiplicity as a function of P 2 hT , which is defined in Eq. 3 below. 'Soft' non-perturbative processes are expected to generate relatively small values of P hT with an approximately Gaussian distribution in P hT [10]. Hard QCD processes are expected to generate large non-Gaussian tails for P hT >1 (GeV/c). They are expected to play an important role in the interpretation of the results reported here, which reach values of P 2 hT up to 3 (GeV/c) 2 .
Transverse-momentum-dependent distributions of charged hadrons in DIS were first measured by the EMC collaboration [11] at CERN, followed by measurements by ZEUS [12] and H1 [13,14] at HERA. These measurements only provided data in a limited dimensional space. Only new-generation experiments provided higher statistics, thereby opening the way to analyse and present the results in several dimensions simultaneously. Recent results were obtained by several fixed-target experiments using various targets and complementary energy regimes, i.e. HERMES [15] at DESY and COMPASS [16] at CERN.
The present paper reports on a new COMPASS measurement of transverse-momentum-dependent multiplicities of charged hadrons and extends the results of our earlier publication on transverse-momentumdependent distributions of charged hadrons [16]. The present measurement enlarges the kinematic coverage in x up to 0.4 instead of 0.12, in Q 2 up to 81 (GeV/c) 2 instead of 10 (GeV/c) 2 and in P 2 hT up to 3 (GeV/c) 2 instead of about 1 (GeV/c) 2 with significantly reduced systematic uncertainties on the normalisation of the P 2 hT -integrated multiplicities [16]. The data reported here represent the most precise results on differential charged hadron multiplicities available at this energy scale. This measurement is unique as its high statistics allows us to analyse the P 2 hT -dependence of charged-hadron multiplicities in four variables simultaneously.
The paper is organised as follows. Section 2 briefly describes the experimental apparatus. Details about the data analysis are given in Sec. 3. The measured charged-hadron multiplicities are presented and compared to previous measurements in Sec. 4. In Sec. 5 fits to the results are presented and discussed. The results are summarised in Sec. 6.

Experimental setup
The set-up of the COMPASS experiment is shortly described in this section. A more detailed description can be found in Ref. [17]. It is a fixed-target experiment, which uses the CERN Super Proton Synchrotron M2 beamline that is able to deliver high-energy hadron and muon beams. The data were collected in 2006 using a naturally polarised µ + beam of 160 GeV/c with a momentum spread of 5%. The intensity was 4 × 10 7 s −1 with a spill length of 4.8 s and a cycle time of 16.8 s. The momentum of each incoming muon was measured before the COMPASS experiment with a precision of 0.3%. The trajectory of each incoming muon was measured before the target in a set of silicon and scintillating fibre detectors. The muons were impinging on a longitudinally polarised solid-state target located inside a superconducting magnet. The target consisted of three cells that were located along the beam one after the other. It was filled with 6 LiD beads immersed in a liquid 3 He/ 4 He mixture. The admixtures of H, 3 He and 7 Li in the target lead to an effective excess of neutrons of about 0.2%. To first approximation, it can be regarded as an isoscalar deuteron target and will be referred to as such in the following. The polarisation of the middle cell (60 cm length) was opposite to that of the two outer cells (30 cm long each), and the polarisation was reversed once per day. In order to obtain spin-independent results, the target polarisation was averaged over by combining the data from all three target cells. Since the data taking in the two target polarisation states was well balanced and remaining polarisation-dependent effects are very small, this procedure ensures that for the data analysis the target can be considered as unpolarised.
The COMPASS two-stage spectrometer was designed to reconstruct scattered muons and produced hadrons in a wide range of momentum and polar angle, where the latter reaches up to 180 mrad. Particle tracking is performed by a variety of tracking detectors that are located before and after the two spectrometer magnets. The direction of the reconstructed tracks at the interaction point is determined with a precision of 0.2 mrad. The momentum resolution is 1.2% in the first spectrometer stage and 0.5% in the second one. The trigger is made by hodoscope systems supplemented by hadron calorimeters. Muons are identified downstream of hadron absorbers.
3 Multiplicity and data analysis 3.1 Multiplicity extraction The differential multiplicity M h for charged hadrons, where h denotes a long-lived charged hadron (π + , π − , K + , K − , p orp), is defined as the ratio between the differential semi-inclusive cross section d 4 σ h and the differential inclusive cross section d 2 σ DIS : Hadron multiplicities are measured in the four-dimensional (x, Q 2 , z, P 2 hT ) space. The bin limits in the four variables are presented in Table 1.
The data used in the present analysis were collected during six weeks in 2006. The data analysis comprises event and hadron selection, the correction for radiative effects, the determination of and the correction for the kinematic and geometric acceptance of the experimental set-up as well as for detector inefficiencies, detector resolutions and bin migration, and the correction for diffractive vector-meson production. Differential hadron multiplicities are evaluated as the ratio of hadron yields d 4 N h in every interval of (x, Q 2 , z, P 2 hT ) and the number of DIS events d 2 N DIS in every interval of (x, Q 2 ) corrected as described above : Here, η DIS and η h denote the correction factors accounting for radiative effects in the inclusive and in the semi-inclusive case, respectively, a h accounts for acceptance effects, and C DIS(h) denotes the correction factor accounting for the diffractive vector-meson contribution in the case of an inclusive (semi-inclusive) measurement. The (x, Q 2 ) dependence is omitted for simplicity as it enters all terms. All corrections are evaluated in the four-dimensional (x, Q 2 , z, P 2 hT ) bins except C DIS , η DIS and η h , which are evaluated only in bins of x and Q 2 . Further kinematic dependences of η h upon z and P 2 hT will be discussed in Sec. 3.2.

Event and hadron selection
The present analysis uses events taken with 'inclusive triggers', i.e. the trigger decision is based on scattered muons only. The selected events are required to have a reconstructed interaction vertex associated with an incident and a scattered muon track. This vertex has to lie inside a fiducial target volume. The incident muon energy is constrained to the range from 140 GeV to 180 GeV. In addition to the kinematic constraints given by the spectrometer acceptance, the selected events are required to have Q 2 >1 (GeV/c) 2 and W > 5 GeV/c 2 . These requirements select the DIS regime and exclude the nucleon resonance region. The relative virtual-photon energy is constrained to the range 0.1 < y < 0.9 to exclude kinematic regions where the momentum resolution degrades and radiative effects are most pronounced [18]. In the range 0.003 < x < 0.4, the total number of inclusive DIS events is 13 × 10 6 , which corresponds to an integrated luminosity of 0.54 fb −1 . The (x, Q 2 ) distribution of this selected 'DIS sample' is shown in Fig. 1, where a strong correlation between x and Q 2 is observed as expected in fixed-target experiments. For a selected DIS event, all reconstructed tracks associated with the primary interaction vertex are considered. Hadron tracks must be detected in detectors located before and after the magnet in the first stage of the spectrometer. The fraction of the virtual-photon energy transferred to a final-state hadron is constrained to 0.2 < z < 0.8. The lower limit excludes the target fragmentation region, while the upper one removes muons wrongly identified as hadrons and excludes the region with larger contributions from diffractive ρ 0 production. This selection yields the 'hadron sample' with a total of 4.3 × 10 6 and 3.4 × 10 6 positively and negatively charged hadrons, respectively.
The corrections for QED higher-order effects are applied on an event-by-event basis taking into account the target composition. They are computed as a function of x and y according to the scheme described in Ref. [19]. For the hadron yields, the correction is calculated by excluding the elastic and quasi-elastic tails. The correction factors η h and η DIS are evaluated in bins of x and Q 2 . They are found to be smaller than 12 % for x < 0.01 and are smaller than 5% elsewhere. An attempt to evaluate the smearing due to radiative effects as a function of z and P 2 hT was performed using a Monte Carlo (MC) simulation, where radiative effects were simulated using the RADGEN generator [20]. A possible impact on the (z, P 2 hT ) dependence of the results due to radiative effects is accounted for in the systematic uncertainties of the P 2 hT -dependence of the multiplicities.

Acceptance correction
The hadron multiplicities must be corrected for geometric and kinematic acceptances of the experimental set-up as well as for detector inefficiencies and resolutions, and for bin migration. The correction for a possible misidentification of electrons as hadrons is included in the acceptance correction. The full correction factor is evaluated using a MC simulation of the muon-deuteron deep inelastic scattering processes. Events are generated using the LEPTO generator [21], where the parton hadronisation mechanism is simulated using the JETSET package [22] with the tuning from Ref. [23]. Secondary hadron interactions are simulated using the FLUKA package [24]. The experimental set-up is simulated using the GEANT3 toolkits [25] and the MC data are reconstructed using the same software that was used for the experimental data [17]. The kinematic distributions of the experimental data are fairly well reproduced by the MC simulation.
In order to minimise a possible dependence on the physics generator used in the simulation and to exclude kinematic regions with large acceptance corrections, a four-dimensional evaluation of the acceptance correction factor a h is performed in narrow kinematic bins. In each (x r , Q 2 r ) kinematic bin, where r denotes the reconstructed values of the variables, the acceptance correction is calculated as the ratio of reconstructed (d 2 N h r ) and generated (d 2 N h g ) hadron yields, where both are evaluated using the simulated DIS sample after reconstruction : An advantage of this definition is that the correction for muon acceptance cancels as it enters both numerator and denominator. The generated values of kinematic variables are used for the generated particles and the reconstructed values of kinematic variables are used for the reconstructed particles. All reconstructed MC events and particles are subject to the same kinematic and geometric selection criteria as the data, while the generated ones are subject to kinematic requirements only. The acceptance correction factor exhibits an almost flat behaviour as a function of z and P 2 hT in most (x, Q 2 ) bins, except at high x for P 2 hT > 1 (GeV/c) 2 , where it remains larger than 0.4. Elsewhere, its average value is above and close to 0.6 for P 2 hT > 0.5 (GeV/c) 2 and is less than or equal to 0.6 for P 2 hT < 0.5 (GeV/c) 2 . As an example, Fig. 2 shows the acceptance as a function of P 2 hT for positively charged hadrons. The two panels show the two z bins between 0.4 and 0.8, with two bins in (x, Q 2 ) in each case. The acceptance correction factors for positively and negatively charged hadrons are found to be very similar, with differences on the level of 0.02-0.04.

Diffractive vector meson contribution
The final-state hadron(s) selected as described above may also originate from diffractive production of vector mesons (ρ 0 , φ , ω) that decay into lighter hadrons (π, K, p) [4,26,27]. This process, which can be described by the fluctuation of the virtual photon into a vector meson that subsequently interacts diffractively with the nucleon through multiple gluon exchange, is different from the interaction of the virtual photon with a single quark in the DIS process. The fraction of selected final-state hadrons originating from diffractive vector-meson decays and their contribution to the SIDIS yields are estimated in each kinematic bin using two Monte Carlo simulations. The first one uses the LEPTO generator to simulate SIDIS events, and the other one uses the HEPGEN generator [28] to simulate diffractively produced ρ 0 6 The COMPASS Collaboration and φ events. Further channels, which are characterised by smaller cross sections, are not taken into account. Events with diffractive dissociation of the target nucleon represent about 25% of those with the nucleon staying intact and are also simulated. The simulation of these events includes nuclear effects, i.e. coherent production and nuclear absorption as described in Ref. [28]. The contribution of pions originating from ρ 0 decay to the hadron sample increases with z, and reaches up to 40 − 50% for z close to 1. For kaons, the contribution from φ decay is concentrated in the z range 0.4 − 0.6, where it reaches up to 15%. The correction factors are separately evaluated for the DIS sample and the hadron sample : Here, f VM DIS denotes the fraction of diffractively produced vector-mesons present in the DIS sample, while f ρ 0 π and f φ K denote the fraction of ρ 0 and φ decay products in the hadron sample, respectively. The fraction of pions, kaons, and protons in the latter sample, which is denoted by F π, K, p , amounts to about 75%, 20% and 5%. The fractions F i (i = π, K, p) and f VM i (i = π, K and VM = ρ 0 , φ ) are evaluated as functions of x, Q 2 , z and P 2 hT . In the following, the general behaviour of some of the above discussed correction factors is illustrated. The correction factor to account for diffractive ρ 0 production in the DIS yield is shown in Fig. 3 as a function of x in the five Q 2 bins. It reaches a maximum value of about 4% in the lowest Q 2 bin. The correction factor for the contribution of diffractively produced ρ 0 mesons to the pion sample, (1 − f ρ 0 π ), is shown in Fig. 4(a) as a function of P 2 hT in the four z bins for the lowest Q 2 bin, where it has the largest value. It reaches a maximum value of about 25% for P 2 hT ∼ 0.12 (GeV/c) 2 in the highest z bins, i.e. 0.6 < z < 0.8, and decreases to few percent at small z. The correction factor for the contribution of diffractively produced φ mesons to the kaon sample, (1 − f φ K ), is shown in Fig. 4(b). In this case, the maximum correction of about 35% is reached at very small P 2 hT in the middle z bin, i.e. 0.4 < z < 0.6. Fig. 3: Correction factor to the DIS yield due to diffractive ρ 0 production as a function of x in the five Q 2 bins.

Systematic uncertainties
The dominant contributions to the systematic uncertainties originate from the uncertainties on the determination of the acceptance correction factor and those of the diffractive vector-meson contribution. The uncertainty on the acceptance calculation is evaluated by varying in the MC simulation both the PDF set and the JETSET parameters describing the hadronisation mechanism. The acceptance correction is estimated for each MC sample and the largest deviation with respect to the values obtained using the MC simulation described in section 3.3 is quoted as a systematic uncertainty. The validity of the correction for the electron contamination is confirmed by comparing the simulated and measured electron distributions for momenta below 8 GeV/c, where electrons are identified using the RICH detector. In order to check a possible dependence on the target cell, in which the event vertex is located, the multiplicities are independently measured from the three target cells. Results from upstream and downstream target cells agree within 2-3%, while the agreement is better than 1% with the middle target cell. These differences are well covered in the acceptance correction uncertainty. A total uncertainty of 5% is estimated for the multiplicities.
The cross section for exclusive production of ρ 0 calculated in HEPGEN is normalised to the phenomenological model of Ref. [29]. The theoretical uncertainty on the predicted cross section in a kinematic region close to COMPASS kinematics amounts to about 30%. This results in an uncertainty on the diffractive vector-meson correction factor, which amounts up to 5 − 6% mainly at small values of x, Q 2 and P 2 hT , and large values of z. Nuclear effects may be caused by the presence of 3 He/ 4 He and 6 Li in the target. The EMC Collaboration has studied in detail such nuclear effects in a similar kinematic range using carbon, copper and tin targets [11]. A z-dependent decrease of 5% was observed for the multiplicities obtained using copper compared to the ones obtained using deuterium. While the effect was larger for tin, no such effect was found for carbon, so that possible nuclear effects in the present experiment are expected to be very small and are hence neglected. When comparing the results obtained from the data taken in six different weeks, no difference is observed. Their numerical values are available on HepData [30] with and without correction for diffractive vectormeson production. It should be noted that a few (x, Q 2 ) kinematic bins are discarded in the lowest (Fig. 5) and the highest (Fig. 8) bins of z because of low statistical precision as well as large acceptance correction factors (Sec. 3.3). The average values of x and Q 2 in the various kinematic bins are evaluated using the DIS sample and are given in Table 2. The results obtained by integrating the multiplicities presented here over P 2 hT are in very good agreement with those of Ref. [26], where the multiplicities of charged pions are measured as a function of z in a restricted momentum range based on an independent analysis of the same data. Multiplicities are larger for positively than for negatively charged hadrons. This difference significantly increases as x increases and shows a weak variation with Q 2 . It is observed to also depend on z and it increases in the range of large z, i.e. z > 0.4, which confirms the observations made in Ref. [26]. Besides their magnitude, the P hT dependence of the multiplicities shows a significant variation with x at fixed Q 2 (as well as with Q 2 at fixed x) for any interval of z. These observations are separately illustrated in Figs. 9 and 10 and discussed in detail in the following.
The comparison between the multiplicities of positively and negatively charged hadron is illustrated as a function of x and Q 2 in Fig. 9 for z = 0.35. On the top row, h + and h − multiplicities are presented at Q 2 1.3 (GeV/c) 2 in the smallest and the largest x bins with average values x = 0.0062 and x = 0.039, respectively. In the right column, h + and h − multiplicities are similarly presented at x 0.04 in the smallest and the largest Q 2 bins with average values 1.4 (GeV/c) 2 and 8.3 (GeV/c) 2 , respectively. At fixed Q 2 , the ratio of h + to h − multiplicities ranges from about 1 in the first x bin to about 1.3 in the last x bin. This increase as a function of x confirms the expectation from valence u-quark dominance, i.e. the dominance of scattering off u-quarks. At fixed x, the ratio of h + to h − multiplicities decreases from 1.3 in the first Q 2 bin to about 1.2 in the last Q 2 bin. While no significant difference is observed in the P 2 hT -dependence of h + and h − multiplicities, the P 2 hT -dependence of the multiplicities is observed to flatten at large values of P 2 hT , where contributions from higher-order QCD processes like QCD Compton and photon-gluon fusion (PGF) are expected to dominate. The data suggest that flattening occurs both as Q 2 increases (at fixed x) and when x decreases (at fixed Q 2 ).
In Figure 10, the comparison between h + and h − multiplicities is illustrated as a function of z. The multiplicities are presented as a function of P 2 hT in the four z intervals in a given (x, Q 2 ) bin with average values x = 0.149 and Q 2 = 9.78 (GeV/c) 2 . A high x bin is chosen, where the difference in the magnitude of the multiplicities is most recognisable. The ratio of h + to h − ranges from about 1.1 in the first z bin to about 2 in the last z bin, reflecting the fact that part of the negative hadrons (K − andp) can not be produced by the favoured fragmentation of a nucleon valence quark, which enhances the expected flavour dependence of TMD-FFs. Another feature of the data is the variation of the P 2 hT -dependence with increasing z for both small and large values of P 2 hT . In particular, the data show a tendency to flatten at large P hT as z decreases, which emphases a significant z-dependence of the hadron transverse momentum with respect to the transverse momentum of the fragmenting quark, p ⊥ .
Another intriguing effect is observed in the kinematic domain 1 (GeV/c) 2 < Q 2 < 1.7 (GeV/c) 2 and 0.6 < z < 0.8, in the range of small P 2 hT . Charged hadron multiplicities do not exhibit an exponential form in P 2 hT in this kinematic region and show an unexpected flat dependence at very small values of P 2 hT . This effect is also present in the earlier published [16] distributions of charged hadrons as a function of P 2 hT . It is illustrated in Fig. 11, which shows the multiplicity of positive hadrons as a function of P 2 hT up to 0.8 (GeV/c) 2 at Q 2 = 1.25 (GeV/c) 2 and x = 0.0062 (left-hand side) and at Q 2 = 4.52 (GeV/c) 2 and x = 0.043 (right-hand side). It should be noted that this particular kinematic region suffers from the highest contribution of the ρ 0 decay products to the charged hadron sample (Fig. 4, blue curve) evaluated using the MC simulation. This effect is further discussed in Sec. 5.
The multiplicities shown in Figs 5-11 agree with the previous measurement of hadron distributions performed by COMPASS [16]. However, as mentioned in Sec. 1, this measurement considerably extends the kinematic range and reduces the statistical and systematic uncertainties, in particular the uncertainties on the normalisation of the P 2 hT -integrated multiplicities.

Comparison with other measurements
The multiplicities presented above are compared in Figs. 12-14 to results from previous semi-inclusive measurements in similar kinematic regions. The experiments are compared in Table 3. In order to compare the present COMPASS results on TMD hadron multiplicities with the corresponding ones by EMC [11], our data sample is reanalysed in bins of z and W 2 according to the binning given in Ref. [11]. The EMC measurements are performed in slightly different kinematic ranges in Q 2 and y, as shown in Tab. 3. While for the measurement described in this paper a deuteron target was used, EMC used proton and deuteron targets and also four different beam energies, which led to four different kinematic ranges. The comparison shown in Fig. 12, where the sum of h + and h − multiplicities is presented as a function of P 2 hT in four W 2 bins in the range 0.2 < z < 0.4, demonstrates good agreement between COMPASS and EMC results. According to the study in Ref. [10], the P 2 hT -dependence of the EMC data could be explained in the simple collinear parton model up to 8 (GeV/c) 2 in P 2 hT . In Figure 13, the multiplicities of positively charged hadrons are compared in the four bins of z to the multiplicities of positively charged pions measured by the HERMES Collaboration [15], where both were corrected for diffractive vector-meson contribution. The measurements by HERMES cover the kinematic range Q 2 > 1 (GeV/c) 2 and 0.023 < x < 0.6. For this comparison, the COMPASS h + multiplicities are integrated over x in the closest possible range 0.02 < x < 0.4 and also over Q 2 . It should be noted that the two experiments cover different ranges in Q 2 . While the highest Q 2 value reached by HERMES is 15 (GeV/c) 2 , COMPASS reaches 81 (GeV/c) 2 . Despite this difference, a reasonable agreement in the magnitude of the measured multiplicities is found for z < 0.6 and small P 2 hT . Most likely due to the differences in kinematic coverage, the agreement between the two sets is rather modest, and the data sets exhibit different dependences upon P 2 hT . In addition, a dip is observed in the HERMES data at very small transverse momenta, i.e. P 2 hT ∼ 0.05 (GeV/c) 2 . This dip, which is not observed in the shown Q 2integrated distribution, appears to be very similar to the trend shown in Fig. 11 by the COMPASS data at low Q 2 .
In Figure 14, the h + multiplicities are compared to the π + semi-inclusive cross section measured by the E00-18 experiment [31] at Jefferson Lab. The measurement by the E00-18 was performed at z = 0.55 and x = 0.32 in the range 2 (GeV/c) 2 < Q 2 < 4 (GeV/c) 2 . The COMPASS results are given at similar (x, z) values, i.e. z = 0.5, x = 0.3, and span the range 7 (GeV/c) 2 < Q 2 <16 (GeV/c) 2 . Similar to the case of the comparison of COMPASS and HERMES data shown in Fig. 13, here the observed different P hT -dependence could be due to the different Q 2 values of the two measurements.   The P hT -dependence of the cross section for semi-inclusive measurements of hadron leptoproduction was empirically reasonably well described by a Gaussian parameterisation for the k T -and p h⊥ -dependence of TMD-PDFs and TMD-FFs in the range of small P hT , i.e. P hT < 1 (GeV/c). This Gaussian parameterisation leads to a P 2 hT -dependence of the multiplicities of the form : where the normalisation coefficient N and the average transverse momentum P 2 hT , i.e. the absolute value of the inverse slope of the exponent in Eq. 8, are functions of x, Q 2 and z.
A fairly good description of SIDIS [1] data was reached with the Gaussian parameterisation without considering neither the z nor the quark flavour dependence of TMD-FFs. Recent semi-inclusive measurements of transverse-momentum-dependent hadron multiplicities [15] and distributions [16] aimed at an extraction of both k 2 ⊥ and p 2 ⊥ . These two observables, however, were found to be too strongly anti-correlated to be disentangled [8,32,33]. In order to extract them, a combined analysis of both the differential transverse-momentum-dependent hadron multiplicities and the spin-independent azimuthal asymmetries in SIDIS may be required. In the following we will discuss separately fits in the region of small P hT and in the full range of P 2 hT accessible by COMPASS, i.e. 0.02 (GeV/c) 2 < P 2 hT < 3 (GeV/c) 2 . The hadron multiplicities presented in Figs. 5-8 are fitted in each (x, Q 2 , z) kinematic bin in the range 0.02 (GeV/c) 2 < P 2 hT < 0.72 (GeV/c) 2 using the single-exponential function given in Eq. 8. Using only statistical uncertainties in the fit, reasonable values of χ 2 per degree of freedom (χ 2 dof ) are obtained in all (x, Q 2 ) bins, except for low values of Q 2 and small values of z, i.e. z < 0.3, where the χ 2 dof values are significantly larger than 3 in most of the x bins. Including the systematic uncertainties in the fit by adding them in quadrature to the statistical ones significantly improves the values of χ 2 dof , whereas the fitted parameters remain unchanged. The z 2 -dependence of P 2 hT obtained from the fits is shown in Fig. 15 for h + in the five Q 2 bins available in a given x bin. A non-linear dependence of P 2 hT on z 2 is observed in the range of small x and Q 2 , in contrast to the range of large x and Q 2 where it becomes linear. In addition, P 2 hT significantly increases with Q 2 at fixed x and z, especially at high z. The h + multiplicities have larger values of P 2 hT than the h − ones at large z, while no significant difference is observed at small z. This conclusion confirms the one made in our previous publication [16], where a detailed study of the kinematic dependence of P 2 hT was presented and discussed. As mentioned earlier in Sec. 4.2, the kinematic region of small Q 2 and large z, i.e. Q 2 < 1.7 (GeV/c) 2 and 0.6 < z < 0.8, shows an intriguing effect in the range of small P 2 hT . As can be seen from Fig. 11, in this range h + and h − multiplicities do not exhibit an exponential form in P 2 hT and show an unexpected flat dependence at very small values of P 2 hT . Figure 16(a) shows the multiplicity of positively charged hadrons as a function of P 2 hT up to 0.8 (GeV/c) 2 at Q 2 = 1.25 (GeV/c) 2 and x = 0.006. While a single-exponential function reasonably describes the P 2 hT -dependence for 0.3 < z < 0.4, the experimental data clearly deviate from this functional form as z increases, with χ 2 dof values increasing from 1.8 in the smallest z bin to 4.6 in the largest one. As an example, Fig. 16(b) shows h + multiplicities at larger Q 2 , i.e. Q 2 = 4.65 (GeV/c) 2 and x = 0.075, where the single-exponential function fits the data well in all z bins.
The measured charged-hadron multiplicities show that in the range of small P hT , i.e. for P hT < 1 (GeV/c) 2 , the simple parameterisation using a single-exponential function describes the P 2 hT -dependence of the re-  Fig. 15: Average transverse momentum P 2 hT , as obtained from the fit of h + multiplicities using the single Gaussian parameterisation, shown as a function of z 2 . The eight panels correspond to the eight x-bins as indicated, where in each panel data points from all five Q 2 bins are shown. Error bars denote to statistical uncertainties. sults quite well for not too large values of Q 2 . For increasing Q 2 , the P 2 hT -dependence of the multiplicities changes as can be seen in Fig. 9. A more complex parameterisation appears to be necessary to fit the data, as shown in Ref. [35] .

The full measured P hT range
Up to now, only one study [10] was performed to describe the full range in P hT using a Gaussian parameterisation for the k T -and p h⊥ -dependence of TMD-PDFs and TMD-FFs in the range P hT < 1 GeV/c, and calculating pQCD higher order collinear contributions in the range P hT > 1 GeV/c. A reasonable description of semi-inclusive hadron multiplicities and cross sections measured by the EMC [11] and ZEUS [12] Collaborations, respectively, was achieved. Below, we attempt to describe the P 2 hTdependence of the above presented charged-hadron multiplicities over the full P hT -range explored by COMPASS, i.e. 0.02 (GeV/c) 2 < P 2 hT < 3 (GeV/c) 2 , using the following two parameterisations : The first function (F 1 ) is defined as the sum of two single-exponential functions (Eq. 9). While N 1 and N 1 denote the normalisation coefficients, α 1 and α 1 denote the inverse slope coefficients of the first and the second exponential function, respectively. All coefficients depend on x, Q 2 and z. Figure 17 shows in a typical (x, Q 2 , z) bin the multiplicities of positively charged hadrons as a function of P 2 hT fitted using F 1 . As described above for Ref. [10], the two exponential functions in our parameterisation F 1 can be attributed to two completely different underlying physics mechanisms that overlap in the region P 2 hT 1 (GeV/c) 2 . Figure 18 shows, as an example, multiplicities of positively charged hadrons as a function of P 2 hT , measured at Q 2 ∼ 1.25 (GeV/c) 2 for two bins of x with average values x = 0.006 and x = 0.016, in the four z bins. Only statistical uncertainties are shown and used in the fit. Values of χ 2 dof of about 1 are obtained in all (x, Q 2 , z) bins, except for a few (6 out of 81) bins, where where values as small as 0.52 and as large as 2.52 are obtained. The normalisation coefficients N 1 and N 1 are found to have a strong variation with x and z and a rather weak variation with Q 2 , reflecting the (x, Q 2 ) dependence of collinear PDFs and the z-dependence of collinear FFs. The inverse slope α 1 has an average value of about 0.23 (GeV/c) 2 for Q 2 < 3 (GeV/c) 2 and about 0.28 (GeV/c) 2 for larger values of Q 2 . Its dependence on z 2 is discussed below using Fig. 19. The inverse slope α 1 has an average value of about 0.6 (GeV/c) 2 and shows a rather weak variation with x and Q 2 .
The so-called Tsallis function F 2 [36], see Eq. 10, describes the two different kinds of power-law behaviour in the two regions of P hT through a single function. The advantage of this function is that it provides both the inverse slope parameter T that characterises the small-P hT range and the exponent 1/(1−q) that parameterises the power-law tail at large P hT . The charged-hadron multiplicities (Figs. [5][6][7][8] are fitted in each (x, Q 2 , z) bin using only statistical uncertainties. Reasonable values of χ 2 dof are obtained in most bins except for 11 out of 81 bins where they are larger than 2 reaching up to 3.65. The exponent parameter q has an average value of about 1.2. The exponent 1/(1 − q) strongly depends on x with a weaker dependence on z and no variation with Q 2 , while its z-dependence is observed to increase with x. The inverse slope parameter T ranges between 0.15 (GeV/c) 2 and 0.4 (GeV/c) 2 and shows a significant non-linear dependence on z 2 over the full z-range.
The inverse slopes P 2 hT , α 1 and T , which were obtained using the fitting functions given in equations 8, 9 and 10 respectively, are presented and compared to each other as a function of z 2 in (x, Q 2 ) bins in Fig. 19. A weak non-linear dependence on z 2 is observed at small x and Q 2 , which becomes more pronounced at larger values of Q 2 . The inverse slope T reproduces the same z 2 -dependence as that of P 2 hT described in Fig. 8. It is observed to be in fair agreement with α 1 except for z > 0.6.
A comparison between the data and the fit function F 1 is shown in Fig. 20(a) in a typical kinematic bin with Q 2 = 2.12 (GeV/c) 2 and x = 0.011. The upper panel shows the multiplicities of positive hadrons as a function of P 2 hT and the corresponding fit function, and the lower panel shows the ratio between the data and the fit. A comparison between the two fitting functions F 1 and F 2 is shown in Fig. 20(b) for the same (x, Q 2 , z) bin. The P 2 hT -dependence of h + multiplicities is equally well described by the two functions F 1 and F 2 , as can be seen from the ratio in Fig. 20(b). The same agreement is obtained for negatively charged hadrons.   Fig. 19: Average transverse momentum obtained from the fit of h + multiplicities using the three fit functions given in Eqs. 8, 9, 10: P 2 hT , α 1 and T as a function of z 2 in (x, Q 2 ) bins. Multiplicities of positively charged hadron as a function of P 2 hT for Q 2 = 2.12 (GeV/c) 2 , x = 0.011 and z = 0.35. The black dotted curve represents the first exponential function of Eq. 9, the blue dashed curved represents the second exponential function of Eq. 9, and the red curve represents the sum. Only statistical uncertainties are shown and used in the fit. Lower panel: The ratio of the experimental points to the fit as a function of P 2 hT . (b) Comparison between the fits obtained using F 1 (Eq. 9) and F 2 (Eq. 10) in the same kinematic bin as in (a).

Summary
We have measured differential multiplicities of charge-separated hadrons in semi-inclusive measurements using muons of 160 GeV/c impinging on an isoscalar (deuteron) target. Using a high-statistics data set collected in 2006, the measurement covers a wide kinematic domain of Q 2 > 1 (GeV/c) 2 , W > 5 GeV/c 2 , 0.003 < x < 0.4, 0.2 < z < 0.8 and 0.02 (GeV/c) 2 < P 2 hT < 3 (GeV/c) 2 . The results are presented as a function of the square of the hadron transverse momentum P 2 hT in three-dimensional bins of x, Q 2 and z, which leads to a total of 4918 experimental data points. The numerical values are available on HepData [30] with and without subtraction of the estimated contribution of diffractive vector-meson production in SIDIS.
The h + multiplicities are only slightly larger than the h − ones in most of the bins, while for large x and z this difference increases. No significant difference between h + and h − is observed in the shape of the P 2 hT -dependence of the multiplicity. Both h + and h − multiplicities are observed to flatten at very small values of P 2 hT in the kinematic region of low x and Q 2 and large z, where contributions from diffractive vector-meson production are the highest. Our results are compared to earlier measurements of hadron multiplicities and cross sections by EMC, HERMES and JLab. Good agreement was found with EMC for W 2 < 150 (GeV/c) 2 , although the EMC data were collected at different beam energies and with different targets. In order to compare with HERMES, we have integrated our multiplicities over the phase space that is common to both experiments. While reasonable agreement is obtained at small z and P 2 hT , differences are observed for large z and P hT where neither magnitudes nor P 2 hT -dependences agree. The π + semi-inclusive cross section measured by the E00-18 experiment at JLab shows fair agreement with COMPASS h + multiplicities, albeit with some discrepancy in the P hT -dependence that might be explained by the difference in the kinematic ranges of the measurements.
In the range of small-P 2 hT , i.e. P 2 hT < 1 (GeV/c) 2 , the measured multiplicities were successfully fitted using a single Gaussian parameterisation. A non-linear z 2 -dependence of the average transverse momentum is observed in the range of small x and Q 2 , which confirms the conclusions of Ref. [16], while it is almost linear for large values of x and Q 2 . In order to fit the multiplicities over the full P hT -range measured by COMPASS, a more complex functional form is required, i.e. either a sum of two Gaussian functions or the so-called Tsallis function. All fits reproduce the data well and their inverse slopes agree well with one another already when using only statistical uncertainties in the fits.

Acknowledgements
We gratefully acknowledge the support of the CERN management and staff and the skill and effort of the technicians of our collaborating institutes. This work was made possible by the financial support of our funding agencies.