CERN-TH-2024-091, MPP-2024-121

e1 e2 e3 e4 e5

11institutetext: Theoretical Physics Department, CERN, 1211 Geneva 23,
Switzerland 22institutetext: Max-Planck-Institut für Physik, Boltzmannstraße
8, 85748 Garching, Germany 33institutetext: Institute for Theoretical Physics,
University of Tübingen, Auf der Morgenstelle 14, 72076 Tübingen, Germany
44institutetext: Physik-Department, Technische Universität München, James-
Franck-Strasse 1, 85748 Garching, Germany

# An event generator for neutrino-induced Deep Inelastic Scattering and
applications to neutrino astronomy

Silvia Ferrario Ravasioaddr1,e1  Rhorry Gauldaddr2,e2  Barbara Jägeraddr3,e3
Alexander Karlbergaddr1,e4  Giulia Zanderighiaddr2,addr4,e5
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(August 28, 2024)

###### Abstract

We extend the recently presented, fully exclusive, next-to-leading-order
accurate event generator for the simulation of massless neutral- and charged-
current deep inelastic scattering (DIS) to the case of incoming neutrinos. The
generator can be used to study neutrino-nucleon interactions at (ultra) high
energies, and is relevant for a range of fixed-target collider experiments and
large-volume neutrino detectors, investigating atmospheric and astrophysical
neutrinos. The matching with multi-purpose event generators such as PYTHIA 8
is performed with the POWHEG method, and accounts for parton showering and
non-perturbative effects such as hadronization. This makes it possible to
investigate higher-order perturbative corrections to realistic observables,
such as the distribution of charged particles. To illustrate the capabilities
of the code we provide predictions for several differential distributions in
fixed-target collisions for neutrino energies up to $1\leavevmode\nobreak\
\mathrm{PeV}$.

††journal: Eur. Phys. J. C

###### Contents

  1. 1 Introduction
  2. 2 Details of the implementation
    1. 2.1 Fixed-target experiments
    2. 2.2 Nucleon targets
    3. 2.3 Variable neutrino flux
    4. 2.4 On the momentum mappings, mass effects and possible extensions to more complex processes
  3. 3 Fixed-order validation
  4. 4 Phenomenological results
    1. 4.1 Particle multiplicities
    2. 4.2 Energy-based distributions
    3. 4.3 Charm production
  5. 5 Conclusions
  6. A Appendix DIS process selection

##  1 Introduction

Neutrinos, together with photons, are the most abundant elementary particles
in the universe. While the properties of photons are extremely well
understood, there are still many outstanding questions regarding neutrinos.
For instance, the origin and nature of neutrino masses (Dirac vs. Majorana
mass) is not understood, nor is their mass hierarchy (normal vs. inverted
ordering). Furthermore, neutrinos provide a portal to beyond the Standard
Model (BSM) physics, making neutrino experiments at the luminosity frontier
sensitive to such BSM interactions (see e.g. Ref. Batell:2009di  for a
review).

Neutrino properties are difficult to measure because they only interact though
the weak force. For this reason, their study often requires large-volume
detectors, which have enabled the discovery of (ultra) high-energy cosmic
neutrinos in 2014, the observation of an astrophysical source of energetic
neutrinos accompanied by gamma-ray emissions in 2018, and the determination of
the oscillation properties of multi-$\mathrm{GeV}$ energy atmospheric
neutrinos (see e.g. Ref. Ackermann:2022rqc  for a review of these and several
recent results). Ongoing experiments such as ANTARES ANTARES:2011hfw , Baikal
BAIKAL:2005qnn , IceCube IceCube:2016zyt , and KM3NeT KM3Net:2016zxf , will
continue to extract information on (ultra) high-energy neutrinos to which
their detectors are exposed. Moreover, a range of proposed next-generation
detectors will facilitate precise measurements of (ultra) high-energy
neutrinos from atmospheric and cosmic sources. This advancement will usher in
a new era of precision, enabling to probe neutrino properties, their
interactions, and fundamental symmetries at the highest possible energies.
Furthermore, this programme will be instrumental to discover and characterize
the astrophysical sources of the most energetic cosmic and gamma-rays.

Additional data opportunities come from high-luminosity experiments. For
example, measurements of neutrino-matter scattering at collider facilities
(e.g. charm production measured by NuTeV NuTeV:2007uwm ) have provided
important information on the hadron structure. Forward-physics facilities such
as SND@LHC SNDLHC:2022ihg ; SNDLHC:2023pun , SHiP SHiP:2015vad , and
FASER$\nu$ FASER:2019dxq ; FASER:2020gpr ; FASER:2023zcr , are already taking
data and the Forward Physics Facility (FPF) is on the horizon for the HL-LHC
Anchordoqui:2021ghd ; Feng:2022inv . A major goal of each of these experiments
is to extract the flavour and energy-dependence of the neutrino flux to which
their detector is exposed. This requires, in addition to a detailed
understanding of the detector, precise knowledge of the expected differential
rates of neutrino-nucleon scattering for varying neutrino flavour and energy.

At the large energies under consideration (multi-$\mathrm{GeV}$ and above),
the scattering rate of neutrinos with matter is dominated by the deep
inelastic scattering (DIS) process. The role of theory in this context is thus
an important one: it provides a well defined (and rigorously tested)
computational framework, that of collinear factorisation Collins:1989gx , to
predict the differential scattering rates of neutrinos. This framework is
reliable provided the exchanged momentum, $Q^{\mu}$, satisfies $|Q^{2}|\gtrsim
m_{p}^{2}$, $m_{p}$ being the proton mass, and can be applied across many
orders of magnitude in neutrino energy. It relies on a combination of
perturbative QCD ingredients, and of the knowledge of the universal partonic
content of the colliding hadrons (as extracted from global analyses of hadron
collider data), see Ref. Ethier:2020way  for a recent review.

This theoretical framework can be straightforwardly applied to the case of
(ultra) high-energy neutrino-nucleon scattering by expressing the differential
cross-section in terms of DIS structure functions (see for example the
discussion in Section II of Cooper-Sarkar:2011jtt ). The structure functions
encapsulate the strong dynamics of the nucleon as struck by an exchanged gauge
boson, and they can be predicted through the convolution of parton
distribution functions (PDFs) with a set of perturbatively calculated
coefficient functions. The simplicity of this approach stems from the fact
that the structure functions provide an inclusive description of all QCD
radiation in the scattering process. On the other hand, it is limited as
predicted cross-sections are differential only in quantities inclusive over
QCD radiation, such as the leptonic momentum transfer $Q^{2}$ and the Bjorken
momentum fraction, $x_{\mathrm{B}}$. The massless hard coefficient functions
that enter into the structure functions have been computed at 3-loops
SanchezGuillen:1990iq ; vanNeerven:1991nn ; Zijlstra:1991qc ; Zijlstra:1992qd
; Zijlstra:1992kj ; vanNeerven:1999ca ; vanNeerven:2000uj ; Moch:1999eb ;
Moch:2004xu ; Vermaseren:2005qc ; Vogt:2006bt ; Moch:2007rq ; Davies:2016ruz ;
Blumlein:2022gpp . Following the structure-function approach, dedicated
theoretical studies of neutrino-nucleon DIS at high energies have appeared
over the years, both at leading-order (LO) Gandhi:1998ri ; Gluck:1998js ;
Cooper-Sarkar:2007zsa ; Connolly:2011vc , next-to-leading order (NLO) Cooper-
Sarkar:2011jtt  and recently at next-to-next-to-leading order (NNLO) in QCD
Bertone:2018dse ; Xie:2023suk . The impact of the physics effects due to
heavy-quark masses, nuclear modifications of PDFs, and resummation of
small-$x$ contributions has been studied in Refs. Bertone:2018dse ;
Xie:2023suk , the role of certain classes of QED effects has been investigated
in Refs. Seckel:1997kk ; Alikhanov:2015kla ; Gauld:2019pgt ; Zhou:2019vxt ;
Zhou:2019frk ; Xie:2023qbn , and effects beyond collinear factorisation have
also been discussed Jalilian-Marian:2003ghc ; Fiore:2005wf ; Block:2013nia ;
Albacete:2015zra ; Goncalves:2015fua ; Arguelles:2015wba .

Predictions obtained in this way provide an important benchmark for
differential DIS cross-sections in terms of QCD-inclusive quantities (e.g.
distributions of $Q^{2}$ and $x_{\mathrm{B}}$), as well as the total cross-
section. However, they do not provide an exclusive description of the
radiation which is generated in the scattering process. This is a significant
limitation for many analyses at current (and future) neutrino experiments
which aim to reconstruct the energy and direction of the incoming neutrino,
and which rely on an accurate description of the properties of final-state
radiation (such as the distribution of electromagnetically charged and neutral
particles) to do so. A step towards overcoming this issue is made in the
current work with the development of an event generator for the simulation of
neutrino-induced massless neutral- and charged-current DIS based on the POWHEG
Nason:2004rx ; Frixione:2007vw  method. The predictions obtained with this
program are accurate at NLO in QCD and can be matched with a multi-purpose
Shower Monte Carlo generator to provide a fully exclusive description of the
scattering process. The implementation is based on the existing generator for
charged-lepton induced DIS processes presented in Banfi:2023mhz , and has been
implemented in the publicly available framework POWHEG-BOX-RES Jezo:2015aia .
The code can be obtained from svn://powhegbox.mib.infn.it/trunk/User-
Processes-RES/DIS.

While this paper was being finalised, an NLO accurate event generator
implementation for lepton-hadron DIS was presented Buonocore:2024pdv . This
implementation is based on the POWHEG-BOX-V2 framework Alioli:2010xd , and has
a particular focus on processes with a heavy lepton, such as a tau neutrino,
and/or a heavy charm quark in the final state. We briefly discuss the
differences between the two codes in 2.4.

The structure of the paper is as follows: in Sec. 2 we summarise the main
details of the process implementation and new features as compared to the
existing generator which describes charged-lepton induced DIS; a validation of
the code for various DIS subprocesses is provided in Sec. 3; in Sec. 4 we
present phenomenological results for several distributions of charged
particles and charmed hadrons for incident neutrino energies of $10^{5}$ and
$10^{6}\leavevmode\nobreak\ \mathrm{GeV}$. Concluding remarks are presented in
Sec. 5. A complete list of all the new features in the code, and how to use
them, is provided in A.

##  2 Details of the implementation

In this section we discuss the extensions needed to augment the POWHEG-BOX-RES
generator for massless neutral- and charged-current DIS, presented in Ref.
Banfi:2023mhz , to allow for the inclusion of initial-state neutrinos and
generic (massive) nuclear targets. The POWHEG-BOX-RES framework combines NLO-
QCD calculations with parton showers (PS) according to the POWHEG method, and
was originally only designed to handle hadron-hadron collisions. One of the
main novelties of Ref. Banfi:2023mhz  was the design of new momentum mappings
that preserve the special kinematics of DIS in the FKS subtraction formalism
Frixione:1995ms ; Frixione:2007vw  as implemented in the POWHEG-BOX-RES
framework.

The original generator of Ref. Banfi:2023mhz  was designed to describe DIS
reactions resulting from the collision of a massless proton with a charged
lepton, relevant to interpret data from, for instance, HERA and the
forthcoming Electron Ion Collider (EIC). It was since extended to also include
polarised beams in Ref. Borsa:2024rmh .

The extension presented here contains three new major features: 1. The
incoming lepton can now be of any species, in particular it can be a neutrino
or a charged lepton; 2. The code can now handle a massive nucleon at rest, of
relevance to fixed-target experiments; 3. A variable flux can be supplied for
the incoming lepton beam. The handling of massive nucleons at rest is
described in Sec. 2.1, and a discussion of how to consistently account for the
nuclear target PDFs can be found in Sec. 2.2. Although in this paper we focus
on phenomenological studies of neutrino beams with fixed energy, we discuss
how to include a variable flux in Sec. 2.3. Finally in Sec. 2.4 we comment on
our momentum mappings and how mass effects are approximately included.

###  2.1 Fixed-target experiments

By default, the POWHEG-BOX-RES can only handle collisions of massless beams.
In this section we therefore describe how to perform fixed-target collisions,
using a set of massless beams. Denoting the energies of two massless colliding
beams in the laboratory frame by $E_{1}$ and $E_{2}$, the POWHEG-BOX builds
the four-momenta of the beam particles as follows:

$\displaystyle k_{\rm beam,1}=$
$\displaystyle\left\\{E_{1},0,0,+E_{1}\right\\},$ $\displaystyle k_{\rm
beam,2}=$ $\displaystyle\left\\{E_{2},0,0,-E_{2}\right\\}.$ (1)

These four-vectors are then used to construct the momenta of the incoming
elementary fermions entering the scattering process.

To account for the collision of a beam of massless particles of energy $E$
with a fixed target nucleon (i.e. proton or neutron) of mass $m$ we extend
this approach by effectively treating the nucleon as massless. In the fixed-
target frame the true momenta are given by the lepton beam momentum, $P_{1}$,
and the fixed target momentum, $P_{2}$,

$\displaystyle P_{1}=$ $\displaystyle\left\\{E,0,0,E\right\\},$ $\displaystyle
P_{2}=$ $\displaystyle\left\\{m,0,0,0\right\\}.$ (2)

From these momenta we obtain a centre-of-mass energy, $E_{\mathrm{CM}}$, via

$E_{\mathrm{CM}}^{2}=(P_{1}+P_{2})^{2}=2mE+m^{2}.$ (3)

We then trivially observe that if we pick $E_{1}=E_{2}=E_{\mathrm{CM}}/2$ in
Eq. (2.1) we can construct a set of massless momenta that coincide with the
centre-of-mass frame of the fixed-target collision. Now consider the boost
from the centre-of-mass frame to the _true_ fixed-target frame. Applying this
boost to our newly constructed massless momenta we can construct massless beam
momenta in Eq. (2.1) where the energies of the beams are set to

$\displaystyle E_{1}=E+m/2,\qquad E_{2}=m/2.$ (4)

Both the massless centre-of-mass and massless fixed-target momenta satisfy
$k_{\rm beam,1}+k_{\rm beam,2}=P_{1}+P_{2}$ by construction, but do not
preserve the mass of $P_{2}$. In practice we expect the massles construction
to be reliable as long as $m/E\ll 1$. The two sets of momenta result in
equivalent predictions, since they are related by a boost, but in practice we
find that using the centre-of-mass momenta is numerically more stable for
ultra-high energy collisions ($E/m\gtrsim 10^{5}-10^{6}$). We provide both
options in the code, as described in A.

We note that when interfacing the events to the parton shower, e.g. PYTHIA
8, the actual mass of the nucleon is restored while retaining the centre-of-
mass energy of the two beams, thereby restoring the correct kinematics.

###  2.2 Nucleon targets

When considering lepton scattering off the nucleons of a bound nucleus, it is
important to differentiate whether the nucleon target is a proton or a
neutron. This distinction is relevant for the eventual matching to the parton
shower, where the quantum numbers of the nucleon remnant must be known. The
selection of the nucleon type in the powheg.input file can be made by setting
the integer `ih2`, as described in A. For the selection of a neutron, we
provide the option to either directly use neutron PDFs, or to instead provide
a set of proton PDFs which the program then internally converts via an isospin
transformation. The latter option has been added because some nuclear PDF
fitting groups (which assume isospin symmetry) provide the nuclear PDFs in the
format of average bound proton PDFs.

Taking as an example the scattering of neutrinos with H2O molecules, the total
cross section is given by

$\displaystyle\sigma^{\text{H}_{2}\text{O}}_{\nu}=2\sigma^{p}_{\nu}+Z\sigma^{p/O}_{\nu}+(A-Z)\sigma^{n/O}_{\nu}\,,$
(5)

where $\sigma^{p}_{\nu}$, $\sigma^{p/O}_{\nu}$, and $\sigma^{n/O}_{\nu}$ are
the cross sections for free protons, bound protons and bound neutrons,
respectively, and $Z=A-Z=8$ for oxygen. In this case one has to perform three
different runs: The first run using free protons, the second using bound
protons, and the third using bound neutrons. For both the bound protons and
neutrons one should use nuclear PDFs. The final showered result is then given
by combining these three runs according to the above equation.

When considering scattering on a single nucleus (such as oxygen), one could
generate events using a PDF which is the appropriate admixture of protons and
neutrons in the target nucleus. This would then require two instances of the
parton shower – one for the proton and one for the neutron – that one selects
event by event with the probability determined by the relative fraction of the
PDFs for protons and neutrons in the nucleus. For an extension of the PYTHIA
8 Monte Carlo event generator that enables the simulation of collisions
between a generic hadron beam on a generic nuclear target see Ref.
Helenius:2024vdj . That work combines the extension of PYTHIA 8 to deal with
heavy ion collisions Bierlich:2018xfw , and the extension to collisions of a
varying hadron beam on a proton target Sjostrand:2021dal .

###  2.3 Variable neutrino flux

By default, we consider a monochromatic incoming lepton flux. To account for
the typical environment of a neutrino-induced DIS process our new
implementation additionally provides an option for a realistic neutrino flux.
The user can implement a realistic flux by modifying the function
pdf_lepton_beam, which is contained in the file lepton_flux.f. If importance
sampling associated with the lepton’s energy fraction is required, the user
can modify the function sample_x_lepton, also contained in the same file. This
function builds the lepton’s energy fraction given a random number.

The correct modeling of such a flux depends on the specific experiment and
goes beyond the scope of this publication. A detailed study for SND@LHC,
FASER$\nu$, and the planned FPF experiments FLArE and FASER$\nu$2, using our
code and framework, will be presented in Ref. RojoDIS .

###  2.4 On the momentum mappings, mass effects and possible extensions to
more complex processes

In Ref. Banfi:2023mhz  we introduced new momentum mappings, focusing on the
fully massless case, and used them to implement a DIS generator in the POWHEG-
BOX-RES framework. A POWHEG-BOX-V2 generator was presented in Ref.
Buonocore:2024pdv , where such mappings have been generalised to account for
an explicit lepton-mass dependence. This mass dependence can be relevant when
studying processes involving $\tau$ leptons for $Q$ values not much higher
than the mass of the $\tau$ lepton, as probed by the FASER$\nu$ and SHiP
experiments. Additionally, the initial-state map of Ref. Buonocore:2024pdv
supports heavy coloured final-state particles. In Ref. Buonocore:2024pdv
there is no dedicated treatment of the collinear singularities associated with
the emissions from a final-state heavy quark. This would have required an
extension of the work of Refs. Barze:2012tt ; Buonocore:2017lry  to the DIS
case. Instead, contributions associated with emissions collinear to a heavy
quark, as well as power-suppressed terms, are included at fixed-order accuracy
as a separate regular contribution, involving potentially large mass
logarithms. Therefore, when the centre-of-mass energy becomes very large
relative to the relevant quark masses - as is the case in (ultra) high-energy
neutrino collisions – the massless QCD calculation, available in both codes,
has to be preferred. Indeed we stress that, even in the massless
approximation, when generating radiation in POWHEG, mass thresholds for the
heavy-quarks are present so that the leading mass-logarithms associated with
collinear final-state emissions are included to all orders. Therefore, in
POWHEG events, radiation with a transverse momentum smaller than the mass of
the emitting quark is vetoed, effectively mimicking a dead cone. Furthermore,
we also stress that even for calculations where final-state quarks or leptons
are treated as massless in the matrix-elements, the generated momenta of the
POWHEG events are reshuffled to include finite masses and that the subsequent
parton shower is fully aware of mass effects, including the correct decays of
$\tau$ leptons.

We also note that, in the massless limit, the maps of Refs. Banfi:2023mhz ;
Buonocore:2024pdv  as well as the handling of final-state radiation are
identical. For initial-state radiation instead, while the kinematic map is the
same, they differ in the definition of the POWHEG hardness parameter away from
the soft and collinear limits. Denoting by $\xi$ and $y$ the energy-fraction
and the cosine of the emission angle and by $\bar{s}$ the centre-of-mass
energy of the underlying Born, the two definitions are given by

$\displaystyle\quad t_{\rm ISR}$
$\displaystyle=\frac{\xi^{2}}{2-\xi(1+y)}\bar{s}(1-y),$ $\displaystyle\text{in
Ref.\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Banfi:2023mhz}{\@@citephrase{(}}{\@@citephrase{)}}}},$ (6)
$\displaystyle\quad t_{\rm ISR}$ $\displaystyle=\frac{\xi^{2}}{2(1-\xi
y)}\bar{s}(1-y),$ $\displaystyle\text{in Ref.\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Buonocore:2024pdv}{\@@citephrase{(}}{\@@citephrase{)}}}}.$
(7)

It is evident that the two definitions are identical in the soft ($\xi\to 0$)
and in the collinear ($y\to 1$) limits. We thus conclude that the two codes
have the same formal accuracy.

The POWHEG-BOX-RES framework, specifically designed to handle hadronic
scattering processes that contain decaying resonances and thus require the
inclusion of radiative corrections not only in the production, but also in the
decay process, is particularly well-suited for extending our approach to other
processes relevant for the phenomenology of hadron-hadron collisions, as well
as including electroweak corrections in DIS. In particular, for processes such
as vector boson fusion or vector boson scattering, that can be modelled as
generalised two-fold DIS processes, the POWHEG-BOX-RES framework is best
suited to handle the two hadronic sub-sectors with a factorised approach. In
this sense, our POWHEG-BOX-RES implementation of the genuine DIS process
provides a stepping stone towards the development of suitable generators for
such more complex hadron-hadron collision processes. It is also more
straightforward to include soft photon emissions connecting the leptonic and
the hadronic sectors of the DIS process in the POWHEG-BOX-RES framework. This
feature will be essential for the inclusion of electroweak corrections in the
generator.

##  3 Fixed-order validation

To validate our new implementation, we perform a comparison with existing
fixed-order predictions for selected DIS processes where a neutrino is
scattering off an oxygen target. Specifically, we compute the quantity

$\displaystyle\sigma^{\text{i/O}}_{\nu}=Z/A\,\sigma^{p/O}_{\nu}+(A-Z)/A\,\sigma^{n/O}_{\nu}\,,$
(8)

which is the per-nucleon cross-section for an (isoscalar) oxygen target.

In our work, we have used the set of nuclear PDFs
`nNNPDF30_nlo_as_0118_p_O16` AbdulKhalek:2022fyi , which is provided in a
variable flavour number scheme ($n_{f}^{\rm max}=5$) and is expressed in terms
of average bound-proton states. We note that top quark contributions to DIS
are expected to be negligible below neutrino energies of about $1$ PeV. If
higher energies are considered, the inclusion of the top-quark contributions
could become relevant for the CC process, see Ref. Garcia:2020jwr . We
generated separately a sample for proton ($p$) and neutron ($n$) targets as
described in Sec. 2.2. The neutron PDF is obtained from the proton one using
isospin relations, as described in A. The central renormalisation $\mu_{R}$
and factorisation scales $\mu_{F}$ are set to the momentum transfer $Q$. Scale
uncertainties are estimated by performing an independent variation of
$\mu_{R}$ and $\mu_{F}$ by a factor $2$ up and down, subject to the constraint
$1/2\leq\mu_{R}/\mu_{F}\leq 2$. We impose a lower cutoff on $Q$ of $Q_{\rm
min}=2.0\leavevmode\nobreak\ \mathrm{GeV}$ which ensures the PDFs and the
strong coupling $\alpha_{\mathrm{s}}$ are never evaluated at scales below
$1.0\leavevmode\nobreak\ \mathrm{GeV}$.

For the masses and widths of the electroweak gauge bosons we start from the
on-shell values given in the PDG ParticleDataGroup:2022pth

$\displaystyle m_{W}^{\rm OS}=80.3770\leavevmode\nobreak\
\mathrm{GeV}\,,\qquad\Gamma_{W}^{\rm OS}=2.085\leavevmode\nobreak\
\mathrm{GeV}\,,$ $\displaystyle m_{Z}^{\rm OS}=91.1876\leavevmode\nobreak\
\mathrm{GeV}\,,\qquad\Gamma_{Z}^{\rm OS}=2.4955\leavevmode\nobreak\
\mathrm{GeV}\,,$ (9)

and convert them to the pole values as described e.g. in Ref. Denner:2019vbn ,
which are then used as input values for the simulations. For the Fermi
constant and the weak mixing angle we use

$G_{F}\\!=\\!1.1663787\\!\times\\!10^{-5}\leavevmode\nobreak\
\mathrm{GeV}^{-2}\\!,\;\;\sin^{2}\theta_{W}=0.2316.$ (10)

The value of electromagnetic coupling $\alpha$ is derived from these
parameters as $\alpha=\sqrt{2}/\pi G_{F}m_{W}^{2}\sin^{2}\theta_{W}$. For the
charged current process this effectively implies the replacement
$\alpha/\sin^{2}\theta_{W}\to G_{F}^{2}m_{W}^{2}$ when evaluating the squared
amplitude. This choice ensures the resummation of the leading universal
electroweak corrections Denner:1991kt . A similar replacement also takes place
for the neutral current process, while additional dependencies on
$\sin^{2}\theta_{W}$ appearing in the squared amplitude are described by our
chosen value of $\sin^{2}\theta_{W}$ (which is fixed to the measured effective
weak mixing angle). This approach provides an accurate normalisation of the
couplings, and ensures that the measured on-shell values of the boson masses
enter the propagators for both the charged and neutral current processes we
are describing.

For the entries of the Cabibbo-Kobayashi-Maskawa matrix we have used

$\displaystyle V_{ud}=V_{cs}=0.97446\,,$ $\displaystyle
V_{us}=V_{cd}=0.22456\,,$ $\displaystyle V_{tb}=1\,,$ (11)

with all other entries zero.

The fixed-order predictions are provided at both NLO and NNLO, and have been
obtained using the implementation from Bertone:2018dse , which relies on APFEL
Bertone:2013vaa  for the computation of the DIS structure functions up to NNLO
vanNeerven:1991nn ; Zijlstra:1991qc ; Zijlstra:1992qd ; Zijlstra:1992kj ;
Moch:1999eb . In each case the same NLO accurate nuclear PDF set specified
above is used. The structure functions have been benchmarked against Hoppet
Salam:2008qg ; Bertone:2024dpm  and the fixed-order predictions have been
cross-checked against predictions from disorder Karlberg:2024hnl .

In the following we denote by LO+PS and NLO+PS predictions at LO and NLO,
respectively, matched to parton shower. For the NLO+PS predictions shown below
we interface our POWHEG-BOX implementation to PYTHIA 8.308 Bierlich:2022pfr ,
with default settings (Monash tune Skands:2014pea ), and we use the simple
shower with fully-local recoil option Cabouat:2017rzi . For the results
presented in this section, QED radiation and hadronization effects are not
included.

We have performed comparisons of cross sections differential with respect to
the DIS variables $Q^{2}$ and $x_{\mathrm{B}}$ with different neutrino
energies for both charged current (CC) and neutral current (NC) processes in
the case of either incoming neutrinos or antineutrinos for the scattering off
an oxygen target at rest, i.e. the reactions $\nu_{e}O\to e^{-}X$,
$\bar{\nu}_{e}O\to e^{+}X$, $\nu_{e}O\to\nu_{e}X$, and
$\bar{\nu}_{e}O\to\bar{\nu}_{e}X$, where $X$ denotes the unresolved hadronic
final state of the DIS reaction. We show explicit results for the selected
processes $\nu_{e}O\to e^{-}X$ and $\nu_{e}O\to\nu_{e}X$ in Fig. 1 and Fig. 2,
respectively. In both cases we consider fixed-target collisions with a
neutrino energy of $E_{\nu}=0.1\leavevmode\nobreak\ \mathrm{PeV}$,
corresponding to a neutrino-nucleon centre-of-mass energy of
$\sqrt{s}=431.74\leavevmode\nobreak\ \mathrm{GeV}$. In Figs. 1 and 2 we show
the differential results with respect to $\ln(Q^{2}/{\rm GeV}^{2})$ (left
panel) and $\ln(x_{\mathrm{B}})$ (right panel) for CC and NC, respectively.
For the LO+PS, NLO+PS and NNLO predictions, we show scale variation
uncertainties, while statistical errors are much smaller and not shown here.

(a)

(b)

Figure 1: Differential cross-section (per-nucleon) for the charged-current
scattering of a neutrino $\nu_{e}$ of energy $E_{\nu}=0.1\leavevmode\nobreak\
\mathrm{PeV}$ on oxygen, with respect to $\ln(Q^{2}/{\rm GeV}^{2})$ (left) and
$\ln(x_{\mathrm{B}})$ (right) at LO+PS (green), NLO+PS (blue), pure NLO
(violet) and NNLO (red). The widths of the bands indicate scale uncertainties
estimated by a 7-point variation of $\mu_{R}$ and $\mu_{F}$ by a factor of two
around the central value $Q$. The lower panels show ratios to the respective
NLO+PS results with $\mu_{R}=\mu_{F}=Q$.

(a)

(b)

Figure 2: Analogous to Fig. 1 for the neutral current process
$\nu_{e}O\to\nu_{e}X$.

Cross-sections with the cut $Q>2\leavevmode\nobreak\ \mathrm{GeV}$ for
$E_{\nu}=0.1$ PeV
---
Process | NLO+PS (pb) | NNLO (pb)
$\nu_{e}O\to e^{-}X$ | $200.68^{+2.87}_{-3.53}\,\text{(scales)}\,^{+2.68}_{-3.29}\,\text{(PDFs)}$ | $197.92^{+1.21}_{-1.02}\,\text{(scales)}$
$\bar{\nu}_{e}O\to e^{+}X$ | $168.32^{+2.73}_{-3.34}\,\text{(scales)}\,^{+2.64}_{-3.34}\,\text{(PDFs)}$ | $165.73^{+1.16}_{-0.99}\,\text{(scales)}$
$\nu_{e}O\to\nu_{e}X$ | $75.97^{+1.25}_{-1.39}\,\text{(scales)}\,^{+0.76}_{-0.91}\,\text{(PDFs)}$ | $74.81^{+0.44}_{-0.41}\,\text{(scales)}$
$\bar{\nu}_{e}O\to\bar{\nu}_{e}X$ | $64.85^{+1.21}_{-1.33}\,\text{(scales)}\,^{+0.78}_{-0.82}\,\text{(PDFs)}$ | $63.75^{+0.42}_{-0.40}\,\text{(scales)}$

Table 1: Total cross-section with the cut $Q>2\leavevmode\nobreak\
\mathrm{GeV}$ for a selection of DIS processes with a (anti-)neutrino of
energy $E_{\nu}=0.1$ PeV at NLO+PS and NNLO accuracy. The quoted uncertainties
are due to scale variation. For the NLO+PS results we also indicate the size
of the PDF uncertainties in the second entry.

Cross-sections with the cut $Q>2\leavevmode\nobreak\ \mathrm{GeV}$ for
$E_{\nu}=1$ PeV
---
Process | NLO+PS (pb) | NNLO (pb)
$\nu_{e}O\to e^{-}X$ | $624.49^{+14.14}_{-16.44}\,\text{(scales)}\,^{+15.26}_{-15.42}\,\text{(PDFs)}$ | $613.42^{+5.02}_{-3.70}\,\text{(scales)}$
$\bar{\nu}_{e}O\to e^{+}X$ | $598.05^{+14.00}_{-16.33}\,\text{(scales)}\,^{+15.81}_{-15.90}\,\text{(PDFs)}$ | $587.09^{+4.99}_{-3.68}\,\text{(scales)}$
$\nu_{e}O\to\nu_{e}X$ | $258.59^{+6.48}_{-7.11}\,\text{(scales)}\,^{+5.67}_{-5.69}\,\text{(PDFs)}$ | $253.61^{+2.06}_{-1.61}\,\text{(scales)}$
$\bar{\nu}_{e}O\to\bar{\nu}_{e}X$ | $248.73^{+6.43}_{-7.07}\,\text{(scales)}\,^{+5.82}_{-5.58}\,\text{(PDFs)}$ | $243.78^{+2.05}_{-1.60}\,\text{(scales)}$

Table 2: Analogous to Tab. 1, now for $E_{\nu}=1$ PeV.

We observe that at low-to-moderate values of $Q^{2}$, within the given scale
uncertainties, the fixed-order NLO predictions agree with the NLO+PS results
and are very similar to the LO+PS results. Obviously, the impact of higher-
order corrections is small on this observable. For the Bjorken variable we
find agreement between the NLO and the NLO+PS results, as expected for this
inclusive quantity. Technically we expect the agreement between NLO and NLO+PS
to be near-perfect, as the shower without QED radiation preserves the lepton
momenta. However, as discussed in Sec. 2.4, the POWHEG-BOX performs a small
momentum reshuffling to account for the finite quark and lepton masses, and
additionally, as was discussed in Sec. 2.1, at event level the nucleon mass is
restored. This reshuffling has a tiny impact on the $Q^{2}$ and
$x_{\mathrm{B}}$ distributions, as was also discussed in Ref. Banfi:2023mhz .

It is worth noticing that the NLO+PS result is not always contained with in
the scale variation band of the LO+PS result. The perturbative uncertainties
of the LO+PS result are not expected to be fully covered by a standard scale
variation, as at this order only $\mu_{F}$ can be varied, while $\mu_{R}$ does
not even enter. On the other hand, we see that the NNLO prediction is fully
contained within the scale variation band of the NLO+PS prediction, thereby
establishing confidence in the reliability of our prediction.

In addition to the differential validation, we also report results for the
per-nucleon cross section, with a cut $Q\geq 2$ GeV, obtained up to NNLO
accuracy in Tab. 1 for $E_{\nu}=0.1$ PeV and Tab. 2 for $E_{\nu}=1$ PeV. The
results are given for a selection of processes and (anti)-neutrino energies.
The central prediction and the uncertainty due to scale variations are shown
in each case. It has been checked that the NLO entries obtained with this
generator (labelled as NLO+PS) reproduce exactly, including scale variations,
the NLO results based on the structure function computation. For that reason
we only show the NLO+PS results. We have additionally reported the
uncertainties due to the nuclear PDFs computed at NLO. Typically these
uncertainties are in the range of $(1-2)\%$ and are similar in size to those
of the scale uncertainties at NLO. Finally, we note that the structure
functions are non-zero below $Q_{\rm min}$ (and hence so is the cross-
section), but the description of this region goes beyond the applicability of
collinear factorisation. Alternative (data-driven) approaches exist to
describe the low-$Q$ region, see for example Bodek:2002vp ; Bodek:2003wd ;
Bodek:2004pc ; Bodek:2010km ; Bodek:2021bde  and, more recently, Ref.
Candido:2023utz .

##  4 Phenomenological results

As highlighted in Sec. 1, a major advantage of the NLO+PS simulation over the
NLO predictions is that they enable a fully exclusive simulation of final-
state radiation while retaining the NLO accuracy of the hard scattering
process. In this section we consider full particle level predictions obtained
with our NLO+PS generator interfaced to PYTHIA 8. We use the same PDFs,
scale settings, and electroweak input parameters specified in Sec. 3, but we
also include QED radiation and hadronization effects in the PYTHIA 8
simulation, which allow us to provide predictions for the production of
hadrons, and to investigate properties of their distributions. We note that
the inclusion of QED corrections can have important consequences for the
description of charged-lepton based observables (see the recent discussion in
Ref. Plestid:2024bva ), and that the leading corrections are naturally
included (and resummed) by the parton shower in the following. Specifically,
we consider fixed-target collisions on oxygen atoms for electron neutrinos
with energies of $0.1$ and $1\leavevmode\nobreak\ \mathrm{PeV}$, which are
primarily relevant for analyses aiming to measure the flux of cosmic
neutrinos.

###  4.1 Particle multiplicities

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3:  Charged particle multiplicity distribution (left) and multiplicity
ratio between charged and neutral particles (right) obtained at NLO+PS (blue)
and LO+PS (green) accuracy for neutrino induced CC DIS, panels (a),(b), and NC
DIS panels (c),(d), on an oxygen target with a neutrino energy of
$E_{\nu}=0.1\leavevmode\nobreak\ \mathrm{PeV}$, and for CC DIS with
$E_{\nu}=1\leavevmode\nobreak\ \mathrm{PeV}$, panels (e),(f). The widths of
the bands indicate scale uncertainties estimated by a 7-point variation of
$\mu_{R}$ and $\mu_{F}$ by a factor of two around the central value $Q$. The
lower panels show ratios to the respective NLO+PS results with
$\mu_{R}=\mu_{F}=Q$.

Water-based detector concepts rely on observing the Cherenkov radiation
pattern generated by charged particles in the detector volume. An accurate
modelling of particle multiplicities in such scattering events is therefore
critical. Charged particle multiplicities, as well as the ratio of charged to
neutral particle multiplicities are shown in Fig. 3 for $\nu_{e}$-induced CC
and NC DIS at $E_{\nu}=0.1\leavevmode\nobreak\ \mathrm{PeV}$, (upper and
middle panels), and CC DIS at $E_{\nu}=1\leavevmode\nobreak\ \mathrm{PeV}$
(lower panels). The multiplicity distribution at
$E_{\nu}=0.1\leavevmode\nobreak\ \mathrm{PeV}$ peaks for a number of charged
particles, $n_{\mathrm{ch}}$, of about 18 in both the CC and NC cases. At
$E_{\nu}=1\leavevmode\nobreak\ \mathrm{PeV}$ the peak is shifted to around
$n_{\mathrm{ch}}=22$.

As a consequence of charge conservation, an odd (even) number of charged
particles is generated in CC neutrino scattering off protons (neutrons).
Furthermore, because of the different flavour composition and associated PDFs
of these two types of target particles, the absolute scattering rate is
different for CC on a proton and on a neutron. The combination of these
effects leads to the observed “oscillatory” behaviour for the
$n_{\mathrm{ch}}$ distributions. This feature is slightly less pronounced at
higher neutrino energies, as the contribution from PDFs at smaller values of
$x$, where the isospin asymmetric contribution of valence quarks is less
important, becomes more relevant. We note that the ratio
$n_{\mathrm{ch}}/n_{\mathrm{neut}}$ peaks at smaller values for the NC
process. Generally, we observe a reduction of scale uncertainty when including
NLO corrections and considerable shape changes induced by NLO effects which
are outside the LO scale uncertainty band, both for the charged particle
multiplicities, as well as the ratios. When considering higher neutrino
energies we notice that the charged particle multiplicity increases, as
expected, and that the NLO corrections are becoming yet more pronounced and
the theoretical uncertainty stemming from scale variation increases.

It is interesting to note that the centre-of-mass energies considered here are
comparable to those of the HERA collider. Our NLO+PS implementation opens up
the opportunity for the re-tuning of event generators such as PYTHIA 8,
which could be relevant given the large impact of NLO+PS corrections on
particle multiplicities.

###  4.2 Energy-based distributions

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Similar to Fig. 3, but for the energy of the leading charged
particle, $E_{1,\rm chg}$, (left) and the mean charged particle energy
$\langle E_{\rm chg}\rangle$ (right).

In Fig. 4 we compare the predictions for the energy of the hardest charged
particle, $E_{1,\rm chg}$, and the mean charged particle energy, $\langle
E_{\rm chg}\rangle$, as predicted at LO+PS and NLO+PS accuracy. We notice that
these energy distributions are genuinely different for the CC and NC cases.
This is due to the fact that in the CC case the outgoing lepton contributes to
both distributions, while this is not the case for NC. For this reason, NLO
corrections turn out to be moderate in the CC case, which is dominated by the
lepton kinematics, but considerable for NC. We note that, generally, for the
determination of $E_{1,\rm chg}$ and $\langle E_{\rm chg}\rangle$ all charged
particles (i.e. hadrons and leptons) are taken into account. If, however, the
outgoing charged lepton is not included in the definition of $E_{1,\rm chg}$
or $\langle E_{\rm chg}\rangle$ in the CC case, it is observed that the
resultant distributions (and the behaviour of the NLO corrections) are similar
to those of the NC case. Like for the case of the particle multiplicity, the
LO scale uncertainty band significantly underestimates the size of higher-
order effects as it does not overlap with the NLO band in the majority of the
phase space. When going to higher energies (plots (e) and (f)), the peaks of
the distributions move accordingly and we find that, as for particle
multiplicities, NLO corrections become more pronounced.

###  4.3 Charm production

(a)

(b)

(c)

(d)

Figure 5:  $D$-meson energy distributions at NLO+PS (blue) and LO+PS (green)
accuracy for neutrino induced CC (left) and NC (right) DIS with a neutrino
energy of $E_{\nu}=0.1\leavevmode\nobreak\ \mathrm{PeV}$, panels (a),(b), and
$E_{\nu}=1\leavevmode\nobreak\ \mathrm{PeV}$, panels (c),(d). The widths of
the bands indicate scale uncertainties estimated by a 7-point variation of
$\mu_{R}$ and $\mu_{F}$ by a factor of two around the central value $Q$. The
lower panels show ratios to the respective NLO+PS results with
$\mu_{R}=\mu_{F}=Q$.

It is also interesting to investigate the effect of QCD corrections on
$D$-meson distributions. This is relevant as, through semi-leptonic decays,
$D$-mesons provide a source of energetic muons which can mimic a starting
track signature similar to that arising from muon-neutrino induced CC. As
discussed in Sec. 2.4, despite being based on a purely massless calculation,
once interfaced to a parton shower, our event generator is well suited to
describe DIS processes involving heavy quarks if their mass is much smaller
than $Q$, as considered in this section. In fact, at the considered neutrino
energies, the typical $Q^{2}$ value which dominates the cross-section is far
in excess of the charm quark mass (i.e. $|Q^{2}|\gg m_{c}^{2}$), as shown in
Fig. 1(a). In such a kinematic regime a massless approach to describing the
scattering process is the appropriate one, and ensures a resummation of the
logarithmically enhanced terms in both the initial and final-state.

We consider here the production of stable $D$-mesons at LO+PS and NLO+PS
accuracy, where the $D$-mesons are produced using the hadronization feature of
PYTHIA 8. In Fig. 5 we present the distribution of the $D$-meson energy,
$E_{D}$, in the CC and NC cases, respectively. We find that in the CC case NLO
corrections are moderate for low energies, but become large for high values of
$E_{D}$, where the cross section peaks.

The CC case is dominated by scattering off $d$\- and $s$-quark distributions,
while NC involves primarily a $c$-PDF, which is generated perturbatively and
has a large factorization scale dependence. For this reason, for NC DIS the
scale uncertainties are larger than in the CC case. These are substantially
reduced at NLO. In each case, the NLO corrections are essential for a
reasonable description of the shape of the energy distribution.

##  5 Conclusions

This work presents a number of extensions to the simulation of neutral- and
charged-current deep inelastic scattering (DIS) Banfi:2023mhz  in the POWHEG-
BOX-RES. First, the code has been extended to accommodate an incoming neutrino
beam. Second, the incoming lepton is no longer required to be monochromatic,
as in standard high-energy DIS experiments. Instead, any incoming lepton flux
can be included. Moreover, an option is provided to straightforwardly account
for the kinematics of fixed-target experiments. Furthermore, more flexible
options for the nuclear targets are now supported.

With the new implementation we have provided sample results for fiducial
cross-sections, standard DIS variables, as well as neutral and charged
particle distributions for various neutrino-induced DIS processes. In our
sample numerical analyses we put a particular focus on the kinematic regime
relevant for the investigation of cosmic neutrinos with the IceCube detector.
We note, however, that our program is not restricted to this application, but
can be employed for the simulation of any neutrino-induced DIS process. In
general, we find that an NLO+PS simulation is necessary to achieve theory
uncertainties below approximately 10%.

The code, along with the new features discussed in this article, is publicly
available via the POWHEG-BOX-RES repository. The reliance on the POWHEG-BOX-
RES framework, which is well-suited for describing complex reactions involving
multiple competing and interfering sub-processes, will enable us to further
improve the description of hadron-collider processes such as vector boson
scattering and vector boson fusion, going beyond what is already available in
POWHEG-BOX-V2. These reactions can be described as (generalized) two-fold DIS
processes, and are highly relevant for the phenomenology of the Large Hadron
Collider. Additionally, our approach paves the way for the simulation of
electroweak corrections in DIS consistently accounting for photon radiation in
the hadronic and leptonic sectors.

## Acknowledgments

We are grateful to Luca Buonocore, Giovanni Limatola, Paolo Nason and
Francesco Tramontano for discussions and to Georg Raffelt for providing useful
references. In particular, we thank Paolo and Francesco for discussions on the
treatment of collisions involving heavy nuclei. We also acknowledge
stimulating discussions with Andrea Banfi during early stages of this work. We
are also indebted to Melissa van Beekveld, Eva Groenendijk, Peter Krack, Juan
Rojo, and Valentina Schutze Sanchez for having triggered this project and
tested a pre-release version of our code. Finally, we are grateful to Alfonso
García Soto for advice and discussions related to experimentally motivated
observables. The work of BJ was supported by the German Research Foundation
(DFG) through the Research Unit FOR 2926. GZ would like to thank CERN for
hospitality while this work was being finalized.

##  Appendix A DIS process selection

In this appendix we summarise inputs that can be used to select the process
and settings in the powheg.input file, which are specific to the DIS case.

####  Lepton beam.

The flavour of the incoming lepton must be specified using `ih1 int`, where
the integer number int is the identifier of the desired lepton in the Particle
Data Group numbering convention ParticleDataGroup:2022pth .

The energy of the lepton beam must be specified using ebeam1 double with a
double-precision number double. By default the code assumes a fixed lepton
energy. To use a variable flux add the option
`fixed_lepton_beam 0`, (see Sec. 2.3 for more details). A variable flux should
be provided in terms of a boost-invariant energy fraction, of the lepton beam
with repect to the maximum energy available.

####  Hadron beam/target.

The selection of the nucleon type in the powheg.input file must be chosen by
setting the value of int in `ih2 int`. We currently support protons and
neutrons. To that end the following options are available:

  1. `ih2 1 #proton target, input proton PDF`

  2. `ih2 2 #neutron target, input proton PDF`

  3. `ih2 22 #neutron target, input neutron PDF`

Depending on the selection for ih2, the PDF specified via the entry lhans2
according to the numbering scheme of the LHAPDF repository Buckley:2014ana  is
interpreted either as a proton or a neutron PDF.111lhans1 must be set to the
same value as lhans2, even if not used, in case the PDF implementation of the
running of the QCD coupling constant (alphas_from_pdf 1) is to be used.

The energy of the hadron is selected via the mandatory entry ebeam2 double. By
default, the code assumes that the hadron beam is massless, with a
longitudinal momentum equal to its energy. For fixed-target collisions, one
has to add the option

    fixed_target 1

In this case the value of the entry for ebeam2 is interpreted as the mass of
the nucleon (i.e. proton or neutron).

####  Hard process selection.

Both CC and NC processes can be simulated within our framework. To select the
desired channel (for a given type of lepton beam, as specified by the value of
ih1), one can use the following option

    channel_type int

with int=3 for CC, and int=4 for NC. In case of charged-lepton induced NC DIS,
the boson exchanged in the $t$-channel has to be specified using `vtype` with

  1. 1.

`vtype 1 # photon exchange only`

  2. 2.

`vtype 2 # Z echange only`

  3. 3.

`vtype 3 # photon+Z exchange`

####  Generation cuts.

The user must specify cuts on the DIS invariants $Q^{2}$, $x_{\mathrm{B}}$ and
$y_{\mathrm{DIS}}=Q^{2}/(x_{\mathrm{B}}S)$, with $S=(P_{1}+P_{2})^{2}$. The
values of Qmin and Qmax are supposed to be provided in units of GeV. For
example, to probe all the available phase space, one should set

    Qmin 1d0
    Qmax 1d8
    xmin 0d0
    xmax 1d0
    ymin 0d0
    ymax 1d0

where Qmax has been set to a value much larger than the center-of-mass energy.
We stress that Qmin=1 GeV is the lowest value accepted by the code, since the
validity of a perturbative QCD approach to describe the cross section is no
longer guaranteed for small $Q^{2}$.

We note that it is possible to fix up to 2 of these variables by setting the
minimum and maximum values equal to each other. In any case, the code will
never generate events outside the physically allowed bounds.

####  Final-state particles masses.

Notice that all particles entering the hard process are treated as massless in
our NLO calculation. As in most POWHEG-BOX implementations, a small
reshuffling of the momenta can be applied when generating events, so as to
give a finite mass to all the final-state massive particles. The mass of the
charged leptons and of the heavy quarks can be specified, e.g. using

    electron_mass 0.51099891d-3
    muon_mass 0.1056583668d0
    tauon_mass 1.77684d0
    charm_mass 1.5d0
    bottom_mass 4.5d0

The numbers chosen in this example correspond to the default values. More
comments on how mass effects are approximately included in our generator are
given in 2.4.

## References

  * (1) B. Batell, M. Pospelov, A. Ritz, Exploring Portals to a Hidden Sector Through Fixed Targets, Phys. Rev. D 80 (2009) 095024.  arXiv:0906.5614, doi:10.1103/PhysRevD.80.095024.
  * (2) M. Ackermann, et al., High-energy and ultra-high-energy neutrinos: A Snowmass white paper, JHEAp 36 (2022) 55–110.  arXiv:2203.08096, doi:10.1016/j.jheap.2022.08.001.
  * (3) M. Ageron, et al., ANTARES: the first undersea neutrino telescope, Nucl. Instrum. Meth. A 656 (2011) 11–38.  arXiv:1104.1607, doi:10.1016/j.nima.2011.06.103.
  * (4) V. Aynutdinov, et al., Search for a diffuse flux of high-energy extraterrestrial neutrinos with the nt200 neutrino telescope, Astropart. Phys. 25 (2006) 140–150.  arXiv:astro-ph/0508675, doi:10.1016/j.astropartphys.2005.12.005.
  * (5) M. G. Aartsen, et al., The IceCube Neutrino Observatory: Instrumentation and Online Systems, JINST 12 (03) (2017) P03012.  arXiv:1612.05093, doi:10.1088/1748-0221/12/03/P03012.
  * (6) S. Adrian-Martinez, et al., Letter of intent for KM3NeT 2.0, J. Phys. G 43 (8) (2016) 084001.  arXiv:1601.07459, doi:10.1088/0954-3899/43/8/084001.
  * (7) D. Mason, et al., Measurement of the Nucleon Strange-Antistrange Asymmetry at Next-to-Leading Order in QCD from NuTeV Dimuon Data, Phys. Rev. Lett. 99 (2007) 192001.  doi:10.1103/PhysRevLett.99.192001.
  * (8) G. Acampora, et al., SND@LHC: the scattering and neutrino detector at the LHC, JINST 19 (05) (2024) P05067.  arXiv:2210.02784, doi:10.1088/1748-0221/19/05/P05067.
  * (9) R. Albanese, et al., Observation of Collider Muon Neutrinos with the SND@LHC Experiment, Phys. Rev. Lett. 131 (3) (2023) 031802.  arXiv:2305.09383, doi:10.1103/PhysRevLett.131.031802.
  * (10) M. Anelli, et al., A facility to Search for Hidden Particles (SHiP) at the CERN SPS (4 2015).  arXiv:1504.04956.
  * (11) H. Abreu, et al., Detecting and Studying High-Energy Collider Neutrinos with FASER at the LHC, Eur. Phys. J. C 80 (1) (2020) 61.  arXiv:1908.02310, doi:10.1140/epjc/s10052-020-7631-5.
  * (12) H. Abreu, et al., Technical Proposal: FASERnu (1 2020).  arXiv:2001.03073.
  * (13) H. Abreu, et al., First Direct Observation of Collider Neutrinos with FASER at the LHC, Phys. Rev. Lett. 131 (3) (2023) 031801.  arXiv:2303.14185, doi:10.1103/PhysRevLett.131.031801.
  * (14) L. A. Anchordoqui, et al., The Forward Physics Facility: Sites, experiments, and physics potential, Phys. Rept. 968 (2022) 1–50.  arXiv:2109.10905, doi:10.1016/j.physrep.2022.04.004.
  * (15) J. L. Feng, et al., The Forward Physics Facility at the High-Luminosity LHC, J. Phys. G 50 (3) (2023) 030501.  arXiv:2203.05090, doi:10.1088/1361-6471/ac865e.
  * (16) J. C. Collins, D. E. Soper, G. F. Sterman, Factorization of Hard Processes in QCD, Adv. Ser. Direct. High Energy Phys. 5 (1989) 1–91.  arXiv:hep-ph/0409313, doi:10.1142/9789814503266_0001.
  * (17) J. J. Ethier, E. R. Nocera, Parton Distributions in Nucleons and Nuclei, Ann. Rev. Nucl. Part. Sci. 70 (2020) 43–76.  arXiv:2001.07722, doi:10.1146/annurev-nucl-011720-042725.
  * (18) A. Cooper-Sarkar, P. Mertsch, S. Sarkar, The high energy neutrino cross-section in the Standard Model and its uncertainty, JHEP 08 (2011) 042.  arXiv:1106.3723, doi:10.1007/JHEP08(2011)042.
  * (19) J. Sanchez Guillen, J. Miramontes, M. Miramontes, G. Parente, O. A. Sampayo, Next-to-leading order analysis of the deep inelastic R = sigma-L / sigma-total, Nucl. Phys. B 353 (1991) 337–345.  doi:10.1016/0550-3213(91)90340-4.
  * (20) W. L. van Neerven, E. B. Zijlstra, Order alpha-s**2 contributions to the deep inelastic Wilson coefficient, Phys. Lett. B 272 (1991) 127–133.  doi:10.1016/0370-2693(91)91024-P.
  * (21) E. B. Zijlstra, W. L. van Neerven, Contribution of the second order gluonic Wilson coefficient to the deep inelastic structure function, Phys. Lett. B 273 (1991) 476–482.  doi:10.1016/0370-2693(91)90301-6.
  * (22) E. B. Zijlstra, W. L. van Neerven, Order alpha-s**2 QCD corrections to the deep inelastic proton structure functions F2 and F(L), Nucl. Phys. B 383 (1992) 525–574.  doi:10.1016/0550-3213(92)90087-R.
  * (23) E. B. Zijlstra, W. L. van Neerven, Order alpha-s**2 correction to the structure function F3 (x, Q**2) in deep inelastic neutrino - hadron scattering, Phys. Lett. B 297 (1992) 377–384.  doi:10.1016/0370-2693(92)91277-G.
  * (24) W. L. van Neerven, A. Vogt, NNLO evolution of deep inelastic structure functions: The Nonsinglet case, Nucl. Phys. B 568 (2000) 263–286.  arXiv:hep-ph/9907472, doi:10.1016/S0550-3213(99)00668-9.
  * (25) W. L. van Neerven, A. Vogt, NNLO evolution of deep inelastic structure functions: The Singlet case, Nucl. Phys. B 588 (2000) 345–373.  arXiv:hep-ph/0006154, doi:10.1016/S0550-3213(00)00480-6.
  * (26) S. Moch, J. A. M. Vermaseren, Deep inelastic structure functions at two loops, Nucl. Phys. B 573 (2000) 853–907.  arXiv:hep-ph/9912355, doi:10.1016/S0550-3213(00)00045-6.
  * (27) S. Moch, J. A. M. Vermaseren, A. Vogt, The Longitudinal structure function at the third order, Phys. Lett. B 606 (2005) 123–129.  arXiv:hep-ph/0411112, doi:10.1016/j.physletb.2004.11.063.
  * (28) J. A. M. Vermaseren, A. Vogt, S. Moch, The Third-order QCD corrections to deep-inelastic scattering by photon exchange, Nucl. Phys. B 724 (2005) 3–182.  arXiv:hep-ph/0504242, doi:10.1016/j.nuclphysb.2005.06.020.
  * (29) A. Vogt, S. Moch, J. Vermaseren, Third-order QCD results on form factors and coefficient functions, Nucl. Phys. B Proc. Suppl. 160 (2006) 44–50.  arXiv:hep-ph/0608307, doi:10.1016/j.nuclphysbps.2006.09.101.
  * (30) S. Moch, M. Rogal, A. Vogt, Differences between charged-current coefficient functions, Nucl. Phys. B 790 (2008) 317–335.  arXiv:0708.3731, doi:10.1016/j.nuclphysb.2007.09.022.
  * (31) J. Davies, A. Vogt, S. Moch, J. A. M. Vermaseren, Non-singlet coefficient functions for charged-current deep-inelastic scattering to the third order in QCD, PoS DIS2016 (2016) 059.  arXiv:1606.08907, doi:10.22323/1.265.0059.
  * (32) J. Blümlein, P. Marquard, C. Schneider, K. Schönwald, The massless three-loop Wilson coefficients for the deep-inelastic structure functions F2, FL, xF3 and g1, JHEP 11 (2022) 156.  arXiv:2208.14325, doi:10.1007/JHEP11(2022)156.
  * (33) R. Gandhi, C. Quigg, M. H. Reno, I. Sarcevic, Neutrino interactions at ultrahigh-energies, Phys. Rev. D 58 (1998) 093009.  arXiv:hep-ph/9807264, doi:10.1103/PhysRevD.58.093009.
  * (34) M. Gluck, S. Kretzer, E. Reya, Dynamical QCD predictions for ultrahigh-energy neutrino cross-sections, Astropart. Phys. 11 (1999) 327–334.  arXiv:astro-ph/9809273, doi:10.1016/S0927-6505(99)00006-7.
  * (35) A. Cooper-Sarkar, S. Sarkar, Predictions for high energy neutrino cross-sections from the ZEUS global PDF fits, JHEP 01 (2008) 075.  arXiv:0710.5303, doi:10.1088/1126-6708/2008/01/075.
  * (36) A. Connolly, R. S. Thorne, D. Waters, Calculation of High Energy Neutrino-Nucleon Cross Sections and Uncertainties Using the MSTW Parton Distribution Functions and Implications for Future Experiments, Phys. Rev. D 83 (2011) 113009.  arXiv:1102.0691, doi:10.1103/PhysRevD.83.113009.
  * (37) V. Bertone, R. Gauld, J. Rojo, Neutrino Telescopes as QCD Microscopes, JHEP 01 (2019) 217.  arXiv:1808.02034, doi:10.1007/JHEP01(2019)217.
  * (38) K. Xie, J. Gao, T. J. Hobbs, D. R. Stump, C. P. Yuan, High-energy neutrino deep inelastic scattering cross sections, Phys. Rev. D 109 (11) (2024) 113001\.  arXiv:2303.13607, doi:10.1103/PhysRevD.109.113001.
  * (39) D. Seckel, Neutrino photon reactions in astrophysics and cosmology, Phys. Rev. Lett. 80 (1998) 900–903.  arXiv:hep-ph/9709290, doi:10.1103/PhysRevLett.80.900.
  * (40) I. Alikhanov, Hidden Glashow resonance in neutrino–nucleus collisions, Phys. Lett. B 756 (2016) 247–253.  arXiv:1503.08817, doi:10.1016/j.physletb.2016.03.009.
  * (41) R. Gauld, Precise predictions for multi-TeV and PeV energy neutrino scattering rates, Phys. Rev. D 100 (9) (2019) 091301.  arXiv:1905.03792, doi:10.1103/PhysRevD.100.091301.
  * (42) B. Zhou, J. F. Beacom, Neutrino-nucleus cross sections for W-boson and trident production, Phys. Rev. D 101 (3) (2020) 036011.  arXiv:1910.08090, doi:10.1103/PhysRevD.101.036011.
  * (43) B. Zhou, J. F. Beacom, W-boson and trident production in TeV–PeV neutrino observatories, Phys. Rev. D 101 (3) (2020) 036010.  arXiv:1910.10720, doi:10.1103/PhysRevD.101.036010.
  * (44) K. Xie, B. Zhou, T. J. Hobbs, The photon content of the neutron, JHEP 04 (2024) 022.  arXiv:2305.10497, doi:10.1007/JHEP04(2024)022.
  * (45) J. Jalilian-Marian, Enhancement and suppression of the neutrino nucleon total cross-section at ultrahigh-energies, Phys. Rev. D 68 (2003) 054005, [Erratum: Phys.Rev.D 70, 079903 (2004)].  arXiv:hep-ph/0301238, doi:10.1103/PhysRevD.68.054005.
  * (46) R. Fiore, L. L. Jenkovszky, A. V. Kotikov, F. Paccanoni, A. Papa, Asymptotic neutrino-nucleon cross section and saturation effects, Phys. Rev. D 73 (2006) 053012.  arXiv:hep-ph/0512259, doi:10.1103/PhysRevD.73.053012.
  * (47) M. M. Block, L. Durand, P. Ha, D. W. McKay, Implications of a Froissart bound saturation of $\gamma$*-p deep inelastic scattering. II. Ultrahigh energy neutrino interactions, Phys. Rev. D 88 (1) (2013) 013003.  arXiv:1302.6127, doi:10.1103/PhysRevD.88.013003.
  * (48) J. L. Albacete, J. I. Illana, A. Soto-Ontoso, Neutrino-nucleon cross section at ultrahigh energy and its astrophysical implications, Phys. Rev. D 92 (1) (2015) 014027.  arXiv:1505.06583, doi:10.1103/PhysRevD.92.014027.
  * (49) V. P. Goncalves, D. R. Gratieri, Investigating the effects of the QCD dynamics in the neutrino absorption by the Earth’s interior at ultrahigh energies, Phys. Rev. D 92 (11) (2015) 113007.  arXiv:1510.03186, doi:10.1103/PhysRevD.92.113007.
  * (50) C. A. Argüelles, F. Halzen, L. Wille, M. Kroll, M. H. Reno, High-energy behavior of photon, neutrino, and proton cross sections, Phys. Rev. D 92 (7) (2015) 074040.  arXiv:1504.06639, doi:10.1103/PhysRevD.92.074040.
  * (51) P. Nason, A New method for combining NLO QCD with shower Monte Carlo algorithms, JHEP 11 (2004) 040.  arXiv:hep-ph/0409146, doi:10.1088/1126-6708/2004/11/040.
  * (52) S. Frixione, P. Nason, C. Oleari, Matching NLO QCD computations with Parton Shower simulations: the POWHEG method, JHEP 11 (2007) 070.  arXiv:0709.2092, doi:10.1088/1126-6708/2007/11/070.
  * (53) A. Banfi, S. Ferrario Ravasio, B. Jäger, A. Karlberg, F. Reichenbach, G. Zanderighi, A POWHEG generator for deep inelastic scattering, JHEP 02 (2024) 023.  arXiv:2309.02127, doi:10.1007/JHEP02(2024)023.
  * (54) T. Ježo, P. Nason, On the Treatment of Resonances in Next-to-Leading Order Calculations Matched to a Parton Shower, JHEP 12 (2015) 065.  arXiv:1509.09071, doi:10.1007/JHEP12(2015)065.
  * (55) L. Buonocore, G. Limatola, P. Nason, F. Tramontano, An event generator for Lepton-Hadron Deep Inelastic Scattering at NLO+PS with POWHEG including mass effects (6 2024).  arXiv:2406.05115.
  * (56) S. Alioli, P. Nason, C. Oleari, E. Re, A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX, JHEP 06 (2010) 043.  arXiv:1002.2581, doi:10.1007/JHEP06(2010)043.
  * (57) S. Frixione, Z. Kunszt, A. Signer, Three jet cross-sections to next-to-leading order, Nucl. Phys. B 467 (1996) 399–442.  arXiv:hep-ph/9512328, doi:10.1016/0550-3213(96)00110-1.
  * (58) I. Borsa, B. Jäger, Parton-shower effects in polarized deep inelastic scattering (4 2024).  arXiv:2404.07702.
  * (59) I. Helenius, M. Utheim, Hadron-ion collisions in Pythia and the vector-meson dominance model for photoproduction (6 2024).  arXiv:2406.10403.
  * (60) C. Bierlich, G. Gustafson, L. Lönnblad, H. Shah, The Angantyr model for Heavy-Ion Collisions in PYTHIA8, JHEP 10 (2018) 134.  arXiv:1806.10820, doi:10.1007/JHEP10(2018)134.
  * (61) T. Sjöstrand, M. Utheim, Hadron interactions for arbitrary energies and species, with applications to cosmic rays, Eur. Phys. J. C 82 (1) (2022) 21.  arXiv:2108.03481, doi:10.1140/epjc/s10052-021-09953-5.
  * (62) M. van Beekveld, S. Ferrario Ravasio, E. Groenendijk, P. Krack, J. Rojo, V. Schutze Sanchez, A Phenomenological Analysis of LHC Neutrino Scattering at NLO Accuracy Matched to Parton Showers (In preparation).
  * (63) L. Barze, G. Montagna, P. Nason, O. Nicrosini, F. Piccinini, Implementation of electroweak corrections in the POWHEG BOX: single W production, JHEP 04 (2012) 037.  arXiv:1202.0465, doi:10.1007/JHEP04(2012)037.
  * (64) L. Buonocore, P. Nason, F. Tramontano, Heavy quark radiation in NLO+PS POWHEG generators, Eur. Phys. J. C 78 (2) (2018) 151.  arXiv:1711.06281, doi:10.1140/epjc/s10052-018-5638-y.
  * (65) R. Abdul Khalek, R. Gauld, T. Giani, E. R. Nocera, T. R. Rabemananjara, J. Rojo, nnnpdf3.0: evidence for a modified partonic structure in heavy nuclei, Eur. Phys. J. C 82 (6) (2022) 507.  arXiv:2201.12363, doi:10.1140/epjc/s10052-022-10417-7.
  * (66) A. Garcia, R. Gauld, A. Heijboer, J. Rojo, Complete predictions for high-energy neutrino propagation in matter, JCAP 09 (2020) 025.  arXiv:2004.04756, doi:10.1088/1475-7516/2020/09/025.
  * (67) R. L. Workman, et al., Review of Particle Physics, PTEP 2022 (2022) 083C01.  doi:10.1093/ptep/ptac097.
  * (68) A. Denner, S. Dittmaier, Electroweak Radiative Corrections for Collider Physics, Phys. Rept. 864 (2020) 1–163.  arXiv:1912.06823, doi:10.1016/j.physrep.2020.04.001.
  * (69) A. Denner, Techniques for calculation of electroweak radiative corrections at the one loop level and results for W physics at LEP-200, Fortsch. Phys. 41 (1993) 307–420.  arXiv:0709.1075, doi:10.1002/prop.2190410402.
  * (70) V. Bertone, S. Carrazza, J. Rojo, APFEL: A PDF Evolution Library with QED corrections, Comput. Phys. Commun. 185 (2014) 1647–1668.  arXiv:1310.1394, doi:10.1016/j.cpc.2014.03.007.
  * (71) G. P. Salam, J. Rojo, A Higher Order Perturbative Parton Evolution Toolkit (HOPPET), Comput. Phys. Commun. 180 (2009) 120–156.  arXiv:0804.3755, doi:10.1016/j.cpc.2008.08.010.
  * (72) V. Bertone, A. Karlberg, Benchmark of deep-inelastic-scattering structure functions at $\mathcal{O}(\alpha_{s}^{3})$ (4 2024).  arXiv:2404.15711.
  * (73) A. Karlberg, disorder: Deep inelastic scattering at high orders (1 2024).  arXiv:2401.16964.
  * (74) C. Bierlich, et al., A comprehensive guide to the physics and usage of PYTHIA 8.3 (3 2022).  arXiv:2203.11601, doi:10.21468/SciPostPhysCodeb.8.
  * (75) P. Skands, S. Carrazza, J. Rojo, Tuning PYTHIA 8.1: the Monash 2013 Tune, Eur. Phys. J. C 74 (8) (2014) 3024.  arXiv:1404.5630, doi:10.1140/epjc/s10052-014-3024-y.
  * (76) B. Cabouat, T. Sjöstrand, Some Dipole Shower Studies, Eur. Phys. J. C 78 (3) (2018) 226.  arXiv:1710.00391, doi:10.1140/epjc/s10052-018-5645-z.
  * (77) A. Bodek, U. K. Yang, Modeling deep inelastic cross-sections in the few GeV region, Nucl. Phys. B Proc. Suppl. 112 (2002) 70–76.  arXiv:hep-ex/0203009, doi:10.1016/S0920-5632(02)01755-3.
  * (78) A. Bodek, U. K. Yang, Modeling neutrino and electron scattering inelastic cross- sections in the few GeV region with effective LO PDFs TV Leading Order, in: 2nd International Workshop on Neutrino-Nucleus Interactions in the Few GeV Region, 2003.  arXiv:hep-ex/0308007.
  * (79) A. Bodek, I. Park, U.-k. Yang, Improved low Q**2 model for neutrino and electron nucleon cross sections in few GeV region, Nucl. Phys. B Proc. Suppl. 139 (2005) 113–118.  arXiv:hep-ph/0411202, doi:10.1016/j.nuclphysbps.2004.11.208.
  * (80) A. Bodek, U.-k. Yang, Axial and Vector Structure Functions for Electron- and Neutrino- Nucleon Scattering Cross Sections at all $Q^{2}$ using Effective Leading order Parton Distribution Functions (11 2010).  arXiv:1011.6592.
  * (81) A. Bodek, U. K. Yang, Y. Xu, Inelastic Axial and Vector Structure Functions for Lepton-Nucleon Scattering 2021 Update (8 2021).  arXiv:2108.09240.
  * (82) A. Candido, A. Garcia, G. Magni, T. Rabemananjara, J. Rojo, R. Stegeman, Neutrino Structure Functions from GeV to EeV Energies, JHEP 05 (2023) 149.  arXiv:2302.08527, doi:10.1007/JHEP05(2023)149.
  * (83) R. Plestid, B. Zhou, Final state radiation from high and ultrahigh energy neutrino interactions (3 2024).  arXiv:2403.07984.
  * (84) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, G. Watt, LHAPDF6: parton density access in the LHC precision era, Eur. Phys. J. C 75 (2015) 132.  arXiv:1412.7420, doi:10.1140/epjc/s10052-015-3318-8.