Matter effects in neutrino visible decay at future long-baseline experiments

Abstract

Neutrino visible decay in the presence of matter is re-evaluated. We study these effects in two future long-baseline experiments where matter effects are relevant: DUNE (1300 km) and a hypothetical beam aimed towards ANDES (7650 km). We find that matter effects are negligible for the visible component of neutrino decay at DUNE, being much more relevant at ANDES. We perform a detailed simulation of DUNE, considering $\nu _\mu$ disappearance and $\nu _e$ appearance channels, for both FHC and RHC modes. The sensitivity to the decay constant $\alpha _3$ can be as low as $2\times 10^{-6}$ eV$^2$ at 90% C.L., depending on the neutrino masses and type of coupling. We also show the impact of neutrino decay in the determination of $\theta _{23}$ and $\delta _\mathrm{CP}$, and find that the best-fit value of $\theta _{23}$ can move from a true value at the lower octant towards the higher octant.

1 Introduction

The neutrino oscillation phenomenon is the only firm evidence of physics beyond the Standard Model. It is then of interest to study at depth if neutrino oscillations could give further information regarding a more complete description of Nature. In this case, one would hope that measurements in neutrino experiments would eventually deviate from the expectations of the standard neutrino oscillation paradigm. Such kind of deviations have been extensively studied in the literature, notable examples being decoherence [1,2,3,4,5,6,7,8,9,10,11,12,13] and non-standard interactions [14,15,16,17,18,19,20,21].

Light neutrino decay constitutes a third way of producing changes in neutrino oscillation experiments. In particular, models where the neutrino decay is due to its coupling to a massless scalar, called a Majoron, have been analyzed in detail by many authors [11, 22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49]. In these kind of models, neutrinos and Majorons can have two types of couplings, scalar ($g_s$) or pseudoscalar ($g_p$):

$$\begin{aligned} \mathcal {L}_\mathrm{int}=\frac{(g_s)_{ij}}{2}{\bar{\nu }}_i\nu _jJ+i\frac{(g_p)_{ij}}{2}{\bar{\nu }}_i\gamma _5\nu _j J. \end{aligned}$$

(1)

The decay widths generated by such couplings are well known and can be found, for instance, in [28, 46]. Relevant constraints to these couplings can be found in the literature [50,51,52,53,54,55].

In general, most works consider that the light neutrinos can decay into an unobservable product, such as a light sterile neutrino or an undetectable active neutrino. This is the so-called invisible decay (ID). However, the possibility of observing the decay into active products, referred to as visible decay (VD), has also been studied in the past [28], and has recently been reconsidered in the context of current long-baseline experiments [46]. The trait of this scenario is an excess of neutrinos of lower energies, which are the decay products of the neutrinos from the original flux. This was followed by an analysis in the framework of future facilities [47], which included matter effects. The question of how to properly include matter effects within the oscillation with decay scenario was also addressed in [56], where they considered the decay of a heavier sterile neutrino.

In this paper we take a different approach in the inclusion of matter effects, following closely the prescription outlined in [28]. The role of matter in the visible component of the decay process is analyzed by comparing two different baselines (corresponding to DUNE [57] and ANDES [58]), which correspond to different matter densities. In addition, we make a systematic study of the sensitivity to the decay parameters at the DUNE experiment, including all oscillation channels and the two modes of operation (FHC and RHC).

This paper goes as follows. On Sect. 2, we develop the theoretical framework to be used for describing neutrino decay in matter. On Sect. 3, we assess the impact of matter effects in neutrino decay within the aforementioned baselines. Finally, on Sect. 4 we discuss the details of our DUNE simulation and statistical analysis and evaluate the sensitivity of this experiment to the decay constant $\alpha _3$. On that Section we also study the impact of neutrino decay on oscillation fits.

2 Neutrino flux including matter effects and decay

The main idea behind our procedure relies on carefully identifying the basis where each process takes place. As is common, one defines the interaction (or flavour) basis as the one where the charged lepton mass matrix is diagonal, which is also the basis where neutrino interaction states ($\nu _e$, $\nu _\mu$, $\nu _\tau$) are produced. However, on this basis the neutrino mass matrix is not diagonal, meaning that the interaction states are a superposition of mass eigenstates. This mass matrix is diagonalized by the PMNS matrix $(U_0)_{\alpha i}$, where $\alpha =e,\,\mu ,\,\tau$, and $i=1,2,3$, which connects both states. We refer to the basis where the neutrino mass matrix is diagonal as the mass basis, and it is in this basis where the decay widths $\Gamma _i$ of each neutrino are defined.

As is well known, when neutrinos travel through matter, the effective Hamiltonian H is not diagonal, neither in the interaction nor mass bases. For our case, in terms of $\alpha _i=E\, \Gamma _i$, we have on the interaction basis:

$$\begin{aligned} H= & {} \frac{1}{2E}\,U_0\left( \begin{array}{ccc} m_1^2 &{} &{} \\ &{} m_2^2-i\,\alpha _2 &{} \\ &{} &{} m_3^3-i\,\alpha _3 \end{array}\right) U_0^\dagger \nonumber \\&\quad +\left( \begin{array}{ccc} \sqrt{2}G_F N_e &{} &{} \\ &{} 0 &{} \\ &{} &{} 0 \end{array}\right) , \end{aligned}$$

(2)

where E is the neutrino energy, $G_F$ is Fermi’s constant, and $N_e$ is the electron density in matter. We have assumed normal ordering with the lightest neutrino being stable.

This structure motivates the introduction of a new basis [59,60,61,62], which we refer to as the matter basis, where the effective Hamiltonian is diagonal. Given the presence of the decay widths, H is diagonalized by non-unitary ${\tilde{U}}_{\alpha I}$ matrices ($I={\tilde{1}},\,{\tilde{2}},\,{\tilde{3}}$):

$$\begin{aligned} {\tilde{U}}^{-1}H\, {\tilde{U}} = H^\mathrm{diag}, \end{aligned}$$

(3)

where $H^\mathrm{diag}$ has complex eigenvalues:

$$\begin{aligned} {\tilde{m}}_I^2-i\,{{\tilde{\alpha }}}_I=2E\,(H^\mathrm{diag})_{II} \end{aligned}$$

(4)

It is important to take into account that, in this new basis, all eigenstates can have a non-zero decay width. For instance, if we start with only $\alpha _3$ different from zero, we can obtain non-vanishing $\tilde{\alpha }_1$, $\tilde{\alpha }_2$ and $\tilde{\alpha }_3$.

Let us now understand how the combination of matter effects and decay affect a neutrino flux. We denote the flux arriving at the detector in absence of flavour transitions by $d\Phi _\alpha ^{(r)}/dE_{\alpha }$, where $\alpha$ refers to the flavour, $r=(+,\,-)$ is the helicity, and $E_\alpha$ is the energy, which are determined when the neutrino is produced. In the presence of oscillations and decay, the neutrino flux of flavour $\beta$, helicity s and energy $E_\beta$, can be calculated with:

$$\begin{aligned} \frac{d \Phi ^{(s)}_\beta }{dE_\beta }=\int P_\mathrm{dec}\left( \nu _\alpha ^{(r)}\rightarrow \nu _\beta ^{(s)}\right) \frac{d\Phi ^{(r)}_\alpha }{dE_\alpha }dE_\alpha . \end{aligned}$$

(5)

The transition function $P_\mathrm{dec}(\nu _\alpha ^{(r)}\rightarrow \nu _\beta ^{(s)})$ is [34, 46]:

$$\begin{aligned} P_\mathrm{dec}\left( \nu _\alpha ^{(r)}\rightarrow \nu _\beta ^{(s)}\right)= & {} \left| \sum _{I=1}^3\left( {\tilde{U}}^{(r)}\right) _{I\alpha }^{-1}\right. \nonumber \\&\left. \times \exp \left[ -i\frac{{\tilde{m}}_I^2 L}{2E_\alpha }\right] \exp \left[ -\frac{\tilde{\alpha }_I L}{2E_\alpha }\right] {\tilde{U}}_{\beta I}^{(s)}\right| ^2 \nonumber \\&\times \,\delta _{rs}\,\delta (E_\alpha -E_\beta ) \nonumber \\&+ P_\mathrm{vis}(E_\alpha ,\,E_\beta ) \end{aligned}$$

(6)

The first term gives the ID contribution, which also includes neutrino oscillations. This term describes the oscillated neutrino flux, taking into account the loss of neutrinos due to the decay. The surviving neutrinos do not change their energy nor helicity, so $\delta (E_\alpha -E_\beta )$ and $\delta _{rs}$ terms need to be added. Notice that the mixing depends on the helicity: ${\tilde{U}}^{(-)}={\tilde{U}}$ and ${\tilde{U}}^{(+)}={\tilde{U}}^*$.

As mentioned in the Introduction, in this work we assume that the neutrinos from the original flux decay into lighter ones. The decay products, for all practical purposes, have the same direction as the initial neutrinos, and reach the far detector at the same time. Thus, they become a secondary contribution to the neutrino flux. This secondary contribution is described by the last term of Eq. (6), and is referred to as VD.

One important thing is that VD does not have a $\delta (E_\alpha -E_\beta )$ function, so the state after the decay can have a different energy compared to the one before the decay. As the effective Hamiltonian depends on the energy, this means that the eigenvalues and mixing (Eqs. (3) and (4)) before the decay shall be different to those after the decay. Thus, we denote the matter basis before the decay with a tilde, and the matter basis after the decay with a hat (for example, ${\tilde{m}}$ vs ${\hat{m}}$).

We refer to the addition of ID and VD contributions as full decay (FD). In contrast, the results without neutrino decay are referred to as standard oscillations (SO).

The matrices ${\tilde{U}}$ and ${\hat{U}}$ relate the interaction eigenstates $\alpha ,\,\beta$ with the matter eigenstates $I,\,J$, which is important to connect the neutrino production and detection with its propagation. In the same way, to connect the neutrino propagation with its decay, we need the matrices ${\tilde{C}}$ and ${\hat{C}}$, that relate the matter eigenstates $I,\,J$ with the mass eigenstates $i,\,j$. We define them in the following way:

$$\begin{aligned} {\tilde{C}}_{Ij}^{(r)}=\sum _{\rho =e,\mu ,\tau }{\tilde{U}}_{\rho I}^{(r)}(U_0)_{\rho j}^{(r)*}&{\hat{C}}_{Ij}^{(s)}=\sum _{\rho =e,\mu ,\tau }{\hat{U}}_{\rho I}^{(s)}(U_0)_{\rho j}^{(s)*} \end{aligned}$$

(7)

With this notation, assuming that the neutrino decays only once in its path, we can use the methods described in [28] to write:

$$\begin{aligned} P_\mathrm{vis}(E_\alpha ,\,E_\beta )= & {} \int d\ell \,\left| \sum _{I={\tilde{1}}}^{{\tilde{3}}}\left( {\tilde{U}}^{(r)}\right) ^{-1}_{I\alpha } \exp \left[ -i\frac{{\tilde{m}}_I^2\, \ell }{2E_\alpha }\right] \right. \nonumber \\&\left. \times \exp \left[ -\frac{\tilde{\alpha }_I\, \ell }{2E_\alpha }\right] \sum _{i=2}^3\sum _{j=1}^{i-1}{\tilde{C}}^{(r)}_{Ii}\right. \nonumber \\&\quad \times \,\sqrt{\frac{d}{dE_\beta }\Gamma _{\nu _i^r\rightarrow \nu _j^s}(E_\alpha )} \nonumber \\&\left. \times \sum _{J={\hat{1}}}^{{\hat{3}}}\left( {\hat{C}}^{(s)}\right) _{jJ}^{-1}\exp \left[ -i\frac{{\hat{m}}_J^2 (L-\ell )}{2E_\beta }\right] \right. \nonumber \\&\times \left. \exp \left[ -\frac{{\hat{\alpha }}_J (L-\ell )}{2E_\beta }\right] {\hat{U}}_{\beta J}^{(s)}\right| ^2 \end{aligned}$$

(8)

The above equation can be understood as follows: first, the source generates a neutrino interaction eigenstate, described by the index $\alpha$. The propagation occurs on the basis where $H(E_\alpha )$ is diagonal, so we need to use ${\tilde{U}}^{-1}$ to rotate into the matter eigenstates, with index I. These states propagate a distance $\ell$, and then decay. The decay, however, is defined on the mass basis in vacuum, so we switch to this basis (index i) using ${\tilde{C}}$. After the decay, one has a mass eigenstate $\nu _j^{(s)}$ that must be propagated a distance $(L-\ell )$. Nevertheless, this state is again not an eigenstate of $H(E_\beta )$, meaning we need to change basis again. For this, we use ${\hat{C}}^{-1}$ to rotate into matter eigenstates (index J). These states propagate the distance $(L-\ell )$, and then interact with the detector. One obtains the final flavour $\beta$ by switching back to the interaction basis using ${\hat{U}}$.

Assuming constant $N_e$, one can evaluate Eq. (8) using [28]. We obtain:

$$\begin{aligned}&P_\mathrm{vis}(E_\alpha ,\,E_\beta )\nonumber \\&\quad = 2\sum _{I={\tilde{1}}}^{{\tilde{3}}}\sum _{J={\hat{1}}}^{{\hat{3}}}\sum _{M={\tilde{1}}}^{{\tilde{3}}}\sum _{N={\hat{1}}}^{{\hat{3}}} \left( {\tilde{U}}^{(r)}\right) _{I \alpha }^{-1}\left( {\tilde{U}}^{(r)}\right) _{M\alpha }^{-1*}\,{\hat{U}}_{\beta J}^{(s)}\,{\hat{U}}_{\beta N}^{(s)*} \nonumber \\&\qquad \times E_\beta \frac{\left[ (E_\beta /E_\alpha )\tilde{\alpha }_{<IM>}-{\hat{\alpha }}_{<JN>}\right] -i\left[ (E_\beta /E_\alpha )\Delta {\tilde{m}}^2_{IM}-\Delta {\hat{m}}^2_{JN}\right] }{\left[ (E_\beta /E_\alpha )\tilde{\alpha }_{<IM>}-{\hat{\alpha }}_{<JN>}\right] ^2+\left[ (E_\beta /E_\alpha )\Delta {\tilde{m}}^2_{IM}-\Delta {\hat{m}}^2_{JN}\right] ^2} \nonumber \\&\qquad \times \left\{ \exp \left[ -i\frac{\Delta {\hat{m}}_{JN}^2 L}{2E_\beta }\right] \exp \left[ -\frac{{\hat{\alpha }}_{<JN>}L}{2E_\beta } \right] \right. \nonumber \\&\qquad \left. -\exp \left[ -i\frac{\Delta {\tilde{m}}_{IM}^2 L}{2E_\alpha }\right] \exp \left[ -\frac{\tilde{\alpha }_{<IM>}L}{2E_\alpha } \right] \right\} \nonumber \\&\qquad \times \sum _{i=2}^3\sum _{j=1}^{i-1}\sum _{m=2}^3\sum _{n=1}^{m-1} {\tilde{C}}_{Ii}^{(r)}{\tilde{C}}_{Mm}^{(r)*}\left( {\hat{C}}^{(s)}\right) _{jJ}^{-1}\left( {\hat{C}}^{(s)}\right) _{nN}^{-1*}\nonumber \\&\qquad \times \sqrt{\frac{d}{dE_\beta }\Gamma _{\nu _i^r\rightarrow \nu _j^s}(E_\alpha )\frac{d}{dE_\beta }\Gamma _{\nu _m^r\rightarrow \nu _n^s}(E_\alpha )} \end{aligned}$$

(9)

where we have generically denoted $\alpha _{<IJ>}=\alpha _I+\alpha _J$. The calculation of $d\Gamma (\nu _\mathrm{ini}^r\rightarrow \nu _\mathrm{fin}^s\,J)/dE_\beta$ on the mass basis has been done before [46].

In this work, we shall again consider only one decay channel $(\nu ^r_3\rightarrow \nu ^s_1 J)$ and only one non-vanishing coupling. In this limit, we find our final expression:

$$\begin{aligned}&P_\mathrm{vis}(E_\alpha ,\,E_\beta )\nonumber \\&\quad = 2\sum _{I=1}^3\sum _{J=1}^{3}\sum _{M=1}^3\sum _{N=1}^3 \left( {\tilde{U}}^{(r)}\right) _{I \alpha }^{-1}\left( {\tilde{U}}^{(r)}\right) _{M\alpha }^{-1*}{\hat{U}}_{\beta J}^{(s)}{\hat{U}}_{\beta N}^{(s)*} \nonumber \\&\qquad \times \frac{\left[ (E_\beta /E_\alpha )\tilde{\alpha }_{<IM>}-{\hat{\alpha }}_{<JN>}\right] -i\left[ (E_\beta /E_\alpha )\Delta {\tilde{m}}^2_{IM}-\Delta {\hat{m}}^2_{JN}\right] }{\left[ (E_\beta /E_\alpha )\tilde{\alpha }_{<IM>}-{\hat{\alpha }}_{<JN>}\right] ^2+\left[ (E_\beta /E_\alpha )\Delta {\tilde{m}}^2_{IM}-\Delta {\hat{m}}^2_{JN}\right] ^2} \nonumber \\&\qquad \times \left\{ \exp \left[ -i\frac{\Delta {\hat{m}}_{JN}^2 L}{2E_\beta }\right] \exp \left[ -\frac{{\hat{\alpha }}_{<JN>}L}{2E_\beta } \right] \right. \nonumber \\&\qquad \left. -\exp \left[ -i\frac{\Delta {\tilde{m}}_{IM}^2 L}{2E_\alpha }\right] \exp \left[ -\frac{\tilde{\alpha }_{<IM>}L}{2E_\alpha } \right] \right\} \nonumber \\&\qquad \times {\tilde{C}}_{I3}^{(r)}{\tilde{C}}_{M3}^{(r)*}\left( {\hat{C}}^{(s)}\right) _{1J}^{-1}\left( {\hat{C}}^{(s)}\right) _{1N}^{-1*}\nonumber \\&\qquad \times \left( \frac{(E_\beta /E_\alpha )\,\alpha _3}{E_\alpha }\right) \left( 1-\frac{m_3^2}{E_\alpha ^2}\right) ^{-1/2}\frac{x_{31}^2}{(x_{31}^2-1)}F^{rs}_g(E_\alpha ,\,E_\beta ) \nonumber \\&\qquad \times \Theta _H(E_\alpha -E_\beta )\,\Theta _H(x_{31}^2E_\beta -E_\alpha ) \end{aligned}$$

(10)

where $g=\{g_s,\,g_p\}$ indicates the non-vanishing coupling, $x_{if}=m_i/m_f>1$, $\Theta _H(x)$ is the Heaviside function, and:

$$\begin{aligned} F^{\pm \pm }_{g_s}(E_\alpha ,\,E_\beta )&= \frac{1}{E_\alpha \,E_\beta }\frac{(E_\alpha +x_{if}E_\beta )^2}{(x_{if}+1)^2}\nonumber \\ F^{\pm \mp }_{g_s}(E_\alpha ,\,E_\beta )&= \frac{(E_\alpha -E_\beta )}{E_\alpha \,E_\beta }\frac{(x_{if}^2E_\beta -E_\alpha )}{(x_{if}+1)^2} \end{aligned}$$

(11a)

$$\begin{aligned} F^{\pm \pm }_{g_p}(E_\alpha ,\,E_\beta )&= \frac{1}{E_\alpha \,E_\beta }\frac{(E_\alpha -x_{if}E_\beta )^2}{(x_{if}-1)^2}\nonumber \\ F^{\pm \mp }_{g_p}(E_\alpha ,\,E_\beta )&= \frac{(E_\alpha -E_\beta )}{E_\alpha \,E_\beta }\frac{(x_{if}^2E_\beta -E_\alpha )}{(x_{if}-1)^2}. \end{aligned}$$

(11b)

The reader should be aware that in the next-to-last line of Eq. (10), we have the plain $\alpha _3$ on the mass basis in vacuum. This represents the neutrino-Majoron couplings.

We point out that in [46] it was shown that the scalar and pseudoscalar couplings give undistinguishable effects when the $\nu _1$ mass, $m_\mathrm{lightest}$, is vanishing. In contrast, when $m_\mathrm{lightest}=0.07$ eV, its largest value allowed by cosmology [63], the choice of coupling leads to a different phenomenology. From now on, we shall write $x_{31}\rightarrow \infty$ when $m_\mathrm{lightest}$ is vanishing, and $x_{31}\rightarrow 1$ when $m_\mathrm{lightest}=0.07$ eV.

An interesting feature of Eq. (10) is that one could have oscillations not only before, but also after the decay, even in the case of only one decay channel. However, we have checked that, for current neutrino beam energies, this would occur for values of $N_e$ ten times larger than those found on Earth. Thus, we do not pursue this effect any further.

3 Impact on $(\Phi \times \sigma )$ at DUNE and ANDES

As a first approach in our analysis, we characterize the spectrum of the flux, weighted by its corresponding cross-section. We use:

$$\begin{aligned} (\Phi \times \sigma )_\beta \equiv \sum _{s}\sigma ^{s,\,\mathrm CC}_\beta (E_\beta )\frac{d\Phi ^{(s)}_\beta }{dE_\beta }, \end{aligned}$$

(12)

where $d\Phi ^{(s)}_\beta /dE_\beta$ was defined in Eq. (5) and $\sigma _\beta ^{s,\,\mathrm CC}$ is the charged-current cross-section for a neutrino with helicity s and flavour $\beta$. This has been evaluated according to the signal channels in the AEDL rules presented in Sect. 4 (Table 1), implying that this parameter encodes both contributions from helicity-conserving and helicity-changing VD.

We are interested in experiments where matter effects can be relevant in neutrino decay. We shall consider DUNE [57], and a future hypothetical experiment based on the Agua Negra Deep Experiment Site (ANDES) underground laboratory [58].

DUNE is one of the most promising experiments in neutrino physics, and is foreseen to begin operation within the next ten years. It comprises a 1300 km baseline, starting at the Long-Baseline Neutrino Facility (LBNF) at Fermilab, and extending to Sanford Underground Research Facility (SURF). This corresponds to an average matter density of $\rho _\mathrm{DUNE}=2.96$ g/cm$^3$. The beam will be generated using a 1.07 MW primary proton beam from the Main Injector, and is expected to be detected using a massive liquid argon time-projection chamber (LArTPC) detector, located deep underground at SURF. Among other goals, the DUNE experiment aims to measure the CP violating phase $\delta _\mathrm{CP}$, determine the neutrino mass ordering, resolve the octant for the atmospheric mixing angle, search for proton decay and detect and measure the $\nu _e$ flux from a core-collapse supernova within our galaxy [57].

The neutrino flux at the Far Detector used in our simulation was provided by the DUNE collaboration [64], in both neutrino mode (Forward Horn Current – FHC) and antineutrino mode (Reverse Horn Current – RHC). We assume $1.47\times 10^{21}$ protons-on-target (POT) per year [57, 64, 65], for 3.5 years on each mode. The 40 kt Far Detector will consist of four LArTPC modules [66], which provide a fine-grained image of neutrino scattering events. The latter information shall be relevant in our full simulation, to be presented on Sect. 4.

On the other hand, ANDES is a proposed underground laboratory in the Southern Hemisphere. It would be built in the deepest part ($\sim$1750 m) of the Agua Negra tunnel, linking Argentina and Chile below the Andes mountain range [58, 67]. Its construction, together with the tunnel, is planned to happen around 2018–2026. The topics of the ANDES scientific program include neutrinos, dark matter, geophysics, biology, among others [68].

With the aim of evaluating Eq. (12) in the context of ANDES, we assume a hypothetical neutrino beam starting at Fermilab. The corresponding neutrino flux is identical to the one we use for DUNE, but properly scaled to match the distance from source to detector. The ANDES baseline would be of order 7650 km [69], which corresponds¹ to an average matter density $\rho _\mathrm{ANDES}=4.7$ g/cm$^3$. Thus, we can expect matter effects in neutrino decay to be much more important than in DUNE.

Fig. 1

$(\Phi \times \sigma )$ at the DUNE baseline, for $\nu _\mu$ disappearance (top) and $\nu _e$ appearance (bottom) channels. The left (right) column shows results considering the FHC (RHC) flux. The black and red curves corresponds to ID and VD, respectively. Solid lines consider matter effects, while dashed lines show an equivalent scenario in vacuum

We fix the neutrino mass differences and mixing parameters at the following values [71]: $s^2_{12}=0.306$, $s^2_{23}=0.441$, $s^2_{13}=0.02166$, $\Delta m^2_{21}=7.5\times 10^{-5}$ eV$^2$ and $\Delta m^2_{31}=2.524\times 10^{-3}$ eV$^2$, for normal ordering. We also set $\delta _\mathrm{CP}=-\pi /2$ and $\alpha _3=4\times 10^{-5}$ eV$^2$, the latter corresponding roughly to a $10\%$ of $\langle E_\alpha \rangle /L$ for DUNE. For VD, we set $x_{31}\rightarrow \infty$, such that the type of coupling becomes irrelevant. For definiteness, from now on we choose only a non-vanishing scalar $(g_s)_{31}$ coupling.

In Fig. 1 we show the expected $(\Phi \times \sigma )$ for $\nu _\mu$ disappearance and $\nu _e$ appearance at DUNE. We separate the ID and VD components of the flux, in order to properly understand the difference between each contribution. For comparison, we show equivalent curves for the case where no matter effects are present. It is important to note that, even though not shown, for the value of $\alpha _3$ we are using we have confirmed that the ID curves are almost identical to those obtained with the SO hypothesis.

For $\nu _\mu$ disappearance (top row), we find that matter effects are almost inexistent for neither ID nor VD components. The ID contribution dominates both the FHC and RHC modes. Moreover, for FHC we have VD between one and two orders of magnitude lower than ID, while for RHC the difference between contributions is somewhat smaller. We can understand why this ID-VD difference in RHC is smaller than the one for FHC in terms of the helicity-changing decay channel. For FHC, the latter decay channel implies that part of the $\nu ^{(-)}$ flux becomes $\nu ^{(+)}$. In contrast, on RHC we have a fraction of $\nu ^{(+)}$ turning into $\nu ^{(-)}$. Since the $\nu ^{(-)}$ have a larger cross-section than the $\nu ^{(+)}$, the relative ID-VD difference in $(\Phi \times \sigma )$ becomes smaller for RHC. This observation is consistent with our results in [46] for T2K, and shall also be relevant in the $\nu _e$ appearance channel.

This result is strictly valid in the $x_{31}\rightarrow \infty$ limit. For larger masses ($x_{31}\rightarrow 1$), the behaviour depends on the coupling. For a scalar coupling, the helicity-flipping decay channel is suppressed, while for a pseudoscalar coupling it is slightly enhanced. This shall have a direct impact on the VD spectrum.

Regardless of these points, the fact remains that ID dominates the flux. Therefore, since for this value of $\alpha _3$ the ID component coincides with SO, we expect that $\nu _\mu$ disappearance shall not be strongly modified by neutrino decay. Then this channel will constrain $\Delta m^2_{32}$ and $\sin ^22\theta _{23}$ reliably, independently of the presence of decay.

For $\nu _e$ appearance, the low value of $\alpha _3$ again implies that the ID spectrum shall be very similar to the one of SO. However, we find that the VD contribution at low energy can dominate the flux, being up to one order of magnitude larger than ID for the RHC flux. Consequently, we can expect a better constraint on $\alpha _3$ by using $\nu _e$ appearance instead of $\nu _\mu$ disappearance data, which is again consistent with the results in [46]. We find that the VD contribution for $\nu _e$ appearance is slightly smaller for FHC than for RHC. As in $\nu _\mu$ disappearance, this can be understood in terms of helicity-flipping decay channels, and the different cross-sections.

Fig. 2

$(\Phi \times \sigma )$ at ANDES, for $\nu _\mu$ disappearance (top) and $\nu _e$ appearance (bottom) channels. The left (right) column shows results considering the FHC (RHC) flux. Curves as in Fig. 1

Another important feature is that matter effects are only relevant for the ID contribution to $\nu _e$ appearance. Given the similarity between ID and SO, these correspond to the typical matter effects for the SO scenario. Thus, a useful result is that, to a very good approximation, one can ignore matter effects in VD completely at this baseline.

In Fig. 2, we show both channels for the ANDES baseline. Here, since the baseline is much larger than in DUNE, matter effects are much more relevant. This time, we use $\alpha _3=8\times 10^{-6}$ eV$^2$, which again corresponds roughly to $10\%$ of $\langle E\rangle /L$ for ANDES.

A first observation is that, qualitatively, the patterns we observed at DUNE are repeated in ANDES. Nevertheless, quantitatively, matter effects make a difference. It is remarkable that even for ID in $\nu _\mu$ disappearance the discrepancy between vacuum and matter effects is noticeable.

The most striking consequence of the inclusion of matter effects is the suppression of the ID contribution to $\nu _e$ appearance at low energy, leaving a dominant VD component, with respect to the vacuum case. We want to emphasize that even though both DUNE and ANDES feature a large VD contribution in $\nu _e$ appearance, the reason for this in each case is different. In addition, VD decay is affected in more strongly than in DUNE, at most being reduced by $\sim 14\%$.

In this channel, the opposite effect happens at larger energy. In vacuum, the difference between ID and VD is smaller than the one in matter. This is due to an enhancement of the ID component in front of matter effects.

Fig. 3

The VD contribution to $(\Phi \times \sigma )$, for $\nu _e$ appearance, at DUNE (top) and ANDES (bottom). The left (right) column shows results considering the FHC (RHC) flux. Solid (dashed) lines include (do not include) matter effects. We show results for $x_{31}\rightarrow 1$ with scalar and pseudoscalar couplings in red and blue, respectively. The green curve considers $x_{31}\rightarrow \infty$, where the two types of coupling give the same effects

For larger neutrino masses ($x_{31}\rightarrow 1$), the size and shape of the VD contribution to $(\Phi \times \sigma )$ varies with the type of coupling, as shown in Fig. 3. For all plots we have the same qualitative behaviour. For FHC, the curve with $x_{31}\rightarrow 1$ and scalar coupling is always above the other curves. In contrast, for RHC, it is the curve with $x_{31}\rightarrow 1$ and pseudoscalar coupling the largest one, with the exception of very small energy, where the curve with $x_{31}\rightarrow \infty$ can have higher values. This can be explained in terms of helicity-changing decays, as was discussed earlier.

An important difference between the $x_{31}\rightarrow 1$ and $x_{31}\rightarrow \infty$ scenarios is that the spectrum of the latter favours lower energies. This makes this situation more susceptible to the low energy thresholds of an experiment.

4 Sensitivity and parameter fits at DUNE

The impact of neutrino decay to the sensitivity to $\alpha _3$ and to parameter fits, at DUNE, has been studied earlier in [47, 48]. Nevertheless, there is still room for further development in these analyses. To begin with, the analysis in [48] only considered ID, so our results with FD are complementary. With respect to [47], which also takes FD into account, we also include the $\nu _\mu$ disappearance channel into our analysis, and display the sensitivity as a function of oscillation parameters. In addition, our approach strictly follows the procedure outlined in [28], while the formulae used in [47] seem to use additional assumptions which we do not have (in particular, the arguments used to build Eqs. (9) and (13) of that work).

4.1 Event generation

In order to calculate the number of events within an energy bin, we follow the same procedure as in [46]. The number of events of flavour $\beta$ in the energy bin i, with helicity s and going through interaction $int=\{CC,\,NC\}$, is obtained from:

$$\begin{aligned} N^{(s),\mathrm int}_{i,\beta }= \int dE_\beta \,K^\mathrm{int}_i(E_\beta )\,\sigma ^{s,\mathrm int}_\beta (E_\beta )\frac{d\Phi ^{(s)}_\beta }{dE_\beta }, \end{aligned}$$

(13)

where $\sigma _\beta ^{s,\mathrm int}(E_\beta )$ is the cross section for process int, and:

$$\begin{aligned} K^\mathrm{int}_i(E_\beta )=\int _{E_\mathrm{i, min}}^{E_\mathrm{i, max}}dE_\mathrm{bin}\,\epsilon ^\mathrm{int}_\beta (E_\mathrm{bin})\, R^\mathrm{int}(E_\mathrm{bin}-E_\beta ). \end{aligned}$$

(14)

The detector efficiency $\epsilon ^\mathrm{int}_\beta (E_\mathrm{bin})$ and resolution function $R^\mathrm{int}(E_\mathrm{bin}-E_\beta )$ are taken from the DUNE Collaboration public files [64]. Both Eqs. (13) and (14) are used to calculate signal and background events, with $\epsilon ^\mathrm{int}_\beta (E_\mathrm{bin})$ and $R^\mathrm{int}(E_\mathrm{bin}-E_\beta )$ properly adjusted, according to the information in [64].

For ID, this calculation is carried out within GLoBES [72, 73]. However, as the neutrino decay effective Hamiltonian is not Hermitian, we could not use the built-in diagonalization algorithms. Thus, for this case, we modified the source probability library, importing our own probability matrix into GLoBES.

Table 1 AEDL rules for $\nu _e^{(-)}$, $\nu _e^{(+)}$ appearance, and $\nu _\mu ^{(-)}$, $\nu _\mu ^{(+)}$ disappearance. We denote $\nu _\alpha =\nu _\alpha ^{(-)}$ and ${\bar{\nu }}_\alpha =\nu _\alpha ^{(+)}$

In the case of VD, we generated the fluxes in Eq. (5) externally, and used these in Eq. (13), in GLoBES, to calculate the event rates. Finally, we modified the channels of each GLoBES rule, such that helicity change was taken into account. For example, given the channel $\nu _\alpha ^{(-)} \rightarrow \nu _\beta ^{(-)}$, we add $\nu _\alpha ^{(-)} \rightarrow \nu _\beta ^{(+)}$.

The full set of rules for each channel, for signal and background, including both ID and VD, is shown in Table 1. The Table separates contributions to signal and background based on the final state. However, they are added in the $\chi ^2$ analyses presented below. Notice that for ID, we include the contributions coming from both of the original $\nu _\mu ^{(-)}$ and $\nu _\mu ^{(+)}$ components of the flux, for both $\nu _\mu$ disappearance and $\nu _e$ appearance channels. This is also done for VD in the $\nu _e$ appearance channel. However, we find it suffices to include only the FHC $\nu _\mu ^{(-)}$ (RHC $\nu _\mu ^{(+)}$) components in VD for $\nu _\mu$ disappearance, with the other component being negligible.

4.2 Analysis and results

We begin by studying the DUNE sensitivity to $\alpha _3$ in the FD scenario. As $\theta _{13}$ is fixed by reactor $\nu _e$ disappearance measurements [74], which are not strongly affected by FD [46], we focus on the effect of varying $\theta _{23}$ and $\delta _\mathrm{CP}$. As in the previous Section, we present results for $x_{31}\rightarrow \infty$, fixing the other oscillation parameters at their best fit values. Data from both $\nu _\mu$ disappearance and $\nu _e$ appearance are taken into account.

To compute $\chi ^2$, we calculate the number of events $N_i$ for each combination of $\theta _{23}$, $\delta _\mathrm{CP}$ and $\alpha _3$, for the energy bin i. $N_i$ is calculated using Eq. (13), adding over helicities and interactions, for signal and backgrounds. This is compared to the events generated by a fixed set of parameters, refered to as true values. We define $\chi ^2$ using:

$$\begin{aligned}&\chi ^2(\theta _{23},\,\delta _\mathrm{CP},\,\alpha _3,\,\theta ^\mathrm{true}_{23},\,\delta ^\mathrm{true}_\mathrm{CP},\,\alpha ^\mathrm{true}_3) \nonumber \\&\quad = \sum _i^\mathrm{bins} \frac{\left( N_i\left( \theta _{23},\delta _\mathrm{CP},\alpha _3 \right) - N_i\left( \theta ^\mathrm{true}_{23},\delta ^\mathrm{true}_\mathrm{CP},\alpha ^\mathrm{true}_3\right) \right) ^2}{N_i\left( \theta ^\mathrm{true}_{23},\delta ^\mathrm{true}_\mathrm{CP},\alpha ^\mathrm{true}_3\right) } \end{aligned}$$

(15)

Fig. 4

Sensitivity to $\alpha _3$ as a function of $\theta ^\mathrm{true}_{23}$. Left: the sensitivity when combining $\nu _\mu$ disappearance and $\nu _e$ appearance, after both FHC and RHC runs. The curves show the sensitivity after marginalizing over $\delta ^\mathrm{true}_\mathrm{CP}$, with the shaded area below giving the sensitivity for fixed values of the phase. Right: sensitivity curves for only FHC or RHC runs, including only the $\nu _e$ appearance channel

Fig. 5

As in Fig. 4, but as a function of $\delta _\mathrm{CP}^\mathrm{true}$

To analyze the sensitivity to $\alpha _3$, we set $\theta _{23}=\theta _{23}^\mathrm{true}$ and $\delta _\mathrm{CP}=\delta _\mathrm{CP}^\mathrm{true}$. The sensitivity to $\alpha _3$ as a function of $\theta _{23}^\mathrm{true}$ is obtained by marginalizing over $\delta _\mathrm{CP}^\mathrm{true}$, and setting $\alpha _3^\mathrm{true}=0$ eV$^2$, that is:

$$\begin{aligned} \left. \chi ^2(\theta ^\mathrm{true}_{23},\,\delta ^\mathrm{true}_\mathrm{CP},\,\alpha _3,\,\theta ^\mathrm{true}_{23},\,\delta ^\mathrm{true}_\mathrm{CP},\,0)\right| _{\min \delta ^\mathrm{true}_\mathrm{CP}} \end{aligned}$$

(16)

Similarly, the sensitivity to $\alpha _3$ as a function of $\delta _\mathrm{CP}^\mathrm{true}$ is obtained with:

$$\begin{aligned} \left. \chi ^2(\theta ^\mathrm{true}_{23},\,\delta ^\mathrm{true}_\mathrm{CP},\,\alpha _3,\,\theta ^\mathrm{true}_{23},\,\delta ^\mathrm{true}_\mathrm{CP},\,0)\right| _{\min \theta ^\mathrm{true}_{23}} \end{aligned}$$

(17)

In Fig. 4, we show the sensitivity to $\alpha _3$, for different values of a given $\theta _{23}^\mathrm{true}$. On the left panel, we use the full potential of DUNE by combining both $\nu _\mu$ disappearance and $\nu _e$ appearance, for both FHC and RHC modes. The sensitivity improves for larger values of $\sin ^2\theta _{23}^\mathrm{true}$, since the VD contribution is proportional to this parameter. This means that for a fixed $\chi ^2$, a larger $\theta ^\mathrm{true}_{23}$ requires a smaller $\alpha _3$. We find that DUNE is sensitive at $3\sigma$ ($5\sigma$) to values of $\alpha _3$ around $(4-7)\times 10^{-6}$ eV$^2$ ($(0.7-1.1)\times 10^{-5}$ eV$^2$). Moreover, the shaded areas below the curves indicate the sensitivity for fixed values of the phase, that is, without marginalizing. We consider the full range of $\delta _\mathrm{CP}^\mathrm{true}$, and see that the uncertainty in this parameter has little impact on the sensitivity to $\alpha _3$.

The right panel of Fig. 4 shows the sensitivity using only the $\nu _e$ appearance channel, after marginalizing, for different modes. An analogous plot for the $\nu _\mu$ disappearance channel would have a much worse sensitivity to $\alpha _3$, which follows from our arguments in Sect. 3. As expected from our earlier discussion, both FHC and RHC curves have the same downward slope.

The $3\sigma$ ($5\sigma$) sensitivity for FHC is around $(1.0-1.6)\times 10^{-5}$ eV$^2$ ($(1.6-2.7)\times 10^{-5}$ eV$^2$). For RHC instead, the sensitivity is around $(5-9)\times 10^{-6}$ eV$^2$ ($(0.8-1.4)\times 10^{-5}$ eV$^2$). Thus, we see that sensitivity is determined principally by the RHC mode. This is easy to understand using the study of $(\Phi \times \sigma )$ in Sect. 3: the RHC flux is dominated by $\nu ^{(+)}$, which have a smaller cross-section than $\nu ^{(-)}$, such that the helicity-changing decays in VD have a greater impact.

The sensitivity to $\alpha _3$ as a function of $\delta ^\mathrm{true}_\mathrm{CP}$ is shown in Fig. 5. On the left panel, we have the same curves and shaded areas as in Fig. 4, but this time marginalizing over $\theta _{23}^\mathrm{true}$. We find the curves to be almost flat, suggesting that the value of $\delta ^\mathrm{true}_\mathrm{CP}$ is not important for the determination of $\alpha _3$.

Nevertheless, the right panel of Fig. 5 shows a very different picture. Here, we find that the sensitivity has a very strong dependence on $\delta ^\mathrm{true}_\mathrm{CP}$, particularly on the FHC mode. This, of course, is due to the influence of $\delta ^\mathrm{true}_\mathrm{CP}$ on the total number of events. For example, for the FHC mode, a positive $\delta ^\mathrm{true}_\mathrm{CP}$ would reduce the events coming from ID, which has a stronger dependence on $\delta ^\mathrm{true}_\mathrm{CP}$ than VD. As a consequence of this reduction, a relatively small $\alpha _3$ is sufficient to generate a VD contribution that can be comparable to the ID component. This makes this point more sensitive to low values of $\alpha _3$. In contrast, a negative $\delta ^\mathrm{true}_\mathrm{CP}$ implies a larger number of events expected from ID, such that a larger $\alpha _3$ is required to reach the same level of sensitivity compared to positive $\delta _\mathrm{CP}^\mathrm{true}$. Thus, the ratio between the ID and VD components has a relevant impact on the sensitivity.

The difference between the largest and smallest values of $\alpha _3$ in a sensitivity curve is most pronounced in the $5\sigma$ case of FHC, where the ID contribution clearly dominates (see Fig. 1), and is modulated by the value of $\delta _\mathrm{CP}^\mathrm{true}$. For the $3\sigma$ case the behaviour is the same, but the difference in $\alpha _3$ is diminished, due to the lesser number of ID events.

Due to CP conjugation, the situation for the RHC flux is opposite: scenarios with positive $\delta ^\mathrm{true}_\mathrm{CP}$ are less sensitive to $\alpha _3$, in comparison to those with negative values. Here, the difference between the largest and smallest $\alpha _3$ is very small, that is, it has a milder dependence on $\delta _\mathrm{CP}^\mathrm{true}$, as the overall number of events is also small and the VD contribution is usually comparable to the ID one. When combining both FHC and RHC, we find that the latter has a stronger pull on the sensitivity. In addition, the curve averages out, leaving an almost flat result.

Fig. 6

Sensitivity to $\alpha _3$ in DUNE, combining $\nu _\mu$ disappearance and $\nu _e$ appearance, after both FHC and RHC runs. We marginalize over $\delta _\mathrm{CP}^\mathrm{true}$ and $\theta _{23}^\mathrm{true}$. The horizontal lines indicate the $3\sigma$ and $5\sigma$ confidence levels

In Fig. 6 we show the sensitivity to $\alpha _3$ for both $x_{31}\rightarrow 1$ and $x_{31}\rightarrow \infty$ scenarios. As we have emphasized earlier, only in the former case it is necessary to distinguish between scalar and pseudoscalar couplings. In the Figure, the pseudoscalar coupling has a much better sensitivity, reaching values of $\alpha _3=3.8\times 10^{-6}$  eV$^2$ ($6.4\times 10^{-6}$ eV$^2$) at $3\sigma \,(5\sigma )$, compared to the scalar coupling, which reaches $\alpha _3=5.2\times 10^{-6}$ eV$^2$ ($8.8\times 10^{-6}$ eV$^2$). This is due to the increased helicity-changing transitions that are typical of this coupling, increasing tensions with RHC expectations. In fact, for the pseudoscalar coupling, we find that the RHC-only curve clearly dominates the overall sensitivity. In contrast, the scalar coupling has less helicity-changing transitions, so the sensitivity does not improve as much. Nevertheless, the $x_{31}\rightarrow 1$ scenario with scalar coupling has a better sensitivity than the $x_{31}\rightarrow \infty$ case, which reaches $\alpha _3=6.1\times 10^{-6}$ eV$^2$ ($1.0\times 10^{-5}$ eV$^2$). The reason for this is that, as we can see in Fig. 3, the VD peak of the $x_{31}\rightarrow \infty$ scenario appears at very low values of energy, which are cut off by the experimental thresholds included in our simulation.

These numbers can be compared to other results in the literature, at $90\%$ C.L. For instance, for FD at T2K and MINOS [46], the best limit is of $\alpha _3<5.6\times 10^{-5}$ eV$^2$. The authors of [48] work in the context of ID at DUNE, and obtain a sensitivity around $1.5\times 10^{-5}$ eV$^2$. In contrast, the authors of [47] take FD, and their best sensitivity is as low as $3.4\times 10^{-6}$ eV$^2$. At this confidence level, our best sensitivity is of $2.0\times 10^{-6}$ eV$^2$, which is comparable to the limit obtained with atmospheric neutrinos with ID, of $\alpha _3<2.2\times 10^{-6}$ eV$^2$ [35].

Fig. 7

Allowed regions, taking SO as a theoretical hypothesis, assuming that the measured signal was generated by SO (gray region), ID (blue curve) or FD (red curve). We include $\nu _\mu$ disappearance and $\nu _e$ appearance channels. The left (right) column takes $\theta _{23}^\mathrm{true}=42.8^\circ$ ($47.2^\circ$), while the upper (lower) row has $\delta ^\mathrm{true}_\mathrm{CP}=90^{\circ }$ ($-90^\circ$). We take $\alpha ^\mathrm{true}_3=4\times 10^{-5}$ eV$^2$. Dashed and solid lines (dark and light regions) correspond to $3\sigma$ and $5\sigma$ confidence levels, respectively, with the dots indicating the best-fit points

For the final part of our analysis, we compare the impact of the SO, ID and FD scenarios in the determination of $\theta _{23}$ and $\delta _\mathrm{CP}$. For each scenario we generate a specific number of events by fixing oscillation and decay parameters in combinations of true values: $\theta _{23}^\mathrm{true}=\{42.8^\circ ,\,47.2^\circ \}$, $\delta ^\mathrm{true}_\mathrm{CP}=\{90^{\circ },\,-90^\circ \}$. For the ID and FD scenarios, we set $\alpha ^\mathrm{true}_3=4\times 10^{-5}$ eV$^2$, as was done on Sect. 3. When performing the fit, we assume SO as a theoretical hypothesis, that is, we evaluate:

$$\begin{aligned} \chi ^2\left( \theta _{23},\,\delta _\mathrm{CP},\,0,\,\theta ^\mathrm{true}_{23},\,\delta ^\mathrm{true}_\mathrm{CP},\,\alpha _3^\mathrm{true}\right) \end{aligned}$$

(18)

where, of course, $\alpha _3^\mathrm{true}=0$ eV$^2$ in the SO scenario. We include both $\nu _\mu$ disappearance and $\nu _e$ appearance, with both FHC and RHC runs.

Results are shown in Fig. 7. We find small differences between SO and ID, which is due to $\alpha _3$ not being large enough to have any effect. The reason for including these curves is mainly for comparison with FD, which is the main focus of this work. If we had used a larger $\alpha _3$, we would have found a stronger impact on the determination of $\theta _{23}$, as was reported in detail in [48].

On the other hand, the FD scenario is much more susceptible to values of $\alpha _3$ of this order of magnitude. This is consistent with Fig. 1, that is, that the additional VD component can dominate the flux at low energy, which in turn is the main reason for the increased sensitivity shown in Fig. 6.

The regions can be understood by comparing the impact of FD on $\nu _\mu$ disappearance and $\nu _e$ appearance. First, since we know from Sect. 3 that $\nu _\mu$ disappearance is not strongly affected by FD, we find the same kind of constraints on the minimum and maximum possible values of $\theta _{23}$, including the octant degeneracy [75,76,77].

The important modification on the fit happen because of the $\nu _e$ appearance channel. As reported in Sect. 3, the inclusion of VD leads to an increase in events in both FHC and RHC modes within our simulated data sample. This fact is the reason why the SO fit in general prefers larger values of $\theta _{23}$, since the leading term in the $\nu _\mu \rightarrow \nu _e$ oscillation probability is proportional to $\sin ^2\theta _{23}$. In fact, true points generated with FD in the lower octant of $\theta _{23}$ can have an SO solution on the higher octant. For true points generated in the higher octant, the $\nu _\mu$ disappearance constraint does not allow $\theta _{23}$ to increase.

Our true points were generated for values of $\delta ^\mathrm{true}_\mathrm{CP}$ where the CP asymmetry is maximal. However, in the fit the SO regions tend to prefer values of $\delta _\mathrm{CP}$ closer to $\pm \pi$, such that the asymmetry is diminished.

5 Conclusions

Starting from the treatment given in [28], we have formulated a description for matter effects in visible neutrino decay. In order to understand these effects, we have implemented a $(\Phi \times \sigma )$ study in two scenarios. On the first one, we use the flux and baseline for DUNE, while on the second one we use the same flux, but consider the corresponding baseline for ANDES. For DUNE, we find that only the ID component of $\nu _e$ appearance can be affected by matter, the effect for all other components, such as VD at $\nu _e$ appearance, can be completely ignored. In contrast, for ANDES, not only do we have an enhanced effect on the ID component due to matter, but also find that the VD component receives a relevant modification. This is especially important for $\nu _e$ appearance, but the effects can also be seen on $\nu _\mu$ disappearance.

Another important part of this work was devoted to the calculation of the sensitivity of DUNE to $\alpha _3$. To this end, we performed as very detailed simulation of DUNE, described in Sect. 4, using the publicly available files of the Collaboration [64]. On that section, we showed the dependence of the sensitivity on $\theta _{23}$ and $\delta _\mathrm{CP}$, using both $\nu _\mu$ disappearance and $\nu _e$ appearance channels, and both FHC and RHC modes. Our final sensitivity to $\alpha _3$ depended on the lightest neutrino mass (encoded on the value of $x_{31}$) and on the type of coupling between the neutrino and the Majoron (scalar or pseudoscalar). For the $x_{31}\rightarrow 1$ scenario, we found the sensitivity of DUNE to be:

$$\begin{aligned} \alpha ^{(s)}_3< 2.8\times 10^{-6}{~\mathrm eV}^2,&\alpha ^{(p)}_3 < 2.0\times 10^{-6}{~\mathrm eV}^2~ \end{aligned}$$

(19)

while for $x_{31}\rightarrow \infty$:

$$\begin{aligned} \alpha _3^{(s,\,p)}<3.2\times 10^{-6}{~\mathrm eV}^2 \end{aligned}$$

(20)

We note that these values of sensitivity are the best ones obtained so far for long-baseline experiments, and are comparable to the current limits set by atmospheric neutrinos using ID [35].

Finally, in order to understand the impact of a non-zero $\alpha _3$ on the determination of oscillation parameters, we performed a fit on $\theta _{23}$ and $\delta _\mathrm{CP}$ assuming SO, with data generated for the FD scenario. We found that the allowed regions would shift towards larger values of $\theta _{23}$, and towards CP-conserving values of $\delta _\mathrm{CP}$.