Constraining visible neutrino decay at KamLAND and JUNO

We study visible neutrino decay at the reactor neutrino experiments KamLAND and, JUNO. Assuming the Majoron model of neutrino decay, we obtain constraints on the couplings between Majoron and neutrino as well as on the lifetime/mass of the most massive neutrino state i.e., $\tau_{3} / m_{3}$ or $\tau_{2} / m_{2}$, respectively, for the normal or the inverted mass orderings. We obtain the constraints on the lifetime $\tau_{2} / m_{2} \geq 1.4 \times 10^{-9}~\rm{s/eV}$ in the inverted mass ordering for both KamLAND and JUNO at 90% CL. In the normal ordering in which the bound can be obtained for JUNO only, the constraint is milder than the inverted ordering case, $\tau_{3} / m_{3} \geq 1.0 \times 10^{-10}~\rm{s/eV}$ at 90% CL. We find that the dependence of lightest neutrino mass ($=m_{\rm{lightest}}$), $m_1 (m_3)$ for the normal (inverted) mass ordering, on the constraints for the different types of couplings (scalar or pseudo-scalar) is rather strong, but the $m_{\rm{lightest}}$ dependence on the lifetime/mass bound is only modest.


Introduction
In the Standard Model (SM) of particle physics, neutrinos are stable particles. This is not only true in the original formulation of the SM, in which neutrinos are massless but also true in practice in the neutrino mass embedded version, the νSM. In the latter, the neutrino has a finite lifetime due to a nonzero mass and the lepton flavor mixing. But, the lifetime is extremely long, > 10 45 s for radiative decay [1,2]. Since such a very long lifetime is practically unmeasurable, neutrinos can be regarded as stable particles in the νSM. Therefore, if neutrino decay is detected it will imply evidence for new physics beyond the SM.
It appears that most of the foregoing analyses of ν 3 lifetime were done under the assumption of invisible decay, namely, the case that decay products are unobservable. See, however, refs. [30,46,[52][53][54] for the analyses with visible neutrino decay. Moreover, the majority of the works devoted to the analyses of neutrino decay so far restrict themselves into the case of NO.
In this paper, we discuss the bound on neutrino lifetime with visible neutrino decay. We consider both mass orderings, NO and IO. To treat visible neutrino decay we must specify the model which allows neutrinos decay, and we use the Majoron model [55][56][57][58][59][60] as a concrete model of visible neutrino decay (see Section 3). To place the bound on neutrino decay lifetime, we analyze the reactor neutrino experiments, KamLAND [61] and JUNO [62]. For the former we use the real data in ref. [61], and for the latter the simulated one assuming the total number of 140,000 events which would be obtained with an exposure of 220 GW·years or somewhat more dependinig on the actual avilability (which is expected to be ∼ 85-90 %) of reactors. 1 Under the visible neutrino decay hypothesis, there appear a few new features in the analysis: • Unlike the case of invisible decay, the decay products include active neutrino states, which we call the "daughter" neutrinos, 2 and they can produce additional events in the detectors; 1 Because of this feature and for a very simplified code, our analysis may be called more properly as the one for the "JUNO-like" setting. 2 For notations of the "parent" and "daughter" neutrinos, see Section 4.1 for the definitions.
• There is a clear difference in the constraints we will obtain between the cases of NO and IO. In the IO, ν 1 and ν 2 decay into ν 3 and ν 3 , which leads to a significant deficit of inverse beta decay events due to the large ν e component in the parent ν 1 and ν 2 mass eigenstates. Whereas in the NO, ν 3 decays into ν 1 and ν 2 as well as ν 1 and ν 2 . Since the parent ν 3 states are much less populated by ν e due to the small value of θ 13 , the effect of decay on the ν e spectrum is only modest. Now, we must spell out our attitude on the astrophysical neutrino bound on neutrino decay. Though the bound is likely to be correct and is probably robust we do not use the lifetime bound as granted in our analysis. The reasons for doing this is twofold: (1) The lifetime bound from the reactor neutrino experiments is completely independent of the bounds obtained by the solar and the supernova data. (2) The analysis to derive the solar neutrino bounds on the Majoron couplings is not simple. Most notably, the antineutrino appearance from the sun is involved, which requires a separate analysis. In a variety of contexts, it does make sense to obtain the laboratory bounds even though the astrophysical bounds are much stronger than the laboratory ones. 3 In our analysis, for simplicity, we turn on only the Majoron couplings g 23 and g 13 . In principle we can turn on all the couplings including g 12 , but the analysis becomes far more complicated. It is also very likely that the qualitative features of the bound obtained for the Majoron couplings remain unchanged in our reduced setting. Therefore, only the following decay modes are allowed in our setting: In this paper, after understanding all the above points, we concentrate on deriving the reactor neutrino bound on the Majoron couplings g 23 and g 13 , and the corresponding τ /m in both the NO and the IO. Thus, we explore systematically for the first time, assuming visible neutrino decay, the constraints that can be imposed on ν 3 lifetime (in the case of NO) and on ν 2 and ν 1 lifetimes (in the case of IO) by using the medium-and long-baseline reactor anti-neutrinos experiments. Yet, we must mention that our analysis is based on the Majoron model, and is done under the assumption of switching off the coupling between ν 1 , ν 2 , and Majoron.

Brief recollection of the existing bounds on neutrino decay
In most of the existing literatures, the bounds on neutrino decay have been calculated for NO, and hence Table 1 contains the bound for the NO which uses the ν 3 /ν 3 decay mode only. The tabulated bounds in Table 1 span the region from a few ×10 −12 to a few × 10 −10 /eV. These bounds, which utilize the artificial neutrino beams, are very loose compared with the solar neutrino bounds [22][23][24][25][26][27][28][29][30]. The latter which is usually quote as the one for ν 2 is: τ 2 /m 2 > ∼ 7.02 × 10 −4 s/eV at 99% C.L. [28].

Phenomenological aspects of visible neutrino decay
To describe visible neutrino decay we use the Majoron model with the following interaction Lagrangian where φ is a Majoron field, and g ij S and g ij PS represent, respectively, scalar and pseudoscalar couplings which are complex in general, with the neutrino mass eigenstate indices i, j = 1, 2, 3. Given the model Lagrangian (3.1), we have the two-body decay modes ν 3 → ν i /ν i + φ (i = 1, 2) in the case of NO, and ν i → ν 3 /ν 3 + φ (i = 1, 2) in the case of IO. As we stated in Section 1, we switch off the decay mode ν 2 → ν 1 /ν 1 + φ (i = 1, 2) in the IO. In this work, we assume that the Majoron is massless.
Phenomenology of neutrino decay depends crucially on the following two factors, • if neutrinos undergo visible or invisible decay, that is if the decay products are experimentally detectable or not 4 , • if the neutrino masses exhibit NO or IO.
We emphasize that the above, seemingly-obvious statements do indeed provide the key to understand the results in this paper. This fact is best summarized in Figure 1 in which the ν e disappearance probabilities in the absence or presence of the decay are plotted as a function of the anti-neutrino energy at a few characteristic distances: The top (L = 1.5 km), middle (L = 52 km) and the bottom (L = 180 km) panels correspond, respectively, to the far detectors in Daya Bay, JUNO, and the KamLAND experiments. The left (right) panels in Figure 1 are for the NO (IO). We note that the probability is shown in Figure 1 is the effective one in the sense that it is defined as the ratio of the ν e flux at the detector with oscillation plus decay effects to the flux without them. The former includes the contribution of daughter neutrinos which exists in the case of visible decay.  The left (right) panel is for the NO (IO). We show the effective probabilities (defined as the ratio of the νe flux arriving at the detector divided by the original flux at the detector in the absence of oscillation for a given neutrino energy) assuming standard oscillations without any decay (labeled "osc. w/o decay" in blue), with visible decay effect but without the daughter neutrino contribution (labeled "decay w/o app" in red) and with visible decay including the daughter neutrino contribution (labeled "decay w app" in black). For the visible decay, we consider gS = gPS = 0.2 for NO and gS = gPS = 0.1 for IO. Note that in some of the plots, the individual curves are too close together to be distinctly seen.
We first observe that the effect of decay is very noticeable in the case of IO (right panels) despite that the assumed magnitude of couplings for the IO case is smaller than that for NO , as can be seen in the right panels of Figure 1. It is because decay of the higher mass eigenstates ν 1 and ν 2 , which occurs copiously in the reactor-produced ν e flux, leads to a much stronger reduction of the survival probability P (ν e → ν e ) than the NO case (see below). The effect can be seen clearly in the right panels of Figure 1 with the Majoron coupling constants g S = g PS = 0.1. In the case of NO (left panels), on the other hand, the effect of decay is small, irrespective of whether the contribution from the daughter neutrinos is included or not. Unlike the case of IO, the reduction of the ν e flux is minor as the parent ν 3 component is small in reactor ν e due to suppression by small |U e3 | 2 = sin 2 θ 13 . Therefore, the effect of visible decay in the case of NO is just to dampen the atmospheric-∆m 2 31 driven neutrino oscillation [51].
In visible decay, an additional effect, a pile-up of events at low energies due to the daughter neutrino contribution, should be observed 5 . For the IO, for our choice of couplings g S = g PS = 0.1 this effect is barely noticeable by eyes in Figure 1 but only for longer baseline as in the case of KamLAND at lower energies as a small difference between the cases without (red curves) and with (black curves) daughter contributions. We confirmed that by using somewhat larger values of couplings the pile-up effect become more prominent but its effect is tiny in any case because of small (∝ |U e3 | 2 ) ν e component in the decay product ν 3 . The effect is negligible for the NO because the decay effect itself is small. We, of course, in any case scenario, NO or IO, always include the daughter neutrino contribution irrespective its importance, on which some comments may follow later.

The oscillation probabilities with neutrino decay
We first recapitulate the formulas of the neutrino oscillation probabilities in the simultaneous presence of flavor oscillations and decay [52,[63][64][65]. We treat the system as in vacuum which is a good approximation for the reactor neutrino experiments.

Neutrino decay: General formula
We start by examining a generic case in which each of the three massive neutrinos oscillate and decay at the same time. When a neutrino of flavour α with energy E α is produced at the distance L = 0, the differential probability that a neutrino of flavour β with energy in the interval E β + dE β is detected at the distance L, can be written as [52,63,65] (4.1) 5 We remark here that in the case of reactor neutrino experiments for the baseline of ∼ O(100) km, the area covered by the decay beam spread is expected to be < ∼ O(10 −2 ) m 2 which is much smaller than the detector sizes of KamLAND/JUNO. Therefore, we assume that in practice all the daughter neutrinos which are produced from the parent neutrinos emitted from the source toward the direction of the detector, reach the detector [63].
In Eq. (4.1), τ i and m i represent ν i 's proper lifetime and mass, respectively, and the indices r and s specify, respectively, parent and daughter neutrino helicities. The matrix element U (s) βi = U βi (U * βi ) corresponds to the case for positive (negative) helicity. The first term of the differential probability in Eq. (4.1) describes the contribution from a parent neutrino of flavor α which survived after propagating a distance L. While the second term which contains A ν r α →ν s β (E α , E β , L ) is the contribution of daughter neutrinos, with energy E β , produced through the decay of a parent neutrino at the distance L (< L) with energy E α , during propagation, with the potential possibility of helicity-flipped contribution.
The first term of Eq. (4.1) is often called the "invisible" contribution. But, the case that we examine in this paper has no invisible decay; the decay products always include active neutrinos and hence are always visible in principle. To prevent confusion we call the first and the second terms of Eq. (4.1) as the "parent" and "daughter" contributions, respectively. We remark that if a neutrino undergoes invisible decay, the differential probability is given by the first term of Eq. (4.1). Then, what is the difference between the invisible decay and the parent contribution of our visible decay? The answer is that in the case of un-observable final states (such as sterile neutrinos), the decay width Γ = 1/τ does not contain information about the final states. Whereas in our case Γ, which is computed with the Majoron model, does contain information of final states, such as the mass of the daughter neutrino. The lightest neutrino mass dependence of the event spectrum will be demonstrated in Section 7.

Parent contribution in visible neutrino decay
Let us calculate the contribution from the parent neutrinos i.e, the first term in Eq. (4.1). It gives the whole contribution in the case of invisible neutrino decay. By the nature of this term, helicity flip cannot be involved in it. Then, after integration over the neutrino energy E α we obtain (i, j = 1, 2, 3) Using the standard parametrization of the flavor mixing matrix [66], and substituting τ 1 , τ 2 → ∞ for the NO, we obtain for ν e → ν e channel 6 P parent ee (NO) = cos 4 θ 12 cos 4 θ 13 + sin 4 θ 12 cos 4 θ 13 + sin 4 And for the inverted mass ordering, on substituting τ 3 → ∞, we get

Daughter contribution in visible neutrino decay
We calculate the contribution of daughter neutrinos, the second term in Eq. (4.1). We assume that the couplings between the neutrinos and the Majoron are real quantities and hence there are no decay-related complex phases. Since we are interested only in the electron antineutrino disappearance probabilities here we drop the helicity indices r and s keeping in mind that the quantities correspond to antineutrinos and that only helicity preserving decays can be observed in reactor experiments. However, it should be noted that the full decay-widths include the sum over helicity-preserving as well as helicity-flipping partial decay-widths. In this work, we assume the CP-violating phase δ CP as well as the Majorana phases to be 0. Thus, we can also ignore the complex conjugation of the flavor matrix elements. The transition amplitude A να→ν β (E α , E β , L ) for the NO where ν 3 decays to ν 1 or ν 2 , is given by [52,63,65] Here E α represents the energy of the parent neutrinos, while E β = p β +m 2 d /(2E β ) represents the energy of the daughter neutrinos. Whereas the transition amplitude for the IO, where ν 2 or ν 1 decay to ν 3 , takes the form Here E α = p α +m 2 p /(2E α ) represents the energy of the parent neutrinos while E β represents the energy of the daughter neutrinos.
In the above equations, Γ ij is the partial decay width and W ij represents the normalized energy distribution function for the daughter neutrino for the decay ν i → ν j . The explicit formulas are given in the Appendix A.
Eqs. (4.5) and (4.6) describe the process in which ν α is produced at L = 0, propagates as the parent neutrino ν p to L , then it decays into the daughter state ν d at this distance L and is detected as ν β at L > L after traversing the distance L − L . In the NO, p = 3 and d = 1, 2, while in the IO p = 1, 2 and d = 3.
Using the Eq. (4.5) we can compute the visible decay term in Eq. (4.1), we find that for the NO, and for the IO,

Sketchy descriptions of KamLAND and JUNO
In this section, we briefly describe the details of KamLAND and JUNO, which are phenomenologically relevant to our work.

KamLAND
The KamLAND (Kamioka Liquid Scintillator Antineutrino Detector) reactor neutrino experiment consists of 1 kton of highly purified liquid scintillator detector based in Japan. KamLAND detects neutrinos coming from 16 nuclear power plants with a range of distances that go from 140 km to 215 km. The average distance corresponds to ∼ 180 km. The experiment ran in the reactor anti-neutrino mode from 2002 to 2012, collecting a total exposure of 4.90 × 10 32 target-proton-years. We consider the data presented in [61] to perform our analysis of neutrino decay. The information regarding backgrounds and systematic uncertainties have also been taken from [61]. The expected advantage of KamLAND over JUNO to study the decay effect is, as we could see in the plot of probabilities in Figure 1, the longer average baseline, which is about 3.4 times the JUNO's baseline, leading to larger decay effects.

JUNO
The JUNO (Jiangmen Underground Neutrino Observatory) experiment [62] is a future neutrino experiment that will be based in China. It is expected to start taking data from the year 2022. JUNO has been designed with the primary goal to measure the neutrino mass ordering but it will also be able to measure the oscillation parameters such as θ 12 , ∆m 2 21 and ∆m 2 31 with much better precision. The detector consists of a 20 kton fiducial mass of liquid scintillator and is located at an average distance of ∼ 53 km from Yangjiang and Taishan nuclear power plants. The remote reactor cores at Daya Bay and Huizhou will also have a small contribution to the total flux arriving at the JUNO detector. In our experimental set-up, we consider the various reactor cores with different thermal powers and baselines as described in Table 2 of [62]. The description of backgrounds and systematic uncertainties have been taken from [62]. The exposure is set such that a total of 140, 000 events are obtained (including the backgrounds).
The signal in both of the experiments is the inverse beta-decay (IBD) events in the energy range ∼ [1. 8,8] MeV, which essentially determines the electron antineutrino disappearance probability for a given, relatively well known IBD reaction cross-sections. The main background in a search for visible neutrino decay in JUNO is the geo-neutrino events at low energies. We consider their contributions similarly as done in [51]. The expected advantage of JUNO over KamLAND to study the neutrino decay is (i) much larger statistics and (ii) better energy resolution which is crucial for the case of NO as we will see later.
To simulate the KamLAND and the JUNO experiments, we have used the GLoBES [67,68] software package. The event rates and statistical-χ 2 calculations have also been performed using GLoBES.

Features of event rates in the presence of decay
In this section, using the observed and the expected event rates at KamLAND and JUNO, respectively, we stress upon the following two features that crucially affect the results in the presence of visible decay.
1. Dependence of event rates on the neutrino mass ordering, and 2. Dependence of event rates on the lightest neutrino mass 7 .
Note that the lightest neutrino is ν 3 in the case of IO, and ν 1 in the case of NO.
In Figs . We show the rates for standard oscillations without any decay (labeled "osc. w/o decay" in blue), visible decay without the daughter neutrino contribution (labeled "decay w/o app" in red) and visible decay including the daughter contribution (labeled "decay w app" in black). For the visible decay, we consider gS = gPS = 0.1. Also shown are the observed KamLAND data indicated by the black solid circles with error bars which are taken from [61].
background geo-neutrinos. In the left and the right panels, we assume 8 m lightest = 10 −3 eV and m lightest = 10 −1 eV, respectively. We assume the following values of the oscillation [69] and decay parameters to generate these event rates. . We show the rates for standard oscillations without any decay (labeled "osc. w/o decay" in blue), visible decay without the daughter neutrino contribution (labeled "decay w/o app" in red) and visible decay including the daughter contribution (labeled "decay w app" in black). For the visible decay, we consider gS = gPS = 0.2. Also shown are the observed KamLAND data indicated by the black solid circles with error bars which are taken from [61].
From Figure 2, as expected from the probabilities shown in the right panels of Figure 1, we see that the effects of decay are significant for IO. This is because, for IO, ν 2 or ν 1 mass eigenstate decays to ν 3 /ν 3 . Thus, the decay is expected to affect the ∆m 2 21 -driven oscillations. For KamLAND and JUNO, these oscillations are much larger in magnitude compared to ∆m 2 31 -driven-oscillations due to the large value of θ 12 . Therefore, decay effects are also large. Furthermore, when one considers the full visible decay including the contribution from daughter neutrinos, a pile-up of events at lower energies is noticeable for m lightest =10 −1 eV (see text below for the m lightest dependence on the decay effect). However, this is still a small effect as the appearance of ν e is suppressed due to the smallness of |U e3 | 2 .
We note that the case considered to generate the results shown in Figure 2, g S = g PS = 0.1 for the IO case, is turned out to be excluded as we will see later.
The decay effects in the event spectrum also depend on the lightest neutrino mass as can be seen by comparing the left and right panels in Figure 2. We first observe that, as long as the results shown in Figure 2 is concerned, for relatively small ( < ∼ 0.1) values of couplings, the lightest neutrino mass dependence comes mainly from the mass dependence in the visible part of the probabilities given in Eqs. (4.3) and (4.4). We remind the readers that the daughter contributions coming from Eqs. (4.7) or (4.8) in the effective probabilities shown in Figure 1 is quite small, which should be reflected in the event number distributions.
However, expressions shown in Eqs. (4.3) and (4.4) are not useful to understand the lightest neutrino mass dependence we can see in Figure 2 since the decay lifetime τ i appears in these equations also depend on neutrino masses. Therefore, we should take a closer look at the expressions of Γ functions given in the Appendix, taking into account that By looking into the expression of decay width Γ functions in Eq. (A.1), we can say that the origin of the lightest neutrino mass dependence comes from two parts: (i) the part which is given by the square of the mass of parent neutrino, m 2 i , a factor common for both helicity flipping and conserving processes, and (ii) the part which is described by the dimensionless functions f (x), h(x) and k(x) shown in the Appendix, which have dependence on the both parent and daughter masses of neutrinos as well as on the helicity of daughter neutrino, if it is flipped or conserved.
Let us first take a look the part (ii) which looks more complicated. For the case g S = g PS , the total rate coming from this part is proportional to the sum of (f (x) + h(x) + k(x))/x as we can see from Eq. (A.1) in Appendix A. We observe that the variation of the lightest neutrino mass have little impact on these functions (mainly due to f (x) which is dominant), and therefore induces little impact on the total rate, at most a factor of ∼ 2-3 for both mass orderings (see Figure 1 in Ref. [52] for the NO where (f (x) + h(x) + k(x))/x is shown as a function of x).
On the other hand, the part (i), mass square of the parent neutrinos, m 2 i , has stronger dependence on the lightest neutrino masses for the both mass orderings. In the IO, for the case of ν 2 → ν 3 decay, the parent mass is m 2 = m 2 lightest + ∆m 2 21 + |∆m 2 31 | which implies that the two different values of lightest neutrino mass lead to m 2 2 ∆m 2 atm = 2.40×10 −3 eV 2 for m lightest = 10 −3 eV, and m 2 2 m 2 lightest = 10 −2 eV 2 for m lightest = 10 −1 eV. As a result, the decay width is an order of magnitude larger in the case of m lightest = 10 −1 eV, leading to the small but visible difference between the left and right panels for the IO in Figure 2. In the NO, the situation is similar but m lightest dependence is not visible in Figure 3 because the impact of decay itself is small due to small value of θ 13 as explained in the end of Section 3. We must remind the readers that these particular dependences of the decay rate on the parent neutrino mass may be model-dependent, which should be kept in mind in interpreting our results.
From Figure 3, for the case of NO, it can be seen that the KamLAND experiment is almost insensitive to the decay of ν 3 to ν 2 /ν 2 or ν 1 /ν 1 . This is expected because the ∆m 2 31 -driven oscillations are averaged out at the baselines relevant to the KamLAND experiment, and therefore, the distortions in the spectrum due to decay cannot be seen. (Though a small pile-up of events due to daughter neutrinos is seen, KamLAND cannot place a useful bound on τ 3 in the NO case, probably because it is not statistically significant.) For the case of the JUNO experiment, we find that for the NO there is a very small effect of decay on the event rates. In the NO, s 2 13 suppressed ν 3 component decays into ν 2 /ν 2 or ν 1 /ν 1 , whose effect is mainly just to dampen the ∆m 2 31 -driven oscillations, as explored previously in Ref. [51].
In the case of IO, because of more than a factor of three longer average baselines of KamLAND which leads to the larger effect of neutrino decay than that for JUNO, Kam-LAND should be able to place a stronger constraint on ν 2 lifetime, if the number of events was similar to that of JUNO. However, this advantage is largely compensated by much lower statistics of KamLAND with the total number of events, 2611 (including backgrounds) [61], which is smaller than those assumed for JUNO by a factor of 54. Therefore, interpretation of the KamLAND bound on ν 2 lifetime, which is only slightly better than JUNO as will be reported in Section 7, must be done with care.

Constraints on visible neutrino decay by KamLAND and JUNO
In this section, we present the results of our analysis to obtain the constraints on visible neutrino decay imposed by the KamLAND data, and by a simulated data of JUNO assuming the total number of events equals to 140,000 which can be obtained by the exposure of ∼ 220 GW·years (total reactor thermal power times running period with ∼ 90% of reactor avilability assuming 100% detection efficiency). We exhibit the obtained constraints by drawing the 90% C.L. exclusion contours in the g S − g PS plane in Figure 4 for IO, and in Figure 5 for NO, for both KamLAND and JUNO. To translate these results to the constraints on the ratio of lifetime τ to the mass m, we show the equal τ /m contours as a function of g S and g PS in Figs. 4 and 5.

Analysis procedure
We now describe the numerical procedure for calculating the χ 2 for excluding decay. For KamLAND, we consider the data presented in [61] while for JUNO, we simulate the "true events rates" assuming that neutrinos undergo only standard oscillations and that no neutrino decay occurs. To simulate the true events rates in the case of JUNO, we take the values of the oscillation parameters used in Section 6.
We then "fit" the observed/simulated data with the calculated event rates assuming the existence of neutrino decay, by varying freely the values of g S and g PS in addition to varying the standard oscillation parameters θ 12 , θ 13 , ∆m 2 21 and |∆m 2 31 | in their currently-allowed 3σ ranges. Note that in the fit we keep the test mass ordering same as the true one. These events rates are called the "test event rates". The binned-χ 2 are calculated using GLoBES including the marginalization over the systematic uncertainties. We also add χ 2 due to the Gaussian priors corresponding to the test oscillation parameters that are varied in the fit. The formula for the total χ 2 is given by Here, n = 17 (n = 200) is the total number of energy bins for KamLAND (JUNO). F i and D i are the theoretical number of events and the observed number of events, respectively, in a given i-th bin. ξ k are the systematic uncertainty parameters with standard deviation σ k ; and θ bf j is the best fit value of a given oscillation parameter θ j (that are varied in the fit) with a 1σ uncertainty σ j . For both KamLAND and JUNO, we consider an overall normalization error of ξ 1 = 5% for signal and 20% for the background events and an energy calibration error of ξ 2 = 3%. For a given choice of the test decay parameters g S and g PS , we select the least χ 2 total that is obtained after marginalizing over all the test oscillation parameters. The ∆χ 2 for a given g S and g PS is obtained through: ∆χ 2 = χ 2 total − χ 2 total, smallest . We show the resulting contours corresponding to ∆χ 2 =4.61 for 2 degree of freedom (DOF) as a function of the test g S and g PS .

KamLAND and JUNO bounds on neutrino decay: the inverted mass ordering
We first discuss the potential of the experiments to exclude visible decay for the case of IO, shown in Figure Figure 1 in [52]. Expressed in terms of τ /m, for IO, JUNO excludes τ /m 1.1 × 10 −9 s/eV, which happened to be the same value as that of KamLAND. We see that inclusion of the daughter neutrino contributions does not affect in any significant way the decay exclusion sensitivity. In all the panels, the two curves with and without the daughter's contributions are nearly coincident.
It is remarkable that despite that the number of events obtained by KamLAND used for our analysis is 54 times smaller than that for JUNO (2611 vs 140,000), both experiments give very similar bounds (sensitivities). As mentioned at the end of the previous section, we understand that this is mainly because of more than 3 times larger average baseline of KamLAND (∼ 180 km) compared to JUNO (∼ 53 km) which can largely compensate the much smaller statistics of KamLAND.
Let us now try to understand qualitatively the dependence of the value of m lightest on the sensitivities to g S and g PS we can see in Figure 4. Since the contributions from daughter neutrinos are small, we just need to pay attention to the decay width Γ in Eq. (A.1), in particular for the helicity conserving case r = s which is dominant.
We first note that in the limit of vanishing m lightest , which corresponds to x → ∞, the functions f (x), h(x) and k(x) given in eq.(A.2) tend to become equal. See also Figure 1 of Ref. [52]. It implies that both couplings, g S and g PS , contribute equally to neutrino decay, which explains why the bounds are nearly symmetric to g S and g PS in the case of m lightest = 10 −3 eV.
On the other hand, as the assumed true value of m lightest is increased (or x is decreased), the function f (x)/x (h(x)/x) is increased (decreased), as can be seen from the first equation in (A.1) and also from Figure 1 of Ref. [52]. It means that the scalar (pseudo-scalar) coupling g S (g PS ) becomes more (less) important for decay. This is the reason why the bound on the scalar (pseudo-scalar) coupling become tighter (milder) at a large value of m lightest = 10 −1 eV independent of the mass orderings.

JUNO bound on neutrino decay: the normal mass ordering
We consider only the JUNO experiment as it was shown previously (see Section 6) that there is little sensitivity to the decay ν 3 → ν 1,2 in KamLAND. 9 JUNO will be able to place a limit on the decay effect because of the much larger statistics and good energy resolution which allows to detect the damping-like effect in the ∆m 2 31 driven oscillation we can see in the left middle panel of Figure 1. From Figure 5, we see that for m lightest = 10 −3 eV, JUNO can exclude neutrino visible decay at 90% C.L. if g S 0.28 and g PS 0.33. For m lightest = 10 −1 eV, JUNO can exclude neutrino decay at 90% C.L. if the mass ordering is normal and g S 0.16 and g PS 2.90. The constraint on g PS is much weaker mainly due to the same reason described in Section 7.3 for the IO case. Expressed in terms of τ /m, we find that for NO, JUNO excludes τ /m 7.5 × 10 −11 s/eV. 10 A comparison between Figs. 4 and 5 indicates that the bounds crucially depend on the choice of the mass ordering: the constraints are milder in the NO by an order of magnitude. In the NO, as in the case of IO, the daughter neutrinos give essentially no contribution to exclusion of decay, leading to almost complete degeneracy of the red and the black curves in Figure 5.   In Table 2, we summarize the constraints obtained by KamLAND and JUNO discussed in this section.

Conclusions
In this work, we have discussed the effect of visible neutrino decay which can be detected by observing the positron energy spectrum due to the IBD reaction in reactor neutrino experiments. Modifications of the spectrum not only in the shape but also in the normalization are important. We have obtained the constraints on the lifetime of higher-mass state neutrinos (ν 3 in the NO, and ν 2 or ν 1 in the IO) in medium and long-baseline reactor experiments.
We have used the Majoron model to calculate the neutrino decay rate. We took the two experimental settings, KamLAND and JUNO. They are the most relevant ones because of the long baseline, ∼ 180 km for KamLAND, and high energy resolution < ∼ 3% expected for JUNO in construction. For KamLAND we use the latest data, and for JUNO we assume the exposure of about 220 GW·years which would produce 1.4 × 10 5 events. In comparison with the constraints on neutrino decay often expressed as the bound on τ /m (τ 3 /m 3 for the NO, and τ 2 /m 2 for the IO) at 90% C.L. for 1 degree of freedom, we have provided the corresponding information in Table 3.  We found that the lifetime bounds depend crucially on whether the neutrino mass ordering is normal or inverted. Roughly speaking, the results we obtained for JUNO shows that for the IO, the bounds are better than the one for the NO approximately by a factor of 20. In looking into closer detail, KamLAND is insensitive to the decay of ν 3 in the case of NO, probably because of insufficient energy resolution to measure small wiggles of the atmospheric-scale high-frequency oscillations. JUNO, on the other hand, can rule out τ /m 1.0 × 10 −10 s/eV for the ν 3 mass eigenstate. We note that this value is quite similar and consistent with the bound of 9.3 ×10 −11 s/eV for the same confidence level (90% C.L.), obtained in [51] where the invisible neutrino decay for JUNO was studied. For the case of IO, the bounds of roughly the same order of magnitude are obtained on the decay of the ν 2 mass eigenstate by both KamLAND and JUNO. We find that τ /m 1.4 × 10 −9 s/eV is ruled out by both KamLAND and JUNO for the ν 2 mass eigenstate.
We have observed that, in each of these cases, there is no significant improvement in the sensitivity by including the daughter neutrino contributions in the analyses. It is because the effect is suppressed through |U e3 | 2 , which is small. However, in the case of IO, a decrease of ν e flux by the decay of ν 2 and ν 1 mass eigenstates produces a clear signature of neutrino decay, yielding a stringent lifetime bound mentioned above. For the NO, neutrino decay acts merely as a damping effect of the fast ∆m 2 31 -driven oscillations both in JUNO and KamLAND, rendering detection of decay effect harder, in particular, for the latter. We also mention that our lifetime bound depends on m lightest through m lightest dependence of the decay rate.

A Auxiliary formulae
We have the auxiliary functions, respectively the decay rate, Γ rs ij of ν j of initial mass state of mass m i and helicity r and final neutrino mass m j and helicity s , and normalized spectrum of daughter distribution, W rs ij (E α , E β ) with E α the energy of initial state and E β the energy of the final neutrino [52,63,65] Γ rs ij = f (x) = .