Explaining the MiniBooNE excess by a decaying sterile neutrino with mass in the 250 MeV range

The MiniBooNE collaboration has reported an excess of $460.5\pm 95.8$ electron-like events ($4.8\sigma$). We propose an explanation of these events in terms of a sterile neutrino decaying into a photon and a light neutrino. The sterile neutrino has a mass around 250 MeV and it is produced from kaon decays in the proton beam target via mixing with the muon or the electron in the range $10^{-11} \lesssim |U_{\ell 4}|^2 \lesssim 10^{-7}$ ($\ell = e,\mu$). The model can be tested by considering the time distribution of the events in MiniBooNE and by looking for single-photon events in running or upcoming neutrino experiments, in particular by the suite of liquid argon detectors in the short-baseline neutrino program at Fermilab.


Introduction
The MiniBooNE collaboration has published evidence for an excess of electron-like events of 381.2 ± 85.2 above their background expectation [1], confirming previous hints present in both, neutrino and anti-neutrino beam modes [2]. The combined excess of 460.5 ± 95.8 events corresponds to a significance of 4.8σ. The collaboration presents the results in the context of ( -) ν µ → ( -) ν e neutrino oscillations, under the hypothesis of a sterile neutrino with a neutrino mass-squared difference ∆m 2 of order 1 eV 2 , motivated by a previous claim from LSND [3]. The interpretation of the above mentioned results in terms of neutrino oscillations with an eV-scale sterile neutrino is in strong conflict with data on ( -) ν e and ( -) ν µ neutrino disappearance at the ∆m 2 ∼ 1 eV 2 scale [4][5][6]. This motivates to look for other new-physics explanations, beyond sterile neutrino oscillations.
In this paper we propose a sterile neutrino in the 150 to 300 MeV mass range, which is produced in the beam target from kaon decay via mixing either with electron or muon neutrinos. Subsequently it decays inside the MiniBooNE detector into a photon and a light neutrino. Since the electromagnetic shower of a photon inside MiniBooNE cannot be distinguished from the one of an electron or positron the photon can explain the observed excess events. We study the energy and angular spectra and predict a specific time distribution of the events. In order to obtain a reasonable fit to the angular distribution, we are driven to heavy neutrino masses around 250 MeV, which can be produced by kaon decays in the beam target. For lighter neutrino decays, the signal is too much forward peaked, inconsistent with MiniBooNE data [7]. Then the heavy neutrinos are only moderately relativistic and therefore our signal has a specific time structure, which provides a testable signature of our model [8]. The required parameters are consistent with all laboratory, astrophysics, and cosmology bounds. Current bounds and sensitivities of the upcoming short-baseline program at Fermilab for N → γν with the heavy neutrino N in the relevant mass range have been discussed in ref. [8]. Our model differs from various previously discussed explanations of the MiniBooNE and LSND anomalies based on the decay of a sterile neutrino. In the explanations of refs. [9,10] and [11,12] the heavy neutrino is produced by ν µ scattering inside the detector and has to decay with a very short lifetime into a photon or an e ± pair, respectively. The photon model from ref. [10] is by now excluded by searches for radiative neutrino decays from kaons by the ISTRA+ experiment [13], see also [14,15]. For other decay scenarios and related work see refs. [6,[16][17][18][19].
The article is structured as follows. In section 2 we introduce the model and in section 3 we describe the calculation of the MiniBooNE signal, including the time, energy, and angular event distributions. We present our χ 2 fit to the data in section 4. The results are discussed in terms of the model parameters in section 5, which includes also a discussion of other constraints on the model and possible tests in existing or upcoming experiments. In section 6 we conclude. Details of the heavy neutrino flux calculation are given in appendix A, in appendix B we discuss the impact of the timing cut on the MiniBooNE fit result.

The model
We consider one heavy Dirac neutrino N with mass m N that mixes with the SM neutrinos, parameterized by the leptonic mixing matrix U . The sub-matrix U i with = e, µ, τ and i = 1, 2, 3 is approximately the PMNS matrix that gives rise to neutrino oscillations, and the matrix elements U 4 allow N to interact with the weak currents and the lepton doublets of the Standard Model. Focusing on the case m N = O(100) MeV, we consider effective four-fermion interactions between the heavy neutrino and the SM particles, which are the mesons and leptons at this energy scale. Of particular importance is the following effective operator: where G F is the Fermi constant, q u and q d are up-type and down-type quarks, respectively, V is the CKM matrix, and is a charged lepton. Fixing the CKM matrix element to V us , the operator in eq. (1) allows us to calculate the branching ratio of the kaon into a lepton = e, µ and the neutrino N . For later use we define the following quantity: which takes into account the mixing of the heavy neutrino and the kinematical factors related to the finite mass of the neutrino [20]. Here, x i = m i /m K and we use Br(K → µν) = 0.636 and Br(K → eν) = 1.6 × 10 −5 . The factor ρ (m N ) is normalized to the branching ratio of K → µν, since we use the kaon induced ( -) ν µ flux in MiniBooNE to derive the heavy neutrino flux in both cases, K → N µ and K → N e, see appendix A.
In order to obtain the decay N → νγ into a light neutrino and a photon we introduce another effective operator to parameterize the possible interaction of N with a photon and light neutrinos via its magnetic moment [9,21] 1 : with the electromagnetic field strength tensor F µν = ∂ µ A ν −∂ ν A µ , the anti-symmetric tensor σ µν = γ µ γ ν − γ ν γ µ , and the unknown energy scale Λ. The operator O N νγ could be created at the loop level, for instance, such that we expect 1/Λ to be a combination of an inverse mass, unknown coupling constants, and a typical loop suppression factor. The operator in eq. (3) allows N to decay via the process N → νγ, with the total width in the rest frame of N given by To predict the energy and angular event spectra in MiniBooNE, we will need the differential decay rates with respect to the photon momentum p γ and the angle θ between the photon and N momenta in the laboratory frame: The minimum value of p γ is in backward direction, p γ,min = (E N − p N )/2, and the maximum value in forward direction, p γ,max = (E N + p N )/2. The phenomenology of the magnetic moment operator from eq. (3) has been studied extensively in ref. [15], see also [8,24,25] for recent considerations. In general this operator provides also a production channel for the heavy neutrinos [14,15]. Comparing with the results of ref. [15] we will see that for decay rates relevant for our scenario, the production via mixing and weak boson mediated kaon decay as described in relation to eq. (2) will be the dominant production mechanism.
The neutrino mixing parameters U 4 allow for various decay modes of N into SM particles via weak boson exchange; depending on its mass into a number of leptons, or also into a 1 The operator in eq. (3) has been chosen as a specific example for a possible decay mechanism, which we use below to study the relevant phenomenology. Other operators (including dimension-6 operators) inducing N → νγ in the case of Majorana neutrinos have been considered e.g., in refs. [22,23]. lepton and one or more mesons, which have been computed e.g. in refs. [26,27]. In the mass range of interest to us, m π < m N < m K , the dominant decay modes are N → ± π ∓ and N → νπ 0 . Using the results of ref. [27] the decay rate can be estimated by Here, g(m π , m lept , m N ) is a dimensionless kinematical function depending on the decay channel [27], "lept" indicates either a light neutrino or a charged lepton of flavour = e, µ, and f h ≈ 130 MeV is the pion decay constant. As we will see below, for large portions of the parameter space for Γ N →νγ and U 4 required to explain the MiniBooNE events, the decays N → leptπ will be sub-leading compared to N → νγ.
Note that a decay width of the scale indicated in eq. (4) corresponds to lifetimes much shorter than milliseconds, and therefore our sterile neutrino decays well before Big Bang nucleo-synthesis and hence does not affect cosmology. Heavy neutrinos in the 100 MeV mass range are at the border of being relevant for supernova cooling arguments. Limits from supernova 1987A on heavy neutrino mixing are avoided in our scenario [28], while the limits due to the magnetic moment operator derived in ref. [15] will be relevant in part of the parameter space able to explain the MiniBooNE excess, see also [29,30].
To summarize, the relevant phenomenology of our model is determined by three independent parameters, which we chose to be the heavy neutrino mass: m N , the mixing with the e or µ flavour: |U 4 | 2 , and the decay width into the photon: Γ N →νγ . We will present the parameter space where the MiniBooNE excess can be explained in terms of those three parameters in section 5 below.

The MiniBooNE excess events
Our analysis proceeds as follows: first we construct the kaon flux at the BNB from the given flux of the muon neutrinos. From the kaon flux we derive the flux of the heavy neutrinos and work out its time structure. Then we calculate the energy and angular spectra of the photon from the heavy neutrino decays inside the detector and inside the time window defined by the MiniBooNE collaboration.
In order to calculate the flux of heavy neutrinos Φ N (p N ) we proceed as follows. We depart from the kaon contribution to the ( -) ν µ fluxes provided by the MiniBooNE collaboration ref. [31]. Assuming that this flux is dominated by the two-body decay K → νµ we reconstruct the initial kaon flux, from which in turn we can calculate the heavy neutrino flux at MiniBooNE by taking into account the modified angular acceptance of the detector due to the non-negligible effect of the heavy neutrino mass on the angular distribution. Details of this procedure are provided in appendix A. Note that the flux Φ N (p N ) obtained in this way depends on the mass of the heavy neutrino, which we keep implicit to simplify notation.

Time spectrum
A heavy neutrino with momentum p N arrives at the detector at distance L after a time with the MiniBooNE baseline L 540 m. The ultra-relativistic light neutrinos all arrive after t 0 1.8 µs, the heavy neutrinos generally arrive later. In order to calculate the time distribution of the events, we first convert the neutrino flux Φ N (p N ) into a function of time, with the Jacobian |dp N /dt| = p N t/(t 2 − t 2 0 ), which follows form eq. (8). In the decay model, an additional momentum dependence appears due to the effect of the Lorentz boost on the decay rate, which leads to a factor m N /p N , see eq. (15) below. Finally, to construct the time spectrum we need to include the time structure of the proton beam, which we approximate with a step-function being non-zero from t = 0 to t = δt = 1.6 µs [32]. Therefore, we obtain the time spectrum T (t) in the following way: We show the time distribution of the decay events inside the detector for a typical heavy neutrino mass in fig. 1. The contribution of the monochromatic peak from the stopped kaon decays is visible in the discontinuous part of the red curve around t = 3 µs. It is important to notice, that the neutrino appearance analysis from MiniBooNE considers only events that occur between t 0 and t 0 + 1.6 µs after each beam spill [33]. The fraction of our heavy neutrino signal inside the analysis window is denoted in blue in the figure, those that arrive after t 0 + δt are too late to be included and are denoted in red. The fraction of the events inside the timing window is 41% (34%) in the neutrino (antineutrino) mode. Therefore, we predict a significant fraction of delayed events. Those could be searched for in the MiniBooNE data. Note that MiniBooNE records events within a time window of about 19.8 µs around each beam spill (cf. ref. [32]). A detailed investigation of events in this time region can be a definite test of our model. Using timing information to test heavy neutrino decays has been suggested previously in ref. [8].

Event numbers, energy and angular spectra
The number of heavy neutrinos that decay inside the detector are obtained by integrating over the heavy neutrino flux φ N , together with the probability P dec that the long lived particles decay within its fiducial volume. Furthermore, a detection efficiency has to be included that is an empirical function of the signal energy, here approximated with the momentum of the photon, and the decays have to occur inside a timing window as discussed above. These considerations are summed up in the following master formula: Here, POT denotes the number of protons on target, which is 12.84 (11.27) × 10 20 for the neutrino (antineutrino) mode. The factor ρ (m N ) has been defined in eq. (2) and it includes the mixing matrix element |U 4 | 2 and the branching ratio of the kaon decays into heavy neutrinos. Br νγ = Γ N →νγ /Γ tot is the branching ratio for the decay N → νγ, with Γ tot being the total decay width of N . In the relevant mass range we have Γ tot ≈ Γ N →νγ + Γ π with Γ π given in eq. (7). Furthermore, A MB = π(5 m) 2 is the effective area of the MiniBooNE detector, andˆ is the MiniBooNE detection efficiency [33] (p γ ) averaged over the photon momentum distribution for a given p N . P dec is the probability that the heavy neutrino decays inside the detector, and w time is a timing-related weight. Using the heavy neutrino arrival time t N from eq. (8) the latter is given by For the decay probability we have where L 1 , L 2 denote the distance of the front and back ends of the detector from the beam production and we assume an effective value of ∆L ≡ L 2 − L 1 = 8 m. Here Γ tot is the heavy neutrino width in the rest frame, a factor m N /E N takes into account the boost into the lab frame of the detector, and L i × E N /p N is the time the neutrino needs to reach the position L i . In eq. (15) we have used an approximation, which holds for the MiniBooNE baseline of L ≈ 540 m when Γ tot 10 −15 MeV. In this approximation the number of events is proportional to We remark that in order to fit the observed shapes of the angular and energy spectra, in principle we have to recast the energy of the photon-induced Cherenkov shower (here assumed to be identical to the momentum of the photon from the heavy neutrino decay) into the energy of a light neutrino, E ν , assuming a quasi-elastic scattering event. It turns out, however, that this recasting yields a relative difference on the percent level, such that we will keep the photon energy as our proxy for the reconstructed neutrino energy in the following, for simplicity.
Using the linear approximation for the decay probability (15), and the differential decay widths from eqs. (5) and (6), we construct the predicted angular spectrum A(z) with z ≡ cos θ and energy spectrum E(E ν ) in the following way: where the neutrino flux depends on the horn mode and includes forward and backward decays of the parent meson. The photon momentum is related to the scattering angle by In our model the ratio of signal events for the two horn polarisations (R pred ) is determined from the corresponding fluxes, and can be calculated from eq. (11) for a given heavy neutrino mass. We show R pred as a function of m N in fig. 2 for the heavy neutrino being produced together with a muon. The ratio is almost identical when the heavy neutrino is produced together with an electron, aside from the fact that larger values for m N are kinematically accessible.

Fit to the data
In order to test our model we perform a fit to both, the angular and energy spectra. would fit the angular and energy information simultaneously by using the 2-dimensional distribution of the data. Unfortunately this information is not available, and therefore we have to fit the energy and angular spectra separately and check if results are consistent. 2 In each case the fit is done fitting simultaneously both the neutrino and anti-neutrino spectra.
Using the results of the previous section, we parameterize our model with two effective parameters, which we chose to be N total and the sterile neutrino mass m N . The predicted number of events in a given bin i of the energy or angular data is given by N ν f ν i (m N ) and Nνfν i (m N ) for the neutrino and anti-neutrino polarization, respectively. Here, where R pred (m N ) is the predicted ratio of events in the neutrino and anti-neutrino modes shown in fig. 2, and f ν i (m N ), fν i (m N ) are the predicted relative contributions for each bin, normalized to 1. They are derived from the corresponding differential spectra given in eqs. (16) and (17).
We define the following χ 2 -function to perform the analysis: Here i labels the angular or energy bins, a labels the background contributions (sum over a is implicit), O ν i and Oν i are the number of events in each bin i for the ν andν mode respectively. B a i andB a i are the different a background contributions in each bin i, b a andb a are the pull parameters that account for their uncertainty σ a andσ a , which are taken from table 1 of ref. [1], see Tab. 1. Possible correlations are not taken into account. Furthermore, a totally uncorrelated systematic uncertainty of 20% is considered in each bin to account for possible spectral shape uncertainties: The results of the angular and energy analyses are shown in figure 3. To be specific, we assume the production mode K → N µ, results for K → N e are very similar. The energy spectrum provides a best fit point at m N = 250 MeV and N total = 640 with closed allowed regions. At 68% confidence level our fit allows the masses to vary between 190 MeV and 295 MeV and normalisations between 425 and 865 events, consistent within 1σ with the observed number of excess events N obs = 460.5 ± 95.8, as indicated by the green band in the plot. Note that this comparison is only indicative, since the data used in our fit (taken from fig. 14 of ref. [1]) uses a lower energy threshold and therefore the number of excess events is somewhat larger, consistent with our best fit value. The best fit point has χ 2 min /dof = 58.1/36 which corresponds to a p value of about 1% (see discussion below). In contrast, the angular spectrum only provides an upper bound on N total which is in some tension with the energy fit.
In order to investigate the quality of the fit we show in figure 4 the predicted energy and angular spectra for m N = 250 MeV and N total fixed to 400, chosen within the 1σ range of the observed value. From the left panel we see that our model explains well the excess events in the energy spectrum. A significant contribution to the χ 2 comes from bins above 1 GeV, where a signal is neither observed nor predicted. This explains the rather low p-value of only 1%. The right panel shows that the angular shape for the anti-neutrino mode is in good agreement with the observations, while the signal is somewhat too much forward peaked in the neutrino mode. From comparing the neutrino-mode spectra in the two panels (blue histograms), the tension between energy and angular fit is apparent. While The angular fit provides an upper limit on N total . In green we show the measured excess of events and its 1σ uncertainty. We assume here heavy neutrino mixing with the muon; results for the electron are very similar.
the energy spectrum would prefer to increase the normalization, this would clearly worsen the prediction in the forward angular bin. Note however, that the largest contribution to the angular χ 2 comes from the three bins around cos θ = 0, including the one with the downward fluctuation. Those bins are difficult to explain by any smooth function. A general discussion of the angular event distribution in decay models can be found in ref. [7]. We stress that to definitely assess the viability of our model a joint energy and angular fit should be performed, including detailed acceptances and efficiencies suitable to our signature. Let us also mention that both the timing cut and the implementation of the angular acceptance is important to predict the angular shape, since both affect mostly the signal from the decay of "slow" neutrinos, which give the main contribution to events with cos θ < 1. In appendix B we show the fit results without imposing the 1.6µs timing cut, which leads to an improved angular fit. Below we proceed under the assumption that our model does provide an acceptable fit to MiniBooNE data. neutrino decay rate Γ N →νγ . In fig. 5 we show the 1 and 2σ contours for those two parameters for the neutrino mass fixed at the best fit point. The straight part on the left side corresponds to the linear approximation for the decay probability, eq. (15), where event numbers are proportional to the product |U 4 | 2 Γ N →νγ . The linear approximation breaks down when the decay length becomes shorter than the MiniBooNE baseline and most of the neutrinos decay before reaching the detector. This leads to the upturn of the allowed region visible in the plots for decay rates Γ N →νγ 10 −15 MeV. This value depends only weakly on m N and defines a minimum value of |U 4 | 2 needed to explain the excess of roughly 2 × 10 −11 . The lower limit on |U 4 | 2 is shown as a function of the heavy neutrino mass in fig. 6. The dark and light orange shaded regions correspond to the 1σ and 2σ range for m N as shown in fig. 3. Note that we do not consider masses below 150 MeV in order to avoid N production due to pion decays. We see that the excess can be explained by a wide range of values for the mixing and for the decay rate. Let us now consider other constraints on those parameters.

NOMAD:
The single photon signature predicted in our model can be searched for in various other neutrino experiments. A rather sensitive search comes from the NOMAD experiment at CERN. An overview of the experiment is given in ref. [37]. An analysis searching for single photon events (motivated by the MiniBooNE observation) yields 78 observed events in forward direction versus 76.6 ± 4.9 ± 1.9 expected, which was interpreted as a null result and an upper bound of 18 events at 90% CL has been set on single photon events [36]. We can interpret this bound as a limit within our model. We use the kaon-produced muon neutrino flux from ref. [38] and construct the heavy neutrino flux as described in appendix A. The number of heavy neutrino decays in the detector is estimated as in eq. (11). The POT is 2.2 × 10 19 and we use a constant reconstruction efficiency of 90%, an analysis efficiency of 8%, and a trigger efficiency of 30% [36]. In order to take into account an analysis cut on the observed energy we consider only heavy neutrinos with  Figure 5: Parameter region in the plane of Γ N →νγ and |U 4 | 2 that is consistent with the observed MiniBooNE excess at 1 and 2σ for the neutrino mass fixed at the best fit point. The left (right) plot assumes that N is produced from a kaon decay with an associated electron (muon) and corresponds to m N = 260 (250) MeV. Also shown are upper limits on |U 4 | 2 from NA62 [34] and E949 [35] and the region excluded by NOMAD from the search in ref. [36], interpreted in our model, as well as the region for Γ N →νγ disfavoured by SN1987A [15]. momentum greater than 1.5 GeV. We do not apply any time window for the events.
The parameter space excluded by NOMAD by requiring that less than 18 events are predicted is shown as the dark gray shaded region in fig. 5. Since the baseline of NOMAD is shorter than MiniBooNE and the neutrino energies are higher, the decay rate for which neutrinos start decaying before reaching the detector is shifted to higher values of Γ N →νγ for NOMAD compared to MiniBooNE and NOMAD excludes the "non-linear" part of the parameter space. In the linear regime for the decay probability the NOMAD bound is always consistent with the value |U 4 | 2 Γ N →νγ required to explain MiniBooNE. We have checked that the predicted number of events in NOMAD is about 5.4 × 10 −3 times the signal events in MiniBooNE, with very little dependence of this number on m N within the interesting range. Therefore, the NOMAD bound limits the available parameter space to the linear regime but does not provide a further constraint on the range of the parameters.
Limits from kaon experiments. Due to the long lifetimes of the heavy neutrinos the vast majority of the produced N decay outside the detectors in most kaon experiments. However, an observable feature of their presence is given by an additional peak in the spectrum of the lepton from the decaying kaon. The NA62 experiment has recently published a search for heavy neutral leptons that are produced in kaon decays. Not observing candidates for kaon decays into heavy neutrinos, they placed upper limits at the 90% CL of around |U 4 | 2 ∼ 10 −7 for = e, µ and heavy neutral leptons with masses between 170 and 448 MeV for = e and between 250 and 373 MeV for = µ [34]. Earlier searches for heavy neu-  → eN (µN ). Excluded parameter space from peak searches in the kaon decay spectra of electron and muon from the NA62 [34] and E949 experiments [35] is shown as gray shaded regions. The regions disfavoured by SN1987A constraints on Γ N →νγ [15] are shown as blue shaded regions.
trinos from the E949 experiment studied the muon spectra from about 10 12 stopped kaon decays. In their analysis, the collaboration derived the still most stringent upper limits at the 90% CL on |U µ4 | 2 down to 10 −9 for heavy neutrinos with masses between 175 and 300 MeV [35]. We show the region excluded by E949 and NA62 in figs. 5 and 6 as a gray area. 3 To summarize sofar, as visible in fig. 6, several orders of magnitude in mixing are available to explain the MiniBooNE excess in this model. For a fixed value of m N , for each value of |U 4 | 2 in the allowed range, the value of the decay rate can be adjusted such that the event numbers are kept constant. Requiring N total = 400 events in MiniBooNE we have approximately where we have used the linear approximation for the decay probability and the fact that then event numbers are proportional to the product |U 4 | 2 Γ N →νγ . The power of the mass dependence has been obtained by fitting the numerical result with a power law and it is rather accurate in the range 150 MeV < m N < 300 MeV.
Constraint from SN1987A. As discussed in ref. [15], a heavy neutrino interacting via the operator (3) may contribute to the cooling rate of a supernova. In order to be consistent with the neutrino observation of SN1987A, too fast cooling has to be avoided. This argument can be used to disfavour certain regions in the parameter space of m N and Γ N →νγ . In the parameter region of our interest those considerations lead to a lower bound on the decay rate of approximately [15] For decay rates fulfilling this bound, the heavy neutrino is sufficiently trapped inside the supernova to avoiding too fast cooling. The bound shown in eq. (22) holds in the relevant mass range for our scenario, up to m N ≈ 320 MeV; heavier neutrinos are gravitationally trapped inside the supernova [40]. The region disfavoured by the bound (22) is shown in figs. 5 and 6 as blue shaded region, where in order to translate the bound into |U 2 4 | as shown in fig. 6 we assume our explanation of the MiniBooNE events, using the relation (21). We see that in the mass range where the SN bound applies, the mixing is limited to 10 −11 |U 2 4 | few × 10 −9 , while for m N > 320 MeV values of |U 2 4 | up to the kaon bounds of order 10 −7 are allowed. Once these constraints from the magnetic moment operator are imposed, the supernova limits on mixing from ref. [28] are satisfied; they disfavour the region m N 100 MeV and |U 2 4 | 10 −8 .
By comparing the decay rate from eq. (21) with the mixing induced decay rate in pions given in eq. (7) we find that Γ π Γ N →νγ for We see from fig. 6 that for the largest allowed mixing angles in the high-mass region this condition may be violated. In this case, the decays N → νπ 0 and N → ± π ∓ can provide an additional observable signature. Note, however, that in the region where the linear approximation breaks down, Γ π Γ N →νγ is satisfied and we can use Γ tot ≈ Γ N →νγ for calculating the decay probability according to eq. (14).

Other searches and tests of the model
The PS191 and E816 experiments: The dedicated PS191 experiment searched for displaced vertices from the decay of heavy neutrinos in the mass range from a few tens of MeV to a few GeV. Not having found such vertices PS191 placed limits on the mass-mixing parameter space [41,42]. It is important to notice that these limits are not applicable in the here considered model, because the decays of the heavy neutrino into a photon and a light neutrino do not produce a visible vertex in the decay volume.
We remark that the experiment observed an excess of electron-like events [43], which was interpreted as electron-neutrino scatterings in the calorimeter, but might be also induced by the photons from the N decay in our model. This finding is backed up by the PS191 successor at BNL, the experiment E816 [44]. Unfortunately the collaborations do not provide the details on the kaon flux, such that we cannot quantify the respective signal strengths in our model.

LSND and KARMEN:
The LSND [3] and KARMEN [45] experiments produce neutrinos from muon decay at rest and therefore heavy neutrinos with masses of 100 MeV will not be produced. The interactions of the 800 MeV proton beam with the target might produce a few slow-moving kaons, which could give rise to a heavy neutrino flux that is small compared to the one at MiniBooNE. Furthermore, the standard search in LSND and KARMEN requires a coincidence signal between a prompt positron and delayed neutron capture from theν e inverse beta decay process, which is rather distinctive from the pure electro-magnetic signal induced by the single photon decay in our model. Therefore, we predict a negligible event rate in those experiments.
T2K, NOνA, and other running neutrino experiments: Modern neutrino detectors, such as the near detectors of NOνA, T2K are generally not expected to confuse a single photon with charged current electron neutrino scattering due to their more sophisticated detectors. Recently the T2K collaboration published results for a search for heavy neutrinos [46] by looking for events with two tracks, for instance from the decays N → µ ± π ∓ or N → ± ∓ . The limits, comparable to those from PS191 and E949, are not applicable to our model. A search for single photon events in T2K has been published recently in [47]. We have roughly estimated the sensitivity of this result to our model and found that the resulting limits are weaker than the ones from NOMAD discussed above.
It is important to realize that the signal of our model mimics neutral current produced π 0 decays where one photon was not reconstructed, which may interfere with the control regions of any analysis and affect results in a non trivial way [19]. An analysis that searches for single photons in the data in all running neutrino experiments may be able to shed light on the MiniBooNE excess. The relative signal strength between experiments is fixed by the fluxes and allows to reject the hypothesis.

ISTRA+:
The ISTRA+ experiment searched for and excluded the process K ± → µ ± N , N → νγ for 30 MeV ≤ m N ≤ 80 MeV [13] and for very short neutrino lifetimes. With about 300 million events on tape, the experiment could in principle be sensitive to heavy neutrinos in the here considered mass range.
The Fermilab short-baseline neutrino program: The short-baseline neutrino (SBN) program at Fermilab consists of three liquid argon detectors in the booster neutrino beam line: SBND, MicroBooNE, and Icarus [48,49], with the MicroBooNE detector already running and producing results. A sensitivity study to heavy neutrino decays, including the photon decay mode has been performed in [8], see also [50]. Liquid argon detectors will be very suitable to search for the signal predicted here, since such detectors can discriminate photons from electrons. The main characteristics of the three detectors are summarized in   [48,49] compared to MiniBooNE. We assume the same POT as quoted in [8]. For MiniBooNE we sum the POT in neutrino and antineutrino modes. The row "Ratio" indicates the ratio of signal events relative to MiniBooNE based on the scaling with the factor in eq. (24). In the row "Events" we give the predicted number of events assuming 400 signal events in MiniBooNE.
tab. 2. Since they are located in the same beam as MiniBooNE we can roughly estimate the expected number of events by scaling with the proportionality factor where V is the detector volume and L the distance of the detector from the neutrino source. Note that the simple scaling with this assumes that the linear approximation for the decay probability holds for all baselines. In the table we give this scaling factor for each experiment relative to MiniBooNE ("Ratio"). Assuming 400 signal events in MiniBooNE, we can predict then the expected number of events by multiplying with this ratio. As is clear from the last row in tab. 2 a significant number of events is predicted for each of the three detectors, under the quoted assumptions on the available POT [8].
Atmospheric and solar neutrinos. The magnetic moment operator (3) can lead to the up-scattering of atmospheric or solar neutrinos to the heavy neutrino, which can give observable signals in IceCube [24] or dark matter detectors [25], respectively. The latter can test heavy neutrinos with mass below ∼ 10 MeV. The sensitivities of IceCube derived in ref. [24] from atmsopheric neutrinos are in the relevant mass range, but are about one order of magnitude too weak in Γ N →νγ to start constraining the parameter space relevant for our MiniBooNE explanation.

Conclusions
We presented a model with a heavy neutrino of mass around 250 MeV that is produced from kaon decays at the beam-target interaction via the mixing |U 4 | 2 ( = e, µ) and decays after traveling over several hundred meters into a light neutrino and a single photon via an effective interaction. We demonstrated that it is possible to account for the event numbers and spectral shape of the electron-like excess in the MiniBooNE data under the assumption that single photons are indistinguishable from single electrons. Some tension appears for the angular distribution of excess events, which are somewhat too much forward peaked.
A quantitative assessment of this tension requires a dedicated analysis of MiniBooNE data including a careful treatment of angular and timing acceptances. The excess events can be explained for a wide range of mixing parameters of roughly 10 −11 |U 4 | 2 10 −7 , consistent with existing bounds, see fig. 6. The model makes clear predictions and can be tested in the following way: • Delayed events in MiniBooNE: Due to the non-negligible mass of the heavy neutrino, we predict a characteristic time structure of the signal with a significant fraction (up to 60%) of events outside the time window corresponding to the time structure of the beam and assuming speed of light for the propagation to the detector. Therefore, the model can be tested by looking for delayed events in MiniBooNE data.
• Single photon events in SBN detectors: We predict a sizable number of single photon events in all three liquid argon detectors of the Fermilab short-baseline neutrino program (SBND, MicroBooNE, ICARUS). These detectors have good photon identification abilities and should be able to confirm or refute our hypothesis.
Assuming that the decay N → νγ is induced by the dimension-5 operator of the magnetic moment type, see eqs. (3,4), the decay rates required to explain MiniBooNE would correspond to a suppression scale Λ of roughly 10 4 TeV Λ 10 7 TeV. If the magnetic moment operator is generated at 1-loop level, we expect generically where g is a coupling constant and M np is the mass scale of some new physics. We see that for moderately small g, M np can be in the TeV range and therefore potentially accessible at the LHC. If N is a Majorana neutrino, there will be a contribution to the light neutrino mass via the type I seesaw mechanism of order m ν m 2 D /m N |U 4 | 2 m N , where m D |U 4 |m N is the Dirac mass of N . In the upper range of the allowed region for |U 4 | 2 , the seesaw contribution to m ν is too large. However, it is interesting to note that for m N 250 MeV and |U 4 | 2 10 −10 the seesaw contribution to m ν is of order 0.025 eV, just of the right order of magnitude for light neutrino masses. Furthermore, the magnetic moment operator from eq. (3) will induce also a contribution to the light neutrino mass via a 1-loop diagram [15], whose size in general depends on the UV completion of the operator (3). Both contributions-from seesaw and magnetic moment operator-can be avoided (or suppressed) if N is a Dirac (or pseudo-Dirac) particle. While our scenario has all ingredients to generate light neutrino masses, we leave it for future work to identify consistent models explaining light neutrino masses and mixing in this framework. For m N comparable to m K and for sufficiently large | p K |, also heavy neutrinos that are emitted backwards with respect to p K can reach the detector. This means, from the kaon flux we construct two heavy neutrino fluxes: one from the forward emitted N with cos θ = +1, and one from the backward emitted N with cos θ = −1. Both, the backward Φ bwd N (p N ) and the forward Φ fwd N (p N ) emitted fluxes are normalized to the original light neutrino flux Φ νµ (p νµ ). The peak in the kaon spectrum from the stopped kaons (shown in figure 7) gives rise to monochromatic heavy neutrinos of energy p N,0 . We add this separately to the analysis and call it the "monochromatic peak".
Geometrical acceptance: We work under the assumption that the kaon momentum is always parallel to the beam line. In the experiment, neutrinos (light or heavy) are not produced with cos θ = ±1, but rather with an angle θ = 0 + δθ, π − δθ, such that | cos θ| = 1 − ε. This deviation ε stems from the angles that are smaller than or equal to the solid angle of the detector, which we approximate with θ D ≈ tan θ D = r/L, where r is the radius of the detector and L is the distance from the source. The maximal acceptance angle of the heavy neutrino in the lab frame is given by here θ rest N is the kaon rest frame decaying angle. The component of the momentum perpendicular to the beam line is not affected by the kaon boost and the parallel one is given by expression (28). For small angles sin θ ∼ θ, cos θ ∼ ±1, the acceptance angle in the rest frame, for the backward and the forward decay, can be easily solved: Since the decay in the rest frame is isotropic, the heavy neutrino flux can be corrected by adding a geometrical factor given by the ratio between the maximum acceptance angles in the kaon rest frame for the heavy and light neutrinos: Here we assume that the angular acceptance for light neutrinos is already included in Φ νµ as provided by the collaboration. For small angles, eq. (32) turns into Note that only the light neutrinos decaying in the forward direction reach the detector, so the heavy neutrino acceptance angle, for both backward and forward directions, have to be compared to the light neutrino one in the forward direction.
We have checked that for the kaon energies at MiniBooNE the small angle approximation (33) works very well. On the other hand for the kaon energies in NOMAD, this approximation does not hold because the kaon momentum can be larger. 4 We have checked that taking the approximated expression for the geometrical factor for the NOMAD prediction gives an extra enhancement, i.e., we are somewhat over-predicting the number of events in NOMAD. This makes the limit somewhat too strong and is therefore conservative in what concerns the consistency with MiniBooNE, and hence we stick to the approximated expression.
The geometrical factors (33) can be expressed as a function of the heavy neutrino momentum performing an inverse boost Solving for p K we obtained where upper (lower) signs apply to the forward (backward) geometrical factors. Note that in the forward decay case p N starts from p N,0 and in the backward decay from 0. Finally, the heavy neutrino flux is given by: Total signal events Figure 8: Allowed regions at 1, 2, and 3σ for the energy (orange regions) and angular (dashed blue curves) in the N total versus m N parameter space, without imposing the 1.6µs timing cut. The left (right) panel assumes heavy neutrino mixing with the electron (muon). The best fit of the energy spectral (angular) fit is indicated with a black (blue) cross. In green we show the measured excess of events and its 1σ uncertainty.

B Impact of the timing cut
The timing cut of 1.6µs after each beam spill discussed in sec. 3.1 has a strong impact on the predicted event spectrum, since it removes events from "slow" heavy neutrinos, which would provide a less forward peaked angular distribution for the photon events. In order to illustrate the importance of the timing cut, we show in this appendix results without requiring arrival within 1.6µs, i.e., we include all events from N decays in the predicted signal. From Fig. 8 we see that in this case also the angular fit shows preference for non-zero signal event numbers, and the 1σ allowed regions overlap between the energy and angular spectral fits. The best fit point for the energy spectrum degrades only marginally from χ 2 min /dof = 58.1/36 with timing cut to 62.8/36 without timing cut in the case of muon mixing. For the electron mixing we obtain an energy spectrum best fit with χ 2 min /dof = 61.9/36. Without the timing cut the angular fit yields χ 2 min /dof = 32.1/18 (30.0/18) for the mixing with the muon (electron), corresponding to a p-value of 1.1% (3.7%). In Fig. 9 we show the resulting energy and angular spectra for an example point in the 1σ overlap region. In comparison with fig. 4 we clearly observe an improved angular fit, while still maintaining a good description of the energy distribution. As discussed in the main text, the formally still rather low p-value is a consequence of the scattered data points with small error bars in the tail of the distributions.