Statistical Evaluation of Experimental Determinations of Neutrino Mass Hierarchy

Neutrino mass hierarchy (MH), i.e. whether the mass of the third generation neutrino ($\nu_{3}$ mass eigenstate) is greater or less than the masses of the first and the second generation neutrinos ($\nu_{1}$ and $\nu_{2}$), is one of the main questions to be answered in the standard model. Besides its fundamental importance to neutrino oscillation physics, the resolution of the neutrino MH plays an important role for the search of neutrinoless double-beta decay, which would determine whether neutrino is a Dirac or Majorana fermion. With the recent discovery of a large value of $\sin^{2}2\theta_{13}$ from Daya Bay [1, 2, 3, 4], T2K [5], MINOS [6], Double Chooz [7], and RENO [8], the stage for addressing the neutrino MH has been set. It became one of the major goals of current and next generation long baseline neutrino experiments (T2K [9], NO$\nu$A [10] and LBNE [11]) and atmospheric neutrino experiments (Super-K [12], MINOS [13], PINGU [14], and INO [15]). Meanwhile, the idea of utilizing a reactor neutrino experiment to determine the MH is also intensively discussed [16, 17, 18, 19, 20].

The objective of this paper is to present appropriate ways to do statistical analysis that will help determine the neutrino mass hierarchy. We start by introducing a few symbols and state the physics problem in terms of a pair of statistical hypotheses. Let $m_{1}$, $m_{2}$ and $m_{3}$ denote the masses of the $\nu_{1}$, $\nu_{2}$ and $\nu_{3}$ mass eigenstate neutrinos, and let $\Delta m_{ij}^{2}\equiv m_{i}^{2}-m_{j}^{2}$ for $i,j=1,2,3$. As reviewed in Ref. [21], it is known that $\Delta m_{21}^{2}>0$ from measurements of solar neutrinos given the definition of mixing angle $\theta_{12}$. Whereas the sign of $\Delta m_{32}^{2}$ is so far unknown, and it’s common to use NH and IH to denote the two hypotheses, the normal hierarchy and the inverted hierarchy, respectively:

A unique feature to the above hypotheses testing problem is that, there are additional, rather strong information regarding the parameter $\Delta m_{32}^{2}$ that need to be taken into account properly. Actually, based on previous experiments, a $68\%$ confidence interval of $M^{2}_{32}\equiv|\Delta m^{2}_{32}|$ is given by $(2.43\pm 0.13)\times 10^{-3}$ eV² [22].

We will mainly address two aspects of the hypotheses testing problem. The first one concerns conducting a test after data has been collected. We discuss a classical testing procedure based on a $\Delta\chi^{2}$ statistic (Eq. 3), or equivalently, the procedure of constructing confidence intervals by inverting the test. As a matter of fact, the classical procedure is derived upon the assumption that the best estimator of $\Delta m^{2}_{32}$ based on experimental data would follow a distribution that is approximately Gaussian. But due to existing constraints on $M^{2}_{32}$, this assumption is far from being satisfied. Consequently, actual levels of the resulting confidence intervals may deviate substantially from their nominal levels, as we demonstrate in Sec. II. Instead, a general way to construct confidence intervals that are true to their nominal levels is the Feldman-Cousins approach [23], which we also illustrate in detail in Sec. II.

Still, there is a fundamental limitation to the use of confidence intervals. Note that in the MH determination problem, one of the most crucial questions is, what is the chance that the MH is indeed NH (or IH) given the observed experimental data? Classical confidence intervals are not meant to answer this question directly, whereas credible intervals reported by a Bayesian procedure is. In Sec. III, we present a Bayesian approach, which effortlessly incorporates prior information on $M^{2}_{32}$ and output the easy-to-understand (posterior) probability of NH and IH to conclude the test. We will emphasize the importance to differentiate the Bayesian credible interval from the classical confidence interval.

The second aspect of the hypotheses testing problem that we address concerns assessment of experiments in their planning stage. It is critical to evaluate the “sensitivity” of a proposed experiment, i.e., its capability to distinguish NH and IH. Since this evaluation is performed before data collection, it has to be based on potential data from the experiment. An existing evaluation method (such as employed in [24, 11, 25]) assumes that the most typical data set under one hypothesis, say NH, happens to have been observed. Such a data set is referred to as the Asimov data set [26]. The method then calculates $\overline{\Delta\chi^{2}}$, which stands for the statistic $\Delta\chi^{2}$ in Eq. 3, with the extra bar indicating its dependence on the Asimov data set. It can be seen that $\overline{\Delta\chi^{2}}$ reflects how much the Asimov data set under NH disagrees with the alternative model, IH. It is then common practice to quantify the amount of disagreement by finding the p-value corresponding to $\overline{\Delta\chi^{2}}$ after comparing it to the quantiles of a chi-square distribution with one degree of freedom (choice of MH). Finally, one minus this p-value is sometimes reported as a quantitative assessment of the experiment. We will show in Sec. II that the comparison of the value of $\overline{\Delta\chi^{2}}$ to the quantiles of a chi-square distribution is not justified, when previous knowledge impose constraints on the range of possible values of the parameter $\Delta m^{2}_{32}$.

As an alternative solution, we adopt a Bayesian framework and develop a set of new metrics for sensitivity to evaluate the potential of experiments to identify the correct hypothesis.

The paper is organized as follows. In Sec. II, we review the steps to construct classical confidence intervals for the parameter $\Delta m^{2}_{32}$. In Sec. III, we describe a Bayesian approach that reports the probability of each hypothesis of MH given observed data set. We further extend this Bayesian method to help assess the sensitivity for future experiments. In Sec. IV, we illustrate the Bayesian approach for a simplified version of the MH problem. In particular, analytical formula of the approximations for the probability of the hypotheses, and those for the sensitivity metrics are provided. Also, a numerical comparison is made between the $\overline{\Delta\chi^{2}}$ based on the Asimov data set and the sensitivity metrics based on the Bayesian approach. Finally, discussions and a summary are presented in Sec. V and Sec. VI, respectively.

In this section, we review a classical statistical procedure of forming confidence intervals. For the problem of determining the neutrino mass hierarchy, we demonstrate that the procedure is valid in one scenario, but fails in another where known constraints on $M^{2}_{32}$ are taken into consideration. In the latter case, the Feldman-Cousins method [23] based on Monte Carlo (MC) simulation is recommended to obtain valid confidence intervals.

Consider a spectrum that consists of $n$ energy bins. Assume that the expected number of counts in each bin is a function of $\Delta m^{2}_{32}$ and a nuisance parameter $\eta$. For simplicity, we denote $\Delta m^{2}_{32}$ by $\theta$. Then for the $i$th bin, let $\mu_{i}(\theta,\eta)$ and $N_{i}$ represent the expected and the observed counts of neutrino induced reactions, respectively. When $\mu_{i}$ is large enough, the distribution of $N_{i}$ can be well approximated by a Gaussian distribution with mean $\mu_{i}$ and standard deviation $\sqrt{\mu_{i}}$.

Once the data $x=\{N_{i},i=1,\ldots,n\}$ are observed, the deviations from the expected values $\{\mu_{i}(\theta,\eta),i=1,\ldots,n\}$ are often calculated to help measure the implausibility of the parameter $(\theta,\eta)$. Specifically, when the systematic uncertainties are omitted, and that certain available knowledge concerning the parameters $\theta$ and $\eta$ are taken into consideration, one useful definition of the deviation is given by

Here, the general notation $\delta w$ represents the standard deviation of a variable $w$. So $\delta N_{i}=\sqrt{\mu_{i}}$, and the corresponding $\chi^{2}_{stat}$ term is called the Pearson’s chi-square. Also, note that $|\theta|=M^{2}_{32}$, and it is taken from [22] that $|\theta_{0}|=2.43\times 10^{-3}$ eV² and $\delta|\theta|=0.13\times 10^{-3}$ eV².

Based on Eq. 2 and a standard procedure discussed in Ref. [22], confidence intervals can be obtained for the parameter of interest $\theta$ ($\Delta m^{2}_{32}$), the sign of which is an indicator of the neutrino MH. First, define $\theta_{\min}$ to be the best fit to the data in the sense that $(\theta_{\min},\eta_{\min})=\arg\min_{\theta,\eta}\chi^{2}(\theta,\eta)$ where the minimum is taken over $\Theta\times H$, the space of all possible values of $(\theta,\eta)$. Here, the general notation $\arg\min_{w}h(w)$ denotes the value of $w$ which corresponds to the minimum of the given function $h$. Note that $\theta_{\min}$ suggested by the observed data set will not be exactly the true value of the parameter $\theta$, and a repetition of the experiment would yield a data set that corresponds to a different $\theta_{\min}$. So instead of reporting only $\theta_{\min}$, it is more rational to report a set of probable values of $\theta$ that fit the observed data not too much worse than that of the best fit, and state how trust worthy the set is. Indeed, for any given $\theta$, let $\eta_{\min}(\theta)=\arg\min_{\eta}\chi^{2}(\theta,\eta)$, and define

then a level $a$ confidence interval based on Eq. 3 is defined to be

where we use the standard set-builder notation $\{h(w):\text{restriction $w$}\}$ to denote a set that is made up of all the points $h(w)$ such that $w$ satisfies the restriction to the right of the colon. The key in constructing Eq. 4 is to specify the correct threshold value $t_{a}$ for a given confidence level $a$. (See the final paragraph of this section for a more detailed description of what confidence level means.) Most commonly examined confidence levels use $a=68.27\%(1\sigma)$, $95.45\%(2\sigma)$, $99.73\%(3\sigma)$, which are often linked to threshold values $t_{a}=$ 1, 4, 9 respectively [22]. Note that these three values are the $68.27\%$, $95.45\%$ and $99.73\%$ quantiles of the chi-square distribution with one degree of freedom, respectively. They are used as threshold values because the parameter space $\Theta$ is of dimension one and that, under certain regularity conditions, $\Delta\chi^{2}_{min}(\theta)$ would follow approximately a chi-square distribution with one degree of freedom when $\theta$ is the true parameter value. This procedure and its extensions to cases where $\theta$ is of higher dimension have been successfully applied in many studies [27, 28, 29, 30, 31, 32, 24, 11, 25, 33] in order to constrain various parameters in the neutrino physics.

Although this procedure has been widely used in analyzing experimental data, note that it is not universally applicable. Its limitations has been addressed by Feldman and Cousins [23]. Below, we illustrate this point through a simple MC simulation study. It will be shown that, in a situation that is similar in nature to the MH determination problem in Eq. 1 where there exist special constraints on the possible values of $\theta$, the aforementioned threshold values based on chi-square approximation could result in bad confidence intervals. That is, the actual coverage probabilities of the intervals strongly disagree with their nominal levels.

In the simulation, we set $n=10$, $\mu_{i}(\theta)=1000+15\cdot\theta$ for $i=1,\ldots,n$. (Here, no nuisance parameter $\eta$ is introduced, and all the expected bin counts are assumed equal for simplicity. Nevertheless, these assumptions are not essential to the purpose of our simulation.) The following two cases are investigated:

• •

  Case I: $\Theta=(-\infty,\infty)$,

• •

  Case II: $\Theta=\{-1,1\}$.

Case I is a typical situation where nothing was known about $\theta$ before the current experiment, whereas case II is designed to imitate the situation where existing measurements of $|\theta|=M^{2}_{32}$ are very accurate at around $2.43\times 10^{-3}$ eV², and we simply denoted this value to $1$ for clarity of presentation. Further, the definition of deviation analogue to Eq. 2 is taken to be $\chi^{2}(\theta)=\sum_{i}\frac{(N_{i}-\mu_{i}(\theta))^{2}}{\mu_{i}(\theta)}$ for case I. For case II, the chi-square definition is $\chi^{2}(\theta)=\sum_{i}\frac{(N_{i}-\mu_{i}(\theta))^{2}}{\mu_{i}(\theta)}+\frac{(|\theta|-|\theta_{0}|)^{2}}{(\delta|\theta|)^{2}}$ with experimental constrains on $|\theta|$. It is then reduced to $\chi^{2}(\theta)=\sum_{i}\frac{(N_{i}-\mu_{i}(\theta))^{2}}{\mu_{i}(\theta)}$ with $\theta$ being only 1 or -1.

Under each case, we set the true value of $\theta$ to be $\theta_{0}=1$, based on which 100,000 MC samples are simulated, denoted by $\{N_{1}^{(j)},\cdots,N_{10}^{(j)}\}$ for $j=1,\ldots,100000$. Then for the $j$th sample, confidence intervals of levels $a=68.27\%$, $95.45\%$, $99.73\%$ are constructed according to Eq. 4 using threshold values 1, 4, 9, respectively. Finally, at each of the three levels, we record the proportion of confidence intervals out of the 100,000 that include the truth $\theta_{0}=1$. The results are reported in the last three columns of Table 1. It can be seen that, in case I, the actual coverage probabilities closely match the nominal levels. However, in case II, the actual coverage probabilities are always higher!

Without too much technical detail, we try to explain the reason why the chi-square procedure produced valid confidence intervals for case I, but not for case II. In general, having observed data $x$ from a parametric model $P(x|\theta)$, a sensible test for a pair of hypotheses, $H_{0}:\theta\in\Theta_{0}$ and $H_{1}:\theta\in\Theta-\Theta_{0}$ (the counterpart of $H_{0}$), is the likelihood ratio test that is based on the test statistic

where $\theta_{0,\min}=\arg\min_{\{\theta\in\Theta_{0}\}}P(x|\theta)$, and $\theta_{\min}=\arg\min_{\{\theta\in\Theta\}}P(x|\theta)$ are the best fit over the null parameter set $\Theta_{0}$ and the full parameter set $\Theta$, respectively. If the observed data $x$ yields a large $\Delta\chi^{2}_{min}$, it means that $\Theta_{0}$ is implausible, which further leads to the rejection of $H_{0}$. Note that, the statistic $\Delta\chi^{2}_{min}(\theta)$ in Eq. 3 is a special case of Eq. 5 with $\Theta_{0}$ consisting of a single point, $\theta$.

In order to determine the correct threshold values in rejecting, or equivalently, in constructing confidence intervals defined by Eq. 4, the distribution/quantiles of $\Delta\chi^{2}_{min}(\theta_{0})$ considering all possible data set needs to be known, when the true parameter value is some $\theta_{0}\in\Theta_{0}$. An important result in statistics, Wilks Theorem [34, 35] states that, under certain regularity conditions, $\Delta\chi^{2}_{min}(\theta_{0})$ follows approximately a chi-square distribution with degree of freedom equal to the difference between the dimension of $\Theta$ and that of $\Theta_{0}$, when the data size is large. (In our problem, the data size is simply $\sum_{i}N_{i}$.) The main regularity conditions are, as we quote [35], “the model is differentiable in $\theta$ and that $\Theta_{0}$ and $\Theta$ are (locally) equal to linear spaces”. Essentially, such conditions imply that $\theta_{\min}$ follows an approximately Gaussian distribution centered at the true $\theta$ value, which eventually implies an approximate chi-square distribution for $\Delta\chi^{2}_{min}(\theta_{0})$.

In case I of our simulation, the best estimation of $\theta$ can be calculated directly from the number of events in each bin: $\theta_{\min}=\big[\sqrt{\sum^{n}_{i=1}N_{i}^{2}/n}-1000\big]/15$   ¹¹ 1 Note that the above $\theta_{\min}$ can be closely approximated by $\big[(\sum^{n}_{i=1}N_{i})/n-1000\big]/15$, which is indeed the exact maximum likelihood estimator for $\theta$ had we assumed that each count $N_{i}$ follows a Poisson distribution with mean $\mu_{i}(\theta)=1000+15\cdot\theta$. .The aforementioned regularity conditions are satisfied in this case, and the distribution of $\theta_{\min}$ and that of $\Delta\chi^{2}_{min}(\theta_{0})$ follow approximately the Gaussian and the chi-square distribution respectively as what Wilks theorem predicts. In the top two panels of Fig. 1, we reconfirm this fact by comparing their histograms based on the 100,000 MC samples (black shaded area) to the probability density function of the Gaussian and the chi-square distribution (blue long dash-dotted line). On the other hand, the full parameter space in case II consists of two isolated points and clearly violates the conditions required by Wilks theorem. Indeed, in case II, the best estimation of $\theta$ is given by

follows a Bernoulli distribution, and $\Delta\chi^{2}_{min}(\theta_{0})$ follows a distribution quite different from a canonical chi-square distribution. Approximations to the actual distributions of $\Delta\chi^{2}_{min}(\theta_{0})$ and $\theta_{min}$ can be obtained from the 100,000 MC samples, and are shown (black shaded area) in the bottom two panels of Fig. 1. Further, an analytical approximation (red dash-dot-dotted line) to the distribution is derived in Appendix A. The analytical calculation implies that, independent of whether the truth $\theta_{0}$ is $1$ or $-1$, the p-value ²² 2 The p-value at $t$ is defined to be the percentage of potential measurements that result in the same or a more extreme value of the test statistic, say $\Delta\chi^{2}_{min}$, than $t$. corresponding to an observed value of $\Delta\chi^{2}_{min}(\theta_{0})$, say $t$, is approximately given by

for any $t>0$; and the p-value is 1 for any $t\leq 0$. Here erf is the Gaussian error function: $\text{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}dt$. (We use the general notation $P(A)$ to denote the probability of an event $A$.)

The discussions above suggest that, when constructing confidence intervals in special cases where conditions of Wilks theorem do not hold (or that the user can not be sure if the conditions hold), the regular threshold values (such as $t_{\alpha}=$ 1, 4, 9 mentioned earlier) should not be taken for granted. Instead, alternative thresholds based on MC or case-specific analytical approximations are needed. We recommend using the MC method with a large MC sample size whenever possible, because unlike other methods, it is guaranteed to produce a valid confidence interval for $\theta$. We hereby review how to produce a valid 1-$\sigma$ (68.27%) confidence interval for $\theta$ using the MC method [23]. This method can easily be generalized to build confidence intervals of any level.

• •

  Having observed data $x=\{N_{1},\cdots,N_{n}\}$, apply the following procedure to every $\theta$ in the parameter space $\Theta$ (fix one $\theta$ at a time):

  1. 1.

     Calculate $\Delta\chi^{2}_{min}(\theta)^{x}$ with Eq. 3 based on the observed data.

  2. 2.

     Simulate a large number of MC samples, say $\{x^{(j)}\}_{j=1}^{T}$, where $x^{(j)}=\{N_{1}^{(j)},\cdots,N_{n}^{(j)}\}$ is generated from the model with true parameter value $\theta$. For $j=1,\ldots,T$, calculate $\Delta\chi^{2}_{min}(\theta)^{(j)}$, that is Eq. 3 based on the $j$th MC sample $x^{(j)}$. This produces an empirical distribution of the statistic $\Delta\chi^{2}_{min}(\theta)$.

  3. 3.

     Calculate the percentage of MC samples such that $\Delta\chi^{2}_{min}(\theta)^{(j)}<\Delta\chi^{2}_{min}(\theta)^{x}$. Then $\theta$ is included in the 1-$\sigma$ confidence interval if and only if the percentage is smaller than 68.27%.

One can easily check that p-values analytically obtained from Eq. 6 for case II (basically the MC method) are consistent with the simulation results listed in Table 1.

On a separate issue that was also emphasized in Ref. [23], classical confidence intervals should not be confused with Bayesian credible intervals. On one hand, the confidence-level of a confidence interval, say $a$, is an evaluation of this interval estimation procedure based on many potential repetitions of the experiment. More specifically, had the experiment been independently repeated $100$ times, applying the estimation procedure to each would result in $100$ intervals, and $a$ represents the proportion of these intervals that we expect to contain the true value of the unknown parameter $\theta$. The level-$a$ confidence interval reported in practice is the result of applying such a procedure to the data observed in the current experiment. On the other hand, a Bayesian credible interval, say of credible-level $b$, is a region in the parameter space such that, given the observed data, it contains the true value of the unknown parameter with probability $b$. In general, an $a$-level confidence interval does not coincide with an $a$-level Bayesian credible interval. In other words, if $C_{a}$ is an $a$-level confidence interval built from the observed data $x$, then it is generally inappropriate to give the interpretation that $P(\theta\in C_{a}|x)$ (the probability of true $\theta$ inside $C_{a}$ given data $x$) is $\alpha$. Nevertheless, in Appendix B, we discuss when confidence intervals approximately match Bayesian credible intervals. In the next section, we present a Bayesian approach to the problem of determining neutrino mass hierarchy.

In this section, we consider a situation where $\theta$ can only take on two possible values, $1$ and $-1$, which correspond to the hypotheses NH and IH respectively. This simplified setting is motivated by the fact that existing measurements of $|\theta|=M^{2}_{32}$ are very accurate at around $2.43\times 10^{-3}$ eV², and we simply denote this value to $1$ for clarity of presentation. It is a special case of the Bayesian treatment in the previous section, where $P(\theta|NH)$ and $P(\theta|IH)$ are assigned degenerate distributions at $1$ and $-1$ respectively. That is, $P(\theta=1|NH)=P(\theta=-1|IH)=1$. Further, there is no nuisance parameter $\eta$. As a result, the expected bin counts will be denoted by $\mu_{i}^{NH}=\mu(1)$ and $\mu_{i}^{IH}=\mu(-1)$ respectively.

Below, we showcase numerical calculations of various sensitivity criteria for this example. In particular, we introduce approximations that are simple functions of a term commonly known as “$\overline{\Delta\chi^{2}}$ of Asimov data set” in the physics literature. According to the definition in Ref. [26], “the Asimov data set” under hypothesis MH is given by $x^{MH}=(\mu_{1}^{MH},\cdots,\mu_{N}^{MH})$, where $\mu_{i}^{MH}=\mu_{i}(\theta_{0}^{MH},\eta_{0}^{MH})$ and $(\theta_{0}^{MH},\eta_{0}^{MH})=\arg\max_{(\theta,\eta)}P(\theta,\eta|MH)$ is the prior mode under MH. In words, the Asimov data set is the most typical data set under the most likely parameter values based on prior knowledge subject to the given model.

Interestingly, $\overline{\Delta\chi^{2}}$ is itself often used as a measure of sensitivity. Here, we’ll contrast the typical usage of $\overline{\Delta\chi^{2}}$ to that of the sensitivity criteria developed in the previous section. More accurate evaluations of these sensitivity criteria are also attainable via MC methods.

Suppose that the proposed experiment will collect enough data such that the expected counts under NH and IH are much larger than the difference between them: $\mu_{i}^{NH}\sim\mu_{i}^{iH}>>|\mu_{i}^{NH}-\mu_{i}^{iH}|$. Using the notations introduced in Sec. II, if the nature is NH, then the observed counts $N_{i}$ can be represented as

where $g_{1},\cdots,g_{n}$ are mutually independent standard Gaussian random variables. Then, the statistic $\Delta\chi^{2}$ of Eq. 11 becomes

Here, the subscript $T=NH$ indicates that the nature is NH. Since $\mu_{i}^{iH}>>|\mu_{i}^{NH}-\mu_{i}^{iH}|$, the summation of the last two terms in Eq. 21 is negligible as it is approximately $\sum_{i=1}^{n}\frac{\mu_{i}^{NH}-\mu_{i}^{IH}}{\mu_{i}^{IH}}\cdot(g_{i}^{2}-1)$ by a Taylor expansion of the last term. Therefore, $\Delta\chi^{2}_{T=NH}$ follows a Gaussian distribution, with mean and standard deviation: