# The Strangest Proton?

Ferran Faura  Nikhef Theory Group, Science Park 105, 1098 XG Amsterdam, The
Netherlands  Shayan Iranipour  DAMTP, University of Cambridge, Wilberforce
Road, Cambridge, CB3 0WA, United Kingdom  Emanuele R. Nocera  Nikhef Theory
Group, Science Park 105, 1098 XG Amsterdam, The Netherlands  Juan Rojo  Nikhef
Theory Group, Science Park 105, 1098 XG Amsterdam, The Netherlands  Department
of Physics and Astronomy, VU Amsterdam, 1081 HV Amsterdam, The Netherlands
Maria Ubiali  DAMTP, University of Cambridge, Wilberforce Road, Cambridge, CB3
0WA, United Kingdom

###### Abstract

We present an improved determination of the strange quark and anti-quark
parton distribution functions of the proton by means of a global QCD analysis
that takes into account a comprehensive set of strangeness-sensitive
measurements: charm-tagged cross sections for fixed-target neutrino-nucleus
deep-inelastic scattering, and cross sections for inclusive gauge-boson
production and $W$-boson production in association with light jets or charm
quarks at hadron colliders. Our analysis is accurate to next-to-next-to-
leading order in perturbative QCD and specifically includes charm-quark mass
corrections to neutrino-nucleus structure functions. We find that a good
overall description of the input dataset can be achieved and that a
strangeness moderately suppressed in comparison to the rest of the light sea
quarks is strongly favored by the global analysis.

††preprint: Nikhef 2020-021

Introduction. An accurate determination of the strange quark and anti-quark
parton distribution functions (PDFs) of the proton Gao  _et al._ (2018);
Kovarik  _et al._ (2019); Ethier and Nocera (2020) is key to carrying out
precision phenomenology at current and future colliders, specifically for
measuring fundamental parameters of the Standard Model such as the mass of the
$W$ boson Aaboud  _et al._ (2018a) and the Weinberg angle Sirunyan  _et al._
(2018), see also Bagnaschi and Vicini (2019). Because of the limited
experimental information available, however, the strange quark and anti-quark
PDFs remain significantly more uncertain than the other light sea quark PDFs.

The strange quark and anti-quark PDFs have been determined from neutrino-
nucleus deep-inelastic scattering (DIS) for a long time, specifically from
measurements of dimuon cross sections, whereby the secondary muon originates
from the decay of a charmed meson, $\nu_{\mu}+N\to\mu+c+X$ with $c\to
D\to\mu+X$ Bazarko  _et al._ (1995); Kayis-Topaksu  _et al._ (2008); Mason
_et al._ (2007); Samoylov  _et al._ (2013). When interpreted in terms of the
ratio between strange and non-strange sea quark PDFs,
$R_{s}\equiv(s+\bar{s})/(\bar{u}+\bar{d})$, these measurements favor values
around $R_{s}\simeq 0.5$ when PDFs are evaluated at values of the momentum
fraction $x=0.023$ and scale $Q=1.6$ GeV. Therefore it came as a surprise when
a QCD analysis of the $W$\- and $Z$-boson rapidity distributions measured by
the ATLAS experiment in proton-proton collisions Aad  _et al._ (2012a)
suggested instead a ratio closer to $R_{s}\simeq 1$. This result was
subsequently corroborated by an analysis based on an increased integrated
luminosity Aaboud  _et al._ (2017). Complementary information on the strange
quark and anti-quark PDFs is provided by $W$-boson production in association
with light jets Sutton (2019) and charm quarks Stirling and Vryonidou (2012),
the latter process being dominated by the partonic scattering $g+s\to W+c$.
However, measurements of these processes performed by ATLAS Aaboud  _et al._
(2018b); Aad  _et al._ (2014) and CMS Chatrchyan  _et al._ (2014); Sirunyan
_et al._ (2019) point to opposite directions: the former to $R_{s}\simeq 1$;
the latter to $R_{s}\simeq 0.5$.

This state of affairs has been a source of controversy in recent accurate
determinations of PDFs. The NNPDF3.1 global analysis Ball  _et al._ (2017)
found that, whereas the ATLAS $W$, $Z$ dataset does favour a larger total
strangeness, its $\chi^{2}$ remains poor, and it has a moderate impact in the
global fit. Conversely, the recent CT18 global analysis Hou  _et al._ (2019)
presented fits with and without the ATLAS measurement of Aaboud  _et al._
(2017), with the resulting PDFs being significantly different. Dedicated
studies of the strange quark and anti-quark PDFs have been presented Lai  _et
al._ (2007); Kusina  _et al._ (2012); Alekhin  _et al._ (2015a, 2018); Cooper-
Sarkar and Wichmann (2018), however these often focus on a restricted set of
processes or datasets, or are based on theoretical and methodological
assumptions that can potentially bias the results.

A global reinterpretation of all of the strangeness-sensitive measurements
within an accurate theoretical and methodological framework appears to be
therefore compelling in order to establish whether or not a proton strangeness
crisis occurs. This paper fulfills this purpose: we present an improved
determination of the strange quark and anti-quark PDFs, accurate to next-to-
next-to-leading order (NNLO) in perturbative QCD, by expanding the NNPDF3.1
analysis Ball  _et al._ (2017) in two respects. First, we take into account an
extended piece of experimental information which is relevant to constrain the
strange quark and anti-quark PDFs. Second, we improve the theoretical
description of dimuon neutrino DIS data, by implementing NNLO charm-quark mass
corrections, and of $W$$+$$c$ production data, by including a theoretical
uncertainty that accounts for the unknown NNLO QCD corrections; we also
explicitly enforce the positivity of the $F_{2}^{c}$ structure function. See
Ball  _et al._ (2015) and references therein for a comprehensive description
of the NNPDF methodology.

Input dataset. The bulk of the dataset included in our analysis corresponds to
the one used in Ball  _et al._ (2018), which is in turn a variant of the
dataset used in the NNPDF3.1 NNLO analysis Ball  _et al._ (2017). It contains
in particular measurements of dimuon neutrino-nucleus DIS cross sections from
the NuTeV experiment Mason  _et al._ (2007), and of inclusive gauge boson
production in proton-(anti)proton collisions from several Tevatron and LHC
experiments Aaltonen  _et al._ (2010); Abazov  _et al._ (2007); Aad  _et al._
(2012b); Aaboud  _et al._ (2017); Khachatryan  _et al._ (2016). These
measurements represented the most constraining source of experimental
information on the strange quark and anti-quark PDFs in the NNPDF3.1 analysis.

We supplement this dataset with a number of new measurements. Concerning
neutrino-nucleus DIS, we include measurements of the ratio of dimuon to
inclusive charged-current cross sections,
$\mathcal{R}_{\mu\mu}(\omega)=\sigma_{\mu\mu}(\omega)/\sigma_{\rm
CC}(\omega)$, from the NOMAD experiment Samoylov  _et al._ (2013). The data is
presented for three kinematic variables $\omega$: the neutrino beam energy
$E_{\nu}$, the momentum fraction $x$, and the square root of the final-state
invariant mass $\sqrt{\hat{s}}$. Because experimental correlations are not
provided amongst measurements in different kinematic variables, only one
measurement can be included in the fit at a time: we select the $n_{\rm
dat}=19$ data points as a function of $E_{\nu}$, and verify that similar
results can be obtained with the $\sqrt{\hat{s}}$ set, see the Supplementary
Material (SM). The kinematic sensitivity of the NOMAD measurements is roughly
$0.03\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$<$}}x\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}0.7$.

Concerning proton-proton collisions, we augment the inclusive gauge boson
production measurement from the ATLAS experiment (7 TeV) Aaboud  _et al._
(2017) with the off-peak and forward rapidity bins (not included in NNPDF3.1)
for a total of $n_{\rm dat}=61$ data points. Furthermore, we include the
$n_{\rm dat}=37$ data points corresponding to the ATLAS (7 TeV) Aad  _et al._
(2014) and CMS (7 TeV and 13 TeV) Chatrchyan  _et al._ (2014); Sirunyan  _et
al._ (2019) $W$$+$$c$ measurements; for ATLAS, we consider the charm-jet
dataset, which is amenable to fixed-order calculations (instead of the
$D$-meson dataset). Finally, we take into account the $n_{\rm dat}=32$ data
points corresponding to the ATLAS $W$+jets measurement (8 TeV) differential in
the transverse momentum of the $W$ boson Aaboud  _et al._ (2018b). Overall,
these LHC processes are sensitive to the proton strangeness in the region
$10^{-3}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}x\mathrel{\hbox to0.0pt{\lower
4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}0.1$. Our analysis
contains a total of $n_{\rm dat}=4096$ data points; experimental correlations
within each dataset are available for all of the new measurements considered
here and are therefore included in our analysis.

Theoretical calculations. The theoretical settings adopted in our analysis
closely follow those described in Ball  _et al._ (2017, 2018) (whereby, in
particular, the charm PDF is fitted), with two improvements. First, the
positivity of the structure function $F_{2}^{c}$ is now enforced with a
procedure similar to that described in Ball  _et al._ (2015) for light quarks.
This additional constraint is required to prevent the fitted charm PDF
becoming unphysically negative once the new datasets are included in the fit.
Second, we incorporate the recently computed NNLO charm-quark mass corrections
Berger  _et al._ (2016); Gao (2018) in the description of the NuTeV and NOMAD
measurements. We do so by multiplying the next-to-leading order (NLO)
theoretical prediction in the FONLL general-mass scheme Forte  _et al._
(2010); Ball  _et al._ (2011a) by a $K$-factor defined as the ratio between
the NNLO result in the fixed-flavor number (FFN) scheme with and without the
charm-mass correction. This approach provides a good approximation of the
exact result, because theoretical predictions in the FFN and FONLL schemes are
very close for the NuTeV and NOMAD kinematics (see the SM). These $K$-factors
are in general smaller than unity, and thus enhance the (anti-)strange quark
PDF when accounted for in the fit. Nuclear corrections to neutrino-nucleus
data are not included: for NuTeV, they were demonstrated to be immaterial Ball
_et al._ (2019); for NOMAD, they mostly cancel out in the ratio
$\mathcal{R}_{\mu\mu}(\omega)$, see the SM.

Additional considerations apply specifically to the theoretical treatment of
the new datasets. The NOMAD observables are evaluated as two-dimensional
integrals over the differential DIS cross sections, e.g. for $E_{\nu}$,

$\sigma_{i}(E_{\nu})=\int_{x_{0}}^{1}\frac{dx}{x}\int_{Q^{2}_{\rm
min}}^{Q^{2}_{\rm
max}(x)}dQ^{2}\frac{d^{2}\sigma_{i}}{dxdQ^{2}}(x,Q^{2},E_{\nu})\,,$ (1)

where $Q^{2}_{\rm max}(x)=2m_{p}E_{\nu}x$ and $x_{0}=Q^{2}_{\rm
min}/(2m_{p}E_{\nu})$, with $m_{p}$ the proton mass. While the NOMAD
measurements are reconstructed for $Q^{2}\geq 1$ GeV2, we assume $Q^{2}_{\rm
min}=m_{c}^{2}$, where $m_{c}=1.51$ GeV is the initial parametrization scale
adopted in our analysis Ball  _et al._ (2017). We explicitly verified that
results are unaffected if $Q^{2}_{\rm min}=1$ GeV2 is chosen instead. The
integrand in Eq. (1) is either the dimuon ($i=\mu\mu$) or the inclusive
charged-current ($i={\rm CC}$) cross section. We compute Eq. (1) and the ratio
$\mathcal{R}_{\mu\mu}$ with APFEL Bertone  _et al._ (2014); our results agree
(at permille level) with those of Gao (2018); Bertone  _et al._ (2019) (after
the correction of a bug in APFEL), see the SM for the details of the numerical
benchmark. Theoretical predictions for inclusive $W$\- and $Z$-boson
production and for $W$-boson production with charm quarks and light jets are
evaluated at NLO using MCFM+APPLgrid Boughezal  _et al._ (2017); Carli  _et
al._ (2010), and are supplemented with NNLO QCD $K$-factors evaluated with
FEWZ Gavin  _et al._ (2011). Because NNLO QCD corrections are unknown for
$W$$+$$c$ production, in this case we accompany the data with an additional
correlated uncertainty, estimated from the 9 point scale variations of the NLO
calculation Abdul Khalek  _et al._ (2019a, b).

Table 1: A list of the PDF fits presented in this work, see text for details on the fit ID, dataset and theory settings. Fit ID | Dataset | Theory
---|---|---
str_base | Same as Ball  _et al._ (2018) | Same as Ball  _et al._ (2018) +
NNLO $K$-fact NuTeV
$F_{2}^{c}$ positivity
str_prior | Same as str_base + | Same as str_base +
ATLAS $W$, $Z$ (full) Aaboud  _et al._ (2017) | correlated theory
ATLAS $W$$+$$c$ Aad  _et al._ (2014) | uncertainty for
CMS $W$$+$$c$ Chatrchyan  _et al._ (2014); Sirunyan  _et al._ (2019) | unknown NNLO QCD
ATLAS $W$+jets Aaboud  _et al._ (2018b) | corrections
str | Same as str_prior + | Same as str_prior +
NOMAD $E_{\nu}$ set Samoylov  _et al._ (2013) | NNLO $K$-fact NOMAD

Procedure. We perform the three fits summarized in Table 1. The first fit
(str_base) is our baseline, and corresponds to the fit of Ball  _et al._
(2018) with the addition of the NNLO charm-mass $K$-factors for the NuTeV data
and of the positivity constraint on $F_{2}^{c}$. See the SM for a comparison
with the NNPDF3.1 PDF set. This fit is then supplemented with all the new LHC
data to obtain the second fit (str_prior), for which we generate $N_{\rm
rep}=850$ Monte Carlo replicas. This second fit is finally supplemented with
the NOMAD data, specifically the set that depends on $E_{\nu}$, to determine
the third fit (str). Bayesian reweighting and unweighting Ball  _et al._
(2012, 2011b) are used in this last step, because they allow one to evaluate
the two-dimensional integral in Eq. (1) only once, a task that would otherwise
be computationally very intensive in a fit. After reweighting, one ends up
with $N_{\rm eff}=105$ effective replicas, from which we construct a set of
$N_{\rm rep}=100$ replicas. Variants of the str fit based on a perturbatively
generated charm PDF and on the inclusion of the alternative NOMAD
$\sqrt{\hat{s}}$ dataset are discussed in the SM.

Results. In Table 2 we summarize the values of the $\chi^{2}$ per data point
obtained from each of the three fits and for the datasets discussed above:
$\chi^{2}_{\rm bs}$ for str_base; $\chi^{2}_{\rm pr}$ for str_prior; and
$\chi^{2}_{\rm str}$ for str. By comparing these values across the three fits,
we observe that the description of the new datasets — which, in particular, is
not optimal for the ATLAS $W$, $Z$ data in the str_base fit and for the NOMAD
data in the str_base and str_prior fits — markedly improves as soon as they
are included in subsequent fits. The largest effect is witnessed by the NOMAD
data, whose $\chi^{2}$ decreases from about 9 in the str_base and str_prior
fits to about 0.6 in the str fit. The $\chi^{2}$ for all of the other datasets
is in general not affected upon addition of the NOMAD data in the str fit,
except for the NuTeV dataset, whose $\chi^{2}$ is further reduced in
comparison to the str_prior fit. We therefore conclude that the dataset is
overall consistent and well described in the str fit. The NOMAD dataset is
illustrated in Fig. 1, where we compare the corresponding theoretical
predictions, obtained from the str_prior and str fits (that is, before and
after the inclusion of the NOMAD data), and the experimental measurements. See
the SM for similar comparisons for other strangeness-sensitive datasets.

Table 2: Values of the $\chi^{2}$ per data point for the strangeness-sensitive datasets discussed in this work obtained from the three fits summarized in Table 1. Values in square brackets correspond to datasets not included in the corresponding fit. Process |  Dataset | $n_{\rm dat}$ | $\chi^{2}_{\rm bs}$ | $\chi^{2}_{\rm pr}$ | $\chi^{2}_{\rm str}$
---|---|---|---|---|---
$\nu$DIS ($\mu\mu$) |  | 76/76/95 | 0.76 | 0.71 | 0.53
|  NuTeV Mason  _et al._ (2007) | 76/76/76 | 0.76 | 0.71 | 0.53
|  NOMAD Samoylov  _et al._ (2013) | —/—/19 | [9.3] | [8.8] | 0.55
$W$, $Z$ (incl.) |  | 391/418/418 | 1.45 | 1.40 | 1.40
|  ATLAS Aaboud  _et al._ (2017) | 34/61/61 | 1.96 | 1.65 | 1.67
$W$$+$$c$ |  | —/37/37 | [0.73] | 0.68 | 0.60
|  CMS Chatrchyan  _et al._ (2014); Sirunyan  _et al._ (2019) | —/15/15 | [1.04] | 0.98 | 0.96
|  ATLAS Aad  _et al._ (2014) | —/22/22 | [0.52] | 0.48 | 0.42
$W$+jets |  ATLAS Aaboud  _et al._ (2018b) | —/32/32 | [1.58] | 1.18 | 1.18
Total |  | 3981/4077/4096 | 1.18 | 1.17 | 1.17
Figure 1: Comparison between the theoretical predictions, obtained from the
str_prior and str fits, and the experimental data for the NOMAD $E_{\nu}$
dataset. The inset displays the ratio to the central value of each measured
data point. Figure 2: The ratio $R_{s}$ at $Q=10$ GeV computed from the
str_base, str_prior and str fits. The inset displays the relative PDF
uncertainty on the ratio, $\delta R_{s}/R_{s}$. Figure 3: A comparison of
$R_{s}$ (top) and $s^{+}$ (bottom) obtained from the str fit at $Q=10$ GeV and
the corresponding central values obtained from other recent PDF
determinations. In the bottom panel we normalise to the central value of str,
and the insets display the relative PDF uncertainties for each case.

We now turn to study the impact of the new datasets included in the fits on
the ratio $R_{s}$. The corresponding assessment for individual PDF flavors is
provided in the SM. In Fig. 2 we display the ratio $R_{s}$, and its relative
PDF uncertainty $\delta R_{s}/R_{s}$, at $Q=10$ GeV for the three fits
discussed above. The impact of the new datasets is clearly visible. Concerning
the central value, collider datasets do not alter its expectation (the results
obtained from the str_base and str_prior fits are almost identical); the NOMAD
dataset, instead, prefers a somewhat more suppressed strange sea for
$x\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$>$}}0.1$. Concerning uncertainties, collider datasets lead to a
reduction of the relative uncertainty on $R_{s}$ of about 4% for
$x\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$<$}}0.1$; the NOMAD dataset, instead, reduces it by about a factor
of two for $x\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}0.1$. Overall, the impact of the new
datasets depends on $x$, and is mostly significant for $x=0.2$, where the
uncertainty on $R_{s}$ is reduced from 20% to 8%. For $x\mathrel{\hbox
to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}0.3$
no experimental constraints are available, therefore the PDF uncertainty blows
up. These conclusions are unaffected if the NOMAD $\sqrt{\hat{s}}$ dataset is
used instead of the $E_{\nu}$ one.

Fig. 3 compares the str fit with the results obtained from the CT18/CT18A Hou
_et al._ (2019) (CT18A is a variant of CT18 that includes the ATLAS $W$, $Z$
data), MMHT14 Harland-Lang  _et al._ (2015), and ABMP16 Alekhin  _et al._
(2017) fits. They all include only a subset of the data listed in Table 2, in
particular: the NuTeV dataset is part of all PDFs; the NOMAD dataset is only
part of ABMP16; and the off-peak and forward ATLAS $W$, $Z$ bins, the
$W$$+$$c$ and the $W$+jets datasets are not part of any of these PDF sets. The
upper (lower) panel of Fig. 3 displays the absolute (normalised) central
values of $R_{s}$ ($s^{+}$), with the insets displaying the relative PDF
uncertainties for each case. Our str determination agrees with the CT18A and
ABMP16 results within uncertainties in the data region; it, however,
overshoots the CT18 and MMHT14 results. Note that the very small PDF
uncertainties of the ABMP16 result should be realistically rescaled by a
tolerance factor $T=\chi^{2}>1$ Gao  _et al._ (2018), which is however not
accounted for in their analysis. With this caveat, our results for $s^{+}$ and
$R_{s}$ are also the most precise, in particular around $x\sim 0.1$, thanks to
the wider dataset (and specifically of NOMAD) utilized to constrain the
strange quark and anti-quark PDFs.

Fig. 4 displays the values of $R_{s}$ for the three fits of Table 1 and the
fits shown in Fig. 3; here $R_{s}$ is evaluated at $Q=1.6$ GeV and $x=0.023$
as in the ATLAS studies Aad  _et al._ (2012a); Aaboud  _et al._ (2017)
claiming a symmetric strange quark sea. Fig. 4 makes it clear the consistent
effect of the new datasets included in our analysis: the value of
$R_{s}=0.78\pm 0.20$ in the str_base fit is made more precise by the LHC
datasets, which reduce its uncertainty by about a third, without altering much
its central value, $R_{s}=0.76\pm 0.12$; then the neutrino-DIS NOMAD dataset
shifts this number towards a lower value by a half-sigma bringing in also a
further moderate reduction of the uncertainty, $R_{s}=0.71\pm 0.10$. A larger
reduction of the uncertainty on $R_{s}$, because of NOMAD data, is obtained if
$R_{s}$ is evaluated at a higher value of $x$, see Fig. 2, and the SM.

Our final str result indicates that the strange sea is neither highly
suppressed ($R_{s}\simeq 0.5$), as suggested by neutrino data, nor fully
symmetric ($R_{s}\simeq 1$), as allegedly preferred by the ATLAS $W$, $Z$
production measurements. Actually $R_{s}$ turns out to lie halfway between
these two limiting scenarios. We find similar conclusisons by studying the
momentum fraction $K_{s}$, see the SM. Our result is in agreement with the
determination of $R_{s}$ obtained from other recent PDF determinations within
uncertainties. We therefore conclude that the result $R_{s}=1.13\pm 0.11$,
reported in Aaboud  _et al._ (2017) from an analysis of HERA and ATLAS $W$,
$Z$ data within the xFitter framework Alekhin  _et al._ (2015b), is not
compatible with ours, possibly because it is affected by a restricted dataset
and/or methodological limitations.

Figure 4: The ratio $R_{s}$ computed at $x=0.023$ and $Q=1.6$ GeV from the
three fits of Table 1 and other recent PDF analyses.

We finally emphasize that, apart from the more extensive dataset, our analysis
differs from all of the other PDF determinations shown in Figs. 3-4 in that
the charm-quark PDF is fitted on the same footing as the other light-quark
PDFs Ball  _et al._ (2016). This feature was demonstrated to improve the
description of DIS and LHC datasets, and in particular to partially relieve
tensions between the NuTeV and the ATLAS $W$, $Z$ datasets Ball  _et al._
(2017). The interplay between the charm- and the strange-quark PDFs is further
addressed in the SM, where we find that revisiting our analysis with a
perturbative charm PDF leads to a worse fit quality while affecting only
marginally our conclusions on $R_{s}$.

Summary. By means of a state-of-the-art global analysis, which combines all
the relevant experimental and theoretical inputs, we have achieved a precise
determination of the strangeness content of the proton. We have demonstrated
the compatibility of all strangeness-sensitive datasets; quantified their
relative impact on the fit; compared our results to other recent global
analyses; and assessed the robustness of our results with respect to
methodological choices. Our str PDF set, available in the LHAPDF format
Buckley  _et al._ (2015) together with its perturbative charm
counterpart,111See http://nnpdf.mi.infn.it/nnpdf3-1strangeness/. represents an
important input for phenomenology. Our determination of the strange and anti-
strange quark PDFs could be further stress tested with more exclusive
processes, e.g., measurements of kaon production in semi-inclusive DIS
(SIDIS). Studies of the strange PDFs based on SIDIS Airapetian  _et al._
(2014); Borsa  _et al._ (2017); Sato  _et al._ (2019) notoriously prefer a
suppressed strangeness, but are also subject to the potential bias coming from
their sensitivity to the fragmentation of the strange quarks into kaons.

In short, our analysis demonstrates that no proton strangeness crisis occurs,
that the strange PDF can be precisely determined, and that the proton is not
as strange as it has often thought to be.

Acknowledgments. We are grateful to Jun Gao for providing us with the NNLO
$K$-factors for the NuTeV and the NOMAD measurements, and for help into the
benchmark of the corresponding NLO calculations. We thank Valerio Bertone for
assistance with the usage of APFEL, and Rabah Abdul-Khalek, Stefano Forte, and
Cameron Voisey for a critical reading of the manuscript. E.R.N. is supported
by the European Commission through the Marie Skłodowska-Curie Action ParDHonS
FFs.TMDs (grant number 752748). J.R. is partially supported by the Netherlands
Organization for Scientific Research (NWO). M.U. and S.I. are partially
supported by the STFC grant ST/L000385/1 and by the Royal Society grant
RGF/EA/180148. The work of M.U. is also funded by the Royal Society grant
DH150088.

## References

  * Gao  _et al._ (2018) J. Gao, L. Harland-Lang, and J. Rojo, Phys. Rept. 742, 1 (2018), arXiv:1709.04922 [hep-ph] .
  * Kovarik  _et al._ (2019) K. Kovarik, P. M. Nadolsky, and D. E. Soper, (2019), arXiv:1905.06957 [hep-ph] .
  * Ethier and Nocera (2020) J. J. Ethier and E. R. Nocera, Ann. Rev. Nucl. Part. Sci. , 1 (2020), arXiv:2001.07722 [hep-ph] .
  * Aaboud  _et al._ (2018a) M. Aaboud _et al._ (ATLAS), Eur. Phys. J. C78, 110 (2018a), arXiv:1701.07240 [hep-ex] .
  * Sirunyan  _et al._ (2018) A. M. Sirunyan _et al._ (CMS), Eur. Phys. J. C78, 701 (2018), arXiv:1806.00863 [hep-ex] .
  * Bagnaschi and Vicini (2019) E. Bagnaschi and A. Vicini, (2019), arXiv:1910.04726 [hep-ph] .
  * Bazarko  _et al._ (1995) A. O. Bazarko _et al._ (CCFR), Z. Phys. C65, 189 (1995), arXiv:hep-ex/9406007 .
  * Kayis-Topaksu  _et al._ (2008) A. Kayis-Topaksu _et al._ (CHORUS), Nucl. Phys. B 798, 1 (2008), arXiv:0804.1869 [hep-ex] .
  * Mason  _et al._ (2007) D. Mason _et al._ , Phys. Rev. Lett. 99, 192001 (2007).
  * Samoylov  _et al._ (2013) O. Samoylov _et al._ (NOMAD), Nucl.Phys. B876, 339 (2013), arXiv:1308.4750 [hep-ex] .
  * Aad  _et al._ (2012a) G. Aad _et al._ (ATLAS), Phys.Rev.Lett. (2012a), arXiv:1203.4051 [hep-ex] .
  * Aaboud  _et al._ (2017) M. Aaboud _et al._ (ATLAS), Eur. Phys. J. C77, 367 (2017), arXiv:1612.03016 [hep-ex] .
  * Sutton (2019) M. Sutton (ATLAS), PoS DIS2019, 034 (2019).
  * Stirling and Vryonidou (2012) W. Stirling and E. Vryonidou, Phys.Rev.Lett. 109, 082002 (2012), arXiv:1203.6781 [hep-ph] .
  * Aaboud  _et al._ (2018b) M. Aaboud _et al._ (ATLAS), JHEP 05, 077 (2018b), arXiv:1711.03296 [hep-ex] .
  * Aad  _et al._ (2014) G. Aad _et al._ (ATLAS), JHEP 1405, 068 (2014), arXiv:1402.6263 [hep-ex] .
  * Chatrchyan  _et al._ (2014) S. Chatrchyan _et al._ (CMS), JHEP 02, 013 (2014), arXiv:1310.1138 [hep-ex] .
  * Sirunyan  _et al._ (2019) A. M. Sirunyan _et al._ (CMS), Eur. Phys. J. C79, 269 (2019), arXiv:1811.10021 [hep-ex] .
  * Ball  _et al._ (2017) R. D. Ball _et al._ (NNPDF), Eur. Phys. J. C77, 663 (2017), arXiv:1706.00428 [hep-ph] .
  * Hou  _et al._ (2019) T.-J. Hou _et al._ , (2019), arXiv:1912.10053 [hep-ph] .
  * Lai  _et al._ (2007) H. L. Lai _et al._ , JHEP 04, 089 (2007), arXiv:hep-ph/0702268 .
  * Kusina  _et al._ (2012) A. Kusina, T. Stavreva, S. Berge, F. Olness, I. Schienbein, K. Kovarik, T. Jezo, J. Yu, and K. Park, Phys. Rev. D 85, 094028 (2012), arXiv:1203.1290 [hep-ph] .
  * Alekhin  _et al._ (2015a) S. Alekhin, J. Blumlein, L. Caminadac, K. Lipka, K. Lohwasser, S. Moch, R. Petti, and R. Placakyte, Phys. Rev. D91, 094002 (2015a), arXiv:1404.6469 [hep-ph] .
  * Alekhin  _et al._ (2018) S. Alekhin, J. Blümlein, and S. Moch, Phys. Lett. B 777, 134 (2018), arXiv:1708.01067 [hep-ph] .
  * Cooper-Sarkar and Wichmann (2018) A. Cooper-Sarkar and K. Wichmann, Phys. Rev. D 98, 014027 (2018), arXiv:1803.00968 [hep-ex] .
  * Ball  _et al._ (2015) R. D. Ball _et al._ (NNPDF), JHEP 04, 040 (2015), arXiv:1410.8849 [hep-ph] .
  * Ball  _et al._ (2018) R. D. Ball, S. Carrazza, L. Del Debbio, S. Forte, Z. Kassabov, J. Rojo, E. Slade, and M. Ubiali (NNPDF), Eur. Phys. J. C78, 408 (2018), arXiv:1802.03398 [hep-ph] .
  * Aaltonen  _et al._ (2010) T. A. Aaltonen _et al._ (CDF), Phys. Lett. B692, 232 (2010), arXiv:0908.3914 [hep-ex] .
  * Abazov  _et al._ (2007) V. M. Abazov _et al._ (D0), Phys. Rev. D76, 012003 (2007), arXiv:hep-ex/0702025 .
  * Aad  _et al._ (2012b) G. Aad _et al._ (ATLAS), Phys.Rev. D85, 072004 (2012b), arXiv:1109.5141 [hep-ex] .
  * Khachatryan  _et al._ (2016) V. Khachatryan _et al._ (CMS), Eur. Phys. J. C76, 469 (2016), arXiv:1603.01803 [hep-ex] .
  * Berger  _et al._ (2016) E. L. Berger, J. Gao, C. S. Li, Z. L. Liu, and H. X. Zhu, Phys. Rev. Lett. 116, 212002 (2016), arXiv:1601.05430 [hep-ph] .
  * Gao (2018) J. Gao, JHEP 02, 026 (2018), arXiv:1710.04258 [hep-ph] .
  * Forte  _et al._ (2010) S. Forte, E. Laenen, P. Nason, and J. Rojo, Nucl. Phys. B834, 116 (2010), arXiv:1001.2312 [hep-ph] .
  * Ball  _et al._ (2011a) R. D. Ball _et al._ (The NNPDF ), Nucl. Phys. B849, 296 (2011a), arXiv:1101.1300 [hep-ph] .
  * Ball  _et al._ (2019) R. D. Ball, E. R. Nocera, and R. L. Pearson (NNPDF), Eur. Phys. J. C79, 282 (2019), arXiv:1812.09074 [hep-ph] .
  * Bertone  _et al._ (2014) V. Bertone, S. Carrazza, and J. Rojo, Comput.Phys.Commun. 185, 1647 (2014), arXiv:1310.1394 [hep-ph] .
  * Bertone  _et al._ (2019) V. Bertone, R. Gauld, and J. Rojo, JHEP 01, 217 (2019), arXiv:1808.02034 [hep-ph] .
  * Boughezal  _et al._ (2017) R. Boughezal, J. M. Campbell, R. K. Ellis, C. Focke, W. Giele, X. Liu, F. Petriello, and C. Williams, Eur. Phys. J. C 77, 7 (2017), arXiv:1605.08011 [hep-ph] .
  * Carli  _et al._ (2010) T. Carli _et al._ , Eur.Phys.J. C66, 503 (2010), arXiv:0911.2985 [hep-ph] .
  * Gavin  _et al._ (2011) R. Gavin, Y. Li, F. Petriello, and S. Quackenbush, Comput. Phys. Commun. 182, 2388 (2011), arXiv:1011.3540 [hep-ph] .
  * Abdul Khalek  _et al._ (2019a) R. Abdul Khalek _et al._ (NNPDF), Eur. Phys. J. C 79, 931 (2019a), arXiv:1906.10698 [hep-ph] .
  * Abdul Khalek  _et al._ (2019b) R. Abdul Khalek _et al._ (NNPDF), Eur. Phys. J. C, 79:838 (2019b), arXiv:1905.04311 [hep-ph] .
  * Ball  _et al._ (2012) R. D. Ball, V. Bertone, F. Cerutti, L. Del Debbio, S. Forte, _et al._ , Nucl.Phys. B855, 608 (2012), arXiv:1108.1758 [hep-ph] .
  * Ball  _et al._ (2011b) R. D. Ball _et al._ (The NNPDF), Nucl. Phys. B849, 112 (2011b), arXiv:1012.0836 [hep-ph] .
  * Harland-Lang  _et al._ (2015) L. A. Harland-Lang, A. D. Martin, P. Motylinski, and R. S. Thorne, Eur. Phys. J. C75, 204 (2015), arXiv:1412.3989 [hep-ph] .
  * Alekhin  _et al._ (2017) S. Alekhin, J. Blümlein, S. Moch, and R. Placakyte, Phys. Rev. D96, 014011 (2017), arXiv:1701.05838 [hep-ph] .
  * Alekhin  _et al._ (2015b) S. Alekhin _et al._ , Eur. Phys. J. C75, 304 (2015b), arXiv:1410.4412 [hep-ph] .
  * Ball  _et al._ (2016) R. D. Ball, V. Bertone, M. Bonvini, S. Carrazza, S. Forte, A. Guffanti, N. P. Hartland, J. Rojo, and L. Rottoli (NNPDF), Eur. Phys. J. C76, 647 (2016), arXiv:1605.06515 [hep-ph] .
  * Buckley  _et al._ (2015) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, _et al._ , Eur.Phys.J. C75, 132 (2015), arXiv:1412.7420 [hep-ph] .
  * Airapetian  _et al._ (2014) A. Airapetian _et al._ (HERMES), Phys. Rev. D89, 097101 (2014), arXiv:1312.7028 [hep-ex] .
  * Borsa  _et al._ (2017) I. Borsa, R. Sassot, and M. Stratmann, Phys. Rev. D96, 094020 (2017), arXiv:1708.01630 [hep-ph] .
  * Sato  _et al._ (2019) N. Sato, C. Andres, J. J. Ethier, and W. Melnitchouk (JAM), (2019), arXiv:1905.03788 [hep-ph] .
  * Abdul Khalek  _et al._ (2020) R. Abdul Khalek, J. J. Ethier, J. Rojo, and G. van Weelden, (2020), arXiv:2006.14629 [hep-ph] .

Supplementary Material

Computation of the NOMAD observables. As explained in the main text, the NOMAD
experiment measured the ratio of dimuon to inclusive charged-current cross
sections, $\mathcal{R}_{\mu\mu}$, as a function of the neutrino beam energy
$E_{\nu}$, the momentum fraction $x$, and the square root of the partonic
center-of-mass energy $\sqrt{\hat{s}}$. The ratio $\mathcal{R}_{\mu\mu}$ as a
function of $E_{\nu}$ is evaluated by means of Eq. (1), where the expression
of the integrand reads

$\frac{d^{2}\sigma_{i}}{dxdQ^{2}}(x,Q^{2},E_{\nu})=\frac{G_{F}^{2}M_{W}^{2}}{4\pi}\frac{1}{(Q^{2}+M_{W}^{2})^{2}}\left[\left(Y_{+}-\frac{2m_{p}^{2}x^{2}y^{2}}{Q^{2}}\right)F_{2}^{i}(x,Q^{2})-y^{2}F_{L}^{i}(x,Q^{2})+Y_{-}xF_{3}^{i}\right]K^{i}\,,$
(2)

with the index $i=\mu\mu,\,{\rm CC}$ denoting either the dimuon or the
inclusive charged-current quantities. They enter, respectively, the numerator
and the denominator of $\mathcal{R}_{\mu\mu}$. The kinematic factors
$Y_{\pm}=1\pm(1-y)^{2}$ are related to the inelasticity $y$, which can in turn
be expressed as $y=Q^{2}/(2m_{p}E_{\nu}x)$; $G_{F}$, $M_{W}$ and $m_{p}$ are
respectively the Fermi constant, the mass of the $W$ boson and the mass of the
proton. The factor $K^{i}$ is either the identity, for $i={\rm CC}$, or the
charm semileptonic branching ratio $B_{\mu}$, for $i=\mu\mu$. In the latter
case we use the $E_{\nu}$-dependent parametrization
$B_{\mu}(E_{\nu})=a(1+b/E_{\nu})^{-1}$, with the values of the parameters $a$
and $b$ determined in Samoylov  _et al._ (2013), $a=0.097\pm 0.003$ and
$b=6.7\pm 1.8$. The corresponding uncertainty is included in the covariance
matrix of the measurement. Both the charm (for $i=\mu\mu$) and the total (for
$i={\rm CC}$) structure functions $F_{p}^{i}$ ($p=2,L,3$) are evaluated with
APFEL Bertone  _et al._ (2014).

For illustrative purposes, in the left panel of Fig. 5 we display the charm
production cross section, Eq. (1) with $i=\mu\mu$, as a function of $E_{\nu}$
in the kinematic range measured by the NOMAD experiment. The cross section is
obtained in the FFN scheme (with $n_{f}=3$) at different perturbative orders
using the NNPDF3.1 NNLO PDF set (consistently with $n_{f}=3$). The inset
displays the ratio to the leading order (LO) calculation. Higher-order
corrections clearly suppress the cross section, in particular as $E_{\nu}$
increases. For instance, in the highest energy bin the NNLO cross section is
about 10% smaller than the LO prediction. The size of the NNLO correction is
comparable or larger than the size of the NLO one, therefore its inclusion is
mandatory to achieve a good description of the NOMAD data.

While the comparison of Fig. 5 is presented in the FFN scheme, all the fits
discussed in this work are based on the FONLL general-mass variable flavor
number scheme. As discussed in Ball  _et al._ (2011a), the theoretical
predictions for neutrino-DIS charm production obtained in either the FFN or in
the FONLL schemes turn out to be very close in the kinematic range accessed by
the NuTeV experiment. We explicitly checked this statement, and that the same
applies also to the NOMAD measurements. To this purpose, we explicitly
computed the relative difference between the FONLL-A and FFN scheme
predictions for the NuTeV and NOMAD datasets based on structure functions
accurate to $\mathcal{O}\left(\alpha_{s}\right)$. We found that differences
were less than 1% in the entire kinematic range for NuTeV, and of about 1.5%
irrespective of the value of $E_{\nu}$ for NOMAD. These differences are well
below the experimental and the PDF uncertainties. We therefore conclude that
using a NNLO K-factor determined in the FFN scheme in fits that otherwise use
the FONLL scheme is unlikely to affect our results.

Figure 5: Left: the integrated dimuon cross section, Eq. (1) ($i=\mu\mu$),
computed as a function of the neutrino beam energy $E_{\nu}$ in the kinematic
range measured by the NOMAD experiment. The cross section is computed in the
FFN scheme ($n_{f}=3$) with the NNPDF3.1 NNLO PDF set (also the $n_{f}=3$
version) for various perturbative orders. The inset displays the ratio to the
LO calculation. Right: the same cross section, accurate to NLO, computed with
APFEL (apfel) and with a code based on Gao (2018) (ref). The ratio of the two
computations is shown in the inset.

We tested the robustness of our computation of Eqs. (1)-(2) in two further
respects. First, we performed a benchmark against the independent computation
based on Gao (2018). In the process, a bug affecting the APFEL calculation of
NLO charged-current coefficient functions was identified and corrected, a fact
that made our computation in excellent agreement with that based on Gao
(2018). Second, we estimated the impact of nuclear corrections on our
predictions. To this purpose, we recomputed them with the recently presented
nNNPDF2.0 NLO Fe nuclear PDF set Abdul Khalek  _et al._ (2020) and we compared
the result with the predictions obtained with the NLO free proton PDF set
consistently determined in Abdul Khalek  _et al._ (2020). The full set of
correlations between Fe and proton PDFs were therefore appropriately taken
into account. The relative difference between the two computations (with and
without nuclear PDF corrections) turned out to range between 3%, in the lowest
$E_{\nu}$ bin, and a fraction of percent, in the bins at the highest
$E_{\nu}$. These differences are smaller than both the data and PDF
uncertainties, therefore neglecting nuclear PDF uncertainties is a justified
approximation which should have a negligible impact on our results. We also
note that the effects of nuclear corrections due to a Fe target were
explicitly studied in Ball  _et al._ (2019) for the NuTeV case and found to be
immaterial in a global fit of PDFs.

Stability with respect to the choice of NOMAD observable. The NOMAD cross
section ratio $\mathcal{R}_{\mu\mu}$ is presented in terms of three kinematic
variables characteristic of the DIS process: the energy of the neutrino beam
$E_{\nu}$, the momentum fraction $x$, and the square root of the invariant
mass of the final state (or partonic center-of-mass energy) $\sqrt{\hat{s}}$.
Because the three datasets are reconstructed from the same underlying
measurement, and the associated experimental correlation matrix is unknown, it
is not possible to include the three of them simultaneously in a fit without
incurring in double counting. The best-fit str presented in this work is based
on the NOMAD cross section ratio as a function of the neutrino energy
$E_{\nu}$. We selected this specific set because, among the three kinematic
variables, $E_{\nu}$ is the one more directly accessed in the experimental
setup.

Process |  Dataset | $n_{\rm dat}$ | $\chi^{2}_{\rm str}$ | $\chi^{2}_{\rm str\\_s\\_hat}$ | $\chi^{2}_{\rm str\\_pch}$
---|---|---|---|---|---
$\nu$DIS ($\mu\mu$) |  | 95/95/95 | 0.53 | 0.52 | 0.63
|  NuTeV Mason  _et al._ (2007) | 76/76/76 | 0.53 | 0.55 | 0.61
|  NOMAD Samoylov  _et al._ (2013) | 19/15/19 | 0.55 | 0.35 | 0.69
$W,Z$ (incl.) |  | 418/418/418 | 1.40 | 1.39 | 1.40
|  ATLAS Aaboud  _et al._ (2017) | 61/61/61 | 1.67 | 1.64 | 1.80
$W$$+$$c$ |  | 37/37/37 | 0.60 | 0.66 | 0.68
|  CMS Chatrchyan  _et al._ (2014); Sirunyan  _et al._ (2019) | 15/15/15 | 0.96 | 1.00 | 1.00
|  ATLAS Aad  _et al._ (2014) | 22/22/22 | 0.42 | 0.43 | 0.46
$W$+jets |  ATLAS Aaboud  _et al._ (2018b) | 32/32/32 | 1.18 | 1.18 | 1.18
Total |  | 4096/4092/4096 | 1.17 | 1.17 | 1.20
Table 3: Same as Table 2 now comparing the results of fits str (fitted charm,
NOMAD $E_{\nu}$ dataset), str_s_hat (fitted charm, NOMAD $\sqrt{\hat{s}}$
dataset), and str_pch (perturbative charm, NOMAD $E_{\nu}$ dataset).

Figure 6: Left: same as Fig. 1 now for the fit based on the NOMAD
$\sqrt{\hat{s}}$ distribution str_s_hat. Right: comparison of the ratio
$R_{s}$ computed at $Q=10$ GeV with fits based on either the NOMAD $E_{\nu}$
(str) or the NOMAD $\sqrt{\hat{s}}$ (str_s_hat) datasets.

Figure 7: Left: same as Fig. 1 now for corresponding perturbative charm fits.
Right: same as the right panel of Fig. 6 now comparing the ratio $R_{s}$ in
the fits str (fitted charm) and str_pch (perturbative charm).

Here we study the stability of our results when the NOMAD $E_{\nu}$ dataset is
replaced by its $\sqrt{\hat{s}}$ counterpart in the fit. To this purpose, we
reweight the str_prior fit with the NOMAD $\sqrt{\hat{s}}$ dataset, to obtain
the str_s_hat fit. After reweighting, one ends up with $N_{\rm eff}=135$
effective replicas (out of $N_{\rm rep}=850$ initial replicas), from which we
construct an ensemble of $N_{\rm rep}=100$ replicas. The number of effective
replicas is similar to that obtained after reweighting with the $E_{\nu}$
dataset. In Table 3 we report the values of the $\chi^{2}$ per data point
computed from the str_s_hat fit for all of the strangeness-sensitive datasets
considered in this work. The format is similar to that used in Table 2, from
which we report, for ease of comparison, the corresponding values of the
$\chi^{2}$ obtained from our best-fit str. We observe that a very similar fit
quality is achieved in the str_s_hat and str fits, not only for the NOMAD
data, but also for all of the other datasets: the differences in the values of
the $\chi^{2}$ between the two fits are smaller than statistical fluctuations.
The impact of the NOMAD $\sqrt{\hat{s}}$ dataset is further displayed in Fig.
6. In the left panel, we compare the data and the corresponding theoretical
predictions computed with the str_base and str_s_hat fits. We observe a
reduction of the relative uncertainty similar to that reported in Fig. 1 for
the NOMAD $E_{\nu}$ dataset. In the right panel, we compare the ratio $R_{s}$
computed from the str and str_s_hat fits: both the central values and the PDF
uncertainties of $R_{s}$ are very similar. This fact confirms the robustness
of our analysis, and the independence of the results upon the choice of the
NOMAD dataset included in the fit.

Figure 8: Comparison between theoretical predictions and experimental data for
some of the measurements used in this work, in addition to those from the
NOMAD experiment presented in Fig. 1. Top: the $W$$+$$c$ lepton rapidity
distributions (separately for $W^{+}$ and $W^{-}$) corresponding to the ATLAS
measurement at 7 TeV Aad  _et al._ (2014). Bottom: the $W$$+$$c$ lepton
rapidity distributions (sum of $W^{+}$ and $W^{-}$) corresponding to the CMS
measurements at 7 TeV Chatrchyan  _et al._ (2014) and 13 TeV Sirunyan  _et
al._ (2019). The insets display the theory to data ratio. Theoretical
predictions are computed with the str_base and str fits.

Figure 9: Same as Fig. 8 now for the $Z$ dilepton rapidity distributions from
the ATLAS measurement of Aaboud  _et al._ (2017) for both the central and
forward selection cuts (top); and the charm dimuon cross sections from the
NuTeV measurement of Mason  _et al._ (2007), for both neutrino and
antineutrino beams (bottom). In the latter case only the ratio to data is
shown; data points are sorted by their ID value, roughly corresponding to
increasing $x$, $Q^{2}$, and $y$ values when the plot is read from left to
right.

Impact of the treatment of the charm PDF. The description of the charm- and
strange-quark PDFs is known to be correlated. For instance, charged-current
weak boson production at colliders is sensitive at LO to both strange and
charm quark PDFs via the $c\bar{s}$ and $s\bar{c}$ initial-state partonic
channels. For this reason, it is interesting to assess the robustness of our
results against the theoretical treatment of the charm-quark PDF in our fits.
To this purpose, we repeated the analysis carried out in the str_prior and str
fits assuming a purely perturbative charm-quark PDF (that is, generated from
the gluon and light quarks via DGLAP evolution). This assumption is used in
other global analyses discussed in this work, see Figs. 3-4. This way we
obtained the corresponding str_prior_pch and str_pch fits. In this case we
generated only $N_{\rm rep}=500$ replicas in the str_prior_pch fit; after
reweighting we are left with $N_{\rm eff}=157$ effective replicas, from which
we constructed an ensemble of $N_{\rm rep}=100$ replicas in the str_pch fit.

We collect the values of the $\chi^{2}$ per data point obtained from the
str_pch fit for the usual datasets in Table 3. For ease of comparison, the
corresponding values, obtained with the str fit, are also reported. We observe
that the perturbative charm fit (str) achieves a better description of the
strangeness-sensitive datasets, and of the global dataset overall, than the
fitted charm fit (str_pch). We note in particular the $\chi^{2}$ values of the
ATLAS $W$, $Z$ and of the total datasets, which increase respectively from
1.67 to 1.80 and from 1.17 to 1.20 when comparing the str and the str_pch
fits. We therefore confirm that fitting charm achieves an improved description
of the experimental data.

To further illustrate the impact of the perturbative charm assumption, in the
left panel of Fig. 7 we reproduce Fig. 1 now with the str_prior_pch and
str_pch fits. In comparison to the fitted charm case, see Fig. 1, the
prediction obtained from the str_prior_pch fit is closer to the NOMAD data and
has smaller uncertainties than its counterpart obtained with the str_prior
fit. The agreement after reweighting is, however, comparably good in the two
cases. Finally, the right panel of Fig. 7 displays the ratio $R_{s}$ obtained
from our fitted and perturbative charm best fits, str and str_pch. While PDF
uncertainties turn out to be very similar in the two cases, the perturbative
charm fit prefers a central value which is systematically larger than the one
obtained from the fitted charm fit. The size of the shift, however, is at most
at the one-sigma level in units of the PDF uncertainty bands, in line with
previous studies Ball  _et al._ (2017, 2016).

Comparison between experimental data and theory predictions. Here we extend
the comparison between theoretical predictions and data to some of the
measurements used in this work, other than those from the NOMAD experiment
presented in Fig. 1. In Figs. 8 and 9 we display: the $W$$+$$c$ lepton
rapidity distributions corresponding to the ATLAS measurement of Aad  _et al._
(2014) (for both $W^{+}$ and $W^{-}$) and to the CMS measurements (sum of
$W^{+}$ and $W^{-}$) of Chatrchyan  _et al._ (2014); Sirunyan  _et al._ (2019)
(respectively at 7 TeV and 13 TeV); the $Z$ dilepton rapidity distributions
from the ATLAS measurement of Aaboud  _et al._ (2017) (for both the central
and forward selection cuts); and the charm dimuon cross sections from the
NuTeV measurement of Mason  _et al._ (2007) (for both neutrino and
antineutrino beams). In the first three cases, the insets display the ratio of
the theory to the central value of the experimental measurement. In the last
case, only this ratio is shown; data points are sorted by their ID value,
roughly corresponding to increasing $x$, $Q^{2}$, and $y$ values when the plot
is read from left to right. Theoretical predictions are evaluated with the
str_base and str fits.

A fair agreement between data and theory is found in all cases, as expected
from the pattern of $\chi^{2}$ values reported in Table 2. However, we clearly
see that the size of the PDF uncertainty relative to the size of the data
uncertainty depend on the dataset. Concerning the ATLAS and CMS $W$$+$$c$
measurements, experimental uncertainties span the range between 10% and 20%,
and are consistently larger than PDF uncertainties. We note that the PDF
uncertainties in the theory predictions are markedly reduced in the str fit in
comparison to str_base fit, as highlighted by the ratios in the insets.
Concerning the ATLAS $Z$ measurement, the total experimental uncertainty is
much smaller, around 2% for the central rapidity bin, and is comparable to the
PDF uncertainty. We therefore expect this measurement to be one of the most
constraining amongst all of the LHC measurements considered in this work.
Interestingly, once the NOMAD dataset is included in the fit, the central
value of the theoretical prediction approaches the central value of the ATLAS
data, and PDF uncertainties are slightly reduced. A similar trend can be
observed for the forward selection data and for the NuTeV measurement. This
behaviour is a further sign of the good overall compatibility of all of the
datasets, and in particular of neutrino DIS and LHC gauge boson production
measurements.

Figure 10: The three PDF sets displayed in Fig. 2 now for other flavor
combinations: the gluon, the total quark singlet, and the charm, up, down, and
strange sea PDFs. The comparison is performed at $Q=100$ GeV. The upper panels
display the PDFs normalised to the str_base fit (ref); the lower panels
display the one-sigma uncertainty, relative to each fit (0). Figure 11: Same
as Fig. 10 for the up, down, and strange valence PDFs. Opposite to Fig. 10,
PDFs (upper panels) and their uncertainties (lower panels) are not normalised,
as suitable for valence-type distributions.

Implications for other PDF flavor combinations. In the main manuscript we
focused on the impact of the new datasets on the ratio $R_{s}$ and the total
quark-antiquark strange distribution $s^{+}$, see Figs. 2-3. We now briefly
discuss the corresponding impact on other PDF flavor combinations. In Fig. 10
we compare the gluon, the quark singlet, and the charm, up, down, and strange
sea PDFs resulting from the str_base, str_prior and str fits at $Q=100$ GeV.
The upper panels display the PDFs normalised to the str_base fit; the lower
panels display the corresponding one-sigma uncertainty relative to each fit.
In Fig. 11 we display a similar comparison now for the absolute PDFs and
uncertainties of the up, down, and strange quark valence distributions. From
these comparisons, we observe that the new datasets have a little impact on
the gluon and the quark singlet PDFs, both on central values and on
uncertainties, as expected. A bigger effect, entirely due to the NOMAD
dataset, is observed instead on the charm PDF: in the str fit, its central
value is suppressed in comparison to the str_prior fit; uncertainties are
reduced by up to a factor 2 for $x\simeq 0.05$. Note that the NOMAD data is
indirectly sensitive to the charm PDF through its interplay with the
$s\bar{c}$ and $\bar{s}c$ contributions to $W$-boson production. Concerning
valence distributions, the new datasets leave central values unchanged, but
reduce the uncertainties in the entire kinematic range. The pattern of
uncertainty reduction is the same for the three quark valence combinations:
the electroweak LHC datasets constrain the distributions at low to mid values
of $x$, $x\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}0.1$, while the NOMAD datasets do so
at larger values of $x$, $x\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}0.1$. The two datasets are therefore
complementary, and concur together to make all valence quark PDFs more
precise. Finally, we explicitly verified that the impact of the new datasets
on the up and down sea quark PDFs is minimal.

Additional metrics to quantify the strangeness suppression. In Fig. 4 we
displayed the ratio $R_{s}$ for a single kinematic point, $x=0.023$ and
$Q=1.6$ GeV, for ease of comparison with the ATLAS studies Aad  _et al._
(2012a); Aaboud  _et al._ (2017). However, as illustrated in Fig. 3, $R_{s}$
is a function of $x$ and $Q$. We therefore repeat the comparison displayed in
Fig. 4 for a different kinematic point, $x=0.13$ and $Q=100$ GeV, see Fig. 12.
In this case one finds values of $R_{s}$ typically smaller than in the
previous case. For example, for the str fit one can compare the value
$R_{s}=0.71\pm 0.10$ for $(x,Q)=(0.023,1.6~{}{\rm GeV})$ with the value
$R_{s}=0.64\pm 0.05$ for $(x,Q)=(0.13,100~{}{\rm GeV})$. We observe however
that the qualitative features displayed in Fig. 4, namely the reduction of PDF
uncertainties between the str_base and the str fits, and the relative size of
predictions obtained from the various PDF sets, are recovered in Fig. 12.

Another measure of the strange to light sea quark suppression is the ratio of
the corresponding momentum fractions

$K_{s}(Q^{2})\equiv\frac{\int_{0}^{1}dx\,x\left(s(x,Q^{2})+\bar{s}(x,Q^{2})\right)}{\int_{0}^{1}dx\,x\left(\bar{u}(x,Q^{2})+\bar{d}(x,Q^{2})\right)}\,.$
(3)

The middle and right panels of Fig. 12 display $K_{s}(Q^{2})$ for $Q=1.6$ GeV
and $Q=100$ GeV, respectively. The qualitative interpretation of this quantity
is consistent with that of $R_{s}$, in particular, PDF uncertainties are
reduced by a factor of two in the str fit with respect to the str_base fit.
The values of $K_{s}$ grow with the scale $Q$, as expected due to DGLAP
evolution effects: for instance, using the str fit, one finds $K_{s}=0.64\pm
0.07$ at $Q=1.6$ GeV, and $K_{s}=0.81\pm 0.04$ at $Q=100$ GeV. We therefore
conclude that our results are consistent upon different choices of metric used
to quantify the strangeness suppression with respect to the rest of the light
sea quarks.

Figure 12: Left panel: same as Fig. 4 now for $x=0.13$ and $Q=100$ GeV. Middle
and right panels: the ratio of the momentum fraction of the strange sea to
that of the up and down quark sea, Eq. (3), for $Q=1.6$ GeV and $Q=100$ GeV.
Figure 13: Comparison between the NNPDF3.1 NNLO fit Ball  _et al._ (2017) and
the baseline fit used in this work, str_base. We show the gluon, the quark
singlet, and the charm, up, down, and strange quark sea PDFs at $Q=100$ GeV,
normalized to the central value of NNPDF3.1 (ref).

Comparison with NNPDF3.1. Our baseline fit str_base differs from NNPDF3.1 Ball
_et al._ (2017) in several respects. These include: the treatment of inclusive
jet production from ATLAS and CMS with NNLO $K$-factors, see Ball  _et al._
(2018); an updated treatment of non-isoscalarity effects and branching
fractions in neutrino-DIS data, see Ball  _et al._ (2019); the inclusion of
the NNLO massive corrections to the NuTeV structure functions; the new
$F_{2}^{c}$ positivity constraint; and the correction of the APFEL bug found
in the benchmark reported in Fig. 5. For completeness, we compare the NNPDF3.1
and str_base parton sets in Fig. 13. Specifically we display the gluon, the
quark singlet, the charm, up, down, and strange quark sea PDFs at $Q=100$ GeV,
normalized to the central value of NNPDF3.1. In comparison to NNPDF3.1, in the
str_base fit we observe: an increase in the central value of the strange PDF
for $x\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$>$}}10^{-3}$ (mostly due to the correction of the APFEL bug); a
similar effect in the case of the charm PDF for $x\mathrel{\hbox
to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$>$}}10^{-2}$ (mostly due to the new $F_{2}^{c}$ positivity
constraint); a moderate rearrangement of the quark flavor separation at medium
and large-$x$; and a harder gluon at large-$x$ (mostly due to the improved
NNLO treatment of jet data). All in all, while the two fits agree within
uncertainties, the improvements introduced in the str_base fit justify its
adoption as the baseline for the present study.