# High-energy neutrino event simulation at NLO in Genie for KM3NeT and other observatories

HHigh-energy neutrino event simulation at NLO inGenie for KM3NeT and other observatories

The KM3NeT Collaboration ‡ ∗ ‡ [email protected] High-energy neutrino astronomy is a key pillar of multi-messenger astronomy and has the po-tential to also advance fundamental neutrino physics. An accurate simulation of the neutrinointeractions is key in any analysis of neutrino telescope data.The currently available generator codes of neutrino interactions use leading order expressions forthe differential cross sections which govern rate and kinematics of the events. These calculationsare affected by well-known large theoretical uncertainties, hence the need of beyond leading orderformalism. Also, a consistent use of modern Parton Density Functions, which are derived usingNLO (or NNLO) frameworks, is required.For this reason and others, the

GENIE event generator of neutrino interactions, which is a standardtool in the few-GeV energy regime, has so far been stated to be valid up to only about 1 TeV.The work reported here consists of a high-energy (up to 10 GeV) extension, which will beavailable in

GENIEv4 . It interfaces

PYTHIA6 for the hadronization and employs a pragmaticmethod to assign the ﬂavour of the struck and outgoing quarks, so that the heavy quark production,including bottom and top, is also correctly simulated.

Corresponding authors:

Alfonso Garcia † , Aart Heijboer . Nikhef Institute for Subatomic Physics, Amsterdam, The Netherlands36th International Cosmic Ray Conference -ICRC2019-July 24th - August 1st, 2019Madison, WI, U.S.A. ∗ for collaboration list see PoS(ICRC2019)1177 † Speaker. c (cid:13) Copyright owned by the author(s) under the terms of the Creative CommonsAttribution-NonCommercial-NoDerivatives 4.0 International License (CC BY-NC-ND 4.0). http://pos.sissa.it/ a r X i v : . [ h e p - e x ] A ug igh-energy neutrino event simulation at NLO in Genie Alfonso Garcia

1. Introduction

Neutrino telescopes have become a fundamental tool in astroparticle physics: the strong ev-idence of high energy cosmic neutrinos [1][2] and the association of a neutrino event to a ﬂaringblazar [3] being key examples. In addition, they can be used to study neutrino properties, such asoscillations [4] or cross sections [5]. A new generation of neutrino telescopes is foreseen for thenext years with KM3NeT [6] and IceCube-Gen2 [7], which will provide unprecedented sensitivityin many analyses. These developments should be accompanied by accurate tools for the simulationof neutrino interactions in or near the detector.Monte Carlo simulations of the neutrino interactions rely on computations of the Deep In-elastic Scattering (DIS) neutrino-nucleon cross section. At the moment, leading-order calculationsof the double-differential cross section underpin the main neutrino event generator codes (such asGENHEN [8] (based on LEPTO [9]), ANIS [10] or NuGen (a modiﬁed version of ANIS). In thelast years, several groups have provided calculations beyond leading order (i.e NLO or NNLO): theCSMS [11] and BRG [12].Such NLO cross-section calculations, however, do not constitute an event generator. A fullevent generator requires sampling of the full double-differential cross section in an efﬁcient way,and the treatment of the outgoing particles, which result from the breakup of the target nucleon.We present a neutrino generator that uses NLO cross sections, and can use the corresponding(also NLO) modern Parton Density Functions (PDFs). The scattering off heavy quarks and thehadronic shower development is simulated by

PYTHIA6 [13]. The resulting generator, named

HEDIS , is implemented as an extension to the

GENIE project [14].

GENIE has become a standardtool for long baseline experiments, which detect neutrinos with energies in the 1-10 GeV range. Inparticular, the (LO) DIS implementation in

GENIE is valid only for energies below about 1 TeV.With the inclusion of

HEDIS the use of [14] can be extended up to EeV energies.In the following sections, we outline the main features of neutrino scattering in the high en-ergy regime and we provide a description of the

HEDIS package. Additionally, we present threedifferent calculations of the neutrino cross sections: CSMS, BRG and the calculation used in theKM3NeT LoI [6].

2. Deep Inelastic Scattering

The neutrino-nucleon scattering can be modelled using the DIS formalism when the trans-ferred momentum is high enough such that the quark content of the target nucleon is resolved.The inclusive cross section of the Charged Current (CC) interaction can be described in terms ofthe Bjorken scaling variables x = Q / M N and y = ( E ν − E lep ) / E ν as (neglecting the mass of thelepton)d σ CC ν − N d x d y = G F M N E ν M W π ( Q + M W ) (cid:2) y x / F ( x , Q ) + ( − y ) F ( x , Q ) + y (cid:0) − y / (cid:1) F ( x , Q ) (cid:3) (2.1)A similar expression can be derived for the Neutral Current (NC) interactions. F , , are thestructure functions, which contain all the information about the structure of the nucleon or, in otherwords, the underlying QCD dynamics. 1 igh-energy neutrino event simulation at NLO in Genie Alfonso Garcia

The structure functions are the fundamental feature of the DIS model and they can be writtenas follows (assuming massless partons and the renormalisation and factorization scales both chosento be Q ) F i (cid:0) x , Q (cid:1) = ∑ j (cid:90) x d zz f j (cid:0) z , Q (cid:1) C i , j (cid:0) x / z , Q (cid:1) , (2.2)where j labels gluon and quarks, f j are the PDFs and C i , j are the coefﬁcient functions.The PDFs describing the quark and gluon densities in the proton are obtained from ﬁts tohadron collider data, applying the DGLAP formalism [15]. They form an input to the cross-sectioncalculation and, in GENIE, are obtained via the LHAPDF [16] package. In this program a wide listof different PDF ﬁts is available. They are stored in lookup tables as function of log x and log Q and some interpolation routines are implemented in order to obtain smooth results. The PDFs areprovided in a limited phase space so extrapolation procedures must be used for outer regions.The coefﬁcient functions that enter into the structure function computation depend on the orderof the calculation [17]. At LO, most of the coefﬁcients vanish or become delta functions, leading toa linear relation between the structure functions and the PDFs. This makes the structure functions(and thus the differential cross section) fast and easy to compute. Most existing neutrino generatorsuse this approach, relying at the same time on NLO PDFs, despite the inconsistency.In HEDIS , NLO structure functions are required, which are computed using numerical DGLAPevolution by an external package called

APFEL [18]. Another package, called

QCDNUM [19], wasalso tested leading, to similar results as shown in Fig 1. As the calculation of the structure func-tions is expensive, they are computed in

HEDIS prior to the event generation and stored in lookuptables as a function of log x and log Q for each interaction channel. For this, functionality from LHAPDF is used, which has similar tables for the PDFs.

Figure 1:

Ratio of the

QCDNUM17 and

APFEL proton structure functions for neutrino CC interactions usingthe HEDIS-CSMS conﬁguration (see Tab. 1). Small discrepancies are found when Q > m charm because theimplementation of the heavy-quark mass effects is slightly different in the two codes. As a result, the application that computes F i is the central feature of the HEDIS package. Itsdesign allows the manipulation of several features, such as order in pQCD with a desired massscheme for heavy quark production [20]. As explained in Sec. 1, three different calculations havebeen tested in this framework. In Tab. 1, their main features are listed.

The total neutrino cross section is important to determine the overall event rate of neutrinos2 igh-energy neutrino event simulation at NLO in Genie

Alfonso Garcia

Conﬁguration HEDIS-CSMS HEDIS-BGR HEDIS-LOIOrder pQCD NLO NLO LOMass scheme ZM-VFNS [21] FONLL [22] -PDF HERAPDF15 [23] NNPDF31sx+LHCb [24] CTEQ6D [25]x limits [ − , ] [ − , ] [ − , ] Q limits [ , ] [ . , ] [ . , ] Table 1:

Main features of the three different conﬁgurations used in

HEDIS . The limits are obtained fromthe

LHAPDF6 lookup tables. and can be used as a benchmark in the comparison with other published calculations to check thevalidity of

HEDIS .The total cross section is computed by integrating the differential expression over the allowedkinematic range in x and y . Currently, this integration is performed using a 2 dimensional grid withlog spacing. More complex algorithms are available, such as VEGAS [26], but the performancesare found to be very similar. Furthermore, the grid approach allows us to store the maximal dif-ferential cross section, which is relevant in the event generation process as explained in Sec. 2.3.Nevertheless, in the future these other integration methods will be available.The boundaries of the kinematical phase space can be derived using the invarinat mass of thehadronic system W and the four-momentum transfer Q . Assuming energy and momentum conser-vation, M N < W < √ s − M lep and E lep − | (cid:126) p lep | < ( Q + M lep ) / E ν < E lep + | (cid:126) p lep | . In practice, thephase space may shrink because the region of validity of the PDFs is smaller (see Sec. 2). Thus, theintegration limits depends on the conﬁguration used to compute tables for the structure functions.At low neutrino energies ( ∼

100 GeV), there is a signiﬁcant contribution to the cross sectionsfrom partons at Q ≡ − Q regime within GENIE should resolve this.At the highest neutrino energies, the PDFs are probed at very low x -values, which are oftenpoorly constrained by experimental data.The total cross section for nuclei is obtained from the appropriate sum of the neutron- andproton cross sections. F i are computed with PDFs from free nucleons, so nuclear-binding effects,such as shadowing [27], are not accounted for. None of the conﬁgurations presented in this workaccount for these nuclear effects. However, they could be introduced in future versions by follow-ing the procedure described in [12], where nuclear PDFs are computed using the EPPS16 globalanalysis [28].Fig. 2 shows the neutrino total cross section for an isoscalar target using the three differentcross-section calculations described in Sec. 1. These predictions have been compared with the In order to match the result from CSMS calculation, the mass of the top quark is not taken into account once 5ﬂavours are activated. The lower limit is found to be different in

LHAPDF5 lookup table ( x PDFmin = − ). The values of F i ’freeze’ below Q PDFmin = .

69 (i.e. F i ( x , Q ) = F i ( x , Q PDFmin ) if Q < Q PDFmin ). igh-energy neutrino event simulation at NLO in Genie Alfonso Garcia results from the

HEDIS package using the conﬁgurations described in Tab. 1. The result matchesat the 1% level in each case, demonstrating that the formalism adopted in

HEDIS is accurate. Theobserved differences between the three calculations at low energies is due to the different Q cut-off applied to the PDFs (see Tab. 1). Besides, BRG predicts a lower CC cross section because theproduction of heavy quarks is more suppressed in the FONLL mass scheme.

10 [GeV] ν E0.70.80.911.11.21.3 C S M S σ / σ CSMS (2001)HEDIS-CSMSBGR (2018)HEDIS-BGRKM3NeT LOI (2016) HEDIS-LOI

CC cross section µ ν

10 [GeV] ν E0.70.80.911.11.21.3 C S M S σ / σ NC cross section µ ν

10 [GeV] ν E0.70.80.911.11.21.3 C S M S σ / σ CC cross section µ ν

10 [GeV] ν E0.70.80.911.11.21.3 C S M S σ / σ NC cross section µ ν

10 [GeV] ν E0.70.80.911.11.21.3 C S M S σ / σ CSMS (2001)HEDIS-CSMSBGR (2018)HEDIS-BGRKM3NeT LOI (2016) HEDIS-LOI

CC cross section µ ν

10 [GeV] ν E0.70.80.911.11.21.3 C S M S σ / σ NC cross section µ ν

10 [GeV] ν E0.70.80.911.11.21.3 C S M S σ / σ CC cross section µ ν

10 [GeV] ν E0.70.80.911.11.21.3 C S M S σ / σ NC cross section µ ν Figure 2:

Ratio of the total cross section (per nucleon) with respect to CSMS for muon-neutrino scatteringwith an isoscalar target via CC (left) and NC (right). The results from three different calculations are shownwith dots. Lines show the results using

HEDIS for the conﬁgurations shown in Tab. 1. The predictions from

HEDIS are computed up to a certain E ν because above that value the bulk of the differential cross sectionlies bellow the x lower limit from the PDFs (obtained from LHAPDF6 ). In case of a charged-current interaction the outgoing lepton is simulated straightforwardlyfrom the kinematic variables sampled from the double differential cross section. The neutrinointeraction will in general result in the breakup of the target nucleon, which produces multipleoutgoing particles, which are more complex to simulate.The simulation of the hadronic ﬁnal state is, in GENIE, performed by

PYTHIA6 , which takesas input a description of the hadronic system just after the Z / W ± q → q (cid:48) exchange, in which aquark q from the nucleon has been transformed in a high-momentum quark q (cid:48) . In the input to the PYTHIA6 hadronization routine, the state is represented by three particles: one is the outgoingquark produced in by the interaction ( q (cid:48) ). The other two are baryons or mesons, which are com-posed by a combination of the quarks from the nucleon remnant and the (sea-quark) companion ofthe hit parton q , ¯ q .The scheme is straightforward to implement at LO, where the ﬂavour of the interacting quark q can be readily sampled from the relative contributions of the PDFs to the cross section, and theﬂavour of q (cid:48) follows from this (off-diagonal CKM matrix elements are included). However, beyondleading order, this is not a simple relation, as e.g. gluons will also contribute to the cross sectionat NLO, but the outgoing quark is not ﬁxed for gluons. As a ﬁrst pragmatic step, we sample theﬂavour of the outgoing quark using the LO structure functions. After x and Q are sampled fromthe NLO cross section, the outgoing quark ﬂavour f will be randomly picked with probability: P f = ∑ j V CKM f j σ LO f ( x , Q ) σ LOtot ( x , Q ) . (2.3)4 igh-energy neutrino event simulation at NLO in Genie Alfonso Garcia

In this manner, we beneﬁt from NLO double-differential cross sections, while only using aleading-order approximation in the assignment of the ﬂavour of the outgoing quark. The PDFsused contain heavy ﬂavour quarks as part of the sea, so that also processes like b → t are included.For the LO ’ﬂavour-picking’, we apply the ’slow rescaling’ [29] prescription to account for thekinematic suppression of the production of t , b and c quarks. In the NLO cross sections, sucheffects are accounted for by the mass-scheme.We note that the ﬂavour of the outgoing quark is typically of limited relevance in neutrinotelescopes as they are not sensitive to microscopic details of the hadronic shower. However, anotable exception is the leptonic decay of heavy ﬂavour, which may produce additional leptons. To simulate neutrino interactions on an event by event basis a sampling method must be usedto select the underlying kinematics of the event once the interaction process has been selected.This is usually done using the acceptance-rejection method. In

HEDIS , for a particular neutrinoenergy, a maximal differential cross section is computed and two independent variables ( x and y )are sampled using the double differential distribution of the chosen process.Because HEDIS is used in a wide energy range, the phase space of the variables is very broadas shown in Fig 3. Therefore, the sampling is performed over log x and log y taking into accountthe corresponding Jacobian. − − − − − − − log7 − − − − − − − y l og d x d y σ d × σ -p Charged Current µ ν [GeV] × =5 ν E − − − − − − − log7 − − − − − − − y l og d x d y σ d × σ -p Charged Current µ ν [GeV] × =5 ν E Figure 3:

Muon-neutrino-proton Charged Current differential cross section for E ν =

100 GeV (left) and E ν = GeV (right). These distributions were obtained using the conﬁguration HEDIS-BGR. The sharpcut-off in the diagonal is due to the Q lower limit of the structure functions (see Tab. 1). In Fig 4, we compare the x and y distribution for the three conﬁgurations (see Tab. 1) usingneutrino-Oxygen CC and NC interactions with a monochromatic neutrino ﬂux at two differentenergies. Overall, the distribution are similar for the three conﬁguration. Nevertheless, the amountof top quarks predicted by BGR differs from the other two conﬁgurations. The main reason is thatthe treatment of the mass for heavy-quark production is more accurate in BGR, leading to a largesuppression at these energies.

3. Conclusions and Outlook

The work presented here describes the main feature of the HEDIS package, which will be partof the next release

GENIEv4 . It will allow to use this Monte Carlo framework to study neutrino5 igh-energy neutrino event simulation at NLO in Genie

Alfonso Garcia − − − − − − log050001000015000200002500030000 e v en t s − − − − − − − − − − log0500010000150002000025000300003500040000 e v en t s Charged Current -O µ ν [GeV] = 10 ν E HEDIS-LOIHEDIS-LOI (charm)HEDIS-CSMSHEDIS-CSMS (charm)HEDIS-BGRHEDIS-BGR (charm) − − − − − − − − log0500010000150002000025000300003500040000 e v en t s Charged Current -O µ ν [GeV] = 10 ν E HEDIS-LOIHEDIS-LOI (top)HEDIS-CSMSHEDIS-CSMS (top)HEDIS-BGRHEDIS-BGR (top) − − − − − − log0500010000150002000025000 e v en t s Figure 4:

Kinematics of muon-netrino-Oxygen Charged Current interactions: log y (left) and log x (right). 10 neutrino interactions were simulated using a monochromatic neutrino ﬂux of E ν = GeV(top) and E ν = GeV (bottom). Each color represents a conﬁguration described in Tab. 1 (nuclear effectsare not accounted for). Dashed lines represent the contribution for a particular quark production. interactions at higher energies. The main feature is the new approach to calculate the structurefunctions, which allows running at higher orders in pQCD. Neutrino observatories will beneﬁtfrom this extension of the software. Also, it provides theoreticians a tool to test their calculationsin an standard way, which will be accessible to experimentalists.In addition, we plan to further develop

HEDIS in the future. In the short term, it will be impor-tant to include in the framework a method to quantify the error of these calculations. Currently, themain uncertainty arrises from the PDFs. Therefore, some extra developments are needed to incor-porate them in

HEDIS . Furthermore, an extension to

PYTHIA8 [30] is foreseen for next versionsof

GENIE , so this new tool will beneﬁt from that. This will allow us to study more sophisticatedapproaches for the hadronization step, such as parton showering.More accurate calculations of the Glashow resonance and order sub-leading processes havebeen developed recently and it would be interesting to include them in this framework [31]. Anotherinteresting upgrade is the extension of this model to the low energy regime where most of the longbaseline experiments work. In order to do so, calculations with state-of-the-art-PDFs must bestudied in the low Q regime. 6 igh-energy neutrino event simulation at NLO in Genie Alfonso Garcia

4. Acknowledgements

We are grateful to Juan Rojo, Rhorry Gauld and Valerio Bertone for their willingness to sharetheir expertise in pQCD. We would like also to thank the Genie Collaboration for supporting thisproject.

References [1] IceCube Collaboration,

Phys. Rev. Lett. (2014) 101101 [ ].[2] ANTARES Collaboration,

Astrophys. J. Lett (2018) 1 [ ].[3] IceCube Collaboration,

Science (2018) 147 [ ].[4] ANTARES Collaboration,

JHEP (2019) 113 [ ].[5] IceCube Collaboration, Nature (2017) 596 [ ].[6] KM3NeT Collaboration,

J. Phys.

G43 (2016) 084001 [ ].[7] IceCube-Gen2 Collaboration,

PoS

FRAPWS2016 (2017) 004 [ ].[8] KM3NeT/ANTARES Collaborations,

Workshop on MC simulation for neutrino telescopes (2018).[9] G. Ingelman et al.,

Comput. Phys. Commun. (1997) 108 [ hep-ph/9605286 ].[10] A. Gazizov et al.,

Comput.Phys.Commun. (2005) 203 [ astro-ph/0406439 ].[11] A. Cooper-Sarkar et al.,

JHEP (2011) 042 [ ].[12] V. Bertone et al., JHEP (2018) 217 [ ].[13] T. Sjostrand et al., JHEP (2006) 2006 [ hep-ph/0603175 ].[14] C. Andreopoulos et al., Nucl. Instrum. Meth.

A614 (2010) 1 [ astro-ph/0406439 ].[15] G. Altarelli et al.,

Nucl. Phys.

B126 (1977) 298.[16] A. Buckley et al.,

Eur. Phys. J.

C75 (2015) 132 [ ].[17] A. Vogt et al.,

Nucl. Phys. Proc. Suppl. (2006) 44 [ hep-ph/0608307 ].[18] V. Bertone et al.,

Comput. Phys. Commun. (2014) 1647 [ ].[19] M. Botje,

Comput. Phys. Commun. (2011) 490 [ ].[20] S. Forte et al.,

Nucl. Phys.

B834 (2010) 116 [ ].[21] J. Collins et al.,

Nucl. Phys.

B278 (1986) 934.[22] R. D. Ball et al.,

Phys. Lett.

B754 (2016) 49 [ ].[23] ZEUS, H1 Collaboration,

Proc. 40th Int. Symp. Mult. Dyn. (2010) [ ].[24] R. Gauld et al.,

Phys. Rev. Lett. (2017) 072001 [ ].[25] J. Pumplin et al.,

JHEP (2002) 2002 [ hep-ph/0201195 ].[26] T. Ohl, Comput. Phys. Commun. (1999) 13 [ hep-ph/9806432 ].[27] L. Frankfurt,

Nuclei, Phys. Rept. (2012) 255 [ ].[28] K. Eskola et al.,

Eur. Phys. J.

C77 (2017) 163 [ ].[29] R. M. Barnett,

Phys. Rev. Lett. (1976) 1163.[30] T. Sjostrand et al., Comput. Phys. Commun. (2008) 852 [ ].[31] R. Gauld,

Nikhef 2019-011 (2019) [ ].].