Visible neutrino decay in the light of appearance and disappearance long-baseline experiments

We investigate the present constraints from MINOS and T2K experiments for the neutrino decay scenario induced by non-diagonal couplings of Majorons to neutrinos. As novelty, on top of the typical invisible decay prescription, we add the contribution of visible decay, where final products can be observed. This new effect depends on the nature of the neutrino-Majoron coupling, which can be of scalar or pseudoscalar type. Using the combination of disappearance data from MINOS and disappearance and appearance data from T2K, for normal ordering, we constrain the decay parameter α ≡ E Γ for the heaviest neutrino, where E and Γ are the neutrino energy and width, respectively. We find that when considering visible decay within appearance data, one can improve current neutrino long-baseline constraints up to α<O10−5eV2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \alpha <\mathcal{O}\left(1{0}^{-5}\right)\ {\mathrm{eV}}^2 $$\end{document}, at 90% C.L., for both kinds of couplings, which is better by one order of magnitude compared to previous bounds.

In the last years, the interplay between oscillations and fast neutrino decay has been studied in a model-independent manner [52,69]. This has been achieved by considering that the decay products of the neutrino decay are invisible to the detector. Such a situation can happen, for example, if the decay product consists of lighter, sterile states [40,[46][47][48][49][50][51][53][54][55][56][57][58][59][60][61][62]. Throughout this work, we refer to this scenario as invisible decay. In this situation, one can bound the neutrino lifetime-to-mass ratio τ /m or, alternatively, the parameter α = E Γ, related to the neutrino width Γ evaluated at the energy E. In terms of these parameters, the two most important constraints to date are the following: • For solar neutrinos, the studies in [60,61] have constrained the lifetime of ν 2 , giving τ 2 /m 2 ≥ 7.2 × 10 −4 s · eV −1 at 99% C.L., equivalent to α 2 < 9.1 × 10 −13 eV 2 .
In contrast, not much work has been done in the direction of visible decay, where the final decay products involve active neutrinos, which can be detected. In oscillation experiments, this manifests itself as an additional contribution to invisible decay. We refer to the addition of both visible and invisible contributions as full decay. Nevertheless, in order to describe the visible contribution, one needs an expression for the differential width, so it is not possible to carry out this analysis in a completely model-independent way.
Models for fast neutrino decay usually involve the coupling of two neutrinos with a massless scalar, called a Majoron [71][72][73][74]. This interaction can proceed through scalar (g s ) or pseudoscalar (g p ) couplings. As we shall see in this paper, both couplings can be related to the α decay constant by g 2 s,p ∼ 16πα/m 2 , where m is the neutrino mass. It is then of our interest to understand how such an interaction affects the results of oscillation experiments, for both possible kinds of coupling.
To this end, we implement a full decay formalism within two experiments: MINOS and T2K. As mentioned before, the ν µ → ν µ disappearance channels of both experiments were used in [59] to study invisible decay. Of course, one of our objectives is to update the results in these channels within the full decay framework. However, this time we also have data from the T2K ν µ → ν e appearance channel, for both neutrinos and antineutrinos. As far as we are aware of, this is the first time that visible decay is considered as modification of a neutrino appearance phenomenon, and we find it has important consequences.
Before we proceed, we need to point out that important constraints on neutrino -Majoron couplings exist [75][76][77][78][79][80]. If one assumes that these can be directly translated as bounds on the decay constant α, and assuming for illustration a neutrino mass m = 0.01 eV, they would imply constraints as low as α < O 10 −25 eV 2 . However, this translation is not always straightforward, and most of the bounds rely on additional assumptions, apart from the Majoron hypothesis. Thus, constraints coming from neutrino decay within oscillation experiments play an important role in transparently understanding the kind of phenomenology these new couplings bring, as they are established under controlled experimental conditions. Therefore, in the following, in addition to updating the bounds in [59], we give a detailed explanation of how the neutrino spectrum of long-baseline experiments is modified, hoping that these insights might be useful in other scenarios.
We present our framework in Section 2, and provide insight on how the neutrino transition operator will be modified by the different kind of coupling. In Section 3, we describe our experimental setup for T2K and MINOS, as well as describe our procedure when carrying out the fit. In Section 4 we present our results for each experiment, and then combine the data. Since the combination does not provide a better fit in comparison to the one from standard oscillations, we extract bounds on our decay parameter α. We conclude in Section 5.
We also include three Appendices. In Appendix A, we present the general, model independent, full decay framework, within the oscillation scenario. In Appendix B we give details of the neutrino -Majoron coupling, and explain how we obtain the formulae used in this work. Finally, in Appendix C we describe how we carry out the MINOS and T2K simulations.

Neutrino Visible Decay in Oscillation Experiments
We consider a neutrino oscillation experiment, where a neutrino flux is directed towards a detector located a large distance away. During their propagation, the neutrinos are subject to an evolution function which takes into account their oscillation and decay, such that the flux arriving at the detector is modified.
The first term in the evolution function includes the standard oscillation contributions, multiplied by an exponential governed by the neutrino decay parameter α i . For stable neutrino mass eigenstates, all α i are zero, and we recover the standard oscillation formula. In the following, we shall refer to this term as the contribution from invisible decay (ID).
If we included this first term only, we would be describing the decay of the mass eigenstates into lighter non-interacting states. However, we contemplate the possibility that these lighter states consist of active neutrinos with lower energy. This is taken into account by the second term, which we refer to as the contribution from visible decay (VD). For the rest of this work, we refer to the addition of the two terms as full decay (FD).
The general formula for the visible decay function G rs αβ (E α , E β ) can be found on Appendix A. In this work we use a simplified version, considering only one decay channel (ν r i → ν s f J, where ν r i and ν s f are on the mass basis). It has the following form: We can understand Eq. (2.3) in the following way. The term (1 − e −Γ i L ) |U αi | 2 , when multiplied by the original flux φ ν r α , gives the amount of ν r i mass eigenstates that decay. When the latter is multiplied by the normalized spectrum, 1 , we obtain the number of ν s f mass eigenstates, with energy between E β and E β + dE β , produced from the decay. The amount of final ν s β flavour eigenstates is then obtained by changing back to the flavour basis with |U βf | 2 .
The full expression for dΓ(ν r i → ν s f J)/dE β is derived from the interaction Lagrangian, which is written in terms of scalar (g s ) or pseudoscalar (g p ) neutrino -Majoron couplings. Details are given in Appendix B. Furthermore, if we take only one non-vanishing coupling (either g s or g p ), the final expression is greatly simplified: with the visible decay function given by Here we have x if = m i /m f > 1, the label g = {g s , g p } indicates the non-vanishing coupling, Θ H (x) is the Heaviside function and we have replaced the total width Γ i → α i /E α . The  cases: For a given x ij and energy, we have that ±± and ±∓ transitions are complementary, which is expected, since F ±± g + F ±∓ g = 1. Studying the behaviour of these functions is important in order to understand how scalar and pseudoscalar couplings affect our analysis of T2K and MINOS. To this end, we show in Figure 1 the functions F rs g , defined in Eq. (2.5), for fixed E α = 3 GeV, and for several values of x if , as a function of E β . On the top (bottom) rows we show results for scalar (pseudoscalar) couplings. The left column shows chirality-conserving transitions (±±), while the right columns shows chirality-changing transitions (±∓).
For large x if (purple curve), we see that both scalar and pseudoscalar couplings have the same behaviour (compare purple curve within the same column). For E β close to E α , we find that ±± processes (left column) are favoured, while for lower energies, ±∓ processes (right column) dominate.
For lower x if (dashed blue curve) we have a different behaviour depending on the coupling. For scalar couplings, we see that ±± transitions dominate all final energies. Meanwhile, pseudoscalar couplings have a mixed behaviour, with a clear preference for ±∓ transitions.
From now on, we assume a normal ordering scenario, with ν 3 unstable and decaying exclusively into ν 1 . We label α 3 → α, m 1 → m light and (g s,p ) 31 → g s,p (see Appendix B). Notice that cosmological bounds [81] require i m i < 0.23 eV. For normal ordering, and taking squared mass differences at their best fit points, for example, this would rule out x 31 < 1.24. This bound is taken into account in our final result.

T2K
The T2K experiment [82] has been running in the latest years, sending a neutrino beam to the Super-Kamiokande (SK) detector [9], located 295 km away from the source. The detection process uses Cerenkov radiation to identify neutrinos, however, this is blind to the charge of the associated lepton. Thus, it is unable to distinguish between neutrinos and antineutrinos, as a consequence, in this work we shall sum the neutrino and antineutrino contributions when calculating the number of events. The experiment has two running modes, called "neutrino" and "antineutrino" runs.
The neutrino run of the T2K experiment consisted in delivering a primarily ν (−) µ beam to the SK detector, with a final luminosity of 7.48 × 10 20 protons on target (POT). Their results in the ν disappearance channel imposed strong bounds on the sin 2 2θ 23 and ∆m 2 32 parameter space [15,82]. In addition, through the observation of ν appearance channel, it provided the first indication of non-zero θ 13 [83].
In the current antineutrino run, a primarily ν (+) µ beam illuminates the SK detector. The luminosity released to the public corresponds to 7.47×10 20 POT, and the combination of both runs has given hints favouring a negative value of δ CP [84,85]. In our analysis we use both datasets.
On the following, we shall focus on the T2K neutrino run. For the antineutrino run, the same analysis applies, taking the CP conjugates of the states.
Given that ν is the largest component in the beam [82], we use the FD mode (visible and invisible contributions), for both appearance and disappearance channels. For the other components of the flux (ν For the ν µ disappearance channel, we include FD for ν and ID for ν µ . As usual, the signal consists of charged-current quasielastic (CCQE) interactions, while the main background includes the charged-current non-quasielastic (CCnQE) and neutral current (NC) interactions.
The ν e appearance channel has many contributions. As mentioned previously, we consider the FD for ν e . In this case, we also include the FD for ν e , but take only the ID contribution for ν e . Both CCQE and CCnQE interactions are considered part of the signal. The background consists of the ν (±) e → ν (±) e processes, as well as NC contributions.
We also include full decay when considering the NC contribution to the events. Notice that the latter shall now depend on the neutrino mixing angles, due to the non-unitarity of the decay mechanism [40].
We give more details on the generation of events at T2K in Appendix C.1.

MINOS
MINOS was a long-baseline neutrino experiment [16] which used two detectors and was exposed to a neutrino beam produced at Fermilab (NuMI beam line [86]). The latter is a two-horn-focused neutrino beam that can be configured in two ways: Forward Horn Current (FHC), to produce a beam optimized for muon neutrinos, and Reverse Horn Current (RHC), for a beam optimized for muon anti-neutrinos. The Near and the Far Detectors are located at 1 km and 735 km from the target, respectively. The data set used in our analysis comes from FHC mode with 10.71 × 10 20 POT. The beam composition was 92.9% of ν e . We use the data set that comprises the charged current (CC) contained-vertex neutrino and anti-neutrino disappearance data [16,87].
Notice that, in contrast to T2K, the MINOS magnetized muon spectrometer does distinguish between neutrinos and antineutrinos. Thus, even though we only use the FHC mode, we do include the ν contribution. Analogously, for ν (+) µ disappearance data, we take into account the FD for ν The details about the MINOS reconstruction data can be found in Appendix C.2.

Statistical Analysis for T2K and MINOS
The relevant parameters in this study, which we shall vary, are s 2 23 , s 2 13 , δ CP , ∆m 2 32 , α and m light . We keep fixed s 2 12 = 0.306 and ∆m 2 21 = 7.5 × 10 −5 eV 2 . To perform the fit for T2K, we use a χ 2 function similar to the one used in [82]. In order to take into account systematic errors, we include normalisation and energy calibration nuisance parameters, n x and t x , respectively, for signal (x = s) and background (x = b). For a given set of oscillation and decay parameters, and for each channel, the χ 2 is minimized with respect to n x and t x , adding appropriate pull factors [88,89]: Here, N β, fit is the sum of expected signal and background events per bin, which also involve the nuisance parameters. Denoting the numerical prediction per bin for signal and background ν where E min (E max ) is the minimum (maximum) energy in the analyzed neutrino spectrum, E bin is the average bin energy, andÊ = 1 2 (E max + E min ) is the average energy of the spectrum.
The other parameters appearing in Eq. (3.1) are N β, obs , which is the observed number of ν (−) β and ν (+) β events [84,85], and σ nx and σ tx are the respective uncertainties in normalisation and tilt, set both equal to 10%.
In the case of MINOS, using a similar notation, we use the following χ 2 function: The parameter σ bin represents the statistical and systematic uncertainty extracted from data set, while σ ns , σ nb are equal to 14.7% and 4%, respectively [14]. Finally, we can also include information from reactor data. For example, the Daya Bay experiment [90] has given the most precise bound on s 2 13 , which is sin 2 2θ 13 = 0.092 ± 0.017 [11]. We have estimated the impact of neutrino decay on this experiment, and find that the ratio of the background-subtracted spectra to prediction assuming no oscillation is modified within the given error bars [91]. We then assume that no reactor disappearance experiment is affected by neutrino decay. Thus, only for our final result constraining the decay parameter α, we include an additional pull factor on our χ 2 function: with s 2 13,reactor = 0.0243 and σ 13,reactor = 0.0026, following the same procedure as in [82]. In our fit, we vary m light and α logarithmically, from 5 × 10 −3 to 10 −1 eV and 2 × 10 −6 to 5 × 10 −4 eV 2 , respectively. The oscillation angles θ 13 and θ 23 are also scanned, with 0 ≤ s 2 13 ≤ 0.05 and 0.32 ≤ s 2 23 ≤ 0.86. The phase δ CP is scanned between −π and π, while the squared mass difference ∆m 2 32 is varied between 1.6 × 10 −3 and 2.8 × 10 −3 eV 2 .

Analysis in T2K
In the following, we explore how the inclusion of the neutrino FD distorts the currently allowed regions for oscillation parameters, obtained from standard oscillations (SO).
In Figure 2, we show the spectrum of neutrino appearance and disappearance events for the neutrino run, for SO, ID and FD. We have fixed s 2 23 = 0.532, s 2 13 = 0.022 and ∆m 2 32 = 2.545 × 10 −3 eV 2 , which correspond to the best-fit parameters for T2K data. We . Spectrum of disappearance (top row) and appearance (bottom row) events, for the neutrino run. The spectrum for SO is shown in red, dotted (dashed) for δ CP = −π/2 (+π/2). The spectrum for ID is shown dashed, in black, and the spectrum for FD is shown solid, in black. For both we set δ CP = +π/2. We show results for scalar (pseudoscalar) couplings on the left (right). T2K data is shown in gray [84].
have set α = 4 × 10 −5 eV 2 , which is roughly 10% of the mean lifetime that T2K should be sensitive to (≈ E/L). We have also set m light = 7 × 10 −2 eV to maximize the difference between the scalar and pseudoscalar scenarios (see Figure 1).
The top panels show neutrino disappearance events. We see that for the chosen value of α, there is no big difference between the SO and FD scenarios. Nevertheless, even though not shown, one finds that for larger values of α the decay scenario erases the oscillation dip. This is due to neutrino decay behaving as a decoherence term in the oscillation formulae.
The bottom panels show neutrino appearance events. Here, the SO scenario is shown for δ CP = ±π/2. Notice that one obtains a larger amount of events when δ CP = −π/2. In contrast, we show ID and FD scenarios only with δ CP = +π/2. A first important feature is that, although the ID contribution is similar to the SO result (for the same value of δ CP ), when comparing the FD and ID spectra it is evident that the VD contribution is sizeable. Thus, this decay scenario has a stronger impact on appearance than on disappearance measurements. A second important feature is that the FD scenario, with δ CP = π/2, can have a similar number of events as the SO result with δ CP = −π/2 at the cost of bringing some distortion to the low-energy (high-energy) part of the spectrum, for scalar (pseudoscalar) couplings.   Figure 2, but for the antineutrino run, with T2K data from [85].
We show in Figure 3 the corresponding spectra for the antineutrino run. For disappearance events (top row), the conclusions are similar to those for Figure 2, that is, for the given value of α, neutrino decay does not significantly modify the spectrum.
On the other hand, as in Figure 2, for the appearance spectrum we show SO with δ CP = ±π/2. On this case, δ CP = −π/2 leads to a smaller number of events compared to δ CP = +π/2. We again show the spectrum for ID and FD with δ CP = +π/2, and again find the VD contribution to be sizeable. However, contrary to the neutrino run, this time the number of events greatly exceeds the SO prediction for δ CP = −π/2. In fact, the ID contribution alone is already too large. This means that the antineutrino run shall be relevant when disentangling the value of δ CP within the decay scenario.
On both Figures we find that the FD results differ for scalar and pseudoscalar couplings. The reason for this is that we have chosen a large m light . In the previous section (Figure 1), we found that for scalar couplings, large values of m light (x 31 → 1) favour ν (±) → ν (±) over ν (±) → ν (∓) transitions. The opposite behaviour is seen in pseudoscalar couplings, which allow a much larger proportion of chirality-changing transitions. This, of course, shall modify results, as ν (−) and ν (+) states have different cross-sections.
Given the argument above, we see that on the neutrino run the peak of the spectrum is smaller for the pseudoscalar than for the scalar coupling. As said previously, pseudoscalar couplings give a larger rate of ν (+) , and since these have a smaller cross-section, the number of events is lower. Furthermore, for the antineutrino run, the pseudoscalar coupling gives a larger rate of ν (−) , which have higher cross-sections, and thus lead to more events.
In Figure 4, we show allowed regions in several subspaces of the parameter space, under the hypothesis of the FD scenario. In all plots, red (blue) curves refer to scalar (pseudoscalar) couplings. Moreover, we show the corresponding regions on the SO scenario in shades of grey.
If we concentrate on the s 2 23 − ∆m 2 32 subspace (top right panel), dominated by disappearance measurements, we notice that the resulting allowed regions for scalar and pseudoscalar couplings are equal, and very similar to those for the SO scenario. This confirms that the FD scenario has little impact on neutrino disappearance events.
The FD scenario has a much larger impact in the s 2 13 − δ CP region (top left panel), where appearance measurements are crucial. We find that the upper bound to s 2 13 is similar to the one for SO, for both scalar and pseudoscalar couplings. However, the lowest possible value of s 2 13 decreases. Moreover, we find that a large region of positive δ CP is ruled out, in agreement with the T2K fit for the SO scenario.
One can get a better insight on the situation in appearance measurements from the other panels, where we show the allowed values of α as a correlation with s 2 13 and m light . On the bottom left panel, we see that for low values of α, the allowed range for s 2 13 is well bounded (consistent with SO). Then, as α increases, the FD formalism starts to influence, and the range for s 2 13 is shifted to smaller values. This means that the lower number of neutrinos from oscillations, due to the smaller s 2 13 , is compensated by the extra neutrinos coming from decay. Finally, in both cases, too large values of α generate too many appearance events, which is incompatible with T2K data, restricting the allowed α to values of O 10 −5 eV 2 . This automatically leaves unaffected the s 2 23 − ∆m 2 32 subspace, which requires much larger values of α to modify significantly the spectrum.
On the other hand, on the bottom right panel, one finds the confidence level contours for α as a function of m light . For scalar (pseudoscalar) couplings, we find that for large values of m light , the contours reach higher (lower) values of α. For pseudoscalar couplings, large mass values favour chirality-changing transitions, which gives tensions when applied to both neutrino and antineutrino runs, as seen in Figures 2-3. Thus, in this case, smaller values of α are compatible with data. Opposed to the latter case, for scalar couplings with large m light , chirality-changing transitions are minimized, so visible decay has a smaller impact, and larger values of α are allowed. In contrast, for smaller m light , both scalar and pseudoscalar curves tend to the same value of α. This means that in this limit both couplings are undistinguishable, as expected from Figure 1 (for example, for x 31 = 100).
Notice that the best-fit points correspond to non-vanishing values of α. As we shall see later, these points are not statistically significant. Thus, from the T2K-only χ 2 analysis, we conclude that the FD solution does not improve the quality of the fit in a statistically significant way. Therefore, from T2K data we will later obtain constraints on α (see Figure 8).

Analysis in MINOS
In Figure 5, we show the spectrum of muon neutrino (top panels) and antineutrino (bottom panels) disappearance events of FHC mode, for SO, ID and FD scenarios. The first (second) column shows the scalar (pseudoscalar) couplings case. The panels follows the convention presented in Section 4.1, where we have fixed the oscillation parameters s 2 23 = 0.41, s 2 13 = 0.0243, ∆m 2 32 = 2.41 × 10 −3 eV 2 , which correspond to the best fit of MINOS to standard oscillation model, and δ CP = 0, because it is not sensitive in MINOS. When the decay is present, we have set α = 3.0 × 10 −4 eV 2 , and m light = 0.05 eV to investigate the difference between the scalar and pseudoscalar scenarios in MINOS (see Figure 1).
It is interesting to compare the ID with the FD case to see the effect of adding visible neutrino decay. As explained in the previous Section, and in Figure 1, when we fixed a large value of m light , this means that for scalar couplings we have a preference for ±± transitions. In contrast, pseudoscalar couplings have a clear preference for ±∓ transitions.  Figure 5. The spectrum of muon and anti-muon neutrinos events is presented for the SO, ID and FD scenarios. We fixed the oscillation parameters in ∆m 2 32 = 2.41 × 10 −3 eV 2 , s 2 23 = 0.41, s 2 13 = 0.0243, and δ CP = 0. We also fixed the decay parameter α = 3.0×10 −4 eV 2 and m light = 0.05 eV. The MINOS data points were taken from Ref. [16,87].
This explain why the visible neutrino decay contribution for muon neutrino spectrum is more significant in the scalar case. This happens because the NuMI beam configuration was composed of more than 90% of muon neutrinos favoring the transition ν µ . By the other side, looking the muon anti-neutrino spectrum, we see that the contribution of visible neutrino decay is higher for pseudoscalar case. Using the same previous reason, this occur because the transition ν is favored for the value of m light chosen. The best fit obtained for ID and FD scenario were for an α different from zero. Then, in MINOS, this effect results in an impact on oscillation parameters compared to the one obtained for SO. This same effect is observed at Ref. [59]. However, if we take the minimum χ 2 for FD and SO, which are 42.04 and 45.62 respectively, and taking into account the data bins and different number of parameters in each scenario, one can demonstrate through the Akaike Information Criterion [92,93] that there is no statistical difference between both  Figure 6. The allowed regions for the cases scalar and pseudoscalar using the disappearance signals from MINOS at FHC mode. The projections shown are to respect the oscillation parameters and the decay parameter. models.
In Figure 6, we show the 1σ and 90% C.L. allowed regions for many combination of the relevant parameters, already computing ν µ . The red (blue) curves refer to scalar (pseudoscalar) couplings. To compare the effect of neutrino decay, we add as well the limit when α → 0, which corresponds to SO scenario (full gray) in Eq. (2.2).
The allowed region in the (s 2 23 −∆m 2 32 ) plane show that the FD scenario has a significant impact when compared with SO scenario. This means that for MINOS, the contribution of events coming from visible decay is not enough to constrain the asymmetry at the probability of ID scenario. Another important observation is that we do not find significant differences between scalar and pseudoscalar couplings except by a region in (s 2 23 − ∆m 2 32 ) plane. We can see it for small values of ∆m 2 32 and higher for s 2 23 . We also present the two-dimensional projections (α − ∆m 2 32 ) and (α − s 2 23 ) in Figure 6. We can observe that for large values of α, the contours increase the range in ∆m 2 32 , and either allow higher values of s 2 23 . The right bottom panel in Figure 6 show the correlation between α and m light . We can conclude by it, that MINOS is not able to distinguish the α constraint for scalar and pseudoscalar couplings.

Combination of T2K and MINOS
In order to combine T2K and MINOS data, for every points in the parameter space we have added the respective χ 2 values for each experiment. We show the result in Figure 7. The panels follow the same conventions as in Figures 4 and 6.
On the top right panel we show the allowed s 2 23 − ∆m 2 32 subspace. We find that the resulting allowed area is not as large as the one in the MINOS-only analysis, in fact, for both scalar and pseudoscalar scenarios, the contours match those for the SO result. This is due to the value of α being strongly constrained by ν e appearance, down to values where MINOS has no sensitivity. This is analogous to the T2K-only case, where we saw that no effects were visible for ν µ disappearance. Therefore, the T2K analysis dominates the fit.
On the top left panel, we show the s 2 13 − δ CP region. The allowed area in the pseudoscalar scenario is very similar to the T2K-only result. For the scalar scenario, however, when compared to the T2K-only analysis, we find that the regions allow even smaller values of s 2 13 . In addition, in this scenario no value of δ CP is ruled out at 90% C.L. The qualitative behaviour observed on the bottom left panel, showing the s 2 13 −α plane, is very similar to the one of Figure 4. An evident difference is that, for both couplings, the best-fit for α is increased to a value of O 10 −5 . This implies that the MINOS data shall still relevant when bounding this parameter. Nevertheless, we shall find these best-fit values not to be statistically relevant.
Finally, on the bottom right panel, we can understand why all values of δ CP are allowed in the scalar coupling scenario. In the SO scenario, if we fix s 2 13 at the reactor best-fit value, positive δ CP implies too few (too many) events in the T2K neutrino (antineutrino) run. If we only consider data from the T2K neutrino run, we would find in both FD scenarios that positive δ CP is incompatible with a vanishing value of α. This means that the lesser number of events in the neutrino run, coming from the oscillation contribution, is compensated by the additional contribution from visible decay. This was previously explained in Section 4.1 (see the lower panels of Figure 2). Nevertheless, the data from the T2K antineutrino run is incompatible with this solution (lower panels of Figure 3), such that positive δ CP is not allowed. This is precisely what happened in Figure 4. However, when adding the MINOS data, we find that the latter pulls the fit towards larger α, for all values of δ CP for both couplings, reintroducing the ambiguity in δ CP for the scalar case. This ambiguity with δ CP does not happen for the pseudoscalar case, as we have seen, the incompatibility of a value of α of O 10 −5 , for δ CP = +π/2, in front of antineutrino data is much more serious. One would expect that, if further T2K antineutrino data exhibits the same preference for SO, the positive δ CP solution would eventually also be excluded for scalar coupling.
In Figure 8 we present show ∆χ 2 as a function of α, for the two scenarios (scalar and pseudoscalar couplings) and three fits (T2K, MINOS and T2K+MINOS) considered. We also add additional curves showing the impact of including reactor data, as explained in Section 3.3.
The addition of reactor data adds a strong pull factor towards non-zero s 2 13 . However, as the solutions for both couplings involve values of s 2 13 compatible with reactor data, the addition of the latter does not modify significantly the results.

Conclusions
In this work we have studied the implications of visible neutrino decay on T2K and MINOS experiments. We have considered that the decay proceeds through the coupling of two neutrinos and a Majoron, which can be either of scalar or pseudoscalar nature. We present our results for both couplings, separately, and describe their effects in detail.
In our analysis we consider neutrino disappearance and appearance channels, both in T2K and MINOS. The results for MINOS point towards a non-zero α as a best-fit solution. Nevertheless, the inclusion of T2K neutrino and antineutrino data excludes this result, being consistent with a vanishing α. Consequently, we put bounds on this parameter.
The power of the T2K data lies on the strong influence that visible neutrino decay has on appearance channels. The additional events from the full decay framework are significant even for values of α much smaller than the characteristic E/L where the neutrino experiment has its largest sensitivity. In contrast, for the neutrino disappearance spectrum, the inclusion of visible and invisible neutrino decay components does not differ from what is obtained when considering only the invisible component.
Our results depend on the type of non-zero coupling. For scalar coupling, we have briefly seen that adding reactor data improves the bounds. If we denote the allowed values of α, given the coupling c, due to the experiment exp, by α (c) exp , the best constraint we find at 90% C.L. is: Notice that these bounds depend crucially on the T2K antineutrino data, which is currently in its first years of data taking. We expect that the bounds shall improve with further data.

A Full Formula for Visible Neutrino Decay
The G rs αβ (E α , E β ) function describing visible neutrino decay can be derived from the original formula of [52], using the methods of [69]. Assuming the final state neutrinos do not oscillate, the result is: αm . Notice that, following Eq. (1.1), the PMNS matrix has been defined without the Majorana phases. In Eq. (A.1), the latter are taken into account within the parameter ω sr ijmn = γ s jj − γ r ii − γ s nn + γ r mm . We also denote: If we only have one unstable neutrino, this expression can be considerably simplified. We chose ν r 3 to be unstable, which means we fix i = m = 3 in Eq. (A.1), so we obtain: Furthermore, if only one decay channel is allowed (for example, ν r 3 → ν s 1 J, obtained by fixing j = n = 1 in Eq. (A.1)), we have: One must note that the previous formulae have been written in the convention m 1 < m 2 < m 3 , which coincides with the normal hierarchy. To obtain results for the inverted hierarchy, one either must write m 1 → m 3 , m 2 → m 1 and m 3 → m 2 , or modify the full PMNS matrix such that: U · diag 1, e i γ 22 2 , e i γ 33

B Formulae for Visible Neutrino Decay Rates
In the following, we shall concentrate on the differential decay width dΓ(ν r i → ν s f J)/dE β , for the general case where several decay channels exist. Our procedure follows [69], where we have a neutrino -Majoron coupling: where the neutrinos are on the mass eigenstates basis. With this, the partial decay rates of a neutrino with energy E α are [69,94]: , and: For our description of neutrino decay, we need the differential decay rate: The Θ(E α , E β ) function fixes the angle between the momenta of the incoming and outgoing neutrinos: which effectively imposes the following bound on E f : On the relativistic limit, this reduces to: Thus, in terms of the Heaviside function Θ H , we have: The neutrino decay amplitudes M(ν f J) can be computed as a function of (g p ) 2 if and (g s ) 2 if , the neutrino masses and the neutrino energies [69,94]. These can be recasted using the partial width of neutrinos, given by Eqs. (B.2). We first define: We can invert the Eq. (B.9) and solve for one coupling, such that either the scalar or the pseudoscalar couplings is written: Using Eqs. (B.10), we can re-write the expression of the amplitude M(ν f J) [69,94] in terms of α ij and one coupling. For example, in terms of (g if ), we have: Alternatively, in terms of (g p ) if , we find: For both cases, we have defined: These results are very useful, as they allow us to simplify the differential decay width. If we go to the limit where one coupling is zero, the other coupling is automatically determined by the value of α if . If, in addition, there exists only one allowed decay channel, then the whole process is governed by α i .
These sort of considerations lead us to the principal formula used in this work, Eq. (2.4).

C Formulae for T2K and MINOS
In order to obtain the number of events within an energy bin, one needs to multiply Eq. (2.1) by the ν s β cross-section and perform the E β integral considering the particular bin width, efficiency and energy resolution function. We calculate the number of events in the energy bin i, with chirality s and going through interaction int, using: where σ s,int β (E β ) denotes the appropriate cross section, and the "detection kernel" K int i (E β ) for bin i is defined: Here, int β (E bin ) denotes the detector efficiency, while the resolution function R(E int β,rec − E β , ς int β ) reflects our capacity to reconstruct the true neutrino energy E β . In both experiments we model this using a gaussian function: (C. 3) The variable E int β,rec is related to the bin energy E bin by a possible energy shift, and ς int β denotes the energy resolution width.

C.1 T2K
For T2K, the neutrino fluxes appearing in Eq. (2.1) are obtained from [82], and the cross sections are taken from GENIE [95]. For ν µ disappearance, we take the shift and resolution width ς int µ from [96]: with equal values for neutrino and antineutrino channels.
For ν e appearance, we use the "migration matrix" from [97], which is equivalent to using: For NC events, we fit the shape of the spectra shown in [84,98]. We set:  [82] as well as relevant plots in [84,85,98].

C.2 MINOS
The φ r να (E α ) described in Eq. (2.1) is the flux at Far Detector for MINOS. It was obtained by the product of the flux at the Near Detector (ND) and the Beam Matrix f F/N (E α , E β ), following: where E α is the bin energy and E β are the adjacent bins of E α . The ND flux φ (ND)r να (E α ) was taken from references [14,99]. To construct the Beam Matrix f F/N (E α , E β ), we used the matrix M (E α , E β ) extracted from reference [100], and also two Gaussian functions G l (G r ) with resolutions widths having energy dependence in the maximum until the second order. So, the f F/N (E α , E β ) can be described as: The functions G l (G r ) is responsible for extrapolate the flux of higher (lower) to lower (higher) energies.
In MINOS, the only considered interaction is CC. The cross sections σ s,CC β are taken from references [101,102]. We have taken E CC µ,rec = E bin , and assume the resolution widths ς CC β (E β ) to be linear in energy. The efficiency CC β of the detector was extracted from reference [14].