# Explaining the MiniBooNE Excess Through a Mixed Model of Oscillation and
Decay

S. Vergani<EMAIL_ADDRESS>University of Cambridge, Cambridge CB3
0HE, United Kingdom  N.W. Kamp<EMAIL_ADDRESS>Dept. of Physics, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA  A. Diaz<EMAIL_ADDRESS>Dept. of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139,
USA  C.A. Argüelles<EMAIL_ADDRESS>Dept. of Physics, Harvard
University, Cambridge, MA 02138, USA  J.M. Conrad<EMAIL_ADDRESS>Dept. of
Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  M.H.
Shaevitz<EMAIL_ADDRESS>Dept. of Physics, Columbia University,
New York, NY, 10027, USA  M.A. Uchida<EMAIL_ADDRESS>University
of Cambridge, Cambridge CB3 0HE, United Kingdom

###### Abstract

The electron-like excess observed by the MiniBooNE experiment is explained
with a model comprising a new low mass state ($\mathcal{O}(1)$ eV)
participating in neutrino oscillations and a new high mass state
($\mathcal{O}(100)$ MeV) that decays to $\nu+\gamma$. Short-baseline
oscillation data sets are used to predict the oscillation parameters. Fitting
the MiniBooNE energy and scattering angle data, there is a narrow joint
allowed region for the decay contribution at 95% CL. The result is a
substantial improvement over the single sterile neutrino oscillation model,
with $\Delta\chi^{2}/dof$ = 19.3/2 for a decay coupling of $2.8\times 10^{-7}$
GeV-1, high mass state of 376 MeV, oscillation mixing angle of $7\times
10^{-4}$ and mass splitting of $1.3$ eV2. This model predicts that no clear
oscillation signature will be observed in the FNAL short baseline program due
to the low signal-level.

## Introduction

For the past 25 years, anomalies have been observed in short-baseline (SBL)
neutrino oscillation experiments. These have been studied under a model called
“3+1” that introduces a new non-interacting, hence “sterile,” state with mass
of $\mathcal{O}(1\,\rm{eV})$, in addition to the three Standard Model (SM)
neutrino states. In this model, $\nu_{\mu}\rightarrow\nu_{e}$ appearance,
$\nu_{e}$ disappearance, and $\nu_{\mu}$ disappearance searches should all
point to neutrino oscillations at $L/E\sim 1$ m/MeV, where $L$ is the distance
a neutrino of energy $E$ travels, with a consistent set of flavor mixing
parameters [1, 2, 3, 4]. However, while individually the data appear to fit
oscillations, global fits find a small probability that all of the relevant
data sets are explained by the same parameters [2, 3], as measured by the
Parameter Goodness of Fit (PGF) test [5, 6]. In particular, appearance data
from MiniBooNE produces large tension between appearance and disappearance in
the 3+1 model. This is because the 3+1 best-fit parameters from the other data
sets yield a poor fit to the lowest energy range of the MiniBooNE anomaly [7].
Therefore, there is significant interest in explanations for MiniBooNE beyond
the 3+1 model; for example, one can consider decays of a sterile neutrinos
into active neutrinos and singlet scalars [8, 9].

The MiniBooNE anomaly is a 4.8$\sigma$ excess of electron-like events observed
in interactions from a predominantly muon neutrino beam in a Cherenkov
detector [10], which cannot distinguish between electromagnetic showers from
electrons and photons. Hence, a favored alternative to the 3+1 model has been
to introduce MeV-scale heavy neutral leptons (HNLs) that decay via
$\mathcal{N}\rightarrow\nu\gamma$ within the detector, where the photon is
then misidentified as an electron [11, 12, 13, 14, 15, 16, 17, 18]; see
Refs.[19, 20, 21, 22, 23, 24, 25, 26, 27] for misidentified di-electron
scenarios. These initial studies of $\mathcal{N}$-decay models describe the
MiniBooNE energy distribution well but omit the 3+1 oscillations predicted
from fits to the other anomalies.

In this work, we explore a combination of the two explanations by fitting the
MiniBooNE energy and angle distributions using a combined model,
3+1+$\mathcal{N}$-decay. The 3+1 oscillation component has been obtained by
fitting SBL data sets other than MiniBooNE appearance. We will show that such
a model explains the data well, identifying a highly limited range for the
four model parameters: the mixing angle, $\sin^{2}2\theta$, and mass
splitting, $\Delta m^{2}$, for the oscillation; and the HNL mass,
$m_{\mathcal{N}}$, and photon coupling, $d$, for the decay.

## Model

The combination of eV-scale and MeV-scale mass states is motivated if the two
are members of a family of $\mathcal{N}_{j}$ where $j=1,2,3$. If the mass
splittings are similar to the quark and charged-lepton sectors, then the
family might also include a keV-scale member [28, 29]. All members may
interact with photons at a weak level through a dipole portal interaction
[17], also known as neutrino magnetic moment [30, 31, 32, 33, 34, 35]. Thus,
the $\mathcal{N}_{1}=\nu_{4}$ can decay, but the lifetime is typically longer
than the age of the Universe [31, 36]. The keV-scale mass state,
$\mathcal{N}_{2}$, would have a lifetime that is too long to be observed
through decay in terrestrial experiments but could explain observed X-ray
lines [28, 37]. Only the $\mathcal{N}\equiv\mathcal{N}_{3}$ would decay on
length scales relevant to SBL experiments. Conversely, only the eV-scale mass
state would have sufficiently small mass splitting with respect to the light
neutrino states [38, 39, 40, 41] to form observable oscillations. In the
3+1+$\mathcal{N}$-decay model, any given SBL experiment may be sensitive to
signatures of 3+1 oscillations only, $\mathcal{N}\rightarrow\nu\gamma$ only,
or both.

Used to Test | References (Flux Type) | Type of Fit
---|---|---
$\bar{\nu}_{e}$ disappearance |  [42, 43, 44, 45, 46] (Reactor) |
$\nu_{e}$ disappearance |  [47, 48, 49] (Source) |
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ appearance |  [50, 51]($\pi/\mu$ DAR) | $\uparrow$
$\nu_{\mu}\rightarrow\nu_{e}$ appearance |  [52]($\pi/\mu$ DIF) | 3+1-only
$\bar{\nu}_{\mu}$ disappearance |  [53, 54, 55, 56] ($\pi/\mu$ DIF) | $\downarrow$
$\nu_{\mu}$ disappearance |  [57, 54, 58, 59] ($\pi/\mu$ DIF) |
$3+1+\mathcal{N}$ |  [10] (MiniBooNE BNB $\nu$) | $\mathcal{N}$
Table 1: Data sets used in this paper. These include reactor, radioactive,
decay-at-rest (DAR) and decay-in-flight (DIF) neutrino sources. For a detailed
overview of the data sets used in each of the fits see Supplemental Table 2.

In our model, the production and decay of $\mathcal{N}_{j}$ is due to a dipole
portal interaction between left-handed neutrinos, photons, and right-handed
HNLs. The $\mathcal{N}_{j}$ are added to the SM Lagrangian using the following
term [15, 17]:

$\begin{split}\mathcal{L}&\supset\mathcal{L}_{SM}+\sum_{j=1}^{3}\big{[}\bar{\mathcal{N}}_{j}(i\not{\partial}-M_{j})\mathcal{N}_{j}\\\
&+\sum_{\alpha\in\\{e,\mu,\tau\\}}(d_{\alpha
j}\bar{\nu}_{i}\sigma_{\mu\nu}F^{\mu\nu}\mathcal{N}_{j}+h.c.)\big{]},\end{split}$
(1)

where the $\nu_{i}$ correspond to the light neutrino mass states and
$F^{\mu\nu}$ is the electromagnetic field strength. The dimension-full
$d_{\alpha j}$ couplings control the strength of the electromagnetic
interactions between neutrino species, namely the strength of process like
$\mathcal{N}_{j}\to\nu_{i}\gamma$. Note that $d_{\alpha j}$ reflects the
effective dipole coupling of $\mathcal{N}_{j}$ to the weak eigenstate
$\nu_{\alpha}$.

This leads to two production mechanisms for $\mathcal{N}_{j}$: coupling to
virtual photons produced in meson decays, such as $\pi^{0}$, and Primakoff
upscattering of active neutrinos to $\mathcal{N}$ as they traverse material.
Feynman diagrams for the two production (left, middle) and decay (right)
processes are shown in Fig. 1.

In our analysis, we considered only production and decay of the third mass
state $\mathcal{N}$. This follows if $d_{\alpha j}$ is the same for all
generations and is found to be sufficiently small, such that upscattering is
rare, because then the small masses of states 1 and 2 lead to lifetimes that
are too long for an SBL experiment to observe decay. In this case, we define
$d\equiv d_{\alpha j}$ as a universal coupling. An alternative explanation if
$d$ is found to be large is that the $d_{\alpha j}$ vary with family member,
suppressing decays of the first and second mass states. In this case, we
define $d\equiv d_{\alpha 3}$, where the coupling of $\mathcal{N}_{3}$ to all
light neutrino species is assumed to be the same. In the decay, the
polarization of $\mathcal{N}$ must be considered [60, 61, 62]. The photon from
a right-handed Dirac $\mathcal{N}$ decay has a $(1-\cos\theta)$ distribution,
where $\theta$ is the angle between the $\mathcal{N}$ and photon momentum
vectors; conversely, a left-handed Dirac $\mathcal{N}$ will decay with a
$(1+\cos\theta)$ distribution for the photons. Production through an
unpolarized virtual, off-shell photon yields an equal combination of right-
handed and left-handed $\mathcal{N}$, leading to an effectively isotropic
photon decay distribution. For the case of upscattering, where the
$\mathcal{N}$ is produced from an interaction with a left-handed muon
neutrino, the outgoing $\mathcal{N}$ will be right-handed and the
$(1-\cos\theta)$ angular distribution of the photons must be considered. All
three new mass states are related via a mixing matrix to the flavor states.
Calling the new sterile flavors $s_{j}$, the mass and flavor states are
related by:

$\nu_{\alpha}=\sum_{i}U_{\alpha i}\nu_{j}+U_{\alpha
3+i}\mathcal{N}_{i},~{}(\alpha=e,\mu,\tau;s_{1},s_{2},s_{3}),$ (2)

where $U_{\alpha j}$ is the extended $6\times 6$ neutral-lepton mixing matrix
[1].

$\pi_{0}$$\nu_{\alpha}$$\gamma$$\mathcal{N}$

$\nu_{\alpha}$$\mathcal{N}$$\gamma$$\mathrm{N}$$\mathrm{X}$

$\mathcal{N}$$\gamma$$\nu_{\alpha}$

Figure 1: $\mathcal{N}$ production from $\pi^{0}$ Dalitz decay (left) and
$\nu$ upscattering (middle), and decay (right).

## Constraints

Three SBL experiments have relevant limits to $\mathcal{N}$ with mass $>10$
MeV. While not appearing directly in the fit, the viable solution must fall
outside of these limits. NOMAD and CHARM-II were high-energy neutrino
experiments with too small $L/E$ to be sensitive to the 3+1 parameters under
discussion. The NOMAD analysis searched directly for photons from HNL decay
[14]. CHARM-II could not differentiate electrons from photons, and the limit
is derived from comparing $\nu_{\mu}$-electron elastic scattering (ES) data to
the SM prediction [63]. At larger $\mathcal{N}$ masses, which are
kinematically inaccessible in the former process, contributions from
$\nu_{\mu}$-nucleon upscattering are also present [22]. However, a detailed
analysis of this process in CHARM-II has not yet been performed. LSND has also
released $\nu_{\mu}$-electron scattering results in agreement with the SM,
placing a limit on $\mathcal{N}$ [17].

Cosmological observations place constraints on additional neutrino species. In
order to alleviate tension with a light sterile neutrino [64, 65, 66, 67, 68,
69, 70], one can invoke either noncanonical cosmological scenarios [71, 72] or
secret neutrino interactions [73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83]. HNL
interactions similar to those studied in this model also play an important
role in cosmology [84], where they may impact Big Bang nucleosynthesis, relax
cosmological bounds on neutrino masses [85], or explain the Hubble tension
[86].

## Oscillation Global Fit

If the maximum energy of the neutrino source is too small to produce the
heavier $\mathcal{N}$ state, then the SBL experiments can observe only 3+1
oscillations. We refer to this collection of SBL experiments that are not
sensitive to $\mathcal{N}$ decay as “3+1-only” experiments. The references for
relevant “3+1-only” experiments used in this analysis are provided in Table 1
(top), including experiments with anomalies of significance from $2\sigma$ to
$4.8\sigma$ and experiments consistent with $\nu$SM oscillations. The
experiments, individually listed, can also be found in Supplemental Table 2.
We use these experiments to determine the oscillation parameters
$\sin^{2}(2\theta)$ and $\Delta m^{2}$ in a 3+1-only fit. Notably, the
“3+1-only” experiments exclude all MiniBooNE results, as we will then use
MiniBooNE to fit for the $\mathcal{N}$ parameters $m_{\mathcal{N}}$ and $d$.

As shown in Table 1 (bottom), we used MiniBooNE data to fit for the
$\mathcal{N}$ parameters $m_{\mathcal{N}}$ and $d$, given the oscillation
parameters from the 3+1-only fit. MiniBooNE has excesses in three appearance
data subsets: neutrino-mode [87], antineutrino mode [88], and with an off-axis
beam [89]. All three cases are compatible with either 3+1 or HNL explanations.
However, the latter two running modes had low statistics and more limited data
releases, so we restricted our fit to the neutrino-mode sample.

For the 3+1-only fit, mixing between heavy neutrinos and the three lightest
mass states is constrained to be small by terrestrial measurements at
accelerators [90]. Further, oscillations involving the two largest mass states
do not contribute to the explaining the anomalies considered in this
work.Therefore, we have explicitly assumed no mixing between the two largest
mass states and the active states. The only relevant squared-mass-splitting
$\Delta m^{2}$ is between the lightest mostly sterile and the mostly active
states, where the masses of the latter are assumed to be degenerate and
negligible. We concentrated on the $4\times 4$ neutral-lepton-mixing submatrix
that relates the lightest mass states to their flavor states. For $U_{\alpha
k}$, where $\alpha$ is the flavor and $k$ is the mass state, the mixing angles
for the appearance and disappearance oscillation signatures are not
independent: $\sin^{2}2\theta_{ee}=4(1-|U_{e4}|^{2})|U_{e4}|^{2}$ (electron
flavor disappearance); $\sin^{2}2\theta_{\mu\mu}=4(1-|U_{\mu 4}|^{2})|U_{\mu
4}|^{2}$ (muon flavor disappearance); and $\sin^{2}2\theta_{\mu e}=4|U_{\mu
4}|^{2}|U_{e4}|^{2}$ (appearance). This implies that the electron and muon
flavor disappearance signals must be consistent with the
$\nu_{\mu}\rightarrow\nu_{e}$ appearance signal, limiting the range of
$\sin^{2}2\theta_{e\mu}$.

This analysis employed the 3+1 global fitting code described in Ref. [3]. The
list of experiments used in the 3+1-only fit can be found in Supplemental
Table 2. Compared to Ref. [3], we have added new disappearance results from
the STEREO experiment [46] and updated PROSPECT results [45]. In this update,
we have not included the results from NEUTRINO-4 [91], since the collaboration
has not provided an appropriate data release. We have also not included the
latest result from IceCube [92, 93], which shows a preferred region at 90% CL
compatible with our light sterile neutrino best-fit point, since the
collaboration has not provided enough information to reproduce the analysis.
To reiterate, the 3+1-only global fit omitted the MiniBooNE neutrino-mode,
antineutrino-mode, and off-axis appearance data sets.

The best-fit parameters are $\Delta m^{2}=1.32\,{\rm eV}^{2}$ and
$\sin^{2}2\theta_{e\mu}=6.9\times 10^{-4}$. In past 3+1 fits, the tension
between the appearance and disappearance data sets [3], as measured using the
PG test, has been very high, with a probability of $4\times 10^{-6}$
(4.5$\sigma$) that the data are explained by the same parameters. Without
MiniBooNE appearance in the fit, the probability increases to $7\times
10^{-3}$ ($2.5\sigma$). Thus, the tension is, in large part, due to the
MiniBooNE appearance data set, which we hypothesize has the additional
component of $\mathcal{N}$-decay, and, hence, poor agreement with 3+1-only.

## MiniBooNE Analysis

In order to fit MiniBooNE data for $\mathcal{N}$ decay, we wrote a Monte Carlo
simulation for the production and decay of $\mathcal{N}$ in the Booster
Neutrino Beam in neutrino mode. Two processes were included for production:
Dalitz-like $\pi^{0}$ decay and Primakoff upscattering $\nu A\to\mathcal{N}A$.
The latter is by far the dominant $\mathcal{N}$-production mode for $10\,{\rm
MeV}<m_{\mathcal{N}}<1000\,{\rm MeV}$. Therefore, we neglected the $\pi^{0}$
decay contribution throughout this study. For the Primakoff mode, we generated
incident $\nu_{\mu}$ and $\nu_{e}$ events from the MiniBooNE neutrino-mode
flux [94]. We then simulated the upscattering rate on both standard upper-
continental crust nuclei [95] and the MiniBooNE CH2 detector medium, using Eq.
A6 in Ref. [17] to calculate the total interaction rate and momentum transfer.
This process produced a sample of right-hand-polarized $\mathcal{N}$ events,
predominately forward peaked due to the $1/t^{2}$ dependence of the
differential cross section.

Figure 2: Preferred regions to explain the MiniBooNE excess in $E_{\nu}^{QE}$
(pink) and $\cos\theta$ (green) as a function of dipole coupling and
$\mathcal{N}$ mass. The black star indicates
$\\{d,~{}m_{\mathcal{N}}\\}=\\{2.8\times
10^{-7}\;\text{GeV}^{-1},~{}376\;\text{MeV}\\}$, which lies in the joint 95%
CL allowed region for both distributions. Constraints from other experiments
are also shown at the 95% CL.

Simulated $\mathcal{N}$ which enter the MiniBooNE detector were forced to
decay into a photon and a neutrino, taking into account polarization, and
weighted by the decay probability. To incorporate the detector efficiency,
$eff$, we performed a linear fit to the reconstructed gamma-ray efficiency as
a function of true energy [96], $eff=(-0.12\,{\rm GeV^{-1}})*E_{true}+0.29$,
which we used to weight the $\mathcal{N}\to\nu\gamma$ events. The true energy
and angle of the photons were smeared independently according to the
resolution given by the MiniBooNE collaboration. More details on the
simulation can be found in Appendix A.

Ideally one would fit the background-subtracted 2-D distribution of visible
energy, $E_{v}$, vs. scattering angle, $\theta$, of the MiniBooNE events to
the prediction from $3+1+\mathcal{N}$-decay. However, the systematic errors
for this distribution have not been released by MiniBooNE. They are only
available for $E_{\nu}^{QE}$, which is a combination of $E_{v}$ and $\theta$
given by [97]

$E_{\nu}^{QE}=\dfrac{2(M_{n}-B)E_{v}-((M_{n}-B)^{2}+M_{e}^{2}-M_{p}^{2})}{2((M_{n}-B)-E_{v}+\sqrt{(E_{v}^{2}-M_{e}^{2})}\cos\theta)},$
(3)

where $M_{n}$, $M_{p}$, and $M_{e}$ are the neutron, proton and electron
masses, and $B$ is the binding energy of carbon. This formula accurately
describes the neutrino energy in the case of two-body charged-current neutrino
scattering with no final state interactions, assuming the neutrinos enter the
detector along the axis from which $\theta$ is measured. Thus, it is
applicable to the oscillation component of the excess. Though $E_{\nu}^{QE}$
has no physical meaning when applied to the photons from $\mathcal{N}$ decay,
the decay kinematics cause most events to be show up at low $E_{\nu}^{QE}$. We
performed a fit to the MiniBooNE excess in $E_{\nu}^{QE}$ using statistical
and systematic uncertainties. We also performed a separate fit to the
scattering angle distribution, although only statistical uncertainty is
available in this case.

To isolate the decay component, we subtracted from the MiniBooNE excess the
predicted contribution of the oscillation component, which was determined from
the $3+1$-only global fit without MiniBooNE data. The remaining excess was fit
to the model for dipole production, decay, and observation in the detector as
described above. Fig. 2 shows confidence regions for both fits in
$\\{d,m_{\mathcal{N}}\\}$ parameter space, computed assuming Wilks’ theorem
with two degrees of freedom is valid for the test statistic
$\chi^{2}(d,m)-\min_{d,m}(\chi^{2}(d,m))$ [98]. We found a region of parameter
space consistent with both distributions at the 95% CL near $d=3\times
10^{-7}$ GeV-1 and $m_{\mathcal{N}}=400$ MeV.

## Results

Fig. 2 shows that the allowed regions from MiniBooNE fits are substantially
lower in $d$ than the NOMAD or LSND limits. The overlapping solution is also
at substantially higher $m_{\mathcal{N}}$ than kinematically accessible by
LSND. Supernova results [17] set limits in a band from approximately
$d=10^{-7}$ to $10^{-11}$ GeV-1, which is below the solution we found for
MiniBooNE. Thus, our preferred region lies in a window of allowed parameters.

We now consider an example HNL decay contribution for $d=2.8\times 10^{-7}$
GeV-1 and $m_{\mathcal{N}}$ = 376 MeV, indicated by the star in Fig. 2. This
corresponds to the best fit to the $E_{\nu}^{QE}$ distribution within the
joint 95% CL allowed region from the $E_{\nu}^{QE}$ and $\cos\theta$ fits.
Table 2 shows the $\chi^{2}$ values for the $3+1$ and $3+1+\mathcal{N}$-decay
fits to both distributions, indicating significant improvement for the
$3+1+\mathcal{N}$-decay model. The global oscillation fit gives tight
constraints requiring $\Delta m^{2}\approx 1.32$ eV2, but allows values of
$\sin^{2}2\theta_{\mu e}\in[3\times 10^{-4},2\times 10^{-3}]$ at the 90% CL.
The same $\mathcal{N}$ decay fit procedure outlined above has been performed
for each end of the allowed $\sin^{2}2\theta_{\mu e}$ range. In each case we
again examined the $\\{d,m_{\mathcal{N}}\\}$ point that best fits the
$E_{\nu}^{QE}$ distribution within the joint 95% CL region. The $\chi^{2}$
values for these fits are also given in Table 2, indicating a preference for a
smaller oscillation contribution in MiniBooNE in order to explain both the
$\cos\theta$ and $E_{\nu}^{QE}$ distributions via $\mathcal{N}\to\nu\gamma$.
Table 2 also gives the $\chi^{2}$ values for the null case, with neither eV-
scale oscillations nor HNL decay.

Parameters | $\chi^{2}/dof$
---|---
($\sin^{2}2\theta$,$d$,$m_{\mathcal{N}}$) | $3+1+\mathcal{N}$ | $3+1$
| ${E_{\nu}^{QE}}$ | ${\cos\theta}$ | ${E_{\nu}^{QE}}$ | ${\cos\theta}$
(0.30, 3.1, 376) | 5.7/8 | 32.1/18 | 30.5/10 | 86.4/20
(0.69, 2.8, 376) | 7.9/8 | 31.4/18 | 27.3/10 | 71.8/20
(2.00, 5.6, 35) | 20.2/8 | 36.7/18 | 27.6/10 | 40.8/20
(0, 0, 0) | 34.1/10 | 99.4/20 | same | same
Table 2: $\chi^{2}$/dof values for $3+1$ and $3+1+\mathcal{N}$-decay models
obtained by comparing expectations to the MiniBooNE excess in $E_{\nu}^{QE}$
and $\cos\theta$. The parameters in column one refer to ($\sin^{2}2\theta_{\mu
e}\times 10^{-3}$, $d\times 10^{-7}$ [GeV${}^{-1}]$, $m_{\mathcal{N}}$ [MeV]).
The mass splitting is 1.32 eV2 in all cases. The null case (no oscillations
and no HNL decay) is also shown in the last row.

Figure 3:  $E_{\nu}^{QE}$ (left) and $\cos\theta$ (right) distributions of the
MiniBooNE excess for a representative point of the $3+1+\cal{N}$-decay model.
The error bars on the energy distribution include systematic and statistical
errors, while for the angular distribution only statistical errors are
included.

Fig. 3 presents the MiniBooNE excess in $E_{\nu}^{QE}$ (left) and $\cos\theta$
(right) compared with the model prediction. This figure includes both the
global fit oscillation contribution for $\Delta m^{2}=1.3$ eV2 and
$\sin^{2}2\theta_{ee}=6.9\times 10^{-7}$, and HNL decay contribution for
$d=2.8\times 10^{-7}$ GeV-1 and $m_{\mathcal{N}}$ = 376 MeV. We re-emphasize
here that the oscillation contribution shown on these plots comes from a
global fit to the 3+1-only model not including MiniBooNE data. Good agreement
is observed for both distributions, especially noting again the lack of
systematic errors for the angular distribution, which dominate over
statistical errors in the MiniBooNE electron-like analysis [10].

Along with energy and angle, it is essential for the model to be consistent
with the recently published timing distribution of the MiniBooNE excess with
respect to the proton-beam bunch. The excess events occur within $\pm 4$ ns of
the observed $\nu_{\mu}$ events. The time of flight depends on the location at
which the HNL is produced. For the preferred $m_{\mathcal{N}}\sim 400\,$MeV,
the majority of $\mathcal{N}$ events come from upscattering within the
MiniBooNE detector followed by $\mathcal{N}$ decay after travel distances of
$\mathcal{O}(50\,{\rm cm})$. This leads to timing delays of $<2\,{\rm ns}$,
well within the MiniBooNE constraint.

## Conclusion

We have explored a 3+1+$\mathcal{N}$-decay model through fits to the 3+1-only
and MiniBooNE-neutrino-mode data sets. The former yields best-fit oscillation
parameters of $\Delta m^{2}=1.3$ eV2 and $\sin^{2}2\theta_{e\mu}=6.9\times
10^{-4}$. The latter narrows the HNL mass to be within $300-400$ MeV and the
dipole coupling strength to be within $2.5-3.5\times 10^{-7}\,\rm{GeV}^{-1}$.
This model produces a $\Delta\chi^{2}$/dof improvement of 19.3/2 (40.3/2)
compared to the global 3+1 scenario for the fit to the MiniBooNE energy
(angular) distribution. This large improvement in $\Delta\chi^{2}$ motivates a
more detailed analysis by MiniBooNE. Ideally, the experiment would perform a
joint fit to the two-dimensional visible energy and angle distribution, using
a full covariance matrix.

Our model also makes very specific predictions for the experiments now running
in Fermilab’s Short-Baseline Neutrino Program [99]. These experiments make use
of liquid-argon time-projection chambers (LArTPCs) that can separate photon
showers from electron showers with $\sim 85\%$ accuracy [100, 101]. Because
the experiments run in the same beamline and are located within $\sim$50 m of
MiniBooNE, the flux is nearly identical. Thus, the ratio of oscillation to HNL
decay contributions for the far detectors — MicroBooNE and ICARUS — will be
very similar to that of the MiniBooNE case presented here, with $\sim 75\%$ of
the excess events predicted to be single photons. The photon signature will
have large backgrounds even in an LArTPC detector, especially from neutral
current $\Delta$ baryon production with decays to one or two photons plus a
neutron, as well as photons from neutrino interactions produced outside the
detector [102]. However, this background rate will be well constrained from
reconstructed $\Delta$-decay events with a proton. Also, the energy-angle
correlation of the photon in the decay, which depends strongly on
$m_{\mathcal{N}}$, can be used for background rejection since this parameter
is predicted to have a narrow range of values. The oscillation rates for
MicroBooNE and ICARUS are predicted to be low. Within statistics, we predict
that no clear oscillation signature will be observed.

Beyond the SBN program, our model can also produce signatures in $\nu_{\mu}$
ES searches at MINER$\nu$A and NO$\nu$A, as well as dedicated searches for
single photons in T2K [103]. Additionally, at neutrino energies of
$\mathcal{O}(100\,\rm{GeV})$ the long decay length produces “double-bang”
morphologies in IceCube and CCFR [63, 104, 105, 106], which would be a smoking
gun signature for our model.

In summary, we have presented a new physics model including neutrino-partners
with masses of $\mathcal{O}$(1 eV) that participate in oscillations and
$\mathcal{O}$(100 MeV) that decay to single photons. This model can
simultaneously explain the MiniBooNE anomaly and relieve tension in the global
experimental picture for 3+1 oscillations. The results indicate very narrow
ranges of HNL decay and oscillation parameters; thus, this is a highly
predictive result that can be further tested by existing experiments in the
near future.

## 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 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. We thank the MiniBooNE Collaboration for useful input, and
Gabriel Collin and William Louis for comments on the draft of this paper.
Finally, we thank Matheus Hostert and Ryan Plestid for useful discussions.

## References

  * Collin et al. [2016] G. Collin, C. Argüelles, J. Conrad, and M. Shaevitz, Phys. Rev. Lett. 117, 221801 (2016), eprint 1607.00011.
  * Dentler et al. [2018] M. Dentler, A. Hernández-Cabezudo, J. Kopp, P. A. N. Machado, M. Maltoni, I. Martinez-Soler, and T. Schwetz, JHEP 08, 010 (2018), eprint 1803.10661.
  * Diaz et al. [2020] A. Diaz, C. Argüelles, G. Collin, J. Conrad, and M. Shaevitz, Phys. Rept. 884, 1 (2020), eprint 1906.00045.
  * Böser et al. [2020] S. Böser, C. Buck, C. Giunti, J. Lesgourgues, L. Ludhova, S. Mertens, A. Schukraft, and M. Wurm, Prog. Part. Nucl. Phys. 111, 103736 (2020), eprint 1906.01739.
  * Maltoni et al. [2002] M. Maltoni, T. Schwetz, and J. W. F. Valle, Phys. Rev. D 65, 093004 (2002), eprint hep-ph/0112103.
  * Maltoni and Schwetz [2003] M. Maltoni and T. Schwetz, Phys. Rev. D 68, 033020 (2003), eprint hep-ph/0304176.
  * Giunti et al. [2013] C. Giunti, M. Laveder, Y. F. Li, and H. W. Long, Phys. Rev. D 88, 073008 (2013), eprint 1308.5288.
  * Dentler et al. [2020] M. Dentler, I. Esteban, J. Kopp, and P. Machado, Phys. Rev. D 101, 115013 (2020), eprint 1911.01427.
  * de Gouvêa et al. [2020] A. de Gouvêa, O. Peres, S. Prakash, and G. Stenico, JHEP 07, 141 (2020), eprint 1911.01447.
  * Aguilar-Arevalo et al. [2021a] A. A. Aguilar-Arevalo et al. (MiniBooNE), Phys. Rev. D 103, 052002 (2021a), eprint 2006.16883.
  * Gninenko [2009] S. Gninenko, Phys. Rev. Lett. 103, 241802 (2009), eprint 0902.3802.
  * McKeen and Pospelov [2010] D. McKeen and M. Pospelov, Phys. Rev. D 82, 113018 (2010), eprint 1011.3046.
  * 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. 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.
  * 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.
  * Fischer et al. [2020] O. Fischer, A. Hernández-Cabezudo, and T. Schwetz, Phys. Rev. D 101, 075045 (2020), eprint 1909.09561.
  * Bertuzzo et al. [2018a] E. Bertuzzo, S. Jana, P. A. N. Machado, and R. Zukanovich Funchal (2018a), eprint 1808.02500.
  * Bertuzzo et al. [2018b] E. Bertuzzo, S. Jana, P. A. N. Machado, and R. Zukanovich Funchal, Phys. Rev. Lett. 121, 241801 (2018b), eprint 1807.09877.
  * Ballett et al. [2019a] P. Ballett, S. Pascoli, and M. Ross-Lonergan, Phys. Rev. D 99, 071701 (2019a), eprint 1808.02915.
  * Argüelles et al. [2019] C. A. Argüelles, M. Hostert, and Y.-D. Tsai, Phys. Rev. Lett. 123, 261801 (2019), eprint 1812.08768.
  * Ballett et al. [2019b] P. Ballett, M. Hostert, and S. Pascoli, Phys. Rev. D99, 091701 (2019b), eprint 1903.07590.
  * Ballett et al. [2020] P. Ballett, M. Hostert, and S. Pascoli, Phys. Rev. D 101, 115025 (2020), eprint 1903.07589.
  * Abdullahi et al. [2020] A. Abdullahi, M. Hostert, and S. Pascoli (2020), eprint 2007.11813.
  * Abdallah et al. [2020a] W. Abdallah, R. Gandhi, and S. Roy (2020a), eprint 2010.06159.
  * Abdallah et al. [2020b] W. Abdallah, R. Gandhi, and S. Roy, JHEP 12, 188 (2020b), eprint 2006.01948.
  * Drewes et al. [2017] M. Drewes et al., JCAP 01, 025 (2017), eprint 1602.04816.
  * Abada et al. [2020] A. Abada, K. N. Abazajian, A. Abdullahi, S. K. Agarwalla, J. A. Aguilar, W. Altmannshofer, C. A. Argüelles, A. B. Balantekin, G. Barenboim, B. Batell, et al., _Opportunities and signatures of non-minimal heavy neutral leptons_ , Letter of Interest submitted to the Snowmass2021 (2020), https://www.snowmass21.org/docs/files/summaries/NF/SNOWMASS21-NF2_NF3-EF9_EF0-RF4_RF6-CF1_CF0-TF8_TF11_Matheus_Hostert-041.pdf.
  * 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).
  * Dvornikov and Studenikin [2004] M. Dvornikov and A. Studenikin, Phys. Rev. D 69, 073001 (2004), eprint hep-ph/0305206.
  * Giunti and Studenikin [2015] C. Giunti and A. Studenikin, Rev. Mod. Phys. 87, 531 (2015), eprint 1403.6344.
  * Brdar et al. [2021] V. Brdar, A. Greljo, J. Kopp, and T. Opferkuch, JCAP 01, 039 (2021), see arXiv version v3 for updated constraint plot summary., eprint 2007.15563.
  * Nieves [1983] J. F. Nieves, Phys. Rev. D 28, 1664 (1983).
  * Abazajian [2017] K. N. Abazajian, Phys. Rept. 711-712, 1 (2017), eprint 1705.01837.
  * Tanabashi et al. [2018] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
  * Esteban et al. [2019] I. Esteban, M. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni, and T. Schwetz, JHEP 01, 106 (2019), eprint 1811.05487.
  * de Salas et al. [2018] P. de Salas, D. Forero, C. Ternes, M. Tortola, and J. W. F. Valle, Phys. Lett. B 782, 633 (2018), eprint 1708.01186.
  * Capozzi et al. [2016] F. Capozzi, E. Lisi, A. Marrone, D. Montanino, and A. Palazzo, Nucl. Phys. B 908, 218 (2016), eprint 1601.07777.
  * Declais et al. [1995] Y. Declais et al., Nucl. Phys. B434, 503 (1995).
  * Ko et al. [2017] Y. J. Ko et al. (NEOS), Phys. Rev. Lett. 118, 121802 (2017), eprint 1610.05134.
  * Alekseev et al. [2018] I. Alekseev et al. (DANSS), Phys. Lett. B787, 56 (2018), eprint 1804.04046.
  * Andriamirado et al. [2021] M. Andriamirado, A. B. Balantekin, H. R. Band, C. D. Bass, D. E. Bergeron, D. Berish, N. S. Bowden, J. P. Brodsky, C. D. Bryan, T. Classen, et al. (PROSPECT Collaboration), Phys. Rev. D 103, 032001 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.032001.
  * Almazán et al. [2020] H. Almazán, L. Bernard, A. Blanchet, A. Bonhomme, C. Buck, P. del Amo Sanchez, I. El Atmani, J. Haser, F. Kandzia, S. Kox, et al. (STEREO Collaboration), Phys. Rev. D 102, 052002 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.102.052002.
  * Conrad and Shaevitz [2012] J. M. Conrad and M. H. Shaevitz, Phys. Rev. D85, 013017 (2012), eprint 1106.5552.
  * Kaether et al. [2010] F. Kaether, W. Hampel, G. Heusser, J. Kiko, and T. Kirsten, Phys. Lett. B685, 47 (2010), eprint 1001.2731.
  * Abdurashitov et al. [2009] J. N. Abdurashitov et al. (SAGE), Phys. Rev. C80, 015807 (2009), eprint 0901.2200.
  * Aguilar-Arevalo et al. [2001] A. Aguilar-Arevalo et al. (LSND), Phys. Rev. D 64, 112007 (2001), eprint hep-ex/0104049.
  * Armbruster et al. [2002] B. Armbruster et al. (KARMEN), Phys. Rev. D65, 112001 (2002), eprint hep-ex/0203021.
  * Astier et al. [2003] P. Astier et al. (NOMAD), Phys. Lett. B570, 19 (2003), eprint hep-ex/0306037.
  * Cheng et al. [2012] G. Cheng et al. (MiniBooNE, SciBooNE), Phys. Rev. D86, 052009 (2012), eprint 1208.0322.
  * Stockdale et al. [1984] I. E. Stockdale et al., Phys. Rev. Lett. 52, 1384 (1984).
  * Adamson et al. [2012] P. Adamson et al. (MINOS), Phys. Rev. Lett. 108, 191801 (2012), eprint 1202.2772.
  * Adamson et al. [2011] P. Adamson et al. (MINOS), Phys. Rev. D84, 071103 (2011), eprint 1108.1509.
  * Mahn et al. [2012] K. B. M. Mahn et al. (SciBooNE, MiniBooNE), Phys. Rev. D85, 032007 (2012), eprint 1106.5685.
  * Dydak et al. [1984] F. Dydak et al., Phys. Lett. 134B, 281 (1984).
  * Adamson et al. [2016] P. Adamson et al. (MINOS), Phys. Rev. Lett. 117, 151803 (2016), eprint 1607.01176.
  * Formaggio et al. [1998] J. A. Formaggio, J. M. Conrad, M. Shaevitz, A. Vaitaitis, and R. Drucker, Phys. Rev. D 57, 7037 (1998).
  * Baha Balantekin and Kayser [2018] A. Baha Balantekin and B. Kayser, Ann. Rev. Nucl. Part. Sci. 68, 313 (2018), eprint 1805.00922.
  * de Gouvea et al. [2021] A. de Gouvea, P. J. Fox, B. J. Kayser, and K. J. Kelly (2021), eprint 2104.05719.
  * 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.
  * Hannestad et al. [2012] S. Hannestad, I. Tamborra, and T. Tram, JCAP 07, 025 (2012), eprint 1204.5861.
  * Lattanzi and Gerbino [2018] M. Lattanzi and M. Gerbino, Front. in Phys. 5, 70 (2018), eprint 1712.07109.
  * Knee et al. [2019] A. M. Knee, D. Contreras, and D. Scott, JCAP 07, 039 (2019), eprint 1812.02102.
  * Berryman [2019] J. M. Berryman, Phys. Rev. D 100, 023540 (2019), eprint 1905.03254.
  * Gariazzo et al. [2019] S. Gariazzo, P. F. de Salas, and S. Pastor, JCAP 07, 014 (2019), eprint 1905.11290.
  * Hagstotz et al. [2020] S. Hagstotz, P. F. de Salas, S. Gariazzo, M. Gerbino, M. Lattanzi, S. Vagnozzi, K. Freese, and S. Pastor (2020), eprint 2003.02289.
  * Adams et al. [2020a] M. Adams, F. Bezrukov, J. Elvin-Poole, J. J. Evans, P. Guzowski, B. O. Fearraigh, and S. Söldner-Rembold, Eur. Phys. J. C 80, 758 (2020a), eprint 2002.07762.
  * Gelmini et al. [2004] G. Gelmini, S. Palomares-Ruiz, and S. Pascoli, Phys. Rev. Lett. 93, 081302 (2004), eprint astro-ph/0403323.
  * Hamann et al. [2011] J. Hamann, S. Hannestad, G. G. Raffelt, and Y. Y. Y. Wong, JCAP 09, 034 (2011), eprint 1108.4136.
  * Hannestad et al. [2014] S. Hannestad, R. S. Hansen, and T. Tram, Phys. Rev. Lett. 112, 031802 (2014), eprint 1310.5926.
  * Dasgupta and Kopp [2014] B. Dasgupta and J. Kopp, Phys. Rev. Lett. 112, 031803 (2014), eprint 1310.6337.
  * Archidiacono et al. [2015] M. Archidiacono, S. Hannestad, R. S. Hansen, and T. Tram, Phys. Rev. D 91, 065021 (2015), eprint 1404.5915.
  * Archidiacono et al. [2016] M. Archidiacono, S. Hannestad, R. S. Hansen, and T. Tram, Phys. Rev. D 93, 045004 (2016), eprint 1508.02504.
  * Saviano et al. [2014] N. Saviano, O. Pisanti, G. Mangano, and A. Mirizzi, Phys. Rev. D 90, 113009 (2014), eprint 1409.1680.
  * Chu et al. [2015] X. Chu, B. Dasgupta, and J. Kopp, JCAP 10, 011 (2015), eprint 1505.02795.
  * Cherry et al. [2016] J. F. Cherry, A. Friedland, and I. M. Shoemaker (2016), eprint 1605.06506.
  * Chu et al. [2018] X. Chu, B. Dasgupta, M. Dentler, J. Kopp, and N. Saviano, JCAP 11, 049 (2018), eprint 1806.10629.
  * Song et al. [2018] N. Song, M. C. Gonzalez-Garcia, and J. Salvado, JCAP 10, 055 (2018), eprint 1805.08218.
  * Farzan [2019] Y. Farzan, Phys. Lett. B 797, 134911 (2019), eprint 1907.04271.
  * Cline [2020] J. M. Cline, Phys. Lett. B 802, 135182 (2020), eprint 1908.02278.
  * Shakya and Wells [2019] B. Shakya and J. D. Wells, JHEP 02, 174 (2019), eprint 1801.02640.
  * Escudero et al. [2020] M. Escudero, J. Lopez-Pavon, N. Rius, and S. Sandner, JHEP 12, 119 (2020), eprint 2007.04994.
  * Berbig et al. [2020] M. Berbig, S. Jana, and A. Trautner, Phys. Rev. D 102, 115008 (2020), eprint 2004.13039.
  * Aguilar-Arevalo et al. [2007] A. A. Aguilar-Arevalo et al. (MiniBooNE), Phys. Rev. Lett. 98, 231801 (2007), eprint 0704.1500.
  * Aguilar-Arevalo et al. [2013] A. A. Aguilar-Arevalo et al. (MiniBooNE), Phys. Rev. Lett. 110, 161801 (2013), eprint 1303.2588.
  * Adamson et al. [2009] P. Adamson et al. (MiniBooNE, MINOS), Phys. Rev. Lett. 102, 211801 (2009), eprint 0809.2447.
  * Abe et al. [2019a] K. Abe et al. (T2K), Phys. Rev. D 100, 052006 (2019a), eprint 1902.07598.
  * Serebrov et al. [2019] A. P. Serebrov et al. (NEUTRINO-4), Pisma Zh. Eksp. Teor. Fiz. 109, 209 (2019), eprint 1809.10561.
  * Aartsen et al. [2020a] M. Aartsen et al. (IceCube), Phys. Rev. Lett. 125, 141801 (2020a), eprint 2005.12942.
  * Aartsen et al. [2020b] M. Aartsen et al. (IceCube), Phys. Rev. D 102, 052009 (2020b), eprint 2005.12943.
  * Aguilar-Arevalo et al. [2009a] A. A. Aguilar-Arevalo, C. E. Anderson, A. O. Bazarko, S. J. Brice, B. C. Brown, L. Bugel, J. Cao, L. Coney, J. M. Conrad, D. C. Cox, et al., Physical Review D 79 (2009a), ISSN 1550-2368, URL http://dx.doi.org/10.1103/PhysRevD.79.072002.
  * Clarke and Washington [1924] F. W. Clarke and H. S. Washington, Tech. Rep. (1924), URL http://pubs.er.usgs.gov/publication/pp127.
  * Pavlovic et al. [2012] Z. Pavlovic, R. Van de Water, and S. Zeller, Tech. Rep. (2012).
  * Aguilar-Arevalo et al. [2021b] A. A. Aguilar-Arevalo, B. C. Brown, J. M. Conrad, R. Dharmapalan, A. Diaz, Z. Djurcic, D. A. Finley, R. Ford, G. T. Garvey, S. Gollapinni, et al. (MiniBooNE Collaboration), Phys. Rev. D 103, 052002 (2021b), URL https://link.aps.org/doi/10.1103/PhysRevD.103.052002.
  * Wilks [1938] S. S. Wilks, The Annals of Mathematical Statistics 9, 60 (1938), URL https://doi.org/10.1214/aoms/1177732360.
  * Antonello et al. [2015] M. Antonello et al. (MicroBooNE, LAr1-ND, ICARUS-WA104) (2015), eprint 1503.01520.
  * Adams et al. [2020b] C. Adams, M. Alrashed, R. An, J. Anthony, J. Asaadi, A. Ashkenazi, S. Balasubramanian, B. Baller, C. Barnes, G. Barr, et al., Journal of Instrumentation 15, P02007–P02007 (2020b), ISSN 1748-0221, URL http://dx.doi.org/10.1088/1748-0221/15/02/P02007.
  * Acciarri et al. [2017] R. Acciarri, C. Adams, J. Asaadi, B. Baller, T. Bolton, C. Bromberg, F. Cavanna, E. Church, D. Edmunds, A. Ereditato, et al. (ArgoNeuT Collaboration), Phys. Rev. D 95, 072005 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.95.072005.
  * MicroBooNE Collaboration [2020] MicroBooNE Collaboration (2020), URL https://microboone.fnal.gov/wp-content/uploads/MICROBOONE-NOTE-1087-PUB.pdf.
  * Abe et al. [2019b] K. Abe et al. (T2K), J. Phys. G 46, 08LT01 (2019b), eprint 1902.03848.
  * Coloma [2019] P. Coloma, Eur. Phys. J. C 79, 748 (2019), eprint 1906.02106.
  * de Barbaro et al. [1992] P. de Barbaro et al. (CCFR), in _26th International Conference on High-energy Physics_ (1992), URL https://lib-extopc.kek.jp/preprints/PDF/1993/9301/9301432.pdf.
  * Budd et al. [1992] H. S. Budd et al. (CCFR), in _Beyond the Standard Model III (Note change of dates from Jun 8-10)_ (1992), URL https://lib-extopc.kek.jp/preprints/PDF/1992/9212/9212531.pdf.
  * Aguilar-Arevalo et al. [2009b] A. A. Aguilar-Arevalo, C. E. Anderson, A. O. Bazarko, S. J. Brice, B. C. Brown, L. Bugel, J. Cao, L. Coney, J. M. Conrad, D. C. Cox, et al. (MiniBooNE Collaboration), Phys. Rev. D 79, 072002 (2009b), URL https://link.aps.org/doi/10.1103/PhysRevD.79.072002.
  * deNiverville et al. [2017] P. deNiverville, C.-Y. Chen, M. Pospelov, and A. Ritz, Phys. Rev. D 95 (2017), eprint 1609.01770.

Supplemental Material

##  Appendix A Further Explanation of the MiniBooNE $\mathcal{N}$ Simulation

### Primakoff Upscattering Simulation

In order to simulate the Primakoff process $\nu A\to\mathcal{N}A$ for a given
dipole $d$ coupling and HNL mass $m_{\mathcal{N}}$, we first generated a
sample of $5\times 10^{5}$ $\nu_{\mu}$ and $\nu_{e}$ events with energies
according to the MiniBooNE neutrino-mode flux from [94] and initial azimuthal
angles generated randomly in $\cos\theta$ within a cone encompassing the
MiniBooNE detector. We chose a scattering location uniformly in column
density. For the purposes of this study, we have taken the detector to be a
sphere of CH2 with a radius of 6.1 m, surrounded by a concentric sphere of air
with radius 9.1 m, intended to represent the detector hall. The rest of the
volume is taken to be standard upper-continental crust. We next select a
nucleus for the scattering event according to atomic abundances and scattering
cross sections. If the scattering event took place in the dirt, we use atomic
abundances from [95], which are reproduced in Table 1. Upscattering events
inside the MiniBooNE detector happen exclusively off of CH2. Events almost
never happen in the air surrounding MiniBooNE due to the low density.
Scattering cross sections are calculated by numerically integrating the
differential cross section, adapted from Ref. [35]:

$\begin{split}\frac{d\sigma}{dt}=\frac{2\alpha
d^{2}}{m}\bigg{[}&F_{1}^{2}(t)\Big{(}\frac{1}{E_{r}}-\frac{1}{E_{\nu}}+m_{\mathcal{N}}^{2}\frac{E_{r}-2E_{\nu}-M}{4E_{\nu}^{2}E_{r}M}+m_{\mathcal{N}}^{4}\frac{E_{r}-M}{8E_{\nu}^{2}E_{r}^{2}M^{2}}\Big{)}\\\
+&\frac{F_{2}^{2}(t)}{4M^{2}}\Big{(}\frac{2M}{E_{\nu}^{2}}((2E_{\nu}-E_{r})^{2}-2E_{r}M)+m_{\mathcal{N}}^{2}\frac{E_{r}-4E_{\nu}}{E_{\nu}^{2}}+\frac{m_{\mathcal{N}}^{4}}{E_{\nu}^{2}E_{r}}\Big{)}\bigg{]},\end{split}$
(A1)

where $\alpha$ is the fine structure constant, $d$ is the dipole coupling,
$E_{\nu}$ is the SM neutrino energy, $m_{\mathcal{N}}$ is the mass of the
heavy neutrino, $M$ is the target mass, $t=-(p_{\mathcal{N}}-p_{\nu})^{2}$ is
the momentum transfer, $E_{r}=-t/2M$ is the target recoil energy, and
$F_{1/2}(t)$ are the charge/magnetic target form factors, respectively. Note
that the term proportional to $E_{r}m_{N}^{4}$ in the $F_{1}$ line only exists
for spinless nuclei, and must be replaced for nonzero spin nuclei [15]. In the
case of coherent scattering off of a nucleus, $F_{1}$ receives a $Z^{2}$
enhancement and is therefore dominant over $F_{2}$; which has therefore been
neglected for this study. Here $F_{1}$ is given by the dipole approximation

$F_{1}^{A}=\frac{1}{(1+(r_{A}^{2}q^{2})/12)^{2}},$ (A2)

where $r_{A}$ is the nuclear radius of the target. In the case of inelastic
scattering off of nucleons, the form factors are calculated by solving the
following system of equations, repeated here from Appendix A of Ref. [17].

$\begin{split}&G^{\\{p,n\\}}_{E}=F_{1}^{\\{p,n\\}}-\frac{Q^{2}}{4m_{\\{p,n\\}}^{2}}F_{2}^{\\{p,n\\}}=\\{G_{D},0\\},\\\
&G^{\\{p,n\\}}_{E}=F_{1}^{\\{p,n\\}}-\frac{Q^{2}}{4M^{2}}F_{2}^{\\{p,n\\}}=\mu_{\\{p,n\\}}G_{D},\\\
&G_{D}=\frac{1}{(1+Q^{2}/0.71~{}\text{GeV}^{2})^{2}},\\\
&\mu_{\\{p,n\\}}=\\{2.793,-1.913\\}.\end{split}$ (A3)

The total cross section for scattering off of a nuclear target $A$ is given by
the incoherent sum of the nuclear and nucleon scattering cross sections,
namely

$\sigma_{A}=Z^{2}\int_{t_{min}}^{t_{max}}\frac{d\sigma_{\nu
A\to\mathcal{N}A}}{dt}dt+Z\int_{t_{min}}^{t_{max}}\frac{d\sigma_{\nu
p\to\mathcal{N}p}}{dt}dt+(A-Z)\int_{t_{min}}^{t_{max}}\frac{d\sigma_{\nu
n\to\mathcal{N}n}}{dt}dt,$ (A4)

where the lower bounds for each integral are given in Appendix C of Ref. [17]
and the upper bounds are calculated by requiring physical scattering angles
$|\cos(\theta)|<1$. From this we calculate the probability of scattering off a
given nucleus $A_{k}$ (with atomic fractional abundance in dirt/air/CH2
$F_{k}$) as

$P_{k}=\frac{F_{k}\sigma_{A_{k}}}{\sum_{k^{\prime}}F_{k^{\prime}}\sigma_{A_{k^{\prime}}}}.$
(A5)

Nucleus | Z | A | Crust Mass Fraction | Crust Atomic Fraction | Nuclear Radius [MeV-1]
---|---|---|---|---|---
O | 8 | 16 | 0.466 | 0.627 | 0.00218
Si | 14 | 28 | 0.277 | 0.213 | 0.00252
Al | 13 | 27 | 0.081 | 0.065 | 0.00247
Fe | 26 | 56 | 0.05 | 0.019 | 0.00301
Ca | 20 | 40 | 0.037 | 0.02 | 0.00281
K | 19 | 39 | 0.027 | 0.015 | 0.00277
Na | 11 | 23 | 0.026 | 0.024 | 0.00241
Mg | 12 | 24 | 0.015 | 0.013 | 0.00247
Ti | 22 | 48 | 0.004 | 0.002 | 0.0029
P | 15 | 31 | 0.001 | 0.001 | 0.00257
SUPPL. TABLE 1: Relevant parameters of the ten most abundant nuclei in the
Earth’s upper crust according to [95]

At this stage, we decide whether each upscattering event occurred coherently
off of the nucleus or inelastically off of a proton or neutron by considering
the relative cross sections. We then pull a random momentum transfer from Eq.
A1. Once we have chosen a value for $t$, the heavy neutrino energy and
scattering angle are fixed [15]:

$\displaystyle E_{\mathcal{N}}=E_{\nu}-E_{r}$ (A6)
$\displaystyle\cos(\theta)=\frac{E_{\nu}-E_{r}-ME_{r}/E_{\nu}-m_{\mathcal{N}}^{2}/2E_{\nu}}{\sqrt{E_{\nu}^{2}+E_{r}^{2}-2E_{\nu}E_{r}-m_{\mathcal{N}}^{2}}}.$
(A7)

The $1/t^{2}$ dependence of Eq. A1 creates a preference for
$E_{\mathcal{N}}\approx E_{\nu}$ and $\cos(\theta)\approx 1$. At this stage,
if the scattering angle of the $\mathcal{N}$ is greater than the solid angle
of the MiniBooNE detector (considering the scattering location), the event is
rejected. If not, we multiply the existing weights by the probability that the
heavy neutrino decays via $\mathcal{N}\to\nu\gamma$ in MiniBooNE, which is
given by the following expression [17]:

$\displaystyle
P_{decay}=\exp(\frac{-L_{enter}}{L_{decay}})-\exp(\frac{-L_{exit}}{L_{decay}})$
(A8) $\displaystyle L_{decay}=4\pi\hbar c\frac{\beta
E_{\mathcal{N}}}{d^{2}(m_{\mathcal{N}}c)^{4}},$ (A9)

where $L_{enter/exit}$ denote the distance between the creation point of the
HNL and the entry/exit point in the MiniBooNE detector, respectively. The
culmination of this simulation chain is a weighted sample of $\mathcal{N}$
events which come from the Primakoff upscattering process and decay in the
MiniBooNE detector. The last remaining step is to calculate the POT required
to get $N=5\times 10^{5}$ upscattering events along the BNB beamline. This
will depend on the dipole coupling and heavy neutrino mass in general, and can
be calculated using

$N=\text{POT}\int_{E_{\nu}=0~{}\text{GeV}}^{10~{}\text{GeV}}dE_{\nu}\phi(E_{\nu})\sum_{\text{targets}~{}A}\sigma_{A}(E_{\nu})n_{A},$
(A10)

where $\phi(E_{\nu})~{}[\nu/\text{POT}]$ is again the neutrino flux in
MiniBooNE [94], $\sigma_{A}~{}[\text{cm}^{2}]$ is the total upscattering cross
section for nuclear target $A$, and $n_{A}[\text{nuclei}/\text{cm}^{2}]$ is
the column density for $A$ along the beamline.

Figure 1 shows the energy/angular fit $95\%$ CL allowed regions in
$\\{d,m_{\mathcal{N}}\\}$ parameter space for the three $\sin^{2}2\theta$
values considered in Table 2. One can see that the overlap region between the
two fits becomes less strict for larger oscillation contributions. This is
partially because larger oscillation contributions give a poorer fit to the
MiniBooNE excess, as shown by the $\Delta\chi^{2}$ plots in Figure 2. Here one
can see that for the largest $\sin^{2}2\theta$, there is no longer a closed
contour in either the $E_{\nu}^{QE}$ or $\cos\theta$ fit at $3\sigma$ CL.
Figure 3 shows each of the $E_{\nu}^{QE}$ and $\cos\theta$ predictions from
Table 2 compared with the MiniBooNE excess. These plots further indicate a
preference for a smaller $eV$-scale oscillation contribution (considering
global best-fit oscillation parameters) in order to fit MiniBooNE.

Finally, we consider the timing delay distribution for HNL decays in the
MiniBooNE detector. The timing delay is defined as the time between HNL
production and decay minus the same time it would take for a speed of light
neutrino to travel the same distance. Figure 4 shows this delay for the three
different oscillation amplitude cases, indicating timing delays small enough
to be consistent with the MiniBooNE excess timing distribution [10].

SUPPL. FIG. 1:  95% CL allowed regions in HNL dipole model parameter space
from energy and angular fits to the MiniBooNE excess, considering an
oscillation contribution with $\Delta m^{2}=1.3$ eV2 and three different
amplitudes: $\sin^{2}2\theta=6.9\times 10^{-4}$ (top left),
$\sin^{2}2\theta=3\times 10^{-4}$ (top right), and $\sin^{2}2\theta=2\times
10^{-3}$ (bottom).

SUPPL. FIG. 2: Maps of
$\Delta\chi^{2}=\chi^{2}(d,m_{\mathcal{N}})-\min_{\\{d,\mathcal{N}\\}}\chi^{2}(d,m_{\mathcal{N}})$
in HNL dipole model parameter space for the energy (left) and angular (right)
fits to the MiniBooNE excess, considering an oscillation contribution with
$\Delta m^{2}=1.3$ eV2 and three different amplitudes:
$\sin^{2}2\theta=6.9\times 10^{-4}$ (top), $\sin^{2}2\theta=3\times 10^{-4}$
(middle), and $\sin^{2}2\theta=2\times 10^{-3}$ (bottom).

SUPPL. FIG. 3: Example HNL dipole model predictions for the energy (left) and
angular (right) distributions compared with the MiniBooNE excess, considering
an oscillation contribution with $\Delta m^{2}=1.3$ eV2 and three different
amplitudes: $\sin^{2}2\theta=6.9\times 10^{-4}$ (top),
$\sin^{2}2\theta=3\times 10^{-4}$ (middle), and $\sin^{2}2\theta=2\times
10^{-3}$ (bottom). The specific $d$ and $m_{\mathcal{N}}$ values shown here
come from Table 2.

SUPPL. FIG. 4: Timing delay distributions for HNL events in the MiniBooNE
detector, as defined in the text. The specific $d$ and $m_{\mathcal{N}}$
values shown here come from Table 2. Three different $eV$-scale oscillation
amplitudes are considered: $\sin^{2}2\theta=6.9\times 10^{-4}$ (top left),
$\sin^{2}2\theta=3\times 10^{-4}$ (top right), and $\sin^{2}2\theta=2\times
10^{-3}$ (bottom).

### Meson Decay Simulation

The $\pi^{0}$s have been simulated using a Sanford-Wang (SW) parametrization
with $\pi^{0}$ distribution calculated as the average between the $\pi^{+}$
and $\pi^{-}$ distributions with coefficients taken from [107]. A rejection
sampling from previous studies [108] has been used to simulate momenta
$p_{\pi^{0}}\in[0,8]\>$ GeV, emission angles with respect to the direction of
the beam $\theta_{\pi^{0}}\in[0,\frac{\pi}{2}]$, and angle
$\phi_{\pi^{0}}\in[0,2\pi]$ of $\pi^{0}$ produced at the source in lab frame.
Four-momentum vectors have been created with the aforementioned triplets
$p_{\pi^{0}}$, $\theta_{\pi^{0}}$, and $\phi_{\pi^{0}}$ according to
$p^{\mu}_{\pi^{0}}=(E_{\pi^{0}},p_{\pi^{0}}[\cos(\phi_{\pi^{0}})\sin(\theta_{\pi^{0}}),\sin(\phi_{\pi^{0}})\sin(\theta_{\pi^{0}}),\cos(\theta_{\pi^{0}})])$.
A decay probability $P_{decay}$ has been associated to each $\pi^{0}$
according to Eq. A8. In this study, $\mathcal{N}$ have been created via the
three body Dalitz-like decay $\pi^{0}\rightarrow\mathcal{N}+\nu+\gamma$.
Assuming $\nu$ to be massless, the energy $E_{\mathcal{N}}$ in the $\pi^{0}$
rest frame has been constrained to be
$E_{\mathcal{N}}\in(m_{\mathcal{N}},\frac{{m_{\pi^{0}}}^{2}+{m_{\mathcal{N}}}^{2}}{2m_{\pi^{0}}})$
and the differential branching ratio in the $\pi^{0}$ rest frame for each
$\mathcal{N}$ has been evaluated using Eq. A.10 in [17]. The simulated
$\mathcal{N}$ have been subsequently boosted to the lab frame of the
$\pi^{0}$. Considering the channel $\mathcal{N}\rightarrow\gamma+\nu$ with a
massless $\nu$, the $\gamma$ has been simulated assuming their energies
$E_{\gamma}=\frac{m_{\mathcal{N}}}{2}$, emission angle
$cos(\theta_{\gamma})\in[-1,1]$, angle $\phi_{\gamma}\in[0,2\pi]$, and
associated four-momentum vectors to be
$p^{\mu}_{\gamma}=(E_{\gamma},p_{\gamma}[\cos(\phi_{\gamma})\sin(\theta_{\gamma}),\sin(\phi_{\gamma})\sin(\theta_{\gamma}),\cos(\theta_{\gamma})])$.
The resulting $\gamma$s have then been boosted to the lab frame of the
$\mathcal{N}$ and multiplied by a normalization constant
$N_{const}=\frac{N_{\pi^{0}}}{N_{sim}}$. $N_{\pi^{0}}=1.17*10^{21}$ is the
number of $\pi^{0}$ produced in neutrino mode over the lifetime of MiniBooNE
for $18.75*10^{20}$ POT [10] and a $\pi^{0}$ multiplicity per POT 0.9098 [17].
$N_{sim}$ is the effective number of $\pi^{0}$ decays actually simulated.

##  Appendix B Further Discussion on the Fits

Experiment | Ref. | Type | In $3+1+\mathcal{N}$ osc fit? | In $3+1+\mathcal{N}$ decay fit?
---|---|---|---|---
Bugey | [42] | $\bar{\nu}_{e}\rightarrow\bar{\nu}_{e}$ | Yes | No
NEOS | [43] | $\bar{\nu}_{e}\rightarrow\bar{\nu}_{e}$ | Yes | No
DANSS | [44] | $\bar{\nu}_{e}\rightarrow\bar{\nu}_{e}$ | Yes | No
PROSPECT | [45] | $\bar{\nu}_{e}\rightarrow\bar{\nu}_{e}$ | Yes | No
STEREO | [46] | $\bar{\nu}_{e}\rightarrow\bar{\nu}_{e}$ | Yes | No
KARMEN/LSND Cross Section | [47] | $\nu_{e}\rightarrow\nu_{e}$ | Yes | No
Gallium | [48, 49] | $\nu_{e}\rightarrow\nu_{e}$ | Yes | No
SciBooNE/MiniBooNE | [53] | $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{\mu}$ | Yes | No
CCFR | [54] | $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{\mu}$ | Yes | No
MINOS | [55, 56] | $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{\mu}$ | Yes | No
SciBooNE/MiniBooNE | [57] | $\nu_{\mu}\rightarrow\nu_{\mu}$ | Yes | No
CCFR | [54] | $\nu_{\mu}\rightarrow\nu_{\mu}$ | Yes | No
CDHS | [58] | $\nu_{\mu}\rightarrow\nu_{\mu}$ | Yes | No
MINOS | [59] | $\nu_{\mu}\rightarrow\nu_{\mu}$ | Yes | No
LSND | [50] | $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ | Yes | No
KARMEN | [51] | $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ | Yes | No
MiniBooNE (BNB) | [88] | $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ | No | No
NOMAD | [52] | $\nu_{\mu}\rightarrow\nu_{e}$ | Yes | No
MiniBooNE (NuMI) | [89] | $\nu_{\mu}\rightarrow\nu_{e}$ | No | No
MiniBooNE (BNB) | [97] | $\nu_{\mu}\rightarrow\nu_{e}$ | No | Yes
SUPPL. TABLE 2: Table of experimental results used in the fits in this paper.

This Appendix provides further information on the fit results presented in
this paper.

Table 2 provides a table of the experimental results used in the fits, with
references. Notation in Column 3 refers to electron flavor disappearance, muon
flavor disappearance, and muon-to-electon flavor appearance for neutrinos and
antineutrinos. The lastt two columns indicate the data sets used in the fits
to the new $3+1+\mathcal{N}$ model. This model has two contributions: the
oscillation contribution, established through a 3+1 oscillation fit to
experiments in Column 3 and the $\mathcal{N}$ decay fit which is established
through a fit to the MiniBooNE 2020 data set only, as indicated in Column 4.
This omits MiniBooNE (BNB) antineutrino data and MiniBooNE (NuMI) data from
the $\mathcal{N}$ decay fit, which have anomalous signals consistent with
MiniBooNE (BNB) neutrino data, but at low statistics and with high
backgrounds. Including these requires additional development of Monte Carlo
beamline simulations, and the modest results will not change the conclusion of
this paper, thus we conclude this is beyond the scope of this article.

Below, we provide the results of the oscillations-only global fits using the
data indicated in Table 2. A thorough description of the fitting process is
provided in Ref. [3]. In addition to the experiments listed in Ref. [3], we
have updated the MiniBooNE (BNB) neutrino appearance data set [97], updated
the PROSPECT antineutrino disappearance data set [45], and introduced the
STEREO antineutrino disappearance data set.

In reporting the tension, we make use of the probability associated with the
Parameter Goodness of Fit Test, or PGF Test, which is a standard in our field.
In this test [6], the global (glob) set data is divided into disappearance
(dis) and appearance (app) sets. One defines the $\chi^{2}$ and degrees of
freedom as:

$\displaystyle\chi^{2}_{PGF}$ $\displaystyle=$
$\displaystyle\chi^{2}_{glob}-(\chi^{2}_{app}+\chi^{2}_{dis}),$ (A11)
$\displaystyle N_{PGF}$ $\displaystyle=$
$\displaystyle(N_{app}+N_{dis})-N_{glob}~{}.$ (A12)

We quote the associated probability as the tension. Table 3 lists the inputs
to the calculation that appears in this paper.

For the 3+1-only fit used in this article, we run the oscillation fit on the
experiments that would not be sensitive to the $\mathcal{N}$ decay. As stated
earlier, this procedure finds a best fit point of $\Delta m^{2}=1.32$ eV2 and
$\sin^{2}2\theta_{e\mu}=6.9\times 10^{-4}$, with an allowed range
$\sin^{2}2\theta_{\mu e}\in[3\times 10^{-4},2\times 10^{-3}]$ at the 90% CL.
The best-fit regions are shown on the leftmost plot in Fig 5. The tension
within this fit is demonstrated by separately fitting the appearance and
disappearance data sets, which are shown in the middle and right plot in 5,
respectively. This tension is found to have a p-value of $7\times 10^{-3}$,
with the relevant parameters summarized in Table 3.

SUPPL. FIG. 5: These plots show the best fit regions for the 3+1+only
oscillation data sets used in this work. The left plot shows the best fit
regions of the global fits, with 90%, 95%, and 99% regions shown in red,
green, and blue, respectively. The middle plot shows the best fit region when
only fitting to the appearance subset of data, while the right plot shows the
best fit region when only fitting to the disappearance subset. The tension
within this data set can be qualitatively seen by the fact that the best fit
regions do not overlap between the appearance-only and disappearance-only data
sets.

For completeness, we also provide a global fit including all experiments
listed in Table 2. These results, which we’ll label “3+1-complete,” are not
used in the preceding analysis. The best fit regions for this $3+1$ global fit
is shown in the left plot of Fig. 6, with the best fit point at
$\sin^{2}(2\theta_{\mu e})=1\times 10^{-3}$ and $\Delta m_{41}^{2}=1.32$ eV2.
As has been seen before [1, 2, 3, 4], this model suffers from a tension
between the appearance and disappearance data sets. The best-fit regions of
the appearance and disappearance data sets are shown in the middle and right
plot of Fig. 6, respectively. The tension found in the global data set is
found to have a p-value of $8\times 10^{-7}$, with the relevant parameters
summarized in Table 3.

SUPPL. FIG. 6: These plots show the best fit regions for the 3+1+complete
oscillation data sets. This serves as an update to the fits provided in [3],
and are not used in the preceding analysis. The left plot shows the best fit
regions of this complete global fit, with 90%, 95%, and 99% regions shown in
red, green, and blue, respectively. The middle plot shows the best fit region
for the appearance subset of data, while the right plot shows the best fit
region for the disappearance subset. As seen before, the global sterile
neutrino oscillation model suffers from substantial tension.

Our analysis fits to the MiniBooNE (BNB) neutrino energy distribution, which
has been presented with statistical and systematic error. Specifically, the
points and errors are taken from Fig. 19 of Ref. [97]. In this case, the data
are the excess events when the MiniBooNE measurement is compared to the
constrained backgrounds. It is important to note that for most neutrino energy
bins, the systematic error dominates, and so is necessary to include in the
uncertainty when performing the fit. The result of the $3+1+\mathcal{N}$
energy fit has been shown in the paper in Fig. 3, left, and here we reproduce
and enlarge the same figure in Fig. 7. Light pink indicates the oscillation
component, with $\Delta m^{2}=1.3$ eV2 and $\sin^{2}2\theta=6.9\times
10^{-4}$. The darker pink regions indicate the HNL decay component from both
coherent and incoherent upscattering production, with $d=2.8\times 10^{-7}$
GeV-1 and $m_{\mathcal{N}}=376$ MeV.

SUPPL. FIG. 7: Enlarged version of Fig. 3, left, in the paper. The prediction
for the oscillation contribution to the $3+1+\mathcal{N}$ component is shown
in light pink. The prediction for the $\mathcal{N}$ decay component, from both
coherent and incoherent upscattering, is shown in darker pink. See text in
this appendix for explanation of the data.

MiniBooNE has not provided a data release for the angular distribution of the
electromagnetic shower in the excess events. The data shown in Fig. 3, right,
is obtained by subtracting the unconstrained backgrounds from the MiniBooNE
measurement shown in Fig. 8 of Ref. [97]. Only statistical errors are provided
by MiniBooNE at present. We note that one would expect the systematic
uncertainty to be substantially larger, as was the case for the neutrino
energy. Therefore, it would be interesting to repeat this analysis considering
a robust treatment of the systematic errors in the angular distribution. Fig.
8 reproduces and enlarges Fig. 3, right, showing the oscillation component in
light green and the HNL decay component from coherent/incoherent upscattering
production in darker green. One can see that the oscillation component does
not address the forward-peak in the MiniBooNE data while the decay component
primarily addresses that peak.

SUPPL. FIG. 8: Enlarged version of Fig. 3, right, in the paper. The prediction for the oscillation contribution to the $3+1+\mathcal{N}$ component is shown in light green. The prediction for the $\mathcal{N}$ decay component, from both coherent and incoherent upscattering, is shown in darker green. See text in this appendix for explanation of the data.  Fit type: | 3+1-only | 3+1-complete
---|---|---
$\chi^{2}_{app}$ | 48 | 79
$N_{app}$ | 2 | 2
$\chi^{2}_{dis}$ | 557 | 557
$N_{dis}$ | 3 | 3
$\chi^{2}_{glob}$ | 615 | 664
$N_{glob}$ | 3 | 3
$\chi^{2}_{PGF}$ | 10 | 28
$N_{PGF}$ | 2 | 2
$p$-value | 7E-03 | 8E-07
$N\sigma$ | 2.5$\sigma$ | 4.8$\sigma$
SUPPL. TABLE 3: PGF test results for 3+1 and our model, where Eqs. A11 and A12
explain how $\chi^{2}_{PGF}$ and $N_{PGF}$ are determined.