A Novel Approach to Study Atmospheric Neutrino Oscillation

We develop a general theoretical framework to analytically disentangle the contributions of the neutrino mass hierarchy, the atmospheric mixing angle, and the CP phase, in neutrino oscillations. To illustrate the usefulness of this framework, especially that it can serve as a complementary tool to neutrino oscillogram in the study of atmospheric neutrino oscillations, we take PINGU as an example and compute muon- and electron-like event rates with event cuts on neutrino energy and zenith angle. Under the assumption of exact momentum measurements of neutrinos with a perfect e-$\mu$ identification and no backgrounds, we find that the PINGU experiment has the potential of resolving the neutrino mass hierarchy and the octant degeneracies within 1-year run, while the measurement of the CP phase is significantly more challenging. Our observation merits a serious study of the detector capability of estimating the neutrino momentum for both muon- and electron-like events.


Introduction
In the last two years, the field of neutrino physics has significantly advanced by constraining the reactor angle θ 13 . The T2K experiment [1] was the first to report a hint of nonzero reactor angle, followed by MINOS [2] and Double CHOOZ [3] which added up to a confidence level above 3 sigma. It was measured accurately by Daya Bay [4] and RENO [5] in March and April 2012, respectively, reaching 7.7 sigma [6] by the October of the same year.
The relatively large reactor angle opens up opportunities [7] for determining the mass hierarchy, the octant of the atmospheric mixing angle, and the CP phase. The first could be achieved with a medium baseline reactor experiment [8,9] and long baseline accelerator experiments [10,11,12,15] could measure all three of them. Atmospheric neutrino experiments [13]- [35] could offer alternative ways to accomplish the same.
Recent studies have focused on magnetized detectors, which can distinguish neutrinos from antineutrinos [15,18,19,21]. Equipped with this capability, a detector of 50-100 Kton (∼ 10 3 tons) scale is enough to distinguish the mass hierarchy. Large volume water-Cherenkov or ice-Cherenkov detectors of tens of Mton (∼ 10 6 tons) scale could offer an alternative. DeepCore [22], the existing in-fill to IceCube can reach down to energies of O(10) GeV and has recently reported the observation of muon neutrino disappearance oscillations [23] and an electron neutrino flux consistent with expectations [24,25], demonstrating the capabilities of an low-energy extension. DeepCore has also some sensitivity to neutrinos from the MSW resonance region [37] around E ν ≈ 5 ∼ 10 GeV. It can however only partially cover it [17,20] and to really exploit it a lower threshold detector would be needed.
There has been extensive interest [30,31,32,33,34,38,39] recently by the IceCube Collaboration and theoretical community to extend the existing IceCube neutrino telescope [26] with an in-fill array called PINGU (Precision Icecube Next Generation Upgrade) [27] that could detect neutrino events of O(1) GeV. Such a detector opens up the opportunity of detecting more patterns of the atmospheric neutrino oscillation behavior, which is diluted at higher energy scale, especially due to matter effects [29]. A large benefit is the expected high event statistics at low energies. Event rates of O(100, 000) per year from atmospheric neutrinos allow for measurements with small statistical uncertainty. During the preparation of this draft, a preliminary experimental study [42] appeared. In Europe a similar detector to PINGU is being considered as part of the Km3NeT project. Our studies can be transferred to this ORCA -(Oscillation Research with Cosmics in the Abyss) [43].
The expectation of a high statistics sample down to 1 GeV scale makes determining the mixing parameters with atmospheric neutrinos very promising [17,20,30,32]. The paper [30,32,33] adopts oscillograms [29] to depict the structure of oscillation resonances [37,44,45] when atmospheric neutrinos travel through the Earth. The ability to determine the mass hierarchy, the octant of the atmospheric mixing angle, and the CP phase is studied. The event numbers and difference between normal hierarchy (NH) and inverted hierarchy (IH) are shown in oscillograms. In [35], the Bayesian approach is explored in a generic way while the Toy Monte Carlo based on an extended unbinned likelihood ratio test statistic is implemented in [36]. When combined with accelerator experiments, the sensitivities on the octant of the atmospheric mixing angle [13,21,30,46] and the CP phase [47] can be enhanced.
In Sec. 2, we first develop a general framework of decomposing the neutrino oscillation probabilities and the event rates in the propagation basis, and apply it to the symmetric Earth matter profile in order to analytically disentangle the effects of the neutrino mass hierarchy, the atmospheric angle, and the CP phase. In Sec. 3, we calculate and display the event rates that can be observed at PINGU. Based on these results, we try to establish the potential of atmospheric neutrino measurement at PINGU in Sec. 4, while its dependence on the input values of the neutrino mass hierarchy, the atmospheric mixing angle and the CP phase can be fully understood in our decomposition formalism. Finally, we summarize our conclusions in Sec. 5. For more details about the basic inputs, including the atmospheric neutrino fluxes, cross sections, effective fiducial volume of PINGU, and the Earth matter profile, as well as the numerical methods of evaluating the neutrino oscillation probabilities through Earth, please refer to Sec. A

Disentangling Parameters in the Propagation Basis
We first develop a general framework in the propagation basis [55,56] for phenomenological study of neutrino oscillation. It can analytically decompose the contributions of the neutrino mass hierarchy, the atmospheric mixing angle, and the CP phase. This decomposition method can serve as a complementary tool to the neutrino oscillogram [29] for the analysis of atmospheric neutrino oscillations, and can apply generally to other types of neutrino oscillation experiments.

Propagation Basis
In the propagation basis, the atmospheric mixing angle θ 23 [55] and the CP phase δ [56] can be disentangled from the other mixing parameters as well as the Earth matter potential. This can be seen from the effective Hamiltonian, represents the matter effect which is proportional to neutrino energy E ν and the matter potential V (x). The mass differences are denoted as, The lepton-flavor mixing matrix U relates the flavor basis (ν α = ν e , ν µ , ν τ ) and the mass eigenstates (m νi = m i , i = 1, 2, 3), ν α = U αi ν i , (2.4) and can be parametrized as U ≡ O 23 (θ a )P δ O 13 (θ r )P † δ O 12 (θ s ), (2.5) where c α ≡ cos θ α and s α ≡ sin θ α . The solar, the atmospheric, and the reactor mixing angles are labelled as, (s, a, r) ≡ (12,23,13) , (2.6) according to how they were measured. For convenience, we denote the three rotation matrices in (2.5) from the left to the right as O 23 , O 13 , and O 12 respectively. The 2-3 mixing matrix O 23 and the rephasing matrix P δ can be extracted out as overall matrices [56], In this way, O 23 and P δ are separated from the neutrino mass hierarchy, which is encoded in the first term inside the square bracket, as well as the matter effect, represented by the second term. In other words, the atmospheric mixing angle θ a and the CP phase δ are disentangled from the remaining mixing parameters, analytically. This is a significant simplification in the analysis of neutrino oscillation, especially the atmospheric neutrino oscillation that suffers from complicated matter profile, as described in Sec. A.2.
To make it explicit, the original Hamiltonian H can be rotated to the equivalent H ′ in the propagation basis through a similar transformation, There are only four mixing parameters involved in H ′ , the two mass squared differences δm 2 s and δm 2 a , the solar mixing angle θ s , and the reactor mixing angle θ r . Correspondingly, we can define a propagation basis (ν ′ i ) [55,56] that is related to the flavor basis (ν α ) and the mass eigenstates (ν i ) as follows: The transformed Hamiltonian H ′ is the effective Hamiltonian defined in the propagation basis. Once the neutrino oscillation amplitudes, are calculated with the Hamiltonian H ′ in the propagation basis, the oscillation amplitudes in the flavor basis, are simply obtained by the unitarity transformation, This makes the formalism much simpler. We find the propagation basis very useful in the phenomenological study of neutrino oscillation. It allows us to analytically factor out θ a [55] and δ [56] from the numerical evaluation of the oscillation amplitudes that involves many factors and can be very complicated. The contributions of the still unknown neutrino mass hierarchy, the octant of the atmospheric mixing angle θ a , and the CP phase δ are now disentangled from each other. A general formalism based on this feature can help to reveal the pictures behind neutrino oscillation phenomena. This is especially important when the three unknown parameters are under close investigations at current and future neutrino experiments.

Oscillation Probabilities
The oscillation probabilities are measured in the flavor basis. It is necessary to explicitly express the flavor basis amplitude matrix S in terms of its counterpart S ′ in the propagation basis by the unitary transformation with O 23 P δ , namely to expand (2.12). According to the definition of the mixing matrix in (2.5) , O 23 P δ can be explicitly written as, The mixing from the propagation to the flavor basis occurs between the second and the third indices. We can expect the first element of S ′ to be unaffected when (2.13) is combined with (2.12) [56], Note that only the elements among e and µ flavors are shown since they are sufficient to derive all the flavor basis oscillation probabilities, from ν e and ν µ (as well as fromν e andν µ , as shown below). Explicitly we find, where R and I gives the real and imaginary parts, respectively. The dependence on the atmospheric mixing angle θ a and the CP phase δ can be clearly seen in the above expressions. The transition probability into ν τ are then obtained by unitarity conditions, while we neglect contributions from tiny components of ν τ andν τ flux in the atmospheric neutrinos [51]. The oscillation probabilities for antineutrinos are then obtained simply as, by reversing the sign of the matter potential in the Hamiltonian (2.1) and the CP phase δ in the neutrino mixing matrix (2.5), which is identical to the parametrization adopted in Review of Particle Physics [57].

Simplifications with Symmetric Matter Profile
The expressions in (2.16) can be significantly simplified in the approximation of the symmetric or reversible matter profile along the baseline, such as those of atmospheric neutrinos in the earth whose matter profile is approximately spherically symmetric as in PREM [48] adopted in our study. It has been known that [58] the oscillation amplitude matrix after experiencing a reversible matter profile is symmetric in the absence of CP violation. This is indeed the case for the oscillation amplitudes through the Earth in the propagation basis, giving, Based on the above observation, the atmospheric neutrino oscillation amplitudes (2.14) can be further simplified, As a convention, we adopt those elements S ′ ij with i ≤ j. It is now manifest that the flavor oscillation amplitudes S eµ and S µe differ only by the CP phase and the expression for S µµ (2.14d) is greatly simplified in (2.20d). The oscillation probabilities now read, Throughout our studies in this report we adopt the expression (2.21) for computing the oscillation probabilities in our numerical calculation, which are exact in the limit of the symmetric earth matter profile PREM [48] and neglecting the depth of the detector beneath the earth surface as compared to the baseline lengths. The oscillation probabilities for antineutrinos P αβ are then computed as in (2.18).

Expansion of Oscillation
Probabilities with respect to x a = cos 2θ a and δm 2 s Although the expressions (2.21) for the oscillation probabilities P αβ , and P αβ via (2.18), are simple enough to perform numerical analysis efficiently, we can obtain further insight by keeping only the leading terms of the following two small parameter of the three neutrino model, x a ≡ cos 2θ a = 1 − sin 2 2θ a = 0.21 +0.06 −0.10 , (2.22a) δm 2 s |δm 2 a | = 0.032 ± 0.002 , (2.22b) whose numerical values are constrained from the data [57,59,60], as summarized below in (4.3). First, by expanding c a and s a in terms of x a : the oscillation probabilities P αβ (2.21) are expanded as, We can clearly identify the linear terms of x a in P eµ and P µe , which are identical, and also in P µµ . In the above expansion, we keep the terms of order x 2 a , which turn out to have significant impacts in the measurement of x a despite the smallness of x 2 a 0.05 at 90% confidence level. Furthermore, we introduce a short-hand notation, cos δ ′ ≡ 2c a s a cos δ ≈ 1 − x 2 a cos δ , sin δ ′ ≡ 2c a s a sin δ ≈ 1 − x 2 a sin δ . (2.25) in (2.24), without expanding the factor 1 − x 2 a , since all the δ-dependence in the transition probabilities (2.21) are functions of 2c a s a cos δ and 2c a s a sin δ. The uncertainty of the δ-measurement should be modulated by the factor 1/ 1 − x 2 a . We find it quite useful to express the oscillation probabilities P αβ in (2.24) and the corresponding antineutrino oscillation probabilities P αβ as, where P αβ with k = 1, · · · , 6 are the coefficients of corresponding terms linear in x a , cos δ ′ , sin δ ′ , x a cos δ ′ , x 2 a , and cos 2 δ ′ , respectively. The magnitude of these coefficients determines the experimental sensitivity of measuring the two mixing parameters, x a and δ. The coefficients P (k) αβ of (2.24) are shown in the following table.
It is clearly seen from (2.27) that P ee = P (ν e → ν e ) has no dependence on θ a and δ, all the other oscillation probabilities have terms P (1) αβ linear in x a , the coefficients of cos δ ′ are the same for P eµ and P µe , those of sin δ ′ have the same magnitude but the opposite sign between P eµ and P µe , while P µµ has no dependence on sin δ ′ . Most of these properties of the oscillation probabilities are expected from theoretical considerations, while they are made explicit in (2.26) and (2.27). The corresponding coefficients for the antineutrino oscillations, P where S ′ ij are the oscillation amplitudes in the propagation basis which are obtained from S ′ ij by reversing the sign of the matter potential a(x): The relation (2.18) between the ν andν oscillation probabilities, P αβ and P αβ , respectively, is simplified significantly in the propagation basis where the matter dependence and the δ dependence of the oscillation amplitudes are factorized. The parameter dependences of the oscillation probabilities P αβ and P αβ are further simplified significantly when we take account of the smallness of the mass squared difference δm 2 s as compared to |δm 2 a |, (2.22b). We note in the propagation basis Hamiltonian (2.8) that if we set δm 2 s ≡ δm 2 12 = 0, then the oscillation occurs only between ν ′ 1 and ν ′ 3 , and hence the transitions between ν ′ 1 and ν ′ 2 , and those between ν ′ 2 and ν ′ 3 should be suppressed, in the propagation basis. In our numerical study of the atmospheric neutrino oscillations in the energy range 2 GeV < E ν < 20 GeV, we find |S ′ 12 | < 0.15, and |S ′ 23 | < 0.06. We therefore obtain the following approximation by dropping all the terms of order (δm 2 s /δm 2 a ) 2 : (2.31c) and (2.27) is further simplified as follows, In this approximation, there are only 6 independent oscillation factors in the propagation basis, of which |S ′ 11 | 2 determines the overall rates P , and the coefficients of the cos δ ′ terms are P µµ . For the x 2 a term, its coefficient |S ′ 22 − S ′ 33 | 2 is also independent. The coefficient P     We examine the energy and zenith angle dependence of these 6 oscillation factors in Fig. 1 and Fig. 2. Shown in Fig. 1(a) and (b) are the E ν dependence of the coefficients P , respectively, for the baseline along five zenith angles, cos θ z = −1, −0.9, −0.8, −0.6 and −0.4. In each panel, the solid and dashed curves are for ν andν oscillations, respectively, shown by the thick lines for NH and by the thin lines for IH. It should be noted that the coefficient |S ′ 11 | in Fig. 1 not only determines P (0) ee , as in (2.16a), but also governs all the coefficients of x a , P µµ in the approximation (2.31). We immediately notice in Fig. 1(a) the absence of the significant oscillation in P µµ (NH), the oscillation curves for the same contributions, solid-thin and dashed-thick lines, show the vacuum-oscillation like pattern. They are consequences of the absence of the MSW resonance in these cases, as explained in Sec. A.4. Conversely, the strong oscillation pattern for P ee (NH) and P ee (IH) in Fig. 1(a) and the significant deviation from the vacuum oscillation pattern for P (0) µµ (NH) and P (0) µµ (IH) in Fig. 1(b) are both consequences of the MSW resonance at E ν ∼ 6 GeV for the earth matter density of ρ ∼ 5g/cm 2 along the baseline with cos θ z < −0.6; see Fig. 12. More generally, we find and vice versa for P (k) αβ (IH). The relative minus signs for k = 2, 3, 4, as compared to the relations (2.28) within the same hierarchy, are consequences of the extra minus sign in S ′ 13 and S ′ 22 − S ′ 33 , when both δm 2 a and a(x) reverse signs in the limit of vanishing δm 2 s . Therefore, if we observe the presence or absence of the MSW resonance effects in ν andν oscillations, we can determine the neutrino mass hierarchy. However, as shown clearly in Figs. 1(a) and (b), the oscillation probabilities of ν in one mass hierarchy are very similar to those ofν in the other mass hierarchy. Therefore, the capability of an atmospheric neutrino detector that cannot distinguish particle charges depend critically on the difference in the flux times cross section products of ν andν, as shown in Fig. 8.
µe and I(S In Fig. 2, we show the coefficients of cos δ ′ = 1 − x 2 a cos δ and sin δ ′ = 1 − x 2 a sin δ, which determines the sensitivity of the neutrino oscillation among the e and µ flavors on the CP phase δ. The real part R(S ′ 12 S ′ * 13 ) shown in Fig. 2(a) governs all the coefficients of cos δ ′ , see (2.32), whereas the imaginary part I(S ′ 12 S ′ * 13 ) in Fig. 2(b) dictates the sin δ ′ coefficients of P  (2.33) between the ν oscillation in IH and theν oscillation in NH, and vice versa, between the ν oscillation in NH and theν oscillation in IH also holds rather well, despite small phase-shifts due to oscillations in δm 2 s . In addition, we note the smallness of their magnitudes, typically at the level of 3% for cos θ z < −0.8, being terms of order δm 2 s /δm 2 a . Note that they are larger at lower energies, E ν 6 GeV. Consequently, the measurement of the CP phase δ may require sensitivity to the ν µ ↔ ν e andν µ ↔ν e oscillations at relatively low energies. Finally, in Fig. 3(a) we show the coefficients of the cross term x a cos δ ′ , which appears only in P (4) µµ and P (4) µµ , and in Fig. 3(b) the coefficients P (5) µµ and P (5) µµ of the quadratic term x 2 a . We should note that P (4) αβ is of the same order as P (2) αβ and P αβ , and also satisfy the same features described in the last paragraph. Generally speaking, the coefficient of the x a cos δ mixing term is small in magnitude as compared to those of x a , P as can be inferred from the |S ′ 11 | 2 plots of Fig. 1(a), especially at high energies of E ν 5 GeV. Therefore, we expect little δ-dependence in the x a measurement. On the other hand, the coefficients of x 2 a , P µµ and P (5) µµ shown in Fig. 3(b), are large in magnitudes and can dominate the terms linear in x a , even for x 2 a ∼ 0.04, especially for the neutrino oscillation in IH and the antineutrino oscillation in NH where 1 2 (1 − |S ′ 11 | 2 ) has very small magnitude at E ν 4 GeV, see the thin-solid and thick-dashed curves in Fig. 1(a). Consequently, their contributions can be significant in the measurement of x a .

Event Rates at a Charge-Blind Detector
We are now ready to study systematically the event rate distributions of atmospheric neutrino observation at PINGU. It cannot distinguish particle charges, but has the capability of resolving high energy µ ± tracks from e ± and/or hadronic showers [23,24,25].
Out from the four components in the atmospheric neutrino flux, namely the fluxes of electron-and muonneutrino/antineutrino, as shown in Fig. 8, the PINGU detector [27] is assumed to observe both electron-like and muon-like events [28], , (3.1b) by summing over contributions from charged-current (CC) e ± and µ ± production events. Since the fiducial volume is universal for both e ± and µ ± channels, as explained in Sec. A.1, it serves as an overall factor. For each flavor, neutrino and antineutrino contribute with the corresponding CC cross sections. It should be noted that we neglect the contributions from tau-neutrino/antineutrino since the tau-neutrino flux is very small [51] and also because the charged-current τ ± production events followed by their pure-leptonic decays contribute mainly to events with low observable energies, which may not contribute much due to the smaller fiducial volume. In addition, the tau neutrino cross section is small compared to the CC electron and muon neutrino cross sections [52]. Contributions from the charged-current τ ± production events as well as neutral-current events will be studied elsewhere.
Since the number of signal events depends on the oscillation probabilities linearly, which have been decomposed into six terms in (2.26), the event rates can also be decomposed accordingly, By combining with the explicit expressions of the decomposed oscillation probabilities in (2.27), the coefficients for electron-like event number rates are, For brevity, the arguments E ν and cos θ z have been omitted. Note that there is no term with x a cos δ ′ , x 2 a or cos 2 δ ′ dependence for electron-like events since P For the muon-like events, we find, Note that for muon-like events, the crossing term x a cos δ ′ has nonvanishing coefficient N µ . Consequently, with muon-like events included, we should observe some correlation between the measurements of the atmospheric mixing angle x a and the CP phase δ, as will be described in Sec. 4. In addition, a nonzero N (6) µ , which is one further order of magnitude smaller than N (4) µ , is kept according to (2.27) just to show its magnitude. By combining everything together, the atmospheric neutrino flux, cross section and effective fiducial volume of Fig. 8 in Sec. A.1, and the oscillation probabilities discussed in Sec. 2, the energy and the zenith angle dependences of the coefficients for muon-and electron-like event rates are shown in Fig. 4, Fig. 5, Fig. 6, and Fig. 7. Note the different scales of the plots, which are adjusted to show the structure of the coefficients.   Fig. 4(a) show significant oscillatory behavior for both NH (thick-solid lines) and IH (thick-dashed lines). However, the huge hierarchy dependences in the ν µ → ν µ oscillation (MSW [37] resonance only for NH) and in theν µ →ν µ oscillation (MSW resonance only for IH) as shown in Fig. 1(b) diminish significantly because of the cancellation between the ν µ andν µ contributions. Because the flux times cross section for ν µ is a factor of about three larger than that forν µ as shown in Fig. 8(b), the hierarchy dependence of the ν µ → ν µ oscillation survives, resulting in the smaller rate for IH at the MSW resonant energy of ∼ 6 GeV at cos θ z − 0.8. Especially at cos θ z −0.9, shown in the top two panels of Fig. 4(a), the nearly maximal resonant oscillation of ν µ → ν µ for NH at E ν ∼ 4 GeV shown by the thick-red curves in the top two panels of Fig. 1(b) gives rise to the significant difference in the muon-like event rate in the 3 ∼ 5 GeV region due to the so-called parametric resonance [44,45]. Because of the large event numbers, there is a possibility that these differences can be identified in experiments and that the neutrino mass hierarchy is determined. We note, however, that the finite energy and angular resolution of real experiments may make it difficult to identify differences which depend strongly on the energy, such as those in the oscillation phase observed at E ν 4 GeV at all cos θ z . On the other hand, the hierarchy dependence of the electron-like event rate N (0) e , shown also by thick-red lines in Fig. 4(b), has little dependence on cos θ z and does not oscillate in E ν . Although both the overall rate and the difference is small, the event is consistently higher for NH than IH in the broad energy range of 2 ∼ 10 GeV, reflecting the MSW and parametric enhancements of the ν µ → ν e oscillation, P , for NH; see (2.32) and Fig. 1(a). Such moderate E ν and cos θ z dependences of the electron-like event rate on the mass hierarchy may allow actual experiments to identify the difference.
Let us now examine the coefficients N are both proportional to 1 2 (1 − |S ′ 11 | 2 ) in the approximation of (2.32), and hence the energy-angular dependences are mild especially at high energy region of 4 ∼ 10 GeV, just like the electron-like event rate N (0) e . This will help experiments to measure x a . Because of the positive sign of N (1) µ , the mass hierarchy determination by using only the muon-like events should be easier for x a < 0 (sin 2 θ a < 0.5) than for x a > 0 (sin 2 θ a > 0.5). The trend can be reversed when the electron-like events are also included because N (1) e has negative sign and has relatively larger magnitude. Likewise, x a will be measured more accurately for NH than for IH when only the muon-like events are studied, whereas the measurement for IH can be significantly improved by including the electron-like events in the analysis. The dependence on the CP phase δ is shown in Fig. 5 and Fig. 6, respectively, for NH and IH. In both figures, the N α which are measured in unit of 1000/GeV as shown in Fig. 4. If we restrict our attention to the higher energy region of E ν > 4 GeV which is less sensitive to the experimental energy-angular smearing effects, the electron-like events in Fig. 5(b) have higher sensitivity to both cos δ [red-solid lines] and sin δ [blue-dashed lines] than the muon-like events in Fig. 5(a). All the four coefficients N α for α = µ and α = e have larger magnitudes at lower energies, E ν 3 GeV, although they oscillate rapidly with E ν . The expected energy resolution of the PINGU detector may smear out those rapid oscillation. However, in a certain cos θ z region the coefficients tend to have a definite sign which may survive after the energy smearing. For instance, let us examine the sin δ ′ measurement by using the coefficients N tends to be negative in the whole region. They tend to oscillate about zero at cos θ z = −1.0, and the sign reverses at cos θ z = −0.6. Therefore, if the angular resolution of experiments can resolve cos θ z = −0.9 (θ z ∼ 154 • ) from cos θ z = −0.6 (θ z ∼ 127 • ), then it might be possible to measure sin δ ′ by using the total number of events including the low energy region. The same applies for the cos δ ′ measurements, for which the coefficient N e . Although a quantitative study with realistic event simulation is beyond the scope of the present paper, probability of using the low energy data for measuring δ may worth serious studies.
In case of IH, the coefficients N (2) α and N α behave as in Fig. 6 (a) and (b), respectively, for muon-like (α = µ) and electron-like (α = e) events. Although the general trend of the oscillation patterns looks very similar between  for NH and Fig. 6 for IH, the magnitude of the coefficients at high energies (E ν > 4 GeV) are significantly smaller for IH than those for NH. We should therefore expect that the δ measurement is more difficult for IH than for NH, when higher energy data are used. On the other hand, the cos θ z dependence of the sign and the magnitude of the coefficients integrated over the low energy region below 4 GeV are similar to those of NH for cos θ z −0.6. Dedicated studies with realistic energy and angular resolution may reveal the possibility of measuring δ even in the IH case.
µ is similar to those of N (2) α and N α for cos δ and sin δ in Fig. 5 and Fig. 6, as expected from an order of magnitude estimates in (2.30). Also, its energy dependence is similar to that of cos δ coefficient N (2) µ . We therefore expect that the uncertainty of the cos δ measurement depends on the sign of x a when the muon-like events are used in the analysis.
The magnitude of the coefficient N (5) µ of x 2 a in Fig. 7(b) is of the same order of magnitude as that of N (1) α for the linear term in Fig. 4. Note that the maximum value of the coefficient N (1) µ of x a in Fig. 4(a) is around 140, 250, 690, and 570 for cos θ z = −1, −0.9, −0.8, and −0.6, respectively, for NH, whereas those of the IH are smaller. On the other hand, the maximal values of the x 2 a coefficient N Fig. 7(b) are 1800, 1870, 1264, and 1050, for the same region. Around the peak of N (5) µ , its magnitude is typically more than one order of magnitude larger than that of N

A Simple χ 2 Analysis
In this section, we examine the potential sensitivity of the PINGU experiment to the three unknown neutrino oscillation parameters at the neutrino event level. In other words, we assume that both the neutrino energy (E ν ) and its momentum direction (cos θ z ) are measured exactly for each event, and by ignoring uncertainties in the neutrino flux, cross sections and effective fiducial volume, the probability to misidentify µ ± and e ± events, as well as backgrounds from τ -decay and neutral current events. Although these assumptions are far from reality, the results are still useful in identifying the maximum information hidden in the data, motivating and directing studies with full detector simulations. Especially, it can demonstrate that the decomposition method in the propagation basis is extremely powerful to reveal the hidden patterns behind the neutrino oscillogram.

χ 2 Function
We introduce a conventional χ 2 technique to investigate experimental sensitivities on the neutrino mass hierarchy, the atmospheric mixing angle θ a and its octant, as well as the CP phase δ. Note that the result of this method, χ 2 min corresponds to the so-called "average experiment" [62] or "Asimov data set" [63]. The uncertainty from statistical fluctuation can be easily estimated to be ∆(χ 2 min ) ≈ 2 ∆(χ 2 min ) [64]. It not only applies to the case of discrete variables such as the neutrino mass hierarchy and the octant, but actually applies generally as long as the binned event number is large enough such that statistical fluctuation can be approximated by Gaussian distribution. Based on these two key parameters, χ 2 min and its variation ∆(χ 2 min , statistical interpretation can be made. The χ 2 function receives contributions from the statistical uncertainty of the event numbers as functions of the neutrino energy E ν and its momentum direction cos θ z , as well as external constraint on neutrino oscillation parameters which has been denoted as χ 2 para , We generate events with the mean values of the parameters in (4.3) for one of the mass hierarchies, except for x a = cos 2 θ a − sin 2 θ a , for which we examine three input values ±0.2 and 0, which are consistent with the present constraints. As for the CP phase δ, we examine four cases, 0, π, and ±π/2. We then use MINUIT [65] to find the minimum of the χ 2 function by varying all the six parameters, δm 2 a , δm 2 s , sin 2 2θ s , sin 2 2θ r , x a , and δ.
Not all the information in atmospheric neutrino mixing pattern can be retrieved after reconstructing the events. The mixing pattern with low energy and small | cos θ z | will be lost due to smearing and detector resolution. To see how this would affect the result, we will apply simple event cuts on neutrino energy E ν and the zenith angle θ z when presenting the results below.
The dependence on the true values of the neutrino oscillation parameters is consistent with those in the previous studies [13,16,33,34], except for the x a -dependence of the hierarchy sensitivity as explained in Sec. 4.2.

The Mass Hierarchy
As shown in Fig. 1, the neutrino mass hierarchy can be determined by observing the MSW resonances due to the Earth matter effect which occurs only for NH in neutrino oscillations and for IH in antineutrino oscillations. Although the differences are partially cancelled for a detector like PINGU which is incapable of distinguishing neutrino from antineutrino, it is still possible to determine the neutrino mass hierarchy with atmospheric neutrino oscillation because of incomplete cancellation. The hierarchy can be defined as, where the χ 2 minimum is obtained by setting the neutrino mass hierarchy to be normal (NH) or inverted (IH).  Table 1: The dependence of hierarchy sensitivity ∆χ 2 MH on the input values of the atmospheric angle's deviation,x a , and the CP phase,δ, with 1-year running of PINGU. The cases of both NH and IH, muon-and electron-like events have been considered with event cut E ν > 6 GeV and cos θ z < −0.4.
Since the CP phase δ and the atmospheric mixing angle θ a have not been pinned down yet, their values would affect the distinguishability of the mass hierarchy. The dependence of the hierarchy sensitivity ∆χ 2 MH on the input values of δ and the parameter x a = cos 2 θ a − sin 2 θ a is summarized in Table 1. Four typical cases ofδ = n 2 π with n = 0, 1, 2, 3 respectively and three possibilities forx a = ±0.2, 0.0 have been shown for both NH in the left and IH in the right. Since muon-like events are easier to be measured, we first show the results with only muon-like events in the upper part and then also include the electron-like events in the lower part. The results in Table 1 are obtained with the event cuts E ν > 6 GeV and cos θ z < −0.4 which will be detailed below. Even with this limited parameter space, the hierarchy sensitivity ∆χ 2 MH is sizable, being larger than 35 for NH and 48 for IH for all cases of input values for δ and θ a , under the assumption of a perfect detector. The smallest value of ∆χ 2 MH in each block has been marked as black bold numbers.
With only muon-like events, the hierarchy sensitivity is comparable between NH and IH. This is because in the considered energy range, E µ > 6 GeV, the zeroth order event rates N (0) µ alternate without preference for neither NH or IH. If the energy cut is lowered down to E ν > 4 GeV, the event rate for IH dominates over the one with NH, rendering larger hierarchy sensitivity for NH. The situation does not change much when the electron-like events are included. This can be seen by comparing the results withx a = 0 where complication for the x a terms are absent.
For all four blocks, the hierarchy sensitivity is larger withδ = 0 • than the result withδ = 180 • . This is because the in the considered energy range E ν > 6 GeV, the coefficient N (2) µ of cos δ ′ dominates. In addition, it is mainly negative for NH and positive for IH, as shown in Fig. 5. With a positive cos δ ′ , the difference between NH and IH is enlarged when N (2) µ cos δ ′ is added to the zeroth order N (0) µ . The dependence on x a is a little more complicated due to the presence of both the linear and quadratic terms. Let us first compare the results withx a = ±0.2 which have the same contribution from the quadratic term. For NH, the hierarchy sensitivity is larger for negativex a . This is because in most part of the considered parameter region, especially in the regions cos θ z −0.7 and cos θ z −0.9, negative x a makes the difference between NH and IH larger. This trend remains when the energy cut is lowered down to E ν > 4 GeV and is enhanced when electron-like events are included. The contribution from the quadratic term N (5) µ is always positive and its hierarchy dependence is the opposite to that of N (0) µ , leading to a negative effect. If only linear term of x a is present, the dependence is monotonically decreasing and is lifted by the quadratic term whenx a vanishes. For IH, the linear term coefficient N (1) µ is much smaller reducing the difference betweenx a = ±0.2 and making the increase atx a = 0 prominent. Note that the dependence onx a is different from the results in [13,16,33,34] due to detector responses.
In the upper-left block of Table 1 for NH with only muon-like events, the dependence on δ is smaller than that on x a . For each row with fixed input value of δ, the variation is around 9 ∼ 12, while for each column with fixed input value of x a , it is around 5 ∼ 6. This property applies for all the other three blocks. In the upper-right block for IH with only muon-like events, the variant in rows is around 9 ∼ 10, but the variant in columns is much smaller being around 2 ∼ 3. Such trends are expected since the magnitude of the observable coefficients of x a , namely N (1) µ , is much larger than those of cos δ ′ and sin δ ′ , N (2) µ and N  Fig. 4. The x a -dependence of χ 2 MH is consistently larger for NH than for IH, because the coefficient N (1) α of x a is larger for NH than for IH, as shown in Fig. 4. It is remarkable that when electron-like events are included in the analysis, the hierarchy distinguishing power increases significantly for x a = −0.2, but not much for x a = +0.2. This is because of the negative sign of N (1) e , shown in Fig. 4(b), which enlarges the hierarchy dependence of the event rate for negative x a .
Although the absolute magnitude of ∆χ 2 MH in Table 1 for a perfect detector without systematic uncertainty do not have much significance, the relative importance of electron-like events and possible impacts of the x a value in the hierarchy determination may want further studies. Note that the neutrino mass hierarchy can be resolved no matter what true values of the atmospheric angle and the CP phase can be, in contrast to the CP-hierarchy and octant-hierarchy degeneracies from which accelerator based neutrino experiments suffer [66].  Table 2: The dependence of the hierarchy sensitivity ∆χ 2 MH on event selection cuts for 1-year running of PINGU. The cases of NH and IH, µ and e+µ-like events respectively have been considered. The input values of x a and δ are taken as those in Table 1 that produces the smallest value of ∆χ 2 MH . Namely, (x a , δ) = (−0.2, 180 • ) for µ-like events with IH, and (x a , δ) = (+0.2, 180 • ) for the others.
Those events at low energy and/or small | cos θ z | will be largely smeared out, and oscillating features may be averaged out. This is because the energy smearing mainly comes from neutrino scattering, which is expected to scale as a linear function δE ∝ E, and statistical fluctuation, which scales as δE ∝ √ E. On the other hand, the neutrino oscillation period shrinks quickly at low energy, approximately as a quadratic function ∆E ∝ E 2 . No oscillation signal can be expected to survive below some energy threshold. For angular resolution, δθ z , it should be roughly a constant in the neutrino frame. When converted to the parameter in phase space, the resolution, δ(cos θ z ) = sin θ z δθ z , is much larger for horizontal events, sin θ z ≈ 1. Hence, these regions may not contribute much in more realistic studies, and can be omitted by simply applying event selection cuts E ν > E cut ν and cos θ z < cos θ cut z . The dependences on E cut ν and cos θ cut z are shown in Table 2 where the input values of δ and x a are chosen corresponding to the smallest value of ∆χ 2 MH in each block of Table 1 respectively. We observe that the results depend strongly on the range of E ν and cos θ z , which can be effectively analyzed in real experiments. For instance, when E cut ν is raised from 6 GeV to 8 GeV, ∆χ 2 MH drops by nearly a factor of 3 ∼ 7, while for E ν > 8 GeV, changing of the zenith angle coverage from cos θ z < −0.4 to cos θ z < −0.6 reduces ∆χ 2 MH by further factor of 2. It is therefore very important to have low energy threshold of the detector and to study the smearing effects in detail.

The Atmospheric Angle and Its Octant
Once the mass hierarchy is determined, χ 2 fit can be performed with the correct mass hierarchy as an input, and the same χ 2 function (4.1) can be used to measure the atmospheric mixing angle and its octant discussed here, as well as the CP phase δ which will be discussed in Sec. 4.4. In order to make the global minimum to be at the input value of x a =x a = x input a , we modify the last term in χ 2 para (4.2) as follows, which keeps the uncertainty in sin 2 2θ a the same as in (4.3) while shifting the mean value to the input, sin 2 2θ a = 1−x 2 a . This avoids small but inessential dependence of the x a measurement on the input values of x a in the region |x a | < 0.2.

the Atmospheric Angle
We compute the minimum of the χ 2 function as a function of the parameter x a , by varying all the other 5 parameters with MINUIT2 [65], to obtain the region, χ 2 min (x a ) ≤ 1, whose half-width is defined as the expected uncertainty of the x a measurement, ∆(x a ) or ∆(x 2 a ). The results for NH and IH are listed in the left and right panels of Table 3, respectively. For each hierarchy, four input valuesδ = 0 • , 90 • , 180 • , 270 • of the CP phase δ and three input values x a = 0.0, ±0.2 for the atmospheric angle θ a have been tested. Both the uncertainty of x a , ∆(x a ), and that of x 2 a , ∆(x 2 a ), are given in Table 3.  Table 3: The dependence of the uncertainty of x a and x 2 a , corresponding to χ 2 min (x a ) ≤ 1 and χ 2 min (x 2 a ) ≤ 1, respectively, on the input valuesx a andδ after 1-year running of PINGU. The case of both muon-and electron-like events for both NH and IH has been considered with event cuts E ν > 6 GeV and cos θ z < −0.4.
Note that the uncertainty ∆(x a ) is larger if the atmospheric mixing angle is maximal, or whenx a = 0 (sin 2 2θ a = 1). This is a consequence of relatively large coefficient N (5) µ of x 2 a , shown in Fig. 7(b), as compared to the coefficients N  Fig. 4 (a) and (b), respectively. The x a -dependence of the muon-like events are expected to be, ignoring the other small terms. The variation of the number of events at x a =x a is then, Since the N for |x a | = 0.2, and it vanishes forx a = 0, the combined effective coefficient of δx a has larger magnitude for nonzerox a . This explains the reduced uncertainty ∆(x a ) for x a = ±0.2. The slight difference between the two mirror casesx a = −0.2 andx a = 0.2 comes from the first term. Since N (1) µ is positive, cancellation happens in the combined effective coefficient whenx a is negative, leading to systematically larger ∆(x a ). There is some slight dependence on the input valueδ of the CP phase, but not sizable.
To make a direct comparison with the external constraint (4.5), the resolution ∆(x 2 a ) is also shown in Table 3. Since sin 2 2θ a ≡ 1 − x 2 a , the uncertainty ∆(sin 2 2θ a ) is exactly ∆(x 2 a ). For x a ≈ |0.2|, it can be roughly estimated as ∆(x 2 a ) ≈ (2x a )∆(x a ) ≈ 0.4∆(x a ) which is typically a factor of 5 ∼ 6 smaller than the uncertainty of 0.03 given in (4.3) and (4.5). If the true value of x a vanishes, its uncertainty should be estimated as ∆(x 2 a ) ≈ [∆(x a )] 2 which is roughly 0.001 for NH, which is smaller than the current uncertainty by a factor of 30. The uncertainties ∆(x a ) and ∆(x 2 a ) are slightly larger for IH due to the smaller coefficients N (1) α as shown in Fig. 4. Summing up, a neutrino telescope like PINGU has a potential to resolve the octant degeneracy of the atmospheric mixing angle θ a , and to reduce the uncertainty of sin 2 2θ a by a factor of 5 to 30 within one year of running. Although our simulation does not take account of energy and zenith angle resolutions, we expect this high potential to be confirmed in more realistic simulations because the coefficients N (1) µ and N (1) e do not oscillate much with E ν or cos θ z .

Octant Sensitivity
If the linear terms N  Table 4 Table 4: The dependence of the octant sensitivity on the input valuesx a andδ after 1-year running of PINGU. The case of both muon-and electron-like events for both NH and IH has been considered with event cuts E ν > 6 GeV and cos θ z < −0.4.
We can see that the octant sensitivity forx a = ±0.1 is much smaller than the one forx a = ±0.2, as expected. With smaller distance between the two mirrors, the difference due to the linear term is much smaller. And we can estimate the significance to scale roughly as ∆χ 2 octant ∝ x 2 a , according to (4.1), as verified by the results in Table 4. Between the two mirrors, the octant sensitivity is always larger for positive x a due to the same sign between N (1) µ and N (5) µ . This makes the effective linear term coefficient in (4.7) larger with positivex a , hence enhances the octant sensitivity. This trend is further enhanced by including the electron-like events which has only a negative linear term coefficient N (1) e but not quadratic term. With positivex a , the event rates become smaller, explaining the further enhancement.
The dependence on the neutrino mass hierarchy is much easier to be understood. For NH, the octant sensitivity is much larger than that for IH, because the linear term coefficients N (1) µ and N (1) e have much larger magnitude for NH. There is small dependence on the CP phase. It comes from the cross term x a cos δ ′ whose coefficient N (4) µ is mainly positive for NH and negative for IH in the considered energy range E ν > 6 GeV, especially around cos θ z ≈ −0.6, as shown in Fig. 7(a). The effective linear term in (4.7) becomes, With nonzerox a , especially whenx a ≈ ±0.2, the quadratic term coefficient N dominates. Cancellation between cosδ ′ N (4) µ and 2x a N (5) µ happens if they have opposite signs, leading to smaller octant sensitivity as shown in Table 4. For NH, the octant sensitivity is larger for cosδ 0 whenx a is negative and cosδ 0 whenx a is positive. It is the opposite for IH. This trend remains when the electron-like events are also included in the analysis, since the cross term coefficient N (4) e for the electron-like event rates is zero, although the sensitivity can be enhanced a lot. Note that the cross term does not have sizable effect on the uncertainty of measuring the atmospheric mixing angle but manifests itself in the octant sensitivity.
Since there is not so much oscillation behavior in N µ , and N µ , these trends in the octant sensitivity ∆χ 2 Octant , obtained with neutrino events without full simulation, should survive the smearing effects from the neutrino scattering and detector resolutions.

Uncertainty of the CP Phase
In this section, we will show the capability of PINGU in measuring the CP phase δ in terms of χ 2 min (δ). The dependence of χ 2 min (δ) on δ comes from fixing δ and fitting the other five parameters to find the minimum. The uncertainty ∆(δ) are then obtained according to the condition χ 2 min (δ) < 1 for each case of the input values,x a = 0, ±0.2, and δ = 0, ±π/2, π, as well as both NH and IH. Since the δ-dependent terms have tiny coefficients, N α , as shown in Fig. 5 and Fig. 6, it is very challenging to determine the CP phase and can only be possible with a much longer time. We check the results after 10-years running of PINGU with both muon-and electron-like events.
As explained in Sec. 2 and Sec. 3, the event rates depend on the CP phase δ only through the terms proportional to cos δ and sin δ, since the coefficient N (6) µ of cos 2 δ in the expansion (3.2) is negligibly small. Therefore, the δ-dependence of χ 2 min (δ) can be approximated as a quadratic function of cos δ and sin δ, where the covariance matrix V is a 2 × 2 real symmetric matrix, and ∆(δ) is the uncertainty on the CP phase. For various inputs, the results from exact χ 2 minimization have been shown in Table 5.
∆(δ) NH IH x a = −0.2x a = 0.0x a = +0.2x a = −0.2x a = 0.0x a = +0.2 Table 5: Dependence of the uncertainty of measuring the CP phase δ on the input valuesx a andδ, after 10-years running of PINGU, for the NH (left 3 columns) and for the IH (right 3 columns). Both muon-and electron-like events in the region of cos θ z < −0.4 are used for E ν > 6 GeV (upper 4 rows) and for E ν > 4 GeV (lower 4 rows).
Since the coefficients N (2) α and N (3) α are much smaller for IH, as shown in Fig. 5 and Fig. 6, the resultant ∆(δ) in Table 5 is expected to be larger. Besides, there is very slight dependence on the input values of the atmospheric mixing angle and the CP phase, in contrast to the case at accelerator based neutrino experiments [66]. Since the δ-dependence of the event rates occurs only through terms of cos δ and sin δ, the χ 2 min (δ) function is a periodic function of δ. It vanishes at δ =δ and has only one period of oscillation in the range of [0, 2π] with peak around |δ −δ| ≈ π. With energy cut E ν > 6 GeV, the maximum value of χ 2 min (δ) is around 25 ∼ 40 for NH, and just 6 ∼ 8 for IH, depending on the input values ofx a andδ. The value increases to around 12 ∼ 16 for IH if we include events down to E ν = 4 GeV. The function χ 2 min (δ) can be approximated by a quadratic function (4.10) of δ −δ up to more than 3σ for all the results quoted in Table 5 with the only exception of E ν > 6 GeV for IH where it holds up only to around 2.5σ.
There is some slightx a -dependence of the uncertainty ∆(δ), which can then be read off directly by expressing the covariance matrix in terms of the decomposition coefficients, N (2) α +x a N (4) α and N (3) α , for an input value of x a =x a as follows: where we neglect the small δ-dependent terms in the denominator of (4.11b). Since N share the same sign in the considered energy range of E ν > 6 GeV, especially around cos θ z ≈ −0.8 ∼ −0.6, as shown in Fig. 5,  Fig. 6, and Fig. 7, the combination N (2) µ + N (4) µxa has larger magnitude whenx a is positive. This explains the reduced uncertainty atx a = +0.2. Note that there is an overall factor 1 −x 2 a due to the modulation cos δ ′ = 1 − x 2 a cos δ and sin δ ′ = 1 − x 2 a sin δ. It can increase the CP phase uncertainty ∆(δ) if the atmospheric mixing angle is not maximal, namely 1 − x 2 a < 1. Theδ dependence of ∆(δ) can be expressed as, for the integration region of 6 GeV < E ν < 20 GeV and cos θ z < −0.4 with 10 years running of PINGU. From the results of (4.13) and (4.14), we find that the uncertainties ∆(sin δ) and ∆(cos δ) are rather large when only the muon-like events of E ν > 6 GeV are used in the analysis. This is especially the case of IH, for which the uncertainties are larger than unity even after 10-years of running. In other words, the uncertainty ∆(δ) can be rather broad. Therefore, the electron-like events are essential to measure the CP phase.
Generally speaking, ∆(cos δ) is slightly smaller than ∆(sin δ) due to larger magnitude of N (2) α , resulting in slightly smaller uncertainties ∆(δ) atδ = ±90 • than those atδ = 0 • , 180 • . In addition, the correlation turns out to be small |ρ| ∼ 0.2 and negative for all the cases. From the combined results in (4.13) and (4.14), we can estimate the smallest and the largest uncertainty of ∆(δ) as a function ofδ, for IH. All the results quoted in Table 5 lie within the above range. The largest differences are found to be about a few degrees for the IH case with E ν > 6 GeV.

Conclusions and Outlook
In this work, we develop a general decomposition formalism in the propagation basis, which is extremely useful for the phenomenological study of neutrino oscillation. It can analytically separate the contributions of the three unknown parameters, namely, the neutrino mass hierarchy, the atmospheric mixing angle, and the CP phase. In this way, the pattern behind χ 2 minimization can be revealed clearly, especially for atmospheric neutrino which experiences very complicated Earth matter profile. Hence, it can serve as a complementary tool to the neutrino oscillogram. The latter is designed for the overall pattern, especially the resonance behaviors, in the atmospheric neutrino oscillations, while our decomposition method can unveil more hidden structures behind the oscillogram. In addition, the decomposition method can apply generally to any type of neutrino oscillation experiment.
To illustrate the powerfulness of this decomposition formalism, we study in detail the ability of PINGU in determining the neutrino mass hierarchy, the atmospheric angle θ a and its octant, as well as the CP phase δ by measuring the oscillation pattern of atmospheric neutrinos. Both muon-and electron-like events have been considered. Our results suggest that PINGU has the potential to determine the mass hierarchy and the octant of the atmospheric mixing angle θ a within one year of operation if the neutrino energy and the zenith angle can be measured accurately in the region E ν = 6 ∼ 20 GeV and cos θ z < −0.4. The uncertainty of measuring the value of θ a can be reduced by a factor of 5 ∼ 30, while the determination of the CP phase δ is significantly more challenging. The dependence on the input values of the neutrino mass hierarchy, the atmospheric mixing angle, and the CP phase can be fully understood in our decomposition formalism. Our findings merit a serious investigation of the physics potential of PINGU with realistic detector response which we expect to be underway within the IceCube/PINGU Collaboration.

Acknowledgements
We thank Naotoshi Okamura for stimulating discussions in the early stage of our investigation. SFG is grateful to Hong-Jian He for kind support and Center for High Energy Physics of Tsinghua University (TUHEP), where part of this work is done, for hospitality. JSPS has provided SFG a generous postdoc fellowship, which is deeply appreciated, to continue the study at KEK. This work is supported in part by Grant-in-Aid for Scientific research (No. 25400287) from JSPS. CR is grateful for support from the Center for Cosmology and AstroParticle Physics (CCAPP) at The Ohio State University, where part of this work was performed. This work is supported in part by the Basic Science Research Program through the National Research Foundation of Korea funded by the Ministry of Education, Science and Technology (2013R1A1A1007068).

A. Atmospheric Neutrino Oscillation
In this appendix, we introduce the input and method that we used to evaluate the atmospheric neutrino oscillation in this study. They include the atmospheric neutrino fluxes, interaction cross sections, and an energy-dependent effective detector volume, the Earth matter density profile and the numerical procedure.

A.1. Atmospheric Neutrino Flux, Cross Sections and Effective Volume
The atmospheric neutrino flux depends on many factors. First, it varies with the geographic location, mainly due to the earth magnetic field at the source regions. In addition, it dependents on the neutrino momentum direction, the zenith and the azimuth angles. Seasonal effects can also modulate the neutrino flux. For our study, we use an annual and azimuth angle averaged neutrino flux computed for the South Pole [51]. The earth magnetic field effect, which introduces the largest modification on neutrino fluxes, depending on the position, are mostly relevant for neutrino energies below the one considered in this study; therefore our results can be easily transferred to detectors at other geographic locations. With these factors taken into consideration, the neutrino flux is a function of neutrino energy and the zenith angle.
The energy dependence of the atmospheric neutrino flux is shown in Fig. 8(a). It can be seen that ν µ andν µ fluxes dominate over the ν e andν e fluxes. They drop very quickly with increasing energy, decreasing by four orders of magnitude from E ν = 1 GeV to E ν = 20 GeV. For both flavors, the antineutrino flux is slightly smaller than the neutrino flux due to an asymmetry between the π + and π − production spectra in cosmic ray air-showers [51].
The atmospheric neutrino flux will be modulated by the neutrino interaction cross section with nuclei. In this study, we use the charged current (CC) cross sections generated by NEUGEN3 [53]. Since the cross sections increase linearly as function of neutrino energy E ν , higher energy event rates can be enhanced, by roughly an order of magnitude. In addition, the neutrino cross section is a factor of 2 − 3 larger than that of the antineutrino. The difference between the neutrino and antineutrino rates becomes even more significant when the atmospheric neutrino flux in Fig. 8(a) is multiplied with cross sections, as shown in Fig. 8(b). This difference between the neutrino and antineutrino event rates is critical for the mass hierarchy determination with a detector incapable of telling leptons (µ and e) from antileptons (μ andē), as will be made explicit in Sec. A.4. The PINGU effective fiducial volume used for our study assume a geometry of 20 additional strings within a radius of 75 m inside the DeepCore volume [54] with an inter-string spacing of approximately 26 m. It was required that at least 20 optical sensors would register a Cherenkov photon. At this level, events are assumed to be reconstructable. At the relevant energies in this study, the effective fiducial volume for electron and muon neutrinos is approximately equal [38]. Under this assumption, the fiducial volume is universal and will not affect the relative amount of electron and muon neutrinos/antineutrino event rates. The fiducial volume increases with energy, further increasing the event rates at higher energies. However, due to the much larger neutrino flux at lower energies, the highest event rates are expected to be detected at lower energies as can be seen in Fig. 8(c).  Fig. 9: Comparison of the effective fiducial volume reported by the PINGU Collaboration [54], and the one adopted in the work of Akhmedov et. al. [30,32].
In Fig. 9, we compare the effective fiducial volume reported by the PINGU Collaboration [54], that we adopt in this study [see Fig. 8], and the one adopted in the pioneering work of Akhmedov et. al. [30,32]. In the left panel, the absolute values of both PINGU (solid curve) and the Akhmedov et. al.'s paper (dashed curve) are shown while the ratio between the one used in Akhmedov et. al.'s work and the PINGU curve is presented in the right panel. We can see that in most of the regions, the one implemented in the Akhmedov et. al.'s paper is larger than the one from the PINGU Collaboration by a factor of around 3. It is only in the low energy end that this ratio decreases.

A.2. Earth matter profile
We use the preliminary Earth reference model (PREM) [48] to describe the Earth's matter density distribution. It represents the standard framework of interpreting the seismological data to determine the Earth matter density, assuming a spherically symmetric Earth without describing the chemical composition. The Earth radius varies from 6353 km to 6384 km due to rotational flattening, which distorts Earths shape from that of an ideal sphere. The variation is relatively small (0.5%) and hence it is a reasonable assumption to use a spherically symmetric Earth with radius R = 6371 km. As the matter potential felt by neutrino depends on the electron density, we need to make assumptions about the chemical composition.
The matter potential felt by neutrino is proportional to the electron density, while V (x) receives an extra minus sign for antineutrino. The parameter G F is the Fermi constant and N e (x) is the electron number density as a function of position, x ≡ d/L(θ z ) with d being the distance traversed by neutrino and L(θ z ) denoting the path length corresponding to the neutrino zenith angle θ z . The electron number density N e depends on the chemical composition of the Earth, which cannot be measured directly. It is approximated as a linear function of the matter density ρ as N e = Y ρ/m p , where m p is the nucleon mass. The coefficient Y is the ratio between electron and nucleon number densities, given by n p /(n p + n n ) where n p is the number density for proton and n n the number density for neutron. Its value can vary significantly from light elements Y ≃ 0.5 to heavier elements, which increasingly have more neutrons per protons (for example Y = 0.466 for F e). The core is expected to be dominated by iron, which takes an 85% share, rendering a smaller Y core = 0.468. For mantle, Oxygen (Y = 0.5), M agnesium (Y = 0.494), and Silicon (Y = 0.498) take 44%, 23%, and 21% of the total weight respectively, leading to a larger Y mantle = 0.497 [50]. The overall uncertainty in the averaged Y core and Y mantle is expected to be small, even if element abundances carry uncertainties of about 10% [50]. This is because the uncertainty in the averaged Y 's is approximately a product of the 10% variation in the element abundances and the variation in Y for individual element which also has around 10% variation. Hence, the estimated matter potential V would receive an uncertainty around 1%. In the current study, we omit this small uncertainty since it would not affect the potential of the atmospheric neutrino experiments at a qualitative level. The Earth's radius is denoted as R while d is the path length of the neutrino traversing the Earth. For the sectioned profile, ten sublayers, from the inside to the outside, have been sliced into n = 2 2,3,3,0,1,1,0,0,0,0 equal steps along the path. See Appendix A.3 for details. Fig. 10 shows the matter density along different paths, denoted by the neutrino zenith angle θ z . For cos θ z = −1, the path runs through the Earth center, corresponding to vertically up-going neutrino, with horizontal neutrinos denoted by cos θ z = 0. The path length is related to the zenith angle as L(θ z ) = 2R| cos θ z | and is symmetric with respect to the central point. For convenience, the horizontal axis in Fig. 10 is defined as x ≡ d/L(θ z ) to make comparison between different paths. It should be noted that there is a big discontinuity at the boundary between the core and mantle. For cos θ z = −1, this boundary rests around |x − 0.5| ≈ 0.27 since the core's size (R core = 3480 km) is almost half of the Earth's radius, R = 6371 km. In addition, the core can be divided into the inner core and the outer one which are separated by a boundary at r = 1221.5 km. The inner/outer core and mantle boundaries are very well known and have an uncertainty of less than 10 km [49]. For the mantle, there are eight sublayers with boundaries at r = (5701, 5771, 5971, 6151, 6346, 6356, 6368) km respectively. Of these eight sublayers, the three outermost layers have constant matter density. Depending on the zenith angle θ z , neutrinos pass these various layers sequentially.

A.3. Numerical Method for Oscillation Amplitude Matrix S ′ in the Propagation Basis
Although the dependences of the oscillation probabilities on the atmospheric angle θ a and the CP phase δ can be expressed analytically, as elaborated in Sec. 2, the other oscillation parameters are still entangled with the matter effect inside the oscillation amplitudes S ′ ij in the propagation basis, as explicitly shown in (2.8). Since the matter profile has a very complicated structure, depending on the neutrino zenith angle θ z , it is necessary to find an accurate and efficient numerical method to evaluate S ′ ij . As an approximation, we first replace the PREM profile within each sublayer, as shown in Fig. 10, by a constant density averaged along the path. Within each constant potential, the oscillation amplitude matrix A i for threeneutrino oscillation can be evaluated exactly [67]. The full amplitude matrix is a sequential matrix product of these individual ones. In this way, we can keep the discontinuity between sublayers, especially the periodic mantle-coremantle structure, which is important for parametric resonances. But the slowly varying behavior within each sublayer is averaged out. As an illustration, the resultant oscillation probabilities of the ν e → ν e channel for neutrino zenith angle cos θ z = −1.0, −0.9, −0.8, respectively, are shown in Fig. 11 with thin lines, in contrast to the exact solutions with thick lines. We can observe that this simple approximation shows very good agreements with the exact solution, especially that the peaks and troughs appear at almost exactly the same energies. In other words, the resonance features are maintained. But their amplitude can differ up to 10% which cannot be ignored. A more precise method is needed for a precision analysis such as χ 2 minimization to obtain the physics potential of atmospheric neutrino oscillation experiments.
To account for the finer structure, we further divide each sublayer into several sections, within each the matter density is approximated by the averaged value along the path. Since the matter density has different slopes within different PREM layers, the number of sections is chosen accordingly as n = 2 i with i = (2, 3, 3, 0, 1, 1, 0) for the density-varying sublayers from the inner core to the outer crust. In this way, accuracy and efficiency can be balanced leading to an optimized program. In principle, with sublayers divided into more number of sections, the result would be closer to the exact solution. This is verified for linear potential by comparing with exact solution [68] at various distances. For all the paths (cos θ z ) along the PREM matter distribution, we confirm that slicing the sublayers into finer sections gives no visible effects on the oscillation probabilities as we check by doubling the number of divisions with n = 2 i+1 . Hence, we call the solution with n = (4, 8, 8, 1, 2, 2, 1) as the exact solution in this study. The difference in the probabilities is found to be in the order of 10 −3 , which can be safely ignored. The oscillation probabilities P ee of the exact solution are shown in Fig. 11 with thick lines.

A.4. Normal Hierarchy v.s. Inverted Hierarchy
The significant discontinuity in the matter density between the Earth's mantle and the core leads to an interesting pattern in the atmospheric neutrino oscillation probabilities. Along trajectories with cos θ z < −0.84, neutrinos experience a large jump in the potential causing parametric resonance [44] which is also known as oscillation length resonance [45]. With periodic matter density profile, the oscillation probability is largely enhanced. In addition, there is MSW resonance in the wide range of cos θ z and E ν [37] when neutrino crosses the mantle region. This makes the oscillation pattern of atmospheric neutrinos very sensitive to the neutrino mass hierarchy. We briefly discuss the difference between NH and IH here. In Fig. 12, the electron neutrino (antineutrino) survival probabilities, P ee (Pēē), along the paths with neutrino zenith angle cos θ z = −1, −0.8 and −0.4 are plotted for both NH in the left and IH in the right panel. For NH, strong oscillation patterns are found in the oscillation probability P ee (thick curves) of neutrino, while the oscillation probability Pēē (thin curves) of antineutrino has much smaller variation. This significant difference is because of the MSW resonance effect which can be demonstrated in the much simpler two-neutrino oscillation case. Under this simplified circumstance, the effective mixing angle can be analytically expressed as, sin 2 θ = sin 2θ/ sin 2 2θ + (cos 2θ − 2EV /δm 2 ) 2 , where θ and δm 2 are the true mixing angle and the true mass squared difference. The resonance happens at cos 2θ = 2EV /δm 2 , leading to a maximal effective mixing angle sin 2θ = 1. With a full treatment of three-neutrino oscillation, the probability resonates around E = cos 2θ ij δm 2 ij /2V , where θ ij and δm 2 ij are the relevant true mixing angle and the true mass squared difference. For the solar mass squared difference δm 2 12 ≡ m 2 2 − m 2 1 , the resonance energy is around 100 MeV, with the typical matter potential in the mantle region, which is below the energy region considered in our studies. For the atmospheric mass squared difference δm 2 13 ≡ m 2 3 − m 2 1 , instead, the resonance occurs around 4-6 GeV, which is within the accessible region of the investigated atmospheric neutrino oscillation in this study. It should be noted that the MSW resonance only occurs with δm 2 13 > 0 (NH) for neutrinos (V > 0), and with δm 2 13 < 0 (IH) for antineutrinos (V < 0), since cos 2θ 13 = 1 − 2 sin 2 θ 13 > 0.
It is much simpler to determine the neutrino mass hierarchy if the detector is capable of distinguishing neutrinos from antineutrinos. With the MSW resonance around 2 ∼ 7 GeV observed in neutrinos rather than antineutrinos, then the mass hierarchy must be normal and vice versa. In other words, the existence or the absence of the MSW resonance can serve as a solid discriminator of the neutrino mass hierarchy. For a detector without the capability to distinguish neutrinos from antineutrinos, the MSW resonance could still be used to determine the mass hierarchy. This is made possible by the differences in the neutrino and antineutrino fluxes as well as their charged current cross sections discussed in Sec. A.1. The hierarchy sensitivity obtained from the residual difference between NH and IH can still be sizable, if large enough event rates are collected by a huge underground detector such as PINGU.