---
abstract: |
    [\
    ]{} Heavy neutrinos with additional interactions have recently been proposed as an explanation to the MiniBooNE excess. These scenarios often rely on marginally boosted particles to explain the excess angular spectrum, thus predicting large rates at higher-energy neutrino-electron scattering experiments. We place new constraints on this class of models based on neutrino-electron scattering sideband measurements performed at MINER$\nu$A and CHARM-II. A simultaneous explanation of the angular and energy distributions of the MiniBooNE excess in terms of heavy neutrinos with light mediators is severely constrained by our analysis. In general, high-energy neutrino-electron scattering experiments provide strong constraints on explanations of the MiniBooNE observation involving light mediators.
author:
- 'Carlos A. Argüelles'
- Matheus Hostert
- 'Yu-Dai Tsai'
bibliography:
- 'main.bib'
title: |
    Testing New Physics Explanations of MiniBooNE Anomaly\
    at Neutrino Scattering Experiments
---

**Introduction –** Non-zero neutrino masses have been established in the last twenty years by measurements of neutrino flavor conversion in natural and human-made sources, including long- and short-baseline experiments. The overwhelming majority of data supports the three-neutrino framework. Within this framework, we have measured the mixing angles that parametrize the relationship between mass and flavor eigenstates to few-percent-level precision [@Esteban:2018azc]. The remaining unknowns are the absolute scale of neutrino masses and their origin, the CP-violating phase, and the mass ordering of the neutrinos. Nevertheless, anomalies in short-baseline accelerator and reactor experiments [@Athanassopoulos:1996jb; @Aguilar:2001ty; @AguilarArevalo:2007it; @Aguilar-Arevalo:2018gpe] challenge this framework and are yet to receive satisfactory explanations. Minimal extensions of the three-neutrino framework to explain the anomalies introduce the so-called sterile neutrino states, which do not participate in Standard Model (SM) interactions in order to agree with measurements of the Z-boson invisible decay width [@ALEPH:2010aa]. Unfortunately, these minimal scenarios are disfavoured as they fail to explain all data [@Collin:2016aqd; @Capozzi:2016vac; @Dentler:2018sju; @Diaz:2019fwt]. This has led the community to explore non-minimal scenarios. Along this direction, it is interesting to study well-motivated neutrino-mass models that can also explain the short-baseline anomalies and are testable in the laboratory. In this work, we will investigate a class of neutrino-mass-related models that have been proposed as an explanation of the anomalous observation of $\nu_e$-like events in MiniBooNE [@Aguilar-Arevalo:2018gpe].

MiniBooNE is a mineral oil Cherenkov detector located in the Booster Neutrino Beam (BNB), at Fermilab [@AguilarArevalo:2008yp; @AguilarArevalo:2008qa]. Using data collected between 2002 and 2017, the experiment has observed an excess of $\nu_e$-like events that is currently in tension with the standard three-neutrino prediction and is beyond statistical doubt at the $4.7 \sigma$ level [@Aguilar-Arevalo:2018gpe]. While it is possible that the excess is fully or partially due to systematic uncertainties or SM backgrounds (see, *e.g.*, [@AguilarArevalo:2008rc; @Aguilar-Arevalo:2012fmn; @Hill:2010zy]), many Beyond the Standard Model (BSM) explanations have been put forth. These new physics (NP) scenarios typically require the existence of new particles, which can: participate in short-baseline oscillations [@Murayama:2000hm; @Strumia:2002fw; @Barenboim:2002ah; @GonzalezGarcia:2003jq; @Barger:2003xm; @Sorel:2003hf; @Barenboim:2004wu; @Zurek:2004vd; @Kaplan:2004dq; @Pas:2005rb; @deGouvea:2006qd; @Schwetz:2007cd; @Farzan:2008zv; @Hollenberg:2009ws; @Nelson:2010hz; @Akhmedov:2010vy; @Diaz:2010ft; @Bai:2015ztj; @Giunti:2015mwa; @Papoulias:2016edm; @Moss:2017pur; @Carena:2017qhd], change the neutrino propagation in matter [@Liao:2016reh; @Liao:2018mbg; @Asaadi:2017bhx; @Doring:2018cob], or be produced in the beam or in the detector and its surroundings [@Gninenko:2009ks; @Gninenko:2010pr; @Dib:2011jh; @McKeen:2010rx; @Masip:2012ke; @Masip:2011qb; @Gninenko:2012rw; @Magill:2018jla]. These models either increase the conversion of muon- to electron-neutrinos or produce electron-neutrino-like signatures in the detector, where in the latter category one typically exploits the fact that the LSND and MiniBooNE are Cherenkov detectors that cannot distinguish between electrons and photons. Although many MiniBooNE explanations lack a connection to other open problems in particle physics, recent models [@Bertuzzo:2018ftf; @Bertuzzo:2018itn; @Ballett:2018ynz; @Ballett:2019cqp; @Ballett:2019pyw] are motivated by neutrino-mass generation via hidden interactions in the heavy-neutrino sector. In particular, a common prediction of these models is the upscattering of a light neutrino into a heavy neutrino, usually with masses in the tens to hundreds of MeV, which subsequently decays into a pair of electrons. To reproduce the MiniBooNE excess angular distribution either the heavy neutrino must have moderate boost factors and the pair of electrons produced need to be collimated [@Bertuzzo:2018itn], or the heavy neutrino two-body decays must be forbidden [@Ballett:2018ynz].

In this article, we introduce new techniques to probe models that rely on the ambiguity between photons and electrons to explain the MiniBooNE observation, using the dark neutrino model from [@Bertuzzo:2018itn; @Bertuzzo:2018ftf] as a benchmark scenario. Our analysis extends to all models with new marginally boosted particles produced in coherent-like neutrino interactions, as they predict large number of events at higher energies [@Gninenko:2009ks; @Gninenko:2010pr; @Dib:2011jh; @McKeen:2010rx; @Masip:2012ke; @Masip:2011qb; @Gninenko:2012rw; @Magill:2018jla]. Thus, our analysis uses high-energy neutrino-electron scattering measurements [@Auerbach:2001wg; @Deniz:2009mu; @Bellini:2011rx; @Park:2013dax; @Valencia:2019mkf; @Park:2015eqa; @Valencia-Rodriguez:2016vkf; @DeWinter:1989zg; @Geiregat:1992zv; @Vilain:1994qy]. This process is currently used to normalize the neutrino fluxes, due to its well-understood cross section, and has been a fertile ground for light NP searches [@Pospelov:2017kep; @Lindner:2018kjo; @Magill:2018tbb]. Here, however, we expand the capability of these measurements to probe BSM-produced photon-like signatures, by developing a new analysis using previously neglected sideband data. Our technique is complementary to recent searches for coherent single-photon topologies [@Abe:2019cer]. Since the upscattering process has a threshold of tens to hundreds of MeV, we focus on two high-energy neutrino experiments:  [@Park:2013dax; @Valencia:2019mkf; @Park:2015eqa; @Valencia-Rodriguez:2016vkf], a scintillator detector in the Neutrinos at the Main Injector (NuMI) beamline at Fermilab, and CHARM-II [@DeWinter:1989zg; @Geiregat:1992zv; @Vilain:1994qy], a segmented calorimeter detector at CERN along the Super Proton Synchrotron (SPS) beamline. These experiments are complementary in the range of neutrino energies they cover and have different background composition. In all cases a relevant sideband measurement exists, allowing us to take advantage of the excellent particle reconstruction capabilities of and the precise measurements at CHARM-II to constrain NP.

![[*Upscattering cross section compared to the quasi-elastic.*]{} The quasi-elastic cross section on Carbon ($6p^+$) is shown as a function of the neutrino energy (solid black line). The coherent (solid blue) and incoherent (dashed blue) scattering NP cross sections are also shown for the benchmark point of [@Bertuzzo:2018itn]. In the background, we show the BNB flux of $\nu_\mu$ at MiniBooNE (light gray), and the NuMI beam neutrino flux at MINER$\nu$A for the LE (light golden) and ME (light blue) runs in neutrino mode.\[fig:cross\_section\]](cross_sections.pdf){width="49.00000%"}

**Model –** We consider a minimal realisation of dark neutrino models [@Bertuzzo:2018ftf; @Bertuzzo:2018itn; @Ballett:2018ynz; @Ballett:2019cqp; @Ballett:2019pyw] that can explain MiniBooNE. This comprises of one Dirac heavy neutrino[^1], $\nu_4$, with its associated flavor state, $\nu_D$. The dark neutrino $\nu_D$ is charged under a new local U$(1)^\prime$ gauge group, which is part of the particle content and gauge structure needed for mass generation. The dark sector is connected to the SM in two ways: kinetic mixing between the new gauge boson and hypercharge, and neutrino mass mixing. We start by specifying the kinetic part of the NP Lagrangian $$\mathscr{L}_{\rm kin} \supset
\;\; \frac{1}{4} \hat{Z}^{\prime}_{\mu \nu} \hat{Z}^{\prime \mu \nu} + \frac{\sin{\chi}}{2} \hat{Z}^{\prime}_{\mu \nu} \hat{B}^{\mu \nu} + \frac{m_{\hat{Z}^\prime}^2}{2} \hat{Z}^{\prime \mu} \hat{Z}^\prime_{\mu},$$ where $\hat{Z}^{\prime \mu}$ stands for the new gauge boson field, $\hat{Z}^{\prime \mu\nu}$ its field strength tensor, and $\hat{B}^{\mu \nu}$ the hypercharge field strength tensor. After usual field redefinitions [@Chun:2010ve], we arrive at the physical states of the theory. Working at leading order in $\chi$ and assuming $m_{Z^\prime}^2/m_{Z}^2$ to be small, we can specify the relevant interaction Lagrangian as $$\mathscr{L}_{\rm int} \supset \;\;g_D \overline{\nu}_D \gamma_\mu \nu_D Z^{\prime \mu}
 + e \varepsilon Z'^{\mu}J^{\rm EM}_{\mu},$$ where $J^{\rm EM}_{\mu}$ is the SM electromagnetic current, $g_D$ is the U$(1)^\prime$ gauge coupling assumed to be $\mathcal{O}(1)$, and $\varepsilon \equiv c_{\rm w} \chi$, with $c_{\rm w}$ being the cosine of the weak angle. Additional terms would be present at higher orders in $\chi$ and mass mixing with the SM $Z$ is also possible, though severely constrained. After electroweak symmetry breaking, $\nu_D$ is a superposition of neutrino mass states. The flavor and mass eigenstates are related via $$\nu_\alpha = \sum^{4}_{i=1} U_{\alpha i}\nu_{i}, \quad (\alpha=e,\mu,\tau,D),$$ where $U$ is a $4\times4$ unitary matrix. It is expected that $|U_{\alpha 4}|$ is small for $\alpha = e, \mu, \tau$, but $|U_{D4}|$ can be of $\mathcal{O}(1)$ [@Parke:2015goa; @Collin:2016aqd].

![*New physics prediction at ME and CHARM-II.* Neutrino-electron scattering data in $dE/dx$ at (top) and in $E\theta^2$ at CHARM-II (bottom). Error bars are too small to be seen. For both experiments, we show the $\nu-e$ signal and total background prediction quoted (after tuning at MINER$\nu$A), as well as the NP prediction (divided by 10 at CHARM-II). The cuts in our analysis our shown as vertical lines. \[fig:NP\_events\]](both_cartoon.pdf){width="50.00000%"}

**MiniBooNE signature and region of interest–** The heavy neutrino is produced from an active flavour state upscattering on a nuclear target $A$, $\nu_\alpha A \to \nu_4 A$. The upscattering cross section is proportional to $\alpha_D \alpha_\textsc{qed}\varepsilon^2 |U_{\alpha 4}|^2$, dominated by $|U_{\mu 4}|$ since all current accelerator neutrino beams are composed mainly of muon neutrinos. This production can happen off the whole nucleus in a coherent way or off individual nucleons. For $m_{Z^\prime} \lesssim 100$ MeV, the production will be mainly coherent, but for heavier masses, such as the ones considered in [@Ballett:2018ynz], incoherent upscattering dominates. In Fig. \[fig:cross\_section\], we show the NP cross section at the benchmark point of [@Bertuzzo:2018itn] and compare it with the quasi-elastic cross section. By superimposing the cross section on the neutrino fluxes of and MiniBooNE, we make it explicit that the larger energies at and CHARM-II are ideal to produce $\nu_4$. Once produced, $\nu_4$ predominantly decays into a neutrino and a dielectron pair, $\nu_4 \to \nu_\alpha e^+ e^-$, either via an on-shell [@Bertuzzo:2018itn] or off-shell [@Ballett:2018ynz] $Z^\prime$ depending on the choice of $m_4$ and $m_{Z^\prime}$. In this work, we restrict our discussion to the $m_4 > m_{Z^\prime}$ case, where the upscattering is mainly coherent and is followed by a chain of prompt two body decays $\nu_4 \to \nu_\alpha (Z^\prime \to e^+ e^-)$. The on-shell $Z^\prime$ is required to decay into an overlapping $e^+e^-$ pair, setting a lower bound on its mass of a few MeV. Experimentally, however, $m_{Z^\prime} > 10$ MeV for $e \epsilon \sim 10^{-4}$ to avoid beam dump constraints [@Bauer:2018onh]. Increasing $m_{Z^\prime}$ increases the ratio of incoherent to coherent events and makes the electron pair less overlapping. Even though we focus on overlapping $e^+e^-$ pairs, we note that a significant fraction of events would appear as well-separated showers or as a pair of showers with large energy asymmetry, similarly to neutral current (NC) $\pi^0$ events. The asymmetric events also contribute to the MiniBooNE excess and offer a different target for searches in $\nu-e$ scattering data.

A fit to the neutrino energy spectrum at MiniBooNE was performed in [@Bertuzzo:2018itn] and is reproduced in [Fig. \[fig:final\_plot\]]{}. We have performed our own fit to the MiniBooNE energy spectrum using the data release from [@Aguilar-Arevalo:2018gpe], and our results agree with [@Bertuzzo:2018itn], when we simulate the signal at MiniBooNE and the analysis cuts in the same way. This fit leads to preferred values of $m_4$ close to 100 MeV and $|U_{\mu 4 }| \sim 10^{-4}$. Unfortunately, this energy-only fit neglects the distribution of the excess events as a function of their angle $\theta$ with respect to the beam. This is important, as the total observed excess contains only $\approx 50\%$ of the events in the most forward bin ($0.8 < \cos{\theta} < 1.0$), with a statistical uncorrelated uncertainty of 5% on this quantity.

As was recently pointed out in [@Jordan:2018qiy], few NP scenarios can reproduce the angular distribution of the MiniBooNE excess. Among these are models where new unstable particles are produced in inelastic collisions in the detector, such as the present case. Here, large $\theta$ can be achieved by tweaking the mass of the heavy neutrino; the signal becomes less forward as $\nu_4$ becomes heavier. To show this, we use our dedicated Monte Carlo (MC) simulation to asses the values of $m_4$ preferred by MiniBooNE data [^2]. For $m_{Z^\prime} = 30$ MeV and $m_4 = 100$, $200$, and $400$ MeV, we find that 98%, 87%, and 70% of the NP events would lie in the most forward bin, respectively. We find the predicted angular distribution to be more forward than [@Bertuzzo:2018itn] due to an improved MiniBooNE simulation; see Supplementary Material for details. This simulation discrepancy is understood and only strengthens our conclusions. Thus the relevant region for the MiniBooNE angular distribution is $m_4 \gtrsim 400$ MeV for $m_{Z^\prime} = 30$ MeV.

[**Our analysis –**]{} Neutrino-electron scattering measurements predicate their cuts in the following core ideas: no hadronic activity near the interaction vertex, small opening angle from the beam, $E_e \theta^2 \lesssim 2 m_e$, and the requirement that the measured energy deposition, $dE/dx$, be consistent with that of a single electron. For the NP events, when the coherent process dominates and the mass of the $Z^\prime$ is small, the first two conditions are often satisfied. However, the requirement of a single-electron-like energy deposition removes a significant fraction of the new-physics induced events. This presents a challenge, as the NP events are mostly overlapping electron pairs and will potentially be removed by the $dE/dx$ cut. In order to circumvent this problem, we perform our analysis not at the final-cut level, but at an intermediate one. This is done differently for CHARM-II and MINER$\nu$A: the CHARM-II experiment provides data as a function of $E_e \theta^2$ without the $dE/dx$ cut, and provides data as a function of the measured $dE/dx$ after analysis cuts on $E_e \theta^2$.

We have developed our own MC simulation for candidate electron pair events in MiniBooNE, MINER$\nu$A and CHARM-II; see the Supplementary Material for more details on detector resolutions, precise signal definition, and resulting distributions. We only consider the coherent part of the cross section to avoid hadronic-activity cuts, which is conservative. We also select only events with small energy asymmetries and small opening electron angles. When required, we assume the mean $dE/dx$ in plastic scintillator to follow the same shape as the NC $\pi^0$ prediction. Our prediction for new physics events for the BP point is show in Fig. \[fig:NP\_events\] on top of the ME and CHARM-II data and MC prediction.

The CHARM-II analysis is mostly based on Fig. 1 of [@Vilain:1994qy]. This sample is shown as a function of $E\theta^2$ and does not have any cuts on $dE/dx$. It contains all events with shower energies between $3$ and $24$ GeV, and our final cut on $E\theta^2$ is fixed at $28$ MeV. For , the event selection is identical for the low-energy (LE) and medium-energy (ME) analyses [@Park:2015eqa; @Valencia:2019mkf]. The minimum shower energy required is $0.8$ GeV in order to remove the $\pi^0$ background and have reliable angular and energy reconstruction. Events are kept only when they meet the following angular separation criterion: $E_e \theta^2 < 3.2\times 10^{-3}~{\rm ~GeV ~rad^2}$. A final cut is applied, ensuring $dE/dx < 4.5~{\rm MeV} / 1.7~{\rm cm}$. The analyses use the data outside the previous $dE/dx$ cut to constrain backgrounds. This sideband is defined by all events with $E_e\theta^2 > 5 \times 10^{-3} {\rm ~GeV ~rad^2}$ and $dE/dx < 20~{\rm MeV}/1.7~{\rm cm}$. Using this sideband measurement, the collaboration tunes their backgrounds by ($0.76$, $0.64$, $1.0$) for ($\nu_e$CCQE, $\nu_\mu$NC, $\nu_\mu$CCQE) processes in the LE mode. Our LE analysis uses the data shown in Fig. 3 of [@Park:2015eqa] where all the cuts are applied except for the final $dE/dx$ cut. In our final event selection, we require that the sum of the energy deposited be more than $4.5$ MeV$/ 1.7$ cm, compatible with an $e^+e^-$ pair and yielding an efficiency of $90\%$.

The recent ME data contains an excess in the region of large $dE/dx$ [@Valencia:2019mkf], where the NP events would lie. However, this excess is attributed to NC $\pi^0$ events, and grows with the shower energy undershooting the rate require to explain the MiniBooNE anomaly. With normalization factors as large as 1.7, the collaboration tunes primarily the NC $\pi^0$ prediction in an energy dependent way. After tuning, the total NC $\pi^0$ sample corresponds to $20\%$ of the total number of events before the $dE/dx$ cut.

To place our limits, we perform a rate-only analysis by means of a Pearson’s $\chi^2$ as test statistic; detailed definition is given in the Supplementary Material. We incorporate uncertainties in background size and flux normalization as nuisance parameters with Gaussian constraint terms. For the neutrino-electron scattering and BSM signal, we allow the normalization to scale proportionally to the same flux uncertainty parameter. The background term also scales with the flux-uncertainty parameter but has an additional nuisance parameter to account for its unknown size. We obtain our constraint as a function of heavy neutrino mass $m_4$, and mixing $|U_{\mu 4}|$ assuming a $\chi^2$ with two degrees of freedom [@Tanabashi:2018oca].

In our nominal LE (ME) analysis, we allow for 10% uncertainty on the flux [@Aliaga:2016oaz], and 30% (40%) uncertainty on the background motivated by the amount of tuning performed on the original backgrounds. Note that the nominal background predictions in the LE (ME) analysis overpredicts (underpredicts) the data before tuning, and that tuning parameters are measured at the 3% (5%) level [@Park:2013dax; @Valencia:2019mkf]. We also perform a background-ignorant analysis in which we assume 100% uncertainty for the background normalization, which changes our conclusions by only less than a factor of two. This emphasizes the robustness of our bound, since the NP typically overshoots the low number of events in the sideband. For the benchmark point of [@Bertuzzo:2018itn], we predict a total signal of 232 (4240) events for LE (ME).

For CHARM-II, the NP signal lies mostly in a region with small $E\theta^2$. Thus, we constrain backgrounds using the data from $28 < E\theta^2 < 60$ MeV rad$^2$. This sideband measurement constrains the normalization of the backgrounds in the signal region at the level of $3\%$. The extrapolation of the shape of the background to the signal region introduces the largest uncertainty in our analysis. For this reason, we raise the uncertainty of the background normalization from $3\%$ to a conservative $10 \%$ when setting the limits. Flux uncertainties are assumed to be $4.7\%$ and $5.2\%$ for neutrino and antineutrino mode [@Allaby:1987bb], respectively, and are applicable to the new-physics signal, $\nu-e$ scattering prediction, and backgrounds. Uncertainties in the $\nu-e$ scattering cross sections are expected to be sub-dominant and are neglected in the analysis [@deGouvea:2006hfo]. For CHARM-II, the NP also yields too many events in the signal region, namely $\approx 2.2\times10^{5}$ events for the benchmark point of [@Bertuzzo:2018itn] in antineutrino mode. If we lower $|U_{\mu4}| = 10^{-4}$ and $m_4 = 100$ MeV, CHARM-II would still have $\approx 3 \times 10^3$ new physics events.

![[*New constraints on dark neutrinos as a MiniBooNE explanation.*]{} The fit to the MiniBooNE energy distribution from [@Bertuzzo:2018itn] is shown as closed yellow (orange) region for one (three) sigma C.L., together with the benchmark point (${\bf\odot}$). Our constraints are shown at 90% C.L. for LE in blue (solid – 30% background normalization uncertainty, dashed – conservative 100% case), for ME in cyan (solid – 40% background normalization uncertainty, dashed – conservative 100% case), and for CHARM-II in red (solid – 3% background normalization from the sideband constraint, dashed – conservative 10% case). Vertical lines show the percentage of excess events at MiniBooNE that lie in the most forward angular bin. Exclusion from heavy neutrino searches is shown as a hatched background. Other relevant assumed parameters are shown above the plot; changing them does not alter our conclusion.\[fig:final\_plot\]](bounds.pdf){width="3.38in"}

[**Results and conclusions –**]{} The resulting limits on dark neutrinos using neutrino-electron scattering experiments are shown in the $|U_{\mu 4}|$ vs $m_4$ plane at 90% confidence level (CL) in [Fig. \[fig:final\_plot\]]{}. The MiniBooNE fit from [@Bertuzzo:2018itn] is shown, together with vertical lines indicating the percentage of events at MiniBooNE that populate the most forward angular bin. We have chosen the same values of $\varepsilon$, $\alpha_D$, and $m_{Z^\prime}$ as used in [@Bertuzzo:2018itn], and shown their benchmark point ($m_4 = 420$ MeV and $|U_{\mu 4}|^2 = 9 \times 10^{-7}$) as a dotted circle. For these parameters, we can conclude that a good angular distribution at MiniBooNE is in large tension with neutrino-electron scattering data. We note that the MiniBooNE event rate scales identically to our signal rate in all the couplings, and the dependence on $m_{Z^\prime}$ is subleading due to the typical momentum transfer to the nucleus, provided $m_{Z^\prime} \lesssim 100$ MeV . This implies that changing the values of these parameters does not modify the overall conclusions of our work. In addition, for this realization of the model, larger $m_{Z^\prime}$ implies larger values of $m_4$, increasing their impact on neutrino-electron scattering data. Our and CHARM-II results are mutually reinforcing given that they impose similar constraints for $m_4 \lesssim 200 $ MeV. For larger masses, the kinematics of the signal becomes less forward and the production thresholds start being important. This explains the upturns visible in our bounds, where we observe it first in and later in CHARM-II as we increase $m_4$, since CHARM-II has higher beam energy.

We emphasize that our analysis is general, and can be adapted to other models. In fact, any MiniBooNE explanation with heavy new particles faces severe constraints from high-energy neutrino-electron scattering data if the signal is free from hadronic activity. This is realised, for instance, in scenarios with heavy neutrinos with dipole interactions [@Gninenko:2009ks; @Gninenko:2010pr; @Dib:2011jh; @McKeen:2010rx; @Masip:2012ke; @Masip:2011qb; @Gninenko:2012rw; @Magill:2018jla]. Our bounds can also be adapted to other scenarios with dark neutrinos and heavy mediators [@Ballett:2018ynz; @Ballett:2019pyw]. For those, however, we do not expect our bounds to constrain the region of parameter space where the MiniBooNE explanation is viable, since most of the signal at MiniBooNE contains hadronic activity which would be visible at and CHARM-II.

In the near future, our new analysis strategy could be used in the up-coming ME results on antineutrino-electron scattering. The NP cross section, being the same for neutrino and antineutrinos, is thus more prominent on top of backgrounds. This class of analyses will also greatly benefit from improved calculations and measurements of coherent $\pi^0$ production and single-photon emitting processes. This is particularly important given the excess seen in the ME analysis. A new result can also be obtained by neutrino-electron scattering measurements at NO$\nu$A, which will sample a different kinematic regime as its off-axis beam peaks at lower energies and expects fewer NC $\pi^0$ events per ton. Beyond neutrino-electron scattering, the BSM signatures we consider could be lurking in current measurements of $\pi^0$ production, *e.g.*, at MINOS [@Adamson:2016hyz] and MINER$\nu$A [@Wolcott:2016hws] [^3], and in analyses like the single photon search performed by T2K [@Abe:2019cer]. Thus, if dark neutrinos are indeed present in current data, our technique will be crucial to confirm it.

To summarize, a variety of measurements are underway to further lay siege to this explanation of the MiniBooNE observation and, simultaneously, start probing testable neutrino mass generation models, as well as other similar NP signatures.

Acknowledgements {#acknowledgements .unnumbered}
================

We thank Janet Conrad, Kareem Farrag, Alberto Gago, Gordan Krnjaic, Trung Le, Pedro Machado, Kevin Mcfarland, and Jorge Morfin for useful discussions, and Jean DeMerit for carefully proofreading our work. The authors would like to thank Fermilab for the hospitality at the initial stages of this project. Also, the authors would like to thank Fermilab Theory Group and the CERN Theory Neutrino Platform for organizing the conference “Physics Opportunities in the Near DUNE Detector Hall,” which was essential to the completion of this work. CAA would especially like to thank Fermilab Center for Neutrino Physics summer visitor program for funding his visit. MH’s work was supported by Conselho Nacional de Ciência e Tecnologia (CNPq). CAA is supported by U.S. National Science Foundation (NSF) grant No. PHY-1801996. This document was prepared by YDT 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.

**Supplementary Material**

Our analysis discussed in the main text is now described in more detail and all assumptions used in our simulations are summarized. We start by discussing our statistical method, and then discuss our Monte Carlo (MC) simulation, stating our signal definitions more precisely. Later, we show a few kinematical distributions from our dedicated MC simulation, including the angular distributions at MiniBooNE used to obtain the vertical lines in [Fig. \[fig:final\_plot\]]{}. In order to furher aid reproducibility of our result, we also make our Monte Carlo events for some parameter choices available on <span style="font-variant:small-caps;">g</span>it<span style="font-variant:small-caps;">h</span>ub [^4].

Statistical Analysis
====================

Our statistical analysis uses the Pearson-$\chi^2$ as a test statistic, where the expected number of events is scaled by nuisance parameters to incorporate systematic uncertainties. Our test statistic reads

$$\begin{aligned}
\chi^2(\vec\theta, \alpha, \beta) = \frac{ (N_{\rm data} - N_\mathrm{pred}(\vec\theta, \alpha, \beta) )^2 }{ N_\mathrm{pred}(\vec\theta, \alpha, \beta)} + \left(\frac{\alpha}{\sigma_\alpha}\right)^2  +
\left(\frac{\beta}{\sigma_\beta}\right)^2,\end{aligned}$$

with the number of predicted events given by $$\begin{aligned}
N_\mathrm{pred}(\vec\theta, \alpha, \beta) = (1+\alpha + \beta) \mu_{\rm MC}^{\rm BKG} + (1+\alpha) \mu_{\rm MC}^{\nu-e} + (1 + \alpha) \mu_{\rm BSM}(\vec\theta),\end{aligned}$$

where $\vec\theta$ are the model parameters, while $\alpha$ and $\beta$ are nuisance parameters that incorporate uncertainties from the overall rate and the background rate, respectively. Here, $N_{\rm data}$ stands for the total rate observed in the experiment, $\mu_{\rm MC}^{\rm BKG}$ the quoted background rates, $\mu_{\rm MC}^{\nu-e}$ the quoted $\nu-e$ events, and $\mu_{\rm BSM}(\vec\theta)$ the predicted number of BSM events calculated using our MC. We discuss the choice of these systematic uncertainties, namely the choice of $\sigma_\alpha$ and $\sigma_\beta$, when describing the simulation of each experiment below. To obtain our results we use the test statistic profiled over the nuisance parameters, namely $\chi^2(\vec\theta) = \min_{(\alpha, \beta)} \left(\chi^2(\vec\theta, \alpha, \beta) \right)$, and use the test-statistic thresholds given in [@pdg].

Simulation Details
==================

  Experiment    Detector Resolution                                                                                                                 Overlapping                                                          Analysis Cuts
  ------------- ----------------------------------------------------------------------------------------------------------------------------------- -------------------------------------------------------------------- -----------------------------------------------------------------------------------------
  MiniBooNE                                                                                                                                                                                                              
                $\sigma_E/E = 12\%/\sqrt{E_e/{\rm GeV}}$ $\sigma_\theta = 4{^\circ}$                                                                $E_{+} > 30$ MeV $E_{-} > 30$ MeV $\Delta \theta_{\pm} < 13^\circ$   N/A
  MINER$\nu$A                                                                                                                                                                                                            
                $\sigma_E/E = 5.9\%/\sqrt{E_e/{\rm GeV}} + 3.4\%$ $\sigma_\theta = 1{^\circ}$                                                       $E_{+} > 30$ MeV $E_{-} > 30$ MeV $\Delta \theta_{\pm} < 8^\circ$    $E_{\rm vis} > 0.8$ GeV $E_{\rm vis} \theta^2 < 3.2$ MeV $Q^2_{\rm rec} < 0.02$ GeV$^2$
  CHARM-II                                                                                                                                                                                                               
                $\sigma_E/E = 9\%/\sqrt{E/{\rm GeV}} +  11\%$ $\sigma_\theta/{\rm mrad} =  \frac{27 (E/{\rm GeV})^2 +14}{\sqrt{E/{\rm GeV}}} + 1$   $E_{+} > 30$ MeV $E_{-} > 30$ MeV $\Delta \theta_{\pm} < 4^\circ$    $3$ GeV $<E_{\rm vis} < 24$ GeV $E_{\rm vis} \theta^2 < 28$ MeV

We generate events distributed according to the upscattering cross section for the process $\nu_\mu A \to \nu_4 A$, where $A$ is a nuclear target. We only discuss upscattering on nuclei, as the number of incoherent scattering on protons is much smaller for the relevant $Z^\prime$ masses; see [Fig. \[fig:cross\_section\]]{}. We then implement the chain of two-body decays: $\nu_4 \to \nu_\mu Z^\prime$ followed by $Z^\prime \to e^+ e^-$. To go from our MC true quantities to the predicted experimental observables, we perform three procedures. First, we smear the energy and angles of the $e^+$ and $e^-$ originating from the decay of the $Z^\prime$ according to detector-dependent Gaussian resolutions. Next, we select all events with that satisfy the $e^+e^-$ overlapping condition given in [Table \[tab:parameters\]]{}. Namely, if the condition is satisfied they are assumed to be reconstructed as a single electromagnetic (EM) shower. This guarantees that the events behave like a photon shower inside the detector [^5]. Finally, for and CHARM-II, these samples are subject to analysis-dependent kinematic cuts to determine if they contribute to the $\nu-e$ scattering sample. Detector resolutions, requirements for the dielectron pair to be overlapping, and analysis-dependent cuts are summarized in [Table \[tab:parameters\]]{}. We now list the experimental parameters used in our simulations for each individual detector.

#### CHARM-II

The CHARM-II experiment is simulated using the CERN West Area Neutrino Facility (WANF) wide band beam  [@Vilain:1998uw]. The total number of protons-on-target (POT) is $2.5 \times 10^{19}$ for the $\nu$ and $\overline{\nu}$ runs combined. We assume glass to be the main detector material, (SiO$_2$), such that we can treat neutrino scattering off an average target with $\langle Z\rangle=11$ and $\langle A \rangle = 20.7$ [@DeWinter:1989zg; @Vilain:1993sf]. The fiducial volume in our analysis is confined to a transverse area of $320$cm$^2$, which corresponds to a fiducial mass of $547$t, and the detection efficiency is taken to be $76\%$; efficiency for $\pi^0$ sample is quoted at $82\%$ [@Vilain:1992wx]. We reproduce the total number of $\nu-e$ scattering events with $3$ GeV $< E_{\rm vis} <24$ GeV, namely $2677+2752$, to within a few percent level when setting the number of POTs in $\nu$ mode to be $1.69$ of that in the $\overline{\nu}$ mode [@Geiregat:1991md]. We assume a flux uncertainty of $\sigma_\alpha = 4.7\%$ for neutrino, and $\sigma_\alpha = 5.2\%$ for antineutrino beam [@Vilain:1992wx]. The background uncertainty is constrained to be $\sigma_\beta = 3\%$ using the data with $E_{\rm vis} \theta^2 > 28$ MeV, where the number of new physics events is negligible.

#### MINER$\nu$A

For our MINER$\nu$A simulation, we use the low-energy (LE) and medium-energy (ME) NuMI neutrino fluxes [@AliagaSoplin:2016shs]. The total number of POT is $3.43\times 10^{20}$ for LE data, and $11.6\times10^{20}$ for ME data. The detector is assumed to be made of CH, with a fiducial mass of $6.10$ tons and detection efficiencies of $73\%$ [@Parke:2015goa; @Valencia:2019mkf]. We assume a flux uncertainty of $\sigma_\alpha = 10\%$ for both the LE and ME modes [@Aliaga:2016oaz]. Due to the tuning performed in the sideband of interest, the uncertainties on the background rate are much larger. For the LE, we take $\sigma_\beta = 30\%$, while for the ME data $\sigma_\beta = 50\%$. Although tuning is significant for the coherent $\pi^0$ production sample, the overall rate of backgrounds in the sideband with large $dE/dx$ does not vary by more than $20\%$ ($40\%$) in the LE (ME) tuning.

#### MiniBooNE

To simulate MiniBooNE, we use the Booster Neutrino Beam (BNB) fluxes from [@AguilarArevalo:2008yp]. Here, we only discuss the neutrino run, although the predictions for the antineutrino run are very similar. We assume a total of $12.84 \times 10^{20}$ POT in neutrino mode. The fiducial mass of the detector is taken as $450$t of CH$_2$. In order to apply detector efficiencies, we compute the reconstructed neutrino energy under the assumption of CCQE scattering $$\begin{aligned}
 E_\nu^{CCQE} = \frac{E_{\rm vis} m_p}{m_p - E_{\rm vis} (1 - \cos{\theta}) },\end{aligned}$$ where $E_{\rm vis} = E_{e^+} + E_{e^-}$ is the total visible energy after smearing. Under this assumption, we can apply the efficiencies provided by the MiniBooNE collaboration [@Aguilar-Arevalo:2012fmn]. Using our MC we can reproduce well the distributions obtained using the MiniBooNE Monte Carlo data release provided for oscillation analyses.

Kinematic Distributions
=======================

As an important check of our calculation and of the explanation of the MiniBooNE excess within the model of interest, we plot the MiniBooNE neutrino data from 2018 [@Aguilar-Arevalo:2018gpe] against our MC prediction in Suppl. Fig. \[fig:MB\_distributions\]. We do this for three different new physics parameter choices. We set $m_{Z^\prime} = 30$ MeV, $\alpha \epsilon^2 = 2\times10^{-10}$ and $\alpha_D = 1/4$ for all points, but vary $|U_{\mu 4}|^2$ and $m_4$ so that the final number of excess events predicted by the model at MiniBooNE equals 334.

![image](Enu_reco.pdf){width="49.00000%"} ![image](Theta_reco.pdf){width="49.00000%"}

To verify that the new physics signal is important in neutrino-electron studies, we also plot kinematical distributions for the benchmark point (BP) discussed in the main text for different detectors. This corresponds to $m_{Z^\prime} = 30$ MeV, $\alpha \epsilon^2 = 2\times10^{-10}$, $\alpha_D = 1/4$, $|U_{\mu 4}|^2 = 9\times10^{-7}$ and $m_4 = 420$ MeV. The interesting variables are the energy asymmetry of the dielectron pair $$|E_{\rm asym}| = \frac{|E_+ - E_-|}{E_+ + E_-},$$ as well as the separation angle $\Delta \theta_{e^+e^-}$ between the two electrons. These variables are plotted in Suppl. Fig. \[fig:other\_distributions\] at MC truth level, before any smearing or selection takes place. We also plot the total reconstructed energy $E_{\rm vis} = E_{e^+} + E_{e^-}$ and the quantity $E_{\rm vis} \theta^2$, where $\theta$ stands for the angle formed by the reconstructed EM shower and the neutrino beam. The visible energy, $E_{\rm vis}$, and angle, $\theta$, are computed after smearing, but before the selection into overlapping pairs takes place.

![image](Easym_v1.pdf){width="49.00000%"} ![image](DeltaTheta_v1.pdf){width="49.00000%"}\
![image](Evis_v1.pdf){width="49.00000%"} ![image](Etheta2_v1.pdf){width="49.00000%"}

[^1]: Models with the decay of Majorana particles will lead to greater tension with the angular distribution at MiniBooNE due to their isotropic nature [@Formaggio:1998zn; @Balantekin:2018ukw].

[^2]: Since the released MiniBooNE data do not provide the correlation between angle and energy, and their associated systematics, an energy-angle fit is not possible.

[^3]: This $\nu_e$CCQE measurement by observes a significant excess of single photon-like showers attributed to diffractive $\pi^0$ events. These are abundant in similar realizations of this NP model [@Ballett:2018ynz].

[^4]: [[github.com/mhostert/DarkNews]{}](https://github.com/mhostert/DarkNews).

[^5]: For MiniBooNE, we also include events that are highly asymmetric in energy, *i.e.*, $E_{\pm} > 30$ MeV and $E_{\mp} < 30$ MeV, where the most energetic shower defines the angle with respect to the beam.
