Measuring the atmospheric neutrino oscillation parameters and constraining the 3+1 neutrino model with ten years of ANTARES data

The ANTARES neutrino telescope has an energy threshold of a few tens of GeV. This allows to study the phenomenon of atmospheric muon neutrino disappearance due to neutrino oscillations. In a similar way, constraints on the 3+1 neutrino model, which foresees the existence of one sterile neutrino, can be inferred. Using data collected by the ANTARES neutrino telescope from 2007 to 2016, a new measurement of $\Delta m^2_{32}$ and $\theta_{23}$ has been performed - which is consistent with world best-fit values - and constraints on the 3+1 neutrino model have been derived.

Neutrino oscillations have been detected by a variety of experiments, studying solar as well as atmospheric neutrinos, but also neutrinos produced from nuclear reactors and particle accelerators. For a comprehensive review see [4].
Atmospheric neutrinos are produced through the interaction of cosmic rays with nuclei in the Earth's atmosphere. Their flux spans many orders of magnitude in energy, from GeV to hundreds of TeV. Being isotropic to first order, it allows to investigate a large range of baselines on the Earth's surface, from ∼10 km of vertically down-going to ∼10 4 km of vertically up-going neutrinos.
In this paper the muon disappearance channel (P νµ→νµ ) is studied. The vacuum survival probability for a muon neutrino of energy E interacting at a distance L from its creation point is given by: where U µi , U µj are elements of the PMNS matrix U , and ∆m 2 ji = m 2 j − m 2 i are the mass splittings between two mass eigenstates. The rightmost term describes the "single ∆m 2 dominance" approximation, relevant in the energy domain considered for this analysis. Here the ν µ survival probability depends only on U µ3 = sin θ 23 cos θ 13 and ∆m 2 32 . For a vertically up-going atmospheric ν µ , the first minimum of the survival probability described in Equation 1 is reached at energies of ∼ 25 GeV. The formalism given in Eq. 1 is further modified by matter effects [5,6,7] as the neutrinos propagate through the Earth. Throughout the paper, oscillation probabilities are calculated with the OscProb package [8] which treats matter effects for an arbitrary number of neutrino families numerically without approximations.
The ANTARES neutrino telescope [9] has been designed and optimised for the exploration of the high-energy Universe by using neutrinos as cosmic probes. However, its energy threshold of about 20 GeV is sufficient, even if at the edge, to be sensitive to the first atmospheric oscillation minimum, making also the study of neutrino oscillations possible. As neutrinos and antineutrinos are indistinguishable on an event-by-event basis in neutrino telescopes, in the following muon (electron) neutrinos are refered to the sum of contributions from both neutrinos and antineutrinos.
A previous analysis of ANTARES data, covering the data acquisition period from 2007 to 2010, represented the first study of this kind performed by a neutrino telescope, and measured the atmospheric neutrino oscillation parameters, ∆m 2 32 and θ 23 [10]. In the present work, data collected during 10 years have been studied with a new analysis chain that also includes a more comprehensive treatment of various systematic effects.
Despite the fact that neutrino oscillation is a well established phenomenon, some observed experimental anomalies, such as the ones reported by the LSND [11] and MiniBooNE [12] collaborations, seem to indicate a deviation from the standard 3-flavour picture. These discrepancies could be partially explained by introducing in the model an additional neutrino state. However, since the number of weakly interacting families of light neutrinos is limited to three by the LEP results [13], the additional neutrinos have to be sterile, i.e., they do not undergo weak interactions. The 3+1 neutrino model foresees the existence of one sterile neutrino in addition to the three standard ones. A choice has to be made, how to extend the mixing matrix U from three to four families. In this analysis the convention from [14] (see "supplementary materials") is adopted: where R ij is the rotation matrix for angle θ ij . If j − i > 1, R ij also contains a CP-violating phase, δ ij . Six new real mixing parameters have to be accounted for: three new mixing angles, θ 14 , θ 24 and θ 34 , a new mass splitting, ∆m 2 41 , and two new phases, δ 14 and δ 24 . In line with other analyses of sterile neutrinos in the muon disappearance channel [14,15,16], θ 14 = 0 is assumed, which also eliminates any dependency on δ 14 . Even though a sterile neutrino does not interact as the active flavours, its presence would still modify the oscillation pattern of the standard neutrinos, due to the fact that the standard neutrino flavours could oscillate into these additional sterile species. In particular, for up-going ν µ in the energy range of 20-100 GeV, non-zero values of U µ4 and U τ 4 with can lead to distortions in their survival probability. This is illustrated in Figure 1 which shows the ν µ survival probability for maximal mixing of θ 23 and different combinations of the mixing parameters θ 24 , θ 34 and δ 24 . If only θ 34 is non-zero, the survival probability of ν µ with respect to the non-sterile hypothesis is only modified close to the first oscillation minimum. The case of both θ 24 and θ 34 being non-zero leads instead to a significant shift of the first oscillation minimum in energy (depending on δ 24 ) and modifies the event rate up to energies of few hundred GeV, easily accessible with ANTARES. The fast wiggles due to ∆m 2 41 = 0.5 eV 2 will be smeared out by detector resolution effects, therefore no sensitivity to this parameter is expected. The surprisingly strong effect of δ 24 on the ν µ survival probability, neglected in all similar analyses so far, is further detailed in the Appendix.
Since the effect of an additional sterile neutrino would be visible in the same energy and zenith range as the ν µ disappearance, the same analysis chain and data sample can be exploited to constrain the 3+1 neutrino model parameters. In this paper, the results of an investigation aiming to constrain the mixing angles θ 24 and θ 34 of the 3+1 neutrino model are also reported.
The paper is organised as follows: in Section 2 the ANTARES neutrino telescope is briefly described and its detection principle is illustrated; the ANTARES data sample as well as the Monte Carlo (MC) chain are presented in Section 3, while the event reconstruction is discussed in Section 4. Section 5 is dedicated to the event selection and the minimisation procedure. The results are presented in Section 6, while conclusions are given in Section 7.

The ANTARES neutrino telescope
The ANTARES neutrino telescope is located in the Mediterranean Sea, 40 km off the coast of Toulon, France, at a mooring depth of about 2475 m. The detector was completed in 2008. ANTARES 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 photomultiplier tube (PMT) from Hamamatsu [17], whose axis points 45 • downwards. All signals from the PMTs that pass a threshold of 0.3 single photoelectrons (hits) are digitised and sent to the shore station [18,19]. The on-shore trigger system [20] performs a hit selection based on causality relations and builds events under the hypothesis that the selected hits originate from Cherenkov radiation induced by relativistic charged particles as they are produced in neutrino interactions close to the ANTARES instrumented volume.
The main sources of optical background registered by the ANTARES PMTs are represented by Cherenkov light from decay products of the radioactive isotope 40 K, naturally present in sea-water, by light emitted through bioluminescence by living organisms, and by energetic atmospheric muons, which can penetrate deeply under the sea and reach the detector from above. The aim of the MC production is to reproduce in the most realistic way the events expected at the detector, as well as the response of the apparatus when recording these events. In order to account for changes of the environmental conditions, as well as for the different operational status of the detector and its components over time, a run-by-run MC approach is applied [21]. A typical run lasts few hours. Several time dependent conditions are taken from real data and applied to the run-by-run MC. First, temporarily or permanently non-operational OMs are masked in the simulation. Secondly, background light conditions, which might vary due to bioluminescence, are measured every 104 ms for each individual OM. These samples are directly used as input for the background light simulation. Thirdly, individual OM efficiencies are considered, as calculated on an approximately weekly basis from 40 K coincidence rates [22]. Finally, the acoustics based position calibration, performed every few minutes, is applied. All these detailed inputs assure an authentic description of the detector response for each individual run. Remaining uncertainties are small and can be handled as global parameters which are discussed below. They are included in the analysis as systematic uncertainties.

ANTARES data and Monte Carlo samples
Neutrino interactions of all flavours have been simulated with the GENHEN [23] package, developed inside the ANTARES Collaboration. It allows to reproduce neutrino interactions in the GeV to multi-PeV energy range. MC neutrino events can be weighted to reproduce different physical expectations. For atmospheric neutrinos with E ν ∈ [20 − 100] GeV, a MC sample almost three hundreds times larger than the data sample is available. The model by Honda et al. [24] for the Fréjus site is used in this work.
Even though the sub-marine location of ANTARES provides a good shielding against atmospheric muons, still a large amount of them will reach the detector. The event generator used in ANTARES to simulate atmospheric muons is MUPAGE [25]; the energy and angular distributions, as well as the multiplicity of muons propagating in sea water are parameterised. The contribution from this background is also evaluated from the data itself.
Particle propagation and Cherenkov light production are simulated using a GEANT-based [26] package [23], which takes into account all relevant physics processes and computes the probability that photons emitted by a particle reach the OM surface, producing a hit. Finally the detector response is simulated, including the digitisation and filtering of hits. At this stage a realistic optical noise is added on each OM for each data acquisition run of the detector, and the time evolution of the detector configuration is accounted for as described above.

Event reconstruction
Charged-current (CC) interactions of muon neutrinos produce a muon propagating through the detector and inducing Cherenkov light. They are identified as track-like events. The event reconstruction and selection used in the analysis have been optimised to select such events. On the other hand, ν e CC interactions, as well as neutral-current interactions (NC) 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 and without muons in the final state. All these events constitute an additional source of background for this study.
Events have been reconstructed using two different algorithms, described in detail in [27,28]. In the following discussion these algorithms will be referred to as method A and method B, respectively. Both are optimised for events induced by GeV-scale ν µ CC interactions. 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 find the best track. Events can have a single-line topology (SL), if all the selected hits have been recorded in the same detector line, or a multi-line topology (ML), when hits belong to OMs of different lines. Method B consists of a chain of fits, aimed to improve at each step the track estimation. Starting from a hit selection, a first prefit, based on a directional scan with a large number of isotropically distributed directions, is performed. The best 9 directions are used as starting points for the final likelihood (log L) fit.
Once the muon track has been reconstructed, its length, L µ , is computed. This is done, for ML events, by projecting back to the track the first and last selected hit. For SL events, since a vertex estimation is not possible due to the lack of azimuth information, the track length is estimated from the z-coordinates of the uppermost and lowermost storey which have recorded the selected hits and taking into account the reconstructed zenith angle.
The muon energy estimation is based on the fact that muons in the few-GeV energy range can be treated as minimum ionising particles, and their energy can be estimated from their track length L µ : where the factor 0.24 GeV/m represents the energy loss of muons in sea water in the energy range of 10-100 GeV [29]. This quantity is used in the following as estimator for the neutrino energy. The energy resolution of fully contained muons is dominated by the spacing of the detector elements and is found to be around 5 GeV. For muons leaving the detector only a lower limit for their energy can be derived, corresponding to their visible length inside the instrumented volume. More details on the muon energy resolution can be found in [10].

Analysis
To achieve the best sensitivity to the measurement of the oscillation parameters, a set of quality criteria has been applied. The selection of ν µ CC events has been optimised by performing a preliminary Monte Carlo (MC) sensitivity study, before applying the whole analysis chain to data. The main parameter on which the selection is based is the reduced χ 2 for method A and the log L for method B. Events reconstructed by method A and passing the corresponding event selection are kept. The events discarded by this procedure are further reconstructed by method B; they are kept in the analysed sample if the corresponding selection criteria are passed. Only events which are reconstructed as up-going are used in the following. A minimum number of five storeys with selected hits is required, in order to minimise the background induced by atmospheric muons.
In Figure 2 the distribution of the MC true neutrino energy, E T , for selected ν µ CC events is shown. For the histogram with the solid line no neutrino oscillations are assumed, while the dashed one refers to a 2-flavour oscillation scenario with maximal mixing and ∆m 2 32 = 2.46 × 10 −3 eV 2 . As can be seen, atmospheric neutrino oscillations affect the expected event distribution for E T 100 GeV. About 7590 well-reconstructed ν µ CC events are expected in a live time of 2830 days when oscillations are neglected. Roughly one half of these events are reconstructed with method A (ML), while methods A (SL) and B both contribute with approximately 25% to this event sample. Further, ∼40 ν e CC events are selected. Oscillations reduce the number of expected events by ∼720 events. This reduction is dominantly seen in the A (SL) sample (∼ 60%) which contains the lowest energetic and most vertical events, while the other two reconstruction methods contribute each about 20%. ν τ CC events reduce this oscillation signal by ∼ 20 events, taking into account the energy-dependent cross section ratio σ(ν τ CC)/σ(ν µ CC) (about 0.5 at 25 GeV), the 17% branching ratio of the muonic τ decay and the resulting soft spectrum of the produced muons. Figure 3 shows the distribution of the reduced χ 2 SL for method A (SL) events where data are compared to simulated atmospheric neutrinos and background atmospheric muons. While the MC reproduces quite well the data in the signal region dominated by the neutrino signal, a disagreement between the MC expectation and data is visible for larger χ 2 SL . Both data and MC follow an exponential law in this region, but with different slopes. For this reason, the number of background atmospheric muons in the signal region has been determined from data itself. The distribution in Figure 3 has been parameterised in the region dominated by atmospheric muons (χ 2 SL > 0.8) with four different exponential fits by varying   the fit range. Each fit has been extrapolated into the signal region, and its corresponding integral has been computed. The mean of these integrals has been used to estimate the number of atmospheric muon background, and its uncertainty has been computed from the errors on the fitted function parameters. Summing up the results of this method for events that have been reconstructed by method A (SL and ML) and method B, and combining the corresponding errors in quadrature, a total background of 740 ± 120 atmospheric muons has been determined. This value is subsequently used as a Gaussian prior mean value and uncertainty in the minimisation procedure. The energy and direction distribution of the atmospheric muon background has been, instead, estimated directly from MC.
After applying the event selection criteria described above on the data sample, a total of 7710 events have been selected, 1950 from method A (SL), 3682 from method A (ML) and 2078 from method B. In Figure 4 the event distribution as a function of the logarithm of the reconstructed energy, log 10 (E reco /GeV), and the cosine of the reconstructed zenith angle, cos θ reco , is shown. The distribution of the MC expectation assuming no neutrino oscillation (left panel) is compared to what is observed in data (right panel). Eight bins in log 10 (E reco /GeV) have been considered, seven from 1.2 to 2.0, plus an additional underflow bin which accounts for all events with log 10 (E reco /GeV) < 1.2; there are 17 bins in cos θ reco , from 0.15 to 1.0, the latter denoting vertically up-going events. The final fit has been performed on the 2-dimensional histograms shown in Figure 4. The fit follows a log-likelihood approach, by minimising the function: where the first sum runs over the histogram bins of log 10 (E reco /GeV) and cos θ reco , N data i,j is the number of events in bin (i,j) and N M C i,j (p,η) the corresponding number of expected MC events in the same bin. This number depends on the set of oscillation parameters,p, that are under investigation, as well as on the set of parameters related to systematic uncertainties,η, as described in the next subsections. The dependency on oscillation parameters is taken into acount for CC interactions of all neutrino flavours which contribute to the final event sample. The second sum runs over the number of nuisance parameters taken into account, < η k > being the assumed prior of the parameter k, and σ η k its uncertainty. The log-likelihood function converges to the standard χ 2 for bins with high statistics. For bins with a small number of entries the log-likelihood is more adequate.
Since the treatment of the systematic uncertainties slightly differs between the standard atmospheric oscillation analysis and the sterile neutrino analysis, they are described separately in the following subsections.

Treatment of systematics for the standard oscillation analysis
The standard oscillation analysis accounts for six sources of systematic uncertainties. Three are related to the atmospheric neutrino flux. A global neutrino normalisation factor, n ν , which is left unconstrained during the fit, accounts for uncertainties on the total number of expected events. A variation ∆γ in the nominal neutrino flux spectral index has been used as additional nuisance parameter. Uncertainties on the neutrino/anti-neutrino flux ratio, ν/ν, and on the flux asymmetry between up-going and horizontal neutrinos, ν up /ν hor , have also been taken into account. These uncertainties [30] have been parametrised by the IceCube Collaboration [31]. Such parameterisations compute 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. The two ratios considered in this analysis have been found to be strongly correlated, thus a unique nuisance parameter is considered in the fit.
An additional source of systematic uncertainty is the limited knowledge of the neutrino interaction model. At the energy of interest for this study, the cross section is dominated by deep inelastic scattering (DIS), with a smaller contribution 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 with gSeaGen [32], which uses GENIE [33] to model neutrino interactions. The dominant systematic is found to be related to the axial mass for CC resonance neutrino production, M A . Its default value is 1.12 ± 0.22 GeV [33]. By varying this parameter by ±1σ, the correction with respect to the expected number of events has been computed as a function of the true neutrino energy and this parameterisation is used in the final fit.
Apart from the oscillation parameters under investigation, ∆m 2 32 and θ 23 , the other oscillation parameters may play a role, but their effect is limited for this study. In particular, θ 13 is left free in the fit but treated with a Gaussian prior at θ 13 = (8.41 ± 0.28) • , which is taken from a global fit [34] as well as the values of the solar neutrino parameters, which are kept fixed: ∆m 2 21 = 7.37 × 10 −5 eV 2 and sin 2 θ 12 = 0.297. Different values of δ CP have been tested at the stage of the MC sensitivity study and found to have no impact on the final result. Therefore δ CP is fixed at zero.
The number of atmospheric muons, N µ , contaminating the neutrino sample, is treated as an additional nuisance parameter. Its value and uncertainty, determined with the data-driven technique, are used as a prior.
Finally, detector and sea water related systematics have been studied as well. 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 [22] which makes the chosen 10% variations a conservative benchmark value, in line with early studies performed on ANTARES OMs [17]. The water absorption length had been measured several times at the ANTARES site [35]. The different measurements, taken at two different wavelengths, vary within about 10%.
The correction to the event rates, obtained by dividing the event rates from the modified MC simulation (r var ) and the one from the nominal MC simulation (r nom ), has been computed as a function of the MC neutrino energy and zenith angle for ν µ CC events, reconstructed as up-going. 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 . Figure 5 shows the distribution of the event ratios as a function of true neutrino energy, together with its parameterisation. 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 values of the fitted parameters A , B , A w and B w are listed in Table 1. The effects of A and A w are taken into account in the minimisation procedure by the global normalisation factor, n ν , which is left unconstrained, while B and B w are covered by the uncertainty of the prior on the spectral index, ∆γ (see Table 2).  Table 1: Fitted values for the parameterisation of the event weight correction with a variation of ±10% from the nominal value of the OM photon detection efficiency and water absorption length.

Treatment of systematics for the sterile oscillation analysis
For the sterile analysis, the flux as well as the cross section related systematic uncertainties are treated in the same way as described in the previous subsection.
Since the effect of a sterile neutrino would modify the oscillation pattern in a similar way as ∆m 2 32 and θ 23 do, these parameters are considered to be one of the sources of systematic uncertainty for this analysis. Both ∆m 2 32 and θ 23 are left unconstrained as recommended in [36]. The other standard oscillation parameters are treated as previously discussed.
As discussed in Section 1, the addition of a sterile neutrino in the model implies six new mixing parameters to be accounted for. The mixing angle θ 14 and its associated phase δ 14 have been fixed at zero, since they mainly affect the ν e channel. The fast oscillations due to ∆m 2 41 0.5 eV 2 are unobservable due to the limited energy resolution of the detector, making ∆m 2 41 not measurable. It has been kept fixed at 0.5 eV 2 . The choice of the neutrino mass hierarchy (NMH) as well as δ 24 are expected to impact the result. Therefore both normal and inverted hierarchy (NH/IH) and various values of δ 24 have been tested during the fit. Furthermore, to ensure the stability of the fit procedure, the atmospheric muon contamination has been fixed at the value found by the standard oscillation analysis. It has been verified that this choice does not lead to better constraints with respect to the case of a free muon contamination.

Results
The minimisation procedure has been done using the ROOT package Minuit2 [37], applied to the function introduced in Equation 5. Results are presented in the following subsections, for the standard oscillation analysis and the sterile oscillation analysis, respectively.

Results for the standard oscillation analysis
In Table 2 the complete list of all the fitted parameters for the standard oscillation analysis is shown, together with their best-fit values and their priors. Due to the high energy threshold of ANTARES this analysis is not sensitive to the NMH. The results hold for both NH and IH. The best-fit value is found for ∆m 2 32 at (2.0 +0.4 −0.3 ) × 10 −3 eV 2 , which is compatible with the current world best-fit value [38]. The mixing angle θ 23 is found to be compatible with maximal mixing within its error. The global normalisation factor for neutrinos, n ν , is found to be 18% lower. This value is within the atmospheric neutrino flux uncertainties and it is compatible with what was reported by other analyses [31]. A nonnegligible pull is found on ν/ν. This parameter seems to compensate for the low value of n ν : this has been derived from an alternative fit, for which all nuisance parameters but n ν have been fixed, to allow a more direct comparison with the result reported in [10]. Under these conditions n ν = 1.04 ± 0.02 is found. Concerning the spectral index correction, ∆γ, no significant distortion from the nominal value is observed. The fitted value for the atmospheric muon contamination shows a strong pull and it is found incidentally close to the MC expectations. For both θ 13 and M A the best fit values and their errors are found at the corresponding prior, which indicates no sensitivity to these parameters. This can be understood as the ν µ survival probability does not depend on sin θ 13 but only on cos θ 13 = 0.99 (see Eq. 1) whereas M A mainly affects neutrinos with energies below the detection threshold of ANTARES. 0.0 ± 1.0 0.0 ± 1.0 Table 2: Priors and fitted values obtained from the minimisation for all the parameters considered in the standard oscillation analysis.
The distribution of the ratio between the reconstructed energy and the cosine of the reconstructed zenith is shown in Figure 6. This ratio is affected by the oscillation phenomenon as can be seen for the lowest values of E reco / cos θ reco . For comparison, also the distribution of MC assuming no neutrino oscillation, as well as the one assuming the world best-fit values [38] are shown. The latter two are calculated with all nuisance parameters at their nominal values. Such a 1D distribution does not carry the full information exploited in the fit, which is performed on the 2D distribution shown in Figure 4. While compatible with world data, ANTARES results seem to prefer a somewhat shallower (or energy shifted) oscillation minimum.  Figure 6: E reco / cos θ reco distribution for data (black), MC without oscillation (red), MC assuming the world best-fit values (blue) [38] and MC assuming best-fit values of this analysis (green). The left plot shows event numbers while the right plot illustrates the event ratio with respect to the MC without oscillations.
In Figure 7 the 90% CL contour obtained in this work, in the plane of sin 2 θ 23 and ∆m 2 32 , is compared to those published by other experiments. The 1D projections, after profiling over the other variable, are shown as well. Confidence level contours have been computed by looping over a fine grid of values in ∆m 2 32 and θ 23 and minimising the negative log-likelihood over all the other parameters. The non-oscillation hypothesis has been tested by performing the minimisation with a fixed null value of the oscillation parameters, and it is discarded with a significance of 4.6σ, compared to 2.3σ in our previous analysis [10].

Results for the sterile oscillation analysis
In Table 3 the complete list of all the fitted parameters for the sterile oscillation analysis for NH and IH is shown, together with their best-fit values and their priors. While θ 24 is found to be compatible with zero, the best fit for θ 34 is found at a non-zero value. This can be understood from the slight preference of the ANTARES data for a shallower oscillation dip (see discussion related to Fig. 6), which can be easily provided by a non-zero value of sin θ 34 (see Fig. 1). The non-sterile hypothesis is found at −2∆ log L = 4.4 which corresponds to a 2-parameter p-value of 11%. The fitted values of ∆m 2 32 and θ 23 are slightly different but consistent with respect to the ones obtained in the standard oscillation analysis. The complex phase δ 24 is found at 180 • . For IH instead the fit prefers δ 24 = 0 • with otherwise identical results, as expected from the degeneracy between NMH and δ 24 (see Appendix). For the other parameters a similar behaviour as for the standard oscillation analysis is observed.
Exclusion contours are built by applying Wilks' theorem. In Figure 8 the resulting 90% and 99% CL exclusion limits have been computed on a 2D grid in the plane of the two matrix elements, namely |U µ4 | 2 = sin 2 θ 24 and |U τ 4 | 2 = sin 2 θ 34 cos 2 θ 24 . The exclusion limit for unconstrained δ 24 , which corresponds to both [NH,δ 24 = 180 • ] or [IH,δ 24 = 0 • ], can be directly compared to the IceCube/DeepCore [16] (IH) limit. Also shown are limits for NH and δ CP = 0 • which allow a direct comparison with the results from IceCube/DeepCore [16] (NH) and Super-Kamiokande [15]. All three experiments find the best fit for |U τ 4 | 2 to differ from zero. Our results exclude regions of the parameter space not yet excluded by other experiments.
The IceCube/DeepCore analysis [16] is limited to events with reconstructed energy lower than 56 GeV, while the distortion on the oscillation pattern possibly produced by the presence of a sterile neutrino Figure 7: Contour at 90% CL in the plane of sin 2 θ 23 and ∆m 2 32 obtained in this work (black line) and compared to the results by other experiments: IceCube/DeepCore (red) [31], Super-Kamiokande (green) [39], NOνA (purple) [40], T2K (blue) [41], and MINOS (light blue) [42]. The lateral plots show the 1D projections on the plane of the two oscillation parameters under study.  would be evident also at higher reconstructed energies. The present analysis includes events with reconstructed energy up to 100 GeV. It has been verified that the ANTARES limits degrade when restricting the analysis to events with E reco < 56 GeV. In this work both of the standard atmospheric oscillation parameters ∆m 2 32 and sin 2 (2θ 23 ) are left unconstrained in line with the IceCube/DeepCore analysis [16].  After profiling over the other variable, the following limits on the two matrix elements can be derived: |U µ4 | 2 < 0.007 (0.13) at 90% (99%) CL, |U τ 4 | 2 < 0.40 (0.68) at 90% (99%) CL.

Conclusions
Ten years of ANTARES data have been analysed to provide a measurement of the atmospheric neutrino oscillation parameters. The analysis chain has been optimised with respect to our previously published study, by combining two track reconstruction algorithms and introducing a more elaborate treatment of various sources of systematic uncertainties. The results, ∆m 2 32 = (2.0 +0.4 −0.3 )×10 −3 eV 2 and θ 23 = (45 +12 −11 ) • , are consistent with what has been published by other experiments. The non-oscillation hypothesis is discarded with a significance of 4.6σ.
Exploiting the same analysis chain and the same data set, a further study has allowed to constrain, for the first time with ANTARES, the parameter space of the 3+1 neutrino model, which foresees the existence of one sterile neutrino.

Appendix : Sterile neutrinos and matter effects
For the analysis presented in this paper, oscillation probabilities are evaluated with the software package OscProb [8]. However, in this Appendix some common approximations are applied, to derive analytical formulae. These are NOT used for the analysis itself but allow to get a better understanding of the interplay between different parameters. The ν µ survival probability in vacuum in the 3 + 1 model can be simplified with the following two hypotheses [43,15]: first, it is assumed, that the first generation decouples completely, i.e. ∆m 2 21 = 0 and θ 12 = θ 13 = 0; second, fast wiggles due to oscillations involving m 4 are assumed to be unobservable, i.e. sin 2 (∆m 2 4i L/4E) = 1/2 for all i. This yields with P µµ the ν µ survival probability in the 3-flavour scheme, i.e. without additional sterile neutrinos. Only |U µ4 | 2 = sin 2 θ 24 can be probed in this scheme, which is applied in most accelerator based ν µ disappearance analyses. However, when analysing atmospheric neutrinos, matter effects cannot be neglected. An analytical formalism is developped in Eqs. 4.13-4.25 of [15]. In Eq. 4.13, a complex phase is present in the non-diagonal term of the matrix, which is neglected, i.e. set to zero, in subsequent steps. If instead this phase is kept, sin 2θ s in Eq. 4.16 acquires an extra term exp(−iδ).
For δ = 0 the original expressions from [15] are reproduced. With A 32 = ∆m 2 32 /E ν and A s = √ 2 2 G F N n (|U µ4 | 2 + |U τ 4 | 2 )/2 (G F the Fermi constant and N n the neutron density) the ν µ survival probability in matter is fully defined and can be written equivalently to Eq. 9 (see also Eq. 4.23 of [15]): which describes well all features shown in Fig. 1. The impact of the CP-phase δ disappears when either |U µ4 | 2 = 0 or |U τ 4 | 2 = 0, which leads to sin 2θ s = 0. Further, δ → δ + π is completely degenerate with changing the mass hierarchy, i.e. swapping the sign of A 32 if either cos 2θ 23 = 0 or cos 2θ s = 0. Deviation from maximal mixing in θ 23 or from the symmetry between |U µ4 | 2 and |U τ 4 | 2 defining θ s breaks this degeneracy. The impact of the neutrino mass hierarchy on the ν µ survival probability in matter had been pointed out already in [44], while the influence of complex phases is also discussed in [45].