††institutetext: University of Jyväskylä, Department of Physics, P.O. Box 35,
FI-40014 University of Jyväskylä, Finland††institutetext: Helsinki Institute
of Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland

# Dimuons from neutrino-nucleus collisions in the semi-inclusive DIS approach

Ilkka Helenius  Hannu Paukkunen and  Sami Yrjänheikki

###### Abstract

We present a next-to-leading order perturbative QCD calculation of dimuon
production in neutrino-nucleus collisions. This process is typically
calculated by assuming it to be proportional to inclusive charm production,
which requires an effective acceptance correction to take the experimental
cuts on the decay-muon kinematics into account. Here, we instead compute the
dimuon production cross section directly as a convolution of semi-inclusive
deep inelastic scattering to produce charmed hadrons, and a decay function
fitted to $e^{+}e^{-}$ data to produce a muon from the charmed hadrons. The
presented approach is in a good agreement with available experimental data and
will serve as a starting point for higher-order QCD calculations without an
external acceptance correction. The uncertainties arising from the decay
function and scale dependence are sizeably smaller than those from the nuclear
parton distribution functions. We also calculate the effective acceptances
within our approach and compare them to those usually used in global fits of
parton distribution functions, finding differences of the order of $10\,\%$,
depending on the kinematics, perturbative order, and applied parton
distributions.

###### Keywords:

Dimuon production, muon-pair production, nuclear parton distribution
functions, QCD, neutrino SIDIS

##  1 Introduction

A staple of Quantum Chromodynamics (QCD) is collinear factorization
Collins:1989gx , which allows for the calculation of hadronic cross sections
by separating the effects of long- and short-distance physics. The long-
distance effects are described by universal, process-independent, parton
distribution functions (PDFs) Kovarik:2019xvh ; Ethier:2020way ;
Klasen:2023uqj  and fragmentation functions (FFs) Metz:2016swz . While non-
perturbative and thus not calculable from perturbative QCD, their scale
dependence is still predicted by the perturbative Dokshitzer-Gribov-Lipatov-
Altarelli-Parisi (DGLAP) evolution equations Gribov:1972ri ; Gribov:1972rt ;
Dokshitzer:1977sg ; Altarelli:1977zs . The short-distance effects, on the
other hand, are perturbatively calculable due to the asymptotic freedom of
QCD.

One of the most poorly known components in the global analysis of free- and
bound-proton (nuclear) PDFs is the strange-quark distribution Hou:2019efy ;
Bailey:2020ooq ; NNPDF:2021njg ; Eskola:2021nhw ; Duwentaster:2022kpv ;
AbdulKhalek:2022fyi . This in part limits e.g. the precision determination of
Standard Model parameters — like the mass of the $W^{\pm}$ boson ATLAS:2017rzl
; Bagnaschi:2019mzi ; Hussein:2019kqx ; ATLAS:2023fsi  or the weak mixing
angle CMS:2018ktx  — from the data collected at the Large Hadron Collider
(LHC). An important process to constrain the strange-quark content of nucleons
is the semi-inclusive dimuon production in charged-current neutrino-nucleus
($\nu_{\mu}N$) deep-inelastic scattering (DIS) CCFR:1994ikl ; NuTeV:2001dfo ;
NuTeV:2007uwm ; CHORUS:2008vjb ; NOMAD:2013hbk . In comparison to the fully
inclusive neutrino DIS, whose compatibility with the neutral-current DIS in
the case of nuclear targets has been discussed quite thoroughly
Schienbein:2009kk ; Paukkunen:2010hb ; Kovarik:2010uv ; Paukkunen:2013grz ;
Muzakka:2022wey , this process produces a charm quark in the partonic final
state, which then hadronizes into a charmed hadron and eventually decays into
a final state containing a secondary decay muon. Due to the elements of the
Cabibbo–Kobayashi–Maskawa (CKM) matrix, $\left|V_{cs}\right|\approx 0.973$,
$\left|V_{cd}\right|\approx 0.225$, and $\left|V_{cb}\right|\approx 0.042$
Workman:2022ynf , the outgoing charm quark couples mostly to the strange-quark
content of the target. While the neutrino data have been taken with heavy
nuclear targets to obtain sufficient statistics, they are still typically used
in global fits of free-proton PDFs Hou:2019efy ; Bailey:2020ooq ;
NNPDF:2021njg . This is taken into account by either “correcting” the data
Hou:2019efy ; Accardi:2021ysh  for the nuclear effects or treating them as
theoretical uncertainties Ball:2018twp ; Bailey:2020ooq ; NNPDF:2021njg  tied
to the corrections obtained from nuclear PDFs. Ultimately, this creates a
correlation between proton and nuclear PDFs. In the future, experiments
employing far-forward neutrinos originating from proton-proton collisions at
the LHC could provide further constraints on the strange-quark distribution
Cruz-Martinez:2023sdv . Indeed, the first collider neutrinos have already been
detected by the FASER FASER:2023zcr  and SND@LHC SNDLHC:2023pun
collaborations.

By including LHC proton-proton ($pp$) data into a global fit of PDFs, one can
hope to reduce the dependence of the strange-quark PDFs on the treatment of
nuclear corrections. While the inclusive $Z$ and $W^{\pm}$ boson production
processes carry a significant contribution from channels involving a strange
quark Kusina:2012vh , they are still subleading in comparison to the
contributions from other production channels involving up and down quarks.
Moreover, at the high interaction scale $Q^{2}\sim 10^{4}\,{\rm GeV^{2}}$
relevant for these processes, a big part of the behaviour of the strange-quark
PDFs is simply dictated by QCD with gluons splitting into $s\bar{s}$-pairs. As
a result, it is not particularly easy to faithfully resolve the non-
perturbative strange-quark density from such inclusive measurements
Eskola:2022rlm  and e.g. any deficiencies in the theoretical description can
manifest themselves as shifts in the strange-quark PDF just because it is much
less constrained than e.g. the up- and down-quark densities Kusina:2020lyz .
The background from other partonic processes can be suppressed by considering
the production of $W^{\pm}$ bosons associated with a charmed particle
Baur:1993zd ; Stirling:2012vh ; ATLAS:2014jkm ; CMS:2013wql ; CMS:2018dxg ;
ATLAS:2023ibp . However, the interaction scale is still high, and the charm
identification increases the experimental uncertainties. Interestingly, the
LHC $pp$ data have displayed some tensions with the neutrino dimuon data, the
former preferring a larger ratio $R_{s}=(s+\bar{s})/(\bar{u}+\bar{d})$ at
$x\approx 0.02$ ATLAS:2012sjl ; ATLAS:2016nqi ; Cooper-Sarkar:2018ufj ;
NNPDF:2017mvq ; Hou:2019efy ; Bailey:2020ooq . However, higher-order QCD
corrections on the dimuon cross sections tend to be negative Berger:2016inr ;
Gao:2017kkx  and appear to somewhat ease the observed tensions Faura:2020oom .
The strange-quark PDF can also be constrained by semi-inclusive kaon
production in DIS HERMES:2013ztj ; Borsa:2017vwy ; Sato:2019yez  on protons,
although in this case the uncertainties in the parton-to-kaon FFs limit the
precision. As a result, semi-inclusive dimuon production still remains an
integral part in constraining the stange-quark PDFs of the proton.

The semi-inclusive dimuon production cross section in neutrino-nucleus DIS is
typically computed by assuming that it is proportional to the inclusive charm
production cross section

$\mathrm{d}\sigma(\nu_{\mu}N\to\mu^{-}\mu^{+}X)\simeq\mathcal{A}\mathcal{B}_{\mu}\mathrm{d}\sigma(\nu_{\mu}N\to\mu^{-}cX),$
(1)

where the acceptance correction $\mathcal{A}$ compensates for the experimental
energy cut on the decay muon originating from the process and
$\mathcal{B}_{\mu}$ is the semileptonic branching ratio averaged over charm
mesons. At leading-order (LO) fixed-flavour number scheme (no heavy-quark
PDFs), the dimuon cross section indeed factorizes as above, but at higher
perturbative orders or in variable flavour-number scheme this is no longer the
case. This poses a potential problem for increasing the accuracy of the dimuon
production calculation as increasing the perturbative order of the inclusive
charm production cross section does not necessarily translate into an
increased accuracy of the entire dimuon process if the acceptance correction
$\mathcal{A}$ is not updated likewise. Indeed, the acceptance correction
$\mathcal{A}$ is typically calculated in a different framework such as the
next-to-leading order (NLO) DISCO Monte Carlo program Kretzer:2001tc . In this
paper, we forgo the assumption in eq. (1) and instead compute the dimuon cross
section directly by implementing the entire chain of semi-inclusive charmed
hadron production and the subsequent muonic decay. The calculation of semi-
inclusive charmed hadron production effectively resums all logarithms
$\log(Q^{2}/m^{2}_{\rm charm})$ that originate from collinear initial- and
final-state QCD radiation into the scale dependence of PDFs and FFs. At this
stage our implementation will be accurate at the next-to-leading order (NLO)
level in perturbative QCD, accounting for the leading charm-mass corrections,
though an extension to the full general-mass variable flavour number scheme
(GM-VFNS) Kretzer:1997pd  will be straightforward. In addition, by using the
recent next-to-NLO (NNLO) results for semi-inclusive DIS (SIDIS) publicly
available from the literature Abele:2021nyo , we are already able to estimate
the expected size of NNLO corrections in our framework. The central idea is,
that in the approach we present, one can systematically improve the
theoretical description without a need for an external acceptance correction.

The paper is organized as follows. In section 2, we introduce the semi-
inclusive calculation of charmed-hadron production in neutrino-nucleus DIS and
in section 3 we discuss in detail the implementation of the charmed-hadron
decay. In section 4, we then present our results for the dimuon cross sections
comparing them with the NuTeV NuTeV:2007uwm  and CCFR CCFR:1994ikl  data,
assess the uncertainties of our approach, and calculate the effective
acceptance and nuclear corrections. Finally, section 5 summarizes the work and
discusses further developments.

##  2 Semi-inclusive charged-current deep inelastic scattering

In the usual inclusive DIS processes, only the momentum of the outgoing lepton
is measured and one integrates over the momenta of all other final-state
particles. While experimentally simpler, this approach gives no information
about the outgoing particles in the hadronic final state. In SIDIS, however,
also the momentum of a particle in the hadronic final state is measured.

Figure 1: Diagrammatic depiction of the (anti)neutrino-SIDIS process
$\nu_{\mu}N\to\mu hX$. Momenta of the particles are indicated in parentheses.

The charged-current SIDIS cross section for the production of a hadron $h$,
depicted in figure 1, is given by

$\begin{split}\frac{\mathrm{d}\sigma(\nu_{\mu}N\to\mu
hX)}{\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z}&=\frac{G_{F}^{2}M_{W}^{4}}{\left(Q^{2}+M_{W}^{2}\right)^{2}}\frac{Q^{2}}{2\pi
xy}\\\ &\phantom{=}\
\times\bigg{[}xy^{2}F_{1}(x,z,Q^{2})+\left(1-y-\frac{xyM^{2}}{s-M^{2}}\right)F_{2}(x,z,Q^{2})\\\
&\phantom{=\times\bigg{[}}\ \pm
xy\left(1-\frac{y}{2}\right)F_{3}(x,z,Q^{2})\bigg{]},\end{split}$ (2)

where $G_{F}\approx 1.166379\times 10^{-5}\,\mathrm{GeV}^{-2}$ is the Fermi
coupling constant, $M_{W}\approx 80.377\,\mathrm{GeV}$ is the mass of the
$W^{\pm}$ boson, and the $+$ sign in the third term corresponds to neutrino
scattering and the $-$ sign to antineutrino scattering. In eq. (2) the usual
DIS kinematical variables

$\displaystyle Q^{2}$ $\displaystyle=-q^{2}=(k-k^{\prime})^{2}\,,$ (3)
$\displaystyle x$ $\displaystyle=\frac{Q^{2}}{2P_{N}\cdot q}\,,$
$\displaystyle y$ $\displaystyle=\frac{P_{N}\cdot q}{P_{N}\cdot k}\,,$

are accompanied by the momentum fraction of the produced hadron

$z=\frac{P_{N}\cdot P_{h}}{P_{N}\cdot q}.$ (4)

Here $P_{N}$, $P_{h}$, $k$, and $k^{\prime}$ are the momenta of the target
nucleon $N$, outgoing hadron $h$, incoming neutrino, and outgoing muon. At
NLO, the zero-mass structure functions are given by

$\displaystyle F_{2}(x,z,Q^{2})$
$\displaystyle=2x\sum_{ff^{\prime}}\lambda_{+}^{ff^{\prime}}\int_{x}^{1}\frac{\mathrm{d}\xi}{\xi}\int_{z}^{1}\frac{\mathrm{d}\xi^{\prime}}{\xi^{\prime}}f(\xi,\mu_{\text{fact}}^{2})D_{f^{\prime}\to
h}(\xi^{\prime},\mu_{\text{frag}}^{2})$ (5)
$\displaystyle\quad\quad\times\left[C_{\mathrm{LO}}^{ff^{\prime}}(\hat{x},\hat{z})+\frac{\alpha_{s}(\mu_{\text{ren}}^{2})}{2\pi}C_{2}^{ff^{\prime}}(\hat{x},\hat{z};\mu_{\text{fact}}^{2},\mu_{\text{frag}}^{2})\right],$
$\displaystyle F_{1}(x,z,Q^{2})$
$\displaystyle=\sum_{ff^{\prime}}\lambda_{+}^{ff^{\prime}}\int_{x}^{1}\frac{\mathrm{d}\xi}{\xi}\int_{z}^{1}\frac{\mathrm{d}\xi^{\prime}}{\xi^{\prime}}f(\xi,\mu_{\text{fact}}^{2})D_{f^{\prime}\to
h}(\xi^{\prime},\mu_{\text{frag}}^{2})$
$\displaystyle\quad\quad\times\left[C_{\mathrm{LO}}^{ff^{\prime}}(\hat{x},\hat{z})+\frac{\alpha_{s}(\mu_{\text{ren}}^{2})}{2\pi}C_{1}^{ff^{\prime}}(\hat{x},\hat{z};\mu_{\text{fact}}^{2},\mu_{\text{frag}}^{2})\right],$
$\displaystyle F_{3}(x,z,Q^{2})$
$\displaystyle=2\sum_{ff^{\prime}}\lambda_{-}^{ff^{\prime}}\int_{x}^{1}\frac{\mathrm{d}\xi}{\xi}\int_{z}^{1}\frac{\mathrm{d}\xi^{\prime}}{\xi^{\prime}}f(\xi,\mu_{\text{fact}}^{2})D_{f^{\prime}\to
h}(\xi^{\prime},\mu_{\text{frag}}^{2})$
$\displaystyle\quad\quad\times\left[C_{\mathrm{LO}}^{ff^{\prime}}(\hat{x},\hat{z})+\frac{\alpha_{s}(\mu_{\text{ren}}^{2})}{2\pi}C_{3}^{ff^{\prime}}(\hat{x},\hat{z};\mu_{\text{fact}}^{2},\mu_{\text{frag}}^{2})\right],$

where $\hat{x}=x/\xi$ and $\hat{z}=z/\xi^{\prime}$. The LO coefficient
function is
$C_{\mathrm{LO}}^{ff^{\prime}}(\hat{x},\hat{z})=\delta(1-\hat{x})\delta(1-\hat{z})$
when $f$ and $f^{\prime}$ are quarks, and zero otherwise. The couplings are
given by

$\lambda_{+}^{ff^{\prime}}=\begin{cases}\left|V_{ff^{\prime}}\right|^{2}&\text{$f$
and $f^{\prime}$ are quarks}\\\ 1&\text{$f$ or $f^{\prime}$ is a
gluon}\end{cases},$ (6)

and

$\lambda_{-}^{ff^{\prime}}=\begin{cases}-\lambda_{+}^{ff^{\prime}}&\text{$f$
or $f^{\prime}$ is an antiquark}\\\
\lambda_{+}^{ff^{\prime}}&\text{otherwise}\end{cases}.$ (7)

The PDFs $f$ are evaluated at the factorization scale $\mu_{\text{fact}}$, FFs
$D_{f\rightarrow h}$ at the fragmentation scale $\mu_{\text{frag}}$, and the
strong coupling constant $\alpha_{s}$ at the renormalization scale
$\mu_{\text{ren}}$. The quark flavors are determined by charge conservation.
For the CKM matrix elements $V_{ff^{\prime}}$, we use the global-fit values
Workman:2022ynf . The massless $\overline{\rm MS}$ NLO coefficient functions
$C_{i}^{ff^{\prime}}$ for $i=1,2,3$ can be found from refs. Furmanski:1981cw
and deFlorian:1997zj , where the latter also explicitly indicates the scale-
dependent logarithms.111Note that there is a typographical error in ref.
Furmanski:1981cw : on the second line of the coefficient $C_{2}^{F(1)}$ in
Appendix II, the term $\frac{3(1+x)^{2}}{1-x}\ln x$ should be
$\frac{3(1+x^{2})}{1-x}\ln x$. The NNLO coefficient functions are currently
known only approximately. The massless quark-to-quark coefficient functions
are known in the next-to-leading power (NLP) approximation to NNLO
Abele:2021nyo  and even to $\mathrm{N}^{3}\mathrm{LO}$ accuracy Abele:2022wuy
. Although these NNLO coefficient functions are approximate, we can still use
them to get a rough estimate of the NNLO corrections.222At the final stages of
preparing this paper, preprints of the full massless NNLO coefficient
functions for unpolarized and polarized SIDIS have appeared Bonino:2024qbh ;
Bonino:2024wgg ; Goyal:2024tmo . These corrections are given for photon-
mediated interactions, so one has to adapt the electromagnetic couplings to
the charged-current case.

While quark masses are neglected in the hard coefficient functions
$C_{i}^{ff^{\prime}}$, we can take most of their effects into account by
replacing $x$ in the right-hand side of eq. (5) by the so-called slow-
rescaling variable PhysRevLett.36.1163 ; Barone:1995ma

$\chi_{q}\equiv x\left(1+\frac{m_{q}^{2}}{Q^{2}}\right)>x\,,$ (8)

where $m_{q}$ is the quark mass in channels where a single heavy quark is
produced. If several heavy quarks are produced in the partonic scattering,
this can be generalized. In practice, the only relevant channels turn out to
be those in which a single charm quark is produced. The charm mass is taken to
be that used in the PDFs. While the slow-rescaling variable induces an
enhancement to $F_{2}$, as this structure function is directly proportional to
$x$, the net effect also depends on the $x$-dependence of the PDF due to the
shift in the momentum argument at which the PDFs are evaluated.

We note that introducing the rescaling variable $\chi_{q}$ is not just a
prescription but it emerges naturally from the GM-VFNS calculations
Nadolsky:2009ge . Indeed, at leading order, the slow-rescaling variable
actually captures all quark-mass effects. Thus any additional mass corrections
are also suppressed by $\alpha_{s}$ and are, in general, of the order of
$\mathcal{O}(\alpha_{s}m^{2}/Q^{2})$. This is a specific feature of charged-
current DIS, as there are additional mass effects already at the leading order
in the neutral-current case. For a more quantitative check, we have computed
the ratio of inclusive (i.e. not SIDIS) charm production with $m_{c}=1.5\,{\rm
GeV}$ and $m_{c}=1.3\,{\rm GeV}$. This ratio has been computed also in ref.
Gao:2017kkx , including full quark-mass corrections in the fixed flavour-
number scheme. Our calculation, with just the slow-rescaling variable, agrees
to within ${5}{\,\%}$ with ref. Gao:2017kkx  indicating that the rescaling
variable already accounts for the majority of quark-mass effects and is an
accurate-enough approach for our present purposes.

##  3 Decay of charmed hadrons

While the charmed hadrons cannot always be directly reconstructed by the
experiments, their specific decay products can be measured. One of the
prominent decay channels for charmed hadrons is $\mu+\text{anything}$, which
produces the secondary decay muon in the dimuon production process. This
process has been measured most recently by NuTeV NuTeV:2007uwm , CCFR
CCFR:1994ikl , and NOMAD NOMAD:2013hbk .

Figure 2: Diagrammatic depiction of the production process $AB\to CX$. Momenta
of the particles are given in parentheses.

###  3.1 General formalism and the decay function

We begin by assuming that the cross section for the production of a charmed
hadron $h$ in collisions between particles $A$ and $B$ with center-of-mass
energy $\sqrt{s}$, followed by its decay to particle $C$, factorizes into a
Breit-Wigner form,

$\begin{split}\sigma(AB\to
CX)=\frac{1}{2s}\int&\mathrm{d}\Pi(P_{C})\,\mathrm{d}\Pi(P_{X_{1}})\,\mathrm{d}\Pi(P_{X_{2}})\\\
&\times(2\pi)^{4}\delta^{(4)}(P_{A}+P_{B}-P_{C}-P_{X_{1}}-P_{X_{2}})\\\
&\times\mathcal{A}_{p}\frac{1}{(P_{h}^{2}-m_{h}^{2})^{2}+m_{h}^{2}\Gamma_{\mathrm{tot}}^{2}}\mathcal{A}_{d}\,,\end{split}$
(9)

see figure 2. Here, $\mathrm{d}\Pi$ denote the phase space elements, $P_{i}$
are the momenta of different (groups of) particles, $m_{h}$ is the mass of the
hadron $h$, and $\Gamma_{\mathrm{tot}}$ denotes its total decay width. The
symbols $\mathcal{A}_{p}$ and $\mathcal{A}_{d}$ denote the parts associated
with the production and decay of the particle $C$, respectively. Here, we
neglect the possible spin correlations between the production and decay, if
the hadron $h$ is not a spin-$0$ particle. In the narrow-width approximation,

$\frac{1}{(P_{h}^{2}-m_{h}^{2})^{2}+m_{h}^{2}\Gamma_{\mathrm{tot}}^{2}}\simeq\frac{\pi}{m_{h}\Gamma_{\mathrm{tot}}}\delta(P_{h}^{2}-m_{h}^{2}).$
(10)

In order to make the momentum $P_{h}$ of the intermediate-state particle
explicit, we insert unity as

$1=\int\frac{\mathrm{d}P_{h}^{2}}{2\pi}\int\frac{\mathrm{d}^{3}\mathbf{P}_{h}}{(2\pi)^{3}2E_{h}}(2\pi)^{4}\delta^{(4)}(P_{h}-P_{C}-P_{X_{2}})\,,$
(11)

where $E_{h}$ is the energy of the hadron $h$. This results in

$\begin{split}\sigma(AB\to
CX)=\frac{1}{2s}&\int\frac{\mathrm{d}^{3}\mathbf{P}_{h}}{(2\pi)^{3}2E_{h}}\\\
&\times\int\mathrm{d}\Pi(P_{X_{1}})\,\mathcal{A}_{p}(2\pi)^{4}\delta^{(4)}(P_{A}+P_{B}-P_{X_{1}}-P_{h})\\\
&\times\frac{1}{2m_{h}\Gamma_{\mathrm{tot}}}\int\mathrm{d}\Pi(P_{X_{2}})\,\mathrm{d}\Pi(P_{C})\mathcal{A}_{d}(2\pi)^{4}\delta^{(4)}(P_{h}-P_{C}-P_{X_{2}}){\big{|}_{P_{h}^{2}=m_{h}^{2}}}\,.\end{split}$
(12)

Identifying the differential production cross section with

$\frac{\mathrm{d}\sigma(AB\to
hX)}{\mathrm{d}\Pi(P_{h})}=\frac{1}{2s}\int\mathrm{d}\Pi(P_{X_{1}})\,\mathcal{A}_{p}(2\pi)^{4}\delta^{(4)}(P_{A}+P_{B}-P_{X_{1}}-P_{h})\,,$
(13)

and the decay width of the hadron

$\Gamma_{h\to
C}=\frac{1}{2m_{h}}\int\mathrm{d}\Pi(P_{X_{2}})\,\mathrm{d}\Pi(P_{C})\mathcal{A}_{d}(2\pi)^{4}\delta^{(4)}(P_{h}-P_{C}-P_{X_{2}}){\big{|}_{P_{h}^{2}=m_{h}^{2}}}\,,$
(14)

the total cross section becomes

$\sigma(AB\to CX)=\int\mathrm{d}\Pi(P_{h})\frac{\mathrm{d}\sigma(AB\to
hX)}{\mathrm{d}\Pi(P_{h})}\frac{\Gamma_{h\to C}}{\Gamma_{\mathrm{tot}}}\,.$
(15)

We now want to introduce a kinematical cut on the energy of the outgoing
particle $C$, ${E_{C}\geq E_{C}^{\mathrm{min}}}$. This is needed for dimuon
production, as the NuTeV and CCFR data include a cut on the decay muon energy.
To take such a cut into account, we have to make eq. (15) differential in the
momentum of $C$. This amounts to replacing the decay width $\Gamma_{h\to C}$
with a more differential quantity. To this end, we define a dimensionless
decay function $d_{h\to C}$,

$d_{h\to
C}(w)=\frac{1}{2(2\pi)^{3}}\int\mathrm{d}\Pi(P_{X_{2}})\mathcal{A}_{d}(2\pi)^{4}\delta^{(4)}(P_{h}-P_{C}-P_{X_{2}}){\big{|}_{P_{h}^{2}=m_{h}^{2}}}\,,$
(16)

where

$\displaystyle w\equiv\frac{P_{C}\cdot P_{h}}{m_{h}^{2}}\,,\quad 0\leq
w\leq\frac{1}{2}\,.$ (17)

To see why $d_{h\to C}$ depends only on $w$, note that $d_{h\to C}$ is a
Lorentz scalar depending on the momenta $P_{h}$ and $P_{C}$ and thus the only
non-trivial invariant is $P_{C}\cdot P_{h}$. As the mass of $C$ will be
negligible in the present case, it makes sense to normalize this with
$P_{h}^{2}=m_{h}^{2}$. Using eq. (16), the decay width in eq. (14) can be
written as

$\Gamma_{h\to
C}=\frac{1}{2m_{h}}\int\frac{\mathrm{d}^{3}\mathbf{P}_{C}}{E_{C}}d_{h\to
C}(w)=\frac{1}{2m_{h}}\int\mathrm{d}|\mathbf{P}_{C}|\,\mathrm{d}\Omega\,\frac{|\mathbf{P}_{C}|^{2}}{E_{C}}
d_{h\to C}(w)\,.$ (18)

The differential decay width is then

$\frac{\mathrm{d}\Gamma_{h\to
C}}{\mathrm{d}|\mathbf{P}_{C}|}=\frac{\pi}{m_{h}}\frac{|\mathbf{P}_{C}|^{2}}{E_{C}}\int\mathrm{d}(\cos\theta)\,d_{h\to
C}(w)\,,$ (19)

where the angular dependence in $w$ is given explicitly by

$w=\frac{1}{m_{h}^{2}}\left(E_{C}E_{h}-|\mathbf{P}_{C}||\mathbf{P}_{h}|\cos\theta\right)\,,$
(20)

and $\theta$ denotes the relative angle between ${\bf P}_{C}$ and ${\bf
P}_{h}$. Integrating eq. (19) over $|\mathbf{P}_{C}|$ and inserting it back to
eq. (15),

$\mathrm{d}\sigma(AB\to
CX)=\frac{\pi}{m_{h}\Gamma_{\mathrm{tot}}}\int\mathrm{d}\Pi(P_{h})\,\frac{\mathrm{d}\sigma(AB\to
hX)}{\mathrm{d}\Pi(P_{h})}\int\mathrm{d}|\mathbf{P}_{C}|\,\frac{|\mathbf{P}_{C}|^{2}}{E_{C}}\int\mathrm{d}(\cos\theta)\,d_{h\to
C}(w)\,.$ (21)

For convenience, we define $\rho\equiv E_{C}/E_{h}$ and, neglecting the mass
of $C$, introduce an energy-dependent branching fraction $B_{h\to
C}(E_{h},E_{C}^{\mathrm{min}})$ as

$B_{h\to
C}(E_{h},E_{C}^{\mathrm{min}})\equiv\frac{\pi}{m_{h}\Gamma_{\mathrm{tot}}}\int\mathrm{d}\rho\,\rho
E_{h}^{2}\int\mathrm{d}(\cos\theta)\,d_{h\to C}(w){\big{|}_{E_{C}=\rho
E_{h}\geq E_{C}^{\mathrm{min}}}}\,,$ (22)

which allows us to write the cross section simply as

$\mathrm{d}\sigma(AB\to
CX)=\int\mathrm{d}\Pi(P_{h})\,\frac{\mathrm{d}\sigma(AB\to
hX)}{\Pi(P_{h})}B_{h\to C}(E_{h},E_{C}^{\mathrm{min}}).$ (23)

Combining eqs. (19) and (22), we see that

$B_{h\to
C}(E_{h},E_{C}^{\mathrm{min}})=\frac{1}{\Gamma_{\mathrm{tot}}}\int\mathrm{d}^{3}\mathbf{P}_{C}\frac{\mathrm{d}\Gamma_{h\to
C}}{\mathrm{d}^{3}\mathbf{P}_{C}}\bigg{|}_{E_{C}\geq E_{C}^{\mathrm{min}}}.$
(24)

When there is no cut, $E_{C}^{\mathrm{min}}=0$, $B_{h\to C}$ becomes
independent of $E_{h}$ and reduces to the usual branching fraction
$\Gamma_{h\to C}/\Gamma_{\mathrm{tot}}$. For a non-zero energy cut, $B_{h\to
C}$ depends on the energy $E_{h}$. Its dependence on $E_{h}$ can be
precomputed in a dense-enough grid, such that the double integrals in eq. (22)
do not need to be evaluated repeatedly. This significantly speeds up the
evaluation of the cross section in eq. (23).

We will apply the formalism introduced above to the case of dimuon production
process $\nu_{\mu}N\to\mu^{-}\mu^{+}X$. Using eq. (23), we can write the
production cross section as

$\frac{\mathrm{d}\sigma(\nu_{\mu}N\to\mu^{-}\mu^{+}X)}{\mathrm{d}x\,\mathrm{d}y}=\sum_{h}\int_{z_{\text{min}}}^{1}\mathrm{d}z\,\frac{\mathrm{d}\sigma(\nu_{\mu}N\to\mu^{-}hX)}{\mathrm{d}z\,\mathrm{d}x\,\mathrm{d}y}B_{h\to\mu}(E_{h}=zyE_{\nu},E_{\mu}^{\text{min}}),$
(25)

where the sum is over
$h\in\left\\{D^{0},D^{+},D_{s},\Lambda^{+}_{c}\right\\}$. The kinematical
variables $x$, $y$, and $z$ were defined in eqs. (3) and (4). The lower bound
of integration $z_{\mathrm{min}}$ is given by

$z=\frac{E_{h}}{yE_{\nu}}\geq\frac{\max\\{E_{\mu}^{\mathrm{min}},m_{h}\\}}{yE_{\nu}}\equiv
z_{\text{min}}\,,$ (26)

where the relation $E_{h}=zyE_{\nu}$ holds in the rest frame of the nuclear
target.

###  3.2 Fitting the decay function

For the dimuon production process, we need to determine $d_{h\to\mu}$. To do
this, we adapt the formalism introduced in the previous section to $D$-meson
production in $e^{+}e^{-}$ collisions. We use data from CLEO CLEO:2006ivk
given for the semileptonic decays $D\to e^{+}\nu_{e}X$. While this technically
corresponds to $C=e^{+}$ and not $C=\mu^{+}$, we assume that the differences
are negligible as both the electron and muon masses are small compared to the
$D$-meson masses.

The CLEO data are collected at the $\psi(3770)$ mass threshold. The
$\psi(3770)$ decays mostly (with a branching fraction of
$\left(93\begin{subarray}{c}+8\\\ -9\end{subarray}\right)\,\%$ Workman:2022ynf
) into a $D$-meson pair $D\overline{D}$, where the relevant $D^{0}$ and
$D^{+}$ mesons are not produced to rest. Instead, they carry momenta
${0.277}{\,{\rm GeV}}$ and ${0.243}{\,{\rm GeV}}$ for the $D^{0}$ and $D^{+}$
meson, respectively. Although the CLEO data are given separately for the
$D^{0}$ and $D^{+}$ mesons, we combine them into one universal decay function
that we use for all charmed hadrons. This is a reasonable assumption since the
decay function depends only on $w$, which by eq. (20) depends on $m_{h}$ and
$\left|\mathbf{P}_{h}\right|$. Both of these values are very similar for
$D^{0}$ and $D^{+}$.

Figure 3: The branching fractions $B_{h\to\mu}$ of eq. (22) as a function of
the hadron energy $E_{h}$ with $E_{\mu}^{\mathrm{min}}={5}{\,{\rm GeV}}$ for
$D^{+}$, $D^{0}$, $D_{s}$, and $\Lambda^{+}_{c}$. For the decay function, we
use eq. (27) and parameters given in eq. (28).

The data are fitted with the parametrization

$d_{h\to\mu}(w)=Nw^{\alpha}(1-\gamma w)^{\beta}\,,$ (27)

for $0\leq w\leq 1/\gamma$ and zero elsewhere. The parameter $\gamma$ is
expected to be around 2, since $w\leq 1/2$. For the fit, we combine the two
CLEO data sets for $D^{0}$ and $D^{+}$, and take the $D$-meson mass and
momentum parameters to be the mean values of the corresponding $D^{0}$ and
$D^{+}$ values: $m_{h}={1.867}{\,{\rm GeV}}$ and
$|\mathbf{P}_{h}|={0.26}{\,{\rm GeV}}$. The best-fit parameters are then

$N=3.156,\quad\alpha=1.009,\quad\beta=1.799,\quad\gamma=2.101,$ (28)

which are obtained by the nonlinear least-squares method. Figure 3 shows the
branching fractions $B_{h\to\mu}$ of eq. (22) using the parametrization eq.
(27) and parameter values in eq. (28). Here, and in the rest of the paper, we
have assumed that the same decay function $d_{h\to\mu}$ applies for all
charmed hadrons. The asymptotic values as $E_{h}\to\infty$ are consistent with
the usual branching fractions Workman:2022ynf  which gives us confidence
concerning this universality assumption of $d_{h\to\mu}$.

Figure 4: Left: The decay function $d_{h\to\mu}$ of eq. (27) as a function of
$w$. The central curve corresponds to the best-fit parameters given in eq.
(28), while the replica curves correspond to all of the replica fit
parameters. Right: Differential decay width of $D\to e^{+}\nu_{e}X$ as a
function of the $D$ meson momentum, with $D^{0}$ and $D^{+}$ data from ref.
CLEO:2006ivk . The $D^{+}$ data points have been slightly shifted horizontally
for improved readibility. The central fit curve corresponds to eq. (19) with
the parameter values given in eq. (28). The band depicts the ${90}{\,\%}$
confidence-level uncertainty taken from the replica fits.

The fit parameters are highly correlated. To assess the uncertainty of the fit
parameters, we employ a strategy similar to that of refs. Forte:2002fg ;
Ball:2008by . From the original CLEO data set, we generate 1000 new replica
data sets so that each data value $D_{i}$ is transformed into

$D_{i}^{\prime}=D_{i}(1+\delta D_{i}R_{i})\,,$

where $R_{i}$ is a random number drawn from the standard normal distribution
and $\delta D_{i}$ is the total relative uncertainty of the data point.333A
more appropriate choice would be to only include statistical uncertainty, but
the CLEO data includes only the total uncertainty. We then fit these replica
data sets with the same procedure as above to obtain a set of replica fits. An
observable can then be evaluated using all replica fits to estimate the
uncertainty resulting from the fitting procedure. In order to combine the
results obtained with the replica fits into an uncertainty band, we can either
take some central subset of all replica results as in ref. AbdulKhalek:2022fyi
or simply compute the standard deviation. We have found that in the present
case the differences between these two approaches are negligible. The left-
hand panel of figure 4 shows the central result for $d_{h\to\mu}$ accompanied
by the 1000 replica fits, and the right-hand panel compares the fit and its 90
% uncertainty with the CLEO data. This $90\,\%$ uncertainty band is derived by
rejecting $5\,\%$ of the upper- and lower-extreme replica sets.

##  4 Results

Having the decay function for $D$ mesons fitted to the CLEO data we can use it
in eq. (25) and calculate dimuon production cross sections in neutrino-nucleus
SIDIS without resorting to the approximate factorization of the acceptance
correction discussed in the context of eq. (1). We will compare these
calculations with data from the NuTeV NuTeV:2007uwm  and CCFR CCFR:1994ikl
experiments. These data correspond to $E_{\mu}^{\rm min}=5\,{\rm GeV}$.
Similar dimuon data exist also from the NOMAD collaboration NOMAD:2013hbk . In
this case the data are provided as a ratio between the fully integrated dimuon
and inclusive DIS processes which is beneficial in terms of cancellation of
experimental systematic uncertainties. However, the minimum allowed $Q^{2}$ is
very low, $Q^{2}>1\,{\rm GeV}^{2}$, and the observable is sensitive to the
region which begins to be outside the regime where our calculation can be
applied. In addition, an acceptance correction has already been applied to
account for the finite experimental cut in the final-state muon energy
$E_{\mu}^{\rm min}$, whereas the entire point of the present paper is to avoid
such external correction. While we could technically set $E_{\mu}^{\rm
min}=0$, this would make us also increasingly sensitive to the small-$z$
instabilities of the time-like $Q^{2}$ evolution of the FFs, and also to the
region of $D\rightarrow e^{+}X$ transition not directly constrained by the
CLEO data. As a result, we only consider here the case with a finite
$E_{\mu}^{\rm min}$.

For the nuclear PDFs, we use the EPPS21 Eskola:2021nhw , nCTEQ15HQ
Duwentaster:2022kpv , and nNNPDF3.0 AbdulKhalek:2022fyi  NLO sets. The charm
mass, which enters the cross-section calculation through the slow-rescaling
variable, is taken to be the same as in the PDF set, i.e.
$1.3\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ for EPPS21 and nCTEQ15HQ, and
$1.51\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ for nNNPDF3.0. For the
fragmentation functions, we use the NLO set kkks08 Kneesch:2007ey  for $D^{0}$
and $D^{+}$, and NLO bkk05 Kniehl:2006mw  for $D_{s}$ and $\Lambda_{c}^{+}$.
Both sets use a value of $1.5\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ for the
charm-quark mass, which is also the parametrization scale in these analyses.
In the case of kkks08, the provided sets include a fit using $e^{+}e^{-}$
collision data only from the OPAL/LEP OPAL:1996ikk ; OPAL:1997edj  experiment
at $\sqrt{s}\approx m_{Z}\approx$91\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$,
one with data only from the Belle Belle:2005mtx  experiment at
$\sqrt{s}\approx$10\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$, and a global fit
including data from several experiments. As was pointed out in ref.
Kneesch:2007ey , some tension was found between the low- and high-energy data
sets, which did not allow for a satisfactory fit when combining these
different data sets. In this study we choose to use only the OPAL data, as
this should be theoretically cleanest as potential ambiguities due to mass
effects ought to be negligible at these high energies.

###  4.1 Uncertainties from the decay function

We begin by quantifying the uncertainties associated with the decay function
described in section 3.2. Figure 5 indicates that the cross sections are very
stable in terms of variations of the parameters and the uncertainty associated
with the fitting procedure is minimal with a relative uncertainty of
${2}{\,\%}$, which is of the same order as the uncertainty in the right-hand
panel of figure 4.

Figure 5: Left: Neutrino dimuon cross section evaluated using EPPS21
Eskola:2021nhw  at $y=0.771$ and $E_{\nu}={247}{\,{\rm GeV}}$. The central
curve uses the parameter values given in eq. (28). The band depicts the
${90}{\,\%}$ confidence-level uncertainty taken from the replica fits. The
theoretical computation is compared against NuTeV data from ref. NuTeV:2007uwm
. The cross section values should be multiplied by $G_{F}^{2}ME_{\nu}/100\pi$.
Right: Same as on the left, but normalized to the central result with EPPS21.

###  4.2 Channel decomposition

Figure 6:  Contributions of different partonic channels to the neutrino (left)
and antineutrino (right) dimuon cross sections. Only the relevant channels are
shown, as the other channels are negligible. The cross section is computed
with EPPS21 Eskola:2021nhw  in the NuTeV kinematics. The cross section values
should be multiplied by $G_{F}^{2}ME_{\nu}/100\pi$.

Figure 7: Left: Contributions of different charmed hadrons to the neutrino
dimuon cross section, computed with EPPS21 Eskola:2021nhw  at $y=0.558$ and
$E_{\nu}=$174.29\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$. The contribution of
$\Lambda_{c}^{+}$ is multiplied by $1.14$. The cross section values should be
multiplied by $G_{F}^{2}ME_{\nu}/100\pi$. Right: The same as on the left, but
as ratios to the total cross section.

Next, we study the sensitivity of the cross sections to different partonic
production channels by plotting some contributions separately in figure 6.
Only those contributions in which there is a final-state charm quark present
are shown, as the contributions from gluons and other quarks fragmenting into
a charmed hadron turned out negligible. As expected from the CKM matrix, the
most relevant channel is the $s\to c$ channel ($\bar{s}\to\bar{c}$ for
antineutrino scattering), where a strange quark is taken from the target and a
charm quark fragments into a charmed hadron. In the case of neutrino
scattering, the $d$-quark contribution becomes dominant at $x\gtrsim 0.2$.
This behavior reflects the fact that the $d$-quark distribution includes a
valence contribution that increases the PDF at large values of $x$ and
compensates for the smaller CKM matrix element. In comparison, the
$\bar{d}$-quark initiated processes in the antineutrino scattering have only a
small contribution to the cross section as the $\bar{d}$-quark PDFs are of the
same order as for $\bar{s}$ quarks. Interestingly, we notice that at NLO the
gluon-initiated processes have a negative contribution to the cross section,
which is related to the definitions of the NLO PDFs and FFs in the
$\overline{\mathrm{MS}}$ scheme. Overall, these results confirm that the
considered cross sections are very sensitive to the $s$\- and $\bar{s}$-quark
PDFs.

The contributions of different charmed hadrons to the full dimuon cross
section are shown in figure 7. We find that the relative contributions from
different hadrons are independent of the kinematics as the relative
contributions reflect the ratio

$\frac{f(c\to h)/m_{h}\Gamma_{h}^{\text{tot}}}{\sum_{h}f(c\to
h)/m_{h}\Gamma_{h}^{\text{tot}}}\,,$ (29)

which can be understood from eq. (22). The fragmentation fractions $f(c\to h)$
for each hadron $h$ can be found e.g. in refs. ZEUS:2005pvv ; Lisovyi:2015uqa
and the total decay widths $\Gamma_{h}^{\text{tot}}$ in ref. Workman:2022ynf .
The contribution of $\Lambda_{c}^{+}$ is multiplied by $1.14$ to account for
other baryonic states ZEUS:2005pvv . The biggest contributions of roughly
$40\text{\,}\mathrm{\char 37\relax}$ are from $D^{0}$ and $D^{+}$ production,
with $D_{s}$ contributing around $10\text{\,}\mathrm{\char 37\relax}$. The
smallest contribution comes from $\Lambda_{c}^{+}$ (and other baryon states)
at about $3\text{\,}\mathrm{\char 37\relax}$. The calculated fractions are
indeed well in line with the expectations of eq. (29).

###  4.3 Scale dependence

To gauge the effects of missing higher orders, we employ the standard approach
of computing observables with different combinations of the renormalization
scale $\mu_{\text{ren}}$, factorization scale $\mu_{\text{fact}}$, and
fragmentation scale $\mu_{\text{frag}}$. These scales are varied for the
combinations
$\mu_{\text{ren}},\mu_{\text{fact}},\mu_{\text{frag}}\in\\{\frac{1}{2}Q,Q,2Q\\}$,
with the restriction
$\frac{1}{2}\mu_{\text{ren}}\leq\mu_{\text{fact}},\mu_{\text{frag}}\leq
2\mu_{\text{ren}}$ to limit the appearance of excessively large scale
logarithms. This results in 17 different scale-choice combinations from which
the resulting error band is derived by taking an envelope. For some NuTeV and
CCFR data points, the scale at which the PDFs and FFs are evaluated becomes
smaller than the parametrization scale of the respective analyses. This is
especially the case when computing the scale dependence with
$\mu_{\text{fact}},\mu_{\text{frag}}=\frac{1}{2}Q$. In these cases, we freeze
the factorization scales to the minimum values of $\mu_{\rm fact}$ specified
by the PDF analysis and do not extrapolate beyond the edge of the grid. As we
use the values for $\alpha_{s}$ from the PDF sets, we freeze also the
renormalization scale $\mu_{\text{ren}}$ to the same value as
$\mu_{\text{fact}}$. The fragmentation functions become scale-independent
below the scale $m_{c}=$1.5\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$.

Figure 8: Scale uncertainties in the neutrino (top) and antineutrino (bottom)
dimuon cross sections, computed with EPPS21 Eskola:2021nhw . Absolute cross
sections are shown on the left and cross sections normalized to the central LO
values on the right. The perturbative order in the figure refers only to the
partonic cross section. In each case, the PDF is evaluated at NLO. The NNLO
(NLP) refers to the approximate NNLO calculation of ref. Abele:2021nyo .
Theoretical predictions are compared with the NuTeV data NuTeV:2007uwm . The
uncertainty band depicts the envelope of all scale choice combinations. The
cross section values should be multiplied by $G_{F}^{2}ME_{\nu}/100\pi$.

The resulting uncertainties from the scale variations for the dimuon
production cross section in neutrino and antineutrino scattering are shown in
figure 8 for LO and NLO, including also the central, approximate NNLO
correction. We have not computed an uncertainty band associated with the
approximate NNLO coefficient functions. For this comparison, figure 8 shows
only representative kinematic bins of $y$ and $E$. As is to be expected,
moving from the LO to the NLO calculation decreases the uncertainties from the
scale variations significantly. The scale uncertainty also grows larger at
smaller $x$, which corresponds to smaller $Q^{2}$. Towards this region, any
calculation based on perturbative QCD becomes less reliable due to the
stronger scale dependence of $\alpha_{s}$ and PDFs. In addition, effects from
the missing mass terms at the matrix-element level could become more relevant.
It is also possible that higher-twist effects start to play a role when
$Q^{2}$ is small enough.

The NLO corrections are rather significant throughout the available kinematic
range but stay below $\sim 30\,\%$ apart from the lowest values of $Q^{2}$. At
small values of $x$, the approximate NNLO corrections are typically small.
Towards larger values of $x$ their relevance appears to increase and in the
case of neutrino scattering the approximate NNLO result tends to lie near the
upper edge of the NLO scale-uncertainty band. In the antineutrino-scattering
case, the relative NNLO correction is larger than in the neutrino-scattering
case due to the different $x$-dependence of the $\bar{d}$-quark distribution
compared to the $d$-quark distribution. Interestingly, the approximate NNLO
correction is here always positive unlike in the case of inclusive charm
production in the fixed-flavour number scheme Gao:2017kkx . However, we stress
that here we have applied only the approximative NNLO coefficient functions
and the situation may still change once the full NNLO coefficient functions
are in place. On the other hand, there is no a priori reason why the NNLO
corrections should be exactly the same in these two different processes,
namely inclusive and semi-inclusive charm production.

###  4.4 Nuclear-PDF uncertainties

We will now compare our NLO calculations with different nuclear PDFs with the
experimental data from NuTeV and CCFR. Figures 9 and 10 compare the NuTeV
NuTeV:2007uwm  cross sections on dimuon production in neutrino-nucleus and
antineutrino-nucleus scatterings with NLO calculations using EPPS21,
nCTEQ15HQ, and nNNPDF3.0 nuclear PDFs. While the calculations with the three
different nuclear PDF sets are in agreement with each other when their
uncertainties are accounted for, some systematic differences can be observed.
At the smallest values of $x$, the EPPS21 calculation gives a lower cross
section than the other sets, whereas at $0.1<x<0.4$, nCTEQ15HQ is below the
others. For $x>0.2$ the nNNPDF3.0 set, in turn, yields the highest cross
section. These differences directly reflect the differences in the general
shapes of the strange-quark distributions, as is seen in figure 11, further
underlying the sensitivity of this observable to the $s$-quark PDF. The
calculations with all three nuclear PDFs describe the NuTeV data well within
the uncertainties, though at $x<0.1$ the data seem to slightly prefer the
lower cross section from EPPS21 compared to the others. We note here that
these data were included in the fit of nNNPDF3.0 nuclear PDFs, where the
factorized acceptance correction was used. In the case of antineutrino
scattering, the cross sections are slightly smaller overall and the fall-off
as $x$ increases is steeper. This is related to having a significant large-$x$
contribution from the valence $d$-quark PDFs in the case of neutrino
scattering, as was shown in figure 6. In general, the systematics of the
relative behaviour of calculations with different nuclear PDFs are very
similar between neutrino and antineutrino scattering and within uncertainties
all sets describe the data.

To study the differences between the three PDF sets in more detail, we show a
representative example of the cross sections normalized to the central EPPS21
value in figure 12. It particularly highlights the larger cross sections
obtained with nNNPDF3.0 at $0.3<x<0.4$, in contrast to EPPS21 and nCTEQ15HQ.
This behavior can again be understood from the strange-quark distributions in
figure 11. While these NuTeV data do not reach $x$ values much above $x=0.3$,
the inclusive neutrino-nucleus data still extend to higher values of $x$ and
particularly in the case of antineutrino scattering, the higher large-$x$
strange quarks of nNNPDF3.0 are well visible also there, see e.g. figure 5 of
ref. Klasen:2023uqj .

Figures 13 and 14 show the analogous results for the older CCFR data for
neutrino and antineutrino scattering, respectively. While the $y$ and
$E_{\nu}$ bins are slightly different, the kinematic reach of these data is
similar to the NuTeV data. In general, the agreement with the data and the
calculations applying different nuclear PDF sets are similar as in the case of
the NuTeV data, though the preferred PDF set at small-$x$ varies for different
bins of $y$: the smallest $y$ data tend to be best described with EPPS21, but
in the largest $y$ bins the larger cross sections provided by nCTEQ15HQ and
nNNPDF3.0 are preferred.

Figure 9: Neutrino dimuon cross sections in the NuTeV kinematics evaluated
using the EPPS21 Eskola:2021nhw , nCTEQ15HQ Duwentaster:2022kpv , and
nNNPDF3.0 AbdulKhalek:2022fyi  PDF sets. The calculations are done at next-to-
leading order. The uncertainty bands depict the PDF uncertainties with a
$90\text{\,}\mathrm{\char 37\relax}$ confidence interval. The theoretical
calculations are compared against NuTeV data from ref. NuTeV:2007uwm . The
cross section values should be multiplied by $G_{F}^{2}ME_{\nu}/100\pi$.
Figure 10: Same as figure 9, but for antineutrino scattering.

Figure 11: Left: The strange-quark distribution
$xs(x,Q^{2}=$10\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}^{2}$)$ for the EPPS21
Eskola:2021nhw , nCTEQ15HQ Duwentaster:2022kpv , and nNNPDF3.0
AbdulKhalek:2022fyi  nuclear PDFs. The band indicates the
$90\text{\,}\mathrm{\char 37\relax}$ confidence-level uncertainties of the
nuclear PDFs. Right: The same distribution as on the left, but normalized to
the central EPPS21 value.

Figure 12: Same as in figure 9 (left) and figure 10 (right), but normalized to
the central EPPS21 Eskola:2021nhw  value. Figure 13: Same as figure 9, but
with CCFR kinematics and data from ref. NuTeV:2001dfo . Figure 14: Same as
figure 13, but with antineutrino scattering.

###  4.5 Nuclear effects

As the NuTeV and CCFR data have been taken with an iron target, nuclear
effects are involved Klasen:2023uqj . In fits of free-proton PDFs, one
approach seen in the literature NuTeV:2005wsg ; Hou:2019efy  (see also ref.
Accardi:2021ysh ) has been to model the nuclear correction as an overall
multiplicative factor depending only on $x$, $f(x)$, instead of applying the
correction on a flavor-by-flavor basis. For example, ref. NuTeV:2005wsg  uses
the so-called SLAC/NMC form,

$f(x)=1.10-0.36x-0.28e^{-21.94x}+2.77x^{14.42},$ (30)

which is found by fitting to charged-lepton neutral-current DIS processes
Seligman:1997fe . It should be noted that this correction was originally
fitted to the ratio $F_{2}^{A}/F_{2}^{D}$ of the structure function $F_{2}$
from nuclear targets (iron and calcium) to deuterium.

We can calculate the isospin-corrected multiplicative nuclear effect $R$ by
taking the ratio between the dimuon cross section computed with nuclear PDFs
and the dimuon cross section computed with isospin-corrected proton PDFs

$R\equiv\frac{\mathrm{d}\sigma(\nu_{\mu}N\to\mu\mu X)|_{\text{nuclear
PDF}}}{\mathrm{d}\sigma(\nu_{\mu}N\to\mu\mu X)|_{\text{isospin PDF}}}.$ (31)

The isospin-corrected proton PDFs $f_{A}^{p}$ simply take the number of
protons ($Z$) and neutrons ($N$) in a nucleus into account via isospin
symmetry and are defined as

$\displaystyle u_{A}^{p}(x,Q^{2})$
$\displaystyle=\frac{Z}{A}u^{p}(x,Q^{2})+\frac{N}{A}d^{p}(x,Q^{2}),$ (32)
$\displaystyle d_{A}^{p}(x,Q^{2})$
$\displaystyle=\frac{Z}{A}d^{p}(x,Q^{2})+\frac{N}{A}u^{p}(x,Q^{2}),$ (33)

where $f^{p}$ are the usual proton PDFs. The definition is analogous for
$\overline{u}_{A}^{p}$ and $\overline{d}_{A}^{p}$. For other flavors and the
gluon, $f_{A}^{p}(x,Q^{2})=f^{p}(x,Q^{2})$.

Figure 15 shows the isospin-corrected nuclear ratio for neutrino and
antineutrino scattering. For each nuclear PDF set, we use the corresponding
baseline proton PDF set. For EPPS21 Eskola:2021nhw , this is CT18ANLO
Hou:2019efy , while nCTEQ15HQ Duwentaster:2022kpv  and nNNPDF3.0
AbdulKhalek:2022fyi  come with their own baseline proton PDFs. While the three
considered nuclear PDF sets mostly agree within their uncertainties, there are
some clear differences as well. The EPPS21 result roughly follows eq. (30),
where the correction is above unity in the region $0.05<x<0.26$
(antishadowing) and below unity elsewhere (shadowing, EMC effect). On the
other hand, the central curve of nCTEQ15HQ is above unity when $x<0.13$, but
below it when $x>0.13$. Finally, nNNPDF3.0 prefers to stay below unity in the
entire $x$-range considered here. There are also some differences between
neutrino and antineutrino scattering particularly in the case of nCTEQ15HQ: at
$x=0.4$, the nuclear effect in the neutrino case is nearly double of that in
the case of antineutrinos. The mutual behaviour of nuclear effects with
different nuclear PDFs follows rather well the relative ordering of nuclear
effects in the strange-quark distributions, see e.g. figure 4 of ref.
Klasen:2023uqj . Overall, we conclude that the nuclear effects in the dimuon
processes come with a significant uncertainty.

Figure 15: Isospin-corrected nuclear ratio, defined in eq. (31), and evaluated
using EPPS21 Eskola:2021nhw , nCTEQ15HQ Duwentaster:2022kpv , and nNNPDF3.0
AbdulKhalek:2022fyi  at NLO. The figure on the left shows the neutrino
scattering case and the figure on the right shows the antineutrino case. The
bands depict the PDF uncertainty at $90\text{\,}\mathrm{\char 37\relax}$
confidence level. As a reference, the SLAC/NMC form Seligman:1997fe  of eq.
(30) is shown as the thicker red curve.

###  4.6 Effective acceptance

Figure 16: Acceptance correction, as defined in eq. (34), in neutrino
scattering and in the NuTeV kinematics evaluated using the EPPS21
Eskola:2021nhw , nCTEQ15HQ Duwentaster:2022kpv , and nNNPDF3.0
AbdulKhalek:2022fyi  PDF sets. The calculation is done at NLO, except for the
case explictly marked as LO. The uncertainty bands correspond to the PDF
uncertainties at $90\text{\,}\mathrm{\char 37\relax}$ confidence level. Our
calculation is compared against the DISCO Monte-Carlo calculation in ref.
Mason:2006qa . Figure 17: Same as in figure 16, but for antineutrino
scattering.

The framework we have adopted allows us to directly implement kinematical cuts
for the final-state muons, but in order to compare our results to the
factorized approach in eq. (1), we calculate the equivalent of the effective
acceptance correction $\mathcal{A}$ from

$\mathcal{A}=\frac{1}{\mathcal{B}_{\mu}}\frac{\mathrm{d}\sigma(\nu_{\mu}N\to\mu\mu
X)}{\mathrm{d}\sigma(\nu_{\mu}N\to\mu cX)},$ (34)

where $\mathcal{B}_{\mu}$ is the average semileptonic branching ratio of charm
mesons. The denominator $\mathrm{d}\sigma(\nu_{\mu}N\to\mu cX)$ corresponds to
the inclusive charm production cross section and, as the numerator, it is
calculated in the variable flavour-number scheme approximating the quark-mass
effects with the slow-rescaling parameter in eq. (8). In addition to channels
with a charm quark explicitly in the final state, it includes contributions in
which the presence of the final-state charm in implicit, e.g. the partonic
channel $\bar{c}\rightarrow\bar{s}$ at LO.

It is not obvious what value one should use for $\mathcal{B}_{\mu}$. Some fits
Hou:2019efy ; NNPDF:2017mvq ; Martin:2009iq  have used a value of $0.099\pm
0.012$, which has been used by the CCFR collaboration CCFR:1994ikl . The PDG
value $0.087\pm 0.005$ is an average of several measurements Workman:2022ynf ,
with more values listed in ref. Mason:2006qa . It should be noted that the PDG
value of $\left|V_{cd}\right|$ is also dependent on the value of
$\mathcal{B}_{\mu}$. Several groups have moved away from the CCFR value for
$\mathcal{B}_{\mu}$, see for example refs. Ball:2018twp ; Harland-Lang:2014zoa
for more discussion. For figures 16 and 17, we follow ref. Harland-
Lang:2014zoa  and use $0.0919\pm$10\text{\,}\mathrm{\char 37\relax}$$, which
is based on FNAL E531 data Bolton:1997pq . While the value of
$\mathcal{B}_{\mu}$ can have a significant effect on the dimuon production
cross section values when using eq. (1), its impact in global fits can be
mitigated by including the uncertainty of $\mathcal{B}_{\mu}$ as a correlated
systematic uncertainty in the global analysis.

In figures 16 and 17, we compare the acceptance values, computed according to
eq. (34), to those from a DISCO NLO Monte-Carlo calculation Mason:2006qa  for
neutrino and antineutrino scattering, respectively. It should be noted that
the $10\text{\,}\mathrm{\char 37\relax}$ normalization uncertainty of
$\mathcal{B}_{\mu}$ in eq. (34) is not included in figures 16 and 17. While
the NLO values of $\mathcal{A}$ from the DISCO Monte-Carlo calculation of ref.
Mason:2006qa  are roughly in line with the values computed here directly,
there are clear systematic differences which are roughly at the $10\,\%$
level. Whereas our calculation shows a general trend of decreasing acceptance
with increasing $x$, the DISCO Monte-Carlo calculation tends to have a modest
rise with increasing $x$. Furthermore, we notice that for small values of $y$
and $E$, our results tend to be above the DISCO calculation, but with
increasing $y$ and $E$ the ordering goes the other way around. Thus it is
clear that a single value of $\mathcal{B}_{\mu}$ is not enough to match our
calculation to the DISCO one. The acceptances calculated with different
nuclear PDFs are in a rather good mutual agreement with all sharing the same
systematic behaviour, but especially in the antineutrino case shown in figure
17 there are larger uncertainties originating from the nuclear PDFs. In
addition to the NLO results, figures 16 and 17 also show the LO acceptance
correction, where both the dimuon production and inclusive charm production
cross sections are calculated at LO. As can be seen, the LO and NLO
acceptances differ quite significantly. Based on these comparisons, we see
that the acceptances in eq. (1) depend both on the nuclear PDFs as well as on
the perturbative order of the calculation.

##  5 Conclusion and outlook

We have presented an NLO calculation of dimuon production in semi-inclusive
deep-inelastic neutrino-nucleus scattering. In comparison to how this process
is usually calculated in PDF fits, our calculation involves scale-dependent
FFs to turn final-state charm quarks into charmed hadrons while at the same
time resumming collinear final-state radiation, and a decay function fitted to
$e^{+}e^{-}$ data to produce muons from the charmed-hadron decays. The major
advantage of our approach is that we do not need to rely on the approximative,
factorized acceptance and branching-ratio corrections to account for the
experimental cuts on the final-state muon energies. Instead, the kinematical
cuts are a natural ingredient of the calculation. The perturbative accuracy of
the presented framework can be systematically improved by including higher-
order corrections to the coefficient functions and employing PDFs and FFs
evolving with higher-order splitting functions. The approach is also
generalizable to other processes where an intermediate-state particle decays
into a given particle.

We assessed different theoretical uncertainties arising from the decay
function, parton distribution functions, and variations in the
renormalization, factorization, and fragmentation scales. While the
uncertainties in the decay-function fit translated into few-percent
differences in the dimuon cross sections, the NLO scale uncertainty was of the
order of $10\,\%$ in the perturbative regime. We also studied the size of the
the approximative NNLO corrections in the coefficient functions and found them
to be comparable to the size of the NLO scale uncertainty. By far, the largest
uncertainty was associated with the strange-quark component of nuclear PDFs,
which could be $\sim 100\,\%$ in the data-constrained region. This indicates
that there is a potential for broader use of these data in global fits of
nuclear PDFs.

We also studied the effective acceptance and nuclear effects which are often
taken as inputs in free-proton fits to these dimuon data. As for the effective
acceptance we found systematic differences of the order of $10\,\%$ between
our calculation and widely applied acceptance corrections derived from the
DISCO Monte-Carlo framework. Moreover, the effective acceptance is clearly
dependent on the perturbative order and — particularly in the case of
antineutrino scattering — also PDF dependent. This indicates that using an
external acceptance correction can, perhaps, lead to some biases in PDF fits.
As for the nuclear correction, the latest globally analyzed nuclear PDFs led
to quite significant mutual differences with qualitatively distinct behaviour
and significant uncertainties.

In the near future, we plan to improve the framework introduced here in
several ways. First, it can be straightforwardly extended to full GM-VFNS,
thereby accounting also for the subleading effects of the charm-quark mass.
Second, the calculations for the full massless NNLO SIDIS coefficient
functions have been very recently completed and thus they will replace the
approximative ones which were still used in the current setup. Third,
similarly to the case of inclusive neutrino DIS, electroweak radiative and
target-mass corrections can be included in the SIDIS cross sections as well.
Eventually, it will be interesting to see whether all this will bring insights
to the observed tensions between dimuon production in neutrino DIS and heavy
gauge boson production at the LHC.

## Acknowledgements

We acknowledge the financial support from the Magnus Ehrnrooth foundation
(S.Y.), the Research Council of Finland Project No. 331545 (I.H.), and the
Center of Excellence in Quark Matter of the Research Council of Finland,
project 346326. The reported work is associated with the European Research
Council project ERC-2018-ADG-835105 YoctoLHC. We acknowledge grants of
computer capacity from the Finnish Grid and Cloud Infrastructure (persistent
identifier urn:nbn:fi:research-infras-2016072533) and the Finnish IT Center
for Science (CSC), under the project jyy2580.

## References

  * (1) 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 [hep-ph/0409313].
  * (2) K. Kovařík, P.M. Nadolsky and D.E. Soper, _Hadronic structure in high-energy collisions_ , _Rev. Mod. Phys._ 92 (2020) 045003 [1905.06957].
  * (3) J.J. Ethier and E.R. Nocera, _Parton Distributions in Nucleons and Nuclei_ , _Ann. Rev. Nucl. Part. Sci._ 70 (2020) 43 [2001.07722].
  * (4) M. Klasen and H. Paukkunen, _Nuclear PDFs After the First Decade of LHC Data_ , 2311.00450.
  * (5) A. Metz and A. Vossen, _Parton Fragmentation Functions_ , _Prog. Part. Nucl. Phys._ 91 (2016) 136 [1607.02521].
  * (6) V.N. Gribov and L.N. Lipatov, _Deep inelastic e p scattering in perturbation theory_ , _Sov. J. Nucl. Phys._ 15 (1972) 438.
  * (7) V.N. Gribov and L.N. Lipatov, _e+ e- pair annihilation and deep inelastic e p scattering in perturbation theory_ , _Sov. J. Nucl. Phys._ 15 (1972) 675.
  * (8) Y.L. Dokshitzer, _Calculation of the Structure Functions for Deep Inelastic Scattering and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics._ , _Sov. Phys. JETP_ 46 (1977) 641.
  * (9) G. Altarelli and G. Parisi, _Asymptotic Freedom in Parton Language_ , _Nucl. Phys. B_ 126 (1977) 298.
  * (10) T.-J. Hou et al., _New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC_ , _Phys. Rev. D_ 103 (2021) 014013 [1912.10053].
  * (11) 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) 341 [2012.04684].
  * (12) NNPDF collaboration, _The path to proton structure at 1% accuracy_ , _Eur. Phys. J. C_ 82 (2022) 428 [2109.02653].
  * (13) K.J. Eskola, P. Paakkinen, H. Paukkunen and C.A. Salgado, _EPPS21: a global QCD analysis of nuclear PDFs_ , _Eur. Phys. J. C_ 82 (2022) 413 [2112.12462].
  * (14) P. Duwentäster, T. Ježo, M. Klasen, K. Kovařík, A. Kusina, K.F. Muzakka et al., _Impact of heavy quark and quarkonium data on nuclear gluon PDFs_ , _Phys. Rev. D_ 105 (2022) 114043 [2204.09982].
  * (15) R. Abdul Khalek, R. Gauld, T. Giani, E.R. Nocera, T.R. Rabemananjara and J. Rojo, _nNNPDF3.0: evidence for a modified partonic structure in heavy nuclei_ , _Eur. Phys. J. C_ 82 (2022) 507 [2201.12363].
  * (16) ATLAS collaboration, _Measurement of the $W$-boson mass in pp collisions at $\sqrt{s}=7$ TeV with the ATLAS detector_, _Eur. Phys. J. C_ 78 (2018) 110 [1701.07240].
  * (17) E. Bagnaschi and A. Vicini, _Parton Density Uncertainties and the Determination of Electroweak Parameters at Hadron Colliders_ , _Phys. Rev. Lett._ 126 (2021) 041801 [1910.04726].
  * (18) M. Hussein, J. Isaacson and J. Huston, _A Study of the Role of the PDF Uncertainty on the LHC $W$-Boson Mass Measurement_, _J. Phys. G_ 46 (2019) 095002 [1905.00110].
  * (19) ATLAS collaboration, _Improved W boson Mass Measurement using 7 TeV Proton-Proton Collisions with the ATLAS Detector_ , .
  * (20) CMS collaboration, _Measurement of the weak mixing angle using the forward-backward asymmetry of Drell-Yan events in pp collisions at 8 TeV_ , _Eur. Phys. J. C_ 78 (2018) 701 [1806.00863].
  * (21) CCFR collaboration, _Determination of the strange quark content of the nucleon from a next-to-leading order QCD analysis of neutrino charm production_ , _Z. Phys. C_ 65 (1995) 189 [hep-ex/9406007].
  * (22) NuTeV collaboration, _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].
  * (23) NuTeV collaboration, _Measurement of the Nucleon Strange-Antistrange Asymmetry at Next-to-Leading Order in QCD from NuTeV Dimuon Data_ , _Phys. Rev. Lett._ 99 (2007) 192001.
  * (24) CHORUS collaboration, _Leading order analysis of neutrino induced dimuon events in the CHORUS experiment_ , _Nucl. Phys. B_ 798 (2008) 1 [0804.1869].
  * (25) NOMAD collaboration, _A Precision Measurement of Charm Dimuon Production in Neutrino Interactions from the NOMAD Experiment_ , _Nucl. Phys. B_ 876 (2013) 339 [1308.4750].
  * (26) I. Schienbein, J.Y. Yu, K. Kovarik, C. Keppel, J.G. Morfin, F. Olness et al., _PDF Nuclear Corrections for Charged and Neutral Current Processes_ , _Phys. Rev. D_ 80 (2009) 094004 [0907.2357].
  * (27) H. Paukkunen and C.A. Salgado, _Compatibility of neutrino DIS data and global analyses of parton distribution functions_ , _JHEP_ 07 (2010) 032 [1004.3140].
  * (28) K. Kovarik, I. Schienbein, F.I. Olness, J.Y. Yu, C. Keppel, J.G. Morfin et al., _Nuclear Corrections in Neutrino-Nucleus DIS and Their Compatibility with Global NPDF Analyses_ , _Phys. Rev. Lett._ 106 (2011) 122301 [1012.0286].
  * (29) H. Paukkunen and C.A. Salgado, _Agreement of Neutrino Deep Inelastic Scattering Data with Global Fits of Parton Distributions_ , _Phys. Rev. Lett._ 110 (2013) 212301 [1302.2001].
  * (30) K.F. Muzakka et al., _Compatibility of neutrino DIS data and its impact on nuclear parton distribution functions_ , _Phys. Rev. D_ 106 (2022) 074004 [2204.13157].
  * (31) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * (32) A. Accardi, T.J. Hobbs, X. Jing and P.M. Nadolsky, _Deuterium scattering experiments in CTEQ global QCD analyses: a comparative investigation_ , _Eur. Phys. J. C_ 81 (2021) 603 [2102.01107].
  * (33) NNPDF collaboration, _Nuclear Uncertainties in the Determination of Proton PDFs_ , _Eur. Phys. J. C_ 79 (2019) 282 [1812.09074].
  * (34) J.M. Cruz-Martinez, M. Fieg, T. Giani, P. Krack, T. Mäkelä, T. Rabemananjara et al., _The LHC as a Neutrino-Ion Collider_ , 2309.09581.
  * (35) FASER collaboration, _First Direct Observation of Collider Neutrinos with FASER at the LHC_ , _Phys. Rev. Lett._ 131 (2023) 031801 [2303.14185].
  * (36) SND@LHC collaboration, _Observation of Collider Muon Neutrinos with the SND@LHC Experiment_ , _Phys. Rev. Lett._ 131 (2023) 031802 [2305.09383].
  * (37) A. Kusina, T. Stavreva, S. Berge, F.I. Olness, I. Schienbein, K. Kovarik et al., _Strange Quark PDFs and Implications for Drell-Yan Boson Production at the LHC_ , _Phys. Rev. D_ 85 (2012) 094028 [1203.1290].
  * (38) K.J. Eskola, P. Paakkinen, H. Paukkunen and C.A. Salgado, _Proton-PDF uncertainties in extracting nuclear PDFs from $W^{\pm}$ production in p+Pb collisions_, _Eur. Phys. J. C_ 82 (2022) 271 [2202.01074].
  * (39) A. Kusina et al., _Impact of LHC vector boson production in heavy ion collisions on strange PDFs_ , _Eur. Phys. J. C_ 80 (2020) 968 [2007.09100].
  * (40) U. Baur, F. Halzen, S. Keller, M.L. Mangano and K. Riesselmann, _The Charm content of $W$ \+ 1 jet events as a probe of the strange quark distribution function_, _Phys. Lett. B_ 318 (1993) 544 [hep-ph/9308370].
  * (41) W.J. Stirling and E. Vryonidou, _Charm production in association with an electroweak gauge boson at the LHC_ , _Phys. Rev. Lett._ 109 (2012) 082002 [1203.6781].
  * (42) ATLAS collaboration, _Measurement of the production of a $W$ boson in association with a charm quark in $pp$ collisions at $\sqrt{s}=$ 7 TeV with the ATLAS detector_, _JHEP_ 05 (2014) 068 [1402.6263].
  * (43) CMS collaboration, _Measurement of Associated W + Charm Production in pp Collisions at $\sqrt{s}$ = 7 TeV_, _JHEP_ 02 (2014) 013 [1310.1138].
  * (44) CMS collaboration, _Measurement of associated production of a W boson and a charm quark in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _Eur. Phys. J. C_ 79 (2019) 269 [1811.10021].
  * (45) ATLAS collaboration, _Measurement of the production of a $W$ boson in association with a charmed hadron in $pp$ collisions at $\sqrt{s}=13\,\mathrm{TeV}$ with the ATLAS detector_, _Phys. Rev. D_ 108 (2023) 032012 [2302.00336].
  * (46) ATLAS collaboration, _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 [1203.4051].
  * (47) ATLAS collaboration, _Precision measurement and interpretation of inclusive $W^{+}$ , $W^{-}$ and $Z/\gamma^{*}$ production cross sections with the ATLAS detector_, _Eur. Phys. J. C_ 77 (2017) 367 [1612.03016].
  * (48) A.M. Cooper-Sarkar and K. Wichmann, _QCD analysis of the ATLAS and CMS $W^{\pm}$ and $Z$ cross-section measurements and implications for the strange sea density_, _Phys. Rev. D_ 98 (2018) 014027 [1803.00968].
  * (49) NNPDF collaboration, _Parton distributions from high-precision collider data_ , _Eur. Phys. J. C_ 77 (2017) 663 [1706.00428].
  * (50) 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) 212002 [1601.05430].
  * (51) J. Gao, _Massive charged-current coefficient functions in deep-inelastic scattering at NNLO and impact on strange-quark distributions_ , _JHEP_ 02 (2018) 026 [1710.04258].
  * (52) F. Faura, S. Iranipour, E.R. Nocera, J. Rojo and M. Ubiali, _The Strangest Proton?_ , _Eur. Phys. J. C_ 80 (2020) 1168 [2009.00014].
  * (53) HERMES collaboration, _Reevaluation of the parton distribution of strange quarks in the nucleon_ , _Phys. Rev. D_ 89 (2014) 097101 [1312.7028].
  * (54) I. Borsa, R. Sassot and M. Stratmann, _Probing the Sea Quark Content of the Proton with One-Particle-Inclusive Processes_ , _Phys. Rev. D_ 96 (2017) 094020 [1708.01630].
  * (55) JAM collaboration, _Strange quark suppression from a simultaneous Monte Carlo analysis of parton distributions and fragmentation functions_ , _Phys. Rev. D_ 101 (2020) 074020 [1905.03788].
  * (56) S. Kretzer, D. Mason and F. Olness, _Differential distributions for NLO analyses of charged current neutrino production of charm quarks_ , _Phys. Rev. D_ 65 (2002) 074010 [hep-ph/0112191].
  * (57) S. Kretzer and I. Schienbein, _Charged current leptoproduction of D mesons in the variable flavor scheme_ , _Phys. Rev. D_ 56 (1997) 1804 [hep-ph/9702296].
  * (58) M. Abele, D. de Florian and W. Vogelsang, _Approximate NNLO QCD corrections to semi-inclusive DIS_ , _Phys. Rev. D_ 104 (2021) 094046 [2109.00847].
  * (59) W. Furmanski and R. Petronzio, _Lepton - Hadron Processes Beyond Leading Order in Quantum Chromodynamics_ , _Z. Phys. C_ 11 (1982) 293.
  * (60) D. de Florian, M. Stratmann and W. Vogelsang, _QCD analysis of unpolarized and polarized Lambda baryon production in leading and next-to-leading order_ , _Phys. Rev. D_ 57 (1998) 5811 [hep-ph/9711387].
  * (61) M. Abele, D. de Florian and W. Vogelsang, _Threshold resummation at NLL3 accuracy and approximate N3LO corrections to semi-inclusive DIS_ , _Phys. Rev. D_ 106 (2022) 014015 [2203.07928].
  * (62) L. Bonino, T. Gehrmann and G. Stagnitto, _Semi-inclusive deep-inelastic scattering at NNLO in QCD_ , 2401.16281.
  * (63) L. Bonino, T. Gehrmann, M. Löchner, K. Schönwald and G. Stagnitto, _Polarized semi-inclusive deep-inelastic scattering at NNLO in QCD_ , 2404.08597.
  * (64) S. Goyal, R.N. Lee, S.-O. Moch, V. Pathak, N. Rana and V. Ravindran, _NNLO QCD corrections to polarized semi-inclusive DIS_ , 2404.09959.
  * (65) R.M. Barnett, _Evidence for new quarks and new currents_ , _Phys. Rev. Lett._ 36 (1976) 1163.
  * (66) V. Barone, M. Genovese, N.N. Nikolaev, E. Predazzi and B.G. Zakharov, _The Strange quark distribution_ , _Z. Phys. C_ 70 (1996) 83 [hep-ph/9505343].
  * (67) P.M. Nadolsky and W.-K. Tung, _Improved Formulation of Global QCD Analysis with Zero-mass Matrix Elements_ , _Phys. Rev. D_ 79 (2009) 113014 [0903.2667].
  * (68) CLEO collaboration, _Absolute Branching Fraction Measurements for $D^{+}$ and $D^{0}$ Inclusive Semileptonic Decays_, _Phys. Rev. Lett._ 97 (2006) 251801 [hep-ex/0604044].
  * (69) 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].
  * (70) NNPDF collaboration, _A Determination of parton distributions with faithful uncertainty estimation_ , _Nucl. Phys. B_ 809 (2009) 1 [0808.1231].
  * (71) T. Kneesch, B.A. Kniehl, G. Kramer and I. Schienbein, _Charmed-meson fragmentation functions with finite-mass corrections_ , _Nucl. Phys. B_ 799 (2008) 34 [0712.0481].
  * (72) B.A. Kniehl and G. Kramer, _Charmed-hadron fragmentation functions from CERN LEP1 revisited_ , _Phys. Rev. D_ 74 (2006) 037502 [hep-ph/0607306].
  * (73) OPAL collaboration, _A Study of charm hadron production in Z0 — $>$ c anti-c and Z0 —$>$ b anti-b decays at LEP_, _Z. Phys. C_ 72 (1996) 1.
  * (74) OPAL collaboration, _Measurement of f(c — $>$ D*+ X), f(b —$>$ D*+ X) and Gamma (c anti-c) / Gamma (hadronic) using D*+- mesons_, _Eur. Phys. J. C_ 1 (1998) 439 [hep-ex/9708021].
  * (75) Belle collaboration, _Charm hadrons from fragmentation and B decays in e+ e- annihilation at s**(1/2) = 10.6-GeV_ , _Phys. Rev. D_ 73 (2006) 032002 [hep-ex/0506068].
  * (76) ZEUS collaboration, _Measurement of charm fragmentation ratios and fractions in photoproduction at HERA_ , _Eur. Phys. J. C_ 44 (2005) 351 [hep-ex/0508019].
  * (77) M. Lisovyi, A. Verbytskyi and O. Zenaiev, _Combined analysis of charm-quark fragmentation-fraction measurements_ , _Eur. Phys. J. C_ 76 (2016) 397 [1509.01061].
  * (78) NuTeV collaboration, _Precise measurement of neutrino and anti-neutrino differential cross sections_ , _Phys. Rev. D_ 74 (2006) 012008 [hep-ex/0509010].
  * (79) W.G. Seligman, _A Next-to-Leading Order QCD Analysis of Neutrino - Iron Structure Functions at the Tevatron_ , Ph.D. thesis, Nevis Labs, Columbia U., 1997.  10.2172/1421736.
  * (80) D.A. Mason, _Measurement of the strange - antistrange asymmetry at NLO in QCD from NuTeV dimuon data_ , Ph.D. thesis, Oregon U., 2006.  10.2172/879078.
  * (81) A.D. Martin, W.J. Stirling, R.S. Thorne and G. Watt, _Parton distributions for the LHC_ , _Eur. Phys. J. C_ 63 (2009) 189 [0901.0002].
  * (82) 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) 204 [1412.3989].
  * (83) T. Bolton, _Determining the CKM parameter V(cd) from neutrino N charm production_ , hep-ex/9708014.


Comparison of systemic PMs determined with (red) and without (blue) the likelihood term for line-of-sight velocities. The x-axis and y-axis show the $\mu_{\alpha,*}$ and $\mu_{\delta}$ component, respectively. The sources of the spectroscopic works are those used in <cit.>, Kirby et al., 2013,  Kirby et al., 2015,  Carlin & Sand, 2018, Torrealba et al., 2019,  Simon et al., 2020, Longeard et al., 2021.

As Figure <ref>.

As Figure <ref>.

See Figure <ref>

See Figure <ref>

§ PLOTS ON ORBITAL HISTORIES

Here we include plots showing the orbital evolution in the past 3 Gyr in the triaxial light MW potential with and without the inclusion of the LMC.

Distance from the MW center as a function of time for the past 6 Gyr (at present time t$=0$) for the  "perturbed" and "isolated" "Light MW" potentials (cyan and orange lines, respectively), i.e. in the case with and without the the infall of a massive LMC; the orbits are determined from the observed, error-free, bulk motions.

As <ref>. For Pisces II the orbits shown are obtained using the systemic motion from the inclusion of the spectroscopic information.

As Fig. <ref>. For Tucana V the orbits shown are obtained using the systemic motion from the inclusion of the spectroscopic information.