JUNO's prospects for determining the neutrino mass ordering

The flagship measurement of the JUNO experiment is the determination of the neutrino mass ordering. Here we revisit its prospects to make this determination by 2030, using the current global knowledge of the relevant neutrino parameters as well as current information on the reactor configuration and the critical parameters of the JUNO detector. We pay particular attention to the non-linear detector energy response. Using the measurement of $\theta_{13}$ from Daya Bay, but without information from other experiments, we estimate the probability of JUNO determining the neutrino mass ordering at $\ge$ 3$\sigma$ to be 31% by 2030. As this probability is particularly sensitive to the true values of the oscillation parameters, especially $\Delta m^2_{21}$, JUNO's improved measurements of $\sin^2 \theta_{12}$, $\Delta m^2_{21}$ and $|\Delta m^2_{ee}|$, obtained after a couple of years of operation, will allow an updated estimate of the probability that JUNO alone can determine the neutrino mass ordering by the end of the decade. Combining JUNO's measurement of $|\Delta m^2_{ee}|$ with other experiments in a global fit will most likely lead to an earlier determination of the mass ordering.


I. INTRODUCTION
After the first observation of the so-called solar neutrino puzzle by the Homestake experiment in the late 60's, it took us about 30 years to establish that neutrino flavor oscillations are prevalent in nature, impacting cosmology, astrophysics as well as nuclear and particle physics. In the last 20 years we have consolidated our understanding of neutrino oscillations both at the experimental as well as the theoretical level. We know now, thanks to a great number of experimental efforts involving solar, atmospheric, accelerator and reactor neutrino oscillation experiments, that neutrino oscillations are genuine three flavor phenomena driven by two independent mass squared differences (∆m 2 21 and ∆m 2 32 ) and three mixing angles (θ 12 , θ 13 and θ 23 ) and possibly a charge-parity violating phase (δ CP ). See [1] for a review of how these parameters are defined.
Today a single class of experiments dominates the precision of the measurement of each of these aforementioned parameters [2][3][4]. In the solar sector (12), ∆m 2 21 is determined to less than 3% mainly by the KamLAND [5] long-baselineν e disappearance reactor experiment, while sin 2 θ 12 is determined by the combination of solar neutrino experiments [6][7][8][9][10][11][12][13][14][15][16][17] to ∼ 4%. In the atmospheric sector (23), ∆m 2 32 (or ∆m 2 31 ) and sin 2 θ 23 are dominantly determined by the ν µ andν µ disappearance accelerator experiments MINOS [18], NOvA [19,20] and T2K [21], with corresponding precision of better than 1.5% and 8%, respectively. The mixing sin 2 θ 13 , that connects the solar and atmospheric sectors, is determined by the shortbaselineν e disappearance reactor experiments Daya Bay [22], RENO [23,24] and Double Chooz [25] to a precision of ∼ 3%. Regarding the CP-phase δ CP there is a small tension in the determination among current experiments T2K and NOvA [2,4,26]. The determination of δ CP remains an open problem that probably will have to be addressed by the next generation of long-baseline neutrino experiments such as DUNE and Hyper-K. There is, nevertheless, an important open question that influences the better determination of some of these parameters: what is the neutrino mass ordering?
If one defines the mass eigenstates ν 1 , ν 2 , ν 3 in terms of decreasing amount of electron neutrino flavor content, then the results of the SNO solar neutrino experiment determined that m 1 < m 2 . However, the available information does not allow us to know the complete ordering yet: both m 1 < m 2 < m 3 (normal ordering, NO) and m 3 < m 1 < m 2 (inverted ordering, IO) are compatible with the current data [2,4,26]. The measurement of the neutrino mass ordering is one of the most pressing and delicate challenges of our times. Besides its direct impact on the precise knowledge of the oscillation parameters, neutrino mass ordering affects the sum of neutrino masses from cosmology, the search for neutrinoless double-β decay and ultimately, our better understanding of the pattern of masses and mixing in the leptonic sector.
The use ofν e from nuclear reactors with a medium-baseline detector to determine the mass ordering, exploring genuine three generation effects as long as sin 2 θ 13 few %, was first proposed in [27]. This idea was further investigated in [28] for a general experiment and more recently in [29], specifically for JUNO. In all three of these papers, different artificial constraints were imposed on the ∆m 2 3i 's when comparing the NO and IO spectra. As we will see, any and all of these artificial constrains increases the difference between the NO and IO, see appendix A for a more detailed discussion. In fact, it was shown in Ref. [30] that what these experiments can precisely measure is the effective combination [31] ∆m 2 ee ≡ cos 2 θ 12 ∆m 2 31 + sin 2 θ 12 ∆m 2 32 , and the sign (+ for NO, − for IO) of a phase (Φ ) that depends on the solar parameters. This subtlety is of crucial importance in correctly assessing the sensitivity to the neutrino mass ordering. The Jiangmen Underground Neutrino Observatory (JUNO) [32], a 20 kton liquid scintillator detector located in the Guangdong Province at about 53 km from the Yangjiang and Taishan nuclear power plants in China, will be the first experiment to implement this idea. This medium-baseline facility offers the unprecedented opportunity to access in a single experiment four of the oscillation parameters: ∆m 2 21 , sin 2 θ 12 , |∆m 2 ee |, sin 2 θ 13 and the sign of phase advance, Φ (L/E), which determines the mass ordering. JUNO aims in the first few years to measure ∆m 2 21 , sin 2 θ 12 and |∆m 2 ee | with a precision 1% to be finally able, after approximately 8 years 1 , to determine the neutrino mass ordering at 3σ confidence level (C.L.).
Many authors have studied the neutrino mass ordering determination at medium-baseline reactor neutrino experiments such as JUNO. In Ref. [33] a Fourier analysis was proposed, but no systematic effects were considered. The effect of energy resolution was investigated in Ref. [34]. The importance of also taking into account non-linear effects in the energy spectrum reconstruction was first pointed out in Ref. [35] and addressed in Ref. [36,37], where limited impact on the mass ordering was observed. Matter effects, geo-neutrino background, energy resolution, energy-scale and spectral shape uncertainties were investigated in [38]. The impact of the energy-scale and flux-shape uncertainties was further explored in [39]. The benefits of a near detector for JUNO was demonstrated in [40] and in [41] the impact of the sub-structures in the reactor antineutrino spectrum, due to Coulomb effects in beta decay, was studied in the light of a near detector under various assumptions of the detector energy resolution. This was further explored in [42]. In [43] the distribution for the test statistics, to address the mass ordering determination, was proven to be normally distributed and this was also applied to quantify the JUNO sensitivity. It was also shown that without statistical fluctuations, the mentioned test statistics is equivalent to the widely adopted ∆χ 2 approach used in sensitivity studies. Finally, the combined sensitivity of JUNO and PINGU was also recently studied by the authors of [44], while a combined sensitivity study of JUNO and T2K or NOvA was performed in [45].
One can appreciate the difficulty in establishing the mass ordering with this setup by noticing that after 8 years (2400 days) of data taking, the difference in the number of events for NO and IO is only a few tens of events per energy bin which is smaller than the statistical uncertainty in each bin. It is clear that this formidable endeavor depends on stringent requirements on the experiment's systematic uncertainties, but also on the actual values of the oscillation parameters as well as on statistical fluctuations. This is why we think it is meaningful to revisit the prospect that JUNO can obtain a 3σ preference of the neutrino mass ordering by 2030. This is the task we undertake in this paper.
Our paper is structured as follows. In Sec. II we describe theν e survival probability in a way that highlights the physics that is relevant for medium-baseline reactor neutrino experiments and how it depends on the oscillation parameters. In Sec. III we explain how we simulate the experiment and show the statistical challenges associated with extracting the mass ordering in a medium-baseline reactor experiment. Sec. IV addresses how the following experimental details affect the determination power of the neutrino mass ordering of the JUNO experiment; (A) reactor distribution and backgrounds, (B) bin to bin flux uncertainties, (C) the number of energy bins used in the analysis, (D) the size of the energy resolution. In Sec. V we show how varying the true values of the neutrino oscillation parameters improves or reduces the prospects for JUNO's determination of the neutrino mass ordering. Sec. VI addresses the effects of the non-linear detector response on the mass ordering determination. In Sec. VII we simulate 60 k experiments consistent with the current best fit values and uncertainties of the oscillation parameters. From this simulation we estimate the probabilities that JUNO can determine the mass ordering at ≥ 3σ with 4, 8 and 16 years of data taking. In Sec. VIII we show how JUNO's measurement of ∆m 2 ee , when combined with other experiments can determine the mass ordering after a few years of data taking. Finally in Sec. IX we draw our conclusions. There are four Appendices: Appendix A addresses the effects of imposing artificial constraints on the ∆m 2 's, Appendix B gives a derivation of the oscillation probability used in this paper, Appendix C compares our analysis with the JUNO collaboration's analysis and Appendix D discusses the current impact of T2K, NOvA and the atmospheric neutrino data on the determination of |∆m 2 ee |.

II. THEν e SURVIVAL PROBABILITY
The neutrino survival probability for reactor experiments in vacuum is given by where the kinematic phases are ∆ ij ≡ ∆m 2 ij L/(4E) and P = sin 2 2θ 12 cos 4 θ 13 sin 2 ∆ 21 . This survival probability was first rewritten, without approximation, in a more useful way for the medium baseline reactor experiments in [30], as where ∆m 2 ee , defined in Eq. (1), is the effective atmospheric ∆m 2 for ν e disappearance, see [31,46]. The mass ordering is determined by the sign in front of Φ , '+' ('-') for NO (IO). The phase advance or retardation Φ is Φ = arctan (cos 2θ 12 tan ∆ 21 ) − ∆ 21 cos 2θ 12 . (4) Note that the survival probability depends only on four of the oscillation parameters, θ 12 , θ 13 , ∆m 2 21 and |∆m 2 ee | and the sign for the mass ordering. The determination of the sign of ∆m 2 21 cos 2θ 12 > 0 by SNO [47] is crucial for this measurement. See Appendix B for more details on the survival probability.
This behavior is illustrated in Fig. 1, using the central and 1σ bands for the solar parameters given in Table I taken from the recent global fit [2]. Here, we show the advance/retardation as a function of E at L = 52.5 km, and in Appendix B also as function of L/E.   Table I. Φ as function of L/E is given in Appendix B.
Matter effects are important in the determination of the best fit values of the solar parameters ∆m 2 21 and sin 2 θ 12 . The size of these effects, which do not satisfy the naive expectations, was first given in a numerical simulation in [48], and later explained in a semianalytical way in [49]. For ∆m 2 21 and sin 2 θ 12 , the sizes of these shifts are -1.1% and 0.20%, respectively. Since we are interested only in sensitivities, we can ignore matter effects in the propagation of neutrinos in this paper and will use here the vacuum expression for the survival probability, Eq. (3). However, in a full analysis of real data matter effects must be included. In the next section we will describe details of our simulation of the JUNO reactor experiment.

III. SIMULATION OF A MEDIUM BASELINE REACTOR EXPERIMENT
For the simulation of JUNO we use the information given in Refs. [32,44] but with the updates of baselines, efficiencies and backgrounds provided in Ref. [50]. In order to simulate the event numbers and to perform the statistical analysis we use the GLoBES software [51,52]. We start with an idealized configuration where all reactors that provide the 26.6 GW th total thermal power are at 52.5 km baseline from the detector. The antineutrinos are mainly created in the fission of four isotopes, 235 U (56.1%), 238 U (7.6%), 239 Pu (30.7%) and 241 Pu (5.6%) [53]. For our simulation we use the Huber-Mueller flux predictions [54,55] for each isotope. Theν e propagate to the JUNO detector and are observed via inverse beta decay ν e + p → e + + n [56]. We assume a liquid scintillator detector with a 20 kton fiducial mass and a running time of 2400 days (8 years @ 82% live time) 2 . We will include the detector energy resolution of 3.0% unless otherwise stated. The aforementioned quantities affect the calculation of the event numbers. The number of events, N i , in the i-th bin corresponding to the reconstructed neutrino energy E i is given by Here, N T is a normalization constant taking into account the exposure time, efficiency, fiducial mass of the detector and reactor-detector distance, φ νe (E) is the antineutrino flux, which relates the reconstructed and true neutrino energies. The energy resolution is given by where the prompt energy E p is given by The variable is the detector energy resolution. In this paper we will use = 3.0% except when discussing the effects of varying this parameter in Sec. IV D, where we also will use 2.9% and 3.1%. In Fig. 2 we have plotted the event spectrum for JUNO using 200 bins for 8 years (2,400 live days) of data taking and 26.6 GW th . In the top panel, the blue and red spectra corresponds to ∆m 2 ee [NO] = 2.530 × 10 −3 eV 2 and ∆m 2 ee [IO] = −2.548 × 10 −3 eV 2 , respectively 3 . The ∆m 2 ee for NO is input whereas the value for IO is chosen so as to minimize the ∆χ 2 = χ 2 min [IO] − χ 2 min [NO] between the two spectra. By construction χ 2 min [NO] = 0, so minimizing ∆χ 2 is equivalent to minimizing χ 2 min [IO]. In the lower panel, we plot the difference in event spectra obtained for NO and IO. Note that this difference is less than 20 events/bin. Also shown is the statistical uncertainty in each oscillated bin (orange band), which for all bins exceeds the difference between the NO and IO event spectra 4 . This figure demonstrates the statistical challenges for JUNO to determine the mass ordering and will be addressed in more detail in Sec. VII. For reference, we also show on the right panel of Fig. 2 the corresponding χ 2 distributions. Throughout this paper we will use dashed (solid) lines for the fit with NO (IO). 2 An exposure of 26.6 GW th for 2400 days (8 years @ 82%) is equivalent to 35.8 GW th for 1800 days (6 years @ 82%) as used in [32]. 3 Note, that the value for ∆m 2 ee [IO] does not correspond to any of the artificial constraints on the atmospheric mass splitting imposed in Refs. [27][28][29], see Appendix A for more details. 4 Caveat: if one halves the bin size in this figure the difference between the NO and IO goes down by a factor of 2 whereas the statistical uncertainty by only √ 2 making the difference more challenging to observe. If one doubles the bin size the statistical uncertainty increases by the √ 2 whereas the difference would increase by 2, this improves the situation except for the fact that at low energy there is some washing out of the difference. Events/bin i , is given, as well as plus/minus statistical uncertainty in each oscillated bin (orange band), ± N NO i ≈ ± N IO i . Note, the difference is always within the statistical uncertainty for that bin.
Note that including systematic uncertainties as well as the real distribution of core-reactor distances and backgrounds will further decrease the difference between the two spectra. But first let us address the simulation details and systematic uncertainties.
To perform the statistical analysis we create a spectrum of fake data N dat i for some set of oscillation parameters. Next we try to reconstruct this spectrum varying the relevant oscillation parameters p. For each set p we calculate a χ 2 function where N i ( p, α) is the predicted number of events 5 for parameters p, α = (α 1 , α 2 , . . .) are the systematic uncertainties with their corresponding standard deviations σ k . χ 2 NL is the penalty for the non-linear detector response and will be discussed in more detail in Sec. VI.
As in Ref. [32], we included systematic uncertainties concerning the flux, the detector efficiency (which are normalizations correlated among all bins, i.e. N i → αN i ) and a binto-bin uncorrelated shape uncertainty. The shape uncertainty is simply introduced as an independent normalization for each bin in reconstructed energy, i.e. N i → α i N i .
In the next section we will discuss in detail how some experimental issues can affect JUNO's ability to determine the neutrino mass ordering 6 . We will concentrate on the impact of the real reactor core distribution, the inclusion of background events, the bin to bin flux uncertainty, the number of equal-size energy bins of data and the detector energy resolution. We leave the discussion of the dependence on the true value of the neutrino 5 The number of events includes the background events extracted from Ref. [32]. 6 For a verification of our simulation, see Appendix C. oscillation parameters, on the non-linearity of the detector energy response and on statistical fluctuations for later sections.

IV. MEAN (OR AVERAGE) DETERMINATION OF THE NEUTRINO MASS ORDERING
In the following subsections we will discuss in which way the following quantities affect the determination power of the neutrino mass ordering of the JUNO experiment: Unless otherwise stated, we generate fake data fixing the neutrino oscillation parameters as in Tab. I and assume the nominal values for the energy resolution, number of data bins and total exposure for JUNO given in Tab TABLE II. Nominal values, as well as lowest and largest values, assumed in this paper for the JUNO energy resolution, systematic uncertainties (b2b=bin to bin and the energy scale bias σ bias ), number of energy data bins and exposure. One year is 300 days of live time.

A. Effect of the Reactor Distribution and Backgrounds
The real position of the reactor cores and background events are expected to impact JUNO's sensitivity. Fig. 3 shows the reduction in ∆χ 2 as one goes from the ideal reactor core-detector disposition (all cores at 52.5 km) with no backgrounds included to the real reactor core-detector baseline distribution given in Table III with all backgrounds taken into account. The blue lines, labeled "ideal, wo BG", are the same as on the right panel of Fig. 2.
There are two types of background events at JUNO: one from remote reactors (Daya Bay and Huizhou) and the other includes accidental events, cosmogenic decays and geo-neutrinos. The first we compute, the latter we take from [32].
Notice the χ 2 min [IO] goes from 14.5 (ideal, wo BG) down to 9.1 (real, all BG), a decrease of more than 5 units. The real core positioning alone causes a reduction in sensitivity of 2.8 and the background events an extra 2.6 (1.8 from DB and HZ). We use the real baseline distribution and include all backgrounds in the rest of this paper.   [50]. The total power is 26.6 GW th . The remote reactors Daya Bay (DB) and Huizhou (HZ) produce background events for the neutrino mass ordering.

B. Effect of bin to bin Flux Uncertainties
There is uncertainty related to the exact shape of the reactorν e flux, inherent to the flux calculation. This uncorrelated bin to bin (b2b) shape uncertainty is included in our analysis by varying each predicted event bin with a certain penalty. The primary purpose of the TAO near detector is to reduce this bin to bin shape uncertainty, see [53].
The effect of this systematic bias is shown in Fig. 4. The lines labeled "stat only" is the same as the one labeled "real, all BG" in Fig. 3. We find χ 2 min [IO] = 8.5, 7.1 and 5.6, respectively, for 1%, 2% and 3%. When the b2b systematic uncertainty is not included, we recall, χ 2 min [IO] = 9.1. So if the shape uncertainty is close to 1% (the nominal value), the sensitivity to the neutrino mass ordering is barely affected. However, for 2% and 3% we see a clear loss in sensitivity. This is because increasing the uncorrelated uncertainty for each bin, makes it easier to shift from a NO spectrum into an IO one and vice versa. We use 1% b2b in the rest of the paper.  [53].

C. Effect of varying the number of Energy Bins
Here we examine the impact on χ 2 of changing the size of the neutrino energy bins in the range [1.

D. Effect of varying the Energy Resolution
Next we consider variations of the detector energy resolution. In this section we assume, that this number can be slightly better and slightly worse than the nominal 3.0%. Small variations of the resolution have large impacts on the determination of the neutrino mass ordering, as shown in Fig. 6. The red line corresponds to the nominal energy resolution of 3.0%. The blue and green lines are obtained for 2.9% and 3.1%, respectively, with corresponding χ 2 min [IO] = 9.7 and 7.5. Clearly χ 2 min [IO] is quite sensitive to the exact value of the resolution that will be achieved by JUNO. Therefore even a small improvement on the energy resolution would have a sizable impact on the determination potential of the neutrino mass ordering. However, it appears challenging for JUNO to reach an energy resolution even slightly better than 3.0%, see [57]. Note the red lines in Fig. 6 (3.0% res.) are also the same as the blue lines in Fig. 4 (1% b2b). We always use 3.0% resolution elsewhere in this paper.

V. EFFECT OF VARYING THE TRUE VALUES OF THE NEUTRINO OSCIL-LATION PARAMETERS
In this section we explore how varying the true values of the neutrino oscillation parameters improves or reduces the prospects for JUNO's determination of the neutrino mass ordering. We first consider the variation of single parameters with the others held fixed and then consider the correlations varying both ∆m 2 21 and sin 2 θ 12 with ∆m 2 ee and sin 2 θ 13 held fixed and vice versa.
We start by creating fake data sets using the upper and lower 1σ bounds obtained in Ref. [2] (see Tab. I), always for one parameter at the time. The result of these analyses is shown in Fig. 7, where in each panel we vary one of the parameters as indicated. Here again solid (dashed) lines are used for IO (NO). As can be seen, changes in any of the oscillation parameters can have large effects on the determination power of the neutrino mass ordering. Especially remarkable is the effect of a smaller ∆m 2 21 , which shifts χ 2 min [IO] from 8.5 to 7.1. Note that the best fit value obtained from the analysis of solar neutrino data from Super-K [15] is even smaller than the one considered here and therefore the determination would then be even more difficult. On the other hand side, a larger value of the solar mass splitting improves significantly the determination of the mass ordering. In this case we obtain χ 2 min [IO] = 10.2. The effect of the other parameters is not as pronounced as in the case of the solar mass splitting, but still appreciable: ∆m 2 ee /sin 2 θ 13 / sin 2 θ 12 , within 1σ of their current best fit value, can move χ  In Fig. 8 we show the correlated variation of the χ 2 min [IO] as a function of (sin 2 θ 12 , ∆m 2 21 ) holding (sin 2 θ 13 , ∆m 2 ee ) fixed as well as a function of (sin 2 θ 13 , ∆m 2 ee ) holding (sin 2 θ 12 , ∆m 2 21 ) fixed. Even varying these parameters within 3σ of their current best fit, there are very significant changes to the χ 2 min [IO] contour plots. This implies that JUNO's prospect for the determination of the neutrino mass ordering could be improved or weakened by Nature's choice for the true values of these oscillation parameters. The values that were used in [37] (also in [32]) are shown by the gray stars in these figures.

VI. NON-LINEAR DETECTOR ENERGY RESPONSE
In a liquid scintillator detector, the true prompt energy, E p , (positron energy plus m e ) is not a linear function of to the visible energy, E vis , in the detector. The main energydependent effects are the intrinsic non-linearity related to the light emitting mechanisms (scintillation and Cherenkov emission) and instrumental non-linearities. The non-linear detector response can be modeled by a four parameter function [57][58][59] which relates the true prompt energy to the visible detector energy according to and the coefficients (a 1 , a 2 , a 3 , a 4 ) are determined by prompt energy calibration techniques. We use the prompt energy scale calibration curve shown in Fig. 1 of Ref. [57], which can be well described by the f NL given in Eq. (11) with a 1 = 1.049,ā 2 = 2.062 × 10 −4 ,ā 3 = 9.624 × 10 −2 ,ā 4 = 1.184 .
Then the true neutrino energy, E, is then constructed by E = E p + ∆M .
To allow for deviations from this calibration, we use in our simulation the reconstructed prompt energy, E p , given by 1 , a 2 , a 3 , a 4 ; E p ) .
Note, with this definition E vis is held fixed as we change the a i 's from their calibration values,ā i . In the simulation, we generate a distribution of E p 's for the true mass ordering and use Eq. (12) to generate a distribution of E p 's for the test mass ordering 7 . The allowed range of the a i 's is constrained by including a penalty term for the derivation of when fitting the simulated spectra to the test mass ordering. Explicitly, we allow the a i 's to vary from their calibration values and then penalize the fit by using the simplified χ 2 NL defined, as in Ref. [39], as where {a i , i = 1, ..., 4} are the best fit of these parameters for the test mass ordering spectra, σ bias is the uncertainty on the energy scale, max Ep indicates that we take only the maximal difference which happens at E ∼ 2.75 MeV (see Fig. 9). We consider the following sizes for the bias, σ bias = 0.0, 0.2, 0.4, 0.7%, as well as no penalty, i.e. χ 2 NL = 0. Using Fig. 2 of Ref. [57], we see that JUNO expects an approximately energy independent systematic uncertainty on the energy scale of about 0.7%, mainly due to instrumental non-linearity and position dependent effects. Therefore, σ bias = 0.7% is our nominal value from here on.
On the left panel of Fig. 9 we show the ratio E p /E p as a function of E for the corresponding f NL coefficients listed in Tab. IV obtained for the best fit to IO of the NO input spectra. On the right panel we see the effect of the uncertainty on the energy scale on χ 2 . In particular, as the uncertainty increases χ 2 min [IO] goes from 8.5 (no NL effect) down to 8.0 (0.2%), 7.5 (0.4%) and 7.2 (0.7%). Note that if we introduce the non-linearity shift with no penalty χ 2 min [IO] = 6.8. Even with the nominal 0.7% bias, this is a significant effect, reducing χ 2 min [IO] by more than 1 unit (8.5 to 7.2) and in this manner further lowering the mass ordering discrimination power.
We also observe that the precision on the determination of |∆m 2 ee | is notably degraded when the non-linearity in the energy scale is included. In addition, the best fit value for  9. On the left panel we show the ratio between the reconstructed prompt energy E p and the true prompt energy E p as a function of the neutrino energy E, for σ bias = 0.2% (yellow), 0.4% (green) and 0.7% (red) and no penalty (magenta) for the best fit to IO for the NO spectra. In blue we show the line for perfect reconstruction (no NL) as a reference. On the right panel we show the changes to χ 2 caused by the addition of the corresponding χ 2 NL .

VII. FLUCTUATIONS ABOUT THE MEAN FOR THE NEUTRINO MASS OR-DERING DETERMINATION
It has been already pointed out that statistical fluctuations are important for JUNO, see for instance Ref. [34] where they estimate the statistical uncertainty on ∆χ 2 by an analytical expression and a Monte Carlo simulation. The calculation was performed just after the first measurement of sin 2 θ 13 by RENO and Daya Bay, under different detector resolution and systematic assumptions. It is timely to reevaluate this here.
We have already shown in Fig. 2 that the difference between the spectra for NO and IO is smaller than the statistical uncertainty in each bin. We consider here the effects of fluctuating the number of events in each bin. We evaluate the impact of this fluctuations on the mass ordering determination by performing a simulation of 60000 JUNO pseudoexperiments for each exposure and obtain the distributions given in Fig. 10. To generate this figure, we create a fake data set {N 0 i , i = 1, ..., N bins } using the neutrino oscillation parameters in Tab Fig. 10. These distributions are Gaussian (as was proven analytically in Ref. [43]) with corresponding central values ∆χ 2 = 3.4, 6.7 and 12.4 and standard deviations 3.4, 4.7 and 6.1, respectively. Our pseudo-experiments reveal that after 8 years in only 31% of the trials JUNO can determine the neutrino mass ordering at the level of 3σ or better. We also find that there is even a non negligible probability (∼8%) to obtain the wrong mass ordering, i.e., ∆χ 2 < 0. For a shorter (longer) exposure of 4 (16) years, 5% (71%) of the pseudo-experiments rule out IO at 3σ or more. In these cases in about 16% (2%) of the trials the IO is preferred.

VIII. COMBINING JUNO WITH THE GLOBAL FIT
In the previous section we have shown that the significant impact of statistical fluctuations on top of the detector systematic effects, can make it very challenging for JUNO by itself to determine at 3σ or more the neutrino mass ordering even after 16 years. However, as was shown in [31], muon disappearance experiments measure 8 ∆m 2 µµ ≡ sin 2 θ 12 ∆m 2 31 + cos 2 θ 12 ∆m 2 32 , whose relationship to |∆m 2 ee | is given by the positive (minus) sign is for NO (IO). Therefore, by using muon disappearance measurements we have a constraint on the allowed |∆m 2 ee |'s for the two mass orderings, i.e. . Of course, the measurement uncertainty on |∆m 2 µµ | must be smaller than this 3.1% difference for this measurement to impact the confidence level at which the false mass ordering is eliminated. The short baseline reactor experiments, Daya Bay and RENO, measure the same |∆m 2 ee | for both orderings with uncertainties much larger than JUNO's uncertainty.
This physics is illustrated in Fig. 11 where we show the allowed region in the plane ∆m 2 21 versus |∆m 2 ee | by JUNO for NO (blue) and IO (red) after 2 years of data taking and the corresponding 1σ CL allowed region by the current global fit constraint on |∆m 2 µµ |. We see that the global fit and JUNO NO regions overlap while the corresponding IO regions do not. This tension between the position of the best fit values of |∆m 2 ee | for IO with respect to NO gives extra leverage to the data combination.
Therefore, combining JUNO's measurement of |∆m 2 ee | with other experiments, in particular T2K and NOvA, expressed by the current global fits, see Refs. [2][3][4], turns out to be very powerful in unraveling the neutrino mass ordering at a high confidence level, as shown in the left panel of Fig. 12 for 2 years of JUNO data. As we can see χ 2 min [IO] combined (green solid line) turns out to be about 16. As a result with only two years of JUNO data taking the mass ordering is determined at better than 3σ in 99% of the trials, see right panel of Fig. 12 FIG. 11. The ellipses are the allowed regions for JUNO in the ∆m 2 21 versus |∆m 2 ee | plane for NO (blue, 2 and 3σ CL) and IO (red, 2 and 3σ CL) after 2 years. The best fit for NO (IO) is depicted by a black star (dot). We assume NO here and the ∆χ 2 's are determined with respect to NO best fit point. We use the neutrino oscillation parameters at the values given in Tab. I and take into account the experimental nominal systematic uncertainties and energy resolution given in Tab. II. We also show, as red (for IO) and blue (for NO) bands, the 1σ CL allowed regions by the current global fit constraint on |∆m 2 µµ |. Note, the not overlap for the allowed regions for IO. Appendix D we discuss the separate contributions from T2K, NOvA and the atmospheric neutrino data (Super-Kamiokande and DeepCore) to the χ 2 distribution for the global fit determination of |∆m 2 ee | and the corresponding impact on the combination with JUNO, for completeness and comparison with Fig. 5 of Ref. [45].
So even though JUNO cannot determine the ordering alone, a couple of years after the start of the experiment, it's precise measurement of |∆m 2 ee | will allow us to know the mass ordering at better than 3σ when the measurement on |∆m 2 µµ | from other neutrino oscillation experiments is combined in a global analysis.

IX. CONCLUSIONS
The neutrino mass ordering is one of the most pressing open questions in neutrino physics. It will be most likely measured at different experiments, using atmospheric neutrinos at ORCA [60,61], PINGU [62][63][64], Hyper-K [65] or DUNE [66] or accelerator neutrinos at T2HK [67] or DUNE [68]. It also is a flagship measurement for the up-coming JUNO experiment. This is why we have examined here in detail the impact of various factors on the determination power of the neutrino mass ordering by JUNO.
We have assumed NO as the true mass ordering, but our general conclusions do not depend on this assumption. In this case the power of discrimination can be encoded on the value of χ 2 min [IO], the larger it is the larger the confidence level one can discriminate between the two mass orderings using JUNO.
We have determined that the real reactor distribution and backgrounds account for a reduction in sensitivity of more than 5 units (i.e. χ 2 min [IO] going from 14.5 to 9.1), the bin to bin flux uncertainty, at its nominal value of 1%, to an extra reduction of 0.6 down to χ 2 min [IO] = 8.5, both assuming 3% energy resolution and 200 energy bins. Note that an improvement on the energy resolution from 3% to 2.9%, a challenging feat to achieve, would represent an increase of χ 2 min [IO] from 8.5 to 9.7. The values of neutrino oscillation parameters that will impact JUNO's measurement are currently known within a few % uncertainty. We have determined the effect of these uncertainties on the mass ordering discrimination. We remark, in particular, the influence of the true value of ∆m 2 21 , a smaller (larger) value than the current bet fit could shift χ 2 min [IO] from 8.5 to 7.1 (10.2). Another important factor is the non-linear energy response of the detector. Assuming a bias of 0.7% we have verified that this would decrease χ 2 min [IO] further from 8.5 to 7.2.
We have also examined the consequence of statistical fluctuations of the data by performing 60 k Monte Carlo simulated JUNO pseudo-experiments. Using them we have determined that after 8 (16) years in only 31% (71%) of the trials JUNO can determined the neutrino mass ordering at 3σ or more. This means that JUNO by itself will have difficulty determining the mass ordering. However, JUNO can still be used for a plethora of different interesting physics analysis [32,[69][70][71][72][73][74][75][76][77][78][79][80]. In particular, JUNO will be able to measure ∆m 2 21 , sin 2 θ 12 and |∆m 2 ee | with unmatched precision. This will be very useful to improve our understanding of the pattern of neutrino oscillations and to guide future experiments.
Finally, this inauspicious prediction is mitigated by combining JUNO's |∆m 2 ee | measurement into the current global fits, in particular the measurement of |∆m 2 µµ |. As we have shown, this combination will most likely result in the determination of the mass ordering at better than 3σ with only two years of JUNO data. Our conclusion for the global fits result is consistent with the results of [45]. So we can predict that in approximately two years after the start of JUNO we will finally know, via global analyses, the order of the neutrino mass spectrum, i.e. whether the lightest neutrino mass eigenstate has the most ν e (ν 1 ) or the least ν e (ν 3 ).   Table I.
Since Ω appears only as cos Ω, one could use Ω = 2|∆ ee | ± Φ as in Eq. (3). The factor 1 − sin 2 2θ 12 sin 2 ∆ 21 in front of cos Ω in Eq. (B7), modulates the amplitude of the θ 13 oscillations as this factor varies from 1 to cos 2θ 12 ≈ 0.4 as ∆ 21 goes from 0 to π/2. So the (· · · ) modulates the amplitude and Φ modulates the phase of the θ 13 oscillations.  [44], using the oscillation parameters and technical details of that reference. Our code, written for this paper, gives the solid lines whereas the results extracted from the above reference are dashed lines, normal (inverted) ordering is in blue (red).
Appendix D: On the contribution to the determination of |∆m 2 ee | from the |∆m 2 µµ | sensitive experiments It is informative to examine the contributions of the |∆m 2 µµ | sensitive experiments included in the global fit to the final determination of |∆m 2 ee |. We will focus here on the major players: T2K, NOvA and the atmospheric neutrino oscillation experiments Super-Kamiokande and DeepCore (ATM). The analyses of T2K, NOvA and ATM data shown in this section correspond to the analyses performed in Ref. [2]. For this purpose we show in Fig. 15  From these plots we see that T2K prefers |∆m 2 ee [NO, IO]| closer to the global fit best fit values, while NOvA (ATM) prefers lower (higher) values. Note that both accelerator neutrino oscillation experiments, however, prefer |∆m 2 ee [IO]| smaller than the value JUNO will prefer (NO assumed true). Since none of the χ 2 distributions are very Gaussian at this point, the combined χ 2 min [IO] is a result of broad distributions pulling for different minima that at JUNO's best fit value for |∆m 2 ee [IO]| contribute to an increase of χ 2 min [IO] of about 7 (NOvA), 3 (T2K) and 5 (ATM) units, resulting on the final power of the combination.
The addition of the atmospheric data, and also to a minor extent of MINOS data (which is compatible with NOvA), to the global fit used in this paper explains the difference of about 4 units in the predicted boost for the determination of the mass ordering we show here with respect to what is predicted in Fig. 5  and Super-K and DeepCore atmospheric data, labeled ATM, (lower panel) to the χ 2 fit of |∆m 2 ee | to NO (dashed lines) and IO (solid lines) included in the global fit (blue) and in the combination of the current global fit with 2 years of JUNO data (green). JUNO fit only is in red.