Main

The standard model of particle physics, extended to include neutrino mass, describes three-flavour eigenstates of neutrinos (νe, νμ, ντ) that are related to three mass eigenstates (ν1, ν2, ν3) by the 3 × 3 complex Pontecorvo–Maki–Nakagawa–Sakata unitary mixing matrix UPMNS (refs. 15,16,17). This mixing, together with non-zero neutrino mass, allows for the phenomenon of neutrino oscillation, in which, during propagation, the flavour content of a neutrino evolves at a rate that depends on neutrino mass-squared splittings (\(\Delta {m}_{ij}^{2}\equiv {m}_{i}^{2}-{m}_{j}^{2}\)) and the UPMNS matrix elements. Apart from these oscillation parameters, the rate depends on neutrino energy Eν and neutrino propagation distance L (baseline). Although experiments studying this process in recent decades have provided insights into the details of neutrino masses and mixings12, many open questions remain.

The mixing matrix UPMNS is typically parameterized in terms of three mixing angles (θ12, θ13, θ23) and at least one complex phase δCP (ref. 12). It is unknown whether sin δCP is non-zero; if it is, neutrinos—and thus leptons—violate charge-parity (CP) symmetry and thereby provide a source of matter–antimatter asymmetry in nature17, which is of great interest given the connection between CP violation and the unexplained matter dominance in the Universe7,8,9,10,11. Separately, oscillation experiments have established that the mass-squared splitting \(\Delta {m}_{32}^{2}\) is roughly 30 times larger in magnitude than \(\Delta {m}_{21}^{2}\), but the sign of the former is at present unknown. That is, ν3 may be heavier or lighter than the ν1–ν2 pair, with these two options termed, respectively, the normal (\(\Delta {m}_{32}^{2} > 0\)) and inverted (\(\Delta {m}_{32}^{2} < 0\)) mass orderings. Knowledge of the mass ordering can constrain experimental searches and theory development in a wide range of physics, including absolute neutrino mass measurements18, neutrinoless double beta decay searches to investigate the nature of neutrino mass19, models of supernova explosion and detection20,21, and the cosmological evolution evidenced in cosmic microwave background and large-scale structure measurements22. For the mixing angles, current data suggest θ23 is near 45°, a notable value hinting at a μ/τ flavour symmetry17. Improved precision on this and other mixing angles is essential for gaining a clearer view of flavour mixing and to probe the validity of the three-flavour model.

Long-baseline accelerator neutrino oscillation experiments are well suited to address the above questions. In these, a high-intensity neutrino beam enriched in muon neutrinos (νμ) or muon antineutrinos (\({\bar{{\rm{\nu }}}}_{{\rm{\mu }}}\)) is produced at a particle accelerator and directed through the crust of Earth towards a massive far detector located hundreds of kilometres away. Note that the word ‘neutrino’ is used to mean both neutrino and antineutrino unless stated otherwise. The far detector measures the event rates of νμ and νe—the latter primarily from νμ → νe oscillation—as a function of neutrino energy, from which the oscillation parameters above can be determined. These experiments use near detectors, sited a short distance from the beam source, in which oscillation effects are negligible and a very high neutrino event rate can be measured. The near detectors provide vital control measurements that substantially mitigate large systematic uncertainties in the initial neutrino flux, neutrino-on-nucleus interaction cross-sections and in some cases detector response (for example, energy reconstruction and event selection efficiencies).

Two such experiments are in operation today, T2K and NOvA. Each experiment uses a narrow-band off-axis beam23,24, whose peak energy is near the first oscillation maximum, \({\sin }^{2}\left(\frac{\Delta {m}_{32}^{2}L}{4E}\right)\approx 1\), at the far detector. Note that natural units, where ħ = c = 1, are used throughout. T2K uses an approximately 0.6 GeV neutrino beam from J-PARC in Tokai, Japan, and the 50-kt Super-Kamiokande water Cherenkov detector for its far detector located 295 km away25. In the United States, an approximately 2 GeV beam of NOvA is produced at Fermilab near Chicago, and the 14-kt tracking calorimeter far detector is located 810 km away in northern Minnesota26. Further details on the designs of NOvA and T2K and on long-baseline experiments can be found in the Methods and refs. 25,26,27.

We report here a combined analysis of the datasets from T2K and NOvA previously analysed independently in refs. 13,14. This combination takes advantage of marked complementarity in the sensitivities of the two experiments to the oscillation parameters. In particular, the νμ → νe oscillation probability is a function of (among other things) both δCP and the neutrino mass ordering, and these two effects must be teased apart.

Figure 1 shows the complementarity between the experiments in a simplified case. Sets of oval curves indicate the energy-integrated total νe and \({\bar{\nu }}_{{\rm{e}}}\) event counts expected in the far detectors under various mass ordering and δCP scenarios, with other oscillation parameters held fixed. The measured event counts in NOvA and T2K are shown as black points with error bars.

Fig. 1: The impact of mass ordering and δCP on event rates.
figure 1

a,b, A bi-event plot that shows experimental sensitivity to neutrino mass ordering and δCP, with panels representing the NOvA (a) and T2K (b) cases. Black points with 1σ Poisson statistical error bars show the total number of νe and \({\bar{\nu }}_{{\rm{e}}}\) candidates selected in the far detectors. The oval parametric curves trace out predicted numbers of events under the normal (blue) or inverted (orange) mass ordering assumption as the parameter δCP varies from −π to π. Four specific δCP values are labelled for reference. All other oscillation parameters are kept fixed in this graphic, set to their most probable values from the joint analysis (Extended Data Table 3).

Full size image

As shown in Fig. 1a, there is stronger separation between the mass ordering ovals for NOvA, because of higher beam energies, but as the NOvA data lie near the overlap of the ellipses, there can be ambiguity as to which ordering is correct and (in a correlated way) which values of δCP are preferred. By contrast, T2K has less sensitivity to the mass ordering, but points with similar values of δCP in each hierarchy sit close to one another, and the data lie closest to \({\delta }_{{\rm{CP}}}=-\,\frac{\pi }{2}\), regardless of mass ordering. Combining these datasets can provide simultaneous mass ordering and δCP information with substantially less ambiguity, maximizing the use of current data and informing data-taking strategies for current and future experiments.

This discussion points to a more general observation that the oscillation parameters of interest represent a highly correlated multidimensional space. The analysis reported here calculates a joint Bayesian posterior, using the likelihoods of the experiments defined over the full parameter space. Moreover, we use the full suite of analysis tools from both experiments: detector response models, neutrino energy estimators, near-detector measurements and systematic uncertainties, all within a unified framework for statistical inference. This level of integration is the first for accelerator neutrino experiments, to our knowledge.

The posterior calculation is based on detailed parameterized models of the neutrino flux, cross-sections and detectors that predict the binned spectra of neutrino events in each of our selected samples as a function of the oscillation parameters and a large number of nuisance parameters mostly related to systematic uncertainties in the models. A likelihood is constructed from Poisson probability terms describing the compatibility between the prediction and the observed data in bins of relevant variables. Prior probabilities are set on all parameters as detailed in the Methods.

Both T2K and NOvA have software that explores the posterior using Markov chain Monte Carlo (MCMC) methods28,29 (ARIA for NOvA30 and MaCh3 for T2K31). By containerizing32 the likelihood and prior portions of the code, we can construct and analyse the joint posterior using either of the original MCMC frameworks, in spite of the very different software environments involved. For each fitting framework, ARIA or MaCh3, the native likelihood and priors of the fitter are calculated directly, whereas the likelihood and priors of other experiments are accessed using the software container. In this way, either framework can be used, providing valuable redundancy and thus cross-checks of all statistical inferences.

Although a single set of oscillation parameters naturally applies to both experiments in the joint posterior, the treatment of the many nuisance parameters related to systematic uncertainties is more subtle. Both measurements of the oscillation parameters at present have statistical uncertainties larger than the systematic uncertainties, but the latter are not negligible. We thoroughly surveyed the flux, cross-section and detector models and their systematic uncertainties to determine whether correlations between the experiments affect the analysis at a significant level. Our conclusions from this effort are summarized in the following paragraphs.

Both T2K and NOvA use beams produced by directing accelerated protons onto graphite targets. The hadrons are charge-selected with magnetic horns: positively charged hadrons decay to produce neutrinos, and negatively charged hadrons produce antineutrinos. Many uncertainties on these beam fluxes stem from processes unrelated between the two experiments, for example, the alignment of beam components. Yet, uncertainties on the rate of hadron production in the graphite targets are substantial, and the underlying physics is the same. However, the two experiments use proton beams of different energies (30 GeV for T2K and 120 GeV for NOvA), and the external datasets used to tune the hadron production models of both experiments are different33,34,35. Moreover, the ultimate impact of flux uncertainties on far-detector predictions in NOvA is much smaller than other uncertainties. We, therefore, conclude that at current experimental exposures, the flux uncertainties of the two experiments need not be correlated.

Given the different detector technologies involved, most detector-related uncertainties are independent between the experiments. Furthermore, the very different energy estimation techniques used, combined with model tuning and uncertainty estimation using in situ calibration samples in each experiment, including for the lepton and neutron energy scales, lead to independence between the two detector uncertainty models. We conclude that there are no significant correlations in the detector models.

For neutrino-on-nucleus cross-sections, the underlying physics is the same; in many cases, the same external datasets are used by both experiments to tune and set prior uncertainties on model parameters. Thus, cross-section model correlations are expected. However, in the specific case of NOvA and T2K, the description of this physics differs markedly. The simulation packages differ36,37, the physical models implemented in them differ in many places, the parameterizations differ almost entirely, and customized tunings are necessary and applied, given the specific energies of the experiments, detector technologies and approaches to systematic uncertainty mitigation.

Proper correlations between experiments could be implemented by starting from a common cross-section model spanning different energy ranges and able to describe both the leptonic and hadronic parts of the final state. A joint description is not yet mature and is one of the focuses of the community in the years to come38. Given the differences in the models, a direct mapping of their parameters was deemed not practical at this time. Instead, we studied whether neglecting these correlations could appreciably affect our measurements of the oscillation parameters. The studies are limited to our current experimental exposures and models and would need re-evaluation if applied to any other context.

First, we assessed whether correlations between single systematic parameters in our models could have a substantial impact on our results. For each of \(\Delta {m}_{32}^{2}\), θ23 and δCP, we identified the systematic parameter in each experiment with the largest impact on the measurement of that oscillation parameter. Then, regardless of whether those two systematic parameters made physical sense to correlate, we performed fits to simulated pseudo-data with the parameters fully correlated, uncorrelated and fully anticorrelated. Details of these studies, including how we identified the most impactful parameters, are shown in the Methods. In summary, we saw no case in which the choice of correlation of individual systematic parameters significantly affected the oscillation parameter measurements.

Checking individual parameters does not rule out effects from a mix of systematic parameter variations that combine to produce a net effect that is larger and possibly more degenerate with oscillation effects, representing a potential worst-case scenario for the analyses. Rather than seeking such a set of variations, we directly identified, or in some cases constructed, single systematic parameters for each experiment that have effects similar to each oscillation parameter of interest. We then adjusted the size of the priors on these ‘nightmare’ parameters such that their impact on the measurements is comparable to that of statistical errors and therefore larger than the net effect of all our regular systematic parameters combined. These nightmare parameters were added to our nominal uncertainty models to create augmented models, allowing us to study a case in which systematic effects are comparable to statistical uncertainty. Next, we constructed simulated pseudo-datasets with the nightmare parameters increased in both experiments by one standard deviation above their prior central values. These simulated pseudo-data were then fit three times using the augmented model: once with the nightmare parameters of the experiments fully correlated (matching the pseudo-data), once fully anticorrelated and finally uncorrelated. We find that the oscillation parameter constraints extracted in the fully correlated and uncorrelated cases have negligible differences. However, the incorrect anticorrelated case yields a large bias. We expect that with even larger systematic uncertainties, differences between the correlated and uncorrelated cases would eventually become relevant. However, this study indicates that we are not in such a regime with the current exposures and systematic uncertainties (see the Methods for further results).

Given that no significant biases are seen from neglecting correlations between actual systematic parameters, and the only bias seen with the nightmare parameters comes not from neglecting a correlation but from adding an incorrect one, we choose in most cases to neglect the correlations between the systematic uncertainties of the two experiments. The one exception relates to the approximately 2% normalization uncertainties on all νe and \({\bar{\nu }}_{{\rm{e}}}\) events described in ref. 39. In this case, the uncertainties are implemented identically by T2K and NOvA, and we have correlated them.

We also perform studies in which the joint fit is tested against pseudo-data constructed with a set of discrete model variations not directly accessible using the nominal uncertainty models of the experiments. This procedure was used in the earlier independent T2K analysis14, and we include in the present analysis those model variations seen as most impactful previously. Similarly, we studied a secondary set of variations based on extrapolating the cross-section model of each experiment to the context of the other experiment. Predefined thresholds were used to establish that no substantive changes in the central values or interval widths of the oscillation parameters were seen under these tests, as described in the Methods. For all tested alternative models, all observed changes in credible intervals were within thresholds (see the Methods for further details). Each experiment continues to investigate improvements in its cross-section models, and the studies described here would warrant repeating for larger data exposures and/or updated theoretical understanding. Continued theoretical and experimental effort in this direction is important.

With the joint likelihood and systematic uncertainty model defined, we use our fitting frameworks to analyse the combined datasets of refs. 13,14, finding consistent results between the two frameworks. Unless stated otherwise, we report results using an external constraint on θ13 (named the ‘reactor constraint’ below) and external constraints on \(\Delta {m}_{21}^{2}\) and θ12. The values used for these constraints correspond to the 2020 Particle Data Group summary values40 and are given in the Methods.

We tested the goodness of fit (Methods) of our model to data using the P-value method41, both overall and for each individual sample in the far detectors. All the P-values are within an acceptable range (>0.05 after the look-elsewhere-effect adjustment described in the Methods). The overall P-value to describe all NOvA and T2K samples is 0.75 for full spectral analysis and 0.40 for rate-only analysis, marginalized over both mass orderings. Similar results were obtained without the reactor constraint and in each mass ordering. Thus, the joint oscillation model simultaneously fits T2K and NOvA data well. The P-values are also consistent with those of previous T2K-only and NOvA-only analyses.

We produce parameter estimations using the highest-posterior-density credible intervals and perform discrete hypothesis tests using the Bayes factor formalism. Conclusions related to CP conservation or violation, \(\Delta {m}_{32}^{2}\), \({\sin }^{2}{\theta }_{23}\) and mass ordering have been tested to be robust under the alternative model variations described previously. For the measured oscillation parameters, we report 1σ (68.27%) credible intervals unless noted.

We find \({\sin }^{2}{\theta }_{23}=0.5{6}_{-0.05}^{+0.03}\) without any assumptions on the ordering of the neutrino masses. The fit weakly prefers the upper octant of θ23 (\({\sin }^{2}{\theta }_{23} > 0.5\)) over the lower octant with a Bayes factor of 3.5. Removing the reactor constraint gives no statistically significant preference for either octant (Bayes factor 1.2 for the lower octant compared with the upper octant). We also find \(\Delta {m}_{32}^{2}=2.4{3}_{-0.03}^{+0.04}\times 1{0}^{-3}\,{{\rm{eV}}}^{2}\) assuming the normal ordering and \(\Delta {m}_{32}^{2}=-\,2.4{8}_{-0.04}^{+0.03}\times 1{0}^{-3}\,{{\rm{eV}}}^{2}\) assuming the inverted ordering. This is at present the smallest experimental uncertainty on \(| \Delta {m}_{32}^{2}| \) (Fig. 2), to our knowledge. This conclusion also applies when the reactor constraint is replaced by a flat prior.

Fig. 2: Experimental measurements of \(| {\boldsymbol{\Delta }}{{\bf{m}}}_{{\bf{32}}}^{{\bf{2}}}| \).
figure 2

The measurements assume the inverted ordering preferred by this analysis. Sources for the results from top to bottom, starting with the second line, are as follows: refs. 13,14,43,44,45,46,47,48,49. The normal ordering case is available in Extended Data Fig. 9.

Full size image

There is no statistically significant preference obtained for either of the mass orderings, with a Bayes factor of 1.3 in favour of the inverted ordering with reactor θ13 constraint and 2.5 without reactor θ13 constraint. Although the two experiments individually prefer the normal ordering, the values of other oscillation parameters are more consistent in the inverted ordering, leading to a different ordering preference in the joint fit, although still not statistically significant. The effect on mass ordering preference when additionally incorporating reactor \(\Delta {m}_{32}^{2}\) measurements is discussed in the Methods.

With no assumption on the true mass ordering, we find the 1σ credible interval on δCP to contain [−0.81π, −0.26π] with the highest posterior probability value being −0.47π. We also find that values of δCP around +π/2, an extremum of sin δCP, are outside our 3σ (99.73%) credible intervals, which also holds for either mass ordering separately. Figure 3 shows the joint fit result compared with the individual measurements of NOvA and T2K in the \({\sin }^{2}{\theta }_{23}-{\delta }_{{\rm{CP}}}\) plane, as well as one-dimensional (1D) uniformly binned posterior probability distributions for both mass ordering cases. Assuming the normal ordering, the joint analysis allows a wide range of δCP values, giving a 3σ credible interval of δCP ∈ [−1.38π, 0.30π]. In the case of the inverted ordering δCP ∈ [−0.92π, −0.04π], excluding 56% of the parameter space, the CP-conserving values of δCP = 0 and π are outside the 3σ credible interval. A consistent picture is seen when analysing the Jarlskog invariant, JCP (ref. 42), which is a parametrization-independent measure of CP violation. The CP-conserving value of JCP = 0 falls outside the 3σ credible interval for the inverted ordering, and the above statements are true whether the prior used is uniform in δCP or sin δCP (Fig. 4). This analysis, therefore, provides evidence for CP violation in the lepton sector if the inverted ordering is assumed to be true. However, we do not see a significant preference at present for either mass ordering. Future mass ordering measurements will, therefore, influence the interpretation of these results. See the Methods for more data projections and comparisons.

Fig. 3: Constraints on \({{\bf{\sin }}}^{2}{{\boldsymbol{\theta }}}_{23}\) and δCP.
figure 3

Marginalized posterior probabilities and 1D or 2D Bayesian credible regions of \({\sin }^{2}{\theta }_{23}\) and δCP in the case of the normal (blue, left side) and inverted (orange, right side) neutrino mass ordering with the reactor constraint applied. Shaded areas correspond to 1σ, 2σ and 3σ credible regions. a,b, The 2D panels of \({\sin }^{2}{\theta }_{23}\) vs δCP (a,b) are overlaid with 1σ credible regions from the T2K-only (dark red) and NOvA-only (dark blue) data fits assuming normal (a) and inverted ordering (b). c–f, The 1D panels show the posterior probabilities of \({\sin }^{2}{\theta }_{23}\) (c) and δCP (d) in the normal ordering, and δCP (e) and \({\sin }^{2}{\theta }_{23}\) (f) in the inverted ordering.

Full size image
Fig. 4: Constraints on the Jarlskog invariant.
figure 4

a–d, Marginalized posterior probabilities of the Jarlskog invariant, JCP, in the case of the normal (blue; a,b) and inverted (orange; c,d) neutrino mass ordering with the reactor constraint applied. The posterior distributions use prior distributions either flat in δCP (a,c) or sin δCP (b,d). Shaded areas show the 1σ, 2σ and 3σ Bayesian credible intervals.

Full size image

Methods

The NOvA experiment

The NOvA experiment measures neutrino oscillations using two detectors of functionally identical construction located along the NuMI neutrino beam50 produced at the Fermi National Accelerator Laboratory (Fermilab).

The smaller 0.3-kt near detector is located on the Fermilab campus 1 km downstream from the neutrino production target, whereas the 14-kt far detector is located 810 km away in northern Minnesota. The detectors themselves are highly segmented tracking calorimeters consisting of long PVC cells filled with a mineral-oil-based liquid scintillator. Each cell measures 6.6 cm × 3.9 cm in cross-section, runs the full height or width of the detector (15.5 m for the far detector and 3.9 m for the near detector) and is instrumented with a wavelength-shifting fibre and avalanche photodiode to detect the scintillation light produced when charged particles pass through the cell. The cells are arranged in a series of layers, each with either horizontal or vertical orientation, with the direction alternating between layers to provide three-dimensional (3D) event reconstruction. This segmented design offers the excellent muon and electron classification needed for tagging the incoming neutrino flavour. In particular, electromagnetic showers at typical NOvA energies are much larger than the detector cell widths and thus are well-imaged and distinct from many potential backgrounds. The detectors of NOvA are centred 14.6 mrad off the central axis of the NuMI beam, yielding a narrow-band neutrino beam peaked at 1.8 GeV.

As is typical for particle physics experiments, NOvA makes use of detailed simulations of beam production, neutrino interaction physics and detector response as part of the analysis. Given the matching near and far detectors, NOvA forms its oscillation-dependent predictions of the far-detector event rates directly from data using the millions of neutrino interactions recorded in the near detector. This near-to-far extrapolation process is carried out as a function of multiple kinematic and event classification variables. Uncertainties from the simulations have substantially reduced impact as they enter the oscillation fit only to the extent that they affect the mapping between expected near and far event rates, not the event rates of the individual detectors themselves. Uncertainties on the simulations are taken as the a priori uncertainties from, for instance, the external model constraints or other external data and are supplemented by additional model uncertainties in which a priori coverage was deemed unsatisfactory.

Far-detector data are fitted to the corresponding far-detector predictions to extract oscillation parameter constraints. These data are separated by beam mode (that is, neutrino- or antineutrino-dominated running) and further into \({{\rm{\nu }}}_{{\rm{\mu }}}/{\bar{{\rm{\nu }}}}_{{\rm{\mu }}}\) charged current and \({\nu }_{{\rm{e}}}/{\bar{\nu }}_{{\rm{e}}}\) charged current candidate samples using a convolutional neural network51 whose inputs are the calibrated event images recorded by the detector cells. Subsequent reconstruction of tracks and showers within each event provides kinematic information such as estimated neutrino energy. Far detector \({{\rm{\nu }}}_{{\rm{\mu }}}/{\bar{{\rm{\nu }}}}_{{\rm{\mu }}}\) samples are analysed in bins of neutrino energy and hadronic energy fraction. The \({\nu }_{{\rm{e}}}/{\bar{\nu }}_{{\rm{e}}}\) samples are analysed in bins related to event containment, event classification score and neutrino energy. More details on the analysis techniques, simulation packages, systematic uncertainties and the overall NOvA experimental design can be found in ref. 13 and the references therein.

The T2K experiment

The T2K experiment is composed of the J-PARC neutrino beam, a near site with multiple detectors and the water Cherenkov detector Super-Kamiokande (SK) as the far detector. Full details of the experiment can be found in ref. 25.

The primary detector at the near site, 280 m from the target, is a magnetized off-axis (centred at 43.6 mrad) tracking detector called ND280. While taking the data used in this analysis, ND280 consisted of a π0 detector followed by a tracker consisting of three time-projection chambers interleaved with two hydrocarbon fine-grained detectors (FGD1 and FGD2), all surrounded by an electromagnetic calorimeter. The stability and direction of the neutrino beam are monitored using the on-axis near detector INGRID.

SK is situated 295 km downstream of the neutrino production target, 43.6 mrad off-axis, and contains 50 kt of water. An inner detector (ID) using 11,129 inward-facing 20-inch photomultiplier tubes (PMTs) detects Cherenkov radiation from charged particles traversing the detector. An optically separated outer detector uses 1,885 outward-facing 8-inch PMTs to reject interactions originating outside the ID volume. SK can discriminate between electrons and muons by their Cherenkov ring profiles.

T2K uses a forward-fitting analysis strategy. First, a model that predicts the event spectra at the near and far detectors is defined and tuned to external experimental data. The predictions are generated by simulating the neutrino flux and cross-section as well as the detector response. The model, with variable parameters, is fit to the ND280 data to obtain tuned values of the parameters with uncertainties. The constrained model resulting from this near-detector fit is then used to make SK predictions, which are fit to the SK data to extract oscillation parameters. Complete details for this analysis, including model details, are in ref. 14.

T2K splits data at the near and far detectors into mutually exclusive samples defined by particle identification in each beam mode. At ND280, events are categorized into 18 samples, nine samples in each of FGD1 and FGD2. In neutrino mode, data with one negatively charged muon is split into three samples in each FGD corresponding to the number of pions (0, 1, or >1). In antineutrino mode, data are first split by whether a negatively or positively charged muon is present, and then divided by the number of pions as in the neutrino-mode data, forming six samples in each FGD. For all samples, the data are fit in a 2D space of the muon momentum and the angle between the muon and the average beam direction. The exclusive samples allow the near-detector fit to better constrain parameters related to different neutrino–nucleus interaction modes. At SK, the data are divided into three samples in neutrino mode: one-ring muon-like, one-ring electron-like and one-ring electron-like with one decay electron; in antineutrino mode, only the one-ring muon-like and one-ring electron-like samples are used. The data are binned in reconstructed neutrino energy. All electron-like samples are additionally binned in a second dimension, the angle between the reconstructed electron direction and the beam direction.

Detector systematic uncertainties are evaluated using a variety of sideband samples and calibrations, covering effects such as particle identification, particle momentum reconstruction, secondary particle interactions and fiducial volume effects.

Correlations in flux modelling

The modelling of the neutrino flux depends on many details relating to the incident proton beam, the hadron production target and the magnetic focusing horns. As these details are specific to each experiment, flux systematic uncertainties due to magnetic field variations, component alignment and other beamline properties are uncorrelated between the experiments.

The only possible correlation identified was the pion and kaon production models and the use of hadron interaction experiments to tune them52,53. In the case of NOvA, the primary data are from the NA49 experiment33, which collected thin-target (slices of the target material) data at 158 GeV c−1, which is then scaled to the NuMI beam energy. The NA61/SHINE experiment, which collected data for T2K, uses some of the same detectors and the same beamline as NA49. NA61/SHINE34,35 collected both thin-target and replica-target (a full-sized target) data for T2K at 31 GeV c−1, the J-PARC beam momentum. Checking the consistency of the NA49 and NA61/SHINE data used is difficult, as the data are collected at different beam energies.

The NOvA experiment primarily uses thin-target NA49 hadron production data to tune the particle multiplicities, reweighting interactions and particle propagation inside the target and other beamline materials. By contrast, T2K uses thin and replica-target data from NA61/SHINE to reweight the multiplicities of particles exiting the target. Given these differences in data collection and tuning methodology, and given that flux uncertainties have a suppressed influence after ND data constraints are considered, there is no expectation of significant correlations between flux systematic parameters for NOvA and T2K in the joint fit.

Correlations in detector modelling

The experiments use different detector technologies as well as strategies for forming data samples, which removes most opportunities for correlation. However, the modelling of particle propagation through the detectors derives from the same underlying physics. This propagation is called secondary interaction (SI), and the case of pion SI is noteworthy, as this process is expected to occur in both experiments, and for T2K, it is an important effect. T2K selects exclusive data samples in which a change in reconstructed pion multiplicity can cause migration between samples. By contrast, NOvA uses inclusive selections, and pion SI has minimal effect on the calorimetric energy estimation at NOvA. Thus, we do not expect significant correlations due to pion SI.

Tests of individual parameter correlations

Neutrino-on-nucleus scattering plays a central part in both experiments, but the modelling of this physics has substantial differences between the two individual analyses. These differences, together with the presence of different nuclear targets, neutrino energies and near-detector strategies, mean that direct estimation of systematic uncertainty correlations in the neutrino scattering models is highly non-trivial. As part of this analysis, we tested how significant inter-experimental systematic uncertainty correlations could be, starting by identifying the most impactful systematic uncertainties of T2K and NOvA and exploring correlations between them.

To determine an impactful systematic parameter, we carry out a fit to pseudo-data generated with all parameters at their prior values from our nominal model. Then, for each parameter in turn, we reweight all steps from the obtained MCMC chain to have a tight (‘shrunk’) prior for that parameter around a different value (‘pulled’) to that used to generate the pseudo-data and study the change in the extracted oscillation parameter intervals. This procedure mocks up the result of an external experiment, providing a strong constraint on each systematic parameter at a different value from that preferred by simulated pseudo-data. This ‘shrink and pull’ study allows for assessing the single-parameter impact on the systematic uncertainty and the estimated credible intervals of the measurement of the individual neutrino oscillation parameters.

First, we identify both the systematic parameters of NOvA and T2K with the largest impact on δCP, \({\sin }^{2}{\theta }_{23}\) and \(\Delta {m}_{32}^{2}\) in the joint fit.

For both experiments, the largest change in δCP credible interval comes from uncertainties on νe and \({\bar{\nu }}_{{\rm{e}}}\) normalizations. As discussed, these uncertainties are implemented identically in both experiments, and we have correlated them in the joint analysis. No additional interaction uncertainties in our models have any significant impact on the resulting credible intervals of δCP.

For \({\sin }^{2}{\theta }_{23}\), all the individual interaction systematic parameters have very small effects, changing the width of the 1σ interval by less than 2% when shrunk by 50% and pulled 1σ away from the nominal value. The largest change in credible interval comes from the uncertainty on the neutron visible energy for NOvA, and the two-particle two-hole (2p2h) C/O cross-section scale for T2K (2p2h C/O cross-section scale allows the 2p2h cross-section on carbon to differ from that for oxygen). For \(\Delta {m}_{32}^{2}\), all the individual interaction parameters have a negligible effect on the resulting \(\Delta {m}_{32}^{2}\) credible intervals. Hence, we widened the list of considered parameters and identified the calorimetric energy scale uncertainty of NOvA and the SK energy scale uncertainty of T2K as the most impactful for \(\Delta {m}_{32}^{2}\).

Second, despite there being no a priori reason to expect correlations between these specific parameters, we test whether or not correlating the most impactful T2K parameter with the most impactful NOvA parameter modifies oscillation parameter constraints in the joint fit in a significant way. We simulate pseudo-data to which we perform a joint fit while treating the T2K and NOvA parameters described above as either uncorrelated, fully correlated or fully anticorrelated. We repeat the study for each pair of the most impactful parameters of T2K and NOvA with respect to δCP, \({\sin }^{2}{\theta }_{23}\) and \(\Delta {m}_{32}^{2}\). In the case of \(\Delta {m}_{32}^{2}\), we further inflate the original SK energy scale uncertainty from 2% to 7% to amplify the effect. Finally, we check the extracted 1σ and 2σ credible regions for any substantial differences between the three correlation configurations. These tests are repeated for three sets of pseudo-data generated with oscillation parameter values that are T2K-like, NOvA-like and NuFit-like54, which are chosen to be close to recent data results from the respective collaborations and are given in Extended Data Table 1.

As an example, Extended Data Fig. 1 shows the results in terms of the posterior probability distributions and credible regions of the parameters of interest from the set of fits with the largest single-parameter impact on \({\sin }^{2}{\theta }_{23}\). We conclude that the choice of correlation between single parameters does not significantly change the oscillation parameter constraints derived from the current version of the joint analysis.

Nightmare parameters

As described in the main text, we study correlations in more extreme situations using the so-called nightmare parameters, which are either artificially constructed parameters or existing parameters with highly inflated uncertainties chosen to be deliberately problematic for the individual analyses. The prior uncertainties of the parameters are set so that they are comparable in impact to the statistical uncertainties on the measurements under study. We carry out this procedure separately for simulated measurements of \(\Delta {m}_{32}^{2}\) and θ23. No nightmare study was carried out for δCP because its total systematic uncertainty compared with the statistical uncertainty is much smaller than for the other two cases.

We construct pseudo-datasets with both the NOvA and T2K nightmare parameters shifted by one standard deviation from their prior values, inducing a systematic bias representing a simultaneous and coordinated shift in both NOvA and T2K data. We fit this pseudo-data while treating the NOvA and T2K nightmare parameters as either fully correlated, uncorrelated or anticorrelated. The results of the nightmare parameters correlation study are presented as 1σ credible 2D regions of \(\Delta {m}_{32}^{2}-{\sin }^{2}{\theta }_{23}\) in Extended Data Fig. 2 for both nightmare scenarios. We conclude that there is no significant difference in treating the nightmare parameters as either fully correlated (matching the pseudo-data) or uncorrelated between the experiments, whereas the incorrect anticorrelated case yields a clear bias. We note that these are not general conclusions but are specific to the T2K and NOvA analysis versions and cumulative beam exposures used here. The construction of the nightmare parameters is also not a unique choice, and other formulations of the parameters could be considered.

Out-of-model variations

As described in the main text, we use a set of discrete changes to the base cross-section model to test the robustness of our analysis. For each test, pseudo-data are generated assuming the specific model variation, and these pseudo-data are then fit either with the default analysis directly, which does not incorporate the model variation (‘out-of-model’ case) or with a modified analysis that has had its nominal event spectra altered to match the spectra expected under the varied model (‘in-model’ case). Between these two cases, we require that the width of each of the extracted oscillation parameter intervals changes by no more than 10% (representing a small ‘error on the error’) and that the centre of the interval does not move by more than 50% of the systematic uncertainty (indicating adequate systematic uncertainty coverage of the tested out-of-model variation). Furthermore, we require that taking the largest changes seen across these studies does not affect the stated conclusions on CP violation or mass ordering determination for the analysis.

Three variations were chosen to perform the out-of-model studies:

  • MINERvA 1π: this model suppresses charged current (CC) and neutral current (NC) resonant pion production at low Q2 to ensure good agreement between the MINERvA data55 and the implementation of the Rein–Seghal model in the GENIE v.2 neutrino interaction simulation software37.

  • Non quasi-elastic (non-QE): in the T2K oscillation analysis14, the ND280 data samples with a muon candidate and zero pion candidates are underpredicted by the pre-fit T2K nominal model by 10% in both FGDs, which the fit accounts for by enhancing the charged current quasi-elastic (CCQE) interaction rate. To check this large freedom does not cause bias, an alternate model is produced, in which this underprediction is attributed to only non-QE processes.

  • Pion SI: the pion SI model in the GEANT4 detector simulation toolkit 56 was replaced with the Salcedo–Oset model57 implemented in the NEUT generator36, tuned to π–A scattering data58.

We also used this process to study what happens when fitting pseudo-data constructed for both experiments using the nominal cross-section model of one or the other experiment (T2K-like and NOvA-like studies).

We show example results here for the MINERvA 1π case. Extended Data Fig. 3a,b shows the effect of this alternative model on event spectra used in the analysis. Note that not all event spectra are uniformly binned. Extended Data Figs. 3c–g and 4 compare the in-model and out-of-model fit results. No failures of our criteria are seen in any of the cases. More generally, no significant bias is seen in this joint fit for any of the model variations studied across any of the three tested sets of oscillation parameter values.

Some more recent T2K analyses45 did see criteria failures when considering an alternative nuclear model, HF-CRPA59, and as a result widened their \(\Delta {m}_{32}^{2}\) intervals. Both NOvA and T2K have independently studied the impact of the HF-CRPA model on the analyses used in this joint result, and we estimate that any potential effects in the context of this joint fit are within the thresholds set for our out-of-model variation tests.

Goodness of fit

The posterior-predictive P-value41 technique is used to determine whether a model provides a good fit to the data it is confronted with. We require that the posterior-predictive P-value to obtain the far-detector data in all samples, given the joint post-fit model, is greater than 0.05. We also check the P-values for individual far-detector samples and require that they are greater than 0.05 after allowing for the look-elsewhere effect, using the Bonferroni correction60. All the P-values from the joint fit are shown in Extended Data Table 2. All the P-values (both total and split sample by sample) are within our acceptable range (>0.05), even without taking the look-elsewhere effect into account. This means that the model used in this joint fit—that is, the systematic models of the individual experiments with a shared oscillation parameter model—fits our data well, even when looking at individual samples. The P-values are consistent with previous T2K-only and NOvA-only analyses. The P-value considering rate and shape for all T2K samples in a T2K-only fit is 0.73, whereas the P-value considering all T2K samples in the joint fit is 0.75. Similarly, the P-values for all NOvA samples are 0.56 (NOvA-only fit) and 0.64 (joint fit).

Example posterior predictions61 of the spectra for the νμ and νe subsamples of both experiments, overlaid over the observed data, are shown in Extended Data Fig. 5.

Priors

The default priors on the oscillation parameters for this analysis are as follows: flat between −π and π in δCP, flat between 0 and 1 in \({\sin }^{2}{\theta }_{23}\), flat in \(\Delta {m}_{32}^{2}\) and Gaussian with μ ± σ = (2.18 ± 0.07) × 10−2 in \({\sin }^{2}{\theta }_{13}\). Where alternate priors are used, this is stated in the text.

This analysis is not sensitive to the oscillation parameters \({\sin }^{2}{\theta }_{12}\) and \(\Delta {m}_{21}^{2}\) beyond existing experimental constraints; their Gaussian priors are set to be \({\sin }^{2}{\theta }_{12}=0.307\pm 0.013\) and \(\Delta {m}_{21}^{2}=(7.53\pm 0.18)\times 1{0}^{-5}\,{{\rm{eV}}}^{2}\). These values, along with a Gaussian prior on \({\sin }^{2}{\theta }_{13}\), when it is used, come from the 2020 version of the Particle Data Group (PDG) summary tables40, which were current at the time of the original analyses. Updates to these constraints in more recent versions of the PDG do not change any conclusions.

As well as the standard prior flat in δCP, we also studied the effect of a prior flat in \(\sin {\delta }_{{\rm{CP}}}\) and saw no significant changes in conclusions.

Moreover, the experiments define priors for all of the systematic parameters in their models. These definitions are detailed in the individual experiment analyses underlying this work.

Highest posterior probability values and 1σ credible intervals

Extended Data Table 3 summarizes the highest posterior probability values and credible intervals measured jointly by NOvA and T2K.

Additional oscillation parameter plots

The main text shows the 1D posterior distributions and credible intervals for the Jarlskog invariant, δCP and \({\sin }^{2}{\theta }_{23}\), as well as 2D distributions and credible regions for the latter two. In this section, we present the 1D distributions and credible intervals for δCP, \({\sin }^{2}{\theta }_{23}\), \({\sin }^{2}2{\theta }_{13}\) and \(| \Delta {m}_{32}^{2}| \), and 2D distributions and credible regions for all pairwise combinations of these parameters. These are shown in Extended Data Figs. 6–8, for the cases of marginalized over both mass orderings, conditional on the normal ordering and conditional on the inverted ordering, respectively. The distributions and intervals are shown in a triangle plot, in which a lower triangular matrix of plots shows the 1D distributions along the diagonal and the 2D distributions in each of the off-diagonal positions.

Reactor \({\boldsymbol{\Delta }}{{\boldsymbol{m}}}_{{\bf{32}}}^{{\bf{2}}}\)

The energy-dependent \({\bar{\nu }}_{{\rm{e}}}\to {\bar{\nu }}_{{\rm{e}}}\) oscillation probability measured by reactor experiments is sensitive to \(| \Delta {m}_{32}^{2}| \), and reactor measurements of this parameter are expected to agree with long-baseline measurements only under the correct mass ordering assumption. Under the incorrect ordering assumption, these two techniques are expected to measure incorrect values that differ from one another by about 2–3% (ref. 62). Thus, comparing \(| \Delta {m}_{32}^{2}| \) measurements from accelerator and reactor experiments under both mass ordering hypotheses can inform mass ordering discrimination. The Daya Bay experiment47 provides the tightest constraints on θ13 and also reports a 2D \({\theta }_{13}-\Delta {m}_{32}^{2}\) likelihood that we can directly incorporate into our joint fit instead of the θ13-only prior discussed elsewhere in this study.

The mass ordering Bayes factor obtained when using this 2D reactor constraint is 1.4 in favour of the normal ordering, in contrast to 1.3 in favour of the inverted ordering when using the θ13-only reactor constraint. This slight pull towards a preference for the normal ordering is expected, given the relative agreement of the Daya Bay and NOvA+T2K \(| \Delta {m}_{32}^{2}| \) measurements shown in Fig. 2 (inverted ordering) and Extended Data Fig. 9a (normal ordering). However, there remains no statistically significant mass ordering preference in this combination.

Additional global comparisons

In Extended Data Fig. 9, results of the analysis using the default priors are compared with other experimental measurements. The statement on \(\Delta {m}_{32}^{2}\) precision is still valid for the normal ordering assumption. As in the case of the \({\sin }^{2}2{\theta }_{13}\) result (Extended Data Fig. 9b,c), the long-baseline measurements (in this comparison, without applying the prior from reactor measurements) are consistent with reactor experiments, with larger consistency in the normal ordering than the inverted ordering. We do not strongly prefer either octant of \({\sin }^{2}{\theta }_{23}\) (Extended Data Fig. 9d,e), which is consistent with other modern experiments. The joint analysis result for δCP (Extended Data Fig. 9f,g) is consistent with all experiments and their combinations, although the uncertainty remains large.