Open-source QCD analysis of nuclear parton distribution functions at NLO and NNLO

TL;DR

This paper introduces new nuclear parton distribution functions at NLO and NNLO, derived from deep inelastic scattering data, using an open-source tool with modifications for nuclear analysis, and compares them to existing sets.

Contribution

It provides the first comprehensive NLO and NNLO nuclear PDFs using an open-source framework with modifications for nuclear data analysis.

Findings

New nPDF sets at NLO and NNLO levels.

Comparison with existing nPDFs shows consistency.

Uncertainty estimates highlight data limitations.

Abstract

We present new sets of nuclear parton distribution functions (nPDFs) at next-to-leading order (NLO) and next-to-next-to-leading order (NNLO). Our analyses are based on deeply inelastic scattering data with charged-lepton and neutrino beams on nuclear targets. In addition, a set of proton baseline PDFs is fitted within the same framework with the same theoretical assumptions. The results of this global QCD analysis are compared to existing nPDF sets and to the fitted cross sections. Also, the uncertainties resulting from the limited constraining power of the included experimental data are presented. The published work is based on an open-source tool, xFitter, which has been modified to be applicable also for a nuclear PDF analysis. The required extensions of the code are discussed as well.

Tables

Table 1. Table 1: Summary of experimental DIS data used to determine proton PDFs. In the last two columns the resulting χ 2 superscript 𝜒 2 {\large{\chi}}^{2} values at NLO and NNLO obtained in our analysis are provided.

Exp.	Data set	Year	Ref.	$N_{\mathrm{dp}}$	${\chi }^2$ NLO	${\chi }^2$ NNLO
BCDMS	F2p 100GeV	1996	Benvenuti et al. (1989)	83	88.88	90.98
	F2p 120GeV			90	69.97	67.75
	F2p 200GeV			79	89.46	85.91
	F2p 280GeV			75	66.97	68.73
HERA 1+2	NCep 920	2015	Abramowicz et al. (2015)	377	455.15	475.14
	NCep 820			70	72.47	73.84
	NCep 575			254	225.24	228.97
	NCep 460			204	223.23	223.95
	NCem			159	233.55	229.42
	CCep			39	42.19	44.41
	CCem			42	65.94	68.99
NMC-97	NCep	1997	Arneodo et al. (1997)	100	124.56	111.64
In total:				1559	1845.99	1909.08

Table 2. Table 2: Summary of experimental DIS data used to determine the nuclear PDFs. In the last two columns the resulting χ 2 superscript 𝜒 2 {\large{\chi}}^{2} values at NLO and NNLO obtained in our analysis are provided.

Nucleus	Exp.	Year	Ref.	$N_{\mathrm{dp}}$	${\chi }^2$ NLO	${\chi }^2$ NNLO
D	NMC 97	1996	Arneodo et al. (1997)	120	124.85	118.66
	EMC 90	1989	Arneodo et al. (1990)	21	29.23	31.73
He/D	HERMES	2002	Airapetian et al. (2002)	7	54.64	37.99
	NMC 95, re.	1995	Amaudruz et al. (1995)	13	12.44	12.98
	SLAC E139	1994	Gomez et al. (1994)	11	7.21	4.68
Li/D	NMC 95	1995	Arneodo et al. (1995)	12	7.06	5.93
Be/D	SLAC E139	1994	Gomez et al. (1994)	10	7.84	7.83
Be/C	NMC 96	1996	Arneodo et al. (1996a)	14	14.80	16.19
C	EMC 90	1989	Arneodo et al. (1990)	17	11.01	10.05
C/D	FNAL E665	1995	Adams et al. (1995)	3	5.12	5.91
	SLAC E139	1994	Gomez et al. (1994)	6	15.12	17.16
	EMC 88	1988	Ashman et al. (1988)	9	4.49	3.50
	NMC 95, re.	1995	Amaudruz et al. (1995)	13	38.08	36.52
C/Li	NMC 95, re.	1995	Amaudruz et al. (1995)	10	17.27	13.90
N/D	HERMES	2002	Airapetian et al. (2002)	1	2.20	0.97
Al/D	SLAC E139	1994	Gomez et al. (1994)	10	11.20	14.22
Al/C	NMC 96	1996	Arneodo et al. (1996a)	14	6.51	6.55
Ca	EMC 90	1989	Arneodo et al. (1990)	19	13.17	12.56
Ca/D	NMC 95, re.	1995	Amaudruz et al. (1995)	12	29.61	31.12
	FNAL E665	1995	Adams et al. (1995)	3	4431	6.01
	SLAC E139	1994	Gomez et al. (1994)	6	8.44	9.34
Ca/Li	NMC 95, re.	1995	Amaudruz et al. (1995)	10	7.36	5.16
Ca/C	NMC 95, re.	1995	Amaudruz et al. (1995)	10	6.47	6.70
	NMC 96	1996	Arneodo et al. (1996a)	14	7.14	6.99
Fe	SLAC E140	1993	Dasu et al. (1994)	2	0.05	0.05
Fe/D	SLAC E139	1994	Gomez et al. (1994)	14	34.08	34.18
Fe/C	NMC 96	1996	Arneodo et al. (1996a)	14	9.82	9.96
$\nu$ Fe	CDHSW	1991	Berge et al. (1991)	464	347.74	365.14
$\overline{\nu }$ Fe	CDHSW	1991	Berge et al. (1991)	462	423.06	398.25
Cu/D	EMC 93	1993	Ashman et al. (1993)	19	18.12	17.45
	EMC 88	1988	Ashman et al. (1988)	9	5.59	7.22
Kr/D	HERMES	2002	Airapetian et al. (2002)	1	2.02	2.02
Ag/D	SLAC E139	1994	Gomez et al. (1994)	6	16.24	18.81
Sn/D	EMC 88	1988	Ashman et al. (1988)	8	14.56	9.24
Sn/C	NMC 96	1996	Arneodo et al. (1996a)	14	12.90	7.61
	NMC 96, $Q^2$ dep.	1996	Arneodo et al. (1996b)	134	94.7	79.85
Xe/D	FNAL E665	1992	Adams et al. (1992)	3	2.13	2.53
Au/D	SLAC E139	1994	Gomez et al. (1994)	11	16.64	19.80
Pb/D	FNAL E665	1995	Adams et al. (1995)	2	12.24	13.32
Pb/C	NMC 96	1996	Arneodo et al. (1996a)	14	9.94	6.77
$\nu$ Pb	CHORUS	2005	Onengut et al. (2006)	405	229.11	243.85
$\overline{\nu }$ Pb	CHORUS	2005	Onengut et al. (2006)	405	361.35	328.28
In total:				2336	2072.29	2014.02

Table 3. Table 3: Values of the NLO fit parameters at the initial scale Q 0 2 = 1.69 ​ GeV 2 superscript subscript 𝑄 0 2 1.69 superscript GeV 2 Q_{0}^{2}=1.69\,\mathrm{GeV^{2}} . (SR) means that the normalization for that particular parton is fixed by the sum rules. A dash symbolizes that this parameter was excluded from the fit. Parameter values for the sea quarks, apart from u ¯ ¯ 𝑢 \bar{u} , were derived from the applied constraints s ¯ = s = d ¯ = u ¯ ¯ 𝑠 𝑠 ¯ 𝑑 ¯ 𝑢 \bar{s}=s=\bar{d}=\bar{u} .

$g$	value	$u_v$	value	$d_v$	value	$\overline{u}$	value
$c_{0,0}^g$	7.0352	$c_{0,0}^{u_v}$	(SR)	$c_{0,0}^{dv}$	(SR)	$c_{0,0}^{\overline{u}}$	(SR)
$c_{1,0}^g$	0.2871	$c_{1,0}^{u_v}$	0.6046	$c_{1,0}^{dv}$	0.7376	$c_{1,0}^{\overline{u}}$	-0.1915
$c_{2,0}^g$	14.243	$c_{2,0}^{u_v}$	3.7064	$c_{2,0}^{dv}$	2.9225	$c_{2,0}^{\overline{u}}$	7.5403
$c_{3,0}^g$	11.459	$c_{3,0}^{u_v}$	4.6595	$c_{3,0}^{dv}$	-0.8736	$c_{3,0}^{\overline{u}}$	8.2448
$c_{4,0}^g$	-	$c_{4,0}^{u_v}$	-	$c_{4,0}^{dv}$	-	$c_{4,0}^{\overline{u}}$	-
$c_{1,1}^g$	-50.064	$c_{1,1}^{u_v}$	-0.0616	$c_{1,1}^{dv}$	-52.218	$c_{1,1}^{\overline{u}}$	-7.4250
$c_{1,2}^g$	-0.0008	$c_{1,2}^{u_v}$	0.4455	$c_{1,2}^{dv}$	-0.1002	$c_{1,2}^{\overline{u}}$	-0.0021
$c_{2,1}^g$	-6.5209	$c_{2,1}^{u_v}$	-11.643	$c_{2,1}^{dv}$	3.1722	$c_{2,1}^{\overline{u}}$	-0.2658
$c_{2,2}^g$	0.2039	$c_{2,2}^{u_v}$	0.0002	$c_{2,2}^{dv}$	0.1336	$c_{2,2}^{\overline{u}}$	-0.2754

Table 4. Table 4: Same as table 3 , but at NNLO.

$g$	value	$u_v$	value	$d_v$	value	$\overline{u}$	value
$c_{0,0}^g$	6.2654	$c_{0,0}^{u_v}$	(SR)	$c_{0,0}^{dv}$	(SR)	$c_{0,0}^{\overline{u}}$	(SR)
$c_{1,0}^g$	0.2712	$c_{1,0}^{u_v}$	0.8060	$c_{1,0}^{dv}$	1.0227	$c_{1,0}^{\overline{u}}$	-0.1162
$c_{2,0}^g$	11.334	$c_{2,0}^{u_v}$	3.6897	$c_{2,0}^{dv}$	4.2717	$c_{2,0}^{\overline{u}}$	7.1632
$c_{3,0}^g$	5.0606	$c_{3,0}^{u_v}$	1.6388	$c_{3,0}^{dv}$	-0.6035	$c_{3,0}^{\overline{u}}$	-
$c_{4,0}^g$	-	$c_{4,0}^{u_v}$	-	$c_{4,0}^{dv}$	-	$c_{4,0}^{\overline{u}}$	2.4190
$c_{1,1}^g$	-2.4627	$c_{1,1}^{u_v}$	-0.1080	$c_{1,1}^{dv}$	-2.8603	$c_{1,1}^{\overline{u}}$	-3.2213
$c_{1,2}^g$	-0.0024	$c_{1,2}^{u_v}$	0.3766	$c_{1,2}^{dv}$	-0.0099	$c_{1,2}^{\overline{u}}$	-0.0123
$c_{2,1}^g$	-1.4764	$c_{2,1}^{u_v}$	74.620	$c_{2,1}^{dv}$	1.1235	$c_{2,1}^{\overline{u}}$	-0.0028
$c_{2,2}^g$	0.3704	$c_{2,2}^{u_v}$	-0.0001	$c_{2,2}^{dv}$	0.2357	$c_{2,2}^{\overline{u}}$	-0.9263

Equations

$$\frac{{\mathrm{d}}^2\sigma }{\mathrm{d}x\mathrm{d}y}=N^l\left[y^{2}x{F}_1^l+(1-y)F_2^l\mp \left(y-\frac{y^2}{2}\right)x{F}_3^l\right],$$

$$N^{l,\mathrm{NC}}=\frac{4\pi {\alpha }^2}{xy{Q}^2},$$

$$N^{\nu ,\mathrm{CC}}=\frac{G_{\mathrm{F}}^2{M}_W^4{Q}^2}{4\pi xy{\left(Q^2+M_W^2\right)}^2},$$

$$F_i\left(x,Q^2\right)=\sum _j{C}_i^j\left(x,{\alpha }_s({\mu }^2),{\mu }^2/Q^2\right)\otimes f_j\left(x,{\mu }^2\right),$$

$$F_1^{\nu }$$

$$F_2^{\nu }$$

$$F_3^{\nu }$$

$$F_1^{\bar{\nu }}$$

$$F_2^{\bar{\nu }}$$

$$F_3^{\bar{\nu }}$$

$$W^2\approx Q^2\left(\frac{1}{x}-1\right).$$

$$Q^2\approx yxs.$$

$$x{f}_i^{p/A}\left(x,Q_0^2\right)=c_0{x}^{c_1}(1-x{)}^{c_2}\left(1+c_{3}x+c_4{x}^2\right)$$

$$c_k\to c_k(A)=c_{k,0}+c_{k,1}\left(1-A^{-c_{k,2}}\right)$$

$$f_i^{N/A}\left(x,Q^2\right)=\frac{Z\cdot f_i^{p/A}+(A-Z)\cdot f_i^{n/A}}{A}.$$

$${\int }_0^1\mathrm{d}x{f}_{u_{\mathrm{v}}}^{p/A}\left(x,Q_0^2\right)=2,$$

$${\int }_0^1\mathrm{d}x{f}_{d_{\mathrm{v}}}^{p/A}\left(x,Q_0^2\right)=1,$$

$${\int }_0^1\mathrm{d}x\sum _{i}x{f}_i^{p/A}\left(x,Q_0^2\right)=1.$$

$${\chi }^2=\sum _i\frac{{\left({\mu }_i-{\hat{m}}_i\right)}^2}{{\Delta }_i^2}+\sum _{\alpha }{b}_{\alpha }^2$$

$${\hat{m}}_i=m_i+\sum _{\alpha }{\Gamma }_{i\alpha }{b}_{\alpha }.$$

$${\chi }^2\left(\text{m},\text{b}\right)=$$

$$+$$

$$H_{ij}=\frac{1}{2}\left(\frac{{\delta }^2{\chi }^2}{\delta y_i\delta y_j}\right),$$

$${\chi }^2={\chi }_0^2+\sum _{i,j}{H}_{ij}{y}_i{y}_j.$$

$$y_i=\sum _j{v}_{ij}\sqrt{\frac{1}{{ϵ}_j}}{z}_j,$$

$$\Delta {\chi }^2={\chi }^2-{\chi }_0^2=\sum _i{z}_i^2.$$

$$a_i=a_0\pm \Delta a_i=a_0\pm \Delta {\chi }^2\sum _j\frac{v_{ij}^2}{{ϵ}_j},$$

$$\displaystyle(\Delta X^{\pm})^{2}=\sum_{i=1}^{n_{\mathrm{param}}}\bigg{\{}{\begin{subarray}{c}\text{max}\\ \text{min}\end{subarray}}\bigg{[}$$

$$\displaystyle X(S_{i}^{-})-X(S_{0}),\,0\bigg{]}\bigg{\}}^{2}\,,$$

$${\hat{F}}_2^A\equiv \frac{1}{2}{F}_2^{p,A}+\frac{1}{2}{F}_2^{n,A}$$

$$F_2^A=\frac{Z}{A}{F}_2^{p,A}+\frac{N}{A}{F}_2^{n,A},$$

$${\hat{F}}_2^A=\beta F_2^A,$$

$$\beta =\frac{A}{2}\left(1+\frac{F_2^{n,A}}{F_2^{p,A}}\right)/\left(Z+N\frac{F_2^{n,A}}{F_2^{p,A}}\right).$$

Full text

1. I Introduction

2. II Theoretical framework

3. II.1 Deeply inelastic scattering

4. II.2 PDF parameterization

5. III Analysis procedure

6. III.1 Fitting procedure

7. III.2 Uncertainty analysis

8. III.3 The fitting framework

9. IV Experimental data

10. IV.1 Charged-lepton DIS data

11. IV.2 Neutrino DIS data

12. IV.3 Isoscalar modifications and experimental uncertainties

13. V Results

14. V.1 Proton baseline

15. V.2 Nuclear PDFs

16. V.3 Comparison to data

17. V.4 Comparisons to other nPDF sets

18. V.4.1 Comparison at NLO

19. V.4.2 Comparison at NNLO

20. VI Summary and outlook

21. A PDF parameters

I Introduction

Phenomenology based on collinear factorization Collins et al. (1989) has proven extremely successful in the era of the Large Hadron Collider (LHC). In this approach the parton distribution functions (PDFs) Kovařík et al. (2019), describing the number distributions of quarks and gluons in the colliding hadrons, are factorized from the hard partonic scattering. The latter can be calculated within perturbative quantum chromodynamics (QCD) Dokshitzer et al. (1991); Brock et al. (1995) but the PDFs have to be determined in a global analysis using experimental data and Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations that provide the scale evolution of the PDFs Lipatov (1975); Gribov and Lipatov (1972); Altarelli and Parisi (1977); Dokshitzer (1977). The most precise constraints for PDFs come from high-energy deeply inelastic scattering (DIS) experiments where the hadron structure is probed with a highly-virtual photon or a massive electroweak boson. The data available from the HERA collider Abramowicz et al. (2015), combined with older fixed-target measurements, provide plenty of data points with a broad kinematic reach that can be used to constrain the proton PDFs in the kinematic region relevant at the LHC Gao et al. (2018). However, the data suitable for nuclear PDF (nPDF) analyses are far more sparse Kovařík et al. (2019). In these analyses the nuclear modifications, first observed in DIS experiments with nuclear targets, are assumed to be nonperturbative and absorbed into the PDFs obeying the same scale evolution equations as free protons Armesto (2006).

The global nPDF analyses rely heavily on nuclear DIS data. Compared to the HERA data available for proton PDF fits, the kinematic reach of the fixed-target nuclear data is quite limited. Such data provide direct constraints for quarks, but the gluon distributions are probed only at higher orders in perturbative QCD (pQCD) and via scale-evolution effects. There is potential for improvement by including neutrino-nucleus DIS data that have additional sensitivity to the flavour decomposition of the PDFs due to the different coupling to up- and down-type quarks. Most of the recent next-to-leading order (NLO) nPDF analyses, e.g. DSSZ de Florian et al. (2012), include data for Drell-Yan (DY) dilepton production which provide additional constraints for antiquarks. Some further constraints for gluons have been obtained from pion-production data in d+Au collisions at the Relativistic Heavy-Ion Collider (RHIC) Adler et al. (2007) which were first used in the EPS08 analysis Eskola et al. (2008), and later also in EPS09 Eskola et al. (2009) and nCTEQ15 Kovarik et al. (2016). The EPPS16 Eskola et al. (2017) was the first analysis to use measurements from the LHC by incorporating data for $Z$ and $W^{\pm}$ bosons Khachatryan et al. (2015, 2016); Aad et al. (2015) and dijet production Chatrchyan et al. (2014) in p+Pb collisions. These provide further constraints for the flavour decomposition and gluon nuclear modifications, but the statistics of the Run I data for these observables was still quite limited. The more precise data from Run II from the LHC will bring more constraints especially for the gluon nPDFs Eskola et al. (2019a). Furthermore, the existing heavy-meson data from LHCb demonstrate promising potential to directly measure the gluon shadowing down to $x\sim 10^{-5}$ Zenaiev et al. (2015); Cacciari et al. (2015); Gauld and Rojo (2017); De Oliveira et al. (2018); Eskola et al. (2019b). A complementary channel to study small-$x$ gluons would be direct photon production at forward rapidities Helenius et al. (2014). Such a measurement with the required forward instrumentation could be performed with the proposed FoCal upgrade in ALICE Peitzmann (2018).

Another way to improve the PDFs is to increase the perturbative precision of the analysis. For proton PDFs the current standard is next-to-next-to-leading order (NNLO) Harland-Lang et al. (2015); Dulat et al. (2016); Abramowicz et al. (2015); Ball et al. (2017); Alekhin et al. (2017) in pQCD. The splitting functions required for the DGLAP evolution have been available at the required NNLO accuracy for a while already Vogt et al. (2004); Moch et al. (2004), but despite the recent progress Boughezal et al. (2017); Currie et al. (2017); Campbell et al. (2017); Gaunt et al. (2015); Czakon et al. (2016); Grazzini et al. (2018); Chen et al. (2019); Caola et al. (2019), the number of processes for which a full NNLO calculation is publicly available is still limited. For nPDFs first works in this direction have been performed recently. The first one, KA15 Khanpour and Atashbar Tehrani (2016), includes data from different fixed-target DIS experiments with a lepton beam. The fit is performed in the zero-mass variable-flavour-number-scheme (ZM-VFNS) where all quarks are assumed massless for energies above their mass threshold. The more recent analysis nNNPDF1.0 Abdul Khalek et al. (2019), applies the NNPDF methodology Forte et al. (2002); Ball et al. (2009) in which the resulting nPDFs are determined by a neural network. Also there, the applied data were restricted to neutral-current DIS, but a more realistic general-mass variable-flavour-number-scheme (GM-VFNS) was applied.

In this work we introduce two new nPDF analyses, performed at NLO and NNLO in pQCD, respectively, which we refer to as TUJU19. The presented work is based on the open-source xFitter package Zenaiev (2016); Bertone et al. (2018) (formerly known as HeraFitter Alekhin et al. (2015)) that has been modified in order to also accommodate data from nuclear collisions and suitable PDF parameterizations. In addition to neutral-current DIS with a lepton beam, we also include charged-current neutrino DIS ($\nu$DIS) data with nuclear targets that are sensitive to the flavour decomposition of nonisoscalar nuclei. For a free proton baseline we fit new PDF sets mainly based on the combined HERA I and II data, providing baseline fits consistent with our assumptions and kinematical cuts made for the nPDF analyses. Furthermore, the required extensions of the code will be published, providing a first open-source tool to analyze nuclear PDFs. In this paper we describe the theoretical framework in section II, then discuss the analysis procedure in section III and the selection of experimental data in section IV. The results of the analysis are presented in section V and the work is then summarized in section VI, where an outlook towards future developments is also presented.

II Theoretical framework

II.1 Deeply inelastic scattering

In this analysis, neutral-current (NC) and charged-current (CC) DIS processes are considered: NC in case of electron(positron)-nucleus ($\mathrm{e}A$) and CC for (anti)neutrino-nucleus ($\nu A$) scattering. For these processes the differential cross section is given by

[TABLE]

where $x$ and $y$ are the standard kinematic variables in DIS. The coupling factor $N^{\,l}$ depends on the scattering type and $F^{\,l}_{1,2,3}$ denote the structure functions Roberts (1994); Sterman at al. (2001); Patrignani et al. (2016) for the scattering of lepton $l$. In equation (1), the index $l$ covers different beams including $l=\nu,\,\bar{\nu},\,e^{+},\,e^{-},\,\mu^{+},\,\mu^{-}$. For the nuclear data used in this work two normalization factors are relevant: one for NC DIS in the case of unpolarized leptons and one for the CC DIS of incoming (anti)neutrinos. In the case of unpolarized leptons the normalization factor $N^{\,l,\,NC}$ for NC DIS is Patrignani et al. (2016)

[TABLE]

and in the case of incoming (anti)neutrinos the normalization factor $N^{\,\nu,\,CC}$ for CC DIS is

[TABLE]

where $Q^{2}$ is the virtuality of the intermediate boson that provides the scale at which the nucleons are probed. For the CC processes $G_{F}$ is the Fermi coupling constant and $M_{W}$ is the mass of the $W^{\pm}$ boson. When combining the structure functions into the differential cross sections in alignment with (1), the sign before $F_{3}$ is positive for $\nu$ and $e^{+}$, and negative for $\bar{\nu}$ and $e^{-}$.

In QCD the structure functions $F_{i}$, as introduced in eq. (1), are related to the scale-dependent parton distribution functions $f_{j}(x,\,Q^{2})$, with $j=g$ or $j=q\,,\bar{q}$, via

[TABLE]

where one typically chooses $\mu=Q$. The symbol $\otimes$ in equation (4) denotes a convolution between the parton distribution functions and the Wilson coefficients $C_{i}^{j}$ (see e.g. Refs. van Neerven and Vogt (2000a, b); Vermaseren et al. (2005) for $C_{2}^{j}$ at NLO and NNLO). For example, the structure functions for neutrinos and antineutrinos are given at leading order (LO) by Tanabashi et al. (2018)

[TABLE]

and

[TABLE]

As can be seen from relations (II.1) and (II.1), there is an added value to the analysis by including the neutrino-nucleus DIS data that have additional sensitivity to flavour decomposition of the PDFs due to the different coupling to up- and down-type quarks.

The factorization of the partonic scattering process and the nonperturbative PDFs Collins et al. (1989), as well as the perturbative treatment are valid at sufficiently high energy scales $(Q^{\,2}\gtrsim 1~{}\text{GeV}^{2})$. In this work we have selected kinematic cuts $Q^{\,2}>3.5\,\mathrm{GeV^{2}}$, the Bjorken variable $x<0.7$, and the invariant mass of the hadronic final state $W^{2}>12\,\mathrm{GeV^{2}}$. The latter can be expressed in terms of the other invariant variables as

[TABLE]

Some of the experimental data sets do not specify the invariant $y$, but when the collision energy is known, it can be derived from the relation

[TABLE]

II.2 PDF parameterization

A global DGLAP-based analysis requires a nonperturbative input for the PDFs at the initial scale of the fit. In this analysis parton distributions of a free proton and of a nucleon bound in a nucleus are parameterized as

[TABLE]

for parton flavour $i=g,\,d_{\mathrm{v}},\,u_{\mathrm{v}},\,\bar{u},\,\bar{d},\,\bar{s}$, at the initial scale $Q_{0}^{2}=1.69\,\mathrm{GeV^{2}}$. This form of PDF parameterization is similar to the functional form used in the HERAPDF2.0 analysis Abramowicz et al. (2015) and is motivated by the fact that the main constraints for the free proton PDF baseline come from the same DIS data. To keep the framework consistent we use the same form for the nuclear PDFs.

The main focus of this work is on the nuclear PDFs for which the fit parameters $c_{k}$ in equation (9) are re-parameterized to be dependent on the nuclear mass number $A$ as

[TABLE]

where $k={0,\dots,4}$. A similar form was also successfully used in the nCTEQ15 analysis Kovarik et al. (2016). At the same time, if $A=1$ the $A$-dependent right-hand part of equation (10) becomes zero and the free proton PDFs are recovered by default. The explicit $A$ dependence of the nuclear PDFs also allows us to make predictions for nuclei which were not part of the actual analysis, but are possibly interesting for future experiments.

As discussed, $xf^{p/A}_{i}$ given in equation (9) defines the parton distribution in a proton bound to a nucleus $A$. In addition there are also neutrons in a nucleus which we denote by $f_{i}^{\,n/A}$. The full PDF for a nucleon inside a nucleus can be obtained by averaging over the number of protons and neutrons in nuclei:

[TABLE]

The PDFs of neutrons are not separately fitted, but are determined from the proton PDFs based on isospin symmetry. In addition to this symmetry, we have assumed $s=\bar{s}$ and $s=\bar{s}=\bar{u}=\bar{d}$ as the included DIS data are not sensitive enough to constrain the strange-quark content, nor the sea quark flavour decomposition. In particular, even though the neutrino DIS data are sensitive to the separation of up- and down-type quarks, the kinematic region covered by the incorporated data ($x\gtrsim 0.01$)is where the valence quarks dominate the cross section.

The heavy quarks are treated within the general-mass variable-flavour-number-scheme (GM-VFNS), see Ref. Accardi et al. (2016) for a recent overview. There are several options for GM-VFNS implemented in xFitter111Also fixed-flavour (FF) mass schemes, like e.g. the ABM scheme Alekhin ; Alekhin and Moch (2011), are available in xFITTER., including (S)ACOT schemes Aivazis et al. (1994); Krämer et al. (2000); Tung et al. (2002); Kretzer et al. (2004), RT and RT optimal schemes Thorne and Roberts (1998); Thorne (2006, 2012), as well as a FONLL scheme Cacciari et al. (1998); Forte et al. (2010) for GM-VFNS. In this work we apply the FONLL-A scheme for the NLO analysis and the FONLL-C at NNLO, implemented in the APFEL package Bertone et al. (2014). The heavy-quark masses are fixed to $m_{\mathrm{charm}}=1.43~{}\mathrm{GeV}$ and $m_{\mathrm{bottom}}=4.50~{}\mathrm{GeV}$. The strong coupling constant $\alpha_{S}$ is set to $\alpha_{\mathrm{S}}(M_{\mathrm{Z}})=0.118$ for both the NLO and the NNLO fits.

For the parton distribution functions $xf^{p/A}_{i}$ as defined in equation (9), we assume the baryon number sum rules and the momentum sum rule satisfied by every nucleon in the nucleus,

[TABLE]

[TABLE]

[TABLE]

Strictly speaking, for nuclear parton distribution functions the sum rules are approximations that might not hold for individual nucleons in a nucleus in general, but which are reasonable at the available level of precision in regard to the experimental uncertainties. In this work, equation (12) is used to fix the normalization of $u_{\mathrm{v}}$ quarks in a proton and equation (13) defines the normalization of $d_{\mathrm{v}}$ quarks in a proton, nucleus per nucleus. The momentum sum rule (14) is used to constrain the normalization of the sea quarks. The remaining unconstrained normalization coefficient $c^{g}_{0}$ in the gluon PDF is treated as a regular free parameter during the fitting procedure. Alternatively, the gluon normalization could have been fixed by the momentum sum rule as done in many earlier analyses, e.g. Ref. Eskola et al. (2017).

As described above, we parameterize and fit the PDFs of a proton in a nucleus, and the neutron PDFs are determined based on SU(2) symmetry. In particular, the distributions of $u$ and $d$ quarks are exchanged: $u^{p}\rightarrow d^{n}$ and $d^{p}\rightarrow u^{n}$, which is valid for valence and sea quarks. For completeness we mention that this interchange requires the validity of charge symmetry, and in Ref. Boros et al. (1999) it has been suggested that some charge symmetry violation (CSV) could take place in the small-$x$ region. However, in the $x$ region covered by the nuclear DIS data we use, such effects should be negligible. Besides the DIS experiments, CSV effects can be studied in experiments measuring asymmetries in $W$ boson production. Further experiments and tests of CSV in parton distributions were suggested in Refs. Saito et al. (2000); Londergan et al. (2005). Similarly, possible isospin symmetry violations have been studied in Refs. Martin et al. (2004); Wang et al. (2016). In this work, however, we assume that the charge and isospin symmetries hold.

III Analysis procedure

III.1 Fitting procedure

The optimal values for the parameters are obtained by minimizing ${\large{\chi}}^{2}$ defined as

[TABLE]

with

[TABLE]

Here, $\mu_{i}$ is the value of the measured data point for a given observable, $\Delta_{i}$ is the uncorrelated experimental error, whereas the sum over correlated systematic errors is given by the term $\sum_{\alpha}b^{2}_{\alpha}$ in equation (15). The theoretical predictions for each data point $i$ are represented by $\hat{m}_{i}$, defined in equation (16). There, $m_{i}$ is the actual theoretical value calculated using DGLAP-evolved PDFs with given parameters $\{c_{k}\}$, $\Gamma_{i\alpha}$ are the correlated errors and $b_{\alpha}$ are the so-called nuisance parameters. A nuisance parameter quantifies the strength of the correlated error source $\alpha$, whereas $\Gamma_{i\alpha}$ quantifies the sensitivity of the $i^{\mathrm{th}}$ measurement to the correlated systematic error source $\alpha$. The quality of the fit can be estimated from the resulting ${\large{\chi}}^{2}/N_{\mathrm{dp}}$ ratio, where $N_{\mathrm{dp}}$ is the number of data points. A value ${\large{\chi}}^{2}/N_{\mathrm{dp}}\approx 1$ indicates that the agreement between the theoretical prediction and the measured observable is on average at the level of the experimental uncertainties.

There are several ways to take into account the correlated and uncorrelated uncertainties and to combine the statistical and uncorrelated systematic uncertainties in the ${\large{\chi}}^{2}$ definition in xFitter used in this analysis. Here we use the following definition

[TABLE]

The variables have been introduced in (15) with $\Gamma^{i}_{\alpha}=\gamma^{i}_{\alpha}\,\mu_{i}$, and $\delta_{i,\mathrm{stat}}$ and $\delta_{i,\mathrm{uncorr}}$ are the relative statistical and uncorrelated systematic uncertainties, respectively. The above form in (17) corresponds to a Poisson-like scaling for the statistical experimental uncertainties, whereas the systematic uncorrelated and correlated uncertainties are scaled linearly. This choice for ${\large{\chi}}^{2}$ is similar to the one used in the HERAPDF2.0 analysis incorporating the combined H1 and ZEUS DIS data Abramowicz et al. (2015), on which our proton PDF baseline is mainly based. For consistency the same form has also been used for the nuclear PDF analysis.

III.2 Uncertainty analysis

The minimization of ${\large{\chi}}^{2}$ provides a central set of PDFs with the parameter values providing the best description of the used data. However, the experimental data always contain several uncertainties of different types, such as statistical, systematic or correlated, like e.g. normalization errors. To study how well the experimental data actually constrain the fitted distributions, a separate error analysis needs to be performed. Such an analysis quantifies how much (according to the given criteria) room there is for the parameters to vary so that the resulting cross sections are still in agreement within the experimental uncertainties. The distribution functions resulting from the uncertainty analysis are typically provided as part of the PDF sets.

There are two established methods which can be used for the error analysis, the Hessian method Pumplin et al. (2001a, b) or the Monte Carlo (MC) method Giele and Keller (1998); Giele et al. (2001); Hartland et al. (2019). The former relies on quadratic approximation of ${\large{\chi}}^{2}$ with respect to the parameters $\{c_{k}\}$ near the minimum. In the latter method, the data is varied within the given uncertainties, and for each variation a PDF replica set is fitted. Therefore this method is less sensitive to the form of the parameterization but numerically more demanding. Also the Lagrange multiplier method Pumplin et al. (2001b); Stump et al. (2001) has been used to study the uncertainties (see e.g. Ref. Martin et al. (2003)), but there the error propagation to an observable becomes more involved. In this QCD analysis the Hessian method is used for the analysis of the uncertainties.

The Hessian error analysis is performed assuming a quadratic expansion of the function ${\large{\chi}}^{2}={\large{\chi}}_{0}^{2}+\Delta{\large{\chi}}^{2}$ around its global minimum. Here, ${\large{\chi}}_{0}^{2}$ is the value of the function at the global minimum (with the best-fit parameters $\{k_{0}\}$) and $\Delta{\large{\chi}}^{2}$ is the displacement from the minimum Pumplin et al. (2001a, b). The Hessian matrix $H$ is constructed from the second derivatives of ${\large{\chi}}^{2}$ at the minimum. The matrix elements $H_{ij}$ are defined as

[TABLE]

with $y_{i}$ being the displacement of the parameter $a_{i}$ from its value $a_{0}$ at the minimum. For the analysed function ${\large{\chi}}^{2}$ one writes

[TABLE]

The Hessian matrix is symmetric and thus has a complete set of orthonormal eigenvectors $v_{ij}$. The eigenvectors and the eigenvalues $\epsilon_{j}$ of the Hessian matrix are used to transform the displacements ${y_{i}}$ into a new set of parameters ${z_{i}}$

[TABLE]

leading to a simplified relation

[TABLE]

This representation has the advantage that the surfaces of constant ${\large{\chi}}^{2}$ are spheres in ${z_{i}}$ space with $\Delta{\large{\chi}}^{2}$ being the squared distance from the minimum. The varied parameters ${a_{i}}$ from which the resulting error sets are defined can then be written as

[TABLE]

where $\Delta{\large{\chi}}^{2}$ defines the tolerance criterion determining the allowed growth of ${\large{\chi}}^{2}$. The relation in equation (22) shows that the parameters which correspond to the eigenvectors of the Hessian matrix with large eigenvalues are well determined since their $\Delta a_{i}$ is small, whereas the weakly determined parameters correspond to small eigenvalues. The uncertainties for a given observable $X$ can be calculated via

[TABLE]

where $X(S_{0})$ is the observable calculated with the central parameter set and the $S_{i}^{\pm}$ correspond to the error sets in the positive and negative direction determined from the diagonalized parameter $z_{i}$.

In an ideal case one would choose the tolerance criterion so that $\Delta{\large{\chi}}^{2}=1$. However, since we consider several different data sets which are not necessarily in a mutual agreement with one another, such a choice would underestimate the underlying uncertainty. In this work the tolerance for $\Delta{\large{\chi}}^{2}$ is based on the statistically motivated method as discussed e.g. in Refs. Hirai et al. (2007); Eskola et al. (2017). For the proton baseline with $13$ free fit parameters it becomes $\Delta{\large{\chi}}^{2}=20$ at the 90% confidence level. This choice has also been validated by comparing to the error bands generated with the MC method, though the $\Delta{\large{\chi}}^{2}$ value preferred by the MC method is quite flavour- and kinematics-dependent. Previous nPDF studies have shown that such a statistically motivated method would not fully cover the experimental uncertainties in the nuclear data Kovarik et al. (2016); Eskola et al. (2017). Thus for the nuclear PDF error analysis we increase the tolerance from the statistically motivated value and choose $\Delta{\large{\chi}}^{2}=50$ for our 16 free parameters.

III.3 The fitting framework

The global analyses of the baseline proton and nuclear PDFs were performed with the xFitter Zenaiev (2016); Bertone et al. (2018) tool. The main goal of the xFitter project is to provide an open-source tool to fit proton PDFs with different theoretical assumptions. A schematic view of the fitting procedure and relations to different external programs are shown in figure 1. Being an open-source tool it is available to everyone and makes the research process fully transparent, which is important in order to establish a common knowledge base and a deep understanding of the opportunities and limitations. The released version covers various options like different PDF parameterization forms, mass schemes, etc. Furthermore, xFitter provides interfaces programmed in Fortran or C++ to the commonly used tools like MINUIT James and Roos (1975); Lazzaro and Moneta (2010), QCDNUM Botje (2011), APPLGrid Carli et al. (2010) or APFEL Bertone et al. (2014), etc. The DGLAP evolution routine and the calculation of DIS cross sections are implemented at NNLO. Further functionalities in regards to the future potential and alternative fitting approaches, including dipole models Luszczak and Kowalski (2014, 2017) and small-$x$ resummation Ball et al. (2018); Abdolmaleki et al. (2018), are available in the released version.

In order to perform a nuclear PDF analysis several modifications of the code were required. First, the PDF parameterization had to be adapted for the purpose of nuclear PDFs. Thus, new parameters $c_{k}(A)$ dependent on the nuclear mass number $A$ as per equation (10) were introduced. In order to reflect the new nPDF parameterization, the form of the steering file, as well as the file containing the initial parameters for MINUIT, and the according interpretation routine were adapted. As the next step, the mass number $A$ and proton number $Z$ of a given nucleus for the nucleon decomposition (cf. eq. (11)) of the up and down quarks were included. The possible combinations of data sets for different mass numbers $A$ and proton numbers $Z$ were kept flexible in order to deal with data for ratios between different nuclei. The information on $A$ and $Z$ depends on the data set and thus needs to be provided inside the data files. Therefore, the form and the routine to read the experimental data files were extended accordingly. Additionally, the overall minimization routine FCN was modified so that the DGLAP evolution can evolve nuclear PDFs covering different combinations of $A$ and $Z$ individually. Next, the calculation of sum rules (equations 12, 13, 14) had to be adapted in order to reflect the flexibility of an $A$-dependent normalization. Additionally, to keep the form of the PDF parameterization flexible, a new numerical integration routine for the calculation of sum rules was implemented.

Besides that, the cross section calculation routine was enhanced for the treatment of various isoscalar modifications, as described in subsection IV.3. Three flags identifying the ’NMC’, ’EMC’ or ’SLAC’ forms of the corrections were implemented. Furthermore, experimental nuclear data is often provided in terms of ratios $\sigma(A_{1})/\sigma(A_{2})$ or $F_{2}(A_{1})/F_{2}(A_{2})$ for two different nuclei $A_{1}$ and $A_{2}$. Thus, we had to extend the xFitter mechanisms for the consideration of these ratios by the implementation of a two-step loop. The underlying PDF flavour decomposition for a proton was modified so that in the case of a nucleus the PDF decomposition is applied for a nucleon of the form (11).

Moreover, charged-current (CC) processes for the neutrino DIS data222The implemented calculation of the reduced CC DIS cross section by using the FONLL scheme in xFITTER was customized to the HERA framework. Specifically, a factor $1/2$, e.g., in equation (8) of Ref. Adloff et al. (2003), is included in the reduced cross section to account for the factor 2 in the normalization of the structure functions used there. For the neutrino-nucleus DIS data, we use the nonreduced differential cross section with the prefactor from Tanabashi et al. (2018), and especially with the structure functions defined in equations (II.1) and (II.1). Thus, for the calculation of the CC DIS cross sections, a division by a factor of $2$ was removed from xFitter in the case of neutrino beams. were incorporated in xFitter according to the differential cross section described in subsection II.1. Finally, the uncertainty analysis routine333In this work the option ’DoBands’ has been used to generate asymmetric error bands, which is based on the ’iterate’ method by John Pumplin Pumplin et al. (2001b). Its advantage is that, if necessary, the iteration routine will add a positive value $X$ to all eigenvalues to force the matrix to be positive definite, which is as close as possible to the actual $\hat{H}$. The positive definiteness of the Hessian matrix relies on the second derivatives, which is a difficult computation and is numerically often approximate. One reason is that the minimized function ${\large{\chi}}^{2}$ is not exact, but given by a second-order polynomial in the space of the fit parameters. Thus, if some fit parameters are not well constrained by the data, higher-order polynomial terms of ${\large{\chi}}^{2}$ might become relevant Eskola et al. (2019a). Another point is that the function ${\large{\chi}}^{2}$ might not be as smooth as necessary due to the limited numerical precision at which the DGLAP equations are solved and due to the finite accuracy of the integrals. Pumplin et al. (2001a, b) was modified so that scaling of the error bands (cf. eq. 22) could also be performed for $\Delta{\large{\chi}}^{2}>1.0$. The modifications described in this section will be published as a part of the package later on.

IV Experimental data

IV.1 Charged-lepton DIS data

The QCD analyses presented here have been performed by including the experimental data from DIS measurements. The free proton baseline was fitted using data from the HERA Abramowicz et al. (2015), BCDMS Benvenuti et al. (1989) and NMC Arneodo et al. (1997) experiments, as listed in table 1. The fixed-target DIS data with lepton and neutrino beams used to determine the nuclear parton distribution functions are summarized in table 2. The kinematic reach of the included experimental DIS data in the $(x,Q^{2})$ plane is shown in figure 2 for the different nuclear targets. The applied kinematic cuts $x\leq 0.7$ and $Q^{2}\geq 3.5\,\mathrm{GeV^{2}}$ are illustrated by the dotted lines in the plots. The number of available data points varies for different nuclei. A large number of data points are available for deuteron (D) and the heavier nuclei of carbon (C), iron (Fe), tin (Sn) and lead (Pb), as shown in the upper panel of figure 2. These data points are provided either in the form of absolute cross sections, or as ratios where D is usually used as the reference (denominator). Also calcium (Ca) has been intensively used in the relevant experiments. For the other nuclei (lower panel in figure 2), only a few data points are available. Therefore the nuclei predominantly present in the included data are expected to be better constrained than the nuclei with fewer data points.

IV.2 Neutrino DIS data

Neutrino data were included in the analysis by using the measured cross sections for neutrino and antineutrino beams. The advantage compared to the isospin-averaged structure functions $F_{2}$ and $F_{3}$ utilized in DSSZ de Florian et al. (2012) is that the sensitivity to the flavour decomposition is retained in the cross sections. Another approach was used in EPPS16 Eskola et al. (2017) where normalized (anti)neutrino cross sections were considered. This increases the sensitivity to the shape of the nuclear modifications. In order to extract complete, i.e. without isospin averaging, information from the incorporated neutrino data sets, the absolute cross sections are exploited here.

The data from the CDHSW $\nu\mathrm{Fe}$ experiment Berge et al. (1991) and the CHORUS $\nu\mathrm{Pb}$ experiment Onengut et al. (2006) have been included in this analysis. In addition there is more neutrino scattering data available, e.g. measured cross sections with an $\mathrm{Fe}$ target by the the NuTeV collaboration Tzanov et al. (2006), and also data from the CCFRR collaboration MacFarlane et al. (1984). The data from the CCFRR experiment were excluded from our analysis for two reasons. First, the quantities $x$ and $Q^{2}$, required for the analysis procedure, were not publicly available for the cross sections. Second, only the averaged structure functions $F_{2}$ and $F_{3}$ for $\nu\mathrm{Fe}$ and $\bar{\nu}\mathrm{Fe}$ were available, which lose the sensitivity to flavour decomposition. In regard to NuTeV data, an early study in Ref. Owens et al. (2007) found that these data could be accommodated together with the CHORUS neutrino DIS data when constraining the $d/u$ ratio but with the applied nuclear corrections some tension with other DIS and DY data were observed. Later on, the analyses in Refs. Schienbein et al. (2008, 2009) found some unresolved tension between the NuTeV neutrino DIS data and lepton-nucleus data. In a following work, a similar tension was also found when taking into account neutrino DIS data from the CHORUS and CCFRR experiments in Ref. Kovarik et al. (2011). Simultaneously, a study presented in Ref. Paukkunen and Salgado (2010) concluded that the tension with other data was specifically due to the NuTeV data at certain energies, whereas CDHSW and CHORUS data were well compatible with the existing lepton-nucleus DIS data. This tension was further studied in Ref. Paukkunen and Salgado (2013) where again the NuTeV data were found to be incompatible with the other considered data. Only by normalizing the differential data with the integrated cross section at each energy bin was an acceptable agreement achieved. Due to the demonstrated tension, we have not included the NuTeV neutrino DIS data in this analysis.

IV.3 Isoscalar modifications and experimental uncertainties

Some experimental analyses of charged-lepton DIS have modified the measured structure functions to achieve isospin symmetry for nonsymmetric nuclei such as iron or lead. According to the relations summarized in Ref. Eskola et al. (2017), an isoscalar structure function of a nucleus with the mass number $A$ is defined as

[TABLE]

with $F_{2}^{p,A}$ and $F_{2}^{n,A}$ representing the structure functions of the bound protons and neutrons. By definition, the isoscalar structure function contains an equal number of protons and neutrons, which holds only for specific nuclei. A general structure function for a nucleus with $Z$ protons and $N=A-Z$ neutrons can be written as

[TABLE]

which is not isoscalar if $Z\neq N\neq A/2$. The relation between the isoscalar structure function (24) and the general structure function (25) is given by

[TABLE]

where

[TABLE]

Furthermore, it is assumed that the ratio $F_{2}^{n,A}/F_{2}^{p,A}$ for any nucleus is unchanged compared to that for unbound nucleons, so that the relation

[TABLE]

can be used in (27) to describe the isoscalar modifications. The ratio $F_{2}^{n}/F_{2}^{p}$ for the isoscalar “correction” is parameterized in a different way by each experiment,

• •

EMC Ashman et al. (1993):

[TABLE]

• •

SLAC Gomez et al. (1994):

[TABLE]

• •

NMC Amaudruz et al. (1992):

[TABLE]

[TABLE]

In this work the general form of the structure function (25) is used to calculate the theoretical predictions. In case isoscalar modifications were applied to the measured quantities, for consistency the same modifications are applied to the obtained theoretical results by using relation (26).

Some of the experiments provide normalization uncertainties on top of the systematic and statistical errors. In this work normalization uncertainties have been treated as correlated errors as discussed in Refs. Stump et al. (2001); Pumplin et al. (2001a, b). The correlated uncertainties are treated as implemented in xFitter, described further in Refs. Stump et al. (2001); Botje (2002); Aaron et al. (2009). The same procedure applies if any overall uncertainties are provided in addition to the point-to-point uncertainties, e.g. for the SLAC data Gomez et al. (1994); Bodek et al. (1983); Dasu et al. (1994).

V Results

V.1 Proton baseline

Analyses of nuclear PDFs have often been performed by using an existing proton PDF set as a baseline for the nuclear modifications. In this work, however, we have fitted the proton PDFs using the same setup as for the nuclear PDFs. This ensures that all assumptions like sum rules, parton flavour decomposition, etc., as well as all parameters like coupling constants and quark masses, and also further settings like e.g. the heavy flavour mass scheme, are applied in a consistent way. Furthermore, this paves the way for a future combined proton and nuclear PDF analysis.

As we use xFitter as our analysis framework, the baseline proton PDFs are derived with a very similar setup as for the HERA2.0 PDFs Abramowicz et al. (2015). However, in addition to the combined HERA DIS data we also include data from other experiments (cf. table 1). Another difference is that we use the parameterization in eq. (9) whereas the HERA2.0 analysis includes additional terms for the gluon at the initial scale of the analysis. The obtained parton distribution functions are compared to the HERA2.0 PDFs Abramowicz et al. (2015) in figure 3 at NLO and in figure 4 at NNLO. As shown in Ref. Abramowicz et al. (2015), the HERA2.0 PDFs are well compatible with other state-of-the-art proton PDFs in the kinematic region considered in this work, and since the main focus of this work is on nuclear PDFs we do not present further comparisons to other proton PDF sets. As expected, because of our use of the same fitting framework with similar data and similar definition of ${\large{\chi}}^{2}$, the agreement with the HERA2.0 PDFs is very good both at NLO and NNLO. The observed difference of the gluon PDFs at small $x$ can be traced back to the different parameterization applied. Since we also include data from experiments other than HERA, we have used a larger $\Delta{\large{\chi}}^{2}$ value, namely $\Delta{\large{\chi}}^{2}=20$ (see section III.2), which results in larger uncertainties than those quoted for the HERA2.0 PDFs.

V.2 Nuclear PDFs

The resulting nuclear PDFs, referred to as TUJU19, are presented in figure 5 at NLO and in figure 6 at NNLO for a few different nuclei, together with the fitted proton baseline PDFs. As the sea-quark nPDFs have been assumed to be flavour independent, i.e. $s=\bar{s}=\bar{u}=\bar{d}$, the $x\bar{u}\left(x,\,Q^{2}\right)$ distribution represents all sea quarks in the plots. Many earlier analyses have assumed that the nuclear modifications for the deuteron are negligible and constructed its PDFs from the free proton ones using isospin symmetry. In this work we, instead, treat the deuteron as a nucleus in the fitting procedure. Small deviations from the proton PDFs are found for the proton in a deuteron, as shown in figures 5 and 6. The deviation from the proton PDFs becomes larger with increasing $A$, and significant effects are found in heavy nuclei such as iron and lead. The optimal parameters according to the chosen parameterizations in eqs. (9) and (10) are listed in appendix A.

The nuclear modifications of the PDFs for the lead nucleus, defined as

[TABLE]

where $f_{i}^{\mathrm{p/Pb}}(x,\,Q^{2})$ and $f_{i}^{\mathrm{p}}(x,\,Q^{2})$ are the PDFs for the bound and the free proton, respectively, are shown in figure 7. The NLO and NNLO modifications are compared at the initial scale of the analysis ($Q^{2}=1.69~{}\text{GeV}^{2}$) and at a higher scale ($Q^{2}=100~{}\text{GeV}^{2}$) after DGLAP evolution. In both cases, the ratio of gluon PDFs shows some low-$x$ shadowing and a rapid rise with increasing $x$. This behaviour is similar to what was observed in the HKN07 analysis Hirai et al. (2007), but in our case the enhancements are moderated at higher scales and a recognizable antishadowing peak develops around $x\sim 0.3$. For the sea quarks the typical nuclear modifications – shadowing, antishadowing and EMC suppression – are visible already at the initial scale. However, especially the magnitude of the small-$x$ shadowing differs at different orders, NNLO favoring a stronger effect. Since for gluons the behaviour is opposite we conclude that these differences arise due to the fact that the sea-quark and gluon evolution are coupled and the applied DIS data is not sensitive enough to fully separate the contributions. At higher scales the sea quark modifications are in better agreement though some difference still persists at large $x$. The valence-quark parameters were allowed to be flavour-dependent also for nPDFs and the resulting nuclear effects indeed become rather different for bound $u_{v}$ and $d_{v}$ distributions. For $d_{v}$ again the typical features including shadowing and antishadowing are well visible but for $u_{v}$ we find that some amount of low-$x$ enhancement is preferred. We note that a similar behaviour was observed in the nCTEQ15 analysis Kovarik et al. (2016), although no neutrino DIS data were included there that would provide additional flavour sensitivity especially for the valence sector. One should keep in mind that the full nPDF for an average nucleon will be the sum of those for protons and neutrons, so the opposite behaviour will cancel out to a certain extent.

The uncertainty bands for the nPDFs provided in this work have been generated with $\Delta{\large{\chi}}^{2}=50$ as described in subsection III.2. The resulting uncertainty bands do, however, depend also on the flexibility of the applied parameterization. Due to the limited sensitivity of the applied data to the gluon and sea quark nPDFs, we had to limit the number of $A$-dependent parameters in order to achieve numerical convergence of the fits. Therefore the provided uncertainty bands for the gluon distribution likely underestimate the true uncertainty to some extent, which should be kept in mind when comparing to previous works. In the future, by adding more data providing further constraints one could consider admitting more parameters and therefore allowing for a larger flexibility of the parameterization.

V.3 Comparison to data

The optimal set of nPDF parameters is derived by minimizing the ${\large{\chi}}^{2}$ as defined in eq. (15) by comparing to the measured data presented in table 2. The resulting cross sections, structure functions and ratios are compared to the data we used in figures 8, 10 - 15 for neutral current DIS processes and in figures 16 and 17 for charged current DIS processes with neutrinos for a subset of that data. An overview of the resulting ${\large{\chi}}^{2}$ values, divided by the number of data points $N_{\mathrm{dp}}$, is shown in figure 9 for NLO and NNLO. Values above ${\large{\chi}}^{2}/N_{\mathrm{dp}}>3.0$ have been truncated in this graph for better representation, but the actual numbers are shown in table 2. Figure 9 demonstrates that the agreement between the theoretical predictions and the experimental measurements varies between different data sets. For example, the agreement with most of the data published by the NMC collaboration is excellent, whereas the agreement with the HERMES data is clearly not optimal. In this particular example one needs to point out that the number of data points from the HERMES experiment is much smaller than the number of NMC data points, so that the contribution to the total ${\large{\chi}}^{2}$ is relatively small for the HERMES data. Apart from the few outlying data sets, the overall agreement is found to be very good, and the total ${\large{\chi}}^{2}/N_{\mathrm{dp}}$ is $0.887$ at NLO and $0.862$ at NNLO. Even though some of the data sets are better described at NLO and some at NNLO, the total ${\large{\chi}}^{2}$ values are very close at the different orders. The good agreement is also apparent in figures 8–17. Interestingly, a very good agreement is also achieved for the neutrino data, even though some earlier studies observed difficulties when incorporating these data in a global nPDF analysis Kovarik et al. (2011). However, as concluded in Ref. Paukkunen and Salgado (2013), this likely because of the tension caused by the NuTeV data which we have not included.

V.4 Comparisons to other nPDF sets

Each nPDF analysis is based on a set of assumptions, e.g. the form of the nonperturbative input at the initial scale, the choice of the proton baseline and the kinematical cuts. Therefore, even when based on the same set of data it is not guaranteed that the results will be equivalent. However, some level of agreement – within the estimated uncertainties – is expected.

V.4.1 Comparison at NLO

In figure 18 we compare our obtained nPDFs to those of other recent NLO nPDF analyses, nCTEQ15 and EPPS16, at our initial scale $Q^{2}=1.69~{}\text{GeV}^{2}$ and at $Q^{2}=100~{}\text{GeV}^{2}$. The comparisons are shown for $g$, $\bar{u}$, $u_{v}$, $d_{v}$ and $V=u_{v}+d_{v}$ in a proton bound in lead. For gluons at the initial scale the agreement is not very good, though still well within the uncertainties. Towards higher scales, however, a much better mutual agreement is observed. For sea quarks (here represented by $\bar{u}$) the agreement with the previous analyses is better already at the initial scale, and at $Q^{2}=100~{}\text{GeV}^{2}$ our result is between EPPS16 and nCTEQ15. When comparing the sea-quark distributions one should keep in mind that each analysis has different assumptions for the sea-quark flavour dependence, i.e. in TUJU19 we have assumed $s=\bar{s}=\bar{u}=\bar{d}$, whereas for nCTEQ15 $s=\bar{s}$ and $\bar{u}=\bar{d}$ are connected by an additional factor, and only $s=\bar{s}$ applies for EPPS16. For valence quarks we find that $u_{v}$ tends to stay below (above) the EPPS16 and nCTEQ15 results at $x\gtrsim 0.03$ ($x\lesssim 0.03$) whereas the opposite behaviour is found for $d_{v}$. This can be explained by the fact that in the case of nuclear data only a combination of $u_{v}$ and $d_{v}$ is probed, and even with the included neutrino data the flavour dependence of valence quarks is not well constrained. Indeed, we find a very good agreement between the three analyses for the sum of valence quarks $V$.

The uncertainty bands in our NLO fit are similar to those obtained in the earlier analyses for sea quarks, but for gluons the resulting uncertainties are somewhat smaller. Since both EPPS16 and nCTEQ15 include additional data with some sensitivity to gluons, we conclude that our reduced uncertainties are likely due to the limited number of free parameters in the gluon nPDFs, and the uncertainty due to the lack of data constraints is underestimated. One should note that nCTEQ15 does not provide error sets for the baseline proton PDFs, which partly explains why their uncertainties for the sea quarks (at the initial scale) and the valence quarks tend to be smaller than in EPPS16 and this analysis. The comparisons were generated by using the LHAPDF6 library Buckley et al. (2015) and the published grids.

In figure 19 we compare the nuclear modification of the PDFs as defined in eq. (32) at $Q^{2}=100~{}\text{GeV}^{2}$. Also comparisons to the DSSZ analysis are included, for which only ratios $R_{i}^{\mathrm{p/Pb}}$ (eq. 32) were available with error bands444No LHAPDF6 grids are available for DSSZ. . In most cases the results are compatible within the estimated uncertainties, though some features stand out. A rather prominent feature of our NLO gluons is the large antishadowing around $x\sim 0.3$. Such a large enhancement is not supported by other analyses which include data sensitive to gluon antishadowing and underlines the need for further data sensitive to such effects. However, the obtained gluon shadowing is in good agreement with EPPS16 and nCTEQ15 results, though with somewhat reduced uncertainty estimates. Only the DSSZ $R_{g}^{\mathrm{p/Pb}}$ with very mild shadowing is outside the uncertainty bands in this region. For the flavour dependence of the valence quarks we find a similar behaviour as nCTEQ15 where some small-$x$ enhancement and large-$x$ suppression were observed for $u_{v}$, along with the opposite behaviour for $d_{v}$. However, when calculating the total valence distribution for a complete nucleon, as shown in figure 18, we find a good agreement with the other analyses.

V.4.2 Comparison at NNLO

The comparison of the TUJU19 NNLO nPDF fit to other NNLO nPDF analyses is shown in figure 20 for nNNPDF1.0 Abdul Khalek et al. (2019) and in figure 21 for KA15 Khanpour and Atashbar Tehrani (2016). The comparison at NNLO is separated into two figures since different information is available for the nNNPDF1.0555The nNNPDF1.0 LHAPDF6 set is provided with the assumption that $u=d$ and that $\bar{u}=\bar{d}=\bar{s}=s$ to comply with the LHAPDF format, i.e. to provide individual quark flavours. For our comparison, we have used the provided PDFs and LHAPDF uncertainties at the $90\%$ confidence level, keeping in mind that only the quark multiplets $\Sigma$, $T_{3}$, $T_{8}$ were determined in the nNNPDF1.0 analysis. and KA15666For the comparison at NNLO the KA15 nPDFs, provided on request by the authors, were also included. The ratios $R^{\mathrm{p/Pb}}$ were not available for this analysis. analyses. In the case of nNNPDF1.0 we consider lead nPDFs for gluon $g$ and the quark singlet $\Sigma=u+\bar{u}+d+\bar{d}+s+\bar{s}$ (as per Abdul Khalek et al. (2019)). The comparisons are shown at the two scales $Q_{0}^{2}=1.69~{}\text{GeV}^{2}$ and $Q^{2}=100~{}\text{GeV}^{2}$ for the distribution functions, and at the higher scale $Q^{2}=100\,\mathrm{GeV}^{2}$ for the ratios $R_{i}^{\mathrm{p/Pb}}$ of PDFs in a proton bound in lead compared to PDFs in a free proton. Even though none of the analyses include data directly sensitive to the gluon distribution, a reasonable behaviour is found for the gluon at the initial scale of this analysis, keeping in mind that our uncertainty bands for the gluon are potentially underestimated, as discussed above. The congruence is improved when going towards higher scales (e.g. $Q^{2}=100~{}\text{GeV}^{2}$). Furthermore, a very good agreement is observed for the quark singlet $\Sigma$ shown in the lower panels of figure 20. The values of $x\Sigma$ are within the error bands at the initial scale, and become even more consistent at the higher scale. The observed deviation in the low-$x$ region ($x<0.0005$) reflects the lack of low-$x$ constraints by the available nuclear DIS data. We also show the ratios $R^{\mathrm{p/Pb}}$ for TUJU19 compared to nNNPDF1.0 at NNLO in figure 20. Again a reasonable shape is found for the gluon nuclear modification, and a very good agreement is visible for the quark singlet, as it is well constrained by the incorporated experimental data. The uncertainties of the nNNPDF1.0 distributions are considerably larger compared to the ones found in our analysis. Indeed, as shown in Ref. Abdul Khalek et al. (2019), the uncertainty bands for a combination of singlet and octet contributions become comparable to those in the earlier works where direct data constraints exist. Still, the nNNPDF1.0 uncertainties grow rapidly towards small-$x$ which can be accounted for by the fact that the applied neural network framework is not as prone to parameterization bias as the traditional Hessian error analysis applied here.

The nuclear parton distribution functions TUJU19 at NNLO compared to the results of the KA15 group are presented in figure 21. Here we consider lead nPDFs for gluons, sea quarks (here $\bar{s}$) and the sum of valence quarks, $V=u_{v}+d_{v}$. Again the comparisons are shown at two scales, at the initial scale of KA15, $Q_{0}^{2}=2.0\,\mathrm{GeV}^{2}$, and at $Q^{2}=100~{}\text{GeV}^{2}$. Considering the fact that neither of the analyses includes data directly sensitive to the gluon distribution, a fair agreement is found for the gluon at the initial scales. However, while the agreement between TUJU19 and nNNPDF1.0 remains at higher scales (fig. 20), the gluon distribution from KA15 falls below the other two (fig. 21) at $Q^{2}=100~{}\text{GeV}^{2}$. The $\bar{s}$ distributions in turn are in a reasonable agreement at higher scales, although at the initial scale the KA15 result is considerably above TUJU19. The total valence distributions from TUJU19 are found to be in very good agreement with those of the KA15 analysis. Apart from the gluon nPDFs at the initial scale, the KA15 uncertainties tend to be very small. This may partly follow from the rather rigid parameterization applied, but it may also be due to the chosen low error tolerance $\Delta{\large{\chi}}^{2}=1$.

VI Summary and outlook

We have presented new sets of nPDFs, which we refer to as TUJU19, at NLO and NNLO. Contrary to most previous analyses, our nPDF sets are based on a proton baseline fitted within the same framework, which guarantees consistency throughout the analysis concerning the series of choices on parameter values, assumptions, constraints and kinematical cuts that need to be made when performing a global analysis. The numerical implementation was embedded in the open-source tool xFitter. The source code with all applied modifications required for the treatment of nuclear PDFs will be published, providing a first open-source tool for a global nuclear PDF analysis. The common framework will also enable a simultaneous proton and nuclear PDF analysis in the future.

For the proton PDF baseline analysis DIS data from the HERA, BCDMS and NMC experiments were included. Similarly, the nPDF analysis was based on measurements from fixed-target neutral-current DIS and data from neutrino-nucleus charged-current DIS experiments. The neutrino DIS data were implemented for the first time in a global pQCD analysis at NNLO. The deuteron was treated as a nucleus without neglecting nuclear effects, instead of constructing it as a pure composition of free proton and neutron PDFs as been assumed in several earlier analyses. The resulting cross sections show very good agreement with the included experimental data, both for neutral-current and charged-current DIS processes, as confirmed by the resulting ${\large{\chi}}^{2}/N_{\mathrm{dp}}<1.0$ for the nuclear part of the analysis. The comparisons to the existing nPDF sets demonstrate a reasonable agreement within the error bands. The obtained results are consistent with the expectation that due to the consideration of DIS data only the nPDFs for valence quarks are well constrained by the experimental data, whereas gluon and sea quarks are constrained only indirectly by the included data and mostly by DGLAP evolution. The resulting nPDFs will be published in the LHAPDF6 format including uncertainties for both, the proton baseline and the nuclear PDF analysis, derived with a Hessian uncertainty analysis.

As this is the first NNLO nPDF fit within the developed framework, only fixed-target DIS data with lepton and neutrino beams were included. In the future we plan to also add data for other observables for which theoretical calculations at NNLO exist. The fixed-target proton-nucleus DY dilepton data would provide further constraints for the sea quark distributions. Furthermore, the $W^{pm}$ and $Z$ boson production data from p+Pb collisions at the LHC are sensitive to the flavour decomposition and could therefore help to disentangle observed differences in valence quark nuclear modifications. Even after these, direct gluon constraints will remain sparse. Recently it has been shown that such constraints could be obtained from the existing data for dijet and charmed-meson production in p+Pb collisions at the LHC. Further in the future, an electron-ion collider would provide precision data for nPDF analyses. In order to get the best information from these data the highest possible perturbative precision will be required, and we think our NNLO analysis is an important step in this context.

Acknowledgments

The authors thank the xFitter development team for providing operative support. In particular Dr. Alexander Glazov, Dr. Oleksandr Zenaiev are acknowledged for sharing their expertise on various aspects of the tool and on the details of uncertainty analysis, kindly supported by Prof. Fred Olness. Furthermore, many thanks to Dr. Valerio Bertone for clarifying details on the implementation of the FONLL scheme and to Dr. Francesco Giuli for sharing the reference input files. The authors acknowledge discussions and scientific exchange with the other nPDF groups. We thank Dr. Aleksander Kusina and Prof. Fred Olness for discussing details in regards to the nCTEQ15 publication, and Dr. Karol Kovarik for sharing his expertise on the nCTEQ framework. Furthermore, the authors acknowledge discussions, technical input, and comments on the manuscript by Dr. Hannu Paukkunen. The authors are grateful to Dr. Hamzeh Khanpour for sharing the KA15 PDFs. We thank Dr. Pia Zurita for sending the DSSZ PDFs. The authors acknowledge Dr. Marco Stratmann for valuable discussions and information. The authors acknowledge support by the state of Baden-Württemberg through bwHPC providing the possibility to run the computational calculations on the high-performance cluster. This work was supported in part by the Bundesministerium für Bildung und Forschung (BMBF) grant 05P18VTCA1, Carl Zeiss Foundation and by the Academy of Finland, project n.o. 308301.

Appendix A PDF parameters

Here we show the parameters obtained for the proton and nuclear parton distribution functions presented in section V. The naming convention corresponds to the PDF parameterization given in equations (9) and (10). Table 3 provides the NLO parameters, while table 4 presents the NNLO ones.

Some of the parameters were deliberately excluded from the fit. In most cases this means that in the initial version of the analysis procedure those parameters were used, but the obtained parameter value turned out to be very close to zero with very large uncertainty. Thus, that parameter was considered as not required. Alternatively, some of the nuclear parameters were never included as free parameters since the best fit criterion for nuclear PDFs, ${\large{\chi}}^{2}\leq 1.0$, could be satisfied by the selected subset of the free parameters.