Synergies and prospects for early resolution of the neutrino mass ordering

The measurement of neutrino mass ordering (MO) is a fundamental element for the understanding of leptonic flavour sector of the Standard Model of Particle Physics. Its determination relies on the precise measurement of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta m^2_{31}$$\end{document}Δm312 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta m^2_{32}$$\end{document}Δm322 using either neutrino vacuum oscillations, such as the ones studied by medium baseline reactor experiments, or matter effect modified oscillations such as those manifesting in long-baseline neutrino beams (LB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document}νB) or atmospheric neutrino experiments. Despite existing MO indication today, a fully resolved MO measurement (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ge 5\sigma$$\end{document}≥5σ) is most likely to await for the next generation of neutrino experiments: JUNO, whose stand-alone sensitivity is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 3\sigma$$\end{document}∼3σ, or LB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document}νB experiments (DUNE and Hyper-Kamiokande). Upcoming atmospheric neutrino experiments are also expected to provide precious information. In this work, we study the possible context for the earliest full MO resolution. A firm resolution is possible even before 2028, exploiting mainly vacuum oscillation, upon the combination of JUNO and the current generation of LB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document}νB experiments (NOvA and T2K). This opportunity is possible thanks to a powerful synergy boosting the overall sensitivity where the sub-percent precision of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta m^2_{32}$$\end{document}Δm322 by LB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document}νB experiments is found to be the leading order term for the MO earliest discovery. We also found that the comparison between matter and vacuum driven oscillation results enables unique discovery potential for physics beyond the Standard Model.


Introduction
This document serves as a supplementary material to the work [1]. It describes some details regarding the analyses performed in this reference.

A. Empirical Reproduction of the χ 2 Function for the LBνB-II Experiments
In this section, we shall detail how we computed the number of events for T2K and NOvA. For a constant matter density, without any approximation, appearance oscillation probability for given baseline L and neutrino energy E, can be expressed [2] as P (ν µ → ν e ) = a ν + b ν cos δ CP + c ν sin δ CP , P (ν µ →ν e ) = aν + bν cos δ CP + cν sin δ CP , (1) where a ν , b ν , c ν , aν, bν and cν are some factors which depend on the mixing parameters (θ 12 , θ 23 , θ 13 , δm 2 21 and ∆m 2 32 ), E, L as well as the matter density. This implies that, even after taking into account the neutrino flux spectra, cross sections, energy resolution, detection efficiencies, and so on, which depend on neutrino energy, and after performing integrations over the true and reconstructed neutrino energies, the expected number of ν e (ν e ) appearance events, N νe (Nν e ), for a given experimental exposure (running time) have also the similar δ CP dependence as, N νe = n 0 + n c cos δ CP + n s sin δ CP , Nν e =n 0 +n c cos δ CP +n s sin δ CP , where n 0 , n c , n s ,n 0 ,n c andn s are some constants which depend not only on mixing parameters but also on experimental setups. Assuming that background (BG) events do not depend (or depend very weakly) on δ CP , the constant terms n 0 andn 0 in Eq. (2) can be divided into the signal contribution and BG one as n 0 = n sig 0 + n BG 0 andn 0 =n sig 0 +n BG 0 , as an approximation.
In Table A1, we provide the numerical values of these coefficients which can reproduce quite well the expected number of events shown in the plane spanned by N obs νe and N obs νe , often called bi-rate plots, found in the presentations by T2K [3] and NOvA [4] at Neutrino 2020 Conference, for their corresponding accumulated data (or exposures). We show in the left panels of Figures A1  and A2, respectively, for T2K and NOvA, the bi-rate plots which were reproduced by using the values given in Table A1. Our results are in excellent agreement with the ones shown by the collaborations [3,4]. 1 The χ 2 function for the appearance channel (AC), for a given LBνB experiment, T2K or NOvA, which is based on the total number of events, is simply defined as follows, for each MO, where N obs νe (N obs νe ) is the number of observed (or to be observed) ν e (ν e ) events, and N theo νe (N theō νe ) are the corresponding theoretically expected numbers (or prediction), and (4) Note that the number of events in Eq. (3) include also background events.  [4,5].
Using the number of events given in Eq. (2) with values of coefficients given in Table A1 we performed a fit to the data recently reported by T2K at Neutrino 2020 Conference [3] just varying sin 2 θ 23 and δ CP and could reproduce rather well the ∆χ 2 presented by T2K in the same conference mentioned above, as shown in the right panel of Figure A1. We have repeated the similar exercises also for NOvA and obtained the results, shown in the right panel of Figure A2, which are reasonably in agreement with what was presented by NOvA at at Neutrino 2020 Conference [4]. In the case of NOvA the agreement is slightly worse as compared to the case of T2K. We believe that this is because, for the results shown in Figure A2, unlike the case of T2K, we did not take into account the θ 23 constraint by NOvA (or we have set χ pull in Eq. (4) equals to zero) as this information was not reported in [4].
We note that for this part of our analysis, we considered only the dependence of sin 2 θ 23 and δ CP and ignore the uncertainties of all the other mixing parameters as we are computing the number of events in an approximated way, as described above, by taking into account only the variation due to sin 2 θ 23 and δ CP with all the other parameters fixed (separately by T2K [3] and NOvA [4] collaborations) to some values which are close to the values given in Table 1 of [1].
In particular, we neglected the uncertainty of ∆m 2 32 in the LBνB AC part analysis when it is combined with JUNO plus LBνB DC part analysis to obtain our final boosted MO sensitivities. Strictly speaking, ∆m 2 32 must be varied simultaneously (in a synchronised way) in the χ 2 defined in Eq. (5) when it is combined with the χ 2 defined in Eq. (14). However, in our analysis, we simply add ∆χ 2 obtained from our simplified LBνB AC simulation which ignored ∆m 2 32 uncertainty, to the JUNO's boosted χ 2 (described in detail in the Appendix C). This can be justified by considering that a variation of ∆m 2 32 of about ∼ 1% imply only a similar magnitude of variations in the appearance oscillation probabilities, which would be significantly smaller than the statistical uncertainties of LBνB-II AC mode, which are expected to reach at most the level or ∼5% or larger even in our future projections for T2K and NOvA.
For the MO resolution sensitivity shown in Figure 2 of [1] and used for our analysis throughout this work, we define the ∆χ 2 (labeled as ∆χ 2 AC LBνB ), as where +(-) sign corresponds to the case where the true MO is normal (inverted), and χ 2 AC LBνB is computed as defined in Eq. (3) but with N obs νe(νe) replaced by the theoretically expected ones for given values of assumed true values of θ 23 and δ CP . In practice, since we do not consider the effect of fluctuation for this part of our analysis, χ 2 AC LBνB min = 0 by construction for true MO. We 2   note that when T2K and NOvA are combined, some enhancement of sensitivities in the positive (negative) δ CP region for NMO (IMO) occur (see light green curves in Figure 6 of [1]. This is because that in these δ CP ranges, T2K and NOvA data can not be simultaneously fitted very well by using the common δ CP for the wrong MO, leading to an increase of ∆χ 2 . For simplicity, for our future projection, we simply increase by a factor of 3 both T2K and NOvA exposures, to the coefficients given in Table A1 for both ν andν channels. This corresponds approximately to 8.0 (6.4)×10 21 POT for T2K ν (ν) mode and to 4.1 (3.8)×10 21 POT for NOvA ν (ν) mode, to reflect roughly the currently considered final exposures for T2K [6] ( 10 × 10 21 POT in total for ν andν) and NOvA [5] ( 3.2 × 10 21 POT each for ν andν). This approach implies that our calculation does not consider future unknown optimisations on the ν (ν) mode running.

B. LBνB Disappearance MO Sensitivity
In the upper panel of Figure A3, we show the 4 curves of survival oscillation probabilities, P (ν µ → ν µ ) and P (ν µ →ν µ ) for NMO and IMO, which were obtained by using the best fitted parameters in NuFit5.0 given in Table 1 of [1] for the baseline corresponds to NOvA(L = 810 km) and with the matter density of ρ = 2.8 g/cm 3 . The NMO and IMO cases are shown, respectively, by blue and red colours whereas the cases for ν andν are shown, respectively, by solid thin and dashed thick curves. We observe that all of these 4 curves coincide very well with each other, so differences are very small. In the lower panel of the same Figure A3, we show the differences of these curves, between ν andν channels for both NMO and IMO, as well as between NMO and IMO for both ν andν, as indicated in the legend. We observe that the differences of these oscillation probabilities are ≤1% for the energy range relevant for NOvA.
Two points can be highlighted. First, the fact that the differences between neutrino and anti-neutrino are quite small implies that the matter effects are very small in these channels, hence determining MO by using matter effects based only on LBνB DC would be almost impossible. And second, the fact that the curves for NMO and IMO agree very well implies that the absolute values of the effective mass squared differences, called ∆m 2 µµ , defined in Eq. (11) in Appendix D, which correspond to NMO and IMO cases, should be similar. Indeed, by using the values given in Table 1 of [1], we obtain ∆m 2 µµ = 2.422(−2.431) × 10 −3 eV 2 for NMO (IMO) exhibiting a small ∼ 0.4% difference. In other words, for each channel, ν andν, there are two degenerate solu-tions, one corresponds to NMO and the other, to IMO, which give in practice the same survival probabilities. We stress that this degeneracy can not be resolved by considering LBνB experiment with DC alone.

C. Analytic Understanding of Synergy between JUNO and LBνB based experiments
In this section, we shall detail the relation between true and false ∆m 2 32 solutions in the case of JUNO and LBνB, as they are different. This difference is indeed exploited as the main numerical quantification behind the ∆χ 2 BOOST term which was schematically illustrated in Figure 3  Theν e →ν e survival probability in vacuum can be expressed as [7] Pν e→νe = 1 − c 4 13 sin 2 2θ 12 sin 2 ∆ 21 − where the notation c ij ≡ cos θ ij and s ij ≡ sin θ ij is used, and ∆ ij ≡ ∆m 2 ij L/4E, L and E are, respectively, the baseline and the neutrino energy, and the effective mass squared difference ∆m 2 ee is given by [8] ∆m 2 ee ≡ c 2 12 ∆m 2 31 + s 2 12 ∆m 2 32 = ∆m 2 32 + c 2 12 ∆m 2 21 , (7) and φ is given by tan φ = c 2 12 sin(2s 2 12 ∆ 21 ) − s 2 12 sin(2c 2 12 ∆ 21 ) c 2 12 sin(2s 2 12 ∆ 21 ) + s 2 12 sin(2c 2 12 ∆ 21 ) , where the approximated value of δm 2 φ can be estimated by choosing the average representative energy of reactor neutrinos (∼4 MeV) as  Table 1 of [1] are shown for NMO (solid and dashed red curves) and IMO (solid and dashed blue curves) as a function of neutrino energy, as indicated in the legend. In the lower panel, the differences of these probabilities are shown in percent.
We found that for a given assumed true value of ∆m 2 32 = 2.411×10 −3 eV 2 (corresponding to NMO), we can reproduce very well the false value of ∆m 2 32 = −2.53 × 10 −3 eV 2 (corresponding to IMO) obtained by a χ 2 fit if we use E = 4.4 MeV in Eqs. (9) and (10). The relation between true and false ∆m 2 32 for JUNO is illustrated by the vertical black dashed and black solid lines in Figure A4 (b) and (d).

C.2 LBνB Relation between True-False ∆m 2 32
For LBνB experiments like T2K and NOvA the L/E are such that |∆ 31 | ∼ |∆ 32 | ∼ π/2. From the disappearance channels ν µ → ν µ andν µ →ν µ , it is possible to measure precisely the effective mass squared difference ∆m 2 µµ whose value is independent of the MO. In terms of fundamental mixing and oscillation parameters, ∆m 2 µµ can be expressed, with very good approximation, as [ where in the last line of the above equation, some simplifications were considered based on the fact that best fitted values of sin 2 θ 13 and sin 2 θ 23 in recent global analysis [9] are similar for both MO solutions. By using the relation given in Eq. (13), for a given assumed true value of ∆m 2 32 (common for all experiments) we obtain the yellow colour bands shown in Figure A4 (b) and (d).

C.3 Boosting Synergy Estimation
The extra synergy for MO determination sensitivity by combining JUNO and LBνB DC can be achieved thanks to the mismatch (or disagreement) of the fitted ∆m 2 32 values for the wrong MO solutions between these two types of experiments. For the correct MO, ∆m 2 32 values measured by different experiments should agree with each other within the experimental uncertainties. But for those values which correspond to the wrong MO do not agree. The difference can be quantified and used to enhance the sensitivity.
Following the procedure described in [10,11], we include to the JUNO analysis the external information on the external information on ∆m 2 32 from LBνB with an additional pull term as where χ 2 JUNO implies the χ 2 function for JUNO alone computed in a similar fashion as in [11], σ(∆m 2 32 ) LBνB is the experimental uncertainty on ∆m 2 32 achieved by LBνB based experiments. As typical values in this paper, we consider 3 cases σ(∆m 2 32 ) LBνB = 1, 0.75 and 0.5%.
In order to take into account the possible fluctuation of the central values of the measured ∆m 2 32LBνB we define the extra boosting ∆χ 2 due to the synergy of JUNO and LBνB based experiments as the difference of χ 2 defined in Eq. (14) for normal and inverted MO as, where +(-) sign corresponds to the case where the true MO is normal (inverted). Note that in our simplified phenomenological approach (based on the future simulated JUNO data), for the case with no fluctuation, by construction, χ 2 NMO (IMO) = 0 for NMO (IMO). Let us try to see how the boosting will be realized by applying our discussion to JUNO and T2K for illustration. In Figure A4 for the cases where the true MO is normal in the panel (a) and inverted in the panel (c) we show by the solid (dashed) black curve ∆χ 2 for JUNO alone case for true (false) MO. The difference of ∆χ 2 between true and false MO is 9 if only JUNO is considered implying that the false MO (indicated by the dashed curves) can be rejected at 3σ. On the other hand, let us assume the case where T2K can determine |∆m 2 32 | with 1% uncertainty, and the corresponding ∆χ 2 curves are given by the solid (dashed) blue curves for true (false) MO in the same plots, rejecting the wrong MO only at 2σ by T2K alone. If we combine JUNO and T2K following the procedure described in this section, the resulting ∆χ 2 are given by the solid (dashed) red curves for true (false) MO, rejecting the wrong MO with more than 4σ for both NMO and IMO.
The large (∼ 10) increase of the combined ∆χ 2 for the wrong MO fit comes from the mismatch of the false |∆m 2 32 | values between JUNO (black dashed line) and T2K (yellow colour bands) shown in the panels (b) and (d) of Figure A4. This is nothing the boosting effect, which can be analytically understood and quantified as follows.
Suppose that we try to perform a χ 2 fit assuming the wrong MO. Let us first assume that σ(∆m 2 32 ) JUNO σ(∆m 2 32 ) LBνB and no fluctuation for simplicity (i.e. where the numbers in the last line were estimated for σ(∆m 2 32 ) LBνB = 1%. The case where δ CP = ±π/2 and ∆χ 2 boost ∼ 9 can be directly compared with more precise results shown in Figure 4(a) of [1], see the blue solid curve at δ true CP = ±π/2 which gives ∆χ 2 boost ∼ 8 which is in rough agreement. The expression in Eq. (16) is in agreement with the one given in Eq. (18) of [10] apart from the term δm 2 φ which is not so large.

D. Full 3 ν versus effective 2ν formulation
In the previous discussions found in [10,11], in order to demonstrate the boosting synergy effect between JUNO and LBνB experiments, the effective mass squared differences ∆m 2 ee and ∆m 2 µµ , defined respectively, in Eqs. (7) and (11) originally found in [8] were used. While we used these parameters in some intermediate steps of our computations, as described in Appendix C, we did not use these parameters explicitly in our combined χ 2 describing the extra synergy between JUNO and LBνB (DC) based experiments defined in Eq. (14), as well as in the final sensitivity plots presented in this paper. The main advantage of using these effective mass squared differences is that no a priori assumptions have to be made about any other parameters not accessible by JUNO, in particular, CP phase (δ CP ), whereas by using them, one must specify explicitly the δ CP value as done in Ref. [10].
In order to check the consistency between our work and previous studies, we have explicitly verified that the results do not depend on the parameters used in   true is fixed to the NuFit5.0 best value shown, respectively, in panels (a) and (c) for NMO and IMO. The extra gain in ∆χ 2 BOOST discrimination numerically originates from the fact that the true |∆m 2 32 | solutions should match between JUNO (solid black vertical line) and LBνB (solid blue vertical circle); hence the false solutions (dashed vertical lines) must differ. Panels (b) and (d) illustrate this origin. The relation between true-false ∆m 2 32 solutions is different and complementarity for JUNO and LBνB experiments. The difference is large (≈ 1.5×δm 2 21 ) for JUNO. Instead, LBνB exhibits a smaller difference that modulates with δ CP . So, the relative difference between ∆m 2 32 false JUNO and ∆m 2 32 false LBνB is maximal (minimal) for the δ CP -conserving ±π (0) value. Hence, ∆χ 2 BOOST depends on δ CP , and ambiguity arises (yellow band) from the a priori different values of δ CP for the true or false solutions. The T2K data (red points) contrasts the precision on |∆m 2 32 | now [12] compared to needed scenarios ≤1.0% scenario (blue points and parabolas). The precision of each contribution is indicated by the parabolas' width, where JUNO is fixed to the nominal value [11]. µµ . Expected MO sensitivity to be obtained by JUNO with external information of ∆m 2 µµ coming from LBνB experiments following the procedure described in [10,11], are shown as a function of the precision of ∆m 2 µµ for cos δ CP = ±1 and 0. This plot is similar to Figure 7 in [10], once upgraded to the latest global data inputs [9]. We observe that they are consistent with each other, if the curves for the δ CP values of 0 • (blue) and 180 • (red) were interchanged, as a result of a typo in the legend of [10]. the analysis and in the presentation of the final results, provided that that comparisons are done properly. In Figure A5, we show ∆χ 2 (JUNO⊕LBνB-DC) computed by using explicitly ∆m µµ (instead of using ∆m 2 32 ) in our χ 2 analysis as done in [10,11], as a function of the precision of ∆m 2 µµ . There is general good agreement with the result shown in Figure 7 of [10], if δ CP curves for 0 • and 180 • were interchanged, as described in Figure A5.