An Improved Measurement of Neutrino Oscillation Parameters by the NOvA Experiment

We present new $\nu_\mu\rightarrow\nu_e$, $\nu_\mu\rightarrow\nu_\mu$, $\overline{\nu}_\mu\rightarrow\overline{\nu}_e$, and $\overline{\nu}_\mu\rightarrow\overline{\nu}_\mu$ oscillation measurements by the NOvA experiment, with a 50% increase in neutrino-mode beam exposure over the previously reported results. The additional data, combined with previously published neutrino and antineutrino data, are all analyzed using improved techniques and simulations. A joint fit to the $\nu_e$, $\nu_\mu$, $\overline{\nu}_e$, and $\overline{\nu}_\mu$ candidate samples within the 3-flavor neutrino oscillation framework continues to yield a best-fit point in the normal mass ordering and the upper octant of the $\theta_{23}$ mixing angle, with $\Delta m^{2}_{32} = (2.41\pm0.07)\times 10^{-3}$ eV$^2$ and $\sin^2\theta_{23} = 0.57^{+0.03}_{-0.04}$. The data disfavor combinations of oscillation parameters that give rise to a large asymmetry in the rates of $\nu_e$ and $\overline{\nu}_e$ appearance. This includes values of the CP-violating phase in the vicinity of $\delta_\text{CP} = \pi/2$ which are excluded by $>3\sigma$ for the inverted mass ordering, and values around $\delta_\text{CP} = 3\pi/2$ in the normal ordering which are disfavored at 2$\sigma$ confidence.

ν e, and ν µ → ν µ oscillation measurements by the NOvA experiment, with a 50% increase in neutrino-mode beam exposure over the previously reported results. The additional data, combined with previously published neutrino and antineutrino data, are all analyzed using improved techniques and simulations. A joint fit to the νe, νµ, ν e, and ν µ candidate samples within the 3-flavor neutrino oscillation framework continues to yield a best-fit point in the normal mass ordering and the upper octant of the θ23 mixing angle, with ∆m 2 32 = (2.41 ± 0.07) × 10 −3 eV 2 and sin 2 θ23 = 0.57 +0.03 −0.04 . The data disfavor combinations of oscillation parameters that give rise to a large asymmetry in the rates of νe and ν e appearance. This includes values of the CP-violating phase in the vicinity of δCP = π/2 which are excluded by > 3 σ for the inverted mass ordering, and values around δCP = 3π/2 in the normal ordering which are disfavored at 2 σ confidence.

I. INTRODUCTION
We report new measurements of neutrino oscillation parameters using neutrino and antineutrino data from the NOvA experiment. The data include a 50% increase in neutrino-mode beam exposure over the previously reported results [1]. We perform a joint fit to ν µ (ν µ ) → ν e (ν e ) and ν µ (ν µ ) → ν µ (ν µ ) oscillations utilizing improvements in the analysis of these data.
Numerous experiments [2][3][4][5][6][7][8][9][10] corroborate the paradigm in which three neutrino mass eigenstates (ν 1 , ν 2 , ν 3 ) mix to form the three flavor eigenstates (ν e , ν µ , ν τ ). The mixing can be expressed by the unitary matrix, U PMNS , named for Pontecorvo, Maki, Nakagawa, and Sakata. U PMNS can be parameterized by three mixing angles (θ 12 , θ 23 , θ 13 ) along with a phase (δ CP ) that, if different from 0 or π, indicates violation of Charge-Parity (CP) symmetry. Neutrino mixing gives rise to oscillations from one flavor state to another, dependent on the mixing parameters and the mass splittings (∆m 2 ij ≡ m 2 i − m 2 j ). Using the definition of ν 1 as having the largest ν e contribution, it has been established that ∆m 2 21 is positive, and therefore, the ν 2 mass eigenstate is heavier than ν 1 . However, the sign of the larger mass splitting, ∆m 2 32 , is unknown. If this term is positive, then the third mass eigenstate is the heaviest, and the mass ordering is la-beled as the Normal Ordering (NO) (also referred to as Normal Hierarchy). The alternative is referred to as Inverted Ordering (IO) (or Inverted Hierarchy). Knowing the mass ordering would constrain models of neutrino masses [11][12][13][14][15] and could aid in the resolution of the Dirac or Majorana nature of the neutrino [16,17].
The mass ordering affects the rates of ν µ → ν e and ν µ → ν e oscillations when neutrinos travel through the Earth as compared to a vacuum. Coherent forward scattering on electrons in the Earth's crust enhances the rate of ν µ → ν e oscillations and suppresses ν µ → ν e for the NO while the enhancement and suppression is reversed for the IO. This matter effect [18] changes the oscillation probabilities for NOvA by ∼ 20%. Depending on the value of δ CP and the mass ordering itself, NOvA may be able to exploit the resulting neutrino-antineutrino asymmetry to measure the sign of ∆m 2 32 and thus determine the mass ordering.
NOvA also has sensitivity to δ CP , which will increase the ν µ → ν e oscillation probability if sin δ CP is positive and suppress oscillations if negative (the effect is reversed for antineutrinos). Additionally, a non-zero value of sin δ CP would identify the neutrino sector as a source of CP violation which is central to some explanations of the matter-antimatter asymmetry observed based on leptogenesis [19][20][21][22][23]. Since a measurement of both the mass ordering and δ CP rely on a comparison of ν e and ν e appearance, certain combinations of δ CP and mass ordering will be degenerate with others for NOvA's oscillation baseline.
Here, we reanalyze the data taken in the antineutrinomode beam from June 29, 2016, to February 26, 2019, with an exposure of 12.5 × 10 20 protons on target (POT) delivered during 321.1 s of integrated beam-pulse time. These data are combined with an increased, and reanalyzed, neutrino-mode beam exposure of 13.6 × 10 20 POT from 555.3 s of integrated beam-pulse time recorded between February 6, 2014, to March 20, 2020. During these periods, the proton source achieved an average power of 650 kW, and a peak hourly-averaged power of 756 kW.
In addition to the increased neutrino-mode beam exposure, this analysis introduces various improvements that will be described in detail in the following sections. There are changes to the underlying neutrino interaction simulation, particle propagation, and detector response models. The reconstruction uses a new clustering algorithm and expands the use of neural networks. Furthermore, the Near-to-Far extrapolation method has been expanded to further constrain the FD predictions, which also reduces the impact of systematic uncertainties on the analysis by up to 9% as compared to the previous method. Finally, we have improved some systematic uncertainties and introduced new ones associated with the above changes.

II. THE NOvA EXPERIMENT AND SIMULATIONS
NOvA observes ν µ (ν µ ) → ν e (ν e ) appearance and ν µ (ν µ ) → ν µ (ν µ ) disappearance oscillations using two functionally-identical tracking calorimeters [25] deployed in Fermilab's NuMI beam [26]. Charged particle tracking is accomplished via PVC cells filled with a mineral oilbased liquid scintillator [27]. The cells are 6.6 cm × 3.9 cm in cross section and are oriented in alternating vertical and horizontal planes to achieve 3D reconstruction. The 290 ton Near Detector (ND) is located 100 m underground and ∼1 km from the production target. The main body of the ND is followed by a muon range stack where the active planes are interleaved with steel plates. The 14 kton Far Detector (FD) is located at Ash River, Minnesota, ∼810 km from the source. Being located on the surface with a modest rock overburden, the FD receives a cosmic-ray flux of 130 kHz. This analysis benefits from an updated simulation of the geometries of the detectors and their surroundings that more accurately reflects the surrounding rock composition and detectors as built.
Neutrino interactions are simulated using a custom model configuration of genie 3.0.6 [49,50] tuned to external and NOvA ND data. 1 In this configuration, charged-current (CC) quasi-elastic (QE) scattering is simulated using the model of Nieves et al. [53], which includes the effects of long-range nucleon correlations calculated according to the Random Phase Approximation (RPA) [53][54][55]. The CCQE axial vector form factor is a zexpansion parameterization tuned to neutrino-deuterium scattering data [56]. CC interactions with two nucleons producing two holes (2p2h) are given by the IFIC València model [57,58]. The initial nuclear state is represented by a local Fermi gas in both the QE and 2p2h models, and by a global relativistic Fermi gas for all other processes. Baryon resonance (RES) and coherent pion production are simulated using the Berger-Sehgal models with final-state mass effects taken into account [59,60]. Deep inelastic scattering (DIS) and non-resonant background below the DIS region are described using the Bodek-Yang model [61] with hadronization simulated by a data-driven parameterization [62] coupled to pythia [63]. Bare nucleon cross-sections for RES, DIS, and nonresonant background processes are tuned by genie to neutrino scattering data. Final state interactions (FSI) are simulated by the genie hN semi-classical intranuclear cascade model in which pion interaction probabilities are assigned according to Oset et al. [64] and pion-nucleon scattering data.
The 2p2h and FSI models in this genie configuration are adjusted to produce a NOvA-specific neutrino interaction model tune. The 2p2h model is fit to ν µ CC inclusive scattering data from the NOvA ND. Inspired by 1 Neutrino interactions in this analysis were inadvertently simulated with event kinematics of GENIE configuration N18 10j 00 000 but integrated rates with configuration N18 10j 02 11a. These two configurations have the same model set and differ only in the tune of the resonant, nonresonant background, and DIS free nucleon cross sections, where the N18 10j 00 000 tune used inclusive neutrino scattering data and the N18 10j 02 11a tune used 1π and 2π production in addition to the inclusive neutrino scattering data [51,52]. The predicted Far Detector event spectra generated using N18 10j 02 11a are consistent with the predictions used in this measurement within the systematic uncertainties.
Gran et al. [65], this 2p2h tune enhances the base model as a function of energy and momentum transfer to the nucleus and is applied to all CC 2p2h interactions for both the neutrino and anti-neutrino beams. The parameters governing π ± and π 0 FSI are adjusted to obtain agreement with π + on 12 C scattering data [66][67][68][69][70][71][72].
The propagation of final state particles through the detectors is simulated by an updated version of geant4 (v10.4) [73], which provides the input for the detector response simulation [74]. In addition, a custom patch to the new version implements an exact calculation of the density effect correction to the Bethe equation using Sternheimer's method [75] as opposed to the approximate parameterization used previously (a 1% or less change to the muon range and energy lost in dead material).
The absolute energy scale for both detectors is calibrated using the minimum ionizing portion of stopping cosmic-ray muon tracks [76]. The calibration procedure is now applied separately to the data in shorter time periods to account for an observed 0.3% decrease in detected light per year.

III. RECONSTRUCTION AND SELECTION
The first stage of reconstruction is to group hits, which are measurements of deposited energy in a cell above a preset threshold, into single-neutrino-interaction events. This clustering, performed based on hit proximity in time and space, now uses a new method that reduces the rate of mis-clustered hits in the high occupancy environment of the ND [77]. Mis-clustering had previously led to differences in data-MC selection efficiency, which are now reduced to the sub-percent level. The other reconstruction techniques remain unchanged from the previous analysis [1].
For each event, initial selections are applied to ensure basic data quality. Additionally, events are required to be sufficiently far from the edges of the detector such that energy is not lost to exiting final-state particles, and so entering background events are not selected as signal. These containment criteria have been re-optimized for this analysis due to changes in the geometry model and hit grouping algorithm, but follow the same outline as described in Ref. [1].
A convolutional neural network, CNN evt [78], is used to classify neutrino event candidates into ν e CC, ν µ CC, NC, or cosmogenic background. The network is trained using simulated calibrated hits that have been clustered into single neutrino interactions, as well as cosmogenic data. Scores from CNN evt are used to create two non-overlapping samples of either inclusive ν µ (ν µ ) CC or ν e (ν e ) CC candidate events. Updates to this algorithm provide improved performance and decreased dependency on calorimetric energy, the dominant source of systematic uncertainty in the results presented here. This is achieved by scaling up or down the energy of all hits while training the CNN. The scale factors used are drawn on an event-by-event basis from a normal distribution with a 1 σ range from 0.9 -1.1 [79]. This training procedure reduced the influence of calibration uncertainties on CNN evt classification decisions to a negligible level. Effective rejection of cosmogenic backgrounds at the FD is paramount due to the significant flux of cosmicray particles it receives. A new CNN, trained to identify cosmogenic backgrounds has been introduced, is applied in parallel to cosmic-identifying boosted decision trees (BDTs). The BDTs have been trained on samples selected to contain signal-like cosmogenic particles. Together the CNN and BDTs reduce the cosmic contamination in the selected samples to < 5%, a total reduction of 6 orders of magnitude, comparable to the previous analysis. For fully contained ν e events, the BDT replaces the previous cosmic rejection method, which directly used reconstructed position and kinematic event information.
Neutrino energy, E ν , is determined using different methods for the ν e and ν µ CC candidate events. The energies of ν e CC candidates are parameterized using a quadratic function determined from a 2D fit to the simulated electromagnetic (EM) and hadronic calorimetric energies (E EM and E had respectively). The two components produce different detector responses and are separated using a third CNN classifier that identifies EM-like hit clusters within the event with the remaining clusters being classified as hadronic [80]. For ν µ CC candidates, E ν is the sum of the muon energy, determined by the track length, and the total calorimetric energy of the hadronic system, E had . The muon is identified with a BDT that utilizes track length, multiple Coulomb scattering, and energy deposition, while the hadronic system is taken as all hits not associated with the muon track.
The selection criteria and energy estimation techniques were developed based on ND beam and FD cosmic data, along with simulated samples prior to inspecting the FD beam data distributions. The algorithms were trained separately on neutrino and antineutrino beam modes due to differences in beam purity and interactions.
The sensitivity of the oscillation fit is enhanced by splitting the fully contained ν e and ν e CC, "core", samples into low and high purity bins, based on the scores output by CNN evt . At the FD, the ν e (ν e ) selection efficiency for signal events in the core sample is 54% (64%) 2 . To further increase the efficiency of the FD sample, a "peripheral" selection is included, consisting of events that fail the containment or cosmic rejection requirements but pass more strict selection criteria on the cosmic BDT and CNN evt . This sample increases the total ν e (ν e ) selection efficiency to 63% (75%) 2 but is included only as an integrated rate in the oscillation fits due to possible energy bias caused by particles leaving the detector. Properties of these subsamples are summarized in Table I.
For ν µ CC candidates, the position and amplitude of the oscillation maximum in the FD energy spectra are strongly dependent on ∆m 2 32 and θ 23 , respectively. To maximize the sensitivity to these parameters, the candidates are divided into four equally populated samples based on the hadronic energy fraction, E frac = E had /E ν , which is correlated with energy resolution and background contamination as summarized in Table I. Sensitivity is further increased by using variably-sized E ν bins for these samples.

IV. NEAR-TO-FAR EXTRAPOLATION
This analysis extracts oscillation parameters using data-driven predictions of the FD spectra largely derived from high statistics measurements in the ND. The ν µ (ν µ ) disappearance and ν e (ν e ) appearance signal spectra in the FD are predicted using the spectra of ν µ (ν µ ) CC candidate events in the ND (Fig. 1a). The procedure begins with reweighting the simulation to obtain agreement with the data in each reconstructed E ν bin of the ND ν µ (ν µ ) CC candidate samples. Predicted rates of NC, ν e CC, and ν e CC interactions in the samples (<0.5% total) are taken directly from the simulation and subtracted. The wrong-sign component of the samples (2.9% and 10.5% in the neutrino and antineutrino beams respectively) is also taken directly from the simulation. The resulting corrected ν µ + ν µ CC reconstructed E ν spectra are transformed to true E ν using the simulation. The spectra are then multiplied by the appropriate farto-near ratios of the simulated samples in bins of true E ν . This step accounts for beam divergence, differences in selection efficiency and acceptance between the two detectors, and the differences in the ν µ and ν e cross sections. Oscillation probabilities are applied to yield the predicted disappearance or appearance signal spectra in true E ν at the FD. Matter effects are included in the oscillation probability calculations, with the Earth's crust density assumed to be uniformly 2.84 g/cm 3 [81]. Finally, the predicted spectra are converted back to reconstructed E ν .
To reduce potential bias and the impact of uncertainties from the neutrino interaction model, the extrapolation to predict the disappearance and appearance signals is performed using variables in addition to E ν . As in the previous analysis, the extrapolations for the disappearance samples are done separately in each reconstructed hadronic energy fraction range (as given in Table I), enabling neutrino interaction processes that occur in different inelasticity regions to be constrained independently. In this analysis, the extrapolations for both disappearance and appearance samples are additionally performed separately in bins of reconstructed transverse momentum, p T , of the final state charged lepton. The smaller transverse extent of the ND leads to lower acceptance at higher p T in the ND than in the FD (Fig. 2), which results in the extrapolated predictions being sensitive to the modeling of the p T -dependence of the neutrino interactions. Extrapolating in bins of p T reduces this sensitivity by enabling the ND data to constrain the p Tdependence. In the ND samples, the p T bins divide each E ν bin into three equal populations for the extrapolation, and the resulting FD predictions are summed over the p T bins for the oscillation fit.
Background spectra at the FD are also predicted using data-driven techniques. Cosmogenic backgrounds in both the appearance and disappearance samples are estimated using FD data collected outside the NuMI beam time window. Beam-induced backgrounds in the appearance samples are primarily CC interactions from the irreducible ν e + ν e component of the beam, with contributions from mis-identified NC and ν µ + ν µ CC interactions. The FD spectra for these backgrounds are predicted using the spectra of ν e (ν e ) CC candidate events in the ND (Fig. 1b). Since the relative event rate between the ND and FD is different for the background components, the relative contribution of the different background components in the data needs to be estimated. In neutrino beam-mode these estimates are datadriven [77,83] while they are taken directly form the simulation in antineutrino beam-mode.

V. SYSTEMATIC UNCERTAINTIES
The impacts of systematic uncertainties are evaluated by varying the simulation via event reweighting or simulating alternative event samples and repeating the extrapolation procedure. Uncertainties associated with the neutrino flux, neutron modeling, and detector calibrations are unchanged from the previous analysis [1].
Detector calibration uncertainties remain dominant and are driven by a 5% uncertainty in the calorimetric energy scale. Additionally, a new time-dependent calibration uncertainty is included to account for any residual differences remaining after performing the calibration over shorter time periods as mentioned previously.
Neutrino interaction model uncertainties are evaluated using the event reweighting framework in genie with additional uncertainties constructed by NOvA as follows. Uncertainties on CCQE RPA, low-Q 2 RES suppression, 2p2h, and non-resonant and incoherent Nπ production are established for the new model set using methods similar to those in Ref. [84]. Pion FSI uncertainties are based on comparisons to π + on 12 C scattering data [66][67][68][69][70][71][72] and prior studies using an alternative neutrino interaction generator [85]. Uncertainties on the ν e (ν e ) CC cross section relative to the ν µ (ν µ ) CC cross section due to radiative corrections and possible second-class currents are unchanged from previous analyses [83].
As in the previous analysis, uncertainties are included that are detector specific or account for differences between the ND and FD: the detector masses, beam exposures, kinematic acceptances, beam-induced pile-up, ν e CC selection in the ND, and cosmogenic backgrounds in the FD. The improved hit clustering algorithm reduces pile-up effects in the ND, decreasing uncertainties for the associated data-MC selection efficiency differences. An uncertainty for kinematic acceptance differences between the detectors was overestimated in the previous analysis and is subdominant in this analysis after correction. Extrapolating in p T bins would have substantially reduced the effect of this uncertainty even if left uncorrected.
Uncertainties arising from the custom light model are assigned based on comparison to a more robust response model that was not fully incorporated into the simulation for this analysis. This model is constrained by a sample of ND proton candidates in addition to the muon sample used previously. Differences in the detector response between the proton and muon samples also provide a data-driven uncertainty on the relative production of Cherenkov and scintillation light in the model.
Quantities affected by lepton reconstruction uncertainties include the muon energy scale and lepton angle. The muon energy scale uncertainty now includes a detector mass uncertainty with a component that is uncorrelated between the detectors, plus a correlated component accounting for the Fermi density effect and muon range differences across models. Extrapolating in p T bins introduces a dependence on the reconstructed lepton angle for which a 10 mrad uncorrelated uncertainty is applied. Figure 3 shows the impact of the systematic uncertainties on the measurement of sin 2 θ 23 , ∆m 2 32 , and δ CP as evaluated at the determined best-fit point. The extrapolation method significantly reduces the impact of the detector correlated beam flux and neutrino interaction model uncertainties. In contrast, energy calibration and uncorrelated uncertainties that reflect ND-FD differences are less constrained by extrapolation. Figure 3 also shows the impact of uncertainties for extrapolation with and without p T bins. Extrapolating in p T bins reduces the interaction model uncertainty by 10-30%, and the total systematic uncertainty by up to 9%. Detector calibration, detector response, and neutron modeling un- certainties that affect the reconstructed energy of the recoiling hadronic system, which is correlated with p T , are more modestly reduced. The extrapolation in bins of p T depends on reconstructed lepton kinematics and results in a marginal increase in the associated uncertainties.

VI. RESULTS
The extrapolated predictions of the FD spectra are recomputed for varying oscillation parameters and compared to data using a Poisson negative log-likelihood ratio, −2 ln L. The best-fit parameters minimize −2 ln L.
The following solar and reactor neutrino experiment constraints are used: ∆m 2 21 = 7.53 × 10 −5 eV 2 , sin 2 θ 12 = 0.307, and sin 2 θ 13 = 0.0210 ± 0.0011 [86]. The parameters ∆m 2 32 , sin 2 θ 23 , and δ CP are varied without constraints while the 64 systematic uncertainties are assigned penalty terms equal to the square of the number of standard deviations by which they vary from their nominal values. The value of sin 2 θ 13 is allowed to float similarly. Feldman-Cousins' unified approach [87,88] is used to determine the confidence intervals for the oscillation parameters. All significances given, or plotted, are FC-corrected values. The fitted parameters not shown are profiled over. Figure 4 shows the energy spectra of the ν µ CC, ν µ CC, ν e CC, and ν e CC candidates recorded at the FD. The distributions are compared to the oscillation best-fit expectations. Table II summarizes the total event counts and estimated compositions of the selected samples. The CC candidate event samples recorded at the FD include 211 (105) observed ν µ (ν µ ) → ν µ (ν µ ) events and 82 (33) ν µ (ν µ ) → ν e (ν e ) candidate events. The latter ν e (ν e ) appearance sample has an estimated background of 26.8 +1.6 −1.7 (14.0 +0.9 −1.0 ). This analysis determines a best-fit in the normal mass ordering and upper θ 23 octant (significance of 1.0 σ and 1.2 σ, respectively), where −2 ln L = 173.55 for 175 degrees of freedom (p-value of 0.705). The data disfavor combinations that lead to a strong asymmetry in the rate of ν e versus ν e appearance; therefore, the inverted mass ordering with δ CP = π/2 is excluded at more than 3 σ and the normal mass ordering with δ CP = 3π/2 is disfavored at 2 σ confidence. However, owing to the degeneracies, the 90% confidence level allowed regions cover all values of δ CP given permutations of mass ordering and octant. Thus, the current data do not exhibit a preference concerning CP conservation versus violation. Table III shows the best-fit parameter values for each choice of θ 23 octant and mass ordering. Figure 5 compares the 90% confidence level contours for ∆m 2 32 and sin 2 θ 23 with those of other experiments [91][92][93][94] 3 . Allowed regions in sin 2 θ 23 and δ CP are shown in Fig. 6 and are compared with a recent best fit from T2K [91] 3 .   TABLE II. Event counts at the FD, both observed and predicted at the best-fit point (see Table III  As shown in Fig. 6a, the T2K best-fit point is in the NO but lies in a region that NOvA disfavors. However, some regions of overlap remain. Figure 6b shows that for IO, the T2K allowed region at 90% confidence level is entirely contained within the corresponding NOvA allowed region. This outcome reflects in part the circumstance that T2K observes a relatively more pronounced asymmetry in ν e versus ν e oscillations. Although each experiment reports a mild preference for NO, it has been suggested that a joint fit of the two experiments might converge on an IO solution [97]. Some authors have also explored the possibility that the differences in the ν µ → ν e and ν µ → ν e rates seen by the experiments are explained by additional non-standard matter effects [98,99].
In conclusion, we have presented improved measurements of oscillation parameters ∆m 2 32 , sin 2 θ 23 , and δ CP , including an expanded data set and enhanced analysis techniques with respect to previous publications. These measurements continue to favor the normal mass ordering and upper octant of sin 2 θ 23 , as well as values of the oscillation parameters that do not lead to a large asymmetry in ν µ → ν e and ν µ → ν e oscillation rates.

VII. 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, RFBR, RMES, RSF, and BASIS Foundation, Russia; CNPq and FAPEG, Brazil; STFC, UKRI, and the Ref. [91], the dataset remains unchanged and the same approach is used. The conclusions drawn from the comparisons of the contours remains unchanged. 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.