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

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}$.

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:

$$k_{\rm beam,1}= \left\{E_{1},0,0,+E_{1}\right\},$$

$$k_{\rm beam,2}= \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}$,

$$P_{1}= \left\{E,0,0,E\right\},$$

$$P_{2}= \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

$$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 H₂O molecules, the total cross section is given by

$$\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

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

$$\quad t_{\rm ISR} =\frac{\xi^{2}}{2(1-\xi y)}\bar{s}(1-y), \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

$$\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

$$m_{W}^{\rm OS}=80.3770\leavevmode\nobreak\ \mathrm{GeV}\,,\qquad\Gamma_{W}^{\rm OS}=2.085\leavevmode\nobreak\ \mathrm{GeV}\,,$$

$$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

$$V_{ud}=V_{cs}=0.97446\,,$$

$$V_{us}=V_{cd}=0.22456\,,$$

$$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.

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.

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. vtype 1 # photon exchange only

2. vtype 2 # Z echange only

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.