A Dark Seesaw Solution to Low Energy Anomalies: MiniBooNE, the muon $(g-2)$, and BaBar

A recent update from MiniBooNE has strengthened the observed $4.8\sigma$ excess of $e$-like events. Motivated by this and other notable deviations from standard model predictions, such as the muon $(g-2)$, we propose a solution to low energy anomalies through a dark neutrino sector. The model is renormalizable and can also explain light neutrino masses with an anomaly-free and dark $U(1)^\prime$ gauge symmetry broken at the GeV scale. Large kinetic mixing leads to s-channel production of heavy neutral leptons at $e^+e^-$ colliders, where we point out and explain a $\gtrsim 2\sigma$ excess observed in the BaBar monophoton data. Our model is also compatible with anomalous $e$-like events seen at old accelerator experiments, as well as with an excess of double vertex signatures observed at CCFR.


I. INTRODUCTION
The discovery of neutrino oscillations [1][2][3], and consequently of neutrino masses and mixing, implies that the Standard Model (SM) of particle physics is incomplete. Many extensions have been proposed to explain the origin of neutrino masses, with the Type-I seesaw mechanism [4][5][6][7][8][9][10][11][12] and its variants being the most well studied. Heavy neutral leptons (HNL) are the hallmark of such models and carry a lepton number violating (LNV) Majorana mass, which, barring theoretical prejudice, could take any value from sub-eV to 10 16 GeV. In recent years, renewed attention has been devoted to the MeV -GeV mass scale, as such states can be searched for in an expanding program of fixed-target, meson decay, and collider experiments [13][14][15][16][17][18][19], having consequences for cosmology and the baryon asymmetry of the Universe [20,21]. Two approaches are typically adopted: one of minimality, in which only new neutral fermions are introduced, e.g. [22], and, more recently, one in which the HNLs are considered part of a richer low energy dark sector [23][24][25][26][27][28][29][30][31][32][33][34][35][36], all the more compelling in view of the large abundance of dark matter in our Universe [37][38][39][40]. It has been pointed out that in the second approach the phenomenology can have unique features, requiring the reevaluation of existing bounds and offering new signatures, especially in the presence of multiple portals to the SM [41]. Such an extension of the SM would leave imprints, not just in neutrino experiments, but also in e.g. dark photon and dark scalar searches. Interestingly, some anomalies are present in these areas.
In this letter, we propose a coherent explanation of several experimental anomalies, generating the correct scale The full Lagrangian is given by where flavor indices are implicit, and we write the kinetic mixing between hypercharge and the U (1) mediator X µ , as well as scalar mixing between the Higgs and Φ explicitly. Here, where We diagonalize the mass matrix with a unitary matrix U , defined in terms of sub-blocks U ≡ U α U N U D L U D R T , such thatν m = Uν f ≡ ν N T contains the light neutrinos ν and the (n + 2d) HNLs N . At tree-level, the mostly-active neutrinos get a mass as in the inverse [52,53] and extended seesaw [54,55] models. At the one-loop-level, however, we find an independent finite contribution proportional to M N [56]. This is the same contribution found in Ref. [31], and is analogous to the minimal radiative inverse seesaw [49][50][51]. These independent tree-and loop-level contributions can have opposite signs, leading to cancellations if M X M N . We exploit this fact to achieve neutrino masses compatible with current data. We neglect loop corrections to other mass parameters in the matrix.
The massive dark photon, scalar, and HNLs only couple to the SM via portal operators. After symmetry breaking, the model has two CP-even scalars, the SM Higgs h , which contains a small Φ component with scalar where λ H is the quartic coupling of the Higgs, and a light mostly-dark ϕ . In the neutral gauge boson mass basis, we have a light Z vector boson that couples predominantly to the dark sector current (J µ D ), as well as to the SM electromagnetic (EM), and neutral current (NC), where we assume m Z g X v ϕ m Z , and define ε ≡ c W χ.

III. LOW ENERGY ANOMALIES
Our aim is to show that the model can explain several low energy anomalies, while simultaneously generating the correct scale for light neutrino masses. Since mixing in the light neutrino sector can be generated by appropriate choices of the M D matrix, we work under the simplifying assumption of a single active neutrino generation, in our case ν µ , denoting the lightest non-zero mass eigenstate by ν 3 . We require that m 3 0.05 eV, compatible with the scale suggested by neutrino oscillation experiments [57]. For concreteness, we pick n = 3 sterile and d = 1 vector-like dark neutrinos, although only the three lightest heavy neutrino mass eigenstates N j , j = 4, 5, 6, will be important for the phenomenology we discuss. The heaviest states N 7 and N 8 have masses of several GeVs, and are mostly-sterile states.
Our proposal is illustrated by two benchmark points (BPs), one exhibiting a left-right symmetry and one without. Their properties are shown in Table I but a detailed definition is left to Appendix A. The left-right symmetry in the dark sector of BP-A (ν c D L ↔ ν D R ) is achieved by setting Y L = Y R , and explains the vanishing entries in Table I. This can be shown to be related to CP conservation.
Let us comment on the generic features of our two BPs. We fix m Z = 1.25 GeV and ε 2 = 4.6 × 10 −4 for the dark photon. The three lightest HNLs all have O(100) MeV masses, and decay via neutrino and kinetic mixing as N i → N i−1 e + e − . Specifically, N 5 will typically decay with cτ 0 5 3 cm, leading to displaced e + e − vertices, while N 6 will decay more promptly, cτ 0 6 3 mm. In the case of N 4 , it can only decay into SM particles, N 4 → νe + e − , making it much longer-lived, cτ 0 4 100 km. In addition, N 4 is mostly sterile, which naturally leads to B(Z → N 4 N 4 ) B(Z → N {4,5,6} N {5,6} ), and explains why cτ 0 6 < cτ 0 5 . For concreteness, we fix m ϕ = 1 GeV, forbidding fast N 6 → N j ϕ decays and respecting perturbativity limits on the dark scalar quartic coupling λ Φ . ∆a µ and BaBar -A discrepancy between the most precise ∆a µ measurement performed by the Muon (g−2) collaboration [43] and theoretical calculations [58][59][60][61][62] stands at more than 3.7σ 2 . In view of the efforts by the Muon (g − 2) collaboration to measure this quantity four times more precisely at FNAL [67], it is timely to reconsider the dark photon solution to the ∆a µ puzzle [40]. Minimal dark photon models are excluded by collider and beam dump searches for Z → + − [68][69][70][71]. If a GeV dark photon decays invisibly, then it is subject to strong limits from monophoton searches at BaBar [44]. This constrains ε 2 10 −6 for m Z < 3 GeV by searching for a are the mixing factors in Z Niνj vertices, and αD = g 2 X /4π. Note that Z → ν3ν3 is negligible for the mixings considered. We refer to the MiniBooNE excess as MB, the BaBar excess as BB, and the accelerator experiments as Acc. The zeroes in BP-A are protected by a left-right symmetry (ΛL = ΛR).
missing-mass resonance produced alongside initial-state radiation (ISR), e + e − → γZ . In models where the Z decays semi-visibly inside the detector, B(Z →vis+ / E) 1, this limit can be relaxed. This was proposed in the context of inelastic DM models in Ref. [72], and later criticized in a more conservative analysis [73] (see also [74][75][76]).
In our model, however, the mechanism put forward in Ref. [72] is improved, as more visible energy is deposited in the detector. For the bound to be relaxed above the central value to explain the ∆a µ anomaly, the detection inefficiency for the Z decay products in ISR events ought to be at most 0.22%. Note that in virtually all ISR monophoton events the produced Z promptly decays into N 5 and/or N 6 states, which subsequently lead to one or more e + e − + / E vertices. Such additional particles are hard to miss in the barrel-like BaBar detector, which operates with a 1.5 T magnetic field. In fact, after produced, all N 6 states decay already inside the drift chamber, while for BP-(A,B), we find that, for a typical 2.5 GeV N 5 energy, (79,92)% of N 5 states decay before the electromagnetic calorimeter (ECAL), followed by (11, 5.7)% inside the ECAL, and (8.0, 2.3)% in the muon detection system. Fully invisible decays are rare and satisfy B(Z → N 4 N 4 ) + B(Z → N 4 N 5 ) × P escape N5 2.2 × 10 −3 for the BPs. Visible decays are also negligible, Pseudo-monophotons at BaBar -The dominant production of dark particles in e + e − colliders is s-channel pair production of HNLs due to the large values of α D ε 2 . In particular, the process where ε B is the final detection and selection efficiency of the monophoton analysis at BaBar, not including the probability P γ N5 that the N 5 states decay inside the ECAL and get reconstructed as a photon. For the ISR analysis, ε ISR 0.2 − 3.5%, depending on M 2 miss . In our pseudo-photon case, however, it is impossible to estimate ε B without a dedicated detector simulation and the machine learning algorithm utilized by BaBar. Nevertheless, we fit our model prediction to data, which will give an estimate of the value of P γ N5 ε B required to explain the excess in the model. Since backgrounds are much larger than our signal above M 2 miss > 50 GeV 2 , our fit uses only the data in 24 GeV 2 < M 2 miss < 50 GeV 2 , where a total of 189 events are observed on top of a prediction of 157 background events. Floating P γ N5 ε B for BP-B, we minimize a binned Poisson likelihood, assigning a 1 (5)% normalization uncertainty on the background model. We find a 2.5σ (2.2σ) preference for 53 signal events. Our best-fit point in BP-B is shown in Fig. 1, where events were selected if θ ee < 10 • , and the boost and azimuthal angle cuts were implemented as in Ref. [44]. This corresponds to a total number of 93 pseudo-monophoton events, before any selection cuts. Finally, since both BPs predict very similar shapes, we can make use of Eq. 5 to find P γ N5 ε B (0.26, 0.10)%. A dedicated analysis at BaBar would be able to determine if such numbers are experimentally justified. We also note that our pseudo-monophoton rate is compatible with constraints on B(Υ(1S) → γ + / E) < 5.6 × 10 −6 at 90% C.L. at BaBar [77], provided the e + e − → γ mis-ID rate is less than (100, 77)% for BP-(A,B).
MiniBooNE excess -MiniBooNE is a mineral oil Cherenkov detector in a predominantly ν µ beam with E ν 800 MeV. Recent results with improved background analysis and larger statistics [78] report an excess of 560.6 ± 119.6 (77.4 ± 28.5) e-like events in ν (ν) mode. Initially designed to search for short-baseline oscillations reported by the LSND experiment [79], MiniBooNE reports a much more significant 4.8σ excess. The large tensions with global datasets in oscillation models [80][81][82] (see also [83][84][85]) prompts new scenarios to explain the excess.
We propose that the MiniBooNE excess arises from the decay products of HNLs produced in ν µ upscattering inside the detector, where H = {C, p + } is the hadronic target. The e + e − pairs with small angular aperture or large energy asymmetry mimic a single EM shower in the Cherenkov detector. This is similar to the upscattering explanation proposed in Ref. [41], but successfully achieves fast HNL decay without infringing upon any bounds 3 . A prediction of our signal on top of MiniBooNE neutrino data is shown in Fig. 2 for our BP-B. In our single generation approximation, the upscattering cross section is proportional to is the mixing factor in the ν 3 N j Z interaction and takes O(10 −7 ) values. The scattering is predominantly electromagnetic via Z exchange, and due to the large values of α D ε 2 , no interference with the SM Z is observed. This, together with the purely vectorial couplings of the Z , explains why the signal prefers to be forward with respect to charge-current quasi-elastic scattering. We note that scattering on protons is dominant, and that the angular spectrum predictions can improve when nuclear effects and higher Q 2 scattering regimes are included. The produced e + e − that contributes to the excess has a small invariant mass, with m ee < m 5,6 − m 4 . If m ee is too large, it contributes to the NC π 0 dataset, where an excess is also observed [95]. We estimate the overall detection and signal selection efficiency for our BPs to be 5%. Although many upscatterings lead to N 6 → (N 5 → N 4 e + e − )e + e − , we do not include these double vertex events as a large fraction of them would be excluded by the MiniBooNE cuts.
Old accelerator anomalies -Many accelerator experiments in the 80s and 90s searched for ν µ → ν e transitions at short-baselines, with some of them observing significant excesses. While a neutrino oscillation interpretation of these results is excluded, they can be explained within our model, where the energy dependence and signal characteristics differ from those of oscillation. The largest deviation was observed by the PS-191 experiment at CERN using a E peak ν ∼ 600 MeV ν µ beam and the fine-grained ECAL component of their detector. They observed an excess of 23 ± 8 e-like events on a background of 12 ± 3 events, amounting to a 3σ significance [45,96]. All excess events contained a scattering vertex, followed by an electromagnetic shower < 16 mm away. A follow-up experiment, E-816 [46], was designed to test the PS-191 anomaly at the Brookhaven National Laboratory (BNL) with a wide-band beam of mean energy E ν 1.5 GeV. E-816 also reported an excess of e-like events with a small vertex-shower separation of < 8.8 mm, although at a lower significance of 2σ due to larger systematic errors [46]. In our model, these excesses can be explained by ν µ upscattering to N 6 , which decays very fast to overlapping or energy-asymmetric e + e − pairs, fitting the exponential drop of events as a function of vertex-shower separation. PS-191 and E-816 observed a larger excess than MiniBooNE, which could be explained by the larger N 6 upscattering rate (BP-B) or solely due to different signal reconstruction (BP-A).
Other experiments with E peak ν 1 GeV reported no excess, namely E-734 [97] and E-776 [98,99]. The stringent cuts against π 0 backgrounds would veto most of our e + e − pairs and weaken the constraint. Another set of bounds come from high energy experiments, such as NO-MAD [100], with E ν 24 GeV, and CCFR [101] and NuTeV [102], both with E ν 140 GeV. Their bounds, although very strong under the oscillation hypothesis, are much weaker for our model due to the log E ν growth of the Z mediated neutrino-nucleus cross-sections in comparison to the linear E ν growth in the SM. Finally, we note an unexplained excess of positron events observed at NOMAD [100] in a sideband sample of events containing showers far from the scattering vertex or that had failed kinematic cuts. Such positrons are predicted in our model as coming from asymmetric e + e − pairs in the late decays of our HNLs.
We also note an intriguing excess reported by CCFR in the search for HNLs produced in scattering [47,48,103,104]. The experiment saw evidence for double-vertex events with 9 NC/NC events over an estimated overlay background of 3±0.2 (stat.)±0.4 (syst.). A double-vertex event was defined as one in which there were "two distinct and separate shower regions", and NC/NC refers to two neutral vertices, as opposed to NC/CC events, wherein a second vertex contained a muon candidate. No excess was observed in the NC/CC, which disfavored standard interpretations with HNLs that have large branching ratios to muons. In our model, only NC/NC events appear, mainly from upscattering into N 6 , which immediately decays into N 5 e + e − , with the subsequent N 5 → N 4 e + e − decays typically happening after a few meters at CCFR energies. This leads to good agreement with the 4 to 14 m vertex-shower separation observed, given the typical N 5 energies of 50 GeV. A naive scaling of the cross sections shows that the normalization is compatible with the rate at MiniBooNE and PS-191.

IV. DISCUSSION AND CONCLUSIONS
Let us remark that our BPs satisfy all existing experimental constraints, including decay-in-flight bounds from PS-191 [96,105]. Searches for peaks in the muon spectrum in π + /K + → µ + N j [106,107] are also satisfied due to strong vetoes against visible energy in the detector, as discussed in Ref. [32]. Intriguingly, the latest results from K + → e + N j searches at NA62 [108] indicate an excess at m N = 346 MeV, with |U ej | 2 1.5 × 10 −9 at 2.2σ (3.6σ) global (local) significance. Our model can accommodate this hint by identifying N 5 with the required HNL and switching on the mixing with the electron neutrinos. To take into account the visible decays of our HNL, the required |U ej | 2 is enhanced by a factor ∼ 2 for ∼ 5 ns lifetimes, as quoted by the experiment. For our BPs, we also expect to see an excesses in the muon sector, depending on the K + → µ + N j efficiency at NA62.
There is some freedom in the choice of the HNL parameters while keeping the same key phenomenological features, e.g. HNL decay length and Z branching ratios. For the dark photon parameters, the situation is more constrained. For instance, lower values of m Z , such as 1 GeV with ε 2 = 3 × 10 −4 are possible, and decrease the required α D |V 3j | 2 couplings to explain MiniBooNE, PS-191, and BaBar by a factor of (1.25) 4 ∼ 2.5. Going much below m Z = 1 GeV leads to more forward angular distribution at MiniBooNE and introduces tension with neutrino-electron scattering constraints [90]. A survey of existing bounds and additional BPs are provided in Appendices B and A.
We also want to highlight the left-right symmetry in BP-A, as in that case the lightest HNL ν 4 has vanishing interactions with the Z , except for the |V 45 | 2 vertex. Incidentally, N 4 could lie at the keV scale, and may be a candidate for non-thermal dark matter [109].
Direct searches for our MiniBooNE explanation can be performed at the Short-Baseline Neutrino program at FNAL [110,111], which comprises three Liquid Argon detectors: SBND, µBooNE, and ICARUS. Specifically, for BP-(A,B) we predict that µBooNE [112] would see a total number of ∼ 760 neutrino upscattering events into N 5 and (0, 2800) events into N 6 , before any efficiencies and for a total N POT = 13.6 × 10 20 . While the former would contain a single e + e − pair, the latter events would constitute double vertex events with 10 cm separation. Around 60% of the total number of events are due to coherent scattering, and leave no visible proton tracks. Dedicated studies of the e + e − invariant mass, as well as searches for the double-vertex events would help discriminate our hypothesis from other dilepton MiniBooNE explanations. Other near detectors to experiments like MINERνA, NOνA, and T2K could also shed light on the model. In particular, the incoherent piece of our prediction may be constrained by photon-like showers in ν e CC quasi-elastic measurements and the coherent piece by photon-like showers at neutrino-electron scattering measurements. Searches for e + e − pairs in decay-in-flight at the ND280, the off-axis near detector of T2K, can also constrain HNLs produced in coherent neutrino upscattering inside the detector [113].
Other direct searches can be performed at the NA62 kaon facility [114]. The decays of 75 GeV/c kaons to K + → + α N i followed by N j → N k e + e − would constitute a background-free signature, similar to the one proposed in Ref. [32]. The new physics events would appear as a displaced e + e − vertex with peaked kinematics, The production rate is controlled by |U µj | 2 , where for BP-(A,B) we predict a total K + → µ + (N 6 → N 5 e + e − ) event rate of (1970,2980) for N K = 2.14 × 10 11 fiducial kaon decays and an overall 4% acceptance [115,116].
The dark photon can be searched for in the ISR events at BaBar, Belle-II [73,117], and BESIII [118] by relaxing the vetoes on additional e + e − pairs in the detector. The large value of ε 2 required for the ∆a µ explanation yields several hundred events at BaBar. Direct N j N k pair production, as well as Higgstrahlung e + e − → ϕ Z , would also appear as multiple displaced e + e − vertices at B-factories, and in the fixed-target experiments NA64 [119,120] and LDMX [121], providing a background-free signature for semi-visible dark photons.
In summary, this letter provides an explanation to some of the most prominent low energy anomalies, including the MiniBooNE excess and the ∆a µ anomaly. The phenomenological signatures we presented are achieved in a renormalizable model which extends the SM by an anomaly-free U (1) gauge symmetry and a dark neutrino sector. The model is able to reproduce the correct scale for the light neutrinos, albeit with some level of fine tuning. Phenomenologically, our scenario only requires a semi-visible GeV-scale dark photon that couples to O(100 MeV) HNLs. We show that the dark photon not only evades sensitive searches for missing mass resonances at BaBar, but can actually explain a mild but continuous excess seen in the data thanks to the pseudo-monophotons from N 5 → N 4 e + e − decays. Due to the large kinetic mixing required by ∆a µ , such events naturally arise from s-channel e + e − collisions producing HNLs. We point out that e-like events from upscattering are better able to explain past anomalies reported by PS-191 and E-816, compared to those from excluded oscillation hypotheses. Also curious is the prediction of O(2 cm) lifetime for N 5 , as it leads to double vertex events at neutrino experiments and is compatible with a significant excess reported by CCFR. The novel interplay between portal couplings and exotic decay signatures in our model offer striking signatures at current and upcoming experiments. Observations of displaced vertices at kaon and neutrino experiments, as well as the decays of a semi-visible dark photon, would provide confirmation of our model. In the main text we have focused only on the phenomenological aspects of our model, giving two BPs that can resolve the low energy anomalies. In this appendix, we offer more details on the model side, giving the vertex factors for each relevant interaction that can be used to compute physical observables. The BPs in the main text were given in terms of a model with a single generation of active neutrinos, n = 3 sterile states and d = 1 dark vector-like fermions. We also present two additional BPs to illustrate the ranges of the HNL masses compatible with the phenomenology discussed. In particular, BP-C indicates the smallest scale of m 6 and m 5 which lead to sufficiently fast N 6 and N 5 decays. With BP-D, we illustrate the features of heavier masses.

ACKNOWLEDGMENTS
Following Eq. 2 in the main text, the full mass matrix is given as, The values for the mass matrix parameters used for our BPs are given in Table III. In the mass basis, HNLs mixing with the different flavors is given in Table IV. To clarify the nature of our neutrino couplings to the neutral bosons, we write the explicit vertices in the neutrino mass basis using the flavor gauge boson basis. To leading order in χ and taking light dark photons: Z µ = Z 0 µ +s W χX µ and Z µ = X µ −s W χZ 0 µ . The interactions are given by whereν m is the mass eigenvector and P L,R = (1 ∓ γ 5 )/2. The vertex factors are defined as We show the relevant vertex factors for dark bosons in our BPs in Table V. For all BPs, we take m Z = 1.25 GeV, m ϕ = 1 GeV and ε 2 = 4.6 × 10 −4 . The mixing sin θ 2 is assumed to be negligible for our BPs. The HNLs with masses above m Z , namely N 7 and N 8 , are mostly in the sterile direction, with |V jk | 2 1, and |U N2j The phenomenology of BP-C is similar to BP-B, with cτ (A4) It is important to note that the bounds from NA62 on both |U e5 | 2 and |U e6 | 2 are weakened due to the fast decays of N 5 and N 6 . For N 5 with lifetimes ∼ 5 ns, the experiment expects a weakening of the bound by a factor ∼ 2 [108] implying an effective |U e5 | 2 NA62 1.5 × 10 −9 . The electron mixing requires two active light neutrinos,ν 2 andν 3 . With our chosen Yukawas,ν 2 is massless and mostly in the ν e direction, withν 3 mostly in the ν µ direction. As we do not consider the full 8 × 8 mass matrix with three active light neutrinos, we do not attempt to reproduce the structure of the PMNS matrix, but note that this can be achieved with appropriate choices of the Yukawa couplings in the active sub-block. The scattering cross-section at MiniBooNE is now proportional to i=2,3 |U µi V ij | 2 α D (eε) 2 |U µ3 | 2 |V 3j | 2 α D (eε) 2 , since the |V 2j | mixings are negligible for masslessν 2 .
Appendix B: Survey of Existing Constraints a. Electroweak precision observables An assessment of the impact of kinetic mixing on electroweak precision observables (EWPO) requires a global fit to collider and low energy data. This was performed in Ref. [122], where a model independent bound on ε was derived. For m Z M Z , the authors find ε 2 EWPO < 7.3 × 10 −4 at 95% C.L, just above our value of ε 2 = 4.6×10 −4 . As a sanity check against more recent data, we also directly compute the oblique parameters S, T , and U [123] to leading order in ε = c W χ and µ ≡ g X v ϕ /M SM Z , neglecting the impact of running in the dark couplings and corrections from dark fermion loops. For all our BPs, these are [124][125][126] S 4s 2 W ε 2 (1 + µ 2 )/α = 0.042, Clearly, this is compatible with the current bounds of T < 0.22 and S < 0.14 at 95% C.L. [57]. The constraints on S can be much stronger when fixing T = 0 and U = 0, which would be mostly driven by the 1 to 2σ discrepancies observed between direct M W measurements and the global best fit point. We plan to return to this issue in future communication [56]. b. Deep-inelastic scattering constraints Recently, Ref. [127] appeared setting new model-independent constraints on dark photons using ep + scattering data from HERA [128]. At 95% C.L., the authors find that ε 2 2.9 × 10 −4 , in tension with our BPs. As the authors discuss, inclusion of other datasets weakens the bound, which signals a mild tension between HERA and other experiments. Finally, we note that a naive rescaling of the constraints on contact interactions performed by the ZEUS collaboration [129], where the probability distribution functions were included in the fit, leads to bounds that are weaker by a factor of ∼ 2 than the ones quoted by Ref. [127]. While HERA is certainly sensitive to our model at some level, we believe that it is not excluding it at the 95% C.L. Nevertheless, a trivial modification to our setup to accommodate such bound is to lower m Z = 1 GeV.
c. Z → invisible Dark fermions can be produced in the decays of SM-like Z bosons via its couplings to the dark current. This is induced by kinetic mixing, and to leading order in χ it is This coupling is relevant in our model since g X ε is not so small. A constraint can be derived from LEP measurements of the Z boson decay width [130] and constrains Γ Z→inv < 2 MeV. The largest new physics decay mode is Z → N j N k , for j, k > 3, which even without requiring Here, the Vij ≡ U * D L i UD L j − U * D R i UD R j are the mixing factors in Z Niνj vertices, and αD = g 2 X /4π. Note that Z → ν3ν3 is negligible for the mixings considered here. We refer to the MiniBooNe excess as MB, the BaBar excess as BB, and the accelerator experiments as Acc. The zeroes in BP-A are protected by a left-right symmetry (ΛL = ΛR).
safely below the current constraints even for the largest α D couplings, as it can be shown that 8 j,k |V jk | 2 = 2. Another relevant process is Z → Z ϕ . Neglecting the final state masses, we find also satisfying the constraints independently of the fate of ϕ and Z in the detector. d. h → invisible Searches for Higgs decays to invisible have been performed by CMS [131] and AT-LAS [132]. Latest preliminary results by ATLAS require that B(h →invisible) < 0.13 at 95% C.L. [133], which for which does not lead to strong constraints on our model given that the largest Yukawa we have is for BP-D, where it is a few 10 −5 . The parameters in the scalar potential are also subject to constraints. We consider h → ϕ ϕ induced by the scalar portal coupling λ ΦH . A direct computation neglecting final state masses leads to, which can be interpreted as a strong constraint on our scalar mixing angle θ . This value is currently an order of magnitude below current sensitivity of K L → π 0 ϕ searches at KOTO or K + → π + ϕ searches at NA62.
Decay to a pair of HNLs can also proceed via the dark Yukawas if scalar mixing is present. Neglecting the final state masses, we find where s 2 θ ≡ sin 2 θ is chosen according to B8 for v ϕ = 500 MeV. This clearly satisfies the bound as it can be shown that 8 j,k |S j,k | 2 = 4(|Y L | 2 + |Y R | 2 ). Similarly, the decay h → Z Z is possible and may appear invisible some fraction of the time. Nevertheless, the tree-level rate is and loop-corrections from fermion loops are also negligible.
K → πϕ The KOTO experiment at J-PARC [134] has set stringent constraints on B(K L → π 0 / E). Initial hints of a signal in the latest unblinding [135] have later been revisited due to unexpected charged kaon backgrounds and signal mis-identification [136,137]. The hinted values led to branching ratios much larger than the SM prediction [138], and prompted several new physics studies [35,[139][140][141][142]. In light of the new backgrounds and given stringent constraints on the scalar mixing found above, we refrain from trying to explain these events. We note that KOTO, as well as NA62, are only sensitive to singlet scalar emission in K → πϕ decays for mixings of order s 2 θ ∼ 10 −7 , which for our small v Φ are already excluded due to h → ϕ ϕ decays. A more exotic scalar sector could be invoked to explain this excess, where a new invisible real scalar S could be hidden under background at NA62 if m S ∼ m π [143], and could lead to signatures at KOTO without violating the Grossman-Nir constraints [144].
e. Meson → invisible We consider the decays of vector meson states due to the vector nature of the dark photon couplings. The best current bounds are at the level of B(J/ψ → inv)< 7.2×10 −4 at BES [145], and B(Υ(1S) → inv)< 3.0 × 10 −4 at BaBar [146]. The branching ratios into HNLs, B(V → N j N k ), may be still be slightly above such values, provided a sufficient number of the produced HNLs decay semi-visibily. In general, the branching ratio for the quarkonium states (V) used throughout our article is where we neglected the final state masses, and took m Z = 1.25 GeV and ε 2 = 4.6 × 10 −4 . The decay constants, f Υ(nS) = 498, 430, 336 MeV for n = 2, 3, 4, were extracted from existing V → e + e − measurements in Ref. [147]. For the fully invisible states N 4 N 4 (as well as for light neutrinos) the mixing factor V jk is sufficiently small to avoid the constraints. Production of N 4 N 5 is the next largest contribution and it still satisfies the most stringent limits from BES, since in all BPs is the probability for N 5 to escape detection. f. Pseudo-monophotons As discussed in the main text, s-channel e + e − collisions at BaBar can lead to pseudo-monophoton events. It is not possible to extract a constraint from this without a dedicated detector simulation, as it relies on experimental details such as the efficiency to reconstruct our e + e − as a photon, and on the specifics of the machine learning algorithm. In the main text, however, we proceeded to understand if it is at all feasible to explain a mild excess observed in the monophoton data. By finding a best-fit value for the normalization of events that are "photon-like", we have asked whether such rate is possible within our model and whether the efficiencies it requires are reasonable. Here, photon-like refers to events where a N 4 N 5 pair is produced in the interaction point, followed by a N 5 → N 4 e + e − decay inside the ECAL. When plotting our prediction, we required that the angle of separation between the electrons be less than θ ee < 10 • , and select events within the angular acceptance of the ECAL detector. For our total pair-production rate, we include HNLs produced in Z mediated s-channel e + e − collisions, where the e + e − → (Z ) * → N i N j cross section was found to be neglecting final state masses and where θ CM is the center of mass angle between N i and the collision axis. We also include a contribution from e + e − → Υ(nS) (n = 2, 3, 4), followed by decay into HNLs. We use σ ee→Υ(2S) ( √ s = 10.02 GeV) 7 nb and σ ee→Υ(3S) ( √ s = 10.36 GeV) 4 nb as well as B11. g. Higgstrahlung Another source of HNL production at e + e − colliders is dark higgstrahlung. Due to the large dark coupling, the process e + e − → ϕ Z is important [76], with a differential cross section , (B13) where θ CM is the center of mass angle between the Z and the collision axis. This process is therefore comparable to direct HNL production. Remaining agnostic about the decay products of ϕ , but requiring the decay Z → N j N k , the ratio between direct and higgstrahlung HNL production in our BPs is For production of N 4 N 5 pairs in ours BPs A and B, this constitutes a ratio of just above 3. Given that ϕ decays promplty and visibly, especially into N 6 states, we do not include this contribution in our monophoton discussion, but emphasize that this offers yet more visible signatures at e + e − colliders.
h. Υ(1S) → invisible + γ Vector meson decays to photon plus missing energy are direct probes of our pseudo-monophoton events. The full process is Υ(2S) → π + π − (Υ(1S) → γ + / E), where the π + π − kinematics can be used to identify the Υ(1S) state. The current limits are quoted in terms of the BR into an invisible pseudoscalar, Υ(1S) → γ + A 0 , and into a pair of invisible fermions, Υ(1S) → γχχ. Most relevant to us are the three-body decay limits taken at the smallest χ masses (m χ → 0), where BaBar [77] constrains which was improved by Belle [148] to all at the 90% C.L. More recently BESIII [149] has set the strongest limits on the two-body process but since the missing mass in this process is fixed, the constraint does not apply to us. When implementing these bounds on our model, we used the BRs in B11. For comparing the Υ(1S) → γ + / E constraints to the BaBar pseudo-monophoton rate, only the BaBar limit is taken into account, as P γ N5 is a detector-dependent quantity, and, under the simplifying assumption that it is constant in energy and angle, it is the same for the two processes.
Since the efficiencies for our pseudo-monophoton events are different than in the s-channel production mode, we use the limits above to obtain an upper-bound on the detector-dependent quantity P γ N5 . Neglecting the sub-dominant N 5 N 5 contribution, BaBar sets a limit of B(Υ(1S) → N 4 N 5 )×P γ N5 < 5.6×10 −6 at 90% C.L, which implies P γ N5 < (18, 10, 1.5, 30)% for BP-(A,B,C,D). Since the probability for N 5 to decay inside the ECAL is known to be (14, 13, 14, 10)% for BP-(A,B,C,D) (E N5 5 GeV for s-channel production), we use the limit on P γ N5 to find the largest allowed e + e − → γ mis-ID rate for our explanation of the BaBar excess not to be excluded. Under the approximation that it is independent of the kinematics, this rate is bounded from above by (100, 77, 12, 100)%. Clearly this allows to explain the monophoton excess while remaining consistent with the Υ(1S) → γ / E. i. NA64 searches The fixed-target NA64 [119,120] experiment also sets stringent limits on invisible dark photons. The 100 GeV electrons can produce dark photons via bremsstrahlung in interactions with the dense beam dump material, eW → eWZ , with W a Tungsten nucleus. The bound from 2.84×10 20 electrons on target is shown in Ref. [119] for invisible dark photons with masses as large as m Z ∼ 0.94 GeV. For the latter dark photon mass, it constrains ε 2 10 −4 at 90% C.L. As discussed in Ref. [120], the semi-visible decays of the dark photon can weaken the bound. For our BPs (m Z = 1.25 GeV), we do not have access to the exact value of the invisible-Z constraint, but under a conservative assumption of linear scaling with m Z , we find that NA64 does not constrain our BPs provided ∼ 45% of Z particles produced are vetoed due to the subsequent semi-visible decays of N 5,6 . For the HNLs produced in the decay of the highest energy dark photons (E N5,6 ∼ 50 GeV), we have a typical decay length cτ of O(10 m) and O(1 m) for N 5 and N 6 , respectively. Note that the dark photons are produced inside the beam dump and decay promptly. Under these conservative assumptions, we find that most HNLs decay before or within the instrumented ECAL of NA64, assumed to be a total of ∼ 10 m. The presence of one or more vertices of e + e − pairs would be vetoed from the invisible-Z search due to visible showers in the ECAL, as well as in the additional veto detectors.
j. Beam dump and decay-in-flight searches HNLs heavier than N 4 in our model are unconstrained by decayin-flight searches due to their short lifetimes. On the other hand, N 4 is longer-lived and faces strong constraints from HNL searches at PS-191 [96,105]. If N 4 has new interactions, such as in BP-B, it decays faster than in the minimal HNL models, and the constraints from decays in-flight are modified (see, e.g., the discussion in Ref. [150]). While N 4 is produced in π, K → µN 4 decays, which are controlled by |U µ4 |, its subsequent N 4 → νe + e − decays in BP-B proceed mainly through Z exchange, which is controlled by |V 43 |. In that case, we require where |U * µ4 U e4 | PS is the bound quoted by the PS-191 experiment. The factor F = 3.17 converts the bound on Dirac to Majorana HNLs, and takes into account that PS-191 assumed only charged-current decays in their analysis.
We note that there are additional production channels for N 4 than in the minimal HNL models. On top of the standard meson decays π, K → N , HNLs can also be produced via kinetic mixing in ρ, ω → N i N j and π 0 , η, η → γN i N j , where the vector meson decays dominate. These channels have been explored in the context of two-component fermionic dark sectors in Refs. [151][152][153][154]. In this context, limits on new dark sector fermions that decay to e + e − + / E have been set using CHARM [155] and NuCal [156] data. With the effective field theory approach of Ref. [154], we see that for cτ Ni 10 cm or smaller such constraints are safely avoided due to the short lifetimes. For the longer-lived N 4 , however, a simplified re-scaling of the constraints from Ref. [154], where g 2 /Λ 4 → |V 45 ||U µ4 |G F (eεg X /2m 2 Z ), and g 2 /Λ 4 → |V 45 ||V 43 |(eεg X /2m 2 Z ) 2 for BPs B and C, shows that the constraints are satisfied. A re-analysis of the PS-191 constraints including these new production mechanisms for N 4 , both from vector meson decays as well as from the secondary decays of N 5 , N 6 could set stronger constraints in our parameter space, but is beyond the scope of this work.
We would like to highlight an event found in PS-191 and shown in Fig. 9 of Ref. [157]: it has two tracks in the initial decay detector which subsequently shower in the ECAL. While it is unlikely to be due to photons, as they would not be recorded in the flash tubes, it could be due to two electrons coming from an N 4 decay. We do not elaborate further on this intriguing event.
k. Peak searches -Searches for a missing mass in π/K → µN j decays set stringent limits on |U µ4 | 2 . In our model, the N 5 and N 6 states can be produced, but would lead to visible signatures inside the two most relevant experiments, namely E949 [106] and NA62 [107,108]. As pointed out in Ref. [32], the small probability to miss additional energy deposition in these experiments together with the stringent vetoes against K → µνγ ( * ) backgrounds would, in fact, veto most of our events. For E949, the detection inefficiency is estimated to be larger than 0.5%, the typical photon inefficiency, and so the constraints would be weakened by factors of 200. A similar argument can be made about NA62, where the e + e − signatures would have to be missed by several different detector components. In this way, our BPs are not excluded by peak searches, in particular BP-B where we rely on a relaxation by a factor of ∼ 10 on the E949 limit on |U µ6 | 2 . Note that the mass of both N 4 and N 5 is always below the mass interval constrained by both E949 and NA62. This is important as N 5 has a larger probability to escape these detectors due to its O(10 m) decay lengths in the laboratory frame.
l. Lepton number violating searches Due to electron mixing, we predict large rates for neutrinoless doublebeta decay in BP-D. In addition to the light neutrinos, heavy N j states also contribute and may dominate. Their contribution contains large uncertainties as the momentum dependence of the nuclear matrix element is important at the O(100) MeV mass scale. Nevertheless, a naive rescaling of the effective mass leads to m ββ ≈ 130 meV for p 2 = (100 MeV) 2 when all matrix elements U ei are real. This is to be compared with the current experimental sensitivity of m ββ < 61 − 165 meV at KamLAND-Zen [158]. We interpret this as a suggestion that, unless strong cancellations due to the Majorana phases are at play, we predict observable rates of neutrinoless double beta within current experimental reach. Loop contributions [51] and a full three active flavor treatment of the mixing matrix are also important and should be studied in more detail. Lepton number violating kaon decays of the type K + → µ + (N j → µ + π − ) are important for generic heavy Majorana neutrinos [13,115,159]. This signature proceeds via charged-current branching ratios of N j , and it is much suppressed in our model whenever j > 4, where B(N j → µ + π − ) ∼ 10 −8 − 10 −7 . For N 4 , such decays can in principle have large branching ratios, but the longlifetimes of N 4 renders the experimental searches insensitive.

Appendix C: Old Accelerators Experiments -Additional Details
We now provide additional details regarding the accelerator experiments.
a. PS-191 The PS-191 detector was made of 5 × 5 mm 2 flash tube chambers interleaved with 3 mm thick iron plates giving the detector a very fine granularity, 3 mm of iron or ∼ 17% of a radiation length. This was used to distinguish photons, whose showers started further from the vertex due to conversions, from electrons, which showered immediately. As shown in Fig. 3 of Ref. [45], most single-shower events started within the first chamber, which corresponded to ∼ 16 mm. The analysis was restricted to events above 400 MeV, to avoid π 0 backgrounds. The initial e-like shower sample contained a total 57 events, which, after cuts on energy and distance between vertex and shower start, left a pion background of 7 ± 3 events.
In our model, the most frequent upscattering events at PS-191 would produce a N 6 , which immediately decays into N 5 e + e − . The electromagnetic (EM) shower created by this decay is then close to the upscattering vertex, and would explain the sharp drop in number of events as a function of the shower-vertex distance. The issue of normalization between the number of events required at PS-191 versus those observed at MiniBooNE, which is between a factor 5 to 7 larger, can be explained by the fact that not all N 6 events in MiniBooNE count as signal. This is due to the additional decay of N 5 , which yields a total of four charged leptons that are very rarely mis-reconstructed as a single EM shower. Therefore, by increasing |V 63 | 2 in comparison to |V 53 | 2 , one can increase the ratio of signal events between PS-191 and Mini-BooNE. This happens in BP-B, where |V 63 | 2 4.6|V 53 | 2 , although it is forbidden in BP-A due to the left-right symmetry. It should be noted that the coherent cross section in Iron is larger than in Carbon, and that the most energy-asymmetric e + e − pairs may be reconstructed as a one track plus one shower events.
b. E-816 The E-816 experiment used the same finegrained ECAL as PS-191. The number of events was quoted as a function of the scattering vertex and the start of the shower, allowing for e/γ differentiation. This was measured in units of 0.25 radiation lengths (∼ 17.6 mm). Any photons converting before this would fake electrons, although the exponential nature of the conversion makes such events unlikely. The experiment searched for ν e excesses in the one track one shower (1T1S) sample. To reduce the π 0 background, showers with E 300 MeV were cut from the analysis. According to their simulations, such cuts eliminated ∼ 70% of π 0 s while only removing ∼ 10% of ν e s. After cuts, the π 0 background dropped to ∼ 1.6% of the ν µ interactions and was of the order of the ν e contamination in the beam. The electron excess was then given by the subtraction of the 1T1S events due to pions and those due to intrinsic ν e background from the remaining 1T1S events. They found an excess of 43 ± 17.8 (stat.) ± 9 (sys.) and quoted a significance of 2.4 ± 0.5 σ.
Similar to PS-191, E-816 would also count upscattering events into N 6 as signal when the e + e − s are overlapping or highly energy-asymmetric. The ratio of ν e -like events to ν µ -like events, R = (ν e + ν e )/(ν µ + ν µ ), observed at E-816, R observed /R expected = 1.6 ± 0.9, is compatible but somewhat smaller than the one at PS-191, R observed /R expected = (2 ± 0.5)/(0.7 ± 0.2). The collaboration attributed this to unknown systematic errors in both experiments.
c. E-734 E-734 at Brookhaven National Laboratory ran with peak energy E peak ν ∼ 1.3 GeV at a baseline of ∼ 96 m, and searched for ν µ → ν e transitions [97]. The experiment utilized a filter program to remove events not containing a single electromagnetic shower within an angular interval θ e < 240 mrad relative to the beam direction, with the remaining events scanned by physicists to remove events with more than one shower or additional hadronic activity. It is interesting to note those events with one shower and an associated upstream vertex were used as a control sample of photons. After a cut on the energy, 0.21 < E e ≤ 5.1 GeV, 873 shower events remained. The main backgrounds were identified to be pion production in NC interactions, charged pion production in inelastic CC processes, and those from ν µ − e scattering. Of particular relevance is their cut on the shower energy of E e < 0.9 GeV, reducing the sample to 653 events. The final sample contained 418 events in the energy range 0.9 < E e ≤ 5.1 GeV.
While the experiment saw no excess, we note that most events in our model would not have passed the more stringent cuts. This is mainly due to the larger energies required by the experiment, but also due to the cuts in energy loss, dE/dx, of the shower. Our events would most likely resemble those of the upstream photon-like sample.
d. E-776 E-776, running with both a narrow-(NBB) [98] and wide-(WBB) [99] band beam of mean energy 1.4 GeV, searched for ν e appearance 1 km from the target. A fine grained ECAL consisting of 90 planes of proportional drift tubes interleaved with 1 in. (∼ 0.25 the radiation length) thick concrete absorbers was utilized. A total of 12.8 × 10 3 events were in the full sample, and 1496 shower events were selected in a scan of the sample. After cuts, which included a requirement that E e > 600 MeV, only 55 events remained. Further cuts on EM shower identification were made, e.g. the number of hits in a cluster and the length of the shower. This left a sample of 38 events. To eliminate the π 0 background, the differences in shower profile of pions and electrons were accounted for -the former being wider and more asymmetric. This cut was quoted to have an efficiency of ∼ 80% for rejecting π 0 s at 1 GeV. The final sample contained 17 electron shower-like events, with the remaining 21 constituting the π 0 s. Accounting for the probability of pion-electron mis-ID gave 9.6 events. The observed 17 events was consistent with the background prediction of 18 ± 4.3 (stat.) ± 3.9 (sys.) events (9.6 ± 3.8 (sys.) from π 0 s and 8.8 ± 1.1 (sys.) from ν e s in the beam), and no excess was reported by the experiment.
Due to the cuts on energy and, in particular, shower profile, a large number of our events would be removed in the analysis, weakening the constraint on our model. e. CCFR CCFR searched for production of HNLs with a magnetized toroidal spectrometer-calorimeter, and studied double vertex events from ν − N interactions. The sample at CCFR was selected using a neutral current (NC) trigger, whose threshold for energy deposition in the calorimeter was 10 GeV. To make sure the primary showers were indeed from an NC vertex, it was required that no muons penetrated past the end of the showers. Subsequent cuts selected events with a secondary shower downstream of the first, and further cuts based on kinematical considerations ensured the primary and secondary showers were separated by an angle relative to the beam of < 100 mrad. The remaining events were categorized by those containing two neutral current vertices (NC/NC), and those with a neutral current vertex followed by a charged current vertex (NC/CC), with the latter being accompanied by a visible muon track in the secondary vertex. For the NC/CC events, cuts on distance between the vertices were made and events with separation > 4λ I were selected, where λ I is the nuclear interaction length and λ I ∼ 16.8 cm in Iron. This left 31 events, which was consistent with the estimated background of 36.8 ± 1.7 ± 3.4. For NC/NC events, the backgrounds depended on the shower separation. For long separations, 14λ I , the major background was due to random overlay events, events in which independent neutrino interactions appeared correlated. There were negligible contributions from neutral hadron punch-throughs, events in which hadrons created in the initial interaction were able to "punch" through to the end of the shower and interact further downstream. In this region, 9 events were seen on a background of 3.0 ± 0.2 ± 0.4. It is noteworthy that the distributions of some kinematical variables, e.g. hadronic shower energy, for the excess were consistent with those from the overlay background. For those events with separation < 14λ I , the dominant background was from neutral hadron punch-throughs produced in the initial nuclear interaction. Studies of this background [104] suggested there were large degrees of uncertainty, preventing the collaboration presenting results of any excess in this region.
Our model provides an explanation of the excess observed in the NC/NC sample with the prompt decays of N 6 → N 5 e + e − , where the N 5 travels several meters before decaying to give the secondary shower. At CCFR energies, 50% of N 5 s decay within 10 m. It is also possible to produce the N 5 s directly in upscattering, giving a signature similar to the above. Alternatively, the first shower can be entirely the hadronic shower from the neutrino interaction vertex with the N 6 surviving long enough to decay at the secondary vertex, this is less likely as most N 6 will decay within 50 cm for our BPs.