{'fraction_non_alphanumeric': 0.0410762331838565, 'fraction_numerical': 0.017219730941704037, 'mean_word_length': 4.294532627865961, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Investigations into the flavor dependence of partonic transverse momentum

Recent experimental data on semi-inclusive deep-inelastic scattering from the HERMES collaboration allow us to discuss for the first time the flavor dependence of unpolarized transverse-momentum dependent distribution and fragmentation functions. We find convincing indications that favored fragmentation functions into pions have smaller average transverse momentum than unfavored functions and fragmentation functions into kaons. We find weaker indications of flavor dependence in the distribution functions.


I. INTRODUCTION
Transverse-momentum-dependent (TMD) parton distribution functions (PDFs) and fragmentation functions (FFs) give a multi-dimensional description of partonic structure in momentum space. They are functions of the longitudinal and transverse momentum of partons, with respect to the reference hadron momentum. As such, they offer richer information compared to standard collinear PDFs and FFs, which depend only on the longitudinal momentum. In the last decade, TMD PDFs and FFs have gained increasing attention especially because of emerging data from experiments on semi-inclusive deep-inelastic scattering (DIS) (for reviews see, e.g., [1][2][3]).
In spite of this progress, we have still little knowledge about the most simple and most common of all TMD PDFs: the "unpolarized" distribution, f a 1 (x, k 2 ⊥ ), i.e., the distribution of partons with flavor a summed over their polarization and averaged over the polarization of the parent hadron. The features of the corresponding collinear standard PDF f a 1 (x) strongly depend on the parton flavor a (see, e.g., Refs. [4][5][6][7][8][9]). It comes natural, therefore, to question whether or not partons of different flavors have different transverse-momentum distributions. Several model calculations predict different transverse-momentum behaviors for different quarks [10][11][12][13][14][15][16], although others do not [17][18][19]. Indications of flavor dependence in TMD PDFs come also from pioneering studies in lattice QCD [20]. Therefore, we believe there are compelling motivations to study the flavor dependence of TMD PDFs.
The measurements recently published by the Hermes collaboration [21] are ideal to address this issue, since they refer to semi-inclusive DIS off different targets (protons and deuterons), with different final-state hadrons (chargeseparated pions and kaons), and with multidimensional binning. This is a landmark achievement in the knowledge of the internal structure of hadrons. Earlier data already gave some indications, but were limited in the variety of targets, or final-state hadrons, or multidimensional coverage (see, e.g., [22][23][24][25][26]).
The Compass collaboration has recently released similar data [27]. The amount of statistics is in this case impressive and the kinematic coverage is in general wider than at Hermes. However, at the moment these data are available only for deuteron targets and for unidentified final charged hadrons. Therefore, we decided not to use these data, although they will certainly play an essential role in the near future.
Dealing with semi-inclusive DIS, we need to consider also fragmentation functions and their transverse-momentum dependence. Also in this case, it is possible that different quark flavors fragment into different hadrons with characteristic transverse-momentum distributions [15,28]. This is another fundamental question that has never been addressed at the phenomenological level.
Since our work represents one of the first explorations on this topic, we adopt here a simplified framework, essentially based on a parton-model picture. We perform a leading-order analysis and neglect any modification that can be induced by QCD evolution, both in the collinear PDFs and FFs as well as in the TMD ones. This approximation is justified by the limited range in Q 2 of data: no difficulty arises in describing them with this simplified framework.
All our assumptions, the notation, and the general theoretical framework are briefly outlined in Sec. II. In Sec. III, we describe our fitting procedure. In Sec. IV, we present our final results, and in Sec. V we draw some conclusions and outlooks.

II. THEORETICAL FRAMEWORK
In one-particle semi-inclusive DIS, a lepton with momentum l scatters to a final state with momentum l off a hadron target N with mass M and momentum P , producing (at least) one hadron h in the final state with mass M h and momentum P h : (l) + N (P ) → (l ) + h(P h ) + X . (1) The space-like momentum transfer is q = l − l , with Q 2 = −q 2 . We introduce the usual invariants The available data refer to hadron multiplicities in semi-inclusive DIS, namely to the differential number of hadrons produced per corresponding inclusive DIS event. In terms of cross sections, we define the multiplicities as where dσ h N is the differential cross section for the semi-inclusive DIS process and dσ DIS is the corresponding inclusive one, and where P hT is the component of P h transverse to q. In the single-photon-exchange approximation, the multiplicities can be written as ratios of structure functions (see [29] for details): where ε = 1 − y − 1 4 γ 2 y 2 1 − y + 1 2 y 2 + 1 4 γ 2 y 2 .
We recall that the notation F XY,Z indicates the response of the hadron target with polarization Y to a lepton beam with polarization X and for the virtual photon exchanged in the polarization state Z. Therefore, the numerator of Eq. (4) involves semi-inclusive DIS processes with only unpolarized beam and target. We remark that the above expressions assume a complete integration over the azimuthal angle of the detected hadron. Acceptance effects may modify these formulae, due to the presence of azimuthal modulations in the cross section, though for the data used here such effects were included in the systematic uncertainties. We consider the limits M 2 /Q 2 1 and P 2 hT /Q 2 1. Within them, the longitudinal structure function F U U,L in the numerator of Eq. (4) can be neglected [30]. In the denominator, the standard inclusive longitudinal structure function F L is non negligible and contains contributions of order α S . However, in our analysis we assume a parton-model picture and we neglect such contributions; hence, consistently we neglect the contribution of F L in the denominator of Eq. (4). It may also be noted that in the transverse-momentum analysis of the data, F L induces a change in normalization that depends on x, but is independent of z and P 2 hT , the kinematic variables most relevant in the fitting procedure. Hence, we do not expect large effects on the resulting parameters.
To express the structure functions in terms of TMD PDFs and FFs, we rely on the factorized formula for semiinclusive DIS at low transverse momenta [31][32][33][34][35][36][37][38][39]: Here, H U U,T is the hard scattering part; f a 1 (x, k 2 ⊥ ; µ 2 ) is the TMD PDF for an unpolarized parton of flavor a in an unpolarized proton, carrying longitudinal momentum fraction x and transverse momentum k ⊥ at the factorization scale µ 2 , which in the following we choose to be equal to Q 2 . D a h 1 (z, P 2 ⊥ ; µ 2 ) is the TMD FF for an unpolarized Diagram describing the relevant momenta involved in a semi-inclusive DIS event: a virtual photon (defining the reference axis) strikes a parton inside a proton. The parton has a transverse momentum k ⊥ (not measured). The struck parton fragments into a hadron, which acquires a further transverse momentum P ⊥ (not measured). The total measured transverse-momentum of the final hadron is P hT . When Q 2 is very large, the longitudinal components are all much larger than the transverse components. In this regime, P hT ≈ zk ⊥ + P ⊥ (see also Ref. [42]).
parton of flavor a fragmenting into an unpolarized hadron h carrying longitudinal momentum fraction z and transverse momentum P ⊥ ; the term Y U U,T is introduced to ensure a matching to the perturbative calculations at high transverse momentum. The expression for F U U,T is known up to at least O(α 2 S ), including the resummation of at least nextto-next-to-leading logarithms of the type log (P 2 hT /Q 2 ). However, we are going to use here only the lowest-order expression, which should still provide a good description at low P 2 hT and in a limited range of Q 2 . Eventually, Eq. (6) simplifies to (see, e.g., Refs. [29,40,41]) where the convolution upon transverse momenta is defined as In Fig. 1, we describe our notation for the transverse momenta (in agreement with the notation suggested by the white paper in Ref. [2]), which is also reproduced below for convenience: The Gaussian ansatz consists in assuming the following functional form for the transverse-momentum dependence of both the TMD PDF f a 1 and the TMD FF D a h 1 in Eq. (7): Due to its simplicity, this ansatz has been widely used in phenomenological studies but with constant widths k 2 ⊥ and P 2 ⊥ . Here, for the first time we introduce an explicit dependence on flavor a for both average transverse momenta k 2 ⊥,a and P 2 ⊥,a h . In principle, there are no reasons to prefer the Gaussian ansatz over other functional forms, and indeed more flexible forms should be investigated in the future. Model calculations typically lead to a non-Gaussian behavior [10,12,[17][18][19]43]. The ansatz is also not compatible with the proper QCD evolution of TMD PDFs: it could be at most applicable at one specific starting scale, but would soon be spoiled by QCD corrections. In our analysis, we completely neglect Q 2 evolution, even in the collinear part of the functions, which we evaluate at Q 2 = 2.4 GeV 2 . We can do this only because the range in Q 2 spanned by the Hermes measurements is not large. From now on, we drop the Q 2 dependence of the involved functions.
The convolution on transverse momenta in Eq. (8) can be solved analytically: where the relation between the three variances is In this way, for each involved flavor a the average square value of the transverse momentum P hT of the detected hadron h can be related to the average square values of the intrinsic transverse momenta k ⊥ and P ⊥ , not directly accessible by experiments. Inserting Eq. (10) in Eq. (7), we simplify the multiplicities as If the distribution functions describe a parton a in a proton target, obviously the above expression is valid for N = p, i.e., for a proton target. We can deduce the corresponding result for a neutron target by assuming isospin symmetry. For a deuteron target, we can assume an incoherent sum of proton and neutron contributions. Under these assumptions the necessary label for the parent hadron on PDFs is omitted and PDFs refer to the ones of the proton. We remark also that each quark flavor is described by a single Gaussian with a specific width. The multiplicity is then a sum of Gaussians and thus no longer a simple Gaussian. The above expression can be used with minor modifications also if we assume that the distribution and fragmentation functions for some flavor are themselves sums of Gaussians. We will in fact adopt such an assumption for the up and down quarks, where we distinguish a valence and a sea contribution, each one having a different Gaussian width. For example, the up contribution to the multiplicities is where f uv 1 = f u 1 − fū 1 , and similarly for the down quark. Previous data obtained in unpolarized Drell-Yan and semi-inclusive DIS processes were compatible with calculations based on a Gaussian ansatz for unpolarized TMD PDFs and TMD FFs with flavor-independent constant widths. In this case, Eq. (12) would display a simple Gaussian behavior in P hT with the same width in every target-hadron combination. However, the Hermes multiplicities display significant differences between proton and deuteron targets, and between pion and kaon final-state hadrons. Hence, they strongly motivate our choice in Eq. (9) for a flavordependent Gaussian ansatz.

B. Assumptions concerning average transverse momenta
As mentioned in the previous section, we introduce different widths for the Gaussian forms of the valence and sea components of up and down TMD PDFs. However, we assume that the Gaussian widths of all sea quarks (ū,d, s ands) are the same (i.e., they have the same average square transverse momenta). State-of-the-art parametrizations of collinear PDFs have a more complex structure and introduce differences between sea quarks of different flavors; we leave this flexibility to future studies.
We choose the following functional form for the average square transverse momentum of flavor a: k2 ⊥,a , α, σ, are free parameters. For sake of simplicity, we keep the same exponents α and σ for all flavors. According to the above assumptions, we have three more parameters: k2 ⊥,a for a = u v , d v , sea. In total, we use five different parameters to describe all TMD PDFs. Since the present data have a limited coverage in x, we found no need of more sophisticated choices.
As for TMD FFs, fragmentation processes in which the fragmenting parton is in the valence content of the detected hadron are usually defined favored. Otherwise the process is classified as unfavored. The biggest difference between the two classes is the number of qq pairs excited from the vacuum in order to produce the detected hadron: favored processes involve the creation of at most one qq pair. If the final hadron is a kaon, we further distinguish a favored process initiated by a strange quark/antiquark from a favored process initiated by an up quark/antiquark.
For simplicity, we assume charge conjugation and isospin symmetries. The latter is often imposed also in the parametrization of collinear FFs [47], but not always [48]. In practice, we consider four different Gaussian shapes: The last assumption is made mainly to keep the number of parameters under control, though it could be argued that unfavored fragmentation into kaons is different from unfavored fragmentation into pions. As for TMD PDFs, also for TMD FFs we introduce a dependence of the average square transverse momentum on the longitudinal momentum fraction z, as done in several models or phenomenological extractions (see, e.g., Refs. [15,28,41,[49][50][51]). We choose the functional form The free parameters β, γ, and δ are equal for all kinds of fragmentation functions. In conclusion, we use seven different parameters to describe all the TMD FFs.

A. Selection of data
The Hermes collaboration collected a total of 2688 data points (336 points for each of the 8 combination of target and final-state hadrons), with the average values of (x, Q 2 ) ranging from about (0.04, 1.25 GeV 2 ) to about (0.4, 9.2 GeV 2 ), 0.1 ≤ z ≤ 0.9, and 0.1 GeV ≤ |P hT | ≤ 1 GeV. The collaboration presented two distinct data sets, including or neglecting vector meson contributions. Here, we use the data set where the vector meson contributions have been subtracted. In all cases, we sum in quadrature statistical and systematic errors and we ignore correlations. We always use the average values of the kinematic variables in each bin.
Our analysis relies on the assumption that the transverse-momentum-integrated multiplicities, m h N (x, z, Q 2 ), are well described by currently available parametrizations of collinear PDFs and FFs. However, this is not always true. In order to identify the range of applicability of the collinear results, we compared the multiplicities as functions of x and z with the leading-order (LO) theoretical predictions obtained using the MSTW08LO PDF set [8] and the DSS LO FF set [48]. In the comparison, we neglected the uncertainties affecting the PDFs but we included the ones  [48]. In all cases, the range 0.1 ≤ z ≤ 0.8 was included.
affecting the FFs, obtaining the latter from the plots in Ref. [52]. They are of the order of 5-10% for light quarks fragmenting into pions, of 10-15% for favored kaon FFs, of 50% for all the other cases, and they are larger at higher z.
In Tab. I, we quote the χ 2 per degree of freedom (χ 2 /d.o.f.) obtained in our comparison. Our results are different from the ones quoted in Tabs. IV and VIII of [48] for a few reasons: i) we used the final Hermes data, in particular the set with x and z binning; ii) we included also the lowest z bin (z < 0.2); iii) we did not include any overall normalization constant; iv) we included the theoretical errors on the extracted fragmentation functions. The comparison shows that in general the theoretical predictions do not describe the Hermes data very well. The agreement is particularly bad for π − and K − . However, this is not surprising because: i) the MSTW set of PDFs does not take into account semi-inclusive DIS data, ii) as mentioned above, the DSS set of FFs [48] was deduced using only a preliminary version of the Hermes multiplicities, iii) the Hermes data give very large contributions to the χ 2 of the global DSS analysis. Nevertheless, in our analysis we decided to restrict the ranges to Q 2 > 1.4 GeV 2 and 0.1 < z < 0.8, i.e., excluding the first bin in Q 2 (equivalent also to the lowest x) and the last bin in z. Inclusion of decays from exclusive vector-mesons markedly degrades the χ 2 of the pion channels and increases the global χ 2 (cf. the first and second line of Tab. I). Hence we will present results for only the fits to vector-meson subtracted multiplicities. We checked that our basic conclusions do not change when using data without vector-meson subtraction.
We also noted that a description of data of comparable quality could be achieved by turning off the Q 2 dependence of both collinear PDFs and FFs, and by computing them at the fixed value of Q 2 = 2.4 GeV 2 (cf. the first and third line of Tab. I). Therefore, we decided to systematically neglect any contribution of QCD evolution and to compute all theoretical quantities at the average value of Q 2 = 2.4 GeV 2 .
When considering also the transverse-momentum dependence, the TMD formalism is valid only when P 2 hT Q 2 . In order not to exclude too many data points, we apply the loose requirement P 2 hT < Q 2 /3. This leads to the exclusion of at most two bins at high P 2 hT and low Q 2 . Finally, we exclude also the data points with the lowest |P hT | (|P hT | < 0.15 GeV). A priori, there is no reason to exclude them, but in our attempts we found them particularly difficult to fit, mainly because they often do not follow the trend of the other data points and at the same time they have small errors. In order to be able to fit them, we need to increase the flexibility of our functional forms. We leave this task to future studies.
In summary, we use a total of 1538 data points, approximately 190 for each of the 8 combinations of target and final-state hadrons, which correspond to about 60% of the total 2688 points measured by the Hermes collaboration.

B. Fitting procedure and uncertainties
The fit and the error analysis were carried out using a similar Monte Carlo approach as in Ref. [53], and taking inspiration from the work of the NNPDF collaboration (see, e.g., [54][55][56]). The approach consists in creating M replicas of the data points. In each replica (denoted by the index r), each data point i is shifted by a Gaussian noise with the same variance as the measurement. Each replica, therefore, represents a possible outcome of an independent experimental measurement, which we denote by m h N,r (x, z, P 2 hT , Q 2 ). The number of replicas is chosen so that the mean and standard deviation of the set of replicas accurately reproduces the original data points. In our case, we have found that 200 replicas are more than sufficient.
The standard minimization procedure is applied to each replica separately, by minimizing the following error The sum runs over the i experimental points, including all species of targets N and final-state hadrons h. The theoretical multiplicities m h N,theo and their error ∆m h N,theo do not depend on Q 2 , as explained in the previous section. They are computed at the fixed value Q 2 = 2.4 GeV 2 using the formula in Eq. (12). However, in each z bin for each replica the value of D a h 1 is independently modified with a Gaussian noise with standard deviation equal to the theoretical error ∆D a h 1 . The latter is estimated from the plots in Ref. [52] and it represents the main source of uncertainty in ∆m h N,theo . Finally, the symbol {p} denotes the vector of fitting parameters. The minimization was carried out using the Minuit code. The final outcome is a set of M different vectors of bestfit parameters, {p 0r }, r = 1, . . . M, with which we can calculate any observable, its mean, and its standard deviation. The distribution of these values needs not to be necessarily Gaussian. In this case, the 1σ confidence interval is different from the 68% interval. The 68% confidence interval can simply be computed for each experimental point by rejecting the largest and the lowest 16% of the M values.
Although the minimization is performed on the function defined in Eq. (20), the agreement of the M replicas with the original data is better expressed in terms of a χ 2 function defined as in Eq. (20) but with the replacement m h N,r → m h N , i.e., with respect to the original data set. If the model is able to give a good description of the data, the distribution of the M values of χ 2 /d.o.f. should be peaked around one. In practice, the rigidity of our functional form leads to higher χ 2 values.

IV. RESULTS
In this section, we describe the results obtained by fitting the Hermes multiplicities with the theoretical formula of Eq. (12) and using the Monte Carlo method outlined in the previous section. We performed different kinds of fits with different assumptions. The first one, conventionally named "default fit," includes all the 1538 data points selected according to the criteria described in Sec. III A. In the second one, we excluded data also for the second lowest Q 2 bin, i.e., by selecting Q 2 > 1.6 GeV 2 . This cut reduces the number of data points to 1274. The third scenario corresponds to neglecting kaons in the final state and taking only multiplicities for pions. The last scenario is a fit of the default selection using a flavor-independent Gaussian ansatz. Before discussing each different scenario, here below we list their common features.
As repeatedly mentioned above, in our analysis we neglected completely the effect of Q 2 evolution, even in the collinear PDFs and FFs, and we evaluated all observables at the fixed value Q 2 = 2.4 GeV 2 .
As for the dependence of the TMD average transverse momentum on x, we noticed that the fit is weakly sensitive to the exponents in Eq. (14). We tried fits with α = σ = 0 and obtained good results. However, in order to stress the fact that present data do not constrain these parameters very well, we decided to assign random values extracted from uniform distributions to both the exponents: we consider α as a random number between 0 and 2 and σ as a random number between −0.3 and 0.1 . Better determinations of these parameters require an extended range in x, together with uncorrelated x and Q 2 binnings. The dependence of the TMD FF average transverse momentum on z is governed by Eq. (19); in this case, we decided to keep all three parameters free.
For each scenario, we performed 200 replicas of the fit. In this section, we present the 68% confidence intervals of the parameters over the replicas, computed by rejecting the largest and the lowest 16% of the replicated parameter values. We quote the values as A ± B, where A is the average of the upper and lower limits of the 68% confidence interval and B is their semi-difference. It is understood that much more information is available by scrutinizing the full set of 200 values for each of them. 2 In Tab. II, we list the 68% confidence intervals of the χ 2 /d.o.f. for the different scenarios, including the global result and the outcome for each target-hadron combination, separately. In Tabs. III and IV, we list the 68% confidence intervals for the five fitting parameters for TMD PDFs and for the seven fitting parameters for TMD FFs, respectively.
In all fits, we observe a strong anticorrelation between the distribution and fragmentation transverse momenta. This is not surprising, since the width of the observed P hT distribution is given by Eq. (11). To better pin down the values of k 2 ⊥,a and P 2 ⊥,a h separately for the various flavors a, it will be essential to include also data from electron-positron annihilations and Drell-Yan processes. In any case, a common feature of all scenarios is that the k2 ⊥,a (namely, the average squared transverse momenta of TMD PDFs at x = 0.1) have average values around 0.3 GeV 2 , while the P 2 ⊥,a h (namely the average square transverse momenta of TMD FFs at z = 0.5) have average values around 0.18 GeV 2 . Moreover, the fits prefer large values of the exponents β and δ for TMD FFs, but with large uncertainties; the parameter γ is usually small.
Here below, we discuss in detail the results for the four different scenarios.

A. Default fit
In this scenario, we consider all 1538 data points selected according to the criteria explained in Sec. III A. The quality of the fit is fairly good. The global χ 2 /d.o.f. is 1.63 ± 0.13. In Fig. 2, the distribution of the χ 2 /d.o.f. over the 200 replicas is shown. Many replicas have χ 2 /d.o.f. > 1.5. This indicates some difficulty to reproduce the data correctly. It is not surprising if we take into account that the description of the collinear multiplicities was already difficult (see Tab. I). It may actually seem contradicting that our fit is able to describe the transverse-momentumdependent multiplicities relatively well. This is probably simply due to the fact that the multidimensional binning has many more data points but with much larger statistical errors.
In Tab. II, we list the 68% confidence intervals of the χ 2 /d.o.f. also for each target-hadron combination N → h, separately. The worst result is for D → K − . This may be a bit surprising, also because p → K − is described very well. However, this may be due to the fact that the collinear description of this channel is poor (see Tab. I). We point out also that the systematic errors in D → K − are significantly smaller than p → K − [21]. The second worst agreement is for p → π + , which is not unexpected since statistical errors are smallest in this channel. The π − channels are described decently, which is at odds with the poor description of their collinear multiplicities (see Tab. I). We do not have a reasonable explanation for this feature yet. Maybe, it could be ascribed to the cuts in P hT that we implemented in our fit.   for one selected bin x ∼ 0.15 and Q 2 ∼ 2.9 GeV 2 (out of the total five x bins we used), and for four different z bins (out of the total seven z bins we used). The lowest P 2 hT bin was excluded from the fit, as explained in Sec.III A. The theoretical band is obtained by rejecting the largest and lowest 16% of the replicas for each P 2 hT bin. The theoretical uncertainty is dominated by the error on the collinear fragmentation functions D 1 (z), which induces an overall normalization uncertainty in each z bin. The different values of the fit parameters in each replica are responsible for the slight differences in the slopes of the upper and lower borders of the bands.
In Tab. III, the values of the average square transverse momenta for TMD PDFs are listed. We note that they can range between 0.13 and 0.57 GeV 2 within the 68% confidence interval.
In the left panel of Fig. 5, we compare the ratio k 2 ⊥,dv / k 2 ⊥,uv vs. k 2 ⊥,sea / k 2 ⊥,uv for 200 replicas. The white box represents the point at the center of each one-dimensional 68% confidence interval of the two ratios. The shaded area represents the two-dimensional 68% confidence region, it contains 68% of the points with the shortest distance from the white box. Since for each flavor the x dependence of the average square transverse momenta is the same (see Eq. (14)), these ratios are x-independent. The dashed lines correspond to the ratios being unity and divide the plane into four quadrants. Most of the replicas are in the upper left quadrant, i.e., we have k 2 ⊥,dv < k 2 ⊥,uv < k 2 ⊥,sea . The white box shows that d v is on average about 20% narrower than u v , which is in turn about 10% narrower than the sea. The crossing of the dashed lines corresponds to a flavor-independent distribution of transverse momenta. This crossing point lies at the limit of the 68% confidence region. In a relevant number of replicas d v can be more than 40% narrower than the u v , and the sea can be more than 30% wider than u v . From this fit, it seems possible that the sea is narrower than u v , but unlikely that d v is wider than u v .
In the right panel of Fig. 5, we compare the ratio P 2 ⊥,unf / P 2 ⊥,fav vs. P 2 ⊥,uK / P 2 ⊥,fav in the same conditions as before. All points are clustered in the upper right quadrant and close to its bisectrix, i.e., we have the stable outcome that P 2 ⊥,fav < P 2 ⊥,unf ∼ P 2 ⊥,uK . The width of unfavored and u → K + fragmentations are about 20% larger than the widht of favored ones. In this scenario, we restrict the Q 2 range compared to the default fit by imposing the cut Q 2 > 1.6 GeV 2 . The set of data is reduced to 1274 points. The mean value of the χ 2 /d.o.f is smaller, since we are fitting less data. Moreover, the disregarded Q 2 bin contains high statistics. As for the default fit, the behavior of transverse momenta over the 200 replicas is summarized in Fig. 6. The exclusion of low-Q 2 data leads to partial differences in the features of the extracted TMD PDFs: the average width of valence quarks slightly increases, while the distribution for sea quarks becomes narrower. In the left panel, most of the replicas are in the lower left quadrant, i.e., we have On average, d v quarks are 15% narrower than than u v quarks, which are in turn more than 20% wider than sea quarks. In a relevant number of replicas d v can be more than 40% narrower than the u v , and the sea can be even 50% narrower than u v . In this scenario, it is unlikely that the sea is wider than u v , but it is possible that d v is wider than u v .
In the right panel, the behavior of transverse momenta in fragmentation processes is qualitatively unchanged with respect to the default fit, apart from the fact that the unfavored Gaussian function becomes now more than 25% larger than the favored one.
The crossing point again indicates no flavor dependence and lies just outside the 68% confidence region for TMD PDFs and completely outside the same region for TMD FFs.
We conclude that the low-Q 2 data, being also characterized by low x, can have a significant impact on the analysis of TMD PDFs, in particular the sea components. More data at low x (but possibly at high Q 2 ) are necessary to better constrain the sea quarks TMD PDFs [2,27,57].

C. Fit with pions only
We also choose to fit data related only to pions in the final state, in order to explore the importance of the kaons data set. In this framework, we are left with two independent fragmentation processes: favored and unfavored ones. Accordingly, the number of fit parameters for TMD FFs reduces from 7 to 5 ( P 2 ⊥,fav , P 2 ⊥,unf , β, δ, γ; see Eqs. (15)- (19)).
The agreement between data and the model is the worst (see Tab. (I)). This is due to at least two reasons. First of all, the fit of collinear multiplicities was poor in all the target-hadron combinations involving pions in the final state. Moreover, the high statistics collected for pions (mostly in the low-Q 2 region) leads to higher values of χ 2 . Fig. 7 shows the behavior of transverse momenta over the 200 replicas. As for TMD PDFs, in the left panel most of the replicas are in the lower part, i.e., we have k 2 ⊥,sea < k 2 ⊥,dv < ∼ k 2 ⊥,uv . On average, d v quarks are equally distributed as u v quarks, which are in turn more than 20% wider than sea quarks. In the default fit sea quarks were wider than valence ones and there was a remarkable difference between u v and d v , not evident in this scenario. In any case, in a relevant number of replicas d v can be more than 15% narrower than the u v , but also more than 10% wider than u v . The sea can be even 50% narrower than u v , but it is also not unlikely that the sea is wider than u v Once again, the crossing point for flavor independence lies at the boundary of the 68% confidence region, due to the difference between the distributions of sea quarks and valence quarks. As in the other scenarios, TMD FFs for unfavored processes are wider than favored ones. The difference is comparable to the default fit, with unfavored functions about 20% larger than favored.
Similar fits have been performed in Refs. [24,26], but using data averaged over z, which renders it particularly difficult to disentangle the distribution and fragmentation contributions. To overcome this problem, both fits included also indirect information from the azimuthal cos φ h dependence. The fit of Ref. [26] obtained a small value for the distribution mean square transverse momenta of up quarks, k 2 ⊥,u = 0.07 ± 0.03 GeV 2 , while the down quark mean transverse momentum was compatible with zero, k 2 ⊥,u = −0.01 ± 0.05 GeV 2 (sea quarks were neglected). The previous fit [24], obtained a somewhat different behavior, with a mean transverse momentum of the up quark compatible with zero and k 2 ⊥,d = 0.11 ± 0.13 GeV 2 . In both fits, the average values of the width of the TMD FFs are compatible with our results, but, contrary to our findings, a slight tendency for the favored FF to be larger than unfavored was found. In any case, we remark that the average kinematics of the experiment taken into consideration in Refs. [24,26] are different from Hermes (see also the discussion in Ref. [58]).
Overall, we conclude that kaon data have an important impact in a flavor-dependent analysis, due to the large role played by strange quarks and antiquarks in kaon multiplicities.

D. Flavor-independent fit
In this scenario, we assume a Gaussian ansatz for unpolarized TMD PDFs and TMD FFs with flavor-independent widths, i.e., we neglect any flavor dependence in Eq. (9): Accordingly, the number of fit parameters reduces to 3 for TMD PDFs ( k2 ⊥ , α, σ) and 4 for TMD FFs ( P 2 ⊥ , β, δ, γ). Their values are summarized in Tab. III and IV. The expression (12) for the multiplicities considerably simplifies and the P hT width is the same for every target-hadron combination: The agreement between data and the flavor-independent model is poorer than in the (flavor-dependent) default fit: the central value of the χ 2 /d.o.f. is 1.73 (see Tab. II). This is not surprising, since we are fitting with the same function data for all the available target-hadron combinations, which display sensibly different behaviors. However, these results do not rule out the flavor-independent ansatz. Fig. 8 clearly shows the anti-correlation between k2 ⊥ and P 2 ⊥ induced by Eq. (23). Similar fits have been performed in the past for semi-inclusive DIS and Drell-Yan processes [41,58], also including the effect of gluon resummation [59][60][61] within the so-called Collins-Soper-Sterman formalism [32,62], which is equivalent to taking into account TMD evolution [35,36,38].
The latter functional form is not much different from the one we obtained. The value of the distribution transverse momentum is slightly smaller than our average value, which is compensated by the fact that the fragmentation transverse momentum is slightly higher. Other fits that explored the z dependence in the Gaussian width of TMD FFs can be found in Refs. [41,49]. Comparison with extractions from Drell-Yan experiments (see, e.g., Refs. [41,[58][59][60][61]) is not straightforward, due to the different kinematic conditions and the difficulty to extrapolate the results obtained in the CSS formalism (see also the discussion in Ref. [65]). The mean square transverse momentum obtained from Gaussian fits without TMD evolution [41,58] is larger than in our case, k 2 ⊥ > ∼ 0.7 GeV 2 .

V. CONCLUSIONS
Using the recently published Hermes data on semi-inclusive DIS multiplicities [21], we explored for the first time the flavor dependence of the transverse momenta of both the unpolarized parton distributions (TMD PDFs) and fragmentation functions (TMD FFs). We adopted a simplified framework based on the parton model and neglecting the effects of QCD evolution. Using a flavor-dependent Gaussian ansatz, we obtained different results for multiplicities in eight different target-hadron combinations. We performed several fits of the data in different scenarios: including all bins as described in Sec. III A (the "default fit"), excluding data with Q 2 ≤ 1.6 GeV 2 (equivalent to excluding partons at low x), selecting only pions in the final state, or neglecting any flavor dependence.
Comparing the default fit and the flavor-independent one, we conclude that the flavor-dependent Gaussian ansatz performs better. The difference between the average χ 2 /d.o.f. in the two cases is not striking but, nonetheless, appreciable. We find convincing indications that the unfavored fragmentation functions have larger average transverse momenta with respect to pion favored fragmentation functions. We get weaker indications of flavor dependence for the TMD PDFs. It is very likely to find fits of the available data with differences of the order of 20% in the the mean square transverse momenta of different flavors. In particular, our default fit shows a tendency for valence down quarks to have a narrower distribution than the one of valence up quarks, which in turn is narrower than the one for sea quarks. These features have a potentially large impact on the polarization-dependent TMD extractions [66][67][68][69][70][71][72], where usual flavor-independent transverse momentum parametrizations are assumed in the fragmentation, as even the normalizations extracted for those TMDs depend, directly or indirectly, on the widths of the polarization-averaged TMD FFs.
Apart from the ratios among different flavors, the absolute values of the mean square transverse momenta are compatible with results quoted in the literature. However, it should be kept in mind that there exist strong anticorrelations between mean squared transverse momenta of distribution and fragmentation functions.
This work is a first step in the exploration of the transverse momentum dependence of partons inside hadrons. First of all, it needs to be updated by implementing evolution equations in the TMD framework [36][37][38][39]65]. Secondly, the data set needs to be enlarged to include the recently released Compass data in a wider kinematical domain [27], and, in the following step, to include also data from e + e − annihilations and Drell-Yan processes. Finally, other functional forms different from the Gaussian ansatz should be explored.