Search for non-standard neutrino interactions with 10 years of ANTARES data

Non-standard interactions of neutrinos arising in many theories beyond the Standard Model can significantly alter matter effects in atmospheric neutrino propagation through the Earth. In this paper, a search for deviations from the prediction of the standard 3-flavour atmospheric neutrino oscillations using the data taken by the ANTARES neutrino telescope is presented. Ten years of atmospheric neutrino data collected from 2007 to 2016, with reconstructed energies in the range from ∼16 GeV to 100 GeV, have been analysed. A log-likelihood ratio test of the dimensionless coefficients εμτ and εττ − εμμ does not provide clear evidence of deviations from standard interactions. For normal neutrino mass ordering, the combined fit of both coefficients yields a value 1.7σ away from the null result. However, the 68% and 95% confidence level intervals for εμτ and εττ − εμμ, respectively, contain the null value. Best fit values, one standard deviation errors and bounds at the 90% confidence level for these coefficients are given for both normal and inverted mass orderings. The constraint on εμτ is among the most stringent to date and it further restrains the strength of possible non-standard interactions in the μ − τ sector.


Introduction
Neutrinos have provided the first hint of physics beyond the Standard Model (BSM) through the discovery of neutrino oscillations [1,2], which imply that at least two of the three neutrinos presently known have non-zero masses (see chapter 14 in ref. [3]).The origin and smallness of the neutrino masses compared to the rest of the Standard Model (SM) particles point to new physics at a very high energy scale [4,5] and is at present the subject of intense research.
Interactions not present in the SM are predicted by a wide variety of BSM theories.Considering the presence of new physics at a high energy scale, the effect at a lower scale can be approached in a model independent way by constructing a sum of operators built from the fields which are relevant at the low energy scale.BSM will then show up at a lower scale via nonrenormalisable operators with dimension five or higher [6,7].In the context of neutrino physics, out of all possible neutrino interactions that can be built within this effective theory approach [8], a subset called Non-Standard Interactions (NSIs) has attracted particular attention [9,10].The NSIs are quantified through dimensionless constants (ε αβ ) that appear in the fourfermion interactions expressed as where G F is the Fermi constant, P X (with X=R or L) denotes the chiral projection operators P R,L = 1 2 (1±γ 5 ), f is a first generation SM fermion (e, u or d-quarks), f belongs to the same weak doublet as f , and α and β denote the neutrino flavours: e, µ or τ .The dimensionless coefficients ε f f X αβ and ε f X αβ quantify the strength of NSIs between the neutrinos of flavour α and β and the fermion f ∈ {e, u or d} (for neutral currents) and f = f ∈ {u, d} (for charged currents).The SM scenario is recovered in the limit ε → 0. While charged current (CC) NSIs affect the production and detection processes of neutrino states at scattering experiments [11], neutral current (NC) NSIs would affect the neutrino propagation by coherent forward scattering.In this paper, we constrain NC NSIs that alter the propagation of neutrinos while travelling long distances through the Earth.
The possible existence of NSIs can be studied in a variety of experiments [10,12,13,14].In particular, atmospheric neutrinos provide an excellent opportunity since NSIs will modify the SM potential that describes the matter effects in neutrino flavour oscillations.This will give rise to additional effects that can produce deviations from the expectations for the standard oscillation phenomenon in matter [10,13,14,15,16].High-energy neutrino telescopes are suitable detectors to search for NSI-induced deviations since they are expected to increase with the distance travelled through matter and with energy.Limits on NSIs from atmospheric neutrinos have already been reported by the Super-Kamiokande Collaboration [17] or using Super-Kamiokande data [18], and by the IceCube Collaboration [19,20,21] or using IceCube data [22,23,24].
In this paper, a search for NSI-induced deviations from the expected neutrino oscillation process using 10 years of ANTARES data [25] is presented.
The article is organised as follows.In section 2, the standard paradigm of neutrino flavour oscillations both in vacuum and in matter, as well as the deviations that the presence of NSIs would introduce, are summarised.In section 3, the ANTARES neutrino telescope is described.In section 4, the ANTARES data sample used for this analysis are presented and the Monte Carlo (MC) simulation of neutrino events, their reconstruction and the final event selection is discussed.The details of the analysis are given in section 5 and the results, in section 6.Finally, the conclusions of this work are gathered in section 7.

Neutrino propagation in matter and NSIs
As shown by a series of experiments with solar, atmospheric, accelerator and reactor neutrinos [3], neutrino flavour eigenstates are different from neutrino mass eigenstates.The flavour eigenstates are involved in neutrino production and annihilation in weak interactions, while the mass eigenstates determine the neutrino propagation.Since a neutrino produced with a given flavour evolves according to its superposition in terms of matter eigenstates, it can later interact as a different flavour, giving rise to neutrino flavour oscillations.
Neutrino oscillation probabilities are governed by the Pontecorvo-Maki-Nakawaga-Sakata mixing matrix (PMNS) [26,27] and the differences of squared masses.In the conventional three-flavour scheme with Dirac neutrinos, the relevant parameters are three mixing angles, θ 12 , θ 13 , θ 23 , two masssquared differences, ∆m 2  21 , ∆m 2  31 , and a CP violating phase, δ CP .Our present knowledge of these parameters can be summarised with the following ranges (see chapter 14 in ref. [3]): The octant of θ 23 , the sign of ∆m 2 23 (the mass ordering) and an accurate value of δ CP are yet to be determined.A number of neutrino experiments are planned to improve these measurements [28,29,30,31,32,33,34,35].
When neutrinos propagate in matter, their evolution is affected by the coherent forward elastic scattering on medium.The overall effect can be described by effective potentials associated to the CC and NC interactions.In the case of neutrinos travelling through the Earth, the only relevant potential is the one stemming from the electron neutrino components interacting with electrons in matter.The potential is given by where n e is the electron number density along the neutrino path [36].The relevant Hamiltonian is The PMNS mixing matrix, U , performs the rotation of the relevant mass matrix M2 = diag(0, ∆m 2 21 , ∆m 2  31 ) in the neutrino flavour space, where diag indicates a diagonal matrix with the specified elements.
The model-independent NSIs are introduced as new potentials, where n f is the fermion number density along the neutrino path.H NSI is added to the SM Hamiltonian in Eq. 3. In the present work, neutrinos are assumed to interact with down quarks which are roughly three times as abundant as electrons, n f = n d ≈ 3 n e .The matrix ε (ε α β , α, β = e, µ, τ ) gives the strength of NSIs.Its diagonal terms, if different from each other, give rise to the violation of leptonic universality, while the off-diagonal terms induce flavour changing neutral currents, which are highly suppressed in the SM.
The hermiticity of the hamiltonian matrix reduces the components of the ε matrix to nine real parameters, the three real diagonal plus the three complex off-diagonal elements.In addition, in oscillation experiments this matrix can only be determined up to a global multiple of the identity matrix, so only two parameters of the diagonal are independent, which are conventionally taken to be ε ee − ε µµ and ε τ τ − ε µµ .
To reduce the number of parameters to be fitted, some assumptions are customarily made depending on the nature of the experiment.Since ANTARES detects atmospheric neutrinos and in this study only track-like events are used, we essentially observe charged current interactions of atmospheric muon neutrinos, which are known to oscillate mainly to tau neutrinos 2 .Therefore, we are mostly sensitive to the NSI parameters ε µτ and ε τ τ − ε µµ .In this study, we assume that all other parameters are zero3 .Furthermore, since our sensitivity to phases is low, we assume ε µτ to be real, although this implies a loss of generality.
NSI are supposed to manifest themselves as sub-dominant effects on top of the main phenomenon of neutrino oscillations.In our case, the atmospheric neutrino oscillation channel ν µ → ν τ gives rise to the disappearance of muon neutrinos with a probability that depends on their energy and baseline.Therefore, a deficit with respect to the non-oscillation hypothesis should be observed in the detector as a function of the neutrino energy and arrival direction.In the range from ∼10 GeV to slightly above ∼100 GeV, the detector is capable of estimating the neutrino energy by measuring the muon range [37].The first oscillation minimum in the ν µ → ν µ channel occurs around 20 GeV, and the ν µ deficit is noticeable all the way up to 100 GeV, which enables ANTARES to measure the atmospheric neutrino oscillation parameters (sin 2 θ 23 , ∆m 2  23 ).As an example of the effect of NSI, Figure 1 shows the change in the oscillation pattern arising from non-zero values of ε µτ and ε τ τ − ε µµ .The difference of the survival probabilities for the case of NSI (P NSI νµ→νµ ) and of standard oscillations (P SI νµ→νµ ) is plotted as a function of the cosine of the zenith angle 4 , cos θ and the neutrino energy, E ν , in Figure 1.The left plot shows the difference for neutrinos and the right plot for anti-neutrinos.The test values of the oscillation parameters are adopted from the global neutrino oscillation fit results in ref. [38].Normal mass ordering is assumed.The NSI test point is chosen at (ε µτ , ε τ τ − ε µµ ) = (0.033, 0.147) 5 .Since the NSI effects are not symmetric for neutrinos and anti-neutrinos, they lead to a small but discernible effect in the combined {ν + ν} event distribution.This is important for the ANTARES detector which measures neutrinos via muons and is charge blind.It can be noticed that the effect of NSIs becomes prominent for higher energies and vertical upgoing events.

The ANTARES Detector
The ANTARES neutrino telescope [25] is located in the Mediterranean Sea, 40 km off the coast of Toulon, France, at a depth of about 2475 m.The detector, which was completed in 2008, is composed of 12 detection lines, each one equipped with 25 storeys of 3 optical modules (OMs), except line 12 with only 20 storeys of OMs, for a total of 885 OMs.The horizontal spacing among the lines is ∼60 m, while the vertical spacing between the storeys is 14.5 m.Each OM hosts a 10-inch diameter photomultiplier tube (PMT) from Hamamatsu [39], whose axis points 45 • downwards.All signals from the PMTs that pass a threshold of 0.3 photoelectrons are digitised (hits) and sent to the shore station.The position of the optical modules is determined with an accuracy of ∼10 cm [40] and the overall time calibration is better than 1 ns [41].The main sources of optical background registered by the ANTARES PMTs are the Cherenkov light induced by the decay products of the radioactive isotope 40 K and the bioluminescence.A detailed description of the ANTARES detector and further information about the data acquisition, trigger and calibration systems can be found in [25].

Data set, simulation, reconstruction and event selection
The data sample used in this work is the same as the one used in the previous ANTARES oscillation analysis [37].In that study, the data recorded by ANTARES from 2007 to 2016 (both years included), corresponding to a total livetime of 2830 days, was analysed to select a final sample of 7710 reconstructed track-like events.With this sample, a muon neutrino disappearance study was performed which provided a measurement of the oscillation parameters sin 2 θ 23 and ∆m 2  23 and established limits on the (3+1) sterile neutrino mixing model.The main features of the simulation, event reconstruction and selection are briefly described in what follows.Further details can be found in refs.[37,42,43].
CC ν µ interactions in seawater produce muons propagating through the detector that induce the emission of Cherenkov photons.They are identified as track-like events.All other types of neutrino interactions give rise to Cherenkov photons in the shower-like topology.The event reconstruction and selection used in the analysis have been optimised to select track-like events.On the other hand, ν e CC interactions and NC interactions of all flavours produce hadronic showers.In the case of ν e CC interactions, an electromagnetic shower is produced as well.Moreover, ν τ CC events can be produced as the result of ν µ → ν τ oscillations with or without muons in the final state.All these events constitute an additional source of background for this study.
Simulated atmospheric neutrinos are generated following the flux calculation of Honda et al. [44].Neutrino interactions surrounding ANTARES and inside the instrumented volume are simulated using the GENHEN software package [45].Atmospheric muon events are simulated by MUPAGE [46].Particle propagation is done using a GEANT-based package [47], which also simulates the Cherenkov light production and propagation taking into account the seawater optical properties [48].Optical background from 40 K decays and bioluminescence in water, as measured from counting rates in data, is added on a run-by-run basis [49], which provides a very accurate simulation of the environmental conditions.The overall geometry of the detector, the angular acceptance and quantum efficiency of the PMTs [39], as well as the response of the electronics is then taken into account to produce simulated digitised signals [50].
Two different muon track reconstruction algorithms are used in this analysis [51,52].In method A, a hit selection based on time and spatial coincidences of hits is applied and a χ 2 -fit is performed in order to extract the best track parameters.Events reconstructed by this method can have a single-line topology, if all the selected hits are recorded in the same detector line, or a multi-line topology, when hits belong to OMs of different lines.Method B involves a first prefit based on a directional scan of isotropically distributed directions, followed by a final log-likelihood fit of the track parameters using the best directions as starting points.
In the 10-100 GeV energy range, muons in water behave as minimum ionising particles, so that their energy is proportional to their pathlength, with a proportionality constant of 0.24 GeV/m (see chapter 34 in ref. [3]).The track length of the muon in the detector, L µ , can be calculated projecting back to the reconstructed track the first and the last selected hits.For singleline events, the track length is estimated from the vertical coordinates of the uppermost and lowermost storeys with selected hits and the reconstructed zenith angle.The reconstructed energy from the muon track length does not take into account the energy of the hadronic shower in the interaction vertex, nor the fact that the muon track might be only partially contained in the detector sensitive volume.Whereas for multi-line events, the threshold energy of the final neutrino sample is about 50 GeV, for nearly-vertical singleline events the energy estimate can be as low as 20 GeV or less, close to the first oscillation minimum.In the 20 GeV energy regime, the median angular resolution is 3 • for single-line and 0.8 • for multi-line events and the energy resolution for muons is found to be (50 ± 22)% [53].
The reconstructed muon energy of fully and partially contained events extends to around 100 GeV, but the true energy of the parent neutrinos goes up to a few TeVs.This broad range of neutrino energy brings enhanced sensitivity to the NSI parameters.
The events are selected using a quality criterion optimised using Monte Carlo simulated events.This quality criterion is based on the goodness-offit provided by each reconstruction algorithm: a reduced-χ 2 for method A and a log-likelihood for method B. Events reconstructed by any of the two methods that fulfill the corresponding quality criterion and that are upgoing (cos θ reco > 0.15) are kept for the analysis.
When these criteria are applied to a Monte Carlo sample of atmospheric neutrinos without oscillations plus atmospheric muons, corresponding to the 2830 days of data, a total of 8044 events are selected.More than 95% of these selected events are muon neutrinos, the background being mainly composed of atmospheric muons misreconstructed as upgoing.When applied to the data sample, the procedure yields a total of 7710 events, 5632 from method A (1950 from single-line and 3682 from multi-line) and 2078 from method B, in agreement with the expectations taking into account oscillations [37].Event samples from both methods, A and B, are selected for the final analysis.

Analysis
The selected events are distributed in a two-dimensional histogram of the reconstructed muon energy, E reco , versus the cosine of the zenith angle, cos θ reco , with a total of 136 bins.The muon energy spectrum is divided into eight bins: a first wide bin between 10 −0.3 = 5 GeV and 10 1.2 = 15.8GeV, plus seven bins equally logarithmically spaced between 10 1.2 and 10 2 GeV.The reconstructed cos θ reco is divided into seventeen bins between 0.15 and 1.0, the latter value corresponding to vertically up-going events.The difference in the number of events expected with NSIs (assuming a trial set of values of model parameters, quoted on the plots) and without NSIs is shown in Figure 2. In order to extract the parameters that best describe the observed energy versus zenith angle distribution, the following test statistic is defined: where n i,j and µ i,j (ō, s) are the number of measured and expected events in bin (i, j).{ō} and {s} are the model parameters, which represent the oscillation and the additional nuisance parameters, respectively (see Table 1).
The second sum runs over the systematic uncertainties s k , where ŝk is the assumed prior on the k-th parameter and σ 2 s k , its uncertainty [54].The expected number of events for each reconstructed (E reco , cos θ reco ) bin are calculated by applying oscillation and systematic parameter weights to un-oscillated events in the Monte Carlo event sample.Neutrino oscillation probabilities including NSI effects are calculated using the OscProb package [55] and the Earth density profile is approximated with 44 radial layers of constant density [56].
The impact of possible systematics effects on the analysis is taken into account introducing eight nuisance parameters (see Table 1 of section 6).The first four ones (∆m 2  31 , θ 23 , θ 13 and δ CP ) are those neutrino oscillation parameters whose uncertainties can influence the result of our analysis, while the last four ones (N ν , ∆γ, ν/ν and ∆M A ) reflect the uncertainties in the atmospheric neutrino flux predictions, the detector response and the neutrino interaction models.
The parameters of the standard oscillations are treated as follows.The solar mass splitting, ∆m 2  21 , and the corresponding mixing angle, θ 12 , are fixed to 7.4 × 10 −5 eV 2 and 33.62 • , respectively [3], since they have a negligible influence on the atmospheric neutrino oscillations.On the other hand, the atmospheric oscillation parameters ∆m 2  31 and θ 23 are treated as nuisance parameters left free during the fit, i.e. without any priors.The reactor angle, θ 13 , is assigned a Gaussian prior with a central value of 8.54 • and an uncertainty of ±0.28 • .Although no impact is found on the final results, the CP violating phase δ CP is fitted without prior.Separate fits are performed for the normal (NO) and inverted (IO) mass ordering hypotheses.All the above values are based on ref. [38].
Our treatment of the systematics stemming from the uncertainties in the neutrino flux predictions, the detector response and the neutrino interaction models follows the one we applied in our latest neutrino oscillation analysis [37].
For the atmospheric neutrino flux, the predictions by Honda et al.
[44] at the Frejus site (latitude 45.1 • N, close to that of ANTARES, 42.5 • N), are used.The systematics coming from the uncertainties in these predictions are considered introducing two nuisance parameters: a global neutrino normalisation factor, N ν , and a possible deviation for the standard value of the spectral index, ∆γ, that takes into account possible changes in the neutrino spectrum due to uncertainties in the primary cosmic ray spectrum.The global normalisation factor is treated as a free parameter without con-straints, while the possible change in the spectral index is introduced with a 5%-width Gaussian prior.Since only event rates are observed, these two parameters also take care of some other sources of systematics, as we explain below.
Uncertainties on the neutrino/anti-neutrino flux ratio, ν/ν, and on the flux asymmetry between upgoing and horizontal neutrinos, ν up /ν hor , have also been taken into account.These uncertainties [57] have been parametrised by the IceCube Collaboration [58] computing a correction on the number of expected events as a function of the neutrino energy, flavour, chirality, direction and the value of the uncertainty on the flux ratio.These two ratios considered are found to be strongly correlated [37], thus a unique nuisance parameter is considered in the fit.
The systematics related to the detector response come from the uncertainties in the OM photon detection efficiencies and the water absorption length.Dedicated MC simulations have been generated with modified OM photon detection efficiencies and a modified water absorption length, assuming a variation of ±10% from the nominal value, but keeping the same wavelength dependence.
The overall OM efficiency can be easily adjusted to the measured coincidence rates from 40 K decays [59] which makes the chosen 10% variations a conservative benchmark value, in line with early studies performed on ANTARES OMs [39].The water absorption length has been measured several times at the ANTARES site [48].The different measurements, taken at two different wavelengths, vary within about 10%.
The change in the event rates, expressed as the ratio of the event rates with the modified MC simulation to the one from the nominal MC simulation, has been computed as a function of the MC neutrino energy and zenith angle for ν µ CC events, reconstructed as upgoing.While no zenith-dependent effect is seen, the energy response of the detector is affected by these variations.The resulting distributions have been fitted, in the energy range 10−10 3 GeV, with a function of the form: where E T is the MC true neutrino energy, A , B are the two fitted parameters describing the effect of the modified OM photon detection efficiencies and E 0 = 100 GeV defines the reference energy for A .The ±10% variation in the OM efficiency translates into a ±20% in the event rate at 100 GeV that decreases linearly in log(E T ) to ±10% at 1 TeV (see Fig. 4 of ref. [37] for further details).The effect of the modified water absorption length is described by the same functional form of Eq. 6 using A w and B w as the corresponding fit parameters.The effects of A and A w are taken into account in the minimisation procedure by the global normalisation factor, N ν , while those of B and B w are covered by the uncertainty of the prior on the spectral index, ∆γ.Therefore, no new nuisance parameter has to be introduced to allow for systematics coming from OM efficiencies and water absorption.
The value of the atmospheric muon background contamination and its uncertainty have been obtained from the data itself following the method used in our oscillation paper [37] and described in the following.In the region of events with high χ 2 , which is dominated by atmospheric muons and excluded from the analysis, the data are fitted by exponential functions and extrapolated to the low χ 2 region, where the neutrino signal is dominant and is used for the analysis.The integral in that region gives an estimate of the atmospheric muons.Different fit ranges in the high χ 2 give different exponential fits.The final number of atmospheric muons is taken to be the mean and the uncertainty is estimated from the errors of the parameters of the fitted functions (for further details see section 5 of ref. [37] and section 6.5.3 of ref. [43]).These values are subsequently used as the Gaussian prior mean and deviation in the minimisation procedure.The energy and direction distribution of the atmospheric muon background has been estimated directly from MC.In the present analysis, we fix the atmospheric muon contamination at the value obtained in the oscillation analysis [37].Although this choice is motivated to ensure a stable fit procedure, it has been verified that this does not yield better constraints than an unconstrained muon normalisation.
Finally, a source of systematic uncertainty is the limited knowledge of the neutrino interaction model.In the energy range relevant for this analysis, the cross section is dominated by deep inelastic scattering (DIS) with smaller contributions from quasi elastic (QE) and resonant (RES) scattering.Uncertainties in the DIS cross section can be incorporated in the global flux normalisation factor N ν as well as in the correction to the spectral index ∆γ.For what concerns the QE and RES processes, dedicated studies have been performed [43] with gSeaGen [60], which uses GENIE [61] to model neutrino interactions.The dominant systematic error is found to be the one related to the axial mass for CC resonance neutrino production, M A .Its default value is 1.12 ± 0.22 GeV [61].By varying this parameter by ±1σ, a correction to the expected number of events as a function of the true neutrino energy has been computed.Although the correction is found to be small (see section 6.2.3 and Fig. 6.4 of ref. [43]) the parameterisation of this correction is used in the final fit.

Results
The best-fit point is obtained by minimising −2 log L in Eq. 5.Although the results are not expected to be strongly sensitive to the neutrino mass ordering (NMO), both possibilities are fitted separately.The values of the fitted parameters at the best-fit point for the NO and IO scenarios are shown in Table 1.The flipping of the signs with the change of the NMO in fits stems from the partial degeneracy between the sign of the ε µτ and the NMO [62].

Parameter
Prior Best-fit value (NO) Best-fit value (IO) N ν none 0.9 ± 0.1 0.9 ± 0.1 ∆γ 0.00 ± 0.05 −0.02 ± 0.04 −0.02 ± 0.04 Table 1: List of free parameters used in the maximisation of the likelihood ratio, their priors and the best-fit values for the case of NO and IO hypotheses.
Concerning the nuisance parameters, the following features can be observed.The values of ∆m 2  31 and θ 23 are in agreement within errors with our results for the standard oscillation analysis [37] and with the world best-fit values [3].The value of θ 13 is very close to its prior, indicating poor sensitivity to this parameter, which is to be expected since the ν µ survival probability depends on cos θ 13 , which is very close to one.
The global normalisation factor, N ν , is 10% lower than 1.This is in line with the result obtained in our standard oscillation analysis [37], where a value 18% lower than 1 was found.Both values are within the atmospheric neutrino flux uncertainties and are compatible with what has been reported by other analyses [58].The high pull in ν/ν is similar to the one obtained in our standard oscillation analysis.In that analysis, the conclusion was that this parameter seems to compensate the low value of N ν , as could be seen fixing all the nuisance parameters but N (see section 6.1 in Ref. [37]).
The parameter ∆M RES A remains very close to its prior, indicating a low sensitivity of the fit to it.As explained in the previous section, this is due to the relative high energy threshold of ANTARES.In Figure 3, the expected number of events using the best-fit values of Table 1 compared to the actual number of events observed by the experiment are shown.The number of events as a function of cos θ reco are displayed for the eight considered energy intervals, which are given in the figure in increasing order of energy from left to right and from top to bottom.The black crosses are the data with their statistical uncertainties and the red histograms are the expectations from the best-fit values obtained by maximising the log-likelihood in Eq. 5.The global minimum is found in the fit of NO.The dashed blue histogram is the expectation for the standard oscillation hypothesis.

ANTARES Preliminary
Exclusion contours in the (ε µτ , ε τ τ − ε µµ ) plane are drawn in the form of confidence level intervals for two degrees of freedom (d.o.f.).For each point of a grid in the (ε µτ , ε τ τ − ε µµ ) plane, the quantity −2 log L is minimised leaving free the rest of parameters.The differences between the values of −2 log L at the absolute minimum and at the grid points are used to build significance contours for two d.o.f., assuming that it obeys the χ 2 -distribution asymptotically [63].In Figure 4, the resulting exclusion limits at various C.L. intervals are shown.Finally, limits on ε µτ and ε τ τ − ε µµ are obtained by profiling over the other variable, as shown in the one-dimensional plots in the top and right panels of Figure 4.The 90% C.L. bounds are: We can summarise the results as follows.The hypothesis of standard interactions, i.e. the compatibility with zero of the combined fit result, is disfavoured with a significance of 1.7σ (1.6σ) for the normal (inverted) mass ordering scenario.However, the 90% C.L. contours do contain the null value.The difference from zero is mainly due to the ε τ τ − ε µµ coefficient.Indeed, while the value obtained for ε µτ profiling over the other coefficient is compatible with zero within its 68% C.L. interval, the 90% C.L. intervals for ε τ τ − ε µµ do not include the null value.However, the 95% C.L. interval for ε τ τ − ε µµ does contain zero.Therefore, we conclude that the results do not show clear evidence of deviations from standard interactions The limits obtained in this analysis are among the most stringent published to date for these NSI parameters, in particular for ε µτ .In Table 2, the 90% CL limits on ε µτ and ε τ τ − ε µµ provided by atmospheric neutrino experiments are summarised.
When comparing the results on Table 2 the following caveats should be considered.Some results were obtained assuming ε µτ to be real, as in our analysis, while some others assume it to be complex and fit both the modulus and the phase.For the latter, the number quoted is the 90% CL upper limit on the modulus.The two cases are indicated by notes (a) and (b) in the fourth column of the table.Although the sensitivity to the phase is not high in general for this sort of experiments and the limits tend to weaken with respect to the assumption of a real ε µτ , a complex ε µτ is the most general hypothesis.Some results are given using NSI couplings to down-quarks while some others are given for effective NSI couplings to electrons, protons and neutrons.The limits of the former should be multiplied by a factor ≈3 to compare to the latter.Limits obtained with the effective coupling and down-quarks couplings are denoted as (c) and (d), respectively.Whenever the authors themselves provide the translation from one case to the other we list both intervals (e.g. for IC DeepCore 2018 [19] translated in ref. [20]).All the limits are given for normal ordering.The table also shows the range of the reconstructed energy of the neutrino sample for each analysis.
In addition to the results based on atmospheric neutrino data, there is a large variety of bounds to NSI parameters from other analyses, some of which are collected in ref. [13].Bounds extracted from the data of one or more long baseline accelerator experiments (NOMAD, MINOS, T2K, NOνA) are given, for instance, in refs.[64,65,66].Bounds based on the data from the COHERENT experiment are given in refs.[67,68].Bounds extracted from global fits to the oscillation data plus the results of the COHERENT experiment are given in refs.[69,70,71].Notice that not all of these bounds can be compared right away to those obtained from atmospheric neutrinos and some of them (e.g.those from the CEνNS data) are model dependent.

Conclusions
The search for neutral current non-standard interactions in neutrino propagation using 10 years of ANTARES data is presented.We do not find clear evidence of deviations from standard interactions.For normal (inverted) neutrino mass ordering, the combined fit of both NSI coefficients yields a value of 1.7σ (1.6σ) away from the null result.However, the 68% and 95% confidence level intervals for ε µτ and ε τ τ − ε µµ , respectively, contain the null point.The resulting bounds of ANTARES obtained from this analysis are among the most restrictive to date.
ANTARES has recorded additional data from 2017 to 2022, which represents a ∼50% increase with respect to that used in the present analysis and are currently being analysed.Furthermore, the next generation of neutrino telescopes in the Mediterranean Sea, KM3NeT-ORCA and -ARCA, now under deployment [72,73,74,75], are expected to further improve these limits.

Figure 2 :
Figure2: Difference of number of events with and without NSI as a function of the reconstructed neutrino energy (log 10 (E reco )) and direction (cos θ reco ).Both panels correspond to a NSI test point, parameterised by ε αβ as quoted on the plots for normal (left) and inverted (right) orderings.

Figure 3 :
Figure 3: Data (black points) and MC simulated events from this analysis, as a function of the reconstructed zenith angle, cos θ reco , for the eight different energy, E reco bins.The blue dashed histogram corresponds to the MC assuming standard oscillations and the (mostly overlapping) red histogram corresponds to the MC for the NSI case using the best-fit values obtained for normal ordering.Note the difference in Y-axis scales in different panels.

Figure 4 :
Figure 4: Contours of 68.3% (solid), 90% (dashed) and 95% (dotted) confidence level on the (ε µτ , ε τ τ − ε µµ ) plane, after 10 years of ANTARES livetime for the NO (left) and the IO (right) cases.The cross indicates the best-fit point obtained in both cases.The lateral plots on both panels show the 1D profile likelihood of the respective NSI parameters under study, when the other parameter is fitted over.The dashed straight lines in the lateral plots indicate the 90% C.L. for one d.o.f.

ε− 2 )Table 2 :
µτ or |ε µτ | (×10 −3 ) ε τ τ − ε µµ (×10 Summary of the 90% CL upper limits on the NSI parameters ε µτ (or |ε µτ |) and ε τ τ − ε µµ from different atmospheric neutrino experiments.The intervals indicate the regions not excluded by the corresponding analysis.For |ε µτ | the single upper limit is given.The label "Pub" indicates that the public data have been analysed by researchers outside the collaboration.The conditions applied are as follows: (a) assuming real ε µτ , (b) modulus of complex ε µτ , (c) for effective NSI couplings and (d) for d-quark NSI couplings.All the results are given for normal ordering.The approximate neutrino energy range used in each analysis is also shown.See the main text for explanations.