Deuteron Uncertainties in the Determination of Proton PDFs

We evaluate the uncertainties due to nuclear effects in global fits of proton parton distribution functions (PDFs) that utilise deep-inelastic scattering and Drell-Yan data on deuterium targets. To do this we use an iterative procedure to determine proton and deuteron PDFs simultaneously, each including the uncertainties in the other. We apply this procedure to determine the nuclear uncertainties in the SLAC, BCDMS, NMC and DYE866/NuSea fixed target deuteron data included in the NNPDF3.1 global fit. We show that the effect of the nuclear uncertainty on the proton PDFs is small, and that the increase in overall uncertainties is insignificant once we correct for nuclear effects.

Parton distribution functions (PDFs) are an essential ingredient in the theoretical predictions of hadronic observables at the LHC [1][2][3]. PDFs for the proton are determined via global QCD fits to a range of experimental data, including those where the proton is not in a free state. In particular, these include deep-inelastic scattering (DIS) and Drell-Yan (DY) fixed target collisions involving deuterium and heavy nuclear targets. In these processes the interaction of the proton is altered due to nuclear effects, and this difference propagates through to the fitted PDFs. Measurements involving deuterium targets still play a significant role in the determination of proton PDFs, in particular to separate the up and down flavours for large momentum fraction x, a region which is especially important for searches for physics beyond the Standard Model. Because of this, deuteron corrections have been extensively studied and have been included in PDF analyses via parametrizations of a nuclear smearing function [4][5][6][7][8], inspired by various deuteron wavefunction models [9][10][11][12][13]. This approach relies on model assumptions, which can ultimately bias the determination of the PDFs in a way which is difficult to quantify. Because the precision of the PDFs is now constrained by the data to a few percent for most quark flavours in a wide kinematic range [14], a faithful estimate [15] of the theoretical uncertainty associated with nuclear effects (potentially of comparable size) is becoming necessary.
In a previous study [16] we showed how theoretical uncertainties due to heavy nuclear targets in DIS and DY measurements can be incorporated into global fits of proton PDFs. Specifically, in the framework of the NNPDF methodology (see [17] and references therein for a comprehensive description), we added to the experimental covariance matrix a theoretical covariance matrix, accounting for the additional uncertainties due to nuclear effects. Two distinct procedures were adopted: in the first, the contribution of the nuclear data to the PDF fit is deweighted by an uncertainty that encompasses both the difference between proton and nuclear PDFs and the uncertainty in the nuclear PDFs; in the second, the difference between proton and nuclear PDFs is used to correct the theoretical predictions, while the deweighting only takes the nuclear PDF Figure 1. A schematic representation of the iterative procedure adopted to determine the uncertainty due to deuteron corrections in proton PDF fits, see text for details. The global dataset is the union of the proton and deuteron datasets. uncertainty into account, and is therefore correspondingly smaller. If the uncertainty in the nuclear PDFs is correctly estimated, and smaller than the shift, the second procedure should give more precise results. The nuclear PDFs were determined as an equally weighted replica average of the DSSZ [18], nCTEQ15 [19], and EPPS16 [20] PDF sets for the relevant heavy nuclei (Cu, Fe and Pb). Despite the fact that these are obtained from a global analysis of experimental data taken in a wide variety of processes, there are sizeable differences between them. This suggests that these three sets might not be sufficiently consistent to determine a precise nuclear correction, but can be used to estimate the uncertainty due to nuclear effects, and indeed the second procedure led to a worse global fit than the first.
In this paper, we extend the approach developed for heavy nuclei in [16] to deuterium data. We focus on the dataset included in the NNPDF3.1 PDF determination [21], which is made up of 3978 data points (see [14] for details). Out of these, 418 data points (about 10% of the whole dataset) come from experiments using deuterium targets, specifically SLAC [22], BCDMS [23], NMC [24], and DYE866/NuSea [25]. The DIS data are in the form of deuteron to proton structure function ratios, F d 2 /F p 2 , for NMC, and of deuteron structure functions, F d 2 , for SLAC and BCDMS; the DY data is in the form of ratios of cross sections for a proton beam on a deuteron target to a proton beam on a proton target, σ DY pd /σ DY pp , for DYE866/NuSea. A significant weakness in our treatment of nuclear effects in heavy nuclei was its dependence on externally determined nuclear PDFs. To avoid this when treating deuterium, we fit our own deuterium PDFs directly from the deuteron data by means of a procedure which is iterated to consistency with the global proton fit. In this way we account simultaneously for the nuclear uncertainties in the deuteron when determining global proton PDFs, and the uncertainties in the proton PDF when determining the deuteron PDF (and thus the nuclear correction). The advantage of the new approach is that the deuteron and proton fits are all performed using a consistent theoretical and methodological fitting framework, and the resulting nuclear corrections and their uncertainties are thus equally reliable and unbiased. It is very similar to the self consistent procedure set out in [15] for the simultaneous determination of PDFs and fragmentation functions using experimental data from semi-inclusive DIS.
The logical structure of the procedure is depicted in Fig. 1. The NNPDF3.1 global dataset is split into two disjoint subsets: 'deuteron data', including the aforementioned datasets (SLAC, BCDMS, NMC and DYE866/NuSea); and 'proton data' including all the other datasets. The global dataset is the union of deuteron and proton datasets. Note that the proton data also includes CHORUS [26], NuTeV [27], and DYE605 [28] data taken in experiments on heavy nuclei. The effect of nuclear corrections on these measurements was studied in [16] and is not considered here as we want to focus exclusively on the nuclear effects in deuterium.
The deuteron data themselves may be split into two sets: 'pure' deuteron data, from DIS on a deuteron target (the SLAC and BCDMS data for F d 2 ), and 'mixed' deuteron data, which also involve protons (the NMC data for the ratio F (1) Here f n is the neutron PDF, determined from the proton PDF by assuming exact isospin invariance (and thus in practice by swapping the up and down PDFs).
The proton data are used to determine a first set of pure proton PDFs {f (k) p : k = 1 · · · N rep }, using the usual NNPDF methodology. The central prediction is f , where the angled brackets denote a simple average over the N rep Monte Carlo replicas. If the deuteron data were all pure deuteron data, in practice only the SLAC and BCDMS datasets, we could simply produce a similar set of pure deuteron PDFs {f (k) d : k = 1 · · · N rep }. These by construction will include the nuclear effects and the size of the nuclear correction would be T s ], with f (0) s determined from the proton PDFs using Eq. (1) averaged over proton PDF replicas. To include the mixed deuteron data, in particular the data from NMC and DYE866/NuSea, we have to be more careful since to evaluate the theoretical predictions for a given deuteron PDF we also need a proton PDF. For this we can use the central value of the pure proton fit f p ]. However we must also include the uncertainty in the proton fit, as part of the theoretical (proton) uncertainty in determining the deuteron PDF from these data: this can be done by computing the theory covariance matrix where i, j run over the data points in the mixed deuteron datasets only. Note that this covariance matrix incorporates correlations between the mixed datasets due to their common dependence on the proton PDF. This theoretical covariance matrix is added to the experimental covariance matrix of the mixed deuteron data when performing the deuteron fit, to take account of the uncertainty of the proton PDF in the determination of the deuteron PDFs. Note that the theory covariance matrix Eq. (2) itself depends on the deuteron PDF, which is what we are trying to determine. However the dependence of the fitted deuteron PDF on the uncertainty in the proton PDF is relatively weak, since it only affects the weight of the mixed data in the fit. Thus to a good approximation we can replace f s , determined from the pure proton fit. For a more accurate determination of the deuteron PDFs, we could then iterate to consistency, performing a second fit to the deuteron data where f (2) is determined from the first fit. It is clear that this iterative process would converge very rapidly.
However our aim here is not so much to determine the deuteron PDF, but rather to use it to determine a theoretical covariance matrix that takes into account the nuclear effects in the deuteron data (both pure and mixed) when using these data in a global fit of the proton PDF. Since the size of the nuclear correction is given by the difference between predictions with deuteron and isosinglet PDFs, this theoretical (deuteron) covariance matrix is Again this covariance matrix incorporates correlations between all the nuclear corrections in the various deuteron datasets, due to their common dependence on the deuteron PDF. To perform a global fit of the proton PDF including nuclear uncertainties in the deuteron data, we can simply add the theory covariance matrix to the experimental covariance matrix of the deuteron datasets and perform the proton fit in the usual way; this yields a set of replicas of proton PDFs {f Since the global proton PDFs will be more precise than the pure proton PDFs we started with, it makes sense once again to iterate; we use our global proton PDFs to determine an improved theoretical proton covariance matrix Eq. (2), repeat the deuteron fit, use this to determine an improved theoretical deuteron covariance matrix Eq. (3), and then use this to perform a new global fit of the proton PDF. This is the iterative procedure shown schematically in Fig. 1. Note that through this procedure the deuteron PDF is also iterated concurrently. We expect the iterations to converge very rapidly to a self consistent set of deuteron and (global) proton PDFs for several reasons: firstly, a small change in the proton PDF makes a small difference to the deuteron correction; secondly, we expect the effect of the deuteron correction on the weight of these data in the global fit to be small; thirdly, the influence of the deuteron data in the global fit is already relatively small (just as the influence of the proton PDF on the deuteron fit is small). Note that the deuteron data are not double counted in this procedure; in the deuteron fit they are used to determine (empirically) the nuclear uncertainty, while in the global fit they influence the central value of the proton PDF directly, but taking into account this nuclear uncertainty. Indeed, the nuclear uncertainty reduces the weight of the deuteron data in the global fit, so they actually count less. As a byproduct, we also determine a set of deuteron PDFs.
This realises the first of the two procedures described in [16], whereby the deuteron datasets are deweighted by the nuclear uncertainty but theoretical predictions are not shifted by a nuclear correction. As in [16], we also implement the second procedure: in this case, the theoretical (deuteron) covariance matrix is defined as while the corrections applied to the theoretical predictions p ] for the deuteron datasets are This procedure is implemented in the same way as the first, with the theoretical predictions for the deuteron data corrected before performing each iteration of the global proton fit. Since, unlike in [16], the corrections are determined empirically and self consistently, we expect the second method to be more precise than the first method; the central values of the theoretical predictions should be a little more accurate, and the uncertainty due to nuclear effects in the deuteron correspondingly a little smaller. The set of fits which we performed, all using NNPDF methodology, are summarised in Table 1. They are all accurate to next-to-next-to-leading order (NNLO) in perturbative QCD, heavy quarks are treated in the FONLL scheme and the charm PDF is parametrised in the same way as the lighter quark PDFs. All PDF sets are made of N rep = 100 Monte Carlo replicas. The baseline fit, 'global-base', is a fit equivalent to the base fit performed in [29]. It is a minor variant of the fit presented in [21]: a bug affecting the computation of theoretical predictions for charged-current DIS cross sections has been corrected; positivity of the F c 2 structure function has been enforced; and NNLO massive corrections [30,31] have been included in the computation of neutrino-DIS structure functions. 'proton-ite0' is the corresponding fit based on the proton data alone. We then perform two iterations of the procedure described above (denoted as iteration 1

Dataset Fit ID Description
Baseline Proton and Deuteron global-base Same as base fit in [29].
Iteration 0 Proton proton-ite0 Same as baseline, but restricted to the proton dataset Iteration 1 Deuteron deuteron-ite1 Same as baseline, but restricted to the deuteron dataset and supplemented with a proton covariance matrix determined from the proton-ite0 fit according to Eq. (2). Proton and Deuteron global-ite1-dw Same as baseline, but supplemented with a deuteron covariance matrix determined from the deuteron-ite1 fit according to Eq. (3).

Iteration 2
Deuteron deuteron-ite2 Same as deuteron-ite1, but with a proton covariance matrix determined from the global-ite1-dw fit. Proton and Deuteron global-ite2-dw Same as global-ite1-dw, but with a deuteron covariance matrix determined from the deuteron-ite2 fit. Proton and Deuteron global-ite2-sh Same as global-ite2-sh, but with a deuteron covariance matrix and shifts determined according to Eqs. (4)-(5). and 2), after which we determine a fit of deuteron PDFs (based only on the deuteron data, and supplemented with a proton covariance matrix), and a global fit of proton PDFs (based on the proton and deuteron data, and supplemented with a deuteron covariance matrix). The deuteron fits 'deuteron-ite1' and 'deuteron-ite2', are performed using exactly the same theoretical and methodological settings as the proton fits, except that the isotriplet PDFs are set to zero since the deuteron is isoscalar. After the first iteration we produce a single global fit of proton PDFs, 'global-ite1-dw', in which the deuteron covariance matrix is evaluated according to Eq. (3).
After the second iteration we produce instead two global fits of proton PDFs: 'global-ite2-dw' in which the deuteron covariance matrix is evaluated with Eq. (3), and 'global-ite2-sh', in which the deuteron covariance matrix is evaluated with Eq. (4), and the theoretical predictions are first corrected according to Eq. (5). Before discussing the results of our fits, we look more closely at the pattern of deuteron corrections, and at the deuteron covariance matrix defined in Eq. (3). As representative examples, results are obtained from proton and deuteron PDFs determined after the first iteration. We explicitly checked that they remain stable after an additional iteration.
In Fig. 2 we display the nuclear correction for the deuteron data obtained from our procedure. Specifically, for each data point i (after kinematic cuts), we show the observables computed with the central deuteron PDF, normalised to the expectation value computed with the central proton Data points are ordered in bins of increasing values of momentum fraction x and energy Q. The deuteron correction generally amounts to a few percent for all of the experiments considered. Uncertainties are rather large and the ratio is mostly compatible with one, except for data points at higher values of momentum fraction x and energy Q, where the correction is negative, as expected from models of nuclear shadowing.
To gain a further idea of the effects to be expected from nuclear corrections in deuteron, in  Fig. 2 is paralleled in Fig. 3, in particular concerning the dependence of the size of the nuclear uncertainties on the bin kinematics for each experiment. Moreover we can now see that the deuteron uncertainties are smaller than the data uncertainties for SLAC and BCDMS, while they are comparable for NMC and DYE866/NuSea. This is largely because the pure deuteron measurements from SLAC and BCDMS are of cross-sections, whereas the more precise mixed measurements from NMC and DYE866/NuSea are of cross-section ratios, for which systematic uncertainties largely cancel.  Finally, we show the experimental correlation matrix, ρ C ij = C ij / C ii C jj , and the sum of the experimental and deuteron correlation matrices, ρ C+S ij = (C ij + S ij )/ (C ii + S ii )(C jj + S jj ), as heat plots in Fig. 4 . The theoretical covariance matrix is computed according to Eq. (3), though the qualitative behaviour of the total correlation matrix is unaltered if Eqs. (4)-(5) are used instead. Our procedure captures the sizeable correlations of the deuteron corrections between different bins of momentum and energy, systematically enhancing bin-by-bin (positive and negative) correlations in the data. As we might expect, nuclear uncertainties are also strongly correlated across the different experiments.
We now turn to discuss the results of the fits collected in Table 1. In Tables 2-3 we display the values of the experimental χ 2 per data point (as defined in Eq. (4) of [16]) for the fits of the deuteron PDFs (based on the deuteron data) and of the proton PDFs (based on the global dataset including both deuteron and proton data), respectively. In Table 3 values are displayed both for separate datasets, and for groups of datasets corresponding to measurements of similar   Table 2. The values of the χ 2 per data point for each dataset included in the fits of deuteron PDFs. observables in the same experiment. Indented datasets are subsets of the preceding non-indented dataset.
In order to examine the convergence of our procedure, we must quantifiy the statistical equivalence between pairs of PDFs obtained from the various fits. To this purpose we display in Fig. 5 the distance (as defined in Eq. (63) of [32]) between the central values of the two iterations of fits based on the deuteron data (deuteron-ite1 and deuteron-ite2), and the corresponding two iterations on the global data (global-ite1-dw and global-ite2-dw). For two PDF sets made of N rep = 100 replicas, a distance of d 1 corresponds to statistically equivalent sets, while a distance of d 10 corresponds to sets that differ by one sigma in units of the corresponding standard deviation. Note that in the left panel of Fig. 5 u andū actually denote the combinations (u + d)/2 and (ū +d)/2, where u = d andū =d by definition. These are the isosinglet combinations determined in the fits to deuteron data.
From the results displayed in Tables 2-3 and in Fig. 5, we can conclude that one iteration is sufficient to achieve stability. The variation of the global χ 2 per data point for fits obtained in subsequent iterations is smaller than statistical fluctuations. This is true both in the case of fits of deuteron PDFs (the global χ 2 per data point is 0.97 and 0.98 for the deuteron-ite1 and deuteron-ite2 fits, respectively), and in the case of global fits of proton PDFs (the global χ 2 per data point is 1.16 in both the global-ite1-dw and global-ite2-dw fits), see Tables 2-3 Table 3. The values of the χ 2 per data point for each dataset included in the global fits of proton PDFs. The deuteron data are at the top of the table.  Table 1 for details. For the deuteron fits, u andū actually denote the combinations (u + d)/2 and (ū +d)/2, where u = d andū =d by definition. Results are displayed as a function of x at a representative scale of the deuteron dataset, Q = 10 GeV. The ReportEngine software [33] was used to generate this figure.
upon iteration are u,ū, d,d in the valence region of the deuteron fit, as might be expected. Next, in Fig. 6 we compare the deuteron PDFs obtained from each iteration (deuteron-ite1 and deuteron-ite2). Specifically we show the average of up and down quark, the average of up and down antiquark, the strange quark, and the gluon distributions: because the deuteron is isoscalar, d = u andd =ū by construction. Again we see that the PDFs hardly change from one iteration to the next, so the procedure has converged. In addition in this plot we compare our NNLO deuteron PDFs with a recent NLO determination of nuclear PDFs based on the NNPDF methodology, nNNPDF2.0 [34]. These were obtained by fitting a range of nuclear and proton data, and assuming a smooth dependence on the mass and atomic numbers A and Z. Due to this assumption, which in effect constrains the deuteron as an interpolation between proton and heavy nuclei, their uncertainties are smaller than our own. However their central values are mostly consistent with ours within uncertainties. A discrepancy of about one sigma, in units of the uncertainty of the deuteron-ite2 fit, is observed for the average of up and down quarks around x ∼ 0.1. Such a discrepancy might be explained in light of the fact that the two determinations are at different orders in perturbation theory. Our determination of the deuteron PDFs is thus very conservative.
To explore the nuclear corrections further, in Fig. 7 we show the ratio F d 2 /F p 2 computed with our deuteron PDFs, at Q = 10 GeV. We see that the correction for nuclear effects in deuteron is only a few per cent over the full range of x, and is negative in the valence region, as expected from nuclear shadowing. However the uncertainty in our determination is as large as the correction. For comparison, we also show the same quantity computed using the nNNPDF2.0 deuteron PDFs [34]: these are NLO, but have a smaller uncertainty since, as explained above, in these fits continuity in A/Z is implicitly assumed, which adds a significant constraint. However the reduction in uncertainty due to this constraint is considerably less in the structure function ratio than it was in the PDFs. We also show the parametric correction used in the MMHT14 fits [6], which has four fitted parameters. Again this has a yet smaller uncertainty, particularly at large x, due to the assumed theoretical constraints of the model. However we note that all three of these estimates are mutually consistent, within uncertainties. The estimate obtained here is clearly the most conservative, particularly outside the valence region, as expected since it is free from any model dependence.
Finally we consider the impact of the deuteron uncertainties in the global fit of proton PDFs. We compare the global fits (to the deuteron and proton data), made without the inclusion of the theoretical covariance matrix (global-base), and then with the inclusion of the theoretical covariance matrix after the second iteration, either with Eq. (3) (global-ite2-dw) or with Eq. (4) and the associated shifts, (5) (global-ite2-sh). From Table 3 we conclude that the nuclear corrections give a small improvement in the overall fit quality (the global χ 2 per data point is reduced from 1.18 in the global-base fit to 1.16 in the global-ite2-dw and global-ite2-sh fits, which corresponds to one standard deviation of the χ 2 distribution), and a significant improvement in the fit quality of the deuteron datasets.
Turning to the PDFs themselves, in Fig. 8 we compare the proton PDFs obtained from these three fits. Here we show only the up and down quark and antiquark PDFs, normalised to the global-base fit, and the corresponding relative uncertainties, since the other quark flavours and the gluon PDFs are only scarcely affected by the nuclear corrections in deuteron. We also show, in Fig. 9 the distances (defined as in Fig. 5) between the baseline fit (global-base) and each of the two global fits with deuteron uncertainties included after the second iteration (global-ite2-dw and global-ite2-sh). Again, a distance of d 10 corresponds to sets that differ by one sigma in units of the corresponding standard deviation.
The effect of the nuclear corrections on the PDF central values is largest in the up antiquark in the valence region: it differs by about half a sigma (d ∼ 5 in Fig. 9) in both the global-ite2-dw and global-ite2-sh fits with respect to the global-base fit. As apparent from Fig. 8, the central value of the up antiquark PDFs is suppressed in the valence region, while that of the down antiquark is enhanced. The effect is seen irrespective of whether the theoretical predictions are shifted. The inclusion of the nuclear uncertainty in the global fits of proton PDFs results in a slight increase in the uncertainties in comparison to the global-base fit, but this increase is rather larger in the global-ite2-dw fit than in the global-ite2-sh fit, where it is scarcely visble. This result, combined with the fact that both these fits have comparable quality (see Table 3), leads us to conclude that the shifted fit is to be preferred. This is as expected, given that the uncertainty due to nuclear corrections has been determined self-consistently, and turns out to be a little smaller in the valence region than the nuclear correction itself (see Fig. 2). This in contrast to the result we found in the case of heavy nuclei [16], for which nuclear uncertainties were instead estimated from independent global determinations of nuclear PDFs. Clearly the   Table 1 for details. Results are displayed as a function of x at a representative scale for the deuteron dataset, Q = 10 GeV. The ReportEngine software [33] was used to generate this figure.
self-consistency of our procedure is advantageous, and should therefore be preferred in the case of deuteron data (and for heavy nuclei whenever it is possible to perform a consistently reliable determination of the nuclear PDFs and their uncertainties).
In summary, we have developed an iterative procedure to incorporate theoretical uncertainties due to nuclear effects self-consistently into global fits of proton PDFs that include DIS and DY data on deuterium targets, without any model dependent assumptions regarding the physics of the nuclear corrections. In the framework of the NNPDF3.1 global analysis we have shown that the effect of the additional uncertainty in the global determination of the proton PDFs is small, and can be reduced further by applying an empirical correction to the theoretical predictions of the deuteron data. Such a fit thus leads to slightly more precise PDFs. A similar procedure might be used to improve the determination of kaon fragmentation functions, and thus the strange and anti-strange proton PDFs, by means of semi-inclusive DIS measurements that are sensitive to both. The PDF sets discussed in this work are available in the LHAPDF format [35] from the authors upon request.