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 flavour 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 flaringblazar [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 modified 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 efficient 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 coefficient functions.The PDFs describing the quark and gluon densities in the proton are obtained from fits 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 fits 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 coefficient functions that enter into the structure function computation depend on the orderof the calculation [17]. At LO, most of the coefficients 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 configuration (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
Configuration 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 configurations 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 configuration used to compute tables for the structure functions.At low neutrino energies ( ∼
100 GeV), there is a significant 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 configurations 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 5flavours 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 configurations 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 configurations 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 final 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 flavour of the interacting quark q can be readily sampled from the relative contributions of the PDFs to the cross section, and theflavour 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 fixed for gluons. As a first pragmatic step, we sample theflavour of the outgoing quark using the LO structure functions. After x and Q are sampled fromthe NLO cross section, the outgoing quark flavour 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 benefit from NLO double-differential cross sections, while only using aleading-order approximation in the assignment of the flavour of the outgoing quark. The PDFsused contain heavy flavour quarks as part of the sea, so that also processes like b → t are included.For the LO ’flavour-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 flavour 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 flavour, 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 configuration 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 configurations (see Tab. 1) usingneutrino-Oxygen CC and NC interactions with a monochromatic neutrino flux at two differentenergies. Overall, the distribution are similar for the three configuration. Nevertheless, the amountof top quarks predicted by BGR differs from the other two configurations. 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 flux of E ν = GeV(top) and E ν = GeV (bottom). Each color represents a configuration 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 benefitfrom 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 benefit 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) [ ].].