Combined sensitivity of JUNO and KM3NeT/ORCA to the neutrino mass ordering

This article presents the potential of a combined analysis of the JUNO and KM3NeT/ORCA experiments to determine the neutrino mass ordering. This combination is particularly interesting as it significantly boosts the potential of either detector, beyond simply adding their neutrino mass ordering sensitivities, by removing a degeneracy in the determination of $\Delta m_{31}^2$ between the two experiments when assuming the wrong ordering. The study is based on the latest projected performances for JUNO, and on simulation tools using a full Monte Carlo approach to the KM3NeT/ORCA response with a careful assessment of its energy systematics. From this analysis, a $5\sigma$ determination of the neutrino mass ordering is expected after 6 years of joint data taking for any value of the oscillation parameters. This sensitivity would be achieved after only 2 years of joint data taking assuming the current global best-fit values for those parameters for normal ordering.

The discovery of neutrino flavor oscillations, implying that neutrinos are massive particles, is so far one of the few observational hints towards physics beyond the Standard Model.As such, it has potentially far-reaching implications in many aspects of fundamental physics and cosmology, from the matter-antimatter asymmetry in the Universe to the naturalness problem of elementary particles (see e.g.Refs.[1][2][3]).Since the first conclusive observations of neutrino oscillations at the turn of the century [4][5][6], a variety of experiments targeting solar, reactor, atmospheric and accelerator neutrinos have achieved an increasingly precise determination of the parameters of the neutrino flavor mixing matrix [7][8][9][10].Despite this tremendous progress, some fundamental properties of neutrinos are yet to be determined, such as their absolute masses, whether they are Majorana particles and therefore are their own anti-particle, the existence and strength of CPviolation in the neutrino sector, and the ordering of the masses (m 1 , m 2 and m 3 ) of the neutrino mass eigenstates (respectively, ν 1 , ν 2 and ν 3 ): either normal ordering (NO, m 1 < m 2 < m 3 ) or inverted ordering (IO, m 3 < m 1 < m 2 ).This last question is a prime experimental goal because its determination would have direct consequences on, e.g., the measurement of leptonic CP violation in future long baseline experiments [11] and the interpretation of results from planned experiments searching for neutrino-less double-beta decay to establish the Dirac vs. Majorana nature of neutrinos [12].The measurement of the neutrino mass ordering (NMO) is on the agenda of several ongoing neutrino experiments in the GeV energy domain that probe long-baseline ν µν e oscillations in Earth matter.Such experiments are sensitive to the atmospheric mass splitting ∆m 2  31 [13].However, none of these experiments alone, either accelerator-based (such as T2K [14] or NOνA [15]) or using atmospheric neutrinos (such as Super-Kamiokande [16] or IceCube [17]), has the capability to unambiguously resolve the NMO (in other words, the sign of ∆m 2  31 ) within the next few years.Even combining all available data, including those from reactor experiments, into global neutrino oscillation fits, has so far yielded only a mild preference for normal ordering, which has faded away again since the inclusion of the latest results of T2K [18] and NOνA [19].
Considering that a high-confidence (> 5σ) determination of the NMO with the next-generation accelerator experiments DUNE [20], T2HK [21] and T2HKK [22] is only envisaged for 2030 or beyond, alternative paths to the NMO measurement are being pursued on a shorter timescale.
JUNO [23][24][25] and KM3NeT/ORCA [26] are the two next-generation neutrino detectors aiming at addressing the NMO measurement within this decade.ORCA, the low-energy branch of the KM3NeT network of water Cherenkov neutrino telescopes, will determine the NMO by probing Earth matter effects on the atmospheric neutrino oscillations in the GeV energy range.JUNO is a medium-baseline (∼ 53 km) reactor neutrino experiment that is sensitive to the NMO through the interplay between the fast oscillations driven by ∆m 2  31 and ∆m 2  32 in the νe disappearance channel, where matter effects play only a small role [27].As reactor νe disappearance measurement is not affected by CP violation [23], the JUNO measurement will be independent of the unknown CP violating phase δ CP .Both detectors are currently under construction and target completion within the first half of this decade.The JUNO detector is planned to be completed in 2022, while ORCA is foreseen to be deployed incrementally until 2025, with 6 out of the total 115 detection lines already installed and taking data [28].
Combining the data from the two experiments is not only motivated by their almost simultaneous timelines, but also by the gain in sensitivity that arises from the complementarity of their approaches to the measurement of the NMO.This boost essentially comes from the expected tension between ORCA and JUNO in the best fit of |∆m 2  31 | from the νe and ν µ disappearance channels when the assumed ordering is wrong.Provided that the measurement uncertainties in each experiment are small enough, the wrong mass ordering can be excluded with a high confidence level from the combination of the two datasets, even if the intrinsic NMO sensitivity of each experiment would not reach that level.This effect was first pointed out in relation with accelerator neutrino experiments [29,30], then reassessed in the context of the combination of a reactor experiment (Daya Bay II, now evolved into JUNO) and an atmospheric neutrino experiment (PINGU [31], a proposed low-energy extension of the IceCube neutrino telescope), showing that a strong boost in NMO sensitivity can indeed be reached with a combined fit [32].The potential of this method was further explored in a combined study with JUNO and PINGU using detailed simulation tools for both experiments [33], leading to the same conclusion.
In this paper, a complete study of the combined sensitivity of JUNO and ORCA to the NMO is presented, based on the same theoretical approach and using up-to-date detector configurations and expected performances.The main features and detection principles of the JUNO experiment are described in Sec. 2, along with the standalone χ 2 analysis used to determine the JUNO-only NMO sensitivity.The same is done for ORCA in Sec. 3, based on the updated detector configuration, and simulation and reconstruction tools used for the latest NMO sensitivity study following Ref.[34].In this case, special attention is paid to the treatment of systematic uncertainties, in particular with the introduction of a systematic error on the measured energy scale at the detector level (not considered in Ref. [33] for PINGU).That systematic effect is shown to degrade the precision of the ∆m 2  31 measurement of ORCA alone, thereby affecting the sensitivity of the combined study.
The JUNO/ORCA combined χ 2 analysis is presented in Sec. 4 for the baseline reactor configuration of JUNO, and the most realistic systematics treatment adopted for ORCA.The enhanced sensitivity achieved with the combined χ 2 analysis over the simple sum of individual χ 2 is also demonstrated.Sec. 5 presents further sensitivity studies exploring the impact of the energy resolution in JUNO and changing the number of reactor cores available towards a more optimistic scenario.The stability of the combined performance versus the true value of θ 23 and ∆m 2  31 is also addressed.The conclusions drawn from these results are presented in Sec. 6.

Jiangmen Underground Neutrino Observatory (JUNO)
The Jiangmen Underground Neutrino Observatory [23][24][25] (JUNO) is a multipurpose experiment being built in the south of China.Among its goals is the determination of the NMO via the precise measurement of νe from the Yangjiang and the Taishan Nuclear Power Plants (NPP) located 53 km away from the detector.The JUNO detector is divided into 3 parts: the Central Detector, the Water Cherenkov Detector and the Top Tracker.The Central Detector is composed of 20 kton of liquid scintillator placed in a 35.4 m diameter acrylic sphere.Around this acrylic sphere, about 18k 20" and 26k 3" photomultiplier tubes (PMT) monitor the liquid scintillator volume to detect neutrino interactions occurring inside, in particular the Inverse Beta Decay (IBD) interactions produced by νe from the NPPs.The IBD interactions are detected in the JUNO Central Detector via the prompt detection of the scintillation light of the positron produced in the interaction and its subsequent annihilation, along with the delayed detection of a 2.2 MeV gamma-ray produced in the neutron capture on hydrogen and subsequent de-excitation of the deuteron.Due to the kinematics of the IBD, most of the available energy of the incident νe is transferred to the positron.Therefore, to do a precise measurement of neutrino oscillations, a good energy resolution to measure the visible energy of the prompt signal is critical, as will be discussed later.The measured visible energy is smaller than the incident νe energy by about 0.8 MeV, due to the mass difference between the initial and final particles (−1.8 MeV) and to the light emitted in the positron annihilation (+1.0 MeV).The acrylic sphere is placed in the center of the Water Cherenkov Detector, a cylindrical ultra-pure water pool (44 m height, 43.5 m diameter), which serves to shield the Central Detector from external radioactivity and to provide a veto for atmospheric muons and for muon-induced background such as cosmogenic nuclei and fast neutrons.The Water Cherenkov Detector and the Top Tracker, located on top of it to precisely track atmospheric muons, compose the Veto System of JUNO.
In addition to the JUNO detector described above, the project also includes the Taishan Antineutrino Observatory (JUNO-TAO) detector [35].This detector will be installed at a distance of 30 m from one of the Taishan's reactors to determine the reactor νe spectrum with a better energy resolution than JUNO, effectively reducing the impact of possible unknown substructures in the reactor neutrino spectra [36] on the measurement of neutrino oscillations.
The precise distance of JUNO to each reactor core of the Yangjiang and Taishan NPPs, provided in Ref. [23], is used in this study rather than just the distance to the NPP complex.In addition to them, the NPPs of Daya-Bay at 215 km and Huizhou at 265 km will also contribute to the total number of detected reactor neutrinos.However, given the much larger distance, the oscillation pattern will not be the same and these neutrinos are part of JUNO's intrinsic background.The Yangjiang NPP is already fully operational, with 6 reactor cores and a total of 17.4 GW of thermal power, as is the Daya-Bay NPP, with a similar total thermal power.The Taishan NPP has already 2 reactor cores operational out of the 4 initially foreseen, totaling a thermal power of 9.2 GW.At present, it is unknown if the remaining 2 reactor cores, which would bring another additional 9.2 GW of thermal power, will be built.The Huizhou NPP is under construction and is expected to be ready by about 2025 [25] with 17.4 GW thermal power.

Modeling JUNO for this study
For this study, the performance of JUNO closely follows that provided in Ref. [23].In particular, a 73% IBD detection efficiency and an energy resolution of 3%/ E/MeV are assumed in the Central Detector which contains 1.5 × 10 33 target protons.Given that the energy resolution is critical for the JUNO sensitivity, the impact of a ±0.5%/E/MeV change in energy resolution is discussed in Sec.5.1.As in Ref. [23], the nominal running time of the experiment is considered to be 1000 effective days every 3 years.
The observed νe spectrum in JUNO will be produced by the interplay between the spectrum produced by the NPPs, the IBD cross-section [37], and the neutrino oscillations which are to be measured.To determine the spectrum produced by the NPPs, the ILL1 νe spectra [38][39][40] are used given that the flux normalization is in better agreement with previous data [41].The fine structure of the spectrum will be precisely determined independently using JUNO-TAO, as mentioned beforehand.Therefore, it is not explicitly included in the spectrum shape.The spectrum is calculated assuming a fission fragments content of 235 U : 239 Pu : 238 U : 241 Pu = 0.564 : 0.304 : 0.076 : 0.056 which is similar to the one from Daya-Bay [42], and using the fission energies for each of these isotopes from Ref. [43].
In addition to reactor neutrinos, the IBD event selection contains some background events.The backgrounds considered in this analysis are taken from Ref. [23], in terms of their rate, shape, and uncertainties.The three dominant components of this background are cosmogenic events, geo-neutrinos, and accidental coincidences mainly from radioactive background, with expected rates of about 1.6, 1.1, and 0.9 events per day, respectively.As a reference, the Daya-Bay and Huizhou NPPs are expected to yield a total of 4.6 events per day in JUNO while the Taishan and Yangjiang NPPs are expected to produce a total of 54.3 events per day, assuming the normal ordering world best-fit [44] oscillation parameters, considering all 4 Taishan NPP reactors operational.
A notable difference from Ref. [23] is our conservative choice to consider, as baseline, only the 2 existing reactors in the Taishan NPP (referred to as "JUNO 8 cores" hereafter) rather than the foreseen 4 reactors (referred to as "JUNO 10 cores" hereafter).This reduces the total number of expected signal neutrinos by about 25% in this study in comparison to Ref. [23].For completeness, the JUNO official baseline with a total of 4 Taishan reactors is also considered and discussed in Sec.5.1.
Although not yet completed, the Huizhou NPP is considered to be active for the whole duration of JUNO in both cases, even if it is possible that JUNO will start before its completion.This is again a conservative assumption given that the Huizhou NPP νe are an intrinsic background to the neutrino oscillation measurements in JUNO.
The expected event distribution as a function of the visible energy is shown in Fig. 1, for 6 years of data with 8 cores and assuming the best-fit oscillation parameters from Ref. [44] for normal ordering.The events corresponding to the expected remaining non-reactor neutrino background are also highlighted in the plot.They are concentrated in the lower energy part of the measured spectrum where the energy resolution of JUNO is not sufficient to see the rapid oscillation pattern.

Sensitivity analysis
In JUNO, the measurement of the neutrino oscillation parameters and, in particular, of the neutrino mass ordering, is done by fitting the measured positron energy spectrum.As shown in Fig. 1, this spectrum exhibits two notable features.The first is a slow oscillatory behavior due to θ 12 and ∆m 2  21 , which causes a large deficit in the number of events in the whole energy range shown in Fig. 1 and that has a minimum at about 2.2 MeV.The second is a rapid oscillatory behavior due to θ 13 and ∆m 2  31 that starts being visible in the figure at about 2 MeV, and for which several periods are shown.If inverted ordering is assumed rather than normal ordering in Fig. 1, the position of the rapid oscillation maxima and minima would change, as the oscillation frequencies producing the pattern would be different [23].This is the signature that JUNO will use to measure the NMO.The sensitivity of JUNO to determine the NMO is calculated using the χ 2 difference between the data being fitted under the two ordering hypotheses.This strategy is also applied for the combined analysis.While this procedure will be further detailed in Sec. 4, it is useful to describe already at this point the χ 2 function used in JUNO.
It is also worth pointing out here that in this study no statistical fluctuations are added to any of the simulated samples, and the expected statistical uncertainties at various timelines are taken into account in the computed χ 2 values.This approach is commonly referred to as "Asimov" approach [45], and the "Asimov" term will be used in this paper to highlight that this approximation is being used.
For this analysis the measured JUNO visible energy spectrum is divided in n = 207 bins of 0.03 MeV between 1.00 MeV and 7.21 MeV.The JUNO χ 2 function used in this analysis has the following form: In this expression, ∆ is a 1 × n matrix whose content is the difference between the observed and expected rates.∆ is defined as ∆ = D − (S + B), where D, S, and B correspond, respectively, to the data, the signal prediction for a given set of oscillation parameters, and the expected background.D is given by the Asimov sample for the true value of the oscillation parameters, while S and B correspond to the expected signal given the test oscillation parameters and the background.∆ T is the transpose matrix of ∆.
The matrix M whose inverse is present in Eq. ( 1) is a n × n covariance matrix.This matrix is calculated as M = M stat + M S + M B .M stat corresponds to the statistical uncertainty in each bin from the total expected number of events (S + B).M S and M B correspond, respectively, to the covariance matrices of the signal and background as described in Ref. [23].
This JUNO-only analysis was validated by comparing the NMO sensitivity with previous results in Refs.[23,33].Cross-checks have been performed using the same set of oscillation parameters and reactor cores as in those references, showing agreement within 0.1-0.5 χ 2 units.

Oscillation Research with Cosmics in the Abyss (ORCA)
The KM3NeT Collaboration is currently building a set of next-generation water Cherenkov telescopes in the depths of the Mediterranean Sea [26].Two tridimensional arrays of PMTs will be deployed: ARCA and ORCA (for Astroparticle and Oscillation Research with Cosmics in the Abyss, respectively).ARCA is a gigaton-scale detector which will mainly search for astrophysical neutrinos in the TeV-PeV energy range.ORCA, subject of this study, is a denser and smaller array (Mton-scale) optimized for oscillation physics with atmospheric neutrinos at energies above 1 GeV.
The ORCA detector will contain 115 detection units, each of them being a vertical line about 200 m long, anchored to the seabed and supporting 18 digital optical modules.These modules are glass spheres that contain 31 PMTs and related electronics.The array of detection units is arranged in a cylindrical shape with an average radius of 115 m and an average distance between the lines of 20 m.On each detection unit, the vertical spacing between the optical modules is around 9 m.The total instrumented volume of ORCA covers about 6.7 × 10 6 m 3 , corresponding to 7.0 Mtons of seawater.
The ORCA PMTs detect the Cherenkov light induced by charged particles originating in neutrino interactions in and around the detector.Such detectable interactions occur mainly through charged current (CC) and neutral current (NC) deep inelastic scattering processes off nucleons in water molecules [46].The pattern and timing of the digitalized PMT output, or hits, recorded by the digital optical modules are used to identify neutrino events and reconstruct their energy and angular direction.The topologies of neutrino-induced events in the energy range of interest for ORCA can be separated into two broad classes.If the final state includes a sufficiently energetic muon, a track-like signature will be produced.This is the case for CC interactions of ν µ /ν µ , and ν τ /ν τ with subsequent muonic decay of the tau lepton.Shower-like events correspond to all other interaction channels, where only hadronic and electromagnetic showers are produced.This includes CC interactions of ν e /ν e and ν τ /ν τ with non-muonic decays, as well as NC interactions of all flavors.

Modeling the ORCA detector for this study
The analysis presented here is based on a detailed Monte Carlo (MC) simulation of the ORCA detector response to atmospheric neutrinos.The generation of neutrino interactions in seawater in the energy range 1-100 GeV is performed using gSeaGen [47], a GENIE [48]-based software developed within the KM3NeT Collaboration.All secondary particles are tracked with the package KM3Sim [49] based on GEANT4 [50] which simulates and propagates the photons induced by the Cherenkov effect, accounting also for light absorption and scattering, and records the signals reaching the PMTs.The optical background due to radioactive decays of 40 K naturally present in seawater is simulated by adding uncorrelated (single-PMT) and correlated (inter-PMT) random noise hits, based on the rates measured with the first deployed detection units [51].The background of atmospheric muons is also simulated.The PMT response, readout, and triggering are simulated using KM3NeT custom software packages.The resulting trigger rate is about 54 kHz for noise events, 50 kHz for atmospheric muons and about 8 mHz for atmospheric neutrinos.The total simulated sample includes more than 15 years of atmospheric neutrinos, 1.4 days of noise events, and 14 days of atmospheric muons, which proves sufficient to probe the background contamination at a percent level [34].
The MC neutrino sample and event selection adopted in this study are identical to those used in the latest ORCA NMO sensitivity study and are extensively described in Ref. [34].From the detected signals, the energy and direction of the events are reconstructed using dedicated algorithms developed for shower-like [52] and track-like [53] event topologies.
A set of preselection cuts are applied, requiring that events are well reconstructed, with an up-going direction (corresponding to a reconstructed zenith angle θ reco > 90 • ), and obey certain criteria of containment in the detector instrumented volume.Events that pass the preselection cuts are then processed by a classification algorithm based on random decision forests for the determination of event topologies (or particle identification -PID) and background rejection.Based on a single score scale η from 0 to 1 provided by the classifier, the sample is then divided into 3 classes of events, also called PID classes: tracks (0.7 < η ≤ 1), intermediate events (0.3 < η ≤ 0.7), and showers (0 ≤ η ≤ 0.3).The two extreme PID classes provide a higher purity level for the genuine track-like and shower-like events.As discussed in Ref. [34], the selection and background suppression cuts are sufficient to reduce pure noise events to a negligible rate and result in an atmospheric muon background contamination of only 3%.The impact of such contamination is expected to be insignificant and atmospheric muons are not included in this study.
The analysis relies on the computation of the expected energy and zenith angle (E, θ) distributions of atmospheric neutrino events for each PID class.Such distributions are obtained with the SWIM package [54], a KM3NeT analysis framework developed for calculating event distributions for a given hypothesis using a full MC approach to model the detector response.
The incoming atmospheric neutrino flux is taken from Ref. [55], for the Gran Sasso site without mountain over the detector, assuming minimum solar activity.The probabilities of neutrino flavor transitions along their path through the Earth are computed with the software OscProb [56], using a radial model of the Earth with 42 concentric shells of constant electron density, for which mass density values are fixed and follow the Preliminary Reference Earth Model [57].The rate of events interacting around the detector is computed for each interaction type ν x ∈ { ν NC} using neutrino-nucleon cross-sections weighted for water molecules as obtained with GENIE.
In order to obtain the expected event distribution as observed by the detector, the SWIM package uses a binned detector response matrix built from the MC sample that maps the events generated with interaction type ν x and true variables (E true , θ true ) into the corresponding reconstructed variables (E reco , θ reco ) and PID class (track, intermediate and shower).Given that there are 8 different interaction types and 3 different PID bins, the global response matrix is a collection of 24 4-dimensional matrices used for the transformation (E true , θ true ) −→ (E reco , θ reco ).These matrices are built using MC-generated events and the outcome of their processing through the reconstruction and classification algorithms, so that the ensemble of matrices account for detection and reconstruction efficiencies, misidentification probabilities and errors on reconstructed variables (including all correlations).This approach is different from the one in Ref. [34] which uses parametrized response functions obtained from the MC distributions.
While this full MC method ensures that all the information on the detector response, including potential correlations between parameters, is taken into account, its accuracy depends on the size of the MC sample.To account for statistical fluctuations in the MC production that could affect the response functions used to build the matrix, the Beeston-Barlow light method [58] has been adopted, as described in the next section.
Fig. 2 depicts the expected neutrino event distributions calculated with SWIM for all PID classes, binned in reconstructed variables (E reco , cos θ reco ) for 6 years of data taking with ORCA, assuming normal ordering and oscillation parameters from Ref. [44].In this study, 20 linear bins are chosen for the reconstructed cosine of the zenith angle while the reconstructed energy is binned logarithmically with 20 bins in the range of 2-80 GeV.

Sensitivity analysis
The ORCA analysis uses the 2D distribution of expected neutrino events in each PID class as a function of the reconstructed neutrino energy and the cosine of the zenith angle.The SWIM package performs the computation and minimization of the test statistic chosen as the log-likelihood ratio between a hypothetical model and data, which for this analysis is an Asimov dataset with an assumed true hypothesis.
The ORCA χ 2 is built as follows: The first sum is the Poisson likelihood of the data d l given the expectation µ l at bin l, where the latter also depends on the nuisance parameters p j .The β parameters are introduced based on the "Beeston-Barlow light method" [58] to account for fluctuations due to finite MC statistics.
The second sum accounts for the Gaussian priors on nuisance parameters p j with mean values p 0 j and variances σ j .The third sum represents the Gaussian priors to the β parameters that are expected to be normally distributed.
The fluctuations β k are assumed to be bin-to-bin uncorrelated and independent of the model parameters.Given this assumption, the values of β k are solved analytically as the solution of ∂χ 2 /∂β k = 0.In addition, the variances k of β k can also be evaluated with a probabilistic model which describes the calculation of the response matrix as a single binomial process [59,60].Finally, both β and its variance can be estimated analytically and used directly in the calculation of the χ 2 without any requirement for additional minimization.The full description of this procedure can also be found in Ref. [54].This implementation results in a ∼ 0.2σ decrease in the sensitivity, reflecting a correction of the overestimation caused by the limited MC sample size.
Two different sets of systematic uncertainties are used in this study (see Tab. 1).The "baseline" scenario corresponds to the standard set of ORCA systematics adopted for oscillation analyses, similar to Ref. [34].Uncertainties related to the incident flux include the spectral index of the neutrino flux energy distribution (free without any constraints) and the flux skew systematics.These skew parameters are introduced to describe the uncertainties in the ratio of different event types, namely ν e /ν e , ν µ /ν µ , and (ν e + νe )/(ν µ + νµ ), while preserving in each case the total associated flux [54].They are constrained with the priors adapted from Ref. [61].A NC normalization systematic is implemented as a scaling factor for the number of NC events.To account for detector-related uncertainties, an energy scale systematic is introduced as a global shift of the neutrino true energy in all detector response functions.This implementation captures Table 1: Baseline (see Ref. [34]) and optimistic (see Ref. [33]) scenarios for the treatment of systematics considered in the ORCA analysis.The cross (×) indicates that the systematic is not included.

Parameter
Baseline  A second scenario is based on the study of Ref. [33], developed for the PINGU detector.That analysis does not apply normalization factors to the PID classes.It uses an overall scaling factor which represents a universal systematic uncertainty on all effective areas (or equivalently, on the combined ν + ν event rate).This effective area scaling, together with an energy scale uncertainty introduced at the flux level, are the only systematics introduced to account, e.g., for potential variations in the detection efficiency of the optical modules.Contrary to the baseline case, these systematics are not introduced at the detector response level and are therefore considered as more optimistic in the rest of our study.The difference between the two approaches is further illustrated in Fig 3, showing the implementation of each set of systematic uncertainties into the workflow of the SWIM framework.The "baseline" systematic set is believed to be more accurate to describe the uncertainties in the ORCA detector.It is therefore used for all presented results, unless when stated explicitly that the "optimistic" systematics are used for the sake of crosschecks and comparisons.
4 Combined JUNO and ORCA analysis

Combination strategy
As described in Secs. 2 and 3, the detectors involved in this combined analysis work in very different conditions.This is true both in terms of their detection techniques and backgrounds, and of the sources and energies of neutrinos relevant for each individual analysis.The only common parameters to both experiments are the neutrino oscillation parameters, which are the core of the present study.
It is also important to note at this point that not all parameters used to describe standard neutrino oscillations have an impact on the results of this analysis.On one hand, JUNO has no sensitivity to either θ 23 or δ CP as this experiment measures νe → νe oscillations which do not depend on those parameters [23].On the other hand, ORCA has negligible sensitivity to θ 12 and ∆m 2  21 as the measured ν µ + νµ oscillations happen at a much smaller L/E than the one required for the development of oscillations with a frequency given by ∆m 2  21 [26].The four oscillation parameters that impact a single experiment are accounted for implicitly in the χ 2 function computation for JUNO and ORCA following the prescription outlined below, while the remaining two oscillation parameters, ∆m 2  31 and θ 13 , have to be considered explicitly in the joint analysis.
In JUNO, for every value of ∆m 2 31 and θ 13 , the χ 2 function is minimized using a grid with 61 uniformly spaced values of sin 2 θ 12 between 0.30225 and 0.31775.The value of ∆m 2  21 is kept fixed at its assumed true value given that JUNO will be able to determine this parameter quickly.Studies have shown that profiling this parameter or keeping it fixed would have an impact smaller than about 0.1 units of χ 2 , which is negligible in the joint analysis.
In ORCA, for every value of ∆m 2 31 and θ 13 , the θ 12 and ∆m 2 21 values are kept fixed to their assumed true values given that ORCA has little sensitivity to those parameters, while θ 23 and δ CP are minimized without constraints.This minimization is performed twice, with the initial value of θ 23 being located in a different octant for each minimization.Only the smallest value is kept as the global minimum of the χ 2 for ORCA.This is done to ensure that the minimizer is not trapped in a possible local minimum.
In order to combine the separate JUNO and ORCA analyses, their obtained χ 2 values at a fixed test value of ∆m 2  31 and θ 13 are calculated and summed.The true value of the oscillation parameters considered in this study are the best-fit values from Ref. [44] obtained "with SK data", unless it is explicitly stated otherwise.For added clarity, those parameters are explicitly shown in Tab. 2. Given that neither JUNO nor ORCA are as sensitive to θ 13 as current reactor neutrino experiments [62][63][64], a prior on that parameter was added to the combined χ 2 from Ref. [44].The full expression used is shown in Eq. ( 3) where ∆m 2  31 and θ 13 are the tested values For each set of true parameters studied, the combined χ 2 from Eq. ( 3) is calculated for each NMO in a 101 × 21 grid in the ∆m 2  31 , sin 2 θ 13 space, called the χ 2 map, centered around the assumed true values of the oscillation parameters and spanning uniformly a ±10% interval in ∆m 2  31 from the central value and a ±6% interval in sin 2 θ 13 from the central value.More explicitly, when assuming true normal ordering with the best-fit values from Ref. [44], the tested values of ∆m 2  31 in the grid will run from −2.78080 × 10 −3 eV 2 to −2.27520 × 10 −3 eV 2 and from 2.27520 × 10 −3 eV 2 to 2.78080 × 10 −3 eV 2 with step of 0.00506 × 10 −3 eV 2 , and those of sin 2 θ 13 in the grid will be from 0.0210278 to 0.0237122 with a step of 0.0001342.It is worth noting that when the true value of the oscillation parameters is changed, as in Sec.5.2, or when assuming inverted ordering, the grid described above is changed so that the central value of the grid corresponds to the true oscillation parameters.
Using the χ 2 map above, for each set of true oscillation parameters tested, the NMO sensitivity is determined by calculating the ) is the minimum value of χ 2 in the χ 2 map in the wrong (true) ordering region of the map.The ∆χ 2 is then converted into a median sensitivity S(σ) = ∆χ 2 [65].The same procedure is also used separately for ORCA and JUNO to obtain the corresponding non-combined sensitivities, computed for each experiment alone.
The ∆χ 2 notation is used above rather than ∆χ 2 , to highlight the fact that an Asimov approximation is being used in the entirety of this paper, therefore, the median sensitivity is always calculated.

Results
Fig. 4 depicts the profile ∆χ 2 scan on the test values of ∆m 2  31 with 6 years of JUNO and ORCA data taking.The four profiles correspond to true normal (top) and inverted (bottom) orderings while fitting the true or wrong ordering.Since the Asimov dataset is used, when assuming the true ordering on the fit both experiments show the same best-fit values at the true ∆m 2 31 and their ∆χ 2 minima yield zero.However, when assuming wrong ordering on the fit, the minima of ∆χ 2 are no longer at zero and the sensitivity to the NMO is obtained from this difference.After 6 years, JUNO will be able to exclude the wrong ordering with the significance of ∼ 2.3σ for either NMO.On the other hand, ORCA is expected to reach a significance of more than 6σ (3σ) for true NO (IO).31 for 6 years of data taking assuming baseline (solid) or optimistic (dashed) systematics.Fig. 4 also shows how the combination of JUNO and ORCA would exceed the NMO sensitivity of each experiment alone.The key advantage of the combination comes from the tension in ∆m 2  31 best fits of the two experiments when assuming the wrong ordering.This tension arises from the fact that each experiment observes neutrino oscillations starting from a different neutrino flavor (ν e for JUNO, ν µ + νµ for ORCA).Due to this difference the effective oscillation frequency will be a different combination of the various ∆m 2 ij for each experiment [29,30].Since the combination requires a single resulting ∆m 2  31 best fit, this tension together with strong constraints in ∆m 2  31 from both experiments, and particularly from JUNO, provides the synergy effect in which the combined ∆χ 2 minimum is enhanced to a higher value than simply adding the ∆χ 2 minima from each experiment.This latter scenario, in which the median sensitivity can be obtained as the square root of the sum, will be referred to as "simple sum" in the following discussion.It is shown only to highlight the benefit from doing the combination between JUNO and ORCA properly.
In Tab. 3, the NMO sensitivities after 6 years of collected data are presented for the combination, each experiment standalone, and the "simple sum" of their sensitivities.The combination reaches 8σ for true NO and 5σ for true IO.This combined sensitivity exceeds the "simple sum" case, which only obtains 7σ for true NO and 4σ for true IO.More importantly, a 5σ significance is achieved for both NMO scenarios within 6 years of combined analysis while each experiment Table 3: Asimov median sensitivity to NMO after 6 years of data taking for each experiment alone, the "simple sum", and the combination of the two experiments, assuming the baseline scenario for systematics.alone, or the "simple sum" of sensitivities, cannot achieve the same performance, sometimes even at significantly longer timescales.
The time evolution of the NMO sensitivity for JUNO, ORCA, and their corresponding combination is presented in Fig. 5 assuming that the two experiments start at the same time.JUNO alone would need 6-10 years of operation to reach 3σ of NMO sensitivity.ORCA has the capability to reach a 5σ significance after 3 years in the case of true NO.However, it would take more than 10 years of exposure to reach 5σ sensitivity in the case of IO.Due to the synergy effect discussed above, the combination would help significantly to reduce the time needed to reach a 5σ NMO sensitivity when compared to ORCA, especially if the neutrino mass ordering is inverted.With the combined analysis, a 5σ significance can be obtained within 2 (6) years in the case of true NO (IO) respectively.
As discussed in Sec.3.2, the ORCA analysis is also performed using a set of systematics similar to those of Ref. [33], as a cross-check for an optimistic approach.Fig. 4 and Fig. 5 show that both the optimistic and the baseline systematics give a very similar ∆χ 2 minimum value and thus yield the same NMO sensitivity for the ORCA-only analysis.However, the optimistic approach provides a much tighter constraint on ∆m 2  31 , as shown in Fig. 4, which causes the combination to reach sensitivities that are 1-2σ higher than in the case of the baseline scenario.This comes from the difference in the implementation of the energy scale systematics.The energy scale applied at the detector response (baseline) is more strongly correlated with ∆m  One of the most challenging design specifications of JUNO is the required energy resolution of the central detector.Reaching a level of about 3%/ E/MeV is essential for JUNO to be able to reach a 3σ sensitivity to determine the neutrino mass ordering by itself.In this sense, if the energy resolution worsens to 3.5%/ E/MeV, the required time to reach a 3σ sensitivity would increase by more than a factor of 2 [23].A significant amount of effort has been made within the JUNO Collaboration to reach this goal of 3%/ E/MeV, and a description of how to get there using a data-driven approach relying on calibration data has been presented in Ref. [66], where a 3.02%/ E/MeV energy resolution has been achieved, with a worsening of this energy resolution to 3.12%/ E/MeV after considering some imperfections in the detector.Nevertheless, it is still extremely interesting to evaluate the sensitivity of the combined NMO analysis to the energy resolution of JUNO.
In the present study a ±0.5%/E/MeV variation of the energy resolution was considered.While this accounts for a larger departure from the JUNO target energy resolution than the one described above, it was chosen to test the robustness of the combination procedure.As shown in Fig. 6, the impact of this variation of the energy resolution in the combined analysis is fairly small in comparison to the impact on the JUNO-only analysis.The reason for this small impact is that the added power to discriminate the neutrino mass orderings in this combination comes mostly from the displacement between the ∆m 2 31 best-fit values obtained by ORCA and JUNO for the wrong ordering assumption rather than from the direct measurement of the neutrino mass ordering in JUNO, as discussed previously.In this scenario, a worsening of the energy resolution would slightly reduce the precision of JUNO to measure ∆m 2  31 , while the best-fit value of ∆m 2  31 for each ordering would not change significantly.Therefore, the tension of the ∆m 2  31 best fit between JUNO and ORCA remains, which preserves the high sensitivity of the analysis.
As discussed in Sec.2.1, there is a possibility that 2 additional reactors could be built at the Taishan NPP, as originally planned.This would double the number of neutrinos produced by that NPP.In this scenario, JUNO by itself would be able to reach 3σ about 3 years earlier,  as shown in Fig. 7.In combination with ORCA however, there is negligible impact to the time required for the combined sensitivity to reach 5σ assuming true normal ordering, at the current best fit value.About 9 months are gained in the inverted ordering scenario which is still a significantly smaller impact than for the standalone JUNO.This behavior is, as in the case of the JUNO energy resolution dependency, due to the fact that the boost obtained from the combination relies on the difference between the JUNO and ORCA best-fit values for ∆m 2  31 in the wrong ordering scenario, rather than due to the precision of each experiment to measure the neutrino mass ordering separately.

Dependence on ∆m 2 31 and θ 23
This section presents the dependence of the analysis on the true value of the oscillation parameters, focusing particularly on θ 23 and ∆m 2 31 .Those parameters are chosen because the true value of θ 23 is known to have a strong influence on the ORCA sensitivity, and because the boost in sensitivity in the combined analysis, as discussed previously, is directly tied to the ∆m 2 31 measurement, and therefore it is essential to ensure that such a boost is valid for any true value of ∆m 2 31 .Fig. 8 shows the dependence of the NMO sensitivity on the true value of θ 23 for 6 years of data taking.As mentioned in Sec.4.1, JUNO has no sensitivity to θ 23 .In the case of ORCA however, the NMO sensitivity potential depends strongly on the true value of θ 23 as this parameter affects the amplitude of the detected oscillation pattern.After 6 years of data taking, ORCA has the sensitivity to reject the wrong ordering with a significance of 3-7σ and only reaches a 5σ sensitivity for true NO with θ 23 in the second octant.The combination curve also follows a similar θ 23 dependence as the ORCA-standalone curve, however thanks to the boost from JUNO, it is shifted to higher sensitivities and the joint fit ensures a 5σ discovery after about 6 years regardless of the true value of θ 23 and of the true NMO.
It is worth noting here that the current global best-fit value of θ 23 is in the upper octant with values of about 49 • for both orderings.This value is the one used in the studies described in Secs.4.2 and 5.1, which explains why in those studies the sensitivity for true NO is always much higher than 5σ after 6 years of data taking.Fig. 9 illustrates the dependence of NMO sensitivity on the true value of ∆m 2 31 .Both JUNO and ORCA standalone sensitivities depict a slight dependence on the true value of ∆m 2  31 with the opposite slope for each experiment.The combination is also quite flat with respect to the true ∆m 2  31 , reaching a significance of 8σ in the case of NO and 5σ in the case of IO.The effect of the boost described previously, relying on the difference between the wrong ordering measurement of ∆m 2  31 , is preserved over the whole ∆m 2 31 range.

Conclusions
This paper presents an evaluation of the sensitivity to the neutrino mass ordering achieved by a combined analysis of the JUNO and KM3NeT/ORCA experiments.It is worth pointing out explicitly that in all cases the combined analysis is more powerful than simply adding the sensitivities for both experiments together.As discussed above, this is due to the tension that arises in the ∆m 2  31 measurement between JUNO and ORCA when the wrong neutrino mass ordering is assumed.
The results show that this combination significantly reduces the time required to reach a 5σ determination of the neutrino mass ordering for any value of the oscillation parameters.In all cases, a 5σ measurement can be obtained within 6 years for the combined analysis, while it could take more than 10 years using only ORCA data, depending on the true ordering.
The gain in time is larger in cases where ORCA alone would require a longer time for reaching a 5σ sensitivity due to the uncertainty on the θ 23 value.In the favorable case of true normal ordering and θ 23 in the second octant, a 5σ NMO determination would be feasible after less than 2 years of data taking with the combined analysis.In this favorable scenario, which also corresponds to the current global best-fit value, the neutrino mass ordering would be determined at least a year ahead of what can be done using only ORCA data.
The boost for the NMO sensitivity obtained by combining JUNO and ORCA presented in this study is in line with what has been presented by previous studies considering the combination of JUNO with the IceCube Upgrade or with PINGU in Refs.[32,33].However, given the differences between PINGU and ORCA, it is important to confirm the result also for the combination of JUNO and ORCA.Of particular interest is the different treatment of the energy scale systematics between this and previous studies.This uncertainty impacts directly the ∆m 2  31 determination with ORCA and thus also the combined result.As shown in this paper, changing the treatment of this systematic uncertainty from an optimistic to a more realistic scenario may significantly affect the power of the combination of JUNO and ORCA.Nevertheless, a 5σ determination of the neutrino mass ordering can be effectively reached even in the ORCA baseline scenario for systematics.
Because the gain in time to reach the determination of the neutrino mass ordering in the combination of JUNO and ORCA does not come exclusively from each experiment's own ability to determine the neutrino mass ordering, the combination is sensitive to systematic uncertainties and detector effects in a different way than either experiment do independently.For instance, even if the JUNO energy resolution is critical for the measurement of the neutrino mass ordering using only JUNO data, it has only a small impact in the combined analysis.Alternatively, changing the ORCA systematics between optimistic and baseline systematics has a small impact on the power of ORCA alone to determine the neutrino mass ordering, however it has a larger impact on the combined analysis.These differences arise from the fact that the combination depends strongly on the measurement of ∆m 2  31 by each experiment rather than simply on their measurements of the neutrino mass ordering directly.In the cases where the gain in time to reach 5σ is of only a few years, this JUNO-ORCA combination is particularly interesting to provide a somewhat independent validation of the result obtained by the ORCA experiment alone, with a different dependency on the systematic uncertainties.

Figure 1 :
Figure 1: Expected event distribution for 6 years of data with JUNO 8 cores as a function of the visible energy of the prompt signal.The current world best-fit [44] oscillation parameters for normal ordering are assumed.The shaded region corresponds to the non-reactor neutrino background events.

Figure 2 :
Figure 2: Expected event distribution for ORCA in 3 PID classes for 6 years of exposure and true normal ordering assumption, with the oscillation parameter values from Ref. [44].

Figure 3 :
Figure 3: Implementation of the two different systematic approaches in the SWIM workflow used in the ORCA analysis.

Figure 4 :
Figure 4: ∆χ 2 profile for only JUNO (red), only ORCA (blue), and the combination of JUNO and ORCA (green) as a function of test values of ∆m 2 31 for 6 years of data taking assuming baseline (solid) or optimistic (dashed) systematics.

Figure 5 :
Figure 5: NMO sensitivity as a function of time for only JUNO (red), only ORCA (blue), and the combination of JUNO and ORCA (green), assuming baseline (solid) or optimistic (dashed) systematics.

Figure 6 :
Figure6: NMO sensitivity as a function of time for only JUNO (red), only ORCA (blue), and the combination of JUNO and ORCA (green), considering a better (dotted) and worse (dashed) energy resolution for JUNO than the nominal one (solid) by ±0.5%/ E/MeV.

Figure 7 :
Figure 7: NMO sensitivity as a function of time for only JUNO (red), only ORCA (blue), and the combination of JUNO and ORCA (green), considering 2 (solid) or 4 (dashed) Taishan NPP reactors, corresponding respectively to 8 or 10 reactor cores at 53 km from JUNO.

Figure 8 :
Figure 8: NMO sensitivity as a function of the true θ 23 value for 6 years of data taking for only JUNO (red), only ORCA (blue), and the combination of JUNO and ORCA (green).The vertical lines indicate the global best-fit values used in this analysis (from Ref. [44]).

Figure 9 :
Figure 9: NMO sensitivity as a function of the true ∆m 231 value for 6 years of data taking for only JUNO (red), only ORCA (blue), and the combination of JUNO and ORCA (green).The vertical lines indicate the global best-fit values used in this analysis (from Ref.[44]).
Deutsche Forschungsgemeinschaft (DFG), Germany; The General Secretariat of Research and Technology (GSRT), Greece; Istituto Nazionale di Fisica Nucleare (INFN), Ministero dell'Università e della Ricerca (MIUR), PRIN 2017 program (Grant NAT-NET 2017W4HA7S) Italy; Ministry of Higher Education Scientific Research and Professional Training, ICTP through Grant [26]tected variations in the parameters affecting the amount of light recorded by the detector, such as the absorption length and the PMT efficiencies, that would not be accounted for in the reconstruction[26].Finally, normalization factors for each PID class are also included to account for any possible systematic effects (in the flux, cross-section, or detector response) that would vary the total number of events in each class.

Table 2 :
[44]al best-fit values for the oscillation parameters (from Ref.[44]) and assumed to be the "true value" in this analysis.Uncertainties are shown for the parameter where a prior based on the global best-fit value was used. of those oscillation parameters, and the last term corresponds to the added prior with sin 2 θ GF 13 being the current global best fit for sin 2 θ 13 and σ sin 2 θ GF