JUNO sensitivity to low energy atmospheric neutrino spectra

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ JUNO sensitivity to low energy atmospheric neutrino spectra © The Author(s) 2021 Published version The JUNO Collaboration The JUNO Collaboration. (2021). JUNO sensitivity to low energy atmospheric neutrino spectra. European Physical Journal C, 81(10), Article 887. https://doi.org/10.1140/epjc/s10052-02109565-z 2021 Eur. Phys. J. C (2021) 81:887 https://doi.org/10.1140/epjc/s10052-021-09565-z Regular Article - Experimental Physics JUNO sensitivity to low energy atmospheric neutrino spectra The JUNO Collaboration Received: 10 June 2021 / Accepted: 19 August 2021 © The Author(s) 2021 Abstract Atmospheric neutrinos are one of the most relevant natural neutrino sources that can be exploited to infer properties about cosmic rays and neutrino oscillations. The Jiangmen Underground Neutrino Observatory (JUNO) experiment, a 20kton liquid scintillator detector with excellent energy resolution is currently under construction in China. JUNO will be able to detect several atmospheric neutrinos per day given the large volume. A study on the JUNO detection and reconstruction capabilities of atmospheric νeand νμfluxes is presented in this paper. In this study, a sample of atmospheric neutrino Monte Carlo events has been generated, starting from theoretical models, and then processed by the detector simulation. The excellent timing resolution of the 3” PMT light detection system of JUNO detector and the much higher light yield for scintillation over Cherenkov allow to measure the time structure of the scintillation light with very high precision. Since νeand νμinteractions produce a slightly different light pattern, the different time evolution of light allows to discriminate the flavor of primary neutrinos. A probabilistic unfolding method has been used, in order to infer the primary neutrino energy spectrum from the detector experimental observables. The simulated spectrum has been reconstructed between 100MeV and 10GeV, showing a great potential of the detector in the atmospheric low energy region. 1 Introduction Atmospheric neutrinos are a naturally occurring neutrino source. They originate from the decays of πand Kproduced in extensive air showers initiated by the interaction of cosmic rays with the Earth’s atmosphere [1–4]. The energy spectrum of primary cosmic rays above 100 MeV can be described by a power law dN dE ∝E−γ, where the spectral index is γ≃2.7forE≤106GeV and γ≃3.0 above that value [5]. ae-mail: [email protected] At energies larger than 5 ×108GeV, the spectrum becomes steeper (γ≃3.2) [6] and it flattens again (γ ≃2.7)when E≥3×109GeV. In the interaction of a single high energy Cosmic Ray with the nuclei of the Earth’s atmosphere, hundredsorthousandsmesonscanbeproduced.Theatmospheric neutrino energy spectrum spans a wide range from the MeV up to the PeV scale and can be roughly described by a power law [7–11]. The spectral index is, in general, steeper than that of primary cosmic rays, since the parent mesons lose a large fraction of their energy before decaying. The spectrum intensity is suppressed at sub-GeV energies reflecting the rigidity cutoff, that describes the shielding provided by the geomagnetic field against the arrival of cosmic rays particles from outside the magnetosphere. Neutrinos originating from muon decays contribute mainly up to a few GeV. The flavor ratio (νμ+¯νμ)/(νe+¯νe) is around two at ∼1GeV and increases as the energy increases, since more muons are likely to reach the Earth’s surface without decaying. At energies above hundreds of GeV, the decay length of πand Kbecomes longer than their path length in the atmosphere, leading to a neutrino flux reduction. At the highest energies, the decay of heavy charmed mesons is expected to dominate the atmospheric neutrino production. Given the very short lifetime of these particles, the associated neutrino flux is commonly referred as “prompt” [12–14]. Since the Earth is mostly transparent to neutrinos below the PeV energy scale, an atmospheric neutrino detector is able to see neutrinos coming from all directions. The distance fromtheproductionpointtothedetectorvariesfrom O(10)to O(104)km, depending on the zenith angle [15]. The angular distribution has a characteristic shape with an increased flux towards the horizontal direction (with respect to the vertical direction), due to the longer path length of parent particles in the atmosphere. In the sub-GeV energy region there is an asymmetry along the East–West axis, which reflects the azimuthal dependence of the rigidity cutoff of the cosmic rays. Atmospheric neutrinos were detected for the first time in the1960s[16,17]. Further measurements led to the discovery 0123456789().: V,-vol 123 887 Page 2 of 16 Eur. Phys. J. C (2021) 81:887 Fig. 1 Present measurements of the atmospheric neutrino energy spectrum, compared with theoretical predictions. Data and models are reported separately for νeand νμ. Figure from [7] of neutrino flavor oscillations in 1998 [18]. Some of the missing pieces in the puzzle of neutrino physics are going to be addressed also by means of atmospheric neutrinos. The field of research is currently very active and several experiments are scheduled in the coming years to answer the unsolved questions. Next-generation detectors for atmospheric neutrino physics plan to significantly improve performances, compared to present ones, by increasing their size and detectiongranularity.Theeffortsare mostlyconcentratedonflavor oscillation physics, pushing the detectors sensitivity for the neutrino mass ordering (MO) and the CP phase δin the neutrino sector. The most prominent examples are DUNE [19], Hyper-Kamiokande [20], INO [21], ORCA [22], and PINGU [23]. In Fig. 1, present measurements of the energy spectrum of atmospheric neutrinos are reported, including predictions from theoretical models. Measurements performed over the last decades, up to present times, are able to cover a very wide range in the neutrino energy, from several hundreds of MeV to several hundreds of TeV. This sector has been explored predominantly by Cherenkov detectors, such as Super-Kamiokande [7] and IceCube [8–10,24]. The Jiangmen Underground Neutrino Observatory (JUNO), currently under construction in China, will be able to detect several atmospheric neutrinos per day. JUNO is going to become the largest liquid scintillator (LS) based detector ever built, having a target LS mass more than one order of magnitude larger than present ones. The large detector mass is one of the key-points for atmospheric neutrino detection, since it is comparable to the largest present water-based detector, Super-Kamiokande. Despite the limited ability of JUNO in tracking single particles after a neutrino interaction, with respect to large Cherenkov detectors, and a slightly reduced accessible statistics, the LS nature of the detector allows more precise measurements towards the low energy region. This sector of the energy spectrum is still not fully covered by present and past experiments. Furthermore, it also corresponds to the region where theoretical models have the largest uncertainties. The atmospheric neutrino flux measurements by means of the JUNO detector allow to investigate the neutrino MO and the θ23 octant. It is possible to pursue also the CP phase δ measurement. In our work, we investigate JUNO’s potential for measuring the atmospheric νeand νμfluxes in the energy range 100MeV–10GeV. 2 JUNO experiment The JUNO experiment [25,26] is a LS neutrino detector currently under construction in a dedicated underground laboratory (about 700 m deep, 1800 m.w.e.) near Kaiping, Jiangmen city, Guandong province (P. R. China). A sketch of the detector is shown in Fig. 2. The central detector (CD) consists of 20kton of LS, contained in a 12 cm thick, highly transparent, acrylic sphere with a diameter of 35.4m. The light produced in the LS is read out by 17,612 20” high quantum efficiency (QE) photomultiplier tubes (PMTs) and 25,600 3” PMTs, providing a total photo-coverage of more than 75%. About 13,000 of the 20” PMTs are Microchannel Plate (MCP) PMTs, developed by the JUNO collaboration and currently being produced by the North Night Vision Technology company. The remaining 5000 20” PMTs consist of the R12860 model produced by Hamamatsu. Both of these PMTs have a photon detection efficiency greater than 27%. For the 20” PMTs, the full waveform will be acquired. Their large photon collection area, however, has the consequence of a large dark noise rate, on average of the order of 30kHz, and a time resolution on single photo-electrons in the range from 1 to 10ns. The additional 3” PMTs, built by the HZC company, are deployed in the 20” PMTs’ lattice structure, in order to reduce any possible systematics due to the loss of linearity in charge reconstruction and to improve the timing measurements [27]. Due to their small area, the 3” PMTswill operate in digitalmode,thus being an independent readout system that can be exploited for cross-calibrating the 20” PMTs energy response. This feature becomes extremely important for high-energy events, where millions of photons are produced. Furthermore, due to the size difference, the Transit Time Spread (TTS) of 3” PMTs is of the order of the nanosecond, while the 20” PMTs one is larger, on average. The acrylic sphere is surrounded by a stainless steel truss structure, which has a diameter of about 40m and constitutes the mechanical support for both the acrylic sphere and the 123 Eur. Phys. J. C (2021) 81:887 Page 3 of 16 887 Fig. 2 Layout of the JUNO detector PMTs. The central detector is submerged in a ∼44m deep water pool (WP) filled with ∼30kton of ultrapure water and instrumented with 2400 20” MCP-PMTs. It acts as an active Cherenkov muon veto and shields the CD against external radioactivity. The walls of the pool are covered with highreflectivity Tyvek film, in order to increase the photon collection and allowing to veto Cosmic Ray muons with >95% efficiency. A Top Tracker (TT) is placed on top of the water pool, to improve the total veto efficiency and the reconstruction of atmospheric muons. The TT consists of three layers of scintillator strip detectors, refurbished after the decommissioning of the OPERA experiment Target Tracker [28]. It has a granularity of 2.6×2.6cm 2and a coverage of about 60% of the WP top surface. The JUNO LS mixture consists of three components: linear alkylbenzene (LAB) as solvent, 2.5 g/l of 2,5diphenyloxazole (PPO) as scintillation fluor and 3 mg/l of 1,4-Bis(2-methylsyryl) benzene (bis-MSB) as wavelength shifter [29]. This mixture ensures an effective light yield of ∼104photons per MeV of deposited energy and an attenuation length greater than 20m for 430nm photons. The designed radio-purity levels of the JUNO LS are O(10−16) g/g for the bulk 238U, 232Th, and 40K contaminants [30]. The calibration of the JUNO CD will be performed by four different systems [31]. An automated calibration unit (ACU) will deploy different radioactive sources along the detector vertical axis. The ACU system is also designed to deploy a laser source, with a photon intensity that can cover a range from hundreds of keV up to O(TeV) equivalent energy. Two Cable Loop Systems (CLS) will instead place sources across two planes. A guide tube (GT) system, installed on the outer circumferenceofthesphere, willprovideinformationregarding non-uniformity at the CD boundary. A Remote Operated Vehicle (ROV) will finally deploy sources in the whole detector volume. Periodical calibration campaigns will ensure to keep the overall energy resolution around 3%/√E/MeV in the MeV energy region, where the analysis for the neutrino MO will focus. Atmospheric neutrinos interacting inside JUNO can produce different final states, depending on the nature of the interaction they undergo. A first distinction can be done between charged-current (CC) and neutral-current (NC) interactions. In the first case, the lepton of the same flavor of the interacting neutrino is produced and therefore the original neutrino information is preserved. In NC interactions, on the contrary, only a hadronic state is visible and the flavor of the interacting neutrino cannot be inferred. In the energy of interest for atmospheric neutrinos, the dominant interaction is the neutrino-nucleon scattering. The most prominent channels are the elastic and quasi-elastic scattering, the resonant production, and the deep inelastic scattering [32]. This last classification concerns only about the final hadronic products, which give information about the original energy of the neutrino, but are not sensitive to the interaction flavor. Instead, the CC or the NC nature of the neutrino interaction implies a fundamental difference in the visible products. Apart from the flavor information, the absence of the flavor-corresponding lepton in the final state for NC interactions means also that the neutrino carries away part of its initial energy, which is not released inside the detector. NC events are therefore expected to be concentrated at lower values of the visible energy, while CC ones dominate at higher energy. 3 Monte Carlo dataset The study of the atmospheric neutrino flux is usually based on the predictions of the expected flux made by Monte Carlo simulations. In this work, we consider the predictions from the latest version of the HKKM model [33],whichwerefer to as HKKM14 hereafter. The model assumes a Cosmic Ray spectrumbasedonBESS[34,35] and AMS-01 [36] measurements. The DPMJET-III [37] and the JAM [38] hadronic interaction models are used for the simulation of the interaction with the Earth’s atmosphere. The HKKM model provides the calculation of the expected atmospheric neutrino flux at different locations, taking into account the latitude and longitude of the detector. The energy range spans from 100MeV up to 10TeV. Solar modulation and the asymmetry in the azimuthal distribution are also considered. In the HKKM parametrization, the neutrino flux is calculated at the source and therefore no oscillation effects are included. The HKKM14 atmospheric neutrino flux prediction, calculated at the JUNO experimental site, is shown in Fig. 3. Hereafter, νeand νμwill be used to label both neutrinos and antineutrinos of electron and muon flavor, respectively. 123 887 Page 4 of 16 Eur. Phys. J. C (2021) 81:887 Fig. 3 Expected atmospheric (neutrino + antineutrino) flux at the JUNO site, for νe(red) and νμ(blue), according to the HKKM14 model [33]. The flux is reported with (full line) and without (dashed line) considering neutrino oscillations In order to get a realistic prediction, neutrino oscillations have been applied to the original flux, including matter effects. The impact of the oscillation effects has been evaluated according to the standard 3-neutrinos mixing scheme [39]. The interaction of atmospheric neutrinos with the JUNO detector hasbeen simulatedby meansof the GENIE Neutrino Monte Carlo Generator [40,41] inside an energy range up to 20GeV. The elemental composition of the neutrino target has been set as the one of the JUNO LS (mainly 12C and 1H, with relative composition of 0.88 and 0.12, respectively). The output of the simulation contains information about the type of interaction that neutrinos undergo, either a CC or a NC one, and the full list of secondary particles and their associated properties (Particle ID, momentum, direction, ...). The contribution of ντCC interactions has been found not to affect theanalysis results byindependentevaluations and havebeen therefore not considered in the present work. Secondary particles produced in the interaction between neutrinos and the JUNO target material have been propagated in the detector by using a GEANT4-based Monte Carlo simulation. The JUNO detector simulation code has been developed within the SNiPER framework [42]. In the detector simulation, severalphysicalprocessesareincluded:electromagneticinteraction, decay, hadronic elastic and inelastic interactions, scintillation (including re-emission), Cherenkov emission, and optical absorption. A detailed optical model, including the optical properties of all the detector materials, is also implemented. The output relevant for the analysis includes the timestamp, the number of photoelectrons, and the position of each PMT hit. A data sample of about 5 ×105νμ+νe events have been generated, hereafter called large data sample, in order to set up the procedure used to reconstruct the atmospheric neutrino spectrum and to understand the detector response over a large statistics of events. A sample of 6500 events have been injected in the simulation as a separate Monte Carlo data sample, corresponding to a detector live time of about 5 years. This smaller data sample is hereafter identified as small data sample. 4 Analysis strategy As a large LS detector, JUNO achieves its best performance oneventswhicharefully-containedwithinthevolume,where a calorimetric measurement can be performed. Partiallycontained events, having some secondaries escaping from the CD active volume, are reconstructed with a worse energy resolution. This analysis therefore targets fully-contained events, to be accounted for reconstruction. This sets an intrinsic upper limit of ∼10GeV on the νμflux, since the highenergy muons produced in a CC interaction always escape the CD volume. For νe(νμ), the “golden” events consist of νe(νμ) fully-contained and CC events and the components to reduce are partially-contained and NC events of all flavor neutrinos. Through-going muons, that may be produced in νμinteractions with materials surrounding the CD, are not considered in this study. 4.1 Fiducial cuts Before applying the analysis selection to isolate νeand νμpopulations, some preliminary cuts are applied to the large neutrino sample, with the aim of removing low-quality events.Afirstcut on theinteractionvertex position isapplied, in order to remove events which release their energy near the edge of the acrylic sphere. These events typically exhibit a loss of linearity between the deposited and the collected energy, because part of the energy is released in the acrylic and water and not in the LS and because the closest PMTs collect a great amount of light and can undergo saturation. A Gaussian smearing with σ= 1m has been applied to the MC interaction point (hereafter called vertex), in order to reproduce the uncertainty on the reconstructed position. We require that Rvertex (i.e. the distance between the vertex and the center of the detector) is less than 16 m to ensure a linear detector response. The precision in the reconstruction vertex at lower energy (in the MeV range) is in general much better than 1m; on the contrary, at the GeV energy scale, secondary particles can deposit their energy on a long track and the events can no longer be considered as point-like. It has been checked that even an error of few meters on the ver123 Eur. Phys. J. C (2021) 81:887 Page 5 of 16 887 tex position does not affect the performance of the selection procedure. As described in Sect. 2, the CD is surrounded by a water Cherenkov detector acting as veto for atmospheric muons. Both muons and secondaries coming from partiallycontained neutrino events can release a certain amount of energy in the WP and produce a large amount of Cherenkov photons. Therefore, in order to remove partially-contained neutrino events and suppress the atmospheric muon background, we require the total number of hits seen by the water pool veto PMTs (NWP hits )to be less than 50, including the contribution from PMT dark noise. This latter term can become important for WP PMTs, because the single-count rate due to WP PMTs dark noise is high (up to several tens of kHz) and the total number of prompt hits from muons Cherenkov light can be small. Hits on WP PMTs are considered in a 200ns time window, which is approximately the time needed by a muon to cross the entire detector. The dark hits contribution is simulated on a statistical basis, assuming a binomial distribution. After applying the fiducial cuts described previously, the simulated large neutrino sample is composed at 97% of fully-contained event. The remaining partiallycontained events are composed at 96% of νμCC interactions. The total efficiency for all νeevents is 68%, for νμevents is 63%. 4.2 Atmospheric muon background The atmospheric muon background consists of the secondary muon flux produced after the interaction of cosmic rays with the atmosphere, in the same way as for neutrinos. The JUNO detector location is about 700m underground, therefore part of the muon radiation is able to penetrate the rock overburden and release energy inside the detector. The energy released by atmospheric muons inside JUNO is comparable with that of particles coming from atmospheric neutrino interactions (hundreds of MeV – several GeV). Muons can mimic the topology of atmospheric neutrino events and can therefore be a source of background. Although the external water Cherenkov veto is designed to reject these events with high efficiency, the atmospheric muon event rate is several orders of magnitude higher than that from atmospheric neutrino interactions. From preliminary calculations, their event rate inside the JUNO CD is around 3–4Hz, corresponding to roughly 105times the atmospheric neutrino event rate, considering that the average energy of atmospheric muons reaching JUNO is 207GeV. The desired acceptance rate for the atmospheric muon background must be therefore at least of the order of 10−5. In order to get a comprehensive picture of the atmospheric muon flux within the framework of this study, a full MC simulation is necessary. Atmospheric muons produce several millions of photons in the JUNO LS and the full detector simulation requires high CPU power and storage. For this purpose, a sample of only 105muon events has been generated, according to the energy and angular distributions evaluated at the JUNO site. The expected muon flux in the detector is calculated within the JUNO Collaboration according to a parametrized model at Earth surface [4] and simulating muons propagation through matter [43]. A detector simulation has been performed. Atmospheric muons in JUNO appear as high-energy tracks which release a large amount of energy both in the WP and in the CD. The fiducial cuts described in Sect. 4.1 require instead a low collected light inside the WP. Hereafter, the readout charge of the event, in terms of the number of PEs collected by CD 20” PMTs, is called NPE. NPE represents the observable used to reconstruct the neutrino energy. The same fiducial cuts have been applied to the muon sample, with the additional request of more than 105NPE, which is the region of interest for the analysis. An acceptance of <2.3×10−5at 90% confidence level is achieved. The accuracy in the estimation of the acceptance will be improved by increasing the Monte Carlo statistics. 4.3 Neutrino flavor identification As mentioned above, νe(νμ) CC interactions are the preferred detection channels, since the corresponding charged leptons have very different behaviours. Electrons lose energy quickly via bremsstrahlung and ionization and even at GeV energies their track length is no more than 1–2 m. On the contrary, muons with energy greater than 1GeV have longer tracks inside the detector volume. Low-energy muons, moreover, may decay inside the scintillator volume and give a delayed energy release from the Michel electron. The above differences make νμCC events more extended in time and space, with respect to νeCC events. The latter component has indeed a much shorter evolution. Hadronic particles are common to all classes of events and make up the visible part of NC events. Hadrons, in general, have a long energy release, because of their interactions and decays. The event time profile can be therefore exploited to discriminate between different classes of events [44]. A highprecisionmeasurement ofthephoton arrival timeis an important requirement. For this reason, the timing information is taken from the data of the 3” PMT system of the JUNO detector, which have a low TTS value. A Gaussian smearing with a typical width σ=1.6ns (taken from preliminary measurements) has been applied to the true Monte Carlo hit time over each 3” PMT. In order to be aligned to a realistic DAQ window, only events inside a 1.2μstime window have been considered. A time residual tres is then defined for each hit onthei-th3”PMTas: 123 887 Page 6 of 16 Eur. Phys. J. C (2021) 81:887 Fig. 4 Distribution of σ(tres)for νμCC (blue), νeCC (red) and NC events (green), for different ranges of NPE ti res =ti hit −n·Ri V c,(1) whereti hit is the hit time on the i-th 3” PMT, nis the refraction index of the JUNO liquid scintillator and Ri Vis the distance between the reconstructed vertex position and the i-th 3” PMT. The time profile of the scintillation light emitted by νeand νμCC events is different and the latter has a more prominent tail; therefore, the RMS of the tres distribution – hereafter called σ(tres)- over the fired 3” PMTs can be used as discrimination parameter. In Fig. 4,theσ(tres)distribution is reported for the three populations: νμCC, νeCC, and NC events. The variable is also reported separately in 4 different intervals of NPE, selected such as to have equal statistics in each of them. The plots in Fig. 4show a good separation between the νeCC and the νμCC component, over the whole energy range. The NC component appears to be overlapped mainly to νμCC events, with a tail also in the νeCC region. The reason is that a large fraction of the hadronic component of the secondaries is made of pions, that either decay to νμ+ μ(π+/π−)ortotwogammas(π0). The first category is almost indistinguishable from νμCC events, while the second one results in electromagnetic showers and resembles the νeCC component. The relative weight of charged and neutral pions in the final state changes across the energy, as well as that of nucleons. This feature is at the origin of the different shape of the σ(tres)distribution for NC events in Fig.4, since each ofthefour bins of NPE correspondstoa different energy interval. Protons and neutrons, moreover, have a time profile similar to the one of muons, because they result in a long-lasting energy release inside the LS. The contribution of NC component, however, becomes less significant at high energy, due to its steeper spectral shape. In NC events, indeed, part of the initial energy of the interacting neutrino is carried away by the neutrino itself and is not deposited inside the detector. Given the different features, two separate selection criteria are used to maximize CC events. In order to separate νeevents, a value of σ(tres)<75ns is required. The cut results in an efficiency for νeevents ≃42%, with respect to the large sample after fiducial cuts, and a residual contamination from νμless than 6%. A requirement of σ(tres)>95ns is required to isolate νμevents. In order to reduce the contribution from NC events at low energy, an 123 Eur. Phys. J. C (2021) 81:887 Page 7 of 16 887 Fig. 5 Distribution of the observable log10(NPE) in the analysis bins, after the νeselection (top) and the νμselection (bottom). The black dots represent the number of selected events in every bin j, with associated statistical error. The filled histograms reproduce the bin composition, in terms of the correct flavor (light blue) and wrong flavor (green) additional requirement of NPE ≥5×105has been set for νμselection, thus limiting the analysis to events with a neutrino energy 400MeV. The efficiency for νμevents is 85% with respect to the large sample after fiducial cuts and the residual νecontamination is less than 20%. The residual NC events are populated both by νeand νμ. The flavor identification procedure has been extensively checked by means of an independent analysis. A variation of the σ(tres)cut has been applied and the resulting efficiency and contamination are in a very good agreement within the statistical and systematic errors reported in this work. In order to test the JUNO performance in reconstructing theatmospheric neutrino flux,weused the small MonteCarlo sample corresponding to ∼5 years of data-taking described in Sect. 3. The energy range considered for the atmospheric νeflux is [−1.00,1.05], expressed in log10(Eν/GeV)units and is divided in seven bins. The corresponding log10(NPE) rangeis[5.0,7.2],dividedinsevenbinsaswell.Similarly,the energy range for νμis [−0.30,1.05], divided in seven bins, andthecorrespondinglog10(NPE)rangeis[5.7,7.2],divided Table 1 Summary of selections flow for νeand νμfluxes, in terms of number of events in the analysis region, applied to the small data sample correspondingto∼5 years ofdatataking. The values are reportedbefore the selections, after the fiducial cuts and after the σ(tres)selection. The residual background is also reported νeνμ Events injected in the simulation 6500 Charge region 1725 1241 Fiducial cuts 1167 773 σ(tres)cut 495 661 Residual background 30 163 in eight bins. The distribution of log10(NPE) is reported in Fig. 5, in the bins used in the analysis. A summary of the small sample population is in Table 1, as a function of the flavor and of the cuts applied, in the NPE regions considered. 4.4 Unfolding The determination of the atmospheric neutrino energy spectrum, starting from the detector experimental observables, is a classical unfolding problem. In this case, the true spectrumis deconvolved from the distribution of the experimental observables, knowing the detector response. In the classical fitting method, on the other hand, the true distribution is extracted from the observables by directly comparing the experimental distribution with the results of a model prediction. The main benefit of the unfolding is that it does not require a particular choice of the spectrum parametrization. In a liquid scintillator detector like JUNO, the main observable for the energy reconstruction is the total number of photoelectrons NPE detected by the 20” PMTs. This value is related to the total energy deposit in the LS and therefore to the neutrino energy. The neutrino energy spectrum Eνis then unfolded from the NPE spectrum. In general, the observable NPE spectrum Ncan be expressed in terms of the primary neutrino spectrum Eas Nj= i Aji Ei,(2) where Aji is the likelihood matrix, which can be estimated by means of a full detector simulation. The relationship in Eq. 2can be inverted by using the unfolding matrix Uij: Ei= j UijNj.(3) The unfolding matrix Uij can be evaluated by means of an iterative Bayesian procedure [45]. In this case, the likelihood matrix Aji can be expressed as the probability P(Nj|Ei)of detecting an event in the j-th bin of the NPE spectrum Nj produced by the interaction of a neutrino in the i-th bin of the 123 887 Page 8 of 16 Eur. Phys. J. C (2021) 81:887 energy spectrum Ei. The values of P(Nj|Ei)are evaluated by means of a Monte Carlo detector simulation and normalized as jAji =1−εi, where εitakes into account the inefficiency in measuring the energy Ei. The wrong-flavor events are also included in Aji. Using Bayes’ theorem, the unfolding matrix Uij can be written as: Uij =P(Ei|Nj)=P(Nj|Ei)P0(Ei) iP(Nj|Ei)P0(Ei).(4) The prior P0(Ei)is the probability for a single event to fall into the i–th energy bin. Once the unfolding matrix is known, afirst estimation of the spectrum can be produced: ˆ Ei= j P(Ei|Nj)Nj.(5) The normalized values of ˆ Eiare used iteratively as the new set of probabilities P(Ei), in order to obtain an updated value of P(Ei|Nj)and therefore of ˆ Ei. The particular choice of the prior and the number of iterations may cause a small bias on the shape of the unfolded spectrum. A small number of iteration may not reflect the information given by the data, while a high number of iteration may amplify statistical fluctuations and distort the spectrum. Since the Bayesian method is strongly data driven, the effect of the particular choice of the prior is in general small, but is still taken into account as a source of systematic uncertainty. The prior should reflect, in principle, the best knowledge of the primary spectrum. The minimum bias is then achieved by adopting the true MC distribution. The strong data–driven nature of the iterative Bayesian method ensures very good results after few iterations. In this work, two iterations have been performed. A soft smoothing is applied to the first value of the probability P1(Ei). As prior distribution, the HKKM14 model has been used. Further details are given in Sect. 4.5. Figure 6shows the likelihood matrix for both νeand νμevents, evaluated according to the binning described in Sect. 4.3 and including the contribution of the background. 4.5 Uncertainties The total uncertainty on the atmospheric neutrino spectrum reconstruction is evaluated in each energy bin, including both contributions from statistics and systematic effects. Statistics The statistical uncertainty is due to the stochastic fluctuations that occur in the data bins. The amount of this fluctuations is visible in Fig. 5, for each observable bin. In order to evaluate their impact in the final unfolded spectrum, 1000 toy data sets have been generated, each time varying the bin content according to a Poisson distribution. The final distribution in each bin of the unfolded spectrum is then fit with a Gaussian function, whose σis quoted as the statistical uncertainty. The statistical contribution ranges from 5% in Fig. 6 Likelihood matrix for νe(top) and νμevents (bottom) the bins with highest statistics up to ∼15% in the highestenergy bins. Selection criteria The selection procedure is in general not intended to produce any bias on the final sample. As explained in Sect. 4.1, fiducial cuts have been used in the unfolding procedure in order to improve the accuracy of the probability evaluations. The energy range of the final reconstructed spectrum is well contained inside the energy range of the Monte Carlo generated events, guaranteeing that the fiducial cuts do not introduce any bias. The neutrino flavor identification based on the time residual selection, on the other hand, could bring some uncertainty in the data bins where the statistics is low: an even small variation in the chosen cut value of σ(tres)could result in a substantially different value of the unfolded flux, due to the wide stochastic fluctuations. The whole analysis has been therefore per123 Eur. Phys. J. C (2021) 81:887 Page 15 of 16 887 Yang37,LeiYang18,XiaoyuYang10,YifanYang10,YifanYang2,Haifeng Yao10, ZafarYasin66,Jiaxuan Ye10,MeiYe10, Ziping Ye31, Ugur Yegin51, Frédéric Yermia47, Peihuai Yi10,NaYin 25, Xiangwei Yin10, Zhengyun You20, Boxiang Yu10, Chiye Yu18, Chunxu Yu33, Hongzhao Yu20,MiaoYu 34, Xianghui Yu33, Zeyuan Yu10, Zezhong Yu10, Chengzhuo Yuan10, Ying Yuan12, Zhenxiong Yuan13,ZiyiYuan 34, Baobiao Yue20, Noman Zafar66, Andre Zambanini51, Vitalii Zavadskyi67, Shan Zeng10, Tingxuan Zeng10, Yuda Zeng20, Liang Zhan10, Aiqiang Zhang13, Feiyang Zhang30, Guoqing Zhang10, Haiqiong Zhang10, Honghao Zhang20, Jiawen Zhang10, Jie Zhang10, Jingbo Zhang21, Jinnan Zhang10, Peng Zhang10, Qingmin Zhang35, Shiqi Zhang20,ShuZhang20,TaoZhang30,XiaomeiZhang10,XuantongZhang10,XueyaoZhang25,YanZhang10,YinhongZhang10, Yiyu Zhang10, Yongpeng Zhang10, Yuanyuan Zhang30, Yumei Zhang20, Zhenyu Zhang34, Zhijian Zhang18, Fengyi Zhao26, Jie Zhao10, Rong Zhao20, Shujun Zhao37, Tianchi Zhao10, Dongqin Zheng19, Hua Zheng18, Minshan Zheng9, Yangheng Zheng14, Weirong Zhong19, Jing Zhou9,LiZhou 10, Nan Zhou22, Shun Zhou10, Tong Zhou10, Xiang Zhou34, Jiang Zhu20, Kejun Zhu10, Zhihang Zhu10, Bo Zhuang10, Honglin Zhuang10, Liang Zong13, Jiaheng Zou10 1Yerevan Physics Institute, Yerevan, Armenia 2Université Libre de Bruxelles, Brussels, Belgium 3Universidade Estadual de Londrina, Londrina, Brazil 4Pontificia Universidade Catolica do Rio de Janeiro, Rio, Brazil 5Pontificia Universidad Católica de Chile, Santiago, Chile 6Universidad Tecnica Federico Santa Maria, Valparaiso, Chile 7Beijing Institute of Spacecraft Environment Engineering, Beijing, China 8Beijing Normal University, Beijing, China 9China Institute of Atomic Energy, Beijing, China 10 Institute of High Energy Physics, Beijing, China 11 North China Electric Power University, Beijing, China 12 School of Physics, Peking University, Beijing, China 13 Tsinghua University, Beijing, China 14 University of Chinese Academy of Sciences, Beijing, China 15 Jilin University, Changchun, China 16 College of Electronic Science and Engineering, National University of Defense Technology, Changsha, China 17 Chongqing University, Chongqing, China 18 Dongguan University of Technology, Dongguan, China 19 Jinan University, Guangzhou, China 20 Sun Yat-Sen University, Guangzhou, China 21 Harbin Institute of Technology, Harbin, China 22 University of Science and Technology of China, Hefei, China 23 The Radiochemistry and Nuclear Chemistry Group in University of South China, Hengyang, China 24 Wuyi University, Jiangmen, China 25 Shandong University, Jinan, China, and Key Laboratory of Particle Physics and Particle Irradiation of Ministry of Education, Shandong University, Qingdao, China 26 Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou, China 27 Nanjing University, Nanjing, China 28 Guangxi University, Nanning, China 29 East China University of Science and Technology, Shanghai, China 30 School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai, China 31 Tsung-Dao Lee Institute, Shanghai Jiao Tong University, Shanghai, China 32 Institute of Hydrogeology and Environmental Geology, Chinese Academy of Geological Sciences, Shijiazhuang, China 33 Nankai University, Tianjin, China 34 Wuhan University, Wuhan, China 35 Xi’an Jiaotong University, Xi’an, China 36 Xiamen University, Xiamen, China 37 School of Physics and Microelectronics, Zhengzhou University, Zhengzhou, China 38 Institute of Physics, National Yang Ming Chiao Tung University, Hsinchu 39 National United University, Miao-Li, Taiwan 40 Department of Physics, National Taiwan University, Taipei 123 887 Page 16 of 16 Eur. Phys. J. C (2021) 81:887 41 Charles University, Faculty of Mathematics and Physics, Prague, Czech Republic 42 University of Jyvaskyla, Department of Physics, Jyvaskyla, Finland 43 IJCLab, Université Paris-Saclay, CNRS/IN2P3, 91405 Orsay, France 44 Univ. Bordeaux, CNRS, CENBG, UMR 5797, F-33170 Gradignan, France 45 IPHC, Université de Strasbourg, CNRS/IN2P3, F-67037 Strasbourg, France 46 Centre de Physique des Particules de Marseille, Marseille, France 47 SUBATECH, Université de Nantes, IMT Atlantique, CNRS-IN2P3, Nantes, France 48 III. Physikalisches Institut B, RWTH Aachen University, Aachen, Germany 49 Institute of Experimental Physics, University of Hamburg, Hamburg, Germany 50 Forschungszentrum Jülich GmbH, Nuclear Physics Institute IKP-2, Jülich, Germany 51 Forschungszentrum Jülich GmbH, Central Institute of Engineering, Electronics and Analytics - Electronic Systems (ZEA-2), Jülich, Germany 52 Institute of Physics and Excellence Cluster PRISMA+, Johannes-Gutenberg Universität Mainz, Mainz, Germany 53 Technische Universität München, München, Germany 54 Eberhard Karls Universität Tübingen, Physikalisches Institut, Tübingen, Germany 55 INFN Catania and Dipartimento di Fisica e Astronomia dell Università di Catania, Catania, Italy 56 Department of Physics and Earth Science, University of Ferrara and INFN Sezione di Ferrara, Ferrara, Italy 57 INFN Sezione di Milano and Dipartimento di Fisica dell Università di Milano, Milano, Italy 58 INFN Milano Bicocca and University of Milano Bicocca, Milano, Italy 59 INFN Milano Bicocca and Politecnico of Milano, Milano, Italy 60 INFN Sezione di Padova, Padova, Italy 61 Dipartimento di Fisica e Astronomia dell’Università di Padova and INFN Sezione di Padova, Padova, Italy 62 INFN Sezione di Perugia and Dipartimento di Chimica, Biologia e Biotecnologie dell’Università di Perugia, Perugia, Italy 63 Laboratori Nazionali di Frascati dell’INFN, Roma, Italy 64 University of Roma Tre and INFN Sezione Roma Tre, Roma, Italy 65 Institute of Electronics and Computer Science, Riga, Latvia 66 Pakistan Institute of Nuclear Science and Technology, Islamabad, Pakistan 67 Joint Institute for Nuclear Research, Dubna, Russia 68 Institute for Nuclear Research of the Russian Academy of Sciences, Moscow, Russia 69 Lomonosov Moscow State University, Moscow, Russia 70 Comenius University Bratislava, Faculty of Mathematics, Physics and Informatics, Bratislava, Slovakia 71 Department of Physics, Faculty of Science, Chulalongkorn University, Bangkok, Thailand 72 National Astronomical Research Institute of Thailand, Chiang Mai, Thailand 73 Suranaree University of Technology, Nakhon Ratchasima, Thailand 74 Department of Physics and Astronomy, University of California, Irvine, California, USA 123