The NOvA Collaboration

# Expanding neutrino oscillation parameter measurements in NOvA using a
Bayesian approach

M. A. Acero  Universidad del Atlantico, Carrera 30 No. 8-49, Puerto Colombia,
Atlantico, Colombia  B. Acharya  University of Mississippi, University,
Mississippi 38677, USA  P. Adamson  Fermi National Accelerator Laboratory,
Batavia, Illinois 60510, USA  N. Anfimov  Joint Institute for Nuclear
Research, Dubna, Moscow region 141980, Russia  A. Antoshkin  Joint Institute
for Nuclear Research, Dubna, Moscow region 141980, Russia  E. Arrieta-Diaz
Universidad del Magdalena, Carrera 32 No 22-08 Santa Marta, Colombia  L.
Asquith  Department of Physics and Astronomy, University of Sussex, Falmer,
Brighton BN1 9QH, United Kingdom  A. Aurisano  Department of Physics,
University of Cincinnati, Cincinnati, Ohio 45221, USA  A. Back  Indiana
University, Bloomington, Indiana 47405, USA  N. Balashov  Joint Institute for
Nuclear Research, Dubna, Moscow region 141980, Russia  P. Baldi  Department of
Physics and Astronomy, University of California at Irvine, Irvine, California
92697, USA  B. A. Bambah  School of Physics, University of Hyderabad,
Hyderabad, 500 046, India  A. Bat  Bandirma Onyedi Eylül University, Faculty
of Engineering and Natural Sciences, Engineering Sciences Department, 10200,
Bandirma, Balıkesir, Turkey  Department of Physics, Erciyes University,
Kayseri 38030, Turkey  K. Bays  School of Physics and Astronomy, University of
Minnesota Twin Cities, Minneapolis, Minnesota 55455, USA  Illinois Institute
of Technology, Chicago IL 60616, USA  R. Bernstein  Fermi National Accelerator
Laboratory, Batavia, Illinois 60510, USA  T. J. C. Bezerra  Department of
Physics and Astronomy, University of Sussex, Falmer, Brighton BN1 9QH, United
Kingdom  V. Bhatnagar  Department of Physics, Panjab University, Chandigarh,
160 014, India  D. Bhattarai  University of Mississippi, University,
Mississippi 38677, USA  B. Bhuyan  Department of Physics, IIT Guwahati,
Guwahati, 781 039, India  J. Bian  Department of Physics and Astronomy,
University of California at Irvine, Irvine, California 92697, USA  School of
Physics and Astronomy, University of Minnesota Twin Cities, Minneapolis,
Minnesota 55455, USA  A. C. Booth  Particle Physics Research Centre,
Department of Physics and Astronomy, Queen Mary University of London, London
E1 4NS, United Kingdom  Department of Physics and Astronomy, University of
Sussex, Falmer, Brighton BN1 9QH, United Kingdom  R. Bowles  Indiana
University, Bloomington, Indiana 47405, USA  B. Brahma  Department of Physics,
IIT Hyderabad, Hyderabad, 502 205, India  C. Bromberg  Department of Physics
and Astronomy, Michigan State University, East Lansing, Michigan 48824, USA
N. Buchanan  Department of Physics, Colorado State University, Fort Collins,
CO 80523-1875, USA  A. Butkevich  Institute for Nuclear Research of Russia,
Academy of Sciences 7a, 60th October Anniversary prospect, Moscow 117312,
Russia  S. Calvez  Department of Physics, Colorado State University, Fort
Collins, CO 80523-1875, USA  T. J. Carroll  Department of Physics, University
of Texas at Austin, Austin, Texas 78712, USA  Department of Physics,
University of Wisconsin-Madison, Madison, Wisconsin 53706, USA  E. Catano-Mur
Department of Physics, William & Mary, Williamsburg, Virginia 23187, USA  J.
P. Cesar  Department of Physics, University of Texas at Austin, Austin, Texas
78712, USA  A. Chatla  School of Physics, University of Hyderabad, Hyderabad,
500 046, India  R. Chirco  Illinois Institute of Technology, Chicago IL 60616,
USA  S. Chaudhary  Department of Physics, IIT Guwahati, Guwahati, 781 039,
India  B. C. Choudhary  Department of Physics and Astrophysics, University of
Delhi, Delhi 110007, India  A. Christensen  Department of Physics, Colorado
State University, Fort Collins, CO 80523-1875, USA  T. E. Coan  Department of
Physics, Southern Methodist University, Dallas, Texas 75275, USA  A.
Cooleybeck  Department of Physics, University of Wisconsin-Madison, Madison,
Wisconsin 53706, USA  L. Cremonesi  Particle Physics Research Centre,
Department of Physics and Astronomy, Queen Mary University of London, London
E1 4NS, United Kingdom  G. S. Davies  University of Mississippi, University,
Mississippi 38677, USA  P. F. Derwent  Fermi National Accelerator Laboratory,
Batavia, Illinois 60510, USA  P. Ding  Fermi National Accelerator Laboratory,
Batavia, Illinois 60510, USA  Z. Djurcic  Argonne National Laboratory,
Argonne, Illinois 60439, USA  M. Dolce  Department of Physics and Astronomy,
Tufts University, Medford, Massachusetts 02155, USA  D. Doyle  Department of
Physics, Colorado State University, Fort Collins, CO 80523-1875, USA  D.
Dueñas Tonguino  Department of Physics, University of Cincinnati, Cincinnati,
Ohio 45221, USA  E. C. Dukes  Department of Physics, University of Virginia,
Charlottesville, Virginia 22904, USA  A. Dye  University of Mississippi,
University, Mississippi 38677, USA  R. Ehrlich  Department of Physics,
University of Virginia, Charlottesville, Virginia 22904, USA  M. Elkins
Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011,
USA  E. Ewart  Indiana University, Bloomington, Indiana 47405, USA  P. Filip
Institute of Physics, The Czech Academy of Sciences, 182 21 Prague, Czech
Republic  J. Franc  Czech Technical University in Prague, Brehova 7, 115 19
Prague 1, Czech Republic  M. J. Frank  Department of Physics, University of
South Alabama, Mobile, Alabama 36688, USA  H. R. Gallagher  Department of
Physics and Astronomy, Tufts University, Medford, Massachusetts 02155, USA  F.
Gao  Department of Physics, University of Pittsburgh, Pittsburgh, Pennsylvania
15260, USA  A. Giri  Department of Physics, IIT Hyderabad, Hyderabad, 502 205,
India  R. A. Gomes  Instituto de Física, Universidade Federal de Goiás,
Goiânia, Goiás, 74690-900, Brazil  M. C. Goodman  Argonne National Laboratory,
Argonne, Illinois 60439, USA  M. Groh  Department of Physics, Colorado State
University, Fort Collins, CO 80523-1875, USA  Indiana University, Bloomington,
Indiana 47405, USA  R. Group  Department of Physics, University of Virginia,
Charlottesville, Virginia 22904, USA  A. Habig  Department of Physics and
Astronomy, University of Minnesota Duluth, Duluth, Minnesota 55812, USA  F.
Hakl  Institute of Computer Science, The Czech Academy of Sciences, 182 07
Prague, Czech Republic  J. Hartnell  Department of Physics and Astronomy,
University of Sussex, Falmer, Brighton BN1 9QH, United Kingdom  R. Hatcher
Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA  M. He
Department of Physics, University of Houston, Houston, Texas 77204, USA  K.
Heller  School of Physics and Astronomy, University of Minnesota Twin Cities,
Minneapolis, Minnesota 55455, USA  V Hewes  Department of Physics, University
of Cincinnati, Cincinnati, Ohio 45221, USA  A. Himmel  Fermi National
Accelerator Laboratory, Batavia, Illinois 60510, USA  B. Jargowsky  Department
of Physics and Astronomy, University of California at Irvine, Irvine,
California 92697, USA  J. Jarosz  Department of Physics, Colorado State
University, Fort Collins, CO 80523-1875, USA  F. Jediny  Czech Technical
University in Prague, Brehova 7, 115 19 Prague 1, Czech Republic  C. Johnson
Department of Physics, Colorado State University, Fort Collins, CO 80523-1875,
USA  M. Judah  Department of Physics, Colorado State University, Fort Collins,
CO 80523-1875, USA  Department of Physics, University of Pittsburgh,
Pittsburgh, Pennsylvania 15260, USA  I. Kakorin  Joint Institute for Nuclear
Research, Dubna, Moscow region 141980, Russia  D. M. Kaplan  Illinois
Institute of Technology, Chicago IL 60616, USA  A. Kalitkina  Joint Institute
for Nuclear Research, Dubna, Moscow region 141980, Russia  J. Kleykamp
University of Mississippi, University, Mississippi 38677, USA  O. Klimov
Joint Institute for Nuclear Research, Dubna, Moscow region 141980, Russia  L.
W. Koerner  Department of Physics, University of Houston, Houston, Texas
77204, USA  L. Kolupaeva  Joint Institute for Nuclear Research, Dubna, Moscow
region 141980, Russia  R. Kralik  Department of Physics and Astronomy,
University of Sussex, Falmer, Brighton BN1 9QH, United Kingdom  A. Kumar
Department of Physics, Panjab University, Chandigarh, 160 014, India  C. D.
Kuruppu  Department of Physics and Astronomy, University of South Carolina,
Columbia, South Carolina 29208, USA  V. Kus  Czech Technical University in
Prague, Brehova 7, 115 19 Prague 1, Czech Republic  T. Lackey  Fermi National
Accelerator Laboratory, Batavia, Illinois 60510, USA  Indiana University,
Bloomington, Indiana 47405, USA  K. Lang  Department of Physics, University of
Texas at Austin, Austin, Texas 78712, USA  P. Lasorak  Department of Physics
and Astronomy, University of Sussex, Falmer, Brighton BN1 9QH, United Kingdom
J. Lesmeister  Department of Physics, University of Houston, Houston, Texas
77204, USA  A. Lister  Department of Physics, University of Wisconsin-Madison,
Madison, Wisconsin 53706, USA  J. Liu  Department of Physics and Astronomy,
University of California at Irvine, Irvine, California 92697, USA  M.
Lokajicek  Institute of Physics, The Czech Academy of Sciences, 182 21 Prague,
Czech Republic  M. MacMahon  Physics and Astronomy Department, University
College London, Gower Street, London WC1E 6BT, United Kingdom  J. A. Lock
Department of Physics and Astronomy, University of Sussex, Falmer, Brighton
BN1 9QH, United Kingdom  S. Magill  Argonne National Laboratory, Argonne,
Illinois 60439, USA  M. Manrique Plata  Indiana University, Bloomington,
Indiana 47405, USA  W. A. Mann  Department of Physics and Astronomy, Tufts
University, Medford, Massachusetts 02155, USA  M. T. Manoharan  Department of
Physics, Cochin University of Science and Technology, Kochi 682 022, India  M.
L. Marshak  School of Physics and Astronomy, University of Minnesota Twin
Cities, Minneapolis, Minnesota 55455, USA  M. Martinez-Casales  Department of
Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA  V.
Matveev  Institute for Nuclear Research of Russia, Academy of Sciences 7a,
60th October Anniversary prospect, Moscow 117312, Russia  B. Mehta  Department
of Physics, Panjab University, Chandigarh, 160 014, India  M. D. Messier
Indiana University, Bloomington, Indiana 47405, USA  H. Meyer  Department of
Mathematics, Statistics, and Physics, Wichita State University, Wichita,
Kansas 67260, USA  T. Miao  Fermi National Accelerator Laboratory, Batavia,
Illinois 60510, USA  V. Mikola  Physics and Astronomy Department, University
College London, Gower Street, London WC1E 6BT, United Kingdom  W. H. Miller
School of Physics and Astronomy, University of Minnesota Twin Cities,
Minneapolis, Minnesota 55455, USA  S. Mishra  Department of Physics, Institute
of Science, Banaras Hindu University, Varanasi, 221 005, India  S. R. Mishra
Department of Physics and Astronomy, University of South Carolina, Columbia,
South Carolina 29208, USA  R. Mohanta  School of Physics, University of
Hyderabad, Hyderabad, 500 046, India  A. Moren  Department of Physics and
Astronomy, University of Minnesota Duluth, Duluth, Minnesota 55812, USA  A.
Morozova  Joint Institute for Nuclear Research, Dubna, Moscow region 141980,
Russia  W. Mu  Fermi National Accelerator Laboratory, Batavia, Illinois 60510,
USA  L. Mualem  California Institute of Technology, Pasadena, California
91125, USA  M. Muether  Department of Mathematics, Statistics, and Physics,
Wichita State University, Wichita, Kansas 67260, USA  K. Mulder  Physics and
Astronomy Department, University College London, Gower Street, London WC1E
6BT, United Kingdom  D. Myers  Department of Physics, University of Texas at
Austin, Austin, Texas 78712, USA  D. Naples  Department of Physics, University
of Pittsburgh, Pittsburgh, Pennsylvania 15260, USA  A. Nath  Department of
Physics, IIT Guwahati, Guwahati, 781 039, India  S. Nelleri  Department of
Physics, Cochin University of Science and Technology, Kochi 682 022, India  J.
K. Nelson  Department of Physics, William & Mary, Williamsburg, Virginia
23187, USA  R. Nichol  Physics and Astronomy Department, University College
London, Gower Street, London WC1E 6BT, United Kingdom  E. Niner  Fermi
National Accelerator Laboratory, Batavia, Illinois 60510, USA  A. Norman
Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA  A.
Norrick  Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA
T. Nosek  Charles University, Faculty of Mathematics and Physics, Institute of
Particle and Nuclear Physics, Prague, Czech Republic  H. Oh  Department of
Physics, University of Cincinnati, Cincinnati, Ohio 45221, USA  A. Olshevskiy
Joint Institute for Nuclear Research, Dubna, Moscow region 141980, Russia  T.
Olson  Department of Physics, University of Houston, Houston, Texas 77204, USA
A. Pal  National Institute of Science Education and Research, Khurda, 752050,
Odisha, India  J. Paley  Fermi National Accelerator Laboratory, Batavia,
Illinois 60510, USA  L. Panda  National Institute of Science Education and
Research, Khurda, 752050, Odisha, India  R. B. Patterson  California Institute
of Technology, Pasadena, California 91125, USA  G. Pawloski  School of Physics
and Astronomy, University of Minnesota Twin Cities, Minneapolis, Minnesota
55455, USA  O. Petrova  Joint Institute for Nuclear Research, Dubna, Moscow
region 141980, Russia  R. Petti  Department of Physics and Astronomy,
University of South Carolina, Columbia, South Carolina 29208, USA  R. K.
Plunkett  Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA
A. Rafique  Argonne National Laboratory, Argonne, Illinois 60439, USA  L. R.
Prais  University of Mississippi, University, Mississippi 38677, USA  V. Raj
California Institute of Technology, Pasadena, California 91125, USA  M.
Rajaoalisoa  Department of Physics, University of Cincinnati, Cincinnati, Ohio
45221, USA  B. Ramson  Fermi National Accelerator Laboratory, Batavia,
Illinois 60510, USA  B. Rebel  Fermi National Accelerator Laboratory, Batavia,
Illinois 60510, USA  Department of Physics, University of Wisconsin-Madison,
Madison, Wisconsin 53706, USA  P. Roy  Department of Mathematics, Statistics,
and Physics, Wichita State University, Wichita, Kansas 67260, USA  O. Samoylov
Joint Institute for Nuclear Research, Dubna, Moscow region 141980, Russia  M.
C. Sanchez  Florida State University, Tallahassee, Florida 32306, USA
Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011,
USA  S. Sánchez Falero  Department of Physics and Astronomy, Iowa State
University, Ames, Iowa 50011, USA  P. Shanahan  Fermi National Accelerator
Laboratory, Batavia, Illinois 60510, USA  P. Sharma  Department of Physics,
Panjab University, Chandigarh, 160 014, India  S. Shukla  Department of
Physics, Institute of Science, Banaras Hindu University, Varanasi, 221 005,
India  A. Sheshukov  Joint Institute for Nuclear Research, Dubna, Moscow
region 141980, Russia  D. K. Singha  School of Physics, University of
Hyderabad, Hyderabad, 500 046, India  W. Shorrock  Department of Physics and
Astronomy, University of Sussex, Falmer, Brighton BN1 9QH, United Kingdom  I.
Singh  Department of Physics and Astrophysics, University of Delhi, Delhi
110007, India  P. Singh  Particle Physics Research Centre, Department of
Physics and Astronomy, Queen Mary University of London, London E1 4NS, United
Kingdom  Department of Physics and Astrophysics, University of Delhi, Delhi
110007, India  V. Singh  Department of Physics, Institute of Science, Banaras
Hindu University, Varanasi, 221 005, India  E. Smith  Indiana University,
Bloomington, Indiana 47405, USA  J. Smolik  Czech Technical University in
Prague, Brehova 7, 115 19 Prague 1, Czech Republic  P. Snopok  Illinois
Institute of Technology, Chicago IL 60616, USA  N. Solomey  Department of
Mathematics, Statistics, and Physics, Wichita State University, Wichita,
Kansas 67260, USA  A. Sousa  Department of Physics, University of Cincinnati,
Cincinnati, Ohio 45221, USA  K. Soustruznik  Charles University, Faculty of
Mathematics and Physics, Institute of Particle and Nuclear Physics, Prague,
Czech Republic  M. Strait  School of Physics and Astronomy, University of
Minnesota Twin Cities, Minneapolis, Minnesota 55455, USA  L. Suter  Fermi
National Accelerator Laboratory, Batavia, Illinois 60510, USA  A. Sutton
Florida State University, Tallahassee, Florida 32306, USA  Department of
Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA
Department of Physics, University of Virginia, Charlottesville, Virginia
22904, USA  K. Sutton  California Institute of Technology, Pasadena,
California 91125, USA  S. Swain  National Institute of Science Education and
Research, Khurda, 752050, Odisha, India  C. Sweeney  Physics and Astronomy
Department, University College London, Gower Street, London WC1E 6BT, United
Kingdom  A. Sztuc  Physics and Astronomy Department, University College
London, Gower Street, London WC1E 6BT, United Kingdom  B. Tapia Oregui
Department of Physics, University of Texas at Austin, Austin, Texas 78712, USA
P. Tas  Charles University, Faculty of Mathematics and Physics, Institute of
Particle and Nuclear Physics, Prague, Czech Republic  T. Thakore  Department
of Physics, University of Cincinnati, Cincinnati, Ohio 45221, USA  J. Thomas
Physics and Astronomy Department, University College London, Gower Street,
London WC1E 6BT, United Kingdom  Department of Physics, University of
Wisconsin-Madison, Madison, Wisconsin 53706, USA  E. Tiras  Department of
Physics, Erciyes University, Kayseri 38030, Turkey  Department of Physics and
Astronomy, Iowa State University, Ames, Iowa 50011, USA  J. Trokan-Tenorio
Department of Physics, William & Mary, Williamsburg, Virginia 23187, USA  Y.
Torun  Illinois Institute of Technology, Chicago IL 60616, USA  J. Urheim
Indiana University, Bloomington, Indiana 47405, USA  P. Vahle  Department of
Physics, William & Mary, Williamsburg, Virginia 23187, USA  Z. Vallari
California Institute of Technology, Pasadena, California 91125, USA  K. J.
Vockerodt  Particle Physics Research Centre, Department of Physics and
Astronomy, Queen Mary University of London, London E1 4NS, United Kingdom  T.
Vrba  Czech Technical University in Prague, Brehova 7, 115 19 Prague 1, Czech
Republic  M. Wallbank  Department of Physics, University of Cincinnati,
Cincinnati, Ohio 45221, USA  T. K. Warburton  Department of Physics and
Astronomy, Iowa State University, Ames, Iowa 50011, USA  M. Wetstein
Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011,
USA  D. Whittington  Department of Physics, Syracuse University, Syracuse NY
13210, USA  D. A. Wickremasinghe  Fermi National Accelerator Laboratory,
Batavia, Illinois 60510, USA  T. Wieber  School of Physics and Astronomy,
University of Minnesota Twin Cities, Minneapolis, Minnesota 55455, USA  J.
Wolcott  Department of Physics and Astronomy, Tufts University, Medford,
Massachusetts 02155, USA  M. Wrobel  Department of Physics, Colorado State
University, Fort Collins, CO 80523-1875, USA  S. Wu  School of Physics and
Astronomy, University of Minnesota Twin Cities, Minneapolis, Minnesota 55455,
USA  W. Wu  Department of Physics and Astronomy, University of California at
Irvine, Irvine, California 92697, USA  Y. Xiao  Department of Physics and
Astronomy, University of California at Irvine, Irvine, California 92697, USA
B. Yaeggy  Department of Physics, University of Cincinnati, Cincinnati, Ohio
45221, USA  A. Yankelevich  Department of Physics and Astronomy, University of
California at Irvine, Irvine, California 92697, USA  K. Yonehara  Fermi
National Accelerator Laboratory, Batavia, Illinois 60510, USA  Y. Yu  Illinois
Institute of Technology, Chicago IL 60616, USA  S. Zadorozhnyy  Institute for
Nuclear Research of Russia, Academy of Sciences 7a, 60th October Anniversary
prospect, Moscow 117312, Russia  J. Zalesak  Institute of Physics, The Czech
Academy of Sciences, 182 21 Prague, Czech Republic  R. Zwaska  Fermi National
Accelerator Laboratory, Batavia, Illinois 60510, USA

###### Abstract

NOvA is a long-baseline neutrino oscillation experiment that measures
oscillations in charged-current $\nu_{\mu}\rightarrow\nu_{\mu}$
(disappearance) and $\nu_{\mu}\rightarrow\nu_{e}$ (appearance) channels, and
their antineutrino counterparts, using neutrinos of energies around 2 GeV over
a distance of 810 km. In this work we reanalyze the dataset first examined in
our previous paper [Phys. Rev. D 106, 032004 (2022)] using an alternative
statistical approach based on Bayesian Markov Chain Monte Carlo. We measure
oscillation parameters consistent with the previous results. We also extend
our inferences to include the first NOvA measurements of the reactor mixing
angle $\theta_{13}$ and the Jarlskog invariant. We use these results to
quantify the strength of our inferences about CP violation, as well as to
examine the effects of constraints from short-baseline measurements of
$\theta_{13}$ using antineutrinos from nuclear reactors when making NOvA
measurements of $\theta_{23}$. Our long-baseline measurement of $\theta_{13}$
is also shown to be consistent with the reactor measurements, supporting the
general applicability and robustness of the PMNS framework for neutrino
oscillations.

††preprint: APS/123-QED

##  I Introduction

In the three-flavor neutrino oscillation paradigm, transitions among the three
flavor eigenstates $\nu_{e}$, $\nu_{\mu}$, $\nu_{\tau}$ are governed by the
matrix elements between these states and the mass eigenstates $\nu_{1}$,
$\nu_{2}$, $\nu_{3}$.111The same is true for their antineutrino counterparts
$\bar{\nu}_{e}$, $\bar{\nu}_{\mu}$, $\bar{\nu}_{\tau}$ and $\bar{\nu}_{1}$,
$\bar{\nu}_{2}$, $\bar{\nu}_{3}$. Throughout this paper the symbol $\nu$ will
refer to both neutrinos and antineutrinos unless otherwise specified. The
elements of this Pontecorvo–Maki–Nakagawa–Sakata (PMNS) matrix have been
constrained by numerous experiments using neutrinos with a variety of energy
spectra over various baselines [1, 2, 3, 4, 5, 6, 7, 8, 9]. Decomposing the
PMNS matrix yields a set of rotation-like “mixing angles” $\theta_{12}$,
$\theta_{13}$, $\theta_{23}$, and a phase $\delta_{\textrm{CP}}$. Contemporary
neutrino oscillation experiments seek to make precision measurements of these
parameters, as well as the differences between the squared mass eigenvalues
($\Delta m^{2}_{ij}\equiv m_{i}^{2}-m_{j}^{2}$). These results have
fundamental implications for models of neutrino mass and lepton flavor [10,
11, 12, 13, 14], as well as models of baryogenesis via charge-parity (CP)
symmetry violation [15, 16, 17, 18, 19]. Where multiple experiments can access
the same parameters using different neutrino flavors or energies, the overall
validity of the three-neutrino framework can also be tested.

Long-baseline (LBL) accelerator neutrino oscillation experiments measure
oscillations in $\nu_{\mu}\rightarrow\nu_{\mu}$ (disappearance) and
$\nu_{\mu}\rightarrow\nu_{e}$ (appearance) channels. These channels constrain
the mixing angles $\theta_{13}$, $\theta_{23}$, the mass-squared splitting
$\Delta m^{2}_{32}$, and the CP-violating phase $\delta_{\textrm{CP}}$.
Current measurements of $\theta_{23}$ are consistent with maximal mixing
($\theta_{23}$ = $\pi/4$), which would suggest a $\mu\text{-}\tau$ symmetry in
their mixing into the $\nu_{3}$ mass eigenstate; a nonmaximal value such as
$\theta_{23}$ $<\pi/4$ (lower octant, LO) or $\theta_{23}$ $>\pi/4$ (upper
octant, UO) would indicate a preferential coupling of $\nu_{\tau}$ or
$\nu_{\mu}$, respectively, with $\nu_{3}$. Current experimental uncertainties
on $\theta_{23}$ are the largest among the mixing angles [20]. LBL oscillation
measurements, where neutrinos traverse significant quantities of matter, are
also impacted by the coherent forward scattering of $\nu_{e}$s on electrons in
the Earth [21]. This modifies the oscillation probabilities
$\mathcal{P}(\nu_{\mu}\rightarrow\nu_{e})$ and
$\mathcal{P}(\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e})$ with opposite signs.
The direction of the resulting change in rate depends on whether $\nu_{3}$ is
the heaviest neutrino state (Normal Ordering, NO) or the lightest (Inverted
Ordering, IO). Current observations from LBL (and other) experiments prefer
the NO, though the strength of the NO hypothesis depends on which measurements
are included [22]. The value of $\delta_{\textrm{CP}}$ is currently
essentially unknown; LBL experiments provide the only constraints, but they
are weak.

NOvA is a long-baseline neutrino oscillation experiment that observes the
$\nu_{\mu}$ disappearance and $\nu_{e}$ appearance channels using neutrinos of
energies around 2 GeV over a distance of 810 km. Previously [23, 24, 25, 26,
27, 28, 29], we have presented NOvA constraints on $\Delta m^{2}_{32}$,
$\sin^{2}\theta_{23}$, and $\delta_{\textrm{CP}}$ using a classical
frequentist approach. However, the Feldman–Cousins technique that is required
in order to obtain correct classical confidence regions for these variables
[30, 31] poses challenges when confronted with highly degenerate sets of
parameters. It also does not allow for post hoc transformations of the
variables considered in the analysis. In this work we present a new analysis
of the dataset from Ref. [29], based on Bayesian Markov Chain Monte Carlo,
which enables us to extend our inferences to include $\theta_{13}$ and the
Jarlskog invariant $J$, for which $J\neq 0$ unambiguously indicates CP
violation. We use these results to examine the implications of assuming short-
baseline, nuclear-reactor antineutrino constraints on $\theta_{13}$ when
making measurements of other oscillation parameters in NOvA.

We also investigate the consistency of the PMNS framework by comparing the
constraint from reactor experiments with our long-baseline measurement of
$\theta_{13}$.

##  II 3-flavor neutrino oscillations in NOvA

In this paper, we reanalyze data collected from an exposure of $13.6\times
10^{20}$ 14 kton-equivalent protons on target (POT) in the neutrino-enriched
beam mode and $12.5\times 10^{20}$ POT in the analogous antineutrino mode. The
dataset, simulations, reconstruction, and estimation of systematic
uncertainties remain unchanged for this analysis. A brief overview of these
components follows, with detailed descriptions available in Ref. [29].
Extensive discussion of the new analysis method, its implementations, and the
resulting inferences are presented in Sec. III.

###  II.1 The NOvA experiment

NOvA observes oscillations using neutrinos from the NuMI beamline [32] at
Fermilab using two functionally identical tracking calorimeter detectors that
differ primarily in size. The 0.3 kton near detector (ND), the smaller of the
detectors, is located 1 km from the neutrino production target, 100 m
underground. The far detector (FD), by contrast, is 14 kton and is located on
the surface at Ash River, Minnesota, 810 km from the target. Both detectors
are built from PVC cells of $3.9\times 6.6\,\mathrm{cm^{2}}$ cross-sectional
area and 3.9 m (ND) or 15.5 m (FD) length, which are arranged in alternating
horizontal and vertical planes and filled with a mineral-oil-based liquid
scintillator. A stack of alternating active planes and steel plates is placed
behind the ND to range out muons while a small rock overburden is placed above
the FD to aid in rejecting the cosmic-ray background. The detectors are placed
14.6 mrad from the central axis of the beam to receive a narrow-band neutrino
flux predominantly between 1 and 3 GeV.

###  II.2 Simulation and selection

We use Geant4 (v10.4) [33, 34] to simulate the production of hadrons from
interactions of the primary proton beam with the target as well as their
transport through the beam optics. These simulations are reweighted using the
Package to Predict the FluX (PPFX) [35] to include constraints from external
hadron production data [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,
49, 50, 51, 52, 53, 54]. Simulated interactions of neutrinos that arise from
decays of those hadrons are generated using GENIE 3.0.6 [55, 56] and modified
by corrections we derive from NOvA ND and external data. NOvA ND data is used
to produce a NOvA-specific tune of the IFIC València model [57, 58] that
describes charged-current neutrino scattering from correlated pairs of
nucleons. The NOvA tune also modifies final-state interactions (simulated with
the GENIE hN full intranuclear cascade model) using pion–nucleon scattering
data [59, 60, 61, 62, 63, 64, 65]. The outgoing final-state particles are
propagated through the detector using Geant4 and a custom NOvA readout
simulation [66].

We use groups of spatially and temporally proximate cells with activity above
threshold to apply basic data quality and containment selection cuts to events
in both data and simulation. Within these groups we assign vertices and
reconstruct likely particle trajectories. Ultimately we divide the events into
$\nu_{\mu}$ charged-current (CC), $\nu_{e}$ CC, neutral-current (NC), or
cosmogenic background categories using NOvA’s convolutional neural network
(CNN)-based classifier [67]. Boosted decision trees (BDTs) are used to further
reject cosmic backgrounds in the FD samples. Both sets of tags are utilized
together to create $\nu_{\mu}$ CC and $\nu_{e}$ CC candidate samples. Fully
contained $\nu_{e}$ candidates at the FD are further divided into high and low
purity samples based on the CNN score in order to enhance the signal-to-
background rejection capability of the fit. To improve the statistical
strength of the fit, we recover an additional sample of “peripheral” events
that fail the containment or cosmic rejection BDT but pass stricter particle
ID requirements. We estimate neutrino energy for $\nu_{\mu}$ CC events using
the muon track length and the total deposited calorimetric energy of the
hadronic system. The energy for $\nu_{e}$ CC events is estimated using a
function of calorimetric energy that takes as input the energy of the event’s
reconstructed trajectories divided into electromagnetic and hadronic
components, as identified by a separate CNN-based classifier [68].

###  II.3 Near-to-far extrapolation and systematics

Predictions for the neutrino event rates at the FD are constrained using the
high-statistics neutrino interactions measured in the ND. These measurements
imply corrections to the ND prediction, which we propagate to the FD by
adjusting for the differing efficiency and flux between the detectors using
simulations. We call this process “extrapolation.” We apply oscillations to
this data-driven prediction when comparing to FD data during inference (Sec.
III.2 below).

The signal spectra for both $\nu_{\mu}$ disappearance and $\nu_{e}$ appearance
channels are predicted using the $\nu_{\mu}$ spectra at the ND. The near-to-
far extrapolation for the $\nu_{\mu}$ disappearance samples at the FD is
performed in quartiles of hadronic energy fraction
($f_{\mathrm{had}}=E_{\mathrm{had}}/E_{\nu}$, with $E_{\nu}$ being the
reconstructed neutrino energy and $E_{\mathrm{had}}$ the reconstructed
hadronic energy).

Figure 1: Reconstructed neutrino energy distribution of selected data events
(black crosses) in FD $\nu_{e}$ CC samples (top) and FD $\nu_{\mu}$ CC samples
(bottom) in neutrino-enriched beam mode (left) and antineutrino-enriched beam
mode (right). The colored bands correspond to the range of 68.3% (darkest),
95.4%, and 99.7% (lightest) of the extrapolated FD spectra produced using the
combinations of the oscillation and systematic parameters sampled by our MCMC
algorithms (illustrating the posterior distributions described further in Sec.
III.2). The FD $\nu_{e}$ samples are divided into bins of low (I) and high
(II) particle ID confidence as well as the peripheral (III) sample discussed
in Sec. II.2. The four $E_{frac}$ $\nu_{\mu}$ subsamples have been combined
together in each of the lower two plots.

Extrapolating in the $f_{\mathrm{had}}$ bins has the effect of grouping events
that share similar hadronic system characteristics so that compatible events
are constrained together in the FD sample, despite the detectors’ somewhat
different acceptances. Performing the oscillation fit in these bins enhances
sensitivity to the oscillation dip as a function of $E_{\nu}$ since we achieve
a finer energy resolution for bins with smaller $f_{\mathrm{had}}$. We further
subdivide the near-to-far extrapolation into three bins of reconstructed
transverse momentum of the outgoing charged lepton $p_{T}$ for both appearance
and disappearance channels. Like the $f_{\mathrm{had}}$ subdivision, this
allows us to better match the constraints from the much smaller ND to the FD,
in this case adjusting for the differing containment of events with leptons
that emerge at large angles relative to the beam direction (i.e., large
$p_{T}$). The predictions in $p_{T}$ bins are summed prior to the oscillation
fit. We constrain the small beam backgrounds in the $\nu_{e}$ appearance
channels with a similar procedure based on ND $\nu_{e}$ candidates after first
decomposing them into NC, $\nu_{\mu}$ CC, and intrinsic $\nu_{e}$ categories
using data-driven constraints. ($\bar{\nu}_{e}$ beam backgrounds are all
constrained together rather than being decomposed this way.) Cosmic
backgrounds for all FD samples are determined from dedicated FD cosmic data
samples. The remaining minor backgrounds are estimated from simulation. More
detail on these procedures can be found in our previous paper [29].

We evaluate the impact of systematic uncertainties on the analysis by
repredicting the sample spectra described above with altered parameters in the
simulation. Uncertainties in the neutrino flux and interaction model are
treated using event reweighting. Uncertainties in the detector calibration and
custom modeling of light in the detectors, on the other hand, must be fully
resimulated. The 67 systematically shifted simulations that arise from these
techniques are extrapolated to the FD using ND data via the same process as
above, which constrains the impact of the uncertainties on the predicted
spectra. Parameters corresponding to these uncertainties are marginalized over
in the oscillation inference in Sec. III.

Figure 1 shows the reconstructed energy spectra of the data observed at the FD
using the $\nu_{e}$ and $\nu_{\mu}$ CC selections described in Sec. II.2,
separated by the neutrino-enriched and antineutrino-enriched beam modes.
Overlaid on these plots are bands of FD predictions produced according to the
extrapolation procedure above. These illustrate the spectra predicted using
the 68.3%, 95.4%, and 99.7% highest probability values of the systematics just
described and of the relevant oscillation parameters, determined using the
MCMC algorithms that will be detailed in the next section.

A Poisson likelihood [20] computed over the bins between the FD data and
predictions such as those shown here forms the likelihood used in Bayes’
theorem below.

##  III Oscillation inferences using Markov Chain Monte Carlo

We derive posterior probability distributions for relevant oscillation
parameters using Bayes’ theorem [69]. Marginalizing away the nuisance
parameters of our model, which include the tens of systematic uncertainties
described in Sec. II.3, is a challenging problem because it requires an
integral over many dimensions. We therefore turn to a Monte Carlo method for
computing the posterior: Markov Chain Monte Carlo (MCMC). In MCMC, we draw a
collection of sequential samples from the posterior with frequency
proportional to the posterior density. Histograms that approximate the
posterior shape (with accuracy governed by the sample count) may be computed
in any variable(s) of interest using these samples. In so doing, any
dimensions not explicitly summed are implicitly marginalized.222For an
accessible introduction to MCMC, the reader is referred to Ref. [70]. An
exhaustive treatment may be found at Ref. [71].

Numerous algorithms for obtaining MCMC samples exist; we have implemented two
for this analysis. The conclusions obtained from them agree with one another.
Though descriptions of both methods are readily found in the literature, there
are certain implementation choices that must be made for each, which we
discuss in Sec. III.1. Sec. III.2 lays out our resulting inferences on the
parameters.

###  III.1 NOvA MCMC implementations

####  III.1.1 “ARIA:” $\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$algorithm

The traditional MCMC algorithm, the
Metropolis–Rosenbluth–Rosenbluth–Teller–Teller
($\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$) method333Historically this method
was known as the “Metropolis” (or “Metropolis–Hastings”) method, referring to
the first author of the seminal papers [72, 73]. We follow more recent
convention and hereafter refer to it by all of the authors’ names. See also
Ref. [74]. is straightforward. We call our implementation ARIA, in honor of
Arianna Rosenbluth, who first implemented the method in machine code [75].
Proceeding from an initial seed in the parameter space, subsequent samples are
selected by proposing a jump to a new set of coordinates, and accepting or
rejecting that proposal according to an acceptance rule [20]. This process is
repeated until a sufficient number of samples have been collected. There is no
explicit stopping criterion. The method also does not specify the distribution
to be used in the proposal algorithm. In our implementation, we use the most
common choice, which is of a multivariate Gaussian. We determine the
characteristic length scale of the Gaussian empirically in order to optimize
sampling efficiency (see App. A). We also “thin” the resulting chains to
reduce autocorrelations among the samples (see App. B). Our ARIA results below
have $5\times 10^{5}$ effective samples.

####  III.1.2 “Stan:” Hamiltonian MCMC

Though the $\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$ method proposes samples
quickly, they are typically highly autocorrelated. Its sampling proposals can
also be inefficient if the posterior is sharply concentrated. Other MCMC
methods have been developed to address these shortcomings, including one
called “Hamiltonian” MCMC inference (HMCMC). We implemented a C++ interface to
the Stan modeling platform [76] to obtain HMCMC samples.

The main difference between HMCMC and $\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$
is how proposals are generated. Rather than proposing randomly, HMCMC views
the posterior surface as a topographical one that can be explored by a
fictitious particle. Samples correspond to trajectories under the influence of
a gravitational potential whose gradient corresponds to that of higher
posterior density. Endowing the particle with an initial momentum that
counterbalances the centripetal force from gravitation produces stable
trajectories that traverse the highest density region of posterior space [77].
HMCMC does this by numerically integrating Hamilton’s equations for the
fictitious particle system with its position $\vec{q}$ (which correspond to
the parameters of interest) and momentum $\vec{p}$ coordinates, and a
Hamiltonian $H=-\log(\text{posterior})$. This approach produces samples that
are nearly uncorrelated without thinning at the expense of additional
computing cycles to compute the gradient of the posterior. (Compare the
$7\times 10^{7}$ Stan samples we obtained, which require no thinning, to the
$5\times 10^{5}$ effective samples mentioned above for ARIA that remain after
thinning.) We find Stan’s default choices of the sampling distribution for the
pseudoparticle kinetic energies and the integration stopping condition to be
sufficient for our needs (see App. C).

Its topographical nature means HMCMC is ill-suited to parameters that assume
only one of a discrete set of values, which would manifest as discontinuities
in the trajectories considered. This presents a difficulty in neutrino
oscillation parameter inference, where the absolute value of $\Delta
m^{2}_{32}$ is known with relatively good precision, but its sign remains an
important unknown. While it is possible to allow HMCMC to explore the entire
range of $\Delta m^{2}_{32}$, we find in practice that this results in poor
exploration, as few trajectories manage to “jump” across the wide disfavored
region at low $|\Delta m^{2}_{32}|$. Instead, at the end of each trajectory
determination, we introduce a separate
$\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$-like step which considers the
possibility of changing the sign of $\Delta m^{2}_{32}$ according to the prior
probability chosen for it (50%, i.e., uniform prior; see Sec. III.1.3 below).
If the acceptance ratio between both mass orderings satisfies the
$\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$ criteria, the proposed sign is
retained; if not, it is reverted to its previous value.

####  III.1.3 Choices of prior

In Bayesian inference, the posterior probabilities are influenced by the
choice of prior probability densities, which encode assumptions made about the
parameters before the data is examined. If the data used for a measurement is
sufficient, its constraint on the posterior will overwhelm the prior, and the
prior choice is unimportant. However, when data is sparse, the choice of prior
may affect the result. The priors we choose differ according to the parameters
being considered.

##### Parameters of interest: $|\Delta m^{2}_{32}|$, $\mathrm{sgn}(\Delta
m^{2}_{32})$, $\sin^{2}\theta_{23}$

We prefer to use “uninformed” priors, which do not favor any particular value,
for the physics parameters we intend to measure directly. In practice, this
usually amounts to a prior uniform in the variable in question. However,
uniformity is not preserved under a change of variable: for instance, a prior
uniform in $\theta_{23}$ is not uniform in the measured variable
$\sin^{2}\theta_{23}$. In our results below, we have studied the impact of
priors uniform both in a particular variable and relevant functions of it, and
we report when the prior choice significantly affects the results.

##### Parameter of interest: $\delta_{\textrm{CP}}$

While $\delta_{\textrm{CP}}$ is intended to have a uniform prior as well, it
receives additional special treatment. Its cyclical nature as the phase of a
complex number results in an infinite set of values having identical
consistency with the data. This can cause MCMC samplers never to converge on a
single value for $\delta_{\textrm{CP}}$ and consequently to sample much more
slowly. To combat this problem, we developed a novel special prior over
$\delta_{\textrm{CP}}$:

$\Pi(\delta_{\textrm{CP}})=\begin{cases}\frac{1}{2}\sin^{2}\left(\frac{1}{4}(\delta_{\textrm{CP}}+\pi)\right),&-1\leq\delta_{\textrm{CP}}/\pi\leq
3\\\ 0,&\text{otherwise}.\end{cases}$ (1)

This function is illustrated in Fig. 2.

Figure 2: Special prior distribution from Eq. 1 used for
$\delta_{\textrm{CP}}$.

This prior forces the samplers to remain near a single phase of
$\delta_{\textrm{CP}}$, $0\leq\delta_{\textrm{CP}}/\pi\leq 2$, as the prior
vanishes outside of [-1, 3]. Because it makes the transition to the vanishing
regions in a differentiable manner, this prior is suitable for use with HMCMC.
Moreover, the sum of the prior’s value at every point
$0\leq\delta_{\textrm{CP}}/\pi\leq 2$ with that of all its complex branch
points outside that interval is a value that does not depend on
$\delta_{\textrm{CP}}$:

$\sum_{n=-\infty}^{\infty}\Pi(\delta_{\textrm{CP}}+2\pi n)=\frac{1}{2}.$ (2)

This implies that the prior behaves identically to a uniform prior when used
in conjunction with the oscillation probability. For the rest of this paper,
whenever we refer to a “uniform” prior in $\delta_{\textrm{CP}}$, we mean this
prior.

When in the subsequent sections we study the effect of choosing a prior
uniform in $\delta_{\textrm{CP}}$ vs. $\sin(\delta_{\textrm{CP}})$, we
reweight the MCMC samples obtained with the prior above, taking the Jacobian
factor $\partial(\sin(\delta_{\textrm{CP}}))=\cos(\delta_{\textrm{CP}})$ as
the weight.

##### Parameter of interest: $\sin^{2}2\theta_{13}$

As discussed further in Sec. III.2 below, we consider two separate cases for
$\sin^{2}2\theta_{13}$: one where its prior is treated identically to that of
$\sin^{2}\theta_{23}$ (see III.1.3), and one where measurements from short-
baseline reactor antineutrino oscillations are applied as a constraint. In the
latter case, we impose a Gaussian prior with standard deviation obtained from
the 2019 world average of reactor measurements [78], analogous to the
treatment in our most recent frequentist result [29]:
$\sin^{2}2\theta_{13}=0.085\pm 0.003$.

##### Systematic uncertainties

For all systematic uncertainties we use a unit Gaussian distribution as their
prior.

###  III.2 Oscillation inferences: results and discussion

In this section we describe our inferences regarding the PMNS neutrino
oscillation parameters, given the data, model, and the Bayesian methodology
described above. Two separate sets of results are obtained. The first uses a
Gaussian prior on $\sin^{2}2\theta_{13}$, imposing a constraint from reactor
antineutrino experiments (Sec. III.2.2-III.2.3). The second uses a prior
uniform in $\sin^{2}2\theta_{13}$, yielding results constrained only by the
NOvA data (Sec. III.2.4). Samples obtained from ARIA and Stan produce
essentially identical distributions in all the variables considered (App. D).

####  III.2.1 Goodness of fit

Figure 3: Posterior predictive $p$-values from real data MCMC samples, with
$\theta_{13}$ constraint from reactor experiments applied. The purple
distribution is a scatterplot of the binned $\chi^{2}$ computed between the
model and real data spectra ($x$-axis), against a similar $\chi^{2}$ between
the model and pseudodata spectra ($y$-axis), both divided by the number of
degrees of freedom (DOF) in the fit, computed for each MCMC sample. (See the
text for more on how the pseudodata spectra are constructed.) The dashed and
solid contours show 1 $\sigma$ intervals from the same posterior-predictive
distributions calculated only for $\nu_{e}$ (dark blue, left solid),
$\bar{\nu}_{e}$ (red, right solid), $\nu_{\mu}$ (light blue, left dashed) and
$\bar{\nu}_{\mu}$ (light red, right dashed) data samples. The posterior
predictive $p$-value is the fraction of the distribution that lies above the
diagonal $\chi^{2}_{\textrm{data}}=\chi^{2}_{\textrm{pseudodata}}$ line
(black). The fact that the purple distribution is centered at $(1,1)$ and is
evenly distributed above and below the the diagonal, with the associated
$p$-value of 0.56, indicates a good fit.

To evaluate the goodness of the fit for this Bayesian analysis we use
posterior predictive $p$-values (PPP) [79]. In a PPP test, the coordinates
from each MCMC sample are used to form a prediction for the observed spectra.
A $\chi^{2}$ statistic is computed between each prediction and the data, which
we denote as $\chi^{2}_{\textrm{data}}$. A second prediction is made for each
sample by applying Poisson fluctuations to the prediction above, and a second
$\chi^{2}_{\textrm{pseudodata}}$ calculated between this pseudodata and the
original prediction. Because the data spectra are unchanged in this process,
$\chi^{2}_{\textrm{data}}$ incorporates only variations in the oscillation
parameters and systematic uncertainties. Conversely, since the pseudodata
distributions have Poisson fluctuations applied but use the same oscillation
parameters and systematic uncertainty pulls as the base model for each MCMC
sample, $\chi^{2}_{\textrm{pseudodata}}$ treats only statistical
uncertainties.

The distribution of these
$(\chi^{2}_{\textrm{data}},\chi^{2}_{\textrm{pseudodata}})$ pairs for the
entire ensemble of spectra considered in Fig. 1 is shown as the purple shading
in Fig. 3. The PPP then consists of the fraction of points in this ensemble
that lie above the $\chi^{2}_{\textrm{data}}=\chi^{2}_{\textrm{pseudodata}}$
line. In the limit of infinite MCMC samples, a model that perfectly describes
the data apart from statistical variations will produce a PPP of 0.5. We
observe that the shaded distribution in Fig. 3 is distributed evenly around 1
unit of $\chi^{2}$ per degree of freedom in both axes, and the PPP we obtain
is 0.56. Both of these imply that the model is a good representation of the
data.

We also find good $p$-values for the $\nu_{\mu}$ and $\bar{\nu}_{\mu}$ samples
considered independently, whose 68.3% credible regions are indicated by the
dashed light blue and light red colored contours in Fig. 3, respectively. By
contrast, the effect of fluctuations in the smaller-statistics samples can be
seen more readily in the analogous $\nu_{e}$ (solid dark blue) and
$\bar{\nu}_{e}$ (solid red) contours. Both contours are much larger than their
$\nu_{\mu}$ counterparts, particularly along the y-axis, which corresponds to
the dimension where statistical uncertainties are considered. In the $\nu_{e}$
contour, the offset downwards from unity along the x-axis, and the
corresponding shift above the diagonal along the y-axis, together suggest a
set of fluctuations that are relatively close to the Asimov (unfluctuated)
model prediction. Each member of the ensemble of pseudodata spectra thus
typically has larger $\chi^{2}$ relative to the Asimov than that of the data.
This relative closeness of the $\nu_{e}$ prediction to the data can also be
seen in the top left panel of Fig. 1. The unusual shape of this contour arises
because the PPP distribution for $\nu_{e}$ consists of two modes superimposed
upon one another, corresponding to the high-probability regions within the NO
that will be discussed below in Fig. 7. On the other hand, the situation is
reversed for the $\bar{\nu}_{e}$ spectra, where the data fluctuations result
in larger deviations from the Asimov prediction than the bulk of the
pseudodata spectra. The fact that combining all four subsamples together
produces a PPP that is closer to 0.5 than any of them are individually is good
evidence that the deviations from PPP of 0.5 in the various subsamples are
predominantly driven by statistical, rather than systematic, effects.444The
reader may find further information on interpretations of posterior predictive
$p$-values in Refs. [80, 81]. We conclude therefore that our set of MCMC
samples reflect PMNS parameters with a good description of the physics
exhibited in our data.

####  III.2.2 PMNS parameter measurements

We produced Bayesian credible regions for the PMNS neutrino oscillation
parameters using the data spectra and the MCMC samplers described in sections
III.1 and II.3, respectively. These credible regions together with the
posterior probability distributions are shown in Figs. 4 and 7. In both cases
the reactor $\theta_{13}$ constraint used is the 2019 PDG’s combination of
extant measurements [78] and applied in the form of a Gaussian prior on
$\sin^{2}2\theta_{13}$, as explained in section III.1.3. In the figures in
this section, credible intervals that show the normal and inverted mass
orderings separately are created by first making one shared credible interval
that spans both. The separate panels then display the relevant regions from
this shared interval that apply to the specified ordering. By constructing
them in this manner, we ensure the NO and IO intervals share a highest-
posterior density point and posterior probability distribution, preserving any
NOvA preference towards one of the mass orderings in the credible region
limits. Similarly, the distributions and intervals labeled “both orderings”
are created from MCMC samples over all values for $\Delta m^{2}_{32}$, summing
together the normal and inverted ordering posteriors before extracting the
credible intervals.

Table 1 shows the highest posterior probability points together with the
credible intervals enclosing 68.3% of the posterior (henceforth called the
$1\sigma$ contour in analogy with Gaussian $p$-values in a frequentist
context). These points are given for all the PMNS oscillation parameters of
interest, split into both, normal, and inverted mass orderings (using the same
methodology regarding the mass ordering as for the figures). For some of the
parameters the $1\sigma$ region spans disjoint areas; we denote this with a
union symbol $\cup$. These high posterior probability regions are in generally
good agreement with the frequentist analysis of the same dataset [29].

Table 1: The highest posterior density points (HPD) together with the 1 $\sigma$ Bayesian credible interval limits for the PMNS parameters of interest marginalized over all the mass ordering (MO) hypotheses. Marginalization over the mass orderings is explained at the beginning of Sec. III.2.2. In these results a Gaussian prior corresponding to reactor constraint on $\sin^{2}2\theta_{13}$ (see sec. III.1.3) is applied. | MO | HPD | $1\,\sigma$
---|---|---|---
$\delta_{\textrm{CP}}$ | Both | $0.91\pi$ | $[0.02\pi,0.31\pi]\cup[0.68\pi,1.67\pi]$
Normal | $0.89\pi$ | $[0.54\pi,1.07\pi]\cup[1.99\pi,0.48\pi]$
Inverted | $1.44\pi$ | $[1.26\pi,1.65\pi]$
$\sin^{2}\theta_{23}$ | Both | $0.56$ | $[0.45,0.49]\cup[0.52,0.59]$
Normal | $0.56$ | $[0.44,0.59]$
Inverted | $0.56$ | $[0.55,0.57]$
$\Delta m^{2}_{32}$ | Normal | $2.39$ | $[2.32,2.46]$
($\times 10^{-3}$) | Inverted | $-2.44$ | $[-2.47,-2.41]$
Figure 4: Binned posterior probability densities (shaded) for
$\sin^{2}\theta_{23}$–$\Delta m^{2}_{32}$, marginalized over both mass
orderings and plotted separately for the normal (top) and the inverted
(bottom) mass orderings (marginalization over the mass orderings is explained
at the beginning of Sec. III.2.2). Contours indicate regions enclosing 68.3%,
95.5%, and 99.7% of the posterior probability, and for convenience are labeled
in terms of Gaussian standard deviations $\sigma$ using the corresponding
traditional $z$-scores (i.e., the number of standard deviations away from the
central value).

Figure 4 shows the $\sin^{2}\theta_{23}$–$\Delta m^{2}_{32}$ plane, where the
denser MCMC samples (darker color) and larger credible regions in the upper
panel relative to the lower indicate a mild preference for the normal
ordering. This conclusion holds in the presence of the entire systematic
uncertainty model discussed in Sec. II.3. However, the preferred regions in
each mass ordering depend in a nontrivial way on the systematic uncertainties,
as illustrated in Fig. 5. The most significant effect is in $\Delta
m^{2}_{32}$, where the systematic effects not only broaden the preferred
region, but also shift the central value to larger absolute magnitudes. The
most important of the uncertainties contributing to this movement is in the
absolute calibration of the calorimetric energy scale. Because this directly
affects reconstructed neutrino energies, it shifts the expected number of
events in the trough of the $\nu_{\mu}$ and $\bar{\nu}_{\mu}$ disappearance
spectra and is thus anticorrelated with $\Delta m^{2}_{32}$, as can be seen in
Fig. 6. Shifting $\Delta m^{2}_{32}$ in this way also moves
$\sin^{2}\theta_{23}$ closer to maximal disappearance
($\sin^{2}{\theta_{23}^{\mathrm{MD}}}\approx 0.51$; the value depends slightly
on $\sin^{2}\theta_{13}$). Because the $\nu_{\mu}$ disappearance spectra are
essentially identical for values reflected across the maximal disappearance
line, the credible regions are nearly symmetric around it (the degeneracy is
broken only by the $\nu_{e}$ appearance spectra, which have less statistical
power). Thus, the credible regions for smaller
$|\sin^{2}\theta_{23}-\sin^{2}{\theta_{23}^{\mathrm{MD}}}|$ appear narrower,
even though the sensitivity is unchanged. This effect is responsible for the
apparently slightly tighter constraint on $\sin^{2}\theta_{23}$ observed in
Fig. 5 under the effect of systematic uncertainies.

Figure 5: Credible interval comparisons for $\sin^{2}\theta_{23}$–$\Delta
m^{2}_{32}$ between statistical-only fits and fits including all the NOvA
systematic parameters. The external constraint on $\theta_{13}$ from the
reactor experiments was applied in these fits. Figure 6: Binned posterior
probability density (shaded) with 1, 2, and 3 $\sigma$ credible intervals in
$|\Delta m^{2}_{32}|$ and the NOvA absolute calibration systematic
uncertainty, where the horizontal axis is measured in units of standard
deviations from the nominal value. The external constraint on $\theta_{13}$
from reactor experiments was applied in this fit.

The situation is different in the $\delta_{\textrm{CP}}$–$\sin^{2}\theta_{23}$
plane, which is shown in Fig. 7. Here, we observe preferences for CP-
nonconserving values of $\delta_{\textrm{CP}}$ (i.e., nonintegral values of
$\delta_{\textrm{CP}}/\pi$) in both normal and inverted orderings, and for the
upper octant ($\sin^{2}\theta_{23}>0.5$) of $\theta_{23}$ (reflected also in
Fig. 4), though the CP-conserving points are only weakly disfavored in NO.
However, in contrast to the conclusions for $\sin^{2}\theta_{23}$–$\Delta
m^{2}_{32}$, these inferences are minimally affected by the presence of
systematic uncertainties, as Fig. 8 makes clear—even if the resolution does
degrade slightly and the credible regions grow. Here the minimal impact of
systematic uncertainties owes primarily to the smaller statistics of the
$\nu_{e}$ and especially $\bar{\nu}_{e}$ appearance samples that drive the
sensitivity to these parameters, which results in statistical uncertainties
dominating the uncertainty budget.

Figure 7: Binned posterior probability density (shaded) with 1, 2, and 3
$\sigma$ credible intervals for $\delta_{\textrm{CP}}$–$\sin^{2}\theta_{23}$,
marginalized over both mass orderings for the normal mass ordering (top, blue)
and the inverted mass ordering (bottom, red). In this fit a Gaussian
$\theta_{13}$ prior from reactor experiments was applied. Figure 8: Credible
interval comparisons for $\delta_{\textrm{CP}}$–$\sin^{2}\theta_{23}$ between
statistical-only fits and fits including all the NOvA systematic parameters,
with external constraint on $\theta_{13}$ from reactor experiments.

To more quantitatively assess the mass ordering and octant preferences, we
give the posterior probabilities inferred for each combination of hypotheses
in Table 2. We also express them in a less prior-dependent way using Bayes
factors. In both cases the evidence for either option is weak.555The reader
unfamiliar with the interpretation of Bayes factors is referred to the
standard treatments of Jeffreys [82] or Kass & Raftery [83]. The analogous
Gaussian $p$-values also correspond to significances of less than 1 $\sigma$.
The interpretations of the octant and the mass ordering hypothesis preferences
are in good agreement with the 2020 frequentist analysis of the same dataset,
which used profiling with Feldman–Cousins corrections instead of
marginalization [29, 31]. This general agreement is also true for the
parameters’ intervals, with small differences expected given the differing
statistical methods used. To examine the CP-conservation situation more
comprehensively, taking the whole PMNS matrix into consideration, we will use
the Jarlskog invariant measure, explored in the next section.

Table 2: Bayes factors (posterior probabilities) in all the $\theta_{23}$ octant and mass ordering hypotheses combinations. A weak preference towards the normal mass ordering and upper octant is observed. Numbers extracted from a fit with the external $\theta_{13}$ constraint. Probabilities summed across rows or columns may differ slightly from the totals due to rounding. In these results the reactor constraint on $\sin^{2}2\theta_{13}$ is applied. | Normal | Inverted | Total
---|---|---|---
| Ordering | Ordering
Upper Octant | 0.71 (0.42) | 0.26 (0.21) | 1.67 (0.63)
Lower Octant | 0.35 (0.26) | 0.13 (0.12) | 0.60 (0.38)
Total | 2.08 (0.68) | 0.48 (0.33) | (1.0)

####  III.2.3 CP violation – Jarlskog invariant

The Jarlskog invariant is a measure of the strength of charge-parity violation
that is independent of how the mixing matrix is parameterized. The neutrino-
mixing variant [84] parallels its development in the quark sector [85]. Under
the three-neutrino-flavor assumption, its definition is:

$\displaystyle J\equiv$
$\displaystyle\cos(\theta_{12})\cos^{2}(\theta_{13})\cos(\theta_{23})\sin(\theta_{12})$
$\displaystyle\times\sin(\theta_{13})\sin(\theta_{23})\sin(\delta_{\mathrm{CP}}),$
(3)

Figure 9: Posterior probability density for the Jarlskog invariant in NO (top
plot) and IO (bottom plot). The top half of each panel shows the posterior
with a prior uniform in $\sin\delta_{\textrm{CP}}$, while the bottom half uses
a prior uniform in $\delta_{\textrm{CP}}$. A CP-conserving line is drawn at
$J=0$. The external constraint in $\theta_{13}$ from the reactor experiments
is used.

where if any of the factors is zero, the invariant vanishes and the neutrino
mixing matrix is CP-conserving. Nonzero values of $J$, on the other hand,
indicate CP violation. We produce this measurement by taking the MCMC chain
with all the oscillation parameters’ values and calculating the Jarlskog
invariant at each MCMC step. As NOvA is not sensitive to the value of
$\theta_{12}$, for $J$ we use a Gaussian prior analogous to the one used for
the reactor constraint (sec. III.1.3) but whose range consists of the PDG’s
2019 average of solar and LBL reactor neutrino measurements:
$\sin^{2}\theta_{12}=0.307\pm 0.013$ [78], the same value used in our previous
results.

As described in section III.1.3, we use uninformed priors for the oscillation
parameters where possible, meaning that the priors are uniform in the variable
that is being shown. Since the Jarlskog invariant is written in terms of
$\sin\delta_{\textrm{CP}}$, the natural formulation is in terms of a prior
uniform in $\sin\delta_{\textrm{CP}}$. At the same time, there are theoretical
considerations that suggest a prior uniform in $\delta_{\textrm{CP}}$ may be
more appropriate [86]. We therefore consider both a prior uniform in
$\sin\delta_{\textrm{CP}}$ and one uniform in $\delta_{\textrm{CP}}$. Priors
on the other oscillation parameters were investigated; these showed very weak
effects on the $J$ posterior due to the far higher sensitivity of NOvA data to
the other PMNS parameters constrained above besides $\delta_{\textrm{CP}}$.

Figure 9 shows the inferred values of the Jarlskog invariant extracted from a
fit to the NOvA data with the external constraint on $\theta_{13}$ from the
reactor experiments, marginalized over the normal (top) and inverted (bottom)
mass orderings. The regions most favored by the data tend towards CP violation
in the NO, although the 1 $\sigma$ intervals lie very close to $J=0$. The IO’s
preferred regions are more markedly distant from CP conservation, particularly
for the prior uniform in $\delta_{\textrm{CP}}$. As expected, these trends
mirror those observed when considering CP-conserving and CP-violating values
of $\delta_{\textrm{CP}}$ in Sec. III.2.2.

A more quantitative way of measuring the NOvA preference for CP conservation
with the Jarlskog invariant is through Bayes factors. These can be calculated
with the Savage–Dickey density ratio method, which computes Bayes factors for
point hypotheses [87, 88]. Using this method we can calculate the Bayes factor
for the CP-conserving value $J=0$, nested under the unconstrained hypothesis
where $J$ can take any value. Bayes factors always depend formally on the
choice of prior, but in many circumstances (such as those discussed above), if
two hypotheses being compared use the same prior, it will cancel. This is not
the case in the Savage–Dickey method. Therefore, we compute the Bayes factor
for $J=0$ under both priors for $\delta_{\textrm{CP}}$ considered above. We
then invert it to obtain the associated Bayes factors for CP violation over CP
conservation, $J\neq 0$.

Table 3: Bayes factors for preference of CP violation over CP conservation, extracted using the Savage–Dickey method at $J=0$. Priors uniform in $\delta_{\textrm{CP}}$ and $\sin\delta_{\textrm{CP}}$ are both shown. The preferences are given for the normal (NO), inverted (IO), and both (BO) mass orderings. Prior | NO | IO | BO
---|---|---|---
Uniform $\sin\delta_{\textrm{CP}}$ | 1.3 | 3.5 | 1.5
Uniform $\delta_{\textrm{CP}}$ | 1.1 | 4.4 | 1.5

Table 3 shows the Bayes factors for CP violation against CP conservation for
normal, inverted, and both mass orderings, calculated for priors uniform in
$\sin\delta_{\textrm{CP}}$ or $\delta_{\textrm{CP}}$. All of these
probabilities point towards a preference for CP nonconservation, although the
preferences are minimal regardless of the $\delta_{\textrm{CP}}$ prior or the
mass ordering (and an analogous frequentist $p$-value would indicate
significances less than 1 $\sigma$, apart from the inverted ordering, where
they range $1.1-1.2$ $\sigma$666We emphasize that these significances cannot
be read off of Fig. 9. The probability density shown there is marginalized
independently at each value of $J$. However, the binary hypothesis test $J=0$
vs. $J\neq 0$ requires a simultaneous marginalization across the whole $J\neq
0$ space. For the latter, a point-hypothesis treatment, such as the the
Savage–Dickey formalism we use here, is necessary.).

Figure 10: Posterior probability density (purple shading)—that is, the
probability that a given ($\nu_{e}$, $\bar{\nu}_{e}$) pair of counts is the
true underlying value in nature (see text)—compared to the measurement and
associated expected statistical variance (black cross) of the total number of
$\bar{\nu}_{e}$ candidates vs. $\nu_{e}$ candidates. The black solid contours
indicate the regions enclosing 68.3% of the posterior. Each panel overlays a
set of ellipses corresponding to the predicted event counts over the range of
possible values of $\delta_{\textrm{CP}}$, with the parameter notated nearby
held at that value (blue for NO, red for IO, with varying shading depending on
the parameter value) and other parameter values as given in Table 4. The
triangular markers show the highest posterior density in the
$\nu_{e}$–$\bar{\nu}_{e}$ space for when restricting to NO (blue) and IO (red)
hypotheses. The posteriors shown here employ a uniform prior over
$\sin^{2}2\theta_{13}$.

####  III.2.4 Using only NOvA constraints on $\theta_{13}$

NOvA’s oscillation measurements simultaneously constrain a combination of
mixing angles ($\theta_{23}$, $\theta_{13}$), the mass-squared splitting
$\Delta m^{2}_{32}$, the CP phase $\delta_{\textrm{CP}}$, and the neutrino
mass ordering. Therefore, applying a strong external 1D constraint on
$\theta_{13}$ from reactor antineutrino experiments—and thereby reducing the
available solution space—increases the sensitivity to other parameters. This
is especially true of the octant of $\theta_{23}$, as we will show below.
However, using an uninformed (uniform) prior in $\sin^{2}2\theta_{13}$ enables
us to directly compare a NOvA-only measurement against that of the reactor
experiments. In so doing we can examine the robustness of the PMNS description
across short-baseline reactor antineutrino experiments and long-baseline
accelerator-based oscillation experiments. We also may study whether our
preferences change in the presence of an external constraint on the data.

(a) $\delta_{\textrm{CP}}$–$\sin^{2}\theta_{23}$

(b) $\sin^{2}\theta_{23}$–$\Delta m^{2}_{32}$

Figure 11: Posterior probability density (shaded) with 1, 2, and 3 $\sigma$
credible intervals for $\delta_{\textrm{CP}}$–$\sin^{2}\theta_{23}$ (left) and
$\sin^{2}\theta_{23}$–$\Delta m^{2}_{32}$ (right), marginalized over both mass
orderings, plotted separately for NO (top, blue) and IO (bottom, red).
Contours extracted from a fit with prior uniform in $\sin^{2}2\theta_{13}$.

Without an external constraint on $\sin^{2}2\theta_{13}$, the relevant degrees
of freedom present a complex space with many intercorrelations among the
parameters and numerous possible solutions. We can use the type of posterior
distributions shown in Fig. 1 to explore these possibilities: for example, by
computing the total number of predicted $\nu_{e}$ and $\bar{\nu}_{e}$
candidates for the parameters of each MCMC sample and comparing the
distribution of these predictions in ($\nu_{e}$, $\bar{\nu}_{e}$) candidate
space to what we observe in the data. This is shown in Fig. 10. To guide our
intuition, we overlay ellipses corresponding to the ranges of predicted number
of events for the possible values of $\delta_{\textrm{CP}}$ at fixed values of
the other parameters, subject to constraint from our data in that the highest
posterior density (HPD) value in each ordering is used for parameters not
explicitly varied. We note several important features. First, the highest-
density (highest-probability) region lies essentially equally along the
overlap region between NO and IO ellipses in all panels, meaning we will find
good solutions for either, given an appropriate value of
$\delta_{\textrm{CP}}$. However, the NO ellipses subtend more of the posterior
region, which will result in a slight overall preference for NO. Second, it is
clear from the left panel that the value of $\sin^{2}2\theta_{13}$ preferred
by the reactor average, $\sin^{2}2\theta_{13}=0.085$, is very close to the
highest-probability region for NOvA. Third, comparing the left and center
panels, there is obvious degeneracy between $\sin^{2}2\theta_{13}$ and
$\sin^{2}\theta_{23}$, since independently varying them produces similar
effects on the predictions. And finally, the right panel demonstrates that the
value of $\Delta m^{2}_{32}$ plays a minor (though not negligible) role in the
appearance channel; its primary constraints arise from the $\nu_{\mu}$
disappearance measurement. These features will be examined more quantitatively
in the distributions that follow. Table 4 shows the highest posterior
probability points together with the 1 $\sigma$ credible intervals for all the
PMNS oscillation parameters of interest, extracted from the fit with a prior
uniform in $\sin^{2}2\theta_{13}$. Similarly to Table 1, the values are split
into both, normal, and inverted mass orderings with any disjoint $1\sigma$
regions denoted with a union symbol $\cup$. As compared to the results of the
fit with the external constraint from reactor experiments, the central values
for $\delta_{\textrm{CP}}$ and $\sin^{2}\theta_{23}$ shown here are the most
different. This is not unexpected since there is also a degeneracy between
these two parameters, and both $\delta_{\textrm{CP}}$ and
$\sin^{2}\theta_{23}$ have multiple areas of high probability.

Table 4: The highest posterior density points (HPD) together with the 1 $\sigma$ Bayesian credible interval limits for the PMNS parameters of interest marginalized over all the mass ordering (MO) hypotheses. Marginalization over the mass orderings is explained at the beginning of Sec. III.2.2. Extracted from a fit with prior uniform in $\sin^{2}2\theta_{13}$. | MO | HPD | $1\,\sigma$
---|---|---|---
$\delta_{\textrm{CP}}$ | Both | $1.53\pi$ | $[0.67\pi,1.88\pi]\cup[0.06\pi,0.12\pi]$
Normal | $0.85\pi$ | $[0.48\pi,1.04\pi]\cup[1.99\pi,0.41\pi]$
Inverted | $1.52\pi$ | $[1.25\pi,1.75\pi]$
$\sin^{2}2\theta_{13}$ | Both | $0.087$ | $[0.071,0.107]$
Normal | $0.084$ | $[0.065,0.108]$
Inverted | $0.094$ | $[0.083,0.106]$
$\sin^{2}\theta_{23}$ | Both | $0.46$ | $[0.43,0.50]\cup[0.53,0.58]$
Normal | $0.46$ | $[0.43,0.59]$
Inverted | $0.46$ | $[0.44,0.48]$
$\Delta m^{2}_{32}$ | Normal | $2.39$ | $[2.33,2.46]$
($\times 10^{-3}$) | Inverted | $-2.44$ | $[-2.48,-2.40]$

Fig. 11 shows our preferred regions in
$\delta_{\textrm{CP}}$–$\sin^{2}\theta_{23}$ and $\sin^{2}\theta_{23}$–$\Delta
m^{2}_{32}$ spaces without the external $\theta_{13}$ constraint. Comparing
with Figs. 4 and 7, it is evident that removing the constraint diminishes the
NOvA sensitivity especially to the octant of $\theta_{23}$ (there is now near
symmetry across $\sin^{2}\theta_{23}=0.5$). Although the octant preference is
weakened, the central values of the PMNS parameters $\sin^{2}\theta_{23}$,
$\Delta m^{2}_{32}$, and especially $\delta_{\textrm{CP}}$ do not otherwise
see any significant change. This is as expected from Fig. 10. The results
without the reactor constraint are therefore fully compatible with the
standard NOvA results with the external constraint on $\sin^{2}2\theta_{13}$
applied. We also observe that although the sensitivity is reduced, and
although that for a given value of $\delta_{\textrm{CP}}$ there is always a
combination of parameters within either mass ordering compatible with the
data, certain $\delta_{\textrm{CP}}$–$\sin^{2}\theta_{23}$–ordering
combinations are still excluded with reasonable confidence, such as
$(\delta_{\textrm{CP}}=\frac{\pi}{2},\mathrm{IO})$ and
$(\delta_{\textrm{CP}}=\frac{3\pi}{2},\sin^{2}\theta_{23}=0.5,\mathrm{NO})$.
Combinations such as these produce strong asymmetries in the ($\nu_{e}$,
$\bar{\nu}_{e}$) counts, and as Fig. 10 makes clear, such parameter values lie
in regions outside the posterior point cloud, and are thus disfavored. Table 5
examines the posterior probabilities for each of the octant and mass-ordering
hypotheses more quantitatively. Compared to Table 2, we note a similar
weakened preference for the upper octant of $\theta_{23}$ for the fit without
the reactor constraint.

Table 5: Bayes factors (posterior probabilities) for all the $\theta_{23}$ octant and mass ordering hypotheses, with a marginal preference towards the normal mass ordering and upper octant. Numbers extracted from a fit with a uniform prior in $\sin^{2}2\theta_{13}$. Probabilities summed across rows or columns may differ slightly from the totals due to rounding. | Normal | Inverted | Total
---|---|---|---
| Ordering | Ordering
Upper Octant | 0.53 (0.35) | 0.20 (0.17) | 1.05 (0.51)
Lower Octant | 0.39 (0.28) | 0.42 (0.20) | 0.95 (0.49)
Total | 1.70 (0.63) | 0.59 (0.37) | (1.0)
Figure 12: Comparison of the 2D $\sin^{2}2\theta_{13}$–$\sin^{2}\theta_{23}$
credible intervals (bottom) and 1D $\sin^{2}2\theta_{13}$ posterior
probability densities (top) between both (purple), normal (light-blue) and
inverted (red) mass orderings with a prior uniform in $\sin^{2}2\theta_{13}$
and $\sin^{2}\theta_{23}$. The reactor experiments’ 1 $\sigma$ interval in
$\sin^{2}2\theta_{13}$ from the PDG 2019 [78] is shown in yellow hatched bar.

In Fig. 12 we study the impact of applying the reactor constraint on our mass-
ordering inference. As expected from Fig. 10, the normal ordering contains
more posterior probability, as evidenced by the larger NO contour in the lower
panel. Moreover, as we see in the marginal posterior distribution shown in the
top panel, the posterior restricted to IO prefers generally larger values of
$\sin^{2}2\theta_{13}$. The associated correlation with $\sin^{2}\theta_{23}$
also pushes $\theta_{23}$ into the lower octant, as seen in the lower panel;
this point will be developed further momentarily. Because the reactor value
(the yellow bar indicates its 1 $\sigma$ range) lies at lower
$\sin^{2}2\theta_{13}$, this results in the slightly stronger preference for
the NO in Table 2 as compared to 5. However, the difference is small, which
indicates that the (mild) NO preference observed in the data is largely
independent of the reactor constraint.

The degeneracy in the measurement between $\sin^{2}\theta_{23}$ and
$\sin^{2}2\theta_{13}$ we noted in Fig. 10 can be studied directly by
examining their joint posterior probability distribution, marginalized over
all other parameters and the mass ordering. We show this in Fig. 13. The
central panel exhibits a clear anticorrelation between the octant of
$\theta_{23}$ and the value of $\sin^{2}2\theta_{13}$, which is expected since
both parameters enter at leading order in the $\nu_{\mu}\rightarrow\nu_{e}$
and $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ oscillation probabilities. Here
the overlap of the reactor measurements (again indicated by the yellow hatched
bar) with our marginal posterior for $\sin^{2}2\theta_{13}$ (right panel)
favors the upper octant over the lower octant when we constrain the results to
specifically the upper or lower octant of $\theta_{23}$. Thus, we see that the
preference for the upper octant of $\theta_{23}$ in Table 2 is an emergent
behavior that arises from the application of the reactor constraint. We also
note that though the marginal posterior for $\sin^{2}\theta_{23}$, in the top
panel, shows a higher posterior density in the lower octant, the total
posterior probability integrated across the upper octant is slightly larger
than the corresponding lower octant probability; but as Table 5 makes clear,
this preference is entirely insignificant.

Figure 13: $\sin^{2}\theta_{23}$–$\sin^{2}2\theta_{13}$ posterior probability
densities with 1, 2, and 3 $\sigma$ credible intervals, marginalized over both
mass orderings. The posterior density was extracted from a fit with a prior
uniform in $\sin^{2}2\theta_{13}$. The reactor experiments’ 1 $\sigma$
interval in $\sin^{2}2\theta_{13}$ from the PDG 2019 [78] is shown in the
yellow hatched bar. The right panel shows the posterior probability for
$\sin^{2}2\theta_{13}$, with its contributions from the upper octant (UO,
transparent, red outline) and lower octant (LO, transparent, blue outline) of
$\theta_{23}$.

We emphasize that reactor neutrino experiments and accelerator neutrino
experiments measure the PMNS $\sin^{2}2\theta_{13}$ by examining different
sectors of neutrino oscillations, over a wide range of baselines. Reactor
neutrino experiments measure the $\bar{\nu}_{e}$ survival probability
$\mathcal{P}(\bar{\nu}_{e}\rightarrow\bar{\nu}_{e})$ with low-energy (few-MeV)
$\bar{\nu}_{e}$s over a short (few-km) baseline. Conversely, accelerator
neutrino experiments simultaneously measure $\nu_{e}$ appearance in a few-GeV
$\nu_{\mu}$ beam, and $\bar{\nu}_{e}$ appearance in a $\bar{\nu}_{\mu}$ beam,
both over a long (hundreds-of-km) baseline. In long-baseline measurements
$\mathcal{P}(\nu_{\mu}\rightarrow\nu_{\mu})$,
$\mathcal{P}(\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{\mu})$,
$\mathcal{P}(\nu_{\mu}\rightarrow\nu_{e})$, and
$\mathcal{P}(\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e})$ are all exploited to
constrain the PMNS oscillation parameters, including $\sin^{2}2\theta_{13}$.
Thus, the consistency observed between long- and short-baseline measurements
lends strong support to the PMNS interpretation of neutrino oscillations.

##  IV Conclusions

We have refit the NOvA dataset of $13.6\times 10^{20}$ POT in neutrino beam
mode and $12.5\times 10^{20}$ POT in antineutrino mode using a Bayesian
statistical approach. NOvA data continues to be consistent with maximal mixing
for $\sin^{2}\theta_{23}$ and the regions of the $\sin^{2}\theta_{23}$–$\Delta
m^{2}_{32}$ and $\delta_{\textrm{CP}}$–$\sin^{2}\theta_{23}$ spaces preferred
by this analysis are consistent with those of the previous analysis done using
a frequentist fit.

With the introduction of the Bayesian analysis, we also expand the neutrino
oscillation parameters measured by NOvA. We report for the first time NOvA
measurements that do not require constraints on $\sin^{2}2\theta_{13}$ from
reactor antineutrino oscillations. Moreover, we are also able to include new
results for $\sin^{2}2\theta_{13}$ and the Jarlskog invariant $J$; these were
impractical to produce under the preceding frequentist method due to the
necessity of Feldman–Cousins corrections. The inferences on $J$ provide a
model-independent measurement of CP violation and indicate that NOvA data has
a weak preference for CP violation, which becomes more pronounced when
assuming the inverted mass ordering.

As NOvA measures a convolution of mixing angles, $\delta_{\textrm{CP}}$, and
mass ordering in the oscillation probabilities, removing the external
constraint on $\theta_{13}$ reduces our constraining power to determine the
octant of $\theta_{23}$, mass ordering, and $\delta_{\textrm{CP}}$. While our
sensitivity is reduced without the external constraint (particularly for the
octant), we note that the conclusions arising from the analysis remain
unchanged.

This measurement of $\sin^{2}2\theta_{13}$ using electron neutrinos and
antineutrinos with energies in the GeV range and propagating for hundreds of
kilometers is fully consistent with measurements performed using few-MeV
electron antineutrinos from nuclear reactors propagating for a few kilometers.
The consistency of results using the PMNS framework across a broad regime of
conditions bolsters its applicability. More stringent tests of CP violation
and the consistency of our $\sin^{2}2\theta_{13}$ measurement with those of
reactors will be possible with increased statistics in upcoming NOvA
measurements.

##  V Acknowledgments

This document was prepared by the NOvA collaboration using the resources of
the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of
Energy, Office of Science, HEP User Facility. Fermilab is managed by Fermi
Research Alliance, LLC (FRA), acting under Contract No. DE-AC02-07CH11359.
This work was supported by the U.S. Department of Energy; the U.S. National
Science Foundation; the Department of Science and Technology, India; the
European Research Council; the MSMT CR, GA UK, Czech Republic; the RAS, MSHE,
and RFBR, Russia; CNPq and FAPEG, Brazil; UKRI, STFC and the Royal Society,
United Kingdom; and the state and University of Minnesota. We are grateful for
the contributions of the staffs of the University of Minnesota at the Ash
River Laboratory, and of Fermilab. For the purpose of open access, the author
has applied a Creative Commons Attribution (CC BY) license to any Author
Accepted Manuscript version arising.

##  Appendix A Determining step sizes for ARIA

As noted in Sec. III.1.1, an $\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$ chain is
derived sample-by-sample using a repeated two-step procedure:

  1. 1.

Proposal: The coordinates of a potential new sample are selected from a
probability distribution centered on the current sample (or initial seed).

  2. 2.

Acceptance: The proposal selected above is either accepted or rejected
according to the rule of detailed balance, i.e., that every step in the chain
be exactly reversible.

If accepted, the proposed coordinates become the next sample. If rejected, the
previous sample is repeated to become the next sample.

The $\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$ algorithm does not specify the
distribution to be used in step 1 above, however. In our implementation, we
use the most common choice, a multivariate Gaussian:

$g(\vec{x}\,^{\prime}|\vec{x})=(2\pi)^{-\frac{N}{2}}(\det\Sigma)^{-\frac{1}{2}}\exp\left(-\frac{1}{2}\left(\vec{x}\,^{\prime}-\vec{x}\,\right)^{T}\Sigma^{-1}\left(\vec{x}\,^{\prime}-\vec{x}\,\right)\right),$
(4)

where $\vec{x}$ represents the current sample coordinates,
$\vec{x}\,^{\prime}$ the proposed next coordinates, and $N$ the dimensionality
of the coordinate space. The matrix $\Sigma$ imposes a length scale on the
“distance” between successive samples, and (especially when it is diagonal)
its elements are usually called the “step sizes” of the sampling for each
degree of freedom. The ideal asymptotic fraction of samples accepted in step
2, $\alpha$, is 23.4% under a wide range of circumstances [89, 90]. Though
this figure is strictly true only for $N\rightarrow\infty$, it has been shown
to hold approximately even for parameter counts as low as $N=5$ [91]. Because
the outcome of step 2 is related to the proposals generated in step 1, we
tuned the values of $\Sigma$ to arrive at $\alpha=23.4\%$.

Our overall heuristic in the tuning procedure is to maintain step sizes that
yield similar autocorrelations (defined rigorously below) across all the
parameters. This results in the most efficient exploration of the parameter
space [71]. We first optimized the step sizes for the parameters of interest,
$\theta_{13}$, $\theta_{23}$, $|\Delta m_{32}^{2}|$, and
$\delta_{\textrm{CP}}$. We constructed a chain that sampled only those
parameters using a unit matrix for $\Sigma$. We computed $\alpha$ for this
chain and scaled the relevant elements of $\Sigma$ in order to arrive at a
tolerable preliminary acceptance rate of about 20%. We then computed the
$k$-lag autocorrelation for each parameter $\theta$, which measures the
average correlation between MCMC sample $n$ and sample $n+k$ across all $n$
[92]:

$r_{k}=\frac{\sum_{n=1}^{N-k}{\left(\theta_{n}-\bar{\theta}\right)\left(\theta_{n+k}-\bar{\theta}\right)}}{\sum_{n=1}^{N}{\left(\theta_{n}-\bar{\theta}\right)^{2}}},$
(5)

where $\theta_{n}$ refers to the value of parameter $\theta$ at step $n$, and
$\bar{\theta}$ is its mean value. These autocorrelations are shown in Fig. 14.
Using these $r_{k}$, we further adjusted the elements of $\Sigma$ so the
oscillation parameters would have similar autocorrelations.

Figure 14: Lag autocorrelation (see Eq. 5) computed using a
$\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$ chain sampling only oscillation
parameters.

To optimize the step sizes for the nuisance parameters (systematic
uncertainties), we constructed a chain sampling only those parameters, again
beginning with $\Sigma$ entries of unity for them. As with the oscillation
parameters, we adjusted $\Sigma$ to ensure none of the nuisance parameters had
significantly different autocorrelations from the others. We then constructed
a new, much longer chain, and subsequently computed a covariance matrix over
the nuisance parameters from it. The Cholesky decomposition of this matrix,
$L$, was used in the next step.

A final chain, this time sampling both oscillation and nuisance parameters,
was constructed using the adjusted $\Sigma$ entries for both. While sampling,
the proposed values for the nuisance parameters were multiplied by the
decomposed covariance scaled by a tunable factor, $\beta L$. Using this chain,
we recomputed the autocorrelations for all the parameters. The elements of
$\Sigma$ were readjusted to obtain similar autocorrelations from the
oscillation parameters, now in the presence of the nuisance parameters. We
also adjusted $\beta$ to yield similar autocorrelations to those of the
oscillation parameters. A final global scale was applied to $\Sigma$ and
$\beta$ to finally arrive at $\alpha=23.4\%$. The autocorrelations for the
oscillation parameters at the end of the tuning procedure are shown in Fig.
15.

Figure 15: Lag autocorrelation (see Eq. 5) computed using a
$\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$ chain after step size tuning with all
parameters included.

##  Appendix B Determining warmup and thinning fractions for ARIA

The property of sample proportionality to posterior density in the
$\textrm{M}\textrm{R}^{2}\textrm{T}^{2}$ method is only guaranteed as
asymptotic behavior. Therefore, it is usually necessary to discard some number
of samples $N_{0}$ at the beginning of the Markov chain while the chain “burns
in” or “warms up.” Moreover, it is usually impossible to find a choice of step
sizes $\Sigma$ in Eq. 4 that proposes samples that are both fully uncorrelated
with the previous one and whose acceptance probabilities are high enough to
not impose severe computing requirements. We “thin” our chains by discarding
all but the $k$th sample to reduce autocorrelations. The resulting fraction
$N_{\mathrm{eff}}=\frac{1}{k}(N-N_{0})$, the “number of effective samples,”
corresponds to the statistical power of the chain.

Because autocorrelations in our analysis are relatively long (see Fig. 15) and
ARIA runs fairly quickly, we produced very long chains of $5\times 10^{6}$
samples each. Thus, when we compared our posterior densities using $N_{0}=0$,
$N_{0}=\\{1,3,5\\}\times 10^{3}$, $N_{0}=\\{1,3,5\\}\times 10^{4}$, and
$N_{0}=\\{1,3\\}\times 10^{5}$, we found them all to be indistinguishable.

Proper thinning results in minimal autocorrelations. Though there is no
universally agreed-upon threshold below which autocorrelations are considered
acceptable, we found the plateau in Fig. 15 of around or less than 1% to be
sufficient for our needs. Therefore we chose to thin at $k=10^{4}$.

Figure 16: Posterior probability density in $\Delta m^{2}_{32}$ for normal
ordering (left) and inverted ordering (right). Samples from Stan are shown as
the shaded distributions (1 $\sigma$, 2 $\sigma$, 3 $\sigma$ regions indicated
by darkest to lightest fill). Samples from ARIA are shown as the red (normal
ordering) or blue (inverted ordering) overlaid lines (solid, dashed, dotted
lines represent the boundaries of the 1 $\sigma$, 2 $\sigma$, 3 $\sigma$
regions, respectively).

Figure 17: Posterior probability density in $\sin^{2}\theta_{23}$ (left) and
$\delta_{\textrm{CP}}$ (right), marginalized over both mass orderings. Line
styles as in Fig. 16.

##  Appendix C Algorithm choices for HMCMC

As noted in Sec. III, HMCMC generates proposals by numerically integrating a
Hamiltonian for a fictitious particle, whose potential arises from treating
the log-posterior in analogy to gravity:

$\displaystyle\frac{d\vec{q}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\partial H}{\partial\vec{p}}$
$\displaystyle=\frac{\partial T}{\partial\vec{p}}$ (6)
$\displaystyle\frac{d\vec{p}}{dt}$ $\displaystyle=-$
$\displaystyle\frac{\partial H}{\partial\vec{q}}$
$\displaystyle=-\frac{\partial T}{\partial\vec{q}}-\frac{\partial
V}{\partial\vec{q}}$

where $T$ and $V$ are the kinetic and potential energies of the system,
respectively.

There are two ingredients of HMCMC left unspecified by the method. In both
cases Stan’s default choices were found to be suitable for our needs. The
first is the distribution of kinetic energies from which $T$ in Eq. 6 is
chosen. Stan’s default is the Euclidean-Gaussian kinetic energy distribution:

$T(\vec{q},\vec{p})=\frac{1}{2}\vec{p}\,^{\text{T}}M^{-1}\vec{p}+\log|M|+\text{const.}$
(7)

Here the mass matrix $M$ (analogous to the effect of mass in gravitation) is a
parameter that is automatically inferred by Stan during its warm-up sampling
by iterative adjustments based on a running covariance over the samples. The
second implementation choice is how long the integrator is allowed to run for
each particular trajectory. Stan uses an algorithm called No-U-Turns (NUTS)
[93], which is a heuristic method that halts integration when two
trajectories, extending in each direction from the starting point along the
initial momentum, begin to converge towards one another. An upper limit of
integrator steps is also supplied as a parameter to Stan; in this analysis, we
find that all our trajectories end within $2^{11}$ steps.

##  Appendix D Equivalence of ARIA and Stan results

We extracted the posterior distributions for all of the results shown in this
paper using both the ARIA and Stan samplers described in Sec. III. The
posteriors were in all cases nearly indistinguishable, with tiny differences
that occasionally caused the boundaries of credible intervals to shift by
single bins. This is illustrated in Figs. 16-17.

## References

  * Fukuda  _et al._ [1998] Y. Fukuda _et al._ (Super-Kamiokande), Phys. Rev. Lett. 81, 1562 (1998), arXiv:hep-ex/9807003 [hep-ex] .
  * Fukuda  _et al._ [2002] S. Fukuda _et al._ (Super-Kamiokande), Phys. Lett. B 539, 179 (2002), arXiv:hep-ex/0205075 [hep-ex] .
  * Ahmad  _et al._ [2002] Q. R. Ahmad _et al._ (SNO), Phys. Rev. Lett. 89, 011301 (2002), arXiv:nucl-ex/0204008 [nucl-ex] .
  * Eguchi  _et al._ [2003] K. Eguchi _et al._ (KamLAND), Phys. Rev. Lett. 90, 021802 (2003), arXiv:hep-ex/0212021 [hep-ex] .
  * Michael  _et al._ [2006] D. G. Michael _et al._ (MINOS), Phys. Rev. Lett. 97, 191801 (2006), arXiv:hep-ex/0607088 [hep-ex] .
  * Abe  _et al._ [2011] K. Abe _et al._ (T2K), Phys. Rev. Lett. 107, 041801 (2011), arXiv:1106.2822 [hep-ex] .
  * Abe  _et al._ [2012] Y. Abe _et al._ (Double Chooz), Phys. Rev. Lett. 108, 131801 (2012), arXiv:1112.6353 [hep-ex] .
  * An  _et al._ [2012] F. P. An _et al._ (Daya Bay), Phys. Rev. Lett. 108, 171803 (2012), arXiv:1203.1669 [hep-ex] .
  * Ahn  _et al._ [2012] J. K. Ahn _et al._ (RENO), Phys. Rev. Lett. 108, 191802 (2012), arXiv:1204.0626 [hep-ex] .
  * Mohapatra and Smirnov [2006] R. N. Mohapatra and A. Y. Smirnov,  _Elementary particle physics. Proceedings, Corfu Summer Institute, CORFU2005, Corfu, Greece, September 4-26, 2005_ , Ann. Rev. Nucl. Part. Sci. 56, 569 (2006), arXiv:hep-ph/0603118 [hep-ph] .
  * Nunokawa  _et al._ [2008] H. Nunokawa, S. J. Parke, and J. W. F. Valle, Prog. Part. Nucl. Phys. 60, 338 (2008), arXiv:0710.0554 [hep-ph] .
  * Altarelli and Feruglio [2010] G. Altarelli and F. Feruglio, Rev. Mod. Phys. 82, 2701 (2010), arXiv:1002.0211 [hep-ph] .
  * King [2015] S. F. King, J. Phys. G 42, 123001 (2015), arXiv:1510.02091 [hep-ph] .
  * Petcov [2018] S. T. Petcov, Eur. Phys. J. C 78, 709 (2018), arXiv:1711.10806 [hep-ph] .
  * Fukugita and Yanagida [1986] M. Fukugita and T. Yanagida, Phys. Lett. B 174, 45 (1986).
  * Buchmuller and Plumacher [1996] W. Buchmuller and M. Plumacher, Phys. Lett. B 389, 73 (1996), arXiv:hep-ph/9608308 .
  * Buchmuller  _et al._ [2005a] W. Buchmuller, P. Di Bari, and M. Plumacher, Annals Phys. 315, 305 (2005a), arXiv:hep-ph/0401240 .
  * Buchmuller  _et al._ [2005b] W. Buchmuller, R. D. Peccei, and T. Yanagida, Ann. Rev. Nucl. Part. Sci. 55, 311 (2005b), arXiv:hep-ph/0502169 .
  * Pilaftsis [1997] A. Pilaftsis, Phys. Rev. D 56, 5431 (1997), arXiv:hep-ph/9707235 .
  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Wolfenstein [1978] L. Wolfenstein, Phys. Rev. D 17, 2369 (1978).
  * Esteban  _et al._ [2020] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, JHEP 09, 178, arXiv:2007.14792 [hep-ph] .
  * Adamson  _et al._ [2016a] P. Adamson _et al._ (NOvA), Phys. Rev. D 93, 051104 (2016a), arXiv:1601.05037 [hep-ex] .
  * Adamson  _et al._ [2016b] P. Adamson _et al._ (NOvA), Phys. Rev. Lett. 116, 151806 (2016b), arXiv:1601.05022 [hep-ex] .
  * Adamson  _et al._ [2017a] P. Adamson _et al._ (NOvA), Phys. Rev. Lett. 118, 151802 (2017a), arXiv:1701.05891 [hep-ex] .
  * Adamson  _et al._ [2017b] P. Adamson _et al._ (NOvA), Phys. Rev. Lett. 118, 231801 (2017b), arXiv:1703.03328 [hep-ex] .
  * Acero  _et al._ [2018] M. A. Acero _et al._ (NOvA), Phys. Rev. D 98, 032012 (2018), arXiv:1806.00096 [hep-ex] .
  * Acero  _et al._ [2019] M. A. Acero _et al._ (NOvA), Phys. Rev. Lett. 123, 151803 (2019), arXiv:1906.04907 [hep-ex] .
  * Acero  _et al._ [2022a] M. A. Acero _et al._ (NOvA), Phys. Rev. D 106, 032004 (2022a), arXiv:2108.08219 [hep-ex] .
  * Feldman and Cousins [1998] G. J. Feldman and R. D. Cousins, Phys. Rev. D 57, 3873 (1998), arXiv:physics/9711021 .
  * Acero  _et al._ [2022b] M. A. Acero _et al._ (NOvA), Submitted to Phys. Rev. D (2022b), arXiv:2207.14353 [hep-ex] .
  * Adamson  _et al._ [2016c] P. Adamson _et al._ , Nucl. Instrum. Methods Phys Res., Sect. A 806, 279 (2016c), arXiv:1507.06690 [hep-ex] .
  * Agostinelli  _et al._ [2003] S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A 506, 250 (2003).
  * Geant4 Collaboration [2017] Geant4 Collaboration, geant4-data.web.cern.ch, https://geant4-data.web.cern.ch/ ReleaseNotes/ReleaseNotes4.10.4.html (2017).
  * Aliaga  _et al._ [2016] L. Aliaga _et al._ (MINERvA), Phys. Rev. D 94, 092005 (2016), [Addendum: Phys. Rev.D95,no.3,039903(2017)], arXiv:1607.00704 [hep-ex] .
  * Paley  _et al._ [2014] J. M. Paley _et al._ (MIPP), Phys. Rev. D90, 032001 (2014), arXiv:1404.5882 [hep-ex] .
  * Alt  _et al._ [2007] C. Alt _et al._ (NA49), Eur. Phys. J. C49, 897 (2007), arXiv:hep-ex/0606028 [hep-ex] .
  * Abgrall  _et al._ [2011] N. Abgrall _et al._ (NA61/SHINE), Phys. Rev. C84, 034604 (2011), arXiv:1102.0983 [hep-ex] .
  * Barton  _et al._ [1983] D. S. Barton _et al._ , Phys. Rev. D27, 2580 (1983).
  * Seun [2007] S. M. Seun,  _Measurement of $\pi-K$ ratios from the NuMI target_, Ph.D. thesis, Harvard U. (2007).
  * Tinti [2010] G. M. Tinti,  _Sterile neutrino oscillations in MINOS and hadron production in pC collisions_ , Ph.D. thesis, Oxford U. (2010).
  * Lebedev [2007] A. V. Lebedev,  _Ratio of pion kaon production in proton carbon interactions_ , Ph.D. thesis, Harvard U. (2007).
  * Baatar  _et al._ [2013] B. Baatar _et al._ (NA49), Eur. Phys. J. C73, 2364 (2013), arXiv:1207.6520 [hep-ex] .
  * Skubic  _et al._ [1978] P. Skubic _et al._ , Phys. Rev. D18, 3115 (1978).
  * Denisov  _et al._ [1973] S. P. Denisov, S. V. Donskov, Yu. P. Gorin, R. N. Krasnokutsky, A. I. Petrukhin, Yu. D. Prokoshkin, and D. A. Stoyanova, Nucl. Phys. B61, 62 (1973).
  * Carroll  _et al._ [1979] A. S. Carroll _et al._ , Phys. Lett. 80B, 319 (1979).
  * Abe  _et al._ [2013] K. Abe _et al._ (T2K), Phys. Rev. D87, 012001 (2013), [Addendum: Phys. Rev. D87, no.1, 019902 (2013)], arXiv:1211.0469 [hep-ex] .
  * Gaisser  _et al._ [1975] T. K. Gaisser, G. B. Yodh, V. D. Barger, and F. Halzen, in  _14th International Cosmic Ray Conference (ICRC 1975) Munich, Germany, August 15-29, 1975_ (1975) pp. 2161–2166.
  * Cronin  _et al._ [1957] J. W. Cronin, R. Cool, and A. Abashian, Phys. Rev. 107, 1121 (1957).
  * Allaby  _et al._ [1969] J. V. Allaby _et al._ (IHEP-CERN), Phys. Lett. 30B, 500 (1969).
  * Longo and Moyer [1962] M. J. Longo and B. J. Moyer, Phys. Rev. 125, 701 (1962).
  * Bobchenko  _et al._ [1979] B. M. Bobchenko _et al._ , Sov. J. Nucl. Phys. 30, 805 (1979), [Yad. Fiz. 30, 1553 (1979)].
  * Fedorov  _et al._ [1978] V. B. Fedorov, Yu. G. Grishuk, M. V. Kosov, G. A. Leksin, N. A. Pivnyuk, S. V. Shevchenko, V. L. Stolin, A. V. Vlasov, and L. S. Vorobev, Sov. J. Nucl. Phys. 27, 222 (1978), [Yad. Fiz.27,413(1978)].
  * Abrams  _et al._ [1970] R. J. Abrams, R. L. Cool, G. Giacomelli, T. F. Kycia, B. A. Leontic, K. K. Li, and D. N. Michael, Phys. Rev. D1, 1917 (1970).
  * Andreopoulos  _et al._ [2010] C. Andreopoulos _et al._ , Nucl. Instrum. Meth. A 614, 87 (2010), arXiv:0905.2517 [hep-ph] .
  * Andreopoulos  _et al._ [2015] C. Andreopoulos, C. Barry, S. Dytman, H. Gallagher, T. Golan, R. Hatcher, G. Perdue, and J. Yarba, The GENIE Neutrino Monte Carlo Generator: Physics and User Manual (2015), arXiv:1510.05494 [hep-ph] .
  * Nieves  _et al._ [2011] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Rev. C 83, 045501 (2011), arXiv:1102.2777 [hep-ph] .
  * Gran  _et al._ [2013] R. Gran, J. Nieves, F. Sanchez, and M. J. Vicente Vacas, Phys. Rev. D 88, 113007 (2013), arXiv:1307.8105 [hep-ph] .
  * Allardyce  _et al._ [1973] B. W. Allardyce _et al._ , Nucl. Phys. A 209, 1 (1973).
  * Saunders  _et al._ [1996] A. Saunders, S. Hoeibraten, J. J. Kraushaar, B. J. Kriss, R. J. Peterson, R. A. Ristinen, J. T. Brack, G. Hofman, E. F. Gibson, and C. L. Morris, Phys. Rev. C 53, 1745 (1996).
  * Meirav  _et al._ [1989] O. Meirav, E. Friedman, R. R. Johnson, R. Olszewski, and P. Weber, Phys. Rev. C 40, 843 (1989).
  * Levenson  _et al._ [1983] S. M. Levenson _et al._ , Phys. Rev. C 28, 326 (1983).
  * Ashery  _et al._ [1981] D. Ashery, I. Navon, G. Azuelos, H. K. Walter, H. J. Pfeiffer, and F. W. Schleputz, Phys. Rev. C 23, 2173 (1981).
  * Ashery  _et al._ [1984] D. Ashery _et al._ , Phys. Rev. C 30, 946 (1984).
  * Pinzon Guerra  _et al._ [2017] E. S. Pinzon Guerra _et al._ (DUET), Phys. Rev. C 95, 045203 (2017), arXiv:1611.05612 [hep-ex] .
  * Aurisano  _et al._ [2015] A. Aurisano, C. Backhouse, R. Hatcher, N. Mayer, J. Musser, R. Patterson, R. Schroeter, and A. Sousa (NOvA),  _Proceedings, 21st International Conference on Computing in High Energy and Nuclear Physics (CHEP 2015): Okinawa, Japan, April 13-17, 2015_ , J. Phys. Conf. Ser. 664, 072002 (2015).
  * Aurisano  _et al._ [2016] A. Aurisano, A. Radovic, D. Rocco, A. Himmel, M. D. Messier, E. Niner, G. Pawloski, F. Psihas, A. Sousa, and P. Vahle, JINST 11 (09), P09001, arXiv:1604.01444 [hep-ex] .
  * Psihas  _et al._ [2019] F. Psihas, E. Niner, M. Groh, R. Murphy, A. Aurisano, A. Himmel, K. Lang, M. D. Messier, A. Radovic, and A. Sousa, Phys. Rev. D 100, 073005 (2019), arXiv:1906.00713 [physics.ins-det] .
  * Mr. Bayes and Mr. Price [1763] Mr. Bayes and Mr. Price, Philosophical Transactions (1683-1775) 53, 370 (1763).
  * Speagle [2020] J. S. Speagle, A conceptual introduction to markov chain monte carlo methods (2020), arXiv:1909.12313 [stat.OT] .
  * Brooks  _et al._ [2011] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng,  _Handbook of Markov Chain Monte Carlo_ (2011) pp. 1–592.
  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
  * Hastings [1970] W. K. Hastings, Biometrika 57, 97 (1970).
  * Gubernatis [2005] J. E. Gubernatis, Phys. Plasmas 12, 057303 (2005).
  * [75] M. Rosenbluth, AIP Conference Proceedings 690, 22.
  * Stan Development Team [2021] Stan Development Team,  _Stan Modeling Language Users Guide and Reference Manual, version 2.26.1_ , Stan Development Team (2021).
  * [77] M. Betancourt, arXiv:1701.02434 [stat.ME] .
  * Tanabashi  _et al._ [date] M. Tanabashi _et al._ (Particle Data Group), Phys. Rev. D 98, 030001 (2018 and 2019 update).
  * Gelman  _et al._ [1996] A. Gelman, X.-L. Meng, and H. Stern, Statistica Sinica 6, 733 (1996).
  * Gelman  _et al._ [2002] A. Gelman, Y. Goegebeur, F. Tuerlinckx, and I. v. Van Mechelen, Journal of the Royal Statistical Society Series C: Applied Statistics 49, 247 (2002).
  * Gelman [2013] A. Gelman, Electronic Journal of Statistics 7, 2595 (2013).
  * Jeffreys [1961] H. Jeffreys,  _Theory of Probability_ , 3rd ed. (Oxford, Oxford, England, 1961).
  * Kass and Raftery [1995] R. E. Kass and A. E. Raftery, Journal of the American Statistical Association 90, 773 (1995).
  * Denton and Parke [2019] P. B. Denton and S. J. Parke, Phys. Rev. D 100, 053004 (2019).
  * Jarlskog [1985] C. Jarlskog, Phys. Rev. Lett. 55, 1039 (1985).
  * Denton and Pestes [2021] P. B. Denton and R. Pestes, JHEP 05, 139, arXiv:2006.09384 [hep-ph] .
  * Dickey [1971] J. M. Dickey, The Annals of Mathematical Statistics 42, 204 (1971).
  * Mulder  _et al._ [2020] J. Mulder, E. J. Wagenmakers, and M. Marsman 10.48550/arxiv.2004.09899 (2020).
  * Gelman  _et al._ [1997] A. Gelman, W. R. Gilks, and G. O. Roberts, The Annals of Applied Probability 7, 110 (1997).
  * Roberts [1998] G. O. Roberts, Stochastics and Stochastic Reports 62, 275 (1998).
  * Roberts and Rosenthal [2001] G. O. Roberts and J. S. Rosenthal, Statistical Science 16, 351 (2001).
  * Croarkin and Tobias [2013] C. Croarkin and P. Tobias, eds.,  _NIST/SEMATECH e-Handbook of Statistical Methods_ (2013).
  * Hoffman and Gelman [2014] M. D. Hoffman and A. Gelman, Journal of Machine Learning Research 15, 1593 (2014).