PDF bias and flavor dependence in TMD distributions

Transverse momentum dependent (TMD) distributions match collinear parton density functions (PDF) in the limit of small transverse distances, which is accounted for by global extractions of TMD distributions. We study the influence of the collinear PDF value and uncertainties on the determination of unpolarized TMD distributions and the description of Drell-Yan (DY) and Z-boson production measurements at low transverse momenta. We take into account, for the first time, flavor-dependent non-perturbative TMD profiles. We carry out a Bayesian analysis to incorporate the propagation of PDF uncertainties into TMD extractions. We find that collinear PDF uncertainties and non-perturbative TMD flavor dependence are both essential to obtain reliable TMD determinations, and should be included in future global analyses.


Introduction
While determinations of the proton's collinear parton distribution functions (PDFs) [1] from fits to experimental data have been carried out for four decades, and PDFs are nowadays known with high precision, the systematic investigation of the transverse momentum dependent (TMD) [2] distribution functions has started much more recently. In the last few years determinations of TMD distributions [3] from fits to experimental data have been performed in the context of low transverse momentum Drell-Yan (DY) production and semi-inclusive deep inelastic scattering (SIDIS) [4][5][6][7][8], small-x deep inelastic scattering [9], nonlinear TMD evolution [10], and parton branching [11]. The results are collected in the public library TMDlib [12,13], similarly to LHAPDF for the case of PDF.
A physical cross section σ with hard scale Q and measured transverse momentum q T Q is described in terms of TMD distributions by a factorization formula with the schematic form (up to power corrections in q T /Q and Λ QCD /Q) where b is the transverse distance Fourier conjugate to transverse momentum, σ ij are perturbatively calculable partonic cross sections, and f i and f j are TMD parton distributions. These depend on the mass and rapidity scales µ and ζ through appropriate evolution equations, involving both perturbative and non-perturbative (NP) components, and are to be determined from fits to experiment.
Current fits of TMD distributions from DY and SIDIS measurements are performed by using not only the TMD factorization and evolution framework but also additional inputs which exploit the relationship between TMD distributions and the collinear PDFs. These relations follow from the operator product expansion (OPE) applied to the TMD operator. The OPE expresses the TMD distributions at small distances b in terms of PDFs via perturbatively calculable coefficients, with power corrections in b, as follows where f i (x, b; µ, ζ) is the TMD distribution, q j (x, µ) is the PDF, C ij are the Wilson coefficient functions, µ and ζ are factorization scales, the subscripts f and f indicate parton flavors, and the last term on the right hand side is the power suppressed correction at small b. An expansion such as eq. (1.2) is used, for instance, in the context of precision studies of the DY transverse momentum q T spectrum at the Large Hadron Collider (LHC) [14] to carry out comparisons of small-q T logarithmic resummations between computer codes based on TMD distributions (e.g., artemide [6], NangaParbat [7]) and computer codes which perform perturbative small-q T resummation without systematically introducing TMD distributions (e.g., DYturbo [15,16], reSolve [17,18], Radish [19,20], SCETlib [21,22], CuTe [23]). To carry out TMD fits an ansatz is made for the NP TMD distribution at large b, where the power corrections to OPE become sizeable. This is of the form [4][5][6][7][8][24][25][26] 3) where f f NP is a function which contains power corrections to eq. (1. 2), behaves as f f NP (x, b) ∼ 1 + O(b 2 ) for b → 0, and is to be fitted to experimental data. This ansatz reproduces eq. (1.2) for b → 0 but allows for modifications at finite b. Current fits determine the NP TMD distributions, once a choice is made for the collinear PDFs in eq. (1.3).
The PDF choice in eqn. (1.3) determines the value of TMD distribution and influences f NP (we address this as PDF bias). Consequently, the extraction based solely on the central value of PDF are incomplete, and PDF uncertainties should be incorporated. The way the PDF error propagates has nontrivial consequences which have not been studied in the literature so far. The purpose of this work is to develop methods for addressing these issues, to critically assess the current status of phenomenological analyses based on TMD factorization and OPE expansion, and to perform the first systematic investigation of the role of PDF bias in TMD determinations.
To this end, we take PDF sets HERA20 [27], NNPDF3.1 [28], CT18 [29] and MSHT20 [30] as representatives of different methodological approaches at next-to-next-to-leading order (NNLO); we use the implementation of the TMD factorization formula as in refs. [4,31] to perform fits of experimental measurements for unpolarized DY and Z-boson production at small q T , from fixedtarget to LHC energies; we devise and implement a Bayesian procedure to propagate the PDF uncertainties, for each PDF set, to TMD extractions. This enables us to present results for the fitted NP TMD distributions which display, for the first time, both experimental and PDF uncertainties. We find that the role of the collinear PDFs is very significant in current TMD phenomenology. We discuss the issues associated with the combination of the different sources of TMD uncertainty, and propose an approach based on the bootstrapping method. An important aspect addressed by our work concerns the physical properties of the NP distributions f f NP in eq. (1.3). The fits based on TMD factorization and evolution which have so far been performed in the literature have assumed flavor-independent parameterizations for f f NP , with the flavor dependence solely in collinear PDFs, motivated by the belief that available data have little sensitivity to TMD flavor structure. In this work we go beyond this assumption and include flavor dependence in the NP TMD distributions. We find that the the inclusion of TMD flavor dependence improves the quality of the fits significantly, leading both to a better agreement of theory with data and to more consistent results among different collinear PDF sets. For each set, we examine the distribution of χ 2 values over PDF replicas. We demonstrate that the spread in the χ 2 distribution among replicas is much reduced with flavor-dependent TMD compared to the flavor-independent case, and that this conclusion applies consistently across collinear sets.
A complete treatment of the problems studied in this paper should eventually involve simultaneous fits of PDF and TMD distributions, requiring much more computational power than is available to us at present. The analysis presented in this paper underlines essential elements which should be included in future TMD studies.
The article is structured as follows. Sec. 2 provides basic theoretical inputs and notation, and presents the flavor dependent model used for f f NP . Sec. 3 describes the data used in the analysis. The statistical methods are reported in sec. 4. Sec. 5 presents the results of the fits and the discussion of their consequences. We give conclusions in sec. 6. We collect further details on the fits in the appendix.

Theory inputs
In this section we summarize the basic theoretical inputs of our analysis: i) TMD factorization formula for the DY cross section and evolution equations for the TMD distributions, and ii) models for the non-perturbative contributions to TMD factorization and evolution formulas.
The new element in this section is the flavor-dependent TMD profile in eq. (2.4): to our knowledge, flavor dependence of the TMD has not been considered before in the context of DY or SIDIS fits including TMD evolution.

DY factorization formula
The TMD factorized expression for the differential cross-section of the DY process h 1 + h 2 → γ * /Z(→ + − ) + X in the low transverse momentum region can be written as [32,33] dσ Here Q 2 is the invariant mass of the vector boson, q T is the transverse component of its momentum (relative to the scattering plane), and Λ QCD is the characteristic hadronic scale. The index f runs over all active quark flavors. The distribution f 1 is the unpolarized TMDPDF, depending on the lightcone momentum fraction x, transverse distance b, mass and rapidity scales µ and ζ; σ 0 is the leading-order DY cross section, H is the perturbative hard factor containing higher-order corrections, and J 0 is the Bessel function. The formula is valid up to power-suppressed corrections in q T /Q and Λ QCD /Q, indicated in the second line of eq. (2.1). The DY factorization formula (2.1) has been rederived in [34][35][36]. We use its implementation as given in refs. [4-6, 24, 25, 31, 37-41].
The evolution of TMD distributions in mass and rapidity is given by the pair of equations where Γ cusp is the anomalous dimension for the cusp of light-like Wilson lines, γ V is the anomalous dimension of the quark vector form factor and D is the Collins-Soper (CS) kernel [42,43]. We will use the solution of the evolution equations according to the ζ-prescription [4,6,24,37]. Using the expansion (1.3), this solution is expressed in terms of PDFs, f f NP and the matching coefficients [44][45][46][47].
The perturbative orders in the strong coupling α s which will be used in the calculations presented in the following sections are summarised in tab. 1 for each element of the cross section. The definition of D resum is given in the next subsection. The resulting cross section corresponds to the logarithmic order NNLL according to the terminology in [14]. Table 1: Summary of the perturbative orders used for each element of the cross section.

Non-perturbative models
The NP content of the DY differential cross section is encoded in two functions: the CS kernel D in eq. (2.2) and the TMD distribution f f NP in eq. (1.3). We write the CS kernel as where the first term D resum on the right hand side is the resummed expression for D at small b [48] with NNLO perturbative coefficients [49,50], while the second term is the non-perturbative model, given in terms of b * (b) = b/ 1 + (b/B NP ) 2 , with parameters c 0 and B NP to be determined from data. Motivated by the results of the fit [4], for simplicity in what follows we will take a fixed value B NP = 2 GeV −1 and leave only c 0 to be fitted to data. The model in eq. (2.3) is linear at b → ∞: this ansatz is partially supported by the studies [51,52] and lattice computations [53][54][55], and it has already been used in ref. [4]. In refs. [31,56] alternative models, with quadratic asymptotic behavior (as in the studies [57][58][59]) and with constant asymptotic behavior (in the spirit of the s-channel picture [60][61][62]), have also been analyzed. Corresponding investigations based on these different asymptotic behaviors are left to future work. The main new feature of the NP treatment in the present analysis concerns the TMD distribution f f NP . We include flavor dependence in the TMD profile, and take the parameterization with λ f 1,2 > 0 and λ 0 > 0. The model is characterized by an exponential asymptotic fall-off at b → ∞ and a Gaussian-like shape at intermediate b [6,58,[63][64][65][66]. While the parameter λ 0 is taken to be universal for all flavors, the parameters λ 1,2 are taken to be flavor dependent. We distinguish u, d,ū,d and sea cases, where sea is used for (s,s, c,c, b,b) flavors. In total we have 11 free parameters. Since the parameters λ 1,2 are almost uncorrelated across flavors, the presence of unnecessary fitting parameters (i.e. not well restricted in a particular setup) does not lead to an overfit.
The numerical implementation is made with artemide [6], which can be found in open-access repository [67]. The PDF values and their evolution are taken from LHAPDF [68]. Artemide is a FORTRAN code, with a Python interface. The evaluation of a single DY-cross-section point includes two convolutions of PDFs, Hankel type integrals and three integrations over the phase space. Additionally, for some measurement one needs to take into account fiducial cuts. The artemide code is optimized for such computations, and it uses various numerical tricks, such as precomputed grids, specially optimized integration algorithms and parallel evaluation. Even so the computation of one value of χ 2 for the data set takes from 30 seconds to a few minutes, depending on the PDF input, NP parameters, experimental fiducial cuts and hardware. Therefore, a single minimisation procedure (made with the help of the iMINUIT package [69]) requires typically a few dozen hours on an average computer.

Complete data set
The complete data set used in this study is given in tab. 2. We restrict the fit to data points in the low transverse momentum region by implementing the cut q T /Q < 0.25 1 . In addition we implement an extra cutting rule for the very precise data typically encountered at LHC. Given a data point p(1 ± σ), with p being the central value and σ its uncorrelated relative uncertainty, corresponding to some values of q T and Q, we include it in the fit only if In other words, if the (uncorrelated) experimental uncertainty of a given data point is smaller than the theoretical uncertainty associated to the expected size of the power corrections, we drop it from the fit. The resulting data set in tab. 2 contains 507 data points, and spans a wide range in mass, from Q = 4 GeV to Q = 150 GeV, and in x, from x ∼ 0.5 · 10 −4 to x ∼ 1. The data roughly split into low-energy and high-energy points. The low-energy subset contains data from fixed target experiments E288 [73], E605 [74] and E772 [75], and PHENIX [76]. The high-energy subset contains data from the neutral current DY (Z/γ-boson) measured at the Tevatron [77][78][79][80][81] and LHC [82][83][84][85][86][87][88][89]. The subsets have similar number of points but those in the high-energy subset are one order of magnitude more precise. The data that we are using are very similar to those included in previous fits [4,5,7] but we consider also the recent CMS Z-boson production data at 13 TeV 2 [86].
The treatment of particular aspects of the measurements such as bin-integration, nuclear modifications and normalization conditions are as in ref. [4] and we refer to it for extra details. Notice that we use the absolute values of the cross-section, whenever available, and perform the full binintegration to accurately incorporate the phase-space effects. The only modification of the data treatment in comparison to [4] is the splitting of the low-energy measurements (namely E288 at E beam = 300 and 400 GeV, E605 and E772) into two independent subsets below and above the Υ-resonance. This allows us to treat the normalization error independently for these energies. The reason for doing this is a possible inconsistency between these energy regions inside the measurements, which leads to tensions. Total 507 309 *Bins with 9 Q 11 are omitted due to the Υ resonance. Table 2: Summary table for the data included in the fit. For each data set we report: reference publication, centre-of-mass energy, coverage in Q and y or x F , possible cuts on the fiducial region, and number of data points that survive the cut of eq. (3.1).

Reduced data set
Given the data set defined above, we observe that the sensitivity to NP parameters of the TMD distributions varies strongly within the data set. For instance, data points with q T ∼ 15-25 GeV have little power in constraining NP parameters, irrespective of their precision, because they are deeply inside the resummation region, where the Hankel integral is dominated by the b 1GeV −1 contribution, and determined by PDF values only. The inclusion of these points in the fit is not harmful but rather time consuming. To speed up the fitting procedure, we have considered a reduced set of data points.
To identify the relevant points for the TMD extraction, we tested the sensitivity of the theory prediction to the variation of NP parameters, and included only the points that are sensitive to the variation. Given a NP parameter p, we compute the cross-section for several values of p distributed in the range p ± δp, and compute the sensitivity coefficient between p and the cross-section, by the formula where ∆σ is the uncorrelated experimental uncertainty of the point. The sensitivity coefficient indicates how strongly a prediction for a data point depends on the parameter p. E.g., s σ,p = 1 means that the variation of p by δp gives rise to a variation in the theory prediction equal to the experimental uncertainty.
We computed the sensitivity coefficient for each point in the data set and for each parameter of our NP ansatz. The values of δp are taken using the Hesse estimation for its uncertainty multiplied by a factor 5. Points with s > 0.4 for at least one parameter were included in the reduced set. This test was run for each of the PDF sets that we studied, and the union of data was taken. This provided us with a very conservative selection of impacting data.
Afterwards, we excluded the sets with too few points, such as LHCb (13 TeV) (1 surviving point) and CDF (run1) (5 points remaining out of 33). Additionally, we excluded the ATLAS (7 TeV) measurement since it has shown an anomalous behaviour and provides significant tension with other similar measurements 3 .
The reduced set contains 309 data points. Most excluded points are from the region characterized by high energy and high q T . The low-energy subset, in contrast, is very sensitive to the NP input and lost only 15 points due to their large error bands. We have checked, by evaluating several random cases, that the minimum of the χ 2 of the reduced set almost coincides (up to 2 digits) with the minimum of χ 2 of the complete set.

Statistical framework
The central point of the present study is the impact of the PDF uncertainty on the extraction of TMD distributions. Therefore, we distinguish two sources of uncertainties: i) the experimental uncertainties and ii) the uncertainty of the collinear PDFs. These have different nature, and could not be easily combined together. In fact, the best approach to account for the PDF uncertainty would be a simultaneous global fit of collinear PDF and TMDPDF, which is challenging. In the present work, we use the strategy of combining PDF and experimental uncertainty based on Bayesian statistics, described in this section.

Definition of the χ 2 -function
Our main tool is the χ 2 -test function, which estimates the goodness of a theory prediction. It is defined as where i and j enumerate points of the data set. For the ith point, t i is the corresponding theoretical prediction, m i is the experimental value, and all the information about the uncertainties for individual data points and the relation between them is given by the covariance matrix V ij . The definition of V ij is crucial for an adequate analysis of the data and interpretation of the results. Its construction distinguishes two types of uncertainties: uncorrelated and correlated. In general, the ith point is given as Here σ i,stat and σ i,unc are the uncorrelated statistical and systematic uncertainties, respectively, that estimate the degree of knowledge of the ith data point regardless of all other data in the set. The correlated uncertainties σ (k) i,corr (with k = 1, ..., n) estimate the relation between the statistical fluctuations of the ith point and all others in the set. The covariance matrix V ij can be written as [90,91] Using this expression in the computation of eq. (4.1) allows us to obtain a reliable estimation of data-to-theory agreement, while taking into account the nature of experimental uncertainties. The minimum of the χ 2 -test function determines the preferred values of the NP parameters λ, which in turn characterize the preferred shape of TMD distributions. The fit procedure consists in the minimization of χ 2 , which is performed with the iMinuit package [69].
In fitting a multivariate function one usually finds several sets of parameters giving numerous local minima of the function. Some of these give mathematically correct but physically unsound results, and therefore a strategy should be implemented to discard them. The MINUIT algorithm can handle the problem of choosing the correct minimum, unless the minima are significantly separated in the parameter space and/or the step (the amount by which the parameters move in one iteration) is too small. In our case, with the χ 2 -function depending on 11 parameters, we have faced the problem of the minimisation algorithm selecting a definitely unphysical minimum. Each occurrence is characterized by values of parameters that (almost) suppress a contribution of some flavor.
In order to disregard such contributions we introduced penalty terms into the minimisation function. The penalty function is defined as This function is non-zero only if the parameters of the u/d andū/d distributions differ from each other by an order of magnitude. So the minimisation was performed for the function

Input PDFs and their uncertainties
Four PDF sets were used in our study, the bold font denoting the shorthand name used to identify them in the rest of this article: • HERA20. The NNLO extraction by the H1 and ZEUS collaborations presented in ref. [27] with Hesse-like error band. The LHAPDF entry is HERAPDF20_NNLO_VAR with id = 61230.
• MSHT20. The NNLO extraction by the MSHT collaboration presented in ref. [30] with Hesse error band. The LHAPDF entry is MSHT20nnlo_as118 with id = 27400.
The comparison of these PDF sets for u-and d-quarks at scale µ = 2 GeV is presented in fig. 1.
In what follows the analyses are done using Bayesian statistics, which requires representing the PDFs as Monte-Carlo (MC) ensembles. As the NNPDF31 set is already given in this form, no further pre-processing is required. The other three distributions have a Hessian definition of uncertainty bands. The corresponding MC ensembles are generated using the prescription given in ref. [92]. Namely, for a distribution f (x) with 68% C.I. for each eigenvector given by f ± i (i = 1, . . . , D, with f + i and f − i the upper and lower bounds, respectively), a MC replica is generated by where R (k) i is a univariate random number. Using this method we generated 1000 MC replicas for HERA20, CT18 and MSHT20, which were used in the analyses. We have checked that the uncertainty bands obtained from the generated MC replicas are almost identical to the original Hesse uncertainty bands. The comparison of uncertainty bands for different PDF sets is also shown in fig.1.

Determination of uncertainties & presentation of the result
Within Bayesian statistics, the propagation of uncertainties to the free parameters is made by fitting each member of the input ensemble. Our input ensemble is generated accounting for two independent sources of uncertainty: EXP: The experimental uncertainty is accounted for by generating pseudo-data. A replica of the pseudo-data is obtained adding Gaussian noise to the values of the data points (and scaling the uncertainties if required). The parameters of the noise are dictated by the correlated and uncorrelated experimental uncertainties. The procedure is described in ref. [90]. We considered such 100 replicas.

PDF:
The uncertainty due to the collinear PDF is accounted for by using each PDF replica as input. We considered 1000 replicas generated as described in sec. 4.2.
The number of pseudo-data replicas is lower because the resulting uncertainty is much smaller than the one coming from the PDF distribution, as demonstrated in the following.
Fitting each member of the ensemble, we end up with the two-dimensional set λ ij , where λ is a 11-dimensional vector of fitting parameters (including λ f 1,2 , λ 0 and c 0 ), i runs over the pseudo-data replicas, and j runs over the PDF replicas. This set is distributed in accordance to the experimental and PDF uncertainties propagated through our fitting procedure. Due to the computational limitations mentioned in subsec. 2.2, calculating the full distribution with 100 × 1000 = 10 5 members is unrealistic. To simplify the task, we consider two distributions. The first one, labeled EXP, when the replicas of the data are fitted with the central PDFs. The second one, called PDF, when the replicas of the PDFs are fitted with the central (original) data. Symbolically, where i = 0 indicates the replica with unmodified data, and j = 0 a replica computed with the central value of the PDF. The PDF and EXP distributions are treated as independent. In the PDF case the distributions are also notably non-Gaussian. Therefore, we estimate the 68%C.I. for a parameter using the bootstrap method. Namely, we compute the [16%, 84%] quantile for a large number of samples and take the average interval. The results are presented as a δa1 δa2 , where a is the mean value of the parameter and δa's are the distances to the 68% C.I. boundary.
A drawback of performing the fit independently for each case is the appearance of several issues when wanting to join the two forms in a single meaningful one. The main problem is that, for some parameters, the distance between the mean values of the distributions is large. This renders impossible any naive description of a joined distribution; for example, the average value of the extracted TMDPDFs would have a large χ 2 , meaning that the corresponding TMDPDFs would not provide an adequate description of the data. Lesser, but not unimportant, problems are the correlations between individual points of TMDPDFs, and the PDF bias of the EXP case. Therefore, finally we present the central value and the 68% C.I., computed as described below.
The central value of our fit is obtained as the TMDPDF computed with the central values of the PDFs and NP parameters. The latter are computed as the weighted average of the PDF and EXP cases, where i = PDF or EXP, and σ i = (δa 1i − δa 2i )/2 is the half-size of the 68% C.I.. In this way, the central value incorporates information from both cases but retains the knowledge of which case gives a better determination of the considered parameter. The final parameters obtained by this procedure for each PDF set are shown in tab. 5. Also, the central distribution so defined returns a reasonable value of χ 2 , which should be similar to the result that we would obtain if we had performed one joined fit. The joined 68% C.I., instead, is computed from the (equal weight) Table 3: Distribution of the values of χ 2 for the central replica over the reduced data set in fits with different PDF inputs.
drawn from the merge of PDF and EXP and take the average interval. For those parameters for which PDF and EXP coincide, the width of the sampled distributions will be narrow, the contribution to the total uncertainty closely following the one of EXP (red band in fig. 3). For those parameters for which PDF and EXP do not significantly overlap, the width of the sampled distributions will be broader, the contribution to the total uncertainty closer, but narrower, to the one from PDF (green band in fig. 3). Overall, the bootstrapping method results in an uncertainty band that resembles the one we would expect if performing the fit with simultaneous replicas of the data and the PDFs. This procedure (rather than computation of average values of parameters) accounts for the correlation between different values of TMDPDFs.

Agreement between data and theory
The individual values of χ 2 for each experiment obtained for the reduced and complete data sets are given in tabs. 3-4, respectively. The tables also report the total χ 2 of the fits and they show that  Table 4: Distribution of the values of χ 2 for the central replica over the TMD data set in fits with different PDF input.
an overall reasonable description of the data is achieved with all PDF sets. Note that the ATLAS 7 TeV data have a generally higher χ 2 similarly to what is found also in other TMD fits [4,5,7,93]. The χ 2 distribution among the PDF and EXP replicas is shown in fig. 2. In general they are consistent with each other and their spread is highly reduced with respect to previous fits that use a flavor independent profile. More details on this are discussed in the appendix. This confirms the relevance of taking into account the flavor dependence of the NP TMD distributions f f NP . For all cases the PDF replicas provide a much larger dispersion of the χ 2 values than the EXP ones. The shapes of χ 2 -distributions are visibly different for different collinear PDFs, despite the resulting uncertainties in cross-section and TMDPDFs being similar. This may be due to the complexity of the corresponding parameterizations. In particular, the NNPDF and CT18 cases have more disturbed shapes of individual replicas (in comparison to HERA20 and MSHT20 cases) and consequently a larger spread of the χ 2 distribution.
The quality achieved in describing the data within the present fit is illustrated in figs. 3-4, for different PDF cases. Given the variety of experiments included in the fit and the number of PDF sets used, we present here only a fraction of all obtained results. We refer the interested reader to   Left panel: the ratio dσ experiment /dσ theory for Z-boson production at 8 TeV measured by the ATLAS experiment with MSHT20. Right panel: the ratio dσ experiment /dσ theory for Z-boson production at 13 TeV at the CMS experiment with NNPDF3.1. The red band is the EXP-uncertainty. The light-green band is the PDF-uncertainty. The blue band is the combined uncertainty. Only the filled bullets are included into the fit. the supplementary material for the complete collection of plots. The left panel of fig. 3 illustrates the Z-boson production measured by ATLAS at √ s = 8 TeV [83], the most precise data set in our analysis. On the right panel of fig. 3, we show the comparison of the theory prediction to data for the Z-boson production measured by CMS at √ s = 13 TeV [86]. Notice that these data were not included in the fit. As an example of the lowest energy measurements considered in the analysis, we present in fig. 4 the DY-process measured by the E288 experiment [73]. In all cases, the PDFuncertainty is larger than the EXP-uncertainty. This is especially pronounced for the high-energy measurements, for which the experimental uncertainty is small, and thus the uncertainty band is dominated by the PDF error. For q T 10 GeV, the EXP-uncertainty becomes negligible. This is the resummation regime, insensitive to NP TMD effects. In a few cases (central rapidity Z-boson production with CT18, and the lowest energy bins for DY process with NNPDF3.1), the PDF and EXP bands do not overlap, illustrating an unrevealed tension in the fitting procedure.   Figure 4: Example of the data description at low energy. Left panel: ratio dσ experiment /dσ theory for the DY process at E288 experiment with 200 GeV beam-energy with CT18. Right panel: ratio dσ experiment /dσ theory for the DY process at E288 experiment with 400 GeV beam-energy with MSHT20. Red band is the EXP-uncertainty. Light-green band is the PDF-uncertainty. The blue band is the combined uncertainty. The filled bullets are included into the fit. The dashed red vertical lines illustrate the cut q T < 0.25Q discussed at the beginning of sec. 3.
To better appreciate the role of the PDF and EXP uncertainty bands in figs. 3-4, we next consider the theoretical uncertainty bands obtained by variation of the perturbative scales. We perform the scale variation according to the ζ prescription approach in [6,37]. This amounts to varying two scales, the factorization scale in the DY cross section formula and the small-b matching scale in the OPE expansion of the solution of TMD evolution equations. (The small-b matching scale of the CS-kernel is present in this approach but its variation is not included in the calculation.) We vary these scales by factors c in the range [0. 5,2], and take the maximum symmetrized deviation. The resulting bands are shown in fig. 5.
We observe that the PDF uncertainty bands in fig. 3 are comparable or larger than the perturbative scale variation bands in fig. 5. This underlines that DY transverse momentum measurements in the TMD region are potentially useful to place constraints on PDFs. For this purpose one needs to employ a theoretical framework capable of describing the low transverse momentum region. One could envisage doing this in a resummation framework, formulated in terms of collinear PDF only, or in a TMD framework, in which a joined fit of both PDFs and TMDPDFs will put extra constraints on the PDFs.

Extracted TMD distributions
The values and error bars of the fitted parameters for the TMD distributions and CS kernel are given in tab. 5 and plotted in fig. 6. There, we report the bands obtained from the estimation of the PDF error (blue), the EXP error (red), and the averaged result (black). The central values for the PDF, EXP and joined cases do not coincide because they are computed with the respective replicas as explained in sec. 4. We observe from fig. 6 that the PDF error (blue) is generally much larger than the EXP error (red). That is, the error due to a single PDF-set uncertainty is always the most significant. This aspect can be relevant for future PDF analyses and should be taken into account by future TMDPDF fits.
We also see from fig. 6 that there is a general agreement among parameters within error bands. However, each PDF case has a few parameters whose values deviate significantly from the rest. These are {λ d 2 , λ s 1,2 } for CT18, {λ u 1 , λd 2 } for NNPDF3.1, {λ u 2 , λū 2 } for HERA20, and λū 1 for MSHT20. This highlights the tension in the corresponding domain between PDF and TMDPDF extractions. This may also be related to the larger values of χ 2 found in flavor-independent fits of TMDPDFs.
The fits in refs. [4,5,7,24] have found deficits in the predicted cross sections compared to measurements, with the deficits being small in the case of high-energy measurements (typically, 1-4%) but more significant for low-energy measurements (typically, 30-40%). In most cases, the discrepancy in the normalization is compensated by the correlated uncertainty of the experiment (e.g., E288 measurements have 25% uncertainty due to the beam luminosity), and thus does not significantly increase the value of χ 2 . Each PDF set requires its own normalization factor, e.g. for the central rapidity Z-boson production measured at ATLAS 8 TeV, the deficits are {3.6, 1.4, 2.1, 6.4}% for {MSHT20, HERA20, NNPDF3.1, CT18} cases. Another example is E288 at 200 GeV with {30, 49, 40, 39}%, correspondingly.
The actual shapes of the extracted TMDPDF are summarized in figs. 7, 8 and 9. We plot the "optimal" [4,24] TMD distribution, that is, the distribution defined according to the ζ-prescription [6,37] as the reference, scaleless TMD distribution. Fig. 7 provides the overall picture as a function of b and x. A more detailed view is offered by the slice in b in fig. 8 and the slice in x in fig. 9, showing the TMDPDF for each PDF set divided by the averaged central value of all PDF cases. The size of the uncertainty varies strongly from the low-b region, where we have a good knowledge

MSHT20
HERA20 NNPDF3.1 CT18 Figure 6: Comparison of the parameter values. Black is the final result. Blue is the value from the fit of the PDF case. Red is the value from the fit of the EXP case.
of the perturbative expansion, to the non-perturbative high-b region ( fig. 8). For b ≥ 2 GeV −1 the relative uncertainty on the TMD distributions is not less than 60-80%. For comparison, figs. 8 and 9 also show the uncertainty band obtained in the fit [4], labelled SV19. One of the main outcomes of the present work is that, compared to previous DY and SIDIS fits such as [4], the TMD uncertainty obtained in this paper is about 4-5 times larger, as a result of the improved analysis framework taking into account the propagation of collinear PDF uncertainties to the TMD extraction and the flavor dependence of the TMD profile. Furthermore, note that the present extraction has a non-zero uncertainty band also at b = 0, which was instead forbidden by construction in all previous studies.
The results presented in this section indicate that, compared to previous DY and SIDIS fits, the approach of the present paper leads to a more reliable estimate of the TMD uncertainties, to a reduced spread in the χ 2 distribution for each PDF set, and to a better agreement between different PDF sets. Nevertheless, we see from figs. 8 and 9 that the TMDPDFs extracted with different PDFs still display significant differences. A similar remark applies to the CS kernel: results for the CS 1.
1.  fig. 10. To further investigate the variation of TMDPDFs with the PDF set, we turn to Mellin moments of the TMDPDFs. We consider the following quantities where we set x min = 10 −5 . One may expect that the dependence on the collinear PDF is reduced, for n = 1 and n = 2, due to momentum sum rules fulfilled by each PDF set. Note however that f NP in eq. (2.4) contains a (weak) dependence on x. At any rate, the investigation of moments is significant because ratios for different flavours can be measured independently using lattice methods, see e.g. refs. [54,94]. The ratios of the first moments for the light quarks with respect to the u-quark are shown in fig. 11. We see that, with respect to point-to-point comparisons, they display a better global agreement. A similar behaviour appears for the ratios of second and third moments. The corresponding plots are given in the supplementary material. We conclude this subsection by noting that the results for unpolarized TMD distributions presented above will also be important for the phenomenology of spin-dependent TMD distributions, where the unpolarized cross-section serves as the normalization for the measurements of polarization asymmetries and angular distributions. For example, the uncertainty band from the fit [4] resulted into 10-15% of the total uncertainty band for the Sivers function [25]. With the results of the present paper one can, on one hand, expect an extra ∼ 50% uncertainty due to the inclusion of the PDF uncertainty. On the other hand, one can also expect a better agreement with experiments, due to a proper flavor dependence in the unpolarized sector.

Predictions for W -boson production
The flavor dependence of the TMD distributions could significantly impact the description of processes mediated by the W -boson [98][99][100]. In this work we have demonstrated that flavor dependence is essential to improve the agreement between theory and data. Including W -boson production data in future studies of unpolarized TMDPDFs will thus be very important.
Since it significantly slows down the fitting procedure, we have not included W -boson production data in this study. Instead, we compare the prediction made with the present extraction to the Parameter MSHT20 HERA20 NNPDF31 CT18  Table 5: The values of the NP parameters obtained in the fit.     Figure 12: Comparison of the differential cross section for W-boson production measured by ATLAS at 7 TeV with the theory prediction using different PDFs and normalized to HERA20 case.
one made in [39]. The computation is performed with the artemide code, appropriately adapted for the computation of the transverse-mass-differential cross-section as described in ref. [39]. The values of the χ 2 for comparison with the data by Tevatron and LHC experiments [85,[95][96][97] are presented in tab. 6. We find a small general improvement compared to ref. [39], which is expected because most of W -boson production data belong to the resummation region. We see that the χ 2 values for the CMS measurement are large in comparison with the other data considered. This may be understood by noting that this measurement is integrated over the full range of the dilepton mass, and the integration range covers low-mass regions in which power corrections in q T /Q are non-negligible. As a result, the theory prediction deviates from the measurement starting from smaller values of q T , producing larger χ 2 -values. This effect is partially canceled in the cross-section ratios W − /W + and Z/W , which show a substantial agreement between theory and the CMS measurements. Additional plots on this are presented in the supplementary material.

Conclusions
In this work we have performed the first quantitative study of the influence of PDF on fits to DY transverse momentum measurements based on TMD factorization. We have referred to this as the PDF bias issue, arising from the fact that an OPE is applied to the TMD operator and an ansatz is made for the TMD distribution in terms of collinear PDFs. We have highlighted the quantitative significance of this issue in unpolarized TMD fits. We have stressed that this is relevant also in polarized TMD analyses, for which unpolarized cross sections are used to normalize angular distributions and asymmetries.
We have carried out a Bayesian procedure to propagate PDF uncertainties to the extraction of TMD parton distributions. We have examined four PDF sets (MSHT20, NNPDF31, CT18, HERA20), representative of different NNLO PDF methodologies, and we have performed a TMD determination from DY data including for the first time both experimental and PDF uncertainties. We have found that the PDF uncertainties are larger than the DY experimental uncertainties for all values of b (or p T ). As a result of the improved analysis framework, we have obtained reliable estimates of TMD uncertainties, with the size of TMD error bands being significantly increased with respect to previous TMD fits which do not perform the full PDF bias analysis.
We have included for the first time flavor dependence in the non-perturbative TMD distribution f f NP . Previous fits include flavor dependence in the collinear PDF only. We have found that flavordependent TMD profiles reduce the spread in χ 2 distributions for each PDF set, improving the agreement between data and theory, and help obtain more consistent results among different PDF sets.
The results of this paper indicate that including collinear PDF uncertainties in TMD extractions and taking into account the flavor dependence of NP TMD distributions are both essential to obtain reliable TMD determinations from DY (and SIDIS) transverse momentum data. Future phenomenological studies, which incorporate these features with more powerful computational and statistical tools than those used here, are warranted.

Acknowledgments
This study was supported by Deutsche Forschungsgemeinschaft (DFG) through the research Unit FOR 2926, "Next Generation pQCD for Hadron Structure: Preparing for the EIC", project number 430915485. I.S. is supported by the Spanish Ministry grant PID2019-106080GB-C21. This project has received funding from the European Union Horizon 2020 research and innovation program under grant agreement Num. 824093 (STRONG-2020). S.L.G. is supported by the Austrian Science Fund FWF under the Doctoral Program W1252-N27 Particles and Interactions.

A Appendix
In this appendix we provide a few more details and checks on the fits performed in sec. 5.
In fig. 2 we have presented our results for the distribution of χ 2 values among the PDF and EXP replicas, and we have observed that the χ 2 spread is much reduced with respect to previous fits in the literature, thanks to the flavor-dependent profile used in the present work for the NP TMD distributions f f NP . Here we illustrate the role of the flavor dependence explicitly by reporting the result for the χ 2 distribution which we obtain by repeating the calculation of fig. 2 but replacing the flavor-dependent f f NP model of eq. (2.4) with the flavor-independent model of ref. [4]. The result is shown in fig. 13 for the case of the NNPDF3.1 PDF set [28].
The result in fig. 13 is to be compared with the third panel in fig. 2. We see that, in contrast to the χ 2 distribution found in fig. 2, the distribution of χ 2 -values over PDF replicas in fig. 13 shows an unsatisfactorily broad shape, with about 64% of the replicas having χ 2 /N pt > 2. We have checked that the unsatisfactory χ 2 -values for the subset of replicas are not due to a single problematic measurement but rather they are common to all data.
We have also checked that the issues described above are not specific to the NNPDF3.1 PDF set [28]. We performed similar tests using HERA20 [27]    [29] (χ 2 0 /N pt = 1.26), and CJ15nlo [104] (χ 2 0 /N pt = 1.82), where χ 2 0 is the χ 2 -value for the central PDF replica. All these PDF sets are characterized by the same issues as NNPDF3.1. This confirms that the essential element at the origin of the difference between the χ 2 distributions in fig. 2 and  fig. 13 is the flavor dependence of the NP TMD distributions f f NP . We next discuss the effect of different PDF replicas on the shape of the predictions for the transverse momentum distribution. One might wonder whether the change in PDF replicas results into an effect primarily on the normalization but not on the q T shape of the predictions. In fig. 14 we illustrate that this is not the case. That is, fig. 14 indicates that the large spread in the χ 2 distribution observed above is due to different PDF replicas inducing different q T -shapes of predictions. The variety of shapes is a consequence of the structure of the convolution within OPE, which correlates the b and x dependences.