Quantifying the sensitivity of oscillation experiments to the neutrino mass ordering

Determining the type of the neutrino mass ordering (normal versus inverted) is one of the most important open questions in neutrino physics. In this paper we clarify the statistical interpretation of sensitivity calculations for this measurement. We employ standard frequentist methods of hypothesis testing in order to precisely define terms like the median sensitivity of an experiment. We consider a test statistic $T$ which in a certain limit will be normal distributed. We show that the median sensitivity in this limit is very close to standard sensitivities based on $\Delta\chi^2$ values from a data set without statistical fluctuations, such as widely used in the literature. Furthermore, we perform an explicit Monte Carlo simulation of the INO, JUNO, LBNE, NOvA, and PINGU experiments in order to verify the validity of the Gaussian limit, and provide a comparison of the expected sensitivities for those experiments.


INTRODUCTION
The ordering of neutrinos masses constitutes one of the major open issues in particle physics. The mass ordering is called "normal" ("inverted") if ∆m 2 31 ≡ m 2 3 − m 2 1 is positive (negative). Here and in the following we use the standard parameterization for the neutrino mass states and PMNS lepton mixing matrix [1]. Finding out which of these two possibilities is realized in Nature has profound implications for the flavor puzzle, as well as phenomenological consequences for cosmology, searches for neutrino mass, and for neutrinoless double-beta decay. Therefore, the determination of the mass ordering is one of the experimental priorities in the field. In particular, with the discovery of a large value of θ 13 [2][3][4][5] an answer within a decade or so is certainly possible and first hints may be obtained even sooner in global fits to the world's neutrino data.
New information is expected to come from long-baseline experiments, like T2K [6] and NOνA [7,8], which look for the appearance of ν e (ν e ) in a beam of ν µ (ν µ ). Proposals for a more long-term time frame include LBNE [9][10][11], LBNO [12], a superbeam based on the ESS [13], and eventually a neutrino factory [14]. Matter effects [15][16][17] will induce characteristic differences between the neutrino and antineutrino channels, which in turn will allow inference of the mass ordering, see e.g., Refs. [18,19] for early references. The fact that a comparison of neutrino and antineutrino channels is performed also implies that the leptonic CP phase δ cannot be ignored and has to be included in the analysis as well. A selective set of recent sensitivity studies for present and future proposed long baseline oscillation experiments can be found in Refs. [20][21][22][23][24][25][26][27][28][29][30][31].
Another possibility to determine the mass ordering arises from observing the energy and zenith angle dependence of atmospheric neutrinos in the GeV range, which will also have the mass ordering information imprinted by matter effects [32][33][34][35][36][37]. The flux of atmospheric neutrinos follows a steep power law with energy and thus the flux in the GeV range is quite small and requires very large detectors. IceCube technology can be adapted to neutrino energies in the GeV range by reducing the spacing of optical modules, eventually leading to the PINGU extension [38] and a similar low-energy modification can also be implemented for neutrino telescopes in the open ocean, called ORCA [39]. Another way to overcome the small neutrino flux is to separate neutrino and antineutrino events using a magnetic field like in the ICal@INO experiment [40,41] (INO for short in the following). Mass ordering sensitivity calculations have been performed for instance in Refs. [42][43][44][45][46][47][48][49][50] for PINGU/ORCA and in Refs. [51][52][53][54][55][56][57][58][59] for INO or similar setups.
Finally, the interference effects between the oscillations driven by ∆m 2 21 and ∆m 2 31 in the disappearance ofν e provide a third potential avenue for this measurement. In particular, this approach has been put forward in the context of reactor neutrinos [60]. JUNO [61,62] will comprise a 20 kt detector at a baseline of about 52 km of several nuclear reactors. A similar project is also discussed within the RENO collaboration [63]. The possibility to use a precision measurement of theν e survival probability at a nuclear reactor to identify the neutrino mass ordering has been considered by a number of authors, e.g., Refs. [49,62,[64][65][66][67][68][69][70][71][72][73][74].
This impressive experimental (and phenomenological) effort has also resulted in a renewed interest in potential issues arising from the statistical interpretation of the resulting data [75,76] (see also [71]), which can be summarized as: Given that the determination of the mass ordering is essentially a binary yes-or-no type question, are the usual techniques relying on a Taylor expansion around a single maximum of the likelihood applicable in this case? The goal of this paper is to answer this question within a frequentist framework for a wide range of experimental situations, including disappearance as well as appearance measurements. The answer we find in this paper can be stated succinctly as: The currently accepted methods yield approximately the expected frequentist coverage for the median experimental outcome; quantitative corrections typically lead to a (slightly) increased sensitivity compared to the standard approach prediction. The methods applied in the following are analogous to the ones from Ref. [77], where similar questions have been addressed for the discovery of θ 13 and CP violation. In the present work we strictly adhere to frequentist methods; Bayesian statistics is used to address the neutrino mass ordering question in Ref. [78], see also Refs. [75,76] for Bayesian considerations.
The outline of our paper is as follows. We first review the principles of hypothesis testing in a frequentist framework in Sec. 2, apply them to the case of the mass ordering, define the sensitivity of the median experiment and discuss the relation to the standard sensitivity based on ∆χ 2 values from the Asimov data set. In Sec. 3 we consider the Gaussian limit, where all relevant quantities, such as sensitivities can be expressed analytically. Details of the derivation can be found in App. A, as well as a discussion of conditions under which the Gaussian approximation is expected to hold. In Sec. 4 we present results from Monte Carlo simulations of the INO, PINGU, JUNO, NOνA, and LBNE experiments. The technical details regarding the simulations are summarized in App. B. We show that for most cases the Gaussian approximation is justified to good accuracy, with the largest deviations observed for NOνA. In Sec. 5 we present a comparison between the sensitivities expected for the different proposals, illustrating how the sensitivities may evolve with date. We summarize in Sec. 6, where we also provide a table which allows to translate the traditional performance indicator for the mass ordering (∆χ 2 without statistical fluctuations) into well defined frequentist sensitivity measures under the Gaussian approximation. We also comment briefly on how our results compare to those in Refs. [75,76].

Frequentist hypothesis testing
Let us start by reviewing the principles of frequentist hypothesis testing, see e.g., Ref. [1]. First we consider the case of so-called "simple hypotheses", where the hypothesis we want to test, H, as well as the alternative hypothesis H do not depend on any free parameters. H is conventionally called null hypothesis. In order to test whether data can reject the null hypothesis H we have to choose a test statistic T . A test statistic is a stochastic variable depending on the data which is chosen in such a way that the more extreme the outcome is considered to be, the larger (or smaller) the value of the test statistic is. Once the distribution of T is known under the assumption of H being true, we decide to reject H at confidence level (CL) 1 − α if the observation is within the α most extreme results, i.e., if T > T α c , where the critical value T α c is defined by with p(T |H) being the probability distribution function of T given that H is true. The probability α is the probability of making an "error of the first kind" (or type-I error rate), i.e., rejecting H although it is true. It is custom to convert this probability into a number of Gaussian standard deviations. In this work we will adopt the convention to use a doublesided Gaussian test for this conversion, such that a hypothesis is rejected if the data is more than nσ away (on either side) from the mean. This leads to the following conversion between nσ and the value of α: 1 α(n) = 2 √ 2π ∞ n dx e −x 2 /2 = erfc n √ 2 ⇔ n = √ 2 erfc −1 (α). (2.2) This definition implies that we identify, for instance, 1σ, 2σ, 3σ with a CL (1 − α) of 68.27%, 95.45%, 99.73%, respectively, which is a common convention in neutrino physics. However, note that nσ is sometimes defined differently, as a one-sided Gaussian limit, see e.g., Eq. (1) of Ref. [79]. This leads to a different conversion between nσ and α, namely which would lead to a CL of 84.14%, 97.73%, 99.87% for 1σ, 2σ, 3σ. In order to quantify how powerful a given test is for rejecting H at a given CL we have to compute the so-called "power" of the test or, equivalently, the probability of making an "error of the second kind" (or type-II error rate). This is the probability β to accept H if it is not true: where now p(T |H ) is the probability distribution function of T assuming that the alternative hypothesis H is true. Obviously, β depends on the CL (1 − α) at which we want to reject H. A small value of β means that the rate for an error of the second kind is small, i.e., the power of the test (which is defined as 1 − β) is large. The case we are interested in here (neutrino mass ordering) is slightly different, since both hypotheses (normal and inverted) may depend on additional parameters θ, a situation which is called "composite hypothesis testing". This is for instance the case of long baseline oscillation experiments, where the value of δ has a large impact on the sensitivities to the neutrino mass ordering. In this case the same approach is valid while keeping a few things in mind: • We can reject the hypothesis H only if we can reject all θ ∈ H. Thus, with This ensures that all θ ∈ H are rejected at confidence level (1 − α) if T > T α c . 2 • The rate of an error of the second kind will now depend on the true parameters in the alternative hypothesis: with T α c defined in Eq. (2.6). It is important to note that in a frequentist framework this cannot be averaged in any way to give some sort of mean rejection power, as this would require an assumption about the distribution of the parameters implemented in Nature (which is only possible in a Bayesian analysis [78]). Sticking to frequentist reasoning, we can either give β as a function of the parameters in the alternative hypothesis, or quote the highest and/or lowest possible values of β within the alternative hypothesis. 2 Here we assume that for given data the value of the observed test statistic T is independent of the true parameter values. This is the case for the statistic T introduced in Eq. (2.10), but it will not be true for instance for the statistic T mentioned in footnote 3.

Application to the neutrino mass ordering
In the search for the neutrino mass ordering, we are faced with two different mutually exclusive hypotheses, namely H NO for normal ordering and H IO for inverted ordering. Both hypotheses will depend on the values of the oscillation parameters (which we collectively denote by θ) within the corresponding ordering. In particular, appearance experiments depend crucially on the CP-violating phase δ. Hence, we have to deal with the situation of composite hypothesis testing as described above. Let us now select a specific test statistic for addressing this problem.
A common test statistic is the χ 2 with n degrees of freedom, which describes the deviation from the expected values of the outcome of a series of measurements x i of the normal distributions N (µ i , σ i ): The further the observations are from the expected values, i.e., the more extreme the outcome, the larger is the χ 2 . If the mean values µ i depend on a set of p parameters θ whose values have to be estimated from the data one usually considers the minimum of the χ 2 with respect to the parameters: According to Wilk's theorem [80] this quantity will follow a χ 2 distribution with n−p degrees of freedom, whereas ∆χ 2 (θ) = χ 2 (θ) − χ 2 min will have a χ 2 distribution with p degrees of freedom. The χ 2 distributions have known properties, and in physics we often encounter situations where data can be well described by this method and the conditions for Wilk's theorem to hold are sufficiently fulfilled, even when individual data points are not strictly normal distributed. In general, however, it is not guaranteed and the actual distribution of those test statistics has to be verified by Monte Carlo simulations [81].
Coming now to the problem of identifying the neutrino mass ordering, one needs to select a test statistic which is well suited to distinguish between the two hypotheses H NO and H IO . Here we will focus on the following test statistic, which is based on a log-likelihood ratio and has been used in the past in the literature: where θ is the set of neutrino oscillation parameters which are confined to a given mass ordering during the minimization. Let us stress that the choice of T is not unique. In principle one is free to chose any test statistic, although some will provide more powerful tests than others. 3 3 In the case of simple hypotheses the Neyman Pearson lemma [82] implies that the test based on the likelihood ratio is most powerful. For composite hypotheses in general no unique most powerful test is known. An alternative choice for a test statistic could be for instance the statistic T (θ) = χ 2 (θ) − χ 2 min , where χ 2 min is the absolute minimum including minimization over the two mass orderings, and θ generically denotes the (continuous) oscillation parameters. This statistic is based on parameter estimation and amounts to testing whether a parameter range for θ remains at a given CL in a given mass ordering. We It is important to note that within a frequentist approach, rejecting one hypothesis at a given α does not automatically imply that the other hypothesis could not also be rejected using the same data. Instead, the only statement we can make is to either reject an ordering or not. The value of T = 0 therefore does not a priori play a crucial role in the analysis. Let us illustrate this point at an example. In the left panel of Fig. 1, we show the distributions of the test statistics T for both mass orderings obtained from the simulation of a particular configuration of the JUNO reactor experiment. Experimental details will be discussed later in Sec. 4.1. In the right panel we show the corresponding critical values T α c for testing both orderings and how they depend on the chosen confidence level 1 − α. The curves for testing the different orderings cross around α = 5.2%, indicated by the dotted lines. This represents the unique confidence level for which the experiment in question will rule out exactly one of the orderings, regardless of the experimental outcome. If, for instance, we would choose to test whether either ordering can be rejected at a confidence level of 90%, then there is a possibility of an experimental outcome T with T 0.1 c,IO < T < T 0.1 c,NO , implying that both orderings could be rejected at the 90% CL. This situation is indicated by the have checked by explicit Monte Carlo simulations that typically the distribution of T is close to a χ 2 distribution with number of d.o.f. corresponding to the non-minimized parameters in the first term (the approximation is excellent for JUNO but somewhat worse for LBL experiments). Sensitivity results for the mass ordering based on T will be reported elsewhere.
dash-dotted line in the right panel of Fig. 1 and applies to the purple region. Thus, in order to claim a discovery of the mass ordering, it will not be sufficient to test one of the orderings. If both orderings were rejected at high confidence, it would mean either having obtained a very unlikely statistical fluctuation, underestimating the experimental errors, or neither ordering being a good description due to some new physics. Conversely, if we would choose α = 0.01 < 0.052 (dashed line in both panels, white region in right panel), then there is the possibility of obtaining T 0.01 c,NO < T < T 0.01 c,IO , meaning that neither ordering can be excluded at the 99% CL.
The CL corresponding to the crossing condition T α c,NO = T α c,IO provides a possible sensitivity measure of a given experiment. We will refer to it as "crossing sensitivity" below. 4 If T α c,NO ≈ −T α c,IO (as it is the case for the example shown in Fig. 1), this is equivalent to testing the sign of T . This test has been discussed also in Ref. [75,76]. From the definition of the sensitivity of an average experiment which we are going to give in the next subsection it will be clear that the crossing sensitivity is rather different from the median sensitivity, which is typically what is intended by "sensitivity" in the existing literature. It should also be noted that the critical values for the different orderings, as well as the crossing of the critical values, in general are not symmetric with respect to T = 0. The fact that Fig. 1 appears to be close to symmetric is a feature of the particular experiment as well of the test statistic T . This would not be the case for instance for the statistic T mentioned in footnote 3. Finally, note that Fig. 1 is only concerned with the critical value of T and its dependence on α. As such, it does not tell us anything about the probability of actually rejecting, for instance, inverted ordering if the normal ordering would be the true one (power of the test). As discussed above, this probability will typically also depend on the actual parameters within the alternative ordering and can therefore not be given a particular value. However, for the crossing point of the critical values, the rejection power for the other ordering is at least 1 − α.

Median sensitivity or the sensitivity of an average experiment
Let us elaborate on how to compare such a statistical analysis to previous sensitivity estimates massively employed in the literature, in particular in the context of long-baseline oscillation experiments. The most common performance indicator used for the mass ordering determination is given by for testing normal ordering, with an analogous definition for inverted ordering. This quantity corresponds to the test statistic T defined in Eq. (2.10) but the data x i are replaced by the predicted observables µ i (θ 0 ) at true parameter values θ 0 . Since no statistical fluctuations are 4 In the case of composite hypotheses, where the distribution of T depends on the true values of some parameters (e.g., the CP phase in the case of long-baseline experiments), we define T α c,NO and T α c,IO in analogy to Eq. (2.6), i.e., we chose the largest or smallest value of T α c (θ), depending on the mass ordering. Hence, the crossing sensitivity is independent of the true values of the parameters. included in this definition it is implicitly assumed that it is representative for an "average" experiment. (This is sometimes referred to as the Asimov data set [79], and T 0 is sometimes denoted as "∆χ 2 " [75].) T 0 is then evaluated assuming a χ 2 distribution with 1 dof in order to quote a "CL with which a given mass ordering can be identified". In the following, we will refer to this as the "standard method" or "standard sensitivity". Note that T 0 by itself is not a statistic, since it does not depend on any random data. The interpretation of assigning a χ 2 distribution to it is not well defined, and is motivated by the intuition based on nested hypothesis testing (which is not applicable for the mass ordering question). In the following we show that actually the relevant limiting distribution for T (but not for T 0 ) is Gaussian, not χ 2 .
The formalism described in section 2 allows a more precise definition of an "average" experiment. One possibility is to calculate the CL (1 − α) at which a false hypothesis can be rejected with a probability of 50%, i.e., with a rate for an error of the second kind of β = 0.5. In other words, the CL (1 − α) for β = 0.5 is the CL at which an experiment will reject the wrong mass ordering with a probability of 50%. We will call the probability α(β = 0.5) the "sensitivity of an average experiment" or "median sensitivity". This is the definition we are going to use in the following for comparing our simulations of the various experiments to the corresponding sensitivities derived from the standard method.
Let us note that the median sensitivity defined in this way is not the only relevant quantity in order to design an experiment, since in practice one would like to be more certain than 50% for being able to reject a wrong hypothesis. Under the Gaussian approximation to be discussed in the next section it is easy to calculate the sensitivity α for any desired β, once the median sensitivity is known.

THE GAUSSIAN CASE FOR THE TEST STATISTIC T
A crucial point in evaluating a statistical test is to know the distribution of the test statistic. In general this has to be estimated by explicit Monte Carlo simulations, an exercise which we are going to report on for a number of experiments later in this paper. However, under certain conditions the distribution of the statistic T defined in Eq. (2.10) can be derived analytically and corresponds to a normal distribution [75]: where N (µ, σ) denotes the normal distribution with mean µ and standard deviation σ and the + (−) sign holds for true NO (IO). 5 In general T NO 0 and T IO 0 may depend on model parameters θ. In that case the distribution of T will depend on the true parameter values and we have to consider the rules for composite hypothesis testing as outlined in section 2. We provide a derivation of Eq. (3.1) in App. A, where we also discuss the conditions that need to be fulfilled for this to hold in some detail. In addition to assumptions similar to the ones necessary for Wilk's theorem to hold, Eq. • we are dealing with simple hypotheses, or consider composite hypotheses at fixed parameter values, or • if close to the respective χ 2 minima the two hypotheses depend on the parameters "in the same way" (a precise definition is given via Eq. (A.21) in the appendix), or • if T 0 is large compared to the number of relevant parameters of the hypotheses.

Simple hypotheses
Let us now study the properties of the hypothesis test for the mass ordering based on the statistic T under the assumption that it indeed follows a normal distribution as in Eq. (3.1). First we consider simple hypotheses, i.e., T 0 does not depend on any free parameters. As we shall see below, this situation applies with good accuracy to the medium-baseline reactor experiment JUNO.
For definiteness we construct a test for H NO ; the one for H IO is obtained analogously. Since large values of the test statistic favor H NO over the alternative hypothesis H IO , we would reject H NO for too small values of T . Hence, we need to find a critical value T α c such that Let us now compute the power p of the test, i.e., the probability p with which we can reject H NO at the CL (1 − α) if the alternative hypothesis H IO is true. As mentioned above, p is related to the rate for an error of the second kind, β, since p = 1 − β. This probability is given by where the last approximation assumes T 0 ≡ T NO 0 ≈ T IO 0 , a situation we are going to encounter for instance in the case of JUNO below. We shown p = 1 − β as a function of T 0 for several values of α in the lower left panel of Fig. 2.
Equation (3.3) (or the lower left panel of Fig. 2) contains all the information needed to quantify the sensitivity of an experiment. In particular, it allows to address the question of how likely it is that the wrong mass ordering will be rejected at a given CL. For example, let us consider an experiment with a median sensitivity of 4σ, which implies T 0 ≈ 14.7. If we now demand that we want to reject the wrong mass ordering with a probability of 90% (β = 0.1), then this experiment will be able to do this only at slightly more than 99% CL. In the right panel of Fig. 2 we show β as a function of α for several fixed values of T 0 using Eq. (3.3). This plot allows a well defined interpretation of the "∆χ 2 " used in the standard method (i.e., T 0 ) under the Gaussian approximation. For a given T 0 and a chosen sensitivity α we can read off the probability with which the experiment will be able to reject the wrong ordering at the (1 − α) CL. Now it is also straight forward to compute the median sensitivity, which we have defined in section 2.3 as the α for which β = 0.5. From Eq. (3.3) we obtain Using our standard convention Eq. (2.2) to convert α into standard deviations the median sensitivity is nσ, with n = √ 2 erfc −1 1 2 erfc T 0 2 (median sensitivity). 3). The latter correspond to the "standard sensitivity" of n = √ T 0 and n = √ T 0 /2 for the crossing sensitivity. The edges of the green and yellow bands are obtained from the conditions on the rate for an error of the second kind β = 1/2 ± 0.6827/2 and β = 1/2 ± 0.9545/2, respectively. β = 1/2 ± 0.9545/2, respectively. They indicate the range of obtained rejection confidence levels which emerge from values of T within 1σ and 2σ from its mean assuming true IO.
Note that if we had used the 1-sided Gaussian rule from Eq. (2.3) to convert the probability Eq. (3.4) we would have obtained n = √ T 0 for the median sensitivity. Indeed, this corresponds exactly to the "standard sensitivity" as defined in section 2.3. 6 We show this case for illustration as dashed curve in Fig. 3. The dashed vertical lines in the right panel of Fig. 2 show explicitly that using this convention we obtain β = 0.5 at nσ exactly for T 0 = n 2 . Note that with our default convention we actually obtain an increase in the sensitivity compared to √ T 0 used in the "standard method". The exponential nature of erfc implies that the difference will not be large, in particular for large T 0 , see Fig. 3. For instance, the values of T 0 corresponding to a median sensitivity of 2σ, 3σ, 4σ according to Eq. (3.5) are 2.86, 7.74, 14.7, respectively, which should be compared to the standard case of T 0 = n 2 , i.e., 4,9,16. In summary, we obtain the first important result of this paper: the sensitivity obtained by using the standard method is very close to the median sensitivity within the Gaussian approximation. 6 We would have obtained the result n = √ T 0 also when using a 2-sided test to calculate α from the distribution of T combined with the 2-sided convention to convert it into standard deviations. Note, however, that for the purpose of rejecting a hypothesis clearly a 1-sided test for T should be used, and therefore we do not consider this possibility further.
Before concluding this section let us also mention the sensitivity defined by the crossing point T NO c = T IO c discussed at the end of section 2.2. This is the sensitivity α for which the critical values are the same for both orderings, which implies that regardless of the outcome of the experiment exactly one of the two hypotheses can be rejected at that CL. In the Gaussian approximation this implies that α = β, i.e., the rates for errors of the first and second kinds are the same. Using Eq. (3.2) and the analog expression for IO we obtain by The corresponding sensitivity is shown as red solid curve in Fig. 3. For this curve we use our default convention to convert α into σ according to Eq. (2.2). If we instead had used the 1sided Gaussian convention from Eq. (2.3) to convert the probability Eq. (3.6) we would have obtained the simple rule n = √ T 0 /2 (dashed red curve). This can be seen also in the right panel of Fig Fig. 3 one can see that for a "typical" experimental outcome the sensitivity will be significantly better than the one given by the crossing condition.

Composite hypotheses
Let us now generalize the discussion to the case where T 0 depends on parameters. This will be typically the situation for long-baseline experiments, where event rates depend significantly on the (unknown) value of the CP phase δ. It is straight forward to apply the rules discussed in section 2 assuming that T = N (T NO 0 (θ), 2 T NO 0 (θ)) for normal ordering and T = N (−T IO 0 (θ), 2 T IO 0 (θ)) for inverted ordering. First we must ensure that we can reject NO for all possible values of θ at (1−α) confidence. Hence, Eq. (3.2) becomes, i.e., we have to choose the smallest possible T α c . Considering T α c from Eq. (3.2) as a function of T 0 , we see that T α c has a minimum at T 0 = 2[erfc −1 (2α)] 2 , and the value at the minimum is −2[erfc −1 (2α)] 2 . This minimum is also visible in Fig. 2 (upper left panel). Hence, we have whereT NO 0 is the minimum of T NO 0 (θ) with respect to the parameters θ.
The expression for the rate for an error of the second kind, Eq. (3.3) will now depend on the true values of θ in the alternative hypothesis: The median sensitivity is obtained by setting β(θ) = 0.5. This leads to the equation T IO 0 (θ) = −(T α c ) min which has to be solved for α. Note that this is a recursive definition, since which case in Eq. (3.8) to be used can only be decided after α is computed. However, it turns out that in situations of interest the first case applies. In this case we have 2 for α corresponding to the median sensitivity. Hence, we obtain the result that is a useful expression for estimating the median sensitivity for composite hypotheses in the Gaussian approximation. We will confirm this later on by comparing it to the full Monte Carlo simulations of long-baseline experiments. Also note the similarity with the expression in case of simple hypotheses (see Eq. (3.4)).
Finally we can also calculate the "crossing sensitivity" by requiring (T α c ) NO min = (T α c ) IO min , for which exactly one hypothesis can be rejected. Again this is a recursive definition, however, ifT NO 0 T IO 0 it turns out that only the second case in Eq. (3.8) is relevant. This leads to where the last relation holds forT 0 ≡T NO 0 ≈T IO 0 , which again is very similar to the case for simple hypotheses, Eq. (3.6).

MONTE CARLO SIMULATIONS OF EXPERIMENTAL SETUPS
Let us now apply the methods presented above to realistic experimental configurations. We have performed Monte Carlo (MC) studies to determine the sensitivity to the neutrino mass ordering for three different types of experiments, each of which obtains their sensitivity through the observation of different phenomena: (a) JUNO [61]: interference (in the vacuum regime) between the solar and atmospheric oscillation amplitudes at a medium baseline reactor neutrino oscillation experiment; (b) PINGU [38] and INO [40]: matter effects in atmospheric neutrino oscillations; (c) NOνA [7] and LBNE [11]: matter effects in a long baseline neutrino beam experiment. In each case we have followed closely the information given in the respective proposals or design reports, and we adopted bench mark setups which under same assumptions reproduce standard sensitivities in the literature reasonably well. The specific details that have been used to simulate each experiment are summarized in App. B.
TABLE I: Sensitivity of the JUNO reactor experiment for 4320 kt GW yr exposure for two different assumptions on the energy resolution. We give the value of the test statistic without statistical fluctuation, T 0 , and the "standard sensitivity" √ T 0 σ. The median sensitivity is calculated according to Eq. (3.4). The "crossing sensitivity" corresponds to the CL where exactly one mass ordering can be rejected regardless of the outcome, which is calculated according to Eq. (3.6).

Medium-baseline reactor experiment: JUNO
For the simulations in this paper we adopt an experimental configuration for the JUNO reactor experiment based on Refs. [61,62,83], following the analysis described in Ref. [49]. A 20 kt liquid scintillator detector is considered at a distance of approximately 52 km from 10 reactors with a total power of 36 GW, with an exposure of 6 years, i.e., 4320 kt GW yr. The energy resolution is assumed to be 3% 1 MeV/E. For further details see App. 1.
The unique feature of this setup is that the sensitivity to the mass ordering is rather insensitive to the true values of the oscillation parameters within their uncertainties. Being aν e disappearance experiment, the survival probability depends neither on θ 23 nor on the CP phase δ, and all the other oscillation parameters are known (or will be known at the time of the data analysis of the experiment) with sufficient precision such that the mass ordering sensitivity is barely affected. Therefore we are effectively very close to the situation of simple hypotheses for this setup. Note that although the mass ordering sensitivity is insensitive to the true values, the χ 2 minimization with respect to oscillation parameters, especially |∆m 2 31 |, is crucial when calculating the value of the test statistic T . In the left panel of Fig. 1 we show the distribution of the test statistic T from a Monte Carlo simulation of 10 5 data sets for our default JUNO configuration. For each true mass ordering we compare those results to the normal distributions expected under the Gaussian approximation, namely N ( are the values of the test statistic without statistical fluctuation (Asimov data set). For the considered setup we find T NO 0 = 10.1 and T IO 0 = 11.1, and we observe excellent agreement of the Gaussian approximation with the Monte Carlo simulation, see also, e.g., Ref. [70].
Hence we can apply the formalism developed in section 3 directly to evaluate the sensitivity of the experiment in terms of T NO 0 and T IO 0 . For instance, Eq. (3.4) gives for the median sensitivity α = 7.3 (4.3) × 10 −4 for testing normal (inverted) ordering, which corresponds to 3.4σ (3.5σ). As discussed in section 3 those numbers are rather close to the "standard sensitivity" based on n = √ T 0 , which would give 3.2σ (3.3σ). For the given values of T NO 0 and T IO 0 we can now use Fig. 2 to obtain the probability to reject an ordering if it is false (i.e., the power of the test) for any desired confidence level (1 − α). The confidence level at which exactly one mass ordering can be rejected (crossing point T NO  from Eq. (3.6) as α = 5.2% or 1.9σ, see also Fig. 1. Those numbers are summarized in Tab. I. There we give also the corresponding results for the same setup but with a slightly worse energy resolution of 3.5% 1 MeV/E, in which case significantly reduced sensitivities are obtained, highlighting once more the importance to achieve excellent energy reconstruction abilities. We have checked that also in this case the distribution of T is very close to the Gaussian approximation.

Atmospheric neutrinos: PINGU and INO
We now move to atmospheric neutrino experiments, which try to determine the mass ordering by looking for the imprint of the matter resonance in the angular and energy distribution of neutrino induced muons. The resonance will occur for neutrinos (antineutrinos) in the case of normal (inverted) ordering. The INO experiment [40] uses a magnetized iron calorimeter which is able to separate neutrino and antineutrino induced events with high efficiency, which provides sensitivity to the mass ordering with an exposure of around 500 kt yr (10 year operation of a 50 kt detector). Alternatively, the PINGU [38] experiment, being a low-energy extension of the IceCube detector, is not able to separate neutrino and antineutrino induced muons on an event-by-event basis. This leads to a dilution of the effect of changing the mass ordering, which has to be compensated by exposures exceeding 10 Mt yr, which can be achieved for a few years of running time. In both cases the ability to reconstruct neutrino energy and direction will be crucial to determining the mass ordering.
Our simulations for the INO and PINGU experiments are based on Refs. [58] and [49], respectively. We summarize the main characteristics of our default setups in Tab. II, further technical details and references are given in App. 2. Let us stress that the sensitivity of this type of experiments crucially depends on experimental parameters such a systematic uncertainties, efficiencies, particle identification, and especially the ability to reconstruct neutrino energy and direction. Those parameters are still not settled, in particular for the PINGU experiment, and final sensitivities may vary by few sigmas, see for instance Refs. [38,48]. Our setups should serve as representative examples in order to study the statistical properties of the resulting sensitivities. While the final numerical answer will depend strongly on to be defined experimental parameters, we do not expect that the statistical behavior will be affected significantly.
In Figs. 4 and 5 we show the distributions of the test statistic T for the INO and PINGU experiments, respectively, obtained from a sample of 10 4 simulated data sets for each mass ordering, using the default setups from Tab. II. We observe good agreement with the Gaus- sian approximation (see also Ref. [46] for a simulation in the context of PINGU). Those results justify the use of the simple expressions from section 3 also for INO and PINGU in order to calculate median sensitivities or rates for errors of the first and second kind. In Fig. 5 we illustrate the dependence of the distributions for PINGU on the true value of θ 23 . From this figure it is clear that the true value of θ 23 plays an important role for the sensitivity to the mass ordering, with better sensitivity for large values of θ 23 (a similar dependence is holds also for INO, see, e.g., Refs. [54,58]). The dependence on other parameters is rather weak (taking into account that, at the time of the experiment, θ 13 will be known even better than today). Let us discuss the θ 23 dependence in more detail for the case of PINGU, where from now on we use the Gaussian approximation. The problem arises when calculating the critical value for the test statistic T in order to reject the null-hypothesis at a given CL. If we follow our rule for composite hypothesis, Eq. (2.6), and minimize (for NO) or maximize (for IO) T α c (θ 23 ) over θ 23 in the range 35 • to 55 • we obtain the black dashed curves in Fig. 6. This is equivalent to using Eq. (3.10). The chosen range for θ 23 corresponds roughly to the 3σ range obtained from current data [84]. However, this may be too conservative, since at the time of the experiment T2K and NOνA will provide a very accurate determination of sin 2 2θ 23 . Hence, θ 23 will be known with good precision up to its octant, see for instance Fig. 5 of Ref. [85]. If we minimize (maximize) T α c (θ 23 ) only over the two discrete values θ true 23 and 90 • − θ true 23 we obtain the blue solid curves in Fig. 6. The green and yellow bands indicate the corresponding 68.27% and 95.45% probability ranges of expected rejection significances. The dotted curves show the corresponding information but using only the true value of θ 23 when calculating T α c . This last case corresponds to the ideal situation of perfectly knowing θ 23 (including its octant), in which case NO and IO become simple hypotheses. The median sensitivity for known θ 23 is not shown in the figure for clarity, but it is very similar to the blue solid curves.
We obtain the pleasant result that all three methods give very similar values for the median sensitivity, ranging from 2σ at θ 23 35 • up to 5σ (6σ) rejection of NO (IO) at θ 23 55 • . Only for the NO test and θ 23 50 • we find that taking the octant degeneracy into account leads to a larger spread of the 68.27% and 95.45% probability ranges for the sensitivity, implying a higher risk of obtaining a rather weak rejection. Actually, this region of parameter space (true IO and θ 23 > 45 • ) is the only one where the octant degeneracy severely affects the sensitivity to the mass ordering [48]. Let us emphasize that the octant degeneracy is always fully taken into account when minimizing the χ 2 . Here we are instead concerned with the dependence of the critical value T α c on θ 23 .

Long-baseline appearance experiments: NOνA and LBNE
Long-baseline neutrino beam experiments try to identify the neutrino mass ordering by exploring the matter effect in the ν µ → ν e appearance channel. Whether the resonance occurs for neutrinos or for antineutrinos will determine the mass ordering. A crucial feature in this case is that the appearance probability, and therefore also the event rates, depend significantly on the unknown value of the CP phase δ. Most likely δ will remain unknown even at the time the mass ordering measurement will be performed, and therefore taking the δ dependence into account is essential. In the nomenclature of sections 2 and 3 we are dealing with composite hypothesis testing. In this work we consider three representative experimental configurations to study the statistical properties of the mass ordering sensitivity, namely NOνA [7], LBNE-10 kt, and LBNE-34 kt [11], which provide increasing sensitivity to the mass ordering. Tab III summarizes their main features, while further details are given in App. 3.
Figs. 7 and 8 show the probability distributions for the test statistic T defined in Eq. (2.10), for the NOνA and LBNE-10 kt setups, respectively. The distributions are shown for both mass orderings, and for different values of δ, as indicated in each panel. Our results are based on a sample of 6 × 10 5 simulations for NOνA and 4 × 10 5 for LBNE-10 kt per value of δ, and we scan δ in steps of 10 • . As can be seen from the figures, both the shape and mean of the distributions present large variations with the value of δ. From the comparison between the two figures it is clear that the NOνA experiment will achieve very limited sensitivity to the mass ordering. On the other hand, for the LBNE-10 kt setup the situation is much better: the overlapping region is reduced, and is only sizable for certain combinations of values of δ in the two mass orderings.
We also note that for NOνA there are clear deviations from the Gaussian shape for the T distributions, while for the LBNE-10 kt experiment they are close to the Gaussian approximation discussed in section 3, namely T = N (±T 0 (θ), 2 T 0 (θ)). For comparison, in Fig. 8 the Gaussian approximation is overlaid on the histograms from the Monte Carlo. Those results are in agreement with the considerations of App. A. As discussed there, one expects that the median of the T distribution should remain around ±T 0 , even if corrections to the shape of the distribution are significant. We have checked that this does indeed hold for NOνA. Furthermore, assuming that there is only one relevant parameter (δ in this case), Eq. (A.24) implies that deviations from Gaussianity can be expected if T 0 ∼ 1, which is the case for NOνA, whereas for T 0 1 (such as for LBNE) one expects close to Gaussian distributions for T .
One can also notice in Figs. 7 and 8 that the shape of the distributions for a given value of δ in one ordering is rather similar to the mirrored image of the distribution corresponding to the other mass ordering and −δ. The reason for this is the well-known fact that the standard mass ordering sensitivity is symmetric between changing the true ordering and δ → −δ, i.e., T NO 0 (δ) ≈ T IO 0 (−δ), see e.g., Figs. 8 and 9 of Ref. [8] and Fig. 4-13 of Ref. [11]. 7 Furthermore, using the formalism in App. A, in particular Eq. (A.24), one can show that also the deviations from the Gaussian distribution will obey the same symmetry. Below we will show that despite the deviations from Gaussianity for NOνA, the final sensitivities obtained from the Monte Carlo will be surprisingly close to the Gaussian expectation. As expected, this will be even more true for LBNE-10 kt.
Due to the strong dependence on the CP phase δ we need to choose the critical value T α c such that the null hypothesis can be rejected at (1 − α) CL for all possible values of δ, see discussion in sections 2 and 3.2. This is illustrated in Fig. 9, which is analogous to Fig. 1 (right panel) for a fixed CL. The continuous (dashed) black curves in Fig. 9 show the values of T α c that lead to the probability of 5% to find a smaller (larger) value of T under the hypothesis of a true normal (inverted) ordering as a function of the true value of δ. The left panel shows the result for NOνA, while the right panel corresponds to LBNE-34 kt. The number of data sets simulated for LBNE-34 kt in this case is 10 5 per value of δ, which is again scanned in steps of 10 • . As discussed in Sec. 2, a composite null hypothesis can only be rejected if we can reject all parameter sets θ ∈ H. In our case, this would imply rejecting the hypothesis for all values of δ. Therefore, in order to guarantee a CL equal to (1 − α), the most conservative value of T α c will have to be chosen. This automatically defines two values T α c (NO) and T α c (IO), which are the values which guarantee that a given hypothesis can be rejected at the 95% CL. These values will generally be different, and are indicated in the figures by the arrows. In Fig. 9 we encounter the two situations already discussed in Sec. 2 (cf. Fig. 1):    Fig. 1, which defines the CL at which exactly one of the hypotheses can be excluded.
Let us now evaluate the rate for an error of the second kind corresponding to a given value of α. After the value of T α c is determined for a given hypothesis and α, we can compute the rate for an error of the second kind, β, as a function of the true value of δ, as discussed in Sec. 2. We show this probability in Fig. 10 for the NOνA and the LBNE-10 kt experiments in the left-and right-hand panels, respectively. To be explicit, we show the probability of accepting normal ordering at 1σ, 2σ, 3σ CL, i.e., α = 32%, 4.55%, 0.27%, (regardless of the value of δ in the NO) although the true ordering is inverted. This probability depends on the true value of δ in the IO, which is shown on the horizontal axis. By doing a cut at β = 0.5 on the left-hand panel (indicated by the dotted line), we can get an idea on the median sensitivity that will be obtained for NOνA: for δ = −90 • it will be around 1σ, while for δ = 90 • it will reach almost the 3σ level. This seems to be roughly consistent with the expected standard sensitivities usually reported in the literature, see for instance Ref. [8]. Similarly, for LBNE-10 kt, we expect that the sensitivity for the median experiment will be around 3σ for δ = −90 • , while for other values of δ we expect it to be much larger. This is also in agreement with the results from Ref. [11], for instance.
Let us now investigate in detail how our median sensitivity compares to the "standard sensitivities" widely used in the literature. In Fig. 11 the solid thick curves show the results for the median sensitivity derived from full MC simulations. The shaded green and yellow bands are analogous to those shown in Fig. 3, and show the range in the number of sigmas with which we expect to be able to reject NO if IO is true in 68.27% and 95.45% of the experiments, respectively. We also show how these results compare to the Gaussian approximation discussed in section 3. The value of the χ 2 is computed without taking statistical fluctuations into account (what is called T 0 in Sec. 2). We then use Eq. (3.10) to compute the confidence level (1 − α) at which the normal ordering can be rejected with a probability of 50% if the inverted ordering is true, as a function of the true value of δ in the IO. Then, for the dot-dashed curves we use a 2-sided Gaussian to convert α into number of σ, i.e., Eq. (2.2), the same prescription is also used for the MC result. We observe good agreement, in particular for LBNE. This indicates that, for the high-statistics data from LBNE, we are very close to the Gaussian limit, whereas from the smaller data sample (and smaller values of T 0 ) in NOνA deviations are visible, but not dramatic. We also show the results using a 1-sided Gaussian, Eq. (2.3), to convert α into number of sigmas, which leads to n = √ T 0 , i.e., the standard sensitivity. This is shown by the dashed lines. As discussed in Sec. 2 we observe that the standard sensitivity slightly under-estimates the true sensitivity. 8 Finally, the dotted horizontal line in Fig. 11 corresponds to the significance of the crossing point T NO c = T IO c defined in Sec. 2, i.e., the confidence level at which exactly one hypothesis can be excluded regardless of the outcome of the experiment. The results are independent of the value of δ, and guarantee that the rate for an error of the second kind β is at most equal to α, unlike for the median experiment where β = 0.5. The results for the crossing point are also consistent with the Gaussian expectation Eq. (3.11).

COMPARISON BETWEEN FACILITIES: FUTURE PROSPECTS
In this section we give a quantitative comparison between the different experiments that have been considered in this paper. We do a careful simulation of all the facilities using the details available in the literature from the different collaborations, see App. B for details. We have checked that our standard sensitivities are in good agreement with the respective proposals or design reports. Nevertheless, we do not explore in which way the assumptions made in the literature towards efficiencies, energy resolution, angular resolution, systematics, etc may affect the results, with the only exception of JUNO, as we explain below. Since we are mainly interested in the statistical method for determining the mass ordering, such analysis is beyond the scope of this paper. Our results will be shown as a function of the date, taking the starting points from the official statements of each collaboration. Obviously, such projections always are subject to large uncertainties. Fig. 12 shows the median sensitivities for the various experiments, i.e., the number of sigmas with which an "average experiment" for each facility can rejected a given mass ordering if it is false. In some sense this is similar to the standard sensitivity of √ T 0 commonly applied in the literature. A different question is answered in Fig. 13, namely: what is the probability that the wrong mass ordering can be rejected at a confidence level of 3σ? The confidence level has been chosen arbitrarily to 3σ, based on the convention that this would correspond to "evidence" that the wrong ordering is false. Below we discuss those 8 Note that traditionally the "standard sensitivity for IO" denotes the case when IO is true and refers to the sensitivity to reject NO. In the language of the present paper we call this a "test for NO". This is also consistent with the formula in the Gaussian approximation, Eq. (3.10), which contains T IO 0 when considering a test for NO. This has to be taken into account when comparing e.g., Fig. 11 (corresponding to a test for NO) to similar curves in the literature. plots in some detail.
In order to keep the number of MC simulations down to a feasible level, we use the Gaussian approximation whenever it is reasonably justified. As we have shown in Sec. 4, this is indeed the case for PINGU, INO, and JUNO. With respect to the LBL experiments, even though we have seen that the agreement with the Gaussian case is actually quite good (see Fig. 11), there are still some deviations, in particular in the case of NOνA. Consequently, in this case we have decided to use the results from the full MC simulation whenever possible. The results for the NOνA experiment are always obtained using MC simulations, while in the case of LBNE-10 kt the results from a full MC are used whenever the number of simulations does not have to exceed 4 × 10 5 (per value of δ). As was mentioned in the caption of Fig. 11, this means that, in order to reach sensitivities above ∼ 4σ (for the median experiment), results from the full MC cannot be used. In these cases, we will compute our results using the Gaussian approximation instead. As mentioned in App. A, the approximation is expected to be quite accurate precisely for large values of T 0 . Finally, for LBNE-34 kt, all the results have to be computed using the Gaussian approximation, since the median sensitivity for this experiment reaches the 4σ bound already for one year of exposure only, even for the most unfavorable values of δ.
For each experiment, we have determined the parameter that has the largest impact on the results, and we draw a band according to it to show the range of sensitivities that should be expected in each case. Therefore, we want to stress that the meaning of each band may be different, depending on the particular experiment that is considered. In the case of long baseline experiments (NOνA, LBNE-10 kt and LBNE-34 kt), the results mainly depend on the value of the CP-violating phase δ. In this case, we do a composite hypothesis test as described in Secs. 2 and 3.2, and we draw the edges of the band using the values of true δ in the true ordering that give the worst and the best results for each setup. Nevertheless, since for these experiments the impact due to the true value of θ 23 is also relevant, we show two results, corresponding to values of θ 23 in the first and second octant. In all cases, the octant degeneracy is fully searched for (see App. 3 for details). In the case of PINGU and INO, the most relevant parameter is θ 23 . We find that, depending on the combination of true ordering and θ 23 the results will be very different. Therefore, in this case we also do a composite hypothesis test, using θ 23 as an extra parameter. Finally, the case of JUNO is somewhat different. In this case, the uncertainties on the oscillation parameters do not have a big impact on the results. Instead, the energy resolution is the parameter which is expected to have the greatest impact, see for instance Ref. [73] for a detailed discussion. Therefore, in this case the width of the band shows the change on the results when the energy resolution is changed between 3% 1 MeV/E and 3.5% 1 MeV/E. For JUNO we do a simple hypothesis test, as described in Sec. 3.1.
The starting dates assumed for each experiment are: 2017 for INO [86], 2019 for PINGU [38] and JUNO [61] and 2022 for LBNE [87]. Note that the official running times for PINGU and JUNO are 5 and 6 years, respectively. For illustrative purposes we extend the time in the plots to 10 years, in order to see how sensitivities would evolve under the adopted assumptions about systematics. For the NOνA experiment, we assume that the nominal luminosity will be achieved by 2014 [8] and we consider 6 years of data taking from that moment on.
From the comparison of Figs. 12 and 13 one can see that, even though the median sensitivity for INO would stay below the 3σ CL, there may be a sizable probability (up to ∼ 40%) that a statistical fluctuation will bring the result up to 3σ. For NOνA, such probability could even go up to a 60%, depending on the combination of θ 23 , δ and the true mass ordering. In the case of LBNE, the dependence on the true value of δ is remarkable, in particular for the power of the test. We clearly observe the superior performance of the 34 kt configuration over the 10 kt one. For 34 kt a 3σ result can be obtained at very high probability for all values of δ, and for some values of δ a much higher rejection significance of the wrong ordering is achieved with high probability.
For the atmospheric neutrino experiments INO and PINGU we show the effect of changing the true value of θ 23 from 40 • to 50 • . The effect is particularly large for PINGU and a true NO. As visible in Fig. 6, for NO the sensitivity changes significantly between 40 • and 50 • , whereas for IO they happen to be similar, as reflected by the width of the bands in Figs. 12 and 13. The reason for this behavior is that for true IO and θ 23 > 45 • the mass ordering sensitivity is reduced due to the octant degeneracy [48]. In the context of PINGU, let us stress that the precise experimental properties (in particular the ability to reconstruct neutrino energy and direction) are still very much under investigation [38]. While we consider our adopted configuration (see Sec. 4.2 and App. 2 for details) as a representative bench mark scenario, the real sensitivity may be easily different by few standard deviations, once the actual reconstruction abilities and other experimental parameters are identified. To lesser extent this applies also to INO.
Let us also mention that in this work we only consider the sensitivity of individual experiments, and did not combine different setups. It has been pointed out in a number of studies that the sensitivity can be significantly boosted in this way [48,49,58,59,85]. We also expect that in this case, if the combined T 0 is sufficiently large, the Gaussian approximation should hold. However, we stress that a detailed investigation of this question is certainly worth pursuing in future work.

DISCUSSION AND SUMMARY
The sensitivity of a statistical test is quantified by reporting two numbers: 1. the confidence level (1 − α) at which we want to reject a given hypothesis, which corresponds to a rate for an error of the first kind, α; and 2. the probability p with which a hypothesis can be rejected at CL (1 − α) if it is false (the power of the test), which is related to the rate for an error of the second kind, β = 1 − p .
In this work we have applied this standard approach to the determination of the type of the neutrino mass ordering. With the help of those concepts it is straight forward to quantify the sensitivity of a given experimental configuration aiming to answer this important question in neutrino physics. We consider a test statistic T (see Eq. (2.10)) in order to perform the test, which is based on the ratio of the likelihood maxima under the two hypotheses normal and inverted ordering. Under certain conditions, see App. A, the statistic T is normal distributed (Gaussian approximation) [75]. In the limit of no statistical fluctuations (Asimov data set) the test statistic T becomes the usual ∆χ 2 (up to a sign) massively used in the literature for sensitivity calculations. In this work we denote this quantity by T 0 (in Ref. [75] it has been denoted by ∆χ 2 ). The sensitivity of an average experiment (in the frequentist sense) can be defined as the confidence level (1 − α) at which a given hypothesis can be rejected  = T IO 0 . The columns show T 0 , the standard sensitivity n = √ T 0 , the median sensitivity (Eqs. (3.4), (3.5)), the crossing sensitivity where exactly one hypothesis is rejected (equivalent to testing the sign of T , Eq. (3.6)), the probability β of accepting a mass ordering at the 3σ CL although it is false (rate for an error of the second kind, Eq. (3.3)), and the range of rejection confidence levels obtained with a probability of 68.27% and 95.45%. We convert CL into standard deviations using a 2-sided Gaussian.
with a probability β = 50% ("median sensitivity"). An important result of our work is the following: The sensitivity obtained by using the standard method of taking the square-root of the ∆χ 2 without statistical fluctuations is very close to the median sensitivity obtained within the Gaussian approximation for the test statistic T .
In section 3 we provide simple formulas, based on the Gaussian approximation, which allow quantification of the sensitivity in terms of error rates of the first and second kind for a given T 0 . For instance, Eqs. (3.3) and (3.9) contain simple expressions for the computation of β for given values of α and T 0 , whereas Eq. (3.10) allows the computation of the median sensitivity in terms of T 0 . In Tab. IV we give a collection of sensitivity measures based on the Gaussian approximation for the three example values T 0 = 9, 16, 25. The columns "std. sens." and "median sens." demonstrate explicitly the statement emphasized above, that the median sensitivity is close to the n = √ T 0 rule. The crossing sensitivity corresponds to the CL at which exactly one of the two hypotheses can be rejected. This is similar to testing the sign of the test statistic T , a test which has been discussed in Ref. [76] and also mentioned in Ref. [75]. By construction, this test gives smaller confidence levels than the median sensitivity and is not necessarily connected to what would be expected from an experiment. We give in the table also the probability for accepting a hypothesis at the 3σ level although it is false (rate for an error of the second kind). The last two columns in the table give the range of obtained rejection significance with a probability of 68.27% and 95.45% (assuming that the experiment would be repeated many times). Those are a few examples of how to apply the equations from section 3. These sensitivity measures provide different information and all serve to quantify the sensitivity of an experiment within a frequentist framework. They can be compared to similar sensitivity measures given in Ref. [75] in a Bayesian context (see, e.g., their Tab. IV).
In the second part of the paper we report on the results from Monte Carlo simulations for several experimental setups which aim to address the neutrino mass ordering: the mediumbaseline reactor experiment JUNO, the atmospheric neutrino experiments INO and PINGU, and the long-baseline beam experiments NOνA and LBNE. In each case we have checked by generating a large number of random data sets how well the Gaussian approximation is satisfied. Our results indicate that the Gaussian approximation is excellent for JUNO, INO, and PINGU. For NOνA the T distributions deviate significantly from Gaussian (strongly dependent on the true value of the CP phase δ), however the Gaussian expressions for the sensitivities still provide a fair approximation to the results of the Monte Carlo. For LBNE the Gaussian approximation is again fulfilled reasonably well. This is in agreement with our analytical considerations on the validity of the Gaussian approximation given in App. A, where we find that for experiments with T 0 large compared to the number of relevant parameters Gaussiantiy should hold. Hence, we expect that the Gaussian approximation should hold to very good accuracy also for experiments with a high sensitivity to the mass ordering, such as for instance a neutrino factory [14,88,89] or the LBNO experiment [12], when explicit Monte Carlo simulations become exceedingly unpractical due to the very large number of data sets needed in order to explore the high confidence levels.
In section 5 we provide a comparison of the sensitivities of the above mentioned facilities using the statistical methods discussed in this paper. Figures 12 and 13 illustrate how the median sensitivity and the probability to reject the wrong mass ordering at 3σ CL for the various experiments, respectively, could evolve as function of time based on official statements of the collaborations. While this type of plots is subject to large error bars on the time axis (typically asymmetric) as well as concerning actual experimental parameters, our results indicate that it is likely that the wrong mass ordering will be excluded at 3σ CL within the next 10 to 15 years. and similar for H . Hereθ α are the parameters at the minimum, which will be different for each hypothesis. In practice often the variances have to be estimated from the data itself, e.g., σ i ≈ σ i ≈ √ x i . In the following we will assume σ i = σ i . Let us note that generalization to correlated data is straight forward. The test statistic T from Eq. (2.10) is then given by T = X 2 (H ) − X 2 (H). In the following we will derive the distribution of T .
The distributions of X(H) and X(H ). Let us assume for definiteness that H is true. First we consider X 2 (H), and we derive the well-known result, that X 2 (H) is distributed as Under H, the y i (θ 0 α ) = n i are N standard normal distributed variables with N (0, 1). Then we have X 2 (H) = min θ i [y i (θ α )] 2 . The minimum condition is Asymptotically the parameter values at the minimumθ α will converge to the true values θ 0 α . Therefore we assume and expand Here and in the following sums run over α, β = 1, . . . , P and i, j, k = 1, . . . , N if not explicitly noted otherwise. Then the minimum condition Eq. (A.4) becomes and we obtain . Then Eq. (A.7) can be written as and The matrix (V iα ) defined in Eq. (A.9) is a rectangular N × P matrix which per construction obeys the orthogonality condition i V iα V iβ = δ αβ . Hence, we can always complete it by N − P columns to a full orthogonal N × N matrix such that k V ik V jk = δ ij and k V ki V kj = δ ij . Then we have where w r ≡ i V ri n i are N − P independent variables distributed as N (0, 1). This shows explicitly that if H is true, X 2 (H) is distributed as a χ 2 with N − P d.o.f. [80].
Let us now derive the distribution of X 2 (H ) under the assumption that H is true. Again we define however, now y i will not be standard normal distributed as N (0, 1), since per assumption x i have mean µ i (θ 0 α ) (and not µ i ). Nevertheless we can assume that y i (θ α ) can be expanded around a fixed reference point θ * α , such that the minimum in the wrong hypothesis,θ α , converges asymptotically towards it. We write and Here n i are N (0, 1) as before, but n i are N (m i , 1). Now the calculation proceeds as before and we arrive at where w r ≡ i V ri n i are now N − P independent normal variables with mean w r = i V ri m i . Then X 2 (H ) has a so-called non-central χ 2 distribution with N − P d.o.f. and a non-centrality parameter ∆ = N r=P +1 w r 2 .
The distribution of the test statistic T . Let us now consider the test statistic T = X 2 (H ) − X 2 (H). Using Eqs. (A.11) and (A.15) we find: The first term in Eq. (A.17) is just a constant, independent of the data. Using the definition of m i in Eq. (A.14) and comparing with Eq. (A.11) one can see that this term is identical to X 2 (H ) but replacing the data x i by the prediction for H at the true values: This is nothing else than the usual "∆χ 2 " between the two hypotheses without statistical fluctuations, compare Eq. (2.11). The second term in Eq. (A.17) is a sum of N standard normal variables, i a i n i . This gives a normal variable with variance i a 2 i . It is easy to show from Eq. (A.17) that the variance is 4T 0 and Eq. (A.17) can be written as with n standard normal. Hence, we find that if the term in Eq. (A.18) can be neglected, T is gaussian distributed with mean T 0 and standard deviation 2 √ T 0 [75]. Consider now the term in Eq. (A.18). Using n i n j = δ ij and the orthonormality of V and V we obtain that the mean value of this term is P − P . Hence, if the number of parameters in the two hypotheses is equal (as it is the case for the neutrino mass ordering) the mean value of T remains T 0 as in the Gaussian approximation, Eq. (A.20). For testing hypotheses with different numbers of parameters the mean value will be shifted from T 0 . However, even if the mean value remains unaffected, the higher moments of the distribution can still be modified. Under which conditions can the term in Eq. (A.18) be neglected?
• Obviously this term is absent if no parameters are estimated from the data, P = P = 0, i.e., for simple hypotheses. This applies in particular, if we compare the two hypotheses for fixed parameters.
• The term in Eq. (A.18) will also vanish for This condition has a geometrical interpretation. Consider the N dimensional space of data. Varying P parameters θ α the predictions µ i (θ α ) describe a P dimensional subspace in the N dimensional space. The operator V V T is a projection operator into the tangential hyperplane to this subspace at the X 2 minimum. This can be seen by considering the definition of V in Eq. (A.9) and of B in Eq. (A.5), which show that V is determined by the derivatives ∂µ i /∂θ α at the minimum. Similar, V V T projects into the P dimensional tangential hyperplane at the minimum corresponding to H . Hence, the condition (A.21) means that the hyperplanes of the two hypotheses have to be parallel at the minima. Obviously this condition can be satisfied only if the dimensions of the hyperplanes are the same, i.e., P = P .
We note also that the condition (A.21) is invariant under a change of parameterization, which amounts to B → BS with S αβ ≡ ∂θ α /∂θ β being a P × P orthogonal matrix describing the variable transformation θ α →θ α . Such a transformation would just change the orthogonal matrix R, but leave the operator V V T invariant. Roughly speaking we can say that sufficiently close to the respective minima θ 0 α and θ * α , the two hypotheses should depend on the parameters "in the same way", where the precise meaning is given by Eq. (A.21).
• Irrespective of the above conditions, we can neglect Eq. (A.18) if its variance is much smaller than the variance of the term in Eq. (A.17), which is given by 4T 0 . Eq. (A.18) is the difference of two χ 2 distributions with P and P d.o.f., respectively. The χ 2 n distribution has a mean and variance of n and 2n, respectively. Hence, we should be able to neglect this term if T 0 P, P , i.e., for high sensitivity experiments.
Example with one parameter. To simplify the situation let us consider the case where just one parameter θ is estimated from the data, for both H and H . The matrix V iα defined in Eq. (A.9) becomes now just a normalized column vector and similar for V i . The term in Eq. (A.18) is now just the difference of the square of two standard normal variables: n 2 − n 2 , with n = i V i n i and n = i V i n i . As mentioned above for the general case, we see that n 2 − n 2 = 0. The variance of this term is obtained as where we have used that n 4 i = 3. We can write (V T V ) 2 = Tr[V V T V V T ] = cos 2 ϕ, where ϕ is the angle between the two hyperplanes (i.e., lines, in this case) for H and H . Hence we find that the variance is zero if |V T V | = 1, i.e., the lines are parallel. And we have a measure to estimate when Eq. (A.20) is valid, namely when the variance of Eq. (A.18) is small compared to the variance of the second term in Eq. (A.17). In the example of one parameter this means Since sin 2 ϕ ≤ 1 we find that for T 0 1 the gaussian approximation is expected to be valid if only one parameter is estimated from data.

Appendix B: Simulation details
In the following, we describe the main details that have been used to simulate the experimental setups considered in this work. Unless stated otherwise the true values for the oscillation parameters have been set to the following values [84], and the χ 2 (or the test statistic T ) has been minimized with respect to them by adding Gaussian penalty terms to the χ 2 with the following 1σ errors: Unless otherwise stated, we assume the true value of θ 23 to be in the first octant. Nevertheless, the region around π/2 − θ 23 would not be disfavored by the penalty term since it is added in terms of sin 2 2θ 23 instead of θ 23 . Therefore, we also look for compatible solutions around ∼ π/2 − θ 23 (the so-called octant degeneracy [90]) and keep the minimum of the χ 2 between the two.

Medium baseline reactor experiment: JUNO
We adopt an experimental configuration for the JUNO experiment based on Refs. [61,62,83], following the analysis described in Ref. [49]. We normalize the number of events such that for the default exposure of 20 kt × 36 GW × 6 yr = 4320 kt GW yr we obtain 10 5 events [61,83]. The energy resolution is assumed to be 3% 1 MeV/E. We perform a χ 2 analysis using 350 bins for the energy spectrum. This number is chosen sufficiently large such that bins are smaller (or of the order of) the energy resolution. We take into account an overall normalization uncertainty of 5% and a linear energy scale uncertainty of 3%. Uncertainties in the oscillation parameters sin 2 θ 13 and sin 2 θ 12 are included as pull parameters in the χ 2 using true values and uncertainties according to Eq. (B.1), while |∆m 2 31 | is left free when fitting the data. For this parameter a dense grid is computed and the minimum is manually searched for. We have updated the analysis from Ref. [49] by taking into account the precise baseline distribution of 12 reactor cores as given in Tab. 1 of Ref. [62] (including also the Daya Bay reactors at 215 and 265 km). This reduces T 0 by about 5 units compared to the idealized situation of a point-like source at 52.47 km (the latter being the power averaged distance of the 10 reactors not including the Daya Bay reactors). Adopting the same assumptions as in Ref. [62] we find for a 4320 kt GW yr exposure T 0 ≈ 11.8, which is in excellent agreement with their results, see red-dashed curve in Fig. 2 (right) of Ref. [62].
Our analysis ignores some possible challenges of the experiment, in particular the effect of a non-linearity in the energy scale uncertainty [70], see also Ref. [62,74]. While such issues have to be addressed in the actual analysis of the experiment, our analysis suffices to discuss the behavior of the relevant test statistic and sensitivity measures.

Atmospheric neutrino experiments: PINGU and INO
For the simulation of the ICal@INO experiment we use the same code as in Ref. [58], where further technical details and references are given. Here we summarize our main assumptions. We assume a muon threshold of 2 GeV and assume that muon charge identification is perfect with an efficiency of 85% above that threshold. As stressed in Refs. [53,54] the energy and direction reconstruction resolutions are crucial parameters for the sensitivity to the mass ordering. We assume here the "high" resolution scenario from Ref. [58], which corresponds to a neutrino energy resolution of σ Eν = 0.1E ν and neutrino angular resolution of σ θν = 10 • , independent of neutrino energy and zenith angle. More realistic resolutions have been published in Ref. [59]. While those results are still preliminary, we take our choice to be representative (maybe slightly optimistic), justified by the fact that we obtain sensitivities to the mass ordering in good agreement with Ref. [59]. With our assumptions we find 242 µ-like events per 50 kt yr exposure assuming no oscillations (sum of neutrino and anti-neutrino events) in the zenith angle range −1 < cos θ < −0.1. We divide the simulated data into 20 bins in reconstructed neutrino energy from 2 GeV to 10 GeV, as well as 20 bins in reconstructed zenith angle from cos θ = −1 to cos θ = −0.1. We then fit the two-dimensional event distribution in the 20 × 20 bins by using the appropriate χ 2 -definition for Poisson distributed data. Our default exposure for INO is a 50 kt detector operated for 10 yr.
For the PINGU simulation we use the same code as in Ref. [49], where technical details can be found. In particular, we adopt the same effective detector mass as a function of neutrino energy, with the threshold around 3 GeV, and the effective mass rises to about 4 Mt at 10 GeV and 7 Mt at 35 GeV. For the reconstruction abilities we assume that neutrino parameters are reconstructed with a resolution of σ Eν = 0.2E ν and σ θν = 0.5/ E ν /GeV. This corresponds to about 13 • (9 • ) angular resolution at E ν = 5 GeV (10 GeV). We stress that those resolutions (as well as other experimental parameters) are far from settled. With our choice we obtain mass ordering sensitivities in good agreement with Ref. [48], which are somewhat more conservative than the official PINGU sensitivities from Ref. [38]. For a 3 yr exposure and θ 23 = 45 • we obtain T 0 ≈ 7.5.
For both, INO and PINGU, we include the following systematic uncertainties: a 20% uncertainty on the over-all normalization of events, and 5% on each of the neutrino/antineutrino event ratio, the ν µ to ν e flux ratio, the zenith-angle dependence, and on the energy dependence of the fluxes. Moreover, in order to make the Monte Carlo simulation feasible we set ∆m 2 21 = 0, which implies that also θ 12 and the CP phase δ disappear from the problem. The validity of this approximation and/or the expected size of δ-induced effects has been studied for instance in Refs. [48,49,58,59]. Typically T 0 varies by roughly 1-2 units as a function of δ, which is small compared to uncertainties related to experimental parameters such as reconstruction abilities. We do not expect that δ and ∆m 2 21 related effects will change the statistical behavior of the test statistic T significantly, as also the results of Ref. [46] seem to indicate. The sensitivity of this type of experiments is largely dependent on the baseline and neutrino energies considered, which may vary widely from one setup to another. In this work we have studied three different setups, NOνA, LBNE-10 kt, LBNE-34 kt.
The first setup considered, NOνA [7,8], has a moderate sensitivity to the mass ordering, estimated to reach at most 3σ (see for instance Refs. [8,91]). The setup consists of a narrow band beam with neutrino energies around 2 GeV, aiming to a 13 kt Totally Active Scintillator Detector (TASD) placed at a baseline of L = 810 km. NOνA has recently started taking data. The beam is expected to reach 700 kW by mid-2014 [91], and by the end of its scheduled running time it will have accumulated a total of 3.6 × 10 21 PoT, equally split between π + and π − focusing modes. The detector performance has been simulated following Refs. [8,92]. Systematic errors are implemented as bin-to-bin correlated normalization uncertainties over the signal and background rates. These have been set to 5% and 10% for the signal and background rates, respectively, for both appearance and disappearance channels.
The second setup considered in this work is the LBNE proposal [10,11]. LBNE would use a wide band beam with an energy around 2-3 GeV and a baseline of L = 1300 km. The first phase of the project (dubbed in this work as LBNE-10 kt) consists of a 10 kt Liquid Argon (LAr) detector placed on surface. In a second stage, dubbed in this work as LBNE-34 kt, the detector mass would be upgraded to 34 kt and placed underground. The longer baseline and higher neutrino energies make this setup more sensitive to the mass ordering: in its first stage is already expected to reach at least a significance between 2.5 − 7σ, depending on the value of δ. The results also depend significantly on the assumptions on systematics and the beam design, see for instance Ref. [11]. In this work, the detector performance has been simulated according to Ref. [10]. Systematic uncertainties have been set at the 5% level for both signal and background rates in the appearance channels, and at the 5% (10%) for the signal (background) rates in the disappearance channels. Tab. V shows the expected total event rates in the appearance channels for each of the long baseline setups considered in this work. It should be noted the difference in statistics between the LBNE-10 kt and LBNE-34 kt, which is not only due to the larger detector mass but also to a different neutrino beam design. The first stage of the project, LBNE-10 kt, is simulated using the fluxes from the October 2012 Conceptual Design Report, Ref. [10], while for the upgraded version, LBNE-34 kt, we consider the fluxes from Ref. [9]. In both cases the beam power is set to 700 kW.
The simulations for the long baseline beam experiments have been performed using GLoBES [93,94]. In order to generate random fluctuations in the number of events, version 1.3 of the MonteCUBES [95] software was used. In addition to the true values and prior uncertainties for the oscillation parameters given in Eq. (B.1), a 2% uncertainty on the matter density is also considered.