Updated Constraints on Non-Standard Interactions from Global Analysis of Oscillation Data

We quantify our present knowledge of the size and flavor structure of non-standard neutrino interactions which affect the matter background in the evolution of solar, atmospheric, reactor and long-baseline accelerator neutrinos as determined by a global analysis of oscillation data - both alone and in combination with the results on coherent neutrino-nucleus scattering from the COHERENT experiment. We consider general neutral current neutrino interactions with quarks whose lepton-flavor structure is independent of the quark type. We study the dependence of the allowed ranges of non-standard interaction coefficients, the status of the LMA-D solution, and the determination of the oscillation parameters on the relative strength of the non-standard couplings to up and down quarks. Generically we find that the conclusions are robust for a broad spectrum of up-to-down strengths, and we identify and quantify the exceptional cases related to couplings whose effect in neutrino propagation in the Earth or in the Sun is severely suppressed. As a result of the study we provide explicit constraints on the effective couplings which parametrize the non-standard Earth matter potential relevant for long-baseline experiments.


Introduction
Experiments measuring the flavor composition of neutrinos produced in the Sun, in the Earth's atmosphere, in nuclear reactors and in particle accelerators have established that lepton flavor is not conserved in neutrino propagation, but it oscillates with a wavelength which depends on distance and energy. This demonstrates beyond doubt that neutrinos are massive and that the mass states are non-trivial admixtures of flavor states [1,2], see Ref. [3] for an overview.
Under the assumption that the Standard Model (SM) is the low energy effective model of a complete high energy theory, neutrino masses emerge naturally as the first observable consequence in the form of the Weinberg operator [4], the only dimension five operator that can be built within the SM particle content. In this framework the next operators with observable consequences at low energies come at dimension six. They include four-fermion terms leading to Non-Standard Interactions (NSI) [5][6][7] between neutrinos and matter (for recent reviews, see [8,9]), both in charge-current interactions (NSI-CC) (ν α γ µ P L β )(f γ µ P f ) (1.1) and in neutral current interactions (NSI-NC) Here α, β are lepton flavor indices, f, f are SM charged fermions and γ µ are the Dirac gamma matrices; P L is the left-handed projection operator while P can be either P L or P R (the right-handed projection operator). These operators are expected to arise generically from the exchange of some mediator state assumed to be heavier that the characteristic momentum transfer in the ν interaction process. Since operators in both Eqs. (1.1) and (1.2) modify the inelastic neutrino scattering cross sections with other SM fermions they can be bounded by precision electroweak data (see for example Refs. [10][11][12]). In general these "scattering" bounds on NSI-CC operators are rather stringent, whereas the bounds on NSI-NC tend to be weaker. On the other hand, the operators in Eq. (1.2) can also modify the forward-coherent scattering (i.e., at zero momentum transfer) of neutrinos as they propagate through matter via so-called Mikheev-Smirnov-Wolfenstein (MSW) mechanism [5,13]. Consequently their effect can be significantly enhanced in oscillation experiments where neutrinos travel large regions of matter, such as is the case for solar and atmospheric neutrinos. Indeed, the global analysis of data from oscillation experiments in the framework of mass induced oscillations in presence of NSI currently provides some of the strongest constraints on the size of the NSI affecting neutrino propagation [14,15].
Of course, for models with a high energy New Physics scale, electroweak gauge invariance generically implies that the NSI-NC parameters are still expected to be subject to tight constraints from charged lepton observables [16,17], leading to no visible effect in oscillations. However, more recently it has been argued that viable gauge models with light mediators (i.e., below the electroweak scale) may lead to observable effects in oscillations without entering in conflict with other bounds [18][19][20][21] (see also Ref. [9] for a discussion). In particular, for light mediators bounds from high-energy neutrino scattering experiments such as CHARM [22] and NuTeV [23] do not apply. In this framework NSI-NC generated by mediators as light as about 10 MeV can only be constrained by their effect in oscillation data and by the recent results on coherent neutrino-nucleus scattering observed for the first time by the COHERENT experiment [24].
In this work we revisit our current knowledge of the size and flavor structure of NSI-NC which affect the matter background in the evolution of solar, atmospheric, reactor and longbaseline (LBL) accelerator neutrinos as determined by a global analysis of oscillation data. This updates and extends the analysis in Ref. [15] where NSI-NC with either up or down quarks were considered. Here we extend our previous study to account for the possibility of NSI with up and down quarks simultaneously, under the simplifying assumption that they carry the same lepton flavor structure. To this aim, in Sec. 2 we briefly summarize the framework of our study and discuss the simplifications used in the analysis of the atmospheric and LBL data on one side and of the solar and KamLAND sector on the other side. In Sec. 3 we present the results of the updated analysis of solar and KamLAND data and quantify the impact of the modified matter potential on the data description, as well as the status of the LMA-D solution [25] in presence of the most general NSI scenario considered here. In Sec. 4 we describe the constraints implied by the analysis of atmospheric, LBL and reactor experiments, and combine them with those arising from the solar+KamLAND data. We show how the complementarity and synergy of the different data sets is important for a robust determination of neutrino masses and mixing in the presence of these general NSI, and we derive the most up-to-date allowed ranges on NSI couplings. Finally in Sec. 5 we further combine the oscillation bounds with those from the COHERENT experiment and in Sec. 6 we summarize our conclusions. We present the details of the analysis of the IceCube results in Appendix A.

Formalism
In this work we will consider NSI affecting neutral-current processes relevant to neutrino propagation in matter. The coefficients accompanying the new operators are usually parametrized in the form: where G F is the Fermi constant, α, β are flavor indices, P ≡ P L , P R and f is a SM charged fermion. In this notation, ε f,P αβ parametrizes the strength of the new interaction with respect to the Fermi constant, ε f,P αβ ∼ O(G X /G F ). If we now assume that the neutrino flavor structure of the interactions is independent of the charged fermion type, we can factorize ε f,P αβ as the product of two terms: where the matrix ε η αβ describes the neutrino part and the coefficients ξ f,P parametrize the coupling to the charged fermions. Under this assumption the Lagrangian in Eq. (2.1) takes the form: As is well known, only vector NSI contribute to the matter potential in neutrino oscillations. It is therefore convenient to define: Ordinary matter is composed of electrons (e), up quarks (u) and down quarks (d). As stated in the introduction, in this work we restrict ourselves to non-standard interactions with quarks, so that only ξ u and ξ d are relevant for neutrino propagation. It is clear that a global rescaling of both ξ u and ξ d by a common factor can be reabsorbed into a rescaling of ε η αβ , so that only the direction in the (ξ u , ξ d ) plane is phenomenologically non-trivial. We parametrize such direction in terms of an angle η: where we have chosen the normalization so that η = arctan(1/2) ≈ 26.6 • corresponds to NSI with up quarks (ξ u = 1, ξ d = 0) while η = arctan(2) ≈ 63.4 • corresponds to NSI with down quarks (ξ u = 0, ξ d = 1). Note that the transformation η → η + π simply results in a sign flip of ξ u and ξ d , hence it is sufficient to consider −π/2 ≤ η ≤ π/2.

Neutrino oscillations in the presence of NSI
In general, the evolution of the neutrino and antineutrino flavor state during propagation is governed by the Hamiltonian: where H vac is the vacuum part which in the flavor basis (ν e , ν µ , ν τ ) reads Here U vac denotes the three-lepton mixing matrix in vacuum [1,26,27]. Following the convention of Ref. [28], we define U vac = R 23 (θ 23 )R 13 (θ 13 )R 12 (θ 12 , δ CP ), where R ij (θ ij ) is a rotation of angle θ ij in the ij plane andR 12 (θ 12 , δ CP ) is a complex rotation by angle θ 12 and phase δ CP . Explicitly: where c ij ≡ cos θ ij and s ij ≡ sin θ ij . This expression differs from the usual one "U " (defined, e.g., in Eq. (1.1) of Ref. [29]) by an overall phase matrix: U vac = P U P * with P = diag(e iδ CP , 1, 1). It is easy to show that in the absence of non-standard interactions such rephasing does not affect the expression of the probabilities and produces therefore no visible effect: in other words, when only standard interactions are considered the physical interpretation of the vacuum parameters (∆m 2 21 , ∆m 2 31 , θ 12 , θ 13 , θ 23 and δ CP ) is exactly the same in both conventions. The advantage of defining U vac as in Eq. (2.8) is that the CPT transformation H vac → −H * vac , whose relevance for the present work will be discussed below, can be implemented exactly (up to an irrelevant multiple of the identity) by the following transformation of the parameters: which does not spoil the commonly assumed restrictions on the range of the vacuum parameters (∆m 2 21 > 0 and 0 ≤ θ ij ≤ π/2). Concerning the matter part H mat of the Hamiltonian which governs neutrino oscillations, if all possible operators in Eq. (2.1) are added to the SM Lagrangian we get: where the "+1" term in the ee entry accounts for the standard contribution, and describes the non-standard part. Here N f (x) is the number density of fermion f as a function of the distance traveled by the neutrino along its trajectory. In Eq. (2.11) we have limited the sum the charged fermions present in ordinary matter, f = e, u, d. Since quarks are always confined inside protons (p) and neutrons (n), it is convenient to define: Taking into account that N u (x) = 2N p (x) + N n (x) and N d (x) = N p (x) + 2N n (x) and also that matter neutrality implies N p (x) = N e (x), Eq. (2.11) becomes: which shows that from the phenomenological point of view the propagation effects of NSI with electrons can be mimicked by NSI with quarks by means of a suitable combination of up-quark and down-quark contributions. Our choice of neglecting ε e αβ in this work does not therefore imply a loss of generality.
Since this matter term can be determined by oscillation experiments only up to an overall multiple of the identity, each ε f αβ matrix introduces 8 new parameters: two differences of the three diagonal real parameters (e.g., ε f ee − ε f µµ and ε f τ τ − ε f µµ ) and three off-diagonal complex parameters (i.e., three additional moduli and three complex phases). Under the assumption that the neutrino flavor structure of the interactions is independent of the charged fermion type, as described in Eq. (2.2), we can write ε p αβ = ε η αβ ξ p and ε n αβ = ε η αβ ξ n , which leads to: so that the phenomenological framework adopted here is characterized by 9 matter parameters: eight related to the matrix ε η αβ plus the direction η in the (ξ p , ξ n ) plane. We finish this section by reminding that as a consequence of the CPT symmetry, neutrino evolution is invariant if the Hamiltonian H ν = H vac + H mat is transformed as H ν → −(H ν ) * . This requires a simultaneous transformation of both the vacuum and the matter terms. The transformation of H vac is described in Eq. (2.9) and involves a change in the octant of θ 12 as well as a change in the neutrino mass ordering (i.e., the sign of ∆m 2 31 ), which is why it has been called "generalized mass ordering degeneracy" in Ref. [28]. As for H mat we need: see Refs. [15,28,30]. As seen in Eqs. (2.11), (2.13) and (2.14) the matrix ε αβ (x) depends on the chemical composition of the medium, which may vary along the neutrino trajectory, so that in general the condition in Eq. (2.15) is fulfilled only in an approximate way. The degeneracy becomes exact in the following two cases: 1 • if the effective NSI coupling to neutrons vanishes, so that ε n αβ = 0 in Eq. (2.13). In terms of fundamental quantities this occurs when ε u αβ = −2ε d αβ , i.e., the NSI couplings are proportional to the electric charge of quarks. In our parametrization this corresponds to η = 0 as shown in Eqs. (2.5) and (2.14); • if the neutron/proton ratio Y n (x) is constant along the entire neutrino propagation path. This is certainly the case for reactor and long-baseline experiments, where only the Earth's mantle is involved, and to a good approximation also for atmospheric neutrinos, since the differences in chemical composition between mantle and core can safely be neglected in the context of NSI [14]. In this case the matrix ε αβ (x) becomes independent of x and can be regarded as a new phenomenological parameter, as we will describe in Sec. 2.2.
Further details on the implications of this degeneracy for different classes of neutrino experiments (solar, atmospheric, etc.) will be provided later in the corresponding section.

Matter potential in atmospheric and long-baseline neutrinos
As discussed in Ref. [14], in the Earth the neutron/proton ratio Y n (x) which characterize the matter chemical composition can be taken to be constant to very good approximation. The PREM model [31] fixes Y n = 1.012 in the Mantle and Y n = 1.137 in the Core, with an average value Y ⊕ n = 1.051 all over the Earth. Setting therefore Y n (x) ≡ Y ⊕ n in Eqs. (2.11) and (2.13) we get ε αβ (x) ≡ ε ⊕ αβ with: If we drop ε e αβ and impose quark-lepton factorization as in Eq. (2.14) we get: In other words, within this approximation the analysis of atmospheric and LBL neutrinos holds for any combination of NSI with up, down or electrons and it can be performed in terms of the effective NSI couplings ε ⊕ αβ , which play the role of phenomenological parameters. In particular, the best-fit value and allowed ranges of ε ⊕ αβ are independent of η, while the bounds on the physical quantities ε η αβ simply scale as (cos η + Y ⊕ n sin η). Moreover, it is immediate to see that for η = arctan(−1/Y ⊕ n ) ≈ −43.6 • the contribution of NSI to the matter potential vanishes, so that no bound on ε η αβ can be derived from atmospheric and LBL data in such case.
Following the approach of Ref. [14], the matter Hamiltonian H mat , given in Eq. (2.10) after setting ε αβ (x) ≡ ε ⊕ αβ , can be parametrized in a way that mimics the structure of the vacuum term (2.7): where R ij (ϕ ij ) is a rotation of angle ϕ ij in the ij plane andR 23 (ϕ 23 , δ NS ) is a complex rotation by angle ϕ 23 and phase δ NS . Note that the two phases α 1 and α 2 included in Q rel are not a feature of neutrino-matter interactions, but rather a relative feature of the vacuum and matter terms. In order to simplify the analysis we neglect ∆m 2 21 and also impose that two eigenvalues of H mat are equal (ε ⊕ = 0). The latter assumption is justified since, as shown in Ref. [32], strong cancellations in the oscillation of atmospheric neutrinos occur when two eigenvalues of H mat are equal, and it is precisely in this situation that the weakest constraints can be placed. Setting ∆m 2 21 → 0 implies that the θ 12 angle and the δ CP phase disappear from the expressions of the oscillation probabilities, and the same happens to the ϕ 23 angle and the δ NS phase in the limit ε ⊕ → 0. Under these approximations the effective NSI couplings ε ⊕ αβ can be parametrized as: With all this the relevant flavor transition probabilities for atmospheric and LBL experiments depend on eight parameters: (∆m 2 31 , θ 13 , θ 23 ) for the vacuum part, (ε ⊕ , ϕ 12 , ϕ 13 ) for the matter part, and (α 1 , α 2 ) as relative phases. Notice that in this case only the relative sign of ∆m 2 31 and ε ⊕ is relevant for atmospheric and LBL neutrino oscillations: this is just a manifestation of the CPT degeneracy described in Eqs. (2.9) and (2.15) once ∆m 2 21 and ε ⊕ are set to zero [14].
As further simplification, in order to keep the fit manageable we assume real NSI, which we implement by choosing α 1 = α 2 = 0 with ϕ ij range −π/2 ≤ ϕ ij ≤ π/2. It is important to note that with these approximations the formalism for atmospheric and LBL data is CP-conserving. We will go back to this point when discussing the experimental results included in the analysis.
In addition to atmospheric and LBL experiments, important information on neutrino oscillation parameters is provided also by reactor experiments with a baseline of about 1 km. Due to the very small amount of matter crossed, both standard and non-standard matter effects are completely irrelevant for these experiments, and the corresponding P ee survival probability depends only on the vacuum parameters. However, in view of the high precision recently attained by both reactor and LBL experiments in the determination of the atmospheric mass-squared difference, combining them without adopting a full 3ν oscillation scheme requires a special care. In Ref. [33] it was shown that, in the limit ∆m 2 21 ∆m 2 31 as indicated by the data, the P µµ probability relevant for LBL-disappearance experiments can be accurately described in terms of a single effective mass parameter ∆m 2 In the rest of this work we will therefore make use of ∆m 2 µµ as the fundamental quantity parametrizing the atmospheric mass-squared difference. For each choice of the vacuum mixing parameters in U vac , the calculations for the various data sets are then performed as follows: • for atmospheric and LBL data we assume ∆m 2 21 = 0 and set ∆m 2 31 = ∆m 2 µµ ; • for reactor neutrinos we keep ∆m 2 21 finite and set ∆m 2 31 = ∆m 2 µµ + r 2 ∆m 2 21 .
In this way the information provided by reactor and long-baseline data on the atmospheric mass scale is consistently combined in spite of the approximation ∆m 2 21 → 0 discussed above. Note that the correlations between solar and reactor neutrinos are properly taken into account in our fit, in particular for what concerns the octant of θ 12 .

Matter potential for solar and KamLAND neutrinos
For the study of propagation of solar and KamLAND neutrinos one can work in the one mass dominance approximation, ∆m 2 31 → ∞ (which effectively means that In this approximation the survival probability P ee can be written as [34,35] P ee = c 4 13 P eff + s 4

13
(2.20) The probability P eff can be calculated in an effective 2 × 2 model described by the Hamiltonian H eff = H eff vac + H eff mat , with: where we have imposed the quark-lepton factorization of Eq. (2.14) and used the parameterization convention of Eq. (2.8) for U vac . The coefficients ε η D and ε η N are related to the original parameters ε η αβ by the following relations: Note that the δ CP phase appearing in Eq. (2.21) could be transferred to Eq. (2.22) without observable consequences by means of a global rephasing. Hence, for each fixed value of η the relevant probabilities for solar and KamLAND neutrinos depend effectively on six quantities: the three real oscillation parameters ∆m 2 21 , θ 12 and θ 13 , one real matter parameter ε η D , and one complex vacuum-matter combination ε η N e −iδ CP . As stated in Sec. 2.2 in this work we will assume real NSI, implemented here by setting δ CP = 0 and considering only real (both positive and negative) values for ε η N . Unlike in the Earth, the matter chemical composition of the Sun varies substantially along the neutrino trajectory, and consequently the potential depends non-trivially on the specific combinations of couplings with up and down quarks -i.e., on the value of η. This implies that the generalized mass-ordering degeneracy is not exact, except for η = 0 (in which case the NSI potential is proportional to the standard MSW potential and an exact inversion of the matter sign is possible). However, as we will see in Sec. 3, the CPT transformation described in Eqs. (2.9) and (2.15) still results in a good fit to the global analysis of oscillation data for a wide range of values of η, and non-oscillation data are needed to break this degeneracy [36,37]. Because of the change in the θ 12 octant implied by Eq. (2.9) and given that the standard LMA solution clearly favors θ 12 < 45 • , this alternative solution is characterized by a value of θ 12 > 45 • . In what follows we will denote it as "LMA-D" [25].

Analysis of solar and KamLAND data
Let us start by presenting the results of the updated analysis of solar and KamLAND experiments in the context of oscillations with the generalized matter potential in Eq. (2.22). For KamLAND we include the separate DS1, DS2, DS3 spectra [38] with reactor fluxes as determined by Daya-Bay [39]. In the analysis of solar neutrino data we consider the total rates from the radiochemical experiments Chlorine [40], Gallex/GNO [41] and SAGE [42], the results for the four phases of Super-Kamiokande [43][44][45][46] (including the 2055 days separate day and night spectra from Ref. [46] of Super-Kamiokande IV), the combined data of the three phases of SNO as presented in Ref. [47], and the results of both Phase-I and Phase-II of Borexino [48][49][50].
We present different projections of the allowed parameter space in Figs. 1-3. In the analysis we have fixed sin 2 θ 13 = 0.022 which is the best-fit value from the global analysis of 3ν oscillations [29,51]. 2 So for each value of η there are four relevant parameters: ∆m 2 21 , sin 2 θ 12 , ε η D , and ε η N . As mentioned above, for simplicity the results are shown for real ε η N . Also strictly speaking the sign of ε η N is not physically observable in oscillation experiments, as it can be reabsorbed into a redefinition of the sign of θ 12 . However, for definiteness we have chosen to present our results in the convention θ 12 ≥ 0, and therefore we consider both positive and negative values of ε η N . Fig. 1 shows the two-dimensional projections on the oscillation parameters (θ 12 , ∆m 2 21 ) for different values of η after marginalizing over the NSI parameters, while Fig. 2 shows the corresponding two-dimensional projections on the matter potential parameters (ε η D , ε η N ) after marginalizing over the oscillation parameters. The one-dimensional ranges for the four parameters as a function of η are shown in Fig. 3.
The first thing to notice in the figures is the presence of the LMA-D solution for a wide range of values of η. This is a consequence of the approximate degeneracy discussed in the previous section. In particular, as expected, for η = 0 the degeneracy is exact and the LMA-D region in Fig. 1 is perfectly symmetric to the LMA one with respect to maximal θ 12 . Looking at the corresponding panels of Fig. 2 we note that the allowed area in the NSI parameter space is composed by two disconnected regions, one containing the SM case (i.e., the point ε η D = ε η N = 0) which corresponds to the "standard" LMA solution in the presence of the modified matter potential, and another which does not include such point igure 1. Two-dimensional projections of the 1σ, 90%, 2σ, 99% and 3σ CL (2 dof) allowed regions from the analysis of solar and KamLAND data in the presence of non-standard matter potential for the oscillation parameters (θ 12 , ∆m 2 21 ) after after marginalizing over the NSI parameters and for θ 13 fixed to sin 2 θ 13 = 0.022. The best-fit point is marked with a star. The results are shown for fixed values of the NSI quark coupling parameter η. For comparison the corresponding allowed regions for the analysis in terms of 3ν oscillations without NSI are shown as black void contours.
igure 2. Two-dimensional projections of the 1σ, 90%, 2σ, 99% and 3σ CL (2 dof) allowed regions from the analysis of solar and KamLAND data in the presence of non-standard matter potential for the matter potential parameters (ε η D , ε η N ), for sin 2 θ 13 = 0.022 and after marginalizing over the oscillation parameters. The best-fit point is marked with a star. The results are shown for fixed values of the NSI quark coupling parameter η. The panels with a scale factor "[×N ]" in their lower-left corner have been "zoomed-out" by such factor with respect to the standard axis ranges, hence the grey square drawn in each panel always corresponds to max |ε η D |, |ε η N | = 2 and has the same size in all the panels. For illustration we also show as shaded green areas the 90% and 3σ CL allowed regions from the analysis of the atmospheric and LBL data. Note that, as a consequence of the periodicity of η, the regions in the first (η = −90 • ) and last (η = +90 • ) panels are identical up to an overall sign flip. disfavored and does not appear at the displayed CL's.
In order to further illustrate the η dependence of the results, it is convenient to introduce the functions χ 2 LMA (η) and χ 2 LMA-D (η) which are obtained by marginalizing the χ 2 for a given value of η over both the oscillation and the matter potential parameters with the constraint θ 12 < 45 • and θ 12 > 45 • , respectively. With this, in the left panel of Fig. 4 we plot the differences χ 2 LMA (η) − χ 2 no-NSI (full lines) and χ 2 LMA-D (η) − χ 2 no-NSI (dashed lines), where χ 2 no-NSI is the minimum χ 2 for standard 3ν oscillations (i.e., without NSI), while in the right panel we plot χ 2 LMA-D (η) − χ 2 LMA (η) which quantifies the relative quality of the LMA and LMA-D solutions. From this plot we can see that even for the analysis of solar and KamLAND data alone (red lines) the LMA-D solution is disfavored at more than 3σ when η −40 • or η 86 • . Generically for such range of η the modified matter potential in the Sun, which in the presence of NSI is determined not only by the density profile but also by the chemical composition, does not allow for a degenerate solution compatible with KamLAND data. In particular, as discussed below, for a fraction of those η values the NSI contribution to the matter potential in the Sun becomes very suppressed and therefore the degeneracy between NSI and octant of θ 12 cannot be realized. In what respects the LMA solution, we notice that it always provides a better fit (or equivalent for η = 0) than the LMA-D solution to solar and KamLAND data, for any value of η. This does not have to be the case in general, and indeed it is no longer so when atmospheric data are also included in the analysis. We will go back to this point in the next section.
From the left panel in Fig. 4 we see that the introduction of NSI can lead to a substantial improvement in the analysis of solar and KamLAND data, resulting in a sizeable decrease of the minimum χ 2 with respect to the standard oscillation scenario. The maximum gain occur for η −64 • and is about 11.2 units in χ 2 (i.e., a 3.3σ effect), although for most of the values of η the inclusion of NSI improves the combined fit to solar and KamLAND by about 2.5σ. This is mainly driven by the well known tension between solar and KamLAND data in the determination of ∆m 2 21 . The phenomenological status of such tension has not changed significantly over the last lustrum, and arises essentially from a combination of two effects: (a) the 8 B measurements performed by SNO, SK and Borexino does not show any evidence of the low energy spectrum turn-up expected in the standard LMA-MSW [5,13] solution for the value of ∆m 2 21 favored by KamLAND, and (b) the observation of a non-vanishing day-night asymmetry in SK, whose size is considerably larger than what predicted for the ∆m 2 21 value indicated of KamLAND. With the data included in the analysis this results into a tension of ∆χ 2 ∼ 7.4 for the standard 3ν oscillations. Such tension can be alleviated in presence of a non-standard matter potential, thus leading to the corresponding decrease in the minimum χ 2 for most values of η -with the exception of the range −70 • η −60 • . Furthermore, as seen in the lower panel in Fig. 3 the allowed range of ∆m 2 21 implied by the combined solar and KamLAND data is pretty much independent of the specific value of η, except again for −70 • η −60 • in which case it can extend well beyond the standard oscillation LMA values.
The special behaviour of the likelihood of solar and KamLAND in the range −70 • η −60 • is a consequence of the fact that for such values the NSI contributions to the matter potential in the Sun approximately cancel. As mentioned in the previous section, the  Figure 3. 90% and 3σ CL (1 dof) allowed ranges from the analysis of solar and KamLAND data in the presence of non-standard neutrino-matter interactions, for the four relevant parameters (the matter potential parameters ε η D and ε η N as well as the oscillation parameters ∆m 2 21 and sin 2 θ 12 ) as a function of the NSI quark coupling parameter η, for sin 2 θ 13 = 0.022. In each panel the three undisplayed parameters have been marginalized.
matter chemical composition of the Sun varies substantially along the neutrino production region, with Y n (x) dropping from about 1/2 in the center to about 1/6 at the border of the solar core. Thus for −70 • η −60 • (corresponding to −2.75 tan η −1.75) the effective NSI couplings ε αβ (x) = ε p αβ + Y n (x)ε n αβ ∝ 1 + Y n (x) tan η → 0 vanish at some point inside the neutrino production region. This means that for such values of η the constraints on the NSI couplings from solar data become very weak, being prevented from disappearing completely only by the gradient of Y n (x). This is visible in the two upper panels in Fig. 3 and in the panels of Fig. 2 with η in such range, where a multiplicative factor 2-8 has to be included to make the regions fit in the same axis range. Indeed for those  values of η the allowed NSI couplings can be so large that their effect in the propagation of long-baseline reactor neutrinos through the Earth becomes sizeable, and can therefore lead to spectral distortions in KamLAND which affect the determination of ∆m 2 21 -hence the "migration" and distortion of the LMA region observed in the corresponding panels in Fig. 1. In particular, it is precisely for η = −64 • for which the "migration" of the KamLAND region leads to the best agreement with the solar determination of ∆m 2 12 , whereas for η = −68 • we find the worst agreement. In any case, looking at the shaded green regions in the corresponding panels of Fig. 2 we can anticipate that the inclusion of atmospheric and LBL oscillation experiments will rule out almost completely such very large NSI values.
As for θ 12 , looking at the relevant panel in Fig. 3 we can see that its determination is pretty much independent of the value of η, however a comparison between colored and void regions in Fig. 1 shows that its allowed range always extends to lower values than in the standard 3ν case without NSI. This is expected since the presence of non-diagonal NSI parametrized by ε η N provides another source of flavor transition, thus leading to a weakening of the lower bound on θ 12 .
We finish this section by noticing that two of the panels in Figs. 1 and 2 correspond to the values of NSI only with f = u (η ≈ 26.6 • ) and only with f = d (η ≈ 63.4 • ) and can be directly compared with the results of our previous global OSC+NSI analysis in Ref. [15]. For illustration we also show in one of the panels the results for η = −44 • which is close to the value for which NSI effects in the Earth matter cancel.

Results of the global oscillation analysis
In addition to the solar and KamLAND data discussed so far, in our global analysis we also consider the following data sets: • atmospheric neutrino data: this sample includes the four phases of Super-Kamiokande (up to 1775 days of SK4 [52]) in the form of the "classical" samples of e-like and µ-like events (70 energy and zenith angle bins), together with the complete set of DeepCore 3-year µ-like events (64 data points) presented in Ref. [53] and publicly released in Ref. [54]. The calculations of the event rates for both detectors are based on the atmospheric neutrino flux calculations described in Ref. [55]. In addition, we also include the results on ν µ -induced upgoing muons reported by IceCube [56][57][58], based on one year of data taking; • long-baseline experiments: we include here the ν µ andν µ disappearance as well as the ν e andν e appearance data in MINOS [59] (39, 14, 5, and 5 data points, respectively), the ν µ andν µ disappearance data in T2K [60] (39 and 55 data points, respectively), and the ν µ disappearance data in NOνA [61] (72 data points). As mentioned in Sec. 2, in order to keep the fit manageable we restrict ourselves to the CP-conserving scenario. At present, the results of the full 3ν oscillation analysis with standard matter potential show a hint of CP violation [29,51], which is mainly driven by the LBL ν e andν e appearance data at T2K [60] and NOνA [61]. Conversely, allowing for CP violation has negligible impact on the determination of the CP-conserving parameters in the analysis of MINOS appearance data and of any LBL disappearance data samples, as well as in our analysis of atmospheric events mentioned above. Hence, to ensure full consistency with our CP-conserving parametrization we have chosen not to include in the present study the data from the ν e andν e appearance channels in NOνA and T2K. This also renders our fit only marginally sensitive to the neutrino mass ordering.
In what follows we will refer to the long-baseline data included here as LBL-CPC; • medium-baseline (MBL) reactor experiments: since these experiments are largely insensitive to matter effects (either standard or non-standard), the results included here coincide with those of the standard 3ν analysis presented in Ref. [51] and illustrated in the black lines of the plot tagged «Synergies: determination of ∆m 2 3 ». Such analysis is based on a reactor-flux-independent approach as described in Ref. [62], and includes the Double-Chooz FD-I/ND and FD-II/ND spectral ratios with 455-day (FD-I), 363-day (FD-II), and 258-day (ND) exposures [63] (56 data points), the Daya-Bay 1230-day EH2/EH1 and EH3/EH1 spectral ratios [64] (70 data points), and the Reno 1500-day FD/ND spectral ratios [65] (26 data points). Figure 5 the two-dimensional projections of the allowed regions in the Earth's matter potential parameters ε ⊕ , ϕ 12 and ϕ 13 (i.e., in the parametrization of Eq. (2.19) with α i = 0) after marginalizing over the oscillation parameters. The green regions show the 90% and 3σ confidence regions (2 dof) from the analysis of atmospheric, LBL-CPC and MBL reactor experiments. Besides the increase in statistics on low-energy atmospheric events provided by the updated Super-Kamiokande and the new DeepCore data samples, the main difference with respect to the analysis in Refs. [14,15] is the inclusion of the bounds on NSI-induced ν µ disappearance provided by IceCube high-energy data as well as the precise information on θ 13 and |∆m 2 31 | from MBL reactor experiments.  Figure 5. Two-dimensional projections of the allowed regions onto the matter potential parameters ε ⊕ , ϕ 12 , and ϕ 13 after marginalization with respect to the undisplayed parameters. The large green regions correspond to the analysis of atmospheric, LBL-CPC, and MBL reactor data at 90% and 3σ CL. For comparison we show in yellow the corresponding results when omitting IceCube and reactor data. The solid colored regions show the 1σ, 90%, 2σ, 99% and 3σ CL allowed regions once solar and KamLAND data are included. The best-fit point is marked with a star.

Let us begin by showing in
To illustrate their impact we show as yellow regions the results obtained when IceCube and reactor data are omitted. For what concerns the projection over the matter potential parameters shown here, we have verified that the difference between the yellow and green regions is mostly driven by IceCube, which restricts the allowed values of the ϕ 12 for |ε ⊕ | ∼ 0.1-1. This can be understood since, for neutrino with energies above O(100 GeV), the vacuum oscillation is very suppressed and the survival probability of atmospheric ν µ arriving at zenith angle Θ ν is dominated by the matter induced transitions  Figure 6. Two-dimensional projections of the allowed regions onto different vacuum parameters after marginalizing over the matter potential parameters (including η) and the undisplayed oscillation parameters. The solid colored regions correspond to the global analysis of all oscillation data, and show the 1σ, 90%, 2σ, 99% and 3σ CL allowed regions; the best-fit point is marked with a star. The black void regions correspond to the analysis with the standard matter potential (i.e., without NSI) and its best-fit point is marked with an empty dot. For comparison, in the left panel we show in red the 90% and 3σ allowed regions including only solar and KamLAND results, while in the right panels we show in green the 90% and 3σ allowed regions excluding solar and KamLAND data, and in yellow the corresponding ones excluding also IceCube and reactor data.
first oscillation maximum for some of the trajectories. Also, the effective parameter ϕ µµ entering in the expression of P µµ depends linearly on ϕ 12 and only quadratically on ϕ 13 , which explains why the bounds on the mixings are stronger for ϕ 12 than for ϕ 13 . As can be seen in Fig. 5, even with the inclusion of IceCube neither upper nor lower bounds on the overall strength of the Earth's matter effects, ε ⊕ , can be derived from the analysis of atmospheric, LBL-CPC and MBL reactor experiments [14,32,67]. 3 This happens because the considered data sample is mainly sensitive to NSI through ν µ disappearance, and lacks robust constraints on matter effects in the ν e sector. As a consequence, when marginalizing over ε ⊕ (as well as over the oscillation parameters) the full flavor projection (ϕ 12 , ϕ 13 ) plane is allowed. On the other hand, once the results of solar and KamLAND experiments (which are sensitive to ν e ) are included in the analysis a bound on ε ⊕ is obtained and the flavor structure of the matter potential in the Earth is significantly constrained.
In Fig. 6 we show the two-dimensional projections of the allowed regions from the global analysis onto different sets of oscillation parameters. These regions are obtained after marginalizing over the undisplayed vacuum parameters as well as the NSI couplings. For comparison we also show as black-contour void regions the corresponding results with the standard matter potential, i.e., in the absence of NSI. As discussed in Sec. 2.2, in the right panels we have chosen to plot the regions in terms of the effective mass-squared difference relevant for ν µ disappearance experiments, ∆m 2 µµ . Notice that, having omitted NOνA and T2K appearance data and also set ∆m 2 21 = 0 in atmospheric and LBL-CPC experiments, the impact of the mass ordering on the results of the fit is greatly reduced.
This figure clearly shows the robustness of the determination of the ∆m 2 21 , |∆m 2 µµ | and θ 23 vacuum oscillation parameters even in the presence of the generalized NSI interactions. This result relies on the complementarity and synergies between the different data sets, which allows to constrain those regions of the parameter space where cancellations between standard and non-standard effects occur in a particular data set. To illustrate this we show as shaded regions the results obtained when some of the data are removed. For example, comparing the solid colored regions with the shaded red ones in the left panel we see how, in the presence of NSI with arbitrary values of η, the precise determination of ∆m 2 21 requires the inclusion of atmospheric, LBL-CPC and MBL reactor data: if these sets are omitted, the huge values of the NSI couplings allowed by solar data for −70 • η −60 • destabilize KamLAND's determination of ∆m 2 21 , as discussed in Sec. 3. The inclusion of these sets also limits the margins for NSI to alleviate the tension between solar and KamLAND data on the preferred ∆m 2 21 value, as can be seen by comparing the full dark-blue and red lines in the left panel of Fig. 4: indeed, in the global analysis the best-fit is achieved for η −44 • , which is precisely when the NSI effects in the Earth matter cancel so that no restriction on NSI contributions to solar and KamLAND data is imposed.
In the same way we see on the right panels that, if the solar and KamLAND data are removed from the fit, the determination of ∆m 2 µµ and θ 23 degrades because of the possible cancellations between NSI and mass oscillation effects in the relevant atmospheric and LBL-CPC probabilities. As NSI lead to energy-independent contributions to the oscillation phase, such cancelations allow for larger values of |∆m 2 µµ |. Comparing the yellow and green regions we see the inclusion MBL reactor experiments, for which NSI effects are irrelevant due to the short baselines involved, is crucial to reduce the degeneracies and provide a NSIindependent measurement of |∆m 2 µµ |. Even so, only the inclusion of solar and KamLAND allows to recover the full sensitivity of atmospheric and LBL-CPC experiments and derive limits on ∆m 2 µµ and θ 23 as robust as the standard ones. The most dramatic implications of NSI for what concerns the determination of the oscillation parameters affect θ 12 . In particular, for generic NSI with arbitrary η the LMA-D solution is still perfectly allowed by the global oscillation analysis, as indicated by the presence of the corresponding region in the left panel in Fig. 6. Turning to Fig. 4 we see that even after including all the oscillation data (dark-blue lines) the LMA-D solution is allowed at 3σ for −38 • η 87 • (as well as in a narrow window around η −65 • ), and indeed for −28 • η 0 • it provides a slightly better global fit than LMA. From Fig. 6 we also see that the lower bound on θ 12 in the presence of NSI is substantially weaker than the standard 3ν case. We had already noticed such reduction in the analysis of solar and KamLAND data for any value of η; here we point out that the cancellation of matter effects in the Earth for η ≈ −43.6 • prevents any improvement of that limit from the addition of Earth-based oscillation experiments.
The bounds on the five relevant NSI couplings (two diagonal differences and three nondiagonal entries) from the global oscillation analysis are displayed in Fig. 7 as a function of η. Concretely, for each value of η we plot as vertical bars the 90% and 3σ allowed ranges (1 dof) after marginalizing with respect to the undisplayed parameters. The left and right panels correspond to the limits for θ 12 within the LMA and LMA-D solution, respectively, both defined with respect to the same common minimum for each given η. For the sake of convenience and comparison with previous results we list in the first columns in Table 1 the 95% CL ranges for NSI with up-quarks only (η ≈ 26.6 • ), down-quarks only (η ≈ 63.4 • ) and couplings proportional to the electric charge (η = 0 • ); in this last case we have introduced an extra √ 5 normalization factor so that the quoted bounds can be directly interpreted in terms of ε p αβ . Let us point out that the sign of each non-diagonal ε η αβ can be flipped away by a suitable change of signs in some of the mixing angles; it is therefore not an intrinsic property of NSI, but rather a relative feature of the vacuum and matter Hamiltonians. Thus, strictly speaking, once the results are marginalized with respect to all the other parameters in the most general parameter space, the oscillation analysis can only provide bounds on |ε η α =β |. However, for definiteness we have chosen to restrict the range of the mixing angles to 0 ≤ θ ij ≤ π/2 and to ascribe the relative vacuum-matter signs to the NSI couplings, so that the ranges of the non-diagonal ε η αβ in Figs. 7 and 8 as well as in Table 1 are given for both signs.
From Fig. 7 and Table 1 we see that the allowed range for all the couplings (except ε η ee − ε η µµ ) obtained marginalizing over both θ 12 octants, which we denote in the table as LMA ⊕ LMA-D, is only slighter wider than what obtained considering only the LMA solution. Conversely, for ε η ee − ε η µµ the allowed range is composed by two disjoint intervals, each one corresponding to a different θ 12 octant. Note that for this coupling the interval associated with the LMA solution is not centered at zero due to the tension between the value of ∆m 2 21 preferred by KamLAND and solar experiments, even after including the bounds from atmospheric and long-baseline data. In general, we find that the allowed ranges for all the couplings do not depend strongly on the value of η as long as η differs enough from the critical value η ≈ −43.6 • . As already explained, at this point non-standard interactions in the Earth cancel out, so that no bound on the NSI parameters can be derived from any Earth-based experiment. This leads to a breakdown of the limits on ε η αβ , since solar data are only sensitive to the ε η D and ε η N combinations and cannot constrain the five NSI couplings simultaneously. In addition to the region around η ≈ −43.6 • , there is also some mild weakening of the bounds on NSI couplings involving ν e for −70 • η −60 • , corresponding to the window where NSI effects in the Sun are suppressed. Apart from these special cases, the bounds quoted in Table 1 are representative of the characteristic sensitivity to the NSI coefficients from present oscillation experiments, which at 95% CL ranges from O(1%) for |ε η µτ | to O(30%) for |ε η eτ | -the exception being, of course, ε η ee − ε η µµ .  Figure 7. 90%, and 3σ CL (1 dof) allowed ranges for the NSI couplings from the global oscillation analysis in the presence of non-standard matter potential as a function of the NSI quark coupling parameter η. In each panel the undisplayed parameters have been marginalized. On the left panels the oscillation parameters have been marginalized within the LMA region while the right panels corresponds to LMA-D solutions. The ranges are defined with respect to the minimum for each η.

Combined analysis of oscillation and COHERENT data
To conclude our study, let us now quantify the impact of adding to our fit the constraints on coherent neutrino-nucleus scattering from the first results of the COHERENT experiment [24]. As discussed in the introduction, while the bounds from oscillation effects apply  Table 1. 2σ allowed ranges for the NSI couplings ε u αβ , ε d αβ and ε p αβ as obtained from the global analysis of oscillation data (left column) and also including COHERENT constraints. The results are obtained after marginalizing over oscillation and the other matter potential parameters either within the LMA only and within both LMA and LMA-D subspaces respectively (this second case is denoted as LMA ⊕ LMA-D).
to models where the NSI are generated by mediators of arbitrarily light masses, for scattering experiments there is a minimum mediator mass below which the contact interaction approximation is not adequate to describe the ν interactions in the detector. For COHER-ENT this threshold is estimated to be O(10 MeV), hence the bounds presented here apply for models for which the mediator responsible for the NSI is heavier than about 10 MeV.
For the statistical analysis of the COHERENT results we follow Ref. [37] and construct χ 2 COH using just the total number of events, according to the expression given in the supplementary material of Ref. [24]. The predicted number of signal events N NSI can be expressed as: where γ is an overall normalization constant, the coefficients f νe = 0.31, f νµ = 0.19, and fν µ = 0.50 are the relative contributions from the three flux components (ν µ ,ν µ and ν e ), and the terms Q 2 wα encode the dependence on the NSI couplings: where we used the effective matrices ε p αβ and ε n αβ defined in Eq. (2.12). In this expression i ∈ {Cs, I} is the sum over the target nuclei, Z i and N i are the corresponding number of protons and neutrons (Z Cs = 55, N Cs = 78 for cesium and Z I = 53, N I = 74 for iodine), and g V p = 1/2 − 2 sin 2 θ W and g V n = −1/2 are the SM vector couplings of the Z boson to protons and neutrons, respectively, with θ W being the weak mixing angle. Note that the neutron/proton ratio in the two target nuclei is very similar, N Cs Z Cs 1.419 for cesium and N I Z I 1.396 for iodine, with an average value Y coh n = 1.407. We can therefore approximate Eq. (5.2) as: After imposing quark-lepton factorization from Eq. (2.14), ε coh αβ can be written as: This is expression formally is identical to Eq. (2.17), except for the numerical value of Y n . This suggests that the analysis of Earth-based oscillation experiments and of coherent scattering data share a number of phenomenological features. In particular, the best-fit value and allowed ranges of ε coh αβ implied by COHERENT are independent of η, while the corresponding bounds on the physical quantities ε η αβ simply scale as (cos η + Y coh n sin η). Also, for η = arctan(−1/Y coh n ) ≈ −35.4 • no bound on ε η αβ can be derived from COHERENT data.
The results of the global analysis of oscillation plus COHERENT data are shown as cyan lines in Fig. 4; the corresponding ranges for the NSI coefficients are shown in Fig. 8 and in the right column of Table 1. As can be seen, the main impact of including COHERENT data is to strongly disfavor the LMA-D solution for a wide range of η. LMA-D is allowed below 3σ only for −38 • η 14 • . This generalizes the results of Ref. [37] to a wider set of NSI-NC with quarks. We also find that the allowed ranges of flavor non-diagonal NSI couplings are moderately reduced. More interestingly, the addition of COHERENT data allows to derive constraints on each of the diagonal parameters separately. This is especially relevant for ε η τ τ for which the bounds become more than an order of magnitude stronger than previous indirect (loop induced) limits [10] for most η values. We notice, however, that COHERENT data are still not strong enough to disfavor the large ranges of NSI allowed by oscillations for η ≈ −43.6 • . Moreover, the cancellation of NSI effects in COHERENT data for η ≈ 35.4 • implies that no separate reconstruction of the diagonal parameters is possible around such value.
We finish by quantifying the results of our analysis in terms of the effective NSI parameters which describe the generalized Earth matter potential and are, therefore, the relevant quantities for the study of long-baseline experiments. The results are shown in Fig. 9 where we plot the dependence of the global χ 2 on each NSI effective couplings after marginalization over all other parameters. 4 Let us point out that, if only the results from Earth-based experiments such as atmospheric, long-baseline and reactor data were included in the analysis, the curves would be independent of η. However, when solar experiments and COHERENT data are also considered the global χ 2 becomes sensitive to the value of η. Given that, what we quantify in Fig. 9 is our present knowledge of the matter potential for neutrino propagation in the Earth for any unknown value of η. Technically this is obtained by marginalizing the results of the global χ 2 with respect to η as well, so that the ∆χ 2 functions plotted in the figure are defined with respect to the absolute minimum for any η (which, as discussed above and shown in Fig. 4, lies close to η ∼ −45 • ). In the upper panels the oscillation parameters have been marginalized within the LMA solution and in the lower ones within the LMA-D solution. Comparing the blue and cyan lines we conclude that COHERENT has a sizeable impact on the results for the LMA-D solution, whereas within the LMA region its present contribution to the determination of the generalized Earth matter potential is marginal.

Summary
In this work we have presented an updated analysis of neutrino oscillation results with the aim of establishing how well we can presently determine the size and flavor structure of NSI-NC which affect the evolution of neutrinos in a matter background. In particular we have extended previous studies by considering NSI with an arbitrary ratio of couplings to up and down quarks (parametrized by an angle η) and a lepton-flavor structure independent of the quark type (parametrized by a matrix ε η αβ ). We have included in our fit all the solar, atmospheric, reactor and accelerator data commonly used for the standard 3ν oscillation analysis, with the only exception of T2K and NOνA appearance data whose recent hints in favor of CP violation are not easily accommodated within the CP-conserving approximation assumed in this work. In addition, we have considered the recent results on coherent neutrino-nucleus scattering from the COHERENT experiment. We have found that: • classes of experiments which are sensitive to NSI only through matter characterized by a limited range of proton/neutron ratio Y n unavoidably exhibit suppression of NSI effects for specific values of η. This is the case for solar data at −70 • η −60 • , for Earth-based (atmospheric, long-baseline, reactor) experiments at η ≈ −44 • , and for COHERENT scattering data η ≈ −35 • . Such cancellations limit the sensitivity to the NSI couplings; • moreover, the interplay between vacuum and matter contributions to the flavor transition probabilities in classes of experiments with limited energy range and/or sensitive only to a specific oscillation channel spoils the accurate determination of the oscillation parameters achieved in the standard 3ν scenario. This is particularly visible in ∆m 2 21 and θ 12 as determined by solar and KamLAND data, as well as in ∆m 2 31 and θ 23 as determined by atmospheric, LBL-CPC and MBL reactor data; • however, both problems can be efficiently resolved by combining together different classes of experiments, so to ensure maximal variety of matter properties, energy ranges, and oscillation channels. In particular, our calculations show that the precise determination of the vacuum parameters is fully recovered (except for θ 12 ) in a joint analysis of solar and Earth-based oscillation experiments, even when arbitrary values of η are considered; • the well-known LMA-D solution, which arises in the presence of of NSI as a consequence of CPT invariance, is allowed at 3σ for −38 • η 87 • from the global analysis of oscillation data. The inclusion of the COHERENT results considerably improves this situation, however even in that case the LMA-D region remains allowed at the 3σ level for −38 • η 14 • .
In addition, we have determined the allowed range of the NSI couplings ε η αβ as a function of the up-to-down coupling η, showing that such constraints are generically robust except for a few specific values of η where cancellations occurs. Finally, in view of the possible implications that generic NSI-NC may have for future Earth-based facilities, we have recast the results of our analysis in terms of the effective NSI parameters ε ⊕ αβ which describe the generalized matter potential in the Earth, and are therefore the relevant quantities for the study of atmospheric or long-baseline experiments.

A Details of the IceCube fit
The number of events measured by the IceCube detector have been provided in a grid with 210 bins [56,58], which depends on the reconstructed neutrino energy (logarithmically spaced in 10 bins ranging from 400 GeV to 20 TeV) and the reconstructed neutrino direction (linearly spaced in 21 bins from cos Θ ν = −1.02 to cos Θ ν = 0.24). To reproduce the number of events of each bin we have computed where φ atm µ,± (E ν , Θ ν ) is the atmospheric muon neutrino flux for neutrinos (+) and antineutrinos (−). Among the different alternatives provided by the IceCube collaboration we have chosen to consider those tagged as "initial", which do not include propagation effects across the Earth. Here A eff i,± (E ν , Θ ν ) is the effective area encoding the detector response to a ν µ with energy E ν and direction Θ ν for the bin 'i'. As effective area we have used the nominal detector response. The quantity P ± µµ (E ν , Θ ν ) is the flavor oscillation probability averaged over the altitude of the neutrino production point, defined as: where κ ± (E ν , Θ ν , h) is the altitude distribution of the flux normalized to one [55], X n (Θ ν ) is the column density along the neutrino trajectory for the nucleon n ∈ {proton, neutron} and σ ± n (E ν ) is the corresponding inclusive cross-section for ν µ . Hence P ± µµ (E ν , Θ ν ) also includes the neutrino flux absorption by the Earth.
In order to reproduce the published fit [58] we need to include in the χ 2 the contribution the systematic uncertainties for every point in the parameter space. Such systematics are included by the collaboration either as a discrete or a continuous nuisance parameter. In our analysis all the systematics are treated as continuous quantities and their effects on the number of events are assumed to be linear. We can divide systematics into two classes: those related to the neutrino flux, and those related to the detector response and the optical properties of the ice. The atmospheric neutrino flux uncertainties are • the normalization (N 0 ) which we assume to be unconstrained; • the tilt of the energy spectrum, which is parameterized by including a factor (E ν /E 0 ) γ with E 0 = 1 TeV, a 5% error on the power law index γ and a central value γ = 0; • the ratio between the pion and the kaon decays contribution to the flux (R π/K ) with a 10% error; • the ratio between the neutrino and the anti-neutrino flux (φ ν /φν) with a 5% error.
The uncertainties associated with the detector response and the ice properties, which are provided by the collaboration in data sets using the same grid as the effective area, are: • the efficiency of IceCube Digital Optical Modules, where as nominal value we have used the table corresponding to 99% efficiency, and as 1σ deviation we have used the table corresponding to 95% efficiency; • the photon scattering in the ice, where the 1σ deviation is defined from the table corresponding to a 10% increase with respect to the nominal response; • the photon absorption in the ice, where the 1σ deviation is defined as a 10% increase in the absorption rate with respect to the nominal response; • the azimuthal anisotropy in the scattering length due to the dust grain shear; here the 1σ deviation is obtained from the data set denoted 'SPICELEA ice model'; • the optical properties of the ice column surrounding each string, where the 1σ deviation is obtained from the data set labelled 'SPICEMIE ice model' which does not include hole ice effects.
For each point in the parameter space the χ 2 [φ atm ] value corresponding to the assumed flux model is calculated from the theoretical predictions and the experimental values by means of a log-likelihood function. The final χ 2 for such point is then chosen by minimizing over all the seven flux models provided by the IceCube collaboration.