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

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$ eV$^2$. 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 O(1 eV), in addition to the three Standard Model (SM) neutrino states. In this model, ν µ → ν e appearance, ν e disappearance, and ν µ disappearance searches should all point to neutrino oscillations at L/E ∼ 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].
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+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θ, and mass splitting, ∆m 2 , for the oscillation; and the HNL mass, m 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 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 N 1 = ν 4 can decay, but the lifetime is typically longer than the age of the Universe [31,36]. The keV-scale mass state, 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 N ≡ 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 arXiv:2105.06470v5 [hep-ph] 27 Sep 2021 Used to Test References (Flux Type) Type of Fit νe disappearance [42][43][44][45][46] (Reactor) νe disappearance [47][48][49] (Source) νµ →νe appearance [50,51](π/µ DAR) ↑ νµ → νe appearance [52](π/µ DIF) 3+1-onlȳ νµ disappearance [53][54][55][56] (π/µ DIF) ↓ νµ disappearance [54,[57][58][59] (π/µ DIF) 3 + 1 + N [10] (MiniBooNE BNB ν) N TABLE I. 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 II. light neutrino states [38][39][40][41] to form observable oscillations. In the 3+1+N -decay model, any given SBL experiment may be sensitive to signatures of 3+1 oscillations only, N → νγ only, or both. In our model, the production and decay of N j is due to a dipole portal interaction between left-handed neutrinos, photons, and right-handed HNLs. The N j are added to the SM Lagrangian using the following term [15,17]: where the ν i correspond to the light neutrino mass states and F µν is the electromagnetic field strength. The dimension-full d αj couplings control the strength of the electromagnetic interactions between neutrino species, namely the strength of process like N j → ν i γ. Note that d αj reflects the effective dipole coupling of N j to the weak eigenstate ν α . This leads to two production mechanisms for N j : coupling to virtual photons produced in meson decays, such as π 0 , and Primakoff upscattering of active neutrinos to 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 N . This follows if d α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 ≡ d αj as a universal coupling. An alternative explanation if d is found to be large is that the d αj vary with family member, suppressing decays of the first and second mass states. In this case, we define d ≡ d α3 , where the coupling of N 3 to all light neutrino species is assumed to be the same. In the decay, the polarization of N must be considered [60][61][62]. The photon from a right-handed Dirac N decay has a (1 − cos θ) distribution, where θ is the angle between the N and photon momentum vectors; conversely, a left-handed Dirac N will decay with a (1 + cos θ) distribution for the photons.
Production through an unpolarized virtual, off-shell photon yields an equal combination of right-handed and lefthanded N , leading to an effectively isotropic photon decay distribution. For the case of upscattering, where the N is produced from an interaction with a left-handed muon neutrino, the outgoing N will be right-handed and the (1 − cos θ) 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: where U αj is the extended 6 × 6 neutral-lepton mixing matrix [1].

CONSTRAINTS
Three SBL experiments have relevant limits to 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 ν µ -electron elastic scattering (ES) data to the SM prediction [63]. At larger N masses, which are kinematically inaccessible in the former process, contributions from ν µ -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 ν µ -electron scattering results in agreement with the SM, placing a limit on N [17].

OSCILLATION GLOBAL FIT
If the maximum energy of the neutrino source is too small to produce the heavier 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 N decay as "3+1-only" experiments. The references for relevant "3+1-only" experiments used in this analysis are provided in Table I (top), including experiments with anomalies of significance from 2σ to 4.8σ and experiments consistent with νSM oscillations. The experiments, individually listed, can also be found in Supplemental Table II. We use these experiments to determine the oscillation parameters sin 2 (2θ) and ∆m 2 in a 3+1only fit. Notably, the "3+1-only" experiments exclude all MiniBooNE results, as we will then use MiniBooNE to fit for the N parameters m N and d.
As shown in Table I (bottom), we used MiniBooNE data to fit for the N parameters m N and d, given the oscillation parameters from the 3+1-only fit. MiniBooNE has excesses in three appearance data subsets: neutrinomode [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-masssplitting ∆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 × 4 neutral-lepton-mixing submatrix that relates the lightest mass states to their flavor states. For U αk , where α 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θ ee = 4(1 − |U e4 | 2 )|U e4 | 2 (electron flavor disappearance); sin 2 2θ µµ = 4(1−|U µ4 | 2 )|U µ4 | 2 (muon flavor disappearance); and sin 2 2θ µe = 4|U µ4 | 2 |U e4 | 2 (appearance). This implies that the electron and muon flavor disappearance signals must be consistent with the ν µ → ν e appearance signal, limiting the range of sin 2 2θ eµ .
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 II. 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 Ice-Cube [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+1only global fit omitted the MiniBooNE neutrino-mode, antineutrino-mode, and off-axis appearance data sets.
The best-fit parameters are ∆m 2 = 1.32 eV 2 and sin 2 2θ eµ = 6.9 × 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 × 10 −6 (4.5σ) that the data are explained by the same parameters. Without MiniBooNE appearance in the fit, the probability increases to 7×10 −3 (2.5σ). Thus, the tension is, in large part, due to the MiniBooNE appearance data set, which we hypothesize has the additional component of N -decay, and, hence, poor agreement with 3+1-only.

MINIBOONE ANALYSIS
In order to fit MiniBooNE data for N decay, we wrote a Monte Carlo simulation for the production and decay of N in the Booster Neutrino Beam in neutrino mode. Two processes were included for production: Dalitz-like π 0 decay and Primakoff upscattering νA → N A. The latter is by far the dominant N -production mode for 10 MeV < m N < 1000 MeV. Therefore, we neglected the π 0 decay contribution throughout this study. For the Primakoff mode, we generated incident ν µ and ν e events from the MiniBooNE neutrino-mode flux [94]. We then simulated the upscattering rate on both standard uppercontinental crust nuclei [95] and the MiniBooNE CH 2 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 N events, predominately forward peaked due to the 1/t 2 dependence of the differential cross section.
Simulated 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, ef f , we performed a linear fit to the reconstructed gamma-ray efficiency as a function of true energy [96], ef f = (−0.12 GeV −1 ) * E true + 0.29, which we used to weight the N → νγ 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, θ, of the MiniBooNE events to the prediction from 3 + 1 + N -decay. However, the systematic errors for this distribution have not been released by MiniBooNE. They are only available for E QE ν , which is a combination of E v and θ given by [97]   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 θ is measured. Thus, it is applicable to the oscillation component of the excess. Though E QE ν has no physical meaning when applied to the photons from N decay, the decay kinematics cause most events to be show up at low E QE ν . We performed a fit to the MiniBooNE excess in E 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 N } parameter space, computed assuming Wilks' theorem with two degrees of freedom is valid for the test statistic [98]. We found a region of parameter space consistent with both distributions at the 95% CL near d = 3 × 10 −7 GeV −1 and m N = 400 MeV. Fig. 2 shows that the allowed regions from MiniBooNE fits are substantially lower in d than the NOMAD or

RESULTS
(left) and cos θ (right) distributions of the MiniBooNE excess for a representative point of the 3 + 1 + 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.
LSND limits. The overlapping solution is also at substantially higher m 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 × 10 −7 GeV −1 and m N = 376 MeV, indicated by the star in Fig. 2. This corresponds to the best fit to the E QE ν distribution within the joint 95% CL allowed region from the E QE ν and cos θ fits. Table II shows the χ 2 values for the 3 + 1 and 3 + 1 + N -decay fits to both distributions, indicating significant improvement for the 3 + 1 + N -decay model. The global oscillation fit gives tight constraints requiring ∆m 2 ≈ 1.32 eV 2 , but allows values of sin 2 2θ µe ∈ [3 × 10 −4 , 2 × 10 −3 ] at the 90% CL. The same N decay fit procedure outlined above has been performed for each end of the allowed sin 2 2θ µe range. In each case we again examined the {d, m N } point that best fits the E QE ν distribution within the joint 95% CL region. The χ 2 values for these fits are also given in Table II, indicating a preference for a smaller oscillation contribution in MiniBooNE in order to explain both the cos θ and E QE ν distributions via N → νγ.  This figure includes both the global fit oscillation contribution for ∆m 2 = 1.3 eV 2 and sin 2 2θ ee = 6.9 × 10 −7 , and HNL decay contribution for d = 2.8 × 10 −7 GeV −1 and m 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 Mini-BooNE 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 ±4 ns of the observed ν µ events. The time of flight depends on the location at which the HNL is produced. For the preferred m N ∼ 400 MeV, the majority of N events come from upscattering within the MiniBooNE detector followed by N decay after travel distances of O(50 cm). This leads to timing delays of < 2 ns, well within the MiniBooNE constraint. CONCLUSION We have explored a 3+1+N -decay model through fits to the 3+1-only and MiniBooNE-neutrino-mode data sets. The former yields best-fit oscillation parameters of ∆m 2 = 1.3 eV 2 and sin 2 2θ eµ = 6.9 × 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 × 10 −7 GeV −1 . This model produces a ∆χ 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 ∆χ 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 ∼ 85% accuracy [100,101]. Because the experiments run in the same beamline and are located within ∼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 ∼ 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 ∆ 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 ∆-decay events with a proton. Also, the energy-angle correlation of the photon in the decay, which depends strongly on m N , can be used for background rejection since this parameter is predicted to have a narrow range of values. The oscillation rates for Mi-croBooNE 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 ν µ ES searches at MINERνA and NOνA, as well as dedicated searches for single photons in T2K [103]. Additionally, at neutrino energies of O(100 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 O(1 eV) that participate in oscillations and 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.

Supplemental Material
Appendix A: Further Explanation of the MiniBooNE N Simulation

Primakoff Upscattering Simulation
In order to simulate the Primakoff process νA → N A for a given dipole d coupling and HNL mass m N , we first generated a sample of 5×10 5 ν µ and ν e events with energies according to the MiniBooNE neutrino-mode flux from [94] and initial azimuthal angles generated randomly in cos θ 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 CH 2 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 I. Upscattering events inside the MiniBooNE detector happen exclusively off of CH 2 . 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]: where α is the fine structure constant, d is the dipole coupling, E ν is the SM neutrino energy, m N is the mass of the heavy neutrino, M is the target mass, t = −(p N − p ν ) 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 4 N 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 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].
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 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(θ)| < 1. From this we calculate the probability of scattering off a given nucleus A k (with atomic fractional abundance in dirt/air/CH 2 F k ) as 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]: The 1/t 2 dependence of Eq. A1 creates a preference for E N ≈ E ν and cos(θ) ≈ 1. At this stage, if the scattering angle of the 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 N → νγ in MiniBooNE, which is given by the following expression [17]: where L enter/exit denote the distance between the creation point of the HNL and the entry/exit point in the Mini-BooNE detector, respectively. The culmination of this simulation chain is a weighted sample of 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 × 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 where φ(E ν ) [ν/POT] is again the neutrino flux in MiniBooNE [94], σ A [cm 2 ] is the total upscattering cross section for nuclear target A, and n A [nuclei/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 N } parameter space for the three sin 2 2θ values considered in Table II. 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 ∆χ 2 plots in Figure 2. Here one can see that for the largest sin 2 2θ, there is no longer a closed contour in either the E QE ν or cos θ fit at 3σ CL. Figure 3 shows each of the E QE ν and cos θ predictions from Table II 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].

Appendix B: Further Discussion on the Fits
This Appendix provides further information on the fit results presented in this paper. Table II 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 + N model. This model has two contributions: the oscillation contribution, established through a 3+1 oscillation fit to experiments in Column 3 and the 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 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 II. 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  Table II. (dis) and appearance (app) sets. One defines the χ 2 and degrees of freedom as: We quote the associated probability as the tension. Table III 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 N decay. As stated earlier, this procedure finds a best fit point of ∆m 2 = 1.32 eV 2 and sin 2 2θ eµ = 6.9 × 10 −4 , with an allowed range sin 2 2θ µe ∈ [3 × 10 −4 , 2 × 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  Table II. Three different eV -scale oscillation amplitudes are considered: sin 2 2θ = 6.9 × 10 −4 (top left), sin 2 2θ = 3 × 10 −4 (top right), and sin 2 2θ = 2 × 10 −3 (bottom). 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. 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 × 10 −3 , with the relevant parameters summarized in Table III. For completeness, we also provide a global fit including all experiments listed in Table II. 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θ µe ) = 1×10 −3 and ∆m 2 41 = 1.32 eV 2 . 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 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. 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 × 10 −7 , with the relevant parameters summarized in Table III.
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 + 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 ∆m 2 = 1.3 eV 2 and sin 2 2θ = 6.9 × 10 −4 . The darker pink regions indicate the HNL decay component from both coherent and incoherent upscattering production, with d = 2.8 × 10 −7 GeV −1 and m N = 376 MeV.
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.