Constraint on Neutrino Decay with Medium-Baseline Reactor Neutrino Oscillation Experiments

The experimental bound on lifetime of nu_3, the neutrino mass eigenstate with the smallest nu_e component, is much weaker than those of nu_1 and nu_2 by many orders of magnitude to which the astrophysical constraints apply. We argue that the future reactor neutrino oscillation experiments with medium-baseline (~ 50 km), such as JUNO or RENO-50, has the best chance of placing the most stringent constraint on nu_3 lifetime among all neutrino experiments which utilize the artificial source neutrinos. Assuming decay into invisible states, we show by a detailed chi^2 analysis that the nu_3 lifetime divided by its mass, tau_3/m_3, can be constrained to be tau_3/m_3>7.5 (5.5) x 10^{-11} s/eV at 95% (99%) C.L. by 100 kt.years exposure by JUNO. It may be further improved to the level comparable to the atmospheric neutrino bound by its longer run. We also discuss to what extent nu_3 decay affects mass-ordering determination and precision measurements of the mixing parameters.


Introduction
Investigations on the possibility of neutrino decay has a long history, see e.g., [1,2]. Since neutrino radiative decay is so tightly constrained [3], decay into invisible final states are more commonly discussed, for example, in the context of majoron models [4,5]. 1 The bound on neutrino lifetime τ depend on (1) whether daughter neutrinos are active or sterile [7], and in the former case (2) which neutrino mass ordering, normal or inverted, is realized. It also depends on (3) whether the neutrinos are Dirac or Majorana particles. However, in any one of these cases whenever the bound exists, its order of magnitude is given by the condition τ /m ∼ L/E for neutrinos with mass m and energy E that traverses distance L [8]. We refer the condition as the kinematic estimate. It is nothing but stating that decay effect is sizeable when traveling time t is comparable to lifetime τ lab of neutrinos in the laboratory (i.e., observer's) frame. 1 To be testable by various means we mention here, one has to arrange the majoron models such that neutrino lifetime is not too long [6].
Astrophysical neutrinos, because of their long path lengths, are the promising sources for yielding the stringent bounds on neutrino lifetime. The solar neutrinos have been used to place bound on the lifetime of ν 2 , which is the dominant component of 8 B neutrinos under the assumption that ν 2 decays into ν 1 (ν 1 ) or into sterile states [8][9][10][11][12]. In most cases the daughter neutrinos were assumed to be unobserved even if the decay into active neutrino is considered. The obtained bound is of the order of the kinematic estimate, τ 2 /m 2 ∼ 10 −4 s/eV. We note that the authors of ref. [11] argued that it is possible to constrain also the lifetime of ν 1 by considering low energy pp and 7 Be solar neutrinos, which have large ν 1 components. Because of lower neutrino energies, the obtained bound for ν 1 is better by about a factor of five than that for ν 2 [11].
Supernova neutrinos are potentially the most powerful source for bound on neutrino decay [13] which could lead to the bound τ /m ∼ 10 5 s/eV according to the kinematic estimate. However, the available data is currently limited to the one which came from SN1987A [14,15]. Moreover, the bound applies only to ν 1 and/or ν 2 which have largē ν e components. Astrophysical neutrinos which have been observed by IceCube [16], in principle, are equally (or more) powerful as supernova neutrinos. But, to place the bound on τ /m we need either identification of the sources, or complete determination of the neutrino flavor ratios [17]. Its typical order is estimated as τ /m ∼ 10 4 s/eV assuming 1 TeV neutrinos from AGN at 100 Mpc [8]. See [18] for recent discussions on the hypothesis of neutrino decays over cosmological distances.
Leptonic decay of mesons leads to a bound on the majoron coupling constant with neutrinos g αβ (α, β = e, µ, τ ) [19] of the order of g 2 ∼ 10 −5 − 10 −4 . For a comprehensive treatment of the majoron coupling bound with the pseudo-scalar as well as the scalar couplings, see e.g., [20] and an update [21]. When translated into ν 2 lifetime assuming decay into active ν 1 , it gives a bound on τ 2 /m 2 in a range comparable with to somewhat weaker than the solar neutrino bound [8]. Using lepton decay channels the authors of [20] also obtained the bound of couplings which include ν τ , α=e,µ,τ |g τ α | 2 < 0.1, which would lead to a less stringent bound on τ 3 /m 3 for active ν 3 decay.
Therefore, so far no stringent bound on ν 3 lifetime appears to exist either from astrophysical or from laboratory neutrinos. Probably, the best way to constrain ν 3 lifetime would be to use neutrino oscillation phenomenon. It is because in this case one can select the channels or region of kinematical phases that are sensitive to ν 3 decay effect. Then, it is natural to use the oscillation driven by ∆m 2 32 ∆m 2 31 , hereafter referred to as atmospheric-scale neutrino oscillation for simplicity. We will see that despite quantum mechanical nature of the phenomenon the kinematic estimate applies. The bound on ν 3 lifetime was obtained [23] by using the data collected by the Super-Kamiokande (SK) atmospheric neutrino observation [24] as well as by the long-baseline accelerator experiments. See [25] for a recent update of the accelerator bound.
In this paper, under the assumption of decay into invisible daughters, we analyze the bound on ν 3 lifetime which will be placed by the future medium-baseline reactor neutrino experiments such as JUNO [26] or RENO-50 [27] via observing the distortion of neutrino energy spectrum. In fact, we argue that they could provide, in principle, the most stringent bound on ν 3 lifetime among all the artificial source neutrino experiments. It is because they can observe the effect of atmospheric-scale oscillations at the baseline around the maximum of the solar-scale oscillation, L OM = 4πE/∆m 2 21 , which is longer than L OM for the atmospheric-scale oscillation by a factor of 30. Assuming five years operation of JUNO, we obtain the bound τ 3 /m 3 > ∼ 7.5 (5.5)×10 −11 s/eV at 95% (99%) CL. See section 2 for details, and for the relationship between our argument and the bound obtainable by using atmospheric neutrino data. The principal objective of the medium-baseline reactor neutrino experiments is to determine the neutrino mass ordering. Then, the immediate question is whether it would be disturbed if possibility of neutrino decay is taken into account in fitting the data, or more drastically, when ν 3 would actually decay. It was also noticed that such experiment has a potential of determining the mixing parameters in a high precision [28,29]. In fact, the recent works with much more elaborate treatment of experimental errors reported an extreme precision of sub-percent level for sin 2 θ 12 and ∆m 2 21 [30,31]. Then, the natural question is whether or to what extent these sensitivities could be affected when possibility of ν 3 decay is turned on. In section 6 we address these questions.

Uniqueness of the medium-baseline reactor neutrino experiments
In this section we try to convince the readers that JUNO/RENO-50 is, in principle, the highest sensitivity experiment among all those which utilize artificial source (or beam) neutrinos in detecting the possible decay effect of ν 3 . In this paper, we consider the case of invisible decay of ν 3 . That is, we assume that the decay products are either some sterile states, or can involve active ν 1 and/or ν 2 state but with significant energy degradation such that daughter neutrinos cannot be observed. The latter possibility necessitates the neutrino mass ordering to be the normal type.
When i-th mass eigenstate neutrino decays with lifetime τ i at rest, the energy E i (more precisely, the energy difference normalized appropriately) of propagating i-th mass eigenstate neutrino can be written as is a Lorentz dilated lifetime. As will be shown in appendix A the decay of ν 3 produces the following two characteristic modifications in the oscillation probabilities in vacuum: 2 • Reduction of the atmospheric-scale oscillation amplitude with the form cos • Decrease of normalization of the probabilities by the amount proportional to 1 − e −Γ 3 L .
Here, L is the distance traveled by neutrinos. To make the point clearer we have used an approximation ∆m 2 31 ≈ ∆m 2 32 ≡ ∆m 2 atm where ∆m 2 ji ≡ m 2 j − m 2 i . The first feature stems from the fact that the effect comes from the interference between the first−second mass eigenstates and the third. The second feature represents the effect of decaying i-th mass eigenstate projected onto the initial and final neutrino flavor states. It is independent of neutrino oscillation and exists even in the limit of vanishing ∆m 2 ji , as it must. This effect is prominent in ν µ (andν µ ) disappearance channel, but is negligible for ν e (andν e ) survival probability. Whereas the first effect is equally important in the all oscillation channels, as will be shown in appendix A.
The decay effect is negligible for baseline L with Γ 3 L 1, and it becomes significant only when Γ 3 L > ∼ 1. In most of the oscillation experiments which use man-made neutrino sources, the baseline is set to the first oscillation maximum of atmospheric-scale oscillation, ∆m 2 atm L/(2E) π. In such setting, it is likely that we obtain the lower bound on the width Γ 3 as where the last equality in parenthesis assumes that the distance L is at the first oscillation maximum. It implies that the kinematic estimate applies. Using (2.4), one can estimate the upper bound on τ 3 /m 3 at In the JUNO or RENO-50 setting, hereafter referred simply as JUNO setting for simplicity, the detector is placed approximately at the maximum of the solar-scale oscillation, ∆m 2 21 L/(2E) π. Then, we can achieve about 30 times longer lifetime bound It should be stressed that the JUNO setting is unique (among artificial neutrino source experiments) in making observation of the atmospheric-scale neutrino oscillation at the distance of the solar-scale oscillation maximum possible. This completes our argument that JUNO is the highest sensitivity experiment among (practically) all the ongoing or proposed artificial-source neutrino experiments as far as the ν 3 lifetime bound is concerned. Of course, the above argument does not necessarily imply that the JUNO bound on ν 3 lifetime must be the severest one achievable by all the neutrino experiments. The likely (and probably unique) exception is the one placed by observation of the atmospheric neutrinos. The naive kinematic estimate for the ν 3 lifetime would lead to the sensitivity to which suggests that the bound by the atmospheric neutrinos is comparable as the one by JUNO. However, in the case of atmospheric neutrino experiments they can observe neutrinos in much wider energy and baseline ranges than those used in the above estimate. Therefore, we expect that the bound on τ 3 /m 3 from the atmospheric neutrinos is tighter than the naive estimate given in (2.7). In fact, the authors of ref. [23], by using the SK atmospheric neutrino data and the others, obtained the bound τ 3 /m 3 > 9.3 × 10 −11 s/eV at 99% CL, which is stronger than the one shown in Eq. (2.7) by about a factor of 3. 3 3 Effect of ν 3 decay on the oscillation probabilities and the observable Theν e survival probability relevant for reactor neutrinos for the baseline L in vacuum is given, under the approximation ∆m 2 31 ≈ ∆m 2 32 ≡ ∆m 2 atm , as P (ν e →ν e ) = 1 − c 4 13 sin 2 2θ 12 sin 2 ∆m 2 where c ij ≡ cos θ ij and s ij ≡ sin θ ij . See appendix A for derivation. For simplicity and as a good approximation, we ignore the matter effect in this work. As we stated in the previous section, there are two types of terms which are affected by the neutrino decay. However, they come with vastly different magnitudes in theν e channel; since the current neutrino data implies s 2 13 0.02, the coefficients of 1 − e −Γ 3 L and cos terms are, respectively, s 4 13 ∼ 5 × 10 −4 and 1 2 sin 2 2θ 13 ∼ 4 × 10 −2 . Therefore, the former oscillation-independentν e attenuation effect should be negligible.
Then, the question is how such decay-affected oscillation probability manifests itself into the observable quantities. To give a feeling to the readers we show in the upper panel of Fig. 1 the energy spectrum of events as a function of positron deposited energy, or visible energy E vis , 4 calculated by convoluting the neutrino flux and the cross section of the inverse β-decay (IBD) (ν e + p → e + + n) reaction. The three curves with different combinations of colors and line types correspond to the following three cases: no decay (black solid line), τ 3 /m 3 = 5 × 10 −11 s/eV (red solid line), and τ 3 /m 3 = 10 −11 s/eV (green dashed line). See appendix B for the procedure to obtain these curves. As we can see from the upper panel in Fig. 1, the decay of ν 3 state tends to average out the fast oscillation (wiggles) driven by ∆m 2 atm . 5 Since decreasing θ 13 and ν 3 decay both reduce the amplitude of atmospheric-scale oscillation, they could potentially be confused with each other. Fortunately, the confusion is precluded by the precision measurement of θ 13 done by the short baseline reactor neutrino experiments. (See the related comment at the end of section 5.) Currently, sin 2 2θ 13 is measured with the uncertainty of 6% [34], while a possibility of reaching the ultimate error of 3% by the end of 2017 is mentioned in [35].
Degradation of the oscillation amplitude can also occur by a finite energy resolution of the detector, which would cause another confusion with the decay effect. To illustrate this point, we show in the lower panel of Fig. 1, the difference between the energy spectra in the absence and in the presence of ν 3 decay for two cases of energy resolution, 3% (red 3 For comparison we note that the kinematic estimate for JUNO, (2.6), is different from our results based on χ 2 analysis by 40% only. It should also be noticed that if a re-analysis in [23] would be done with the currently accumulated SK data, the bound should become even tighter. It is also more stringent by a factor of 1.7 than our JUNO five-years bound to be obtained in section 5. 4 The visible energy Evis is approximately related to neutrino energy E as Evis E − (mn − mp) + me, where mn, mp, me, are, respectively the mass of neutron, proton and electron. 5 In this connection we note that quantum decoherence [32] has the similar effect on energy spectrum ofνe, and the relevant phenomenology is discussed with the JUNO setting in ref. [33].  solid curve) and 6% (blue dashed curve) at 1 MeV. Roughly speaking, doubling the energy resolution causes an amplitude attenuation of the signal by a factor of three. We observe that the detectability of the decay effect is strongly dependent on the energy resolution. Fortunately, the 3% energy resolution is to be reached by JUNO and RENO-50 [26,27]. To our understanding, the remaining questions which need to be addressed are as follows: • As shown in Fig. 1, the decay effect is to reduce the atmospheric-scale oscillation amplitude. This is the oscillation to be utilized to determine the mass ordering in JUNO. Then, the obvious question is to what extent the mass ordering determination could be disturbed if the effect of ν 3 decay is taken into account.
• Another relevant question would be to what extent the sensitivities to mixing parameter measurement in JUNO could be disturbed by ν 3 decay.
We will discuss these issues in section 6 from the following two different viewpoints (assumptions): (1) There is no decay in the input data set, but we consider the decay effect in the output (fit), and (2) ν 3 actually decays with the lifetime which is marginally consistent with the current and the expected JUNO bound.

Analysis method
For definiteness, throughout the paper (including Fig. 1), we assume that the true (input) values of the oscillation parameters are given as follows, which are taken from the best fitted values of the one of the recent global analysis [36], and assume the normal mass ordering (∆m 2 31 > 0), unless otherwise stated. In our statistical analysis, we formally divide the χ 2 function in three terms as as done in [30,31]. The first term, χ 2 stat , is computed by taking the limit of infinite number of bins which is justified because of the large number of events expected at 20 kt detector, 1.4 × 10 5 for 5 years of operation, as follows [30,31], where dN obs /dE vis is the event distributions of the observed (simulated) signal, and ξ i is the flux normalization parameters for reactor neutrinos as well as for geoneutrinos to be varied freely subject to the pull term in χ 2 sys (see below) and we integrate up to E max vis = 8 MeV. In our analysis, following [31], we include the contributions not only from Yangjiang and Taishan reactor complexes at L = 52.5 km (approximated by a single reactor with the thermal power of 35.8 GW) but also the contributions coming from the reactor complexes located at Daya Bay (L = 215 km) and Huizhou (L = 265 km), as well as geoneutrinos. See appendix B for details.
The second term in (4.2) takes into account the current uncertainties of the standard oscillation parameters and is given by wherex i and x fit i (i = 1-4) denote, respectively, the assumed true (input) and fitted values with For the values of σ(x i ), we take the current 1 sigma uncertainties determined by the global fit [36], σ(sin 2 θ 12 ) = 4.1%, σ(∆m 2 21 ) = 2.4%, σ(sin 2 θ 13 ) = 4.6% and σ(∆m 2 31 ) = 1.9%. The last term in (4.2), χ 2 sys , takes into account the contributions of two kind of experimental systematic uncertainties we consider where σ ξreac , σ ξ U and σ ξ Th describe, respectively, the normalization uncertainties for reactor neutrinos, uranium and thorium induced geoneutrinos. For them we take σ ξreac = 3% [37,38] for reactor neutrinos and σ ξ U = σ ξ Th = 20% for geoneutrinos following ref. [31]. In addition to the normalization uncertainties, we also consider the uncertainty of the energy resolution by including the last term in (4.5) with σ η (see below). For simplicity, we assume 100% detection efficiency, det = 1. The effect of uncertainty in the detection efficiency (IBD selection, fiducial volume, etc) is thus neglected, since it is difficult to reliably estimate. However, those errors, being energy independent, should hardly affect the sensitivity to shape-modulating effects caused by the neutrino decay as well as by the differing mass orderings.
With regard to the energy resolution, we only consider the stochastic term, i.e., σ E /E = 0.03(1 + η)/ E/MeV contribution where η is the parameter accounting for the energy resolution variation. 6 The non-stochastic term(s) are expected to be the most relevant, even dominant, at energies > 3 MeV. However, this is a complex experimental matter impossible to anticipate at this stage, so its impact is simply neglected here. Some deterioration of the sensitivity should be expected, if those terms were considered. Because of the illustrated dependence of the energy resolution on the decay sensitivity, its uncertainty is taken into account by using a pull term with σ η = 10%. 7 We feel it appropriate to consider no other systematic error on the energy scale at this stage.
The χ 2 is computed in the following way: In order to derive the expected bound on the ν 3 decay lifetime, for our input data, we consider the JUNO setup to compute dN obs /dE vis assuming no decay (τ 3 /m 3 = ∞) using the oscillation parameters given in (4.1). Then we try to fit such dN obs /dE vis by minimizing the χ 2 function varying freely θ 12 , ∆m 2 21 , θ 13 , ∆m 2 31 , ξ reac , ξ U , ξ Th , η, and τ 3 /m 3 subject to the pull terms in (4.4) and (4.5). In this way we can estimate the sensitivity to τ 3 /m 3 and derive the bound on it, and at the same time, can determine the allowed regions of other parameters as well.
Using the χ 2 function, we will determine the allowed range of τ 3 /m 3 by the condition ∆χ 2 ≡ χ 2 − χ 2 min = 2.71, 3.84 and 6.63 (1, 4 and 9), (4.6) at 90%, 95% and 99% (1, 2 and 3 σ) CL, respectively, for one degree of freedom. For the case where we show the allowed regions in the plane spanned by any combination of two out of nine parameters, we use the condition ∆χ 2 = 2.3, 6.18 and 11.83, respectively, for 1, 2 and 3 σ CL for two degree of freedom. We note that χ 2 min = 0 by construction because we do not take into account the statistical fluctuation in simulating the artificial data.

The bound on τ 3 /m 3
In this section, we derive the bound on τ 3 /m 3 assuming that ν 3 decays into invisible states, whereas ν 1 and ν 2 are stable. For definiteness, we assume the normal mass ordering. But, given the expression of P (ν e →ν e ) in (3.1), we expect that our results are valid also for the case of inverted mass ordering.
As our standard setup, we assume 5 years running of the JUNO detector with the fiducial volume of 20 kt which is placed at L = 52.5 km away from an effective single reactor with the thermal power of 35.8 GW. We assume the detector's running with 100% efficiency, det = 1. But, if det < 1 the running time must be scaled to [5/ det ] years to obtain the same results. Since JUNO may take data for a longer period, we also consider the case of exposure for 15 years. See appendix B for details of our calculation. . ∆χ 2 ≡ χ 2 − χ 2 min is shown, by the red (blue) curves for 5 (15) years of data taking, as a function of the fitted value of τ 3 /m 3 calculated for the JUNO detector placed at L = 52.5 km from a reactor with 35.8 GW thermal power, assuming 5 years of exposure and 100% detection efficiency. We have taken that the true (input) value of τ 3 /m 3 is infinite (stable ν 3 ). The solid curves correspond to the results obtained by using our full χ 2 defined in (4.2) whereas the dashed ones correspond to the case without assuming systematic errors. The contributions from the reactors at Daya Bay and Huizhou as well as those from geoneutrinos are taken into account. The bound comes from the SK atmospheric neutrinos plus long-baseline oscillation experiment obtained in [23] is also indicated by the vertical black dashed line.
In Fig. 2 we show ∆χ 2 ≡ χ 2 − χ 2 min as a function of the fitted value of τ 3 /m 3 where the input (true) value of τ 3 /m 3 is assumed to be infinity, i.e., ν 3 is stable, with oscillation parameters given in (4.1). All the other eight parameters are marginalized in the fit. The case of exposure for 5 (15) years is shown by the blue (red) solid curve. To exhibit the effect of the systematic error onto the bound, we also show by the dashed curves ∆χ 2 for the case without all the systematic errors.
From our result shown by the solid blue curve in Fig. 2, we can conclude that if the data is consistent with no-decay hypothesis, the range can be excluded at 95 (99%) % CL by 5 years exposure by JUNO.
As mentioned in section 2, this bound is about a factor of 1.7 weaker than the bound τ 3 /m 3 > 9.3×10 −11 s/eV at 99% CL based on the data of atmospheric neutrinos (as well as of long-baseline oscillation experiments) [23]. The bound is indicated by the vertical black dashed line in Fig. 2. By considering an extended running of 15 years, the JUNO bound can be improved to 11 (8.5) ×10 −11 s/eV at 95 (99%) % CL, which is barely comparable to the atmospheric neutrino bound. A comparison between the bound of the order of τ 3 /m 3 > ∼ 10 −10 s/eV obtained here and in [23] to the one deduced by using lepton decay channel [20] is described in [23]. They obtained the bound on ν 3 − ν s −majoron coupling as |g s3 | 2 < ∼ 10 −4 (m s /1 eV) −2 at 90% CL assuming m 3 m s , where m s denotes the sterile neutrino mass.
We also considered the hypothetical situation where θ 13 were unknown at the time of JUNO running (not shown). We did it by removing the pull term for θ 13 from our χ 2 , and found that the bound on τ 3 /m 3 would become worse by about a factor of three. It is because the ν 3 decay and the effect coming from the uncertainty in θ 13 can be confused with each other, as we mentioned in Sec. 3.
6 Impact of ν 3 decay on the determination of the mass ordering and the oscillation parameters 6.1 Impact of ν 3 decay on the mass ordering determination Now, we address the question of effect of ν 3 decay onto the mass ordering determination in JUNO/RENO-50, as promised at the end of section 3. We use the same analysis tool as used to obtain the bound on ν 3 lifetime in the previous section. The input mass ordering is always taken to be normal in this section.   . ∆χ 2 ≡ χ 2 − χ 2 min is shown as a function of the fitted value of ∆m 2 ee . The solid blue and red curve correspond, respectively, to the case where the fit is performed assuming the normal (right) and inverted (wrong) mass ordering (MO) for the standard oscillation. The blue (red) dotted and dashed curves corresponds, respectively, to the case where the input value of τ 3 /m 3 is ∞ and 10 −10 s/eV for the case of normal (inverted) mass ordering.
To know the effects of ν 3 decay on the resolution capability of the mass ordering in JUNO, we compute ∆χ 2 ≡ χ 2 − χ 2 min by taking both the normal (right) and in-verted (wrong) mass ordering. 8 In Fig. 3, ∆χ 2 is plotted as a function of ∆m 2 ee obtained by marginalizing all the other parameters. For the abscissa in Fig. 3 we use ∆m 2 ee ≡ |c 2 12 ∆m 2 31 + s 2 12 ∆m 2 32 |, which we believe to be the appropriate variable to discuss resolution of the mass ordering in medium baseline reactor experiments [39]. It is proposed as the effective atmospheric ∆m 2 determined byν e disappearance experiment in vacuum [40], which agrees in a good approximation with the one measured by the reactor θ 13 experiment [41].
The blue and the red curves are for the normal (input) and the inverted (wrong) mass orderings, respectively. The three line types, which are common to the both mass orderings, correspond respectively to: (i) Solid curve: no ν decay in both the input and the fit, the case of standard oscillation, (ii) Dotted curve: no ν decay in the input but ν 3 decay is allowed with prior-unconstrained lifetime in the fit, and (iii) Dashed curve: ν 3 decay with lifetime τ 3 /m 3 = 10 −10 s/eV is assumed in the input and allowed in the fit.
The global features of the results can be summarized as: • From (i) to (ii) there is only minor change (up to ∼ 1) in ∆χ 2 . It means that allowing the decay only in the fit little affects the output.
• From (i) to (iii) there are appreciable changes in ∆χ 2 in a manner which depend very much on the mass ordering. In the case of normal (right) mass ordering the change in ∆χ 2 upon turning on ν 3 decay both in the input and output is modest, slightly opening up the Gaussian parabola.
• If ∆χ 2 difference between the right (input normal) and the wrong mass orderings without decay is denoted as ∆χ 2 no-decay , ∆χ 2 difference at the minima across the different mass orderings becomes ∆χ 2 no-decay − 5 when the ν 3 decay is turned on with lifetime comparable to the current and the JUNO bound, τ 3 /m 3 = 10 −10 s/eV. Thus, we have observed that the ν 3 decay has a big impact on mass ordering resolution in JUNO, significantly worsening the sensitivity.
The readers must be noticed that we carefully avoided to make a quantitative statement on how large is ∆χ 2 no-decay , but restricted ourselves into the change due to the decay effect. It is likely that our procedure to simulate the distortion of the event number energy distribution is insufficient to reliably extract the absolute confidence level for the mass ordering determination. In particular, we do not take into account the uncertainties related to the non-linearity of the energy measurement, whose control would be the key to the success in the mass ordering measurement in JUNO. 9 On the other hand, the omitted error in computing χ 2 may not affect our discussion of the effect of decay in a significant way, because it affects the spectrum as simple reduction of the oscillation amplitude. 10 8 At the end of this section 6.1, we will place a long clarifying remark on how to interpret our ∆χ 2 . 9 This issue was raised in [39], and the crucial importance of controlling the energy-scale errors on mass ordering determination was illuminated with experimental perspectives by the authors of ref. [42]. 10 In addition, there is a subtle issues related to statistical treatment of mass ordering determination. In the usual case one can attribute to ∆χ 2 min the meaning of 1σ significance. The special feature that the fitting parameter is a discrete variable requires a different treatment to evaluate correctly the CL for the mass ordering determination [31,43,44]. Since we do not intend to elaborate this point, and because of 6.2 Impact of decay on the determination of the oscillation parameters We briefly discuss in this section possible effects of ν 3 decay on determination of the standard oscillation parameters in JUNO as well as on the values of output uncertainties. Detailed features of parameter correlations with and without ν 3 decay, including the correlations with the systematic errors, will be discussed in appendix C. As in the previous section 6.1, we consider the three cases, (i) standard oscillation, (ii) decay effect only in the fit and (iii) decay effect both in the input and in the fit. For the finite input value of lifetime we always use τ 3 /m 3 = 10 −10 s/eV.
To present our results in a clear way, we always use the following line symbols throughout figures 4 (this section), 5 (appendix C.1), and 6 (appendix C.2): For the case (ii) we show the 1, 2 and 3σ allowed regions, respectively, by the filled blue, yellow and green colors. For the case (i) the allowed regions are delimited by the black solid curves, and for (iii) by the black dotted curves.
In Fig. 4 we show the allowed regions of the parameters in (a) sin 2 θ 12 -∆m 2 21 (left panel), and (b) sin 2 θ 13 -∆m 2 31 (right panel) space. We observe in the left panel that there is no visible effect of the decay on the determination of the parameters in the 1-2 sector of the MNS matrix, sin 2 θ 12 and ∆m 2 21 . Whereas in the right panel, it is revealed that accuracy in measurement of ∆m 2 31 (of sin 2 θ 13 also but only slightly) is affected when ν 3 decays.  To obtain a more comprehensive view of effect of decay and to represent the change in the sensitivities in a more quantitative way, we present in Table 1 the fractional uncertainties (in %) of all the parameters determinable in the three cases (i-iii) in the JUNO setting. In doing so we do not restrict to the oscillation parameters, but also include the ones which describe the systematic errors. From Table 1 and Fig.4 we observe: the crude nature of our ∆χ 2 construction, we recommend the readers to use the information presented in Fig. 3 only to have a feeling on how ν3 decay affects the determination of the mass ordering in JUNO. Fractional errors in percent (at 1σ) of the oscillation as well as the parameters for systematic uncertainties, before (prior) and after the fit to 5 years of JUNO data with 100% detection efficiency, for the case where the true mass ordering is normal. They are determined by the condition ∆χ 2 = 1 for one degree of freedom. We consider three cases (i) standard oscillation fit without decay (ii) decay sensitivity fit (with true τ 3 /m 3 = ∞) and (iii) decay measurement fit (with true τ 3 /m 3 = 10 −10 s/eV • Comparison between the columns (i) and (ii) indicates that if the decay effect is considered only in the fit, its impact is virtually absent except for θ 13 and the energy resolution uncertainty parameter η.
• Comparison between the columns (i) and (iii) tells us that if the decay effect is considered for both in the input and the fit, there is appreciable impact of the size of approximately 10%-30% but only for sin 2 θ 13 , ∆m 2 31 and 1+η.
As a whole, we conclude that the impact of the decay on the determination of the mass and mixing parameters is rather small.

Conclusions
In this paper, we have investigated the question of how strong a constraint on decay lifetime of the massive neutrino state ν 3 can be placed by the medium-baseline (L ∼ 50 km) reactor neutrino experiments, JUNO or RENO-50, which we referred simply as JUNO in most part in the text. Assuming decay into invisible states and τ 1 , τ 2 τ 3 , we found that the bound τ 3 /m 3 > 7.5 (5.5) × 10 −11 s/eV at 95% (99%) CL can be obtained by JUNO with its five years exposure with 100% efficiency.
In fact, there is a simple reason why the JUNO setting can offer the chance of deriving the most stringent bound on the neutrino lifetime among all the experiments which utilize neutrinos from the artificial sources. Given that such experiments are designed to have sensitivities at around the first oscillation maximum (in exploring the atmospheric-scale neutrino oscillations), the kinematic estimate implies τ /m ∼ L/E 2π/∆m 2 atm . However, if there is such an experiment that can explore atmospheric-scale oscillations at the distance of solar-scale oscillation maximum, the bound on τ /m would become severer by a factor of ∆m 2 atm /∆m 2 21 30. This is what JUNO does. The bound on τ 3 /m 3 we obtained for 5 years exposure of JUNO has the same order of magnitude as (but is somewhat weaker than) the one obtained with the atmospheric neutrino data. We have examined the question of whether the 15 years running of JUNO can tighten up the bound to the level of the current atmospheric neutrino bound, and obtained an affirmative answer. It is nice to see that the comparably strong ν 3 lifetime bound can be deduced by the two completely different experiments, one observing reactor neutrinos at the medium baseline, and the other measuring the atmospheric neutrinos.
We have also discussed to what extent the ν 3 decay could affect the determination of the mass ordering as well as the precision parameter measurement by JUNO. We considered the two cases (in the notation in section 6), (ii) the decay effect is absent in the input data set but allowed in the output (fit), and (iii) ν 3 decay in input with the lifetime τ 3 /m 3 = 10 −10 s/eV, which is marginally consistent with the current and our JUNO bound. For the mass ordering determination, we found the significant impact of the decay (reduction of ∆χ 2 by ∼ five units) but only for the case (iii). With regard to influence of decay to the measurement of the oscillation parameters, we found that the impact is rather small except for ∆m 2 31 again only for the case (iii), in which the uncertainty of ∆m 2 31 would become ∼ 30% larger.

A Neutrino oscillation probabilities in the presence of neutrino decay in vacuum
The neutrino evolution in vacuum with neutrino decay is governed by the usual Schrödinger equation by replacing the neutrino energy E i by the one in (2.1). The S matrix whose elements describe neutrino flavor transition in vacuum as ν α (x) = S αβ (x)ν β (0) is given by where x is the distance traveled by neutrino and the notations for U αi , m i , E, and Γ i are the same as in section 2. Then, the general expression of the three-flavor expression of oscillation probability with neutrino decay in vacuum is given by [45] The first term in (A.2), the oscillation probability in the absence of neutrino decay, is given by the familiar expression, If the approximations Γ 1 x ∼ Γ 2 x 1, and Γ 3 x ∼ 1 holds we can simply (A.2). The disappearance oscillation probability takes the particularly simple form as where in the last line we have used the approximation ∆ 31 ≈ ∆ 32 ≡ ∆ atm . 11 In the appearance channels (α = β), the oscillation probability also simplifies: Under the approximation ∆ 31 ≈ ∆ 32 ≡ ∆ atm and because of anti-symmetry of Im[U αi U * βi U * αj U βj ] (= Im[C αβij ]) under i ↔ j, CP violating term in the decay-width dependent part of the oscillation probabilities vanishes.
Some remarks are in order: • In ν µ disappearance channel the coefficients of 1 − e −Γ 3 x and cos ∆ atm terms are given by |U µ3 | 4 = s 4 23 c 4 13 0.24 and 2|U µ3 | 2 1 − |U µ3 | 2 0.50, respectively. Therefore, the oscillation independent ν µ depletion effect may be as important as the effect of diminishing amplitude in the atmospheric-scale oscillations.
• In ν e appearance channel the corresponding coefficients of the two Γ 3 -dependent terms are given by |U e3 | 2 |U µ3 | 2 = s 2 23 c 2 13 s 2 13 1.1×10 −2 and twice of that 2.2×10 −2 , respectively. They are comparable with each other, and are similar in magnitude as the main oscillation term in vacuum.

B Event energy spectrum
The distribution of the number of events coming from the inverse β-decay (IBD) reaction, ν e + p → e + + n, as a function of the visible energy is given by, where n p is the number of target (free protons), t exp is the exposure, det is the detection efficiency, dφ i (E)/dE is the differential fluxes of reactor neutrinos and geoneutrions, dσ(E ν , E e )/dE e is the IBD cross section, P i (ν e →ν e ; L i , E) is theν e survival probabilities and R(E e , E vis ) is the Gaussian resolution function (see below).
For simplicity, we set det = 1, and as a reasonable approximation for our purpose, we ignore the neutron recoil in the IBD reaction and simply assume that neutrino energy, E, and the positron energy, E e , is related as E e = E − (m n − m p ) E − 1.3 MeV. Due to the finite energy resolution, the event distribution can not be obtained as a function of E e (true positron energy) but as a function of the reconstructed or so called visible energy, E vis , which is approximately related to neutrino energy as E vis E − (m n − m p ) + m e , after taking into account the energy resolution (see the text below). Regarding the cross section, dσ(E ν , E e )/dE e , we use the one found in [46].
For the JUNO detector, we assume the same proton fraction 11% as the Daya Bay detectors [47] which implies ∼ 1.44 × 10 33 free protons for 20 kt. The differential flux of reactor neutrino dφ(E)/dE can be computed as, where P th is the thermal power of the reactor, E 210 MeV is the average energy released by per fission computed by taking into account the ratios of the fuel compositions of the reactor (see below).
We can replace, in a good approximation, the reactor complex consisting of 6 and 4 reactors, respectively, at Yangjiang and Taishan sites by a single reactor with the thermal power of 35.8 GW placed at the baseline L = 52.5 km from the JUNO detector. We also include the contributions from the far reactor complexes at Daya Bay (with the baseline of 215 km) and Huizhou (with the baseline of 265 km) sites, which contribute, respectively, about 3% and 2% in terms of the total number of events.
For the reactor spectra S(E) we use the convenient analytic expressions found in [38] with the typical fuel compositions of the reactors, 235 U: 239 Pu: 238 U: 241 Pu = 0.59: 0.28: 0.07: 0.06, obtained by taking time average of Fig. 21 of [48]. We note that S(E) is nothing but the number of neutrinos being emitted per fission per energy (MeV).
Furthermore, we also include, for completeness, geoneutrinos coming from the decays of U and Th inside the earth in the same way as done in [31], despite that it is not important for our main purpose. Assuming the input (true) geoneutrinos fluxes of φ(U) = 4.0 × 10 6 cm −2 s −1 and φ(Th) = 3.7 × 10 6 cm −2 s −1 , the expected U and Th geoneutrinos induced events are, respectively, ∼ 1.9 × 10 3 and ∼ 5.4 × 10 2 for 5 years of exposure at JUNO detector. For simplicity, we consider only the averaged standard oscillation for geoneutrinos, and ignore the decay effect for them. Since the presence of geoneutrinos has minor impact on the standard oscillation study at JUNO, their decay effect is even minor and we believe that this is a fairly good approximation. In fact we have verified that the presence of geoneutrinos has virtually no impact on the sensitivity to the decay effect. In the fit, we let the fluxes of geoneutrinos vary freely subject to the pull terms in Eq. (4.5) with σ ξ U = σ ξ Th = 20%.
R(E e , E vis ) is the function which takes into account the finite energy resolution of the detector and is given by where the energy resolution is assumed to be [49], where η is introduced to take into account the uncertainty in the energy resolution. The energy resolution in (B.4) amount to consider only the stochastic term. The expected total number of events at JUNO for the 5 years of exposure with 100% detection efficiency is 1.41 × 10 5 . The individual contributions from medium-baseline reactor sites (at Yangjiang and Taishan), far reactor sites (at Daya Bay and Huizhou), geoneutrinos are, respectively, 1.39×10 5 , 6.70×10 3 and 2.40×10 3 .

C Correlations among the oscillation, systematic and decay parameters
Here, we discuss the correlations among the mass and the mixing parameters, systematic uncertainty parameters, and the decay parameter τ 3 /m 3 . For this purpose we have examined all possible correlations among the mixing parameters and the uncertainty parameters. 12 The calculation is done under the same conditions (running time etc.) as assumed for calculating the lifetime bound given in Fig. 2.
The general features of the correlations among the oscillation parameters, systematic uncertainty parameters and the decay parameter (m 3 /τ 3 ) can be summarized as follows: • The effect of ν 3 decay is visible only in the contours which involve the energy resolution σ E (1 + η), θ 13 and ∆m 2 31 , as implied by the results shown in Table 1. 13 • The correlations between τ 3 /m 3 and σ E (1+η), θ 13 exist, but in easily understandable way.
There arises an extended allowed region up to sin 2 θ13 0.044 at 2σ C.L. with a larger flux normalization by about 5%. In Fig. 5, we show the correlations in space spanned by (a) sin 2 θ 13 -τ 3 /m 3 (left panel) and (b) σ E (1 + η) -τ 3 /m 3 (right panel). We notice that there are correlations, i.e., the effect of decay, both in the left and right panels, but only in near the lower end of the allowed region of the lifetime. From the left panel of Fig. 5 we learn that the decay effect (diminishing the amplitude of the wiggles) can be compensated to some extent by making θ 13 larger, which produces a new allowed region toward the large θ 13 direction.
Corresponding to the newly emerged region, there also arises a new allowed region in the σ E (1 + η) -τ 3 /m 3 space, as seen in the right panel of Fig. 5. It is located at near the top of the peninsula, a slightly distorted region toward small σ E (1 + η) direction in the right panel of Fig. 5. The shape reflect a better energy resolution in the newly allowed region. We note that due to the above sin 2 θ 13 -σ E (1 + η) correlation, there always exist a slightly extended allowed region in any one of the correlation plots which involve sin 2 θ 13 or σ E (1 + η) toward the direction of larger θ 13 and smaller σ E (1 + η).

C.2 Correlations of parameters: miscellaneous cases
We now show, for completeness, the allowed regions spanned by all the remaining 32 combinations of parameters considered in this work (which are described explicitly in the footnote 12), except for those already shown in Figs. 4 and 5. The meanings of the filled colors and line types are the same as in Figs. 4 and 5.
As we can see from these plots that, in general, there are no significant differences among the allowed regions shown by the filled colors (no decay in input but allowed in the fit), solid black curve (standard oscillation fit without decay) and black dashed curves (decay effect both in the input and in the fit). Therefore, as a whole, the impact of the decay is rather small. Furthermore, we note that there is no significant newly induced correlations due to the decay effect.