Measurement of the atmospheric ν e and ν μ energy spectra with the ANTARES neutrino telescope

This letter presents a combined measurement of the energy spectra of atmospheric ν e and ν μ in the energy range between ∼ 100 GeV and ∼ 50 TeV with the ANTARES neutrino telescope. The analysis uses 3012 days of detector livetime in the period 2007–2017


Introduction
Atmospheric neutrinos are secondary particles produced by cosmic ray (CR) primaries interacting in the Earth's atmosphere.Due to the need of very large detectors, only a few measurements of the differential flux exist, namely from the AMANDA [1,2], IceCube [3,4,5,6,7] and Super-Kamiokande [8] Collaborations, and a historical measurement from the Frejus Collaboration [9].The ANTARES Collaboration has reported a measurement of the atmospheric ν µ energy spectrum in [10].
Different theoretical frameworks are available to estimate atmospheric neutrino fluxes [11,12,13,14].At energies from 100 GeV to 1 PeV, the main source of ν µ are semi-leptonic and three-body decays of charged kaons, while the contributions from pion and muon decays dominate below 100 GeV.This conventional neutrino flux tends towards a power law Φ , where γ CR is the spectral index of the primary CRs.
Above 100 GeV and up to some tens of TeV, atmospheric ν e 's come mostly from decays of neutral and charged kaons, and have the same spectral index of conventional ν µ .Below 100 GeV, ν e are predominantly produced by muon decays.The ν µ /ν e flux ratio is ∼2 in the GeV range and increases with energy, reaching a factor of ∼20 at 1 TeV.
At high energies, equal fluxes of ν µ and ν e are produced by the decays of charged and neutral D-mesons.Because of the very short lifetime of these mesons, the resulting flux is called prompt neutrino flux [15,16] and its energy spectrum, Φ p ν ∝ E −γ CR ν , follows the primary spectrum up to very high energies.The transition from the region in the spectrum dominated by conventional neutrinos to prompt neutrinos is expected to occur at E ν ∼ 1 PeV for ν µ and around E ν ∼30 TeV for ν e .As a rule of thumb, the primary CR energy is about 20 times higher than the energy of the secondary neutrino.Uncertainties on the conventional flux models at neutrino energies above 1 TeV are mainly due to a poor knowledge of primary CR energy spectrum and composition, and of hadronic interactions, in particular of strange quark production mechanisms [11].For a recent, detailed description of the hadronic interactions leading to inclusive lepton fluxes, refer to [17].
Finally, ν τ production in the atmosphere is rare: this is dominated by the decay D + s → τ + ν τ , followed by τ decay.As oscillation effects for atmospheric ν's are negligible above ∼100 GeV, the contribution from tau neutrinos is not considered in this analysis.
This letter describes a strategy to select shower-like and starting track events ( §2) over the background of atmospheric muons ( §3).The distributions of observed events are unfolded ( §4) to obtain the energy spectra of both atmospheric ν µ and ν e , taking into account the detector acceptance ( §5).The results are compared with those obtained by other experiments ( §6).The ANTARES telescope is not able to distinguish between neutrino and antineutrino events.Hence, the unfolded spectra are the sum of ν e +ν e and of ν µ +ν µ , averaged over the zenith region 90 • -180 • .

The ANTARES detector and neutrino reconstruction algorithms
The ANTARES telescope [18] is a deep-sea Cherenkov neutrino detector, located 40 km offshore Toulon, France, in the Mediterranean Sea.The detector comprises a three dimensional array of 885 optical modules [19], each one housing a 10-inch photomultiplier tube (PMT).The optical modules are distributed along 12 vertical strings anchored to the sea floor at distances of about 70 m from each other and at a depth of about 2500 m.The detection of light from upward going charged particles is optimised with the PMTs facing 45 • downward.Particles above the Cherenkov threshold induce a coherent radiation emitted in a cone with a characteristic angle θ C 42 • in water.For high-energy muons (E µ > 1 TeV), the contribution of the energy losses due to radiative processes increases linearly with E µ , and the resulting electromagnetic showers produce additional light.Completed in 2008, the telescope aims primarily at the detection of neutrino-induced through-going muons.
The signals induced in the PMT by detected photons are referred to as hits [20].The position, time, and collected charge of the hits are used to reconstruct the direction and energy of events induced by neutrino interactions and atmospheric muons.Trigger conditions based on combinations of local coincidences are applied to identify signals due to physics events over the environmental light background due to 40 K decays and bioluminescence [21].For astronomy studies, atmospheric muons and atmospheric neutrinos constitute the main source of background.
This analysis focuses on events induced by neutrinos whose interaction vertices are contained inside (or close to) the instrumented detector volume.These events include: • ν e charged current (CC) interactions producing electromagnetic and hadronic cascades, and neutral current (NC) interactions of neutrinos of all flavors inducing hadronic cascades.Due to the radiation and nuclear interaction lengths in water, the cascades extend up to a maximum distance of ∼10 m from the interaction vertex, much shorter than the distance between detector strings.These events are thus almost point-like at the scale of the detector and are referred to as shower-like events in the following.
• ν µ CC interactions, with a hadronic cascade near the vertex and a starting muon.Most of these muons are minimum ionising particles, and can travel in water about 4 m per GeV of energy, inducing Cherenkov light over large distances with respect to the interaction vertex position.These events with a cascade and a track are referred to as starting track events in the following.
The track reconstruction algorithm used in off-line ANTARES analyses is called AAFit [22] and it is based on a likelihood fit that exploits a detailed parametrisation of the probability density function for the time of the hits.The algorithm provides the track direction with its estimated angular uncertainty; the number of hits used for the reconstruction; and a quality parameter, referred to as Λ.The track direction and the Λ quality parameter are used to remove the largest fraction of atmospheric muons, which are downward going [23].
An auxiliary algorithm denoted as GridFit [24] searches for tracks in 500 different directions covering the full solid angle.The number of hits compatible with a muon track coming from each direction is evaluated and a likelihood fit is performed.The GridFit method defines quality parameters for reconstructed tracks that are usually used in ANTARES studies to improve the rejection of downward going atmospheric muon events.
Shower-like events occurring in the proximity of the instrumented volume are reconstructed with a dedicated algorithm, TANTRA [25]; this likelihood-based method allows to reconstruct the vertex coordinates, the neutrino direction, and the neutrino energy.A parameter associated to the quality of the reconstructed event, M est , is also provided.
All ANTARES analyses follow a blinding policy to avoid possible biases.The cuts and the selection criteria are studied and optimised on a sample of Monte Carlo (MC) simulated events and only at the end of the full selection chain, these cuts are applied to data.A small sample containing 10% of the real data uniformly distributed over livetime is used to verify the agreement with MC events along the selection.
The simulation chain [26] starts with the generation of the event and comprises the generation of Cherenkov light, the inclusion of the environmental optical background extracted from real data, and the digitisation of the PMT signals following a run-by-run strategy.This strategy accounts for seasonal variations related to biological activities and for inefficiencies due to the ageing of the PMTs and to biofouling [27].
At the end of the full simulation chain, a set of MC files is available for each run of real data, stored in the same format.Simulated files are processed with the same reconstruction algorithms and analysis procedures used for the corresponding data.Monte Carlo neutrino events have been generated in the energy range 10 ≤ E ν ≤ 10 8 GeV, separately for ν e , ν µ and their antineutrinos, and for CC and NC processes.Details on the simulation chain, hadronic model for cross sections, interaction kinematics and parton distribution functions are given in [26].The same MC sample can be differently weighted to reproduce the conventional atmospheric neutrinos, the prompt neutrinos and theoretical astrophysical signals.In the present letter, the atmospheric ν e and ν µ fluxes are represented with the same models used in [28], namely, the conventional component follows the spectrum described in [12], extrapolated at higher energy as in [5], and the prompt contribution as calculated in [15].The MC statistics for the atmospheric neutrino sample corresponds to more than two orders of magnitude than for real data.
Finally, for each data run, a file with simulated atmospheric muons (CRµ) is produced with the MUPAGE package [29,30]; in this case, the equivalent MC livetime corresponds to 1/3 of the real run livetime.

Event selection: signal and background
Data collected from 2007 until the end of 2017 have been used.Only runs without high bioluminescence level have been selected.The total livetime corresponds to 3012 days.The background is almost entirely due to CRµ's: after trigger and reconstruction, the expected signal-to-background rate is ∼10 −4 .The background suppression is organised in three different steps.
An initial preselection of shower-like and starting track events combines information from both the track (AAFit) and shower (TANTRA) reconstruction algorithms, according to the following requirements: i) the direction of the event as reconstructed by AAFit must be upward going (i.e., zenith angle > 90 • ), to reject the largest fraction of CRµ's; ii) the TANTRA's reconstructed neutrino interaction vertex must be contained in a cylindrical volume of axial radius of 300 m and height of 500 m, centred at the centre-of-gravity of the detector modules; iii) the TANTRA estimated angular uncertainty on the event direction must be < 30 • and the quality parameter M est < 1000, to remove poorly reconstructed events.Finally, upward going tracks with AAFit quality parameter Λ > −5.2 are discarded to remove muons generated from neutrino interactions outside the detector volume that could have survived the above cuts.These through going events have already been used in the previous measurement of the ν µ spectrum [10].After this preselection, the MC signal is reduced by a factor of two with respect to the trigger and reconstruction level, with ∼350 survived CRµ's for each atmospheric neutrino candidate.
The second step, following [10,23,31], uses the AAFit quality parameter, Λ.The best compromise to suppress the largest percentage of background while keeping a large enough fraction of signal events is obtained by removing events with Λ ≤ −5.7.After this cut, 25% of the signal survives, with about 30 remaining background events for each atmospheric neutrino.Table 1 summarises the number of events passing the preselection and the Λ cut for each MC sample.The last row shows the events in the data sample, after unblinding.
The final classification of events as signal or background is performed with a Boosted Decision Tree (BDT), defined on a multidimensional parameter space.A BDT is an algorithm that belongs to the family of supervised machine learning techniques.To build the classification function, training samples are necessary.CRµ events generated with MUPAGE constitute the background sample; CC+NC interactions of atmospheric ν e are used as signal.This choice is motivated by the fact that the ν e flavor produces the cleanest case of shower-like events and it is the most difficult channel to measure in neutrino telescopes.
For each CRµ or ν e event, the classifier is trained using the following 15 quantities.From the TANTRA shower algorithm, the reconstructed 1) zenith angle and 2) azimuth angle in the local reference frame; 3) interaction vertex coordinates; 4) quality parameter estimator, M est ; 5) number of detector lines with at least one hit; 6) total number of hits used to reconstruct the event; 7) angular resolution associated to a shower-like event.From the AAFit track-like algorithm, the reconstructed 8) zenith angle and 9) azimuth angle in the local reference frame; 10) track length inside the detector volume; 11) quality parameter estimator, Λ; 12) angular resolution associated to a track-like event.From the GridFit track-like algorithm, 13) the quality parameter; 14) the CRµ veto parameter, a likelihood variable based on time sequence and charge of the hits in different storeys of the detector, causally-connected under the assumption of a downward going, minimum ionizing particle; 15) the number of on-time hits, which assumes that the photons are produced at the Cherenkov angle and arrive at the PMT unscattered.
A ranking of the BDT input variables is derived by counting how often each variable is used to split decision tree nodes, and by weighting each split occurrence by its squared separation gain and by the number of events in the node [32].None of the variables is found to be significantly dominant; the variable with the highest ranking is the TANTRA zenith angle (1) with score 0.12, followed with the GridFit quality (13) with score 0.10; in the last two positions, the estimators of the angular resolution from TANTRA (8), with score 0.04, and that from AAFit (12), with score 0.023.
As shown in Fig. 1, the BDT output is an excellent discriminator between events in the atmospheric ν e and background CRµ samples.The BDT distribution obtained from events induced by atmospheric ν µ CC+NC interactions is also included in the plot.As expected, this distribution resembles that of the ν e signal.The BDT condition >0.33 that removes all CRµ's present in runby-run MC events is chosen as selection criteria.An extrapolation of the BDT distribution tail, assuming a Gaussian shape, yields a conservative extrapolation of (at most) ∼3 background events in the final sample, which are considered in the following.
The last column of Table 1 shows the number of events in the final sample, after BDT cut.As  [12] and prompt neutrinos [15].The magenta histogram is the expected contribution from diffuse cosmic ν's, as parameterised in [28].The orange histogram is the sum of all MC contributions and the black crosses are real data (3012 days livetime), after unblinding.The BDT cut value is denoted with a black arrow.
expected, the neutrino sample is still dominated by atmospheric ν µ producing starting tracks: only ∼10% of the selected events originate from ν e .At 1 TeV, the expected flux ratio is Φ νµ /Φ νe ∼20.

Unfolding procedure and detector acceptance
In order to derive from data the ν e and ν µ energy spectra, an unfolding method is used.The two true distributions are deconvolved from the experimentally measured one, based on the best knowledge of the detector and on assumptions made on the interaction rates of the different neutrino flavors.
In counting experiments, events are grouped into certain regions of phase-space, called bins.The main observable quantities in neutrino telescopes are the neutrino direction and energy, which are measured only with finite precision due to inevitable detector effects.Consequently, an event may be assigned to a wrong bin.In addition, due to the presence of background, only a fraction of the events observed in a given bin originates from the searched signal.
The outcome of the unfolding procedure, folded with the detector acceptance and livetime, results in a spectrum that allows a direct comparison with other experiments.Two major classes of unfolding methods exist: algorithms based on matrix inversion or singular value decomposition, such as the TUnfold [33] algorithm used in this analysis; algorithms based on iterative methods or on the use of Bayes' theorem [34].A Bayesian approach has been used, e.g., by the Super-Kamiokande experiment [8] and in our previous measurement of the ν µ energy spectrum using through-going muons [10].For an overview of the commonly used unfolding algorithms, see also [35,36].
The TUnfold algorithm [33] is a widely tested and validated algorithm in the context of highenergy physics and it can handle one or more background sources.The algorithm allows to estimate the number of events in m bins of a true distribution x j , given an observed distribution of y i in n bins: where each bin has a background contribution b i .A ij is a matrix of probabilities describing the migrations from bin j to any of the n bins.The method, interfaced to the ROOT analysis package [37], uses a least square method with Tikhonov regularisation [38] and an optional constraint to fix the total number of events.The least square minimisation requires a number of degrees of freedom such that n − m > 0, meaning that the data y i have to be measured in finer bins than extracted by the unfolding procedure.
The energy estimated by the TANTRA reconstruction algorithm, E reco , is used to construct the distribution of y i .Events in Fig. 1 with BDT> 0.33 are atmospheric ν µ or ν e , with a contamination of less than a few per mill from CRµ and a ∼1% fraction of cosmic neutrinos; both samples are considered as background.The unfolding method requires a (n × m) matrix for the ν µ and ν e energies, with n bins of E reco and m bins of true energy E ν .Monte Carlo samples allow the construction of: • A e ij , a (6 × 3) matrix obtained with the simulated samples of ν e CC+NC interactions; • A µ ij , a (15 × 5) matrix obtained with the simulated samples of ν µ CC+NC interactions.
The chosen number of bins in (E reco , E ν ) for the two samples provides the highest stability in terms of unfolding results applied on MC samples with the same number of events as real data.In the unfolding procedure, the use of E reco is limited to energies between ∼100 GeV and ∼50 TeV.The lower bound is determined by the fact that our reconstruction algorithm cannot reliably reconstruct neutrino energies below 100 GeV.Above 50 TeV, the event statistics are significantly reduced by the requirement of the containment of the interaction vertex within, or near to, the instrumented volume.In addition, cosmic neutrinos, whose flux suffers large uncertainties, start to be the dominant "background".
Figure 2 shows the distribution of E reco for the ν µ sample (blue) with the bin size used for the construction of the A µ ij matrix.For completeness, the distribution of E reco for the ν e sample (red) is superimposed, although the distribution used for the construction of the A e ij matrix has a different binning.The expected contribution in the sample of cosmic ν's of all flavors is also shown.
Concerning the background terms, it includes an extrapolated contribution of 3 track-like CRµ events.Based on the behaviour of atmospheric muons before the BDT cut, the b CRµ i terms are assumed to affect only the ν µ sample, and uniformly in the E reco range.The background from the cosmic neutrino flux (terms b c i ), assuming equipartition (ν e : ν µ : ν τ ) = (1 : 1 : 1), contributes about equivalently to the ν µ and ν e samples, following the E reco distribution shown in Fig. 2.
The unfolding procedure is divided into two parts, whose structure is the same: when focused on the ν µ distribution, the background corresponds to b CRµ ) is equal to the total number of events.The events b µ i and b e i are assumed to be produced with the fluxes given in [12], as in the default Monte Carlo simulation, with free normalization.Possible variations in their spectral indexes are accounted for in the treatment of systematic effects (see §5).
Table 2 presents the information on the unfolded energy in 5 (3) bins for the ν µ (ν e ) sample.The first two columns contain the energy range of the corresponding bin and the weighted central value of the neutrino energy bin, calculated taking into account the steep decrease of the energy spectrum and the detector response.The third column shows the unfolded number of data events as obtained by the algorithm.

The unfolded energy spectrum
To transform the unfolded number of events, N evt , given in Table 2 into a differential energy flux in the proper units (GeV −1 cm −2 s −1 sr −1 ), the following steps are required: i) divide each bin by the livetime of 3012 days, obtaining the event rate integrated in the log 10 of the neutrino energy over the bin; ii) divide by the width of the bin (0.54 for ν µ and 0.9 for ν e ); then, transform the dN evt d log 10 Eν distribution into the dN evt dEν one; iii) divide by the integrated value of the observation solid angle, i.e., 2π sr; iv) divide by the detector effective area, A ef f (E ν ), averaged over the zenith angle of the atmospheric neutrino distribution.
The effective area is the figure of merit for a neutrino telescope, representing the size of a 100% efficient hypothetical target that the detector offers to a certain simulated neutrino flux.It is calculated as where N sel (E ν ) and N gen (E ν ) are, respectively, the number of selected and generated events of a given neutrino energy E ν in the generation volume V gen ; ρ and N A are the matter density and the Avogadro's number; σ(E ν ) is the neutrino cross section; P Earth (E ν ) is the probability of the neutrino to traverse the Earth without being absorbed.Above 100 GeV, there are no corrections needed for oscillation effects.Figure 3 shows the effective area obtained from the selection of events described in this work.
The fourth column of  GeV ; the number of unfolded events assigned to the bin, N evt ; the differential flux (times E 2 ν ) computed in the centre of the bin, E 2 ν Φν , in units of GeV cm −2 s −1 sr −1 ; the statistical error; and the total systematic uncertainty.
The result of the unfolding process depends on the MC simulation via the construction of the response matrix.In turn, the simulation depends on a number of parameters with associated uncertainties.The effects inducing systematic uncertainties on the measurement of the ν µ flux using through-going events have been extensively described in [10]; these also affect the present analysis.The impact is estimated by varying each time only one of the following parameters, and producing dedicated simulation datasets: • Overall sensitivity of the optical modules, changed by +10% and −10%.This includes the uncertainty on the conversion of a photon into a photoelectron as well as the angular dependence of the light collection efficiency of each optical module.
• The uncertainties on water properties, by scaling up and down by 10% the absorption length of light in water with respect to the nominal value.
• The uncertainties related to the neutrino fluxes used in the default response matrix of the unfolding procedures, including a slope change of ±0.1 in the spectral index, independently for ν e and ν µ .
Each modified MC sample was then used as pseudo-data to construct a new response matrix, used for unfolding.The deviation in the content of each E ν bin from the spectrum obtained with the default response matrix, A e ij or A µ ij , corresponds to the systematic uncertainty associated with the parameter variation.For each energy bin, the total uncertainty is computed as the quadratic sum of each contribution, and the resulting value is reported in the last column of Table 2.

Results and conclusions
Figure 4 shows the (ν e +ν e ) and (ν µ +ν µ ) fluxes measured in this work, together with the results from previous experiments.Our unfolded atmospheric neutrino spectra, whose statistical errors are largely dominant over the systematic ones, are about 20% below the most recent computations using the SIBYLL-2.3chadronic interaction model [17,39].
The measurement of the electron neutrino flux at high energy is challenging, because very large detectors are needed to collect sufficient statistics, and due to large systematic uncertainties.Each measurement of IceCube-DeepCore [4] and IceCube [5] rely on about 200 interacting ν e in the polar ice medium.The present measurement is performed in seawater, under completely different environmental conditions and systematic uncertainties, yielding consistent results with the ones obtained in polar ice.During a livetime of 3012 days, about 130 ν e interactions have been reconstructed within the instrumented ANTARES volume.The statistics of the ν e sample is not sufficient to test models above a few tens of TeV, where a significant cosmic flux is present and the transition from the conventional to the prompt flux is expected.Below 100 GeV, the PMT density of the ANTARES detector is insufficient to reconstruct a significant number of events.
Concerning the unfolded ν µ flux, our previous measurement [10] with a sample of through-going events generated by neutrino interactions external to the instrumented volume almost superimposed the SIBYLL-2.3cmodel.The present measurement is based on a totally independent data sample, provided by neutrinos whose reconstructed interaction vertex is inside (or nearby) the instrumented volume of the detector.The present unfolded flux is more close to that of IceCube with 40 lines [6] and is about 25% below both the flux reported in our previous measurement and the one reported by IceCube using 59 strings [7], although consistent within errors.Measured energy spectra of the atmospheric ν e and ν µ using shower-like and starting track events in the ANTARES neutrino telescope (black).The measurements by other experiments (Frejus [9], AMANDA-II [2], IceCube [6,7,4,5], and Super-Kamiokande [8]), as well as the previous ν µ flux measurement using a different ANTARES data sample [10], are also reported.The vertical error bars include all statistical and systematic uncertainties.

Figure 1 :
Figure 1: BDT output for events passing the preselection + Λ cut.The histograms correspond to different MC samples: training CRµ (green), training atmospheric νe (red), atmospheric νµ (blue).The νµ events are not used for BDT training.The green line corresponds to a Gaussian extrapolation of the CRµ histogram.Both νe and νµ fluxes include conventional[12] and prompt neutrinos[15].The magenta histogram is the expected contribution from diffuse cosmic ν's, as parameterised in[28].The orange histogram is the sum of all MC contributions and the black crosses are real data (3012 days livetime), after unblinding.The BDT cut value is denoted with a black arrow.

i , b c i / 2
and to the ν e fraction b e i .When focused on the ν e distribution, the background corresponds to b c i /2 and the muon neutrinos b µ i .The algorithm assumes that i (b µ i + b e i + b CRµ i + b c i

Figure 2 :
Figure 2: Distribution of E reco with the binning used for the construction of the response matrix for the ν µ sample (blue histogram).The red histogram, with the same binning, refers to the ν e sample.The magenta histogram is the expected contribution from a cosmic neutrino flux, as estimated in [28], while the orange histogram includes the sum of all MC contributions.The black crosses correspond to real data.Events in the shaded region are used for unfolding.

Figure 3 :
Figure3: Effective area of the ANTARES neutrino telescope for the events with a vertex inside the instrumented volume and selected by the analysis cuts described in this work: ν e CC+NC (red line), ν µ CC+NC (blue line).The black solid line is the sum of all the interaction channels and neutrino flavors.

Figure 4 :
Figure4: Measured energy spectra of the atmospheric ν e and ν µ using shower-like and starting track events in the ANTARES neutrino telescope (black).The measurements by other experiments (Frejus[9], AMANDA-II[2], IceCube[6,7,4,5], and Super-Kamiokande[8]), as well as the previous ν µ flux measurement using a different ANTARES data sample[10], are also reported.The vertical error bars include all statistical and systematic uncertainties.

Table 1 :
[28]er of events in different Monte Carlo samples surviving the preselection and the cut on track quality parameter Λ (second column) and the final BDT cut (third column).The last row shows the number of data events in 3012 days of livetime.The expectation for cosmic neutrinos is computed assuming the diffuse flux presented in[28].
Table 2 presents the differential flux obtained with the overall procedure.The reported statistical error is determined by the TUnfold method.∆ log E ν log E ν N evt

Table 2 :
The table shows, from the first to last column: ∆ log Eν = log 10