Atmospheric Neutrino Oscillations and New Physics

We study the robustness of the determination of the neutrino masses and mixing from the analysis of atmospheric and K2K data under the presence of different forms of phenomenologically allowed new physics in the nu_mu--nu_tau sector. We focus on vector and tensor-like new physics interactions which allow us to treat, in a model independent way, effects due to the violation of the equivalence principle, violations of the Lorentz invariance both CPT conserving and CPT violating, non-universal couplings to a torsion field and non-standard neutrino interactions with matter. We perform a global analysis of the full atmospheric data from SKI together with long baseline K2K data in the presence of nu_mu ->nu_tau transitions driven by neutrino masses and mixing together with sub-dominant effects due to these forms of new physics. We show that within the present degree of experimental precision, the extracted values of masses and mixing are robust under those effects and we derive the upper bounds on the possible strength of these new interactions in the nu_mu--nu_tau sector.


I. INTRODUCTION
Neutrino oscillations are entering a new era in which the observations from underground experiments obtained with neutrino beams provided to us by Nature -either from the Sun or from the interactions of cosmic rays in the upper atmosphere -are being confirmed by experiments using "man-made" neutrinos from accelerators and nuclear reactors [1].
For atmospheric neutrinos, Super-Kamiokande (SK) high statistics data [2,3] established beyond doubt that the observed deficit in the µ-like atmospheric events is due to the neutrinos arriving in the detector at large zenith angles, and it is best explained by ν µ oscillations. This evidence was also confirmed by other atmospheric experiments such as MACRO [4] and Soudan 2 [5].
The KEK to Kamioka long-baseline neutrino oscillation experiment (K2K) uses an accelerator-produced neutrino beam mostly consisting of ν µ with a mean energy of 1.3 GeV and a neutrino flight distance of 250 km to probe the same oscillations that were explored with atmospheric neutrinos. Their results [6] show that both the number of observed neutrino events and the observed energy spectrum are consistent with neutrino oscillations, with oscillation parameters consistent with the ones suggested by atmospheric neutrinos.
Oscillations are not the only possible mechanism for atmospheric ν µ → ν τ flavour transitions. They can also be generated by a variety of forms of nonstandard neutrino interactions or properties. In general these alternative mechanisms share a common feature: they require the existence of an interaction (other than the neutrino mass terms) that can mix neutrino flavours [7]. Among others this effect can arise due to violations of the equivalence principle (VEP) [8,9,10], non-standard neutrino interactions with matter [11], neutrino couplings to space-time torsion fields [12], violations of Lorentz invariance (VLI) [13,14] and of CPT symmetry [15,16]. From the point of view of neutrino oscillation phenomenology, the most relevant feature of these scenarios is that, in general, they imply a departure from the E −1 energy dependence of the oscillation wavelength [17].
Prior to the highest-statistics SK data, some of these scenarios could provide a good description -alternative to ∆m 2 neutrino oscillations -of the atmospheric neutrino phenomenology [18,19]. However, with more precise data, and in particular with the expansion of the energy range covered by atmospheric neutrino data due to the inclusion of the upward-going muons, these alternative scenarios became disfavoured as leading mechanism to explain the observations [20,21,22]. The results from K2K experiment further singled out oscillations as the dominant mechanism of ν µ ↔ ν τ transitions [23]. 1 The question arises, however, to what point the possible presence of these forms of new 1 Recently [24] SK collaboration has presented a reanalysis of the SK1 data in terms of the reconstructed L/E which allowed them to slightly improve the discrimination between oscillations and alternative mechanisms. Unfortunately, to reproduce such analysis for the subdominant effects discussed here is not possible outside the collaboration.
physics, even if sub-dominant, may affect the derived ranges of masses and mixing from the oscillation analysis of the atmospheric and K2K data. Or in other words, to what level our present determination of the neutrino masses and mixing is robust under the presence of phenomenologically allowed new physics effects. In this paper we address this question by performing a global analysis of the atmospheric and K2K data with ν µ → ν τ transitions driven by neutrino masses and mixing in the presence of some generic forms of new physics. In particular we consider new physics interactions which are vector-like, or tensor-like (scalar interactions cannot be distinguish from oscillations). This allow us to treat, in a model independent way, effects due to the violation of the equivalence principle, violations of the Lorentz invariance both CPT conserving and CPT violating, non-universal couplings to a torsion field and non-standard neutrino interactions in matter. In Sec. II we present the formalism adopted and the data set used. In Sec. III we show the results of our analysis. Conclusions are given in Sec. IV. The technical details of our new statistical analysis of the atmospheric data are described in the appendix.

II. FORMALISM
In what follows we consider some new physics (NP) scenarios which induce new sources of lepton flavour mixing in addition to the "standard" ∆m 2 oscillations (∆m 2 -OSC). We concentrate on flavour mixing mechanisms for which the propagation of neutrinos (+) and antineutrinos (−) is governed by the following Hamiltonian [16]: where ∆m 2 is the mass-squared difference between the two neutrino mass eigenstates, σ ± n accounts for a possible relative sign of the NP effects between neutrinos and antineutrinos and ∆δ n parametrizes the size of the NP terms. The matrices U θ and U ξn,±ηn are given by: where we have also accounted for possible non-vanishing relative phases η n . For concreteness we will focus on NP effects which are induced by tensor-like and vector-like interactions. We denote by tensor-like interactions those with n = 1 leading to a contribution to the oscillation wavelength which grows linearly with the neutrino energy. For example, Eq. (1) can describe the evolution of ν µ and ν τ 's of different masses in the presence of violation of the equivalence principle (VEP) due to non universal coupling of the neutrinos, γ 1 = γ 2 (ν 1 and ν 2 being related to ν µ and ν τ by a rotation θ G ), to the local gravitational potential φ [8,9]. 2 Phenomenology of neutrino oscillations induced or modified by VEP has been widely studied in the literature [25].
In this case For constant potential φ, this mechanism is phenomenologically equivalent to the breakdown of Lorentz invariance induced by different asymptotic values of the velocity of the neutrinos, v 1 = v 2 , with ν 1 and ν 2 being related to ν µ and ν τ by a rotation θ v [13,14]. In this case We denote by vector-like interactions those with n = 0 which induce an energy independent contribution to the oscillation wavelength. This may arise, for instance, from a non-universal coupling of the neutrinos, k 1 = k 2 (ν 1 and ν 2 being related to the ν µ and ν τ by a rotation θ Q ), to a space-time torsion field Q [12], so Violation of CPT due to Lorentz-violating effects also lead to an energy independent contribution to the oscillation wavelength [15,16] with where b i are the eigenvalues of the Lorentz violating CPT-odd operatorν α L b αβ µ γ µ ν β L and θ v is the rotation angle between the corresponding neutrino eigenstates and the flavour eigenstates [16].
In all these scenarios, if the NP strength is constant along the neutrino trajectory, the expression of P νµ→νµ takes the form [16]: where the correction to the ∆m 2 -OSC wavelength, R, and to the global mixing angle, Θ, verify R cos 2Θ = cos 2θ + n R n cos 2ξ n , R sin 2Θ = |sin 2θ + n R n sin 2ξ n e iηn | , 2 VEP for massive neutrinos due to quantum effects discussed in Ref. [10] can also be parametrized as Eq.
with R n being the ratio between ∆m 2 -induced and the NP-induced contributions to the oscillation wavelength For Pν µ→νµ the same expressions hold with the exchange σ + n → σ − n and η n → −η n . For the sake of simplicity, in what follows we concentrate in scenarios with one NP source characterized by a unique n. In this case sin 2 2Θ = 1 R 2 sin 2 2θ + R 2 n sin 2 2ξ n + 2R n sin 2θ sin 2ξ n cos η n , R = 1 + R 2 n + 2R n (cos 2θ cos 2ξ n + sin 2θ sin 2ξ n cos η n ) .
In Fig. 1 we illustrate the effect of the presence of the NP in the atmospheric neutrino events distributions for ∆m 2 -OSC plus sub-dominant CPT-even tensor-like and vector-like NP effects, for some characteristic values of the NP-parameters. In both cases R n is a growing function of E and the NP effects become relevant in the higher energy samples, in particular for upward going muons.
Furthermore, the relevant survival probabilities P νµ→νµ and Pν µ→νµ are not affected by a change in the overall sign of the Hamiltonian, as well as change in the global phase of its non-diagonal components. Therefore, we also have: • θ → θ + π/2 and ξ n → ξ n + π/2, • θ → −θ and ξ n → −ξ n , • η n → −η n . The above set of symmetries allows us to define the ranges of variation of the five parameters as follows: Thus in the general case we cover the mixing parameter space by using, for instance, 0 ≤ sin 2 θ ≤ 1 and 0 ≤ sin 2 2ξ n ≤ 1. For the case of real relative phase, η n ∈ {0, π}, one can absorb the two values of η n into the sign of ξ n . In this case we drop (e) and replace (d) by: and use instead −1 ≤ sin 2ξ n ≤ 1.
Finally we notice that the above derivation is valid for a given sign of σ + n . Keeping the convention of ∆m 2 > 0 and ∆δ n > 0 the survival probability for the opposite sign can be obtained by the exchange In addition, we also consider the special case of vector-like NP due to non-standard neutrino-matter interactions (NSI) [11,19]. In this case the effective Lagrangian describing the evolution of the ν µ -ν τ system can be written as [19,26] where N f (r) is the number density of the fermion f along the path r of the neutrinos propagating in the Earth, and ε ± and ε ′ ± parametrize the deviation from standard neutrino interactions: is the difference between the ν τ +f and the ν µ +f elastic forward scattering amplitudes. The corresponding amplitudes for antineutrinos are given by ε − and ε ′ − . For simplicity we assume that the interactions for neutrinos and antineutrinos are the same, which implies ε ′ Thus the NSI Hamiltonian contains three real parameters, which can be chosen to be ε ′ , |ε| and arg(ε). NSI and their interplay with the oscillations have also been studied in different contexts: among others, in relation to supernova neutrinos [27], to the solar neutrino problem [28], to the LSND results oscillation results [29], to neutrinoless double beta decay [30], and to present and future laboratory neutrinos [31].
Formally Eq. (16) can be seen as a special case of Eq. (1) with n = 0, σ − 0 = −σ + 0 , and Technically the main difference is that NSI only affect the evolution of neutrinos in the Earth, and their strength changes along the neutrino trajectory. Consequently the flavour transition probability cannot be simply read from Eq. (7) and its evaluation requires the numerical solution of the neutrino evolution in the Earth matter. In our calculations we use PREM [33] for the Earth's density profile and a chemical composition with proton/nucleon ratio Y p = 0.497 in the mantle and 0.468 in the core. In what follows for the sake of concreteness we set our normalization on these parameters by considering that the relevant neutrino interaction in the Earth occurs only with down-type quarks.
Concerning the data samples used in the analysis, for atmospheric neutrinos we include in our analysis all the contained events as well as the the upward-going neutrino-induced muon fluxes from the latest 1489 SK data set [2]. This amounts for a total of 55 data points. Details of our new statistical analysis, introduced here for the first time, are presented in the appendix.
For the analysis of K2K we include the data on the normalization and shape of the spectrum of single-ring µ-like events as a function of the reconstructed neutrino energy [6]. We bin the data in five 0.5 GeV bins with 0 < E rec < 2.5 GeV plus one bin containing all events above 2.5 GeV. Details of our analysis of the K2K data were presented in Ref. [32] and will not be repeated here. Let us just comment that together with statistical uncertainties we also account for the systematic uncertainties associated with the determination of the neutrino energy spectrum in the near detector, the model dependence of the amount of nQE contamination, the near/far extrapolation and the overall flux normalization.

III. RESULTS AND DISCUSSION
We now describe the results of our χ 2 analysis of the standard ∆m 2 -OSC+NP scenarios. As discussed in the previous section the analysis contains five parameters: ∆m 2 , θ, ∆δ n , ξ n and η n (or ε ′ , |ε| and arg(ε) for NSI). Our results are summarized in Figs. 2, 3, 4 and 5, where we show different projections of the allowed 5-dimensional parameter space.
In Figs. 2 and 3 we plot two-dimensional projections of the allowed parameter region for the analysis of atmospheric and K2K data in presence of ν µ → ν τ oscillations and different NP effects parametrized in the form Eq. (1). The corresponding results for the case of NSI are presented in Fig. 3. The regions in each panel are obtained after marginalization of χ 2 with respect to the three undisplayed parameters and they are defined at 90%, 95%, 99% and 3σ CL for 2 d.o.f. (∆χ 2 = 4.61, 5.99, 9.21 and 11.83, respectively).
The left panels in Figs. 2 and 3 show the projection of the allowed region on the oscillation parameters ∆m 2 -sin 2 θ plane. The best fit point is marked with a star. For the sake of comparison we also show the lines corresponding to the contours of the allowed regions for pure ∆m 2 -OSC and mark with a triangle the position of the best fit point. The results are shown for the chosen relative sign σ + n = +1. For σ + n = −1 the corresponding region would be obtained by sin 2 θ → 1 − sin 2 θ.
The regions on the right panels of Fig. 2 and 3 show the allowed values for the parameters characterizing the strength and mixing of the NP. The full regions correspond to arbitrary values of the relative phase η n (or equivalently to complex ε parameter for the NSI case) while the lines show the results for real relative phase η n ∈ {0, π} (which for NSI corresponds to ε real and either positive or negative, respectively). For this second case we show the allowed region for −1 ≤ sin 2ξ n ≤ 1 where for σ + n = +1 the sector with −1 ≤ sin 2ξ n ≤ 0 and 0 ≤ sin 2ξ n ≤ 1 correspond to η n = π and η n = 0, respectively, while the opposite The results are shown for the chosen relative sign σ + n = +1; for σ + n = −1 the corresponding region would be obtained by sin 2 θ → 1 − sin 2 θ. The regions on the right panels show the allowed values for the parameters characterizing the strength and mixing of the NP. The full regions corresponds to arbitrary values of the phase η n while the lines correspond to the case η n ∈ {0, π}. holds for σ + n = −1. As discussed in the previous section, for the case of arbitrary phase η n the full mixing parameter space can be covered by 0 ≤ ξ n ≤ π/4, which translates into the symmetry of the allowed region around ξ n = 0.
Several comments are in order. First, from the figures we see that the best fit point for the global ∆m 2 -OSC+NP scenarios is always very near the best fit point of pure ∆m 2 -OSC In other words, the data does not show any evidence of presence of NP even as a subdominant effect. Second, in agreement with SK analysis [2], we find that with the inclusion of the three-dimensional atmospheric fluxes and improved cross sections as well as with the reanalyzed data points from SK, the best fit point and corresponding allowed regions from the atmospheric+K2K neutrino analysis is shifted to slightly lower values of ∆m 2 compared to our previous analysis corresponding to the same data set [32]. Third, the figures illustrate the robustness of the allowed ranges of mass and mixing derived from the analysis of atmospheric and K2K data under the presence of these generic NP effects. Fourth the analysis allow us to derive well-defined upper bounds on the NP strength. From Fig. 2 we see that the bounds on the NP strength parameter ∆δ n tightens for larger values of ξ n , being this effect stronger for NP effects leading to sub-dominant oscillations with stronger energy dependence. In particular, for n = 1 the bound on ∆δ n for sin 2 2ξ n = 1 is a factor ∼ 50 stronger than that for ξ n = 0, while for n = 0 the variation of the bound on ∆δ n with ξ n is at most a factor ∼ 3. This behaviour can be qualitatively understood by studying the modification of the oscillation probability at the best fit point of oscillations, ∆m 2 best = 2.2 × 10 −3 eV 2 and sin 2 2θ best = 1, due to NP effects: ∆P ≡ P ∆m 2 +NP µµ − P ∆m 2 µµ P ∆m 2 µµ = tan 2 φ b − 1 + 2R n,b sin 2ξ n cos η n + R 2 n,b sin 2 2ξ n 1 + 2R n,b sin 2ξ n cos η n + R 2 where φ b = 2.8 (L/10 3 km) (GeV/E) is the ∆m 2 oscillation phase at the best fit point and R n,b = 0.91 × 10 21 (∆δ n /GeV 1−n ) (E/GeV) n+1 is the ratio of NP to the standard oscillation contributions evaluated at the best fit point of oscillations. From Eq. (19) we find that as long as φ b is small the dependence of ∆P on R n,b (and consequently on ∆δ n ) is stronger for larger values of | sin 2ξ n |, which explains the tightening of the bound on ∆δ n . This behaviour was found in Ref. [20] for the case with n = 1.
However, it is worth noticing that the characteristic value of φ b for which NP effects are relevant depends on n since as n increases the effect becomes important only for higher values of E (see Fig. 1). As a consequence, the characteristic φ b for n = 0 is larger than for n = 1. Numerical inspection of Eq. (19) also shows that the variation of the dependence of ∆P on R n,b with sin 2ξ n decreases as φ b increases. This explains the milder dependence of the bound on ∆δ n with the mixing angle sin 2ξ n for n = 0 as compared with the n = 1 case. The figure also illustrates that imposing that the Hamiltonian is real does not substantially affect these conclusions.
The same arguments apply to the results for NSI in Fig. 3. In particular one sees that, as expected, the results for NSI are very similar to to those derived for the n = 0 CPT-odd scenario with the identification in Eq. (II), ∆ρ NSI ∼ 7 × 10 −22 F GeV.
More quantitative conclusions on the robustness of the derived ranges for the oscillation parameters and on the bounds on the NP strength can be obtained from Figs. 4 and 5 where we plot the marginalized ∆χ 2 as a function of the oscillation parameters, ∆m 2 and sin 2 θ, and of the NP strength parameters, for different NP scenarios as labeled in the figures.
In Table I we list the 3σ allowed ranges for ∆m 2 and sin 2 θ. We read that the derived ranges are robust under the presence of these generic forms of NP whose only effect is slightly enlarging the allowed range of ∆m 2 by 15%, and the lower bound on sin 2 2θ by 7% at the 3σ level. We have verified that these conclusions hold for other scenarios characterized by different values of n.
In terms of specific forms of NP the bounds on ∆δ n imply that ATM+K2K limit the possible VLI in the ν µ -ν τ sector via CPT-even effects to |δv| ≤ 8.1 × 10 −25 (1.6 × 10 −24 ) (20) and the possible VEP is constrained to The figure is shown for σ + n = +1. As described in the previous section the results hold for σ + n = −1 with the exchange sin 2 θ → 1 − sin 2 θ (see discussion around Eq. (15)). The individual 3σ bounds in Table I can be read from the corresponding panel with the condition ∆χ 2 ≤ 9.

IV. CONCLUSIONS
In this work we have studied the robustness of our present determination of the neutrino masses and mixing from the analysis of the atmospheric and K2K data under the presence of some new physics effects in the ν µ -ν τ sector. In particular, we have performed a global analysis to atmospheric and K2K data for scenarios where vector-like or tensor-like new physics interactions affect the neutrino evolution together with the standard ∆m 2 -mixing effect.
We have concluded that the data does not show any evidence of these new physics effects even at the sub-dominant level. As a consequence the derived range of oscillation parameters is robust under the presence of those unknown effects. The quantification of this statement is shown in Figs. 4 and 5 and in Table I, from which we read that inclusion of these new physics effects can at most enlarge the allowed range of ∆m 2 by 15% and relax the lower bound on sin 2 2θ by 7% at the 3σ level.
From the analysis we have also derived upper bounds on the strength of the new physics effects in the ν µ -ν τ sector. In particular we show in Eqs. (20) and (22) the bound on the possible violation of Lorentz Invariance via CPT-even and CPT-odd effects in the neutrino evolution respectively. The constraint on the violation of the equivalence principle (VEP) due to non universal coupling of the neutrinos to gravitational potential is given in Eq. (21), while bounds on non-universal couplings of the neutrino to a torsion field are displayed in Eq. (23). The constraints on non-standard neutrino interactions with matter are shown in Eq. (24).  range of sin 2 θ corresponds to σ + n = +1. For σ + n = −1 the corresponding range would be obtained by sin 2 θ → 1 − sin 2 θ.

Acknowledgments
We are particularly indebted to M. Honda for providing us with their new 3-dimensional atmospheric neutrino fluxes. This work was supported in part by the National Science Foundation grant PHY0098527. MCG-G is also supported by Spanish Grants No FPA-2001-3031 and CTIDIB/2002/24.

APPENDIX A: STATISTICAL TREATMENT OF ATMOSPHERIC DATA
We summarize here our updated statistical analysis of the atmospheric data. For convenience we have adopted the pull method used previously by the SK Collaboration (see for instance Refs. [35,36] for details on their latest analysis) and by the Bari group [23,37]. There are however some technical differences which we describe next.
The basic idea of the pull method consists in parametrizing the systematic errors and the theoretical uncertainties in terms of a set of variables {ξ i }, called pulls, which are then treated on the same footing as the other parameters of the model. The χ 2 function can be decomposed into the sum of two parts: where ω denotes the parameters of the model, χ 2 data is the usual term describing the deviation of the experimental results from their theoretical predictions, and the extra term χ 2 pulls provides proper penalties to account for deviations of the systematics and the theoretical inputs from their nominal value. It is convenient to define the pulls in such a way that for each source of systematics or theoretical input i the value ξ i = 0 corresponds to the "expected value" reported by the collaboration or predicted by the theory, and ξ i = ±1 corresponds to a 1σ deviation.
For the Super-Kamiokande experiment χ 2 pulls ( ξ) can be properly written as a positive quadratic function of ξ i . The interpretations of the pulls is particularly transparent if the sources of uncertainties are selected to be independent of each other. In this case the pulls are uncorrelated and the expression of χ 2 pulls is very simple: In its original formulation, the set of pulls selected by Super-Kamiokande [35] did not verify this condition and a correlation matrix between the selected pulls (provided by SK collaboration from their MC simulation) had to be included. In our analysis, however, we have identified the dominant independent sources of systematic uncertainties in SK analysis, and we use them as the basis for our pulls. We have characterized the theoretical and systematic uncertainties in terms of 18 independent sources of error: 4 to parametrize the theoretical uncertainties associated to the atmospheric fluxes (which we describe in Sec. A 1), 6 for the theoretical uncertainties in the interaction cross sections (given in Sec. A 2) and 8 sources of experimental systematic errors (described in Sec. A 3). To the point to which the comparison is possible, this seems close to the approach followed by Super-Kamiokande in their latest analysis [36].
The form of χ 2 data depends on the expected distribution of the experimental results. Under the standard approximation of Gaussian distribution, we have: where R th n (R ex n ) is the ratio between the expected (observed) number of events and the theoretical Monte Carlo for the case of no oscillations. Note that the dependence of χ 2 data on both the parameters ω and the pulls ξ is entirely through R th n ( ω, ξ). In the pull approach, ω and ξ play a very similar role, and in principle should be treated in the same way. However, for the Super-Kamiokande experiment the bounds on ξ implied by χ 2 pulls are in general significantly stronger than those implied by χ 2 data , and it is therefore a good approximation to retain the dependence of χ 2 data on ξ only to the lowest orders. This is done by expanding R th n ( ω, ξ) in powers of ξ i up to the first order: It is easy to prove [37] that under the approximation (A4) Eq. (A1) is mathematically equivalent to the usual covariance definition of the χ 2 which we employed before in Refs. [32,39]. Thus the small difference in the results are not due to the different statistical treatment, but to differences either in the input parameters or in the updated values used for the systematic and theoretical uncertainties.
Furthermore within the present precision one can safely neglect the dependence of π i n on the neutrino parameters ω. With this approximation, we can write: where we have introduced the function χ 2 ( ω) = min {ξ i } χ 2 ( ω, ξ). It is clear from Eq. (A5) that in the present approach the systematic and theoretical uncertainties are completely characterized by the set of quantities {π i n }, which describe the strength of the "coupling" between the pull ξ i and the observable R th n . In the rest of this section we will discuss in detail how we have parametrized and taken into account the various sources of uncertainty and list the derived values for {π i n }.

Flux uncertainties
Flux uncertainties are theoretical uncertainties arising from our limited knowledge of the atmospheric neutrino fluxes. Following Refs. [35,36] we have parametrized them in terms of four pulls: ξ flux norm , ξ flux tilt , ξ flux ratio and ξ flux zenith . The corresponding coefficients π i n are listed in Table II.
• ξ flux norm is the pull associated to the total normalization error, which we set to σ flux norm =20% [38].
• ξ flux tilt is a"tilt" factor which parametrizes possible deviations of the energy dependence of the atmospheric fluxes from the simple power law. Following Refs. [35,36], we define: and assume an uncertainty on the factor δ, σ δ = 5% [35,36]. Also in analogy with Refs. [35,36] we have chosen E 0 = 2 GeV. We then calculate numerically the coefficients π tilt n as follows: we compute the expected number of events for a given bin N n using Φ δ (E) for the central value of δ and for δ ± σ δ and obtain the corresponding coupling π tilt n as the relative change in N n . The results reported in the second column of Table II are obtained neglecting the effect of oscillations. However we have verified that when the dependence of the π tilt n on the neutrino oscillation parameters is properly taken into account we find very similar results.
• ξ flux zenith describes the uncertainty on the zenith angle dependence, which we assume energy independent. As in Ref. [35] we parametrize the coupling for this pull for the  bin n of a given sample as π zenith n = 5% cos θ n . This means that this uncertainty can induce an error in the up/down asymmetry of events which we conservatively take to be 5%. In Ref. [35] the assumed up/down uncertainty was smaller (2.5%) and a separate zenith-pull was introduced for the horizontal-to-vertical ratio uncertainty of 2%. We have verified that within the present precision both parametrizations of the uncertainties in the zenith angle distribution give very similar results.

Cross-section uncertainties
Cross section uncertainties are theoretical uncertainties associated to our ignorance on the neutrino interaction cross section. In our calculation we follow the standard approach and consider separately the contributions to the interaction cross section from the exclusive channels of lower multiplicity: quasi-elastic scattering (QE), and single pion production (1π), and include all additional channels as part of the deep inelastic (DIS) cross section (also refer to as multi-pion). We neglect for simplicity coherent scattering on oxygen and neutralcurrent interactions, which contribute only marginally to the considered data samples.
We assume that each of these three contributions to the cross sections are subject to different sources of uncertainties which allow to consider the corresponding pulls as independent. For each type of neutrino interactions we introduce two pulls: • ξ QE norm , ξ 1π norm , ξ DIS norm describe the total normalization errors. We conservatively assume σ σ QE norm = 15% and σ σ 1π norm = 15%. For the normalization error of the DIS cross section we estimate σ σ DIS norm = 15% for contained events and σ σ DIS norm = 10% for upward-going muons from the spread of theoretical predictions arising from the use of different sets of nucleon structure functions. The relevant coefficients π i n are listed in Table III. They are obtained computing the relative change in the number of expected events in a given data sample arising from the use of either σ i or σ i ± σ σ i norm for each of the three contributions to the cross section.  • ξ QE ratio , ξ 1π ratio , ξ DIS ratio parametrize the uncertainty of the σ i,νµ /σ i,νe ratios. This error is relevant only for contained events, and it is much smaller than the total normalization uncertainty. The numbers listed in Table III are obtained from Ref. [35].

Systematic uncertainties
The systematics uncertainties of the Super-Kamiokande experiment are derived from Tables 9.2,9.3,9.4 and 9.5 of Ref. [35]. We include in our calculations the following sources of systematics: • ξ sys hadron is the pull for the uncertainty associated with the simulation of hadronic interactions; • ξ sys µ/e is the pull for the errors in the particle identification procedure; • ξ sys ring is the pull for the uncertainty coming from the ring-counting procedure; • ξ sys f-vol is the pull for the uncertainty in the fiducial volume determination; • ξ sys E-cal is the pull for the uncertainty in the energy calibration; • ξ sys PC-nrm is the pull for the relative normalization between partially-contained and fullycontained events.
• ξ sys track is the pull for the uncertainty in the track reconstruction of upgoing muons; • ξ sys up-eff is the pull for the uncertainty in the detection efficiency of upgoing muons and the stopping-thrugoing separation.