Invisible neutrino decay : First vs second oscillation maximum

We study the physics potential of the long-baseline experiments T2HK, T2HKK and ESS$\nu$SB in the context of invisible neutrino decay. We consider normal mass ordering and assume that the state $\nu_{3}$ as unstable, decaying into sterile states during the flight and obtain constraints on the neutrino decay lifetime ($\tau_3$). We find that T2HK, T2HKK and ESS$\nu$SB are sensitive to the decay-rate of $\nu_{3}$ for $\tau_{3}/m_{3} \leq 2.72\times10^{-11}$s/eV, $\tau_{3}/m_{3} \leq 4.36\times10^{-11}$s/eV and $\tau_{3}/m_{3} \leq 2.43\times10^{-11}$s/eV respectively at 3$\sigma$ C.L. We compare and contrast the sensitivities of the three experiments and specially investigate the role played by the mixing angle $\theta_{23}$. It is seen that for experiments with flux peak near the second oscillation maxima, the poorer sensitivity to $\theta_{23}$ results in weaker constraints on the decay lifetime. Although, T2HKK has one detector close to the second oscillation maxima, having another detector at the first oscillation maxima results in superior sensitivity to decay. In addition, we find a synergy between the two baselines of the T2HKK experiment which helps in giving a better sensitivity for $\theta_{23}$ in the higher octant. We discuss the octant sensitivity in presence of decay and show that there is an enhancement in sensitivity which occurs due to the contribution from the survival probability $P_{\mu\mu}$ which is more pronounced for the experiments at the second oscillation maxima. We also obtain the combined sensitivity of T2HK+ESS$\nu$SB and T2HKK+ESS$\nu$SB as $\tau_{3}/m_{3} \leq 4.36\times10^{-11}$s/eV and $\tau_{3}/m_{3} \leq 5.53\times10^{-11}$s/eV respectively at 3$\sigma$ C.L.


Introduction
Neutrinos are one of the most fascinating particles in nature. Observation of oscillation of neutrinos signifies evidence of neutrino mass. The very existence of neutrino mass is the reason of neutrinos being so special since in the Standard Model (SM) of particle physics there is no explanation of neutrino mass. The standard formalism of neutrino oscillation requires three mixing angles θ 12 , θ 13 and θ 23 , two mass-squared differences ∆m 2 21 and ∆m 2 31 and one CP phase δ CP . The parameters θ 12 , θ 13 , ∆m 2 21 and the absolute value of ∆m 2 31 have been measured very precisely by oscillation experiments. However, the octant of θ 23 , the sign of ∆m 2 31 (mass-ordering) and δ CP are not yet determined with certainty. Global analysis results [1][2][3] indicate towards normal hierarchy of ∆m 2 31 and θ 23 in the higher octant. The global analysis also disfavours δ CP = π 2 with more than 3σ CL. Many experiments are planned for the near future with the goal of determining these quantities. Some examples of future experiments are DUNE [4][5][6][7][8], T2HK/T2HKK [9,10], ESSνSB [11], JUNO [12], INO [13], PINGU [14], KM3Net-ORCA [15] etc. It is a well known fact that if there is a new physics beyond the SM, then that can result in modification of the oscillation probabilities in these experiments. So, these experiments can be sensitive to the new physics. Also, the new physics can affect the sensitivity of these experiments to the standard parameters. Invisible neutrino decay during the flight is one such new physics idea.
The idea of neutrino decay was first proposed to explain the solar neutrino problem in the very early days [16]. Later, neutrino oscillation with decay solutions were studied as an explanation of the depletion of solar neutrinos [17][18][19][20][21][22][23][24]. These assumed ν 2 to be the unstable state and were able to put bound on the lifetime of ν 2 . The bound from the solar neutrino data is τ 2 /m 2 > 8.5 × 10 −7 s/eV [23]. A recent study analysed low energy solar neutrino data and put bounds on both τ 2 and τ 1 [25]. Supernova neutrinos also give bounds on τ 1 and τ 2 . SN1987A data gives the bound of τ /m > 10 5 s/eV [26]. τ 1 and τ 2 can also be constrained from high resolution multi-ton Xenon detector. Recently, in ref. [27], the authors show that such a detector can give very strong bounds τ 1 /m 1 > ∼ 3 × 10 −2 s/eV and τ 2 /m 2 > ∼ 8 × 10 −3 s/eV at 2σ level using solar neutrinos. Atmospheric and long-baseline experiments give bounds on the ν 3 lifetime. A neutrino decay solution (without any oscillation) was proposed in ref. [28] to explain the atmospheric neutrino problem but this solution fitted the data poorly. Ref. [29,30] considered neutrino decay and neutrino mixing together. This was successful to reproduce the L/E distribution of the Super-Kamiokande (SK) data. However, when zenith angle dependence was used instead of L/E distribution, this model was found to give poorer fit to the SK data [31]. Ref. [29,31] assumed ∆m 2 > 0.1 eV 2 to comply with the K-decay bounds [29]. Therefore ∆m 2 dependent terms were averaged out. These constraints can be relaxed if the unstable state decays to some invisible state with which it has no mixing. There are two scenarios in the literature. The first kept ∆m 2 unconstrained and it explicitly appeared in the probabilities [32]. This fits the SK data with a best-fit value of non-zero decay parameter and ∆m 2 ∼ 0.003 eV 2 . Ref. [33] considered ∆m 2 << 10 −4 eV 2 . Here, the probability does not contain ∆m 2 explicitly. This was able to fit the SK data. However, independent analysis by SK collaboration showed that this scenario gives a poorer fit to the data than only oscillation [34]. Global analysis of atmospheric and long-baseline experiments were performed in ref. [35]. Only oscillation gave best-fit to the SK data and the fit for the oscillation plus decay was not bad. But addition of LBL data from MINOS reduced the fit quality. This analysis put bound on τ 3 /m 3 ≥ 2.9 × 10 −10 s/eV at the 90 % C.L. Ref. [36] studied the oscillation plus decay scenario with unconstrained ∆m 2 for MINOS and T2K data and it found τ 3 /m 3 > 2.8 × 10 −12 s/eV at 90 % C.L. Most of these analyses are done using two generation approximation without matter effect. Authors of ref. [37] performed a complete three generation study of the oscillation plus decay scenario assuming matter effect for the NOvA and T2K preliminary data. This study put a bound on the lifetime on the τ 3 /m 3 as τ 3 /m 3 ≥ 1.50 × 10 −12 s/eV at 3σ. Recently, ref. [38] addressed the IceCUBE track and cascade tension using invisible neutrino decay. They showed that an unstable neutrino with τ /m = 100 s/eV solves the track vs cascade tension.
There have been many studies in literature which discuss the potential of future experiments to the invisible decay. Medium baseline reactor neutrino experiment like JUNO can give a bound on τ 3 /m 3 > 7.5 × 10 −11 s/eV (95 % C.L.) [39]. Decay of ultrahigh energy astrophysical neutrinos can give constraints on decay [40][41][42][43]. Recently, ref [44] showed that IceCUBE can probe τ /m upto 10 s/eV for both mass-orderings for 100 TeV neutrinos coming from a source at a distance of 1 Gpc. The ref. [45] studied the invisible neutrino decay in the context of DUNE. It showed that using charged-current electron type and muon type events, DUNE after 5+5 years of running can put a bound τ 3 /m 3 > 4.50 × 10 −11 s/eV at 90 % confidence level. More recently ref. [46] performed a multi-channel analysis using charged-current, neutral current and tau-channel events and showed that DUNE's sensitivity would be τ 3 /m 3 > 5.2 × 10 −11 s/eV at 90 % confidence level. For the potential of MOMENT experiment to probe invisible neutrino decay see ref [47]. There are also studies on the expected sensitivities from the future atmospheric experiments. See ref. [48,49] for INO and ref. [50] for KM3Net-ORCA.
Invisible neutrino decay can happen for both Dirac and Majorana neutrinos. If neutrinos are Dirac, there can be a coupling between the neutrinos and a light scalar boson [17,51]. This gives the decay channel ν j →ν iR + χ, whereν iR is a right-handed singlet and χ is an iso-singlet scalar. If neutrinos are Majorana particles, neutrino can couple with a Majoron and a sterile neutrino via a pseudo-scalar coupling [52,53]. This gives ν j → ν s + J. LEP data on the Z-decay to invisible particles constraints the Majoron to be dominantly singlet [54]. Neutrinos can also decay to another active state [55][56][57].This is called the visible neutrino decay scenario. This type of decay can happen in the following way. ν j →ν i + J or ν j → ν i + J. If neutrinos are Majorana, the decay product can be observed in the detector. Ref. [58][59][60] discusses the visible decay in the context of longbaseline experiments. In ref. [61], the authors discuss the visible neutrino decay for reactor experiments KAMLAND and JUNO. For visible decay of the astrophysical neutrinos at IceCUBE, see ref [62].
In this paper we study the constraints on invisible neutrino decay which can come from future planned/proposed long baseline experiments -T2HK/T2HKK [9,10] and ESSνSB [11]. We perform a full three flavour study using matter effect and obtain the sensitivity to τ 3 /m 3 for these experiments. The salient feature of the T2HKK and ESSνSB experiment is that they are both designed to have energy peak near the second-oscillation maximum. Since the second oscillation maximum occurs at a lower energy for a particular baseline the effect of decay is expected to be more at the second oscillation maximum. We examine this aspect and delve into the detail of whether the experiments at the second oscillation maximum stand to gain in sensitivity in presence of decay. We also check if the determination of θ 23 can get affected if we assume decay in the data, while the fit does not not assume any decay to be present. In particular, we investigate how the experiments at first and second oscillation maximum fare in this respect and what are the important factors on which the measurement of θ 23 can depend in presence of decay. We also explore how the octant sensitivity of these experiments change in presence of neutrino decay.
The paper is organized in the following way. The next section discusses the neutrino oscillation in presence of invisible neutrino decay. In Section 3 we give our experimental and numerical details, in section 4 we present our results and in section 5, we finally draw our conclusions.

Neutrino oscillation in presence of invisible decay of neutrino
In this section we discuss the propagation of neutrinos in presence of invisible neutrino decay. We assume that the ν 3 is unstable and it decays into a sterile neutrino and a singlet scalar (ν 3 →ν 4 + J) with lifetime τ 3 . In this case we can extend the mass and flavour bases by ν i ν 4 T and ν α ν s T , where i = 1, 2, 3 and α = ν e , ν µ , ν τ . They are related by the following unitary relation.
U is the standard PMNS matrix describing the standard 3 neutrino oscillation. We assume normal hierarchy and m 4 to be the least massive state. We assume that the decay eigenstates and the mass-eigenstates are same. Under these assumptions, we can write the neutrino evolution in presence of matter in the following way. and Here, A is the matter potential, G F is the Fermi constant, E is the energy and n e density of electrons in the earth. We define α 3 = m 3 /τ 3 as the decay rate of the ν 3 state. Since we are considering the decay of ν 3 only, from here onward, we use α instead of α 3 for the ν 3 state. The probability of getting a neutrino in the flavour state ν b for the initial flavour state of ν a is given by Here, a,b denote the flavour states e,µ, τ . The effect of decay comes as the exp(−αL/E) factor in the probability. So an experiment is sensitive to the values of α where, α ∼ E/L. It is clear that, smaller values of α and hence a higher sensitivity can be achieved for longer baseline and for a particular baseline the sensitivity is more for a lower energy. Thus sensitivity is expected to be more at the second oscillation maxima as compared to the first oscillation maxima.

Experimental and simulation details
In this section we describe various experiments and the specifications used in our analysis with the summary of the experiments presented in Tab. 1. First we give the brief descriptions of the experiments followed by details of our numerical simulations.

T2HK
T2HK [9] is a proposed upgradation plan of the currently running T2K experiment in Japan. The neutrinos will be generated at Tokai by an upgraded version of the J-PARC beam. Currently J-PARC gives a beam power of 470 KW, but before T2HK becomes operational, the beam power will be increased beyond 1.3 MW. Although near detector of T2HK is yet to be finalized, there are several ideas. Some of these are like upgradation of the current near detector ND280, building a water Čerenkov detector similar to the far detector but in a smaller scale. The far detector will be the upgradation of the currently running Super-Kamiokande (SK) [64] detector. SK is situated at 295 km away from Tokai at the Kamioka village and it is slightly (2.5 • ) off-axis to the beam-axis. Thus this will give a narrow beam centered around 0.56 GeV which is at the first oscillation maximum of the neutrino oscillation. The upgraded SK, called Hyper-Kamiokande (HK), will consist of two identical cylindrical tanks placed upright. Each tank will contain pure water which will have a fiducial mass of 187 kt each thus giving a total fiducial mass of 374 kt. The new detector will also have better resolution and efficiency compared to SK.

T2HKK
One of the main aim of the future long-baseline experiments is to measure the CP-violation, and T2HK is no exception. T2HK has a very good sensitivity to the CP-violation if the mass-hierarchy is known. However, the CP-violation sensitivity of T2HK drops if masshierarchy is not known beforehand. This is due to a degeneracy between mass-hierarchy and δ CP [65]. In order to solve this problem, there is a proposal to shift one of the water tanks of HK to Korea at about a distance of 1100 km away from the source at Tokai. This proposal is called the T2HKK [10] (the second K is for Korea.). This will have more matter effect. Thus the mass-hierarchy sensitivity will increase significantly breaking the CP-hierarchy degeneracy. However, the second detector will be at the second oscillation maximum and as the baseline is larger and the flux will also decrease.
The Korean detector site is not yet final and there are many possibilities [10]. However, all these sites are in the southern part of the Korean peninsula and lie within a range of 1-3 • off-axis angle with the J-PARC beam line.

ESSνB
ESSνB [11] is another future long-baseline super-beam experiment. The neutrino beam will be generated from the ESS facility at Lund, Sweden. This will have a 2 GeV proton beam with 5 MW beam power. The neutrino beam created by this will have a peak around 0.25 GeV. The far detector will be situated at a distance of about 540 km away from Lund at a mine in Garpenberg. The far detector is proposed to be 500 kt water Čerenkov detector. At 540 km, the second oscillation maximum occurs at 0.35 GeV. Thus the peak of the neutrino beam will be close to the second oscillation maximum for ESSνB.

Simulation details
We have used Global Long-Baseline Experiments Simulator (GLoBES) [66,67] for simulating the long-baseline experiments. We present our results in terms of statistical χ 2 , where, corresponds to the simulated data and (N ) test i is the number of events predicted by the theoretical model. The effect of systematic errors are included by the "pull" method through the "pull" variables ξ. We have incorporated signal normalization error, background normalization error, signal "tilt" error & background "tilt" error in our analysis. Incorporating the errors, the signal and background events can be written as, where, s(b) denotes signal(background). c norm i (c i tilt ) denotes the change in number of events by the variation of the "pull" variable ξ norm (ξ tilt ). The signal and background normalization errors are presented in Tab. 2. In the above equation E i represents the mean reconstructed energy of the i th bin. The maximum and minimum energy of the energy range are E min and E max respectively. The mean energy is given by,Ē = (E max + E min )/2.
Channel T2HK(295 km) T2HK(1100 km) ESSνSB The detector type proposed for T2HK, T2HKK and ESSνSB are water Čerenkov detectors. T2HK & T2HKK have been simulated based on [10] and the ESSνSB experiment has been simulated according to [11]. The detector of ESSνSB is based on the MEMPHYS detector [68], which is a water Čerenkov detector. The neutrino (antineutrino) appearance channels for signal events are ν e (ν e ) and the background events in this channel are neutral current background, intrinsic beam background, mis-identified muons and wrong-sign signal errors. The neutrino (antineutrino) disappearance channels are ν µ (ν µ ), the background channels are neutral current background and wrong-sign muons. The signal and background normalization uncertainties have been considered as presented in Tab. 2. Additionally, a 5% (10%) signal (background) "tilt" error have been considered for T2HK & T2HKK and 0.1% (0.1%) in case of ESSνSB.

Oscillation parameters
True Values Marginalization Range For the calculation of probability in presence of invisible neutrino decay we modified the probability code of the GLoBES. We used θ 12 = 33.82 • , θ 13 = 8.6 • , sin 2 θ 23 = 0.447(0.56) for lower (higher) octant, δ CP = −90 • ,∆m 2 21 = 7.39 × 10 −5 eV 2 and ∆m 2 31 = 2.52 × 10 −3 eV 2 for generating the simulated data (true values) as shown in the Tab. 3. These are consistent with the current ranges allowed by the latest global fits [1]. For the statistical studies we have marginalized over δ CP and θ 23 in their allowed ranges. We kept other standard parameters fixed as we found they have little effects on the marginalization. We also assumed the mass-hierarchy is known and it is normal hierarchy.

Results
In this section, we present the results of our study. First, we give the plots at the probability level in presence of decay for T2HK, T2HKK and ESSνB. These figures can be helpful in understanding many of the results presented in the next sections. Then we give the sensitivity to the decay rate α for each experiment. We also discuss a particular synergy which we observed between the two detectors in T2HKK experiment -one at the first oscillation maximum and one at the second oscillation maximum. We also present the results on the effect of decay on the measurement of θ 23 and octant sensitivity for these experiments

Probability at T2HK/T2HKK and ESSνSB baselines
In fig.1, we present the appearance and disappearance probabilities as a function of energy for the three different baselines under investigation. We give the plots for the no-decay case (α = 0) and also for a representative value of α = 3 × 10 −5 eV 2 . The choice of this α is motivated from the current bounds as obtained from T2K and NOνA data analysis [37]. In both cases we give the plots for two values of θ 23 -one in lower octant (40 • ) and one in the higher octant (51 • ). The probabilities for other values of θ 23 within this range will lie between these two curves. The figures show the interplay between θ 23 and α for the various baselines. The left (right) panel shows ν µ → ν e (ν µ → ν µ ) probabilities. The black (blue) solid and red (green) dashed curves are for α = 0 (α = 3 × 10 −5 eV 2 ) and θ 23 in the lower octant and higher octant respectively. The top, middle and the bottom panel are for the baselines 295 km, 1100 km and 540 km respectively.
In the left panels of fig. 1, we see that for all the baselines, at the oscillation maxima, P µe for the the no decay case is higher for a fixed value of θ 23 . As decay is introduced, the probabilities reduce. On the other hand, for a fixed value of α, lower octant gives smaller probabilities compared to the higher octant. Thus the peak appearance probability can get reduced due to increase in the value of α and/or decrease in the value of the mixing angle θ 23 .
The top panel is for T2HK and here the first oscillation maximum is at ∼0.6 GeV which is also where the flux peaks. For the baseline and energies involved the effect of decay for T2HK is small for the sample value of α considered in the plot. However, we see that the higher and lower octant bands are well separated for T2HK.
Next we focus on the middle and the bottom panels which are relevant baselines for T2HKK (second detector) and ESSνSB. For these experiments the flux peaks near the second oscillation maxima -∼0.7 GeV for the baselines 1100 km and ∼0.35 GeV 540 km. Since these are higher baselines and/or lower energies the effect of decay is more pronounced. The HO and LO bands are much closer in these cases and a reduction in probability due to non-zero values of α can be compensated by increasing θ 23 .
In the right panels of fig. 1 we show the ν µ → ν µ disappearance probability P µµ which plays a crucial role in determining the value of θ 23 in long baseline experiments. In this case, the curves with α = 0 are overlapping irrespective of the octant of θ 23 which shows that the P µµ channel do not have any octant sensitivity in absence of neutrino decay, being dependent on sin 2 2θ 23 . So, the red-dashed and black-solid curves (no-decay) are well separated from the blue-solid and green-dashed curves which reveals the effect of decay. This feature is true for all the three baselines considered in our study. In case of the 295 km baseline the blue curve is closer to the black and red curves compared to the 1100 km and 540 km cases since the effect of decay is less for the 295 km baseline for the α value chosen. Also, it has to be noted that, the otherwise octant degenerate P µµ channel begins to have octant sensitivity once α = 0, resulting in a gap between the blue and green curves.
This can be simply understood if we consider the expression for two-generation survival probability in vacuum In absence of decay the e − αL E term is 1 and the probability depends on sin 2 2θ 23 . But in presence of decay octant sensitivity ensues due to the e − αL E factor. For our later analysis, it is important to understand the dependence of the probability on θ 23 at the survival probability maxima (SPMAX) and minima (SPMIN). At SPMIN, though the probability value is small, the flux peaks near this energy since it corresponds to the oscillation maxima. In figure 2 we show the behaviour of the probability for the 540 km baseline as a function of θ 23 for energies corresponding to SPMAX and SPMIN. Similar behaviour is also true for the two other baselines -295 km and 1100 km. It is seen that at SPMAX, the no-decay probability is always higher than the decay probability. The no-decay probability has no dependence on θ 23 while the decay probability reduces with increasing θ 23 monotonically. This can be easily understood from the expression 4.1. At SPMAX sin 2 ∆m 2 31 L 4E = 0 and the remaining term decreases monotonically with increasing θ 23 for non-zero α. At SPMIN, the θ 23 dependence is more complicated. For lower octant the decay probability is larger as compared to the no-decay probability and the probability reduces with increasing θ 23 till a certain value (depending on α), after which in higher octant it increases with θ 23 . For higher octant, no-decay probability is higher than the decay probability. As we will see, these features of the P µµ probability will play a role in the χ 2 fit.

Sensitivity to the decay
In this section we study the capability of the T2HK, T2HKK and ESSνB experiments in constraining invisible neutrino decay. We assume the neutrinos to be stable while generating the simulated data. In the fit we have marginalized over |∆m 2 31 , θ 23 and δ CP in the range given in Tab. 3. The fig. 3 shows the value of χ 2 for different values of α taken in the fit (test). A higher value of χ 2 is more disfavoured compared to a lower values of χ 2 .
In the left panel we assumed sin 2 θ 23 (true) = 0.447 (lower octant) while generating the data. In the right panel we assumed sin 2 θ 23 (true) = 0.56 (higher octant). The red solid, blue dashed dotted and black dashed curves are for T2HK, T2HKK and ESSνB respectively. Naively, the sensitivity to decay depends on the L/E. From Tab. 1, we see that the second detector of T2HKK has highest L/E and the detector at 295 km has lowest L/E (See Tab. 1). So, we expect T2HKK to have the best sensitivity and T2HK to have the worst, but we observe that ESSνSB has a lower sensitivity than T2HK. Here, marginalization over θ 23 plays an important role. Determination of θ 23 is governed mainly by the P µµ channel. In absence of decay there is no sensitivity to θ 23 at the maxima of survival probability which corresponds to sin 2 ∆m 2 31 L/4E = 0. The maximum precision of θ 23 comes at the oscillation maximum or the survival probability minimum (SPMIN) where, sin 2 ∆m 2 31 L/4E = 1 and can be expressed as, In presence of decay, the survival probability maxima also acquires some sensitivity to θ 23 which increases with increasing value of the decay constant α for a fixed baseline. As seen in the earlier section, for T2HK the octant bands are well separated while the effect of decay is not very significant. Hence marginalization over θ 23 does not play much role. For ESSνSB on the other hand, the octant bands are not so well separated and the effect of decay is more. As a result, marginalization over θ 23 tends to make the probabilities closer by shifting θ 23 to a different value and the sensitivity reduces. The direction of shift depends on the initial value of θ 23 . For SPMIN, decay in the fit, can give same probability as no-decay in data for an increased θ 23 for both lower and higher octant. This can be seen by drawing a horizontal line in fig. 2. For SPMAX, on the other hand, if we fit decay with no-decay in data, then shift of θ 23 towards lower values make the probabilities closer, giving a lower χ 2 . Thus, depending on which energy bins contribute most, the θ 23 in the fit comes at a higher or lower value, reducing the χ 2 . We have checked that if we keep θ 23 as fixed then ESSνSB gives a better sensitivity.
In case of T2HKK, other than a detector close to the second oscillation maximum, there is also a detector at the first oscillation maximum which can measure θ 23 with a good precision. So, the combination of the two detectors, one at the first oscillation maximum and other at the second oscillation maximum gives the best sensitivity. This will be discussed further in the next section. There is also another interesting feature for ESSνSB, which is that for true octant as higher octant we observe a sudden fall of sensitivity for α above α ∼ 3 × 10 −5 eV 2 . We will see in the next section that this is a unique feature of second oscillation maxima and we will discuss this in detail in the coming section. In this section we discuss a synergy between the two baselines L1 = 295 km and L2 = 1100 km of T2HKK in constraining neutrino decay. The first one corresponds to the Japanese detector (JD) and the second one corresponds to the Korean detector (KD). It can happen that the χ 2 obtained by combining two or more experiments can be different from the naive sum of their individual χ 2 s. We call such enhancement of χ 2 as synergy between the experiments. In fig. 4, we present the sensitivity to the α for one of the HK detectors situated at L1 (dark red solid) and L2 (cyan dashed). We also show their combined sensitivity ( blue dashed-dotted) which is nothing but the sensitivity of T2HKK. The grey dotted curve shows the naive summation of the χ 2 s i.e χ 2 L1 + χ 2 L2 . The left panel is for sin 2 θ 23 = 0.447 and the right panel is for sin 2 θ 23 = 0.56. We see that although there is no significant synergy for the left panel, there is some significant synergy for the right panel for α greater than 2.268 × 10 −5 eV 2 , i.e., the χ 2 L1+L2 of T2HKK is larger than the

Synergy between 295 km and 1100 km baselines
This synergy can be understood from the fig. 5 and fig. 6. These figures give the χ 2 for a given value of θ 23 in the fit with a fixed value of α in the fit. The decay constant α = 0 while simulating the data and |∆m 2 31 | and δ CP are marginalized for all the cases. Therefore global minima of the curves in these figures give the χ 2 values for the corresponding α (test) in fig. 4. Fig. 5 is for the case when sin 2 θ 23 = 0.447 and fig.6 is for sin 2 θ 23 = 0.56 in the simulated data. The dark red solid curves assume α = 0 in the fit and green dashed curves assume α = 3 × 10 −5 eV 2 in the fit. First let us consider the fig. 5. Here, sin 2 θ 23 (true) = 0.447. We see that for L1, the global minimum, for a test value of α = 3 × 10 −5 eV 2 in the fit, shifts slightly towards a higher sin 2 θ 23 but stays in the same octant as the α = 0. The shift towards a higher value can be understood by looking at the curve for P µµ . It is seen that at the oscillation maxima (i.e SPMIN), the no-decay probability (black) curve is lower than that with decay (blue curve) for lower octant. Hence when we try to fit data with no-decay with decay in theory, the probabilities can come closer by increasing θ 23 slightly. Similar feature is also observed at L2 in the right panel though. the shift is almost un-noticeable. Sensitivity to the parameter α for the combination L1 and L2 will be given by the global minimum of the sum of χ 2 s of two green dashed curves of fig. 5. As the global minima of these two curves approximately coincide (42 • for L2 and 42.5 • for L1), the global minimum of the sum of their χ 2 s will be approximately equal to the sum of χ 2 s at the global minima of L1 and L2. Therefore, the sensitivity to the combination of L1 and L2 will be approximately equal to the sum of the sensitivities for L1 and L2. Thus there is not significant synergy for this case.
Next we consider the fig. 6. Here, sin 2 θ 23 = 0.56. We see that for L1, like previous case, the global minimum for α = 3 × 10 −5 eV 2 (test) shifts a little towards higher value since at SPMIN, no-decay probability is larger, it can be matched by increased value of θ 23 (cf. ??).
However, for L2, we see that the position of the global minimum drastically changes and now it is in the opposite octant, i.e., lower octant. Again, like the previous case, the sensitivity to the combination of L1 and L2 will be determined by the sum of χ 2 s of two green dashed curves. But here, the positions of the global minima for L1 and L2 are completely different. Therefore, the position of the global minimum of the sum of χ 2 s of L1 and L2 curves is non-trivial and is not equal to the sum of the χ 2 s at the global minima of L1 and L2 curves. Hence, we see the synergy in this case. The shift of the global minimum can be attributed to the interplay between the appearance and disappearance channels.
In order to understand this interplay in fig. 7 we have plotted the appearance and disappearance χ 2 s for two different values of the decay constant α for L2. From the left panel of the figure it is seen that for true sin 2 θ 23 in the lower octant, the minima for both appearance and disappearance channels stay in the correct octant. The octant sensitivity is seen to come from both disappearance and appearance channel. However, for appearance channels, the octant bands are not very widely separated for experiments close to second oscillation maxima, and with increase in α the effect of disappearance channel becomes more pronounced and there is a sharp rise in the octant sensitivity coming from the disappearance channel.
The right panel depicts the situation for the higher octant. In this case when the nodecay in data is fitted with decay, the minima for the appearance channel stays in the same octant, shifting slightly towards higher θ 23 . This behaviour can be understood from the appearance probability since for this HO and no-decay gives a higher probability and hence if we fit with decay, the θ 23 tends to shift to a larger value to get closer to the data. As the  Figure 7.
The χ 2 as a function of θ 23 . The left (right) panels are for true sin 2 θ 23 = 0.447 (sin 2 θ 23 = 0.56) respectively. The red and blue curve are for appearance and disappearance channels respectively. While the total contribution from both the channels is represented by green curves. The solid curves are for α(test) = 3 × 10 −5 eV 2 and the dashed curves for α(test) = 10 −5 eV 2 in the fit. The data have been generated for α = 0. decay constant increases, the shift in the value of θ 23 to get closer to the no-decay curve will be higher.
On the other hand the disappearance minima is seen to shift to the lower octant. For a lower value of the decay constant, the difference between the disappearance χ 2 minimum in the two octants is not very high. On the other hand, the appearance χ 2 rises steeply in the opposite octant and when the appearance χ 2 is added to the disappearance χ 2 , the overall minima comes in the correct octant. But for a higher value of the decay constant, the minima in the wrong octant is much lower than the minima in the correct octant and thus the overall minima comes in the wrong octant, being driven by disappearance channel. The reason for disappearance channel preferring lower octant can be understood from the probability figure 1, by noticing that for peak flux energies of around 0.6 GeV, the redcurves (no-decay, HO) are closer to the blue curves (decay, lower octant). Hence, when data is generated in the higher octant with no decay, disappearance channel prefers the lower octant if α = 0 in the fit, the shift being higher as the decay constant increases. This effect is not present at the detector tuned for the first oscillation maximum at 295 km. From the top-left panel of Fig. 1, we see that the octant bands are wide enough such that the appearance channel gives high octant sensitivity and therefore the appearance channel always dominates and we do not get any false octant solution.
The above discussion can be now used to explain the kink in the sensitivity of ESSνSB ( fig. 3 right panel). For ESSνSB there is only a single detector. For lower values of α, the disappearance channel is still weak and the appearance channel dominates and we get the octant in the higher octant. However, as the α becomes more than some critical value, the disappearance channel begins to dominate and the octant flips to the wrong side. Thus the χ 2 is abruptly increased, giving a lower sensitivity. Figure 8. The ∆χ 2 as a function of θ 23 (test) assuming sin 2 θ 23 (true) = 0.447. The left, middle and the right panels are for T2HK, T2HKK and ESSνB respectively. The dark red solid curves are for standard cases and the green dashed curves are for the cases where α is assumed to be 3 × 10 −5 eV 2 in the simulated data but decay is not considered in the fit.

Role of decay in the θ 23 determination
In this section, we will see how presence of invisible neutrino decay can affect the measurement of θ 23 in T2HK, T2HKK and ESSνB experiments. Fig. 8 gives ∆χ 2 as a function of θ 23 (test). Here, we assumedsin 2 θ 23 = 0.447 in the data. The left, middle and right panels are for T2HK, T2HKK and ESSνB respectively. In the fit, we marginalized over, |∆m 2 31 |, δ CP and kept α fixed at zero. For dark-red solid curves, we assumed stable neutrino in the data and for green dashed curves, we assumed α = 3 × 10 −5 eV 2 in the data.
We see that for the three cases, the best-fit values for data generated for α = 3 × 10 −5 eV 2 and α = 0 are different. We notice that for true α = 3 × 10 −5 eV 2 the best-fit of θ 23 is shifted towards lower values for T2HK and T2HKK. However, for ESSνSB, the shift is in the opposite direction. The shift towards lower values for T2HK and T2HKK is governed by the behaviour of P µµ at the oscillation maxima i.e at SPMIN, where the flux peaks. The data is generated for decay which gives a higher probability than no-decay in the lower octant. Therefore when data is fitted with α = 0, a reduced value of θ 23 gives a better fit since the probability increases with decreasing θ 23 .
The reason for different behaviour for ESSνSB is more complicated. Here, resolution plays a big role. The resolution of the detector smears the imprint of probability over all energy bins. As a result, the feature near the oscillation minima is lost due to the larger bin width of 0.1 GeV. A bin from 0.3 GeV to 0.4 GeV contains events from the survival probability minima as well as the features of the probability around 0.3 GeV where the probabilities corresponding to decay and no decay curves show opposite behaviour than at ∼ 0.35 GeV. As an upshot if one plots the events then the no-decay events are higher than the events with decay in most bins. When the simulated event spectrum is fitted without decay, θ 23 is shifted towards higher value to bring the probability down. The effect of smearing affects the experiments close to second oscillation maxima more since the variation with probability is sharper compared to the first oscillation maxima. Fig. 9 is similar to fig. 8, but here, we assume sin 2 θ 23 = 0.56 in the true data. We see Figure 9. The ∆χ 2 as a function of θ 23 (test) assuming sin 2 θ 23 (true) = 0.56. The left, middle and the right panels are for T2HK, T2HKK and ESSνB respectively. The dark red solid curves are for standard cases and the green dashed curves are for the cases where α is assumed to be 3 × 10 −5 eV 2 in the simulated data but decay is not considered in the fit.
like the previous figure, the best-fit values of θ 23 is changed but in this case, for all three experiments, the shifts of θ 23 are towards a lower value. This can be explained since, for θ 23 in the higher octant the probability for non-zero α is smaller than the probability with α = 0 for both SPMIN and SPMAX and hence no-decay gives larger number of events. For θ 23 in higher-octant, lowering θ 23 reduces the probability and hence gives a better fit when we fit decay in data with no-decay. In Fig. 10, we show the octant sensitivities as a function of the CP phase δ CP for experiments which are tuned to the second oscillation maxima in presence and absence of decay. The top panels are for ESSνSB and the bottom panels are for a hypothetical T2HKK experiment whose both detectors are at 1100 km. The left panels are for no decay and the right panels are for the case when α = 0. The red lines are for appearance channels, the blue lines are for disappearance channels and the green lines are for both appearance and disappearance channels together. In absence of decay disappearance channel does not have any octant sensitivity. However when appearance χ 2 is added to disappearance χ 2 , then the appearance probability being dependent on sin 2 θ 23 rises monotonically in the opposite octant and gives a large octant-sensitive contribution. In presence of decay, the disappearance channel itself has octant sensitivity, which enhances the overall octant sensitivity.

Combined Analysis
Here in Fig. 11, we have shown sensitivities of different experiments to constrain the decay parameter α(test) with (left) and without (right) decay in 'data' assuming the higher octant as the true octant. The plots in the left panel show how precisely these experiments can measure α and hence for a given true value how they can exclude the no decay scenarios. On the other hand, plots in the right panel show the constraints these experiments can put on invisible neutrino decay. The results are shown for all the three experiments ESSνSB, T2HK, T2HKK and the combination of ESSνSB with both T2HK and T2HKK. In the Fig. 11, the black dashed line shows the sensitivity of ESSνSB, the blue dashed-dotted line  is for T2HKK, while the red solid line represent the same for T2HK. The pink and cyan dashed line stands for the combination of ESSνSB with T2HK and T2HKK respectively.
It is observed from the left panel where we have assumed decay in the simulated data, that T2HKK has the highest precision out of the three and at 3σ it can precisely measure α in the range 1.278×10 −5 < α < 4.614×10 −5 eV 2 for the given true value of α = 3×10 −5 eV 2 . For the same true value of α T2HK gives better precision than ESSνSB and at 3σ the range is 3.408 × 10 −6 < α < 5.524 × 10 −5 eV 2 and 4.945 × 10 −6 < α < 7.303 × 10 −5 eV 2 respectively. Combining T2HKK with ESSνSB improves the bounds further and at 3σ the range is 1.597 × 10 −5 < α < 4.231 × 10 −5 eV 2 . Although the combination of T2HK and ESSνSB improves the bound, still it comes to be similar to that of T2HKK. In the right panel of Fig. 11, we observe that by combining two experiments it is possible to improve the constraint on α. Specially, the combination of T2HKK and ESSνSB gives the highest limit on α and at 3σ the constraint is α ≤ 1.193 × 10 −5 eV 2 . Similarly, T2HK+ ESSνSB also improves the constrain and the derived limit in this case is very much equal to that of T2HKK.
In Fig. 12 we present two complementary plots to Fig. 11. If neutrino decays in nature in addition to oscillation, then how it affects the measurements of θ 23 is shown in the left panel of Fig. 12. We have shown the results in the θ 23 (test)-α(test) parameter space at 95% C.L. We observe that the parameter space shrinks noticeably when we combine two experiments. Again in the right panel, we show the allowed region in θ 23 (test)-α(test) parameter space assuming no decay scenario in the simulated data. Here also, the combination of experiments gives better results that the individual experiments.

Summary & Conclusions
We have examined the physics potential of future long-baseline experiments T2HK/T2HKK and ESSνSB in the context of invisible neutrino decay. We performed our study for the case of normal hierarchy where ν 3 is the highest mass state and this is assumed to be unstable, decaying into lighter sterile states. We compared and contrasted the sensitivities to decay of these experiments with special emphasis on the location of the experiments at first and second oscillation maximum and the various factors that can affect the sensitivities.
The effect of decay appears as exp (−αL/E) where, α is the decay constant, L is the distance traveled and E is the energy. Hence, it is expected that second oscillation maxima experiments having higher baselines and/or lower energy than the first effect of decay depends on the L/E factor and second oscillation maximum occurs for lower energy it can be expected that effect of decay can be more pronounced for second oscillation maxima. However, the experimental sensitivities for the decay rate also depend on other vital factors like neutrino flux, resolution of the detector and the sensitivity of the experiment to resolve θ 23 . We expounded in detail which factors play important roles in determining the sensitivity to decay at the first and second oscillation maximum. Among the three experiments T2HK is designed to have its flux peak at the first oscillation maxima, T2HKK has one detector at the first oscillation maximum and one near the second oscillation maximum, while ESSνSB is an experiment at the second oscillation maximum. T2HK has a lower L/E as compared to the other experiments and has less sensitivity to decay. On the other hand being at the first oscillation maxima it has a better sensitivity to θ 23 and its octant. In comparison, ESSνSB has a higher L/E, but because of lower energy resolution as well as poor sensitivity to θ 23 , it has a reduced sensitivity to the decay constant α. The two detector set up of T2HKK helps since the detector at the first oscillation maximum helps in determining θ 23 while the second helps in better constraining α. We have studied the sensitivity of these experiments to the decay rate α for true θ 23 in both lower and higher octant.
We have also presented the octant sensitivity as a function of δ CP for the three experiments. We have obtained an interesting result that in presence of decay, the overall octant sensitivity is enhanced. This can be attributed to the octant sensitive contribution coming from the disappearance channel in presence of decay. Since T2HK is at the first oscillation maxima, the appearance channel already gives a good sensitivity and the effect of disappearance channel is not very significant. But since ESSνSB is at the second oscillation maximum the octant sensitivity coming from appearance channel is not very high and the disappearance channel plays a consequential role in enhancing the octant sensitivity in presence of decay. Similar feature is also observed for the detector at 1100 km baseline with flux peak near the second oscillation maxima.
The octant sensitivity coming from the disappearance channel increases with the decay constant α and also leads to an interesting interplay with the appearance channel in determining the overall minimum for experiments at second oscillation maxima. At the first oscillation maxima on the other hand, the appearance channel dominates. As a result of this different behaviour, for T2HKK, which has one detector at first and one at second oscillation maxima, a synergy is observed for true θ 23 in the higher octant leading to a χ 2 greater than the naive sum of the χ 2 s at two different baselines. Since at second oscillation maximum the appearance channel does not have very good octant sensitivity the χ 2 minimum for the second detector comes in the wrong octant driven by the disappearance channel. The χ 2 minima for the detector at the first oscillation maximum on the other hand stays in the correct octant (the wrong octant solution being disfavoured by the appearance channel). As an upshot, the global minima while combining the two baselines come in a position different from the individual minima resulting in the synergy. On the other hand, for ESSνSB, the first detector is not present and the leaning of the disappearance channel towards the wrong octant, weaken the sensitivity after a certain value of the decay constant α.
Given the correlation between α and θ 23 , the determination of θ 23 can also get affected, if there is decay in nature but it is ignored in the fit resulting in a wrong determination of θ 23 . We have examined this issue in detail and found that the shift of θ 23 from its true value, can be in either direction, depending on the octant of true θ 23 and whether the contribution is coming from bins centered at the maxima or minima of the survival probability, which in turn depends on the smearing and bin-width. This effect is more for the second oscillation maxima experiments since the probability band is narrower.
Apart from the sensitivity study for individual experiments we have also performed sensitivity studies for the two combinations of T2HK+ESSνSB and T2HKK+ESSνSB. We found that the sensitivity increases significantly if we combine the experiments. The sensitivity to the decay-rate attained by the experiments T2HK, T2HKK, ESSνSB and their combinations of T2HK+ESSνSB, T2HKK+ESSνSB have been summarized in the  Table 4. The sensitivity to the decay rate α and the neutrino decay lifetime at 3σ CL for the standalone experiments T2HK, T2HKK, ESSνSB and also for the combinations of T2HK+ESSνSB, T2HKK+ESSνSB.
We have also done a precision study assuming α = 3 × 10 −5 eV 2 in the data and found that the precision is maximum for the combination of T2HKK+ESSνSB while standalone T2HKK has the best precision. This also reiterates the fact that combination of first and second oscillation maximum can give the best sensitivity to decay. Although combination of the experiments improve the sensitivity to α, no significant improvement is observed for θ 23 .
In conclusion, we studied the sensitivity to decay at experiments at the first and second oscillation maxima and discussed the salient features. Overall the combination of first and second oscillation maxima is best for decay since at the first oscillation maximum the precision of θ 23 is more and at the second oscillation maximum the sensitivity to decay could be more because of enhanced baseline (for same energy) and/or lower energy (for same baseline). Our study underscores the importance of the disappearance channel in giving enhanced octant sensitivity for second oscillation maxima experiments in presence of decay and the importance of better energy resolution in obtaining a higher sensitivity to decay.
Note Added: While this work was being written a paper [70] came on arXiv, which also considered neutrino decay in the context of ESSνSB experiment. Our study contains a comparative analysis of T2HK, T2HKK and ESSνSB experiment.