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

S. Aiello, A. Albert, M. Alshamsi, S. Alves Garre, Z. Aly, A. Ambrosone, F.

Ameli, M. Andre, G. Androulakis, M. Anghinolfi, et al.

► To cite this version:

S. Aiello, A. Albert, M. Alshamsi, S. Alves Garre, Z. Aly, et al. Combined sensitivity of JUNO and KM3NeT/ORCA to the neutrino mass ordering. Journal of High Energy Physics, 2022, 03, pp.55. 10.1007/JHEP03(2022)055. hal-03376783

HAL Id: hal-03376783

https://hal.science/hal-03376783

Submitted on 15 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Combined sensitivity of JUNO and KM3NeT/ORCA to the neutrino mass ordering

## KM3NeT collaboration

S. Aiello,ª A. Albert,ᵇba M. Alshamsi,ᶜ S. Alves Garre,ᵈ Z. Aly,ᵉ A. Ambrosone,ᵉ,g

F. Ameli,ʰ M. Andre,ᵢ G. Androulakis,ⱼ,<fnref n="1" /> M. Anghinolfi,ᵏ M. Anguita,ᶜ M. Ardid,ᵇ

S. Ardid,ᵇ M. Aublin,ᶜ C. Bagatelas,ⱼ B. Baret,ᶜ S. Basegmez du Pree,ⁿ

M. Bendahman,ᶜ,'o F. Benfenati,ᶜ,q E. Berbee,ⁿ A.M. van den Berg,ᵈ V. Bertin,ᵉ

S. Biagi,ᶜ M. Boettcher,ᵗ M. Bou Cabo,ᵈ J. Boumaaza,'o M. Bouta,ᵈ M. Bouwhuis,ⁿ

C. Bozza,ᵈ H.Brânzaş,ₓ R. Bruijn,ⁿ,y J. Brunner,ᵉ R. Bruno,ᵃ E. Buis,ᵈ

R. Buompane,ᵉ,aa J. Busto,ᵉ B. Caiffi,ᵏ D. Calvo,ᵈ S. Campion,ᵃb,ḥ A. Capone,ᵃb,ḥ

V. Carretero,ᵈ P. Castaldi,ᵉ,ac S. Celli,ᵃb,ḥ M. Chabab,ᵃd N. Chau,ᶜ,<fnref n="2" /> A. Chen,ᵃe

S. Cherubini,ᶜ,af V. Chiarella,ag T. Chiarusi,ᵇ M. Circella,ah R. Cocimano,ᶜ

J.A.B. Coelho,ᶜ A. Coleiro,ᶜ M. Colomer-Molla,ᶜ,d R. Coniglione,ᶜ P. Coyle,ᵉ

A. Creusot,ᶜ A. Cruz,ai G. Cuttone,ᶜ R. Dallier,aj B. De Martino,ᵉ I. Di Palma,ᵃb,ḥ

A.F. Díaz,ᵈ D. Diego-Tortosa,ᵇ C. Distefano,ᶜ A. Domi,ⁿ,y C. Donzaud,ᶜ D. Dornic,ᵉ

M. Dörr,ak D. Drouhin,ᵇba T. Eberl,al A. Eddyamoui,'o T. van Eeden,ⁿ D. van Eijk,ⁿ

I. El Bojaddaini,ᵈ A. Enzenhöfer,ᵉ V. Espinosa,ᵇ P. Fermani,ᵃb,ḥ G. Ferrara,ᶜ,af

M. D. Filipović,am F. Filippini,ᵉ,q L.A. Fusco,ᵉ T. Gal,al J. García Méndez,ᵃm

A. Garcia Soto,ᵈ F. Garufi,ᵉ,g Y. Gatelet,ᶜ C. Gatius,ⁿ N. Geißelbrecht,al

L. Gialanella,ᵉ,aa E. Giorgio,ᶜ S.R. Gozzini,ᵈ R. Gracia,ⁿ K. Graf,al G. Grella,an

D. Guderian,ᵇᵇ C. Guidi,ᵏ,ao B. Guillon,ᵃp M. Gutiérrez,aq J. Haefner,al S. Hallmann,al

H. Hamdaoui,'o H. van Haren,ar A. Heijboer,ⁿ A. Hekalo,ak L. Hennig,al

J.J. Hernández-Rey,ᵈ J. Hofestädt,al F. Huang,ᵉ W. Idrissi Ibnsalih,ᵉ,aa

G. Illuminati,ᵉ,ᶜ C.W. James,ai D. Janezashvili,as M. de Jong,ⁿ,at P. de Jong,ⁿ,y

B.J. Jung,ⁿ P. Kalaczyński,au O. Kalekin,al U.F. Katz,al N.R. Khan Chowdhury,ᵈ

G. Kistauri,as F. van der Knaap,ᵈ P. Kooijman,ᵧ,bc A. Kouchner,ᶜ,av V. Kulikovskiy,ᵏ

M. Labalme,ᵃp R. Lahmann,al M. Lamoureux,ᶜ,<fnref n="3" /> G. Larosa,ᶜ C. Lastoria,ᵉ A. Lazo,ᵈ

R. Le Breton,ᶜ S. Le Stum,ᵉ G. Lehaut,ᵃp O. Leonardi,ᶜ F. Leone,ᶜ,af E. Leonora,ª

<footnote n="1">¹Deceased.</footnote>

<footnote n="2">²Corresponding author.</footnote>

<footnote n="3">³Also at: Dipartimento di Fisica, INFN Sezione di Padova and Università di Padova, I-35131, Padova, Italy.</footnote>
N. Lessing,ª G. Levi,ª,q M. Lincetto,ª M. Lindsey Clark,ª C. Lipreau,ª,j

C. Llorens Alvarez,ª F. Longhitano,ª D. Lopez-Coto,ª,q A. Lygda,ª L. Maderer,ª

J. Majumdar,ª J. Mańczak,ª A. Margiotta,ª,q A. Marinelli,ª C. Markou,ª L. Martin,ª,j

J.A. Martínez-Mora,ª A. Martini,ª F. Marzaioli,ª,aa S. Mastroianni,ª K.W. Melis,ª

G. Miele,ª,g P. Migliozzi,ª E. Migneco,ª P. Mijakowski,ªu L.S. Miranda,ªw

C.M. Mollo,ª M. Moser,ªl A. Moussa,ª R. Muller,ª M. Musumeci,ª S. L. Nauta,ª

S. Navas,ªq C.A. Nicolau,ª B. Nkosi,ªe B. Ó Fearraigh,ª,n,y M. O'Sullivan,ªi

M. Organokov,ª A. Orlando,ª J. Palacios González,ª G. Papalashvili,ªs R. Papaleo,ª

A. M. Păun,ª G.E. Păvălaş,ª C. Pellegrino,ª,bd M. Perrin-Terrin,ª V. Pestel,ª

P. Piattelli,ª C. Pieterse,ª O. Pisanti,ª,g C. Poirè,ª V. Popa,ª T. Pradier,ª I. Probst,ªl

S. Pulvirenti,ª G. Quéméner,ªp N. Randazzo,ª S. Razzaque,ªw D. Real,ª S. Reck,ªl

G. Riccobene,ª A. Romanov,ª,ao A. Rovelli,ª F. Salesa Greus,ª D.F.E. Samtleben,ª,n,at

A. Sánchez Losa,ª,h,d M. Sanguineti,ª,k,ao D. Santonocito,ª P. Sapienza,ª J. Schnabel,ªl

M.F. Schneider,ªl J. Schumann,ªl H. M. Schutte,ª J. Seneca,ª I. Sgura,ªh

R. Shanidze,ªs A. Sharma,ªx A. Sinopoulou,ª B. Spisso,ªn,f M. Spurio,ª,p,q

D. Stavropoulos,ª S.M. Stellacci,ªn,f M. Taiuti,ª,k,ao Y. Tayalati,ª H. Thiersen,ª

S. Tingay,ªi S. Tsagkli,ª V. Tsourapis,ª E. Tzamariudaki,ª D. Tzanetatos,ª

V. Van Elewyck,ª,av,2 G. Vasileiadis,ªy F. Versari,ª,p,q D. Vivolo,ª,f,aa G. de Wasseige,ª

J. Wilms,ªz R. Wojaczyński,ªu E. de Wolf,ª,n,y T. Yousfi,ª S. Zavatarelli,ª

A. Zegarelli,ª,b,h D. Zito,ª J.D. Zornoza,ª J. Zúñigaªd and N. Zywuckaªt

## JUNO collaboration members

S. Ahmad,ªn J.P.A.M. de André,ª,h,2 E. Baussan,ªh C. Bordereau,ª,m,bq A. Cabrera,ª,f

C. Cerna,ªq G. Donchenko,ªbl E. A. Doroshkevich,ªi M. Dracos,ªbh F. Druillole,ªq

C. Jollet,ªq L.N. Kalousis,ªbh P. Kampmann,ªbe K. Kouzakov,ªbl A. Lokhov,ª,bl

B.K. Lubsandorzhiev,ªi S.B. Lubsandorzhiev,ªi A. Meregaglia,ªq L. Miramonti,ªg

F. Perrot,ªq L.F. Piñeres Rico,ªbh A. Popov,ªbl R. Rasheed,ªq M. Settimo,ªo

K. Stankevich,ªbl H. Steiger,ªk,bp M.R. Stock,ªbp A. Studenikin,ªbl A. Triossi,ªbh

W. Trzaska,ªbr M. Vialkov,ªbl B. Wonsak,ªbj J. Wurtzªbh and F. Yermiaªbo

<footnote>ªINFN, Sezione di Catania, Via Santa Sofia 64, Catania, 95123 Italy</footnote>

<footnote>ᵇUniversité de Strasbourg, CNRS, IPHC UMR 7178, F-67000 Strasbourg, France</footnote>

<footnote>ᶜUniversité de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France</footnote>

<footnote>ᵈIFIC — Instituto de Física Corpuscular (CSIC — Universitat de València), c/Catedrático José Beltrán, 2, 46980 Paterna, Valencia, Spain</footnote>

<footnote>ᵉAix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France</footnote>

<footnote>fINFN, Sezione di Napoli, Complesso Universitario di Monte S. Angelo, Via Cintia ed. G, Napoli, 80126 Italy</footnote>

<footnote n="9">⁹Università di Napoli “Federico II”, Dip. Scienze Fisiche “E. Pancini”, Complesso Universitario di Monte S. Angelo, Via Cintia ed. G, Napoli, 80126 Italy</footnote>

<footnote>hINFN, Sezione di Roma, Piazzale Aldo Moro 2, Roma, 00185 Italy</footnote>

<footnote>iUniversitat Politècnica de Catalunya, Laboratori d'Aplicacions Bioacústiques, Centre Tecnològic de Vilanova i la Geltrú, Avda. Rambla Exposició, s/n, Vilanova i la Geltrú, 08800 Spain</footnote>

<footnote>jNCSR Demokritos, Institute of Nuclear and Particle Physics, Ag. Paraskevi Attikis, Athens, 15310 Greece</footnote>
kINFN, Sezione di Genova, Via Dodecaneso 33, Genova, 16146 Italy

lUniversity of Granada, Dept. of Computer Architecture and Technology/CITIC, 18071 Granada, Spain

mUniversitat Politècnica de València, Instituto de Investigación para la Gestión Integrada de las Zonas Costeras, C/ Paraním, 1, Gandia, 46730 Spain

nNikhef, National Institute for Subatomic Physics, PO Box 41882, Amsterdam, 1009 DB, Netherlands

oUniversity Mohammed V in Rabat, Faculty of Sciences, 4 av. Ibn Battouta, B.P. 1014, R.P. 10000 Rabat, Morocco

pINFN, Sezione di Bologna, v.le C. Berti-Pichat, 6/2, Bologna, 40127 Italy

qUniversità di Bologna, Dipartimento di Fisica e Astronomia, v.le C. Berti-Pichat, 6/2, Bologna, 40127 Italy

rKVI-CART University of Groningen, Groningen, Netherlands

sINFN, Laboratori Nazionali del Sud, Via S. Sofia 62, Catania, 95123 Italy

tNorth-West University, Centre for Space Research, Private Bag X6001, Potchefstroom, 2520 South Africa

uInstituto Español de Oceanografía, Unidad Mixta IEO-UPV, C/ Paraním, 1, Gandia, 46730 Spain

vUniversity Mohammed I, Faculty of Sciences, BV Mohammed VI, B.P. 717, R.P. 60000 Oujda, Morocco

wUniversità di Salerno e INFN Gruppo Collegato di Salerno, Dipartimento di Matematica, Via Giovanni Paolo II 132, Fisciano, 84084 Italy

xISS, Atomistilor 409, Măgurele, RO-077125 Romania

yUniversity of Amsterdam, Institute of Physics/IHEF, PO Box 94216, Amsterdam, 1090 GE, Netherlands

zTNO, Technical Sciences, PO Box 155, Delft, 2600 AD, Netherlands

aaUniversità degli Studi della Campania “Luigi Vanvitelli”, Dipartimento di Matematica e Fisica, viale Lincoln 5, Caserta, 81100 Italy

abUniversità La Sapienza, Dipartimento di Fisica, Piazzale Aldo Moro 2, Roma, 00185 Italy

acUniversità di Bologna, Dipartimento di Ingegneria dell'Energia Elettrica e dell'Informazione "Guglielmo Marconi", Via dell'Università 50, Cesena, 47521 Italia

adCadi Ayyad University, Physics Department, Faculty of Science Semlalia, Av. My Abdellah, P.O.B. 2390, Marrakech, 40000 Morocco

aeUniversity of the Witwatersrand, School of Physics, Private Bag 3, Johannesburg, Wits 2050 South Africa

afUniversità di Catania, Dipartimento di Fisica e Astronomia "Ettore Majorana", Via Santa Sofia 64, Catania, 95123 Italy

agINFN, LNF, Via Enrico Fermi, 40, Frascati, 00044 Italy

ahINFN, Sezione di Bari, via Orabona, 4, Bari, 70125 Italy

aiInternational Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia

ajSubatech, IMT Atlantique, IN2P3-CNRS, Université de Nantes, 4 rue Alfred Kastler — La Chantrerie, Nantes, BP 20722 44307 France

akUniversity Würzburg, Emil-Fischer-Straße 31, Würzburg, 97074 Germany

alFriedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Straße 1, 91058 Erlangen, Germany
am Western Sydney University, School of Computing, Engineering and Mathematics,
Locked Bag 1797, Penrith, NSW 2751 Australia

an Università di Salerno e INFN Gruppo Collegato di Salerno, Dipartimento di Fisica, Via Giovanni Paolo II 132, Fisciano, 84084 Italy

ao Università di Genova, Via Dodecaneso 33, Genova, 16146 Italy

ap Normandie Univ, ENSICAEN, UNICAEN, CNRS/IN2P3, LPC Caen, LPCCAEN, 6 boulevard Maréchal Juin, Caen, 14050 France

aq University of Granada, Dpto. de Física Teórica y del Cosmos & C.A.F.P.E.,
18071 Granada, Spain

ar NIOZ (Royal Netherlands Institute for Sea Research), PO Box 59, Den Burg, Texel, 1790 AB, Netherlands

as Tbilisi State University, Department of Physics, 3, Chavchavadze Ave., Tbilisi, 0179 Georgia

at Leiden University, Leiden Institute of Physics, PO Box 9504, Leiden, 2300 RA, Netherlands

au National Centre for Nuclear Research, 02-093 Warsaw, Poland

av Institut Universitaire de France, 1 rue Descartes, Paris, 75005 France

aw University of Johannesburg, Department Physics, PO Box 524, Auckland Park, 2006 South Africa

ax Università di Pisa, Dipartimento di Fisica, Largo Bruno Pontecorvo 3, Pisa, 56127 Italy

ay Laboratoire Univers et Particules de Montpellier, Place Eugène Bataillon — CC 72, Montpellier
Cédex 05, 34095 France

az Friedrich-Alexander-Universität Erlangen-Nürnberg, Remeis Sternwarte, Sternwartstraße 7,
96049 Bamberg, Germany

ba Université de Haute Alsace, rue des Frères Lumière, 68093 Mulhouse Cedex, France

bb University of Münster, Institut für Kernphysik, Wilhelm-Klemm-Str. 9, Münster, 48149 Germany

bc Utrecht University, Department of Physics and Astronomy, PO Box 80000,
Utrecht, 3508 TA, Netherlands

bd INFN, CNAF, v.le C. Berti-Pichat, 6/2, Bologna, 40127 Italy

be Forschungszentrum Jülich GmbH, Nuclear Physics Institute IKP-2, Jülich, Germany

bf IJCLab, Université Paris-Saclay, CNRS/IN2P3, 91405 Orsay, France

bg INFN Sezione di Milano and Dipartimento di Fisica dell'Università di Milano, Milano, Italy

bh IPHC, Université de Strasbourg, CNRS/IN2P3, F-67037 Strasbourg, France

bi Institute for Nuclear Research of the Russian Academy of Sciences, Moscow, Russia

bj Institute of Experimental Physics, University of Hamburg, Hamburg, Germany

bk Institute of Physics, Johannes-Gutenberg Universität Mainz, Mainz, Germany

bl Lomonosov Moscow State University, Moscow, Russia

bm National United University, Miao-Li

bn Pakistan Institute of Nuclear Science and Technology, Islamabad, Pakistan

bo SUBATECH, Université de Nantes, IMT Atlantique, CNRS-IN2P3, Nantes, France

bp Technische Universität München, München, Germany

bq Univ. Bordeaux, CNRS, CENBG, UMR 5797, F-33170 Gradignan, France

br University of Jyvaskyla, Department of Physics, Jyvaskyla, Finland

E-mail: nchau@apc.in2p3.fr, elewyck@apc.in2p3.fr, jpandre@iphc.cnrs.fr
**ABSTRACT:** 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.

**KEYWORDS:** Neutrino Detectors and Telescopes (experiments), Oscillation

ARXIV ePRINT: 2108.06293
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Jiangmen Underground Neutrino Observatory (JUNO)</td><td>3</td></tr><tr><td>2.1 Modeling JUNO for this study</td><td>4</td></tr><tr><td>2.2 Sensitivity analysis</td><td>5</td></tr><tr><td>3 Oscillation Research with Cosmics in the Abyss (ORCA)</td><td>7</td></tr><tr><td>3.1 Modeling the ORCA detector for this study</td><td>8</td></tr><tr><td>3.2 Sensitivity analysis</td><td>10</td></tr><tr><td>4 Combined JUNO and ORCA analysis</td><td>13</td></tr><tr><td>4.1 Combination strategy</td><td>13</td></tr><tr><td>4.2 Results</td><td>14</td></tr><tr><td>5 Further sensitivity studies</td><td>17</td></tr><tr><td>5.1 Impact of energy resolution in JUNO and 10 reactor cores scenario</td><td>17</td></tr><tr><td>5.2 Dependence on Δm<sup>2</sup><sub>31</sub> and θ<sub>23</sub></td><td>19</td></tr><tr><td>6 Conclusions</td><td>21</td></tr></table>

## 1 Introduction

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-3]). Since the first conclusive observations of neutrino oscillations at the turn of the century [4-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-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 CP-violation in the neutrino sector, and the ordering of the masses ($m_1, m_2$ and $m_3$) of the neutrino mass eigenstates (respectively, $\nu_1, \nu_2$ and $\nu_3$): either normal ordering (NO, $m_1 < m_2 < m_3$) or inverted ordering (IO, $m_3 < m_1 < m_2$).

The question of the mass ordering is a prime experimental goal because its determination will have direct consequences on other branches of neutrino physics, including the
measurement of leptonic CP violation in future long baseline experiments [11], the extraction of cosmological limits on the absolute neutrino mass [12, 13], and the interpretation of results from planned experiments searching for neutrino-less double-beta decay to establish the Dirac vs. Majorana nature of neutrinos [14]. In a broader context, establishing the neutrino mass ordering is also important to help unveil the symmetries hidden behind the flavour mixing, which are connected to the fundamental mechanism giving mass to neutrinos, and to constrain the related models [15, 16].

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 $\nu_\mu-\nu_e$ oscillations in Earth matter. Such experiments are sensitive to the atmospheric mass splitting $\Delta m_{31}^2$ [17]. However, none of these experiments alone, either accelerator-based (such as T2K [18] or NOvA [19]) or using atmospheric neutrinos (such as Super-Kamiokande [20] or IceCube [21]), has the capability to definitively resolve the NMO (in other words, the sign of $\Delta m_{31}^2$) 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 [22] and NOvA [23]. Considering that a high-confidence ($> 5\sigma$) determination of the NMO with the next-generation accelerator experiments DUNE [24], T2HK [25] and T2HKK [26] is only envisaged for 2030 or beyond, alternative paths to the NMO measurement are being pursued on a shorter timescale.

JUNO [27–29] and KM3NeT/ORCA [30] 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 ($\sim 53$ km) reactor neutrino experiment that is sensitive to the NMO through the interplay between the fast oscillations driven by $\Delta m_{31}^2$ and $\Delta m_{32}^2$ in the $\bar{\nu}_e$ disappearance channel, where matter effects play only a small role [31]. As reactor $\bar{\nu}_e$ disappearance measurement is not affected by CP violation [27], the JUNO measurement will be independent of the unknown CP violating phase $\delta_{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 [32].

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. As further discussed in section 4.2, this boost essentially comes from the expected tension between ORCA and JUNO in the best fit of $|\Delta m_{31}^2|$ when the assumed ordering is wrong. Such tension arises from the fact that the two experiments use different neutrino oscillation channels to extract the value of the $|\Delta m_{31}^2|$ parameter: $\bar{\nu}_e$ disappearance for JUNO, and mainly $\nu_\mu(+\bar{\nu}_\mu)$ disappearance for ORCA. Because these channels are sensitive to different combinations of the squared-mass splittings, the resulting value of $|\Delta m_{31}^2|$ obtained with the wrong ordering assumption is different in the two types of experiments. This difference turns out to be
typically larger than the respective uncertainties expected in ORCA and JUNO. The wrong mass ordering can therefore 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 [33, 34], 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 [35], 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 [36]. The potential of this method was further explored in a combined study with JUNO and PINGU using detailed simulation tools for both experiments [37], 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 section 2, along with the standalone $\chi^2$ analysis used to determine the JUNO-only NMO sensitivity. The same is done for ORCA in section 3, based on the updated detector configuration, and simulation and reconstruction tools used for the latest NMO sensitivity study following ref. [38]. 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. [37] for PINGU). That systematic effect is shown to degrade the precision of the $\Delta m^2_{31}$ measurement of ORCA alone, thereby affecting the sensitivity of the combined study.

The JUNO/ORCA combined $\chi^2$ analysis is presented in section 4 for the baseline reactor configuration of JUNO, and the most realistic systematics treatment adopted for ORCA. The enhanced sensitivity achieved with the combined $\chi^2$ analysis over the simple sum of individual $\chi^2$ is also demonstrated. Section 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 $\theta_{23}$ and $\Delta m^2_{31}$ is also addressed. The conclusions drawn from these results are presented in section 6.

## 2 Jiangmen Underground Neutrino Observatory (JUNO)

The Jiangmen Underground Neutrino Observatory [27–29] (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 $\bar{\nu}_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 $\bar{\nu}_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 $\bar{\nu}_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 $\bar{\nu}_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 [39]. This detector will be installed at a distance of 30 m from one of the Taishan's reactors to determine the reactor $\bar{\nu}_e$ spectrum with a better energy resolution than JUNO, effectively reducing the impact of possible unknown substructures in the reactor neutrino spectra [40] on the measurement of neutrino oscillations.

The precise distance of JUNO to each reactor core of the Yangjiang and Taishan NPPs, provided in ref. [27], 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 [29] with 17.4 GW thermal power.

## 2.1 Modeling JUNO for this study

For this study, the performance of JUNO closely follows that provided in ref. [27]. In particular, a 73% IBD detection efficiency and an energy resolution of $3\%/\sqrt{E/\text{MeV}}$ are assumed in the Central Detector which contains $1.5 \times 10^{33}$ target protons. Given that the energy resolution is critical for the JUNO sensitivity, the impact of a $\pm 0.5\%/\sqrt{E/\text{MeV}}$ change in energy resolution is discussed in section 5.1. As in ref. [27], the nominal running time of the experiment is considered to be 1000 effective days every 3 years.

The observed $\bar{\nu}_e$ spectrum in JUNO will be produced by the interplay between the spectrum produced by the NPPs, the IBD cross-section [41], and the neutrino oscillations
which are to be measured. To determine the spectrum produced by the NPPs, the ILL<fnref n="1" />
$\bar{\nu}_e$ spectra [42–44] are used. The ILL-Vogel flux model was chosen for this study given
its flux normalization better matches the latest reactor data [45], since previous studies
found consistent results for all tested spectra, as discussed in ref. [27]. The fine structure of
the spectrum will be precisely determined independently using JUNO-TAO, as mentioned
beforehand.<fnref n="2" /> 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 [46], and using the fission
energies for each of these isotopes from ref. [47].

In addition to reactor neutrinos, the IBD event selection contains some background
events. The backgrounds considered in this analysis are taken from ref. [27], 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 [48] oscillation parameters,
considering all 4 Taishan NPP reactors operational.

A notable difference from ref. [27] 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. [27]. For completeness, the JUNO official baseline with a total of 4 Taishan reactors
is also considered and discussed in section 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 $\bar{\nu}_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 figure 1,
for 6 years of data with 8 cores and assuming the best-fit oscillation parameters from
ref. [48] 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.

## 2.2 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.

<footnote n="1">¹Institut Laue-Langevin.</footnote>

<footnote n="2">²As it is discussed in ref. [39], the JUNO and JUNO-TAO detectors do not see the same fission fraction from the NPPs. While it is important to take into account this difference, it was estimated that using the JUNO-TAO spectrum as input corresponds to about 0.9% bin-to-bin uncertainty in the $\bar{\nu}_e$ spectrum. This value is smaller than the 1% bin-to-bin uncertainty used in the study presented here.</footnote>
<figure><img src="image_2.png" /><figcaption>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 [48] oscillation parameters for normal ordering are assumed. The shaded region corresponds to the non-reactor neutrino background events. The red line corresponds to the neutrino flux from just the Daya-Bay and Huizhou NPPs, which are not in the correct distance from JUNO to help measure the NMO. Statistical error bars on the total number of events for each bin are also shown. Right: difference in number of events expected between normal and inverted ordering scenarios assuming current world best-fit [48] oscillation parameters.</figcaption></figure>

As shown in figure 1, this spectrum exhibits two notable features. The first is a slow oscillatory behavior due to $\theta_{12}$ and $\Delta m^2_{21}$, which causes a large deficit in the number of events in the whole energy range shown in figure 1 and that has a minimum at about 2.2 MeV. The second is a rapid oscillatory behavior due to $\theta_{13}$ and $\Delta 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 figure 1, the position of the rapid oscillation maxima and minima would change, as the oscillation frequencies producing the pattern would be different [27, 49]. This is the signature that JUNO will use to measure the NMO. The sensitivity of JUNO to determine the NMO is calculated using the $\chi^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 section 4, it is useful to describe already at this point the $\chi^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 $\chi^2$ values. This approach is commonly referred to as “Asimov” approach [50], 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 $\chi^2$ function used in this analysis has the following form:

$$ \chi^2 = \Delta \times \mathcal{M}^{-1} \times \Delta^T. \qquad (2.1) $$

In this expression, $\Delta$ is a $1 \times n$ matrix whose content is the difference between the observed and expected rates. $\Delta$ is defined as $\Delta = D - (S + B)$, where $D$, $S$, and $B$ correspond, respectively, to the data, the signal prediction for a given set of oscillation
<table><thead><tr><td>Parameter</td><td>Uncertainty</td></tr></thead><tbody><tr><td>Correlated reactor error</td><td>2.0%</td></tr><tr><td>Uncorrelated reactor error</td><td>0.8%</td></tr><tr><td>Reactor spectrum</td><td>1.0%</td></tr><tr><td>Detector response</td><td>1.0%</td></tr></tbody></table>

<table><caption>**Table 1.** Error envelope used to build covariance matrices for signal (left) and background (right) on JUNO analysis. For details see ref. [27].</caption><thead><tr><th rowspan="2">Background</th><th rowspan="2">Rate<br>(/day)</th><th colspan="2">Uncertainty on</th></tr><tr><th>rate</th><th>shape</th></tr></thead><tbody><tr><td>Cosmogenic</td><td>1.6</td><td>20%</td><td>10%</td></tr><tr><td>Geo-neutrino</td><td>1.1</td><td>30%</td><td>5%</td></tr><tr><td>Accidental</td><td>0.9</td><td>1%</td><td>0%</td></tr><tr><td>Fast neutrons</td><td>0.1</td><td>100%</td><td>20%</td></tr><tr><td><sup>13</sup>C(α, n)<sup>16</sup>O</td><td>0.05</td><td>50%</td><td>50%</td></tr></tbody></table>

parameters for all NPPs, and the expected non-reactor 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. $\Delta^T$ is the transpose matrix of $\Delta$.

The matrix $\mathcal{M}$ whose inverse is present in eq. (2.1) is a $n \times n$ covariance matrix. This matrix is calculated as $\mathcal{M} = \mathcal{M}_{stat} + \mathcal{M}_S + \mathcal{M}_B$. $\mathcal{M}_{stat}$ corresponds to the statistical uncertainty in each bin from the total expected number of events ($S+B$). $\mathcal{M}_S$ and $\mathcal{M}_B$ correspond, respectively, to the covariance matrices of the signal and background. These matrices are created by the sum of the covariance matrices associated with each systematic error shown in table 1 obtained from ref. [27].

The $\chi^2$ defined above is then calculated using different values of the oscillation parameters to estimate the $S$, while keeping $D$ and $B$ fixed at the true assumed value. It is worth noting that given that the $\chi^2$ defined above already accounts for uncertainties from the various systematic errors via the covariance matrices, there is no fitting of the $\chi^2$ with separate nuisance parameters as had been done in ref. [27]. This JUNO-only analysis was validated by comparing the NMO sensitivity with previous results in refs. [27, 37]. 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 $\chi^2$ units.

### 3 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 [30]. 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 \times 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 [51]. 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 $\nu_\mu/\bar{\nu}_\mu$, and $\nu_\tau/\bar{\nu}_\tau$ 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 $\nu_e/\bar{\nu}_e$ and $\nu_\tau/\bar{\nu}_\tau$ with non-muonic decays, as well as NC interactions of all flavors.

## 3.1 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 [52], a GENIE [53]-based software developed within the KM3NeT Collaboration. All secondary particles are tracked with the package KM3Sim [54] based on GEANT4 [55] 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 [56]. 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 [38].

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. [38]. From the detected signals, the energy and direction of the events are reconstructed using dedicated algorithms developed for shower-like [57] and track-like [58] 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 $\theta_{reco} > 90^\circ$, conventionally defined as the angle between the direction to the neutrino source and the upward direction at the detector), 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 $\eta$ 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 < \eta \le 1$), intermediate events ($0.3 < \eta \le 0.7$), and showers ($0 \le \eta \le 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. [38], 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, \theta$) distributions of atmospheric neutrino events for each PID class. Such distributions are obtained with the SWIM package [59], 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. [60], 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 [61], 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 [62]. The rate of events interacting around the detector is computed for each interaction type $\nu_x \in \{\bar{\nu}_e \text{ CC}, \bar{\nu}_\mu \text{ CC}, \bar{\nu}_\tau \text{ CC}, \bar{\nu}_\text{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 $\nu_x$ and true variables ($E_{\text{true}}, \theta_{\text{true}}$) into the corresponding reconstructed variables ($E_{\text{reco}}, \theta_{\text{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_{\text{true}}, \theta_{\text{true}}$) $\rightarrow$ ($E_{\text{reco}}, \theta_{\text{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 accounts for detection and reconstruction efficiencies, misidentification probabilities and errors on reconstructed variables (including all correlations). This approach is different from the one in ref. [38] 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 [63] has been adopted, as described in the next section.

Figure 2 depicts the expected neutrino event distributions calculated with SWIM for all PID classes, binned in reconstructed variables ($E_{\text{reco}}, \cos \theta_{\text{reco}}$) for 6 years of data taking with ORCA, assuming normal ordering and oscillation parameters from ref. [48]. In this
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Left: 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. [48]. Right: difference in expected event distribution between normal and inverted ordering scenarios assuming oscillation parameter values from ref. [48].</figcaption></figure>

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.

## 3.2 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 $\chi^2$ is built as follows:

$$ \chi^2 = -2 \sum_l \left( d_l - \beta_l \mu_l - d_l \ln \frac{d_l}{\beta_l \mu_l} \right) + \sum_j \frac{(p_j - p_j^0)^2}{\sigma_j^2} + \sum_k \frac{(1 - \beta_k)^2}{\sigma_k^2}. \quad (3.1) $$

The first sum is the Poisson likelihood of the data $d_l$ given the expectation $\mu_l$ at bin $l$, where the latter also depends on the nuisance parameters $p_j$. The $\beta$ parameters are introduced based on the “Beeston-Barlow light method” [63] 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_j^0$ and variances $\sigma_j$. The third sum represents the Gaussian priors to the $\beta$ parameters that are expected to be normally distributed.

The fluctuations $\beta_k$ are assumed to be bin-to-bin uncorrelated and independent of the model parameters. Given this assumption, the values of $\beta_k$ are solved analytically as the solution of $\partial\chi^2/\partial\beta_k = 0$. In addition, the variances $\sigma_k$ of $\beta_k$ can also be evaluated with a probabilistic model which describes the calculation of the response matrix as a single binomial process [64, 65]. Finally, both $\beta$ and its variance can be estimated analytically and used directly in the calculation of the $\chi^2$ without any requirement for additional minimization. The full description of this procedure can also be found in ref. [59]. This implementation results in a $\sim 0.2\sigma$ 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 table 2). The “baseline” scenario corresponds to the standard set of ORCA systematics adopted for oscillation analyses, similar to ref. [38]. Uncertainties related to the incident flux include the spectral index of the neutrino flux energy distribution (free without any constraints) and the flux ratio systematics. These parameters are introduced to describe the uncertainties in the ratio of different event types, namely $\nu_e/\bar{\nu}_e$, $\nu_\mu/\bar{\nu}_\mu$, and $(\nu_e + \bar{\nu}_e)/(\nu_\mu + \bar{\nu}_\mu)$, while preserving in each case the total associated flux [59]. They are constrained with the priors adapted from ref. [66]. 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 the effect of undetected 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 [30]. 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.

A second scenario is based on the study of ref. [37], developed for the PINGU detector.<fnref n="3" /> That analysis does not apply normalization factors to the PID classes. It uses an overall

<footnote n="3">³This approach is different from the one adopted in some recent studies performed with the IceCube DeepCore Detector, which include uncertainties related to the PMT optical efficiency and bulk ice properties [67, 68].</footnote>
<table><caption>**Table 2.** Baseline (see ref. [38]) and optimistic (see ref. [37]) scenarios for the treatment of systematics considered in the ORCA analysis. The cross (×) indicates that the systematic is not included.</caption><thead><tr><th>Parameter</th><th>Baseline scenario</th><th>Optimistic scenario</th></tr></thead><tbody><tr><td>Flux spectral index</td><td></td><td>free</td></tr><tr><td>Flux ν<sub>e</sub>/ν&#773;<sub>e</sub> ratio</td><td></td><td>7% prior</td></tr><tr><td>Flux ν<sub>μ</sub>/ν&#773;<sub>μ</sub> ratio</td><td></td><td>5% prior</td></tr><tr><td>Flux (ν<sub>e</sub> + ν&#773;<sub>e</sub>)/(ν<sub>μ</sub> + ν&#773;<sub>μ</sub>) ratio</td><td></td><td>2% prior</td></tr><tr><td>NC normalization</td><td></td><td>10% prior</td></tr><tr><td>Detector energy scale</td><td>5% prior</td><td>×</td></tr><tr><td>PID-class norm. factors</td><td>free</td><td>×</td></tr><tr><td>Effective area scale</td><td>×</td><td>10% prior</td></tr><tr><td>Flux energy scale</td><td>×</td><td>10% prior</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Implementation of the two different systematic approaches in the SWIM workflow used in the ORCA analysis.</figcaption></figure>

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 figure 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 cross-checks and comparisons.

# 4 Combined JUNO and ORCA analysis

## 4.1 Combination strategy

As described in sections 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 $\theta_{23}$ or $\delta_{CP}$ as this experiment measures $\bar{\nu}_e \rightarrow \bar{\nu}_e$ oscillations which do not depend on those parameters [27]. On the other hand, ORCA has negligible sensitivity to $\theta_{12}$ and $\Delta m^2_{21}$ as the measured $\nu_\mu + \bar{\nu}_\mu$ oscillations happen at a much smaller $L/E$ than the one required for the development of oscillations with a frequency given by $\Delta m^2_{21}$ [30]. The four oscillation parameters that impact a single experiment are accounted for implicitly in the $\chi^2$ function computation for JUNO and ORCA following the prescription outlined below, while the remaining two oscillation parameters, $\Delta m^2_{31}$ and $\theta_{13}$, have to be considered explicitly in the joint analysis.

In JUNO, for every value of $\Delta m^2_{31}$ and $\theta_{13}$, the $\chi^2$ function is profiled using a grid with 61 uniformly spaced values of $\sin^2 \theta_{12}$ between 0.30225 and 0.31775. Profiling the $\chi^2$ over $\Delta m^2_{21}$ was shown to change the value of the $\chi^2$ by about 0.1 units or less, which is negligible in the joint analysis. Because of that fact, the value of $\Delta m^2_{21}$ is kept fixed at its assumed true value.

In ORCA, for every value of $\Delta m^2_{31}$ and $\theta_{13}$, the $\theta_{12}$ and $\Delta m^2_{21}$ values are kept fixed to their assumed true values given that ORCA has little sensitivity to those parameters, while $\theta_{23}$ and $\delta_{CP}$ are minimized without constraints. This minimization is performed twice, with the initial value of $\theta_{23}$ being located in a different octant for each minimization. Only the smallest value is kept as the global minimum of the $\chi^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 $\chi^2$ values at a fixed test value of $\Delta m^2_{31}$ and $\theta_{13}$ are calculated and summed. The true value of the oscillation parameters considered in this study are the best-fit values from ref. [48] obtained “with SK data”, unless it is explicitly stated otherwise. For added clarity, those parameters are explicitly shown in table 3. Given that neither JUNO nor ORCA are as sensitive to $\theta_{13}$ as current reactor neutrino experiments [69–71], a prior on that parameter was added to
the combined $\chi^2$ from ref. [48]. The full expression used is shown in eq. (4.1) where $\Delta m_{31}^2$ and $\theta_{13}$ are the tested values of those oscillation parameters, and the last term corresponds to the added prior with $\sin^2 \theta_{13}^{\text{GF}}$ being the current global best fit for $\sin^2 \theta_{13}$ and $\sigma_{\sin^2 \theta_{13}^{\text{GF}}}$ its uncertainty.

$$ \chi^2(\Delta m_{31}^2, \theta_{13}) = \chi_{\text{JUNO}}^2(\Delta m_{31}^2, \theta_{13}) + \chi_{\text{ORCA}}^2(\Delta m_{31}^2, \theta_{13}) + \frac{(\sin^2 \theta_{13} - \sin^2 \theta_{13}^{\text{GF}})^2}{\sigma_{\sin^2 \theta_{13}^{\text{GF}}}^2}. \quad (4.1) $$

For each set of true parameters studied, the combined $\chi^2$ from eq. (4.1) is calculated for each NMO in a 101 × 21 grid in the $(\Delta m_{31}^2, \sin^2 \theta_{13})$ space, called the $\chi^2$ map, centered around the assumed true values of the oscillation parameters and spanning uniformly a ±10% interval in $\Delta m_{31}^2$ from the central value and a ±6% interval in $\sin^2 \theta_{13}$ from the central value. More explicitly, when assuming true normal ordering with the best-fit values from ref. [48], the tested values of $\Delta m_{31}^2$ in the grid will run from $-2.78080 \times 10^{-3} \text{ eV}^2$ to $-2.27520 \times 10^{-3} \text{ eV}^2$ and from $2.27520 \times 10^{-3} \text{ eV}^2$ to $2.78080 \times 10^{-3} \text{ eV}^2$ with step of $0.00506 \times 10^{-3} \text{ eV}^2$, and those of $\sin^2 \theta_{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 section 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 $\chi^2$ map above, for each set of true oscillation parameters tested, the NMO sensitivity is determined by calculating the $\overline{\Delta\chi^2} = \chi_{WO}^2 - \chi_{TO}^2$, where $\chi_{WO}^2 (\chi_{TO}^2)$ is the minimum value of $\chi^2$ in the $\chi^2$ map in the wrong (true) ordering region of the map. The $\overline{\Delta\chi^2}$ is then converted into a median sensitivity $S(\sigma) = \sqrt{\overline{\Delta\chi^2}}$ [72]. The same procedure is also used separately for ORCA and JUNO to obtain the corresponding non-combined sensitivities, computed for each experiment alone.

The $\overline{\Delta\chi^2}$ notation is used above rather than $\Delta\chi^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.

Several studies have reported that this approximation may not be valid for neutrino oscillation studies [73, 74]. Therefore, a Feldman-Cousins procedure should ultimately be applied when interpreting the experimental results in the future. In ref. [59], the Asimov data-set approach was compared with the full Feldman-Cousins correction and observed that the Asimov approach led to an underestimation of the order of $1\sigma$ of the sensitivity for ORCA alone. A full statistical treatment of the JUNO-ORCA combination, beyond the median Asimov sensitivity, may unveil similar effects than those quantified in ref. [59], but this is beyond the scope of the sensitivity study presented here.

## 4.2 Results

Figure 4 depicts the profile $\overline{\Delta\chi^2}$ scan on the test values of $\Delta m_{31}^2$ 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
<table><caption>**Table 3.** Global best-fit values for the oscillation parameters (from ref. [48]) 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.</caption><thead><tr><th>Parameter</th><th>Normal Ordering</th><th>Inverted Ordering</th></tr></thead><tbody><tr><td>sin<sup>2</sup> θ<sub>23</sub></td><td>0.563</td><td>0.565</td></tr><tr><td>sin<sup>2</sup> θ<sub>13</sub></td><td>0.02237<sup>+0.00066</sup><sub>-0.00065</sub></td><td>0.02259 ± 0.00065</td></tr><tr><td>Δm<sup>2</sup><sub>31</sub></td><td>2.528 × 10<sup>-3</sup> eV<sup>2</sup></td><td>-2.435 × 10<sup>-3</sup> eV<sup>2</sup></td></tr><tr><td>δ<sub>CP</sub></td><td>221°</td><td>282°</td></tr><tr><td>sin<sup>2</sup> θ<sub>12</sub></td><td>0.310</td><td></td></tr><tr><td>Δm<sup>2</sup><sub>21</sub></td><td>7.39 × 10<sup>-5</sup> eV<sup>2</sup></td><td></td></tr></tbody></table>

values at the true $\Delta m_{31}^2$ and their $\overline{\Delta\chi^2}$ minima yield zero. However, when assuming wrong ordering on the fit, the minima of $\overline{\Delta\chi^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).

Figure 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 $\Delta m_{31}^2$ 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 ($\bar{\nu}_e$ for JUNO, $\nu_\mu + \bar{\nu}_\mu$ for ORCA). Due to this difference the effective oscillation frequency will be a different combination of the various $\Delta m_{ij}^2$ for each experiment [33, 34]. Since the combination requires a single resulting $\Delta m_{31}^2$ best fit, this tension together with strong constraints in $\Delta m_{31}^2$ from both experiments, and particularly from JUNO, provides the synergy effect in which the combined $\overline{\Delta\chi^2}$ minimum is enhanced to a higher value than simply adding the $\overline{\Delta\chi^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 table 4, 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 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 figure 5 assuming that the two experiments start at the same time. JUNO alone would need 10 years of operation to reach 3σ of NMO sensitivity, in the 8 reactor cores scenario. 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
<figure><img src="image_5.png" /><figcaption>**Figure 4.** $\overline{\Delta\chi^2}$ profile for only JUNO (red), only ORCA (blue), and the combination of JUNO and ORCA (green) as a function of test values of $\Delta m_{31}^2$ for 6 years of data taking assuming baseline (solid) or optimistic (dashed) systematics.</figcaption></figure>

<table><caption>**Table 4.** 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.</caption><thead><tr><td>True NMO</td><td>JUNO, 8 cores</td><td>ORCA</td><td>Simple Sum</td><td>Combination</td></tr></thead><tbody><tr><td>NO</td><td>2.3σ</td><td>6.5σ</td><td>6.9σ</td><td>7.8σ</td></tr><tr><td>IO</td><td>2.4σ</td><td>3.6σ</td><td>4.3σ</td><td>5.1σ</td></tr></tbody></table>

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 section 3.2, the ORCA analysis is also performed using a set of systematics similar to those of ref. [37], as a cross-check for an optimistic approach. Figure 4 and figure 5 show that both the optimistic and the baseline systematics give a very similar
<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

$\Delta\chi^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 $\Delta m_{31}^2$, as shown in figure 4, which causes the combination to reach sensitivities that are 1–2$\sigma$ 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 implementation at flux level (optimistic) is indeed highly correlated to a change in normalization of the flux because of its power-law behavior. On the contrary, the energy scale applied as a shift in all detector response functions (baseline) produces an effect similar to shifting the observed position of the first oscillation maximum, which mimics a variation of $\Delta m_{31}^2$. As a result, the energy scale applied at the detector response (baseline) is more strongly correlated with $\Delta m_{31}^2$ compared to the energy scale at the unoscillated flux (optimistic).

# 5 Further sensitivity studies

## 5.1 Impact of energy resolution in JUNO and 10 reactor cores scenario

One of the most challenging design specifications of JUNO is the required energy resolution of the central detector. Reaching a level of about 3%/$\sqrt{E/\text{MeV}}$ is essential for JUNO to be able to reach a 3$\sigma$ sensitivity to determine the neutrino mass ordering by itself. In this sense, if the energy resolution worsens to 3.5%/$\sqrt{E/\text{MeV}}$, the required time to reach a 3$\sigma$ sensitivity would increase by more than a factor of 2 [27]. A significant amount of effort has been made within the JUNO Collaboration to reach this goal of 3%/$\sqrt{E/\text{MeV}}$, and a description of how to get there using a data-driven approach relying on calibration data has been presented in ref. [75], where a 3.02%/$\sqrt{E/\text{MeV}}$ energy resolution has been achieved, with a worsening of this energy resolution to 3.12%/$\sqrt{E/\text{MeV}}$ after considering
<figure><img src="image_7.png" /><figcaption>Figure 6. 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 $\pm 0.5\%/\sqrt{E/\text{MeV}}$.</figcaption></figure>

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 $\pm 0.5\%/\sqrt{E/\text{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 figure 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 $\Delta m_{31}^2$ 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 $\Delta m_{31}^2$, while the best-fit value of $\Delta m_{31}^2$ for each ordering would not change significantly. Therefore, the tension of the $\Delta m_{31}^2$ best fit between JUNO and ORCA remains, which preserves the high sensitivity of the analysis.

As discussed in section 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\sigma$ about 3 years earlier, as shown in figure 7. In combination with ORCA however, there is negligible impact to the time required for the combined sensitivity to reach $5\sigma$ 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
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

ORCA best-fit values for $\Delta m_{31}^2$ in the wrong ordering scenario, rather than due to the precision of each experiment to measure the neutrino mass ordering separately.

## 5.2 Dependence on $\Delta m_{31}^2$ and $\theta_{23}$

This section presents the dependence of the analysis on the true value of the oscillation parameters, focusing particularly on $\theta_{23}$ and $\Delta m_{31}^2$. Those parameters are chosen because the true value of $\theta_{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 $\Delta m_{31}^2$ measurement, and therefore it is essential to ensure that such a boost is valid for any true value of $\Delta m_{31}^2$.

Figure 8 shows the dependence of the NMO sensitivity on the true value of $\theta_{23}$ for 6 years of data taking. As mentioned in section 4.1, JUNO has no sensitivity to $\theta_{23}$. In the case of ORCA however, the NMO sensitivity potential depends strongly on the true value of $\theta_{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σ depending on the true value of $\theta_{23}$ and true NMO, and only reaches a 5σ sensitivity for true NO with $\theta_{23}$ in the second octant. The combination curve also follows a similar $\theta_{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 $\theta_{23}$ and of the true NMO.

It is worth noting here that the current global best-fit value of $\theta_{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 sections 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.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** NMO sensitivity as a function of the true $\theta_{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. [48]).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** NMO sensitivity as a function of the true $\Delta m^2_{31}$ 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. [48]).</figcaption></figure>

Figure 9 illustrates the dependence of NMO sensitivity on the true value of $\Delta m^2_{31}$. Both JUNO and ORCA standalone sensitivities depict a slight dependence on the true value of $\Delta m^2_{31}$ with the opposite slope for each experiment. The combination is also quite flat with respect to the true $\Delta m^2_{31}$, reaching a significance of $8\sigma$ in the case of NO and $5\sigma$ in the case of IO. The effect of the boost described previously, relying on the difference between the wrong ordering measurement of $\Delta m^2_{31}$, is preserved over the whole $\Delta m^2_{31}$ range.
# 6 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 $\Delta m_{31}^2$ 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\sigma$ determination of the neutrino mass ordering for any value of the oscillation parameters. In all cases, a $5\sigma$ 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\sigma$ sensitivity due to the uncertainty on the $\theta_{23}$ value. In the favorable case of true normal ordering and $\theta_{23}$ in the second octant, a $5\sigma$ 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. [36, 37]. 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 $\Delta m_{31}^2$ 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\sigma$ 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 $\Delta m_{31}^2$ 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\sigma$ 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.

## Acknowledgments

The authors acknowledge the financial support of the funding agencies: Agence Nationale de la Recherche (contract ANR-15-CE31-0020), Centre National de la Recherche Scientifique (CNRS), Commission Européenne (FEDER fund and Marie Curie Program), Institut Universitaire de France (IUF), LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001), Shota Rustaveli National Science Foundation of Georgia (SRNSFG, FR-18-1268), Georgia; 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 AF-13, Morocco; Nederlandse organisatie voor Wetenschappelijk Onderzoek (NWO), the Netherlands; The National Science Centre, Poland (2015/18/E/ST2/00758); National Authority for Scientific Research (ANCS), Romania; Ministerio de Ciencia, Innovación, Investigación y Universidades (MCIU): Programa Estatal de Generación de Conocimiento (refs. PGC2018-096663-B-C41, -A-C42, -B-C43, -B-C44) (MCIU/FEDER), Generalitat Valenciana: Prometeo (PROMETEO/2020/019), Grisolía (ref. GRISOLIA/2018/119) and GenT (refs. CIDEGENT/2018/034, /2019/043, /2020/049) programs, Junta de Andalucía (ref. A-FQM-053-UGR18), La Caixa Foundation (ref. LCF/BQ/IN17/11620019), EU: MSC program (ref. 101025085), Spain.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

## References

[1] M.C. Gonzalez-Garcia and Y. Nir, *Neutrino masses and mixing: evidence and implications*, *Rev. Mod. Phys.* **75** (2003) 345 [hep-ph/0202058] [INSPIRE].

[2] M. Fukugita and T. Yanagida, *Physics of neutrinos and applications in astrophysics*, Springer, Berlin Germany (2003).

[3] R.N. Mohapatra and P.B. Pal, *Massive neutrinos in physics and astrophysics*, World Scientific, Singapore (2004).

[4] SUPER-KAMIOKANDE collaboration, *Evidence for oscillation of atmospheric neutrinos*, *Phys. Rev. Lett.* **81** (1998) 1562 [hep-ex/9807003] [INSPIRE].

[5] SNO collaboration, *Direct evidence for neutrino flavor transformation from neutral current interactions in the Sudbury Neutrino Observatory*, *Phys. Rev. Lett.* **89** (2002) 011301 [nucl-ex/0204008] [INSPIRE].

[6] KAMLAND collaboration, *First results from KamLAND: Evidence for reactor anti-neutrino disappearance*, *Phys. Rev. Lett.* **90** (2003) 021802 [hep-ex/0212021] [INSPIRE].
[7] PARTICLE DATA GROUP collaboration, Review of particle physics, PTEP 2020 (2020) 083C01 [INSPIRE].

[8] I. Esteban, M.C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, The fate of hints: updated global analysis of three-flavor neutrino oscillations, JHEP 09 (2020) 178 [arXiv:2007.14792] [INSPIRE].

[9] P.F. de Salas et al., 2020 global reassessment of the neutrino oscillation picture, JHEP 02 (2021) 071 [arXiv:2006.11237] [INSPIRE].

[10] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri and A. Palazzo, Global constraints on absolute neutrino masses and their ordering, Phys. Rev. D 95 (2017) 096014 [Addendum ibid. 101 (2020) 116013] [arXiv:2003.08511] [INSPIRE].

[11] V. Barger, D. Marfatia and K. Whisnant, Breaking eight fold degeneracies in neutrino CP-violation, mixing, and mass hierarchy, Phys. Rev. D 65 (2002) 073023 [hep-ph/0112119] [INSPIRE].

[12] S. Hannestad and T. Schwetz, Cosmology and the neutrino mass ordering, JCAP 11 (2016) 035 [arXiv:1606.04691] [INSPIRE].

[13] E. Di Valentino, S. Gariazzo and O. Mena, Most constraining cosmological neutrino mass bounds, Phys. Rev. D 104 (2021) 083504 [arXiv:2106.15267] [INSPIRE].

[14] M.J. Dolinski, A.W.P. Poon and W. Rodejohann, Neutrinoless double-beta decay: status and prospects, Ann. Rev. Nucl. Part. Sci. 69 (2019) 219 [arXiv:1902.04097] [INSPIRE].

[15] C.H. Albright and M.-C. Chen, Model predictions for neutrino oscillation parameters, Phys. Rev. D 74 (2006) 113006 [hep-ph/0608137] [INSPIRE].

[16] S.F. King, Models of neutrino mass, mixing and CP-violation, J. Phys. G 42 (2015) 123001 [arXiv:1510.02091] [INSPIRE].

[17] P.F. De Salas, S. Gariazzo, O. Mena, C.A. Ternes and M. Tórtola, Neutrino mass ordering from oscillations and beyond: 2018 status and future prospects, Front. Astron. Space Sci. 5 (2018) 36 [arXiv:1806.11051] [INSPIRE].

[18] T2K collaboration, Search for electron antineutrino appearance in a long-baseline muon antineutrino beam, Phys. Rev. Lett. 124 (2020) 161802 [arXiv:1911.07283] [INSPIRE].

[19] NOvA collaboration, First measurement of neutrino oscillation parameters using neutrinos and antineutrinos by NOvA, Phys. Rev. Lett. 123 (2019) 151803 [arXiv:1906.04907] [INSPIRE].

[20] SUPER-KAMIOKANDE collaboration, Atmospheric neutrino oscillation analysis with external constraints in Super-Kamiokande I-IV, Phys. Rev. D 97 (2018) 072001 [arXiv:1710.09126] [INSPIRE].

[21] ICECUBE collaboration, Determining neutrino oscillation parameters from atmospheric muon neutrino disappearance with three years of IceCube DeepCore data, Phys. Rev. D 91 (2015) 072004 [arXiv:1410.7227] [INSPIRE].

[22] T2K collaboration, Improved constraints on neutrino mixing from the T2K experiment with $3.13 \times 10^{21}$ protons on target, Phys. Rev. D 103 (2021) 112008 [arXiv:2101.03779] [INSPIRE].

[23] NOvA collaboration, Recent three-flavor neutrino oscillation results from the NOvA experiment, J. Phys. Conf. Ser. 1690 (2020) 012172 [INSPIRE].
[24] DUNE collaboration, Long-Baseline Neutrino Facility (LBNF) and Deep Underground Neutrino Experiment (DUNE): conceptual design report, volume 2: the physics program for DUNE at LBNF, arXiv:1512.06148 [INSPIRE].

[25] HYPER-KAMIOKANDE PROTO- collaboration, Physics potential of a long-baseline neutrino oscillation experiment using a J-PARC neutrino beam and Hyper-Kamiokande, PTEP 2015 (2015) 053C02 [arXiv:1502.05199] [INSPIRE].

[26] HYPER-KAMIOKANDE collaboration, Physics potentials with the second Hyper-Kamiokande detector in Korea, PTEP 2018 (2018) 063C01 [arXiv:1611.06118] [INSPIRE].

[27] JUNO collaboration, Neutrino physics with JUNO, J. Phys. G 43 (2016) 030401 [arXiv:1507.05613] [INSPIRE].

[28] JUNO collaboration, JUNO conceptual design report, arXiv:1508.07166 [INSPIRE].

[29] JUNO collaboration, JUNO physics and detector, arXiv:2104.02565 [INSPIRE].

[30] KM3NET collaboration, Letter of intent for KM3NeT 2.0, J. Phys. G 43 (2016) 084001 [arXiv:1601.07459] [INSPIRE].

[31] Y.-F. Li, Y. Wang and Z.-z. Xing, Terrestrial matter effects on reactor antineutrino oscillations at JUNO or RENO-50: how small is small?, Chin. Phys. C 40 (2016) 091001 [arXiv:1605.00900] [INSPIRE].

[32] KM3NET collaboration, First neutrino oscillation measurement in KM3NeT/ORCA, PoS ICRC2021 (2021) 1123 [INSPIRE].

[33] H. Nunokawa, S.J. Parke and R. Zukanovich Funchal, Another possible way to determine the neutrino mass hierarchy, Phys. Rev. D 72 (2005) 013009 [hep-ph/0503283] [INSPIRE].

[34] A. de Gouvêa, J. Jenkins and B. Kayser, Neutrino mass hierarchy, vacuum oscillations, and vanishing $|U(e3)|$, Phys. Rev. D 71 (2005) 113009 [hep-ph/0503079] [INSPIRE].

[35] ICECUBE collaboration, PINGU: a vision for neutrino and particle physics at the South Pole, J. Phys. G 44 (2017) 054006 [arXiv:1607.02671] [INSPIRE].

[36] M. Blennow and T. Schwetz, Determination of the neutrino mass ordering by combining PINGU and Daya Bay II, JHEP 09 (2013) 089 [arXiv:1306.3988] [INSPIRE].

[37] ICECUBE-GEN2 collaboration, Combined sensitivity to the neutrino mass ordering with JUNO, the IceCube Upgrade, and PINGU, Phys. Rev. D 101 (2020) 032006 [arXiv:1911.06745] [INSPIRE].

[38] KM3NET collaboration, Determining the neutrino mass ordering and oscillation parameters with KM3NeT/ORCA, Eur. Phys. J. C 82 (2022) 26 [arXiv:2103.09885] [INSPIRE].

[39] JUNO collaboration, TAO conceptual design report: a precision measurement of the reactor antineutrino spectrum with sub-percent energy resolution, arXiv:2005.08745 [INSPIRE].

[40] D.A. Dwyer and T.J. Langford, Spectral structure of electron antineutrinos from nuclear reactors, Phys. Rev. Lett. 114 (2015) 012502 [arXiv:1407.1281] [INSPIRE].

[41] P. Vogel and J.F. Beacom, Angular distribution of neutron inverse beta decay, $\bar{\nu}_e + p \rightarrow e^+ + n$, Phys. Rev. D 60 (1999) 053003 [hep-ph/9903554] [INSPIRE].

[42] F. Von Feilitzsch, A.A. Hahn and K. Schreckenbach, Experimental beta spectra from Pu-239 and U-235 thermal neutron fission products and their correlated anti-neutrinos spectra, Phys. Lett. B 118 (1982) 162 [INSPIRE].
[43] K. Schreckenbach, G. Colvin, W. Gelletly and F. Von Feilitzsch, *Determination of the anti-neutrino spectrum from U-235 thermal neutron fission products up to 9.5-MeV*, Phys. Lett. B **160** (1985) 325 [INSPIRE].

[44] A.A. Hahn, K. Schreckenbach, G. Colvin, B. Krusche, W. Gelletly and F. Von Feilitzsch, *Anti-neutrino spectra from $^{241}$Pu and $^{239}$Pu thermal neutron fission products*, Phys. Lett. B **218** (1989) 365 [INSPIRE].

[45] DAYA BAY collaboration, *Improved measurement of the reactor antineutrino flux at daya bay*, Phys. Rev. D **100** (2019) 052004 [arXiv:1808.10836] [INSPIRE].

[46] DAYA BAY collaboration, *Spectral measurement of electron antineutrino oscillation amplitude and frequency at Daya Bay*, Phys. Rev. Lett. **112** (2014) 061801 [arXiv:1310.6732] [INSPIRE].

[47] V. Kopeikin, L. Mikaelyan and V. Sinev, *Reactor as a source of antineutrinos: thermal fission energy*, Phys. Atom. Nucl. **67** (2004) 1892 [hep-ph/0410100] [INSPIRE].

[48] I. Esteban, M.C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, *Global analysis of three-flavour neutrino oscillations: synergies and tensions in the determination of $\theta_{23}$, $\delta_{CP}$, and the mass ordering*, JHEP **01** (2019) 106 [arXiv:1811.05487] [INSPIRE].

[49] S.T. Petcov and M. Piai, *The LMA MSW solution of the solar neutrino problem, inverted neutrino mass hierarchy and reactor neutrino experiments*, Phys. Lett. B **533** (2002) 94 [hep-ph/0112074] [INSPIRE].

[50] G. Cowan, K. Cranmer, E. Gross and O. Vitells, *Asymptotic formulae for likelihood-based tests of new physics*, Eur. Phys. J. C **71** (2011) 1554 [Erratum ibid. **73** (2013) 2501] [arXiv:1007.1727] [INSPIRE].

[51] J.A. Formaggio and G.P. Zeller, *From eV to EeV: neutrino cross sections across energy scales*, Rev. Mod. Phys. **84** (2012) 1307 [arXiv:1305.7513] [INSPIRE].

[52] KM3NET collaboration, *gSeaGen: the KM3NeT GENIE-based code for neutrino telescopes*, Comput. Phys. Commun. **256** (2020) 107477 [arXiv:2003.14040] [INSPIRE].

[53] C. Andreopoulos et al., *The GENIE neutrino monte carlo generator*, Nucl. Instrum. Meth. A **614** (2010) 87 [arXiv:0905.2517] [INSPIRE].

[54] A.G. Tsirigotis, A. Leisos and S.E. Tzamarias, *HOU Reconstruction & Simulation (HOURS): a complete simulation and reconstruction package for very large volume underwater neutrino telescopes*, Nucl. Instrum. Meth. A **626-627** (2011) S185 [INSPIRE].

[55] GEANT4 collaboration, *GEANT4 — A simulation toolkit*, Nucl. Instrum. Meth. A **506** (2003) 250 [INSPIRE].

[56] KM3NET collaboration, *Dependence of atmospheric muon flux on seawater depth measured with the first KM3NeT detection units: The KM3NeT Collaboration*, Eur. Phys. J. C **80** (2020) 99 [arXiv:1906.02704] [INSPIRE].

[57] J. Hofestädt, *Measuring the neutrino mass hierarchy with the future KM3NeT/ORCA detector*, Ph.D. thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Erlangen, Germany (2017).

[58] L. Quinn, *Neutrino Mass Hierarchy Determination with KM3NeT/ORCA*, Ph.D. thesis, Aix-Marseille University, France (2018).
[59] S. Bourret, Neutrino oscillations and Earth tomography with KM3NeT-ORCA, Ph.D. thesis, APC, Paris, France (2018).

[60] M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara and S. Midorikawa, Atmospheric neutrino flux calculation using the NRLMSISE-00 atmospheric model, Phys. Rev. D 92 (2015) 023004 [arXiv:1502.03916] [INSPIRE].

[61] J.A.B. Coelho et al., OscProb Neutrino Oscillation Calculator, https://github.com/joaoabcoelho/OscProb.

[62] A.M. Dziewonski and D.L. Anderson, Preliminary reference Earth model, Phys. Earth Planet. Interiors 25 (1981) 297.

[63] R.J. Barlow and C. Beeston, Fitting using finite Monte Carlo samples, Comput. Phys. Commun. 77 (1993) 219 [INSPIRE].

[64] D. Casadei, Estimating the selection efficiency, 2012 JINST 7 P08021 [arXiv:0908.0130] [INSPIRE].

[65] M. Paterno, Calculating efficiencies and their uncertainties, [INSPIRE].

[66] G.D. Barr, T.K. Gaisser, S. Robbins and T. Stanev, Uncertainties in atmospheric neutrino fluxes, Phys. Rev. D 74 (2006) 094009 [astro-ph/0611266] [INSPIRE].

[67] ICECUBE collaboration, Measurement of atmospheric tau neutrino appearance with IceCube DeepCore, Phys. Rev. D 99 (2019) 032007 [arXiv:1901.05366] [INSPIRE].

[68] (ICECUBE COLLABORATION)*, ICECUBE collaboration, All-flavor constraints on nonstandard neutrino interactions and generalized matter potential with three years of IceCube DeepCore data, Phys. Rev. D 104 (2021) 072006 [arXiv:2106.07755] [INSPIRE].

[69] DAYA BAY collaboration, Measurement of the electron antineutrino oscillation with 1958 days of operation at Daya Bay, Phys. Rev. Lett. 121 (2018) 241805 [arXiv:1809.02261] [INSPIRE].

[70] RENO collaboration, Measurement of reactor antineutrino oscillation amplitude and frequency at RENO, Phys. Rev. Lett. 121 (2018) 201801 [arXiv:1806.00248] [INSPIRE].

[71] DOUBLE CHOOZ collaboration, Double CHOOZ $\theta_{13}$ measurement via total neutron capture detection, Nature Phys. 16 (2020) 558 [arXiv:1901.09445] [INSPIRE].

[72] S.S. Wilks, The large-sample distribution of the likelihood ratio for testing composite hypotheses, Annals Math. Statist. 9 (1938) 60 [INSPIRE].

[73] E. Ciuffoli, J. Evslin and X. Zhang, Confidence in a neutrino mass hierarchy determination, JHEP 01 (2014) 095 [arXiv:1305.5150] [INSPIRE].

[74] X. Qian, A. Tan, W. Wang, J.J. Ling, R.D. McKeown and C. Zhang, Statistical evaluation of experimental determinations of neutrino mass hierarchy, Phys. Rev. D 86 (2012) 113011 [arXiv:1210.3651] [INSPIRE].

[75] JUNO collaboration, Calibration strategy of the JUNO experiment, JHEP 03 (2021) 004 [arXiv:2011.06405] [INSPIRE].