# Dipole-Coupled Neutrissimo Explanations of the MiniBooNE Excess Including
Constraints from MINERvA Data

N.W. Kamp  Dept. of Physics, Massachusetts Institute of Technology, Cambridge,
MA 02139, USA  M. Hostert  University of Minnesota, Minneapolis, MN 55455, USA
Perimeter Institute for Theoretical Physics, Waterloo, ON N2J 2W9, Canada  A.
Schneider  Dept. of Physics, Massachusetts Institute of Technology, Cambridge,
MA 02139, USA  S. Vergani  University of Cambridge, Cambridge CB3 0HE, United
Kingdom  C.A. Argüelles  Dept. of Physics, Harvard University, Cambridge, MA
02138, USA  J.M. Conrad  Dept. of Physics, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  M.H. Shaevitz  Dept. of Physics,
Columbia University, New York, NY, 10027, USA  M.A. Uchida  University of
Cambridge, Cambridge CB3 0HE, United Kingdom

###### Abstract

We revisit models of heavy neutral leptons (neutrissimos) with transition
magnetic moments as explanations of the $4.8\sigma$ excess of electron-like
events at MiniBooNE. We first re-examine the preferred regions in the model
parameter space to explain MiniBooNE, considering also potential contributions
from oscillations due to an eV-scale sterile neutrino. We then derive
constraints on the model using neutrino-electron elastic scattering data from
MINERvA. To carry out these analyses, we have developed a detailed Monte Carlo
simulation of neutrissimo interactions within the MiniBooNE and MINERvA
detectors using the LeptonInjector framework. This simulation allows for a
significantly more robust evaluation of the neutrissimo model compared to
previous studies in the literature–a necessary step in order to begin making
definitive statements about beyond the Standard Model explanations of the
MiniBooNE excess. We find that MINERvA rules out a large region of parameter
space, but allowed solutions exist at the $2\sigma$ confidence level. A
dedicated MINERvA analysis would likely be able to probe the entire region of
preference of MiniBooNE in this model.

##  I Introduction

Past, current, and future neutrino experiments offer some of the most
promising avenues for observing physics beyond the Standard Model (BSM). The
discovery of neutrino oscillations, and thus nonzero neutrino mass, is itself
an indication of BSM physics [1, 2], spurring a decades-long oscillation
experimental program spanning orders of magnitude in energy and length scales
[3, 4, 5]. A nearly-consistent three-neutrino mixing paradigm has emerged from
this program; however, anomalous results have also been observed. Two striking
examples are the excess of inverse-beta-decay events at the LSND detector at
the Los Alamos Neutrino Science Center [6] and the excess of electron-like
events observed by the MiniBooNE (MB) experiment at the Fermilab Booster
Neutrino Beam (BNB) [7]. Determining the nature of these excesses is an active
frontier in neutrino physics [8] Historically, both excesses have been
interpreted within the context of a $3+1$ model, in which one introduces an
eV-scale sterile neutrino facilitating short-baseline
$\nu_{\mu}(\overline{\nu}_{\mu})\to\nu_{e}(\overline{\nu}_{e})$ oscillations.
However, models beyond the vanilla $3+1$ scenario may better accommodate these
anomalies within the global experimental landscape, including cosmology and
other neutrino oscillation experiments.

One such model considers MeV-scale heavy neutral leptons (HNLs), or
“neutrissimos” [9, 10], with transition magnetic moments, popularly known as
the dipole portal to HNLs [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24]. Most recently, Ref. [23] showed that the MB anomaly could be fitted with
such a dipole portal extended to include a mass-mixed eV-scale sterile
neutrino. The advantages of this model are twofold: (1) it provides a better
description of the low-energy, forward-angle part of the MB excess compared to
the $3+1$ model alone, and (2) reducing the oscillation-based contribution to
the MB excess alleviates tension in global fits to the $3+1$ model while
retaining an explanation to the LSND anomaly.

For the main result of this paper, we derive constraints on the dipole portal
coupling using existing experimental results from the MINERvA collaboration.
Specifically, MINERvA has performed measurements of the neutrino-electron
elastic scattering ($\nu-e$) rate for the purpose of constraining the Neutrino
Main Injector (NuMI) low-energy (LE) and medium-energy (ME) neutrino fluxes
[25, 26, 27]. Photons produced in the decay of the dipole-coupled neutrino
would mimic the single electromagnetic shower morphology of electron-
scattering (ES) events and would therefore enter as a photon-like (large
$dE/dx$) background in the MINERvA analysis. In addition to deriving new
constraints from MINERvA, we improve upon the MiniBooNE dipole model analysis
performed in Ref. [23].

The analyses described in this article have been carried out using a novel
simulation developed within the LeptonInjector framework [28]. This tool
allows a more robust description of neutrissimo interactions in detector
subsystems as well as the position and kinematics of observable final state
particles. This is an important improvement over previous treatments in the
literature, as an accurate simulation of exotic BSM physics scenarios in
different neutrino detectors will be vital in determining the nature of the
MiniBooNE excess.

The rest of this article is organized as follows. In section II, we review the
dipole-portal sterile neutrino model in more detail. In section III, we
introduce the novel LeptonInjector-based simulation developed for the studies
presented here. In section IV, we refine the preferred regions in dipole
parameter space which explain the energy and angular distributions of the
MiniBooNE excess. In section V, we calculate constraints in dipole parameter
space derived from the MINERvA ES analysis. In section VI, we discuss existing
and projected constraints on the dipole model from current and future neutrino
experiments. Finally, in section VII we discuss the outlook of this mixed
model of oscillation and decay as a solution to the MiniBooNE anomaly in light
of our derived MINERvA constraints.

##  II The Dipole Portal

Figure 1: $\mathcal{N}$ production from $\nu$ upscattering (left) and
$\mathcal{N}$ decay (right).

We introduce a right-handed neutrino $\mathcal{N}$ that couples to the left-
handed neutrino fields via a transition magnetic moment. At the effective
operator level, we have the dimension-six dipole operators,

$\displaystyle\mathscr{L}$
$\displaystyle\supset\frac{1}{\Lambda^{2}}\overline{L_{\alpha}}\widetilde{H}\sigma^{\mu\nu}\mathcal{N}_{R}\left(C_{B}^{\alpha}\,B_{\mu\nu}+C_{W}^{\alpha}\,W^{a}_{\mu\nu}\sigma_{a}\right)+\text{
h.c.},$ (1)

where $\alpha$ is a flavor index, $\widetilde{H}=i\sigma_{2}H^{*}$,
$C_{B}^{\alpha}$, and $C_{W}^{\alpha}$ are Wilson coefficients and $\Lambda$
is the new physics scale. After electroweak symmetry breaking, the dipole
operator gives rise to the electromagnetic transition magnetic moment of
neutrinos,

$\displaystyle\mathscr{L}$ $\displaystyle\supset
d_{\alpha\mathcal{N}}\,\overline{\nu_{\alpha}}\sigma_{\mu\nu}F^{\mu\nu}\mathcal{N}_{R}+\text{
h.c.},$ (2)

where $\nu_{\alpha}$ corresponds to the neutrino Weak eigenstates and
$F^{\mu\nu}$ to the electromagnetic field strength. The dipole parameter is
defined as
$d_{\alpha\mathcal{N}}=(v_{h}/\sqrt{2})(C_{W}^{\alpha}C_{B}^{\alpha}+C_{W}^{\alpha}s_{W})/\Lambda^{2}$,
where $v_{h}$ is the vacuum expectation value of the SM Higgs. The other
transition moments mediated by the $W$ and $Z$ bosons will also be present,
but their low-energy effects are further suppressed by $G_{F}$, and therefore
negligible in our region of interest.

The upscattering signature we are interested in at MiniBooNE is initiated by
muon neutrinos and antineutrinos, so in this work, we always consider the
second-generation coupling $d_{\mu\mathcal{N}}$. In UV-completions of this
operator from one-loop diagrams, the size of the dipole coupling is typically
proportional to $d_{\alpha\mathcal{N}}\sim(1/16\pi^{2})(m_{\beta}/m_{X}^{2})$,
where $m_{\beta}$ is the mass of some charged particle and $X$ is some heavy,
charged scalar, for example. Under the assumption of flavor-conserving
interactions between SM neutrinos and the new physics, one would take
$\beta=\alpha$ and conclude that the transition magnetic moment of the third-
generation neutrinos is much larger. While this need not necessarily be the
case, it still provides enough motivation for us to consider the third-
generation coupling $d_{\tau\mathcal{N}}$, including the case

$\frac{d_{\tau\mathcal{N}}}{d_{\mu\mathcal{N}}}=\frac{m_{\tau}}{m_{\mu}},$ (3)

and neglecting the first-generation couplings altogether.

We note that large neutrino magnetic moments typically imply large Dirac
masses for light neutrinos, $m_{D}\overline{\nu_{L}}\mathcal{N}_{R}$ [29, 30,
31]. In models with heavy neutrinos, this presents two challenges: i) neutrino
masses, schematically given by $m_{\nu}\sim m_{D}^{2}/M$, may be too large,
and ii) $m_{D}$ will generate large mixing between active and heavy neutrinos,
$U_{\alpha\mathcal{N}}\sim m_{D}/M$. Here, $M$ stands for the heavy neutrino
mass scale. The first challenge is easily overcome in models like the inverse-
seesaw where lepton number is approximately conserved and right-handed
neutrino fields combine into pseudo-Dirac $\mathcal{N}$ particles. In that
case, $m_{\nu}$ is controlled by the mass splitting between the Majorana
neutrinos, which may be parametrically small. This is also the case preferred
by the short-baseline phenomenology discussed below, since the single-photons
produced in the decays of Dirac HNLs are less forward. The second challenge,
however, is not so easy to overcome. The mixing between active and heavy
neutrinos remains large even in inverse-seesaw models and the parameter space
in that case is strongly constrained by laboratory limits on
$U_{\alpha\mathcal{N}}$. For instance, decay-in-flight signatures at neutrino
experiments, where $\mathcal{N}$ is copiously produced in meson decays at the
target, can set limits as strong as
$|U_{\alpha\mathcal{N}}|^{2}<\mathcal{O}(10^{-11})$ [32].

There are several models in the literature that can suppress the Dirac mass in
comparison with the magnetic moment of neutrinos [33, 34, 35, 36, 37, 38, 39],
but only to a certain extent and not without fine tuning. We proceed assuming
that $m_{D}$ is sufficiently small so as not to impact the phenomenology, but
note that depending on the amount of fine tuning, bounds on the mixing angles
$|U_{\alpha\mathcal{N}}|$ would also need to be considered.

In performing a fit to MiniBooNE data, we also include a sterile neutrino,
$\nu_{s}$, which mixes with light neutrinos but does not have transition
magnetic moments with active neutrinos. This sterile will be responsible for
short-baseline oscillations with $\Delta m_{41}^{2}$ of $\mathcal{O}(1$
eV${}^{2})$. We provide two examples, the case of a global fit to short-
baseline data excluding MiniBooNE [23] and the case of a joint fit to only
MiniBooNE and the recent MicroBooNE CCQE-like analysis [40]. It is also
important that $\nu_{s}$ and the heaviest neutrino do not mix, as otherwise
$\mathcal{N}_{R}$ would also mix with light neutrinos via $\nu_{s}$. In
summary, our spectrum is defined as

$\nu_{i}\sim\sum_{\alpha=\\{e,\mu,\tau,s\\}}U_{\alpha i}\nu_{\alpha},\text{
for {\hbox{i\leq 4}} and }\nu_{5}\equiv\mathcal{N}$ (4)

where the approximate symbol means that any additional term contains very
small mixing elements. In addition, this mixing would result in a
corresponding dipole coupling $d_{\alpha
s}\overline{\nu}_{\alpha}\sigma_{\mu\nu}F^{\mu\nu}\nu_{s}$, which is strongly
constrained by Big-Bang Nucleosynthesis and stellar cooling. While the former
may be modified à la secret interactions, the latter is significantly more
robust and constrains $d_{\mu s}<2\times 10^{-12}\mu_{B}$ [41].

The relevant interactions for this work are shown in Fig. 1, which include
Primakoff upscattering off of a nuclear target $\nu A\to\mathcal{N}A$ (left)
and the radiative decay of the heavier neutrino
$\mathcal{N}\to\nu_{\alpha}\gamma$ (right). Due to the photon propagator, the
scattering process is dominated by coherent exchange with the nucleus, except
at the largest $m_{\mathcal{N}}$ values. We include helicity-flipping
upscattering, where the helicity of $\mathcal{N}$ is opposite to that of
$\nu_{\alpha}$, and helicity-conserving, where they are the same. The latter
is suppressed by the typical energy of the process,
$m_{\mathcal{N}}^{2}/E^{2}$, and is only relevant at the largest masses. To
give an example, the cross section for an incoming neutrino with energy
$E_{\nu_{\alpha}}=1\;{\rm GeV}$ upscattering off of a carbon nucleus into a
neutrissimo with mass $m_{\mathcal{N}}=50\;{\rm MeV}$ and dipole coupling
$d_{\alpha N}=10^{-6}\;{\rm GeV}^{-1}$ is $\sigma\approx 5\times
10^{-39}\;{\rm cm}^{2}$ [18].

The rest frame decay width of the $\mathcal{N}$ is given by the incoherent sum
over all outgoing flavors,

$\Gamma_{\mathcal{N}\to\nu\gamma}=\left(\sum_{\alpha}|d_{\alpha\mathcal{N}}|^{2}\right)\frac{m_{\mathcal{N}}^{3}}{4\pi}.$
(5)

Throughout this work, we will assume that either $\mathcal{N}$ or light
neutrinos are Dirac particles, such that the differential decay rate is
proportional to $(1\pm\cos{\theta})$ [19]. For the region of dipole parameter
space preferred by MiniBooNE [23], the HNL has a lab frame decay length $L\sim
1-10\;{\rm m}$ for typical MINERvA neutrino energies and $L\leq 1\;{\rm m}$
for typical MiniBooNE neutrino energies. Thus, Primakoff upscattering near or
within each detector is the most relevant production mechanism in this region
of parameter space. For smaller HNL masses, decay lengths become longer and
upscattering in the dirt along the BNB and NuMI beamlines becomes more
important for MiniBooNE and MINERvA, respectively.

##  III Simulating Neutrissimos

Figure 2: Top) schematic representation of HNL production via upscattering and
the subsequent HNL decay within the MiniBooNE (left) and MINERvA (right)
detectors. Detector images have been adapted from Ref. [42] and Ref. [43].
Bottom) example upscattering rates within two of the MINERvA nuclear targets
as simulated using LeptonInjector. The coherent enhancement of the
upscattering cross sections leads to a larger rate in the high-$Z$ components.

In order to describe neutrissimo interactions in the MiniBooNE and MINERvA
detectors, we have developed a custom simulation based in the LeptonInjector
framework [28]. The simulation begins by injecting neutrinos in MiniBooNE and
MINERvA according the BNB and NuMI beam profiles, respectively. For each
neutrino, a flight path is randomly selected within a cone surrounding the
detector. This path is then used to calculate intersections with different
components of the beamline, including bedrock between the target and detector
as well as various detector subsystems. An upscattering location is sampled
along this flight path according to the cross section in each traversed
material. The final state kinematics of the produced HNL are sampled according
to the $d\sigma/dy$ distribution, where
$y\equiv(E_{\nu}-E_{\mathcal{N}})/E_{\nu}$ in the lab frame [41]. See the
Appendix A for more details on the upscattering cross section used in this
analysis. Next, a decay location is sampled along the flight path of the HNL.
The final state kinematics of decay photon are sampled according to
$d\Gamma/d\cos\theta\propto(1-\cos\theta)$. We trace the flight path of the
photon through different materials until it converts to an $e^{+}e^{-}$ pair,
which are assumed to reconstruct as a single electromagnetic shower within the
MiniBooNE and MINERvA detectors.

The above procedure is easily generalizable to other BSM scenarios and
detector configurations. The user needs only to provide (1) the flux of
initial state particles (neutrinos in this case), (2) the total and
differential cross sections for the relevant processes, and (3) the relevant
detector geometry. The flux and cross section can be input either as analytic
expressions or splines; LeptonInjector will interpolate in the latter case.
The detector geometry is set using a simple configuration file in which the
user instantiates any number of volumes. LeptonInjector supports most simple
volumes as well as general extruded polygons–for example, the nuclear targets
in MINERvA shown in Fig. 2, which are subsections of hexagonal prisms).
LeptonInjector will then generate user-specified initial state particles,
sample their interaction locations, and store the kinematics and weights of
the final state particles. The simulation is set up to be as efficient as
possible, such that most generated events create an observable final state
within the fiducial volume and are appropriately down-weighted. In the
analysis presented here, we have tracked final state particles until the
production of an observable $e^{+}e^{-}$ pair. This strategy can be adapted
for other BSM scenarios which involve multiple interactions before the
production of an observable final state.

The simulation performed for the analysis presented here gives a robust
estimation of the $e^{+}e^{-}$ event rate within the fiducial volume of each
experiment. Further, it provides an accurate kinematic description of the
photons which survive fiducialization (namely, the photon energy and angle
with respect to the beamline). This is vital, as it will be shown that the
region in dipole parameter space preferred by the MiniBooNE excess is highly
sensitive to the kinematics of the final state photon. Such an effect has been
appreciated by previous studies in the literature [11, 24], but a
comprehensive fit over the full dipole parameter space was not attempted until
Ref. [23]; the simulation described above allows us to refine the fit from
Ref. [23].

The simulation is also vital in order to properly assess the ability of
MINERvA to constrain the dipole solution to the MiniBooNE excess. By
leveraging our simulation, we will show that the reconstruction efficiency for
radiative HNL decays in MINERvA varies by 3-4 orders of magnitude across the
parameter space. Previous calculations in the literature have assumed a
constant 10% reconstruction efficiency for radiative HNL decays in MINERvA
[44], and have thus incorrectly concluded that MINERvA data rules out the
dipole-portal MiniBooNE solution. Further, they do not simulate the complex
sub-components of the MINERvA detector–this is an important step, as
upscattering in the high-$Z$ nuclear targets just upstream of the fiducial
volume can contribute significantly to the single shower sample for short-
lived HNLs.

Thus, the LeptonInjector-based simulation performed for this study is
essential to evaluate the status of the dipole-coupled neutrissimos as an
explanation for the MiniBooNE excess. As the community turns toward more
exotic BSM explanations of the MiniBooNE excess, it is imperative that these
models are evaluated within the context of realistic detector descriptions.
Our simulation framework is the prefect tool for such a task; though the
version used in this article was developed specifically to study neutrissimo
interactions in MiniBooNE and MINERvA, it can be adapted easily to accommodate
additional BSM scenarios and neutrino detectors.

##  IV Neutrissimos at MiniBooNE

The MiniBooNE detector uses Cherenkov light to detect final state particles
produced in neutrino interactions. As electrons and photons both produce
electromagnetic showers which show up as distorted Cherenkov rings in the
detector, the two particles are indistinguishable in MiniBooNE. Thus, photons
from the dipole model could contribute to the MiniBooNE electron-like excess.
In the case of a nonzero effective dipole coupling, $d_{\mu\mathcal{N}}$, muon
neutrinos from Fermilab’s Booster Neutrino Beam will undergo Primakoff
upscattering into HNL states. This can happen via photon exchange with nuclear
targets both within the dirt between the BNB target and the detector and
within the CH2 detector volume itself.

###  IV.1 Simulation

As described in section III, we use LeptonInjector to simulate the production
and decay of HNLs in MiniBooNE [28]. A schematic depiction of this process is
shown in the top panel of Fig. 2. Muon neutrinos are injected according to the
BNB flux and allowed to upscatter to HNLs along the 541 m baseline between the
BNB target and the detector or within the 818-ton CH2 detector itself. We
describe the MiniBooNE detector as a sphere of CH2 with a total radius of 6.1
m and a fiducial radius of 5 m. The MiniBooNE detector sits within a sphere of
air with a radius of 9 m, meant to represent the detector hall. To simulate
the BNB, muon neutrinos generated in LeptonInjector propagate through 541 m of
dirt before reaching the MiniBooNE detector. The neutrino upscattering and
subsequent HNL decay are simulated according to the proceudre outlined in
section III. We parameterize the post-fiducialization photon detection
efficiency as a linearly decreasing function of the true photon kinetic energy
[23]. We also impose a reconstruction threshold on the photon kinetic energy
of 140 MeV [45]. We also independently smear the reconstructed visible energy
and scattering angle of each photon according to the resolution of each as a
function of true photon kinetic energy. The energy (angular) resolution comes
from a power-law (quadratic) fit to simulated single electromagnetic shower
events in MiniBooNE [46], with typical values of 10% (3∘) for photons from
this model, consistent with figures reported by the MiniBooNE collaboration
[47].

###  IV.2 Analysis Methodology

Figure 3: The 2$\sigma$ and 3$\sigma$ CL preferred regions to explain the
MiniBooNE anomaly in mass-coupling parameter space for a dipole-coupled heavy
neutral lepton. The pink (green) curves correspond to results from fitting the
$E_{\nu}^{\rm QE}$ ($\cos\theta$) distribution. Dipole model fits are
performed after subtracting the oscillation component from a global fit to a
$3+1$ model excluding MiniBooNE data [23]. In general, the energy and angular
distributions prefer different regions of parameter space, though overlap
exists at the $2-3\sigma$ level.

Figure 4: The $E_{\nu}^{\rm QE}$ (top) and $\cos\theta$ distributions at the
example dipole model hypothesis indicated by the black star in Fig. 3. The
darker contribution in each stacked histogram corresponds to the dipole model
prediction, while the lighter contribution corresponds to the oscillation
contribution. The background-subtracted MiniBooNE excess is indicated by the
black data points, with solid and dashed error bars indicating statistical and
statistical+systematic errors, respectively.

The dipole model has previously been studied as a potential solution to the
MiniBooNE anomaly [11, 12, 13, 14, 15, 16, 17, 18, 22, 23, 24]. Specifically,
we expand upon the study performed in Ref. [23], which examined a mixed model
consisting of an eV-scale $\nu_{4}$ facilitating short-baseline
$\nu_{\mu}\to\nu_{e}$ oscillations as well as an MeV-scale HNL decaying to a
photon via the dipole portal mechanism. Ref. [23] found a preference for a
dipole-coupled HNL with $m_{\mathcal{N}}\sim 400\;{\rm MeV}$ and
$d_{\mu\mathcal{N}}\sim 3\times 10^{-7}\;{\rm GeV}^{-1}$.

We make a number of improvements to this analysis. First, we make use of a
more robust nuclear electromagnetic form factor, implementing a data-driven
Fourier-Bessel function parametrization [48, 49, 50] with the data files made
available in Ref. [51]. This should be compared with the simpler dipole
parameterization, which overestimates the differential cross sections at
larger momentum exchange (see Appendix A for more details). This has an impact
on the allowed regions at larger heavy neutrino masses, as the form factor
used in this study drops off much more quickly at larger $Q^{2}$. It also
reduces the contribution from the dipole model at large scattering angles,
thereby making it difficult to explain the back-scattered lepton angular
distribution of the MiniBooNE excess. This effect has been pointed out in
previous studies of the dipole model in MiniBooNE [16]. In this work, we
perform a more detailed analysis of the dipole parameter space to determine
whether solutions exist which can accommodate both the energy and angular
distributions of the excess. The statistical treatment for each distribution
has been improved–we now consider correlated systematic errors in the
reconstructed $E_{\nu}^{\rm QE}$ distribution from the provided covariance
matrix (after constraining with the covariance matrix for MiniBooNE’s
$\nu_{\mu}$ dataset). We also introduce an uncorrelated systematic error of
13% in the $\cos\theta$ distribution, consistent with that in the
$E_{\nu}^{\rm QE}$ distribution.

We perform fits to the excess only in neutrino-mode data, as MiniBooNE has
collected about an order of magnitude more events in this beam configuration
compared to their antineutrino-mode data. We use the simulated photon events
from the above procedure to perform two different spectral analyses across
dipole parameter space: one in the $E_{\nu}^{\rm QE}$ distribution and one in
the $\cos\theta$ distribution. In both cases, we calculate a $\chi^{2}$ test
statistic comparing the dipole model prediction to the remaining excess after
subtracting off the oscillation contribution from the MiniBooNE-less global
fit reported in Ref. [23]. In the $E_{\nu}^{\rm QE}$ fit, we use the electron-
like channel fractional covariance matrix provided by the MiniBooNE
collaboration after constraining with the covariance matrix in the muon-like
channel. No systematic errors are provided for the $\cos\theta$ distribution;
therefore, as mentioned above, we consider an uncorrelated fractional
systematic error of 13% in each bin of the $\cos\theta$ prediction, consistent
with the level in the $E_{\nu}^{\rm QE}$ channel. Confidence regions are drawn
using a $\Delta\chi^{2}$ test statistic, assuming Wilks’ theorem with two
degrees of freedom [52].

###  IV.3 Results

Figure 5: 1$\sigma$ CL preferred and 2$\sigma$ CL allowed regions with regard
to the MiniBooNE anomaly in mass-coupling parameter space for a dipole-coupled
heavy neutral lepton. The pink (green) curves correspond to results from
fitting the $E_{\nu}^{\rm QE}$ ($\cos\theta$) distribution. Dipole model fits
are performed after subtracting the oscillation component from the combined
MiniBooNE + MicroBooNE CCQE-like $3+1$ fit [40]. Mild preference for a dipole-
coupled heavy neutrino with $m_{\mathcal{N}}\lesssim 100\;{\rm MeV}$ is found
at the 1$\sigma$ level.

Figure 6: Similar to Fig. 4. The $E_{\nu}^{\rm QE}$ (top) and $\cos\theta$
(bottom) distributions at the example dipole model hypothesis indicated by the
black star in Fig. 5. The oscillation component comes from the combined
MiniBooNE + MicroBooNE CCQE-like $3+1$ fit [40].

Figure 7: Similar to Fig. 3, but considering the
$\mathcal{N}\to\nu_{\tau}\gamma$ decay channel. On the left, we show the case
with coupling $d_{\tau\mathcal{N}}=(m_{\tau}/m_{\mu})d_{\mu\mathcal{N}}$ and
on the right we consider $d_{\tau\mathcal{N}}=1\times 10^{-5}\;{\rm GeV}^{-1}$
instead. Compared to Fig. 3, preferred regions in the lower HNL mass region
move to lower $d_{\mu\mathcal{N}}$ values.

The result from the fit procedure described above is shown in Figure 3. One
can see that it is difficult to explain the $E_{\nu}^{\rm QE}$ and
$\cos\theta$ distributions through the same dipole-coupled HNL, as the two
distributions prefer different regions of dipole parameter space. The
preferred regions overlap at the $2\sigma$ CL, though some of this overlap
region is in tension with constraints derived from the NOMAD single-photon
analysis [53]. As stated above, the difference between this result and the
result in Ref. [23] is driven mainly by the updated form factor. This reduces
the Primakoff upscattering rate at large scattering angles, requiring the fit
to push to larger $\mathcal{N}$ masses and dipole couplings to explain this
part of the MiniBooNE excess. That being said, there is a region in parameter
space for $d\sim 10^{-6}\;{\rm GeV}^{-1}$ and $m_{\mathcal{N}}\sim 0.5\;{\rm
GeV}$ which is (1) within the $2\sigma$ allowed region for the $E_{\nu}^{\rm
QE}$ distribution, (2) within the $2\sigma$ allowed region for the
$\cos\theta$ distribution, and (3) not ruled out by NOMAD’s single-photon
search [53]. In Fig. 4 we show $E_{\nu}^{\rm QE}$ and $\cos\theta$
distributions for an example parameter point in this region, located at
$d=1.25\times 10^{-6}\;{\rm GeV}^{-1}$ and $m_{\mathcal{N}}=0.47\;{\rm GeV}$.
One can see that this model can describe most of the excess except for the
region at $\cos\theta\lesssim 0$.

This situation might be improved when considering alternative oscillation
scenarios. The above fits assumed the MiniBooNE-less global-fit result, which
found a best-fit solution at $\\{\Delta m^{2}\approx 1.3\;{\rm
eV}^{2},\sin^{2}(2\theta_{\mu e})\approx 6.9\times 10^{-4}\\}$ [23]. We now
consider an alternative sterile neutrino hypothesis: the result from the
recent MiniBooNE + MicroBooNE CCQE-like $3+1$ combined fit performed by the
MiniBooNE collaboration [40].

The MiniBooNE + MicroBooNE CCQE-like $3+1$ combined analysis found a best-fit
solution at $\\{\Delta m^{2}\approx 0.2\;{\rm eV}^{2},\sin^{2}(2\theta_{\mu
e})\approx 0.03\\}$ [40]. This introduces a much larger $\nu_{\mu}\to\nu_{e}$
oscillation component in MiniBooNE. Thus, the dipole model is primarily driven
to explain the lowest energy and most forward-angle portion of the excess. A
mild preference for a dipole-coupled heavy neutral lepton is found at the
$1\sigma$ level. As the dipole model is no longer required to explain the
broad-angle portion of the MiniBooNE excess, the angular fit is able to
accommodate a large range of heavy neutrino masses while the energy fit
prefers lower heavy neutrino masses at $m_{\mathcal{N}}\lesssim 100\;{\rm
MeV}$. The preferred regions in dipole parameter space under this oscillation
hypothesis are shown in Fig. 5. As a benchmark point, we consider a solution
at $d=1.7\times 10^{-7}\;{\rm GeV}^{-1}$ and $m_{\mathcal{N}}=0.08\;{\rm
GeV}$. As shown in Fig. 6, this benchmark point can reasonably describe the
$E_{\nu}^{\rm QE}$ and $\cos\theta$ distributions of the MiniBooNE excess.

One can also consider a nonzero transition magnetic moment coupling between
the $\mathcal{N}$ and the $\nu_{\tau}$ flavor eigenstate. This would open up
the decay channel $\mathcal{N}\to\nu_{\tau}\gamma$, increasing the decay width
by the ratio
$(|d_{\mu\mathcal{N}}|^{2}+|d_{\mu\mathcal{N}}|^{2})/|d_{\mu\mathcal{N}}|^{2}$.
This will have a more pronounced impact on the fit in the lower HNL mass
region of parameter space, as lifetimes in the higher HNL mass region are
sufficiently short such that introducing another decay channel does not
appreciably change the phenomenology.

As discussed in section II, in some UV completions of the dipole model, a
natural scaling given by
$d_{\tau\mathcal{N}}/d_{\mu\mathcal{N}}=m_{\tau}/m_{\mu}$. The resulting
preferred regions in dipole model parameter space under this assumption are
shown on the left panel of Fig. 7. One can see that, compared with Fig. 3,
solutions explaining the $E_{\nu}^{\rm QE}$ distribution have opened at lower
$d_{\mu\mathcal{N}}$ couplings for $m_{\mathcal{N}}\lesssim 100\;{\rm MeV}$.
We also examine the effect of large tau coupling $d_{\tau\mathcal{N}}=1\times
10^{-5}\;{\rm GeV}^{-1}$, which is meant to capture the extent of flexibility
introduced into the dipole model when allowing for nonzero
$d_{\tau\mathcal{N}}$. The preferred regions for this case are shown in the
right panel of Fig. 7. For both cases, we consider an oscillation contribution
given by the MiniBooNE-less global fit.

##  V Neutrissimos at MINERvA

Neutrino upscattering can also occur in the MINERvA detector. We choose to
study MINERvA for two main reasons: i) the NuMI beam is a higher-energy beam
in comparison with the BNB. This is specially true for the medium-energy (ME)
NuMI configuration, where $\langle E_{\nu}\rangle\simeq 7$ GeV, but it is
still the case for the low-energy (LE) configuration, where $\langle
E_{\nu}\rangle\simeq 3$ GeV. This allows us to probe HNLs of larger masses.
ii) it is one of the few accelerator experiments in the few GeV region to have
a dedicated neutrino-electron ($\nu-e$) scattering analysis. While
measurements of this channel have been performed with greater precision at
experiments like CHARM [54] and CHARM-II [55], LSND [56], reactors [57],
Borexino [58, 59], and Super-Kamiokande [60], they are not as well suited for
the study of the MiniBooNE explanations considered here, where HNLs have
hundreds of MeV in mass. With the exception of CHARM and CHARM-II, the
previous experiments operate at energies below the HNL production threshold
and are therefore not sensitive to our region of interest. While we could also
consider CHARM and CHARM-II, we note that they observe larger neutrino-induced
backgrounds thanks to the faster growth of the SM cross section with respect
to the dipole one. In addition, the HNLs would be produced with a larger boost
factor, and therefore escape more often.

Figure 8:  Top) The predicted $E_{\rm sh}\theta_{\rm sh}^{2}$ distribution
before detector smearing and signal selection for three choices of model
parameters at MINERvA. Bottom) The signal selection efficiency of our analysis
cuts, excluding the $dE/dx$ cut, as a function of the HNL mass.

Figure 9: The $dE/dX$ distribution of selected events for the three MINERvA
analyses. From top to bottom: LE FHC, ME FHC, and ME RHC. The latter has the
largest sensitivity due to the smaller backgrounds. All distributions are
shown post-MINERvA tune, except for ME RHC, where it is shown before tuning.

In this work, we will consider three existing measurements of the neutrino-
electron cross section by MINERvA [25, 26, 27]. The first was performed in the
LE configuration of the NuMI beam operating with a forward-horn current (FHC),
optimizing the number of neutrinos produced. The last two were performed in
the ME configuration, one in FHC and the other in reverse-horn current (RHC)
mode, the latter optimizing the number of antineutrinos. The ME RHC
measurement, also the most recent, is particularly sensitive due to the
smaller antineutrino- and neutrino-induced backgrounds. Unlike the dipole
cross section, antineutrino-nucleus weak cross sections are smaller than
neutrino-nucleus cross sections.

###  V.1 Simulation

Target |  $z$-location (cm) |  $z$-extent (cm) | Mass (kg)
---|---|---|---
1-Fe | 13.6 | 2.567 | 370
1-Pb | 13.6 | 2.578 | 317
2-Fe | 31.3 | 2.563 | 370
2-Pb | 31.3 | 2.581 | 317
3-Fe | 53.4 | 2.573 | 197
3-Pb | 53.4 | 2.563 | 141
3-C | 53.4 | 7.620 | 194
Water | 89.5 | 18.06 | 530
4-Pb | 125.6 | 0.795 | 263
5-Fe | 138.9 | 1.289 | 186
5-Pb | 138.9 | 1.317 | 162
Table 1: Specifications of the MINERvA nuclear targets as implemented in
LeptonInjector [28] for this analysis. Each nuclear target is defined as a
hexagonal prism with an apothem of 92 cm. Z positions and extents of each
nuclear target have been taken from Table 4 of Ref. [43]. The coordinate
system is defined such that $z=0$ corresponds to the front of the MINERvA
detector. We have confirmed that the fiducial mass (bounded by an 85 cm
apothem hexagon) of each nuclear target subcomponent matches the fiducial mass
quoted in Table 4 of Ref. [43]. The last column of this table refers to the
mass of each nuclear target subcomponent within the 92 cm apothem hexagonal
prism.

Figure 10: Left) MINERvA constraints in the dipole parameter space at $95\%$
C.L. Solid lines show our nominal limits assuming a $\eta_{\rm bkg}=30\%$
Gaussian systematic uncertainty on the background normalization, and dashed
ones show the constraints assuming an inflated uncertainty of $\eta_{\rm
bkg}=100\%$ on the background. Regions of preference to explain MiniBooNE in
the minimal dipole model are also shown as filled contours at $95\%$ C.L.
Right) Contours of constant $N_{\mathcal{N}}^{\rm fid}$, the total number of
new-physics photons that convert inside the fiducial volume, overlaid on top
of the same parameter space.

As described in section III, we use LeptonInjector [28] to simulate $\nu
A\to\mathcal{N}A$ upscattering inside as well as outside the MINERvA detector.
A schematic depiction of this process is shown in the top panel of Fig. 2. We
include upscattering in the upstream dirt, in the surrounding air, in the
nuclear target planes (detailed below), in the plastic scintillator, as well
as in the outermost electromagnetic calorimeter. While most of these
components are not part of the fiducial volume for the $\nu-e$ analysis, they
can significantly contribute to the signal rate due to the displaced decays of
the HNLs. For long-lived HNLs, upscattering in the dirt dominates the signal
rate, followed by the nuclear target planes, which contain high-density
materials like $\isotope[208]{Pb}$ and $\isotope[56]{Fe}$. Detailed modeling
of the detector geometry and material composition is necessary to correctly
predict the contributions from these different upscattering sites. We define
each nuclear target to be a hexagonal prism with apothem 92 cm, matching that
of the tracker region [61], and z extent given by Table 4 of Ref. [43]. The
six nuclear target planes are detailed in Table 1. In the bottom panel of Fig.
2, we show the positional distribution of the upscattering rate within two of
the MINERvA nuclear targets as simulated using LeptonInjector [28]. We also
consider upscattering within the electromagnetic calorimeter, defined as a
hexagonal prism surrounding the inner detector with 107 cm apothem, and within
the steel veto shield $\sim 1$ m in front of MINERvA [61]. The fiducial volume
of MINERvA for the $\nu-e$ analyses is assumed to be approximately the same
for both the LE and ME analyses and is defined as a hexagon of $81$ cm apothem
with $\sim 2.8$ m z extent inside the plastic-scintillator.

After production, we track the HNL’s path through the detector and force a
decay to occur before the end of the fiducial volume; we then down-weight the
event by the probability of decaying within the considered region. Each decay
produces a photon for which we physically sample a pair-production location.
Events that do not pair-produce within the fiducial volume are removed. This
procedure accounts for events where the HNL decays outside the fiducial
volume, but the photon conversion happens inside of it. This effect is
important for short-lived HNLs since the rate of HNLs produced in the high-
density lead planes can significantly contribute to the signal rate even
though they are not contained in the fiducial volume.

The neutrino fluxes for the LE mode have been taken from Ref. [62] and for the
ME they have been digitized from Ref. [63]. The total exposures used in the
three $\nu-e$ analyses are $3.43\times 10^{20}$ POT for LE-FHC, $1.16\times
10^{21}$ POT for ME-FHC, and $1.22\times 10^{21}$ POT for ME-RHC.

Figure 11: Same as Figure 10 but for the case where HNLs have a larger tau-
neutrino dipole following an approximate scaling of UV completions of the
operator in Equation 1, $d_{\tau\mathcal{N}}=m_{\tau}/m_{\mu}\times
d_{\mu\mathcal{N}}$. The HNL is shorter-lived due to the additional decay
$\mathcal{N}\to\nu_{\tau}\gamma$.

To reduce neutrino-induced backgrounds, MINERvA applies an extensive list of
selection cuts. To properly estimate the resulting efficiency of these cuts in
our HNL signal, it is important to correctly model the reconstruction of the
energy and angle of the single photons. In the absence of a full detector
simulation, we proceed to approximate the detector energy and angular
resolutions as Gaussian functions. For the energy resolution, we take
$\sigma_{E}/E=5.9\%/(\sqrt{E/{\rm GeV}})+3.4\%$ [25], while for the angular
resolution, we take an energy-independent angular resolution of
$\sigma_{\theta}=0.7^{\circ}$, assumed to be isotropic in the shower’s
azimuthal angle 111This is only an approximation, as the MINERvA detector is
not azimuthally symmetric. Nevertheless, the differences in resolution in the
X and Y planes are small [61, 64] and neglected here.. The angular resolution
is implemented by sampling a polar angle $\delta\theta$ from a Gaussian
distribution of standard deviation $\sigma_{\theta}$, rotating the photon by
$\delta\theta$ with respect to its momentum, assigning it an azimuthal angle
$\varphi$ from the uniform distribution $[0,2\pi]$, and finally rotating the
photon back to the laboratory frame by its original polar angle $\theta_{\rm
true}$.

###  V.2 Event selection

Now we discuss the most important signal selection cuts. The analyses [25, 26]
make use of a long list of signal selection cuts, designed to suppress as many
neutrino-nucleus scattering backgrounds as possible. The most worrisome
backgrounds include $\pi^{0}$ production and $\nu_{e}$CC scattering. The
former is particularly important for our radiative decay signal, as it can
give rise to coherent single-photon-like signatures. Cuts related to the
shower radius and transversal as well as longitudinal profiles are not
implemented in our analysis but are expected to have large acceptance due to
our signal being a true single photon (as opposed to two photons from
$\pi^{0}$ or from the $e^{+}e^{-}$ pairs considered in the new physics model
of Ref. [65]). The series of cuts are illustrated in Figure 8, where we show
the acceptance of the cuts as a function of the HNL mass. The selection
acceptance is largely independent of the dipole coupling.

We start with the cut on the reconstructed shower energy, $E_{\rm sh}>0.8$
GeV. To suppress $\nu_{e}$CC backgrounds, a cut on the reconstructed momentum
exchange under the hypothesis of neutrino-nucleon quasi-elastic scattering is
also implemented. It is defined a s $Q^{2}_{\rm QE-reco}=2m_{n}(E_{\nu}^{\rm
QE-reco}-E_{\rm sh})$, with

$E_{\nu}^{\rm QE-reco}=\frac{m_{n}E_{\rm sh}-m_{e}^{2}/2}{m_{n}-E_{\rm
sh}+p_{e}\cos{\theta_{e}}},$ (6)

and the analysis requires $Q^{2}_{\rm QE-reco}<0.02$ GeV2. We note that for
$\nu-e$ scattering this cut can be understood by the following relation,

$Q^{2}_{\rm QE-reco}\simeq(E_{\rm sh}\theta_{\rm
sh})^{2}<\frac{m_{e}E_{\nu}}{2},$ (7)

where we dropped higher-order terms in electron mass and $\theta_{\rm sh}$.
For neutrino-nucleus upscattering in the forward direction (small
$\theta_{\mathcal{N}}$), two-body kinematics for an infinitely heavy nucleus
gives $Q^{2}_{\rm QE-
reco}\simeq(E_{\mathcal{N}}\theta_{\mathcal{N}})^{2}<(E_{\nu}m_{A}-m_{\mathcal{N}}^{2})/2$,
which is a much looser constraint. In addition, the decay of $\mathcal{N}$
introduces even more spread in the angular distribution, so we can already
expect the cut on $Q^{2}_{\rm QE-reco}$ to be very important.

Figure 12: Same as Figure 10 but for the case where HNLs have a larger and
fixed tau-neutrino dipole, $d_{\tau\mathcal{N}}=(100\text{ TeV})^{-1}$. The
HNL is shorter-lived due to the additional decay
$\mathcal{N}\to\nu_{\tau}\gamma$.

The most stringent cut in the analysis; however, is in $E_{\rm sh}\theta_{\rm
sh}^{2}$, required it to be $<3.2$ MeV rad2. The acceptance of this selection
varies from under $10^{-4}$ at high masses to approximately $30\%$ at the
lowest masses. It is larger for Majorana than Dirac HNLs, as Dirac HNLs are
more likely to produce backward-going photons. Finally, the last cut we
implement is the cut on the mean $dE/dx$ of first four scintillator planes,
requiring $dE/dx>4.5$ MeV$/1.7$ cm. Figure 9 shows the $dE/dx$ distributions
of the SM background and HNL decay events, as well as a dashed line
representing the $dE/dx$ cut. The acceptance of this cut is the least-
understood aspect of our analysis as we do not have access to a full detector
simulation. The shape of our $dE/dx>4.5$ MeV$/1.7$ cm distribution is assumed
to be identical to that of the coherent and diffractive $\pi^{0}$ backgrounds
shown in Figure 9. In this approximation, we find this cut has an acceptance
of $93\%$ and $91\%$ for FHC and RHC modes, respectively.

###  V.3 Results

The final $95\%$ C.L. constraints on the dipole model are shown in Figure 10.
The MiniBooNE regions of preference are also shown for comparison. At the
lowest values of $m_{\mathcal{N}}$ where the HNLs are long lived, our
constraints are less sensitive to the MiniBooNE best-fit region than at higher
masses due to the larger HNL boost factors at MINERvA as well as the smaller
fiducial volume when compared to MiniBooNE. For lifetimes longer than
$c\tau^{0}\sim 100$ cm, the event rate in both experiments is dominated by
dirt upscattering. For shorter lifetimes, the event rate is dominated by
upscattering within the detector itself.

One can see that the constraints from MINERvA begin to rule out disfavor
regions of parameter space preferred by MiniBooNE. However, the strongest
MINERvA $2\sigma$ CL limits presented, which come from the ME RHC measurement
and assume 30% uncertainty on the background normalization, do not rule out
the intersection of the $2\sigma$ CL preferred regions from the MiniBooNE
$E_{\nu}^{\rm QE}$ and $\cos\theta$ distributions. This is because the
$\mathcal{N}\to\nu\gamma$ acceptance in the MINERvA ES analysis decreases
rapidly for larger HNL masses, as shown in Figure 8. In the right panel of
Figure 10, we show contours of constant event rate from dipole-coupled HNL
decays in MINERvA. A dedicated single-shower analysis improving the acceptance
for larger HNL masses would likely be sensitive to the entire region of
parameter space preferred by MiniBooNE.

In the left panel of Figure 10, we also show conservative constraints on this
model assuming 100% uncertainty on the background normalization. This is meant
to address the large scale factors (up to factors of $\sim 2$) which have been
applied to the high $dE/dx$ backgrounds in the official MINERvA analysis [27].
These scale factors come from a tuning procedure in kinematic sideband
regions, a process that could potentially wash out any signal from neutrissimo
decays. An optimal analysis would perform a joint fit to both neutrissimo
decays and SM high $dE/dx$ backgrounds to derive constraints (and potentially
allowed regions) on the neutrissimo model presented here; however, such an
analysis is out of the scope of this paper. We also note that the MINERvA
analysis does not include single-photon backgrounds such as radiative
$\Delta(1232)$ decays and coherent single photons. These components are
expected to be small in the energy region of $E_{\rm sh}>800$ MeV [66, 67],
but their inclusion can only make our limits stronger.

We also point out Ref. [68], in which the MINERvA collaboration investigated
an excess in the high $dE/dx$ sideband region of a $\nu_{e}$ charged-current
quasielastic scattering sample. Using topological variables related to the
shower structure, MINERvA concluded that the excess looked more like coherent
or diffractive $\pi^{0}$ production than single photons. This might suggest
that the scale factors in Ref. [27] could also be attributed to additional
$\pi^{0}$ events. However, the analysis presented here suggests that MINERvA
may have unique sensitivity to a neutrissimo-based explanation of the
MiniBooNE excess, thus motivating a more careful separation of one and two
photon events in the high $dE/dx$ region of the MINERvA elastic scattering
samples.

We also derive constraints considering nonzero $d_{\tau\mathcal{N}}$ in Figure
11 and Figure 12. As expected, this impacts the constraints most significantly
at lower HNL masses. For both
$d_{\tau\mathcal{N}}=(m_{\tau}/m_{\mu})d_{\mu\mathcal{N}}$ and
$d_{\tau\mathcal{N}}=(100\;{\rm TeV})^{-1}$ The MINERvA constraints rule out a
large chunk in the middle of the region preferred by the MiniBooNE
$E_{\nu}^{\rm QE}$ distribution. The constraints do not change for
$m_{\mathcal{N}}\gtrsim 200\;{\rm MeV}$, thus the $2\sigma$ overlap between
between the MiniBooNE $E_{\nu}^{\rm QE}$ and $\cos\theta$ distributions
remains valid.

##  VI Discussion

A number of other existing and planned neutrino experiments are sensitive to
an MeV-scale dipole-coupled HNL. Super-Kamiokande can look for single photon
decays from atmospheric neutrinos which upscatter into HNLs within the Earth
[69]. Similarly, one can look for single-photon decays from neutrinos that
upscatter into HNLs within the Earth and propagate to a large-scale
terrestrial detector such as Super-Kamiokande or Borexino [70]. Constraints
from these searches are sensitive to longer-lived HNLs, with typical masses
$\lesssim 10\;{\rm MeV}$ ($\lesssim 100\;{\rm MeV}$) in the solar
(atmospheric) case. These constraints no longer apply for the two cases of
nonzero $d_{\tau\mathcal{N}}$ which we consider in this work, as the HNL
lifetime will be too short to reach the detector. Observed neutrino
interactions from Supernova 1987A can also be used to derive constraints on
the dipole model, as significant upscattering would enhance the stellar
cooling rate, decreasing the neutrino flux observed on Earth [18]. These
constraints require the HNL to be sufficiently long-lived that it can escape
the stellar environment; thus, bounds from Supernova 1987A also do not apply
when we consider nonzero $d_{\tau\mathcal{N}}$.

We also show constraints derived from the NOMAD search for neutrino-induced
single photons, recast as bounds in this parameter space in Ref. [53]. We
rescale them, however, by an overall factor of $\sqrt{2}$ so as to reflect the
decay rate we derived in Equation 5. In addition, we note that these limits
have been obtained with a much less sophisticated simulation than the ones
performed here and that the signal was derived using only events in the pre-
shower detector of NOMAD. In Ref. [53], it is suggested that stronger limits
could be obtained by considering upscattering locations beyond the pre-shower
detector and extended detector volumes where the HNL decay could take place.

We also include the limit imposed by CHARM-II as derived in Ref. [71];
however, we note that it was obtained with a simplified procedure. The
experimental precision on the total neutrino-electron scattering cross section
was used to set limits on the total neutrino-electron upscattering ($\nu
e\to\mathcal{N}e$) cross sections. This is a reasonable assumption at low
values of $m_{\mathcal{N}}$, but potentially breaks down at values close to
the threshold due to differences in kinematics. A similar constraint can be
set using the LSND elastic scattering measurement [18], though it is, in
general, less sensitive than the CHARM-II measurement. A robust re-evaluation
of dipole model constraints from these electron scattering measurements is out
of the scope of this paper. In addition to scattering on electrons, CHARM-II
can provide new limits in the region of interest by considering coherent
neutrino-nucleus upscattering followed by HNL decays into single photons [65].
The sideband with large $dE/dX$ and large values of $E\theta^{2}$ can be used
to set limits, as proposed in [65], although we do not expect them to be as
sensitive due to larger backgrounds and larger boosts.

One can also derive constraints in the dipole-coupled HNLs from LEP through
the $e^{+}e^{-}\to\mathcal{N}\nu_{\ell}$, which can proceed through either the
$\gamma$ or $Z$ mediators [72]. However, these constraints require a strong
enhancement of the mixing between the HNL and SM neutrino; as we consider such
a mixing to be negligible in this model, we do note include constraints from
LEP in our results.

We now discuss the potential for future constraints on dipole-coupled HNLs
from planned measurements. Just like the MINERvA constraint derived in this
work, a neutrino elastic scattering measurement from the NO$\nu$A experiment
would be sensitive to the dipole model [73]. This is especially true of the
DUNE experiment, which has the potential to make a high-statistics neutrino-
electron scattering measurement [74]. This would be particularly advantageous
for the THEIA@DUNE configuration [75] due to its low threshold and large
volume. Dedicated searches at neutrino experiments can further improve
sensitivity to this model. As discussed above, a MINERvA single-shower
analysis without a stringent $E\theta^{2}$ cut could set much stronger
constraints. Experiments which measure CE$\nu$NS, such as COHERENT, NUCLEUS,
and Coherent CAPTAIN-Mills, would also be sensitive to the dipole model by
looking for the coincidence of nuclear recoil from Primakoff upscattering and
a single photon from the HNL decay [76]. These experiments would be most
sensitive to lower mass HNLs with $m_{\mathcal{N}}\lesssim 10\;{\rm MeV}$ due
to the lower energy of typical neutrino sources for CE$\nu$NS experiments.
Existing and upcoming short baseline neutrino experiments, including
MicroBooNE and SBND, have the potential to be sensitive to the neutrissimo
model presented in this paper through a dedicated search for single photon
events [18]. Additionally, neutrino telescopes like IceCube and KM3NeT could
perform searches for events with a double-bang topology from the upscattering
and decay of the HNL, reaching sensitivities of $d_{\mu\mathcal{N}}\sim
10^{-7}\;{\rm GeV}^{-1}$ for $m_{\mathcal{N}}\lesssim 1\;{\rm GeV}$ [71].
High-energy astrophysical tau neutrino observatories such as TAMBO might also
be sensitive to HNL decays from $\nu_{\tau}$ upscattering [77]. Projections
for DUNE [78, 79] estimate that, in the absence of backgrounds, a search for
events with a double-bang morphology could reach $d_{\mu\mathcal{N}}$ values
as low as $\mathcal{O}(10^{-8}\text{ GeV}^{-1})$. Finally, nuclear emulsion
and liquid argon detectors at a future LHC Forward Physics Facility (FPF) will
also be sensitive to transition magnetic moments between HNLs and SM neutrinos
[80].

##  VII Conclusion

In this work, we have explored a mixed model comprising an eV-scale sterile
neutrino and an MeV-scale dipole-coupled HNL. The former facilities
oscillations at short baselines, while the latter introduces the interactions
shown in Fig. 1. The dipole-coupled HNL provides an alternative explanation of
the MiniBooNE excess to the eV-scale sterile neutrino. Thus one can remove
MiniBooNE from global 3+1 fits, reducing tension between appearance and
disappearance experiments while retaining an explanation of the LSND anomaly
[23]. We take the result of the MiniBooNE-less 3+1 global fit as the
oscillation contribution to the MiniBooNE excess and attribute the remaining
excess to decays of the dipole-coupled HNL. We find that spectral fits to the
$E_{\nu}^{\rm QE}$ and $\cos\theta$ distributions prefer different regions of
parameter space in general, though solutions exist which are compatible with
both distributions at the $2\sigma$ confidence level.

We have also derived constraints on the dipole-coupled HNL model using a
MINERvA neutrino-electron elastic scattering measurements [61, 26, 27]. We
find that the most sensitive $\nu-e$ scattering constraints are those obtained
with the NuMI medium-energy mode in antineutrino-enhanced beam configuration.
The constraints from antineutrino-mode are especially strong due to a
reduction in backgrounds at high $dE/dx$, where we expect HNL decays to
contribute. As shown in Fig. 10, MINERvA can exclude large regions of
parameter space preferred by MiniBooNE, but it does not fully exclude it.
There are still allowed MiniBooNE regions at the $2\sigma$ confidence level.
The first is at small $m_{\mathcal{N}}$ values, where HNLs are long-lived and
MINERvA’s small fiducial volume and larger energies reduce the sensitivity.
The second is at larger HNL masses, where the stringent $E\theta^{2}$ cuts
reject most new physics events where decay photons tend to have larger
$E\theta^{2}$. We note that a dedicated search at MINERvA using the same
fiducial volume could significantly improve the signal efficiency in this
large-mass region, and would likely have much better sensitivity, and
potentially probe the entire MiniBooNE-preferred region, as shown in the right
panel of Fig. 10. Nevertheless, as it stands, this mixed model of oscillations
and decay is not dead yet.

## Acknowledgements

MHS is supported by NSF grant PHY-1707971. NSF grant PHY-1801996 supported
CAA, JMC, AD, and NWK for this work. Additionally, CAA is supported by the
Faculty of Arts and Sciences of Harvard University and the Alfred P. Sloan
Foundation. NWK is supported by the NSF Graduate Research Fellowship under
Grant No. 1745302. MAU is supported by the Department of Physics at the
University of Cambridge and SV is supported by the STFC. MH was supported by
Perimeter Institute for Theoretical Physics. Research at Perimeter Institute
is supported by the Government of Canada through the Department of Innovation,
Science and Economic Development and by the Province of Ontario through the
Ministry of Research, Innovation, and Science.

##  Appendix A Upscattering cross section

The cross section for $\mathcal{N}$ production by neutrino upscattering,
$\nu_{\alpha}+A\to\mathcal{N}+A$, has been computed several times in the
literature [81, 82, 41]. We note, however, that all expressions we could find
do not take into account the polarization of the outgoing HNL. This effect is
only important when $m_{\mathcal{N}}/E$ becomes appreciably large and we find
it to be a marginal effect in our calculations.

In the massless limit, a beam of left-handed polarized neutrinos will always
upscatter to right-handed polarized HNLs, assuming the process takes place
purely via the transition magnetic moment. This follows from the chiral
structure of the vertex, $\overline{\nu_{L}}\sigma^{\mu\nu}\mathcal{N}_{R}$.
However, in the massive case, spin and helicity are not equivalent, and both
helicity states of $\mathcal{N}$ can be produced. The helicity-flipping
channel, $\nu_{h}\to\mathcal{N}_{-h}$, typically dominates, while the
helicity-conserving case, $\nu_{h}\to\mathcal{N}_{h}$, will be suppressed by
powers of $m_{\mathcal{N}}/E$, vanishing in the massless limit. Here $h=\pm 1$
denotes the particle’s helicity and $E$ is the typical energy scale of the
scattering process.

We have calculated both terms using the DarkNews code [83], and show our
results in Figure 14. We show a comparison of the upscattering cross section
on Carbon-12 for a few choices of HNL masses for both coherent and proton-
elastic scattering regimes. Scattering on neutrons proceeds only via the
neutron magnetic moment and is much smaller. We only include the proton-
elastic contribution for MiniBooNE, where the proton would be invisible. This
is a conservative approach when deriving the MINERvA limits.

We also show the ratio between helicity-flipping and helicity-conserving
upscattering events inside the MiniBooNE detector as a function of
$m_{\mathcal{N}}$ in Figure 15. The helicity-conserving part is a small
correction, except at the very largest HNL masses, where the rate is
significantly smaller due to the large energy threshold for upscattering.

Figure 13: Comparison of the different nuclear form factors for 12C commonly
used in the literature. The dipole form factor (dashed blue) significantly
overestimates the form factor at large values of the momentum exchange $Q$. In
this paper, we use the Fourier-Bessel parametrization (solid green) for nuclei
for which nuclear data is available, otherwise, we implement the Fermi-
Symmetrized (FS) Woods-Saxon parametrization (dotted yellow).

Figure 14: Comparison between helicity-flipping and helicity-conserving cross
sections for coherent neutrino upscattering on Carbon (left) and free protons
(right) for multiple values of the HNL mass $m_{\mathcal{N}}$. Figure 15: The
ratio between helicity-conserving and helicity-flipping upscattering events on
Carbon at MiniBooNE.

Several model-independent nuclear form factor parametrizations can be found in
the literature. One of them is the Fourier-Bessel parametrization, which
models the charge density in the nucleus as a series of Bessel functions with
a radial cut-off of $R$.

$F^{\rm
FB}(Q)=N\times\frac{\sin(QR)}{QR}\sum_{n}\frac{(-1)^{n}a_{n}}{n^{2}\pi^{2}-Q^{2}},$
(8)

where $1/N=\sum_{n}(-1)^{n}a_{n}/n^{2}\pi^{2}$ is a normalization factor,
ensuring $F(0)=1$. The coefficients $a_{n}$ can be obtained from experimental
data, which is available for a series of common nuclei [48, 49, 50]. We make
use of the machine-readable files provided by Ref. [51].

For nuclei where the nuclear data cannot be found, we implement a Fermi-
symmetrized Woods-Saxon form factor,

$F^{\rm FS-WS}(Q)=\frac{3\pi
a}{r_{0}^{2}+\pi^{2}a^{2}}\frac{a\pi\text{cotanh}({\pi
Qa})\sin({Qr_{0}})-r_{0}\cos({Qr_{0}})}{Qr_{0}\sinh{(\pi Qa)}},$ (9)

where $a=0.523$ fm, $r_{0}=1.03\times A^{1/3}$ fm. These form factors
correctly describe the finite nuclear radius and lead to a strong suppression
of coherent scattering for $Q\gtrsim 200$ GeV. They should be contrasted with
the simpler dipole parametrization

$F^{\rm Dip}(Q)=\frac{1}{1+\frac{Q^{2}}{M_{D}^{2}}}.$ (10)

used in Ref. [23], with $M_{D}=1.18+0.83*A^{1/3}$, and the Helmz form factor

$F^{\rm Helmz}(Q)\frac{3|j_{1}(QR)|}{QR}e^{-Q^{2}s^{2}/2},$ (11)

with $a=0.523$ fm, $s=0.9$ fm, and $R=3.9$ fm.

We provide a comparison of the aforementioned nuclear form factors in Figure
13. It is evident that the dipole parametrization overestimates the cross
section at large values of momentum exchange $Q$. The more sophisticated form
factors used in this work produce more forward angular distributions at
MiniBooNE than what was found in the previous study of Ref. [23].

## References

  * Ahmad et al. [2002] Q. R. Ahmad et al. (SNO), Phys. Rev. Lett. 89, 011302 (2002), eprint nucl-ex/0204009.
  * Ashie et al. [2005] Y. Ashie et al. (Super-Kamiokande), Phys. Rev. D 71, 112005 (2005), eprint hep-ex/0501064.
  * de Salas et al. [2021] P. F. de Salas, D. V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. A. Ternes, M. Tórtola, and J. W. F. Valle, JHEP 02, 071 (2021), eprint 2006.11237.
  * Esteban et al. [2020] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, JHEP 09, 178 (2020), eprint 2007.14792.
  * Capozzi et al. [2021] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri, and A. Palazzo, Phys. Rev. D 104, 083031 (2021), eprint 2107.00532.
  * Aguilar-Arevalo et al. [2001] A. A. Aguilar-Arevalo et al. (LSND), Phys. Rev. D 64, 112007 (2001), eprint hep-ex/0104049.
  * Aguilar-Arevalo et al. [2021] A. A. Aguilar-Arevalo et al. (MiniBooNE), Phys. Rev. D 103, 052002 (2021), eprint 2006.16883.
  * Acero et al. [2022] M. A. Acero et al., in _2022 Snowmass Summer Study_ (2022), eprint 2203.07323.
  * Loinaz et al. [2004] W. Loinaz, N. Okamura, S. Rayyan, T. Takeuchi, and L. C. R. Wijewardhana, Phys. Rev. D 70, 113004 (2004), eprint hep-ph/0403306.
  * Adams et al. [2009] T. Adams et al. (NuSOnG), Int. J. Mod. Phys. A 24, 671 (2009), eprint 0803.0354.
  * Gninenko [2009] S. N. Gninenko, Phys. Rev. Lett. 103, 241802 (2009), eprint 0902.3802.
  * Gninenko [2011] S. N. Gninenko, Phys. Rev. D 83, 015015 (2011), eprint 1009.5536.
  * Dib et al. [2011] C. Dib, J. C. Helo, S. Kovalenko, and I. Schmidt, Phys. Rev. D 84, 071301 (2011), eprint 1105.4664.
  * Gninenko [2012] S. N. Gninenko, Phys. Lett. B 710, 86 (2012), eprint 1201.5194.
  * Masip et al. [2013] M. Masip, P. Masjuan, and D. Meloni, JHEP 01, 106 (2013), eprint 1210.1519.
  * Radionov [2013] A. Radionov, Phys. Rev. D 88, 015016 (2013), eprint 1303.4587.
  * Ballett et al. [2017] P. Ballett, S. Pascoli, and M. Ross-Lonergan, JHEP 04, 102 (2017), eprint 1610.08512.
  * Magill et al. [2018] G. Magill, R. Plestid, M. Pospelov, and Y.-D. Tsai, Phys. Rev. D 98, 115015 (2018), eprint 1803.03262.
  * Balantekin et al. [2019] A. B. Balantekin, A. de Gouvêa, and B. Kayser, Phys. Lett. B 789, 488 (2019), eprint 1808.10518.
  * Balaji et al. [2020a] S. Balaji, M. Ramirez-Quezada, and Y.-L. Zhou, JHEP 04, 178 (2020a), eprint 1910.08558.
  * Balaji et al. [2020b] S. Balaji, M. Ramirez-Quezada, and Y.-L. Zhou, JHEP 12, 090 (2020b), eprint 2008.12795.
  * Fischer et al. [2020] O. Fischer, A. Hernández-Cabezudo, and T. Schwetz, Phys. Rev. D 101, 075045 (2020), eprint 1909.09561.
  * Vergani et al. [2021] S. Vergani, N. W. Kamp, A. Diaz, C. A. Argüelles, J. M. Conrad, M. H. Shaevitz, and M. A. Uchida, Phys. Rev. D 104, 095005 (2021), eprint 2105.06470.
  * Alvarez-Ruso and Saul-Sala [2021] L. Alvarez-Ruso and E. Saul-Sala (2021), eprint 2111.02504.
  * Park et al. [2016] J. Park et al. (MINERvA), Phys. Rev. D 93, 112007 (2016), eprint 1512.07699.
  * Valencia et al. [2019] E. Valencia et al. (MINERvA), Phys. Rev. D 100, 092001 (2019), eprint 1906.00111.
  * Zazueta et al. [2023] L. Zazueta et al. (MINERvA), Phys. Rev. D 107, 012001 (2023), eprint 2209.05540.
  * Abbasi et al. [2021] R. Abbasi et al. (IceCube), Comput. Phys. Commun. 266, 108018 (2021), eprint 2012.10449.
  * Fujikawa and Shrock [1980] K. Fujikawa and R. Shrock, Phys. Rev. Lett. 45, 963 (1980).
  * Pal and Wolfenstein [1982] P. B. Pal and L. Wolfenstein, Phys. Rev. D 25, 766 (1982).
  * Shrock [1982] R. E. Shrock, Nucl. Phys. B 206, 359 (1982).
  * Argüelles et al. [2022] C. A. Argüelles, N. Foppiani, and M. Hostert, Phys. Rev. D 105, 095006 (2022), eprint 2109.03831.
  * Voloshin [1988] M. B. Voloshin, Sov. J. Nucl. Phys. 48, 512 (1988).
  * Barbieri and Mohapatra [1989] R. Barbieri and R. N. Mohapatra, Phys. Lett. B 218, 225 (1989).
  * Babu and Mohapatra [1990] K. S. Babu and R. N. Mohapatra, Phys. Rev. Lett. 64, 1705 (1990).
  * Babu and Mohapatra [1989] K. S. Babu and R. N. Mohapatra, Phys. Rev. Lett. 63, 228 (1989).
  * Leurer and Marcus [1990] M. Leurer and N. Marcus, Phys. Lett. B 237, 81 (1990).
  * Lindner et al. [2017] M. Lindner, B. Radovčić, and J. Welter, JHEP 07, 139 (2017), eprint 1706.02555.
  * Babu et al. [2020] K. S. Babu, S. Jana, and M. Lindner, JHEP 10, 040 (2020), eprint 2007.04291.
  * Aguilar-Arevalo et al. [2022] A. A. Aguilar-Arevalo et al. (MiniBooNE) (2022), eprint 2201.01724.
  * Brdar et al. [2021a] V. Brdar, A. Greljo, J. Kopp, and T. Opferkuch, JCAP 01, 039 (2021a), eprint 2007.15563.
  * Aguilar-Arevalo et al. [2009] A. A. Aguilar-Arevalo et al. (MiniBooNE), Nucl. Instrum. Meth. A 599, 28 (2009), eprint 0806.4201.
  * Aliaga et al. [2014] L. Aliaga et al. (MINERvA), Nucl. Instrum. Meth. A 743, 130 (2014), eprint 1305.5199.
  * Brdar et al. [2021b] V. Brdar, O. Fischer, and A. Y. Smirnov, Phys. Rev. D 103, 075008 (2021b), eprint 2007.14411.
  * Aguilar-Arevalo et al. [2012] A. A. Aguilar-Arevalo et al. (MiniBooNE) (2012), eprint 1207.4809.
  * [46] Private communication with MiniBooNE Collaboration.
  * Shaevitz [2008] M. H. Shaevitz (MiniBooNE), J. Phys. Conf. Ser. 120, 052003 (2008).
  * Fricke et al. [1995] G. Fricke, C. Bernhardt, K. Heilig, L. A. Schaller, L. Schellenberg, E. B. Shera, and C. W. de Jager, Atom. Data Nucl. Data Tabl. 60, 177 (1995).
  * De Vries et al. [1987] H. De Vries, C. W. De Jager, and C. De Vries, Atom. Data Nucl. Data Tabl. 36, 495 (1987).
  * De Jager et al. [1974] C. W. De Jager, H. De Vries, and C. De Vries, Atom. Data Nucl. Data Tabl. 14, 479 (1974), [Erratum: Atom.Data Nucl.Data Tabl. 16, 580–580 (1975)].
  * [51] URL {http://discovery.phys.virginia.edu/research/groups/ncd/index.html}.
  * Wilks [1938] S. S. Wilks, The Annals of Mathematical Statistics 9, 60 (1938), URL https://doi.org/10.1214/aoms/1177732360.
  * Gninenko and Krasnikov [1999] S. N. Gninenko and N. V. Krasnikov, Phys. Lett. B 450, 165 (1999), eprint hep-ph/9808370.
  * Dorenbosch et al. [1989] J. Dorenbosch et al. (CHARM), Z. Phys. C 41, 567 (1989), [Erratum: Z.Phys.C 51, 142 (1991)].
  * Vilain et al. [1994] P. Vilain et al. (CHARM-II), Phys. Lett. B 335, 246 (1994).
  * Auerbach et al. [2001] L. B. Auerbach et al. (LSND), Phys. Rev. D 63, 112001 (2001), eprint hep-ex/0101039.
  * Deniz et al. [2010] M. Deniz et al. (TEXONO), Phys. Rev. D 81, 072001 (2010), eprint 0911.1597.
  * Bellini et al. [2011] G. Bellini et al., Phys. Rev. Lett. 107, 141302 (2011), eprint 1104.1816.
  * Agostini et al. [2019] M. Agostini et al. (Borexino), Phys. Rev. D 100, 082004 (2019), eprint 1707.09279.
  * Fukuda et al. [2001] S. Fukuda et al. (Super-Kamiokande), Phys. Rev. Lett. 86, 5651 (2001), eprint hep-ex/0103032.
  * Park [2013] J. Park, Ph.D. thesis, U. Rochester (2013).
  * Aliaga et al. [2016] L. Aliaga et al. (MINERvA), Phys. Rev. D 94, 092005 (2016), [Addendum: Phys.Rev.D 95, 039903 (2017)], eprint 1607.00704.
  * Bashyal [2021] A. Bashyal, Ph.D. thesis, Oregon State U. (2021).
  * Valencia-Rodriguez [2016] E. Valencia-Rodriguez, Ph.D. thesis, Guanajuato U. (2016).
  * Argüelles et al. [2019] C. A. Argüelles, M. Hostert, and Y.-D. Tsai, Phys. Rev. Lett. 123, 261801 (2019), eprint 1812.08768.
  * Rein and Sehgal [1981] D. Rein and L. M. Sehgal, Phys. Lett. B 104, 394 (1981), [Erratum: Phys.Lett.B 106, 513 (1981)].
  * Wang et al. [2015] E. Wang, L. Alvarez-Ruso, and J. Nieves, Phys. Lett. B 740, 16 (2015), eprint 1407.6060.
  * Wolcott et al. [2016] J. Wolcott et al. (MINERvA), Phys. Rev. Lett. 117, 111801 (2016), eprint 1604.01728.
  * Gustafson et al. [2022] R. A. Gustafson, R. Plestid, and I. M. Shoemaker (2022), eprint 2205.02234.
  * Plestid [2021] R. Plestid, Phys. Rev. D 104, 075027 (2021), eprint 2010.04193.
  * Coloma et al. [2017] P. Coloma, P. A. N. Machado, I. Martinez-Soler, and I. M. Shoemaker, Phys. Rev. Lett. 119, 201804 (2017), eprint 1707.08573.
  * Delgado et al. [2022] F. Delgado, L. Duarte, J. Jones-Perez, C. Manrique-Chavil, and S. Peña (2022), eprint 2205.13550.
  * Bian [2017] J. Bian, in _Meeting of the APS Division of Particles and Fields_ (2017), eprint 1710.03428.
  * Marshall et al. [2020] C. M. Marshall, K. S. McFarland, and C. Wilkinson, Phys. Rev. D 101, 032002 (2020), eprint 1910.10996.
  * Askins et al. [2020] M. Askins et al. (Theia), Eur. Phys. J. C 80, 416 (2020), eprint 1911.03501.
  * Bolton et al. [2021] P. D. Bolton, F. F. Deppisch, K. Fridell, J. Harz, C. Hati, and S. Kulkarni (2021), eprint 2110.02233.
  * Romero-Wolf et al. [2020] A. Romero-Wolf et al., in _Latin American Strategy Forum for Research Infrastructure_ (2020), eprint 2002.06475.
  * Atkinson et al. [2022] M. Atkinson, P. Coloma, I. Martinez-Soler, N. Rocco, and I. M. Shoemaker, JHEP 04, 174 (2022), eprint 2105.09357.
  * Schwetz et al. [2020] T. Schwetz, A. Zhou, and J.-Y. Zhu, JHEP 21, 200 (2020), eprint 2105.09699.
  * Ismail et al. [2022] A. Ismail, S. Jana, and R. M. Abraham, Phys. Rev. D 105, 055008 (2022), eprint 2109.05032.
  * Vogel and Engel [1989] P. Vogel and J. Engel, Phys. Rev. D 39, 3378 (1989).
  * Harnik et al. [2012] R. Harnik, J. Kopp, and P. A. N. Machado, JCAP 07, 026 (2012), eprint 1202.6073.
  * Abdullahi, A. and Hoefken, J. and Hostert, M. and Massaro, D. and Pascoli, S. [2019] Abdullahi, A. and Hoefken, J. and Hostert, M. and Massaro, D. and Pascoli, S., _DarkNews_ , https://github.com/mhostert/DarkNews (2019).