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

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 νμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _\mu $$\end{document} disappearance and νe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _e$$\end{document} appearance channels, for both FHC and RHC modes. The sensitivity to the decay constant α3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _3$$\end{document} can be as low as 2×10-6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\times 10^{-6}$$\end{document} eV2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document} 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 θ23\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{23}$$\end{document} and δCP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _\mathrm{CP}$$\end{document}, and find that the best-fit value of θ23\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{23}$$\end{document} can move from a true value at the lower octant towards the higher octant.

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,. In these kind of models, neutrinos and Majorons can have two types of couplings, scalar (g s ) or pseudoscalar (g p ): a e-mail: jones.j@pucp.edu.pe 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 experi-ment to the decay constant α 3 . On that Section we also study the impact of neutrino decay on oscillation fits.

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 (ν e , ν μ , ν τ ) 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 ) αi , where α = e, μ, τ, 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 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 α i = E i , we have on the interaction basis: 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Ũ α I matrices (I =1,2,3): where H diag has complex eigenvalues: 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 α 3 different from zero, we can obtain non-vanishingα 1 ,α 2 andα 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 (r ) α /d E α , where α refers to the flavour, r = (+, −) is the helicity, and E α is the energy, which are determined when the neutrino is produced. In the presence of oscillations and decay, the neutrino flux of flavour β, helicity s and energy E β , can be calculated with: The transition function P dec (ν (r ) α → ν (s) β ) is [34,46]: 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 δ(E α − E β ) and δ rs terms need to be added. Notice that the mixing depends on the helicity:Ũ (−) =Ũ andŨ (+) =Ũ * . 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 δ(E α − E β ) 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,m vsm).
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Ũ andÛ relate the interaction eigenstates α, β 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 matricesC andĈ, that relate the matter eigenstates I, J with the mass eigenstates i, j. We define them in the following way: With this notation, assuming that the neutrino decays only once in its path, we can use the methods described in [28] to write: The above equation can be understood as follows: first, the source generates a neutrino interaction eigenstate, described by the index α. The propagation occurs on the basis where H (E α ) is diagonal, so we need to useŨ −1 to rotate into the matter eigenstates, with index I . These states propagate a distance , and then decay. The decay, however, is defined on the mass basis in vacuum, so we switch to this basis (index i) usingC. After the decay, one has a mass eigenstate ν (s) j that must be propagated a distance (L − ). Nevertheless, this state is again not an eigenstate of H (E β ), meaning we need to change basis again. For this, we useĈ −1 to rotate into matter eigenstates (index J ). These states propagate the distance (L − ), and then interact with the detector. One obtains the final flavour β by switching back to the interaction basis usingÛ .
Assuming constant N e , one can evaluate Eq. (8) using [28]. We obtain: where we have generically denoted α <I J> = α I + α J . The calculation of d (ν r ini → ν s fin J )/d E β on the mass basis has been done before [46].
In this work, we shall again consider only one decay channel (ν r 3 → ν s 1 J ) and only one non-vanishing coupling. In this limit, we find our final expression: where g = {g s , g p } indicates the non-vanishing coupling, is the Heaviside function, and: The reader should be aware that in the next-to-last line of Eq. (10), we have the plain α 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 ν 1 mass, m lightest , is vanishing. In contrast, when m 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 → ∞ when m lightest is vanishing, and x 31 → 1 when m 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.

Impact on ( × σ ) 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: where is the charged-current cross-section for a neutrino with helicity s and flavour β. 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 ρ 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 δ CP , determine the neutrino mass ordering, resolve the octant for the atmospheric mixing angle, search for proton decay and detect and measure the ν 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×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 (∼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 1 to an average matter density ρ ANDES = 4.7 g/cm 3 . Thus, we can expect matter effects in neutrino decay to be much more important than in DUNE.
In Fig. 1 we show the expected ( × σ ) for ν μ disappearance and ν 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 α 3 we are using we have confirmed that the ID curves are almost identical to those obtained with the SO hypothesis.
For ν μ 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 ν (−) flux becomes ν (+) . In contrast, on RHC we have a fraction of ν (+) turning into ν (−) . Since the ν (−) have a larger cross-section than the ν (+) , the relative ID-VD difference in ( × σ ) becomes smaller for RHC. This observation is consistent with our results in [46] for T2K, and shall also be relevant in the ν e appearance channel.
This result is strictly valid in the x 31 → ∞ limit. For larger masses (x 31 → 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 α 3 the ID component coincides with SO, we expect that ν μ disappearance shall not be strongly modified by neutrino decay. Then this channel will constrain m 2 32 and sin 2 2θ 23 reliably, independently of the presence of decay.
For ν e appearance, the low value of α 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 α 3 by using ν e appearance instead of ν μ disappearance data, which is again consistent with the results in [46]. We find that the VD contribution for ν e appearance is slightly smaller for FHC than for RHC. As in ν μ disappearance, this can be understood in terms of helicity-flipping decay channels, and the different cross-sections.
Another important feature is that matter effects are only relevant for the ID contribution to ν 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 α 3 = 8 × 10 −6 eV 2 , which again corresponds roughly to 10% of E /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 ν μ 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 ν 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 ν 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 ∼ 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.
For larger neutrino masses (x 31 → 1), the size and shape of the VD contribution to ( × σ ) 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 → 1 and scalar coupling is always above the other curves. In contrast, for RHC, it is the curve with x 31 → 1 and pseudoscalar coupling the largest one, with the exception of very small energy, where the curve with x 31 → ∞ 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 → 1 and x 31 → ∞ 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.

Sensitivity and parameter fits at DUNE
The impact of neutrino decay to the sensitivity to α 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 con-sidered ID, so our results with FD are complementary. With respect to [47], which also takes FD into account, we also include the ν μ 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).

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 β in the energy bin i, with helicity s and going through interaction int = {CC, N C}, is obtained from: where σ s,int β (E β ) is the cross section for process int, and: The detector efficiency int β (E bin ) and resolution function R int (E bin − E β ) are taken from the DUNE Collaboration public files [64]. Both Eqs. (13) and (14) are used to calculate signal and background events, with int β (E bin ) and R int (E bin − E β ) 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.
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 ν (−) 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 χ 2 analyses presented below. Notice that for ID, we include the contribu-tions coming from both of the original ν (−) μ and ν (+) μ components of the flux, for both ν μ disappearance and ν e appearance channels. This is also done for VD in the ν e appearance channel. However, we find it suffices to include only the FHC ν (−) μ (RHC ν (+) μ ) components in VD for ν μ disappearance, with the other component being negligible.

Analysis and results
We begin by studying the DUNE sensitivity to α 3 in the FD scenario. As θ 13 is fixed by reactor ν e disappearance measurements [74], which are not strongly affected by FD [46], we focus on the effect of varying θ 23 and δ CP . As in the previous Section, we present results for x 31 → ∞, fixing the other oscillation parameters at their best fit values. Data from both ν μ disappearance and ν e appearance are taken into account.
To compute χ 2 , we calculate the number of events N i for each combination of θ 23 , δ CP and α 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 χ 2 using: To analyze the sensitivity to α 3 , we set θ 23 = θ true 23 and δ CP = δ true CP . The sensitivity to α 3 as a function of θ true 23 is obtained by marginalizing over δ true CP , and setting α true 3 = 0 eV 2 , that is: Similarly, the sensitivity to α 3 as a function of δ true CP is obtained with: In Fig. 4, we show the sensitivity to α 3 , for different values of a given θ true 23 . On the left panel, we use the full potential of DUNE by combining both ν μ disappearance and ν e appearance, for both FHC and RHC modes. The sensitivity improves for larger values of sin 2 θ true 23 , since the VD contribution is proportional to this parameter. This means that for a fixed χ 2 , a larger θ true 23 requires a smaller α 3 . We find that DUNE is sensitive at 3σ (5σ ) to values of α 3 around (4 − 7) × 10 −6 eV 2 ((0.7 − 1.1) × 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 δ true CP , and see that the uncertainty in this parameter has little impact on the sensitivity to α 3 .
The right panel of Fig. 4 shows the sensitivity using only the ν e appearance channel, after marginalizing, for different modes. An analogous plot for the ν μ disappearance channel would have a much worse sensitivity to α 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 sensitivity to α 3 as a function of δ true 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 θ true 23 . We find the curves to be almost flat, suggesting that the value of δ true CP is not important for the determination of α 3 .   Fig. 4, but as a function of δ true CP Nevertheless, the right panel of Fig. 5 shows a very different picture. Here, we find that the sensitivity has a very strong dependence on δ true CP , particularly on the FHC mode. This, of course, is due to the influence of δ true CP on the total number of events. For example, for the FHC mode, a positive δ true CP would reduce the events coming from ID, which has a stronger dependence on δ true CP than VD. As a con-sequence of this reduction, a relatively small α 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 α 3 . In contrast, a negative δ true CP implies a larger number of events expected from ID, such that a larger α 3 is required to reach the same level of sensitivity compared to positive δ true CP . Thus, the ratio between The difference between the largest and smallest values of α 3 in a sensitivity curve is most pronounced in the 5σ case of FHC, where the ID contribution clearly dominates (see Fig. 1), and is modulated by the value of δ true CP . For the 3σ case the behaviour is the same, but the difference in α 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 δ true CP are less sensitive to α 3 , in comparison to those with negative values. Here, the difference between the largest and smallest α 3 is very small, that is, it has a milder dependence on δ true CP , 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.
In Fig. 6 we show the sensitivity to α 3 for both x 31 → 1 and x 31 → ∞ 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 α 3 = 3.8 × 10 −6 eV 2 (6.4 × 10 −6 eV 2 ) at 3σ (5σ ), compared to the scalar coupling, which reaches α 3 = 5.2 × 10 −6 eV 2 (8.8 × 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 → 1 scenario with scalar coupling has a better sensitivity than the x 31 → ∞ case, which reaches α 3 = 6.1 × 10 −6 eV 2 (1.0 × 10 −5 eV 2 ). The reason for this is that, as we can see in Fig. 3, the VD peak of the x 31 → ∞ 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 α 3 < 5.6 × 10 −5 eV 2 . The authors of [48] work in the context of ID at DUNE, and obtain a sensitivity around 1.5×10 −5 eV 2 . In contrast, the authors of [47] take FD, and their best sensitivity is as low as 3.4 × 10 −6 eV 2 . At this confidence level, our best sensitivity is of 2.0×10 −6 eV 2 , which is comparable to the limit obtained with atmospheric neutrinos with ID, of α 3 < 2.2 × 10 −6 eV 2 [35].
For the final part of our analysis, we compare the impact of the SO, ID and FD scenarios in the determination of θ 23 and δ CP . For each scenario we generate a specific number of events by fixing oscillation and decay parameters in combinations of true values: θ true For the ID and FD scenarios, we set α true 3 = 4 × 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: where, of course, α true 3 = 0 eV 2 in the SO scenario. We include both ν μ disappearance and ν 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 α 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 α 3 , we would have found a stronger impact on the determination of θ 23 , as was reported in detail in [48].
On the other hand, the FD scenario is much more susceptible to values of α 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 ν μ disappearance and ν e appearance. First, since we know from Sect. 3 that ν μ disappearance is not strongly affected by FD, we find the same kind of constraints on the minimum and maximum possible values of θ 23 , including the octant degeneracy [75][76][77].
The important modification on the fit happen because of the ν e appearance channel. As reported in Sect. 3, the inclu- 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 ν μ disappearance and ν e appearance channels. The left (right) column takes θ true 23 = 42.8 • (47.2 • ), while the upper (lower) row has δ true CP = 90 • (−90 • ). We take α true 3 = 4 × 10 −5 eV 2 . Dashed and solid lines (dark and light regions) correspond to 3σ and 5σ confidence levels, respectively, with the dots indicating the best-fit points sion 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 θ 23 , since the leading term in the ν μ → ν e oscillation probability is proportional to sin 2 θ 23 . In fact, true points generated with FD in the lower octant of θ 23 can have an SO solution on the higher octant. For true points generated in the higher octant, the ν μ disappearance constraint does not allow θ 23 to increase.
Our true points were generated for values of δ true CP where the CP asymmetry is maximal. However, in the fit the SO regions tend to prefer values of δ CP closer to ±π , such that the asymmetry is diminished.

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 ( × σ ) 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 ν e appearance can be affected by matter, the effect for all other components, such as VD at ν 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 ν e appearance, but the effects can also be seen on ν μ disappearance.
Another important part of this work was devoted to the calculation of the sensitivity of DUNE to α 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 θ 23 and δ CP , using both ν μ disappearance and ν e appearance channels, and both FHC and RHC modes. Our final sensitivity to α 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 → 1 scenario, we found the sensitivity of DUNE to be: α (s) 3 < 2.8 × 10 −6 eV 2 , α ( p) 3 < 2.0 × 10 −6 eV 2 (19) while for x 31 → ∞: α (s, p) 3 < 3.2 × 10 −6 eV 2 (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 α 3 on the determination of oscillation parameters, we performed a fit on θ 23 and δ CP assuming SO, with data generated for the FD scenario. We found that the allowed regions would shift towards larger values of θ 23 , and towards CP-conserving values of δ CP .