Measurement of neutrino-induced neutral-current coherent π 0 production in the NOvA near detector

The cross section of neutrino-induced neutral-current coherent π 0 production on a carbon-dominated target is measured in the NOvA near detector. This measurement uses a narrow-band neutrino beam with an average neutrino energy of 2.7 GeV, which is of interest to ongoing and future long-baseline neutrino oscillation experiments. The measured, flux-averaged cross section is σ ¼ 13 . 8 (cid:2) 0 . 9 ð stat Þ(cid:2) 2 . 3 ð syst Þ × 10 − 40 cm 2 = nucleus, consistent with model prediction. This result is the most precise measurement of neutral-current coherent π 0 production in the few-GeV neutrino energy region.

The cross section of neutrino-induced neutral-current coherent π 0 production on a carbon-dominated target is measured in the NOvA near detector. This measurement uses a narrow-band neutrino beam with an average neutrino energy of 2.7 GeV, which is of interest to ongoing and future long-baseline neutrino oscillation experiments. The measured, flux-averaged cross section is σ ¼ 13.8 AE 0.9ðstatÞAE 2.3ðsystÞ × 10 −40 cm 2 =nucleus, consistent with model prediction. This result is the most precise measurement of neutral-current coherent π 0 production in the few-GeV neutrino energy region. DOI: 10.1103/PhysRevD.102.012004 Neutrinos can interact coherently with target nuclei and produce outgoing pions via either neutral-current (NC) or charged-current (CC) interactions. In the case of an NC interaction, a π 0 is produced: Coherent interactions are characterized by very small momentum transfer to the target nucleus with no exchange of quantum numbers, while the target nucleus remains in its ground state. The characteristic signal topology of NC coherent π 0 production is a single, forward-going π 0 , with no other hadrons in the final state. There are two major motivations for measuring the NC coherent π 0 cross section. First, coherent π 0 production is a contribution to the background of long-baseline ν μ → ν e oscillation measurements. In some neutrino detectors, the photons from π 0 decay are reconstructed as electromagnetic showers, which are often difficult to separate from the showers induced by electrons. An NC π 0 event can be misidentified as a ν e -CC signal event if the two photon showers are not spatially separated or if one shower is undetected. Knowledge of coherent π 0 production provides Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP 3 . a constraint on the size of this background. Second, coherent pion production provides insight into the weak hadronic current structure and serves as a test of the partially conserved axial current (PCAC) hypothesis [1]. Models based upon PCAC relate the coherent pion production to the pion-nucleus elastic scattering cross section at the Q 2 ¼ 0 limit and extrapolate to low but nonzero Q 2 values. Such models include the Rein-Sehgal model [2,3] in the GENIE neutrino generator [4] used for this analysis. Further improvement of PCAC models was made by Berger et al. (Berger-Sehgal model) [5] and others [6][7][8][9][10]. The PCAC models are known to perform well in their intended multi-GeV energy ranges. Their performance in NC interactions at neutrino energies of a few GeV, however, remains to be established. There is another class of models, referred to as microscopic models, that do not rely on PCAC. These models are built using pion production amplitudes at the nucleon level [11][12][13][14][15][16][17] and are expected to be more reliable than PCAC-based models at neutrino energies below 1 GeV, where the Δ resonance dominates the weak pion production. They typically miss contributions from higher resonances above Δ at higher neutrino energy.
The NC coherent π 0 cross section contributes roughly 1% of the inclusive neutrino interactions in the few-GeV neutrino energy region, much smaller than other π 0 production modes. This challenging situation requires the extraction of a small signal with large backgrounds. The backgrounds arise mainly from NC-induced baryon resonance (RES) interactions and π 0 production from NC deep-inelastic scattering (DIS) interactions, where only a single π 0 is reconstructed. Diffractive (DFR) π 0 production, where neutrinos scatter off free protons (hydrogen) with small momentum transfer and produce π 0 's, also contributes to the background. Unlike coherent interactions, a recoil proton is sometimes visible in DFR π 0 production, which makes the DFR event topology potentially different from the coherent signal. The visibility of the recoil proton in a detector depends upon the proton's kinetic energy (T p ), which is determined by T p ¼ jtj=2m p , where jtj is the square of four-momentum transfer to the nucleus and m p is the proton mass.
The coherent process is best identified by a low value of jtj (jtj ≲ ℏ 2 =R 2 , where R is the nuclear radius). However, in NC interactions jtj cannot be determined, because the outgoing neutrino momentum cannot be measured. Alternatively, distinct characteristics of the coherent process can be used to separate coherent from background π 0 production. First, the coherent π 0 production has no other particles in the final state and little vertex activity, while background processes often produce additional nucleons or pions and have larger energy depositions near the neutrino interaction vertex. Second, the π 0 's from coherent production are distinctly forward going. A region with enhanced coherent signal in the 2D space of reconstructed π 0 energy and scattering angle can be defined. The coherent signal is measured as an excess of data events over the background prediction in this region.
This paper reports a measurement of NC coherent π 0 cross section using the fine-grained sampling of neutrino interactions in a predominantly carbon tracking medium afforded by the NOvA near detector (ND) [27] exposed to the off-axis flux of the NuMI beam [28] at Fermilab. The flux-averaged cross section is defined as where N data and N bkg are the number of selected data and simulation-predicted background events, respectively, ϵ is the efficiency of the coherent signal selection calculated from simulation, N target is the number of target nuclei in the detector fiducial volume, and ϕ is the integrated neutrino flux. The NOvA ND consists of 193 metric tons of a fully active tracking calorimeter constructed from polyvinyl chloride (PVC) cells filled with liquid scintillator. The liquid scintillator is 62% of the fiducial mass. The target nuclei for neutrino interactions are predominantly carbon (66.7% by mass), chlorine (16.1%), and hydrogen (10.8%), with small contributions from titanium (3.2%), oxygen (3.0%), and other nuclei. Each cell is 3.9 cm wide, 6.6 cm deep, and 3.9 m in length. Cells are arranged in planes alternating between horizontal and vertical orientations to provide three-dimensional reconstruction of neutrino interactions. The fully active volume of the detector is 12.8 m in length, consisting of 192 contiguous PVC planes with 96 cells each. Each plane is approximately 0.18 radiation lengths. Downstream of the fully active volume is a muon range stack with ten layers of 10-cm-thick steel plates interleaved with pairs of one vertical and one horizontal scintillator plane to enhance muon containment. Scintillation light generated by charged particles passing through a cell is captured by a wavelength-shifting fiber connected to a Hamamatsu avalanche photodiode (APD) [29] at the end of the cell. The APD signals are continuously digitized, and those above a preset threshold are recorded with associated time and charge.
The NuMI neutrino beam is produced by colliding 120 GeV protons from the main injector accelerator on a 1.2-m-long graphite target. Charged hadrons produced in the target are focused by two magnetic horns downstream of the target to select positive mesons which then decay into neutrinos in a 675-m-long decay pipe. This analysis uses data corresponding to 3.72 × 10 20 protons on target (POT). The neutrino beam is simulated by FLUKA [30] and the FLUGG [31] interface to GEANT4 [32]. External thin-target hadron production measurements are used to correct and constrain the neutrino flux via the PPFX package developed for the NuMI beam by the MINERvA Collaboration [33].
The NOvA ND is 1 km from the neutrino source, 100 m underground, and on average 14.6 mrad away from the central axis of the neutrino beam. The neutrino flux seen in the NOvA ND is a narrow-band beam peaked at 1.9 GeV, with 68% of the flux between 1.1 and 2.8 GeV and a mean of 2.7 GeV due to the high-energy tail. The neutrino beam in the 0-120 GeV energy region is predominantly ν μ (91%), with a small contamination from ν e (1%) and antineutrinos (8%). In this measurement, the effect of antineutrinos in the flux is accounted for using the simulation to give a solely neutrino-induced result. The predicted integrated neutrino flux from 0 to 120 GeV in the detector volume used in this analysis is ϕ ν ¼ 123.2 AE 11.6 neutrinos=cm 2 =10 10 POT.
Neutrino interactions in the detector are simulated by the GENIE 2.10.4 neutrino event generator [4] except DFR. The Rein-Sehgal PCAC-based model is used to simulate the coherent process. To simulate NC RES and DIS events, the two major background contributions, the Rein-Sehgal model for baryon-resonance production [34] and the Bodek-Yang model [35], are used. The only DFR model implemented in GENIE is the Rein model [36]. However, the Rein model is valid only for the hadronic invariant mass W > 2 GeV region, which is too high for the energy range of NOvA. To simulate the DFR background for this measurement, events are first generated by the Rein model in GENIE (v2.12.2) and then reweighted to an estimation based upon the PCAC calculation by Kopeliovich et al., which includes both DFR and non-DFR contributions [37,38]. In this estimation, Kopeliovich's prediction of dσ dðjtj−jtj min Þ for inclusive νp → νpπ 0 is fitted with GENIE (without DFR) and an exponential term, where jtj min is the minimum possible value of jtj. The exponential term extracted is considered to be the maximum possible contribution from DFR. The DFR events are simulated independently from RES and DIS, with the interference between them neglected, the effect of which is estimated and taken as systematic uncertainty. More details on the DFR modeling can be found in the Appendix. The cross section and selection efficiency of the DFR model used in this measurement are also provided in the Appendix, so that alternative models may be applied to estimate the impact on this measurement.
The nuclear model used in the simulation is the Bodek-Richie relativistic Fermi gas model with short-range nucleon-nucleon correlations [39,40]. Final-state interactions of hadrons inside the nucleus are simulated in GENIE using an effective intranuclear cascade model [4]. GEANT4 [32] is used to simulate the detector's response to the final-state particles from neutrino interactions. The propagation of photons produced by the simulated energy depositions, the response of the APDs, and the digitization of the resulting waveform is accomplished with a custom simulation package.
In both data and simulation, the recorded cell signals (hits) in the NOvA detector are first collected into groups by their space and time information. Each collection of hits is assumed to come from a single neutrino interaction. The intersection of the particle paths found in the collection using a Hough transform [41] are taken as seeds to find the interaction vertex. Hits are further clustered into "prongs" with defined start points and directions emanating from the vertex. Each prong contains hits attributed to one particle.
The events selected by this analysis are required to have exactly two reconstructed prongs contained in the fully active volume of the detector, both identified as electromagneticlike showers by log-likelihood functions based upon dE=dx information in both the longitudinal and transverse directions of the prongs [42,43]. A convolutional neural network trained for NC-CC separation [44] is used to reject CC events. The energy of the prong is calculated as the sum of the calibrated energy deposited in each cell. The invariant mass is calculated from the momenta and opening angle of the reconstructed prongs assuming both are photons, as where E γ1 and E γ2 are the energies of the two prongs and θ γγ is the opening angle between them. The energy scales of data and simulated events are tuned independently so that the mass peaks match the π 0 mass (134.977 MeV=c 2 ) [45].
Only events with reconstructed π 0 mass between 85 and 185 MeV=c 2 are selected to reduce backgrounds. The momenta of the two reconstructed prongs are summed up to obtain the reconstructed momentum of the π 0 . As shown in Fig. 1, the selected events are high-purity NC π 0 's (90%), including both coherent signal and background arising from NC RES and DIS, with small contributions from DFR π 0 production and other interactions.
The background events may have extra energy, especially in the vertex region, but not enough to be reconstructed as prongs. To better control the background, the NC π 0 sample is further divided into two subsamples using kinematic variables: the ratio of the calorimetric energy included in the reconstructed π 0 to the total energy in the event (E π 0 =E Tot ) and the energy in the vertex region defined as the first eight planes from the reconstructed interaction vertex (E Vtx ). The signal-enhanced sample is defined as events with most of their energy in the two photon prongs (E π 0 =E Tot > 0.9) and low vertex energy (E Vtx < 0.3 GeV) to include most of the coherent signal and reduce background. The rest of the events are defined as a control sample, dominated by π 0 's produced by RES and DIS interactions. The signal and control sample selection is shown in Fig. 2.
The control sample data are used to constrain the background prediction. The simulated distributions of RES and DIS events in the π 0 energy and angle (cos θ with respect to the average beam direction) 2D space are used as templates and scaled to fit the control sample data. RES and DIS have distinct π 0 energy and angle distributions, and together they account for approximately 90% of the total background. The fitting parameters are the normalization factors of the templates. The other background components are kept fixed in the fit. The fit results in an increase of the selected RES background by 17.5 AE 6.2% and a decrease in the DIS background by 43.1 AE 13.8%. The two fitting parameters are strongly anticorrelated. The fit result is applied as a renormalization to the background in the signal sample. It also provides a constraint on the systematic sources affecting backgrounds, which will be discussed later. The energy and angle of the π 0 's in the control sample and the signal sample with the renormalized backgrounds are shown in Figs. 3 and 4. There are notable discrepancies between the signal sample data and simulation, especially in the π 0 angular distribution (Fig. 4, right). The θ π 0 spectrum in the data favors production at angles closer to the beam direction than does the simulation, suggesting that the extrapolation from the Q 2 ¼ 0 PCAC approximation to nonzero Q 2 values in the Rein-Sehgal model needs refinement. Similar discrepancies in pion angular distributions have been reported by the MINERvA experiment in recent measurements of chargedcurrent coherent pion production [46,47]. Further study of systematic uncertainties is ongoing to quantitatively address the discrepancies.
A coherent region in the 2D π 0 energy and angle space is defined as those bins with > 15% predicted coherent π 0 signal purity (Fig. 5, left). The selection is intentionally set loosely to reduce potential systematic uncertainties caused by the discrepancies in the π 0 kinematic distributions mentioned previously. The invariant mass of the signal sample events is shown in Fig. 5, right. The signal selection efficiency is 4.1% according to simulation. Figure 6 shows the selection efficiency as a function of Q 2 along with the GENIE predicted signal Q 2 shape. Alternative coherent models may be applied to estimate the impact on this measurement with the selection efficiency provided.
The normalized background in this coherent region is subtracted from data to obtain the number of measured    Left: Ratio of coherent π 0 signal to total simulated events in the signal sample in the 2D space of π 0 energy and cos θ. The region inside the lines is the coherent region defined as bins with > 15% of total simulated events being coherent π 0 . Right: π 0 invariant mass of the signal sample events in the coherent region as described by the left plot with the background normalized by the control sample data. Vertical lines with arrows show the range of invariant masses accepted into the analysis. signal events. The number of simulated signal events is then normalized to the number extracted from the data. The calculation is iterated until the resulting changes in the estimated signal and background populations become negligible. The outcome of this procedure is the coherent signal content, estimated to be 977 AE 67 (stat) events. Neutrino-and antineutrino-induced coherent π 0 's are indistinguishable in this measurement. GENIE predicts 94% of the signal being neutrino induced. This percentage is used to correct the measurement to solely neutrino induced.
The systematic uncertainties for this analysis arise from the calorimetric energy scale, background modeling, coherent signal modeling, detector response to photon showers, detector simulation, particles entering the detector from external sources, and the simulation of neutrino flux. Datadriven methods are used wherever possible to establish the uncertainties.
The calorimetric energy scale is constrained to within 1% by the π 0 invariant mass distributions of simulation and data which corresponds to a 3.4% uncertainty on the cross section measurement. The background-related uncertainty is constrained by the control sample data through the template fit method. The variations that can arise with the template fit to background are estimated by varying the background-modeling parameters within their AE1σ ranges as assigned by GENIE and then repeating the template fit. The uncertainty from each background-modeling parameter is defined as the maximum change in the measured signal events. To estimate the uncertainty from DFR modeling, the template fit is repeated with DFR added as an additional template with its normalization allowed to float. The fit gives DFR normalization factor 0.0 þ 0.8, which favors no DFR contribution. −100% uncertainty on DFR is assigned based upon this study. Another potential uncertainty source from DFR is the interference between DFR and other pion production channels on hydrogen.
To estimate this effect, an additional systematic variation is created by reweighting all the νp → νpπ 0 events on hydrogen to the Kopeliovich model prediction, including both DFR and non-DFR. The vertex energy cut used to define the signal sample and control sample is subject to nuclear effects which are not well modeled by GENIE. To check the impact, the control sample is redefined by applying the same vertex energy (E Vtx < 0.3 GeV) as the signal sample so that the effect of potential mismodeling of vertex energy cancels out. The resulting difference in the measurement from the nominal is added to the systematic uncertainty from GENIE background modeling.
The uncertainty in the coherent signal modeling results in an uncertainty of the efficiency correction. This effect is evaluated by varying the modeling parameters in the Rein-Sehgal model: axial mass (M A , AE50%) and nuclear radius (R 0 , AE20%) [2][3][4]. To check the effect of the discrepancies in π 0 kinematic distributions on the total cross-section measurement, a test is performed by reweighting the simulated signal to data and comparing to the result obtained before reweighting. A 1% difference is found, which is negligible compared to the signal modeling uncertainty assigned. Bremsstrahlung showers induced by energetic muons from external sources provide a data-driven constraint on the simulation of detector response to photon showers. Those bremsstrahlung showers are identified and the muons are removed to create a single photon control sample in the data and simulation [48]. The sample is subject to the same selection cuts as the π 0 photons, and the uncertainty is evaluated as the 1% difference between the data and simulation in selection efficiency. Lastly, the neutrino flux uncertainty comes from beam focusing and hadron production with external thintarget hadron production data constraints applied [33]. The systematic sources and uncertainties are summarized in Table I. The dominant sources are background modeling and flux uncertainties. The total systematic uncertainty is estimated to be 16.6%.
The flux-averaged cross section of NC coherent π 0 production in this measurement is calculated using Eq. (2). The measured cross section is σ ¼ 13.8 AE FIG. 6. Selection efficiency of NC coherent π 0 signal as a function of Q 2 . The GENIE predicted signal Q 2 shape is shown in gray with arbitrary normalization.
where A i is the atomic number of each chemical element (excluding hydrogen) and n i n Tot is its fraction to the total number of nuclei in the fiducial volume. The factor A 2=3 i is an approximate cross-section scaling between different nuclei in accordance with the Berger-Sehgal model [5]. Other models may differ in the prediction of A dependence of coherent pion production. Figure 7 and Table II show this measurement together with other measurements and the GENIE prediction. All measurements in Fig. 7 are scaled to a carbon target by a scale factor of ðA=12Þ 2=3 for the purpose of comparison. The flux-averaged NC coherent π 0 cross section of this work is in agreement with the cross-section prediction of the Rein-Sehgal model (GENIE implementation), although some discrepancies in the π 0 kinematic distributions are observed. This result is the most precise measurement of NC coherent π 0 production in the few-GeV neutrino energy Neutrino Energy (GeV) 10 20 30 /Nucleus) NOvA Aachen-Padova [18] Gargamelle [19] CHARM [21] SKAT [22] NOMAD [25] MINOS [26] Neutrino Energy (GeV) TABLE II. Summary of NC coherent π 0 measurements. The effective atomic number (A) and the average neutrino energy (hE ν i) are shown for each experiment. The results are reported as total cross section per nucleus, crosssection ratios to inclusive ν μ -CC or to the prediction of Rein-Sehgal model.  [23] and SciBooNE [24]) could be considerably different from the GENIE implementation used by the NOvA measurement. A comparison of the Rein-Sehgal predictions of CC coherent in different generators can be found in Ref. [46]. region and the first such measurement on a carbondominated target in this energy range. It benefits both current and future long-baseline neutrino oscillation experiments in background prediction with reduced uncertainty.

APPENDIX: DIFFRACTIVE PION PRODUCTION
NC DFR pion production on free protons (hydrogen) is a background process to the coherent signal. It produces a forward-going pion with small momentum transfer to the recoil proton and becomes indistinguishable from coherent if the recoil proton is undetected. The recoil protons, when detected, could create additional prongs, increase the vertex energy, or decrease the ratio of E π 0 =E Tot , causing the DFR events failing the selection cuts as a result. The acceptance of DFR, therefore, depends upon the kinetic energy of the recoil proton (T p ), which is related to jtj by T p ¼ jtj=2m p . The DFR selection efficiency in this measurement is shown in Fig. 8 as a function of jtj. It is notable that the selection efficiency decreases with jtj, since the proton energy increases with jtj, and the overall efficiency (1.7%) is considerably lower than the coherent signal (4.1%).
DFR is simulated by the Rein model in GENIE 2.12.2 and is predicted to be about 20% of the coherent cross section on carbon in the few-GeV energy region (Fig. 9). However, the Rein model is intended for the hadronic invariant mass W > 2 GeV region. In the W < 2 GeV region, the interference between DFR and RES or non-RES pion productions makes the performance of the Rein model questionable. Alternatively, the DFR cross section can be estimated by a similar method as in Ref. [47] from the calculation of inclusive νp → νpπ 0 by Kopeliovich et al., which is based upon the PCAC hypothesis and includes both DFR and non-DFR contributions [37,38]. In Fig. 10, the left shows Kopeliovich's predictions of νp → νpπ 0 in dσ djtj at NOvA's average neutrino energy (2.7 GeV). This prediction is compared with the GENIE 2.12.2 prediction of νp → νpπ 0 without DFR, and an enhancement is observed in the low-jtj region. A similar enhancement can be observed at low Q 2 (Fig. 10, right). The DFR contribution to νp → νpπ 0 can be quantified by fitting Kopeliovich's prediction of dσ dðjtj−jtj min Þ with GENIE without DFR plus an exponential term A Ã expð−Bðjtj − jtj min ÞÞ, where A and B are fitting parameters (Fig. 11, left). dσ dðjtj−jtj min Þ is used instead of dσ djtj , since the DFR dσ dðjtj−jtj min Þ follows an exponential form, while dσ djtj deviates from an exponential at low jtj because of the jtj min suppression.
The exponential term extracted from the fit is considered as the maximum possible cross section of DFR, since it may include other contributions to the low-jtj enhancement in the Kopeliovich prediction in addition to DFR. It shows a slightly softer shape in jtj − jtj min than the Rein model prediction (Fig. 11, right). An integral over the exponential gives the estimation of the total DFR cross section at E ν ¼ 2.7 GeV: 3.04 × 10 −40 cm 2 =proton. For the measurement reported in this paper, the DFR background events are first simulated by the Rein model in GENIE 2.12.2 and then reweighted to the estimation from Kopeliovich in both normalization and shape as a function of jtj − jtj min . This reweighting makes a 1% difference in the coherent signal measurement. The above method simulates DFR independently from other pion production channels on hydrogen simulated by GENIE. The interference between DFR and non-DFR channels, however, could potentially affect the rate and shape of both DFR and non-DFR pion productions, and the effect on the measurement reported in the paper needs to be discussed. Since the Kopeliovich model includes both DFR and non-DFR contributions in a coherent way, this effect can be taken into account by simulating all the νp → νpπ 0 events on hydrogen using the Kopeliovich model. This is achieved by reweighting the GENIE simulated νp → νpπ 0 events on hydrogen to the prediction of Kopeliovich as a 2D function of jtj and Q 2 . The background template fit is repeated with the hydrogen contribution fixed. The effect on the measurement is a 2.6% difference from the nominal result, which is considered as an additional systematic uncertainty contribution from DFR. GeV fitted with the GENIE 2.12.2 prediction of the same channel without DFR and an exponential term. The fitted exponential term is considered to be the maximum estimation of DFR dσ dðjtj−jtj min Þ . Right: A shape comparison between the fitted exponential term from Kopeliovich and the term extracted from the Rein model prediction for DFR. A slightly softer shape is observed for the one from Kopeliovich.