JUNO sensitivity to low energy atmospheric neutrino spectra

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 20 kton 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 $\nu_e$ and $\nu_\mu$ 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 $\nu_e$ and $\nu_\mu$ 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 100 MeV and 10 GeV, showing a great potential of the detector in the atmospheric low energy region.


Introduction
Atmospheric neutrinos are a naturally occurring neutrino source. They originate from the decays of π and K produced in extensive air showers initiated by the interaction of cosmic rays with the Earth's atmosphere [1][2][3][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.7 for E ≤ 10 6 GeV and γ 3.0 above that value [5]. At energies larger than 5 × 10 8 GeV, the spectrum becomes steeper (γ 3.2) [6] and it flattens again (γ 2.7) when E ≥ 3 × 10 9 GeV. In the interaction of a single high energy Cosmic Ray with the nuclei of the Earth's atmosphere, hundreds or thousands mesons can be produced. The atmospheric 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][8][9][10][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 ∼1 GeV 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 K becomes 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. Figure 1: Present measurements of the atmospheric neutrino energy spectrum, compared with theoretical predictions. Data and models are reported separately for ν e and ν µ . Figure from [7].
Given the very short lifetime of these particles, the associated neutrino flux is commonly referred as "prompt" [12][13][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 from the production point to the detector varies from O (10) to O(10 4 ) 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 the 1960s [16,17]. Further measurements led to the discovery 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 detection granularity. The efforts are mostly concentrated on flavor 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 Figure 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][9][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 ν e and ν µ fluxes in the energy range 100 MeV -10 GeV.

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 Figure 2. The central detector (CD) consists of 20 kton of LS, contained in a 12 cm thick, highly transparent, acrylic sphere with a diameter of 35.4 m. The light produced in the LS is read out by 17612 20" high quantum efficiency (QE) photomultiplier tubes (PMTs) and 25600 3" PMTs, providing a total photo-coverage of more than 75%. About 13000 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 30 kHz, and a time resolution on single photo-electrons in the range from 1 ns to 10 ns. 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" PMTs will operate in digital mode, 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 40 m and constitutes the mechanical support for both the acrylic sphere and the PMTs. The central detector is submerged in a ∼44 m deep water pool (WP) filled with ∼30 kton 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.6 cm 2 and 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,5-Diphenyloxazole (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 ∼ 10 4 photons per MeV of deposited energy and an attenuation length greater than 20 m for 430 nm photons. The designed radio-purity levels of the JUNO LS are O(10 −16 ) g/g for the bulk 238 U, 232 Th, and 40 K 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 circumference of the sphere, will provide information regarding nonuniformity 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.

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], which we refer to as HKKM14 hereafter. The model assumes a Cosmic Ray spectrum based on BESS [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 100 MeV up to 10 TeV. 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 Figure 3. Hereafter, ν e and ν µ will be used to label both neutrinos and antineutrinos of electron and muon flavor, respectively. 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 has been simulated by means of the GENIE Neutrino Monte Carlo Generator [40,41] inside an energy range up to 20 GeV. The elemental composition of the neutrino target has been set as the one of the JUNO LS (mainly 12 C and 1 H, 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 the analysis results by independent evaluations and have been 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, several physical processes are included: electromagnetic interaction, 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 · 10 5 ν µ + ν 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.

Analysis strategy
As a large LS detector, JUNO achieves its best performance on events which are Fully-Contained within the volume, where a calorimetric measurement can be performed. Partially-Contained 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 ∼10 GeV 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 partiallycontained 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.

Fiducial cuts
Before applying the analysis selection to isolate ν e and ν µ populations, some preliminary cuts are applied to the large neutrino sample, with the aim of removing low-quality events. A first cut on the interaction vertex position is applied, 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 σ = 1 m has been applied to the MC interaction point (hereafter called vertex), in order to reproduce the uncertainty on the reconstructed position. We require that R vertex (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 1 m; 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 vertex position does not affect the performance of the selection procedure.
As described in Section 2, the CD is surrounded by a water Cherenkov detector acting as veto for atmospheric muons. Both muons and secondaries coming from partially-contained 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 (N W P 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 200 ns 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 partially-contained events are composed at 96% of ν µ CC interactions. The total efficiency for all ν e events is 68%, for ν µ events is 63%.

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 700 m 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 -4 Hz, corresponding to roughly 10 5 times the atmospheric neutrino event rate, considering that the average energy of atmospheric muons reaching JUNO is 207 GeV. 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 10 5 muon 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 [43] and simulating muons propagation through matter [44]. 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 Section 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 10 5 NPE, which is the region of interest for the analysis. An acceptance of < 2.3 · 10 −5 at 90% confidence level is achieved. The accuracy in the estimation of the acceptance will be improved by increasing the Monte Carlo statistics.

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 meters. On the contrary, muons with energy greater than 1 GeV 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 ν e CC 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 [45]. A high-precision measurement of the photon arrival time is 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.6 ns (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 µs time window have been considered. A time residual t res is then defined for each hit on the i-th 3" PMT as: where t i hit is the hit time on the i-th 3" PMT, n is the refraction index of the JUNO liquid scintillator and R i V is the distance between the reconstructed vertex position and the i-th 3" PMT. The time profile of the scintillation light emitted by ν e and ν µ CC events is different and the latter has a more prominent tail; therefore, the RMS of the t res distribution -hereafter called σ(t res ) -over the fired 3" PMTs can be used as discrimination parameter. In Figure 4, the σ(t res ) distribution is reported for the three populations: ν µ CC, ν e CC, 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 Figure 4 show a good separation between the ν e CC 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 ν e CC region. The reason is that a large fraction of the hadronic component of the secondaries is made of pions, that either decay to ν µ + µ (π + /π − ) or to two gammas (π 0 ). The first category is almost indistinguishable from ν µ CC events, while the second one results in electromagnetic showers and resembles the ν e CC 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 σ(t res ) distribution for NC events in Figure 4, since each of the four bins of NPE corresponds to a 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 ν e events, a value of σ(t res ) < 75 ns is required. The cut results in an efficiency for ν e events 42%, with respect to the large sample after fiducial cuts, and a residual contamination from ν µ less than 6%. A requirement of σ(t res ) > 95 ns is required to isolate ν µ events. In order to reduce the contribution from NC events at low energy, an additional requirement of NPE ≥ 5 × 10 5 has been set for ν µ selection, thus limiting the analysis to events with a neutrino energy 400 MeV. The efficiency for ν µ events is 85% with respect to the large sample after fiducial cuts and the residual ν e contamination is less than 20%. The residual NC events are populated both by ν e and ν µ .
The flavor identification procedure has been extensively checked by means of an independent analysis. A variation of the σ(t res ) 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 the atmospheric neutrino flux, we used the small Monte Carlo sample corresponding to ∼5 years of data-taking described in Section 3. The energy range considered for the atmospheric ν e flux is [−1.00, 1.05], expressed in log 10 (E ν /GeV) units and is divided in seven bins. The corresponding log 10 (NPE) range is [5.0, 7.2], divided in  seven bins as well. Similarly, the energy range for ν µ is [−0.30, 1.05], divided in seven bins, and the corresponding log 10 (NPE) range is [5.7, 7.2], divided in eight bins. The distribution of log 10 (NPE) is reported in Figure 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. Table 1: Summary of selections flow for ν e and ν µ fluxes, in terms of number of events in the analysis region, applied to the small data sample corresponding to ∼5 years of data taking. The values are reported before the selections, after the fiducial cuts and after the σ(t res ) selection. The residual background is also reported.

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 spectrum is 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 N can be expressed in terms of the primary neutrino spectrum E as where A ji is the likelihood matrix, which can be estimated by means of a full detector simulation. The relationship in Eq. 2 can be inverted by using the unfolding matrix U ij : The unfolding matrix U ij can be evaluated by means of an iterative Bayesian procedure [46]. In this case, the likelihood matrix A ji can be expressed as the probability P (N j |E i ) of detecting an event in the j-th bin of the NPE spectrum N j produced by the interaction of a neutrino in the i-th bin of the energy spectrum E i . The values of P (N j |E i ) are evaluated by means of a Monte Carlo detector simulation and normalized as j A ji = 1−ε i , where ε i takes into account the inefficiency in measuring the energy E i . The wrong-flavor events are also included in A ji . Using Bayes' theorem, the unfolding matrix U ij can be written as: The prior P 0 (E i ) is the probability for a single event to fall into the i-th energy bin. Once the unfolding matrix is known, a first estimation of the spectrum can be produced: The normalized values ofÊ i are used iteratively as the new set of probabilities P (E i ), in order to obtain an updated value of P (E i |N j ) and therefore ofÊ i . 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 P 1 (E i ). As prior distribution, the HKKM14 model has been used. Further details are given in Section 4.5. Figure 6 shows the likelihood matrix for both ν e and ν µ events, evaluated according to the binning described in Section 4.3 and including the contribution of the background.

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 Figure 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 the bins with highest statistics up to ∼15% in the highest-energy bins.
Selection criteria The selection procedure is in general not intended to produce any bias on the final sample. As explained in Section 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 σ(t res ) could result in a substantially different value of the unfolded flux, due to the wide stochastic fluctuations. The whole analysis has been therefore performed by varying the nominal cut value of σ(t res ) by 1 ns steps in a [ -5 ns, +5 ns ] time window. The differences in the unfolded flux are relevant in the bins with less statistics, for the reasons explained above. The total contribution to the bin uncertainty is evaluated as the standard deviation of the flux values distribution in each bin.
Flavor oscillation The current uncertainties on the global fit oscillation parameters are reported in Table 2, which are assumed to be Gaussian. A toy MC has been used to generate 1000 data sets, randomly varying the oscillation parameters within the experimental uncertainties, including the mass ordering and assuming no correlation. The final distribution in the unfolded flux is fit in each bin with a Gaussian function. Since the resulting dispersion is rather small in every bin, the total per-bin uncertainty contribution is quoted as the displacement of the distribution fit peak with respect to the nominal flux value. The total contribution from oscillation parameters uncertainty is estimated to be below 1% on the entire spectrum. The only exception is the first bin of the ν µ spectrum, where oscillation effects are not negligible, which results into an uncertainty corresponding to a σ of 1.2%.

Cross-section
The uncertainties on neutrino cross-section impact directly on the number of observed events. In the MC simulation process, as described in Section 3, neutrino interactions are managed by the GENIE software. The full list of uncertainty sources considered by GENIE is provided in [41]. A comprehensive handling of the whole list is not trivial, since it requires the simultaneous calculation of modified interaction probabilities in a wide parameter space. In this study, the evaluation of the cross-section uncertainty is based on experimental measurements provided by the T2K Collaboration [47][48][49], extrapolated from the associated data releases. Assuming the uncertainty on the measured cross section values to be Gaussian, the related visible spectrum is modified accordingly, within 1σ interval. The propagated uncertainty on the unfolded flux is evaluated by unfolding 1000 toy MC data sets, with NPE bin contents altered according to random variations of the cross section parameters. The unfolded spectra distributions are then fit in each bin with a Gaussian function, whose σ is quoted as the related uncertainty contribution. The uncertainty in the neutrino cross-section values has a large impact in the final reconstructed flux, up to 20%.
Unfolding procedure Although the iterative Bayesian unfolding method is data-driven, the particular MC sample may have an influence on the final result. This means that the initial estimation of the likelihood matrix may have an intrinsic bias, as well as the choice of the prior. The relative impact should be small, but it can have an impact in the unfolding bins with low statistics. The net effect cannot be exactly computed, but a reliable estimation can be achieved by unfolding modified data sets, generated by assuming a primary MC distribution reasonably far from the one used to evaluate the probabilities. The modified spectra are produced from the original MC by means of a re-weighting procedure. The new spectrum can be expressed in the i-th unfolding bin as: where α acts on the absolute normalization and γ on the shape of the primary spectrum. These two parameters are considered to range in the following intervals: ±0.05 for α and ±0.2 for γ.
The size of variation corresponds approximately to a 1σ uncertainty interval in the predicted spectra. In Figure 7 the comparison between each toy data sample and the corresponding unfolding result is reported, together with the fractional deviation between the input and the unfolded result. The deviation is below 1% and turns out to be slightly higher in the case of maximum variation of α and γ in the bins with lower statistics. The conditional probabilities used in the unfolding procedure have been carefully evaluated by using different methods for both ν e and ν µ samples. The relative deviations obtained for different energy bins and different NPE bins are reported in Table  3, as an example, for the ν e sample. The effect on the obtained spectra turns out to be negligible.
The contributions of each uncertainty source are reported in Figure 8, for each unfolding bin. The total uncertainty reported is calculated as the sum in quadrature of all contributions. The neutrino cross-section uncertainty represents the dominant contribution over the whole unfolded spectrum. The statistical uncertainty has also an important weight in high-energy bins. The total  Figure 7: Top panels: modified unfolded spectra (markers), together with the corresponding input (dashed lines). Bottom panels: relative deviation. Left: ν e spectra; right: ν µ spectra. Φ M C represents the nominal flux model [33], which is reported as the black dashed line. Four sets of modified spectra are plotted, whose α and γ values are reported in the figures. flux uncertainty ranges from a minimum value of 10-15% in the O(1 GeV) energy region, up to a 20-25% in the edge bins.

Results and discussion
The unfolded ν e and ν µ energy spectra are shown In Figure 9. The binning is described in Section 4.4. The predicted HKKM14 flux [33] is also reported, both at the source and including oscillation effects along the baseline. The oscillation-induced flux deficit in the ν µ flux below 10 GeV is clearly visible.
JUNO is able to reconstruct the energy spectrum of atmospheric neutrinos in the energy range [100 MeV -10 GeV], usually referred to as the "low-energy" region. This work, although based on simulated data only, shows the good capabilities of a large LS based detector like JUNO to measure the atmospheric neutrino flux. The energy region considered is already populated by other measurements, however some discrepancies still remain. JUNO can provide additional information about this interesting energy region, helping models in constraining their predictions. The quoted uncertainty is competitive with present experimental results and shows a margin of improvement by the increase of exposure time. Although JUNO's design is not optimized for atmospheric neutrino physics, the extremely good performances in the atmospheric neutrino energy reconstruction can be fully exploited for the measurement of the energy spectrum. Moreover, atmospheric neutrinos are a natural source which will be fully accessible from the beginning of data taking.

Conclusions
The JUNO detector has been designed from the beginning as a state-of-the-art detector for neutrino physics. The large dimensions of the detector, as well as its dense instrumentation, pave the way to an entire series of measurements, in a multi-purpose approach. The atmospheric neutrino flux is a natural source that can be observed, from the very beginning of data-taking. Although the detector design is not optimized for this class of events, the large active volume and the fine energy resolution allow to reconstruct the energy spectrum with a competitive precision, especially in the low-energy region.
In this work, a large set of MC events has been generated to evaluate the detector performances. A smaller set has been used to simulate the real data. Thanks to the timing performances of JUNO, the flavor of primary neutrinos can be separated with a limited residual contamination. A rejection power of the order of 10 5 has been applied to reduce atmospheric muon background. The atmospheric neutrino energy spectrum has been reconstructed in the energy range [100 MeV -10 GeV], separately for ν e and ν µ , assuming a ∼5 years detector livetime. The reconstructed spectra lie inside an interesting energy region, where previous measurements show some discrepancies. The results obtained show the good performance of JUNO in detecting the atmospheric neutrino flux in the low energy region, where theoretical models have large uncertainties. The inferred information can provide a fruitful input to constrain flux predictions, which are essential to evaluate the impact of atmospheric neutrinos in the search of rare events. The JUNO Collaboration