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

TL;DR

This paper explores constraints on neutrino decay involving Majorons using MINOS and T2K data, introducing visible decay effects that enhance the sensitivity of long-baseline experiments to decay parameters.

Contribution

It presents the first analysis including visible decay contributions in neutrino decay constraints from long-baseline experiments, improving bounds on decay parameters by an order of magnitude.

Findings

Visible decay improves constraints on neutrino decay parameters.

Decay parameter α constrained to less than 10^{-5} eV^2 at 90% C.L.

Results apply to both scalar and pseudoscalar Majoron couplings.

Abstract

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 $\alpha\equiv E\,\Gamma$ for the heaviest neutrino, where $E$ and $\Gamma$ 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 $\alpha< \mathcal{O}(10^{-5})$ eV$^2$, at 90\% C.L., for both kinds of couplings, which is better by one order of magnitude compared to previous bounds.

Equations

$$U_{\rm PMNS}=\left(\begin{array}[]{ccc}c_{12}c_{13}&s_{12}c_{13}&s_{13}e^{-i\delta_{\rm CP}}\\ -s_{12}c_{23}-c_{12}s_{23}s_{13}e^{i\delta_{\rm CP}}&c_{12}c_{23}-s_{12}s_{23}s_{13}e^{i\delta_{\rm CP}}&s_{23}c_{13}\\ s_{12}s_{23}-c_{12}c_{23}s_{13}e^{i\delta_{\rm CP}}&-c_{12}s_{23}-s_{12}c_{23}s_{13}e^{i\delta_{\rm CP}}&c_{23}c_{13}\end{array}\right),$$

$$\frac{d{N}_{\beta }^s}{d{E}_{\beta }}=\sum _{\alpha ,r}\int d{E}_{\alpha }{P}_{\alpha \beta }^{rs}(E_{\alpha },E_{\beta }){\varphi }_{{\nu }_{\alpha }^r}(E_{\alpha }).$$

$$P_{\alpha \beta }^{rs}(E_{\alpha },E_{\beta })$$

$$G_{\alpha \beta }^{rs}(E_{\alpha },E_{\beta })=\left(1-e^{-{\Gamma }_{i}L}\right){\left|U_{\alpha i}^{(r)}\right|}^2{\left|U_{\beta f}^{(s)}\right|}^2\frac{1}{{\Gamma }_i}\frac{d}{d{E}_{\beta }}\Gamma ({\nu }_i^r\to {\nu }_f^{s}J).$$

$$G_{\alpha \beta }^{rs}(E_{\alpha },E_{\beta })={\left|U_{\alpha i}^{(r)}\right|}^2{\left|U_{\beta f}^{(s)}\right|}^2\left(\frac{1-e^{-{\alpha }_{i}L/E_{\alpha }}}{E_{\alpha }}\right)F_g^{\prime rs}(E_{\alpha },E_{\beta }),$$

$$F_g^{\prime rs}(E_{\alpha },E_{\beta })=\frac{x_{if}^2}{(x_{if}^2-1)}{F}_g^{rs}(E_{\alpha },E_{\beta })\times {\Theta }_H(E_{\alpha }-E_{\beta }){\Theta }_H(x_{if}^2{E}_{\beta }-E_{\alpha }).$$

$$F_{g_s}^{\pm \pm }(E_{\alpha },E_{\beta })$$

$$F_{g_p}^{\pm \pm }(E_{\alpha },E_{\beta })$$

$$\begin{gathered} {\chi }^2(s_{23}^2,s_{13}^2,{\delta }_{\mathrm{CP}},\Delta m_{32}^2,\alpha ,m_{\mathrm{light}})= \\ \sum _{\beta =e,\mu }\underset{\{n_x,t_x\}}{\min }\left(2\sum _{\mathrm{bins}}\left(N_{\beta ,\mathrm{fit}}-N_{\beta ,\mathrm{obs}}\log N_{\beta ,\mathrm{fit}}\right)+{\left\{\frac{n_s^2}{{\sigma }_{ns}^2}+\frac{n_b^2}{{\sigma }_{nb}^2}+\frac{t_s^2}{{\sigma }_{ts}^2}+\frac{t_b^2}{{\sigma }_{tb}^2}\right\}}_{\beta }\right). \end{gathered}$$

$$N_{\beta ,\mathrm{fit}}=\sum _{\begin{array}{c}x=s,b\\ r=+,-\end{array}}{\left(1+n_x+t_x\frac{E_{\mathrm{bin}}-\hat{E}}{E_{\max }-E_{\min }}\right)}_{\beta }{N}_{\beta ,x}^r,$$

$${\chi }^2=\sum _{\begin{array}{c}\mathrm{bins}\\ r=+,-\end{array}}{\left(\frac{(1+n_s)N_{\mu ,s}^r+(1+n_b)N_{\mu ,b}^r-N_{\mu ,\mathrm{obs}}}{{\sigma }_{\mathrm{bin}}}\right)}^2+\frac{n_s^2}{{\sigma }_{ns}^2}+\frac{n_b^2}{{\sigma }_{nb}^2}.$$

$${\chi }_{\mathrm{reactor}}^2={\left(\frac{s_{13}^2-s_{13,\mathrm{reactor}}^2}{{\sigma }_{13,\mathrm{reactor}}}\right)}^2,$$

$${\alpha }_{T2K}^{(s)}<6.3\times 10^{-5}{\mathrm{eV}}^2,$$

$${\alpha }_{TK2+\mathrm{MINOS}}^{(s)}<7.8\times 10^{-5}{\mathrm{eV}}^2,$$

$$G_{\alpha \beta }^{rs}(E_{\alpha },E_{\beta })$$

$${\alpha }_{⟨im⟩}$$

$${\Delta }_{ijmn}^{sr\pm }$$

$$G_{\alpha\beta}^{rs}(E_{\alpha},\,E_{\beta})=\left(\frac{1-e^{-\alpha_{3}L/E_{\alpha}}}{\alpha_{3}/E_{\alpha}}\right)\left|U_{\alpha 3}^{(r)}\right|^{2}\Bigg{\{}\sum_{j=1}^{2}\left|U_{\beta j}^{(s)}\right|^{2}\frac{d}{dE_{\beta}}\Gamma(\nu_{3}^{r}\to\nu_{j}^{s}J)\\ +\left.2\Re e\left[U_{\beta 2}^{(s)}U_{\beta 1}^{(s)*}\exp\left(i\frac{\omega^{sr}_{3231}}{2}\right)\right]\sqrt{\frac{d}{dE_{\beta}}\Gamma(\nu_{3}^{r}\to\nu_{1}^{s}J)\frac{d}{dE_{\beta}}\Gamma(\nu_{3}^{r}\to\nu_{2}^{s}J})\right\}~{}.$$

$$G_{\alpha \beta }^{rs}(E_{\alpha },E_{\beta })=\left(\frac{1-e^{-{\alpha }_{3}L/E_{\alpha }}}{{\alpha }_3/E_{\alpha }}\right){\left|U_{\alpha 3}^{(r)}\right|}^2{\left|U_{\beta 1}^{(s)}\right|}^2\frac{d}{d{E}_{\beta }}\Gamma ({\nu }_3^r\to {\nu }_1^{s}J).$$

$$U\cdot{\rm diag}\left(1,\,e^{i\tfrac{\gamma_{22}}{2}},\,e^{i\tfrac{\gamma_{33}}{2}}\right)\to U\cdot\left(\begin{array}[]{ccc}0&1&0\\ 0&0&1\\ 1&0&0\end{array}\right)\cdot{\rm diag}\left(e^{i\tfrac{\gamma_{33}}{2}},\,1,\,e^{i\tfrac{\gamma_{22}}{2}}\right)~{}.$$

$${\mathcal{L}}_{\mathrm{int}}=\frac{(g_s{)}_{ij}}{2}{\bar{\nu }}_i{\nu }_{j}J+i\frac{(g_p{)}_{ij}}{2}{\bar{\nu }}_i{\gamma }_5{\nu }_{j}J,$$

$$\Gamma ({\nu }_i^{(\pm )}\to {\nu }_f^{(\pm )}J)$$

$$\Gamma ({\nu }_i^{(\pm )}\to {\nu }_f^{(\mp )}J)$$

$$f(x)$$

$$g(x)$$

$$k(x)$$

$$\frac{d}{d{E}_{\beta }}\Gamma ({\nu }_i^r\to {\nu }_f^{s}J)=\frac{m_i{m}_f}{4\pi E_{\alpha }^2}{\left(1-\frac{m_i^2}{E_{\alpha }^2}\right)}^{-1/2}{\left|\mathcal{M}({\nu }_i^r\to {\nu }_f^{s}J)\right|}^2\Theta (E_{\alpha },E_{\beta }).$$

$$\cos \theta =\frac{2E_{\alpha }{E}_{\beta }-(m_i^2+m_f^2)}{2|{\vec{p}}_i||{\vec{p}}_f|},$$

$$\frac{E_{\alpha }}{2}\left(1+\frac{m_i^2}{m_f^2}\right)-\frac{|{\vec{p}}_i|}{2}\left(1-\frac{m_i^2}{m_f^2}\right)\le E_{\beta }\le \frac{E_{\alpha }}{2}\left(1+\frac{m_i^2}{m_f^2}\right)+\frac{|{\vec{p}}_i|}{2}\left(1-\frac{m_i^2}{m_f^2}\right).$$

$$\frac{E_{\alpha }}{x_{if}^2}\le E_{\beta }\le E_{\alpha }.$$

$$\Theta (E_{\alpha },E_{\beta })={\Theta }_H(E_{\alpha }-E_{\beta }){\Theta }_H(x_{if}^2{E}_{\beta }-E_{\alpha }).$$

$${\alpha }_{if}=E_{\alpha }\left(\Gamma ({\nu }_i\to {\nu }_{f}J)+\Gamma ({\nu }_i\to {\bar{\nu }}_{f}J)\right),$$

$$(g_p{)}_{if}^2$$

Full text

1 Introduction

Experiments carried out in the last two decades have categorically confirmed the neutrino oscillation phenomenon as the mechanism responsible for the flavour transitions observed in the neutrino fluxes coming from the sun Cleveland:1998nv ; Abdurashitov:1999zd ; Altmann:2000ft ; Fukuda:2001nj ; Ahmad:2002ka ; Davis:2003xx , the atmosphere Fukuda:1998mi ; Ambrosio:2001je ; Kajita:2016vhj , reactors Araki:2004mb ; An:2012eh ; Abe:2011fz ; Ahn:2012nd and accelerators Adamson:2007gu ; Abe:2012gx ; Adamson:2013whj ; Adamson:2017qqn . These observations confirm that at least two neutrinos are massive, such that two independent mass differences $\Delta m^{2}_{ij}=m^{2}_{i}-m^{2}_{j}\;(i,j=1,\,2,\,3)$ are non-zero GonzalezGarcia:2002dz ; McDonald:2004rx ; Kajita:2016vhj . In addition, interaction and mass eigenstates are related through the PMNS lepton mixing matrix Pontecorvo:1957cp ; Maki:1962mu , defined in terms of three mixing angles and one CP-violation phase $\delta_{\rm CP}$:

[TABLE]

where $s_{ij}=\sin\theta_{ij}$ and $c_{ij}=\cos\theta_{ij}$.

Regardless of the success of the oscillation phenomenon, there is still room for new physics affecting neutrino phenomenology. In particular, several works have studied non-standard interactions GonzalezGarcia:1998hj ; Bergmann:2000gp ; Guzzo:2004ue ; Gago:2001si ; Gago:2001xg ; Fogli:2007tx ; Ohlsson:2012kf ; Esmaili:2013fva , decoherence in oscillations Gago:2000qc ; Gago:2002na ; Fogli:2003th ; Morgan:2004vv ; Lisi:2000zt ; Hooper:2004xr ; Farzan:2008zv ; Oliveira:2010zzd ; robertoPhd ; Oliveira:2013nua ; Berryman:2014yoa and fast neutrino decay Frieman:1987as ; Raghavan:1987uh ; Berezhiani:1991vk ; Berezhiani:1992ry ; Berezhiani:1993iy ; Barger:1999bg ; Beacom:2002cb ; Joshipura:2002fb ; Bandyopadhyay:2002qg ; Ando:2004qe ; Fogli:2004gy ; PalomaresRuiz:2005vf ; GonzalezGarcia:2008ru ; Maltoni:2008jr ; Baerwald:2012kc ; Meloni:2006gv ; Das:2010sd ; Dorame:2013lka ; Gomes:2014yua ; Berryman:2014qha ; Berryman:2014yoa ; Picoreti:2015ika ; Abrahao:2015rba ; Bustamante:2016ciw , among other new physics effects Pantaleone:1992ha ; Bustamante:2010nq ; Barger:2000iv ; Colladay:1998fq ; Esmaili:2014ota .

In the last years, the interplay between oscillations and fast neutrino decay has been studied in a model-independent manner Lindner:2001fx ; PalomaresRuiz:2005vf . 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 Barger:1999bg ; Beacom:2002cb ; Joshipura:2002fb ; Bandyopadhyay:2002qg ; Ando:2004qe ; Fogli:2004gy ; GonzalezGarcia:2008ru ; Maltoni:2008jr ; Baerwald:2012kc ; Meloni:2006gv ; Das:2010sd ; Dorame:2013lka ; Gomes:2014yua ; Berryman:2014qha ; Berryman:2014yoa ; Picoreti:2015ika ; Abrahao:2015rba . Throughout this work, we refer to this scenario as invisible decay. In this situation, one can bound the neutrino lifetime-to-mass ratio $\tau/m$ or, alternatively, the parameter $\alpha=E\,\Gamma$, related to the neutrino width $\Gamma$ 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 Berryman:2014qha ; Picoreti:2015ika have constrained the lifetime of $\nu_{2}$, giving $\tau_{2}/m_{2}\geq 7.2\times 10^{-4}~{}{\rm s}\cdot{\rm eV}^{-1}$ at $99\%$ C.L., equivalent to $\alpha_{2}<9.1\times 10^{-13}$ eV2.

• •

In GonzalezGarcia:2008ru , invisible decay was taken into account within atmospheric neutrino experiments. The authors combined this data with old data from MINOS Weber:2008zzc , bounding the $\nu_{3}$ lifetime. The study gave $\tau_{3}/m_{3}>3.0\times 10^{-10}~{}{\rm s}\cdot{\rm eV}^{-1}$ at $90\%$ C.L., equivalent to $\alpha_{3}<2.2\times 10^{-6}$ eV2.

Although not as competitive as GonzalezGarcia:2008ru , one should also take into account the work in Gomes:2014yua . Here, MINOS and T2K $\nu_{\mu}$ disappearance data was used to constrain invisible decay, giving a limit of $\tau_{3}/m_{3}\geq 2.8\times 10^{-12}~{}{\rm s}\cdot{\rm eV}^{-1}$, at $90\%$ C.L.. This corresponds to $\alpha_{3}<2.4\times 10^{-4}$ eV2.

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 Chikashige:1980qk ; Gelmini:1983ea ; Schechter:1981cv ; Dias:2005jm . 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 $\alpha$ decay constant by $g^{2}_{s,p}\sim 16\pi\alpha/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 $\nu_{\mu}\to\nu_{\mu}$ disappearance channels of both experiments were used in Gomes:2014yua 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 $\nu_{\mu}\to\nu_{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 Kachelriess:2000qc ; Pasquini:2015fjv ; Gando:2012pj ; Agostini:2015nwa ; Hannestad:2005ex ; Archidiacono:2013dua . If one assumes that these can be directly translated as bounds on the decay constant $\alpha$, and assuming for illustration a neutrino mass $m=0.01$ eV, they would imply constraints as low as $\alpha<\mathcal{O}\left({10^{-25}}\right)$ eV2. 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 Gomes:2014yua , 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 $\alpha$. 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.

2 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.

To calculate the flux of $\nu_{\beta}$ and $\bar{\nu}_{\beta}$ arriving at the detector, with energy $E_{\beta}$, due to the oscillation and decay of unstable $\nu_{\alpha}$ and $\bar{\nu}_{\alpha}$ ($\alpha,\beta=e,\mu$), we use:

[TABLE]

Here, $\phi_{\nu_{\alpha}^{r}}(E_{\alpha})$ is the original flux of $\nu^{r}_{\alpha}$ with energy $E_{\alpha}$, and $P_{\alpha\beta}^{rs}(E_{\alpha},\,E_{\beta})$ describes the neutrino oscillation and decay process to $\nu^{s}_{\beta}$. The $r,s=\{(+),\,(-)\}$ indices refer to the neutrino chirality (i.e. $\nu^{(-)}$ for $\nu$, and $\nu^{(+)}$ for $\bar{\nu}$). The evolution function $P_{\alpha\beta}^{rs}(E_{\alpha},\,E_{\beta})$ is obtained by adapting a relevant formula in PalomaresRuiz:2005vf :

[TABLE]

where $U^{(-)}=U_{\rm PMNS}$ and $U^{(+)}=U^{*}_{\rm PMNS}$. The Lorentz-invariant variable $\alpha_{i}$ is related to the total neutrino decay width $\Gamma_{i}$, such that $\alpha_{i}=E_{\alpha}\,\Gamma_{i}$.

The first term in the evolution function includes the standard oscillation contributions, multiplied by an exponential governed by the neutrino decay parameter $\alpha_{i}$. For stable neutrino mass eigenstates, all $\alpha_{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_{\alpha\beta}^{rs}(E_{\alpha},\,E_{\beta})$ can be found on Appendix A. In this work we use a simplified version, considering only one decay channel ($\nu^{r}_{i}\to\nu^{s}_{f}J$, where $\nu^{r}_{i}$ and $\nu^{s}_{f}$ are on the mass basis). It has the following form:

[TABLE]

We can understand Eq. (2.3) in the following way. The term $(1-e^{-\Gamma_{i}L})\left|U_{\alpha i}\right|^{2}$, when multiplied by the original flux $\phi_{\nu_{\alpha}^{r}}$, gives the amount of $\nu^{r}_{i}$ mass eigenstates that decay. When the latter is multiplied by the normalized spectrum, $\tfrac{1}{\Gamma_{i}}\tfrac{d}{dE_{\beta}}\Gamma(\nu_{i}^{r}\to\nu_{f}^{s}\,J)$, we obtain the number of $\nu_{f}^{s}$ mass eigenstates, with energy between $E_{\beta}$ and $E_{\beta}+dE_{\beta}$, produced from the decay. The amount of final $\nu_{\beta}^{s}$ flavour eigenstates is then obtained by changing back to the flavour basis with $\left|U_{\beta f}\right|^{2}$.

The full expression for $d\Gamma(\nu_{i}^{r}\to\nu_{f}^{s}\,J)/dE_{\beta}$ 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:

[TABLE]

with the visible decay function given by

[TABLE]

Here we have $x_{if}=m_{i}/m_{f}>1$, the label $g=\{g_{s},\,g_{p}\}$ indicates the non-vanishing coupling, $\Theta_{H}(x)$ is the Heaviside function and we have replaced the total width $\Gamma_{i}\to\alpha_{i}/E_{\alpha}$. The $F^{rs}_{g}(E_{\alpha},\,E_{\beta})$ functions have chirality conserving ($r=s$) and chirality changing ($r\neq s$) cases:

[TABLE]

For a given $x_{ij}$ and energy, we have that $\pm\pm$ and $\pm\mp$ transitions are complementary, which is expected, since $F_{g}^{\pm\pm}+F_{g}^{\pm\mp}=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^{\prime rs}_{g}$, defined in Eq. (2.5), for fixed $E_{\alpha}=3$ GeV, and for several values of $x_{if}$, as a function of $E_{\beta}$. On the top (bottom) rows we show results for scalar (pseudoscalar) couplings. The left column shows chirality-conserving transitions ($\pm\pm$), while the right columns shows chirality-changing transitions ($\pm\mp$).

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_{\beta}$ close to $E_{\alpha}$, we find that $\pm\pm$ processes (left column) are favoured, while for lower energies, $\pm\mp$ 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 $\pm\pm$ transitions dominate all final energies. Meanwhile, pseudoscalar couplings have a mixed behaviour, with a clear preference for $\pm\mp$ transitions.

From now on, we assume a normal ordering scenario, with $\nu_{3}$ unstable and decaying exclusively into $\nu_{1}$. We label $\alpha_{3}\to\alpha$, $m_{1}\to m_{\rm light}$ and $(g_{s,p})_{31}\to g_{s,p}$ (see Appendix B). Notice that cosmological bounds Ade:2013zuv require $\sum_{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.

3 Experimental Setup and Fitting Procedure

3.1 T2K

The T2K experiment Abe:2015awa has been running in the latest years, sending a neutrino beam to the Super-Kamiokande (SK) detector Kajita:2016vhj , 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 $\nu_{\mu}^{(-)}$ beam to the SK detector, with a final luminosity of $7.48\times 10^{20}$ protons on target (POT). Their results in the $\nu^{(-)}_{\mu}\to\nu_{\mu}^{(-)}$ disappearance channel imposed strong bounds on the $\sin^{2}2\theta_{23}$ and $\Delta m^{2}_{32}$ parameter space Abe:2012gx ; Abe:2015awa . In addition, through the observation of $\nu_{\mu}^{(-)}\to\nu_{e}^{(-)}$ appearance channel, it provided the first indication of non-zero $\theta_{13}$ Abe:2011sj .

In the current antineutrino run, a primarily $\nu_{\mu}^{(+)}$ beam illuminates the SK detector. The luminosity released to the public corresponds to $7.47\times 10^{20}$ POT, and the combination of both runs has given hints favouring a negative value of $\delta_{\rm CP}$ ICHEP-T2K-2016 ; NOW-T2K-2016 . 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 $\nu_{\mu}^{(-)}$ is the largest component in the beam Abe:2015awa , we use the FD mode (visible and invisible contributions), for both appearance and disappearance channels. For the other components of the flux ($\nu_{\mu}^{(+)}$, $\nu_{e}^{(-)}$ and $\nu_{e}^{(+)}$), we only include the FD when the visible part is not negligible.

For the $\nu_{\mu}$ disappearance channel, we include FD for $\nu_{\mu}^{(-)}\to\nu_{\mu}^{(\pm)}$ and ID for $\nu_{\mu}^{(+)}\to\nu_{\mu}^{(+)}$. 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 $\nu_{e}$ appearance channel has many contributions. As mentioned previously, we consider the FD for $\nu^{(-)}_{\mu}\to\nu^{(\pm)}_{e}$. In this case, we also include the FD for $\nu^{(+)}_{\mu}\to\nu^{(\pm)}_{e}$, but take only the ID contribution for $\nu^{(\pm)}_{e}\to\nu^{(\pm)}_{e}$. Both CCQE and CCnQE interactions are considered part of the signal. The background consists of the $\nu^{(\pm)}_{e}\to\nu^{(\pm)}_{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 Berryman:2014yoa .

We give more details on the generation of events at T2K in Appendix C.1.

3.2 MINOS

MINOS was a long-baseline neutrino experiment Adamson:2013whj which used two detectors and was exposed to a neutrino beam produced at Fermilab (NuMI beam line Adamson:2015dkw ). 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\times 10^{20}$ POT. The beam composition was $92.9\%$ of $\nu_{\mu}^{(-)}$, $5.8\%$ of $\nu_{\mu}^{(+)}$, and $1.3\%$ of $\nu_{e}^{(\pm)}$. We use the data set that comprises the charged current (CC) contained-vertex neutrino and anti-neutrino disappearance data Adamson:2013whj ; Cao:2014eca .

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 $\nu_{\mu}^{(-)}$ and $\nu_{\mu}^{(+)}$ beam components in our analysis, separately.

In order to probe the decay framework against the $\nu_{\mu}^{(-)}$ disappearance data we consider the FD for the $\nu_{\mu}^{(-)}\rightarrow\nu_{\mu}^{(-)}$ channel and only the VD $\nu_{\mu}^{(+)}\rightarrow\nu_{\mu}^{(-)}$ contribution. Analogously, for $\nu_{\mu}^{(+)}$ disappearance data, we take into account the FD for $\nu_{\mu}^{(+)}\rightarrow\nu_{\mu}^{(+)}$ and the VD $\nu_{\mu}^{(-)}\rightarrow\nu_{\mu}^{(+)}$.

The details about the MINOS reconstruction data can be found in Appendix C.2.

3.3 Statistical Analysis for T2K and MINOS

The relevant parameters in this study, which we shall vary, are $s^{2}_{23}$, $s^{2}_{13}$, $\delta_{\rm CP}$, $\Delta m^{2}_{32}$, $\alpha$ and $m_{\rm light}$. We keep fixed $s^{2}_{12}=0.306$ and $\Delta m^{2}_{21}=7.5\times 10^{-5}$ eV2.

To perform the fit for T2K, we use a $\chi^{2}$ function similar to the one used in Abe:2015awa . 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 $\chi^{2}$ is minimized with respect to $n_{x}$ and $t_{x}$, adding appropriate pull factors Huber:2002mx ; Fogli:2002pt :

[TABLE]

Here, $N_{\beta,\,\rm 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 $\nu_{\beta}^{(\pm)}$ events as $N^{(\pm)}_{\beta,\,s}$ and $N^{(\pm)}_{\beta,\,b}$, respectively, we have

[TABLE]

where $E_{\rm min}$ ($E_{\rm max}$) is the minimum (maximum) energy in the analyzed neutrino spectrum, $E_{\rm bin}$ is the average bin energy, and $\hat{E}=\tfrac{1}{2}(E_{\rm max}+E_{\rm min})$ is the average energy of the spectrum.

The other parameters appearing in Eq. (3.1) are $N_{\beta,\,\rm obs}$, which is the observed number of $\nu_{\beta}^{(-)}$ and $\nu_{\beta}^{(+)}$ events ICHEP-T2K-2016 ; NOW-T2K-2016 , and $\sigma_{nx}$ and $\sigma_{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 $\chi^{2}$ function:

[TABLE]

The parameter $\sigma_{\rm bin}$ represents the statistical and systematic uncertainty extracted from data set, while $\sigma_{ns}$, $\sigma_{nb}$ are equal to $14.7\%$ and $4\%$, respectively Adamson:2007gu .

Finally, we can also include information from reactor data. For example, the Daya Bay experiment An:2015nua has given the most precise bound on $s^{2}_{13}$, which is $\sin^{2}2\theta_{13}=0.092\pm 0.017$ An:2012eh . 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 An:2016ses . We then assume that no reactor disappearance experiment is affected by neutrino decay. Thus, only for our final result constraining the decay parameter $\alpha$, we include an additional pull factor on our $\chi^{2}$ function:

[TABLE]

with $s_{13,\rm reactor}^{2}=0.0243$ and $\sigma_{13,\rm reactor}=0.0026$, following the same procedure as in Abe:2015awa .

In our fit, we vary $m_{\rm light}$ and $\alpha$ logarithmically, from $5\times 10^{-3}$ to $10^{-1}$ eV and $2\times 10^{-6}$ to $5\times 10^{-4}$ eV2, respectively. The oscillation angles $\theta_{13}$ and $\theta_{23}$ are also scanned, with $0\leq s_{13}^{2}\leq 0.05$ and $0.32\leq s_{23}^{2}\leq 0.86$. The phase $\delta_{\rm CP}$ is scanned between $-\pi$ and $\pi$, while the squared mass difference $\Delta m^{2}_{32}$ is varied between $1.6\times 10^{-3}$ and $2.8\times 10^{-3}$ eV2.

4 Impact on Oscillation Parameter Regions

4.1 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 $\Delta m^{2}_{32}=2.545\times 10^{-3}$ eV2, which correspond to the best-fit parameters for T2K data. We have set $\alpha=4\times 10^{-5}$ eV2, which is roughly $10\%$ of the mean lifetime that T2K should be sensitive to ($\approx E/L$). We have also set $m_{\rm light}=7\times 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 $\alpha$, there is no big difference between the SO and FD scenarios. Nevertheless, even though not shown, one finds that for larger values of $\alpha$ 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 $\delta_{\rm CP}=\pm\pi/2$. Notice that one obtains a larger amount of events when $\delta_{\rm CP}=-\pi/2$. In contrast, we show ID and FD scenarios only with $\delta_{\rm CP}=+\pi/2$. A first important feature is that, although the ID contribution is similar to the SO result (for the same value of $\delta_{\rm 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 $\delta_{\rm CP}=\pi/2$, can have a similar number of events as the SO result with $\delta_{\rm CP}=-\pi/2$ at the cost of bringing some distortion to the low-energy (high-energy) part of the spectrum, for scalar (pseudoscalar) couplings.

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 $\alpha$, neutrino decay does not significantly modify the spectrum.

On the other hand, as in Figure 2, for the appearance spectrum we show SO with $\delta_{\rm CP}=\pm\pi/2$. On this case, $\delta_{\rm CP}=-\pi/2$ leads to a smaller number of events compared to $\delta_{\rm CP}=+\pi/2$. We again show the spectrum for ID and FD with $\delta_{\rm CP}=+\pi/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 $\delta_{\rm CP}=-\pi/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 $\delta_{\rm 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_{\rm light}$. In the previous section (Figure 1), we found that for scalar couplings, large values of $m_{\rm light}$ ($x_{31}\to 1$) favour $\nu^{(\pm)}\to\nu^{(\pm)}$ over $\nu^{(\pm)}\to\nu^{(\mp)}$ 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 $\nu^{(-)}$ and $\nu^{(+)}$ 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 $\nu^{(+)}$, 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 $\nu^{(-)}$, 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_{23}^{2}-\Delta 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_{13}^{2}-\delta_{\rm 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 $\delta_{\rm 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 $\alpha$ as a correlation with $s^{2}_{13}$ and $m_{\rm light}$.

On the bottom left panel, we see that for low values of $\alpha$, the allowed range for $s^{2}_{13}$ is well bounded (consistent with SO). Then, as $\alpha$ 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_{13}^{2}$, is compensated by the extra neutrinos coming from decay. Finally, in both cases, too large values of $\alpha$ generate too many appearance events, which is incompatible with T2K data, restricting the allowed $\alpha$ to values of $\mathcal{O}\left({10^{-5}}\right)$ eV2. This automatically leaves unaffected the $s_{23}^{2}-\Delta m^{2}_{32}$ subspace, which requires much larger values of $\alpha$ to modify significantly the spectrum.

On the other hand, on the bottom right panel, one finds the confidence level contours for $\alpha$ as a function of $m_{\rm light}$. For scalar (pseudoscalar) couplings, we find that for large values of $m_{\rm light}$, the contours reach higher (lower) values of $\alpha$. 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 $\alpha$ are compatible with data. Opposed to the latter case, for scalar couplings with large $m_{\rm light}$, chirality-changing transitions are minimized, so visible decay has a smaller impact, and larger values of $\alpha$ are allowed. In contrast, for smaller $m_{\rm light}$, both scalar and pseudoscalar curves tend to the same value of $\alpha$. 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 $\alpha$. As we shall see later, these points are not statistically significant. Thus, from the T2K-only $\chi^{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 $\alpha$ (see Figure 8).

4.2 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$, $\Delta m^{2}_{32}=2.41\times 10^{-3}~{}\text{eV}^{2}$, which correspond to the best fit of MINOS to standard oscillation model, and $\delta_{\rm CP}=0$, because it is not sensitive in MINOS. When the decay is present, we have set $\alpha=3.0\times 10^{-4}~{}{\rm eV^{2}}$, and $m_{\text{light}}=0.05~{}\text{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_{\text{light}}$, this means that for scalar couplings we have a preference for $\pm\pm$ transitions. In contrast, pseudoscalar couplings have a clear preference for $\pm\mp$ transitions. 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 $\nu_{\mu}^{(-)}\rightarrow\nu_{\mu}^{(-)}$. 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 $\nu_{\mu}^{(-)}\rightarrow\nu_{\mu}^{(+)}$ is favored for the value of $m_{\text{light}}$ chosen.

The best fit obtained for ID and FD scenario were for an $\alpha$ 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. Gomes:2014yua . However, if we take the minimum $\chi^{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 Akaike1 ; Akaike2 that there is no statistical difference between both models.

In Figure 6, we show the 1$\sigma$ and 90$\%$ C.L. allowed regions for many combination of the relevant parameters, already computing $\nu^{(-)}_{\mu}$ combined with $\nu^{(+)}_{\mu}$. The red (blue) curves refer to scalar (pseudoscalar) couplings. To compare the effect of neutrino decay, we add as well the limit when $\alpha\rightarrow 0$, which corresponds to SO scenario (full gray) in Eq. (2.2).

The allowed region in the ($s^{2}_{23}-\Delta 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}-\Delta m^{2}_{32}$) plane. We can see it for small values of $\Delta m^{2}_{32}$ and higher for $s^{2}_{23}$.

We also present the two-dimensional projections ($\alpha-\Delta m^{2}_{32}$) and ($\alpha-s^{2}_{23}$) in Figure 6. We can observe that for large values of $\alpha$, the contours increase the range in $\Delta m^{2}_{32}$, and either allow higher values of $s^{2}_{23}$.

The right bottom panel in Figure 6 show the correlation between $\alpha$ and $m_{\text{light}}$. We can conclude by it, that MINOS is not able to distinguish the $\alpha$ constraint for scalar and pseudoscalar couplings.

4.3 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 $\chi^{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}-\Delta 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 $\alpha$ being strongly constrained by $\nu_{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 $\nu_{\mu}$ disappearance. Therefore, the T2K analysis dominates the fit.

On the top left panel, we show the $s^{2}_{13}-\delta_{\rm 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 $\delta_{\rm CP}$ is ruled out at $90\%$ C.L.

The qualitative behaviour observed on the bottom left panel, showing the $s^{2}_{13}-\alpha$ plane, is very similar to the one of Figure 4. An evident difference is that, for both couplings, the best-fit for $\alpha$ is increased to a value of $\mathcal{O}\left({10^{-5}}\right)$. 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 $\delta_{\rm CP}$ are allowed in the scalar coupling scenario. In the SO scenario, if we fix $s_{13}^{2}$ at the reactor best-fit value, positive $\delta_{\rm 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 $\delta_{\rm CP}$ is incompatible with a vanishing value of $\alpha$. 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 $\delta_{\rm 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 $\alpha$, for all values of $\delta_{\rm CP}$ for both couplings, reintroducing the ambiguity in $\delta_{\rm CP}$ for the scalar case. This ambiguity with $\delta_{\rm CP}$ does not happen for the pseudoscalar case, as we have seen, the incompatibility of a value of $\alpha$ of $\mathcal{O}\left({10^{-5}}\right)$, for $\delta_{\rm CP}=+\pi/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 $\delta_{\rm CP}$ solution would eventually also be excluded for scalar coupling.

In Figure 8 we present show $\Delta\chi^{2}$ as a function of $\alpha$, 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.

As we can see, the MINOS solution for non-zero $\alpha$ is excluded by T2K data, which is compatible with a null value for this parameter. Therefore, we place bounds on $\alpha$. The T2K-only $90\%$ C.L. bounds are $\alpha\lesssim 5.6\times 10^{-5}$ eV2 ($\alpha\lesssim 6.3\times 10^{-5}$ eV2) for the pseudoscalar (scalar) scenario. As expected, the constraint on this parameter for the T2K+MINOS combination is slightly weaker: $6.9\times 10^{-5}$ eV2 ($7.8\times 10^{-5}$ eV2) for the pseudoscalar (scalar) case.

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.

5 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 $\alpha$ as a best-fit solution. Nevertheless, the inclusion of T2K neutrino and antineutrino data excludes this result, being consistent with a vanishing $\alpha$. 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 $\alpha$ 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 $\alpha$, given the coupling $c$, due to the experiment $exp$, by $\alpha^{(c)}_{exp}$, the best constraint we find at $90\%$ C.L. is:

[TABLE]

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.

Acknowledgements.

A.M.G. and J.J.P. acknowledge funding by the Dirección de Gestión de la Investigación at PUCP, through grant DGI-2015-3-0026. They would also like to thank S. Sanchez for discussions on the use of GENIE, as well as O. A. Díaz in the Dirección de Tecnologías de Información at PUCP, for implementing our code within the LEGION system. A.M.G is thankful for the Unicamp hospitality and for the support of the FAPESP funding grant 2015/20129-3. A.M.G and O.L.G.P are also thankful to the Oficina de Internacionalización de la Investigación at PUCP. A.L.G.G. is supported by CAPES, and R.A.G. would like to thanks CNPq grants 461147/2014-9 and 311980/2013-9. O.L.G.P. is thankful for the PUCP hospitality and partial support under the Programa Latino Americano de Física (PLAF) from Brazilian Physics Society. He is also thankful for the support of FAPESP funding grant 2014/19164-6 and CNPq research fellowship  307269/2013-2 and 304715/2016-6 and partial funding from ICTP.

Appendix A Full Formula for Visible Neutrino Decay

The $G_{\alpha\beta}^{rs}(E_{\alpha},\,E_{\beta})$ function describing visible neutrino decay can be derived from the original formula of PalomaresRuiz:2005vf , using the methods of Lindner:2001fx . Assuming the final state neutrinos do not oscillate, the result is:

[TABLE]

Here, we have $K_{jinm}^{\beta\alpha sr}=U^{(s)}_{\beta j}U^{(r)*}_{\alpha i}U^{(s)*}_{\beta n}U^{(r)}_{\alpha m}$. Notice that, following Eq. (1.1), the PMNS matrix has been defined without the Majorana phases. In Eq. (A), the latter are taken into account within the parameter $\omega_{ijmn}^{sr}=\gamma^{s}_{jj}-\gamma^{r}_{ii}-\gamma^{s}_{nn}+\gamma^{r}_{mm}$. We also denote:

[TABLE]

If we only have one unstable neutrino, this expression can be considerably simplified. We chose $\nu_{3}^{r}$ to be unstable, which means we fix $i=m=3$ in Eq. (A), so we obtain:

[TABLE]

Furthermore, if only one decay channel is allowed (for example, $\nu_{3}^{r}\to\nu_{1}^{s}J$, obtained by fixing $j=n=1$ in Eq. (A)), we have:

[TABLE]

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}\to m_{3}$, $m_{2}\to m_{1}$ and $m_{3}\to m_{2}$, or modify the full PMNS matrix such that:

[TABLE]

Appendix B Formulae for Visible Neutrino Decay Rates

In the following, we shall concentrate on the differential decay width $d\Gamma(\nu_{i}^{r}\to\nu_{f}^{s}J)/dE_{\beta}$, for the general case where several decay channels exist. Our procedure follows Lindner:2001fx , where we have a neutrino - Majoron coupling:

[TABLE]

where the neutrinos are on the mass eigenstates basis.

With this, the partial decay rates of a neutrino with energy $E_{\alpha}$ are Lindner:2001fx ; Kim:1990km :

[TABLE]

where $x_{if}=m_{i}/m_{f}>1$, and:

[TABLE]

For our description of neutrino decay, we need the differential decay rate:

[TABLE]

The $\Theta(E_{\alpha},\,E_{\beta})$ function fixes the angle between the momenta of the incoming and outgoing neutrinos:

[TABLE]

which effectively imposes the following bound on $E_{f}$:

[TABLE]

On the relativistic limit, this reduces to:

[TABLE]

Thus, in terms of the Heaviside function $\Theta_{H}$, we have:

[TABLE]

The neutrino decay amplitudes $\mathcal{M}(\nu_{i}^{(\pm)}\to\nu_{f}^{(\pm)}J)$ can be computed as a function of $(g_{p})^{2}_{if}$ and $(g_{s})^{2}_{if}$, the neutrino masses and the neutrino energies Lindner:2001fx ; Kim:1990km . These can be recasted using the partial width of neutrinos, given by Eqs. (B.2). We first define:

[TABLE]

such that $\alpha_{i}=E_{\alpha}\Gamma_{i}=\sum_{f}\alpha_{if}$. We can invert the Eq. (B.9) and solve for one coupling, such that either the scalar or the pseudoscalar couplings is written:

[TABLE]

Using Eqs. (B.10), we can re-write the expression of the amplitude $\mathcal{M}(\nu_{i}^{(\pm)}\to\nu_{f}^{(\pm)}J)$ Lindner:2001fx ; Kim:1990km in terms of $\alpha_{ij}$ and one coupling. For example, in terms of $(g_{if})$, we have:

[TABLE]

Alternatively, in terms of $(g_{p})_{if}$, we find:

[TABLE]

For both cases, we have defined:

[TABLE]

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 $\alpha_{if}$. If, in addition, there exists only one allowed decay channel, then the whole process is governed by $\alpha_{i}$.

These sort of considerations lead us to the principal formula used in this work, Eq. (2.4).

Appendix 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 $\nu_{\beta}^{s}$ cross-section and perform the $E_{\beta}$ 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:

[TABLE]

where $\sigma_{\beta}^{s,\rm int}(E_{\beta})$ denotes the appropriate cross section, and the “detection kernel” $K^{\rm int}_{i}(E_{\beta})$ for bin $i$ is defined:

[TABLE]

Here, $\epsilon^{\rm int}_{\beta}(E_{\rm bin})$ denotes the detector efficiency, while the resolution function $R(E^{\rm int}_{\beta,\rm rec}-E_{\beta},\varsigma_{\beta}^{\rm int})$ reflects our capacity to reconstruct the true neutrino energy $E_{\beta}$. In both experiments we model this using a gaussian function:

[TABLE]

The variable $E^{\rm int}_{\beta,\rm rec}$ is related to the bin energy $E_{\rm bin}$ by a possible energy shift, and $\varsigma^{\rm int}_{\beta}$ denotes the energy resolution width.

C.1 T2K

For T2K, the neutrino fluxes appearing in Eq. (2.1) are obtained from Abe:2015awa , and the cross sections are taken from GENIE Andreopoulos:2009rq .

For $\nu_{\mu}$ disappearance, we take the shift and resolution width $\varsigma_{\mu}^{\rm int}$ from Machado:2011ar :

[TABLE]

with equal values for neutrino and antineutrino channels.

For $\nu_{e}$ appearance, we use the “migration matrix” from Machado:2013kya , which is equivalent to using:

[TABLE]

For NC events, we fit the shape of the spectra shown in ICHEP-T2K-2016 ; Salzgeber:2015gua . We set:

[TABLE]

Finally, the detector efficiency functions are obtained by fitting the number of events on each bin in Figures 20, 28 and 41 in Abe:2015awa as well as relevant plots in ICHEP-T2K-2016 ; Salzgeber:2015gua ; NOW-T2K-2016 .

C.2 MINOS

The $\phi^{r}_{\nu_{\alpha}}(E_{\alpha})$ 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_{\alpha},E_{\beta}})$, following:

[TABLE]

where $E_{\alpha}$ is the bin energy and $E_{\beta}$ are the adjacent bins of $E_{\alpha}$. The ND flux $\phi^{(\text{ND})_{r}}_{\nu_{\alpha}}(E_{\alpha})$ was taken from references Adamson:2007gu ; Adamson:2009ju . To construct the Beam Matrix $f^{F/N}(E_{\alpha},E_{\beta})$, we used the matrix $M(E_{\alpha},E_{\beta})$ extracted from reference JETP-MINOS-saoulidou-2007 , 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_{\alpha},E_{\beta}})$ can be described as:

[TABLE]

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 $\sigma^{s,\rm CC}_{\beta}$ are taken from references Kyberd:2012iz ; Hewett:2012ns . We have taken $E^{\rm CC}_{\mu,\rm rec}=E_{\rm bin}$, and assume the resolution widths $\varsigma^{\rm CC}_{\beta}(E_{\beta})$ to be linear in energy. The efficiency $\epsilon^{CC}_{\beta}$ of the detector was extracted from reference Adamson:2007gu .