Constraints on Neutrino Lifetime from the Sudbury Neutrino Observatory

TL;DR

This paper uses data from the Sudbury Neutrino Observatory to set new lower bounds on the lifetime of neutrino mass state bc2, constraining possible neutrino decay models through energy-dependent survival probability analysis.

Contribution

It provides the first comprehensive fit of a neutrino decay model to all three phases of SNO bc8B solar neutrino data, improving constraints on neutrino lifetime.

Findings

Neutrino lifetime bc2 > 8.08 ^{-5} s/eV at 90% confidence

Combined analysis yields bc2 > 1.04 ^{-3} s/eV at 99% confidence

Results place stringent limits on neutrino decay scenarios.

Abstract

The long baseline between the Earth and the Sun makes solar neutrinos an excellent test beam for exploring possible neutrino decay. The signature of such decay would be an energy-dependent distortion of the traditional survival probability which can be fit for using well-developed and high precision analysis methods. Here a model including neutrino decay is fit to all three phases of $^8$B solar neutrino data taken by the Sudbury Neutrino Observatory. This fit constrains the lifetime of neutrino mass state $\nu_2$ to be ${>8.08\times10^{-5}}$ s/eV at $90\%$ confidence. An analysis combining this SNO result with those from other solar neutrino experiments results in a combined limit for the lifetime of mass state $\nu_2$ of ${>1.04\times10^{-3}}$ s/eV at $99\%$ confidence.

Tables

Table 1. Table 1: Reproduced here are the mixing parameters used in this analysis taken from Daya Bay  An  et al. ( 2017 ) , Super-Kamiokande  Wendell  et al. ( 2010 ) , and KamLAND  Abe  et al. ( 2008 ) results.

Parameter	Value	Ref
$\Delta m_{21}^2$	${7.58}_{-0.13}^{+0.14}{\left(stat\right)}_{-0.15}^{+0.15}\left(syst\right)\times {10}^{-5}$ eV2	Abe et al. (2008)
${\tan }^2{\theta }_{12}$	${0.56}_{-0.07}^{+0.10}{\left(stat\right)}_{-0.06}^{+0.10}\left(syst\right)$	Abe et al. (2008)
$\left|\Delta m_{32}^2\right|$	$2.45\pm 0.06\left(stat\right)\pm 0.06\left(syst\right)\times {10}^{-3}$ eV2	An et al. (2017)
${\sin }^{2}2{\theta }_{13}$	$0.0841\pm 0.0072\left(stat\right)\pm 0.0019\left(syst\right)$	An et al. (2017)
${\sin }^2{\theta }_{23}$	${0.5}_{-0.062}^{+0.058}$	Wendell et al. (2010)

Equations

$$|{\nu }_{\alpha }⟩=\sum _i{U}_{\alpha i}^{\ast }|{\nu }_i⟩.$$

$$⟨{\nu }_i|H_0|{\nu }_j⟩=\frac{1}{2E}\left[\begin{array}{ccc}{m}_1^2& 0& 0\\ 0& m_2^2& 0\\ 0& 0& m_3^2\end{array}\right].$$

$$P_{\alpha \beta }={\left|\text{innerproduct}{\nu }_{\beta }(t){\nu }_{\alpha }\right|}^2={\left|\sum _i{U}_{\alpha i}^{\ast }{U}_{\beta i}{e}^{-i{m}_i^{2}L/(2E)}\right|}^2.$$

$$V_e=\sqrt{2}{G}_F{n}_e.$$

$$H_{MSW}=H_0+U\left[\begin{array}{ccc}{V}_e& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]U^{†}.$$

$${\left|\text{innerproduct}{\nu }_{mi}(V_e(0))\nu (0)\right|}^2={\left|\text{innerproduct}{\nu }_{mi}(V_e(t))\nu (t)\right|}^2$$

$${\varphi }_i={\left|\text{innerproduct}{\nu }_{mi}(V_e){\nu }_e\right|}^2.$$

$$k_i=\frac{{\tau }_i}{m_i}.$$

$${\psi }_i\approx e^{-L/(E{k}_i)}{\varphi }_i=e^{-L/(E{k}_i)}{\left|\text{innerproduct}{\nu }_{mi}(V_e){\nu }_e\right|}^2$$

$$\begin{array}[]{rcl}P_{ee}&=&\sum_{i}\psi_{i}|U_{ie}|^{2}\\ P_{ea}&=&\sum_{i}\psi_{i}|U_{i\mu}|^{2}+\psi_{i}|U_{i\tau}|^{2}.\end{array}$$

$${\Phi }_T={\left(\overline{P_{ee}}+\overline{P_{ea}}\frac{{\sigma }_a}{{\sigma }_e}\right)}^{-1}{\Phi }_e.$$

$$-\ln (\mathcal{L})=\frac{({\Phi }_T-{\Phi }_{SSM}{)}^2}{2({\sigma }_T^2+{\sigma }_{SSM}^2)}$$

Full text

I Introduction

Nuclear reactions in the core of the Sun produce electron flavor neutrinos at rates which can be predicted by solar models. Neutrinos produced in the solar 8B reaction propagate to the Earth and are detected as electron flavor neutrinos with a probability, $P_{ee}$, of roughly $1/3$, with the remainder converted to $\nu_{\mu}$ or $\nu_{\tau}$. Analysis of data from the Sudbury Neutrino Observatory (SNO) Aharmim et al. (2013a) and Super-Kamiokande (SK) Abe et al. (2011a) has shown the origin of this $1/3$ survival probability to be due to mixing of neutrino states with finite mass that are distinct from the flavor states in which neutrinos are produced and interact. Many other experiments Cleveland et al. (1998); Abdurashitov et al. (1999); Anselmann et al. (1992); Abe et al. (2008); Bellini et al. (2014) have made precision measurements of solar neutrino fluxes probing the rich physics of neutrino mixing, and are consistent with this conclusion. With the discovery of finite neutrino mass comes the possibility that neutrinos may be unstable and could decay to some lighter particle.

Neutrino decay was first explored as a possible explanation for the less-than-unity survival probability of electron flavor neutrinos Bahcall et al. (1972) (the so-called “solar neutrino problem”). The flavor-tagging Ahmad et al. (2001) and flavor-neutral Ahmed et al. (2004) detection channels of the SNO detector unambiguously demonstrated by a flavor-independent measurement of the neutrino flux that the solar neutrino problem was not due to neutrino decay. Even though neutrino decay is now known not to be the dominant effect behind the solar neutrino problem, solar neutrinos make an excellent test beam for investigating neutrino decay as a second-order effect.

Astrophysical and cosmological observations provide strong constraints on radiative decay of neutrinos (see references in the Particle Data Group (PDG) review Patrignani et al. (2016)), however constraints on nonradiative decay are much weaker. The energy ranges of solar neutrinos and the baseline between the Sun and the Earth make solar neutrinos a strong candidate for setting constraints on nonradiative decays Beacom and Bell (2002). As such, we consider only nonradiative decays where any final states would not be detected as active neutrinos Beacom and Bell (2002); Berryman et al. (2015); Picoreti et al. (2016); Choubey et al. (2000); Joshipura et al. (2002); Bandyopadhyay et al. (2003); Raghavan et al. (1988). The signal of nonradiative neutrino decay is an energy-dependent disappearance of neutrino flux that can be extracted with a statistical fit to solar neutrino data.

Previous analyses of neutrino lifetime Berryman et al. (2015); Picoreti et al. (2016) utilizing published SNO fits Aharmim et al. (2013b) were limited because the polynomial survival probabilities from Aharmim et al. (2013b) do not well capture the shape distortion of neutrino decay. Additionally, the previously published fits assumed the total flux was conserved, i.e. they inferred $P_{ea}$, the probability of detecting a solar neutrino as a $\nu_{\mu}$ or $\nu_{\tau}$ neutrino at Earth, from the constraint $P_{ee}+P_{ea}=1$ which does not apply in a decaying scenario. Both points are addressed in this analysis by implementing and fitting to a model including neutrino decay that independently calculates $P_{ee}$ and $P_{ea}$.

Precise measurements of neutrino mixing parameters from KamLAND Abe et al. (2008) and Daya Bay An et al. (2017) along with improved theoretical predictions for the 8B flux Serenelli et al. (2009) reduce the uncertainty in the underlying solar neutrino model. Using these constraints, a dedicated fit is performed to SNO data where the electron neutrino survival probabilities, $P_{ee}$ and $P_{ea}$, are calculated directly as an energy-dependent modification to the standard Mikheyev - Smirnov - Wolfenstein (MSW) Wolfenstein (1978); Mikheyev and Smirnov (1985) survival probability.

This paper is organized as follows. In Section II we discuss the SNO detector. Section III reviews the theoretical basis of the measurement. Section IV presents the analysis technique, a likelihood fit of the solar neutrino signal that includes a neutrino decay component. The results are presented in Section V, and Section VI concludes.

II The SNO Detector

The Sudbury Neutrino Observatory was a heavy water Cherenkov detector located at a depth of 2100 m (5890 m.w.e.) in Vale’s Creighton mine near Sudbury, Ontario. The detector utilized an active volume of 1000 metric tons of heavy water (D2O) contained within a 12 m diameter spherical acrylic vessel (AV). The AV was suspended in a volume of ultrapure light water (H2O) which acted as shielding from radioactive backgrounds. This ultrapure water buffer contained 9456 8-inch photomultiplier tubes (PMTs) attached to a 17.8 m diameter geodesic structure (PSUP). These PMTs recorded Cherenkov light produced by energetic particles in the active volume. The effective coverage of the PMTs was increased to 55% Jelley et al. (2009) by placing each PMT inside a non-imaging reflective light concentrator. A schematic diagram of the detector is shown in Figure 1.

SNO was sensitive to three neutrino interaction channels:

[TABLE]

The neutral current interaction (NC) couples to neutrinos of all flavors equally and allowed an unambiguous measurement of the total active neutrino flux. The charged current (CC) and elastic scattering (ES) interactions couple exclusively (CC) or preferentially (ES) to the electron flavor neutrino, which allowed the solar electron neutrino survival probability to be measured.

SNO operated in three phases, which differed in sensitivity to neutrons, and hence to the NC interaction. Phase I was the baseline detector described above in which neutrons were detected via the 6.25 MeV $\gamma$-ray released after capturing on deuterons. Phase II increased the neutron capture efficiency using the higher capture cross section of 35Cl by adding NaCl to the D2O. In addition to the increased cross section, the neutron capture on 35Cl resulted in a cascade of $\gamma$-rays summing to a higher energy of 8.6 MeV, better separating this signal from radioactive backgrounds. Phase III added a Neutral Current Detector (NCD) array inside the active volume for an independent measure of neutron production inside the detector. These NCDs were high purity nickel tubes containing 3He gas, and they were instrumented to utilize the 3He as a proportional counter for thermal neutrons Amsbaugh et al. (2007). For Phase III only there are two sources of detector data: the PMT array data as in Phase I and Phase II and the NCD array data. As these datasets are treated differently in analyses, the PMT data from Phase III will be referred to simply as Phase III with the NCD data being Phase IIIb. A combined analysis of Phase I and II data led to the first low energy measurement of the electron neutrino survival probability Aharmim et al. (2010). That analysis was later extended to incorporate Phase III data Aharmim et al. (2013b), and the analysis described in this paper was based on the analysis described in Aharmim et al. (2013b).

SNO developed a highly detailed microphysical simulation of the detector called SNOMAN Boger et al. (2000). This software could be configured to exactly reflect the experimental conditions at any particular time (for example, the values of the trigger settings during a particular run), allowing accurate Monte Carlo reproduction of the data. Monte Carlo simulations of the various signal and background events generated with statistics equivalent to many years of livetime were used extensively in this analysis. For a detailed description of this simulation package, see Aharmim et al. (2013b).

III Neutrino Decay for 8B Solar Neutrinos

This section provides the theoretical background for the analysis. We begin by reviewing ordinary solar neutrino oscillation and the MSW effect before introducing the effects of possible neutrino decay.

III.1 Neutrino oscillation

Neutrinos are produced and interact in the flavor basis, $\ket{\nu_{\alpha}}$ where $\alpha={e,\mu,\tau}$, however these are not eigenstates of the vacuum Hamiltonian, whose eigenstates (the eigenstates with definite mass) we denote as $\ket{\nu_{i}}$ where $i={1,2,3}$. The flavor basis is related to the mass basis by the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix $U_{\alpha i}$ as follows:

[TABLE]

The free evolution of these states is most easily represented in the mass basis, as these are eigenstates of the Hamiltonian:

[TABLE]

The survival probability for flavor state $\ket{\nu_{\alpha}}$ to be detected as flavor state $\ket{\nu_{\beta}}$ at some later time after free evolution for a distance L is therefore

[TABLE]

III.2 The MSW effect

The MSW Wolfenstein (1978); Mikheyev and Smirnov (1985) effect proposes that the coherent forward scattering of electron flavor neutrinos off of electrons in a material adds a potential energy, $V_{e}$, to electron flavor neutrinos which depends on the local electron density, $n_{e}$:

[TABLE]

Written in the mass basis the Hamiltonian including this effect, $H_{MSW}$, then takes the form

[TABLE]

Notably the evolution of the states in the presence of matter is now much more complicated since the eigenstates now depend on the electron density.

It is useful to introduce the matter mass basis, $\ket{\nu_{mi}(V_{e})}$, consisting of eigenstates of the Hamiltonian $H_{MSW}$ at a particular electron potential $V_{e}$. Note that if the electron density is zero, $H_{MSW}$ reduces to $H_{0}$. Therefore $\ket{\nu_{mi}(0)}\rightarrow\ket{\nu_{i}}$. In many cases the variation of $V_{e}$ is slow enough that the evolution is adiabatic, meaning some state $\ket{\nu(t)}$ has a constant probability to be one of the instantaneous matter mass states at a later time.

[TABLE]

An adiabatic approximation is made in this analysis as with previous SNO analyses Aharmim et al. (2013b). Knowing where in the Sun a neutrino is produced (or more precisely the electron density at the production point), one can calculate the eigenstate composition for as long as the adiabatic condition is satisfied. Once the neutrino reaches the solar radius, vacuum propagation dominates. As vacuum propagation does not change the mass state composition of a state, the neutrinos that arrive at Earth have the same mass state composition as those exiting the Sun. Due to the large distance between the Earth and the Sun, these mass state fluxes can be assumed to be incoherent once they arrive at Earth, and any regeneration of coherence in the Earth is ignored. Therefore, the arrival probability $\phi_{i}$ of neutrino mass state $\nu_{i}$ at Earth due to electron neutrinos produced at an electron potential $V_{e}$ in the Sun in the presence of the MSW effect can be calculated as

[TABLE]

The analytic expression for this value is non-trivial and in practice $H_{MSW}$ is numerically diagonalized in the flavor basis to find $\bra{\nu_{mi}(V_{e})}$ at a particular $V_{e}$ value and compute this projection.

III.3 Modeling a neutrino decay signal

The flux of a particular mass state, $i$, could have some lifetime associated with it, $\tau_{i}$, representing the decay of neutrinos of that mass state. Since the actual neutrino masses are currently unknown, the lifetime may be represented by an effective parameter, $k_{i}$, scaled by the mass of the state:

[TABLE]

Since the Earth-Sun distance is quite large compared to the solar radius, any decay within the sun will be ignored, and decay is only considered while propagating in vacuum from the Sun to the Earth. Here, we consider nonradiative decay to some non-active channel Beacom and Bell (2002), which manifests as disappearance of a mass state. Therefore, the arrival probability, $\psi_{i}$, of a neutrino mass state at Earth in the presence of neutrino decay can be given as

[TABLE]

where $L$ is the radius of the Earth’s orbit (1 AU) and $E$ is the energy of the neutrino. Survival probabilities for electron and muon/tau flavor neutrinos may then be recovered using the PMNS matrix in the usual way:

[TABLE]

III.4 Decay of 8B solar neutrinos

Figure 2 shows the fraction of mass state $\nu_{2}$ in the total neutrino flux as a function of energy. Considering the cross section weighted 8B neutrino energy spectrum, one finds that less than $4\%$ of the detected flux is not mass state $\nu_{2}$. As such, SNO data is dominated by $\nu_{2}$ neutrinos. This analysis is insensitive to decay of mass states $\nu_{1}$ or $\nu_{3}$, and the lifetimes $k_{1}$ and $k_{3}$ are assumed to be infinite.

The signal to be fit is therefore an energy-dependent flux disappearance due to the decay of mass state $\nu_{2}$ neutrinos. This energy dependence is distinct from the MSW effect, allowing an energy-dependent likelihood fit to distinguish between them. In the formalism presented here, decay of mass state $\nu_{2}$ is entirely described by the lifetime parameter $k_{2}$. Examples of $P_{ee}$ for various values of $k_{2}$ are shown Figure 3.

IV Analysis

We performed a likelihood fit over all three phases of SNO data for a finite neutrino lifetime, $k_{2}$, as defined in the previous section. This analysis built on the 3-phase SNO analysis Aharmim et al. (2013b) and the methods are briefly summarized here for completeness but can be found in detail in the previous publication. For each fit, many parameters were floated with constraints. These parameters include background rates, neutrino mixing parameters, and the nominal 8B flux. Systematic uncertainties found not to be strongly correlated with the solar neutrino signal were handled with a shift and refit procedure. For the final result, a likelihood profile for the parameter $k_{2}$ was generated and used to set a lower bound for that parameter. See the following sections for more detail.

IV.1 Data selection

Data selection proceeds in a number of steps. The data are organized in time periods called runs, and the first step is to select runs with nominal detector conditions. This analysis uses the same run list developed for the full analysis of all three phases of the SNO data Aharmim et al. (2013b).

There is also an event-level selection within each run. These cuts remove instrumental backgrounds, muons, and muon followers from the dataset. Again, for this analysis we use the same reconstruction corrections, data cleaning, and high-level cuts used in Aharmim et al. (2013b) for identifying physics events.

We define a region of interest for the analysis in terms of effective recoil electron kinetic energy $T_{\textrm{eff}}$ and radial position $r$, requiring $r<$ 5.5 m, and $3.5\textrm{ MeV}<T_{\textrm{eff}}<20\textrm{ MeV}$. Phase III data is included with a higher range of energies, $6.0\textrm{ MeV}<T_{\textrm{eff}}<20\textrm{ MeV}$, as in previous SNO analyses Aharmim et al. (2013b), since low energy backgrounds were not as well understood in that phase.

IV.2 Blindness

The data from all three phases of SNO were reblinded during the development of the analysis. The fit itself was developed on a statistical ensemble of Monte Carlo datasets. Once the analysis was finalized, the data were unblinded in two stages. The fit was first run on a one-third statistical subsample, to verify that it behaved as expected on real data, before proceeding to fit the full dataset.

IV.3 Fit

We developed a binned likelihood fit that combines all three phases of SNO data. For Phases I and II, we perform a fit in four observable quantities: energy, volume-weighted radius ($\rho=r^{3}/r_{AV}^{3}$), solar angle, and isotropy ($\beta_{14}$). For Phase III data, we perform a fit in three observable quantities: energy, radius, and solar angle. To incorproate Phase IIIb data, we use a constraint from the earlier pulse shape analysis Aharmim et al. (2013b) that determined the number of NCD events that could be attributed to neutrino interactions. For each of these components, the binning of the observable quantities used was that in Aharmim et al. (2013b).

For each class of signal and background events in a phase, a probability distribution function (PDF) with the correct dimensions for that phase is produced using Monte Carlo events. The likelihood of the data being described by a weighted sum of the PDFs for each class of signal and background is maximized by minimizing the negative logarithm of this likelihood with MINUIT James and Roos (1975). The construction of this likelihood function is identical to what is described in the SNO 3-phase analysis Aharmim et al. (2013b) with one exception: the polynomial survival probability from previous SNO analyses is replaced with the survival probability parameterized by the physical quantities described in Section III.

IV.4 Solar Signal

The following sections discuss the inputs to modeling the flux of 8B solar neutrinos as detected by SNO.

IV.4.1 Standard Solar Model

The neutrino model implemented here uses the radial distribution of electron density and radial distribution of the 8B neutrino flux calculated in the BS05(OP) Standard Solar Model (SSM) Bahcall et al. (2005). Uncertainties in these values are not quoted in the original source and are therefore not considered in this fit. These predictions are expected to be uncorrelated with $\nu_{2}$ decay as they are not determined with neutrino measurements.

As earth-bound measurements of the solar neutrino flux would be biased by neutrino decay, a theoretical prediction for the 8B flux is required. Serenelli’s most recent prediction Serenelli et al. (2009) yields a 8B flux of $5.88\times 10^{6}$ cm*-2s-1* with $11\%$ uncertainty which is used as a prior in this fit. For reference the flux from BS05(OP) Bahcall et al. (2005) is $5.69\times 10^{6}$ cm*-2s-1*.

IV.4.2 Neutrino Mixing

Neutrino mixing parameters taken from KamLAND Abe et al. (2008) and Daya Bay An et al. (2017) are reproduced in Table 1. Parameters from KamLAND and Daya Bay were used in this analysis to avoid biasing the result by using values correlated with previous SNO analyses. As these measurements were done with neutrinos produced on Earth, they are expected to be uncorrelated with effects of $\nu_{2}$ decay given existing constraints on neutrino decay. The current limit on $k_{2}$ constrains it to be $>7.2\times 10^{-4}$ s/eV Picoreti et al. (2016) which means at length scales comparable to the diameter of the Earth, the maximum flux fraction lost by $\nu_{2}$ decay for a $10$ MeV neutrino is given by $1-e^{-2R_{earth}/(Ek_{2})}\approx 6\times 10^{-6}$. Such a small fractional loss would have negligible impact on values quoted for mixing parameters. These parameters and their central values are used as priors and floated during the fit.

IV.5 Backgrounds

Besides instrumental backgrounds, which can be easily removed with cuts based on event topology, the main sources of background events are radioactive backgrounds and atmospheric neutrino interactions. A summary of the sources of these and other backgrounds is given in this section.

For Phases I and II, radioactive decays of 214Bi (from uranium and radon chains) and 208Tl (from thorium chains) produce both $\beta$-particles and $\gamma$-rays with high enough energies to pass event selection criteria. In Phase III the lower energy bound was high enough to exclude these backgrounds. The inner D2O, acrylic vessel, and outer H2O volumes are treated as separate sources of radioactive decays due to differing levels of contamination. The PMT array is another source of radioactivity and, despite its increased distance from the fiducial volume, is the dominant source of low-energy backgrounds.

Relevant to all three phases, $\gamma$-rays above $2.2$ MeV may photodisintegrate deuterium resulting in a neutron background. Radon daughters present on the acrylic vessel since construction result in additional neutron backgrounds from ($\alpha$,$n$) reactions on carbon and oxygen in the acrylic. In Phase II the addition of NaCl resulted in a 24Na background from neutron captures on 23Na. 24Na decay produces a $\gamma$-ray with high enough energy to photodisintegrate a deuteron, increasing the neutron background in Phase II. In Phase III the addition of the NCD array inside the acrylic vessel brought additional radioactive backgrounds. Primarily this resulted in an increase of photodisintegration events throughout the detector. Two NCDs with higher levels of radioactivity were treated separately in the analysis.

Additional backgrounds include solar $hep$ neutrinos and atmospheric neutrinos. The $hep$ neutrinos have a higher endpoint than 8B neutrinos, however the predicted flux is approximately a thousand times less Bahcall et al. (2005). The flux of $hep$ neutrinos is fixed to the standard solar model rate in this analysis. Atmospheric neutrinos also have a relatively low flux, and the rate is fixed to results from previous SNO analyses Aharmim et al. (2010).

Finally, there is a class of instrumental background that tends to reconstruct on the acrylic vessel. For Phase III these instrumentals are easily cut in event selection as they were well separated from physics events in the $\beta_{14}$ parameter. Near the lower energy threshold in Phase I and Phase II these events were not as well separated in $\beta_{14}$ resulting in some contamination Aharmim et al. (2010), and this event class was therefore included in the fit for Phase I and Phase II.

For further details on how backgrounds were included in the fit and which in-situ and ex-situ constraints were used, see Appendix B in Aharmim et al. (2013b).

IV.6 Systematics

Parameters that shift, rescale, or affect the resolution of observables used in the fit are treated as systematic uncertainties. Other systematic uncertainties include: parameters that control the shape of the analytic PDF for PMT $\beta-\gamma$ backgrounds, photodisintegration efficiency, and neutron capture efficiencies.

The neutron capture efficiency was found to be strongly correlated with the neutrino parameters and is floated in the fit to correctly account for correlations with the final results.

Less correlated parameters that are well constrained by the data, such as the parameters for the analytic PDF for PMT $\beta-\gamma$ backgrounds, are scanned as an initial step. Each of these systematic parameters is scanned independently with other systematic parameters held fixed while profiling out all floated parameters. This scan produces a likelihood profile, which is fit by an asymmetric Gaussian to determine the central value and uncertainty of the scanned parameter. After a parameter is scanned, its central value and uncertainty is updated to the fit result before scanning the next parameter. This process is repeated until the central values for each parameter stabilize to ensure the global minimum is found. The final central values are retained and fixed during MINUIT minimization, and the impact of their uncertainty on the uncertainty of floated parameters is evaluated with the shift-and-refit procedure described below.

The least correlated parameters that are not well constrained by the data are fixed to predetermined nominal values during the fit, and the impact of their uncertainty on each floated parameter is evaluated with a shift-and-refit method.

The shift-and-refit method draws $10^{6}$ sets of systematic parameters from their respective asymmetric Gaussian distributions. The fit is then re-run many times with the systematic parameters fixed to each of the generated sets. This produces distributions of fitted values for the parameters floated in the fit. The widths of these distributions are taken to represent the systematic uncertainty on the floated parameters.

For a full listing of the systematic uncertainties and how they were handled, see Appendix B in Aharmim et al. (2013b).

IV.7 Bias and pull testing

Significant testing was done on Monte Carlo datasets to ensure the statistical robustness of the fit. For all tests in this section, the solar signal was generated with an assumed $k_{2}$ value of $10^{-4}$ s/eV to test the sensitivity near existing limits, and all other parameters were chosen by randomly sampling the prior distributions for each fake dataset.

The fit was run on each dataset including only the signals and backgrounds for that stage. The fitted values of all floated parameters were recorded and used to produce bias and pull distributions. No significant bias was found, and pull widths were found to be consistent with expectations.

V Results

Figure 4 shows the likelihood profile of $k_{2}$ both with the systematic parameters fixed to central values and with the systematic uncertainties included. The likelihood profile incorporating systematic uncertainties is generated by assuming the shape of the likelihood profile does not change as the systematic parameters vary, but rather simply shifts according to the shift in the fitted value of $k_{2}$ from the shift-and-refit method. Therefore, the systematic uncertainties are included by shifting the fixed systematic profile by each shift in the shift-and-refit distribution for $k_{2}$ and averaging the likelihood at each point.

A shallow minimum at $3.45^{+5.50}_{-1.68}\times 10^{-4}$ s/eV is found, however the upper uncertainty is consistent with infinite lifetime at confidences greater than $85\%$, meaning this analysis is not a significant measurement of neutrino decay. Using Wilks’ theorem Wilks (1938), a lower bound for $k_{2}$ can be set at $k_{2}>8.08\times 10^{-5}$ s/eV at $90\%$ confidence.

V.1 Comparison to previous SNO analyses

The best fit 8B neutrino flux from this analysis is $6.08^{+0.47}_{-0.47}$(stat.)${}^{+0.21}_{-0.22}$(syst.)$\times 10^{6}$ cm*-2s-1* and has slight tension with results of the previous SNO 3-phase analysis: $5.23^{+0.16}_{-0.16}\times 10^{6}$ cm*-2s-1* Aharmim et al. (2013b) where statistical and systematic uncertainties have been combined. With $k_{2}$ fixed to infinite lifetime this analysis results in a 8B neutrino flux of $5.22^{+0.16}_{-0.16}\times 10^{6}$ cm*-2s-1*, again with systematic and statistical uncertainties combined, that is in very good agreement with previous results. The uncertainty with $k_{2}$ allowed to float is much larger due to the additional freedom of neutrino decay in the model and the fact that the lifetime is strongly anti-correlated with the flux. These two parameters are not degenerate only because the effect of neutrino decay is energy-dependent, and this fit to the neutrino energy spectrum can capture that effect. To that end we expect the uncertainty on the 8B neutrino flux from this analysis to be larger than previous analyses.

V.2 Combined analysis results

Any experiment measuring a solar flux can be compared to a standard solar model to constrain neutrino lifetimes. Likelihood profiles of $k_{2}$ generated for other solar experiments can be combined with the profile from this analysis to arrive at a global limit. Particularly, experiments sensitive to lower energy solar neutrinos, such as the $pp$ or 7Be solar neutrinos, can provide strong constraints on $k_{2}$ as the $L/E$ for these neutrinos is greater than for the 8B neutrinos.

To incorporate the results from other experiments, the measured flux reported by an experiment assuming a flux of only electron neutrinos (i.e. $P_{ee}=1$), $\Phi_{e}$, is converted to an total inferred flux, $\Phi_{T}$, by way of a neutrino model that predicts $\overline{P_{ee}}$ and $\overline{P_{ea}}$, the average survival probabilities for that flux, and the relative cross sections, $\sigma_{a}/\sigma_{e}$, where $\sigma_{e}$ is the cross section for electron-flavor neutrinos and $\sigma_{a}$ the cross section for all other neutrino flavors (i.e. $\nu_{\mu}$ and $\nu_{\tau}$)

[TABLE]

The neutrino decay model described in Section III is used, and the averaging is done over the flux-appropriate standard solar model production regions (electron density) Bahcall et al. (2005) and experiment-appropriate cross section weighted neutrino spectra Bahcall et al. (2005).

This $\Phi_{T}$ can be directly compared to standard solar model predictions with the following likelihood term

[TABLE]

where $\sigma_{T}$ and $\sigma_{SSM}$ are the uncertainties on the inferred flux, $\Phi_{T}$, and standard solar model flux, $\Phi_{SSM}$. The mass state $\nu_{1}$ lifetime, $k_{1}$, is a free parameter in the fit and profiled over in producing the final limit on $k_{2}$, as lower energy solar neutrinos may contain significant fractions of $\nu_{1}$. The mass state $\nu_{3}$ lifetime, $k_{3}$, remains fixed to infinity, as all solar neutrinos contain negligible amounts of $\nu_{3}$. The neutrino mixing parameters are constrained as described in Section IV.4.2 and allowed to float.

Following this methodology, a profile for $k_{2}$ is generated using: Super-K Abe et al. (2016), KamLAND Abe et al. (2011b) and Borexino Agostini et al. (2017) 8B results; Borexino Bellini et al. (2011) and KamLAND Gando et al. (2015) 7Be results; the combined gallium interaction rate from GNO, GALLEX, and SAGE Abdurashitov et al. (2009); and the chlorine interaction rate from Homestake Cleveland et al. (1998). For both chlorine and gallium, the predicted interaction rate is computed following the procedure in Section V of Abdurashitov et al. (2009), but using the neutrino model and mixing constraints used elsewhere in this paper. These predicted rates were compared to the measured rates with likelihood terms analogous to Equation 12.

The final profile, combined with this analysis of SNO data, is shown in Figure 5 and constrains $k_{2}$ to be ${>1.04\times 10^{-3}}$ s/eV at $99\%$ confidence.

VI Conclusion

Neutrinos are known to have mass, allowing for potential decays to lighter states. However, analyses of solar neutrino data assuming the MSW solution to the solar neutrino problem are consistent with a non-decaying scenario. By analyzing the entire SNO dataset, using a model that predicts the survival probability of electron-type solar neutrinos allowing for the decay of mass state $\nu_{2}$, we were able to set a limit on the lifetime of neutrino mass state $\nu_{2}$: ${k_{2}>8.08\times 10^{-5}}$ s/eV at $90\%$ confidence. Combining this with measurements from other solar experiments results in a new best limit of ${k_{2}>1.04\times 10^{-3}}$ s/eV at $99\%$ confidence. The improvement from the previous limit, $k_{2}>7.2\times 10^{-4}$ s/eV at $99\%$ confidence Picoreti et al. (2016), was driven by the inclusion of additional solar flux measurements and updated analyses of previously considered experiments.

Acknowledgements.

This research was supported by: Canada: Natural Sciences and Engineering Research Council, Industry Canada, National Research Council, Northern Ontario Heritage Fund, Atomic Energy of Canada, Ltd., Ontario Power Generation, High Performance Computing Virtual Laboratory, Canada Foundation for Innovation, Canada Research Chairs program; US: Department of Energy Office of Nuclear Physics, National Energy Research Scientific Computing Center, Alfred P. Sloan Foundation, National Science Foundation, the Queen’s Breakthrough Fund, Department of Energy National Nuclear Security Administration through the Nuclear Science and Security Consortium; UK: Science and Technology Facilities Council (formerly Particle Physics and Astronomy Research Council); Portugal: Fundação para a Ciência e a Tecnologia. This research used the Savio computational cluster resource provided by the Berkeley Research Computing program at the University of California, Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor for Research, and Chief Information Officer). We thank the SNO technical staff for their strong contributions. We thank INCO (now Vale, Ltd.) for hosting this project in their Creighton mine.