Neutrino Mass Ordering at DUNE: an Extra $\nu$-Bonus

We study the possibility of extracting the neutrino mass ordering at the future Deep Underground Neutrino Experiment using atmospheric neutrinos, which will be available before the muon neutrino beam starts being perational. The large statistics of the atmospheric muon neutrino and antineutrino samples at the far detector, together with the baselines of thousands of kilometers that these atmospheric (anti)neutrinos travel, provide the ideal ingredients to extract the neutrino mass ordering via matter effects in the neutrino propagation through the Earth. Crucially, muon capture by Argon provides excellent charge-tagging, allowing to disentangle the neutrino and antineutrino signature. This is a critical extra benefit of having a Liquid Argon Time Projection Chamber as far detector, that could render a $4\sigma$ extraction of the mass ordering after ten years of exposure.


I. INTRODUCTION
Neutrino oscillation experiments imply the first departure from the Standard Model (SM) of Particle Physics, as they have found overwhelming evidence for the existence of neutrino masses. Despite the accuracy they provide on the neutrino oscillation parameters -which is of the order of the percent level [1] -the sign of the atmospheric mass splitting, ∆m 2 31 , and the value of the CP violating phase δ remain both unknown. The sign of ∆m 2 31 originates two possible scenarios, normal (NO) or inverted ordering (IO) [2,3]. The sensitivity to the neutrino mass spectrum at oscillation experiments is mostly coming from the presence of matter effects [4][5][6][7][8][9][10][11][12][13] in the neutrino and antineutrino propagation. In the normal (inverted) mass ordering scenario, the neutrino flavor transition probabilities will get enhanced (suppressed), while in the case of antineutrino propagation the opposite happens and the antineutrino flavor transition probabilities will get suppressed (enhanced) in the normal (inverted) mass ordering scenario. At long-baseline accelerator experiments, matter effects, and consequently the sensitivity to the mass ordering, increase with the baseline, while these effects will be negligible at short-baseline and medium-baseline experiments. Despite this, when extracting both the mass ordering and the CP violating phase from results of long-baseline facilities, knowledge on the mixing angle θ 13 in vacuum is required. Short-and medium-baseline experiments at reactors have been fundamental to establish strong constraints on such angle. Despite the fact that the neutrino mass ordering remains unknown, current oscillation data are mildly favouring the normal ordering scenario. The authors of Refs. [1,3] have reported a global preference for normal ordering at the level of 2.7σ from all long-baseline accelerator and * chternes@ific.uv.es † gariazzo@ific.uv.es ‡ rashajmu@alumni.uv.es § omena@ific.uv.es ¶ sorel@ific.uv.es * * mariam@ific.uv.es short-baseline reactor data (i.e. T2K, NOνA, K2K, MINOS, Daya Bay, RENO and Double Chooz).
The future long baseline facility DUNE (Deep Underground Neutrino Experiment) [14][15][16] aims to extract the sign of the atmospheric mass splitting and the CP violating phase δ through the golden channels ν µ → ν e and ν µ →ν e , the same channels exploited by the current T2K [17,18] and NOνA [19,20] experiments. However, both quantities can also be extracted using atmospheric neutrino beams 1 . Indeed, the idea of using atmospheric neutrino fluxes to distinguish the type of mass ordering is well-known in the literature since a long time [23,24]. These pioneer studies focused mostly on muon calorimeter detectors, such as MONOLITH [25], MINOS [26] or INO 2 in which the muon charge can be determined, see also Refs. [27][28][29][30][31][32][33][34][35][36][37][38][39][40]. Furthermore, in the absence of a charged current event by event final muon charge discrimination, the addition of the atmospheric neutrino oscillation data to the analysis performed in Ref. [1] improves the preference for normal ordering to the level of 3.4σ, mostly due to the Super-Kamiokande ν µ → ν e data sample [41], where the separation among ν e andν e events is done statistically.
In this manuscript, we exploit the atmospheric neutrino signatures at the DUNE detector, a Liquid Argon Time Projection Chamber (LArTPC). Despite this detection technology, in the absence of a magnetic field, does not allow for a charge identification of the final lepton state, one can make use of a particular event topology available in Argon detectors: muon capture. This bonus process will provide a clean measurement of the muon charge, that will considerably improve the capabilities of DUNE to perform mass ordering measurements with atmospheric neutrinos. Notice that the advantage is twofold, as (i) measurements of the mass ordering could be available before the beam starts, and (ii) the combination with the beam information will notably enhance the expected sensitivity reach. We shall show that muon capture events could greatly enhance the sensitivity to the mass ordering from atmospheric neutrinos only. For an earlier, and preliminary, appraisal of the neutrino mass ordering sensitivity in DUNE using atmospheric neutrinos, including a statistical discrimination between neutrinos and antineutrinos, see Ref. [58]. The latter work was largely based on previous studies in the framework of the LBNE project, see [59].
The structure of the paper is as follows: In Sec. II, we describe the oscillations of atmospheric neutrinos and the matter effects they undergo. Next, in Sec. III, we discuss the simulation of the neutrino event rates at the DUNE far detector and how the muon capture comes into play. Sec. IV contains the description of the statistical method and the main results obtained in this study. Our final remarks are presented in Sec. V.

II. MATTER EFFECTS AND ATMOSPHERIC NEUTRINOS
In atmospheric neutrino experiments, the size of matter effects is given by the effective mixing angle θ 13 in matter, which leads to the golden channel transitions ν µ → ν e , ν e → ν µ ,ν µ →ν e andν e →ν µ and reads, within the simple two-flavor mixing framework, as sin 2 2θ m 13 = sin 2 2θ 13 where the minus (plus) sign refers to neutrinos (antineutrinos). The matter potential is given by A = 2 √ 2G F N e E and N e is the electron number density in the Earth interior. Consequently, matter effects will enhance (deplete) the neutrino (antineutrino) oscillation probabilities P (ν µ → ν e ) and P (ν e → ν µ ) (P (ν µ →ν e ) and P (ν e →ν µ )) if the mass ordering is normal. When the resonance condition is satisfied, matter effects are expected to have their largest contribution. In the case of atmospheric neutrinos, which travel distances of several thousand of kilometers, and for ∆m 2 31 ∼ 2.5 × 10 −3 eV 2 [1], the resonance condition will take place at neutrino energies ∼ 3−8 GeV, depending on the precise value of N e in the Earth's interior.
Matter effects are also present in the muon disappearance channels P (ν µ → ν µ ) and P (ν µ →ν µ ), relevant for both long-baseline and atmospheric neutrino beams. In the simplified case of a constant matter density, the disappearance probability at terrestrial baselines 3 is given by − sin 2 θ m 13 sin 2 2θ 23 × sin 2 1.27 where θ m 13 is that of Eq. (4) The muon survival probabilities will be suppressed (enhanced) if the ordering is normal (inverted), so the effect 3 For an expansion with solar mixing effects included, see Ref. [60].
is opposite to the one present in the ν e → ν µ oscillation channel. Therefore, when dealing with atmospheric neutrino beams, since there is an irreducible muon neutrino background from ν e → ν µ oscillations, the size of the matter effects will be reduced. The distance L traveled through the Earth by these atmospheric neutrino beams is fixed by their arrival zenith angle θ z (with cos θ z = 1 for vertical down-going neutrinos and cos θ z = −1 for Left panels: survival probability P (νµ → νµ) as a function of the neutrino energy E and the cosine of the zenith angle, cos θz, for normal (inverted) ordering in the top (bottom) line. Right panels: same as in the left panels, but for the antineutrino channel.
vertical up-going neutrinos), with R ⊕ the Earth's radius and h 15 km the neutrino production distance from the Earth's surface. The dependence of the survival probabilities P (ν µ → ν µ ) and P (ν µ →ν µ ) on the neutrino energy E and the cosine of the zenith angle, cos θ z , is shown in the left and right panels of Fig. 1 for normal and inverted ordering (top and bottom figures). Notice that, in the case of normal ordering, the resonance takes place at the aforementioned energies (3-8 GeV) for almost vertical up-going neutrinos, −1 < cos θ z < −0.8, while for the inverted ordering, such a resonant enhancement in the transition probabilities will take place in the antineutrino channel instead. Therefore, even if both the angular and the energy resolution of the detector should be optimal, the key ingredient to disentangle matter effects (and, ultimately, the neutrino mass ordering) is to have a detector with muon charge tagging, generally achieved with a magnetized detector. However, as we shall shortly see, LArTPCs allow for such a possibility without the need of a magnetic field.

III. ATMOSPHERIC NEUTRINO EVENTS IN DUNE: MUON CAPTURE IN ARGON
Our statistical analyses will deal with three possible fully contained event samples at atmospheric neutrino detectors: µ − -like events that undergo muon capture (N cap i,j,µ ), the rest of the muons and all of the antimuons that undergo muon decay (N rest i,j,µ ), and e-like events (N i,j,e ) 4 .
Let us start with the µ − -like contained events produced by the interactions of atmospheric up-going neutrinos in the LArTPC DUNE detector. In a LArTPC, both ionization charge [61,62] and scintillation light [63] information can be used to infer the neutrino/antineutrino content in a muon neutrino beam. This is possible by exploiting the signature of µ − capture on Argon nuclei, only available for contained events. In argon, the effective µ − lifetime resulting from the competing decay and nuclear capture processes is given by: where τ cap is the lifetime of the capture process, Q = 0.988 is the Huff correction factor [64] and τ free = 2197.0 ns [65] is the muon lifetime in vacuum. The resulting µ − capture fraction is then given by: The most precise determination of the µ − lifetime in argon was obtained in [66], τ = (616.9±6.7) ns, resulting in This measurement is fully compatible with the earlier measurement of τ = (606 ± 29) ns in [64] and the preliminary result from LArIAT of τ = (626 ± 48) ns [67]. In our analysis, we use the central value and uncertainty in Eq. (8). For our sensitivity estimates, we also assume a 100% efficiency for tagging Michel electrons and positrons from µ ± decays at rest, as done in [59]. Any tagging inefficiency would cause decay events to be mis-interpreted as capture events, and should therefore be avoided for optimal muon neutrino/antineutrino separation. We consider this approximation to be sufficient for the purposes of this feasibility study. Efficiency estimates using detailed DUNE simulations are not yet publicly available. Still, early data from ICARUS [68] and LArIAT [67] have already shown that the Michel electron tagging efficiency can reach values close to unity in LArTPC detectors using either charge or light information. In any case, any Michel electron tagging inefficiency smaller than (1− cap ) 28% will have a sub-dominant contribution to the mixing of muon neutrino and muon antineutrino stopping samples in our analysis, compared to the effect of µ − tracks that do not capture and decay.
Therefore, it appears possible to select a statistically significant, highly pure, sample of µ − -like atmospheric neutrino interactions, with an identification efficiency of cap , as given in Eq. (8). The number of muon-like contained events in the i-th neutrino energy (E r ) and j-th cosine of the zenith angle (c r,ν ) bin (both reconstructed quantities) reads as Er,i dE r,ν cr,ν,j+1 where dφ ν 's are the atmospheric neutrino differential fluxes, σ CC is the CC neutrino cross sections in Argon, N T is the number of available targets, V det is the total volume of the detector, V µ is the effective detector volume and t is the exposure time. Finally, R µ e (E r,ν , E ν ) and R µ θ (θ r,ν , θ ν ) account for the energy and angular smearing.
The µ − -like contained events that undergo muon capture are given by while the remaining muon-like events are given by In Fig. 2 we show an example for the expected number of events, fixing the value the oscillation parameters to the ones in Tab. I and the atmospheric mixing angle to sin 2 θ 23 = 0.547. In the first two panels, we plot the capture and decay events separately, respectively, while showing the combination of both samples in the right panel. Note that when trying to reconstruct the right panel there can be more degeneracies among parameters in the analysis than when fitting the two sets independently and, therefore, one expects to obtain stronger results. This is indeed the case, as we will see below.
In the case of electrons, the number of e-like events in the i-th and j-th bin in (E r , c r,ν ) reads as Er,i dE r,ν cr,ν,j+1 As previously stated, we just consider one electronlike event sample N i,j,e which is computed as the sum of N i,j,e − and N i,j,e + .
Regarding the atmospheric electron and muon (anti) neutrino fluxes, for the differential fluxes dφν α dEν dΩ that appear in Eqs. (9) and (12), we use the results from Ref. [69], albeit very similar numbers would have been obtained using the fluxes from Refs. [70][71][72]. We shall comment in the following section on the errors on these atmospheric neutrino fluxes, that have been properly added to other sources of systematic uncertainties in our numerical studies.
The cross sections for muon and electron (anti)neutrino interactions on Argon nuclei in the 0-10 GeV neutrino energy range have been simulated by means of the GENIE Monte Carlo neutrino event generator [73]. GENIE is extensively used by the neutrino physics community and by the DUNE Collaboration in particular. As our cross section model, we use the total charged-current (anti)neutrino cross sections provided by GENIE version 2.12.10 on 40 Ar nuclei (18 protons and 22 neutrons). The model accounts for a comprehensive list of interaction processes, including quasi-elastic scattering, baryon resonance production, coherent pion production in neutrino-nucleus scattering and deep inelastic scattering. Nuclear effects affecting total cross sections are included. Final state hadronic interactions occurring within the Argon target nucleus are not simulated, but indirectly accounted for via our assumed energy and angular resolution functions.
To compute the effective volume fraction V µ /V det in Eq. (9) for contained muon events, we have approximated the DUNE detector to be made of four independent modules with approximately 13 kton of LAr active mass each, each of them assumed to have an elliptical cylindrical shape of 12 m height and major and minor axis of a = 29 m and b = 7.25 m, respectively. For the calculation of the effective volume we have taken into account the muon range in Argon R µ (E µ ), which depends on the lepton energy. Conservatively, we have also computed the number of µ + -like events restricting ourselves to the contained topology. This assumption eases the comparison with respect to the case in which no flavor tagging is available and ensures good energy reconstruction for the full muon like event sample.
As for the energy and angular smearing inherent to reconstruction processes and final state hadronic interactions within Argon nuclei, R µ e (E r,ν , E ν ) and R µ θ (θ r,ν , θ ν ) in Eq. (9), and R e e (E r,ν , E ν ) and R e θ (θ r,ν , θ ν ) in Eq. (12), are taken to be gaussian functions. The assumed gaus- sian widths σ E /E and σ θ for ν e ,ν e , ν µ andν µ chargedcurrent interactions on argon are shown in Fig. 3. We use the dashed curves in the figure to parametrize the resolutions as functions of the neutrino energy E ν , according to: The numerical values for the parameters in Eq. 13 are reported in Appendix A. The resolution functions were obtained via fast Monte-Carlo simulations as follows, similarly to what was done in [58,59]. First, large samples of mono-energetic neutrino-argon interactions are simulated with GENIE, for the various neutrino flavors (ν e ,ν e , ν µ andν µ ) and for the relevant neutrino energy range 0.5-8 GeV. The GENIE simulation includes nuclear effects. Second, for each event, each final state particle exiting the nucleus has its kinetic energy and angular direction smeared according to the assumptions described in [59]. The relative energy resolutions are taken to be 1%/ √ E e + 1%, 3% and 30%/ √ E had for electrons, muons and hadrons, respectively, where E e and E had are expressed in GeV. The absolute angular resolutions are taken to be 1 • , 1 • and 10 • for the same three final state particle categories and for all energies. Third, the incoming neutrino energy and direction of each interaction is reconstructed as follows: where K r,l and K r,h are the reconstructed charged lepton and hadron kinetic energies, m l is the charged lepton mass, the sum h is intended over all final state hadrons, and p r,ν ≡ p r,l + h p r,h is the reconstructed 3-momentum of the incoming neutrino, where the true neutrino direction is defined along the z axis. Fourth, histograms of the reconstructed neutrino energy and direction are obtained for each (neutrino flavor, neutrino energy) simulated data sample. Fifth, σ E /E ν and σ θ for each data sample are obtained from a gaussian fit to the energy histogram and from the mean of the angle histogram, respectively. The resolution functions are shown in Fig. 13 for each sample via marker symbols. Sixth, the energy dependence of the resolutions functions is parametrized according to Eq. (13).
The behavior of the resolution functions in Fig. 13 can be easily understood. The main effect is that both σ E /E and σ θ improve noticeably as the neutrino energy increases. For σ E /E, this is a direct consequence of the relative energy resolutions assumed for electrons and (especially) hadrons, improving as the particle energies increase. For σ θ , this is due to the Fermi momentum of the target nucleon, whose angular smearing effect is more important at low neutrino energies. A second, smaller, effect can also be appreciated in Fig. 13, namely that antineutrino resolutions are slightly better than neutrino ones, both for σ E /E and σ θ . On the one hand, this is due to the fact that the average inelasticity (or energy fraction carried away by final state hadrons) is somewhat lower in antineutrino interactions [74], and on the other because hadron resolutions are substantially worse than charged lepton ones. An even smaller difference can be appreciated between the relative energy resolutions of electron and muon antineutrinos of the same energy. In this case, electron antineutrino energy resolutions are slightly better, because of the better assumed accuracy in reconstructing electron energy (1%/ √ E e +1%) compared to muon energy (3%).
Our energy resolution assumptions in Fig. 3 are similar to the ones in [58,59], that use similar methodologies and assumptions. They are qualitatively similar also to the ones obtained in more recent studies, see Refs. [74,75]. On the other hand, we are not aware of other neutrino angular resolutions studies in LArTPCs to compare our  findings with.

IV. ANALYSES AND RESULTS
Here, we describe the statistical analysis and how we extract the sensitivity to the neutrino mass ordering. In order to emphasize the impact of the muon capture in Argon, we present two possible analyses. The first case will assume that no charge identification is possible. Then, we will focus on the extra bonus that the muon capture in Argon process provides.
In the following, we define a fiducial mass ordering, true ordering (TO), in order to generate mock data. Then, we try to reconstruct the event rates using the wrong ordering (WO) assumption. Although there is some preference for normal neutrino mass ordering, as previously stated, we shall study also the case of inverted ordering as TO.
We use Eqs. (9) and (12) to generate our mock data, using the oscillation parameters from Tab. I and assuming a 400 kt·yr exposure. We will present our results as a function of the atmospheric angle θ 23 . Therefore, there is no fixed value for this angle in the table. Notice that, since our main sensitivity comes from the ν µ → ν µ channel, the effects of the CP violating phase δ are negligible and therefore we set δ = 0, finding very similar results for other values of the CP phase.
Next, we try to reconstruct the event rates following the two methods mentioned above. Before presenting our results, let us discuss our treatment of systematic uncertainties. We consider several sources of systematic uncertainties in our analyses, coming from the fact that we do not have a perfect knowledge of the atmospheric flux and detector response. In particular, we include an overall rate normalization error accounting for both flux normalization and detector efficiency uncertainties, an error on the ν/ν atmospheric flux ratio, and an error on the ν µ /ν e atmospheric flux ratio. We follow Ref. [59] and assume a 15%, 5%, and 2% gaussian error on these three quantities, respectively. We have verified that adding a systematic on the spectral index of the neutrino flux would have a negligible effect. As explained in the previous section, we also add an uncertainty on cap , see Eq. (8). Apart from the systematic uncertainties, we also marginalize over the oscillation parameters ∆m 2 31 , sin 2 θ 13 and sin 2 θ 23 within their current 3σ ranges for both orderings, namely |∆m 2 31 | ∈ [2.31, 2.60] × 10 −3 eV 2 , sin 2 θ 13 ∈ [0.0196, 0.0244] and sin 2 θ 23 ∈ [0.455, 0.599]. It is well known that the solar parameters do not have big effects in atmospheric neutrino oscillations, hence, they are fixed to their best-fit values throughout the analysis.
Method A: Analysis without muon capture tagging In this case, muons and antimuons cannot be distinguished. We therefore build a χ 2 function in the following way: A (sin 2 θ true 23 ) = min sys,∆m 2 31 ,θ13,θ23 We use a Poissonian χ 2 , which for muons is is the sum of the muon and antimuon contributions. The same formula applies to χ 2 e − +e + , with the replacement µ → e. The results of our analysis with method A are shown as red curves in Fig. 4. Note that the sensitivity ranges between 1.5 and 3.5σ approximately when normal ordering is the TO (solid lines) and between 1.5 and 2σ for true inverted ordering (dashed line).

Method B: Analysis with muon capture tagging
In this other strategy, we use use muon capture to distinguish ∼72% of the muons from antimuons. Therefore, this time our χ 2 function contains three terms, namely χ 2 B (sin 2 θ true 23 ) = min sys,∆m 2 31 ,θ13,θ23 The electron term is the same as for method A, while the other two terms, corresponding to the events with muon capture (cap) and all other events (rest), are given by where X ∈ {cap, rest}, see Eqs. (10) and (11) . The results of the analysis with muon capture are shown in Fig. 4 by the blue curves. As before, true normal ordering is shown as a solid line, while the case of true inverted ordering is represented by a dashed line. The grey band in the figure represents the current 1σ allowed region for sin 2 θ 23 . Note how the sensitivity to the mass ordering is now at the 2.5-4σ level, implying an important improvement with respect to the results obtained with method A. In particular, for the current best-fit point [1] we find that, using atmospheric neutrinos with muon capture, DUNE could measure the neutrino mass ordering at the 3.5σ level. Our method B results can also be compared with the results in the DUNE Conceptual Design Report [58], where a similar sensitivity reach and dependence on sin 2 θ 23 were obtained. Compared to [58], however, our results more clearly highlight the importance of the muon capture tag.

V. CONCLUSIONS
We have explored the advantages of muon capture on Argon nuclei, a process that improves the sensitivity to the neutrino mass ordering using atmospheric neutrino events at the Liquid Argon Time Projection Chamber DUNE far detector. This is a very relevant result, since it comes without any extra cost. Furthermore, it can be combined with DUNE beam neutrino results, allowing for an enhancement in the total sensitivity to the mass ordering determination. It is important to notice that our results are applicable to any experiment using Argon. In the case of accelerator-based neutrinos, where significant ν µ contamination exists in theν µ beam, statistical neutrino/antineutrino separation based on muon capture could also be used to enhance DUNE oscillation sensitivities.   Figure 3 shows our estimated neutrino energy σ E /E ν (top) and neutrino angle σ θ (bottom) resolutions as a function of neutrino energy E ν , for charged-current neutrino interactions on argon. The resolutions are parametrized according to Eq. 13. The A and B parameters describing the relative neutrino energy resolution (in percent), and the C and D parameters describing the neutrino angle resolution (in degrees), are given in Tabs. II and III, respectively. The parameters are given separately for each neutrino flavor: ν e ,ν e , ν µ andν µ . The parameters in Eq. 13 assume that E ν is expressed in GeV.