text:
# A combined treatment of neutrino decay and neutrino oscillations

## 1 Introduction

The most favorable alternative for fast neutrino decay is to introduce an effective decay Lagrangian which couples the neutrino fields to a massless boson carrying lepton number, i.e., a complex scalar field or a Majoron field . One possibility for such a Lagrangian for the case of Majoron decay is111We are not interested in $`\nu _i\nu _i`$ transitions described by $`g_{ii}0`$, since they are forbidden for kinematical reasons. Furthermore, because the $`g_{ij}`$’s are usually assumed to be small, only the first order processes in the Lagrangian $`_{\mathrm{int}}`$ will be relevant to our discussion.

$$_{\mathrm{int}}=\underset{ij}{\underset{i}{}\underset{j}{}}g_{ij}\overline{\nu _{j,L}^c}\nu _{i,L}J,$$
(1)
where the $`\nu `$’s are Majorana mass eigenfields and $`J`$ is a Majoron field.

Neutrino decay has been studied as an alternative to neutrino oscillations, especially for atmospheric  and solar  neutrinos. So far, either decay only or some sequential combination of decay and oscillations has been considered (e.g., MSW-mediated solar neutrino decay ). The values of individual $`g_{ij}`$’s or combinations were restricted, for instance, from double beta decay, pion and kaon decays, or supernovae . For pure neutrino decay without oscillations the constraints on the coupling constants have been used to derive conditions for other parameters. For example, for atmospheric neutrinos $`\mathrm{\Delta }m^20.73\mathrm{eV}^2`$ was found in Ref.. Since this is incompatible with the known (or usually assumed) mass hierarchy of active neutrinos, decay into a sterile neutrino was suggested. However, massive neutrinos allow neutrino decay as well as neutrino oscillations. In this paper, we present a combined model-independent treatment of neutrino decay and neutrino oscillations. We refer to two different cases:

Decay into neutrinos or other particles which are not observable (e.g., sterile decoupled neutrinos222In this paper, sterile decoupled neutrinos refers to neutrinos which can be described by a block diagonal mixing matrix such as
$$U=\left(\begin{array}{cc}U_{3\times 3}& 0\\ 0& U_{n_S\times n_S}\end{array}\right),$$
(2)  where $`n_S`$ is the number of sterile neutrinos, $`U_{3\times 3}`$ the active neutrino mixing matrix, and $`U_{n_S\times n_S}`$ the sterile neutrino mixing matrix.).

Decay into neutrinos which are, in principle, observable (e.g., active neutrinos).

Our paper is organized as follows: In Sec. 2, we will introduce a Majoron decay model for illustration, in order to show two important properties: First, as already mentioned before (e.g., in Ref.), we may be able to detect active decay products. Second, we will show that in such a model, decay products can, in principle, oscillate. In Sec. 3, we will develop a model-independent operator framework for the calculation of transition probabilities and describe its properties. Furthermore, we will show that the phase relationships given by the $`S`$-matrix approach of quantum field theory are satisfied. The transition probability for invisible decay as a generalization of the neutrino oscillation formula will be derived in Sec. 4. As an example, we will apply this to atmospheric neutrinos. In the following section, Sec. 5, this formalism will be used for visible decay in the approximation of almost stable decay products. Again, we will make an application to atmospheric neutrinos. In Sec. 6, we will, for illustration, construct a minimal model in which decay products oscillate. Finally, in Sec. 7, we will present a summary as well as our conclusions.

## 2 Majoron models

Neutrino decay through emission of a massless Nambu–Goldstone boson is currently among the most favorable alternatives for neutrino decay scenarios. Therefore, we will use it as an introductory example for our framework.

### 2.1 Dynamics of Majoron decay

For quantitative estimates we will investigate decay by Majoron emission described by an interaction Lagrangian with Yukawa (scalar or pseudoscalar) couplings (cf., Ref.),

$$_{\mathrm{int}}=\frac{g_1}{2}\overline{\nu }_1\nu _2J+\frac{g_2}{2}\overline{\nu }_1i\gamma _5\nu _2J.$$
(3)
Here $`\nu _1`$ and $`\nu _2`$ are assumed to be Majorana neutrinos. Since the weak interaction couples only to left-handed neutrinos and right-handed antineutrinos, and Majorana particles are identical to their anti-particles, Majorana neutrinos and Majorana antineutrinos can only be distinguished by their spin. The matrix elements in the observer’s rest frame are given by

$`\left|(\nu _2\nu _1)\right|^2`$ $`=`$ $`{\displaystyle \frac{g_1^2}{4}}(A+2)+{\displaystyle \frac{g_2^2}{4}}(A2),`$ (4)
$`\left|(\nu _2\overline{\nu }_1)\right|^2`$ $`=`$ $`{\displaystyle \frac{g_1^2+g_2^2}{4}}\left({\displaystyle \frac{m_1^2+m_2^2}{m_1m_2}}A\right)={\displaystyle \frac{g_1^2+g_2^2}{4}}\left({\displaystyle \frac{1}{x}}+xA\right),`$ (5)
where we have with the assumed mass hierarchy $`m_2>m_1`$ and

$`A`$ $``$ $`{\displaystyle \frac{m_1E_2}{m_2E_1}}+{\displaystyle \frac{m_2E_1}{m_1E_2}}={\displaystyle \frac{1}{x}}{\displaystyle \frac{E_2}{E_1}}+x{\displaystyle \frac{E_1}{E_2}},`$
$`x`$ $``$ $`{\displaystyle \frac{m_2}{m_1}}>1.`$
The differential decay rate for the decay $`\nu _2\nu _1`$ is

$$\frac{d\mathrm{\Gamma }_{21}}{dE_1}=\frac{m_1m_2}{4\pi }\frac{1}{E_2|𝐩_2|}|(E_1)|^2,$$
(6)
with the kinematics in spherical coordinates ($`z`$-axis in the direction of propagation)

$$E_2E_1=|𝐩_2𝐩_1|=\sqrt{𝐩_1^2+𝐩_2^22|𝐩_1||𝐩_2|\mathrm{cos}\theta }.$$
(7)
The condition $`\mathrm{cos}\theta 1`$ implies for the allowed energy range of $`E_1`$ that

$$\frac{E_2}{2}\left(1+\frac{1}{x^2}\right)\frac{|𝐩_2|}{2}\left(1\frac{1}{x^2}\right)E_1\frac{E_2}{2}\left(1+\frac{1}{x^2}\right)+\frac{|𝐩_2|}{2}\left(1\frac{1}{x^2}\right).$$
(8)
For relativistic neutrinos, i.e., $`m_2E_2`$, this reduces to

$$\frac{E_2}{x^2}E_1E_2.$$
(9)
We immediately see that for $`x1`$, i.e., $`\mathrm{\Delta }m^2m^2`$, the energy shift by decay is very small. For an energy resolution

$$\mathrm{\Delta }E_{\mathrm{res}}E_2\left(1\frac{1}{x^2}\right)$$
(10)
the decay products are detected in the same energy bin as the original particles.

Now, we can integrate the differential decay rate given in Eq. (6) over the range determined by Eq. (8). Using Eqs. (4), (5), and (2.1) as well as assuming relativistic neutrinos, i.e., $`m_2E_2`$, we obtain for the total decay rates

$`\mathrm{\Gamma }(\nu _2\nu _1)={\displaystyle \frac{m_1m_2}{16\pi E_2}}[`$ $`g_1^2`$ $`\left({\displaystyle \frac{x}{2}}+2+{\displaystyle \frac{2}{x}}\mathrm{ln}x{\displaystyle \frac{2}{x^2}}{\displaystyle \frac{1}{2x^3}}\right)`$
$`+`$ $`g_2^2`$ $`({\displaystyle \frac{x}{2}}2+{\displaystyle \frac{2}{x}}\mathrm{ln}x+{\displaystyle \frac{2}{x^2}}{\displaystyle \frac{1}{2x^3}})]`$ (11)
and

$$\mathrm{\Gamma }(\nu _2\overline{\nu }_1)=\frac{m_1m_2}{16\pi E_2}(g_1^2+g_2^2)\left(\frac{x}{2}\frac{2}{x}\mathrm{ln}x\frac{1}{2x^3}\right).$$
(12)
### 2.2 Visibility of decay products

Solving Eq. (7) for $`\mathrm{cos}\theta `$ shows that $`\mathrm{cos}\theta `$ as a function of $`E_1`$ always has a minimum at $`E_{1,\mathrm{min}}=\frac{2E_2m_1^2}{m_1^2+m_2^2}`$ (cf., Fig. 1). We can use this to constrain $`\theta `$. For relativistic neutrinos, i.e., $`m_2E_2`$, we obtain

$$\theta _{\mathrm{max}}=\frac{m_1}{2E_2}\left(x^21\right).$$
(13)
Let us assume $`m_1`$ and $`m_2`$ to be comparable by orders of magnitude, i.e., $`x1`$, which is plausible for active neutrinos for the currently assumed mass squared differences. For $`d`$, the distance to the detector, the area $`F`$ covered by a neutrino beam at the detector due to a beam-spread by decay is approximately given by

$$F\pi (\theta _{\mathrm{max}}d)^2=\pi d^2\left(\frac{m_1}{2E_2}\right)^2\left(x^21\right)^2=𝒪\left(\left(\frac{dm_1}{E_2}\right)^2\right)$$
(14)
for $`x=𝒪(1)`$. Table 1 shows typical values for several cases, where neutrino beams and not point sources are assumed. For a beam-spread at the detector (area $`F`$) much smaller than a usual detection area ($`D`$), $`FD`$, everything will be visible, i.e., in the cases of atmospheric, accelerator, and reactor neutrinos.

Let us now introduce a function $`\eta `$ which describes the fraction of the decay products that will arrive at the detector by geometry, i.e., does not escape from detection by changing direction. We will first of all consider the case of a neutrino beam, and discuss the case of a point source later. We introduce this function in a quite general form applicable to any decay model such as a model involving more than two decay products. Thus, $`\eta `$ is a function of intrinsic decay properties such as energy and mass of the original particle $`E_i`$, $`m_i`$, and of the decay products $`E_j`$, $`m_j`$, as well as geometric properties such as length of the baseline $`L`$, decay position $`l`$, and detection area perpendicular to the beam direction $`D`$ (see Fig. 2), i.e.,

$$\eta =\eta (E_i,m_i,E_j,m_j,L,l,D)\eta _{ij}(L,l,D).$$
(15)
Note that the intrinsic properties of the beam and the detector can be independently folded with the transition probability, since they are not directly affected by decay.

Since we assume that we cannot detect particles other than the neutrino decay product, we have to integrate over all momenta and energies except for the angle $`\theta `$ (cf., Fig. 2) and the energy of the decay product $`E^{}`$. In general, these two parameters are to be integrated over the relevant detection region and the energy bin:

$$\eta _{ij}(L,l,D)=\frac{1}{\mathrm{\Gamma }_{ij}^{\mathrm{tot}}}\left(\underset{E_{\mathrm{min}}}{\overset{E_{\mathrm{max}}}{}}\underset{(\mathrm{cos}\theta )_D}{\overset{1}{}}\left|\frac{d^2\mathrm{\Gamma }_{ij}}{d\mathrm{cos}\theta dE^{}}\right|𝑑E^{}d\mathrm{cos}\theta \right).$$
(16)
Here $`[E_{\mathrm{min}},E_{\mathrm{max}}]`$ is the energy bin range and $`(\mathrm{cos}\theta )_D`$ is the cosine of the maximal angle, such that the decay product will still hit the detector. Since we only have one independent parameter in Majoron decay, the double-differential decay rate $`\frac{d^2\mathrm{\Gamma }_{ij}}{d\mathrm{cos}\theta dE^{}}`$ will reduce to $`\frac{d\mathrm{\Gamma }}{dE^{}}`$ times a $`\delta `$-distribution determined by Eq. (7). Assuming that the energy bin width $`\mathrm{\Delta }E_{\mathrm{bin}}`$ and the angle $`\theta _D`$ are very small, i.e., the differential decay rate is approximately constant within the target region and the energy bin range, we obtain for $`\eta `$

$$\eta _{ij}(L,l,D)\frac{1}{\mathrm{\Gamma }_{ij}^{\mathrm{tot}}}\left|\frac{d^2\mathrm{\Gamma }_{ij}}{d\mathrm{cos}\theta dE^{}}\right|_{\mathrm{cos}\theta =1,E^{}=\overline{E}_{\mathrm{bin}}}\mathrm{\Delta }E_{\mathrm{bin}}\mathrm{\Delta }(\mathrm{cos}\theta )_D,$$
(17)
where $`\overline{E}_{\mathrm{bin}}`$ is the mean energy of the target bin considered, $`\mathrm{\Delta }E_{\mathrm{bin}}`$ the energy bin width of that bin, and $`\mathrm{\Delta }(\mathrm{cos}\theta )_D`$ the cosine range, i.e., $`\mathrm{\Delta }(\mathrm{cos}\theta )_D=1(\mathrm{cos}\theta )_D`$, of the angle covered by the detector. From geometry (Fig. 2) we can determine $`\mathrm{\Delta }(\mathrm{cos}\theta )_D`$ for the decay position $`l`$ and a spherical detection area $`D`$:

$$2\pi (Ll)^2\mathrm{\Delta }(\mathrm{cos}\theta )_DD.$$
(18)
Table 1 shows that for solar neutrinos the spread of a beam at the detector is normally larger than the detection area. Figure 3 shows a numerical evaluation of $`\eta _{21}`$ for solar neutrino sample data and a very large detection area for illustration. One can see that in all cases the $`1/(Ll)^2`$-dependence dominates in the range of about $`610^{10}\mathrm{m}`$ to $`910^{10}\mathrm{m}`$ over the effects of the couplings or helicities, so that we obtain the same qualitative behavior. However, the Sun is a radially symmetric ($`4\pi `$) source. It is not a point source, because the production area is large compared to the distance to the Earth. Equation (13) indicates that for decay the maximal angle of deviation of the direction of flight ($`\theta 10^7`$) is much smaller than the angle under which the Sun’s core, where the neutrinos are produced, is observed from the Earth ($`\theta 10^3`$). Hence, equally distributed small deviations of the direction of flight will already be washed out by the flux profile. Besides that, after decay, the profile will still remain radially symmetric. Thus, for the case of the Sun, the function $`\eta _{21}`$ describing a beam can only serve as an example to show the qualitative behavior, and we conclude that $`\eta ^{\mathrm{Sun}}1`$.

Supernova neutrinos are different from solar neutrinos, since they come from a far-distant point source, and hence from a well-defined direction. One might argue that it is impossible to distinguish neutrinos produced by decay in the direction of flight or by decay in slightly different paths due to the radially symmetric flux profile, since decay scatters the neutrinos back into the direction of the detector. However, these two cases may be distinguishable for certain cases or sets of parameters because of the time resolution (e.g., different traveling times or different oscillation lengths due to different path lengths). Hence, supernova neutrinos might be the only nontrivial example of partly visible Majoron decay, i.e., $`\eta <1`$. In general, any far-distant beam source will show similar effects for Majoron decay.

Eventually, we conclude that Majoron decay products are in most cases visible to a good approximation, i.e., $`\eta _{ij}(L,l,D)1`$. In our calculations, we will generalize this to $`\eta _{ij}(L,l,D)\eta _{ij}(L,D)`$. Hence, $`\eta `$ will be assumed not to depend on the distance $`l`$. Moreover, for supernova neutrinos it should be a good approximation to use a step function for $`\eta _{ij}(L,l,D)`$ according to Fig. 3.

### 2.3 Neutrino oscillations of decay products

In this subsection, we demonstrate that neutrino decay products in Majoron decay may oscillate. Since neutrino oscillations are a consequence of phase coherence, several issues connected to this subject need to be addressed for a combined treatment of neutrino decay and neutrino oscillations. We will do that in detail in Subsec. 3.2.

In a charged current weak interaction, it is normally assumed that the phase coherence within a superposition of states as an in state is destroyed. This means that such a process produces a superposition of states with random relative phase shifts, i.e., the states within the superposition cannot interfere. Hence, cross sections referring to the individual states have to be summed incoherently, i.e., by summation of probabilities and not amplitudes (cf., Property 4 in Subsec. 3.2). This comes from kinematical features related to the mass differences of accompanying leptons , i.e., the coherence length is relatively small due to the large mass differences.

In a neutral current weak interaction, the out state can be a coherent superposition of states. One example is $`Z^0`$ decay, producing a coherent superposition of different neutrino and antineutrino mass eigenstates

$$|\nu _Z\frac{1}{\sqrt{3}}\left(|\overline{\nu }_1|\nu _1+|\overline{\nu }_2|\nu _2+|\overline{\nu }_3|\nu _3\right)$$
(19)
by an interaction Lagrangian

$$_{\mathrm{int}}=i\frac{g}{2}\left(\underset{i=1}{\overset{3}{}}\overline{\nu }_{iL}\gamma ^\mu \nu _{iL}\right)Z_\mu ^0.$$
(20)
Majoron decay resembles a neutral current weak interaction in the sense of a Lagrangian similar to Eq. (20). It couples a superposition of different mass eigenstates to a single particle, so that all eigenstates and their wave packets, respectively, experience the same energy shift $`\mathrm{\Delta }E`$ transferred to the massless boson. From the wave packet treatment of neutrino oscillations we know that within a certain coherence length, the mass differences of neutrinos do not destroy the coherence of the mass eigenstates, i.e., neutrino oscillations. Hence, since in Majoron decay the energies of the mass eigenstates are all shifted by the same amount $`\mathrm{\Delta }E`$ in decay, the mass differences will not be able to destroy the coherence of out states in this case either333Later in this paper, a more detailed proof of this analogy applied to our framework will be given in Property 9 in Subsec. 3.2.. Thus, the decay products can oscillate.

Using a Majoron-like decay interaction Lagrangian, we know that to first order of the $`S`$-matrix expansion we can write

$$d\mathrm{\Gamma }(inout)\left|out\left|\underset{\mathrm{}}{\overset{\mathrm{}}{}}d^4x\underset{ij}{\underset{i}{}\underset{j}{}}g_{ij}\overline{\nu }_j\nu _iJ\right|in\right|^2.$$
(21)
This expression implies that the interaction Lagrangian destroys the incoming superposition completely, and creates a superposition of outgoing mass eigenstates with coefficients determined by the coupling constants and the coefficients of the incoming superposition. The outgoing states have a fixed relative initial phase such as the outgoing mass eigenstates in $`Z^0`$ decay . In order to observe neutrino oscillations of decay products, we need to have at least two different mass eigenstates in the out state superposition. We immediately see from Eq. (21) that this implies that at least two of the $`g_{ij}`$’s, say $`g_{i_1j_1}`$ and $`g_{i_2j_2}`$, have to be non-zero for the corresponding states $`|\nu _{i_1}`$ and $`|\nu _{i_2}`$ with non-vanishing in state amplitudes $`\nu _{i_1}`$ and $`\nu _{i_2}`$.

## 3 Operators and properties

In this section, we will introduce effective decay operators for a model-independent combined treatment of neutrino decay and neutrino oscillations. Decay can be described by decay rates $`\mathrm{\Gamma }_{ij}g_{ij}^2`$, which are calculated from the corresponding Feynman diagrams. We will use such an approach to effectively integrate neutrino decay into neutrino oscillation scenarios. For this we do not need details of the specific model as long as neutrino decay can be described by a generic effective Lagrangian destroying one neutrino mass eigenstate and creating another one, i.e., $`_{\mathrm{int}}g_{ij}\overline{\nu }_j\nu _i`$. This also includes multi-particle decay into decay products other than neutrinos, where all momenta and energies of the undetected decay products are to be integrated over in the decay rates444For decay models involving more than one outgoing neutrino we could use a similar framework with slightly different properties. For Majoron-like models higher-order processes are assumed to be suppressed, because the coupling constants are considered to be small.. We will show that our operator framework satisfies the usual properties of the $`S`$-matrix approach and implements the correct phase relations. In certain cases, we will demonstrate this for Majoron decay.

### 3.1 Operators

Let us introduce decay and propagation operators in terms of creation and annihilation operators, $`\widehat{a}^{}`$ and $`\widehat{a}`$, respectively, as well as calculation rules for the transition probabilities. The symbols introduced in the operators will be explained thereafter. In the next subsection, we will show how the operators satisfy the expected properties. We define three operators:

###### Definition 1 (Disappearance operator)

$`𝒟_{}`$ is the transition operator generally known as the “decay operator”. Effectively, it returns the amplitude for an undecayed state remaining undecayed after traveling a distance $`l`$ along the baseline $`L`$:

$$𝒟_{}(l)=\underset{i}{}\mathrm{exp}\left(\frac{\alpha _il}{2E_i}\right)\widehat{a}_i^{}\widehat{a}_i.$$
(22)
###### Definition 2 (Appearance operator)

$`𝒟_+`$ is the differential transition operator, which destroys an in state and creates an out state in $`[l,l+dl]`$ along the baseline $`L`$, i.e., a new state “appears”:

$$𝒟_+(l,L)=\underset{ij}{\underset{i}{}\underset{j}{}}\sqrt{\frac{\alpha _{ij}}{E_i}}\sqrt{\eta _{ij}(L,l,D)}\widehat{a}_j^{}\widehat{a}_ie^{i\xi }.$$
(23)
The phase $`\xi `$ is a random phase taking into account the phase shift by additional (not measured) particles produced in the decay, such as Majorons (cf., Property 7 in Subsec. 3.2). The probability density, which is the square of the amplitude, will have to be integrated over $`l`$.

###### Definition 3 (Propagation operator)

$``$ is the operator propagating a state a distance $`l`$ along the baseline $`L`$555Here we are using the non-relativistic operator, which is different from the relativistic one only by an overall phase factor. This phase factor will cancel in the calculation of transition probabilities. For justification of this approximation see Ref..:

$$(l)=\underset{i}{}\mathrm{exp}\left(iE_il\right)\widehat{a}_i^{}\widehat{a}_i.$$
(24)
In these definitions, $`\eta _{ij}(L,l,D)`$ is the geometrical function introduced in Eq. (15) and defined in Eq. (16). It describes the fraction of decay products which will still arrive at the detector, i.e., is not redirected from the detector by decay. The decay rate is defined as $`\alpha _{ij}\frac{m_i}{\tau _{ij}^0}`$ for $`ij`$ decay, and $`\tau _{ij}^0`$ refers to the rest frame lifetime for that decay channel. Time dilation by the factor $`\gamma _i=E_i/m_i`$ implies that

$$\mathrm{\Gamma }_{ij}^{\mathrm{observer}}=\mathrm{\Gamma }_{ij}^0\gamma _i^1=\frac{1}{\tau _{ij}^0}\frac{m_i}{E_i}=\frac{\alpha _{ij}}{E_i}.$$
(25)
The rate $`\alpha _i_j\alpha _{ij}`$ is the overall decay rate of the state $`i`$, so that $`B\frac{\alpha _{ij}}{\alpha _i}`$ is the branching ratio for $`ij`$ decay. The factor of two in the exponent of $`𝒟_{}`$ and the square root in $`𝒟_+`$ comes from the fact that, besides some phase factors, amplitudes effectively behave like square roots of particles. However, as we will see in the next subsection, the above operators together with the definitions of the transition probabilities also implement the correct phase relations. Let us now define how to calculate the transition probabilities. This is to be done by successive application of decay and propagation operators:

###### Definition 4 (Calculation of transition probabilities)

The transition probability $`P(\nu _\alpha \nu _\beta )^n=P_{\alpha \beta }^n`$ by exactly $`n`$ intermediate decays is for $`n=0`$ given by

$$P_{\alpha \beta }^0=\left|\nu _\beta |(L)𝒟_{}(L)|\nu _\alpha \right|^2,$$
(26)
and for $`n>0`$ given by

$`P_{\alpha \beta }^n`$ $`=`$ $`{\displaystyle \underset{l_1=0}{\overset{L}{}}}\mathrm{}{\displaystyle \underset{l_n=l_{n1}}{\overset{L}{}}}|\nu _\beta |(Ll)𝒟_{}(Ll)`$ (27)
$`\times `$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}\left\{𝒟_+(l_i,L)(l_i)𝒟_{}(l_i)\right\}|\nu _\alpha |^2{\displaystyle }_{i=1}^ndl_i.`$
Here $`l_{i=1}^nl_i`$. The order of $`𝒟_{}`$ and $``$ is arbitrary, since these operators commute as we will see below, in Eq. (32).

In many cases, the total transition probability can be calculated as

$$P_{\alpha \beta }=\underset{n=0}{\overset{N_{\mathrm{max}}}{}}P_{\alpha \beta }^n,$$
(28)
where $`N_{\mathrm{max}}`$ denotes the maximal number of decays possible, depending on the decay model. The individual terms in Eq. (28) add up only if the particles produced by a different number of decays cannot be distinguished, such as by energy resolution, energy threshold or signature (flavor or spin). Depending on the problem, the terms may have to be split appropriately in general666For instance, we may even have two different $`P_{\alpha \beta }^2`$ coming from $`\nu _i\overline{\nu }_j\nu _k`$ and $`\nu _i\nu _j\nu _k`$ described by the appropriate decay rates. Another example is energy binning, where one $`P_{\alpha \beta }^i`$ may have to be split into fractions of different energy ranges to be calculated from the differential decay rate, such as Eq. (6) for Majoron decay..

Finally, we can introduce two simplified definitions for invisible as well as stable or almost stable visible decay products:

###### Definition 5 (Invisible decay)

For invisible decay products we can, in principle, only observe the undecayed particles themselves. Hence, we do not admit appearance operators, i.e.,

$$P_{\alpha \beta }^{\mathrm{invisible}}=P_{\alpha \beta }^0.$$
(29)
###### Definition 6 (Approximations for visible decay)

Considering maximal
one decay or lifetimes longer than $`\frac{Lm}{E}`$, we can neglect repeated decay terms ($`n>1`$) and approximate the total transition probability in Eq. (28) by

$$P_{\alpha \beta }P_{\alpha \beta }^0+P_{\alpha \beta }^1=P_{\alpha \beta }^{\mathrm{invisible}}+P_{\alpha \beta }^{\mathrm{appearance}},$$
(30)
where $`P_{\alpha \beta }^{\mathrm{appearance}}`$ describes the “appearance” of new particles by decay. From Eq. (27) we can read off the first correction term $`P_{\alpha \beta }^1`$ for visible decay in addition to $`P_{\alpha \beta }^{\mathrm{invisible}}`$:

$$P_{\alpha \beta }^1=\underset{l=0}{\overset{L}{}}\left|\nu _\beta |(Ll)𝒟_{}(Ll)𝒟_+(l,L)(l)𝒟_{}(l)|\nu _\alpha \right|^2𝑑l.$$
(31)
### 3.2 Properties

In this subsection, we want to make the most important properties of our operator framework more transparent. Since we deal with decay of states instead of decay of particles, we need to show that this describes particle decay correctly while preserving interference properties for states.

###### Property 1 (Creation and annihilation operators)

The effective creation and annihilation operators, $`\widehat{a}^{}`$ and $`\widehat{a}`$, respectively, are consistent with the $`S`$-matrix approach.

The operator $`\widehat{N}_i\widehat{a}_i^{}\widehat{a}_i`$ in Eqs. (22) and (24) is as usual the occupation number operator, whereas the operator $`\widehat{T}_{ij}\widehat{a}_j^{}\widehat{a}_i`$ in Eq. (23) is an effective transition operator coming from inserting the fermionic neutrino field expansions into an equation like Eq. (21). Here $`i`$ and $`j`$ refer to all quantum numbers, such as spin and momentum. Since we assume that conceptual issues, such as detection in a different energy bin or decay into an antiparticle, are dealt with by splitting Eq. (28) appropriately (cf., Def. 4), we do not need to take care of different quantum numbers here.

###### Property 2 (Commutation relations)

The operators $`𝒟_{}`$ and $``$ satisfy the commutation relation

$$[𝒟_{},]=0,$$
(32)
whereas in general

$$[𝒟_{},𝒟_+]0,[𝒟_+,]0.$$
(33)
These relations can be shown by direct calculation of the commutators, using the fermionic anticommutation relations $`\{\widehat{a}_i,\widehat{a}_j^{}\}=\delta _{ij}`$.

One can also derive Eq. (32) by introducing the notion of a complex mass square, i.e.,

$$\stackrel{~}{m}_i^2m_i^2i\alpha _i.$$

For $`\stackrel{~}{}(l)=_i\mathrm{exp}\left(i\stackrel{~}{E}_il\right)\widehat{a}_i^{}\widehat{a}_i`$, where $`\stackrel{~}{E}_ip+\frac{\stackrel{~}{m}_i^2}{2E_i}`$, we can read off from Eqs. (22) and (24) that

$$𝒟_{}(l)(l)=\stackrel{~}{}(l)=(l)𝒟_{}(l).$$

This means that the disappearance and propagation operators can be combined by the introduction of complex mass squares. Therefore, the commutator reduces to the commutator of two scalars, i.e., the exponentials of the real and imaginary parts of the mass square.

###### Property 3 (Decay into different channels)

Our description for the decay of states corresponds to a description for the decay of particles, which is given by the differential equation system

$$\frac{d}{dt}N_i=\underset{j,ji}{}\mathrm{\Gamma }_{ij}N_i+\underset{j,ji}{}\mathrm{\Gamma }_{ji}N_j.$$
(34)
Here $`\mathrm{\Gamma }_{ij}`$ is the rest frame decay rate for the particle decay $`ij`$ and $`N_i`$ is the number of particles $`i`$. From this we obtain for the disappearance of a particle $`N_{i,0}=1`$, $`N_{j,0}=0`$ for all $`ji`$, and $`\mathrm{\Gamma }_i_{j,ji}\mathrm{\Gamma }_{ij}`$

$$N_i(t)=e^{\mathrm{\Gamma }_it},$$
(35)
and for the appearance of a stable decay product $`N_{i,0}=0`$, $`N_{j,0}=1`$ for all $`ji`$, and $`\mathrm{\Gamma }_i=0`$

$$N_i(t)=\underset{j,ji}{}\frac{\mathrm{\Gamma }_{ji}}{\mathrm{\Gamma }_j}\left(1e^{\mathrm{\Gamma }_jt}\right).$$
(36)
Neglecting propagation operators777In this case, we are not interested in neutrino oscillation probabilities and phase factors will cancel anyway when squaring the probability amplitudes., we obtain from Eq. (26) for the disappearance probability of a mass eigenstate

$$P_{ii}=\left|\nu _i|𝒟_{}(l)|\nu _i\right|^2=e^{\frac{\alpha _i}{E_i}L},$$
(37)
which is identical to the disappearance probability in Eq. (35) if $`\mathrm{\Gamma }_i=\frac{\alpha _i}{E_i}`$, which we have in the observer’s rest frame.

Neglecting propagation operators, we obtain from Eq. (31) for the appearance of a stable decay product ($`\alpha _j=0`$) for exactly one decay

$$P_{ij}=_{l=0}^L\left|\nu _j|𝒟_+(l,L)𝒟_{}(l)|\nu _i\right|^2𝑑l=\eta \frac{\alpha _{ij}}{\alpha _i}\left(1e^{\frac{\alpha _i}{E_i}L}\right),$$
(38)
which is identical to the appearance probability in Eq. (36) if $`\eta =1`$ (everything detected) and $`\mathrm{\Gamma }_i=\frac{\alpha _i}{E_i}`$, which we again have in the observer’s rest frame.

In general, Eqs. (27) and (34) obey the same structure, since transitions among states are described by differential rates as in a Markov chain model.

###### Property 4 (Coherent and incoherent summation of amplitudes)

In quantum field theory, coherent and incoherent summation are often distinguished in the calculation of Feynman diagrams . The notion coherent refers to the summation of amplitudes

$$d\mathrm{\Gamma }\left|\underset{i}{}A_i\right|^2$$

and incoherent to the summation of squares of amplitudes, i.e., summation of probabilities

$$d\mathrm{\Gamma }\underset{i}{}\left|A_i\right|^2.$$

Basically, coherent summation applies to processes which are, in principle, indistinguishable and incoherent summation to processes which are, in principle, distinguishable. Our formalism implements the proper application of coherent or incoherent summation.

To show this, we split the statement into more specific parts in what follows.

###### Property 5 (Incoherent summation and repeated decay)

Repeated decay can be described by the incoherent summation in Eq. (28).

The differential decay rate in terms of the $`S`$-matrix is

$$d\mathrm{\Gamma }(inout)\left|out\left|S\right|in\right|^2,$$
(39)
where $`S`$ can be expanded as usually as

$$S=\underset{n=0}{\overset{\mathrm{}}{}}\frac{i^n}{n!}\underset{\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{}\underset{\mathrm{}}{\overset{\mathrm{}}{}}d^4x_1\mathrm{}d^4x_nT\left\{_{\mathrm{int}}(x_1)\mathrm{}_{\mathrm{int}}(x_n)\right\}.$$
(40)
Thus, the in and out states select the corresponding $`S`$-matrix expansion terms, which will evaluate to non-zero results. For a different number of decays we obtain different decay products characterized by different quantum numbers. Therefore, different $`S`$-matrix expansion terms will be selected in Eq. (39). This means that Feynman diagrams of different orders cannot interfere, even if we do not distinguish them by the measured decay products (often only the neutrinos but no other particles, such as Majorons, are measured). Hence, decay implies the principal possibility to measure the decay products and, thus, in terms of quantum mechanics acts as a measurement. Since repeated decay corresponds to higher order Feynman diagrams in quantum field theory, incoherent summation properly describes it.

###### Property 6 (On-shell propagation in between decays)

It is shown in Ref. that Feynman propagators for virtual particles reduce to on-shell propagators for macroscopic distances888In this reference, this was explicitely shown for the creation and detection processes.. Thus, individual decays in repeated decay can be treated independently without interference. More precisely, Feynman diagrams of the same order do not produce interference terms, because they can (in principle) be distinguished due to the large separation of vertices.

Equation (28) together with Eq. (27) imply that decays are treated separately with decay rates calculated as for independent processes.

###### Property 7 (Phase coherence between in and out neutrinos)

Since at least a third unmeasured particle, e.g., the Majoron, is involved, the overall phase coherence between in and out neutrinos is destroyed.

The random phase $`\xi `$ in Eq. (23) takes care of this. Since it does not appear in $`P_{\alpha \beta }^0`$ and cancels in $`P_{\alpha \beta }^1`$, it is only relevant for repeated decays for $`n>1`$. It also ensures that terms from different individual decays do not produce interference terms, i.e., Feynman diagrams of the same order are summed incoherently in this case (cf., Property 6).

###### Property 8 (Decay model and causality)

The decay model used in Eqs. (26) and (27) properly implements causality.

The operator $`𝒟_{}`$ gives the probability amplitude for a state remaining undecayed until the position $`l`$. The operator $`𝒟_+`$ gives the conditional probability amplitude for decay in the interval $`[l,l+dl]`$ for an undecayed state. Successive alternating applications of $`𝒟_{}`$ and $`𝒟_+`$ give the (conditional) probability amplitude for the decision path. This probability needs to be integrated to obtain the average over all possible decay positions $`l`$ in a causal order, which is taken into account by the integration limits, i.e., a particle cannot decay before it even exists. We are integrating probability amplitudes for different decay positions incoherently, taking into account that undecayed particles cannot interfere with decay products.

###### Property 9 (Neutrino oscillations of decay products)

The appearance operator $`𝒟_+`$ gives the proper relative phase in a superposition of mass eigenstates as an out state. Kinematics induced by different masses does not result in destruction of coherence among the superimposed states such as in charged current weak interactions.

The decay rate is to first order of the $`S`$-matrix expansion as usually defined as

$$\mathrm{\Gamma }(inout)=\frac{1}{T}\left|out|d^4x_{\mathrm{int}}|in\right|^2.$$
(41)
Inserting a Majoron-like interaction Lagrangian (e.g., Eq. (1)) with the corresponding field expansions yields for the decay rate from one flavor $`|\nu _\alpha =_iU_{\alpha i}^{}|\nu _i`$ into another one $`|\nu _\beta =_jU_{\beta j}^{}|\nu _j`$

$$\mathrm{\Gamma }_{\alpha \beta }=\frac{1}{T}\left|\nu _\beta |d^4x_{\mathrm{int}}|\nu _\alpha \right|^2=\frac{1}{T}\left|\underset{i}{}\underset{j}{}U_{\alpha i}^{}U_{\beta j}F_{ij}\right|^2,$$
(42)
and from one mass eigenstate $`|\nu _i`$ into another one $`|\nu _j`$

$$\mathrm{\Gamma }_{ij}=\frac{1}{T}\left|\nu _j|d^4x_{\mathrm{int}}|\nu _i\right|^2=\frac{1}{T}\left|F_{ij}\right|^2.$$
(43)
Here the $`F_{ij}`$’s refer to transition functions depending on the properties of the states $`i`$ and $`j`$ as well as the Lagrangian and the used field expansions. Without loss of generality, for Dirac field expansions of fermions and a Klein–Gordon field expansion of the boson, we obtain for Lagrangians of the Majoron decay type

$`F_{ij}`$ $`=`$ $`\nu _j|{\displaystyle d^4x_{\mathrm{int}}}|\nu _i`$ (44)
$`=`$ $`N_{ij}{\displaystyle d^4xe^{ip_Jx}𝒜\left(\overline{u^j}(𝐩_2)e^{ip_2x}\right)g_{ij}\left(u^i(𝐩_1)e^{ip_1x}\right)},`$
where $`N_{ij}`$ are normalization factors depending on energy, mass, and volume and $`𝒜`$ and $``$ are some operators depending on the Lagrangian, e.g., combinations of charge conjugation and chiral operators. By introducing $`\varphi _{ij}`$ as the phase factor of $`F_{ij}=\left|F_{ij}\right|e^{i\varphi _{ij}}`$ determined by Eq. (44), we can read off the following relation from Eq. (43):

$$\sqrt{\mathrm{\Gamma }_{ij}}e^{i\varphi _{ij}}=\frac{1}{\sqrt{T}}\left|F_{ij}\right|e^{i\varphi _{ij}}=\frac{1}{\sqrt{T}}F_{ij}.$$
(45)
As we will show below, for relativistic neutrinos $`\varphi _{ij}\varphi `$, i.e., $`\varphi _{ij}`$ is independent of the indices $`i`$ and $`j`$.

Without loss of generality, in our framework the rest frame transition probability between two flavor eigenstates for $`\eta _{ij}=1`$ (i.e., everything is detected) is given by

$$P_{\alpha \beta }=\left|\nu _\beta |𝒟_+|\nu _\alpha \right|^2T=\left|\nu _\beta |\underset{i}{}\underset{j}{}\sqrt{\mathrm{\Gamma }_{ij}}\widehat{a}_j^{}\widehat{a}_i|\nu _\alpha \right|^2T.$$
(46)
This can be obtained from Eq. (27) for an infinitesimally short baseline ($`dl=T`$) by using $`𝒟_{}1`$ to zeroth order in the expansion of the exponentials as well as $`\mathrm{\Gamma }_{ij}=\frac{\alpha _{ij}}{E_i}`$. Applying Eq. (45) and the expansions of the flavor eigenstates, we finally obtain from Eq. (46)

$$\mathrm{\Gamma }_{\alpha \beta }=\frac{P_{\alpha \beta }}{T}=\left|\underset{i}{}\underset{j}{}U_{\alpha i}^{}U_{\beta j}\frac{1}{\sqrt{T}}|F_{ij}|\right|^2=\frac{1}{T}\left|\underset{i}{}\underset{j}{}U_{\alpha i}^{}U_{\beta j}|F_{ij}|\right|^2.$$
(47)
This is equivalent to Eq. (42) if $`\varphi _{ij}=\varphi `$ is independent of the indices $`i`$ and $`j`$. In this case, the appearance operator is identical to the expression derived from the $`S`$-matrix expansion and the different masses do not destroy phase coherence in a superposition of states.

It remains to be shown that $`\varphi _{ij}`$ is independent of the indices $`i`$ and $`j`$ for relativistic neutrinos. Since $`\varphi _{ij}=\mathrm{arg}F_{ij}`$, any real factor in $`F_{ij}`$ such as, for instance, a normalization factor, does not affect $`\varphi _{ij}`$. We conclude from Eq. (44) that for real coupling constants the $`i`$ and $`j`$ dependence of $`\varphi _{ij}`$ can only be influenced by the exponentials and the spinors. By spatial integration, the exponentials are transformed into $`\delta `$-distributions implying momentum conservation and, hence, cannot affect the phase. Without loss of generality, the spinors can be written in the Dirac representation as

$$u_{}(𝐩)=N_1\left(\begin{array}{c}1\\ 0\\ N_2p^3\\ N_2(p^1+ip^2)\end{array}\right),u_{}(𝐩)=N_1\left(\begin{array}{c}0\\ 1\\ N_2(p^1ip^2)\\ N_2p^3\end{array}\right),$$
(48)
where $`N_1`$ and $`N_2`$ are (real) normalization factors. For relativistic neutrinos we can assume the three momentum $`𝐩`$ to be identical for all mass eigenstates and absorb the mass dependence of the four-momentum in a correction of the energy $`E_ip+\frac{m_i^2}{2p}`$. Then, they do not affect $`\varphi _{ij}`$ in this approximation, because the phases of these spinors only depend on the common three-momentum. Hence, $`\varphi _{ij}`$ is independent of the indices $`i`$ and $`j`$.

## 4 Invisible decay

In this section, we will treat invisible decay with our operator framework, i.e., decay into particles, such as sterile decoupled neutrinos, which are in principle unobservable. In fact, this is the most often assumed neutrino decay scenario.

### 4.1 Transition probability for invisible decay

The transition probability for invisible decay is calculated according to Eqs. (26) and (29). Since the operator $`𝒟_{}`$ only gives real factors, the calculation for relativistic neutrinos is quite straightforward and similar to the one for the general neutrino oscillation formula. We use $`|\nu _\alpha =_iU_{\alpha i}^{}|\nu _i`$ and define the following quantities:

$`J_{ij}^{\alpha \beta }`$ $``$ $`U_{\alpha i}U_{\alpha j}^{}U_{\beta i}^{}U_{\beta j},`$ (49)
$`\mathrm{\Delta }_{ij}`$ $``$ $`{\displaystyle \frac{\mathrm{\Delta }m_{ij}^2L}{4E}}{\displaystyle \frac{\left(m_i^2m_j^2\right)L}{4E}},`$ (50)
$`\mathrm{\Gamma }_{ij}`$ $``$ $`{\displaystyle \frac{m_iL}{2\tau _{0,i}E_i}}+{\displaystyle \frac{m_jL}{2\tau _{0,j}E_j}}\left({\displaystyle \frac{m_i}{\tau _{0,i}}}+{\displaystyle \frac{m_j}{\tau _{0,j}}}\right){\displaystyle \frac{L}{2E}}\left(\alpha _i+\alpha _j\right){\displaystyle \frac{L}{2E}}.`$ (51)
Here $`E_iE`$ has been used in the decay terms, because for decay the zeroth order terms in $`E_ip+\frac{m_i^2}{2p}`$ do not cancel, and therefore, the higher order terms can be neglected999For neutrino oscillations the zeroth order terms cancel due to multiplication with their complex conjugates in the calculation of the transition probability.. For more details on this kind of approximations see Ref..

After some calculations, not to be shown here, we finally obtain

$`P_{\alpha \beta }^{\mathrm{invisible}}`$ $`=`$ $`\underset{P_{\mathrm{CP}\mathrm{conserving}}}{\underset{}{\underset{P_{\mathrm{pure}\mathrm{decay}}}{\underset{}{{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}\mathrm{}J_{ij}^{\alpha \beta }e^{\mathrm{\Gamma }_{ij}}}}4\underset{i>j}{{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}}\mathrm{}J_{ij}^{\alpha \beta }\mathrm{sin}^2\mathrm{\Delta }_{ij}e^{\mathrm{\Gamma }_{ij}}}}`$ (52)
$``$ $`\underset{P_{\mathrm{CP}\mathrm{violating}}}{\underset{}{2\underset{i>j}{{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}}\mathrm{}J_{ij}^{\alpha \beta }\mathrm{sin}2\mathrm{\Delta }_{ij}e^{\mathrm{\Gamma }_{ij}}}}.`$
In our definition, pure decay admits neutrino flavor mixing, but no neutrino oscillations ($`\mathrm{\Delta }_{ij}=0`$ for all indices $`i`$ and $`j`$).

### 4.2 Example: Application to atmospheric neutrinos

In Ref., neutrino decay was considered as an alternative to atmospheric neutrino oscillations101010In this paper, we use the standard parameterization for the active neutrino mixing matrix
$$U_{3\times 3}=\left(\begin{array}{ccc}c_{12}c_{13}& s_{12}c_{13}& s_{13}e^{i\delta }\\ s_{12}c_{23}c_{12}s_{23}s_{13}e^{i\delta }& c_{12}c_{23}s_{12}s_{23}s_{13}e^{i\delta }& s_{23}c_{13}\\ s_{12}s_{23}c_{12}c_{23}s_{13}e^{i\delta }& c_{12}s_{23}s_{12}c_{23}s_{13}e^{i\delta }& c_{23}c_{13}\end{array}\right),$$
where $`s_{ij}\mathrm{sin}\theta _{ij}`$ and $`c_{ij}\mathrm{cos}\theta _{ij}`$ and where $`\theta _{ij}`$ are the vacuum mixing angles. . Assuming $`\mathrm{\Delta }m_{21}^2=0`$, $`\theta _{12}=\theta _{13}=0`$, and the CP-phase $`\delta =0`$, we obtain from Eq. (52) for decay of $`\nu _2`$ into a sterile decoupled neutrino with the rate $`\alpha \alpha _2`$ the same transition probabilities as in Ref.: For pure decay ($`\mathrm{\Delta }m_{32}^2=0`$) we have

$$P_{\mu \mu }=\left(c_{23}^2e^{\frac{\alpha L}{2E}}+s_{23}^2\right)^2,$$
(53)
and for decay and oscillations

$$P_{\mu \mu }=s_{23}^4+c_{23}^4e^{\frac{\alpha L}{E}}+2c_{23}^2s_{23}^2\mathrm{cos}\left(\frac{\mathrm{\Delta }m_{32}^2L}{2E}\right)e^{\frac{\alpha L}{2E}}.$$
(54)
The authors of Ref. fitted Super-Kamiokande data in the case of decay and oscillations with negligible $`\mathrm{\Delta }m_{32}^2`$ finding $`\alpha =\frac{1}{63}\frac{\mathrm{GeV}}{\mathrm{km}}`$ and $`\mathrm{cos}^2\theta \mathrm{cos}^2\theta _{23}=0.30`$. Pure decay without neutrino flavor mixing, i.e., $`\mathrm{cos}^2\theta =1`$, cannot describe the atmospheric neutrino data . Figure 4 shows the survival probability $`P_{\mu \mu }`$ as a function of the decay rate $`\alpha `$ and the sensitivity $`L/E`$. One can easily see that the transition from pure oscillation to pure decay is, in principle, determined by the curve $`\frac{\alpha L}{E}=𝒪(1)`$.

## 5 Visible decay

In this section, we will assume that the neutrino decay products are, in principle, observable. We know that we may either detect them separately from the undecayed neutrinos by different signatures or energies, or indistinguishably to the undecayed particles with the same signature in the same energy bin111111This strongly depends on the detector properties. For example, for solar neutrino decay into antineutrinos with detection in Borexino or Super-Kamiokande it was discussed in Ref... We postulated in Def. 4 that these conceptual differences manifest themselves in a problem-dependent splitting of Eq. (28) into appropriate parts.

### 5.1 Transition probability for visible decay

We use Eqs. (30) and (31) to first approximation for $`P_{\alpha \beta }^{\mathrm{appearance}}`$, i.e., the case of maximal one decay or sufficiently small decay rates. For the calculation of the transition probability we make the following assumptions:

* Stable decay products or $`\alpha _i\frac{E_i}{L}`$ for all $`i`$.
* The fraction $`\eta _{ij}(L,l,D)`$ of decayed particles not escaping detection by a change of direction does not depend on the decay position $`l`$, i.e., $`\eta _{ij}(L,l,D)=\eta _{ij}(L,D)\eta _{ij}`$ is a (real) constant in $`l`$. In fact, we observed that $`\eta _{ij}(L,l,D)1`$ in most cases of Majoron decay (cf., Subsec. 2.2).
* The masses and decay rates are non-degenerate121212The calculation of integrals changes for degenerate values of the $`\alpha _i`$’s and $`m_i`$’s (see calculation in Appendix A). Nevertheless, we will calculate limits with degenerate values later on..

For $`\mathrm{\Gamma }_{ij}`$ and $`\mathrm{\Delta }_{ij}`$ we use the definitions in Eqs. (50) and (51). Applying Eqs. (30) and (31) and assuming that we cannot distinguish decay products from undecayed particles, i.e., we add $`P_{\alpha \beta }^0`$ and $`P_{\alpha \beta }^1`$, we have for $`P_{\alpha \beta }`$

$`P_{\alpha \beta }`$ $``$ $`P_{\alpha \beta }^0+P_{\alpha \beta }^1=P_{\alpha \beta }^{\mathrm{invisible}}+P_{\alpha \beta }^{\mathrm{appearance}}`$ (55)
$`=`$ $`P_{\alpha \beta }^0+{\displaystyle \underset{l=0}{\overset{L}{}}}\left|\nu _\beta |(Ll)𝒟_{}(Ll)𝒟_+(l,L)(l)𝒟_{}(l)|\nu _\alpha \right|^2𝑑l.`$
#### The invisible decay term.

$`P_{\alpha \beta }^{\mathrm{invisible}}`$ is given by Eq. (52), i.e.,

$`P_{\alpha \beta }^{\mathrm{invisible}}`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}\mathrm{}J_{ij}^{\alpha \beta }e^{\mathrm{\Gamma }_{ij}}`$ (56)
$``$ $`4\underset{i>j}{{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}}\mathrm{}J_{ij}^{\alpha \beta }\mathrm{sin}^2\mathrm{\Delta }_{ij}e^{\mathrm{\Gamma }_{ij}}2\underset{i>j}{{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}}\mathrm{}J_{ij}^{\alpha \beta }\mathrm{sin}2\mathrm{\Delta }_{ij}e^{\mathrm{\Gamma }_{ij}}.`$
#### The appearance term.

$`P_{\alpha \beta }^{\mathrm{appearance}}`$ evaluates after some algebra (cf., Appendix A) to

$`P_{\alpha \beta }^{\mathrm{appearance}}`$ $`=`$ $`\underset{ij}{{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}}\underset{kl}{{\displaystyle \underset{k}{}}{\displaystyle \underset{l}{}}}\sqrt{\eta _{ij}\eta _{kl}}{\displaystyle \frac{L}{E}}{\displaystyle \frac{\sqrt{\alpha _{ij}\alpha _{kl}}}{\left(\mathrm{\Gamma }_{jl}\mathrm{\Gamma }_{ik}\right)^2+4\left(\mathrm{\Delta }_{ij}+\mathrm{\Delta }_{lk}\right)^2}}`$ (57)
$`\times `$ $`\{\mathrm{}(K_{ijkl}^{\alpha \beta })[(\mathrm{\Gamma }_{jl}\mathrm{\Gamma }_{ik})(e^{\mathrm{\Gamma }_{ik}}\mathrm{cos}(2\mathrm{\Delta }_{ki})e^{\mathrm{\Gamma }_{jl}}\mathrm{cos}(2\mathrm{\Delta }_{lj}))`$
$``$ $`2(\mathrm{\Delta }_{ij}+\mathrm{\Delta }_{lk})(e^{\mathrm{\Gamma }_{ik}}\mathrm{sin}(2\mathrm{\Delta }_{ki})e^{\mathrm{\Gamma }_{jl}}\mathrm{sin}(2\mathrm{\Delta }_{lj}))]`$
$``$ $`\mathrm{}(K_{ijkl}^{\alpha \beta })[(\mathrm{\Gamma }_{jl}\mathrm{\Gamma }_{ik})(e^{\mathrm{\Gamma }_{ik}}\mathrm{sin}(2\mathrm{\Delta }_{ki})e^{\mathrm{\Gamma }_{jl}}\mathrm{sin}(2\mathrm{\Delta }_{lj}))`$
$`+`$ $`2(\mathrm{\Delta }_{ij}+\mathrm{\Delta }_{lk})(e^{\mathrm{\Gamma }_{ik}}\mathrm{cos}(2\mathrm{\Delta }_{ki})e^{\mathrm{\Gamma }_{jl}}\mathrm{cos}(2\mathrm{\Delta }_{lj}))]\},`$
where $`K_{ijkl}^{\alpha \beta }U_{\alpha i}^{}U_{\beta j}U_{\alpha k}U_{\beta l}^{}`$ is a generalization of $`J_{ij}^{\alpha \beta }`$.

### 5.2 Limiting cases

Let us now consider some special cases for the general expression for the appearance term in Eq. (57).

#### No oscillations.

If we ignore neutrino oscillations, i.e., for degenerate masses, destruction of coherence in decay, or extremely small $`\mathrm{\Delta }m^2`$’s, we will obtain in the limit $`\mathrm{\Delta }m^20`$

$$P_{\alpha \beta }^{\mathrm{appearance}}\underset{ij}{\underset{i}{}\underset{j}{}}\underset{kl}{\underset{k}{}\underset{l}{}}\sqrt{\eta _{ij}\eta _{kl}}\frac{2\sqrt{\alpha _{ij}\alpha _{kl}}}{\alpha _j+\alpha _l\alpha _i\alpha _k}\mathrm{}(K_{ijkl}^{\alpha \beta })\left(e^{\mathrm{\Gamma }_{ik}}e^{\mathrm{\Gamma }_{jl}}\right).$$
(58)
#### No interference.

For random phases within superpositions of mass eigenstates, i.e.,

$$|\nu _\alpha =a|\nu _1+be^{i\xi _1}|\nu _2+ce^{i\xi _2}|\nu _3$$
(59)
with $`a`$, $`b`$, and $`c`$ some coefficients and $`\xi _1`$ and $`\xi _2`$ random phases, the interference terms in Eq. (57) vanish. We can include this by introducing $`\delta _{ik}\delta _{jl}`$ in the sum of Eq. (57) and we finally obtain

$$P_{\alpha \beta }^{\mathrm{appearance}}\underset{ij}{\underset{i}{}\underset{j}{}}\eta _{ij}\frac{\alpha _{ij}}{\alpha _j\alpha _i}\mathrm{}(K_{ijij}^{\alpha \beta })\left(e^{\frac{\alpha _iL}{E}}e^{\frac{\alpha _jL}{E}}\right).$$
(60)
#### Only one possible decay channel.

Assume that there is only one possible decay channel, i.e., $`\alpha _{ij}\alpha _i\alpha 0`$ and all other $`\alpha _{kl}=0`$. Then, we obtain the simplified formula

$$P_{\alpha \beta }^{\mathrm{appearance}}\eta _{ij}\mathrm{}(K_{ijij}^{\alpha \beta })\left(1e^{\frac{\alpha L}{E}}\right).$$
(61)
### 5.3 Example: Application to atmospheric neutrinos

Here we will extend the example from Subsec. 4.2, where Eq. (54) gives the oscillation formula for decay of $`\nu _2`$ into an invisible decay product such as a sterile neutrino. Let us now assume that $`\nu _2`$ decays into the active neutrino $`\nu _3`$. According to Eq.(55), we are now looking for $`P_{\mu \mu }^{\mathrm{appearance}}`$ in addition to $`P_{\mu \mu }^{\mathrm{invisible}}`$, which can be calculated from Eq. (61) in this limit. If we assume that both undecayed and decayed neutrinos are detected indistinguishably, we will obtain

$$P_{\mu \mu }=P_{\mu \mu }^{\mathrm{invisible}}+P_{\mu \mu }^{\mathrm{appearance}}$$
(62)
with $`P_{\mu \mu }^{\mathrm{invisible}}`$ in Eq. (54) and

$$P_{\mu \mu }^{\mathrm{appearance}}=\eta _{23}s_{23}^2c_{23}^2\left(1e^{\frac{\alpha L}{E}}\right).$$
(63)
Here $`\alpha `$ is the decay rate of $`\nu _2\nu _3`$ decay and $`E`$ is the energy of the undecayed neutrino. Figure 5 shows the different parts of the survival probability for the parameter values given in the figure caption. One can easily see the effect of the increase of the survival probability for $`\frac{L}{E}\frac{1}{\alpha }`$.

## 6 An example for neutrino oscillations of decay products

In this subsection, we will construct a minimal example for observing neutrino oscillations among decay products.

### 6.1 The scenario

Since we want to create a superposition of out states in decay, we have to assume that at least two of the $`g_{ij}`$’s, i.e., $`\alpha _{ij}`$’s, in the interaction Lagrangian are non-zero. We define a mass hierarchy $`m_3>m_2>m_1`$ and only allow the decay chain $`\nu _3\nu _2\nu _1`$, i.e., $`\alpha _{32}\alpha _30`$, $`\alpha _{21}\alpha _20`$, and all the other $`\alpha _{ij}`$’s are equal to zero. In order to demonstrate the physical effects we choose suitable parameters. Hence, we define the $`\mathrm{\Delta }m^2`$’s quite close to each other, i.e., $`\mathrm{\Delta }m_{21}^2=510^4\mathrm{eV}^2`$, $`\mathrm{\Delta }m_{32}^2=3.310^3\mathrm{eV}^2`$, maximal mixing, i.e., $`|\nu _e=\frac{1}{\sqrt{3}}|\nu _1+\frac{1}{\sqrt{3}}|\nu _2+\frac{1}{\sqrt{3}}|\nu _3`$, and the $`\alpha `$’s close to each other, i.e., $`\alpha _2=\frac{1}{1000}\frac{\mathrm{GeV}}{\mathrm{km}}`$ and $`\alpha _3=\frac{1}{1100}\frac{\mathrm{GeV}}{\mathrm{km}}`$. Furthermore, we assume that $`\eta _{ij}1`$ for all $`i`$, $`j`$ and the CP-phase $`\delta =0`$.

### 6.2 Transition probabilities

Let us now look at the disappearance probability of electron neutrinos $`P_{ee}`$. From Eq. (52) we obtain the invisible decay disappearance probability as

$`P_{ee}^{\mathrm{invisible}}`$ $`=`$ $`J_{11}^{ee}+J_{22}^{ee}e^{\alpha _2\frac{L}{E}}+J_{33}^{ee}e^{\alpha _3\frac{L}{E}}`$ (64)
$`+`$ $`2J_{12}^{ee}e^{\alpha _2\frac{L}{E}}\left(12\mathrm{sin}^2\left({\displaystyle \frac{\mathrm{\Delta }m_{21}^2}{4}}{\displaystyle \frac{L}{E}}\right)\right)`$
$`+`$ $`2J_{13}^{ee}e^{\alpha _3\frac{L}{E}}\left(12\mathrm{sin}^2\left({\displaystyle \frac{\mathrm{\Delta }m_{31}^2}{4}}{\displaystyle \frac{L}{E}}\right)\right)`$
$`+`$ $`2J_{23}^{ee}e^{\frac{\alpha _2+\alpha _3}{2}\frac{L}{E}}\left(12\mathrm{sin}^2\left({\displaystyle \frac{\mathrm{\Delta }m_{32}^2}{4}}{\displaystyle \frac{L}{E}}\right)\right).`$
The appearance term is obtained from Eq. (57). Since in this case $`K_{ijkl}^{ee}=U_{ei}^{}U_{ej}U_{ek}U_{el}^{}=\frac{1}{9}`$ and $`\mathrm{}K=0`$, the terms in the sum of Eq. (57) are equal for simultaneous $`ik`$ and $`jl`$ index exchanges. In order to make the result physically more transparent, we recognize that only a limited number of the terms in the sum are non-vanishing and split the formula into three parts:

$$P_{ee}^{\mathrm{appearance}}=P_{ee}^{\mathrm{app},1}+P_{ee}^{\mathrm{app},2}+P_{ee}^{\mathrm{int}}.$$
(65)
Here $`P_{ee}^{\mathrm{app},1}`$ refers to the production of new $`\nu _1`$’s by decay,

$$P_{ee}^{\mathrm{app},1}=K_{2121}^{ee}\left(1e^{\alpha _2\frac{L}{E}}\right),$$
(66)
$`P_{ee}^{\mathrm{app},2}`$ to the production of new $`\nu _2`$’s by decay,

$$P_{ee}^{\mathrm{app},2}=K_{3232}^{ee}\frac{\alpha _3}{\alpha _2\alpha _3}\left(e^{\alpha _3\frac{L}{E}}e^{\alpha _2\frac{L}{E}}\right),$$
(67)
and $`P_{ee}^{\mathrm{int}}`$ to an interference term describing neutrino oscillations among the decay products $`\nu _1`$ and $`\nu _2`$ (here $`K_{3221}^{ee}=K_{2132}^{ee}`$ is used),

$`P_{ee}^{\mathrm{int}}`$ $`=`$ $`4K_{3221}^{ee}{\displaystyle \frac{\sqrt{\alpha _3\alpha _2}}{\alpha _3^2+\left(\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{32}^2\right)^2}}`$
$`\times `$ $`\{\alpha _3[e^{\frac{\alpha _2}{2}\frac{L}{E}}\mathrm{cos}\left({\displaystyle \frac{\mathrm{\Delta }m_{21}^2}{2}}{\displaystyle \frac{L}{E}}\right)e^{\frac{\alpha _2+\alpha _3}{2}\frac{L}{E}}\mathrm{cos}\left({\displaystyle \frac{\mathrm{\Delta }m_{32}^2}{2}}{\displaystyle \frac{L}{E}}\right)]`$
$`+`$ $`(\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{32}^2)[e^{\frac{\alpha _2}{2}\frac{L}{E}}\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Delta }m_{21}^2}{2}}{\displaystyle \frac{L}{E}}\right)e^{\frac{\alpha _2+\alpha _3}{2}\frac{L}{E}}\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Delta }m_{32}^2}{2}}{\displaystyle \frac{L}{E}}\right)]\}.`$
### 6.3 Results and analysis

Figure 6 shows the different parts of the probability $`P_{ee}^{\mathrm{appearance}}`$ for the scenario data. The appearance of $`\nu _2`$ as decay product is described by $`P_{ee}^{\mathrm{app},2}`$. It grows in the beginning, because it is dominated by $`\nu _3\nu _2`$ decay. Since $`\alpha _2>\alpha _3`$, the $`\nu _2`$ is decaying faster than produced. Therefore, $`P_{ee}^{\mathrm{app},2}`$ falls again for large $`\frac{L}{E}`$. The appearance of $`\nu _1`$ is determined by $`P_{ee}^{\mathrm{app},1}`$. Since $`\nu _1`$ is defined to be stable, $`P_{ee}^{\mathrm{app},1}`$ grows until it reaches its asymptotic value $`lim_x\mathrm{}P_{ee}^{\mathrm{app},1}=K_{2121}^{ee}=\frac{1}{9}`$.

Describing neutrino oscillations of decay products, i.e., $`\nu _1`$ and $`\nu _2`$, the interference term $`P_{ee}^{\mathrm{int}}`$ basically oscillates as a beat with the two frequencies determined by $`\mathrm{\Delta }m_{32}^2\pm \mathrm{\Delta }m_{21}^2`$. The $`\mathrm{\Delta }m_{32}^2`$-dependency comes from the phase propagation before decay and the $`\mathrm{\Delta }m_{21}^2`$-dependency from the phase propagation after decay.

Figure 7 shows the different parts of the total survival probability $`P_{ee}`$. The invisible decay probability (disappearance) is the ordinary oscillation probability damped by decay. Note that $`P_{ee}^{\mathrm{appearance}}`$ and $`P_{ee}^{\mathrm{invisible}}`$ can be only sensibly added if the decay products cannot be distinguished from the original ones in the detection process. In this case, the detection rates will be enhanced by a fixed amount at large $`L/E`$. Neutrino oscillations of undecayed particles will also vanish in this region. For small $`L/E`$ there is a correction due to $`\nu _2`$ abundance by decay and in addition due to the interference term.

Neutrino oscillations of decay products, induced by the oscillations of the interference term, can become unobservable due to suppression by the mixing angles or the decay rates different by some orders of magnitude. The first point is quite obvious, since the interference term is in general proportional to

$$K_{ijkl}^{\alpha \beta }=U_{\alpha i}^{}U_{\beta j}U_{\alpha k}U_{\beta l}^{},\text{where}ij,kl,(i,j)(k,l).$$
(69)
This requires that there are coefficients of at least three different states involved as a product, i.e., all of them need to be quite large. The second point is physically obvious, but harder to see in our equations. Since the interference describing neutrino oscillations of decay products can only be observed if two states are produced simultaneously in decay, the couplings to both of them need to be quite large. Large couplings imply large decay rates for both transitions. If one of them is too large, the corresponding state will have decayed before oscillations can be observed. Hence, the decay rates need to be large enough and of the same order of magnitude in order to obtain reasonable results. A measure for the magnitude of the interference effects is the absolute maximum of $`|P_{ee}^{\mathrm{int}}|`$. In Fig. 8, we show the result of a numerical evaluation of this function for different values of the decay rates $`\alpha _3`$ and $`\alpha _2`$. One can easily verify that for too small or too different decay rates, this function gives small values, i.e., we do not see any interference effects.

### 6.4 Application to limiting cases

In Subsec. 5.2, we studied the appearance probability in certain limits. Let us now apply these to our scenario to see what happens if certain terms are assumed to be negligible.

#### No oscillations.

In this limit, all $`\mathrm{\Delta }m^2`$’s are equal to zero (cf., Eq. (58)). Here we can simply use

$$P_{ee}^{\mathrm{appearance}}=P_{ee}^{\mathrm{app},1}+P_{ee}^{\mathrm{app},2}+P_{ee}^{\mathrm{int}}|_{\mathrm{\Delta }m_{ij}^20}.$$
(70)
#### No interference.

This limit is described by Eq. (60), which is equivalent to neglecting the interference term. We obtain

$$P_{ee}^{\mathrm{appearance}}=P_{ee}^{\mathrm{app},1}+P_{ee}^{\mathrm{app},2}.$$
(71)
Figure 9 shows the different limits of the appearance term compared to the original one. In all cases, we obtain similar values for small and large $`L/E`$, since for $`L/E(\mathrm{\Delta }m^2)^1`$ the oscillation terms are ineffective and for $`L/E\alpha ^1`$ they are suppressed by decay. In the limit of no oscillations, there is no oscillating behavior of the curve. Nevertheless, there are interference terms coming from the coherent coupling of the Lagrangian to different mass eigenstates. Hence, in this case we obtain an enhancement in the appearance probability. In the limit of no interference, we obtain a curve similar to the original one, but without oscillations, which are described by the interference term.

## 7 Summary and conclusions

In this paper, we have introduced a combined treatment of neutrino decay and neutrino oscillations in a formal model-independent operator framework. We have shown that it can be easily applied to special cases such as atmospheric neutrinos. We have also seen that for Majoron decay active neutrinos as decay products have to be taken into account in the detection rates and may even be indistinguishable from the original undecayed particles by the detector. In addition, in Majoron-like decay models neutrino decay products may oscillate.

We conclude that the visibility of decay products in Majoron decay models may crucially affect the detection rates depending on the properties of the detector. For other decay models only a certain fraction of the decay products may arrive at the detector due to kinematics. Moreover, neutrino oscillations of decay products are observable for relatively large neutrino flavor mixing only.

One can show that for the combination of neutrino decay and neutrino oscillations, the dimension of the parameter space is equal to $`\frac{3}{2}n(n1)`$, where $`\frac{1}{2}n(n1)`$ parameters come from decay and $`n(n1)`$ from oscillations. Here $`n`$ denotes the number of neutrino flavors. Hence, the number of parameters is increased when neutrino decay is included. Thus, especially in multi-neutrino scenarios, effects not consistent with the conventional neutrino oscillation framework could be described by a combined neutrino decay and neutrino oscillation approach.

As a realistic scenario, we could, for example, use a four-neutrino framework with a decoupled sterile neutrino $`\nu _S`$ (cf., Eq. (2)) and a mass hierarchy $`m_3>m_2>m_1>m_S`$ (cf., Fig. 10). In such a scenario, we have up to six decay rates in addition to the usual mixing parameters. From other conditions such as the observability of supernova neutrinos, we then would need to constrain the lifetimes for individual $`\nu _i`$’s (e.g., $`\nu _1`$ almost stable). Since the maximal number of possible decays is $`N_{\mathrm{max}}=3`$ in this scenario, the exact transition probability $`P_{\alpha \beta }`$ without lifetime constraints would be determined by $`P_{\alpha \beta }^0`$, $`P_{\alpha \beta }^1`$, $`P_{\alpha \beta }^2`$, and $`P_{\alpha \beta }^3`$ in Eq. (28).

Since non-zero neutrino masses imply both neutrino decay and neutrino oscillations, neutrino decay may become interesting not only as an alternative to neutrino oscillations but also as a part of the whole neutrino scenario. Even small $`g_{ij}`$’s may affect detection rates of long traveling neutrinos such as supernova or primordial ones. Flux measurements (e.g., of supernova neutrinos) cannot always be used to extract neutrino mass limits reliably, because of different path lengths of undecayed and decayed neutrinos. Thus, dispersion may be mimiced by decay.

As a final conclusion, our operator formalism is a generalization for a combined treatment of neutrino decay and neutrino oscillations could be of interest to those, who are interested in analyzing neutrino data in a more general way. In addition, it may also be applied in other oscillation scenarios outside of neutrino physics, where coherence is not destroyed in decay.

We would like to thank Evgeny Akhmedov for useful discussions and Jörn Kersten for proof-reading the manuscript. This work was supported by the Swedish Foundation for International Cooperation in Research and Higher Education (STINT), the Wenner-Gren Foundations, the “Studienstiftung des deutschen Volkes” (German National Merit Foundation), and the “Sonderforschungsbereich 375 für Astro-Teilchenphysik der Deutschen Forschungsgemeinschaft”.

## Appendix A Calculation of the appearance term

We are using Eq. (31) as an approximation for the appearance term:

$`P_{\alpha \beta }^{\mathrm{appearance}}`$ $``$ $`{\displaystyle \underset{l=0}{\overset{L}{}}}\left|\underset{A_{\alpha \beta }}{\underset{}{\nu _\beta |(Ll)𝒟_{}(Ll)𝒟_+(l,L)(l)𝒟_{}(l)|\nu _\alpha }}\right|^2𝑑l=`$ (72)
$`=`$ $`{\displaystyle \underset{l=0}{\overset{L}{}}}\left|A_{\alpha \beta }\right|^2𝑑l,`$
where $`A_{\alpha \beta }`$ is the (differential) transition amplitude. Applying Eqs. (22), (23), (24), and $`|\nu _\alpha =_iU_{\alpha i}^{}|\nu _i`$, we obtain for the transition amplitude131313Here we already neglect the random phase factor $`e^{i\xi }`$ from the operator $`𝒟_+`$, since it will cancel anyway by taking its absolute value squared in the transition probability.

$$A_{\alpha \beta }=\underset{ij}{\underset{i}{}\underset{j}{}}U_{\alpha i}^{}U_{\beta j}e^{iE_j(Ll)}e^{\frac{\alpha _j(Ll)}{2E_j}}\sqrt{\frac{\alpha _{ij}}{E_i}}\sqrt{\eta _{ij}}e^{iE_il}e^{\frac{\alpha _il}{2E_i}}.$$
(73)
In this equation, we already have evaluated the creation and annihilation operators for readability. Now we use the approximation for relativistic neutrinos $`E_ip+\frac{m_i^2}{2p}p+\frac{m_i^2}{2E}`$. For neutrino oscillations the three-momentum $`p`$ only gives a common phase factor, which will cancel. For neutrino decay, the exponentials $`e^{\frac{\alpha _il}{2E_i}}`$ and $`e^{\frac{\alpha _j(Ll)}{2E_j}}`$ are real, and thus $`E_ipE`$ in the lowest order nontrivial approximation. Using the relativistic approximations as well as the definition $`\mathrm{\Delta }m_{ab}^2m_a^2m_b^2`$ then yields

$$A_{\alpha \beta }=\underset{ij}{\underset{i}{}\underset{j}{}}\underset{l\mathrm{independent}}{\underset{}{U_{\alpha i}^{}U_{\beta j}\sqrt{\eta _{ij}}\sqrt{\frac{\alpha _{ij}}{E}}e^{\frac{\alpha _jL}{2E}}e^{i\frac{m_j^2}{2E}L}}}\underset{l\mathrm{dependent}}{\underset{}{e^{\frac{(\alpha _j\alpha _i)l}{2E}}e^{i\frac{\mathrm{\Delta }m_{ji}^2}{2E}l}}}.$$
(74)
For the appearance probability we obtain

$`P_{\alpha \beta }^{\mathrm{appearance}}`$ $`=`$ $`{\displaystyle \underset{l=0}{\overset{L}{}}}\left|A_{\alpha \beta }\right|^2𝑑l={\displaystyle \underset{l=0}{\overset{L}{}}}A_{\alpha \beta }A_{\alpha \beta }^{}𝑑l=\underset{ij}{{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}}\underset{kl}{{\displaystyle \underset{k}{}}{\displaystyle \underset{l}{}}}\underset{K_{ijkl}^{\alpha \beta }}{\underset{}{U_{\alpha i}^{}U_{\beta j}U_{\alpha k}U_{\beta l}^{}}}`$
$`\times `$ $`\sqrt{\eta _{ij}\eta _{kl}}\sqrt{{\displaystyle \frac{\alpha _{ij}}{E}}}\sqrt{{\displaystyle \frac{\alpha _{kl}}{E}}}e^{\frac{\alpha _j+\alpha _l}{2E}L}e^{i\frac{\mathrm{\Delta }m_{jl}^2}{2E}L}{\displaystyle \underset{l=0}{\overset{L}{}}}e^{\frac{\alpha _j+\alpha _l\alpha _i\alpha _k}{2E}l}e^{i\frac{\mathrm{\Delta }m_{ji}^2\mathrm{\Delta }m_{lk}^2}{2E}l}𝑑l,`$
where the $`\eta _{ij}`$’s are assumed to be real. For further evaluation we use the integral ($`a`$, $`b`$ real)

$$\underset{0}{\overset{L}{}}e^{(a+bi)l}𝑑l=\frac{1}{a+bi}\left(e^{(a+bi)L}1\right)=\frac{abi}{a^2+b^2}\left(e^{(a+bi)L}1\right)$$
(76)
for $`a+bi0`$. Comparison with Eq. (LABEL:appp1) shows that the latter condition implies

$$\alpha _j+\alpha _l\alpha _i\alpha _k+i\left(m_j^2m_i^2m_l^2+m_k^2\right)0.$$
(77)
Since $`ijkl`$ by the summation rules, this can only happen for $`(i=k)(j=l)(\alpha _i=\alpha _j)`$ for non-degenerate $`m`$’s. Here we assume the $`\alpha `$’s to be non-degenerate for simplicity, so that we can apply Eq. (76). Using the same abbreviations as for invisible decay in Eqs. (50) and (51), we can continue evaluating Eq. (LABEL:appp1) with Eq. (76) by identifying $`a=\frac{\alpha _j+\alpha _l\alpha _i\alpha _k}{2E}=\mathrm{\Gamma }_{jl}\mathrm{\Gamma }_{ik}`$ and $`b=\frac{\mathrm{\Delta }m_{ji}^2\mathrm{\Delta }m_{lk}^2}{2E}=2\left(\mathrm{\Delta }_{ji}+\mathrm{\Delta }_{kl}\right)`$

$`P_{\alpha \beta }^{\mathrm{appearance}}`$ $`=`$ $`\underset{ij}{{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}}\underset{kl}{{\displaystyle \underset{k}{}}{\displaystyle \underset{l}{}}}K_{ijkl}^{\alpha \beta }\sqrt{\eta _{ij}\eta _{kl}}\sqrt{\alpha _{ij}\alpha _{kl}}{\displaystyle \frac{L}{E}}e^{\frac{\alpha _j+\alpha _l}{2E}L}e^{i\frac{\mathrm{\Delta }m_{jl}^2}{2E}L}`$ (78)
$`\times `$ $`{\displaystyle \frac{\mathrm{\Gamma }_{jl}\mathrm{\Gamma }_{ik}2i\left(\mathrm{\Delta }_{ji}+\mathrm{\Delta }_{kl}\right)}{\left(\mathrm{\Gamma }_{jl}\mathrm{\Gamma }_{ik}\right)^2+4\left(\mathrm{\Delta }_{ji}+\mathrm{\Delta }_{kl}\right)^2}}\left(e^{\left(\alpha _j+\alpha _l\alpha _i\alpha _k+i\left(\mathrm{\Delta }m_{ji}^2\mathrm{\Delta }m_{lk}^2\right)\right)\frac{L}{2E}}1\right)`$
$`=`$ $`\underset{ij}{{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}}\underset{kl}{{\displaystyle \underset{k}{}}{\displaystyle \underset{l}{}}}K_{ijkl}^{\alpha \beta }\sqrt{\eta _{ij}\eta _{kl}}\sqrt{\alpha _{ij}\alpha _{kl}}{\displaystyle \frac{L}{E}}{\displaystyle \frac{\mathrm{\Gamma }_{jl}\mathrm{\Gamma }_{ik}2i\left(\mathrm{\Delta }_{ji}+\mathrm{\Delta }_{kl}\right)}{\left(\mathrm{\Gamma }_{jl}\mathrm{\Gamma }_{ik}\right)^2+4\left(\mathrm{\Delta }_{ji}+\mathrm{\Delta }_{kl}\right)^2}}`$
$`\times `$ $`\left(e^{\mathrm{\Gamma }_{ik}}e^{2i\mathrm{\Delta }_{ki}}e^{\mathrm{\Gamma }_{jl}}e^{2i\mathrm{\Delta }_{lj}}\right).`$
Because $`P_{\alpha \beta }^{\mathrm{appearance}}`$ is real, we have $`P_{\alpha \beta }^{\mathrm{appearance}}=\mathrm{}P_{\alpha \beta }^{\mathrm{appearance}}`$. For $`\varphi `$ real, $`y`$ and $`z`$ complex, it can be shown that

$`\mathrm{}(ze^{i\varphi })`$ $`=`$ $`\mathrm{}(z)\mathrm{cos}\varphi \mathrm{}(z)\mathrm{sin}\varphi ,`$ (79)
$`\mathrm{}(iye^{i\varphi })`$ $`=`$ $`\mathrm{}(y)\mathrm{cos}\varphi \mathrm{}(y)\mathrm{sin}\varphi .`$ (80)
We can use this to find

$`\mathrm{}\left(K(c+di)e^{2i\mathrm{\Delta }}\right)`$ $`=`$ $`\mathrm{}(K)c\mathrm{cos}2\mathrm{\Delta }\mathrm{}(K)d\mathrm{sin}2\mathrm{\Delta }`$ (81)
$``$ $`\mathrm{}(K)c\mathrm{cos}2\mathrm{\Delta }\mathrm{}(K)d\mathrm{sin}2\mathrm{\Delta }`$
for $`c`$, $`d`$ real. In comparison with Eq. (78), we identify $`c=\mathrm{\Gamma }_{jl}\mathrm{\Gamma }_{ik}`$ and $`d=2(\mathrm{\Delta }_{ij}+\mathrm{\Delta }_{lk})`$. Applying this formula and re-grouping in $`\mathrm{}K`$ and $`\mathrm{}K`$ yields the final result in Eq. (57).