Measuring Oscillations with a Million Atmospheric Neutrinos

After two decades of measurements, neutrino physics is now advancing into the precision era. With the long-baseline experiments designed to tackle current open questions, a new query arises: can atmospheric neutrino experiments also play a role? To that end, we analyze the expected sensitivity of current and near-future water(ice)-Cherenkov atmospheric neutrino experiments in the context of standard three-flavor neutrino oscillations. In this first in-depth combined atmospheric neutrino analysis, we analyze the current shared systematic uncertainties arising from the common flux and neutrino-water interactions. We then implement the systematic uncertainties of each experiment in detail and develop the atmospheric neutrino simulations for Super-Kamiokande, with and without neutron-tagging capabilities, IceCube Upgrade, ORCA, and Hyper-Kamiokande detectors. We carefully review the synergies and features of these experiments to examine the potential of a joint analysis of these atmospheric neutrino data in resolving the θ 23 octant at 99% confidence level, and determining the neutrino mass ordering above 5 σ by 2030. Additionally, we assess the capability to constrain θ 13 and the CP -violating phase ( δ CP ) in the leptonic sector independently from reactor and accelerator neutrino data. A combination of the atmospheric neutrino measurements will enhance the sensitivity to a greater extent than the simple sum of individual experiment results reaching more than 3 σ for some values of δ CP . These results will provide vital information for next-generation accelerator neutrino oscillation experiments such as DUNE and Hyper-Kamiokande.


I. INTRODUCTION
Atmospheric neutrinos have played a key role in discovering and understanding neutrino oscillations [1,2].The first hint of resilient signatures of deviations from the Standard Model in neutrinos came from a deficit of muon neutrinos in IMB [3] and Kamiokande [4], which were later confirmed by Super-Kamiokande (SuperK) [5].These anomalies are now known to be due to neutrinos having nonzero, small masses, and the flavor states and mass states being misaligned [6].This misalignment produces a characteristic neutrino flavor change, often called oscillations due to its periodic behavior in vacuum, and this flavor change depends on the ratio of the distance from source to detector, L -called baseline -and the neutrino energy, E ν .
One of the reasons why atmospheric neutrinos are good discovery experiments is that they cover ten orders of magnitude in the ratio of baseline to energy, L/E ν .The baseline L covers ranges from 15 km to 12, 700 km and the neutrino energy ranges from O(10 −2 ) GeV to O(10 5 ) GeV.This results in a coverage of L/E ν from O(10 −4 ) km/GeV to O(10 6 ) km/GeV.This broad range of L/E ν and the fact that atmospheric neutrinos go through the largest amounts of matter [7] have allowed them to be used for measuring neutrino oscillation parameters [8,9] and placing stringent bounds on new neutrino states [10,11], non-standard neutrino interactions [12][13][14][15][16][17][18][19][20], space-time symmetries [21], and other physics beyond the Standard Model [22,23]; see Figure 1 for an artistic illustration of atmospheric neutrinos and the detectors used to observe them.
In the last two decades, the neutrino community has developed a vast program using accelerator, reactor, and solar neutrinos to measure their evolution, leading to the conclusion that atmospheric measurements will not contribute to precision neutrino physics.In this article, we change that paradigm, showing that atmospheric neutrinos will improve our knowledge of some of the largest unknowns describing the evolution of neutrino states.See Figure 2 for a comparison between our results, the present status, and the future predictions for the next generation of neutrino experiments.By combining the atmospheric measurements, we aim to achieve the most accurate determination of atmospheric parameters (∆m 2  31 and θ 23 ).Additionally, this combined approach will yield valuable insights into θ 13 due to the Earth matter effect, as we will elaborate upon later.Lastly, the atmospheric measurements are expected to offer a level of precision comparable to current results in determining δ cp .
As we will demonstrate in this article, data from current and soon-to-operate atmospheric neutrino experiments can be combined to address some of the most pressing questions in neutrino physics.The questions that we will discuss in this article can be organized into three categories: determining the neutrino oscillation parameters, establishing the neutrino mass spectra, and measuring the CP-phase in the lepton sector.The neutrino oscillation parameters are encoded in the so-called Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, which relates the neutrino weak and mass eigenstates [32].The precise determination of the lepton mixing parameters is crucial for understanding the neutrino evolution, and it can also be the first indication of a hidden flavor symmetry [33,34].Measuring a large CP -violation in the neutrino sector can also be an explanation to the baryon asymmetry of the early universe via a sphaleron process [35].Finally, the determination of the neutrino mass spectrum in the next few years will significantly impact experiments whose goal is to determine the absolute scale of the neutrino masses [36], to discriminate between the Dirac or Majorana nature of the neutrino masses [37,38], and even to understand the evolution of the universe [39].
In the rest of the section, we will provide a concise summary of the main results from this study.In atmospheric neutrino oscillations, the dominant mixing angle is θ 23 , and the relevant mass-squareddifference is ∆m 2  31 .The atmospheric mixing angle is currently the least precisely determined of the mixing angles and is known to be close to maximal, i.e., sin 2 2θ 23 ≈ 1.However, current experimental data points towards deviations from maximality, but is un-FIG.2. Comparison between the present and future projected sensitivities for the oscillation parameters.This figure showcases the power of atmospheric combined neutrino experiments, as they have smaller or comparable errors on these parameters than accelerator neutrino experiments also shown in this figure.Oscillation parameters that can be measured by atmospheric experiments are arranged in the vertical axis, while their measured precision is quantified in the horizontal axis.The present 1-sigma region allowed by fit to data from T2K [24] (green), NOvA [25] (light orange), SuperK atmospheric neutrinos [26] (light red), Deep Core atmospheric neutrinos [27] (turquoise), the reactor experiments (grey) and the projected sensitivity of Hyper-Kamiokande's accelerator program (blue) for 2030 [28] is compared with the expected 1-sigma region from a combined atmospheric neutrino analysis (this work in red-violet).The sensitivity from DUNE is excluded due to the low significance achievable by this experiment for 2030 assuming a 2029 starting date [29].
able to resolve the octant, i.e., whether θ 23 is smaller or greater than π/4.In terms of the neutrino flavor structure, the maximality and octant question can be rephrased as understanding the relative contribution of tau and muon flavor in the second mass state, where a maximal angle implies equal amounts, the first octant more muon, and the second octant less muon.Therefore, the large muon neutrino component of the atmospheric flux can provide significant knowledge on that parameter.In this article, we will demonstrate that by combining the SuperK, IceCube FIG. 3. Neutrino mass ordering sensitivity as a function of years in operation.The cyan (orange) band shows the sensitivity for rejecting the wrong ordering hypothesis for true normal (inverted) ordering, assuming fixed sin 2 θ13 = 0.022.The width of the bands covers the allowed values for sin 2 θ23 from 0.45 to 0.6.The black dot corresponds to the last reported SuperK neutrino mass ordering analysis [30].For comparison, we also include the prediction of the next-generation longbaseline neutrino experiments that are supposed to start in mid 2027 in the case of HyperK and 2029 in the case of DUNE [31].
Upgrade, ORCA and Hyper-Kamiokande (HyperK) measurements (see Section VI), we will reach a half percent-level precision on ∆m 2  31 and ∼ 2% in the case of sin 2 θ 23 : see Figure 2. Furthermore, the atmospheric measurements will allow us to discriminate the wrong octant at more than 3σ by 2030 assuming the current best-fit value of θ 23 : see Figure 18.As we will discuss in Section VIII, those measurements are primarily limited by statistics.In Appendix C we checked that our main results and conclusions are independent of the benchmark scenario considered.
The neutrino mass spectra enter the neutrino oscillation expressions through the sign of the masssquare differences.The neutrino mass states are conventionally labeled by their decreasing number of electron-neutrino components, where ν 1 is the state that has the most electron neutrino and ν 3 the least; in symbols, |U e1 | > |U e2 | > |U e3 |.From solar neutrino oscillations it is known that the second mass state, ν 2 is heavier than the first one, ν 1 ; however, its not known if ν 3 is heavier (normal ordering, NO) or lighter (abnormal or inverted ordering, IO) than the other two states.This is referred to as the neutrino ordering problem and is one of the main objectives of JUNO [40] and next-generation neutrino experiments DUNE [41] and Hyper-Kamiokande [28].Currently, the combination of neutrino data favors normal ordering mildly, although no conclusive measurement of this has been achieved to date.In this article, we show that the combination of neutrino experiments discussed in this work will determine the neutrino ordering at more than 5σ [42].Atmospheric neutrino experiments can reach 4σ sensitivity before the next generation of neutrino accelerator experiments starts taking data, Figure 3.The sensitivity to this parameter is mainly affected by statistics and the miss-reconstruction of ν τ interaction represented by the ν τ cross section as discussed in Section VIII.In Section III, we summarize the present status of the neutrino cross section measurements along with the most relevant experimental measurement that will happen in the next year and will contribute to increasing the sensitivity over that parameter.
Finally, we turn to the question of CP -violation in the leptonic sector, which has been previously discussed in the context of atmospheric neutrino experiments [43].The combination of the experiments considered in this work will provide a measurement of the CP -violating phase that will shrink the allowed region in a factor five (Figure 2), assuming the preferred value of T2K [44] and SuperK [30,45] measurements, being able to exclude more than half of the allowed parameter space at more than 90% confidence level for any value of δ CP : see Section VII.The capacity to measure the CP -phase is dominated by SuperK and HyperK due to their large low-energy neutrino efficiency and neutrino-antineutrino separation ability.The improved neutrino-antineutrino separation is due to a recent upgrade of SuperK, where the detector is doped with gadolinium (SKGd).See Section VI for a detailed description of the detector response and the simulation used in this analysis.The combination of a precision measurement of the oscillation parameters by the neutrino telescopes and an improve SuperK detector allows us to bring new, and complementary information on the CP -phase compared to that obtained by long-baseline experiments.This is of particular interest since measurements by current long-baseline experiments in the continental United States and Japan are in mild tension.This work shows, for the first time, that atmospheric experiments have the potential to weigh in on this tension.Finally, our work implies that before the operation of the next-generation neutrino detectors -DUNE, Hyper-Kamiokande, and IceCube-Gen2we will have two independent measurements of the CP -phase: one from the combination of accelerator neutrinos (i.e., T2K and NOvA), and another one from the combination of atmospheric neutrinos.
The combination of the SuperK, IceCube Upgrade, ORCA, and HyperK has a twofold purpose of solving open questions in neutrino physics and providing initial input for the next generation of neutrino experiments.In this work, we have developed for the first time the necessary tools to perform such a combined analysis along with the most realistic publicly available simulations for each experiment involved.All of this allows us to make the first in-depth analysis of those three experiments taking into account a detailed description and implementation of the detector responses, and their common systematic uncertainties.The work performed here is well-beyond what has currently been done in any prior global analysis of atmospheric neutrinos.The result is the realistic projection of the sensitivity of atmospheric neutrinos to the remaining mixing parameters and the identification of the uncertainties limiting these measurements, thus establishing a competitive and independent approach from accelerators to the neutrino physics precision era.
The rest of this article is organized as follows.In Section II, we outline the primary characteristics of the atmospheric neutrino flux and discuss associated uncertainties.Section III delves into the interaction of these neutrinos with water, while Section IV focuses on key aspects of reconstructing these interactions.Atmospheric neutrinos needs to cross the Earth before to reach the detector, the main aspect of the neutrino evolution are summarized in Section V.In this analysis, we have incorporated four experiments: Super-Kamiokande, IceCube Upgrade, ORCA and Hyper-Kamiokande.The scope of measurements conducted by each of these experiments are elucidated in Section VI.Section VII contains a description of the main results of this analysis, while Section VIII engages in a comprehensive discussion of their implications.We present out conclusion in Section IX.

II. ATMOSPHERIC NEUTRINO FLUXES
In this section, we will elucidate the key elements of the atmospheric neutrino flux that play a pivotal role in determining oscillation parameters.Additionally, we will examine the primary sources of uncertainty influencing the measurement of this flux and propose a parametrization that incorporates them.
Atmospheric neutrinos are produced by the collision of cosmic rays with Earth's atmosphere.The primary spectrum of cosmic rays spans from MeV to EeV energies [50,51], and is composed of free protons (∼ 80%) and bound nuclei (∼ 20%).Their interaction with nuclei in the atmosphere initiates hadronic showers on average about 20 km above the surface, producing copious amounts of mesons.Neutrinos are produced predominantly from the decay of muons, pions, and kaons, which dominate the muon-neutrino flux below 10, 100, and 10 6 GeV, respectively [52,53].With ν µ , ν e are also produced in the atmosphere, and above the TeV scale there also exists a ν τ flux from the charmed mesons decays [54,55].At the lowest energies, when decay of the mesons is prompt, the spectrum follows the cosmic-ray spectrum, but it softens by approximately one unit in spectral index as the mesons start interacting in the air [56].In Figure 4 (lower panel), we show the measurement of the total neutrino flux carried out by different experiments from ∼ 100 MeV to ∼ 100 GeV.We added the prediction from Honda et al. [49], which is used in this work as a benchmark scenario for the neutrino flux.We use the NuFlux [57] package to interpolate those tables for the energy and directions relevant to this analysis.The three experiments considered in our analysis measure the flux at energies below ∼ 100 GeV.The effective volume for SuperK, IceCube Upgrade, and ORCA as a function of the neutrino energy is shown in Figure 4 (top panel).
The zenith distribution of the neutrino flux at the detector shows an enhancement for the horizontal directions due to the longer paths that mesons have to travel before hitting Earth.An example of this effect is shown in Figure 5 (right) for cos(θ zen ) = 0 for both flavor components of the flux at E = 10 GeV, although the same effect also happens at larger energies.For energies where the mesons and muons have decayed before reaching Earth, there is still a horizontal enhancement due to the spherical geometry [58] of the volume where the neutrinos are produced.The absorption of the mesons and muons by the Earth also contributes to modifying the flavor composition of the initial flux.If all the parent particles are able to decay, we can expect that (ν e + ν e )/(ν µ + ν µ ) ∼ 1/2.As the energy increases and muons hit the Earth, losing their energy, this ratio decreases in the meantime.
In the sub-GeV energy range, the flux shows additional anisotropies due to the interaction of the charged mesons with Earth's magnetic field.At low momentum, the mesons produced in the interaction of the cosmic rays get trapped by the magnetic field and can decay into neutrinos, enhancing the flux at lower energies.The trajectory of the primary cosmic rays that reach Earth also gets modified by the mag-netic effects inducing an east-west asymmetry [59] and a dependence of the flux with the location of the Earth [49].Further, the interaction of the meson flux with Earth's magnetic field modifies the neutrino to antineutrino ratio.At lower energies, the multiple scatterings of the mesons wash out the differences between both fluxes.As the energy increases, the geomagnetic effect becomes less important on the meson fluxes, and the neutrino production is dominated by the primary flux mesons.
The uncertainties in the calculation of the atmospheric neutrino flux from cosmic-ray interactions described above arise from four factors: incident cosmic-ray flux, the hadronic interaction model, the atmospheric air density profile, and the magnetic effect at low energies.Precision measurements of the Earth's atmospheric density have been performed by the NASA Atmospheric Infrared Sounder (AIRS) on board the Aqua satellite [60].This data has been used to compute the atmospheric neutrino fluxes and study the effects of seasonal variations [61][62][63].The effects of seasonal variations on the atmospheric muon neutrino flux are less than 10% for neutrinos below a TeV [62] and average out significantly for data sets that span multiple years.Additionally, cosmic-ray spectral measurements in the relevant energy range have been recently performed by AMS.Hadronic interaction models play a significant role in the uncertainty of atmospheric neutrino calculations; see [53,64] for reviews on this topic.The energy range relevant for our analyses is well covered by accelerator data predominantly by the NA49 and NA61 experiments at CERN, though measurements by HARP, PHENIX and STAR are relevant on the lower and higher energy bands considered in this work [53,64].These sources of uncertainty can be translated to bands on the atmospheric neutrino fluxes.As recently estimated [53], at below ≈ 100 GeV, the normalization uncertainty of ν µ + νµ and ν e + νe is below 10% , while at energies below 10 GeV, the precision on the ratio of muon-to-electron neutrinos is known with a at an error below 2%.Additionally, at energies below ≈ 10 GeV, the ratio of neutrinos to antineutrinos for muons is known with a precision between 1% and 5%, while for electrons this ratio is known to a precision of 10%.
To include all the uncertainties mentioned before, we use a similar parametrization of the flux as in Ref. [64,65] where f α (E, cos ζ) are the Honda table's values in-  For easy comparison, we use a systematic uncertainty budget that is consistent with recent Su-perK [26] and IceCube [66] analyses.However, the choice of uncertainty parameterization used in this article is conservative given the discussion above, and is expected to be further improved by further measurements.In Figure 5, we show the 1σ range of the energy (left) and zenith (right) uncertainties used in our analysis.

III. NEUTRINO WATER CROSS SECTION
The experiments in this work share, in addition to the neutrino flux, the neutrinos cross sections in water.These experiments cover a wide range of neutrino energies in which neutrinos may interact via the exchange of charged or neutral currents.The relevance of distinct interaction channels differs from one experiment to another, since they measure the atmospheric flux at distinct energy scales.Therefore, we need to study the different interaction channels as they may affect the determination of oscillation parameters in a different way in each experiment.In this section, we will provide a concise overview of the essential aspects of the neutrino cross-section that bear relevance to the determination of oscillation parameters.
Charged-current interactions produce a charged lepton that shares flavor with the original neutrino and are divided into three major contributions shown in Figure 6.
• Charged current quasi-elastic (CCQE): These interactions dominate in the lower energy region, below 2 GeV, and scatter off one of the bound nucleons, exchanging a W ± boson, emitting the charged lepton partner of the interacting neutrino.The outgoing nucleon is either a proton, for neutrinos, or a neutron, for antineutrinos.These interactions are most relevant in FIG. 6. Charged-current νµ cross section per nucleon as a function of energy split by its main contributions, QE, RES, and DIS.The cross section models used are from GENIE pre-computed splines, shaded regions correspond to the 1σ priors assumed in Table II, and data points correspond to the most relevant inclusive cross section measurements for this work.
the sub-GeV region of the SuperK atmospheric neutrino dataset and the lowest energy bins for IceCube Upgrade and ORCA, where the sensitivity to the CP -phase resides.Additionally, this channel provides a clear link between the matter-antimatter character of the incoming neutrino and the presence of a neutron in the final state, which becomes relevant when introducing neutron-tagging capabilities in SuperK.
• Resonance production (CC RES): At slightly higher energies up to 4 GeV, neutrinos can excite an entire nucleon, producing a baryon resonance that in turn, quickly decays into a nucleon and single or multiple mesons.The ∆(1232) baryon resonance dominates this channel, producing a single pion in the final state.All three experiments are sensitive to this channel.Similarly to CCQE, in these interactions, antineutrinos tend to produce more neutrons than protons in the final state.Further, in singlepion production, neutrinos are linked to π + and antineutrinos to π − .In turn, oxygen atoms in water quickly absorb π − before decaying, providing a potential signature to separate neutri-nos from antineutrinos from the reconstruction of π + decay products, namely Michel electrons.The SuperK detector is sensitive enough to use these features, which in this energy region, provide sensitivity to the neutrino mass ordering from the core and mantle resonances between 2 GeV and 10 GeV from MSW effect of neutrinos propagating through Earth (Figure 9).Beyond resonance production, single pion final states can also be produced from neutrinos coherently scattering the whole nucleus (Coh π).
• Deep inelastic scattering (CC DIS): At energies above 4 GeV, neutrinos can scatter off a single quark inside the nucleon, producing the corresponding charged lepton plus a hadronic shower in the final state.In this type of event, the flavor reconstruction may get confused unless the outgoing charged lepton has sufficient momentum to be distinguished from the hadronic showers.
Given the large mass of the τ lepton, this is the only channel allowed for charged-current tau neutrino interactions with a threshold of 3.5 GeV.This channel dominates the neutrino interactions measured at IceCube Upgrade and ORCA, and thus the sensitivity to the neutrino mass ordering through the second Earth's matter-effect resonance and to the atmospheric squared mass difference.
On the other hand, neutral-current interactions do not produce an accompanying charged lepton but a neutrino and, therefore, do not carry any information about the flavor of the incoming neutrino in water-Cherenkov detectors.Thus, these interactions are a background for the neutrino flavor oscillation analyses.For neutral currents, only analog resonance and deep-inelastic scattering channels can be reconstructed in water-Cherenkov detectors if the produced mesons or hadronic showers are sufficiently energetic, or if the products decay to other detectable particles.An example of these interactions is the production of a π 0 decaying to a pair of photons, which are a source of background for charged-current electron-like (e-like) events below 1 GeV in SuperK.
Besides, some products of primary interactions undergo secondary interactions within the nuclear media, making the final state of the neutrino interaction more complex and difficult to reconstruct and classify.Secondary interactions become more prominent at higher energies and obscure some distinct signatures expected from primary interactions, like the content of pions and neutrons, weakening the flavor identification as to the separation between neutrinos and antineutrinos in the case of SuperK.
Over the last four decades, several experiments have measured neutrino cross sections with different targets, and over a wide range of energies and channels relevant to this work.The first measurements were done in hydrogen and deuterium bubble chambers [67][68][69].These experiments provided precise measurements but lacked the complexity and features of heavier target nuclei required in current and future neutrino oscillation experiments.The structure and kinematics of nucleons as well as the re-scattering within the nucleus can potentially alter the neutrino interaction cross section and outcome, affecting the measurement of the oscillation parameters [70][71][72].As neutrino physics developed, neutrino cross section measurements on the relevant targets were necessary to account for these nuclear effects.
One of the remaining issues is to achieve the same precision for electron neutrinos and antineutrinos, for which the data is much scarcer since the beams are mainly made of muon neutrinos.
In addition to the aforementioned Minerνa and T2K, there are other experiments with water targets in operation, such as ANNIE [103] and NINJA [104,105].These experiments, together with near future detectors like Hyper-Kamiokande's Intermediate Water-Cherenkov Detector (IWCD) [28] and the upgrade of the current T2K's near detector ND280 [106], will play a crucial role in reducing the uncertainties of both, muon and electron neutrino and antineutrino cross sections as well as differential cross sections at energies from hundreds of MeV to a few GeV on a water target.Additionally, the upgrade of ND280 will measure the ratio between carbon and oxygen to better precision, enabling a more precise extrapolation of cross section measurements on carbon and hydrocarbon targets.Moreover, the tagging of neutrons from neutrino interactions is proven to be an effective tool for separating neutrinos from antineutrinos.It is expected that the ANNIE and IWCD experiments will be Gd-doped detectors exposed to high-intensity neutrino beams, thus providing valu-able input to reduce the current neutron production uncertainties.
In addition to more experimental measurements, the development of more precise nuclear and interaction models is crucial.In this context, two topics will be of great importance in the coming years for neutrino physics.First, the model-independent reassessment of bubble chamber data and, second, the deeper study of nucleon form factors and resonance production from updated models and lattice QCD calculations [107][108][109].The latter work would be especially relevant for the CCQE and RES channels.
For charged-current tau neutrino data is much more limited; it comes primarily from the DONuT [110] and OPERA experiments [111] and, more recently, from tau appearance measurements in the atmospheric neutrino flux from IceCube and SuperK [112,113].These measurements found the cross section to be in good agreement with the expectations within 10% of the expected value.Tau appearance does not play a significant role in extracting the parameters of interest in this article and thus we do not discuss them further.
The systematic uncertainties assumed in this work are summarized in Table II and follow the ranges assumed by the official oscillation analyses from each of the experiments considered in this work.During the upcoming data-taking period of IceCube Upgrade, ORCA, and SuperK, it is expected that the neutrino cross section uncertainties will further improve with new data from experiments in operation and theoretical development.These measurements and theoretical developments are motivated by the upcoming next-generation neutrino detectors, DUNE and Hyper-Kamiokande.

IV. PHYSICS OF WATER CHERENKOV DETECTORS
Due to their large active and instrumented volumes, water-or ice-Cherenkov detectors have proven to be a very effective and successful technology for measuring the properties of neutrinos.
In these kinds of detectors, neutrinos interact with the nuclei in water molecules producing several particles depending on the primary interaction channel and subsequent secondary interactions within the nuclear media.Charged particles with momenta above a given threshold determined by the refractive index of the medium emit Cherenkov radiation detected with photomultiplier tubes (PMTs).Ice and water have sligtly different refractive indexes [121], and thus, Cherenkov energy thresholds; namely 0.13 GeV and 1.2 GeV in Antartic ice and 0.16 GeV and 1.4 GeV in water for muons and protons, respectively.Furthermore, both media exhibit well-understood light-propagating properties, enabling reliable event reconstruction.This is clearly shown by all three experiments considered in this work.On the one hand, thanks to its large photo-coverage and the exhaustive control of the water conditions, SuperK can reconstruct very low-energy events (MeV-scale), which gives a comprehensive reconstruction of atmospheric neutrino events at energies below 1 GeV.On the other hand, IceCube Upgrade and ORCA experiments are capable of precisely parametrizing the optical properties of Antarctic ice and Mediterranean seawater, so it is possible to control such large volumes of unprocessed media and use them as particle physics detectors with energies as small as 1 GeV.In this context, the ORCA experiment provides a better directional reconstruction than IceCube where there is more scattering of photons due to the presence of air bubbles trapped in the ice.As charged particles propagate through a medium, they polarize the medium around them.When these particles move faster than light in that medium, photons from the polarization, unable to keep up with the particles' pace, form a characteristic cone of light around the direction of motion, with a characteristic opening angle depending on the medium's refractive index.The detection of this radiation allows the reconstruction of the direction and production vertex of the charged particle.Additionally, being proportional to the number of photons emitted, the momentum of the particle is inferred from the charge collected by the PMTs.
In the context of neutrinos, this means WC detectors can reliably reconstruct the promptly charged leptons that originate from a neutrino chargedcurrent interaction, namely e ± and µ ± from electron and muon (anti)neutrinos, respectively.Moreover, in large-photo-coverage experiments, some particles are identifies based on the ring patterns; electrons produce diffuse ring patterns due to the electromagnetic showers produced as they propagate through water (e-like or showers), whereas muons, being more massive, produce rings with sharper edges (µ-like or tracks).Other particles such as high-energy photons and charged pions can be detected and reconstructed similarly, the former with a pattern practically indistinguishable from showers and the latter from tracks.All these features extend the capabilities of neutrino WC detectors, allowing for a complete reconstruction of the particles in the final state of the neutrino interaction as well as the decay products from heavier short-lived particles and other charged particles produced in secondary interactions.
There are, however, differences in the reconstruction among the detectors considered.In the case of SuperK, the Cherenkov cone gets projected to the densely instrumented inner surface, which registers all the signals produced in the event.In the case of IceCube and ORCA, strings of PMTs are scattered throughout the volume, enabling the measurement of the energy loss of particles as they travel through the ice (IceCube) and seawater (ORCA) at each step of the track.See Refs.[122][123][124] for recent discussions on reconstruction on these type of detectors.
An additional yet relevant point is the detection of neutrons in WC detectors.Currently, in the context of atmospheric neutrinos, this only applies to SuperK.Despite having a null electric charge, neutrons get captured by hydrogen atoms of water, forming an excited state of deuterium.Its de-excitation emits a 2.2 MeV photon which is detected with an efficiency of ∼20%.Despite the modest tagging efficiency, this capability adds more information to the reconstructed final state of the interaction, improving the oscillation analysis sensitivity, which we will discuss in Section VI B.
Additionally, given the proven relevant of neutron tagging [125], SuperK is being upgraded to dissolve a gadolinium (Gd) salt in the water to improve the detection of neutrons through their capture on this element [126].Gadolinium is the stable isotope with the largest thermal neutron absorption cross section, which, together with the emission of an 8 MeV cascade of photons, as shown in Figure 7, enhances the neutron detection efficiency up to ∼80%.

V. NEUTRINO OSCILLATIONS
In the 3ν scenario, neutrino evolution is described by six parameters: two mass-squared differences (∆m 2  31 and ∆m 2 21 ), three mixing angles (θ 12 , θ 13 , and θ 23 ), and a complex phase that parameterizes the violation of the CP -symmetry in the lepton sector.The so-called solar parameters (∆m 2  21 and θ 12 ) have been measured by solar experiments and Kam-LAND [127][128][129][130] looking for the disappearance of the electron neutrinos or anti-neutrinos.Reactor experiments with a baseline of O(10 0 km) and using a configuration with a near and a far detector has determined θ 13 , becoming the most precisely measured parameter to date [131][132][133].Finally, the atmospheric parameters (∆m 2  31 and θ 23 ) and the CP -violation phase require high-energy beams and, therefore are constrained by long-baseline experiments [24,[134][135][136][137].Using a ν µ flux with energies at or near the GeV scale, those experiments search for the disappearance of muon neutrinos (ν µ → ν µ ) and for the appearance of electron neutrinos (ν µ → ν e ).In the case of muon-disappearance channel, since the matter effects are suppressed by sin 4 θ 23 [138], the oscillation probability at leading order [139] can be written as where ∆m showing that we can constraint ∆m 2 µµ and sin 2 2θ 23 by looking at the muon distribution in the detector.
In the appearance channel the matter effects are more important.To describe long-baseline experiments (LBL) sensitivity using this channel, we use an expression valid in a constant matter neutrino evolution [140][141][142] given by , where V is the matter potential.The dependence of this channel on the mass ordering are proportional to the matter effects.
The ν e appearance also depends on sin 2 θ 23 bringing the possibility to resolve between both θ 23 octants.Also, as we see from equation Equation ( 4), using this channel we have a dependence on θ 13 .Moreover, the possibility of long-baseline experiments in running in neutrino and anti-neutrino modes (or in the atmospheric case, having a mix beam of neutrinos and antineutrinos) enables these experiments to resolve the differences in the oscillation of both modes.The comparison of the neutrino and antineutrino oscillation patterns can be translated into a measurement of δ CP .
The latest results of the global analyses indicate that some of those parameters can be constrained to the percent level [9,143,144], although there are still several sizeable uncertainties in 3ν mixing scenario.Among the less constrained parameters, we have θ 23 , for which values above and below maximal mixing are allowed within 1σ; for the mass ordering, we have a 2σ preference for normal ordering (NO) [145]; and δ CP , for which just a small region around π/2 is excluded by T2K and SuperK at 3σ CL.Atmospheric neutrinos can contribute to narrowing down those uncertainties.

A. Neutrino evolution through the Earth
The large range of baselines covered by atmospheric neutrinos that extend from O(10 km) to O(10 4 km), and the vast energy range where the flux can be measured from O(10 −2 GeV) to O(10 5 GeV) ensures the access to a vast neutrino oscillation phenomenology, as we can see in Figures 8 and 9.In the sub-GeV region, and for baselines ∼ 1000 km, the neutrino evolution is dominated by ∆m 2  21 .The finite energy resolution makes the experiments inaccessible to the oscillation driven by ∆m 2  31 and ∆m 2 32 .The average over those two oscillatory terms enhances the effects of δ CP over neutrino evolution [146,147].The asymmetry between neutrino and anti-neutrino oscillation comes through the Jarlskog invariant [148,149] defined as , where we have factorized the dependence with δ CP .In vacuum, the CP -violation term is given by the product of the three oscillation wavelengths [150] and the Jarlskog invariant; namely where ∆ ij = δm 2 ij L/4E.The average over the two largest mass-splittings suppresses the P CP term by ∼ 1/2 [151].Figure 10 shows the electron-and muonappearance probability for two values of δ CP = 0 and π.In the case of CP -conservation, both probabilities are the same, as seen in the top and bottom panels of the figure.As δ CP takes a different value, both appearance channels separate, increasing the sensitivity over that parameter [152].The CP -violation results results in a different normalization of the probability and a shift in the oscillation phase: these effects are quite broad in energy and more relevant specifically in the sub-GeV energy range.The matter effects depend on the neutrino trajectory along the Earth and create a dependence on the CP -effects with the neutrino direction.
As the neutrino energy rises, the matter effects become more important.At the GeV energy scale and for trajectories crossing the mantle, there is an enhancement of the effective mixing angle θ13 due to the coherent-forward elastic scattering of the neutrinos with the electrons in the Earth, the socalled MSW effect [7,153].The matter-modified mixing angle is given by sin 2 θ13 = sin 2θ 13 (cos θ 13 − 2EV /∆m 2 31 ) 2 + sin 2 2θ 13 .
(6) For energies around 6 GeV and densities around 5g/cm 3 , sin 2 θ13 becomes maximal, giving rise to an enhancement of the flavor conversion: see in Figure 11.The location of the resonance is controlled by sin 2 θ 13 for a given value of ∆m 2  31 , providing sensitivity to this angle.As sin 2 θ 13 becomes larger, the resonance moves to lower energies and densities.The opposite happens if sin 2 θ 13 becomes smaller, we will need larger energies and densities to meet the resonant condition.For very small values of sin 2 θ 13 , the oscillation length at the resonance becomes larger than the size of the Earth [138], making impossible to get a large flavor conversion with atmospheric neutrinos, as is the case with the blue line in Figure 11.
Beyond the MSW effect, at energies around 1 GeV, the oscillation is enhanced for some trajectories crossing Earth's core and mantle due to the matter distribution, the so-called parametric resonance [65,[154][155][156][157][158][159][160][161][162][163][164][165].Both types of resonance happen for neutrinos if the mass ordering is normal and for anti-neutrinos in the case of inverted ordering, as shown in Figures 8 and 9.The differences in the flux and the cross section for both fermions bring the possibility to measure the neutrino mass ordering using atmospheric neutrinos.
In the multi-GeV scale, the neutrino oscillation length gets longer, and the neutrino oscillation is dominated by ∆m 2  31 and θ 23 .The first oscillation minimum for P (ν µ → ν µ ) happens at E ∼ 20 GeV for baselines that cross the Earth: see in Figure 8.The energies where that oscillation minimum happens depend on |∆m 2  31 |, and the oscillation amplitude is controlled by sin 2 2θ 23 , as shown in Figure 12 [166].It is important to notice that the octant of θ 23 can be measured with atmospheric neutrinos in two ways as shown in Figure 12: similarly to LBL experiments through the electron appearance channel as it is proportional to sin 2 θ 23 , Equation ( 4), and through muon disappearance as matter effects break its dependence with sin 2θ 23 due to the enhancement of sin 2 θ13 .
In this analysis, we numerically solve the neutrino evolution using nuSQuIDS [167].For the Earth matter density, we used the Preliminary Reference Earth Model [168] (PREM), which divides the Earth into eleven concentric layers, and for each of the layers, the density is given by a polynomial function that depends on its distance to the center of the Earth.As a benchmark scenario for the mixing parameters, we followed the Nu-Fit results [143].For the solar parameters, we used ∆m 2 21 = 7.42 × 10 −5 eV 2 and sin 2 θ 12 = 0.304 that are fixed in the whole analysis.The reactor angle is fixed at sin 2 θ 13 = 0.022 unless otherwise specified.For the atmospheric parameters, we assumed ∆m 2 31 = 2.5 × 10 −3 eV 2 and sin 2 θ 23 = 0.572, and for the CP -violation phase δ CP = 234 o .This is summarized in Table III.

VI. EXPERIMENTAL ANALYSES
To analyze the sensitivity to the neutrino oscillation parameters, we produce an Asimov dataset for each experiment, and perform a combined fit to the Monte Carlo simulation assuming different values of the oscillation parameters sampled from two 4-dimensional (∆m 2  31 , θ 23 , θ 13 , and δ CP ) grid of points, one for each neutrino ordering.To reduce the impact of stochastic uncertainties from the simulations, we produce large-exposure Monte Carlo datasets for each experiment.We bin the events in observable quantities, which differ from experiment to experiment and are discussed in the following subsection.We then construct the following test statistic, χ 2 , to compare data with prediction.This is given by Equation ( 7), where µ i and O i are the expected number of events in bin i th respectively, and f ij is the fractional change in the number of events in bin i th due to the j th systematic source of uncertainty which takes value η j and has ηj nominal value with σ j error size.The first term in Equation ( 7) is the negative times two of the log-likelihood ratio between the Asimov dataset and the MC at a given point in the oscil-lation parameter space assuming Poisson statistics, and the last term takes into account the penalty from each source of systematic uncertainty, for which gaussianity is assumed.The Asimov dataset is fitted against the MC using a binned χ 2 method assuming the number of events per bin follows a Poisson statistic.For introducing the effect of systematic uncertainties, the number of entries in each bin is FIG.12. Electron appearance (left) and muon disappearace (right) probabilities for cos(θzen) = −0.85.We show impact of sin 2 θ13 in both oscillation channels.The fast oscillations have been smeared assuming a Gaussian uncertainty of 5% E/GeV.re-weighted accordingly and an additional penalty term is introduced in Equation ( 7), see Ref. [169] for details.
Minimizing Equation ( 7) requires solving the system of equations defined by ∇χ 2 = ⃗ 0 over all systematic uncertainties.This procedure brings the Asimov dataset and MC into the best agreement allowed by the size of systematic uncertainties.Usually, this process requires extensive CPU resources as it relies on the numerical computation of partial derivatives with respect to all systematic sources, which has to be repeated for all points in the defined grid of oscillation parameters to be tested while managing large simulation files.This problem can be alleviated by analytically computing the Jacobian of Equation ( 7), at the cost of implementing in the code the partial derivatives of the expected number of events in each bin for each systematic source, as shown in Equation (8).Such implementation improves the convergence time to a minimum by almost two orders of magnitude in this analysis, and provides more robust minimization.For the minimization of χ 2 we use the SciPy Python package [170].
, where The events are divided into the samples according to the definitions of each experiment and binned by their reconstructed zenith angle and energy as explained in Section VI A, Section VI B, Section VI D, and Section VI E.
The f ij are introduced by computing them on run-time before the fit on an event-by-event basis or read-out from publicly available tables.Unlike prior global analyses of atmospheric neutrinos, flux and cross section systematics are common to all experiments considered, while detector systematics are only applied to each experiment.
As mentioned earlier, in this work, we consider four experiments: SuperK, IceCube Upgrade, ORCA, and HyperK.Additionally, we consider the three distinct phases of operation within SuperK: SuperK without H-neutron tagging (SuperK), SuperK with H-neutron tagging (SK-Htag), and SuperK with gadolinium (SKGd).The phases of SuperK are treated in the analysis as three independent data-taking experiments, with uncorrelated detector systematics.The first phase of SuperK covers the first three runs of the experiment, from SuperK-I to SuperK-III; SK-Htag covers the detector from SuperK-IV to SuperK-V with neutron tagging capabilities on hydrogen; finally, SKGd covers the projected running time of SuperK with gadolinium enabling improved neutron tagging.In terms of exposure, we assume the full data-taking periods through SuperK-I to SuperK-V as reported in [30].For SKGd, we project five years of operation with the final concentration of Gd dissolved in water, 0.2%, starting in 2025.
For the soon-to-be-deployed IceCube Upgrade and ORCA, we conservatively foresee five and three years of operation starting in 2025 and 2027, respectively.For HyperK, we assume 2.5 years of operation starting in mid 2027 as foreseen by the collaboration.
Thus, the combined analysis in the report assumes the running of SuperK until the last reported exposure time, the SKGd and IceCube Upgrade datataking period extending from 2025 until 2030, and ORCA from 2027 to 2030.Additionally, HyperK is projected to be completed by mid-2027, ensuring a large amount of statistics by 2030, and thus its included in the fit.The situation is less certain for DUNE, as its final volume has not been defined, for this reason it is not included in the fit.Nonetheless, its expected performance is described in Section VI F. According to the current status of construction, this fit conservatively considers the neutrino oscillation measurement picture these experiments will provide by the end of the decade.

A. Super-Kamiokande
SuperK is a 50-kton cylindrical water-Cherenkov detector located in Kamioka, Japan.The detector is instrumented with more than 11,000 20-inch PMTs in the inner surface and facing inwards.The mountain above shields the experiment from most of the cosmicray muons and the large photo-coverage enables the measurement of low-energy atmospheric neutrinos up to 100 MeV [171].
SuperK started its operation in 1996 and has largely contributed to the current knowledge of neutrinos; particularly, it contributed to the discovery that neutrinos have a definite mass by measuring oscillations in atmospheric neutrinos [172].Over the years, the experiment went through various phases, from SuperK-I to the current SuperK-VI [30].From 2020, the experiment is undergoing a major upgrade going from an ultra-pure to a Gd-doped water Cherenkov detector [173], which will improve its capabilities by enabling highly efficient neutron tagging on gadolinium.
From the latter and according to to [45], SuperK has been able to constrain the values of the atmospheric parameters, ∆m Recently, due to the knowledge of the detector, the SuperK collaboration was able to extend the fiducial volume by 20%, from the usual 22.5 kton to 27 kton in all running periods [30], enlarging the sample size and thus its sensitivity.We assume this improved fiducial volume in our analysis.
For this work, we develop a detailed Monte Carlo simulation to predict the atmospheric neutrino event rate at SuperK.First neutrino events interacting with a water target are generated using the GENIE event generator [187].Then, we implement software emulating the SuperK reconstruction procedures, efficiencies, and resolutions based on publicly available information.In the end, we produce a simulation equivalent to 300 years for each phase of the experiment following the event categories defined in each case.
The official atmospheric neutrino analysis from SuperK separates events into three categories depending on their morphology: fully contained (FC), partially contained (PC), and upward-going muons (Up-µ).In this work, we only simulate the first two as the latter has very mild impact on the sensitivity and is statistically dominated by IceCube Upgrade and ORCA.
The reconstruction is focused on the FC events as they carry the most information and are the most abundant in the energy region sensitive to the oscillation parameters, from 0.1 GeV to 20 GeV.Fully contained events get categorized into Single-Ring or Multi-Ring depending on the number of rings, Sub-GeV or Multi-GeV attending to their reconstructed visible energy, and e-like or µ-like in terms of the reconstructed ID of the most energetic ring.In addition, depending on other reconstructed variables, events get further divided to distinguish neutrinos from antineutrinos, and charged-current from neutralcurrent interactions.Single-Ring Sub-GeV samples are divided into 0-decay-electrons and > 0-decayelectrons for e-like, and 0-decay-electrons, 1-decayelectrons and > 1-decay-electrons for µ-like, where the decay-electrons are the number of delayed electrons reconstructed from the decay of charged pions and muons in the event.
An additional Single-Ring Sub-GeV π 0 -like sample is defined for e-like events passing a π 0 cut, accounting for a fraction of neutral current events producing a π 0 which decays to a pair of photons, but only one of them is reconstructed.For these events, another sample is defined in the case that both photons are reconstructed, Sub-GeV 2-Ring π 0 -like.Single-Ring Multi-GeV e-like events follow the same criteria as their Sub-GeV counterparts, and no further classification is done for µ-like events.Finally, for Multi-Ring Multi-GeV e-like events, a two-step classification is done.The first one tags most of the neutral-current events into the Multi-Ring Other sample.The remaining events are divided into neutrino and antineutrino samples.
As expected from Section V, Sub-GeV samples are those most sensitive to the CP -violating phase.The effect from different values of δ CP in the most relevant samples is shown in Figure 13.Additionally, in Figure 13 we also show the power of Multi-GeV e-like samples to discern between normal and inverted orderings.Despite their name, ν-like samples have still more neutrinos than antineutrinos; this contamination of neutrinos comes due to the modest neutrino-antineutrino separation power and the smaller ν cross section.Therefore, the effects of different values of δ CP or mass ordering scenarios are diluted.More details about the SuperK simulation can be found in Appendix A.
Detector systematics have a secondary effect on the sensitivity to the oscillation parameters, except δ CP , as compared with those of the flux or the neutrinonucleon cross section.To realistically asses the sensitivity of the experiment, we implement a subset of the most relevant SuperK detector systematic uncertainties as detailed as possible within the reach of our simulation.Our implemented systematics follow [45] and are summarized in Table IV for each of the three major detector phases aforementioned.

B. Super-Kamiokande with neutron tagging
After the electronics upgrade in 2008, SuperK was able to lower the signal threshold enough to be sensitive to the 2.2 MeV gamma from the de-excitation of deuterium after neutron capture on hydrogen.Despite being very weak, this signal is reconstructed with an efficiency of ∼20%.
Further, in 2018, the upgrade of the SuperK detector started to make it compatible with dissolving Gd.At 0.2% concentration by mass of Gd sulfate in water, 90% of the thermal neutrons are captured by Gd with a half-life of ∼30 µs, emitting an 8 MeV γ cascade.This is detected by SuperK with an efficiency of 90%, having a final neutron tagging efficiency of ≈80%.Currently, SKGd is running at a third of the goal gadolinium concentration, which is expected to be achieved in the next years [188].
Neutrinos produce on average fewer neutrons than antineutrinos; then in addition to the usual cut in the number of electrons from muon decays for single-ring samples, a cut is applied in the number of tagged neutrons (0 neutrons or ≥1 neutrons) to improve the neutrino-antineutrino separation.This establishes new sample definitions for the analysis of SuperK-IV data.
The SuperK atmospheric neutrino analysis includes this neutron information by defining new event samples to improve the separation between neutrinos and antineutrinos.The previously explained Single-Ring e-like events with 0 decay-electrons are divided into samples with 0 tagged neutrons or one or more tagged neutrons, and similarly, for µ-like events as described in [189].
This improved event classification enhances the sensitivity to those oscillation parameters behaving differently for neutrinos and antineutrinos, namely the CP -phase and the neutrino mass ordering.These effects are shown in Figure 14, for the case of Gdneutron tagging.

C. Hyper-Kamiokande
HyperK is a next-generation water-Cherenkov neutrino and proton-decay experiment in Japan, ex-FIG.13.Impact of different values of δCP (left) and both neutrino mass orderings (right) for the most sensitive SuperK samples.On the left, R µ/e (δCP ; 0) are the ratio of the number of events in Sub-GeV µ-like and e-like samples, with 1 and 0 decay-electrons respectively, for a given value of δCP compared with the case of δCP = 0. On the right, we show the ratio between inverted (NIO) and normal (NNO) orderings for the number of events in SuperK Single-Ring and Multi-Ring Multi-GeV e-like samples.The rest of parameters follow Table III  pected to take over the current SuperK detector in 2027.It will be a scaled-up version of SuperK consisting of a cylinder of 74 m in diameter and 60 m in height, with a total mass of 258 kt [28].The detector design is very similar to that of SuperK, with inner and outer detectors, the fiducial volume is 187 kt, a factor 8.4 larger than that of SuperK.The inner detector will be instrumented with approximately 20,000 20-inch PMTs and additional ∼1,000 multi-PMT modules [190,191].The photo coverage will be 20%, which is half of SuperK, but the performance is expected better due to improved photosensors.

FIG. 14.
Impact of different values of δCP (left) and both neutrino mass orderings (right) for the most relevant SKGd samples.On the left, R µ/e (δCP ; 0) are the ratio of the number of events in Sub-GeV µ-like and e-like samples, with 0 Gd-tagged neutrons, and 1 and 0 decay-electrons respectively, for a given value of δCP compared with the case of CP -conservation.On the right, we show the ratio between inverted (NIO) and normal (NNO) orderings for the number of events in SKGd Single-Ring Multi-GeV e-like samples.The rest of parameters follow Table III.
Given the similarities between SuperK and Hy-perK, we used our SK-Htag simulation as the baseline Monte Carlo for HyperK.This assumption provides a conservative estimate of the detector sensitivity as preliminary studies have already shown that neutron tagging on hydrogen will be more efficient in HyperK than in SuperK, see Ref. [192,193].
In the analysis, the HyperK simulated data-set is divided into the same samples as the SK-Htag, using two times the binning in zenith angle to account for the larger statistics.Following the same logic, we conservatively assume the same detector systematic errors as in for SuperK-IV, Table IV.

D. IceCube and the IceCube Upgrade
IceCube is a one km 3 ice-Cherenkov detector located at the South Pole [194], consisting of 5160 optical sensors (DOM) deployed at depths between 1450 m and 2450 m in the Antarctic glacier.The detector can observe neutrino interactions with energies beyond ∼ 10 GeV.The low-energy part of the atmospheric spectra is measured by IceCube-DeepCore (DeepCore), a denser sub-array of strings in the inner part of the detector of ∼ 10 Mton.For energies around ∼ 10 GeV, two event morphologies can be differentiated in the detector: Tracks, which are gen-erated by the propagation of muons through the ice, and Cascades, produced by the propagation of electrons, taus, hadronic cascades, or electromagnetic cascades.
As described in Section V, IceCube has contributed to the measurement of ∆m 2  31 and sin 2 θ 23 .Using the first three years of data [66,195], the sensitivity obtained is ∆m 2 31 = 2.55 +0.12 −0.11 × 10 −3 eV 2 and sin 2 θ 23 = 0.58 +0.04 −0.13 .These results are systematiclimited, further evaluation of the detector response together with a larger data sample predicts a large improvement in the precision over those parameters.
IceCube is planning an upgrade [196][197][198] that will consist of the deployment of additional strings to increase the total volume, leading to a large sample size with an increased energy range that will extend to lower energies.The upgrade will also increase the number of strings in the volume surrounding DeepCore, reducing the separation between the optical sensors and lowering the energy threshold to ∼ 1 GeV.The new range of energies will increase the total amount of events observed thanks to the soft energy spectra of the atmospheric neutrino flux, as we can see in Figure 15.The cascade sample dominates the low-energy sample due to the larger detector response at low energies.
The upgraded detector will be able to accurately explore the atmospheric parameters (∆m 2  31 , sin 2 θ 23 ).In Figure 15, we show the impact that several values of the mixing parameters have on the event distribution for cos θ ∈ [−1, −0.8].In two upper figures, we see the impact that ∆m 2 31 will have on the cascade (left) and track (right) distributions.As we see from that figure, the main sensitivity comes from tracks with reconstructed energy around ∼ 20 GeV, which coincides with the first minimum in the muondisappearance channel.As we see in Figure 15, tracks are sensitive to sin 2 2θ 23 and therefore can discriminate between maximal mixing and not, but have a very mild sensitivity to the octant.
The new range of energies accessible by the detector will bring the possibility of increasing the sensitivity over other parameters.As mentioned in the Section V, neutrinos are sensitive to mass ordering at the GeV due to Earth's matter effects.In Figure 15, we see the event distribution for both mass orderings.The most significant difference is obtained for tracks with better angular resolution and energies below 10 GeV.At the GeV scale, enhancing the effective θ13 leads to a significant effect in IceCube.At the bottom of Figure 15, we see how two extreme values for θ 13 will modify the cascade and track distribution compared to the best-fit value measured in reactor experiments.The largest sensitivity comes from tracks with energies below ∼ 20 GeV.
In the case of the IceCube upgrade, there is still not a dedicated analysis of the most relevant detector systematics uncertainties and their impact.For that reason, we have used the same systematics as in the search for ν τ carried out by DeepCore [66] extended to the lower energies of the upgrade.The detector uncertainties account for the different properties of the ice as optical absorption and scattering, the overall DOM efficiency, and the DOM response to lateral and head-on light.In Table V, we have a list of all the detector systematics used in this analysis and the 1σ range.

E. KM3Net/ORCA
A new water-Cherenkov neutrino telescope is under construction in the Mediterranean sea, KM3NeT [199].This experiment will be composed of two detectors: ARCA, a cubic kilometer water Cherenkov detector that can observe high-energy neutrinos from the astrophysical sources in the northern sky, and ORCA, a megaton scale experiment that will be able to explore the neutrino properties by measuring the atmospheric neutrino flux.In this article, we focus on the latter.
For this work, we have developed an independent Monte Carlo for ORCA based on the IceCube Upgrade simulation that takes advantage of the event simulation carried out by the IceCube collaboration.In the simulation of ORCA, we only consider events with reconstructed energy in the range of 1.85 GeV to 53 GeV and keep the true neutrino variables from the IceCube simulation.Then, the reconstructed energy, reconstructed zenith, and event morphology are computed and assigned following the distributions in [200].The Monte Carlo event weights for ORCA are translated from those of IceCube following the ratio of effective volumes between both experiments, as shown in Equation (9).
While the IceCube events are classified only into tracks and cascades, the ORCA collaboration reports a more sophisticated event classification including a third (intermediate) class.To implement this, we reassign morphologies based on the results reported by the ORCA Collaboration in [200].Further details of the development of the Monte Carlo simulation can be found in Appendix B.
For the systematics treatment, we employ a similar set of systematics with the IceCube analysis.We use a same set of neutrino cross section-related systematics; for the detector uncertainties, we assume the same electronics behaviors including DOM efficiency rates, and we adapt the ice-related uncertainties to that of the water in the case of ORCA.
In Figure 16, we show the event distribution in reconstructed energy for both orderings (top) and for three different values of sin θ 23 = {0.3,0.58, 0.7}.Worth noting is that intermediate events dominate the sample for the whole energy range.The new event morphology improves the "purity" in the cascade and track samples.The neutrino mass ordering changes the number of cascades with energies around the  Besides the uncertainties related to the flux and the cross section, the measurement carried out in ORCA will be also affected by the uncertainties in the detector response.We consider a free normalization for each of the event morphologies used in the analysis and a 5% error on the energy scale, Table VI.Additionally, in line with the IceCube methodology for addressing detector uncertainties, we have incorporated a set of systematic factors to consider the absorption and scattering of photons in water, as well as the response of the photomultiplier.

F. DUNE
Liquid argon time projection chamber (LArTPC) detectors have demonstrated a good capacity in the reconstruction of sub-GeV neutrinos.The energy an direction of the incoming neutrino can be reconstructed by detecting the tracks of all charged particles produced after the neutrino interaction and identifying them by their topology and energy loss.DUNE is planning to use a 20 kton detector based on this technology to reconstruct beam neutrinos, although they can also measure the atmospheric neutrino flux.
The reconstruction of the sub-GeV atmospheric neutrinos give us access to determine the CP-phase.Following [146], the event distribution in the reconstructed energy and direction has been estimated by a simulation of the neutrino-argon interaction using an event generator.The uncertainties considered for the outgoing protons are 10% in energy and 10 • in the direction.In the case of the leptons, we have used 5% for the energy and 5 • for the direction.Due to the LArTPC capabilities in identifying low energy charged particles, the events are classify by the number of outgoing visible protons.This allows an statistical separation between neutrinos, which dominates the fraction of the with one proton, and anti-neutrinos, which dominate the 0-proton sample.
Considering a 20 kton detector, taking data for 1 year, we get a sensitivity over δ CP around 1.5σ, Fig. 17.In this analysis, we have included only the uncertainties related to the flux and the detector response.Another set of uncertainties related to the neutrino-argon cross section will also affect this sensitivity.The low number of events expected by the end of this decade, and the large uncertainties, makes that the DUNE measurement of the atmospheric neutrino flux cannot contribute significantly to the determination of δ CP by the end of this decade.Therefore, we have decided not included in our analysis.

VII. RESULTS
In our analysis, we explore the combined sensitivity that the present and near-future generation of atmospheric neutrino experiments will have in determining the oscillation parameters for the 3ν mixing scenario by the end of this decade.Through the simulation of the different phases of SuperK, Ice-Cube Upgrade, and ORCA, and the inclusion of all the previously discussed systematic uncertaintiesflux, cross section, and detector-, we investigate the synergies between the three experiments for the sensitivity to the ∆m 2  31 , θ 23 , δ CP and θ 13 oscillation parameters.
Unless otherwise stated, our study of the neutrino oscillation sensitivity presented here assume the values and treatment shown in Table III.
For baselines comparable to Earth's diameter, neutrino oscillations are sizeable at E ν ∼ 100 GeV and lower.The location of the first oscillation minimum occurs at energies around 20 GeV and depends on the value of ∆m 2  31 .The amplitude of this oscillation is modulated by sin 2 θ 23 .This energy region will be measured with large sample sizes by the IceCube Upgrade and ORCA.SuperK also has some sensitivity to this region from the Multi-GeV Multi-Ring samples, but with smaller sample sizes.The track sample dominates the sensitivity to ∆m 2  31 thanks to better angular resolution.Since the atmospheric neutrino flux is dominated by ν µ , tracks are mainly sensitive to P µµ .
On the other hand, as seen in Section V, the muondisappearance channel is sensitive to sin 2 2θ 23 , so it can only distinguish whether sin 2 θ 23 is maximal mixing or not.To distinguish between both octants, we need to consider the appearance channel, proportional to sin 2 θ 23 .Water detectors have a better angular resolution compared to ice detectors due to the larger number of direct photons arriving at the PMTs.For that reason, ORCA, SuperK and HyperK show a better precision for sin 2 θ 23 .The combined analysis of the three experiments shows that ∆m 2   31   can be measured at ∼ 0.6% and the octant of sin 2 θ 23 can be resolved at more than 3σ in the assumed scenario.Figure 18 shows that ∆m 2  31 is dominated by IceCube Upgrade due to its better energy resolution: see Section VI D.
At the GeV scale, neutrinos crossing the mantle undergo a large flavor oscillation due to the MSW resonance at ∼ 6 GeV.As described in Section V, the resonance happens due to the effective enhancement of the θ 13 mixing angle.This affects neutrinos if the mass ordering is normal, and anti-neutrinos in the inverted scenario.Although atmospheric experiments cannot differentiate between neutrinos and anti-neutrinos on an event-by-event basis, the dif-FIG.18. Two-dimensional 90% confidence level regions for sin 2 θ23 and ∆m 2  31 for SuperK (dashed dark blue), IceCube Upgrade (dashed orange), ORCA (dashed green), HyperK (dashed cyan), and combined analysis (solid violet-red).FIG.19.Sensitivty to sin 2 θ13 for SuperK (dashed dark blue), IceCube Upgrade (dashed orange), ORCA (dashed green), HyperK (dashed cyan), and combined analysis (solid violet-red).
ferences in the flux and the cross section allow for statistical separation.Since this is a very localized process in energy and direction, the events with better angular and energy resolution will contribute more to identifying it, and therefore they will contribute to the neutrino mass ordering sensitivity.This is what we have observed in Figure 15, where tracks show a larger modification under the mass ordering.For the case of HyperK, SuperK and ORCA, we see in Figure 16, Figure 13, and Figure 14 that most of the sensitivity to the ordering comes from e-like and cascades samples in the relevant energy region.In Figure 3, we show the sensitivity to the ordering as a function of operation time by combining the three experiments.It would be possible to obtain a 6σ identification of the ordering after five years of SuperK with Gd in addition to its current exposure, five years of IceCube Upgrade, and three years of ORCA even with the obtained 90% C.L. range over sin 2 θ 23 .
The effective enhancement on θ 13 also provides the option to measure this parameter.Although this parameter has been measured with high precision in reactor experiments [131,132,201] looking for the disappearance of ν e and more recently in LBL [44,202] using the appearance channel, atmospheric neutrinos bring us a complementary way to measure it via the matter effect in Earth.Also, this will constitute the first observation of the MSW effect on Earth.This parameter is measured by the three experiments by separate, and the combined analysis will reach an uncertainty smaller than 20%: see Figure 19.
Finally, the CP -phase is the least constrained parameter, for which almost the entire range is allowed at 3σ; T2K [24] and NOvA [25] are the only experiments have shown some sensitivity.So far, the only atmospheric measurement that shows some sensitivity over δ CP comes from SuperK [30], which excludes values around δ CP = π/5 with a significance slightly above 2σ.This parameter is the main target for the next-generation neutrino experiments.Using the presented atmospheric experiments, the sensitivity over some parameter values can be increased up to 99% C.L.: see Figure 20.Still, the sensitivity is dominated by SuperK and HyperK, but IceCube and ORCA can get a 1σ significance thanks to their low-energy measurements.The δ CP is the parameter that benefits the most from the combined analysis due to the large sample sizes from the IceCube Upgrade and ORCA constraining additionally the flux uncertainties which reduce the sensitivity to this parameter.
Previous studies such as the ones given in Refs.[9,143,203] do not consider the correlation between systematic uncertainties as we have done in this work.
The improved methodology, provides enhanced capacity to determine δ CP and θ 23 , and has lesser effect in other parameters.Throughought this work the lines labelled 'Combined Fit' and 'Trivial χ 2 Sum' name our analysis and a simplified uncorrelated work, respectively.The difference between these two lines showcases the impact of the correlated systematics.
The excluded region for δ CP depends on its true value because it predicts very different event distributions for electron and muon samples in SuperK and HyperK, and also, it modifies the cascade (and intermediate events) distribution in the case of IceCube Upgrade (ORCA).Since almost the entire space is allowed, we have explored what fraction of δ CP can be excluded at some confidence levels and for different true values of δ CP : see Figure 21.The most favorable scenario would correspond to δ CP = 0 or π, where 60% of the whole space can be excluded at 90% C.L. and 30% of the space is explored with a significantly larger than 3σ.
In Table VII, we summarize the projected 1σ regions for all four oscillation parameters considered free in this combined analysis and with the true values from Table III.FIG. 21.Excluded fraction of δCP as a function of the true value of δCP .Lines correspond to exclusion at 2σ (blue), 99% (orange), and 3σ (green) for the combined fit assuming normal ordering and fixed sin 2 θ13 = 0.022.40% of the parameter space can be excluded at 99% confidence level for the CP-conserving scenario.For lower confidence levels, the small deviations in the event number expectations can contribute to increase the sensitive region making the sensitivity more uniform.

VIII. DISCUSSION AND OUTLOOK
Atmospheric neutrinos have played a very important role in determining the oscillation parameters since they were observed by SuperK [172].As highlighted in this article, the large range of energies and baselines covered by the flux crossing Earth leads to a great sensitivity in the determination of the atmospheric parameters (∆m 2  31 and θ 23 ).The matter effect at the GeV scale on neutrinos crossing the core and the mantle also brings the possibility of measuring the mass ordering and the so-called reactor angle (θ 13 ).Finally, at the GeV scale and below, the deviations in the appearance channel and statistical differentiation between neutrinos and anti-neutrinos allow the constraint of the CP -violation phase.
The results obtained in our analysis indicate that soon, atmospheric neutrinos can provide a complementary role in constraining the oscillation parameters and improving, in some cases, the precision obtained by long-baseline experiments.For instance, in the case of mass ordering, the combined analysis of LBL has a 2σ preference for inverted ordering due to the tension in δ CP [204], while the latest results from atmospheric neutrinos show a preference for normal ordering with a significance larger than 2σ [30].The results presented in this study demonstrate that atmospheric neutrinos alone are poised to achieve a significance exceeding 6σ.This suggests that these parameters will be measured with remarkable precision by the conclusion of this decade, as illustrated in Figure 3.
The measurement of the atmospheric parameters is currently dominated by LBL experiments.In the case of ∆m 2  31 , the current sensitivity is dominated by T2K [205] and NOvA [202], and it is known with ∼ 1.1% [143] according to latest global analyses.The sensitivity to this parameter can be largely improved up to ∼ 0.5% using just atmospheric neutrinos.The latest results on sin 2 θ 23 indicate a preference for the upper-octant at a significance smaller than 1σ, and maximal mixing is excluded with less than 2σ significance.If we assume the current best-fit value of the global analyses, atmospheric neutrinos can achieve a 3σ sensitivity over maximal mixing, ruling out the lower octant with a higher significance than the LBL experiments, as shown in Figure 22.These results predict that atmospheric neutrinos will rule out the wrong octant of θ 23 with a large significance by the end of the decade.Most of the sensitivity over those two parameters comes from the region above ∼ 10 GeV dominated by the large sample size of IceCube Upgrade and ORCA measurements and, to a lesser extent, SuperK and HyperK.In that region, the event distribution is dominated by ν µ that interact via DIS.Therefore, flux and cross section systematic uncertainties have a mild effect on the sensitivity as shown in Figure 23.
Until now, the CP -violation phase has been explored with a low significance and, only a small region around δ CP = π/2 is excluded at more than 3σ [143].By the end of the decade, a significant fraction of the parameter space will be excluded by atmospheric neutrinos, as shown in Figure 21.That will also help in resolving the actual tension between T2K and NOvA in the determination of δ CP and sin 2 θ 23 [204] thanks to the possibility to differentiate between the two octants of sin 2 θ 23 , as shown in Figure 22.We undertake a systematic study to examine how potential future improvements in various systematic uncertainties will enhance sensitivity regarding that specific parameter.The sensitivity to the CP -phase is heavily impacted by the uncertainties of the flux, detector and, to a lesser extent, cross section, see Figure 23.The flavor ratio and normalization below 1 GeV of the flux, and those related to the CCQE cross section, are the uncertainties reducing the sensitivity to this parameter.Detector systematics have a similar quantitative impact in the sensitivity as those of the flux.An independent and complementary to long-baseline experiments measurement of δ CP is of utmost importance to boost the precise picture of the 3-flavor neutrino mixing.Ancillary measurements of the low-energy cosmic-ray flux and hadron production in proton and 4  2 He scattering with nitrogen and oxygen nuclei, as well as development of more complete models, would narrow these uncertainties with a significant enhancement in the sensitivity to δ CP from atmospheric neutrinos.Further, as discussed in Section III, the current and projected experiments for measuring the neutrino cross section below 1 GeV needed for next-generation accelerator neutrino experiments will provide valuable input for improving the measurement of the CP -phase with atmospheric neutrinos.
Finally, the measurement of θ 13 is dominated by reactor experiments [131,132,201], but lately, LBL experiments have measured this parameter with a considerable precision [44,202].We have explored the possibility of measuring θ 13 in the atmospheric neutrino flux using the matter effects at the GeV scale.The results indicate that by the end of the decade, it will be possible to reach a 20% precision with the atmospheric neutrino flux, see Fig. 19.Although the precision is not comparable to the reactor measurements, it will certainly be on the same order as LBL measurement, Fig. 15.Exploring the mixing parameters at different energy scales might be a convenient way to search for new physics [206][207][208][209][210][211].
Besides, the resulting scenario also enables atmospheric neutrinos to play a prominent role in the measurement of ν τ cross section.Tau neutrinos are not expected from the unoscillated atmospheric flux, but they are measured in the detectors due to neutrino oscillations and strongly depend on the associated oscillation parameters.This way, such a combined analysis will provide very valuable input for the lower end of charged-current ν τ cross section, see the discussion in Refs.[212,213].

IX. CONCLUSION
In this article, we explored the sensitivity of current and soon-to-operate water(ice)-Cherenkov atmospheric neutrino detectors -namely, IceCube Upgrade, ORCA, SuperK, and HyperK -to determine neutrino oscillation parameters.To simulate these experiments we have developed dedicated Monte Carlo simulations and reproduce with good fidelity their experimental results.We incorporate more than 80 sources of systematic uncertainties and treat the correlated uncertainties between experiments, namely those associated to the common cross section and flux.Through a comprehensive study, we motivate a combined data-fit from these experiments by showing the few-percent level precision that it would provide to the measurement of the remaining oscillation parameters -in particular, θ 23 , ∆m 2  31 -and the neutrino mass ordering, as well as providing a constraint on the value of the CP -phase independent from long-baseline neutrino experiments.Furthermore, the tools and the analysis presented in this work comprise the crucial first step to properly include atmospheric neutrinos in global fits.
Additionally, we identify the synergies among experiments and the common systematic uncertainties diminishing the sensitivity for each parameter.We have an extended discussion of the sources of these uncertainties; in particular, regarding the flux and cross section inputs.It is worth noting that while our analysis uses conservative estimations of the flux and cross section uncertainties, both are expected to be improved through new measurements and further theoretical developments.In other words, this is a motivation for further ancillary measurements and methods, since an improved set of aforementioned estimations will greatly benefit the sensitivity of the analyzed experiments, and of particular relevance for the sensitivity to δ CP .
On the detector side, our analysis is also conservative.The reconstructions used for all the experiments considered in this work use traditional reconstruction techniques.These traditional techniques have been shown to underperform compared to new machine-learning-enhanced reconstruction methods, which show greater accuracy and improved execution time, see, e.g., [123,214].Additionally, in the case of the IceCube Upgrade, the reconstruction used in this work does not take full advantage of the nextgeneration sensors which have improved light collection.Improvements are also expected for ORCA as the detector development proceeds.In summary, the detector systematics assumed in this work can be taken as a conservative baseline that is expected to improve as this decade unfolds.
Finally, we emphasize that the results from a combined fit of atmospheric neutrinos would provide a very valuable input for the next-generation neutrino physics program towards the precise measurement of the CP -violating phase in the lepton sector.This combined analysis nurtures itself from more than 40 years of global expertise and measurements of neutrino-water interactions paving the road for next-generation neutrino-water experiments such as Hyper-Kamiokande and IceCube-Gen2.

FIG. 4 .
FIG. 4. Measurements of the atmospheric neutrino flux as a function of the energy.The total neutrino flux measured by different experiments [46-48] together with the energy range covered by the four experiments (SuperK, IceCube Upgrade, ORCA and HyperK) considered in this analysis is shown in the lower panel.On that figure, we have also included the flux prediction from HKKM2014 model [49].On the top panel, we have the effective volume for the three experiments as a function of the neutrino energy.

FIG. 5 .
FIG. 5. Energy (left) and zenith (right) distribution of the atmospheric neutrino flux.The flux is based on Honda et al. [49] calculation.Around the prediction for each flavor component of the flux at the detector, we show the 1σ uncertainty band.

FIG. 10 .
FIG.10.Electron (top) and muon (bottom) appearance probabilities for cos(θzen) = −0.85.We show the impact of δCP in both oscillation channels.The fast oscillations have been smeared assuming a Gaussian uncertainty of 5% E/GeV.

FIG. 11 .
FIG. 11.Electron appearance probability for different values of sin 2 θ13 .The neutrino direction is fixed to cos(θzen) = −0.85.The fast oscillations that happen for E ∼ GeV have been smeared assuming a Gaussian uncertainty of 5% E/GeV.

FIG. 15 .
FIG.15.Event distribution for different set of values of the mixing parameters for cos θ ∈ [−1, −0.8] as a function of the reconstructed neutrino energy (Er).We use the IceCube upgrade MC to evaluate the neutrino reconstructed energy[198].The cascade distribution is shown on the left and tracks on the right.The statistical 1σ error is shown as a band around the event distribution.

FIG. 17 .
FIG. 17. DUNE sensitivity to δCP .We consider 1 year of data taking and two modules of 10 ktonnes each.

FIG. 16 .
FIG. 16.Event distribution as a function of the reconstructed energy and direction for ORCA after 3 years of data taking.Lines correspond to event distributions (with oscillations) for track, cascade, and intermediate morphology classes respectively.

FIG. 22 .
FIG.22.Present sensitivity to the 3ν mixing parameters Comparison between the latest results from reactor and LBL measurements, and the future prediction from the combine analysis of SuperK, IceCube upgrade and ORCA.

FIG. 23 .
FIG.23.From top to bottom, the impact of each atmospheric neutrino flux (left) and cross section (right) systematic uncertainties on the sensitivity to δCP , sin 2 θ23 and ∆m2  31 .For clarity, solid lines represent the most relevant systematics, while dashed lines represent the rest.

TABLE II .
Summary of neutrino-water interactions systematic uncertainties used in this work.

TABLE IV .
Summary of SuperK detector systematic uncertainties used in this work.

TABLE VI .
Summary of ORCA detector systematic uncertainties used in this work.atmospheric resonance, Figure 16 (top).Also, the normalization of the cascade sample contributes to the measurement of sin 2 θ 23 .The same happens for tracks with energies around 5 GeV.For energies above ∼ 10 GeV, tracks are sensitive to sin 2 2θ 23 , what helps in the separation of θ 23 being maximal mixing or not.