Exploring invisible neutrino decay at ESSnuSB

We explore invisible neutrino decay in which a heavy active neutrino state decays into a light sterile neutrino state and present a comparative analysis of two baseline options, 540 km and 360 km, for the ESSnuSB experimental setup. Our analysis shows that ESSnuSB can put a bound on the decay parameter τ3/m3 = 2.64 (1.68) × 10−11 s/eV for the baseline option of 360 (540) km at 3σ. The expected bound obtained for 360 km is slightly better than the corresponding one of DUNE for a charged current (CC) analysis. Furthermore, we show that the capability of ESSnuSB to discover decay, and to measure the decay parameter precisely, is better for the baseline option of 540 km than that of 360 km. Regarding effects of decay in δCP measurements, we find that in general the CP violation discovery potential is better in the presence of decay. The change in CP precision is significant if one assumes decay in data but no decay in theory.


Introduction
The phenomenon of neutrino oscillations in the standard three-flavor scenario can be expressed by three mixing angles θ 12 , θ 23 , and θ 13 , two mass-squared differences ∆m 2 21 = m 2 2 − m 2 1 and ∆m 2 31 = m 2 3 − m 2 1 , and one Dirac-type CP-violating phase δ CP . During the past decades, data from solar, atmospheric, accelerator, and reactor neutrino experiments have successfully been able to determine the values of the parameters θ 12 , θ 13 , ∆m 2 21 , and |∆m 2 31 | to an excellent precision. The parameters, which are unknown at this moment, are: (i) the mass ordering of neutrinos, i.e., ∆m 2 31 > 0 known as normal ordering (NO) or ∆m 2 31 < 0 known as inverted ordering (IO), (ii) the octant of the mixing angle θ 23 , i.e., θ 23 > 45 • known as the higher octant (HO) or θ 23 < 45 • known as the lower octant (LO), and (iii) the true value of δ CP and its precision. At this moment, data from ongoing experiments provide a hint towards the true ordering as NO, the true octant as HO, and δ CP = −90 • [1]. There are many dedicated experiments to establish these hints concretely.
Apart from standard neutrino oscillation physics, there are many new physics scenarios which can be probed at neutrino oscillation experiments. Invisible neutrino decay is an example of one such scenario [2]. In invisible neutrino decay, a heavy neutrino state decays into a light neutrino state, which is sterile and therefore invisible. 1 Theoretically, for Dirac neutrinos, this scenario can arise if there exists a coupling between the neutrinos and a light scalar boson [3]. In this case, decay can be defined as ν j → ν iR + χ, where ν iR is a right-handed singlet and χ is an iso-singlet scalar. For Majorana neutrinos having a pseudo-scalar coupling with a Majoron J and a sterile neutrino ν s , it is possible that ν j → ν s +J [4]. Irrespective of the model, in the presence of neutrino decay, the Hamiltonian

JHEP05(2021)133
for neutrino propagation is modified, and therefore, one can measure the decay parameter as well as the effect of decay on the measurement of the standard oscillation parameters in a neutrino oscillation experiment. In principle, all three neutrino states ν 1 , ν 2 , and ν 3 can decay invisibly. The decay due to ν 2 is severely constrained from solar neutrino data [5] and the constraints from supernova SN1987A [6] is applicable to decays of ν 2 and ν 1 [7]. In ref. [8], decays of supernova neutrinos have been investigated in great theoretical detail. Decay due to ν 3 can be measured in present and future accelerator, atmospheric, and reactor neutrino experiments. Recently, a lot of work has been performed in this direction. Studies of invisible decay in the accelerator neutrino experiments T2K [9], NOνA [10], MINOS [11], DUNE [12], and MOMENT [13] can be found in refs. [14][15][16][17][18]. For the study of invisible neutrino decay in the ongoing atmospheric neutrino experiment Super-Kamiokande (SK) [19] and the future atmospheric neutrino experiment INO [20], see refs. [21][22][23], for the study with atmospheric neutrino data of the future ultra-high energy neutrino experiment KM3NeT-ORCA [24], see ref. [25], and for the medium baseline reactor neutrino experiment JUNO [26], see ref. [27].
In this paper, we study the scenario of invisible decay of the neutrino state ν 3 in the future long-baseline experiment ESSnuSB [28,29]. The primary aim of the ESSnuSB experiment is to measure the Dirac CP-violating phase δ CP with high precision at the second oscillation maximum [30]. Currently, there are two possible baseline options for ESSnuSB under consideration, which are 540 km and 360 km. In the present work, we consider both baseline options to estimate the sensitivity to invisible neutrino decay at ESSnuSB. The topics, which we address in this work, are the following: (i) the capability of ESSnuSB to put bound on the decay parameter assuming there is no decay in Nature, (ii) the capability of ESSnuSB to discover neutrino decay assuming that neutrinos decay in Nature, (iii) how precisely ESSnuSB can measure the decay parameter if there exists neutrino decay in Nature, and (iv) the effect of neutrino decay on the measurement of δ CP . Since the sensitivity of ESSnuSB to measure the neutrino mass ordering, the octant of θ 23 , and the precision of θ 23 is weak [31][32][33][34], we will not address such potential measurements in this work. This paper is organized as follows. In section 2, we will discuss how the phenomenon of neutrino oscillations in the standard three-flavor neutrino oscillation scenario is altered in the presence of invisible neutrino decay. In section 3, we will present the experimental setup of ESSnuSB along with the simulation details. In section 4, we will present our results, and finally in section 5, we will summarize and conclude.

Effects of invisible neutrino decay in neutrino oscillations
In the standard three-flavor neutrino oscillation framework, the evolution equation can be written as

JHEP05(2021)133
where U is the leptonic mixing matrix, E is the neutrino energy, and V = √ 2G F n e is the effective matter potential with G F being the Fermi coupling constant and n e the electron density of matter (along the neutrino trajectory). The sign of V is positive for neutrinos and negative for antineutrinos. Note that the Hamiltonian in eq. (2.1) is Hermitian and can be diagonalized by a unitary transformation. Assuming that the neutrino state ν 3 decays into a sterile state, which effectively means that ∆m 2 31 → ∆m 2 31 − iα 3 [2], the evolution equation is modified to where α 3 ≡ m 3 /τ 3 with τ 3 being the rest-frame lifetime of the neutrino state ν 3 having mass m 3 . We assume that the invisible neutrino state is lighter than the lightest active neutrino state, and therefore, decay is possible for both NO and IO. It is interesting to note that the Hamiltonian with decay is no longer Hermitian and therefore cannot be diagonalized by a unitary transformation. In this case, we follow the prescription given in ref. [35] to numerically diagonalize this Hamiltonian to calculate the neutrino oscillation probabilities. For ESSnuSB, the impact of decay comes in both the appearance channel (ν µ → ν e ) and the disappearance channel (ν µ → ν µ ). In vacuum, using the approximations where s ij ≡ sin θ ij , c ij ≡ cos θ ij , ∆ atm ≡ ∆m 2 31 L/(2E), and Γ 3 ≡ α 3 L/E = m 3 L/(τ 3 E). Note that as eq. (2.3) is derived using the approximation ∆m 2 21 ∼ 0, there is no δ CP term in this equation. This is because the parameter δ CP appears in the appearance channel probability from the interference term between ∆m 2 21 and ∆m 2 31 and we set ∆m 2 21 = 0 in the above for simplicity. From eqs. (2.3) and (2.4), we understand that invisible neutrino decay leads to a depletion in the number of events in both appearance and disappearance channels. We also note that the decay parameter appears together with the quantity L/E. Since the neutrino energy E is same for both baseline options for ESSnuSB, which are L = 540 km and L = 360 km, we expect to observe different effects of decay for different baseline lengths.
We must point out that while the approximate neutrino oscillation probabilities given in eqs. (  full three-flavor effects. We have explicitly checked that the approximate expression of the probability for the disappearance channel given in eq. (2.4) matches with the exact numerical results to a very good accuracy. On the other hand, for the appearance channel, the approximate expression given in eq. (2.3) has a mismatch with the exact probability, since as was pointed out before, we have neglected ∆m 2 21 in this approximation. We will observe in section 4 that the sensitivity of ESSnuSB to invisible neutrino decay comes essentially from the disappearance channel, while the significance of the appearance channel is marginal. Therefore, in the rest of this section, we will only look into the disappearance channel, using the approximate probability given in eq. (2.4) in order to understand the impact of decay.
Since the main aim of the current work is to determine how well ESSnuSB can distinguish the standard case from the decay case, we present in figure 1 the probability difference where P std µµ is the neutrino oscillation probability in the standard case without decay using the mixing angle θ 23 = 49.2 • and P decay µµ is the probability as a function of the effective mixing angle θ 23 with invisible decay using the decay parameter τ 3 /m 3 = 1.0 × 10 −11 s/eV. We divide the probability difference ∆P µµ in eq. (2.5) as follows (2.9) In figure 1, each of these terms are shown for three different values of θ 23 , i.e., 41 • (dotted curves), 45 • (solid curves), and 49 • (dashed curves). The purpose of this figure is mainly to show the effect of θ 23 on the difference between the oscillation probability for no decay and decay. Note that ∆P 2 and ∆P 3 are zero in the standard case with no decay, whereas ∆P 1 does not depend on the decay parameter. Therefore, since ∆P 2 and ∆P 3 are equal to zero if Γ 3 = 0 (i.e., there is no decay), the expressions for ∆P 2 and ∆P 3 are independent of θ 23 , but not of θ 23 . We will use this figure in the later sections to explain our main results on the sensitivity of ESSnuSB to the decay parameter. We can already note a few points from the figure. The term ∆P 1 depends on sin 2 (∆ atm /2) and is hence observed to be oscillating with a small amplitude as a function of L/E, since it is proportional to the difference 4c 2 13 s 2 , where s 23 and s 23 are taken to be very close to each other. The reason is that it does not depend on the decay parameter and the small non-zero oscillatory value it obtains (roughly 3 %) comes from the fact that θ 23 = 49.2 • , while θ 23 is assumed to be 41 • , 45 • , and 49 • , for the three cases, respectively. On the other hand, we see that the terms ∆P 2 and ∆P 3 have rather non-trivial behaviors and values. The term ∆P 2 is independent of ∆ atm , and hence is non-oscillatory. It increases almost linearly with L/E as ∆P 2 c 4 13 s 4 23 α 3 L/E for Γ 3 small. The spread in ∆P 2 with s 23 can be observed in the figure. In contrast, the term ∆P 3 is observed to be oscillating [with a relatively large amplitude proportional to 2c 2 13 s 2 23 (1 − c 2 13 s 2 23 )] as a function of L/E, since it depends on cos ∆ atm . We show by the two grey vertical lines the L/E corresponding to the 360 km and 540 km baseline options, respectively, where we take E = 0.36 GeV, which is the energy where the ESSnuSB flux peaks. Note from the figure that for the 360 km baseline option, ∆P 3 has a peak, whereas for the 540 km baseline option, it has a trough. Contrast this with the behavior of ∆P 2 , which is increasing almost linearly and is positive for both the 360 km and 540 km baseline options. This means that while ∆P 2 and ∆P 3 add up for the 360 km baseline option, they cancel each other for the 540 km one. We will see that this has far reaching consequences when it comes to marginalization over θ 23 .

Experimental setup and simulation details of ESSnuSB
We use the software GLoBES [38,39] to simulate the sensitivity of the ESSnuSB experiment. We consider a water-Cherenkov far detector [40] of fiducial volume 507 kt located at distance of either 540 km or 360 km from the neutrino source. We also consider an identical near detector located at a distance of 0.5 km having volume of 0.1 kt. For the neutrino source, we consider protons of 2.5 GeV originating from a beam of 5 MW capable of delivering 2.7 × 10 23 protons on target per year. We assume a total run-time of 10 Table 2. Best-fit values of the neutrino oscillation parameters and the corresponding 3σ allowed values for the standard three-flavor scenario. The values have been adopted from ref. [1].
divided into 5 years in neutrino mode and 5 years in antineutrino mode. We consider correlated systematics between far and near detectors with the errors as given in ref. [41] and we list them in table 1 for convenience.
We estimate the statistical χ 2 function using where n is the number of energy bins, N true is the number of true events, and N test is the number of test events and incorporate the systematics by the method of pulls. Unless otherwise mentioned, in our simulation, we generate the data with the best-fit values from the global analysis of the world neutrino data as obtained by NuFIT v5.0 [1] and we present them in table 2. In the fit, we minimize over the parameters θ 13 , θ 23 , and δ CP in their current 3σ ranges, as given in table 2. We keep the parameters θ 12 , ∆m 2 21 , and ∆m 2 31 fixed to their best-fit values in the fit. Throughout our analysis, we assume that the mass ordering of the neutrinos is known and is NO.

Simulation results
We discuss the sensitivity of ESSnuSB in presence of invisible neutrino decay. Our strategy is as follows. First, we show how the neutrino oscillation probabilities in the ν µ → ν e and ν µ → ν µ channels are modified due to the presence of decay for both baseline options of ESSnuSB. We also show the event rates without decay, to demonstrate the relevant energy values from where the sensitivity stems. Then, we study the capability of ESSnuSB to put bound on the decay parameter and compare our results with other experimental setups. We also present the discovery potential of ESSnuSB to observe neutrino decay. Next, we study how precisely ESSnuSB can measure the decay parameter if decay exists in Nature. Finally, we discuss the effects of decay in CP violation discovery and CP precision measurements at ESSnuSB.

Discussion at the probability and event level
In figure 2, we present the neutrino oscillation probabilities and neutrino event rates of ESSnuSB as functions of neutrino energy E. The left panel is for the appearance channel, whereas the right panel is for the disappearance channel. In each panel, the blue curve (histogram) corresponds to the probability (event rate) for L = 540 km and the red curve (histogram) corresponds to the probability (event rate) for L = 360 km. The solid curves correspond to the case when there is no decay. The histograms are plotted for the no decay scenario. To generate these curves, we use the best-fit values of the standard neutrino oscillation parameters given in table 2. To understand the effect of decay in the neutrino oscillation probabilities, we consider two values of the decay parameter, i.e., τ 3 /m 3 = 10 −11 s/eV and τ 3 /m 3 = 10 −12 s/eV, and they are shown by dashed and dotted curves, respectively. In the appearance channel (the left panel of figure 2), we realize that for L = 540 km, most of the sensitivity is due to the second oscillation maximum. In fact, for L = 360 km, the sensitivity comes from both the first and second oscillation maxima.
In the disappearance channel (the right panel of figure 2), we note that for L = 540 km, the JHEP05(2021)133 10 -11 10 -10 10 -9  sensitivity actually comes from the second oscillation minimum, whereas for L = 360 km, it stems from the second oscillation maximum. From the figure, we also observe that the effect of decay is non-negligible in both the appearance and disappearance channels. Between τ 3 /m 3 = 10 −11 s/eV and τ 3 /m 3 = 10 −12 s/eV, the separation of the no decay curve and the decay curve is larger for τ 3 /m 3 = 10 −12 s/eV, reflecting the fact that the change in the probability is much larger for a shorter lifetime of the neutrino state. Note that the number of events are higher for L = 360 km as compared to L = 540 km, since the number events are proportional to 1/L 2 .

Sensitivity and discovery potential in presence of invisible neutrino decay
In figure 3, we plot the capability of ESSnuSB to put bound on the decay parameter (or the sensitivity χ 2 ) in the left panel and its potential to discover decay (or the discovery χ 2 ) in the right panel.
In the left panel, we have not considered decay in data, whereas in the right panel, we have not considered decay in theory. The black dotted horizontal lines represent the value of the χ 2 corresponding to 90 % C.L and 3σ. In each panel, the blue curve represents the sensitivity for L = 540 km and the red curve represents the sensitivity for L = 360 km. From the left panel, we observe that the bound on the decay parameter is 1.68×10 −11 s/eV for L = 540 km and 2.64 × 10 −11 s/eV for L = 360 km at 3σ. On the other hand, ESSnuSB can discover decay at 3σ if the true value of the decay parameter is 2.56 × 10 −11 s/eV for L = 360 km and 3.22 × 10 −11 s/eV for L = 540 km. We note that the sensitivity χ 2 is better for the 360 km baseline option of ESSnuSB, whereas the discovery χ 2 is better for the 540 km baseline option of ESSnuSB. To understand this result, we calculate the contribution to the χ 2 function from the individual appearance and disappearance channels for both 540 km and 360 km and for both sensitivity and discovery. We find that the most of the contribution comes from the disappearance channel, which is plotted in figure 4.
In figure 4, as was the case in figure 3, the left panel is for the sensitivity χ 2 and the right panel is for the discovery χ 2 . The blue curves correspond to L = 540 km and red curves correspond to L = 360 km. In each panel, the solid curves correspond to the case when the χ 2 is minimized over the parameter θ 23 in theory and the dashed curves correspond to the case when the parameter θ 23 is kept fixed to its best-fit value in theory. From the panels, we note that when θ 23 is fixed, the 540 km baseline option of ESSnuSB is better than the 360 km baseline option, for both sensitivity and discovery. However, when the χ 2 is minimized over the parameter θ 23 , we find that the L = 360 km option is better for sensitivity, whereas L = 540 km is better for discovery. The reason is that the effect of minimizing over θ 23 is significant for the sensitivity χ 2 for the L = 540 km option. In order to understand this, let us first see why the L = 540 km option is better for both sensitivity and discovery when θ 23 is kept fixed.
In figure 5, we display as a function of E the difference between the number of events without decay and the number of events with decay divided by the number of events with decay, i.e., ∆N/N D (left panel), and the difference between the number of events without decay and the number of events with decay divided by the number of events without decay, i.e., ∆N/N WD (right panel) for the disappearance channel. We generate these panels for the best-fit values of the neutrino oscillation parameters given in table 2    presented curves are for different values of θ 23 , corresponding to either the true value of θ 23 or the value of θ 23 where the χ 2 minimum is occurring, considering both decay and no decay. The value of the decay parameter is τ 3 /m 3 = 10 −11 s/eV. From figure 2, we understand that for the disappearance channel the sensitivity comes from E ∼ 0.4 GeV. Therefore, in these panels, we focus on the values of the probability around E ∼ 0.4 GeV. For L = 540 km, we note that the separation between the blue solid curve (which is the true point for the sensitivity χ 2 ) and the red solid curve (which is the test point for the sensitivity χ 2 for fixed θ 23 ) is larger than the separation between the blue solid curve and the blue dashed curve (which is the test point for the sensitivity χ 2 when θ 23 is minimized). On the other hand, for L = 360 km, the separation between the blue solid curve and the red solid curve is similar to the separation between the blue solid curve and the blue dashed

JHEP05(2021)133
curve. For this reason, the effect of θ 23 marginalization is larger for the sensitivity χ 2 for L = 540 km. Similarly, for the discovery χ 2 , one can see that the separation between the red solid curve (which is the true point for the discovery χ 2 ) and the blue solid curve (which is the test point for the discovery χ 2 for fixed θ 23 ) and the separation between the red solid curve and the red dashed curve (which is the test point for the discovery χ 2 when θ 23 is minimized) is very similar for both L = 540 km and L = 360 km. This is why the discovery χ 2 is not affected by the θ 23 minimization. To see this even more clearly, let us go back to our approximate expression for the survival probability given in eq. (2.4) and the related figure, i.e., figure 1. For the sensitivity study, the data are generated for no decay and we assume θ 23 = 49.2 • . This is then fitted with a theory with decay, while allowing θ 23 to vary within its current 3σ bound. The fit, of course, tries to minimize the difference between the data and the theory -this essentially means that ∆P µµ is reduced to the best possible case in the fit. If we now look at the three terms plotted in figure 1 for the three test values of θ 23 , we note the following. The θ 23 dependence is mostly coming from ∆P 2 and for θ 23 = 49 • , the theory is worst as compared to the data. The difference between data and theory (i.e., ∆P µµ ) starts to reduce as we reduce θ 23 in the fit. Note that the spread in ∆P µµ due to ∆P 2 coming from θ 23 is larger for the 540 km baseline option as compared to the 360 km one. Therefore, it is expected that the effect of θ 23 marginalization will be larger for the 540 km baseline option. Indeed, this is what we have observed in figure 1. In addition, there is yet another point to note that brings a difference between the 360 km and 540 km baseline options. If we look at ∆P 3 , we see that while this term itself is largely independent of θ 23 , it has a very interesting correlation with respect to ∆P 2 , and that makes it relevant in θ 23 marginalization. The term ∆P 3 is oscillatory and we can see from figure 1 that while it has a peak for the 360 km baseline option, it has a trough for the 540 km baseline option. This results in partial cancellation between ∆P 2 and ∆P 3 for the 540 km option. Since ∆P 2 depends on θ 23 , the extent of this cancellation depends on the value of θ 23 and as a result the χ 2 for the 540 km option becomes lower as a result of marginalization over θ 23 .
A similar argument can be used to understand how θ 23 marginalization affects the χ 2 for discovery and why for that case too the effect is larger for the 540 km baseline option. However, there is a major difference. The data are now generated for decay, while the fit is performed for standard oscillations. In figure 1, this would mean that ∆P µµ would flip sign for all three terms, but this is not the main issue. The more important difference is that now we have standard oscillations in the fit. This means that ∆P 2 and ∆P 3 are absent in the fit and only ∆P 1 can be manipulated to minimize the χ 2 . Since ∆P 1 has only a small dependence on θ 23 , the effect of marginalization on the discovery χ 2 is smaller. The difference in θ 23 marginalization between the 360 km and 540 km baseline options can also be seen from the same figure. We observe that for the 360 km baseline option, ∆P 1 is almost zero, and hence independent of θ 23 . Therefore, there is a small effect of marginalization on the discovery χ 2 for this case. For the 540 km baseline option, the effect is larger, since for this case, there is a peak for ∆P 1 . However, even for this baseline option, the θ 23 marginalization effect is significantly smaller as compared to the sensitivity case for the reason mentioned above.  In figure 7, we study the impact of θ 23 marginalization further and show the individual χ 2 for the appearance (red curves) and disappearance (blue curves) channels as functions of θ 23 assumed in the test. The dashed curves are for the case of discovery, whereas the solid curves are relevant for sensitivity. The left panel is displayed for the 540 km baseline option, while the right panel for the 360 km one. The data for all cases are generated at 49.2 • . The effect of the appearance channel is negligible, as was pointed out earlier, so we concentrate only on the disappearance channel. We note that for the 540 km baseline option, the best fit for the sensitivity case comes in the lower octant, while the value of θ 23 = 49.2 • is strongly disfavored. For the discovery case, the effect is less dependent on θ 23 , but the effect of marginalization exists. This is consistent with our discussion in the previous paragraphs in connection with ∆P 2 in figure 1.
In figure 8, we show the combined χ 2 coming from both disappearance and appearance channels as a function of θ 23 in the test. The left panel is displayed for the 540 km baseline option, while the right panel for the 360 km one. The data are generated for θ 23 = 49.2 • . We show this for three different cases -the red solid curves show the χ 2 for the case when we have standard oscillations in both data and theory. Therefore, these curves show the octant sensitivity for the 360 km and 540 km baseline options. We observe that the octant sensitivity of the 360 km baseline option is better than that of the 540 km one. Indeed, the octant sensitivity of the 540 km baseline option is rather poor. The blue dot-dashed curves show the χ 2 for the case when we have no decay in data and decay in theory. This is the case for the sensitivity study. The low θ 23 octant sensitivity of the 540 km baseline option results in the true octant getting disfavored at a very high significance level in the sensitivity case. With decay switched on in the theory, the fit prefers to choose the wrong octant to lower the χ 2 coming from the inclusion of the decay parameter in the fit. in the previous paragraphs. For the 360 km baseline option, the fit also prefers the wrong θ 23 octant for the blue dot-dashed curve, but the χ 2 difference between the two octants is lower than for the 540 km case.
Let us now briefly compare the bound on the decay parameter expected for the ESS-nuSB experiment, with the corresponding values for other accelerator, atmospheric, and reactor neutrino experiments. In table 3, we list the (expected) bounds on τ 3 /m 3 from different experiments along with what has been obtained in this work.
In this table, the bounds for T2K + NOνA, T2K + MINOS, and SK+MINOS are obtained using real data, whereas for the other experiments, the bounds correspond to numerical simulations obtained by different groups. To understand under which assumptions the bounds are obtained, we list the input parameters for different experiments in table 4. For the bounds, which were obtained from real data, we provide the references of the data for the relevant experiments as well as the parameters that are minimized during the fits. For the bounds, which were estimated using simulated data, we provide the details of the fiducial volume, the run-time, and the true values of the oscillation parameters as well as the parameters that are minimized in the fits for the relevant experiments. As the bounds for T2K + NOνA, T2K + MINOS, and SK+MINOS are obtained using real data and all relevant parameters are minimized in the fits, these bounds are robust. For MO-MENT, ESSnuSB, DUNE, JUNO, INO, and KM3NeT-ORCA, the bounds mainly depend on the true values of the oscillation parameters, the fiducial volume, and the run-time.
For different values of these parameters, the bounds can be different. In particular, we note that except for θ 23 and δ CP , the true values of the oscillation parameters are very similar for all experiments. Since the sensitivity mainly comes from the disappearance channel, the bounds are not expected to change much with respect to different values of δ CP . However, different true values of θ 23    From table 3, we note that the bounds obtained from the ongoing atmospheric neutrino experiment SK along with MINOS, future atmospheric experiment INO, and the atmospheric data of the future ultra-high energy neutrino experiment KM3NeT-ORCA are one order of magnitude stronger than the bounds obtained from future accelerator and reactor experiments. On the other hand, the bound obtained from the currently running accelerator experiments T2K and NOνA, along with MINOS, are one order of magnitude lower than the expected bounds from future accelerator and reactor experiments. Among ESSnuSB, DUNE, MOMENT, and JUNO, the expected bound for the JUNO experiment is the strongest and the one for MOMENT is weakest. The sensitivity of ESSnuSB is slightly better than that of DUNE obtained from a charged current analysis, for the baseline option of 360 km and worse for the baseline option of 540 km, but better than the one of MOMENT. Note that a stronger bound for DUNE could possibly be obtained in a combined analysis of both charged current and neutral current events.  Table 4. Details of the input parameters, which are used to estimate the bounds on the decay parameter τ 3 /m 3 .

Precision measurement of the decay parameter
In this section, we discuss the capability of ESSnuSB to measure the value of the decay parameter τ 3 /m 3 , if there exists invisible decay in Nature. In figure 9, we plot the precision

CP sensitivity in presence invisible neutrino decay
In this section, we discuss the CP violation discovery and the CP precision sensitivity of ESSnuSB in presence of invisible neutrino decay. The CP violation discovery potential of an experiment is defined as its capability to distinguish a true value of δ CP from 0 • and 180 • , whereas the CP precision capability of an experiment is defined as its capability to exclude all values of δ CP other than the true value. In the upper panels of figure 10, we present the CP violation discovery ∆χ 2 as a function of the true δ CP , and in the lower panels, we present the 1σ precision as a function of the true δ CP . In each row, the left panel is for L = 540 km and the right panel is for L = 360 km. In all panels, the red curve corresponds to the standard neutrino oscillation scenario, where there is no decay in both data and theory. To study the effect of decay in the measurement of δ CP , we consider two scenarios: (i) the presence of decay in both data and theory (blue curve) and (ii) the presence of decay in data but not in theory (green curve). Scenario (i) addresses the question how the CP sensitivity of ESSnuSB will be altered if there exists decay in Nature and scenario (ii) addresses the question what happens if we try to fit a theory without decay to a data set that includes decay. Scenario (ii) is usually the case because when a data set from a neutrino oscillation experiment is available, it will first be fitted using the standard neutrino oscillation framework. The value of the decay parameter is τ 3 /m 3 = 1.0 × 10 −11 s/eV and to generate the blue curves, we minimize τ 3 /m 3 in the 3σ range as obtained from figure 9. When decay is present in both data and theory, the CP violation sensitivity is better than for the standard scenario around δ CP = ±90 • . The effect is larger for L = 540 km compared to L = 360 km. For CP precision, the sensitivity  is weaker for the true δ CP values of ±90 • and stronger for the true δ CP values of 0 • and 180 • as compared to the CP precision in the standard three-flavor neutrino oscillation scenario. This deviation is larger for L = 360 km than for L = 540 km. Further, the CP violation sensitivity due to standard neutrino oscillations and the sensitivity due to the presence of decay in data but not in theory are almost identical for δ CP = 90 • . However, for δ CP = −90 • , the CP violation sensitivity due to the presence of decay in data is better than for the standard case. These conclusions are true for both baseline length options of ESSnuSB. In this case, the value of χ 2 min is 45.4 for L = 360 km and 82.2 for L = 540 km. Nevertheless, when decay is present in theory but not in data, the sensitivity to precision is significantly different compared to the standard scenario. This can be qualitatively understood in the following way. When there is decay or no decay in both data and theory, the χ 2 minimum for CP precision is always zero and it appears with the true value of δ CP . However, if there is decay in data and not in theory, the χ 2 minimum is non-zero and it can appear with a different value of δ CP other than the true value of δ CP and this can change the sensitivity significantly. This can be viewed as a signature of the presence of decay that, if there exists decay in data and one tries to fit a theory without decay, the CP precision will be significantly different.

Summary and conclusions
Invisible decay is defined as the decay of a heavy neutrino state into a lighter sterile neutrino state. Due to the presence of invisible neutrino decay, the evolution equation of neutrino propagation is altered, and therefore, it is possible to probe invisible decay in neutrino oscillation experiments. The bound on the decay parameters due to decay of ν 1 and ν 2 comes from the solar and supernova neutrino experiments and the decay due to ν 3 can be probed in current and future accelerator, atmospheric, ultra-high energy, and reactor neutrino experiments. In this work, we have studied the physics sensitivity of the ESS-nuSB experimental setup in presence of invisible decay. The primary aim of the ESSnuSB experiment is to measure the leptonic CP-violating phase δ CP at the second oscillation maximum. In the current work, we have considered two baseline options of ESSnuSB, which are 540 km and 360 km, respectively, and studied (i) the capability to put bounds on the decay parameter, (ii) the capability to discover invisible decay, (iii) the bounds on the decay parameter obtained with other experiments, (iv) the potential to measure a particular value of decay parameter, and (v) the effect of decay on the measurement of δ CP . In our work, we have shown that the capability of ESSnuSB to put bounds on the decay parameter or the sensitivity χ 2 is better for the baseline option of 360 km, while its capability to discover decay with a particular value of decay parameter or the discovery χ 2 is better for the baseline option of 540 km. The sensitivity χ 2 is worse for the baseline option of 540 km due to the effect of θ 23 on decay. Our results show that the bound obtained for ESSnuSB with the baseline option of 360 km is better than the one of DUNE (CC) and the bound obtained with the baseline option of 540 km is worse than the one of DUNE (CC), but better than the one of MOMENT. Note that the robustness of the bounds obtained from different experiments depends on experimental specifications, true value of the oscillation parameters, etc. Therefore, the conclusion that has been derived in this work can change based on different assumptions upon which the bounds have been estimated. Assuming true values of the decay parameter outside the bounds obtained for ESSnuSB, we have shown that the capability of ESSnuSB to precisely measure a value of the decay parameter is better for L = 540 km than for L = 340 km and as the lifetime of decay increases, precision becomes worse. Regarding the measurement of δ CP , we have found that in presence of decay in both data and theory, the CP violation sensitivity is better for δ CP = ±90 • , whereas the CP precision capability is stronger for δ CP = 0 • and 180 • and weaker for δ CP = ±90 • compared to the standard three-flavor neutrino oscillation scenario. For CP violation sensitivity, the deviation is larger for L = 540 km, whereas for CP precision capability, the deviation is larger for L = 360 km. Further, if we try to fit data, which include invisible neutrino decay, with a theory that does not incorporate neutrino decay, then there exist a visible difference in the CP violation sensitivity due to effect of decay around δ CP = −90 • , whereas the CP precision capability is largely affected as compared to the standard case. The significant change in the CP precision capability can be viewed as a signature of the presence of decay when one tries to fit data that contain decay with a theory that does not include decay. In summary, our results have shown that ESSnuSB provides a good opportunity to study invisible neutrino decay and the primary goal of the ESSnuSB to measure δ CP can be affected due to the presence of such decay.