New constraints on oscillation parameters from $\nu_e$ appearance and $\nu_\mu$ disappearance in the NOvA experiment

We present updated results from the NOvA experiment for $\nu_\mu\rightarrow\nu_\mu$ and $\nu_\mu\rightarrow\nu_e$ oscillations from an exposure of $8.85\times10^{20}$ protons on target, which represents an increase of 46% compared to our previous publication. The results utilize significant improvements in both the simulations and analysis of the data. A joint fit to the data for $\nu_\mu$ disappearance and $\nu_e$ appearance gives the best fit point as normal mass hierarchy, $\Delta m^2_{32} = 2.44\times 10^{-3}{{\rm eV}^2}/c^4$, $\sin^2\theta_{23} = 0.56$, and $\delta_{CP} = 1.21\pi$. The 68.3% confidence intervals in the normal mass hierarchy are $\Delta m^2_{32} \in [2.37,2.52]\times 10^{-3}{{\rm eV}^2}/c^4$, $\sin^2\theta_{23} \in [0.43,0.51] \cup [0.52,0.60]$, and $\delta_{CP} \in [0,0.12\pi] \cup [0.91\pi,2\pi]$. The inverted mass hierarchy is disfavored at the 95% confidence level for all choices of the other oscillation parameters.


INTRODUCTION
Joint fits of ν µ → ν µ disappearance and ν µ → ν e appearance oscillations in long-baseline neutrino oscillation experiments can provide information on four of the standard neutrino model parameters, |∆m 2  32 |, θ 23 , δ CP , and the mass hierarchy, when augmented by measurements of the other three parameters, ∆m 2  21 , θ 12 , and θ 13 , from other experiments [1].Of the four parameters, the first pair are currently most sensitively measured by ν µ → ν µ oscillations and the second pair are most sensitively measured by ν µ → ν e oscillations.However, the precision with which ν µ → ν e oscillations can measure the second pair of parameters depends on the precision of the measurement of θ 23 since that oscillation probability is largely proportional to sin 2 2θ 13 sin 2 θ 23 .
The determination of the neutrino mass hierarchy is important both for grand unified models [3] and for the interpretation of neutrinoless double beta decay experiments [4].In long-baseline neutrino experiments, it is measured by observing the effect of coherent forward neutrino scattering from electrons in the earth, which enhances ν µ → ν e oscillations for the normal mass hierarchy (NH), ∆m 2  32 > 0, and suppresses them for the inverted mass hierarchy (IH), ∆m 2 32 < 0. For the baselines of current experiments and for fixed baseline length to energy ratio, the magnitude of this effect is approximately proportional to the length of the baseline.
The amount of CP violation in the lepton sector is proportional to | sin δ CP |.For δ CP in the range 0 to 2π, ν µ → ν e oscillations are enhanced for δ CP > π and suppressed for δ CP < π, with maximal enhancement at δ CP = 3π/2 and maximal suppression at δ CP = π/2.
In addition to the NOvA results [5], previous joint fits of ν µ → ν µ and ν µ → ν e oscillations in long-baseline experiments have been reported by the MINOS [6] and T2K [7] experiments.
The data reported here correspond to the equivalent of 8.85×10 20 protons on target (POT) in the full NOvA Far Detector with a beam line set to focus positively charged mesons, which greatly enhances the neutrino to antineutrino ratio.This represents a 46% increase in neutrino flux since our last publication [5].These data were taken between February 6, 2014 and February 20, 2017.
Significant improvements have been made to both the simulations and data analysis.The key updates to the simulations include a new data-driven neutrino flux model, an improved treatment of multi-nucleon interactions, and an improved light model including Cherenkov radiation in the scintillator.The main improvements in the ν µ disappearance data analysis are the use of a deeplearning event classifier and the separation of selected events into different samples based on their energy res-olution.The main improvement for the ν e appearance data analysis is the addition of a signal-rich sample that expands the active volume considered.

II. THE NOvA EXPERIMENT
NOvA [8] is a two-detector, long-baseline neutrino oscillation experiment that samples the Fermilab NuMI neutrino beam [9] approximately 1 km from the source using a Near Detector (ND) and observes the oscillated beam 810 km downstream with a Far Detector (FD) near Ash River, MN.The detectors are functionally identical, scintillating tracker-calorimeters consisting of layered reflective polyvinyl chloride cells filled with a liquid scintillator comprised primarily of mineral oil with a 5% pseudocumene admixture.These cells are organized into planes alternating in vertical and horizontal orientation.The net composition of the detectors is 63% active material by mass.Light produced within a cell is collected using a loop of wavelength-shifting optical fiber, which is connected to an avalanche photodiode (APD).
The FD cells are 3.9 × 6.6 cm in cross section, with the 6.6 cm dimension along the beam direction, and 15.5 m long [10].The FD contains 896 planes, leading to a total mass of 14 kt.The majority of ND cells are identical to those of the FD apart from being shorter (3.9 m long instead of 15.5 m).To improve muon containment, the downstream end of the ND is a "muon catcher" composed of a stack of sets of planes in which a pair of one vertically-oriented and one horizontally-oriented scintillator plane is interleaved with one 10 cm-thick plane of steel.There are 11 pairs of scintillator planes separated by 10 steel planes in this sequence.The vertical planes in this section are 2.6 m high.The ND consists of 214 planes for a total mass of 290 ton.
The FD sits 14.6 mrad away from the central axis of the NuMI beam.This off-axis location results in a neutrino flux with a narrow-band energy spectrum centered around 1.9 GeV in the FD.Such a spectrum emphasizes ν µ → ν e oscillations at this baseline and reduces backgrounds from higher energy neutral current events.The ND sees a line source and so it receives a much larger spread in off-axis angles than the FD does.The ND is positioned at the same average off-axis angle as the FD to maximize the similarity between the neutrino energy spectrum at its location and that expected at the FD in the absence of oscillations.
The beam is pulsed at an average rate of 0.75 Hz.All of the APD signals above threshold from a large time window around each 10 µs beam spill are retained.Because the FD is located on the Earth's surface, it is exposed to a substantial cosmic ray flux, which is only partially mitigated by its overburden of 1.2 m of concrete plus 15 cm of barite.Therefore, we also use cosmic data taken from 420 µs surrounding the beam spill within beam triggers to obtain a direct measure of the cosmic background in the FD.Separate periodic minimum-bias triggers of the  same length as the beam trigger allow us to collect highstatistics cosmic data for algorithm training and calibration purposes.As the ND is 100 m underground, the cosmic ray flux there is negligible.

III. SIMULATIONS
To assist in calibrating our detectors, determining our analysis criteria, and inferring extracted physical parameters, we rely on predictions generated by a comprehensive simulation suite, which proceeds in stages.We begin by using geant4 [11] and a detailed model of the beamline geometry to simulate the production of hadrons arising from the collision of the 120 GeV primary proton beam with the graphite target [12], as well as their subsequent focusing and decay into neutrinos.The resultant neutrino flux is corrected according to constraints on the hadron spectrum from thin-target hadroproduction data using the ppfx tools developed for the NuMI beam by the MINERvA collaboration [13].The correction applied to the underlying model used in the simulation (FTFP BERT) is in the order of 7-10% for both, the ν µ and ν e flux predictions.The uncertainties are in the order of 8% in the peak.Table I shows simulated predictions of the beam composition at the Near and Far Detectors in the absence of oscillations; the ND predicted spectra from 0-20 GeV are shown in Fig. 1.The predicted flux is then used as input to genie [14,15], which simulates neutrino reactions in the variety of materials of which our detectors and their surroundings are composed.We alter its default interaction model as described below.Finally, we use a detailed model of our detectors with a combination of geant4 and custom software to simulate the detectors' photon response to particles outgoing from individual predicted neutrino reactions, including both scintillation and Cherenkov radiation in the active detector materials, as well as the light transport, collection, and digitization processes.The overall energy scales of both detectors are calibrated using the minimum-ionizing portions of stopping cosmic ray muon tracks.
As in our previous results [5,16,17], we augment genie's default configuration by enabling its semi-empirical model for Meson Exchange Current (MEC) interactions [18] to account for the likely presence of interactions of neutrinos with nucleon-nucleon pairs in both chargedand neutral-current reactions.However, in this analysis we no longer reweight the momentum transfer distributions produced by this model, preferring instead to allow fits to the FD data to profile [19] over the substantially improved systematic uncertainty treatment for this component of the model, as described in Sec.V.In our central-value prediction we simply increase the rate of MEC interactions by 20% as suggested by fits to the sample of ND ν µ CC candidate events in our ND data.In addition, we now reweight the output of the default model for quasielastic production to treat the expected effect of long-range nuclear charge screening according to the Random Phase Approximation (RPA) calculations of J. Nieves and collaborators [20,21].Lastly, we continue to reduce the rate of ν µ CC nonresonant single pion production with invariant hadronic mass W < 1.7 GeV to 41% of genie's nominal value [22].

IV. DATA ANALYSIS
In order to infer the oscillation parameters from our data, we compare the spectra observed at the FD with our predictions under various oscillation hypotheses.This process consists of three steps.First, we develop selections to retain ν e and ν µ charged-current (CC) events and to reject neutral-current (NC) events and cosmogenic activity.Second, we apply the relevant subset of these selections (excluding, e.g., cosmic rejection criteria) to samples observed at the ND, where both ν µ disappearance and ν e appearance are negligible, to constrain our prediction for the selected sample composition.Finally, we combine the constrained prediction from the previous step with the predicted ratio of the FD and ND spectra, which accounts for geometric differences between the detectors, the beam dispersion, and the effect of oscillations.The result is used in fits to the neutrino energy spectra of the candidates observed at the FD.The following sections discuss how this procedure unfolds for each of the two analyses separately.

Event selection
Isolation of samples of candidate events begins with cells whose APD responses are above threshold, known as hits; those neighboring each other in space and time are clustered to produce candidate neutrino events [23,24].We pass hits in event candidates that survive basic quality cuts in timing (relative to the 10 µs beam spill), containment, and contiguity into a deep-learning classifier known as the Convolutional Visual Network (CVN) [25].CVN applies a series of linear operations, trained over simulated beam and cosmic data event samples, which extract complex, abstract visual features from each event, in a scheme based on techniques from computer vision [26,27].The final step of the classifier is a multilayer perceptron [28,29] that maps the learned features onto a set of normalized classification scores, which range over beam neutrino event hypotheses (ν e CC, ν µ CC, ν τ CC, and NC) and a cosmogenic hypothesis.We retain events whose CVN score for the ν µ CC hypothesis exceeds a tuned threshold.
To identify the muon in such events, tracks produced by a Kalman filter algorithm [30][31][32] are scored by a knearest neighbor classifier [33] over the following variables: likelihoods in dE/dx and scattering constructed from single-particle hypotheses, total track length, and the fraction of planes along the track consistent with having minimum-ionizing-like dE/dx.The most muon-like of these tracks is taken to be the muon candidate.Events that have no sufficiently muon-like track are rejected.We also discard events where any clusters of activity extend to the edges of the detector or where any track besides the muon candidate penetrates into the muon catcher in the ND.To avoid being considered as cosmogenic, FD events must furthermore be deemed sufficiently signallike by a boosted decision tree (BDT) [34] trained over simulation and cosmic data that considers the positions, directions, and lengths of tracks, as well as the fraction of the event's total hit count associated with the track and the CVN score for the cosmic hypothesis.According to our simulation, the FD selection efficiency for our basic quality and containment cuts, relative to all true ν µ CC events within a fiducial volume, is 41.3%; the efficiency of the CVN and PID constraints applied to the quality-andcontainment sample is 78.1%.The final selected sample is 92.7% ν µ CC.The predicted composition of the sample at various stages in the selection is given in Table II.

Energy estimation and analysis binning
We reconstruct each event's neutrino energy E ν using a function of the muon candidate and hadronic remnant energies, which are estimated separately.The muon candidate energy E µ is determined from the range of the track, calibrated to true muon energy in our simulation.
We estimate the energy of the hadronic component with a mapping of observed non-muon energy to true non-muon energy also calibrated with the simulation [35].The resulting neutrino energy resolution over the whole sample is 9.1% at the FD (11.8% at the ND due to the lower active fraction of the muon catcher) for ν µ CC events.
The precision with which we can measure sin 2 2θ 23 and ∆m 2  32 depends on the ν µ energy resolution, particularly for events near the disappearance maximum, about 1.6 GeV at the NOvA baseline.Accordingly, we optimize the binning in two ways to get the best effective use of our energy resolution.First, we employ a variable neutrino energy binning with finer bins near the disappearance maximum and coarser bins elsewhere.And, second, we further divide the event populations in each energy bin into four populations in reconstructed hadronic energy fraction, (E ν − E µ )/E ν , which correspond to regions of different neutrino energy resolution [36].These divisions are chosen such that the FD populations are of equal size in the unoscillated simulation; however, the boundaries show little sensitivity to the choice of oscillation parameters.Grouping in this manner has the additional advantage of isolating most background cosmic and beam NC events (those typically mistaken for signal events with energetic hadronic systems) along with events of worst energy resolution into a separate quartile from the three quartiles containing the signal events with better resolution.The average ν µ energy resolution in the FD across the whole energy spectrum is estimated to be 6.2%, 8.2%, 10.3%, and 12.3% for each quartile, respectively.
Figure 2 shows a comparison of the reconstructed neutrino energy for the selected ν µ CC events in the ND with simulation shown area-normalized to the data.The means of the distributions agree to within 10 MeV (0.6%).Normalizing the prediction by area removes a 1.3% normalization difference between the data and the simulation and suppresses 10-20% absolute normalization uncertainties due primarily to our knowledge of the neutrino flux and normalization offsets from cross-section uncertainties.The remaining uncertainties arise from shape differences.The full set of uncertainties that are used to compute the error band is described in Sec.V. Figure 3 shows the corresponding distributions divided into the quartiles.
Reconstructed Neutrino Energy (GeV) Comparison of the reconstructed neutrino energy for selected νµ CC events (black dots) in the ND with area-normalized simulation (red line).Shading represents the bin-to-bin systematic uncertainties.The gray area, which is nearly indistinguishable from the lower figure boundary, shows the simulated background.

Constraints from the Near Detector data
As in our previous work [16], we obtain a data-driven estimate for the true neutrino energy spectrum using our observed ND data.To do so, we reweight the simulation in each reconstructed neutrino energy bin to obtain agreement with the ND data, thus correcting the differences observed in Fig. 3.After subtracting the expected background, which is minimal, we pass the resulting reconstructed neutrino energy spectrum through the migration matrix between reconstructed and true neutrino energies predicted by our ND simulation.The corrected prediction is then multiplied by the predicted bin-by-bin ratios of the FD and ND true energy spectra, which includes the effects of differing detector geometries and acceptances, beam divergence, and three-flavor oscillations, to obtain an expected FD true energy spectrum.The latter is finally converted back to reconstructed energy by way of the analogous FD migration matrix.This constrained signal prediction is summed together with the cosmic prediction, whose reconstructed energy distribution is extracted using events in the minimum-bias trigger passing all the selection criteria and normalized using the 420 µs window around the beam bunch, and a simulationbased beam background prediction to compare to the observed FD data.In the current analysis, this extrapolation procedure is performed within each hadronic energy fraction range separately so that neutrino reaction types that favor different regions of the elastic-to-inelastic continuum (and thereby have typically different neutrino energy resolution) can be constrained independently.We find the total number of events in each of the four quartiles, in order from lowest to highest inelasticity, to be adjusted by +12%, −13%, −13%, and +4% relative to the nominal simulation by this method.

Event selection
We employ the same hit finding and time clustering as in the ν µ disappearance analysis, and select events whose ν e CC score under the same CVN algorithm exceeds a tuned selection cut.To further purify the sam-ple of ν e CC candidates, we reconstruct events as follows.First, we build three-dimensional event vertices using the intersection of lines constructed from Hough transforms applied to each two-dimensional detector view separately [37,38].Hits in the same view falling roughly along common directions emanating from these vertices are further grouped into "prongs," which are then matched between views based on their extent and energy deposition [39,40].We use these prongs to remove events where the energy of the event is distributed largely transverse to the neutrino beam direction; our simulation and our large sample of cosmic data taken from minimum-bias triggers indicate such events are typically cosmogenic.We further reject events where the prongs fail containment criteria, where extremely long tracks indicate obvious muons, where there are too many hits for proper reconstruction, or where another event in close proximity in both time and space approaches the top of the detector.To combat background events from cosmogenic photon showers entering through the back of the detector, where the overburden is thinner, we also cut events which appear to be pointing toward Fermilab rather than away from it.These events are distinguished by having the number of planes without hits in the portion of the event closest to Fermilab exceeding the number in the portion farthest from Fermilab, the reverse of the expectation for an electromagnetic shower coming from the neutrino beam direction.Events surviving these selections form our "core" sample in both detectors.The predicted composition of the FD sample at various stages in this selection is given in Table III.
We also construct a second, "peripheral" sample of FD events by considering events that have high scores for the CVN ν e hypothesis but which fail the cosmic rejection or containment criteria.These are subjected to a more focused BDT (distinct from the one mentioned in Sec.IV A) trained over the variables used for the containment and cosmic rejection cuts.The containment variables include the closest distance to the top of the detector and the closest distance to any other face of the detector.Variables distinguishing cosmogenic from beam-induced activity include the transverse momentum fraction of the event and the number of hits in the event.Simulation and our cosmic data sample indicate that events in the signal-like regions of both this BDT and CVN are likely to be signal and not the result of externally entering activity and are therefore retained.Distributions for the peripheral sample illustrating the predicted beam and cosmic response in this BDT and the CVN ν e score, as well as comparing the BDT distribution in data and simulation, are given in Fig. 4. Because events on the periphery of the detector are not guaranteed to be fully contained, peripheral events are summed together into a single bin instead of dividing them by the neutrino energy estimate as is done for the core sample.The FD event counts at two stages of the peripheral selection are noted in Table IV.
The ND event sample is predicted to consist of 42% beam ν e , 30% NC background, and 28% ν µ CC background.These predictions include the effect of the datadriven constraints described in Sec.IV B 3. The simulated FD efficiency for the basic quality and containment cuts used in the combined core and peripheral selections relative to all true ν e CC events within a fiducial volume is 92.6%.The remaining core selections, i.e., CVN and cosmic rejection, retain 58.8% of the true ν e CC events in the quality-and-containment population.With the addition of the peripheral sample under the combined CVN+BDT criteria, this figure rises to 67.4%.Improvements to the selection criteria generate an increase of 6.8% in effective exposure [41] relative to our previous results, while the efficiency gain due to the addition of the peripheral sample yields a further increase of 17.4%.

Energy estimation and binning
To estimate the neutrino energy in ν e candidate events, we construct a second-order polynomial in two variables: the sum of the calibrated hit energies from prongs identified as electromagnetic activity and the sum of the ener-gies of hits in the event not within those prongs.The coefficients of this polynomial are fit to minimize the predicted neutrino energy residuals in selected simulated ν e CC events.Whether a prong is considered electromagnetic or not is determined by a deep learning single particle classifier that utilizes both information from the prong itself and the full event [42].This results in an estimator with 11% resolution for both appearance signal and beam background ν e CC events in both detectors.
The expected appearance signal has a narrow peak at the ν µ disappearance maximum, about 1.6 GeV.Additionally, in this analysis, NC and cosmogenic backgrounds concentrate at low reconstructed energies, and beam ν e backgrounds dominate at high energies.Based on these considerations, figure-of-merit calculations based on simulation suggest we limit the neutrino energies we consider to be between 1 and 4 GeV for the FD core sample and 1-4.5 GeV for the peripheral sample.The corresponding core or peripheral range is used for the ND sample when applying the data constraint detailed in Sec.IV B 3. Each of these is further subdivided into three ranges in the CVN classifier output so as to concentrate the sample of highest purity together.The peripheral event sample is treated as a fourth bin.

Near Detector data constraints
The procedure for using the ND data in the ν e analysis is similar to that used for ν µ , extended to account for the particular natures of the signal and beam background components.Appeared electron neutrinos arise from oscillated beam muon neutrinos, so the ν µ -selected candidates in the ND are used to correct the expected ν e appearance signal with the same procedure detailed in Sec.IV A 3. Additionally, the ν µ -selected events are used to verify the ν e selection efficiency.From the ν µ data and simulated samples, we create two subsets where the reconstructed muon track is replaced by a simulated electron shower with the same energy and direction [43].The ν e selection criteria are applied to these electroninserted samples, and the efficiencies for identifying neutrino events in data and simulation, relative to a loose preselection, are found to match within 2%.
As there is no signal and cosmogenic activity is negligible at the ND, the ν e CC candidates at the ND consist entirely of beam background events, originating from CC reactions of the intrinsic ν e component in the beam and mis-identified NC and ν µ CC events.As in our last result [5], we use a combination of data-driven methods to "decompose" the ν e -selected data into these three categories and constrain them independently.We examine low-and high-energy ν µ CC samples at the ND in order to adjust the yields of the parent hadrons that decay into both ν e and ν µ , which constrains the ν e beam background.We also use the observed distributions of time-delayed electrons from stopping µ decay in each analysis bin to constrain the ratio of ν µ CC and NC interactions.The peripheral sample boundary is chosen at the red line: the majority of signal events lie above and to the right and the sample has little cosmogenic contamination there, while events to the left and below are predominantly cosmogenic and are rejected.Right: comparison of the observed distribution (black points) of the BDT variable for peripheral events with the prediction (stacked histogram).sulting decomposition of the selected ν e candidate sample at the ND therefore agrees with the data distribution by construction.The nominal and constrained predictions are shown compared to the data distribution in Fig. 5.The corrections to the beam ν e , NC and ν µ CC components are extrapolated to the FD core sample using the bin-by-bin ratios of the FD and ND reconstructed energy spectra, for each of the three CVN ranges.The predicted beam backgrounds in the FD peripheral sample are corrected according to the results of the extrapolation for the highest CVN bin in the core sample (see Fig. 5).The sum of the final beam-induced background prediction and the extrapolated signal for given oscillation parameters is added to the measured cosmic-induced backgrounds to compare to the observed FD data.

V. SYSTEMATIC UNCERTAINTIES
We evaluate the effect of potential systematic uncertainties on our results by reweighting or generating new simulated event samples for each source of uncertainty and repeating the entire measurement, including the extraction of signal and background yields, the computation of migration matrices, and the calculation of the ratios of FD to ND expectations using each modified simulation sample and applying our constraint procedures.
The effect of each of these uncertainties on the pre- The effect of the decomposition and constraint procedure on the predicted ND candidate νe spectrum; the stacked histogram shows corrected backgrounds (from bottom, beam νe, νµ CC, NC).The three panels show the results for each of the CVN classifier bins, ranging left to right from lower to higher purity.Predictions for each background class prior to correction are given by the dashed lines.The overall corrections to the normalizations of the yields by category are: beam νe CC, +3.0%; NC, +17.0%; and νµ CC, +18.9%.dicted yields of selected ν e CC candidate events is contained in Table V.We estimate the effects on the extracted oscillation parameters sin 2 θ 23 , ∆m 2  32 and δ CP in the joint fit to be as given in Table VI.These are negligibly different from a ν µ -only fit.The largest effects on this analysis stem from uncertainty in our calibrations and energy scales, in the crosssection and final-state interaction (FSI) models in genie, and in the impact of imperfectly simulated event pileup from the neutrino beam on reconstruction and selection efficiencies at the ND.
Calibration and energy scale To evaluate the uncertainty from calibrations and energy scales, which can affect the two detectors differently, we group these uncertainties into absolute (fully positively correlated between TABLE VI.Sources of uncertainty and their estimated average impact on the oscillation parameters in the joint fit.This impact is quantified using the increase in the one-dimensional 68% C.L. interval, relative to the size of the interval when only statistical uncertainty is included in the fit.Simulated data were used and oscillated with the same parameters as in Table V.Given the asymmetry of the sin 2 θ23 interval with respect to its best fit value, only the change in the upper edge is included.The total systematic uncertainty is calculated by adding the individual components in quadrature.detectors) and relative (anticorrelated or uncorrelated) components.Both absolute and relative muon energy scale uncertainties are < 1% based on a combination of thorough accounting of our detectors' material composition and an examination of the parameters in the Bethe formula for stopping power and the energy-loss model of geant4.The overall energy response uncertainty, on the other hand, is driven by uncertainty in our overall calorimetric energy calibration.To investigate the response, we compare simulated and measured data distributions of numerous channels including the energy deposits of muons originating from cosmogenic-and beam-related activity, the energy spectra of electrons arising from the decay of stopped muons, the invariant mass spectrum of neutral pion decays into photons, and the proton energy scales in ND quasielastic-like events.The uncertainty we use is guided by the channel exhibiting the largest differences, the proton energy scale, at 5%.We take this 5% uncertainty as both an absolute energy uncertainty, correlated between the two detectors, and a separate 5% relative uncertainty, since there are not sufficient quasielastic-like events to perform this check at the FD.Cross sections and FSI Estimates for the majority of the cross section and FSI uncertainties that we consider are obtained using the event reweighting framework in genie [15].However, ongoing effort in the neutrino cross section community and the NOvA ND data suggest some modifications are necessary.First, we apply additional uncertainty to the energy-and momentumtransfer-dependence of CC quasielastic (CCQE) scattering due to long-range nuclear correlations [44] according to the prescription in Ref. [21].Second, as the detailed nature of MEC interactions is not well understood, we construct uncertainties for the neutrino energy

Quartile 2
Events / 0.1 GeV FIG. 6.Comparison of the reconstructed energy spectra of selected νµ CC candidates in FD data (black dots) and best-fit prediction (purple).The sample is split into four reconstructed hadronic energy fraction quartiles labeled 1 through 4, where 1 (4) has the best (worst) energy resolution.The majority of the total background (gray, upper) including the cosmogenic subcomponent (blue, lower) lies in the fourth quartile.
dependence, energy-transfer dependence, and final-state nucleon-nucleon pair composition based on a survey of available theoretical treatments [45][46][47].The normalization of the MEC component is recomputed under each of these uncertainties using the same fit procedure used to arrive at the 20% scale factor for the central value prediction.Third, it is now believed that the inflated value of the axial mass in quasielastic scattering (M QE A ) obtained in recent neutrino-nucleus scattering experiments relative to the light liquid bubble chamber measurements is due to nuclear effects that we are now treating explicitly with the foregoing [48].We thus reduce genie's uncertainty for M QE A to ±5% (a conservative estimate of the bubble chamber range [49,50]) from its default of +25% −15% , while retaining genie's central value M QE A = 0.99 GeV/c 2 .Fourth, we increase the uncertainty applied to nonresonant pion production with three or more pions and invariant hadronic mass of W < 3 GeV to 50% to match the default for 1-and 2-pion cases, based on data-simulation disagreements observed in the ND data.Fifth, and finally, we introduce two separate 2% uncertainties on the ratio of ν e CC and ν µ CC cross sections: one to account for potential differences between them due to radiative corrections, and one to consider the possibility of secondclass currents in CCQE events [7,51].
To validate the uncertainties assigned by genie to the NC backgrounds in our analyses, we performed a study within the ν µ CC candidate sample in the ND that measured the rates of neutrons that were produced at the ends of tracks and subsequently recaptured, emitting photons.This study was done by investigating timedelayed activity consistent with a neutron capture, taking into account the tail of the Michel electron time spectrum.The neutron rate is different for the mostly µ − identified in ν µ CC reactions versus the mostly π ± in NC.This study suggested that the NC cross-section uncertainties provided by GENIE, combined together with the calibration uncertainties mentioned previously, account for any differences between data and simulation.Therefore we no longer include the ad hoc 100% additional uncertainty on NC backgrounds used in previous results [5,16].
Normalization We quantify the uncertainty arising from potential imperfections in the simulation of beaminduced pileup in the ND by overlaying a single extra simulated event onto samples of both simulated and data events.We then examine the selection efficiency of this extra event and assign the 3.5% difference between the data and simulation samples as a conservative uncertainty on the normalization of the ND rate.These are added in quadrature with much smaller uncertainties in the detector mass and the total beam exposure to yield an overall normalization systematic.
Other Other contributions to our systematic uncertainty budget are associated with the improved ppfx flux prediction and potential differences between the acceptances of the ND ν µ selection criteria and the FD ν e sample into which the ND corrections are extrapolated in the ν e analysis.Also substantially reduced are the uncertainties in the light response model used for detector simulation.Previous fits of the parameters in the Birks model for scintillator quenching with a secondorder term [52], using proton tracks in candidate ND ν µ CC quasielastic-like events in data, obtained values inconsistent with other measurements of Birks quenching in liquid scintillator [53,54].Previous results therefore used a variation with the other measurements' values to compute an uncertainty.With the addition of Cherenkov light in scintillator to our detector model, however, we find a best fit at the same values preferred by other experiments.To quantify any residual uncertainty in the light model, in this analysis we take alternate predictions where we alter the scintillation and Cherenkov photon yields in the model within the tolerance of agreement with the ND data while holding the muon response fixed (since it is set by our calibration procedure).

VI. RESULTS
We performed a blind analysis in which the FD data were analyzed only after all aspects of the analysis had been specified.An independent implementation of the methods described in Secs.IV-V for incorporating the Near Detector data constraint and assessing the impact of systematic uncertainties, as well as extracting oscillation parameters via likelihood fitting, was used to check the analysis presented in this paper.It produced results consistent with those shown in the following sections.

A. νµ disappearance data
After selection, 126 ν µ CC candidates are observed in the FD.In the absence of oscillations, we would have expected 720.3 +67.4  −47.0 (syst.)ν µ CC candidates based on the extrapolation from the Near Detector, including an expected background of 5.8 misidentified cosmic rays and 3.4 misidentified neutrino events of other types.
Figure 6 shows the observed energy spectrum in each quartile and the corresponding best fit predictions.As noted earlier, most of the predicted background appears in the fourth (worst resolution) quartile.Figure 7 shows Reconstructed Neutrino Energy (GeV) the data of Fig. 6 summed over all of the quartiles.The neutrino energy spectrum exhibits a sharp dip at about 1.6 GeV.Essentially, sin 2 2θ 23 corresponds to the depth of the dip and ∆m 2  32 corresponds to its location.Both of these measurements are sensitive to the energy resolution, so we expect the best measurement in the quartile with best energy resolution.

B. νe appearance data
After selection we observe 66 ν e CC candidate events in the FD including an expected background of 20.3 ± 2.0 (syst.)events.The composition of the expected background is estimated to be 7.3 beam ν e CC events, 6.4 NC events, 1.3 ν µ CC events, 0.4 ν τ CC events, and 4.9 cosmic rays.
Figure 8 shows the distribution of these events as a function of the reconstructed neutrino energy for the three CVN classifier bins and for the peripheral sample, along with the expected background contributions and the best fit predictions.To give some context to the number of observed ν e events, Fig. 9 shows the number of events expected for the best fit values of ∆m 2  32 and sin 2 θ 23 as a function of δ CP , for the two possible mass hierarchies.

C. Joint fit results
We have performed a simultaneous fit to the binned data shown in Figs. 6 and 8. Systematic uncertainties are incorporated into the fit as nuisance parameters with Gaussian penalty terms.Where systematic uncertainties are common between the two data sets, the nuisance parameters associated with the effect are correlated appropriately.In making these fits and in the con-POT-equiv tours and significance levels that follow, we used the following values for physics parameters measured by other experiments [1]: ∆m 2  21 = (7.53± 0.18) × 10 −5 eV 2 /c 4 , sin 2 θ 12 = 0.307 +0.013 −0.012 , sin 2 θ 13 = 0.0210 ± 0.0011.We use a matter density computed for the average depth of the NuMI beam in the Earth's crust for the NOvA baseline of 810 km using the CRUST2.0model [55], ρ = 2.84 g/cm 3 .

Best fits
Table VII gives the parameter values at the best fit point in each relevant mass hierarchy and θ 23 octant combination.The top line shows the overall best fit, which occurs in the normal mass hierarchy and the upper θ 23 octant; the middle line shows best fit in the lower θ 23 octant for the normal mass hierarchy, which is only slightly less significant; and the bottom line shows the best fit in the inverted mass hierarchy, which is disfavored largely because it predicts fewer ν e appearance events than are observed.The column labeled ∆χ 2 represents the difference in χ 2 between the fit and the overall best fit, where χ 2 in this case is −2lnL with L being the likelihood function calculated using Poisson statistics plus Gaussian penalty terms for the systematic uncertainties.There are no best fit values in the inverted mass hierarchy and lower θ 23 octant because the likelihood has no local maximum in this hierarchy-octant region, as will become clear in Fig. 14.The χ 2 for the overall best fit is 84.6 for 72 degrees of freedom.
The precision measurements of sin 2 θ 23 and ∆m 2   32   come from the ν µ disappearance data.A fit to these data alone gives essentially the same values for these parameters in the normal mass hierarchy.However, the best joint ν µ -ν e fit pulls the value of |∆m 2 32 | up by 0.04 × 10 −3 eV 2 /c 4 from the ν µ disappearance only fit in the inverted mass hierarchy.All of the contours and significance levels that follow are constructed following the unified approach of Feldman and Cousins [56], profiling over unspecified physics parameters and systematic uncertainties.
Figure 10 shows the 1, 2, and 3 σ two-dimensional contours for ∆m 2  32 and sin 2 θ 23 , separately for each mass hierarchy.Figure 11 shows a comparison of 90% confidence level contours for these parameters in the normal mass hierarchy for NOvA, T2K [7], MINOS [6], IceCube [57], and Super-Kamiokande [58].All of the experiments have results consistent with maximal mixing.Note that the range 0.4 to 0.6 in sin 2 θ 23 corresponds to the range 0.96 to 1.00 in sin 2 2θ 23 , which is the variable directly measured in ν µ → ν µ oscillations.Figure 12 shows the analogous contours to those of Fig. 10  vs. sin 2 θ23 parameter space consistent with the νe appearance and the νµ disappearance data at various levels of significance.The top panel corresponds to normal mass hierarchy and the bottom panel to inverted hierarchy.The color intensity indicates the confidence level at which particular parameter combinations are allowed.Figures 13,14, and 15 show the significance with which values of |∆m 2 32 |, sin 2 θ 23 , and δ CP are disfavored in the two mass hierarchies, respectively.The results in Fig. 14 differ from the ones previously reported [16] in that the disfavoring of maximal mixing (θ 23 = π/4) has changed from 2.6 standard deviations (σ) to 0.8 σ in the present results.This change was caused by three changes, each of which moved θ 23 closer to maximal mixing.The largest effect was due to new simulations and calibrations.The two smaller effects were from new selection and analysis procedures and from the additional 2.80 × 10 20 POT of data included here.The additional data taken by itself favored maximal disappearance.In Fig. 15 two curves are shown in the normal mass hierarchy, one for each of the θ 23 octants, corresponding to the near degeneracy shown in Fig. 14.Only one curve is shown for the inverted mass hierarchy since there is only one minimum, which occurs in the upper octant.The point of minimum significance in the inverted mass hierarchy differs among the figures because, although the ∆χ 2 's are identical (see Table VII), the translation of ∆χ 2 to significance depends on which oscillation parameters are profiled.
Table VIII shows the 1 σ confidence intervals for ∆m FIG.11.Comparison of measured 90% confidence level contours for ∆m 2  32 vs. sin 2 θ23 for this result (black line; best-fit value, black point), T2K [7] (green dashed), MINOS [6] (red dashed), IceCube [57] (blue dotted), and Super-Kamiokande [58] (purple dash-dotted).intervals in the inverted mass hierarchy.Finally, we have calculated the significance level for the rejection of the inverted hierarchy using the same procedure as in the above contours and confidence intervals, namely by profiling over all the other physics parameters and the systematic uncertainties.Frequentist coverage was checked following the suggestion of Berger and Boos [59].The entire inverted mass hierarchy region is disfa- Significance at which each value of δCP is disfavored in the normal (blue, lower) or inverted (red, upper) mass hierarchy.The normal mass hierarchy is divided into upper (solid) and lower (dashed) θ23 octants corresponding to the near degeneracy in sin 2 θ23.vored at the 95% confidence level.

FIG. 1 .
FIG.1.Predicted composition of the NuMI beam at the ND with the horns focusing positively charged hadrons.Curves from top to bottom: νµ, νµ, νe, νe.TableIgives the fractional composition for each neutrino flavor integrated from 1-5 GeV.
FIG. 2.Comparison of the reconstructed neutrino energy for selected νµ CC events (black dots) in the ND with area-normalized simulation (red line).Shading represents the bin-to-bin systematic uncertainties.The gray area, which is nearly indistinguishable from the lower figure boundary, shows the simulated background.

10 FIG. 3 .
FIG.3.Comparison of selected νµ CC candidates (black dots) in the ND data to the prediction (red histograms) in the hadronic energy fraction quartiles, where the prediction is absolutely normalized to the data by exposure.The expected background contributions (gray) are smaller in the quartiles with better resolution.The shaded band represents the quadrature sum of all systematic uncertainties.These distributions are the input to the extrapolation procedure described in the text.
FIG. 5.The effect of the decomposition and constraint procedure on the predicted ND candidate νe spectrum; the stacked histogram shows corrected backgrounds (from bottom, beam νe, νµ CC, NC).The three panels show the results for each of the CVN classifier bins, ranging left to right from lower to higher purity.Predictions for each background class prior to correction are given by the dashed lines.The overall corrections to the normalizations of the yields by category are: beam νe CC, +3.0%; NC, +17.0%; and νµ CC, +18.9%.
Figures 13,14, and 15 show the significance with which values of |∆m2  32 |, sin 2 θ 23 , and δ CP are disfavored in the two mass hierarchies, respectively.The results in Fig.14differ from the ones previously reported[16] in that the disfavoring of maximal mixing (θ 23 = π/4) has changed from 2.6 standard deviations (σ) to 0.8 σ in the present results.This change was caused by three changes, each of which moved θ 23 closer to maximal mixing.The largest effect was due to new simulations and calibrations.The two smaller effects were from new selection and analysis procedures and from the additional 2.80 × 10 20 POT of data included here.The additional data taken by itself favored maximal disappearance.In Fig.15two curves are shown in the normal mass hierarchy, one for each of the θ 23 octants, corresponding to the near degeneracy shown in Fig.14.Only one curve is shown for the inverted mass hierarchy since there is only one minimum, which occurs in the upper octant.The point of minimum significance in the inverted mass hierarchy differs among the figures because, although the ∆χ 2 's are identical (see TableVII), the translation of ∆χ 2 to significance depends on which oscillation parameters are profiled.TableVIIIshows the 1 σ confidence intervals for ∆m 2 32 ,

FIG. 12 .
FIG.12.Regions of sin 2 θ23 vs. δCP parameter space consistent with the νe appearance and the νµ disappearance data.The top panel corresponds to normal mass hierarchy (∆m 2 32 > 0) and the bottom panel to inverted hierarchy (∆m2  32 < 0).The color intensity indicates the confidence level at which particular parameter combinations are allowed.

sin 2 θ
23 , and δ CP in the normal mass hierarchy, corresponding to Figs.13-15.There are no 1 σ confidence

FIG. 14 .
FIG.14.Significance at which each value of sin 2 θ23 is disfavored in the normal (blue, lower) or inverted (red, upper) mass hierarchy.The vertical dotted line indicates the point of maximal mixing.
FIG. 15.Significance at which each value of δCP is disfavored in the normal (blue, lower) or inverted (red, upper) mass hierarchy.The normal mass hierarchy is divided into upper (solid) and lower (dashed) θ23 octants corresponding to the near degeneracy in sin 2 θ23.
Table I gives the fractional composition for each neutrino flavor integrated from 1-5 GeV.

TABLE I .
Predicted beam flux composition in the 1 to 5 GeV neutrino energy region in the absence of oscillations.

TABLE II .
Predicted composition of the νµ CC candidate sample in the FD, in event counts, at various stages in the selection process.Oscillation parameters used in the prediction are the best fit values from Sec. VI.
The peripheral sample is a signal-rich subset of νe FD candidates that failed the core cosmic rejection or containment criteria (see text).Left: The two-dimensional BDT-CVN space used in the definition of the peripheral sample.The predicted distribution of νe appearance signal events (boxes) is shown superimposed on the predicted purity in each bin (shaded color).

TABLE III
. Predicted composition of the core νe CC candidate sample at the FD, in event counts, at various stages in the selection process.Oscillation parameters used in the prediction are the best fit values from Sec. VI.These figures do not include the effect of the extrapolation procedure described in Sec.IV B 3. TABLE IV.Predicted composition of the peripheral νe CC candidate sample, in event counts, at two stages in the selection process.Here "basic quality" refers to events that pass beam and detector data quality cuts but fail the core sample containment criteria.Parameters are as in TableIII.
Comparison of the neutrino energy spectra of selected νe CC candidates in the FD data (black dots) with the best fit prediction (purple lines) in the three CVN classifier bins and the peripheral sample.The total expected background (gray, upper) and the cosmic component of it (blue, lower) are shown as shaded areas.The events in the peripheral bin have energies between 1 and 4.5 GeV.Total number of νe CC candidate events observed in the FD (gray) compared to the prediction (color) as a function of δCP.The color lines correspond to the best fit values of sin 2 θ23 and ∆m2  32from TableVII, with the upper two curves (blue) representing two octants in the normal mass hierarchy (∆m232 > 0) and the lower curve (red) the inverted hierarchy (∆m232 < 0).The color bands correspond to 0.43 ≤ sin 2 θ23 ≤ 0.60.All other parameters are held fixed at the best-fit values.

TABLE VII .
Best fit values.See text for further explanation.
Significance at which each value of |∆m 2 32 | is disfavored in the normal (blue, lower) or inverted (red, upper) mass hierarchy.
TABLE VIII. 1 σ confidence intervals for physics parameters in the normal mass hierarchy.