††institutetext: 1INPAC, Shanghai Key Laboratory for Particle Physics and
Cosmology, School of Physics and Astronomy, Shanghai Jiao-Tong University,
Shanghai 200240, China††institutetext: 2Key Laboratory for Particle
Astrophysics and Cosmology (MOE), Shanghai 200240, China††institutetext:
3Center for High Energy Physics, Peking University, Beijing 100871, China

#  Machine learning of log-likelihood functions in global analysis of parton
distributions

DianYu Liu1, ChuanLe Sun1, Jun Gao1,2,3<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Modern analysis on parton distribution functions (PDFs) requires calculations
of the log-likelihood functions from thousands of experimental data points,
and scans of multi-dimensional parameter space with tens of degrees of
freedom. In conventional analysis the Hessian approximation has been widely
used for the estimation of the PDF uncertainties. The Lagrange Multiplier (LM)
scan while being a more faithful method is less used due to computational
limitations, and is the main focus of this study. We propose to use Neural
Networks (NNs) and machine learning techniques to model the profile of the
log-likelihood functions or cross sections for multi-dimensional parameter
space in order to overcome those limitations which work beyond the quadratic
approximations and meanwhile ensures efficient scans of the full parameter
space. We demonstrate the efficiency of the new approach in the framework of
the CT18 global analysis of PDFs by constructing NNs for various target
functions, and performing LM scans on PDFs and cross sections at hadron
colliders. We further study the impact of the NOMAD dimuon data on
constraining PDFs with the new approach, and find enhanced strange-quark
distributions and reduced PDF uncertainties. Moreover, we show how the
approach can be used to constrain new physics beyond the Standard Model (BSM)
by a joint fit of both PDFs and Wilson coefficients of operators in the SM
effective field theory.

###### Keywords:

PDFs, QCD, Machine Learning

##  1 Introduction

Precise understanding of the parton structure of the proton is a central topic
of QCD Gao:2017yyd ; Kovarik:2019xvh . The parton structure can be described
by parton distribution functions (PDFs), which represent distributions of
momentum fractions of the proton carried by quarks and gluons, for instance in
the case of QCD collinear factorization Collins:1989gx . They are usually
determined by fitting to a variety of experimental data, such as data from
proton-proton collision, proton-antiproton collision, electron-proton
collision, and neutrino–nucleus scatterings. Besides, there have also been
recent developments on calculating PDFs from first principles based on the
large momentum effective theory Ji:2020ect  and lattice QCD simulations
LatticeParton:2018gjr .

Especially, PDFs play important roles in LHC studies. For example, PDF
uncertainties represent one of the dominant uncertainties in measurements of
the Higgs boson couplings LHCHiggsCrossSectionWorkingGroup:2016ypw . Better
control of PDF uncertainties are necessary in direct searches for new heavy
resonances Beenakker:2015rna  and indirect searches for new physics beyond the
SM Alioli:2017jdo . Furthermore, PDF uncertainties also have a large impact on
precision measurements of the SM parameters including the strong coupling
constant ATLAS:2020mee , the weak mixing angle and the W boson mass
Kaur:2019ndj ; ATLAS:2017rzl .

Modern analysis of PDFs requires calculations of the log-likelihood functions
from thousands of experimental data points, and scans of multi-dimensional
parameter space with tens of degrees of freedom. There are several groups
providing regular updates of PDFs via global fits, see Refs. Hou:2019efy ;
Bailey:2020ooq ; Ball:2021leu ; Alekhin:2018pai ; H1:2017bml ; Jimenez-
Delgado:2014twa ; Park:2021kgf ; ATLAS:2021vod  for recent results on PDF
determinations. The difference between those PDF sets is mainly due to the
choice of the experimental data sets, the theoretical calculations used, and
the parametrization form of PDFs.

PDF uncertainties can be determined with three methods: the Hessian
Pumplin:2001ct ; Martin:2009iq , Monte Carlo (MC) Forte:2002fg , and Lagrange
Multiplier (LM) Pumplin:2000vx ; Stump:2001gu  method. There also exist
recently developed approaches, meta analysis Gao:2013bia , ePump
Schmidt:2018hvu  and L2 sensitivity Wang:2018heo , on accessing impacts of
experimental data on PDFs based on the Hessian method. In the Hessian method,
the log-likelihood function ($\chi^{2}$) of a global fit is approximated with
a quadratic form of the PDF parameters at the neighborhood of the global
minimum. The uncertainties are thus determined through error PDFs along
eigenvector directions, constructed by requiring the increase of the total
$\chi^{2}$ of 1 or of a certain tolerance. In the MC method, one can obtain
the PDF uncertainties from an ensemble of PDF replicas which are fitted to an
ensemble of “pseudo-data”. Those pseudo-data are generated from the
probability distributions related to the original experimental data sets. On
another hand, for the LM method, PDF uncertainties of an observable can be
determined from the profiled $\chi^{2}$ as a function of the observable,
without relying on any assumptions about the behavior of the $\chi^{2}$ at the
neighborhood of the global minimum. This means PDF uncertainties estimated
from the LM method are more robust than those from the Hessian method.
However, the LM method requires a detailed scan of the PDF parameter space for
every observable studied, which is usually time consuming.

This drawback can be overcome with the help of machine learning (ML). ML has
been widely used in studies of high-energy physics in recent years. In many
cases, ML is used for classifications such as particle identification and
event selection in experimental data analysis Guest:2018yhq . Neural networks
(NNs) are also helpful in regression problems, for example, applications of
NNs in the study of PDFs have been pioneered by the NNPDF collaboration
NNPDF:2014otw . Dependence of PDFs on the momentum fraction are parametrized
using NNs, which ensures a great flexibility Forte:2020yip . On another hand,
dependence of the $\chi^{2}$ or any physics quantity, such as the cross
section, on PDFs is complex in general. NNs offer an opportunity to relate
physics quantities to PDFs efficiently. One can build NNs with PDFs as input
variables to model their PDF dependence. Compared with traditional methods,
NNs can greatly improve efficiencies on generating predictions for those
physics quantities.

With above motivations, in this paper we propose a new approach with which PDF
uncertainties can be calculated efficiently using the LM method with the
assistance of NNs. It takes three steps to achieve this goal. First, we
construct and train NNs to model the $\chi^{2}$ of each individual data set
used in the global fit with PDFs. Second, we construct and train other NNs to
associate the physics quantity to be studied with PDFs. Finally, we can
perform LM scans to determine PDF uncertainties in a robust way. The speed of
LM scans can be improved by several orders of magnitude due to the
introduction of the NNs. We demonstrate above idea in the framework of CT18
NNLO global analysis Hou:2019efy  and beyond. We show how the new approach can
help to understand various PDF uncertainties and the interplay between
different data sets in the global fit. Moreover, we explore several directions
beyond CT18 as will be explained below.

Only a few data sets in the CT18 global fit are sensitive to the strange-quark
PDFs. The dimuon production in neutrino scatterings provides an opportunity to
directly constrain strange-quark distributions in the nucleon. In recent NOMAD
measurements NOMAD:2013hbk , a sample of about $9\times 10^{6}$ events of
inclusive charged-current deep-inelastic scattering (CCDIS), together with
about 15344 events of dimuon production, is collected. The large statistics
lead to a better control on various systematic errors and also an improvement
in statistical uncertainties. We include the NOMAD data in the global fit and
evaluate the impact on the PDFs using the aforementioned approach.

The High Luminosity LHC (HL-LHC) is supposed to accumulate an integrated
luminosity of 3000 fb-1 for ATLAS and CMS and of 300 fb-1 for LHCb
Dainese:2019rgk . We take two of those HL-LHC pseudo-data sets constructed in
Ref. AbdulKhalek:2018rok ; AbdulKhalek:2019mps , the high-mass Drell-Yan data
and the forward W/Z production data, and evaluate their impacts on PDFs. Our
projection shows they can largely improve separations of different flavors,
especially for sea quarks.

In the searches for new physics beyond the SM from scatterings involving
nucleons, for instance at HERA or LHC, one key problem is on the degeneracy of
PDF variations and the new physics contributions, especially in cases when
similar measurements are used in both the global fit of PDFs and in the
searches of new physics. Ideally a joint global fit including both PDFs and
model parameters of the new physics should be performed, see Refs.
Carrazza:2019sec ; Greljo:2021kvv ; Madigan:2021uho ; CMS:2021yzl ;
Iranipour:2022iak  for examples. We demonstrate successful application of our
approach in such scenario by a simultaneous fit of both PDFs and the Wilson
coefficient of lepton-quark contact interactions in the SM effective field
theory (SMEFT).

The rest of this paper is organized as follows. In Section 2, we describe the
basic setup of our approach, including architectures of the NNs, PDF
parametrizations and experimental data sets considered in the global fit. In
Section 3, we discuss performances of the approach and show that the accuracy
of approximations with NNs are far sufficient for phenomenological studies. In
Section 4, we explain the method of LM scans and discuss several features of
the CT18 analysis based on the new approach. In Section 5, we study the impact
of the NOMAD measurements and of the two pseudo-data of HL-LHC on PDFs, and
show a joint fit with both PDFs and new physics contributions. Finally, we
conclude in Section 6.

##  2 Setup of the Neural Network program

In this section, we give a brief introduction to the setup of our NNs,
including the architectures, the input variables, and the target functions. We
further explain the training processes from the generation of samples to the
minimization of the loss function.

###  2.1 Basic setup of NNs

The general structure of NNs includes three parts: the input layer, several
hidden layers and the output layer. Each of these layers contains a collection
of nodes termed by perceptrons. There exist various implementations of NNs,
and we use Keras keras  in this work. From the NNs built by Keras, PDFs as
inputs are associated with either $\chi^{2}$ or physics quantities as outputs.
The log-likelihood function $\chi^{2}$ quantifies agreements between theory
predictions and experimental measurements for each data set and is calculated
according to Hou:2019efy . The physics quantities considered include cross
sections of several benchmark processes at the LHC, and PDFs or their ratios
at different $Q$ values. An example of the architecture of our NNs is shown in
Fig. 1, in which the inputs are PDFs at an initial scale and the outputs are
the $\chi^{2}$ of the fit to an experimental data set. In this figure, the
PDFs $f_{i}(x,Q)$ are evaluated at an initial scale of $Q=1.295$ GeV with $x$
selected among 14 different values, and $i\in\\{g,u,d,\bar{u},\bar{d},s\\}$
runs over all parton flavors. We always assume $s=\bar{s}$ at the initial
scale. They altogether form the input layer with 84 nodes {$I_{1}$, $I_{2}$
$\ldots$ $I_{84}$}. In addition, the differences in setups between NNs for
different target functions are shown in Table 1. The choice on the
architecture is based on the observation that cross sections or evolved PDFs
are in general non-linear functions of the PDF parameters. The $\chi^{2}$ is
positive defined and is a sum of various individual terms that depend on cross
sections quadratically, and thus can be approximated by a more complicated
architecture as prescribed. We include more details on the construction of our
NNs in Appendix A.

Figure 1: An example of the architecture of NNs in this work, taking
$\chi^{2}$ as the target function.

Target |  No. of hidden layers |  No. of nodes for each hidden layer |  Activation functions for each layer |  No. of total params
---|---|---|---|---
$\chi^{2}$ | 2 | 60,40 |  $\tanh$, $(x^{2}+2)$, linear | 7581
$\sigma$, $f_{i}(x,Q)$, $f_{i}(x,Q)$/$f_{j}(x,Q)$ | 1 | 40 |  $\tanh$, linear | 3441

Table 1:  The architecture of NNs in this paper. Structure is set up for
either $\chi^{2}$ or other quantities.

To construct the $l_{th}$ layer of a NN, we define

$b_{i}^{(l)}=\left\\{\begin{array}[]{ll}\sum_{j}w_{ij}^{(l)}I_{j},&(l=1),\\\
\\\ \sum_{j}w_{ij}^{(l)}h_{j}^{(l-1)},&(l>1),\end{array}\right.$ (1)

where $b_{i}^{(l)}$ is the value before the activation of the $i_{th}$ node in
the $l_{th}$ layer, $w_{ij}^{(l)}$ is the weight matrix connecting the
$(l-1)_{th}$ layer to the $l_{th}$ layer, $I_{j}$ is the value of the $j_{th}$
node in the input layer, and $h_{j}^{(l-1)}$ is the value of the $j_{th}$ node
in the $(l-1)_{th}$ layer. The value of the $i_{th}$ node in the $l_{th}$
layer is then obtained by applying the activation function $t^{(l)}$ on
$b_{i}^{(l)}$:

$h_{i}^{(l)}=t^{(l)}(b_{i}^{(l)}).$ (2)

This procedure iterates over all hidden layers, and in the end we obtain a
single value for the output layer. The activation functions used include the
conventional choices of linear, and $\tanh$ types, as well as a customized one
of quadratic form, depending on the target functions and layers. Note we
constrain elements of the weight matrix of the output layer to be positive for
the NN associated with the $\chi^{2}$ since it is positive definite. Elements
in the weight matrix, $w^{(l)}_{ij}$, are trained to minimize the so-called
loss function, which is defined as

$d_{loss}=\frac{1}{n}\sum_{k=1}^{n}\left(A^{k}_{{\rm NN}}(w_{ij})-A^{k}_{{\rm
TR}}\right)^{2},$ (3)

where $n$ is the total number of events in the training sample and $A^{k}_{\rm
TR}$ and $A^{k}_{\rm NN}$ are the truth of the target function and the
prediction from NNs for the $k_{th}$ event.

###  2.2 PDF parametrization form

The parametrization form of PDFs used at the initial scale $Q_{0}$ is

$f_{i}\left(x,Q_{0}\right)=a_{0}x^{a_{1}-1}(1-x)^{a_{2}}P_{i}\left(y;a_{3},a_{4},\ldots\right),$
(4)

where $\\{a_{1},a_{2},\ldots\\}$ are free parameters, and the behavior of
$x^{a_{1}}$ at $x\to 0$ and $(1-x)^{a_{2}}$ at $x\to 1$ is guided by Regge
theory and spectator counting rules respectively.
$P_{i}\left(y;a_{3},a_{4},\ldots\right)$ is a polynomial dependent on
$y\equiv\sqrt{x}~{}(y\equiv 1-(1-\sqrt{x})^{a_{3}})$ for valence quark and
gluon PDFs (light-quark sea PDFs). Parametrization forms used here are the
same as in the CT18 NNLO analysis Hou:2019efy .

For the valence-quark ($u_{v}$ and $d_{v}$) PDF,

$\displaystyle f_{v}\left(x,Q_{0}\right)=$ $\displaystyle
a_{0}x^{a_{1}-1}(1-x)^{a_{2}}P_{v}(y),$ (5) $\displaystyle P_{v}(y)=$
$\displaystyle\sinh\left[a_{3}\right](1-y)^{4}+\sinh\left[a_{4}\right]4y(1-y)^{3}+\sinh\left[a_{5}\right]6y^{2}(1-y)^{2}$
$\displaystyle+(1+\frac{1}{2}a_{1})4y^{3}(1-y)+y^{4}.$

For the gluon PDF,

$\displaystyle f_{g}\left(x,Q_{0}\right)=$ $\displaystyle
a_{0}x^{a_{1}-1}(1-x)^{a_{2}}P_{g}(y),$ (6) $\displaystyle P_{g}(y)=$
$\displaystyle\sinh\left[a_{3}\right](1-y)^{3}+\sinh\left[a_{4}\right]3y(1-y)^{2}+(3+2a_{1})y^{2}(1-y)+y^{3}.$

For the sea quark ($\bar{u}$, $\bar{d}$ and $s\equiv\bar{s}$) PDF,

$\displaystyle f_{\bar{q}}\left(x,Q_{0}\right)=$ $\displaystyle
a_{0}x^{a_{1}-1}(1-x)^{a_{2}}P_{\bar{q}}(y),$ (7) $\displaystyle
P_{\bar{q}}(y)=$
$\displaystyle(1-y)^{5}+a_{4}5y(1-y)^{4}+a_{5}10y^{2}(1-y)^{3}+a_{6}10y^{3}(1-y)^{2}$
$\displaystyle+a_{7}5y^{4}(1-y)+a_{8}y^{5}.$

In all, we have 8 free parameters for valence quarks after applying the
valence sum rules and letting $a_{1}$ be equal for $u_{v}$ and $d_{v}$. We
have 15 free parameters for sea quarks after fixing some of those $a_{i}$ or
letting them be equal for different flavors Hou:2019efy . We are left with 5
free parameters for gluon after applying the momentum sum rule. The total
number of free PDF parameters is 28.

###  2.3 Targets and samples

In this paper, we associated PDFs with $\chi^{2}$ and other physics quantities
through our NNs. Details of these target functions are described in the
following:

$\bullet$ The individual $\chi^{2}$ of each data set in an NNLO global
analysis of PDFs. We use the same 39 experimental data sets as in CT18 NNLO
global analysis. These experimental data sets are summarized in Table 2. The
theoretical calculations used are explained in the CT18 paper Hou:2019efy . We
take those calculations from CT18 except for minor updates on NNLO K-factors
of several data sets. The global $\chi^{2}$ is simply a sum of the 39
individual $\chi^{2}$.

$\bullet$ The cross sections of Higgs boson pair (top-quark pair with a Higgs
boson) production in proton-proton collisions at center of mass energy
$\sqrt{s}=$ 13 TeV or 100 TeV. They are computed at leading (next-to-leading)
order in QCD using MG5_aMC@NLO Alwall:2014hca  and AMCfast Bertone:2014zva  to
provide an interface with APPLgrid Carli:2010rw . We choose these two
processes for demonstrations, and any scattering cross sections at hadron
collisions can be included in a similar way.

$\bullet$ The PDFs and PDF ratios at various $x$ and $Q$ values. They are
obtained using HOPPET Salam:2008qg  with DGLAP evolutions at NNLO.

ID | Experimental data set | $N_{pt}$ | ID | Experimental data set | $N_{pt}$
---|---|---|---|---|---
160 |  HERA I+II 1 $fb^{-1}$, H1 and ZEUS NC and CC reduced cross sec. comb. H1:2015ubc  | 1120 | 101 | BCDMS $F_{2}^{p}$ BCDMS:1989qop  | 337
102 | BCDMS $F_{2}^{d}$ BCDMS:1989ggw  | 250 | 104 | NMC $F_{2}^{d}/F_{2}^{p}$ NewMuon:1996fwh  | 123
108 | CDHSW $F_{2}^{p}$ Berge:1989hr  | 85 | 109 | CDHSW $x_{B}F_{3}^{p}$ Berge:1989hr  | 96
110 | CCFR $F_{2}^{p}$ CCFRNuTeV:2000qwc  | 69 | 111 | CCFR $x_{B}F_{3}^{p}$ Seligman:1997mc  | 86
124 | NuTeV $\nu\mu\mu$ SIDIS osti_879078  | 38 | 125 | NuTeV $\bar{\nu}\mu\mu$ SIDIS osti_879078  | 33
126 | CCFR $\nu\mu\mu$ SIDIS NuTeV:2001dfo  | 40 | 127 | CCFR $\bar{\nu}\mu\mu$ SIDIS NuTeV:2001dfo  | 38
145 | H1 $\sigma_{r}^{b}$ H1:2004esl  | 10 | 147 | Combined HERA charm production H1:2012xnw  | 47
169 | H1 $F_{L}$ H1:2010fzx  | 9 | 201 | E605 Drell-Yan process Moreno:1990sf  | 119
203 | E866 Drell-Yan process $\sigma_{pd}/(2\sigma_{pp})$ NuSea:2001idv  | 15 | 204 | E866 Drell-Yan process $Q^{3}d^{2}\sigma_{pp}/(dQdx_{F})$ NuSea:2003qoe  | 184
225 | CDF Run-1 lepton $A_{ch}$, $p_{Tl}>$ 25 GeV CDF:1998uzn  | 11 | 227 | CDF Run-2 electron $A_{ch}$, $p_{Tl}>$ 25 GeV CDF:2005cgc  | 11
234 | D$\varnothing$ Run-2 muon $A_{ch}$, $p_{Tl}>$ 20 GeV D0:2007pcy  | 9 | 260 | D$\varnothing$ Z rapidity D0:2007djv  | 28
261 | CDF Run-2 Z rapidity CDF:2010vek  | 29 | 266 | CMS 7 TeV 4.7 fb-1, moun $A_{ch}$, $p_{Tl}>$ 35 GeV CMS:2013pzl  | 11
267 | CMS 7 TeV 840 fb-1, electron $A_{ch}$, $p_{Tl}>$ 35 GeV CMS:2012ivw  | 11 | 268 | ATLAS 7 TeV 35 $pb^{-1}$ W/Z cross section, $A_{ch}$ ATLAS:2011qdp  | 41
281 | D$\varnothing$ Run-2 9.7 $fb^{-1}$ electron $A_{ch}$, $p_{Tl}>$ 25 GeV D0:2014kma  | 13 | 504 | CDF Run-2 inclusive jet production CDF:2008hmn  | 72
514 | D$\varnothing$ Run-2 inclusive jet production D0:2008nou  | 110 | 245 |  LHCb 7 TeV 1.0 fb-1 W/Z forward rapidity cross sec. LHCb:2015okr  | 33
246 |  LHCb 8 TeV 2.0 fb-1 $Z\to e^{-}e^{+}$ forward rapidity cross sec. LHCb:2015kwa  | 17 | 249 |  CMS 8 TeV 18.8 fb-1 muon charge asymmetry $A_{ch}$ CMS:2016qqr  | 11
250 | LHCb 8 TeV 2.0 fb-1 W/Z cross sec. LHCb:2015mad  | 34 | 253 | ATLAS 8 TeV 20.3 fb-1, Z $p_{T}$ cross sec. ATLAS:2015iiu  | 27
542 |  CMS 7 TeV 5 fb-1, single incl. jet cross sec., R = 0.7 (extended in y) CMS:2014nvq  | 158 | 544 |  ATLAS 7 TeV 4.5 fb-1, single incl. jet cross sec., R = 0.6 ATLAS:2014riz  | 140
545 |  CMS 8 TeV 19.7 fb-1, single incl. jet cross sec., R = 0.7, (extended in y) CMS:2016lna  | 185 | 573 |  CMS 8 TeV 19.7 fb-1, $t\bar{t}$ norm. double-diff. top $p_{T}$ and y cross sec. CMS:2017iqf  | 16
580 |  ATLAS 8 TeV 20.3 fb-1, $t\bar{t}$ $p_{T}^{t}$ and $m_{t\bar{t}}$ abs. spectrum ATLAS:2015lsn  | 15 |  |  |

Table 2: Experimental data sets involved in the global fit Hou:2019efy .

We first generate randomly a training sample consisting of 6000 replicas of
PDFs and another test sample of 2000 replicas to prevent from over training.
Details about the generation of the replicas of PDFs can be found in Appendix
A. We compute all the target functions ($\chi^{2}$ or physics quantities) for
each of the replicas, which can be time consuming depending on whether the
fast interpolation approaches, like APPLgrid or FastNLO, are used or not.
However, we only need to perform these heavy calculations once for all.
Afterwards we construct a NN for each of the target function considered with
the architectures shown in Table 1. We train each NN for about 10 hours,
depending slightly on the architecture, on a single CPU-core (2.4 GHz)
according to the loss function defined in Eq. (3). Thus for all $\chi^{2}$ of
the 39 individual data sets that takes about 390 core-hours in total for the
training process. We found a very good performance of the resulting NNs
without much tuning on the training process for all target functions studied,
which will be reported in the next section. In a later stage for the
evaluation of the target functions with arbitrary PDF parameters, we can
simply use the optimized NNs rather than direct calculations. Comparison
between computational cost of the NNs and the direct computations are
summarized in Appendix A where substantial improvements in the speed from the
NNs are observed.

##  3 Validation of NNs

In this section, we perform several comparisons between the truths and the
predictions from our NNs before we apply them to further phenomenological
studies. We emphasize that the entire NN approach we discussed so far and in
the following is bound to the CT18 parametrization form, especially with the
CT18 PDF set. All PDF replicas for training and testing are sampled from the
CT18 PDFs. A first attempt of generalization to other parametrization forms or
even independent of PDF parametrization shows promising results, and is
detailed in Appendix A. It should be noted that the NNs should be retrained in
general if the underlying PDF parametrization changes.

###  3.1 $\chi^{2}$ of the global fit

In Fig. 2, we show the predictions to truths ratios of $\chi^{2}$ for three
experimental data sets: measurements of the proton structure function by
BCDMS, measurements of inclusive DIS reduced cross sections at HERA and
measurements of the inclusive jet cross sections at $\sqrt{s}=$ 7 TeV by CMS.
The ratio of total $\chi^{2}$ for the full data set is also shown in the
lower-right panel. The ratios are calculated for the PDFs from the
aforementioned training sample and test sample of NNs as well as the CT18 NNLO
PDFs. The CT18 NNLO PDFs consist of a central PDF set and 56 error PDFs in a
total of 28 Hessian eigenvector directions. The horizontal axis represents the
truths of $\chi^{2}$. Each mark corresponds to a PDF set from these three
samples of PDF sets. The green squares and the blue circles represent the
ratios corresponding to the PDFs from the training sample and test sample
respectively, and the purple triangles represent the ratios corresponding to
the PDFs from CT18 NNLO. We find good agreement between the training and test
samples, although the NN produces greater deviation than the original CT18
NNLO PDFs. We find that the predictions and the truths in general agree within
1 per mille for each data set. For the total $\chi^{2}$, the deviation is
within 0.6 per mille.

((a))

((b))

((c))

((d))

Figure 2: The predictions to truths ratios of $\chi^{2}$ for experimental data
sets for measurements of the proton structure function by BCDMS, measurements
of inclusive DIS reduced cross sections at HERA and measurements of the
inclusive jet cross sections at $\sqrt{s}=$ 7 TeV by CMS as well as for the
full data set.

We define $\Delta\chi^{2}$ as the difference between a certain $\chi^{2}$
value and its value at the best fit, which is conventionally used in the
determination of PDF uncertainties. In Fig. 3, differences between predicted
and true $\Delta\chi^{2}$ denoted as
$\delta(\Delta\chi^{2}_{\alpha})\equiv\Delta\chi^{2}_{\alpha,pre}-\Delta\chi^{2}_{\alpha,tru}$
are demonstrated for each data set, where $\alpha$ represents the PDF set used
in the comparison. Here we choose a sample of PDF sets consisting of the 56
Hessian error PDFs in CT18 NNLO set, which are represented by the marks
distributed along the vertical direction. We also show similar results for the
total $\chi^{2}$. It can be seen that, for each individual data set the
$\delta(\Delta\chi^{2}_{\alpha})$ is at most 1 unit, and the
$\delta(\Delta\chi^{2}_{\alpha})$ for the full data set are within 2 units.
The extent of $\delta(\Delta\chi^{2}_{\alpha})$ for HERA inclusive DIS data
set and for total $\chi^{2}$ is slightly larger than other experimental data
sets. It should be noticed that the number of data points of HERA inclusive
DIS data set and the full data set is 1120 and 3671 respectively. Besides, the
$\Delta\chi^{2}_{\alpha}$ of the full data set for most of the 56 Hessian
error PDFs is about 100 units. This indicates the relative deviation of NNs
predictions is below 2% for $\Delta\chi^{2}$.

Figure 3:  The $\delta(\Delta\chi^{2}_{\alpha})$ corresponding to the 56 error
PDFs of CT18 NNLO are represented by the 56 marks distributed along the
vertical direction for each individual data set and for the full data set.

Furthermore, we compare the $\Delta\chi^{2}$ for the full data set along the
28 eigenvector directions of CT18 NNLO PDFs by scans of PDF parameters. A
variable $d$ is introduced to measure the distance that PDF parameters go
along the direction of a certain eigenvector. The variation of PDF parameters
for the scan along the $j_{th}$ eigenvector direction can be written as:

$a_{i}^{j,scan}(d)=\left\\{\begin{array}[]{ll}d(a_{i}^{2j-1}-a_{i}^{0})+a_{i}^{0},&(d>0),\\\
\\\ d(a_{i}^{0}-a_{i}^{2j})+a_{i}^{0},&(d<0),\end{array}\right.$ (8)

where $i$ represents the index of the PDF parameters, {$a_{i}^{0}$} represents
PDF parameters for the central PDF of CT18 NNLO, {$a_{i}^{2j-1}$} and
{$a_{i}^{2j}$} represent PDF parameters of the two error PDFs in the $j_{th}$
eigenvector direction of CT18 NNLO. The total $\chi^{2}$ are computed using
the new set of PDF parameters. We define
$\Delta{\chi^{2}}\equiv\chi^{2}(d)-\chi^{2}(d=0)$, and compare the truths of
$\Delta\chi^{2}$ and the predictions from NNs as a function of $d$ for a few
selected eigenvector directions in Fig. 4.

((a))

((b))

((c))

((d))

((e))

((f))

Figure 4: The variation of $\Delta\chi^{2}$ with $d$ along the 1st, 6th, 11th,
16th, 22th and 28th CT18NNLO Hessian eigenvector directions.

We find that the NNs can describe well the dependence of the $\Delta\chi^{2}$
on PDF parameters in all Hessian eigenvector directions. The predictions and
the truths in general agree within 2% in all directions, which agrees with
Fig. 3. We also notice that $\Delta\chi^{2}$ has a sizable deviation from
quadratic shape in some directions. The NNs can reproduce well the asymmetric
and non-quadratic behavior of $\Delta\chi^{2}$, which is one of the main
advantages as comparing to the traditional Hessian method. It is justified to
say the deviation of $\chi^{2}$ due to the NNs approximation is negligible for
the PDF parameter space of interest.

###  3.2 Physics quantities

In Fig. 5, we show the ratios of the predictions to the truths for the cross
section of top-quark pair with a Higgs boson production in proton-proton
collisions at a center of mass energy $\sqrt{s}$ = 13 TeV, and for the PDF
ratio $d/u$($x=0.3$, $Q=100$ GeV). The ratios are calculated for the PDFs from
the training sample and test sample of NNs as well as the CT18 NNLO PDFs. We
find good agreements between the distributions of marks for training sample
and for test sample. Deviations for these two physics quantities are in
general within 0.15 per mille and 0.2 per mille, respectively. The performance
of the NNs for these two physics quantities is better than that for
$\chi^{2}$, which is because the dependence of $\chi^{2}$ on PDFs is more
complex than the cases for cross sections or PDF ratios. The dependence of
these physics quantities on PDFs is even close to linear.

((a))

((b))

Figure 5: The ratios of the predictions from NNs to the truths for the cross
section of top-quark pair with a Higgs boson production in proton-proton
collisions at a center of mass energy $\sqrt{s}$ = 13 TeV and the PDF ratio
$d/u$($x=0.3$, $Q=100$ GeV).

We further summarize the relative difference between the predictions from NNs
and the truths for various physics quantities in Fig. 6. For each physics
quantity, 57 marks distributed along the vertical direction correspond to the
results from the 57 CT18 NNLO PDFs. The results for the cross sections of
Higgs boson pair production and top-quark pair with a Higgs boson production
in proton-proton collisions at center of mass energy $\sqrt{s}=$ 13 TeV or 100
TeV are shown in this figure. We find the predictions from NNs and the truths
for these cross sections agree within 0.2 per mille. In addition, the results
for cross sections of $pp\to hh$ and $pp\to ht\bar{t}$ with high invariant
mass $m_{hh}>$ 2.5 TeV or $m_{ht\bar{t}}>$ 2.5 TeV are also shown in Fig. 6,
and the relative difference between the predictions and truths in general
agree within 0.3 per mille. Comparisons for strangeness ratio
$R_{s}\equiv\dfrac{s(x,Q)+\bar{s}(x,Q)}{\bar{u}(x,Q)+\bar{d}(x,Q)}$ and PDF
ratios $d/u$ and $\bar{d}/\bar{u}$ at various $x$ and $Q$ are also shown in
this figure, and the predictions and the truths agree within 1 per mille. We
also show the results for PDF values for $g$ and $s$-quark at various $x$ and
$Q$, and the deviations between the predictions and the truths are within 0.75
per mille.

Figure 6: The relative difference between the predictions from NNs and the
truths for various physics quantities. For each physics quantity, 57 marks
distributed along the vertical direction correspond to the results from 57
PDFs of CT18 NNLO.

##  4 Lagrange Multiplier scans

LM scan is a robust method to estimate PDF uncertainties, which was originally
developed in Refs. Pumplin:2000vx ; Stump:2001gu . In this method a physics
quantity $X(\\{a_{i}\\})$ is introduced to the global fit as a Lagrange
multiplier. Then the new function that needs to be minimized in the global fit
becomes

$\Psi\left(\lambda,\left\\{a_{i}\right\\}\right)\equiv\chi^{2}\left(\left\\{a_{i}\right\\}\right)+\lambda{X}\left(\left\\{a_{i}\right\\}\right),$
(9)

where $\lambda$ is a specified constant. For each value of $\lambda$, one can
determine a set of {$a_{i}$}, $X(\\{a_{i}\\})$ and $\chi^{2}$ by minimizing
$\Psi$. Here the $\chi^{2}$ corresponds to the minimum of a constrained fit
with $X\\{a_{i}\\}$ fixed to the corresponding value. Specially, the central
value of $X(\\{a_{i}\\})$ and the global minimum of $\chi^{2}$,
$\chi^{2}_{min}$, can be determined by setting $\lambda=0$. A parametrically
defined curve ($X$, $\chi^{2}$) can be determined by repeating the
minimization for many values of $\lambda$. This means the $\chi^{2}$ of the
global fit depends on the value of $X(\\{a_{i}\\})$ and can be represented as
$\chi^{2}=\chi^{2}_{min}+\Delta\chi^{2}$. The PDF uncertainty of
$X(\\{a_{i}\\})$ can be determined by requiring $\Delta\chi^{2}+P=T$, here $T$
is the so-called “tolerance factor”. We assume that 90% CL region corresponds
to $T=100$. The penalty term $P$, called Tier-2 penalty, is introduced to
ensure the tolerance will be reached as soon as any data set shows
disagreement at 90% CL. The detailed definition of the penalty term can be
found in Refs. Gao:2017yyd ; Dulat:2015mca .

In comparison, we briefly describe the calculation of PDF uncertainties in the
framework of the Hessian method. Given the physics quantity $X(\\{a_{i}\\})$,
the asymmetric PDF uncertainties can be calculated as Lai:2010vv

$\begin{array}[]{l}\delta^{+}X=\sqrt{\sum_{i=1}^{N_{d}}\left[\max\left(X_{2i-1}-X_{0},X_{2i}-X_{0},0\right)\right]^{2}},\\\
\delta^{-}X=\sqrt{\sum_{i=1}^{N_{d}}\left[\max\left(X_{0}-X_{2i-1},X_{0}-X_{2i},0\right)\right]^{2}},\end{array}$
(10)

where $X_{0}$ represents the value of the physics quantity with the central
PDF of the Hessian set, $X_{2i-1}$ ($X_{2i}$) represents the value of the
physics quantity with the error PDF of the Hessian set in the positive
(negative) direction of the $i_{th}$ eigenvector in the $N_{d}$-dimensional
PDF parameter space.

###  4.1 LM scans on PDFs

We first study PDF values and ratios with LM scans based on the aforementioned
NNs approximation of $\chi^{2}$ and physics quantities. The results are shown
in Fig. 7. The black and the red solid lines represent $\Delta\chi^{2}$ and
$\Delta\chi^{2}+P$ respectively. The dot and the dash lines indicate the
contributions to $\Delta\chi^{2}$ from individual data sets. The blue and the
green vertical dot-dash lines indicate the uncertainties at 90% CL determined
with the LM method by requiring $\Delta\chi^{2}+P$ = 100 and with the Hessian
method from the published CT18 NNLO PDFs, respectively. Among the generic
features of the scans, it can be seen that the profile of the total
$\Delta\chi^{2}$ and individual $\Delta\chi^{2}$ show almost a quadratic
dependence on the variable at the neighborhood of the global minimum, which is
a requirement of the Hessian method. Some individual data sets prefer PDF
values or ratios that differ significantly from those at the global minimum.
Besides, the HERA inclusive DIS data play important roles in all cases, which
can be understood as due to the high experimental precision and the large
number of data points. The penalty term also gives strong constraints on some
PDF values or ratios. In addition, there are some slight differences between
the uncertainties determined with the Hessian method and the LM scans, which
is expected.

((a))

((b))

((c))

((d))

Figure 7: LM scans on the $g$ ($x=0.01$ GeV, $Q=125$ GeV), $s$ ($x=0.3$,
$Q=100$ GeV), $d/u$ ($x=0.3$, $Q=100$ GeV) and $R_{s}$ ($x=0.023$, $Q=1.5$
GeV). The black and the red solid lines represent $\Delta\chi^{2}$ and
$\Delta\chi^{2}+P$ respectively. The dot and the dash lines indicate the
contributions to $\Delta\chi^{2}$ from individual data sets. The blue and the
green vertical dot-dash lines indicate the uncertainties at 90% CL determined
with the LM method by requiring $\Delta\chi^{2}+P$ = 100 and with the Hessian
method from the published CT18 NNLO PDFs, respectively.

In the upper-left panel of Fig. 7, we show the results of LM scans on the
gluon PDF at $Q=125$ GeV and $x=0.01$. We find that the HERA inclusive DIS
data and the LHC jet data give the leading constraints. In addition, the CDF
inclusive jet data (Exp. ID = 504) prefers a smaller value of the gluon PDF.
At the global minimum, the $\chi^{2}$ for the CDF inclusive jet data is
elevated by about 20 units. The Hessian method gives a smaller PDF uncertainty
than the estimation based on the LM scans.

In the upper-right panel, we show the results of LM scans on the $s$-quark PDF
at $Q=100$ GeV and $x=0.3$. In this panel, the NuTeV and the CCFR dimuon data
together with the HERA inclusive DIS data give the dominant constraints. These
experimental data are consistent with the global fit. A marked deviation from
the quadratic shape can be observed in the profile of the $\Delta\chi^{2}$. In
this case, a notable difference in uncertainties manifests between the LM
method and the Hessian method, and the LM method should give the more reliable
result.

In the bottom-left panel of Fig. 7, we show the results of LM scans on the PDF
ratio $d/u$ at $Q=100$ GeV and $x=0.3$. The $d/u$ ratio is dominantly
constrained by the LHCb W and Z boson production and the fixed target
experiments BCDMS and NMC. Contrasted with previous situations, the LM method
gives a smaller PDF uncertainty for the $d/u$ ratio.

In the bottom-right panel, we show the results of LM scans on the strangeness
ratio $R_{s}$ at $Q=1.5$ GeV and $x=0.023$. We find that the NuTeV and the
CCFR dimuon data and HERA inclusive DIS data give the dominant constraints. It
is also worthy noting that the NuTeV dimuon data with anti-neutrinos (Exp. ID
= 125) prefers $R_{s}\approx 0.25$ which is smaller than the best fit value
from the global fit $R_{s}$ = 0.52 and results in a large penalty term. At the
global minimum, the $\chi^{2}$ for the NuTeV dimuon data with anti-neutrinos
is elevated by about 7 units. The LM method predicts
$R_{s}=0.52^{+0.35}_{-0.36}$ at 90% CL that has smaller uncertainties than
$R_{s}=0.52^{+0.39}_{-0.41}$ from the Hessian method.

Above scans have also been performed in the CT18 analysis Hou:2019efy , and
our results are consistent with those, which further proves the validity of
our approach. After demonstrating the great efficiency and validity of our
approach on LM scans, we are now ready to perform a systematic study on the
PDF values and ratios for all flavors at a series of $x$ values spreading over
a wide range.

In Fig. 8, we compare the PDF uncertainties at 68% CL at $Q$ = 1.295 GeV
between the LM method and the Hessian uncertainties from the CT18 NNLO PDFs.
The blue and the red solid lines represent the central values of the CT18 NNLO
and the PDFs determined with the aforementioned NNs approximation of
$\chi^{2}$ respectively. The blue and the red hatched areas represent the
uncertainties determined with the Hessian method and the LM method
respectively. The results are normalized to the central value of CT18 NNLO
PDFs. We find good agreements of both the uncertainties and the central values
between the two methods. A notable difference, however, can be seen for $u$,
$d$, $\bar{u}$, $\bar{d}$ and $s$-quark at small-$x$ ($\lesssim 10^{-4}$), as
well as for $d$, $\bar{u}$ and $s$-quark at large-$x$ ($\gtrsim 0.4$). This
indicates a failure of the quadratic approximation in these regions. The
uncertainties from the LM method can be either larger or smaller than the
uncertainties from the Hessian method depending on the flavor and the $x$
value.

In the lower-right panel, we find that $s$-quark PDFs have large uncertainties
for both $x\lesssim 0.001$ and $x\gtrsim 0.4$. This is because the large-$x$
and small-$x$ behavior of the $s$-quark are mostly constrained by the
extrapolation of the PDF parametrization. The error band of $s$-quark of CT18
NNLO PDFs covers negative PDF values at $x\gtrsim 0.4$. This unphysical
behavior implies a limitation of the Hessian method. On the contrary, the
error band determined with the LM method is bounded above zero in all regions.

((a))

((b))

((c))

((d))

((e))

((f))

Figure 8: The parton distribution functions at $Q=1.295$ GeV for $u$, $d$,
$\bar{u}$, $\bar{d}$, $g$ and $s$. The blue and the red solid lines represent
the central values of the CT18 NNLO and the PDFs determined with the
aforementioned NNs approximation of $\chi^{2}$ respectively. The blue and the
red hatched areas represent the uncertainties at 68% CL determined with the
Hessian method and the LM scans respectively.

We also perform the LM scans on the general PDF ratios that is defined as

$R_{f}\equiv\frac{f_{i}(x_{1},Q)}{f_{j}(x_{2},Q)}.$ (11)

The relative uncertainties of $R_{f}$ are calculated at $Q=1.295$ GeV with
$x_{1}$ and $x_{2}$ selected among 12 values from 3$\times 10^{-5}$ to 0.6
listed in Table 3, and $i,j\in\\{g,u,\bar{u},d,\bar{d},s\\}$ runs over all
parton flavors. The results at 90% CL are shown in Fig. 9 (a) and (b) for the
Hessian and LM method respectively. The x and the y axis indicate the
numerator and the denominator, and color code represents the relative
uncertainties of the ratio $R_{f}$. By comparison of the two panels, we find
good agreements between the uncertainties determined with the LM method and
the Hessian method in most regions. Similar to Fig. 8, there are notable
differences at small-$x_{1}$, especially for those with sea quarks in the
numerator.

| 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
$x$ | $3\times 10^{-5}$ | $7\times 10^{-5}$ | $3\times 10^{-4}$ | $7\times 10^{-4}$ | $3\times 10^{-3}$ | $7\times 10^{-3}$
| 7 | 8 | 9 | 10 | 11 | 12
$x$ | 0.01 | 0.02 | 0.06 | 0.1 | 0.2 | 0.6
Table 3: The $x$ values selected for the calculation of the uncertainties of
$R_{f}$.

((a)) Hessian

((b)) LM

Figure 9: The relative uncertainties of $R_{f}=f_{i}(x_{1},Q)/f_{j}(x_{2},Q)$
determined with the Hessian method and the LM method at 90% CL are shown in
panel (a) and (b) respectively. The color code represents the relative
uncertainties of the ratio $R_{f}$. The relative uncertainties are calculated
at $Q=1.295$ GeV with $x_{1}$ and $x_{2}$ selected among 12 values from
3$\times 10^{-5}$ to 0.6 listed in Table 3, and
$i,j\in\\{g,u,\bar{u},d,\bar{d},s\\}$ runs over all parton flavors.

###  4.2 LM scans on cross sections

Higgs bosons are produced dominantly through gluon fusions at the LHC. The
inclusive gluon-fusion cross-section has been calculated to next-to-next-to-
next-to-leading order in QCD Anastasiou:2015vya , which further reduces the
scale variations and makes the PDF uncertainties even more important. Besides,
the production of Higgs boson pair and top-quark pair associated with a Higgs
boson are of equal importance for studies of the Higgs boson self-coupling and
the top-quark Yukawa coupling.

In Fig. 10 we show the results of LM scans on $\sigma_{pp\to hh}$ and
$\sigma_{pp\to ht\bar{t}}$ at $\sqrt{s}$ = 13 TeV or 100 TeV. For
$\sigma_{pp\to hh}$ at $\sqrt{s}$ = 13 TeV, in the upper-left panel, the
behaviors of $\chi^{2}$ are very much similar to that shown in the upper-left
panel of Fig. 7 for the gluon PDF. That is because the cross section of $pp\to
hh$ at 13 TeV is strongly correlated with the gluon PDF at $x\sim 0.02$.
Constraints from HERA inclusive DIS data, BCDMS proton and deuterium data, CMS
8 TeV jet data and ATLAS 8 TeV Z $p_{T}$ data stand out as expected. In
addition, the BCDMS proton data and ATLAS 8 TeV Z $p_{T}$ data both prefer a
larger cross section contrasted with the BCDMS deuterium data which prefers a
smaller value. For $\sigma_{pp\to hh}$ at $\sqrt{s}$ = 100 TeV, in the upper-
right panel, the constraints are distributed among more data sets and are
related to PDFs at small-$x$.

The cross section of $pp\to ht\bar{t}$ mainly depends on the gluon, $u$-quark
and $d$-quark PDFs. For $\sigma_{pp\to ht\bar{t}}$ at $\sqrt{s}$ = 13 TeV, in
the lower-left panel of Fig. 10, similar behaviors of the $\chi^{2}$ as
$\sigma_{pp\to hh}$ are observed. At $\sqrt{s}$ = 100 TeV, the constraints
from HERA inclusive DIS data predominate. In addition, constraints from NuTeV
dimuon data, CMS jet data and BCDMS proton data also play important roles.

((a))

((b))

((c))

((d))

Figure 10: LM scans on the $\sigma_{pp\to hh}$ and $\sigma_{pp\to ht\bar{t}}$
at $\sqrt{s}=$ 13 TeV or 100 TeV.

###  4.3 Study on impact of individual data sets

In order to assess the contribution from an individual experimental data set,
we remove one data set at a time, and repeat the LM scans on physics
quantities with the rest of the data sets. Difference between the fit with and
without the data set can be an assessment of its contribution.

The results for $R_{s}$ at $x=0.023$ and $Q=1.5$ GeV are shown in Fig. 11.
After the removal of each data set, we find that $R_{s}$ value and its
uncertainty are only changed slightly, as represented by those error bars
comparing to the uncertainty from LM scans with the full data set represented
by the gray band. The subtraction of a single data set shows largest effects
for NuTeV dimuon production data (Exp. ID = 124, 125), CCFR dimuon production
data (Exp. ID = 126, 127), E866 Drell-Yan data (Exp. ID = 204) and HERA
inclusive DIS data (Exp. ID = 160). In addition, the NuTeV dimuon production
data and HERA inclusive DIS data prefer a smaller $R_{s}$ contrasted with E866
Drell-Yan data which prefers a larger value, that is consistent with the
bottom-right panel of Fig. 7.

Figure 11: The results of LM scans on the $R_{s}$ ($x=0.023$, $Q=1.5$ GeV)
with data subtracted. The horizontal axis represents the experimental data set
removed from the LM scans. The blue mark and the red error bar respectively
indicate the central value and uncertainties at 90% CL determined with the LM
method with the rest of the data sets. The green hatched area and the gray
band represent the uncertainties at 90% CL determined with the Hessian method
and the LM method with the full data set respectively.

In Fig. 12 we show the results for $\bar{d}/\bar{u}$ at $x=0.3$ and $Q=100$
GeV. The E866 Drell-Yan ratio data (Exp. ID = 203) gives the dominant
constraints. The fit without E866 Drell-Yan ratio data predicts a result of
$\bar{d}/\bar{u}=1.26^{+0.82}_{-0.59}$, while the fit with the full data set
expects $\bar{d}/\bar{u}=1.28^{+0.20}_{-0.33}$. After the inclusion of E866
Drell-Yan ratio data, the uncertainties of $\bar{d}/\bar{u}$ are reduced by
almost 60%. That is because the penalty term of E866 Drell-Yan ratio data
provides a strong constraint on $\bar{d}/\bar{u}$. In addition, constraints
from NMC deuteron data (Exp. ID = 104) and HERA inclusive DIS data also play
important roles.

Figure 12: The same as Fig.11, but for the results of LM scans on the
$\bar{d}/\bar{u}$ ( $x=0.3$ and $Q=100$ GeV).

In Fig. 13 we show the results for $\sigma_{pp\to hh}$ at $\sqrt{s}$ = 13 TeV.
The constraints from HERA inclusive DIS data predominate as expected. In
addition to that, constraints from BCDMS proton and deuterium data (Exp. ID =
101, 102) and ATLAS 8 TeV Z $p_{T}$ data (Exp. ID = 253) also play important
roles. The fit without HERA inclusive DIS data expects $\sigma_{pp\to
hh}=0.0129^{+0.0007}_{-0.0009}$ pb, while the fit with the full data set gives
$\sigma_{pp\to hh}=0.0131^{+0.0005}_{-0.0007}$ pb. An upward shift of about
$2\times 10^{-4}$ pb is observed when we incorporate HERA inclusive DIS data,
and the uncertainties of $\sigma_{pp\to hh}$ are reduced by almost 20%. In
addition, the HERA inclusive DIS data and BCDMS deuterium data both prefer a
larger $\sigma_{pp\to hh}$ contrasted with BCDMS proton data and ATLAS Z
$p_{T}$ data which prefer a smaller value, that is consistent with the upper-
left panel of Fig. 10.

Figure 13: The same as Fig.11, but for the results of LM scans on the
$\sigma_{pp\to hh}$ at $\sqrt{s}$ = 13 TeV.

###  4.4 Two-Dimensional LM scans

Besides the PDF uncertainties, it is possible to quantify other statistical
estimators such as the correlation between two physics quantities with two-
dimensional LM (2-D LM) scans. That can be achieved by adding a second physics
quantity into Eq. (9). The new function that needs to minimized in the global
fit becomes

$\Psi\left(\lambda_{1},\lambda_{2},\left\\{a_{i}\right\\}\right)\equiv\chi^{2}\left(\left\\{a_{i}\right\\}\right)+\lambda_{1}{X_{1}}\left(\left\\{a_{i}\right\\}\right)+\lambda_{2}{X_{2}}\left(\left\\{a_{i}\right\\}\right),$
(12)

where $\lambda_{1}$ and $\lambda_{2}$ are specified constants, and
$X_{1}(\\{a_{i}\\})$ and $X_{2}(\\{a_{i}\\})$ represent the two physics
quantities of interest. Similar to Eq. (9), the constrained minimum of
$\chi^{2}$ from the global fit depends on $X_{1}$ and $X_{2}$, and can be
written as $\chi^{2}=\chi^{2}_{min}+\Delta\chi^{2}$, where
$\chi^{2}_{min}=\chi^{2}(\lambda_{1}=0,\lambda_{2}=0)$. The contour of
$\Delta\chi^{2}+P$ in the plane of $X_{1}$ vs. $X_{2}$ can be an assessment of
the correlation between $X_{1}$ and $X_{2}$.

As examples in Fig. 14 we show contours for $\bar{d}/\bar{u}$ ($x=0.3$,
$Q=100$ GeV) vs. $R_{s}$ ($x=0.023$, $Q=1.5$ GeV) and $\sigma_{pp\to hh}$
($\sqrt{s}$ = 13 TeV) vs. $\sigma_{pp\to hh}$ ($\sqrt{s}$ = 100 TeV)
determined with the 2-D LM scans. In the left panel, a weak correlation
between strangeness ratio $R_{s}$ and $\bar{u}/\bar{d}$ ratio is observed.
That is because the two quantities are dominantly constrained by different
experimental data sets. At small $\Delta\chi^{2}+P$ the contour shows an
elliptic shape. When $\Delta\chi^{2}+P$ gets larger, the shape of the contour
becomes irregular due to the increase of penalty term contributions. On the
contrary, the right panel demonstrates a strong correlation between
$\sigma_{pp\to hh}$ ($\sqrt{s}$ = 13 TeV) and $\sigma_{pp\to hh}$ ($\sqrt{s}$
= 100 TeV) since both processes are sensitive to gluon PDFs and constrained by
the relevant experimental data sets.

((a))

((b))

Figure 14: Contour plot of $\Delta\chi^{2}$ plus Tier-2 penalty term on the
plane of $\bar{d}/\bar{u}$ ($x=0.3$, $Q=100$ GeV) vs. $R_{s}$ ($x=0.023$,
$Q=1.5$ GeV) and $\sigma_{pp\to hh}$ ($\sqrt{s}$ = 13 TeV) vs. $\sigma_{pp\to
hh}$ ($\sqrt{s}$ = 100 TeV).

##  5 Applications

In this section, we evaluate the impact of the NOMAD measurements and of two
pseudo-data sets of HL-LHC on PDFs based on the new approach. In addition, we
study constraints on the new physics with a joint fit of both PDFs and the
Wilson coefficient of lepton-quark contact interactions in the framework of
the SMEFT.

###  5.1 Constraint from NOMAD data

The charm-quark production in CCDIS process provides a unique sensitivity to
the strange-quark distribution in the nucleon, with a clean signal of two
muons with opposite charges in the final state. Recently, NOMAD collaboration
reported a measurement of dimuon production in the neutrino-iron scattering
experiment NOMAD:2013hbk . A sample of about $9\times 10^{6}$ inclusive CCDIS
events, including 15344 dimuon events, is collected, providing a reduced
statistical uncertainty. Observables are taken to be the ratios of dimuon to
inclusive cross-sections, which provides a large cancellation of the common
systematic uncertainties presented in both the numerator and the denominator.
Final results are distributed among three differential variables: the
reconstructed neutrino energy $E_{\nu}$, the Bjorken $x$ and the partonic
center of mass energy $\sqrt{\hat{s}}$. By the supplement of data from NOMAD,
the improvement in the constraint on $s$-quark PDFs are studied in this
section using the same NNs approach on $\chi^{2}$ mentioned in previous
sections.

On the theoretical side, structure functions in S-ACOT-$\chi$ general mass
scheme up to NNLO are constructed, so that a full consideration of the charm-
quark mass is included Berger:2016inr ; Gao:2017kkx ; Gao:2021fle .
Predictions of inclusive CCDIS and open charm production cross-sections are
made from these constructions, and dimuon cross sections are derived by
further applying the inclusive decay branching ratio of charm quark to muon.
The significant uncertainties of the decay branching ratio contribute as one
of the dominant systematic errors on the dimuon cross sections, which are
summarized in Appendix C.

In Fig. 15 we show comparison of NOMAD data and our predictions at both NLO
and NNLO, as well as the Hessian PDF uncertainties at 68% CL for distributions
over $E_{\nu}$ or $x$. The PDF uncertainties can be as large as 10% in most
regions. This directly comes from the large uncertainties of the predictions
of dimuon cross-sections, and can be further traced back to the poor knowledge
about $s$-quark PDFs. In both distributions, most of the data points are
consistent with our predictions, while a significant deviation can be found in
the last two points of the distribution over Bjorken $x$. That can be due to
the modeling of heavy nuclear corrections used in the experimental analysis.
We will discard those two data points when including NOMAD data in our later
global fit. It is also noted that the inclusion of NOMAD data to the global
fit can improve the consistency with almost no cost of tension with the other
data sets Alekhin:2014sya ; Faura:2020oom . Most of the data lie above our NLO
predictions of central values. Given this fact, an increased $s$-quark PDF is
expected after the inclusion of NOMAD data, and this increase gets larger due
to the negative corrections from NNLO.

((a))

((b))

Figure 15:  NLO (blue line) and NNLO (red line) predictions for ratios of
dimuon to CCDIS inclusive differential cross-sections with respect to neutrino
energy (panel a) and Bjorken $x$ (panel. b). The blue hatched areas represent
the Hessian PDF uncertainties of the NLO predictions at 68% CL. NOMAD data are
also shown with statistical uncertainties and the combination of both
statistical and systematic uncertainties.

((a))

((b))

Figure 16: PDF induced correlations between theory predictions for different
experimental bins, for NOMAD distribution in neutrino energy (a) and in
Bjorken-$x$ (b), calculated with CT18 NNLO Hessian PDF set. Numbers in the
axis represent the center of each bin, and numbers in the table represent the
correlation cosine for each pair of bins.

As mentioned earlier, NOMAD presents measurements on three distributions. The
different sensitivities of distributions over $E_{\nu}$ or $x$ are illustrated
in Fig. 16. It shows the PDF induced correlations among bins of each
distributions calculated with CT18 NNLO Hessian PDFs. We find that for
$E_{\nu}$ distribution, all data points are strongly correlated, and similar
results are found for $\sqrt{\hat{s}}$ distribution which is not shown here.
Both of them only impose constraints on the overall normalization of the
$s$-quark distribution. Thus their constraints are diluted due to the
systematic errors on the inclusive branching ratio of charm quark to muon
(0.094$\pm$0.01). On the other hand, the correlation pattern is nontrivial for
$x$ distribution which imposes further constraints on the shape of $s$-quark
PDFs. We can not simply combine all these distributions from NOMAD data due to
the lack of public statistical correlation between these distributions. Hence
in the following, only the $x$ distribution is included in our global
analysis.

In Fig. 17, we compare $u$, $\bar{u}$, $\bar{d}$ and $s$-quark PDFs at
$Q=1.295$ GeV from fits with and without the inclusion of NOMAD data. The PDF
uncertainties are shown through hatched areas with relevant colors. NOMAD data
are taken from the distribution over Bjorken $x$ excluding the last two
points, with predictions calculated up to NNLO in QCD. Predictions for data
sets 124-127 (dimuon measurements from NuTeV and CCFR) in the global fit are
replaced with their NNLO versions when including NOMAD data, in order to match
on the theoretical precision. Note in the fit without NOMAD data the
predictions for data sets 124-127 are evaluated at NLO similar to those in
CT18. All PDFs are normalized to the central value without NOMAD data in Fig.
17. In the upper-left panel, almost no change occurs in the region $x\gtrsim
0.1$ of $u$-quark PDF, and a negligible downward shift smaller than 2% can be
seen for $x\lesssim 0.05$. Slight downward shifts on both central value and
uncertainty region can also be observed in the $\bar{u}$ (upper-right panel)
and the $\bar{d}$-quark (lower-left panel) PDFs. The downward shifts observed
are required to stabilize the $W$ and $Z$ production cross sections at
collider experiments. The improvement in the constraints on $u$, $\bar{u}$ and
$\bar{d}$-quark PDFs are smaller than about 3%. This insensitivity of $u$,
$\bar{u}$ and $d$-quark PDFs to NOMAD data is an indication of the CKM
suppression in the charm-quark production. The constraint on $s$-quark PDF is,
however, markedly improved around $x=0.05$. In the region of $x\sim 0.05$, the
$s$-quark PDF achieves a factor of two better precision when NOMAD data are
incorporated. This is because NOMAD data peak at neutrino energy
$E_{\nu}\approx 30$ GeV, which implies a sensitivity to kinematic region with
Bjorken $x\sim 1/(1+2M_{nucleon}E_{\nu}/Q^{2})\sim 0.03$ at $Q=1.295$ GeV. An
upward shift of more than 15% is also observed in most regions. It is indeed a
manifestation of the trend of prediction-data comparison shown in Fig. 15.

Both ABM and NNPDF groups considered the impact of NOMAD data Alekhin:2014sya
; Faura:2020oom . As to the analysis of ABM group, an at most 5% downward
shift is reported near region $x\approx 0.05$ at scale $Q=3$ GeV when NOMAD
data are incorporated into the fit with only NuTeV/CCFR data (data sets
124-127 in this paper) Alekhin:2014sya . More data sets are considered in the
work of NNPDF group Faura:2020oom . With the analysis performed there, NOMAD
data together with ATLAS W/Z data sets ATLAS:2012sjl ; ATLAS:2016nqi
contribute to a marked enhancement of $s$-quark PDF at $Q=10$ GeV compared
with CT18 data sets. It is noted that, between these two kinds of data sets,
ATLAS W/Z data sets are already reported to give a larger $s$-quark PDF
compared with CT18 data sets Hou:2019efy , and the work of NNPDF group further
demonstrated that ATLAS W/Z data sets prefer a larger $s$-quark PDF compared
with NOMAD data. Finally, both the two groups and our analysis indicate strong
constraints on $s$-quark PDF in the region near $x\approx 0.05$, given the
incorporation of NOMAD data.

((a))

((b))

((c))

((d))

Figure 17:  The parton distribution functions at $Q=1.295$ GeV for $u$,
$\bar{u}$, $\bar{d}$, and $s$. The red and the blue solid lines represent the
central values with and without NNLO NOMAD data respectively, and the red and
the blue hatched areas represent the respect uncertainties at 68% CL. When
NOMAD (NOM) data are incorporated, 124-127 data sets (NUT) are replaced with
their NNLO version. Central values are normalized to the NLO version.

We also compare the sensitivities to $s$-quark PDF between NOMAD data and the
other data sets in Fig. 18, in which we show LM scans on $s$-quark PDF and
$R_{s}$. This comparison is set at a scale of $Q=1.5$ GeV and $x=0.1$ in panel
(a). In this panel, NOMAD data predominate over the other experimental data
sets. When $x$ gets smaller to be 0.023 as in panel (b), NuTeV and CCFR
neutrino DIS experiments become more important but still the NOMAD data show
the most prominence. In panel (a), the fit without NOMAD data predicts
$R_{s}(x=0.1,Q=1.5\rm{GeV})=0.40^{+0.35}_{-0.20}$ at 90% CL, while fit
including NOMAD data expects
$R_{s}(x=0.1,Q=1.5\rm{GeV})=0.54^{+0.24}_{-0.06}$, giving improved constraints
by a factor of two. In panel (b), the corresponding values are
$R_{s}(x=0.023,Q=1.5\rm{GeV})=0.53^{+0.33}_{-0.38}$ and
$R_{s}(x=0.023,Q=1.5\rm{GeV})=0.70^{+0.40}_{-0.17}$, respectively. NOMAD data
hence give about 20% reduction on PDF uncertainties. It is also noted that a
slight tension exists between NOMAD data and data from the other two neutrino
DIS experiments, i.e., NuTeV and CCFR, in both panels. The latter two
experiments both prefer smaller $R_{s}$s contrasted with NOMAD data which
prefer a larger one. Further investigations on the interplay of the three
experiments and of different theories are included in Appendix C. Moreover,
the ATLAS W/Z data ATLAS:2012sjl , which prefer an especially larger
$R_{s}(x=0.023,Q=1.38\rm{GeV})\sim 1$, show an even stronger tension with
these two neutrino DIS experiments. NOMAD data, however, compromise between
these two extremes. This conclusion is also observed in the analysis of
Faura:2020oom . Meanwhile, a similar result of
$R_{s}(x=0.023,Q=1.6\rm{GeV})=0.71\pm 0.1$ is obtained in that work once NOMAD
data are included.

On the other hand, we let the scale increase to be $Q=100$ GeV in panel (c)
and panel (d). The case with $x=0.3$ shows more sensitive than that with
$x=0.002$. In panel (d), it can be seen that NuTeV and CCFR data become
comparable with NOMAD data. No significant shift in the central value is found
when we incorporate NOMAD data, but an almost 30% better constraints on
$s$-quark PDF is achieved. In panel (c), the sensitivity of NOMAD data becomes
worse due to the favor of large-$x$ at this scale, and collider data now play
important roles. Only improvement of a few percent in the constraint on
$s$-quark PDF can be obtained.

((a))

((b))

((c))

((d))

Figure 18:  LM scans on the $R_{s}$ at $Q=1.5$ GeV and $x=0.1$ or 0.023 (upper
panels), and LM scans on the $s$-quark at $Q=100$ GeV and $x=0.002$ or 0.3
(lower panels). The blue and the green vertical solid (dot-dash) lines
represent the central values (uncertainties) with and without NOMAD data
respectively.

###  5.2 Impact of High-luminosity LHC

LHC data play important roles on constraining PDFs as shown in Table 2. And
the upgrade of the LHC, the HL-LHC, is expected to accumulate a total
integrated luminosity of $\mathcal{L}=3000$ fb-1 for ATLAS and CMS and $300$
fb-1 for LHCb. In this section, we take two of those HL-LHC pseudo-data sets
constructed in Ref. AbdulKhalek:2018rok , and evaluate their impact on PDFs
within the framework of CT18 based on our new approach.

The HL-LHC pseudo-data are generated for processes of Drell-Yan production
with high dilepton invariant mass and W and Z boson production in the forward
region. Details of these pseudo-data are described as follows:

$\bullet$ The distribution of dilepton invariant mass $d\sigma(pp\to
l^{+}l^{-})/dm_{ll}$ of high-mass Drell-Yan process at $\sqrt{s}=$14 TeV,
covered by the ATLAS experiment, is generated according to the following
requirements: $p_{T}^{l_{1(2)}}\geq$ 40 (30) GeV, $|\eta^{l}|\leq$ 2.5, and
$m_{ll}\geq$ 116 GeV. The total number of data points is 21. The binning and
the systematic uncertainties are determined from Refs. ATLAS:2016gic ;
AbdulKhalek:2018rok .

$\bullet$ The distributions for W and Z boson production in the forward region
at $\sqrt{s}=$14 TeV, covered by the LHCb experiment, are generated according
to the following cuts: $p_{T}^{l}\geq$ 20 GeV, 2.0 $\leq\eta^{l}\leq$ 4.5. An
additional requirement for Z production is that 60 GeV $\leq m_{ll}\leq$ 120
GeV. The total number of data points is 90. The binning and the systematic
uncertainties are determined from Refs. LHCb:2015mad ; AbdulKhalek:2018rok .

We include those pseudo-data in the CT18 global fit and quantify their impact
on PDFs. In Fig. 19 we show a comparison of the PDFs with and without HL-LHC
pseudo-data, together with the published Hessian set of CT18. All results are
normalized to the central value of CT18. The PDF uncertainties are shown
through hatched areas with relevant colors. In Fig. 19, a significant
reduction in PDF uncertainties can be found in all cases once including the
pseudo-data, especially for sea quarks. In the upper-left panel, the PDF
uncertainties are reduced by almost a factor of 2, from about 30% to about
15%, at small-$x$. Similar improvements can also be observed in $d$, $\bar{u}$
and $s$-quark PDFs. That is because HL-LHC pseudo-data contribute a great
improvement in statistics, and cover the kinematic regions where PDFs are not
determined well. Specifically, the process of high-mass Drell-Yan is directly
sensitive to sea quarks at large-$x$, and the process of forward W/Z
production constrains the $s$-quark PDF at both small-$x$ and large-$x$. In
the lower-right panel, we find that the HL-LHC gives about 30% reduction on
PDF uncertainties of $s$-quark in the regions of $x\sim 0.01$ and $x\sim 0.1$.
This result highlights the importance of the process of forward W/Z
production.

((a))

((b))

((c))

((d))

Figure 19: The parton distribution functions at $Q=1.295$ GeV for $u$, $d$,
$\bar{u}$ and $s$. The red and the green solid lines represent the central
values without and with HL-LHC pseudo-data respectively, and the red and the
green hatched areas represent the respective uncertainties at 68% CL. The
results are normalized to the central value of CT18 NNLO (blue solid line).

We show the results of LM scans on $R_{s}$ and $u/d$ ratio in Fig. 20. We find
measurements of high-mass Drell-Yan process and forward W/Z production process
at HL-LHC give strong constraints on $R_{s}$ and $u/d$ ratio at both small-$x$
and large-$x$. The PDF uncertainties of $R_{s}$ and $u/d$ are significantly
reduced after the inclusion of pseudo-data.

((a))

((b))

((c))

((d))

Figure 20: LM scans on the $R_{s}$ at $Q=1.5$ GeV and $x=0.023$ or 0.1 (upper
panels), and LM scans on the $d/u$ at $Q=100$ GeV and $x=0.002$ or 0.3 (lower
panels). The blue and the green vertical solid (dot-dash) lines represent the
central values (uncertainties) with and without HL-LHC pseudo-data
respectively.

In the upper-left panel of Fig. 20 for $R_{s}$ at $x=0.023$ and $Q=1.5$ GeV,
the constraints from HL-LHC pseudo-data predominate as expected. In addition
to that, constraints from NuTeV dimuon, CCFR dimuon and HERA inclusive DIS
data also play important roles. The fit without pseudo-data predicts a result
of $R_{s}$ = $0.53^{+0.33}_{-0.38}$ at 90% CL, while fit including pseudo-data
gives $R_{s}$ = $0.54^{+0.22}_{-0.19}$. After the inclusion of pseudo-data,
the PDF uncertainties are reduced by almost 50%. As $x$ increases to 0.1 in
the upper-right panel, HL-LHC forward W/Z data becomes more important. Fit
without pseudo-data gives a result of $R_{s}$ = $0.40^{+0.35}_{-0.20}$, while
fit including pseudo-data gives $R_{s}$ = $0.39^{+0.16}_{-0.16}$.

For $d/u$ at $x=0.002$ and $Q=100$ GeV, in the lower-left panel, the most
strong constraints originate from HL-LHC pseudo-data together with LHC W and Z
boson data and the fixed target experiments E866 and NMC. The results of $d/u$
are $0.946^{+0.038}_{-0.032}$ and $0.954^{+0.026}_{-0.032}$ corresponding to
fit without and with pseudo-data respectively. After the inclusion of pseudo-
data, the PDF uncertainties are reduced by almost 30%. In the lower-right
panel for $x=0.3$ and $Q=100$ GeV, pseudo-data predominate over the other
experimental data sets. Fit without and with pseudo-data give $d/u$ =
$0.376^{+0.028}_{-0.021}$, and $d/u$ = $0.385^{+0.020}_{-0.016}$ respectively.
PDF uncertainties are reduced by almost 25% in this case. Both of the two HL-
LHC processes prefer a larger $d/u$, and their inclusion leads to an increase
of the central value.

In Fig. 21, we show the results for the general PDF ratio $R_{f}$ as defined
in Eq. (11). The uncertainties of $R_{f}$ are also determined with the LM
method. In panel (a), we show the relative uncertainties of $R_{f}$ at 90% CL,
$\Delta R_{f}$, from the fit with inclusion of HL-LHC pseudo-data. To compare
with the results from the fit without HL-LHC pseudo-data, a reduction factor
of $\Delta R_{f}$,

$y_{red}=2\dfrac{\Delta R_{f}^{base}-\Delta R_{f}^{base+HL}}{\Delta
R_{f}^{base}+\Delta R_{f}^{base+HL}},$ (13)

is shown in panel (b). We find that the relative uncertainties of $R_{f}$ have
a noticeable reduction in general. For the case of the $u$ and $\bar{u}$-quark
PDFs as the numerator, we find that the HL-LHC gives about 80% reduction on
relative uncertainties in the region of $x_{1}\lesssim 0.001$, which is
because the HL-LHC pseudo-data give strong constraints on $u$ and
$\bar{u}$-quark PDFs in this region as shown in Fig. 19. In addition, for the
case of the $s$-quark PDFs as the numerator, we find that the HL-LHC gives
about 50% reduction on relative uncertainties in the region of $x_{1}\sim
0.01$. However, for the $f_{d}(x_{1},Q)/f_{u}(x_{2},Q)$,
$f_{d}(x_{1},Q)/f_{\bar{u}}(x_{2},Q)$, $f_{\bar{d}}(x_{1},Q)/f_{u}(x_{2},Q)$
and $f_{\bar{d}}(x_{1},Q)/f_{\bar{u}}(x_{2},Q)$, we find that the reduction
factors on relative uncertainties are minor in the region of $x_{1}\lesssim
1\times 10^{-3}$ and $x_{2}\lesssim 1\times 10^{-3}$. That is because the
correlation between $u$-quark PDFs and $d$-quark PDFs at small-$x$ that
originates from the parametrization form of PDFs. Besides, for the ratios of
gluon PDFs $f_{g}(x_{1},Q)/f_{g}(x_{2},Q)$, the uncertainties are reduced by
only a few percent, which is expected due to the weak correlations between the
two HL-LHC processes and gluon PDFs.

((a)) CT18NNLO + HL-LHC

((b)) Reduction

Figure 21: The relative uncertainties of $R_{f}=f_{i}(x_{1},Q)/f_{j}(x_{2},Q)$
determined with the LM method at 90% CL are shown in panel (a), where
$Q=1.295$ GeV. In panel (b) we show the reduction factors on the relative
uncertainties of $R_{f}$ between with and without HL-LHC pseudo-data.

###  5.3 Constraint on new physics with the global fit

PDFs and their uncertainties play important roles in the indirect searches for
new physics beyond the SM. In this case, the scale of the new physics can be
well beyond the typical scale of hard scatterings, and its effects can be
formally described in the framework of the SMEFT. PDFs are determined by
fitting to a variety of experimental data under the assumption of the SM. This
leads to a problem that the degeneracy of PDF variations and the new physics
contributions cannot be identified. Therefore, to assess and furthermore
constrain the new physics, a joint global fit including both PDFs and model
parameters of new physics should be performed. In this paper, we only consider
one dimension-six operator, namely the lepton-quark contact interactions, to
model the BSM effects in the SMEFT framework,

$\displaystyle\mathcal{L}_{\mathrm{SMEFT}}$
$\displaystyle=\mathcal{L}_{\mathrm{SM}}+\sum_{i,j}\dfrac{c_{ij}}{\Lambda^{2}}(\bar{q_{i}}\gamma_{\mu}q_{i})(\bar{l_{j}}\gamma^{\mu}l_{j})$
$\displaystyle=\mathcal{L}_{\mathrm{SM}}+\dfrac{\tilde{c}}{\Lambda^{2}}\sum_{i,j}e_{q_{i}}e_{l_{j}}(\bar{q_{i}}\gamma_{\mu}q_{i})(\bar{l_{j}}\gamma^{\mu}l_{j}),$
(14)

where $c_{ij}$ is the Wilson coefficient, $l_{j}$ and $q_{i}$ represent fields
of charged leptons and quarks of flavor $j$ and $i$ respectively, and
$e_{q_{i}(l_{j})}$ are the corresponding electric charges. We assume the new
interactions being vector-current type and have a flavor structure similar to
the QED coupling for simplicity. Thus the contributions from the new physics
are parametrized by a single variable of the effective Wilson coeﬀicient
$\tilde{c}$ that is normalized to the QED coupling.

In the case of data sets of the CT18 global fit, the DIS and the Drell-Yan
processes receive contributions from this operator. Processes with relatively
large $Q^{2}$ are especially sensitive to BSM effects, where $Q$ is the
momentum transfer. Most of the data of Drell-Yan process included in the CT18
analysis are collected near the Z-pole region, which is less sensitive to new
physics. Hence, we only consider the HERA DIS process due to its large
$Q^{2}$. The amplitude of SM contributions from QED interactions is
proportional to $1/Q^{2}$, and the amplitude of the BSM contributions is
proportional to $\tilde{c}/\Lambda^{2}$ 111The weak interactions from $Z$
boson induce a different energy dependence of $1/(Q^{2}+M_{Z}^{2})$ which we
neglect here for simplicity.. Hence, the total cross section including the BSM
effects can be written as:

$\sigma_{\mathrm{total}}=(1+\frac{\tilde{c}}{\Lambda^{2}}Q^{2})^{2}\times\sigma_{\mathrm{DIS}}.$
(15)

A new NNs is built by adding the parameter $\tilde{c}/\Lambda^{2}$ into the
input layer. An association between the 29 variables
{$a_{i},\tilde{c}/\Lambda^{2}$} and $\chi^{2}$ is constructed. With the new
NNs, $\chi^{2}$ is recalculated and the results of LM scans on
$\tilde{c}/\Lambda^{2}$ are shown in Fig. 22. HERA inclusive DIS data give the
dominant constraints as expected. The LM scans predict a result of
$\tilde{c}/\Lambda^{2}=0.56^{+9.16}_{-9.16}$ $\mathrm{TeV^{-2}}$ at 90% CL,
which is consistent with the SM. The interplay between PDFs and BSM effects in
the framework of the SMEFT has been studied in previous works Carrazza:2019sec
; Greljo:2021kvv ; Iranipour:2022iak . A simultaneous determination of the
PDFs and BSM effects from DIS data based on the NNPDF framework was presented
in Ref. Carrazza:2019sec . The Wilson coefficients of the lepton-quark contact
interactions (i.e. $l$-$u$, $l$-$d$, $l$-$s$ and $l$-$c$ contact interactions)
are constrained by the HERA inclusive DIS data. The most stringent bounds are
obtained for $u$-quark, followed by $d$-quark, and then $c$-quark and
$s$-quark. The constraint on the Wilson coefficients for $u$-quark converted
to $\tilde{c}/\Lambda^{2}$ is [-6.5 $\mathrm{TeV^{-2}}$, 39.2
$\mathrm{TeV^{-2}}$] at 90% CL. In Ref. Greljo:2021kvv , the BSM effects are
constrained by the high-mass Drell-Yan data. The result converted to
$\tilde{c}/\Lambda^{2}$ is [-6.5 $\mathrm{TeV^{-2}}$, 57.8
$\mathrm{TeV^{-2}}$] at 95% CL.

Figure 22: LM scans on $\tilde{c}/\Lambda^{2}$. The blue vertical solid line
represents the central value of $\tilde{c}/\Lambda^{2}$, and the blue vertical
dot-dash lines represent the uncertainties at 90% CL.

In Fig. 23, we compare $u$, $d$, $\bar{u}$ and $g$ PDFs at $Q$ = 1.295 GeV
determined by fitting with and without the new physics contributions. The PDF
uncertainties are shown through hatched areas with relevant colors. The PDFs
from two fits are almost indistinguishable for both central value and the
uncertainties. In addition, we find that the central value is slightly changed
if the $\tilde{c}/\Lambda^{2}$ is fixed at -8.60 or 9.72 TeV-2, as represented
by the two black solid lines. Specifically, in the upper-left panel, a shift
as large as 2% can be observed in the region of $x\sim 0.02$. Similar shifts
can also be observed in panel (b) and panel (c). Besides, in the lower-right
panel, a shift as large as 10% can be observed at both small-$x$ ($\sim
10^{-4}$) and large-$x$ ($\sim 0.6$). These shifts on PDFs are required to
compensate for the contributions from the new physics on DIS cross sections.
Our approach can be extended to include more EFT operators from new physics
which we leave for future studies.

((a))

((b))

((c))

((d))

Figure 23: The parton distribution functions at $Q=1.295$ GeV for $u$, $d$,
$\bar{u}$ and $g$. The blue and the red solid lines represent the central
values determined by fitting without and with the new physics contributions
respectively, and the blue and the red hatched areas represent the respective
uncertainties at 68% CL. The black solid lines represent the PDFs when
$\tilde{c}/\Lambda^{2}$ is fixed at -8.60 or 9.72 TeV-2.

##  6 Conclusion

Better understanding on parton distributions is essential for precision
physics at hadron colliders, as well as for study of QCD. Nowadays the
analysis of PDFs requires calculations of the log-likelihood functions
$\chi^{2}$ from thousands of experimental data points, and scans of multi-
dimensional parameter space with tens of degrees of freedom. Such analyses
will benefit from development of new methods and improvement of computing
efficiencies, for instance by various interpolation approaches. In this paper
we propose a new approach of using Neural Networks and machine learning
techniques to model the dependence of the $\chi^{2}$ or any physics quantities
on the PDFs. We demonstrate the high accuracy of our approach through detailed
comparisons in the PDF parameter space of interest, taking the CT18 NNLO
analysis as an example. Importantly, compared with direct calculations the
computational cost on calculating $\chi^{2}$ are reduced by several orders of
magnitude. The improvement ensures efficient scans of the full PDF parameter
space and is desirable for the determination of PDF uncertainties.

Based on our NNs, we perform a series of LM scans to reevaluate PDF
uncertainties in the CT18 NNLO analysis, and to understand the interplay
between different data sets. The LM method is generally more reliable through
a scan of the $\chi^{2}$ along the trajectory of constrained minimum of the
physics quantity studied. Our new approach renders such extensive scans almost
costless and ensures the possibility of detailed comparisons of PDF
uncertainties determined from the LM and Hessian method. We first perform LM
scans on PDF values and ratios at various $x$ and $Q$ values, and find the
results from the LM method and the Hessian method agree well in general.
However, a notable difference can be observed in the small and the large-$x$
regions. Since the quadratic approximation fails in the region where PDF
uncertainties are large, and the results from the LM method are more reliable.
Besides, we perform LM scans on the production cross sections of the Higgs
boson pair and the top-quark pair in association with a Higgs boson, at the
LHC or future colliders, as well as two dimensional scans on a pair of PDFs or
cross sections. Furthermore, using LM scans we study the impact of individual
data sets in the CT18 NNLO analysis by subtracting and adding back one data
set at a time.

We show further applications of our approach on several extensions of the CT18
NNLO analysis. Especially, we study the impact of the NOMAD dimuon data on
constraining the strange-quark PDFs. Theoretical predictions are calculated in
the S-ACOT-$\chi$ general mass scheme up to NNLO, based on which the NNs are
constructed and LM scans are performed. We find that the NOMAD data place
stringent constraints on the strange-quark PDFs at intermediate and large-$x$
regions. At $x\sim 0.05$, for example, the PDF uncertainties of the strange
quark are reduced by almost a factor of 2. An upward shift of more than 15% in
the strange-quark PDF as well as slight downward shift in the $u$ and
$d$-quark PDFs are also observed in most regions. We show the interplay of the
NOMAD data and other data sets in the CT18 by detailed LM scans on $R_{s}$ and
$s$-quark PDFs at different scales and $x$ values. The global fit with NOMAD
data predicts $R_{s}(x=0.023,Q=1.5\rm{GeV})=0.70^{+0.40}_{-0.17}$ at 90% CL
and a slight tension between NOMAD and NuTeV data is observed. We also present
a series of variant fits for clarifications on the impact of different theory
predictions and of different choices of the decay branching ratio of the charm
quark.

Afterwards, we study the impact of two HL-LHC pseudo-data constructed in Ref.
AbdulKhalek:2018rok , including the high-mass Drell-Yan data and the forward
W/Z production data. We find potentially large reduction on PDF uncertainties
of the sea quarks. These results highlight the importance of HL-LHC
measurements. Besides, we performed a joint fit on both PDFs and effects of
new physics beyond the SM. We take the lepton-quark contact interactions as an
example that are described by high dimensional operators in the SMEFT. We
determine the effective Wilson coefficient to be
$\tilde{c}/\Lambda^{2}=0.56^{+9.16}_{-9.16}$ TeV-2 at 90% CL as mostly
constrained by the HERA inclusive DIS data. Foreseen extensions of the study
would be to include more SMEFT operators in the joint fit that is under
investigation.

###### Acknowledgements.

This work was sponsored by the National Natural Science Foundation of China
under the Grant No. 11875189 and No.11835005. JG would like to thank members
of CTEQ-TEA collaboration for helpful discussions and proofreading of the
manuscript.

##  Appendix A More on the Neural Network approach

In this appendix we collect various details of the NN approach, including on
the architectures and parametrization dependence, the generation of training
and test samples, and the performances in terms of computational cost. One
important feature of our approach is to use directly PDF values as inputs to
the NNs rather than the PDF parameter themselves. That ensures a great
flexibility of the functional space since we can select PDF values at an
arbitrary number of $x$ points. In our current study with CT18 parametrization
form, we select the $x$ grid consisting of 14 points for each PDF flavor with
their values shown in Table 4. They are selected randomly with the only
criteria being distributed evenly in $\ln x$.

We explain briefly on the mathematical model behind our NNs. The true
dependence of our target function, for instance, the $\chi^{2}$, on the PDF
parameters is uniquely determined by the theory and experimental data used in
the global fit, which we denote as $A_{\rm TR}$. On another hand if we
exchange the PDF parameters by the PDF values at discrete $x$ points, the
mapping is not unique since we have input PDF values far more than the number
of PDF parameters. Thus we arrive at a bunch of possible functions
$\\{A^{*}_{\rm TR}\\}$ depending explicitly on $\\{I_{k}\\}$ which is the PDF
value at the $k_{th}$ node, satisfying

$\chi^{2}_{truth}=A_{{\rm TR}}(a)=A^{*}_{{\rm TR}}(\\{I_{k}(a)\\}).$ (16)

The purpose of our NNs is to construct an explicit function of $\\{I_{k}\\}$
depending on a set of tuneable parameters $t_{\beta}$. By the training
procedure we update $A_{\rm NN}$ iteratively until it converges to the
neighborhood of one of the truth function $A^{*}_{\rm TR}$ with a choice on
the parameters $\hat{t}_{\beta}$. Finally we arrive at our approximation to
the $\chi^{2}$ dependence on the PDF parameters as

$\chi^{2}_{pred}=A_{{\rm NN}}(\\{I_{k}(a)\\};\,\\{\hat{t}_{\beta}\\}).$ (17)

As from above one expects that the outcome NN (or equivalently the solution
$\hat{t}_{\beta}$) in general depends on the parametrization form of PDFs.
However, in practice one can approximate either PDFs or cross sections in
terms of interpolated functions on a dense $x$-grid with sufficient accuracy,
as implemented successfully in APPLgrid Carli:2010rw , FastNLO Kluge:2006xs ,
and FastKernal Ball:2021leu . In that sense there may exist an almost
universal solution for different parametrization forms if one start with a
sufficiently large number of PDF inputs. We leave that for future
investigations.

((a))

((b))

((c))

((d))

Figure 24: The predictions to truths ratios of $\chi^{2}$ for measurement of
the inclusive jet cross sections at $\sqrt{s}=$ 7 TeV by ATLAS when training
the NNs to individual PDF parametrization forms including CT18, MSHT14 and
NNPDF4.0, or an ensemble of PDF replicas of the three.

| 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
$x$ | 3.30$\times 10^{-5}$ | 3.73$\times 10^{-4}$ | 3.50$\times 10^{-3}$ | 1.24$\times 10^{-2}$ | 2.48$\times 10^{-2}$ | 4.34$\times 10^{-2}$ | 8.59$\times 10^{-2}$
| 8 | 9 | 10 | 11 | 12 | 13 | 14
$x$ | 0.118 | 0.167 | 0.206 | 0.302 | 0.406 | 0.637 | 0.831

Table 4: The $x$ values we choose for the training and test samples of NNs.

The 8000 PDF replicas used for training and test are generated through a
randomly sampling of the PDF parameters defined in Eq. (4) with the help of
CT18 NNLO Hessian PDF set. Each replica or PDF parameters ${a^{rep}_{i}}$ is
determined by 28 randomly distributed variables ${r_{j}}$, namely

$a_{i}^{rep}=a_{i}^{0}+\sum_{j=1}^{28}r_{j}(a_{i}^{2j-1}-a_{i}^{2j})/2,$ (18)

where {$a_{i}^{0}$} represent the $i_{th}$ PDF parameter of the central PDF of
CT18 NNLO, {$a_{i}^{2j-1}$} and {$a_{i}^{2j}$} represent the $i_{th}$ PDF
parameter of the error PDFs in the plus and the minus direction of the
$i_{th}$ eigenvector respectively. For each $r_{j}$ we use a Gaussian sampling
with mean value 0 and variance $1/\sqrt{28}$ which ensures coverage of the PDF
parameter space with average increase of global $\chi^{2}$ of a few hundred
units comparing to CT18 best fit as shown in Fig. 2. We note that performances
of the trained NNs are not sensitive to the choice of training samples as far
as we are within or close to the uncertainty range of CT18.

We further test performance of our NN approach with alternative PDF
parametrization forms taking the target function of $\chi^{2}$ of the ATLAS 7
TeV jet data as an example. We have chosen MMHT2014 Harland-Lang:2014zoa ,
NNPDF3.1 NNPDF:2017mvq  and NNPDF4.0 Ball:2021leu  NNLO PDFs with 4000 MC PDF
replicas each generated from the corresponding Hessian PDF sets with LHAPDF6
Buckley:2014ana  222We have not used the native MC replicas of NNPDF since the
numbers of replicas are limited to be 1000 in that case.. We use the same
architecture as used for CT18 except for extensions to include 9 PDF flavors,
namely with $\bar{s}$, $c$, and $b$-quark PDFs in addition. The NNs have been
trained and tested for each individual parametrization form with the
corresponding MC replicas. We find very good performance of the NNs in cases
of MMHT2014 and NNPDF4.0 as shown in Fig. 24, similar to the case of CT18.
Interestingly, we find performance of the same architecture is much better for
the parametrization form of NNPDF4.0 than NNPDF3.1, possibly due to the smooth
conditions applied in NNPDF4.0 Ball:2021leu . We also try to train the NNs
with an ensemble of PDF replicas, 12000 replicas in total, from CT18, MMHT14
and NNPDF4.0. The accuracy of the trained NNs is only marginally worse than
the NNs trained to individual parametrization forms. That hints the
possibility of a universal NN to accommodate for a variety of smooth PDF
parametrizations, as discussed at earlier this section.

Finally we summarize the performances of our NNs in terms of computational
cost in Table 5 comparing to the traditional approaches. Note that we have not
included the time cost for the process of training of the NNs since we do not
need to repeat it in later scans of the PDF parameters. In Table 5 the numbers
indicate the time cost on a single CPU-core (2.4 GHz) of calculating the
target functions for a single point in the PDF parameter space. For $\chi^{2}$
the cost includes those for the calculations of the needed cross sections
(taking 10 points per data set as an example) and for the multiplications with
covariance matrix. In the conventional approach the computing efficiency
varies significantly, e.g., for the $\chi^{2}$, depending on the number of
data points, the perturbative order of the theory calculations, and
importantly whether the fast interpolation algorithms are used or not. Thus
included numbers only represent typical average cost in the CT18 NNLO analysis
for a direct calculation or using fast interpolations (shown in parenthesis).
The fast interpolation method for calculations of a single cross section
involves more PDF values on a dense grid and thus is slower than the NN
approaches. Nevertheless, the NN approaches lead to significant improvement in
general and ensure efficient scans of the PDF parameter space with much less
cost. The NNs was programmed with PYTHON2.7, and we expect further reduction
of the computational cost if transferred into more efficient programming
languages like Fortran/C++. We are planning to provide open source access for
the NN framework used together with trained NNs for various target functions
of CT18 in the near future.

method target cost | $\chi^{2}$ | $\sigma$ | $f(x,Q)$
---|---|---|---
NNs | 0.70 ms | 0.41 ms | 0.37 ms
traditional |  $10^{7}$(200) ms |  $10^{6}$(20) ms | 20(2) ms
Table 5: Comparison on computational cost between NNs and traditional methods.
Numbers in parenthesis represent cases if fast interpolations on PDFs are
used.

##  Appendix B Hessian PDF set

We further generate a Hessian PDF set based on the $\chi^{2}$ profile obtained
with the NN approaches. The Hessian error matrix on the PDF parameters is
calculated using a numeric method of finite difference. We use an iterative
algorithm on diagonalization of the Hessian matrix that is developed in Ref.
Pumplin:2000vx ; Pumplin:2001ct  and used in later CTEQ analyses. The
iterative procedure greatly improves the performance of Hessian approximation
in the case of large number of free parameters (28 here) and in the existence
of flat directions. Once all orthogonal eigenvectors are determined, two error
PDFs are generated for each eigenvector by scanning along the plus and the
minus directions and looking for solutions with $\Delta\chi^{2}+P=100$ (for
90% CL).

((a))

((b))

((c))

((d))

((e))

((f))

Figure 25: The parton distribution functions at $Q=1.295$ GeV for $u$, $d$,
$\bar{u}$, $\bar{d}$, $g$ and $s$. The blue and the red solid lines represent
the central values of CT18NNLO and NNs respectively. The blue and the red
hatched areas represent the uncertainties of CT18NNLO and final PDFs in this
paper at 68% CL respectively.

We compare the PDF uncertainties at 68% CL from the Hessian PDF set to the
published CT18 NNLO PDFs in Fig. 25. We find very good agreements between
predictions of the two Hessian PDF sets in general. However, some notable
differences can be seen for $d$-valence and gluon PDFs at large-$x$ ($\sim
0.4$) as well as for sea quarks at $x\lesssim 10^{-3}$. There are two reasons
that lead to the differences in the new Hessian set and the CT18 set. First as
mentioned earlier in the global fit presented in this paper the NNLO K-factors
used for predictions of the Drell-Yan data have been updated comparing to
those used in the CT18 analysis. Besides, when calculating the Hessian error
matrix numerically we use a step size of $\Delta\chi^{2}=10$ on sampling of
the PDF parameters while a value of $\sim 1$ is used in the CT18 analysis. The
dependence on choices of this step size reflect one intrinsic uncertainty of
the Hessian approaches Pumplin:2000vx .

##  Appendix C Variant fits with NOMAD data

In this appendix we present a series of global fit with inclusion of the NOMAD
data and with different theories or different choices of decay branching
ratios of charm quark to muon. In Sect. 5.1 when comparing to the global fit
without NOMAD data, we use NLO cross sections (but with NNLO PDFs) for NuTeV
and CCFR dimuon data to be consistent with the CT18 analysis. Thus the changes
observed after including the NOMAD data can be due to both the NOMAD data or
the changes of theories for the other two dimuon data. Now we further consider
different choices of the theory predictions, namely either calculated at NLO
or NNLO in QCD, to disentangle their effects.

Furthermore, to compare with dimuon data, one has to convert the cross
sections of charm-quark production to production of dimuon which relies on the
input of inclusive semileptonic branching ratio Br($c\to\mu$). Since NOMAD
dimuon data extend down to $E_{\nu}\sim$ 6 GeV the energy dependence of
Br($c\to\mu$) is taken into account in NOMAD analysis, with a parametrization
form

$Br(c\to\mu)=\dfrac{a}{1+b/E_{\nu}},$ (19)

where $a$ and $b$ are free parameters. In the NOMAD paper it suggests values
of $a=0.094\pm 0.010$ and $b=6.6\pm 3.9$ GeV as measured by the E531
experiment FermilabE531:1988hzw . The uncertainties on parameters $a$ and $b$
will propagate into the unfolded charm-quark cross sections and are treated as
additional correlated systematic errors that are summarized in Table 6 for
distribution in Bjorken-$x$. Other correlated systematic uncertainties for
NOMAD data can be found in Ref. NOMAD:2013hbk . On the other hand, for NuTeV
and CCFR data, since the neutrino energies are sufficiently high, a constant
value of Br($c\to\mu$)=$0.099\pm 0.010$ has been suggested osti_879078  and is
used in the CT18 analysis. The central value is slightly higher than the
parameter $a$ used in our nominal fit of NOMAD data. Thus we perform variant
fits using $a=0.099\pm 0.010$ for NOMAD data to further investigate the impact
of this overall normalization on the outcome PDFs.

In Fig. 26 we compare the strange-quark PDFs at 1.295 GeV from all variant
fits. We show the PDF uncertainties at 68% CL from LM scans for fits with and
without NOMAD data and using NNLO predictions from dimuon production
consistently. That can be compared with Fig. 17 where NLO predictions are used
in fit without NOMAD data. We also present central PDFs obtained with NLO
predictions or with higher branching ratio for NOMAD data. We find including
the NNLO corrections leads to a moderate increase of the strange-quark PDF,
which is in consistent with the conclusions in Ref. Hou:2019efy . The
inclusion of NOMAD data results in about 20% enhancement of the strange-quark
PDF at $x$ around 0.05 and a significant reduction of the PDF uncertainties.
Changing to $a=0.099\pm 0.010$ for NOMAD only induces a minor reduction of the
strange-quark PDF.

We further summarize the total or individual $\chi^{2}$ of all variant fits
together with predictions on $R_{s}(x=0.023,Q=1.5{\rm GeV})$ with
uncertainties at 68% CL in Table 7. By comparison with the $\chi^{2}$ we find
the NNLO predictions in general lead to a slightly worse fit with increase on
$\chi^{2}$ of a few units. However, in all cases the global fit can describe
well various dimuon data as can be seen from the $\chi^{2}$ per number of
degree of freedoms. When comparing fits in the last two rows we find using a
consistent branching ratio in different data sets results in a better fit and
reduced PDF uncertainties.

$x_{\rm Bj}$ | Bin center |  $\sigma_{\mu\mu}/\sigma_{cc}\pm\delta^{stat}\pm\delta^{syst}$ ($10^{-3}$) |  $\delta^{a}$, % |  $\delta^{b}$, %
---|---|---|---|---
0.0000 - 0.0336 | 0.0267 |  13.383 $\pm$ 0.441 $\pm$ 0.289 | 10.6 | 5.3
0.0336 - 0.0511 | 0.0440 |  11.245 $\pm$ 0.380 $\pm$ 0.210 | 10.6 | 6.8
0.0511 - 0.0672 | 0.0598 |  9.991 $\pm$ 0.347 $\pm$ 0.201 | 10.6 | 7.7
0.0672 - 0.0836 | 0.0756 |  9.141 $\pm$ 0.324 $\pm$ 0.189 | 10.6 | 8.3
0.0836 - 0.1000 | 0.0917 |  8.198 $\pm$ 0.297 $\pm$ 0.169 | 10.6 | 8.8
0.1000 - 0.1246 | 0.1122 |  7.176 $\pm$ 0.225 $\pm$ 0.144 | 10.6 | 9.0
0.1246 - 0.1535 | 0.1389 |  6.229 $\pm$ 0.195 $\pm$ 0.118 | 10.6 | 9.4
0.1535 - 0.1870 | 0.1699 |  5.427 $\pm$ 0.171 $\pm$ 0.106 | 10.6 | 9.6
0.1870 - 0.2277 | 0.2066 |  4.837 $\pm$ 0.151 $\pm$ 0.093 | 10.6 | 9.9
0.2277 - 0.2800 | 0.2524 |  4.235 $\pm$ 0.133 $\pm$ 0.083 | 10.6 | 10.0
0.2800 - 0.3590 | 0.3165 |  3.595 $\pm$ 0.113 $\pm$ 0.072 | 10.6 | 10.0
0.3590 - 0.4583 | 0.4036 |  2.955 $\pm$ 0.111 $\pm$ 0.062 | 10.6 | 10.1
0.4583 - 0.5838 | 0.5116 |  2.355 $\pm$ 0.120 $\pm$ 0.055 | 10.6 | 9.9
0.5838 - 0.7500 | 0.6465 |  1.607 $\pm$ 0.150 $\pm$ 0.047 | 10.6 | 9.4
Table 6: NOMAD measurements on the Bjorken-$x$ distribution of dimuon to inclusive CC cross section ratio, including the binning, central values, statistical and total systematic uncertainties. The last two columns show the additional correlated systematic uncertainties in percentages, if converting back to production cross sections of charm-quark, due to input parameter $a$ and $b$ respectively. These additional errors are derived based on theoretical cross sections at NLO with CT18 NNLO PDFs.  Figure 26: Strange-quark PDF at $Q=1.295$ GeV from LM scans in global fits with various conditions. The PDF uncertainties are shown for 68% CL. data sets |  $\chi^{2}_{total}$ (3671/3683) |  $\chi^{2}_{nomad}$ (12) |  $\chi^{2}_{124}$ (38) |  $\chi^{2}_{125}$ (33) |  $\chi^{2}_{126}$ (40) |  $\chi^{2}_{127}$ (38) |  $R_{s}$(0.023,1.5GeV)
---|---|---|---|---|---|---|---
NUT (NLO) | 4272.64 | - | 18.83 | 39.55 | 29.89 | 19.42 | $\mathrm{0.518^{+0.349}_{-0.363}}$
NUT (NNLO) | 4268.77 | - | 21.44 | 32.84 | 34.06 | 22.54 | $\mathrm{0.616^{+0.441}_{-0.377}}$
NUT (NLO) + NOM (NLO) | 4286.28 | 8.39 | 24.81 | 41.95 | 29.30 | 18.79 | $\mathrm{0.593^{+0.256}_{-0.155}}$
NUT (NNLO) + NOM (NNLO) | 4291.47 | 14.47 | 28.13 | 34.26 | 34.21 | 22.15 | $\mathrm{0.695^{+0.384}_{-0.169}}$
B r(c$\to\mu$)=0.099 | 4289.61 | 13.79 | 27.38 | 34.06 | 34.10 | 22.13 | $\mathrm{0.685^{+0.290}_{-0.174}}$
Table 7: Total $\chi^{2}$ and individual $\chi^{2}$ of dimuon data for global
fits with various conditions. Numbers in the first row indicate the total
number of data points. The last column includes predictions on
$R_{s}(x=0.023,Q=1.5{\rm GeV})$ with uncertainties at 90% CL.

## References

  * (1) J. Gao, L. Harland-Lang, and J. Rojo, The Structure of the Proton in the LHC Precision Era, Phys. Rept. 742 (2018) 1–121, [arXiv:1709.04922].
  * (2) K. Kovařík, P. M. Nadolsky, and D. E. Soper, Hadronic structure in high-energy collisions, Rev. Mod. Phys. 92 (2020), no. 4 045003, [arXiv:1905.06957].
  * (3) J. C. Collins, D. E. Soper, and G. F. Sterman, Factorization of Hard Processes in QCD, Adv. Ser. Direct. High Energy Phys. 5 (1989) 1–91, [hep-ph/0409313].
  * (4) X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao, Large-momentum effective theory, Rev. Mod. Phys. 93 (2021), no. 3 035005, [arXiv:2004.03543].
  * (5) Lattice Parton Collaboration, Y.-S. Liu et al., Unpolarized isovector quark distribution function from lattice QCD: A systematic analysis of renormalization and matching, Phys. Rev. D 101 (2020), no. 3 034020, [arXiv:1807.06566].
  * (6) LHC Higgs Cross Section Working Group Collaboration, D. de Florian et al., Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector, arXiv:1610.07922.
  * (7) W. Beenakker, C. Borschensky, M. Krämer, A. Kulesza, E. Laenen, S. Marzani, and J. Rojo, NLO+NLL squark and gluino production cross-sections with threshold-improved parton distributions, Eur. Phys. J. C 76 (2016), no. 2 53, [arXiv:1510.00375].
  * (8) S. Alioli, M. Farina, D. Pappadopulo, and J. T. Ruderman, Precision Probes of QCD at High Energies, JHEP 07 (2017) 097, [arXiv:1706.03068].
  * (9) ATLAS Collaboration, Determination of the strong coupling constant and test of asymptotic freedom from Transverse Energy-Energy Correlations in multijet events at $\sqrt{s}=13$ TeV with the ATLAS detector, .
  * (10) CMS, ATLAS Collaboration, S. Kaur, Standard Model: Electroweak Physics with CMS and ATLAS at 13 TeV, J. Phys. Conf. Ser. 1258 (2019), no. 1 012015.
  * (11) ATLAS Collaboration, M. Aaboud et al., Measurement of the $W$-boson mass in pp collisions at $\sqrt{s}=7$ TeV with the ATLAS detector, Eur. Phys. J. C 78 (2018), no. 2 110, [arXiv:1701.07240]. [Erratum: Eur.Phys.J.C 78, 898 (2018)].
  * (12) T.-J. Hou et al., New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC, Phys. Rev. D 103 (2021), no. 1 014013, [arXiv:1912.10053].
  * (13) S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin, and R. S. Thorne, Parton distributions from LHC, HERA, Tevatron and fixed target data: MSHT20 PDFs, Eur. Phys. J. C 81 (2021), no. 4 341, [arXiv:2012.04684].
  * (14) R. D. Ball et al., The path to proton structure at 1% accuracy: NNPDF Collaboration, Eur. Phys. J. C 82 (2022), no. 5 428, [arXiv:2109.02653].
  * (15) S. Alekhin, J. Blümlein, and S. Moch, NLO PDFs from the ABMP16 fit, Eur. Phys. J. C 78 (2018), no. 6 477, [arXiv:1803.07537].
  * (16) H1 Collaboration, V. Andreev et al., Determination of the strong coupling constant $\alpha_{s}(m_{Z})$ in next-to-next-to-leading order QCD using H1 jet cross section measurements, Eur. Phys. J. C 77 (2017), no. 11 791, [arXiv:1709.07251].
  * (17) P. Jimenez-Delgado and E. Reya, Delineating parton distributions and the strong coupling, Phys. Rev. D 89 (2014), no. 7 074049, [arXiv:1403.1852].
  * (18) S. Park, A. Accardi, X. Jing, and J. F. Owens, CJ15 global PDF analysis with new electroweak data from the STAR and SeaQuest experiments, in 28th International Workshop on Deep Inelastic Scattering and Related Subjects, 8, 2021.  arXiv:2108.05786.
  * (19) ATLAS Collaboration, G. Aad et al., Determination of the parton distribution functions of the proton using diverse ATLAS data from $pp$ collisions at $\sqrt{s}=7$, 8 and 13 TeV, arXiv:2112.11266.
  * (20) J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai, and W. K. Tung, Uncertainties of predictions from parton distribution functions. 2. The Hessian method, Phys. Rev. D 65 (2001) 014013, [hep-ph/0101032].
  * (21) A. D. Martin, W. J. Stirling, R. S. Thorne, and G. Watt, Parton distributions for the LHC, Eur. Phys. J. C 63 (2009) 189–285, [arXiv:0901.0002].
  * (22) S. Forte, L. Garrido, J. I. Latorre, and A. Piccione, Neural network parametrization of deep inelastic structure functions, JHEP 05 (2002) 062, [hep-ph/0204232].
  * (23) J. Pumplin, D. R. Stump, and W. K. Tung, Multivariate fitting and the error matrix in global analysis of data, Phys. Rev. D 65 (2001) 014011, [hep-ph/0008191].
  * (24) D. Stump, J. Pumplin, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai, and W. K. Tung, Uncertainties of predictions from parton distribution functions. 1. The Lagrange multiplier method, Phys. Rev. D 65 (2001) 014012, [hep-ph/0101051].
  * (25) J. Gao and P. Nadolsky, A meta-analysis of parton distribution functions, JHEP 07 (2014) 035, [arXiv:1401.0013].
  * (26) C. Schmidt, J. Pumplin, C. P. Yuan, and P. Yuan, Updating and optimizing error parton distribution function sets in the Hessian approach, Phys. Rev. D 98 (2018), no. 9 094005, [arXiv:1806.07950].
  * (27) B.-T. Wang, T. J. Hobbs, S. Doyle, J. Gao, T.-J. Hou, P. M. Nadolsky, and F. I. Olness, Mapping the sensitivity of hadronic experiments to nucleon structure, Phys. Rev. D 98 (2018), no. 9 094030, [arXiv:1803.02777].
  * (28) D. Guest, K. Cranmer, and D. Whiteson, Deep Learning and its Application to LHC Physics, Ann. Rev. Nucl. Part. Sci. 68 (2018) 161–181, [arXiv:1806.11484].
  * (29) NNPDF Collaboration, R. D. Ball et al., Parton distributions for the LHC Run II, JHEP 04 (2015) 040, [arXiv:1410.8849].
  * (30) S. Forte and S. Carrazza, Parton distribution functions, arXiv:2008.12305.
  * (31) NOMAD Collaboration, O. Samoylov et al., A Precision Measurement of Charm Dimuon Production in Neutrino Interactions from the NOMAD Experiment, Nucl. Phys. B 876 (2013) 339–375, [arXiv:1308.4750].
  * (32) A. Dainese, M. Mangano, A. B. Meyer, A. Nisati, G. Salam, and M. A. Vesterinen, eds., Report on the Physics at the HL-LHC,and Perspectives for the HE-LHC, vol. 7/2019 of CERN Yellow Reports: Monographs.  CERN, Geneva, Switzerland, 2019.
  * (33) R. Abdul Khalek, S. Bailey, J. Gao, L. Harland-Lang, and J. Rojo, Towards Ultimate Parton Distributions at the High-Luminosity LHC, Eur. Phys. J. C 78 (2018), no. 11 962, [arXiv:1810.03639].
  * (34) R. Abdul Khalek, S. Bailey, J. Gao, L. Harland-Lang, and J. Rojo, Probing Proton Structure at the Large Hadron electron Collider, SciPost Phys. 7 (2019), no. 4 051, [arXiv:1906.10127].
  * (35) S. Carrazza, C. Degrande, S. Iranipour, J. Rojo, and M. Ubiali, Can New Physics hide inside the proton?, Phys. Rev. Lett. 123 (2019), no. 13 132001, [arXiv:1905.05215].
  * (36) A. Greljo, S. Iranipour, Z. Kassabov, M. Madigan, J. Moore, J. Rojo, M. Ubiali, and C. Voisey, Parton distributions in the SMEFT from high-energy Drell-Yan tails, JHEP 07 (2021) 122, [arXiv:2104.02723].
  * (37) M. Madigan and J. Moore, Parton distributions in the SMEFT from high-energy Drell-Yan tails, in European Physical Society Conference on High Energy Physics 2021, 10, 2021.  arXiv:2110.13204.
  * (38) CMS Collaboration, A. Tumasyan et al., Measurement and QCD analysis of double-differential inclusive jet cross sections in proton-proton collisions at $\sqrt{s}$ = 13 TeV, arXiv:2111.10431.
  * (39) S. Iranipour and M. Ubiali, A new generation of simultaneous fits to LHC data using deep learning, JHEP 05 (2022) 032, [arXiv:2201.07240].
  * (40) “Keras.” https://keras.io/.
  * (41) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP 07 (2014) 079, [arXiv:1405.0301].
  * (42) V. Bertone, R. Frederix, S. Frixione, J. Rojo, and M. Sutton, aMCfast: automation of fast NLO computations for PDF fits, JHEP 08 (2014) 166, [arXiv:1406.7693].
  * (43) T. Carli, D. Clements, A. Cooper-Sarkar, C. Gwenlan, G. P. Salam, F. Siegert, P. Starovoitov, and M. Sutton, A posteriori inclusion of parton density functions in NLO QCD final-state calculations at hadron colliders: The APPLGRID Project, Eur. Phys. J. C 66 (2010) 503–524, [arXiv:0911.2985].
  * (44) G. P. Salam and J. Rojo, A Higher Order Perturbative Parton Evolution Toolkit (HOPPET), Comput. Phys. Commun. 180 (2009) 120–156, [arXiv:0804.3755].
  * (45) H1, ZEUS Collaboration, H. Abramowicz et al., Combination of measurements of inclusive deep inelastic ${e^{\pm}p}$ scattering cross sections and QCD analysis of HERA data, Eur. Phys. J. C 75 (2015), no. 12 580, [arXiv:1506.06042].
  * (46) BCDMS Collaboration, A. C. Benvenuti et al., A High Statistics Measurement of the Proton Structure Functions F(2) (x, Q**2) and R from Deep Inelastic Muon Scattering at High Q**2, Phys. Lett. B 223 (1989) 485–489.
  * (47) BCDMS Collaboration, A. C. Benvenuti et al., A High Statistics Measurement of the Deuteron Structure Functions F2 (X, $Q^{2}$) and R From Deep Inelastic Muon Scattering at High $Q^{2}$, Phys. Lett. B 237 (1990) 592–598.
  * (48) New Muon Collaboration, M. Arneodo et al., Measurement of the proton and deuteron structure functions, F2(p) and F2(d), and of the ratio sigma-L / sigma-T, Nucl. Phys. B 483 (1997) 3–43, [hep-ph/9610231].
  * (49) J. P. Berge et al., A Measurement of Differential Cross-Sections and Nucleon Structure Functions in Charged Current Neutrino Interactions on Iron, Z. Phys. C 49 (1991) 187–224.
  * (50) CCFR/NuTeV Collaboration, U.-K. Yang et al., Measurements of $F_{2}$ and $xF^{\nu}_{3}-xF^{\bar{\nu}}_{3}$ from CCFR $\nu_{\mu}-$Fe and $\bar{\nu}_{\mu}-$Fe data in a physics model independent way, Phys. Rev. Lett. 86 (2001) 2742–2745, [hep-ex/0009041].
  * (51) W. G. Seligman et al., Improved determination of alpha(s) from neutrino nucleon scattering, Phys. Rev. Lett. 79 (1997) 1213–1216, [hep-ex/9701017].
  * (52) D. A. Mason, Measurement of the strange - antistrange asymmetry at nlo in qcd from nutev dimuon data, Physics of Elementary Particles and Fields (2006).
  * (53) NuTeV Collaboration, M. Goncharov et al., Precise Measurement of Dimuon Production Cross-Sections in $\nu_{\mu}$ Fe and $\bar{\nu}_{\mu}$ Fe Deep Inelastic Scattering at the Tevatron., Phys. Rev. D 64 (2001) 112006, [hep-ex/0102049].
  * (54) H1 Collaboration, A. Aktas et al., Measurement of F2($c\bar{c}$) and F2($b\bar{b}$) at high $Q^{2}$ using the H1 vertex detector at HERA, Eur. Phys. J. C 40 (2005) 349–359, [hep-ex/0411046].
  * (55) H1, ZEUS Collaboration, H. Abramowicz et al., Combination and QCD Analysis of Charm Production Cross Section Measurements in Deep-Inelastic ep Scattering at HERA, Eur. Phys. J. C 73 (2013), no. 2 2311, [arXiv:1211.1182].
  * (56) H1 Collaboration, F. D. Aaron et al., Measurement of the Inclusive e\pmp Scattering Cross Section at High Inelasticity y and of the Structure Function $F_{L}$, Eur. Phys. J. C 71 (2011) 1579, [arXiv:1012.4355].
  * (57) G. Moreno et al., Dimuon production in proton - copper collisions at $\sqrt{s}$ = 38.8-GeV, Phys. Rev. D 43 (1991) 2815–2836.
  * (58) NuSea Collaboration, R. S. Towell et al., Improved measurement of the anti-d / anti-u asymmetry in the nucleon sea, Phys. Rev. D 64 (2001) 052002, [hep-ex/0103030].
  * (59) NuSea Collaboration, J. C. Webb et al., Absolute Drell-Yan dimuon cross-sections in 800 GeV / c pp and pd collisions, hep-ex/0302019.
  * (60) CDF Collaboration, F. Abe et al., Measurement of the Lepton Charge Asymmetry in $W$ Boson Decays Produced in $p\bar{p}$ Collisions, Phys. Rev. Lett. 81 (1998) 5754–5759, [hep-ex/9809001].
  * (61) CDF Collaboration, D. Acosta et al., Measurement of the forward-backward charge asymmetry from $W\to e\nu$ production in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV, Phys. Rev. D 71 (2005) 051104, [hep-ex/0501023].
  * (62) D0 Collaboration, V. M. Abazov et al., Measurement of the muon charge asymmetry from $W$ boson decays, Phys. Rev. D 77 (2008) 011106, [arXiv:0709.4254].
  * (63) D0 Collaboration, V. M. Abazov et al., Measurement of the Shape of the Boson Rapidity Distribution for $p\bar{p}\to Z/\gamma^{*}\to e^{+}e^{-}+X$ Events Produced at $\sqrt{s}$ of 1.96-TeV, Phys. Rev. D 76 (2007) 012003, [hep-ex/0702025].
  * (64) CDF Collaboration, T. A. Aaltonen et al., Measurement of $d\sigma/dy$ of Drell-Yan $e^{+}e^{-}$ pairs in the $Z$ Mass Region from $p\bar{p}$ Collisions at $\sqrt{s}=1.96$ TeV, Phys. Lett. B 692 (2010) 232–239, [arXiv:0908.3914].
  * (65) CMS Collaboration, S. Chatrchyan et al., Measurement of the Muon Charge Asymmetry in Inclusive $pp\to W+X$ Production at $\sqrt{s}=$ 7 TeV and an Improved Determination of Light Parton Distribution Functions, Phys. Rev. D 90 (2014), no. 3 032004, [arXiv:1312.6283].
  * (66) CMS Collaboration, S. Chatrchyan et al., Measurement of the Electron Charge Asymmetry in Inclusive $W$ Production in $pp$ Collisions at $\sqrt{s}=7$ TeV, Phys. Rev. Lett. 109 (2012) 111806, [arXiv:1206.2598].
  * (67) ATLAS Collaboration, G. Aad et al., Measurement of the inclusive $W^{\pm}$ and Z/gamma cross sections in the electron and muon decay channels in $pp$ collisions at $\sqrt{s}=7$ TeV with the ATLAS detector, Phys. Rev. D 85 (2012) 072004, [arXiv:1109.5141].
  * (68) D0 Collaboration, V. M. Abazov et al., Measurement of the electron charge asymmetry in $\boldsymbol{p\bar{p}\rightarrow W+X\rightarrow e\nu+X}$ decays in $\boldsymbol{p\bar{p}}$ collisions at $\boldsymbol{\sqrt{s}=1.96}$ TeV, Phys. Rev. D 91 (2015), no. 3 032007, [arXiv:1412.2862]. [Erratum: Phys.Rev.D 91, 079901 (2015)].
  * (69) CDF Collaboration, T. Aaltonen et al., Measurement of the Inclusive Jet Cross Section at the Fermilab Tevatron p anti-p Collider Using a Cone-Based Jet Algorithm, Phys. Rev. D 78 (2008) 052006, [arXiv:0807.2204]. [Erratum: Phys.Rev.D 79, 119902 (2009)].
  * (70) D0 Collaboration, V. M. Abazov et al., Measurement of the inclusive jet cross-section in $p\bar{p}$ collisions at $s^{(1/2)}$ =1.96-TeV, Phys. Rev. Lett. 101 (2008) 062001, [arXiv:0802.2400].
  * (71) LHCb Collaboration, R. Aaij et al., Measurement of the forward $Z$ boson production cross-section in $pp$ collisions at $\sqrt{s}=7$ TeV, JHEP 08 (2015) 039, [arXiv:1505.07024].
  * (72) LHCb Collaboration, R. Aaij et al., Measurement of forward $\rm Z\rightarrow e^{+}e^{-}$ production at $\sqrt{s}=8$ TeV, JHEP 05 (2015) 109, [arXiv:1503.00963].
  * (73) CMS Collaboration, V. Khachatryan et al., Measurement of the differential cross section and charge asymmetry for inclusive $\mathrm{p}\mathrm{p}\rightarrow\mathrm{W}^{\pm}+X$ production at ${\sqrt{s}}=8$ TeV, Eur. Phys. J. C 76 (2016), no. 8 469, [arXiv:1603.01803].
  * (74) LHCb Collaboration, R. Aaij et al., Measurement of forward W and Z boson production in $pp$ collisions at $\sqrt{s}=8$ TeV, JHEP 01 (2016) 155, [arXiv:1511.08039].
  * (75) ATLAS Collaboration, G. Aad et al., Measurement of the transverse momentum and $\phi^{*}_{\eta}$ distributions of Drell–Yan lepton pairs in proton–proton collisions at $\sqrt{s}=8$ TeV with the ATLAS detector, Eur. Phys. J. C 76 (2016), no. 5 291, [arXiv:1512.02192].
  * (76) CMS Collaboration, S. Chatrchyan et al., Measurement of the Ratio of Inclusive Jet Cross Sections using the Anti-$k_{T}$ Algorithm with Radius Parameters R=0.5 and 0.7 in pp Collisions at $\sqrt{s}=7$ TeV, Phys. Rev. D 90 (2014), no. 7 072006, [arXiv:1406.0324].
  * (77) ATLAS Collaboration, G. Aad et al., Measurement of the inclusive jet cross-section in proton-proton collisions at $\sqrt{s}=7$ TeV using 4.5 $fb^{-1}$ of data with the ATLAS detector, JHEP 02 (2015) 153, [arXiv:1410.8857]. [Erratum: JHEP 09, 141 (2015)].
  * (78) CMS Collaboration, V. Khachatryan et al., Measurement and QCD analysis of double-differential inclusive jet cross sections in pp collisions at $\sqrt{s}=8$ TeV and cross section ratios to 2.76 and 7 TeV, JHEP 03 (2017) 156, [arXiv:1609.05331].
  * (79) CMS Collaboration, A. M. Sirunyan et al., Measurement of double-differential cross sections for top quark pair production in pp collisions at $\sqrt{s}=8$ TeV and impact on parton distribution functions, Eur. Phys. J. C 77 (2017), no. 7 459, [arXiv:1703.01630].
  * (80) ATLAS Collaboration, G. Aad et al., Measurements of top-quark pair differential cross-sections in the lepton+jets channel in $pp$ collisions at $\sqrt{s}=8$ TeV using the ATLAS detector, Eur. Phys. J. C 76 (2016), no. 10 538, [arXiv:1511.04716].
  * (81) S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump, and C. P. Yuan, New parton distribution functions from a global analysis of quantum chromodynamics, Phys. Rev. D 93 (2016), no. 3 033006, [arXiv:1506.07443].
  * (82) H.-L. Lai, M. Guzzi, J. Huston, Z. Li, P. M. Nadolsky, J. Pumplin, and C. P. Yuan, New parton distributions for collider physics, Phys. Rev. D 82 (2010) 074024, [arXiv:1007.2241].
  * (83) C. Anastasiou, C. Duhr, F. Dulat, F. Herzog, and B. Mistlberger, Higgs Boson Gluon-Fusion Production in QCD at Three Loops, Phys. Rev. Lett. 114 (2015) 212001, [arXiv:1503.06056].
  * (84) E. L. Berger, J. Gao, C. S. Li, Z. L. Liu, and H. X. Zhu, Charm-Quark Production in Deep-Inelastic Neutrino Scattering at Next-to-Next-to-Leading Order in QCD, Phys. Rev. Lett. 116 (2016), no. 21 212002, [arXiv:1601.05430].
  * (85) J. Gao, Massive charged-current coefficient functions in deep-inelastic scattering at NNLO and impact on strange-quark distributions, JHEP 02 (2018) 026, [arXiv:1710.04258].
  * (86) J. Gao, T. J. Hobbs, P. M. Nadolsky, C. Sun, and C. P. Yuan, A general heavy-flavor mass scheme for charged-current DIS at NNLO and beyond, arXiv:2107.00460.
  * (87) S. Alekhin, J. Blumlein, L. Caminadac, K. Lipka, K. Lohwasser, S. Moch, R. Petti, and R. Placakyte, Determination of Strange Sea Quark Distributions from Fixed-target and Collider Data, Phys. Rev. D 91 (2015), no. 9 094002, [arXiv:1404.6469].
  * (88) F. Faura, S. Iranipour, E. R. Nocera, J. Rojo, and M. Ubiali, The Strangest Proton?, Eur. Phys. J. C 80 (2020), no. 12 1168, [arXiv:2009.00014].
  * (89) ATLAS Collaboration, G. Aad et al., Determination of the strange quark density of the proton from ATLAS measurements of the $W\to\ell\nu$ and $Z\to\ell\ell$ cross sections, Phys. Rev. Lett. 109 (2012) 012001, [arXiv:1203.4051].
  * (90) ATLAS Collaboration, M. Aaboud et al., Precision measurement and interpretation of inclusive $W^{+}$ , $W^{-}$ and $Z/\gamma^{*}$ production cross sections with the ATLAS detector, Eur. Phys. J. C 77 (2017), no. 6 367, [arXiv:1612.03016].
  * (91) ATLAS Collaboration, G. Aad et al., Measurement of the double-differential high-mass Drell-Yan cross section in pp collisions at $\sqrt{s}=8$ TeV with the ATLAS detector, JHEP 08 (2016) 009, [arXiv:1606.01736].
  * (92) T. Kluge, K. Rabbertz, and M. Wobisch, FastNLO: Fast pQCD calculations for PDF fits, in 14th International Workshop on Deep Inelastic Scattering, pp. 483–486, 9, 2006.  hep-ph/0609285.
  * (93) L. A. Harland-Lang, A. D. Martin, P. Motylinski, and R. S. Thorne, Parton distributions in the LHC era: MMHT 2014 PDFs, Eur. Phys. J. C 75 (2015), no. 5 204, [arXiv:1412.3989].
  * (94) NNPDF Collaboration, R. D. Ball et al., Parton distributions from high-precision collider data, Eur. Phys. J. C 77 (2017), no. 10 663, [arXiv:1706.00428].
  * (95) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, and G. Watt, LHAPDF6: parton density access in the LHC precision era, Eur. Phys. J. C 75 (2015) 132, [arXiv:1412.7420].
  * (96) Fermilab E531 Collaboration, N. Ushida et al., Cross-sections for Neutrino Production of Charmed Particles, Phys. Lett. B 206 (1988) 375–379.