All-flavor constraints on nonstandard neutrino interactions and generalized matter potential with three years of IceCube DeepCore data

We report constraints on nonstandard neutrino interactions (NSI) from the observation of atmospheric neutrinos with IceCube, limiting all individual coupling strengths from a single dataset. Furthermore, IceCube is the first experiment to constrain flavor-violating and nonuniversal couplings simultaneously. Hypothetical NSI are generically expected to arise due to the exchange of a new heavy mediator particle. Neutrinos propagating in matter scatter off fermions in the forward direction with negligible momentum transfer. Hence the study of the matter effect on neutrinos propagating in the Earth is sensitive to NSI independently of the energy scale of new physics. We present constraints on NSI obtained with an all-flavor event sample of atmospheric neutrinos based on three years of IceCube DeepCore data. The analysis uses neutrinos arriving from all directions, with reconstructed energies between 5.6 GeV and 100 GeV. We report constraints on the individual NSI coupling strengths considered singly, allowing for complex phases in the case of flavor-violating couplings. This demonstrates that IceCube is sensitive to the full NSI flavor structure at a level competitive with limits from the global analysis of all other experiments. In addition, we investigate a generalized matter potential, whose overall scale and flavor structure are also constrained.


I. INTRODUCTION
The collective evidence for flavor transitions of neutrinos propagating in vacuum and various types of matter conclusively demonstrates that at least two of the three known active neutrinos have mass [1][2][3][4]. Since these transitions typically exhibit an oscillatory pattern in the ratio of neutrino propagation distance to energy, they are also referred to as "neutrino oscillations" [5,6]. Since the neutrino mass scale is many orders of magnitude smaller than that of charged fermions, qualitatively different mechanisms for generating neutrino masses than in the Standard Model (SM) have been proposed [7].
In ordinary matter, f; f 0 ¼ e, u, d indicate charged leptons and quarks of the first generation. The coefficients ϵ fC αβ and ϵ ff 0 C αβ are the effective NC and CC NSI coupling strengths, respectively, normalized to Fermi's coupling constant G F [25]. Hence, the SM is recovered in the limit ϵ αβ → 0. The only NSI coupling strengths relevant to neutrinos propagating in matter with negligible incoherent interactions are given by NSI couplings with α ≠ β represent new sources of flavor violation (FV), whereas those with α ¼ β accommodate new flavor-diagonal interactions, which could give rise to flavor nonuniversality (NU). What sets neutrino oscillation experiments apart from other experiments is their unique capability to probe BSM scenarios responsible for NSI independently of the new physics energy scale Λ [19]. Detailed global analyses of available neutrino oscillation data (e.g., [26][27][28][29]) allowing for NSI have so far shown no statistically significant evidence for BSM interactions and have thus been used to place limits on NSI in a model-independent manner [30]. Coupling strengths up to ϵ ∼ Oð0.1Þ at a 90% confidence level continue to be allowed.
In this paper, we present a new search for NC NSI using atmospheric neutrinos 1 interacting in the IceCube DeepCore detector. Matter effects are expected for these neutrino trajectories, since their oscillation baselines range up to the order of the diameter of the Earth. Compared to our previous study [31], the present analysis is based on an extended event selection that includes neutrinos of all flavors, with reconstructed energies reaching up to 100 GeV. To obtain a high purity sample, the overwhelming background of atmospheric muons is reduced by approximately eight orders of magnitude through a series of containment and quality selection criteria [32]. Furthermore, whereas our earlier analysis constrained real-valued flavor-violating NSI in the μ-τ sector via the disappearance of atmospheric muon neutrinos, we now constrain multiple, potentially complex-valued, NSI couplings, each through its simultaneous effects in all oscillation channels. In addition, we test a more general NSI flavor structure within a CP-conserving framework, proceeding in close analogy to the analysis of atmospheric and long-baseline accelerator neutrino experiments in recent global NSI fits [30,33,34]. In a rigorous statistical approach, NSI hypotheses are tested by comparing Monte Carlo (MC) expectation to observation.

II. NEUTRINO FLAVOR TRANSITIONS IN EARTH MATTER WITH NONSTANDARD INTERACTIONS A. Evolution equation
While the flavor evolution of ultrarelativistic neutrinos in vacuum depends solely on neutrino energy, masssquared differences and the leptonic mixing matrix [5,6,35,36], evolution in matter introduces a potential position dependence [37,38]. The energy-independent interaction Hamiltonian 2 is identified with the matrix of effective neutrino potentials for coherent forward scattering in the flavor basis. For SM interactions in an unpolarized medium, the interaction Hamiltonian is with the standard matter potential V CC ðxÞ ¼ ffiffi ffi 2 p G F N e ðxÞ [39,40], where N e ðxÞ is the local electron number density. 3 This potential is responsible for the matter effects on neutrino propagation in the Sun and the Earth: the Mikheyev-Smirnov-Wolfenstein (MSW) effect and resonance [20,42] and, in the case of the Earth, parametric enhancement [43][44][45]. See, e.g., [46] for a review.
NSI contributions lead to a straightforward generalization of Eq. (4). Accounting for the hermiticity of the Hamiltonian, H mat can be well approximated by where a term proportional to ϵ ⊕ μμ · 1 was subtracted to reduce the dimensionality without observable consequences. This interaction Hamiltonian makes use of the constant effective NSI couplings, to electrons, protons and neutrons (e, p, n) in Earth matter, the latter including the nearly constant relative 1 We denote both neutrinos and antineutrinos as "neutrinos," unless a distinction is necessary. 2 While for neutrinos, the overall evolution is governed by the sum of the vacuum and matter Hamiltonian H ν ðxÞ ¼ ½H vac þ H mat ðxÞ, the evolution of antineutrinos follows HνðxÞ ¼ ½H vac − H mat ðxÞ Ã [33]. 3 The corresponding Hamiltonian can be cast into the necessary NC form by means of a Fierz transformation [41]. neutron-to-electron number density of the Earth, Y ⊕ n ≡ hN n ðxÞ=N e ðxÞi ≈ 1.051 [30]. The Hamiltonian is described by eight real NSI parameters (five amplitudes and three phases): This will in the following be referred to as "standard parametrization."

B. Generalized matter potential
In addition to the above description of NSI, our analysis is carried out in the alternative parametrization [33,47] with the three matrices on the right side defined as [30] D mat ðxÞ ¼ V CC ðxÞdiagðϵ ⊕ ; ϵ 0 ⊕ ; 0Þ; ð10Þ Here R 12 ðφ 12 Þ and R 13 ðφ 13 Þ correspond to real rotations through the angles φ 12 and φ 13 in the 1-2 and 1-3 planes, respectively, whereasR 23 ðφ 23 ; δ NS Þ denotes a complex rotation through the angle φ 23 and the phase δ NS . Since IceCube DeepCore is sensitive mainly to muon neutrino disappearance and existing data from atmospheric neutrino experiments has little sensitivity to CP-violating effects [30], the dimensionality of this parametrization can be reduced while approximately retaining model independence. We set ϵ 0 ⊕ ¼ 0 [48,49], set the phases α 1;2 ¼ 0, and disregard φ 23 and δ NS as unphysical [33,34] (see the Appendix A for a complete justification). As a result, H mat is real-valued and has three free parameters: Any given point in the three-dimensional parameter space of this "generalized matter potential" (GMP) uniquely corresponds to a point in the standard parametrization described in Sec. II A. When the vacuum Hamiltonian H vac is included in the CP-conserving framework by setting δ CP ¼ 0, we can retain the usual minimal parameter ranges for the standard Pontecorvo-Maki-Nakagawa-Sakata (PMNS) [5,35] mixing parameters and neutrino masssquared differences [50] by choosing the ranges of the matter potential rotation angles as −π=2 ≤ φ ij ≤ π=2.

C. NSI effects on the oscillation probability
In our calculation we assume an atmospheric neutrino production height of h ¼ 20 km above the surface [51]. The zenith angle ϑ then geometrically fixes the oscillation baseline [52] ranging from "upgoing," Earth-crossing (cos ϑ ¼ −1, d ≈ 1.3e4 km) trajectories to "downgoing" We approximate the Earth's matter density profile using twelve concentric uniform-density layers adopted from the preliminary reference Earth model (PREM) [53], with matter densities between about 3 g=cm 3 and 13 g=cm 3 . We take the relative electron-to-nucleon number density Y c e ¼ 0.466 for the Earth's inner and outer core; for the mantle we choose Y m e ¼ 0.496. The (nominal) values for the PMNS mixing parameters and the neutrino mass-squared differences are taken from a global fit to neutrino oscillation data [54,55], except for δ CP ¼ 0, to which the analysis is insensitive.
In Fig. 1 the oscillation probability P μτ is shown for three different NSI parameters as a function of the neutrino energy 2 GeV ≤ E ν ≤ 1000 GeV for an inclined trajectory that only crosses the Earth's mantle. The chosen zenith angle of cosðϑÞ ¼ −0.75 corresponds to a baseline L ≈ 9.6e3 km. We show the corresponding standard interactions (SI) oscillation probability as a reference in each figure. Approximations employed in the discussions below are just for illustrative purposes. All oscillation probabilities underlying the analysis in this paper are obtained by solving the full three-neutrino evolution equation [56].
The top panel of Fig. 1 shows the oscillation probability P μτ that results from varying the GMP parameter ϵ ⊕ while restricting the matter potential to the ee matrix element, i.e., φ ij ¼ 0. All of the cases correspond to a rescaling of the SM matter potential by the factor . The center and bottom panels of Fig. 1 show P μτ when the only source of NSI is ϵ ⊕ ττ − ϵ ⊕ μμ and ϵ ⊕ μτ , respectively, in the latter assuming a real-valued NSI coupling for simplicity. For the more general case of complex-valued ϵ ⊕ μτ , cf. Eq. (8), the value of the complex phase δ μτ affects the impact of the magnitude jϵ ⊕ μτ j on the oscillation probabilities in the μ-τ sector. For example, their leading-order perturbative expansions [57] reveal that a purely imaginary coupling ϵ ⊕ μτ ¼ ijϵ ⊕ μτ j (corresponding to δ μτ ¼ 90°; 270°) is expected to result in less sensitivity at the probability level.
All oscillation probabilities P αβ that are not given in Fig. 1 show similar or subdominant effects and can be found in Appendix B alongside a more detailed phenomenological discussion.

III. EVENT SELECTION WITH ICECUBE DEEPCORE
The in-ice array of the IceCube Neutrino Observatory, located at the Geographic South Pole, consists of 5160 individual photosensors. These form a cubic-kilometervolume detector for Cherenkov emission of charged particles propagating through the Antarctic ice shield [58], allowing for the detection of neutrino interactions. The subarray at the center of IceCube, DeepCore, has an effective mass of approximately 10 Mton and is instrumented approximately five times more densely with respect to the standard IceCube array in order to lower IceCube's neutrino detection energy threshold to a few GeV [59]. In event topologies, we differentiate between "tracklike" extended light depositions along muon trajectories caused by ν μ CC events and other, "cascadelike" events that predominantly consist of electromagnetic and hadronic showers.
The sample of events used in this work was collected in IceCube-DeepCore between April 2012 and May 2015. The event selection criteria only differ at the final selection level from those previously presented in [60,32] (sample "B" therein). Starting from triggered events passing the DeepCore online filter [59], the selection applies coincidence and containment criteria, using the surrounding IceCube modules as active veto. This reduces both the rates of atmospheric μ AE and noise events by approximately eight orders of magnitude, leading to a sample with a purity of approximately 95% in atmospheric neutrinos and antineutrinos. To enhance sensitivity to NSI effects, we do not impose containment of the reconstructed stopping position of the event, and we keep all events whose reconstructed energy lies below 100 GeV while adopting the same lower bound of 5.6 GeV as previous analyses. All reconstructed zenith angles are accepted. We observe 47855 events in the sample, corresponding to an increase of ∼15% compared to sample "B." We employ identical methods as therein for obtaining the expected event distribution in all observable parameters from simulation, including detailed models of photon generation and propagation. The expected sample composition according to reaction channels at the best fit point within one of our NSI hypotheses is shown in Table I, split up into the "cascadelike" and "tracklike" morphological categories that also constitute a binning dimension of the event histograms (cf. Sec. IV B). For each event type, neutrinos dominate in the sample by a factor of about two to three times over antineutrinos, predominantly due to their larger cross sections [51]. Furthermore, muon neutrinos contribute a significantly higher fraction of events (CC: ∼60%) than do electron neutrinos (CC: ∼23%) or, in particular, tau neutrinos (CC: ∼3%). Separation into "tracklike" and "cascadelike" events is based on the likelihood ratio obtained by reconstructing each event under two hypotheses: That of a cascade and track, expected for ν μ CC interactions, and that of a single cascade, expected for all other interactions. The distribution of selected events in energy, direction and morphological category is shown in Fig. 2 (Sec. V).
Apart from the SI case (ϵ ⊕ ¼ 1), the no interactions case (vacuum, ϵ ⊕ ¼ 0) is highlighted as a dashed green line. The blue lines denoting negative values are mostly covered by the dark red lines. Center panel:

A. NSI hypotheses
The IceCube DeepCore event sample introduced in Sec. III is interpreted assuming Schrödinger-like evolution of three active neutrinos and six different matter Hamiltonians. Each represents a distinct NSI hypothesis, as summarized in Table II.
The five phenomenological NSI parameters defined in Eqs. (7) and (8) are assumed to be nonzero "one-by-one" with the four remaining parameters fixed to zero in each case. Here we follow the convention in the literature of constraining NSI coupling strengths by allowing only one to be nonzero at a time [19]. This method is necessarily model dependent; the most generally applicable constraints result from accounting for the correlations between all couplings. These correlations can lead to partial cancellations and thereby to weakened constraints compared to those resulting from assuming one coupling at a time [48]. Nevertheless, we take this approach in the first part of this NSI search, not least because there are several theoretical NSI models that accommodate the possibility of the existence of only a single or a small number of sizable coupling strengths relevant to neutrino propagation [17,61,62]. Unlike this work, none of the previous analyses using IceCube data [31,[63][64][65] have performed measurements of e-μ nonuniversality or of complex couplings.
Testing the generalized matter potential in Eq. (9) with three nonzero parameters ϵ ⊕ , φ 12 , and φ 13 has a reduced model dependence compared to the one-by-one fits.
Moreover, probing all dimensions of this parameter space simultaneously is computationally feasible within our frequentist statistical framework due to the reduced dimensionality with respect to the parametrization in Eqs. (7) and (8).

B. Statistical approach
We perform a χ 2 fit to a histogram of the observed events, binned in the reconstructed ν ð−Þ energy E reco , cosinezenith cosðϑ reco Þ, and event type. For E reco , we employ eight bins covering the range from 10 0.75 ¼ 5.6 GeV to 10 1.75 ¼ 56.2 GeV that are uniformly spaced in log 10 ðE reco =GeVÞ, extended by one bin reaching up to 100 GeV. For cosðϑ reco Þ, we divide the range from −1 to 1 into eight uniformly spaced bins. The third histogram dimension is divided into the two flavor classification categories introduced in Sec. III, namely cascadelike and tracklike events. In total, there are N bins Each hypothesis from Sec. IVA is fit to the threedimensional event histogram through the minimization of a modified Pearson's χ 2 function defined as in [32,60] (see Appendix C for more detail). All fits are performed by finding the global χ 2 minimum χ 2 min ≡ min χ 2 in the multidimensional space of NSI and nuisance parameters. The d ¼ 1, 2, 3-dimensional space defined by the respective NSI parameters 4 is furthermore mapped-out using a dense grid of the same dimension, consisting of N g points fgg i¼1;…;N g ∈ C d . From the difference between the χ 2 values resulting from minimizing with NSI parameter values fixed to the single grid points, fχ 2 min ðg i Þg i¼1;…;N g , and the global χ 2 min , Δχ 2 profiles are obtained.
With d ¼ 1 for flavor-diagonal parameters, d ¼ 2 for flavorviolating complex parameters, and d ¼ 3 for the GMP parametrization.
The sampling grids for all six hypotheses are specified in Table II. In each case, the number of points N g is of the order of 10 2d .
Due to the computational infeasibility 5 of a Feldman-Cousins approach [66], we derive d-dimensional frequentist confidence regions by applying Wilks' theorem to a given Δχ 2 profile, i.e., by assuming that it behaves as a χ 2 distribution with d degrees of freedom [67]. In the case d ¼ 1, these confidence regions correspond to confidence intervals on the sampled NSI parameters. When d ¼ 2 or d ¼ 3, we determine the confidence regions and intervals in all d ¼ 1 and d ¼ 2 parameter subsets from the projections of the original, higher-dimensional Δχ 2 profile.
In order to prevent the lower-dimensional projections from getting biased due to the discrete nature of the samples in the NSI parameters to be optimized, the following routine is employed. For each point in the NSI parameters onto which the high-dimensional Δχ 2 profile is to be projected, we search for local minima on the (one-or two-dimensional) grid spanning the space of NSI parameters that have to be optimized. Each local minimum that is detected is used as a seed for an additional local minimization process. The best fit among the set of minimization outcomes is recorded and employed in the projection.

C. Nuisance parameters
A total of 15 nuisance parameters are optimized in addition to each considered set of NSI parameters. This implies that the χ 2 is a function of between 16 and 18 fit parameters, depending on the fit hypothesis. Table III gives a list of the nuisance parameters found to be relevant throughout MC studies, grouped according to their origin. Each parameter is specified together with its prior constraints, where applicable, as well as its allowed fit range. Both the choice of prior and fit range include our understanding of the behavior of the respective parameter. In addition, the fit ranges are restricted to avoid unphysical parameter space.
We account for seven uncertainties related to the intrinsic flux of atmospheric neutrinos and their detection cross sections, where the unconstrained effective livetime represents several uncertainties related to the overall normalization of the observed event count.
Of the six vacuum Hamiltonian parameters, we only let θ 23 and Δm 2 32 vary, without imposing any prior constraints. The remaining parameters have small impact on the event sample under study and are fixed to θ 12 ¼ 33.62°, [54,55], and δ CP ¼ 0. The detector related uncertainties include optical properties of the deep glacial ice and the photosensors' efficiency of detecting Cherenkov photons-both overall and depending on their angle of incidence.
The normalization of the atmospheric muon background distribution, given as a fraction of the total size of the event sample, is also included as an unconstrained nuisance parameter.
A more detailed interpretation of all nuisance parameters can be found in [32]. This also includes nuisance parameters that were found to be negligible, such as the upwardgoing vs horizontal flux of electron neutrinos and local ice properties. Table IV gives an overview of the outcomes of the six separate NSI fits discussed in Sec. IVA. The outcome of fitting SI is shown in addition in order to set the null hypothesis which is nested within all NSI hypotheses. All fits are performed within the parameter space of the normal ordering, i.e., Δm 2 32 > 0. Depending on the NSI hypothesis under consideration in the respective fit, this choice does not a priori result in any loss of generality of the derived NSI constraints. We return to the mass ordering question below in the context of each set of fit results.

V. RESULTS
All outcomes are characterized by a goodness of fit in the range of 19% to 22%. The goodness of a given fit hypothesis is not determined from the Δχ 2 value in TABLE III. Nuisance parameters employed by all NSI fits, as well as their associated Gaussian priors and fit ranges. For a given parameter with a prior, the range is specified as a number of standard deviations (σ) from the prior's nominal value. See text for the interpretation of all parameters.

Parameter
Prior Fit range Oscillation: Optical efficiency, head-on (a.u.) Atmospheric muons: 35] Table IV but by comparing the observed value of χ 2 min to the test-statistic distribution resulting from fitting the same hypothesis to a large number of statistically independent pseudoexperiments generated assuming SI. No nuisance parameter with an external constraint is found to experience a statistical pull from its best fit value beyond 1.1σ, no matter which fit hypothesis is chosen (see Appendix D for more detail).
For computational reasons, the compatibility of the SI and NSI best fit hypothesis is tested using Δχ 2 instead of pseudoexperiments. In all cases, the best fit SI hypothesis is statistically compatible with the best fit NSI hypothesis: the strongest disfavoring of the SI hypothesis is observed for the assumption of e-μ nonuniversality, at approximately p ¼ 0.3. Our best fit values of the vacuum Hamiltonian parameters Δm 2 32 and θ 23 under the SI hypothesis are compatible with the constraints found in the dedicated IceCube DeepCore analyses of Refs. [32,60]. In addition, the best fit values of Δm 2 32 and θ 23 under the various NSI hypotheses are within 2.5% and 4%, respectively, of the values obtained assuming the SI hypothesis.
In Fig. 2, we show E reco slices of the observed event counts as a function of cosðϑ reco Þ ≤ 0 for the two event Histograms of observed cascadelike events (top row) and tracklike events (bottom row) as a function of cosðϑ reco Þ for different slices in E reco (indicated at the top of each panel), together with the MC expectation under the generalized matter potential fit outcome, labeled as "best fit (GMP)." For display purposes, the eight lowest reconstructed energy bins have been merged into four, and only the upgoing region cosðϑ reco Þ ≤ 0 is shown, where the largest NSI effects are expected. Also shown are the expected event distributions for one particular μ-τ nonuniversality realization (ϵ ⊕ ττ − ϵ ⊕ μμ ¼ 0.10), one μ-τ flavor-violation realization (ϵ ⊕ μτ ¼ 0.050), and one e-μ flavorviolation realization (ϵ ⊕ eμ ¼ −0.30). In each of these three example NSI scenarios, all nuisance parameters are set to their respective global best fit values within the corresponding NSI parameter space. classes (rows). In the figure, we have condensed the eight lowest energy bins into four slices, each of which covers two energy bins of the original binning used in the analysis. We also show the best fit of the generalized matter potential hypothesis, as well as three signal hypotheses with nonzero For the latter three, all nuisance parameters are set to the values obtained by the respective global best fit, with the NSI coupling strengths given in Table IV (see Appendix D for the detailed nuisance  parameter values). Thus, the induced event count differences follow solely from choosing different NSI parameter values compared to those that fit the data best; the event distributions at all three best fit points would be barely distinguishable by eye from the fit of the generalized matter potential in Fig. 2. These differences in count are what is observable of the imprints of NSI on the oscillation probability after superimposing the expected event distributions in E reco , cosðϑ reco Þ, and event classification of both the atmospheric μ AE background and the effective-area weighted oscillated fluxes of neutrinos and antineutrinos of all flavors; see the thirteen sample components in Table I. A. One-by-one fits 1. Flavor-nonuniversal NSI Figure 3 shows observed Δχ 2 profiles as a function of the two differences of the flavor-diagonal NSI coupling strengths, namely ϵ ⊕ ee − ϵ ⊕ μμ and ϵ ⊕ ττ − ϵ ⊕ μμ . In each case, all other coupling strengths are fixed to zero.
The shaded bands give the experimental sensitivity by showing the symmetrical central 68.3% and 90% confidence intervals of the Δχ 2 distributions obtained in fits to pseudoexperiments for the generation of which SI were assumed.
Horizontal dashed lines denote the 68.3th, 90th, and 99.7th percentiles of a χ 2 distribution with one degree of freedom. The vanishing impact from θ 12 and δ CP causes the sign of 1 þ ϵ ⊕ ee − ϵ ⊕ μμ to be fully degenerate with the mass ordering [68]. This applies similarly to μ-τ nonuniversality and jϵ ⊕ μτ j flavor violation. We therefore do not need to test solutions within the inverted ordering explicitly. When interpreted in terms of standard matter effects, the Δχ 2 profile asymmetry about ϵ ⊕ ee − ϵ ⊕ μμ ¼ −1 suggests that the data slightly favors the normal ordering, corresponding to the point ϵ ⊕ ee − ϵ ⊕ μμ ¼ 0, over the inverted ordering, corresponding to the point ϵ ⊕ ee − ϵ ⊕ μμ ¼ −2, at the level of Δχ 2 ≈ 0.5. 6 A more detailed discussion of the profile characteristics and their causes can be found in Appendix E.
μ-τ nonuniversality.-From the right panel of Fig. 3, we find that the observed event sample is fully compatible with NSI that are μ-τ flavor universal, that is, ϵ ⊕ ττ − ϵ ⊕ μμ ¼ 0. In contrast to the observed Δχ 2 profile under the hypothesis of e-μ nonuniversality, here the test statistic keeps increasing for large values of jϵ ⊕ ττ − ϵ ⊕ μμ j, which allows for stringent FIG. 3. Observed Δχ 2 profiles as a function of the effective NSI flavor nonuniversality parameters ϵ ⊕ ee − ϵ ⊕ μμ (left) and ϵ ⊕ ττ − ϵ ⊕ μμ (right), together with the central 68.3% and 90% confidence intervals of the experimental sensitivity shown as shaded bands. See text for details. 6 See [69] for a statistically rigorous study of the neutrino mass ordering in the absence of NSI with two related IceCube DeepCore event samples.
constraints with values of ϵ ⊕ ττ − ϵ ⊕ μμ outside the interval ½−0.041; 0.042 excluded by the data at 90% CL. The sensitivity to this type of NSI stems almost exclusively from its impact on the ν μ andν μ survival probabilities, cf. Sec. II C. We find that the summation over neutrinos and antineutrinos in general does not lead to significant cancellations of the respective NSI signatures.

Flavor-violating NSI
The central panel of the each of the three plots in Fig. 4 shows the observed 90% CL contour (Δχ 2 ≈ 4.61) in the NSI magnitude and complex phase from the fit of a given flavor-violating NSI coupling strength. The projection of the two-dimensional Δχ 2 profile onto the magnitude of the coupling strength is depicted on top, and that onto the complex phase on the right. Lines and shaded bands have the same meanings as in Fig. 3. Note that the SI case is located at the origin. The appropriate entry for Δχ 2 SI in Table IV provides the maximal projected Δχ 2 at which any value of the complex phase is disfavored, due to the projection method and the vanishing amplitude rendering the complex phase unphysical. Fig. 4 (left), the magnitude of e-μ flavor-violating NSI is compatible with zero (or SI) at a significance level of approximately p ¼ 0.3 and an upper bound of jϵ ⊕ eμ j ≤ 0.146 (90% CL) is obtained when the full range 0°≤ δ eμ ≤ 360°is considered. A stronger constraint on the magnitude follows when δ eμ is only allowed to take more disfavored values of 160°AE90°. In this limited range of δ eμ , the NSI magnitude that best fits the data is zero. This explains the "plateau" in the projection onto δ eμ with Δχ 2 ¼ Δχ 2 SI ≈ 1.1 (compare Table IV). A somewhat stronger exclusion of real negative values of ϵ ⊕ eμ (δ eμ ≈ 180°) with respect to the expectation from pseudoexperiments for SI is observed. Within the inverted ordering parameter space, the 90% CL exclusion contour shifts to larger values of jϵ ⊕ eμ j by approximately 10%; no change is observed for the one-dimensional jϵ ⊕ eμ j interval allowed at 90% when the inverted ordering is adopted. Furthermore, the one-dimensional allowed intervals (at any CL) for both parameters obtained under the assumption of Δm 2 32 > 0 also apply to the scenario in which the mass ordering is considered as a nuisance parameter. For the observed e-τ flavor-violation constraints switching from the normal ordering to the inverted ordering parameter space similarly has negligible impact. e-τ flavor violation.-Compared to e-μ flavor violation, we find both qualitatively and quantitatively similar bounds on e-τ flavor violation; see Fig. 4 (middle). The 90% CL upper bound on the NSI magnitude from optimizing over 0°≤ δ eτ ≤ 360°is slightly larger, jϵ ⊕ eτ j ≤ 0.173, and the best fit is well compatible with the SI hypothesis. For δ eτ values in an approximately 90°range around 200°, a somewhat more constraining bound results for the magnitude. The NSI hypotheses corresponding to this limited range of δ eτ best fit the data when jϵ ⊕ eτ j ¼ 0, leading to a plateau in the projection onto δ eτ with Δχ 2 ¼ Δχ 2 SI ≈ 0.5.

e-μ flavor violation.-From
μ-τ flavor violation.-The right side of Fig. 4 suggests that the selected event sample has significantly better sensitivity to μ-τ flavor violation than to flavor violation in the electron sector considered in Sec. II C. When the full δ μτ range is allowed, jϵ ⊕ μτ j ≤ 0.0232 at 90% CL. We find the strongest bounds on the NSI magnitude for real NSI, i.e., for δ μτ ¼ 0°and δ μτ ¼ 180°. The data sample shows almost vanishing sensitivity to δ μτ . This observation is in agreement with Table IV,  maximum of the Δχ 2 projection onto δ μτ . Hypotheses with δ μτ ≈ 125°and δ μτ ≈ 235°result in the weakest constraints on the magnitude, not those with δ μτ ¼ 90°or δ μτ ¼ 270°, for which the contribution of the magnitude jϵ ⊕ μτ j to the oscillation probabilities in the μ-τ sector is reduced, see Sec. II C. Pseudoexperiments suggest that such a deviation is characteristic of considering the joint ν μ CC andν μ CC event distribution. Indeed, the strength of cancellations between the effects of jϵ ⊕ μτ j on neutrino and antineutrino channels in the medium-to high-energy regime depends on the value of δ μτ .

Summary and experiment comparison
Table V compiles a summary of the constraints (at 90% CL) placed by this analysis on the NSI flavor nonuniversality and flavor-violation parameters. Both SI and flavor-universal NSI (in the case of flavor-diagonal couplings) are compatible with each best fit NSI hypothesis. None of the complex phases are constrained at 90% CL.
For comparison with existing measurements, in Fig. 5 we restrict the flavor-violating coupling strengths to the real plane, defined by δ αβ ¼ 0°; 180°, and show the 90% CL intervals for the real-valued signed coupling strengths ϵ ⊕ αβ . The lower limits (δ αβ ¼ 180°) on ϵ ⊕ eμ and ϵ ⊕ eτ are stronger than their upper limits (δ αβ ¼ 0°). The latter reproduce the constraints on the NSI magnitudes jϵ ⊕ eμ j and jϵ ⊕ eτ j that are found under the hypotheses of complex coupling strengths in Table V. In the case of ϵ ⊕ μτ , the upper limit is slightly stronger than the lower limit, −0.0165 ≤ ϵ ⊕ μτ ≤ 0.0130. This range is fully compatible with, and smaller than, that reported by our previous study, −0.020 ≤ ϵ ⊕ μτ ≤ 0.024 [31]. 7 Neither of the magnitudes of the upper or the lower limit reproduces the limit on jϵ ⊕ μτ j in Table V because the sensitivity of the event sample to jϵ ⊕ μτ j is weakest for a complex coupling strength (cf. Fig. 4).
Data from a number of other neutrino experiments has been used to set limits on the NSI coupling strengths ϵ uV;dV αβ , which we have rescaled for consistency with the definition of the effective coupling strengths for Earth matter, Eq. (6). Figure 5 contains results reporting one-dimensional 90% CL intervals, almost all of which are based on oneby-one fits similar to those discussed in Sec. VA. Among these are limits on ϵ ⊕ ττ − ϵ ⊕ μμ and ϵ ⊕ μτ (with correlations) obtained from atmospheric neutrino data collected by Super-Kamiokande [70], as well as limits on ϵ ⊕ μτ from long-baseline accelerator ν ð−Þ μ disappearance data from MINOS [71] and high-energy atmospheric ν ð−Þ μ disappearance data from IceCube [64] (labeled "IC 2017"), respectively. Furthermore, we show the limits on flavor-violating coupling strengths reported by an analysis of the published timing (or flavor) data from coherent elastic neutrinonucleus scattering (CEνNS) at COHERENT [72]. Here, the assumed underlying NSI model based on the exchange of a Z 0 mediator with M Z 0 ∼ O ð10 MeVÞ dictates ϵ u αβ ¼ ϵ d αβ , so that no cancellations between NSI with different quark flavors occur (see [73] for a comprehensive analysis).
While CEνNS only yields constraints that are valid for a new physics energy scale above Oð10 MeVÞ, it is sensitive TABLE V. Summary of 90% CL constraints on NSI nonuniversality and flavor-violation parameters obtained by the one-byone fits in this study, as well as on the parameters of the generalized matter potential, whose fit is discussed in Sec. V B. Δm 2 32 > 0 is assumed everywhere, but does not introduce a loss of generality (see text for details).

Hypothesis
Parameter Allowed interval (90% CL) FIG. 5. Summary of the one-by-one constraints at 90% CL on real NSI nonuniversality and flavor-violation parameters obtained in this study (labeled as "IC DC 2021") compared to previous limits [30,31,64,65,[70][71][72]. Constraints on the magnitudes of complex NSI parameters are given for the respective phase restricted to δ αβ ¼ 0°, 180°. See text for details. 7 After translating from NSI with down quarks to the effective coupling strengths in Eq. (6).
to the individual flavor-diagonal coupling strengths ϵ uV;dV ee and ϵ uV;dV μμ (not depicted in Fig. 5)-in contrast to neutrino oscillation experiments. Similarly, our results are not directly comparable to NSI limits set in collider experiments as these commonly depend strongly on the underlying model and new physics energy scale [74]. Figure 5 additionally allows gauging the impact of the increased event statistics and the inclusion of higher-energy events in our sample compared to a study with public IceCube DeepCore data in [65], labeled "IC DC 2020 (public)." The widths of the 90% CL intervals are smaller by between ∼25% (for ϵ ⊕ ττ − ϵ ⊕ μμ , ϵ ⊕ eτ , ϵ ⊕ μτ ) and ∼50% (ϵ ⊕ eμ ).
We also display the limits derived in a combined analysis of global neutrino oscillation datasets with negligible sensitivity to CP-violating effects [30] ("global 2018"). The global analysis only assumes the coupling strengths ϵ uV;dV αβ to be nonzero; their flavor structure is taken to be independent of the quark type. The fact that correlations between NSI couplings with different flavor indices are fully taken into account explains why these constraints are no more stringent than those found in this study.

B. Generalized matter potential
Our final fit to data employs the generalized matter potential that is characterized by the three intrinsic NSI FIG. 6. Observed 68.3%, 90%, and 99.7% confidence regions for parameters ϵ ⊕ , φ 12 , and φ 13 , together with each parameter's projected one-dimensional Δχ 2 profile. The color in each of the three large panels encodes the local value of the projected twodimensional Δχ 2 profile. The best fit point for each pair of parameters is indicated by a cross. The SI/flavor-universal NSI hypothesis, indicated by the dash-dotted lines, is located at ϵ ⊕ ¼ 1, φ 12 ¼ 0, φ 13 ¼ 0. See text for details. This is the first time the GMP overall scale and flavor structure are constrained simultaneously using IceCube DeepCore data. parameters ðϵ ⊕ ; φ 12 ; φ 13 Þ. Figure 6 shows the resulting constraints, by means of the projected one-and twodimensional Δχ 2 profiles. In terms of the five standard NSI parameters, the indicated best fit, also given in Table IV, corresponds to It is weakly favored over the hypothesis of SI (or flavoruniversal NSI) by Δχ 2 ¼ 2.2, corresponding to p ¼ 0.1, cf. Table IV. This difference cannot be directly derived from any of the projections in Fig. 6, as none of them explicitly show the corresponding grid points ðϵ ⊕ ¼ AE1; The one-dimensional projections yield the following 90% confidence intervals (optimized over the two remaining matter potential parameters and all nuisance parameters in each case): −9°≤ φ 12 ≤ 8°, −14°≤ φ 13 ≤ 9°, and the union of intervals ½−1.2; −0.3 ∪ ½0.2; 1.4 for ϵ ⊕ . The fact that φ 12 and φ 13 are allowed to vary does not have a significant weakening effect on the bounds on ϵ ⊕ at 90% CL, nor does it change the overall shape of its Δχ 2 profile (compare ϵ ⊕ ee − ϵ ⊕ μμ in Fig. 3). The two-dimensional Δχ 2 projection onto ðϵ ⊕ ; φ 12 Þ demonstrates that jϵ ⊕ j ≥ 0.05 is excluded at a significance greater than 99.7% when jφ 12 j ≥ 10°, for any value of φ 13 . Similarly, the projection onto ðϵ ⊕ ; φ 13 Þ implies that jϵ ⊕ j ≥ 0.1 is excluded at a significance greater than 99.7% when jφ 13 j ≥ 20°, for any value of φ 12 . For smaller values of jφ 12 j and jφ 13 j, no 99.7% bound on ϵ ⊕ is obtained.
Conversely, in the projection onto ðφ 12 ; φ 13 Þ, constraints can be set at 90% CL. However, the maximal significance of excluding any particular pair of values of the matter rotation angles cannot exceed the Δχ 2 value of the vacuum hypothesis, which renders both parameters unphysical. Combined with the lacking bound on ϵ ⊕ at the 99.7% CL this results in the "crosslike" shape formed by the corresponding contours in the two upper Δχ 2 projections in Fig. 6.
Finally, we point out that these constraints do not suffer from a loss of generality due to the normal ordering assumption in the fit, for the reasons given in Sec. II A.

VI. CONCLUSION
We have presented a comprehensive study of nonstandard interactions in the propagation of atmospheric neutrinos observed with IceCube DeepCore within the general framework of three flavor neutrino oscillations. Instead of exclusively focusing on NSI in the μ-τ sector, as was done in our previous analysis [31], we have taken an extended approach that tests all five effective flavornonuniversal and flavor-violating NSI coupling strengths for Earth matter individually. In particular, this includes studies of NSI involving the electron flavor, which are not common targets of atmospheric neutrino experiments. All our measurements yield results that are statistically compatible with SM neutrino interactions, i.e., neutrino oscillations with standard matter effects.
The sample of 47855 events with reconstructed energies between 5.6 GeV and 100 GeV was created from three years of data taken with IceCube DeepCore and contains significant contributions from the interactions of neutrinos and antineutrinos of all flavors. One-by-one NSI parameter fits to this sample result in limits (quoted at 90% CL) of similar power with respect to existing global limits on the magnitudes of all five NSI parameters observable by atmospheric neutrino experiments. Those that apply to μ-τ nonuniversality and flavor-violation strengths are of the order of 10 −2 and are as, or more, stringent than limits obtained with other oscillation experiments or other IceCube (DeepCore) event samples. Weaker Oð1Þ constraints apply to e-μ nonuniversality, or, when reinterpreted in terms of SM interactions, to the strength of the Earth's standard matter potential.
With a separate fit we have investigated a more general flavor structure of the Earth's matter potential, in a manner similar to recent global NSI fits [30,33,34]. The adopted parametrization naturally includes NSI hypotheses that lead to cancellations of the induced matter effects in the survival probabilities of atmospheric muon neutrinos and antineutrinos. Within this framework, we have shown that the event sample allows for simultaneous constraints of the overall strength of the matter potential and its flavor structure at 90% CL, whereas no constraint emerges at 99.7% CL.
Because of the vanishing momentum transfer in the coherent forward scattering processes that generate the neutrino matter potential, our constraints apply independently of the new physics energy scale responsible for NSI. This distinguishes our measurements from those performed at experiments investigating coherent neutrino scattering, deep inelastic neutrino scattering, or at high-energy colliders.
Future versions of this analysis may profit from enhanced minimization approaches, as the computational limitations of this analysis are closely connected to the challenges of minimizing a high-dimensional parameter space with a large number of local minima. For upcoming NSI measurements with IceCube and its low-energy extension DeepCore, a significant increase in event statistics and an extended energy range compared to the analysis presented in this paper are expected. Furthermore, the imminent IceCube Upgrade [75] will increase the detection efficiency and improve the reconstruction capabilities for atmospheric neutrinos with respect to DeepCore, and lower the energy threshold to allow high-statistics measurements with ∼1 GeV atmospheric neutrinos. It will thus facilitate the determination of the overall strength of the Earth's matter potential and improve IceCube's ability to distinguish NSI from standard matter effects [76]. mass ordering degeneracy" 8 once NSI are introduced, because H mat ðxÞ → −½H mat ðxÞ Ã can be implemented by [68] In the Earth, where the effective NSI couplings have little variation along the neutrino trajectory, cf. Eq. (6), the degeneracy is almost exact. When H mat is only described by ðϵ ⊕ ; φ 12 ; φ 13 Þ, it is therefore sufficient to restrict Δm 2 31ð32Þ > 0 and test both signs of ϵ ⊕ . The two choices ðϵ ⊕ ¼ AE1; φ 12 ¼ 0; φ 13 ¼ 0Þ correspond to neutrino propagation with SI given the normal ordering ("þ") and the inverted ordering ("−"), respectively.

APPENDIX B: NSI PHENOMENOLOGY AT THE PROBABILITY LEVEL
At the GeV-scale energies considered here, all transitions involving ν e are suppressed in vacuum compared to those not involving ν e . For energies above a few GeV, ν e → ν e;μ;τ transitions are driven by the mixing angle θ 13 and the "atmospheric" mass-squared difference Δm 2 32 , with negligible corrections due to the "solar" mass-squared difference Δm 2 21 . For Δm 2 32 > 0-the baseline assumption in the example cases in this paper-SM matter effects in general lead to an enhancement of the transitions involving ν e , while a negative matter potential in general leads to their suppression. 9 In Figs. 7-11 oscillation probabilities P αβ are shown for different NSI parameters. As all neutrino flavors are considered in this study, no individual oscillation channel can be singled out a priori. However, at the energies considered in this paper, the tau neutrino fluxes generated in the atmosphere are negligible [77], resulting in a restriction to α ∈ ðe; μÞ. Also, neutrino absorption is not relevant below the TeV scale [52]. In all cases, in the absence of intrinsic CP violation (i.e., δ CP ¼ 0 or δ CP ¼ π and real NSI coupling strengths) the Earth's symmetric matter potential with respect to the midpoint of any given trajectory, VðxÞ ¼ VðL − xÞ, implies that P αβ ¼ P βα (apart from negligible short-scale corrections due to h prod ≠ −d det ) [78].
1. Nonzero single NSI parameters a. Nonzero ϵ ⊕ (rescaled SM matter potential) The oscillation probabilities P αβ shown in Fig. 7 result from varying ϵ ⊕ ∈ ½−5; 5 while restricting the matter potential to the ee entry, i.e., φ ij ¼ 0, yielding effects corresponding to a rescaling of the SM matter potential by the factor V CC ðxÞ→V 0 ðxÞ¼ϵ ⊕ V CC ðxÞ¼ð1þϵ ⊕ ee −ϵ ⊕ μμ ÞV CC ðxÞ. Typically, the ν e disappearance probability, 1 − P ee , in vacuum remains small in the limit Δm 2 21 → 0. In contrast to this, when ϵ ⊕ > 0, the resonance condition can be satisfied (given Δm 2 32 > 0 and θ 13 < π=4). In this case, the effective 1-3 mixing angle in matter becomes maximal. A complete disappearance of ν e can therefore be observed in principle at a resonance energy E R which is inversely proportional to the average value of the slowly changing rescaled matter potential along the neutrino trajectory, ϵ ⊕ hV CC i. For the Earth's mantle and SI, E SI R ≈ 6 GeV. The example trajectory in Fig. 7 is chosen such that the oscillation phase leads to a nearly complete disappearance of ν e at this energy. The transition probabilities to ν μ and ν τ are nearly identical since the 2-3 mixing is close to maximal. A complete disappearance is also observed for a value of ϵ ⊕ ≈ 3 and E R ≈ 2 GeV. Negative values of ϵ ⊕ together with Δm 2 32 > 0 do not give rise to a similar enhancement. Consequently, there are no significant transitions ν e → ν μ;τ for ϵ ⊕ < 0. Instead, the antineutrino transitionsν e → ν e;μ;τ are then subject to the matter effects detailed above. Figure 7 further demonstrates that ν e decouples from the evolution and that the transitions ν μ → ν μ;τ proceed as in vacuum for sufficiently high energy, E ν ≳ 20 GeV for the considered trajectory-irrespective of the value of ϵ ⊕ . At sufficiently low energies of a few GeV, the simple twoneutrino picture no longer applies, and rather complex corrections due to 1-3 mixing need to be taken into account. For their discussion see, for example, Refs. [46,79].
In case ϵ ⊕ ττ − ϵ ⊕ μμ is the only source of NSI-cf. Fig. 8 for variations ϵ ⊕ ττ − ϵ ⊕ μμ ∈ ½−0.20; 0.20-the nonuniversality gives rise to an effective potential in the decoupled μ-τ system that was introduced in the previous section. As detailed in [63], the 2-3 mixing in matter is modified according to the standard MSW mechanism, but with a potential V 0 ðxÞ ¼ ðϵ ⊕ ττ − ϵ ⊕ μμ ÞV CC ðxÞ. For a given sign of the nonuniversality, whether the resonance occurs in the neutrino or the antineutrino channel depends on the octant of θ 23 . Since the 2-3 mixing in vacuum is nearly maximal, the introduction of nonzero ϵ ⊕ ττ − ϵ ⊕ μμ in general leads to a reduction of the mixing. The main observable consequence of ϵ ⊕ ττ − ϵ ⊕ μμ is therefore the increased survival probability of both atmospheric ν μ 's andν μ 's across the broad range of energies at which the μ-τ system is decoupled.
In contrast to the ϵ ⊕ ττ − ϵ ⊕ μμ -only case, when ϵ ⊕ μτ is the only nonzero NSI coupling strength-cf. Fig. 9-, the offdiagonal elements V CC ðxÞϵ ⊕ðÃÞ μτ of the two-neutrino interaction Hamiltonian result in qualitatively different effects on the neutrino evolution [63]. A resonance occurs for neutrinos when ϵ ⊕ μτ < 0 and for antineutrinos when ϵ ⊕ μτ > 0, independent of the octant of θ 23 . Resonances at E R ≳ 60 GeV are observed. Since the corresponding oscillation phases are small, the survival probability P μμ becomes nearly maximally enhanced at high energies when ϵ ⊕ μτ < 0. At energies sufficiently far below the resonance, ϵ ⊕ μτ results in a shift in energy of the oscillation pattern in the μ-τ system [63]. When ϵ ⊕ μτ > 0, a shift to higher energies appears for neutrinos, and a shift to lower energies for antineutrinos; the effects are reversed for ϵ ⊕ μτ < 0. At high energies, the two-neutrino survival probability of both ν μ andν μ is reduced compared to the vacuum value [63].
c. Subdominant single NSI parameters: ϵ ⊕ eμ and ϵ ⊕ eτ Similar to ϵ ⊕ ee − ϵ ⊕ μμ , the flavor-violating couplings involving the electron flavor, ϵ ⊕ eμ and ϵ ⊕ eτ , typically are not in the focus of atmospheric neutrino studies, partly due to their weaker impact on the disappearance probabilities of ν ð−Þ μ : It has been shown perturbatively that flavorviolating couplings involving the electron flavor contribute to disappearance probabilities of ν ð−Þ μ only at second order, far away from the 1-3 MSW resonance regime [57]. They enter the oscillation probabilities involving the electron flavor at the second order, lower by one order compared to the four remaining couplings [80].
The range of oscillation probabilities induced by values of ϵ ⊕ eμ ∈ ½−0.30; 0.30 are depicted in Fig. 10. One prominent effect of this is the suppression of the trajectory dependent SI ν e resonance around E ν ≳ 6 GeV for large absolute coupling FIG. 8. Same as Fig. 7 values. Varying ϵ ⊕ eτ within the same range as ϵ ⊕ eμ results in very similar patterns given the exchange of the flavor indices μ ↔ τ [80]. Hence, only ϵ ⊕ eμ results in modifications of the atmospheric oscillation channels involving ν μ across the full range of energies. Characteristically, at the energies shown here it manifests itself in the disappearance of ν μ andν μ and the appearance of ν e andν e . In contrast, ϵ ⊕ eτ induces the conversion ν ð−Þ e ↔ ν ð−Þ τ at high energies. For detailed phenomenological and numerical discussions of the oscillation-probability impact of ϵ ⊕ eμ and ϵ ⊕ eτ at the GeV energy scale (in the context of future long-baseline and atmospheric neutrino experiments) see, e.g., Refs. [81][82][83].

Arbitrary NSI flavor structure
In the generalized matter potential parametrized by Eq. (9), ϵ ⊕ plays the role of the strength of the matter potential and drives the overall sizes of the coupling strengths. This is evident from the fact that all elements of H mat are ∝ ϵ ⊕ , cf. Appendix A.
Once H mat is allowed to take an arbitrary flavor structure, atmospheric neutrino oscillation probabilities are not, in general, treatable analytically, prompting the implementation of well motivated constraints on the parameter space to yield a point of reference for a phenomenological discussion. As discussed in Refs. [48,49,84], in specific regimes of neutrino propagation the three-neutrino evolution in the presence of NSI can be reduced to an analytically treatable effective two-neutrino system, which is rotated with respect to the flavor basis. The specific case investigated is when all NSI parameters are zero except ϵ ⊕ ee − ϵ ⊕ μμ , ϵ ⊕ ττ − ϵ ⊕ μμ , and ϵ ⊕ eτ . Here, two identical eigenvalues result in the "atmospheric parabola" relation which is able to accommodate two-flavor vacuumlike (cf. Sec. II A) ν ð−Þ μ disappearance at high energy, independent of the overall sizes of the involved NSI coupling strengths. The relations given in Appendix A demonstrate that Eq. (B1) is satisfied for φ 12 ¼ 0, in which case the flavor-violating coupling strengths involving the μ flavor are zero, ϵ ⊕ eμ ¼ ϵ ⊕ μτ ¼ 0. Hence, as a point of reference, we show the oscillation probabilities obtained for different values of φ 13 ∈ ½−20°; 20° while keeping the overall strength of the matter potential and the 1-2 matter rotation angle fixed at ϵ ⊕ ¼ 1 and φ 12 ¼ 0, respectively, in Fig. 11. In contrast to the behavior resulting from this, the case φ 13 ¼ 0 and φ 12 ≠ 0 results in high-energy ν The modified Pearson's χ 2 function used here is defined as [32,60] Here, n obs i is the observed number of events in bin i and n exp i is the combined expectation due to ν ð−Þ and background events in the same bin. The expectation n exp i depends on the values of the hypothesis parameters of interest and on the values of several nuisance parameters (cf. Sec. IV C). Its ν ð−Þ contribution is retrieved by reweighting a large sample of simulated events, with an effective livetime that exceeds that of the observed event sample by one order of magnitude. The variance of the expectation, ðσ exp i Þ 2 , is given by the sum of the variance σ 2 i;ν of the expected number of ν ð−Þ events and the variance σ 2 i;bkg of the expected number of background events in the bin (cf. Table I). Generally, the two variances on the right hand side of Eq. (C2) are found to be of similar size, except for some bins in the downgoing region with cosðϑ reco Þ > 0.5, where the uncertainty of the background expectation dominates. This is also the only region of the histogram for which the total variance of the expectation, ðσ exp i Þ 2 , can reach a similar size as the Poisson variance n exp i . Sensitivity to NSI on the contrary predominantly originates from the upgoing region, cosðϑ reco Þ < 0.
The second sum contributing to χ 2 in Eq. (C1) is taken over only those N prior nuisance parameters that are subject to external Gaussian constraints: a deviation Δp j of the jth such parameter from its nominal value is penalized according to the parameter's prior standard deviation σ p j as ðΔp j Þ 2 =σ 2 p j .

APPENDIX D: NUISANCE PARAMETER PULLS
Ten of the 15 nuisance parameters that are optimized together with each considered set of NSI parameters have a Gaussian prior associated, as was introduced in Sec. IV C. The statistical pulls on the best fit values of these nuisance parameters show little variance between the single fit hypotheses (see Table II), showing the small impact on the expected signal of the different best fit NSI parameter hypotheses. In addition, the statistical pulls on the nuisance parameter fit values are within 1.1σ in all of the fits (see Fig. 12), which is expected in case of correctly chosen nuisance parameter priors and ranges. All best fit values of nuisance parameters with no Gaussian prior are well within their allowed ranges listed in Table III. FIG. 11. Same as Fig. 7, but for different realizations of the matter rotation angle φ 13 , with −20°≤ φ 13 ≤ 20°, keeping ϵ ⊕ ¼ 0 and φ 12 ¼ 0 fixed. The blue dashed lines show the probabilities obtained for φ 13 < 0, while the red dashed lines show those obtained for φ 13 > 0. Darker shades represent larger jφ 13 j. See text for details.