Two-Higgs doublet solution to the LSND, MiniBooNE and muon $g-2$ anomalies

We show that one of the simplest extensions of the Standard Model, the addition of a second Higgs doublet, when combined with a dark sector singlet scalar, allows us to: $i)$ explain the long-standing anomalies in the Liquid Scintillator Neutrino Detector (LSND) and MiniBooNE (MB) while maintaining compatibility with the null result from KARMEN, $ii)$ obtain, in the process, a portal to the dark sector, and $iii)$ comfortably account for the observed value of the muon $g-2$. Three singlet neutrinos allow for an understanding of observed neutrino mass-squared differences via a Type I seesaw, with two of the lighter states participating in the interaction in both LSND and MB. We obtain very good fits to energy and angular distributions in both experiments. We explain features of the solution presented here and discuss the constraints that our model must satisfy. We also mention prospects for future tests of its particle content.


I. INTRODUCTION
The rock-like stability of the Standard Model (SM) [1] has provided a powerfully reliable framework for both theoretical and experimental progress in particle physics over many decades. It cannot, however be denied that over this period, the consistent agreement of experimental data (in particular from colliders) with its predictions has also been a source of frustration. This is especially so since there are undeniably strong qualitative reasons, coupled with physical evidence, to expect that there must be physics beyond the ambit of the SM. These reasons, and this evidence, include a) Dark Matter (DM) [2][3][4][5][6], b) the observed matter and anti-matter asymmetry in our Universe [1,7,8], c) the existence of small but non-zero neutrino mass differences [9][10][11][12], with masses widely different in magnitude from those of the charged leptons and quarks, and d) the existence, unsupported by compelling physical reasons, of three families of quarks and leptons with mixings and a large mass hierarchy.
Parallelly, albeit on a relatively smaller scale, extremely important experimental efforts in non-collider settings have supplemented and buttressed the search for new physics. It has gradually become evident that the landscape here is less bleak, and at present one can point to several experiments which report statistically significant discrepancies with respect to the predictions of the SM. Some anomalous results which have garnered attention and spurred significant activity in an effort to understand their origin are: a) excesses in electron events in short baseline neutrino experiments, which are now in tension with muon neutrino disappearance data [13] if interpreted as oscillation effects involving a sterile neutrino; b) observed discrepancies in the values of the anomalous magnetic moment of the muon [14] and the electron [15]; c) a significant excess in the signal versus background expectation in the KOTO experiment [16] which searches for the decay of a neutral kaon to a neutral pion and a neutrino pair; d) discrepancies with SM predictions in observables related to B-decays [17]; and finally, e) anomalies in the decay of excited states of Beryllium [18].
Our focus in this work is on a subset of results in category a) above. Specifically, we address the Liquid Scintillator Neutrino Detector (LSND) excess (e.g., Ref. [19]) and the MiniBooNE (MB) Low Energy Excess (LEE) (e.g., Ref. [20]). In addition to having appreciable statistical significance, they have withstood scrutiny by both theoretical and experimental communities over a period of time. It is thus possible that these results in particular indicate genuine pointers to new physics, as opposed to un-understood backgrounds or detector-specific effects 1 . The solution proposed here also helps resolve the discrepancy between the measured (see, e.g., Ref. [21]) and theoretically predicted (e.g., Ref. [22]) values of the anomalous magnetic moment of the muon.
We show that one of the simplest possible extensions of the SM, the addition of a second Higgs doublet, when acting as a portal to the dark sector, connects and provides an understanding of all three discrepant results mentioned above. Its function as a portal is achieved via its mixing with a dark (i .e., SM singlet) scalar. This mixing in the scalar sector allows heavier dark neutrinos coupled to the singlet scalar become part of the link between the SM and the dark sector. The dark neutrinos play two additional roles: a) they participate in the interaction that we use to explain the excess events in LSND and MB; b) they help generate neutrino masses 1 With regard to LSND and MB, which share many similarities in their overall physics goals and parameter reach (e.g., the ratio of oscillation length versus energy of the neutrino beam), we note that such attribution to their results requires two distinct "mundane" explanations, given that they differ very significantly in backgrounds and systematic errors. via a seesaw mechanism. This lends synergy and economy to the model, the specifics of which we give below. It provides excellent fits to both energy and angular event distributions at LSND and MB. Our paper is organized as follows: Section II briefly gives the specifics of the MB and LSND anomalies and has a brief discussion of the observed discrepancy in the value of the muon g − 2. Section III describes i) the Lagrangian of our model and its particle content, ii) how the couplings of the additional scalars to fermions arise, and iii) the generation of neutrino masses. Section IV focusses on the interaction we use to explain the MB and LSND excesses. Section V gives our results and provides an accompanying discussion of their important features. Section VI discusses the constraints on our model. Section VII provides a concluding discussion, and indicates possible future tests of the model.

II. THE MB, LSND AND MUON g − 2 ANOMALIES
A. Event excesses in MB and LSND Two low energy neutrino experiments, MB (see [20] and references therein) and LSND (see [19], and references therein), have observed electron-like event excesses. Over time, it has become evident that the results of both cannot easily be explained within the ambit of the SM.
MB, based at Fermilab, uses muon neutrino and antineutrino beams produced by 8 GeV protons impinging upon a beryllium target. The neutrino fluxes peak at around 600 MeV (ν µ ) and around 400 MeV (ν µ ). The detector consists of a 40-foot diameter sphere containing 818 tons of pure mineral oil (CH 2 ) and is located 541 m from the target. Starting in 2002, the MB experiment has up to 2019 collected a total of 11.27 × 10 20 Protons on Target (POT) in anti-neutrino mode and 18.75 × 10 20 POT in neutrino mode. Electron-like event excesses of 560.6 ± 119.6 in the neutrino mode, and 79.3 ± 28.6 in the anti-neutrino mode, with an overall significance of 4.8σ have been established in the neutrino energy range 200 MeV< E QE ν < 1250 MeV. Most of the excess is confined to the range 100 MeV < E vis < 700 MeV in visible energy, with a somewhat forward angular distribution, and is referred to as the MB LEE. We note a) that all major backgrounds are constrained by in-situ measurements, and b) that MB, due to being a mineral oil Cerenkov light detector, cannot distinguish photons from electrons in the final state. Additionally, under certain conditions, MB could also mis-identify an e + e − pair as a single electron or positron.
LSND was a detector with 167 tons of mineral oil, doped with scintillator. It employed neutrino and antineutrino beams originating from π − DIF as well as µ decay-at-rest (DAR). The principal detection interaction was the inverse beta decay process,ν e + p → e + + n. The detector looked for Cherenkov and scintillation light associated with the e + and the correlated and delayed scintillation light from the neutron capture on hydrogen, producing a 2.2 MeV γ. The experiment observed 87.9 ± 22.4 ± 6.0 such events above expectations at a significance of 3.8σ, over its run span from 1993 to 1998 at the Los Alamos Accelerator National Laboratory. For reasons similar to those at MB, LSND lacked the capability to discriminate a photon signal from those of an e + , e − or an e + e − pair.
In addition, we mention the KARMEN experiment [23], which, like LSND and MB, employed a mineral oil detection medium, but was less than a third of the size of LSND. It had similar incoming proton energy and efficiencies. Unlike LSND, it saw no evidence of an excess.
There have been numerous attempts to understand both of these excesses. A widely discussed resolution involves the presence of sterile neutrinos with mass-squared values of ∼ 1 − 10 eV 2 , oscillating to SM neutrinos, leading toν e and ν e appearance [24]. It is partially supported by deficits in ν e events in radioactive source experiments and inν e reactor flux measurements as well as results from the reactor experiments. However, this explanation for LSND and MB excesses has had to contend with gradually increasing tension with disappearance experiments and is also disfavoured by cosmological data. For recent global analyses, a full set of references and more detailed discussions of these issues, the reader is referred to Refs. [25][26][27][28][29][30][31].
The tightening of constraints and parameter space for the sterile-active hypothesis has, in turn, led to a large number of proposals to explain one or both of the LSND and MB excesses via new physics [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49]. Many of these scenarios also face a significant number of constraints. For a discussion of these and for related references, we refer the reader to Refs. [50][51][52][53]. It is, however, fair to say that at the present time, the search for a compelling and simultaneous explanation of both the LSND and MB anomalies remains a challenge [54].
B. The muon g − 2 anomaly The Lande g factor, and its deviation from the tree level value of 2, represents one of the most precisely measured quantities in the SM. It thus is also an excellent probe for new physics. Currently there exists a longstanding and statistically significant discrepancy between its measurement [21,55] and the theoretically predicted value, which involves contributions from quantum electrodynamics, quantum chromodynamics and electroweak theory [14,22,56,57]. Specifically, ∆a µ = a meas µ − a theory µ = (2.74 ± 0.73) × 10 −9 . (1) There have been many proposals for new physics which provide possible explanations for this discrepancy (for reviews and a full list of references, see [14,22,56,57].). Our attempt in this work, details of which are provided in the sections to follow, is related to a class of possible solutions suggested by several authors [58][59][60][61][62][63][64][65][66][67][68][69] involving a light scalar with a mass in the sub-GeV range and a relatively weak coupling to muons.

III. THE MODEL
We extend the scalar sector of the SM by incorporating a second Higgs doublet, i .e., the widely studied two Higgs doublet model (2HDM) [70,71] in addition to a dark singlet real scalar 2 φ h . In addition, three right-handed neutrinos help generate neutrino masses via the seesaw mechanism and participate in the interaction described in the next section.
Our model resembles the approach taken in Refs. [47,48]. In particular, it is essentially a more economical version of the model in Ref. [48], without an additional U (1).
are the mass eigenstates, ii) H 0 1 ≈ h is the SM-like Higgs in the alignment limit (i .e., λ 6 ∼ 0) assumed here, and iii) m 2 h λ 1 v 2 . The masses of the extra CP-even physical Higgs states (H, h ) are given by Also, the charged and CP-odd Higgs masses, respectively, are given by In the Higgs basis the relevant Lagrangian L can be written as follows where and Y k ,Ỹ k are the Yukawa couplings in the (Φ 1 , Φ 2 ) basis. We note that X k ij andX k ij are independent Yukawa matrices. The fermion masses receive contributions only from X k ij , since in the Higgs basis only φ h acquires a non-zero VEV while φ H = 0 = φ h , leading to X k = M k /v, where M k are the fermion mass matrices. In this basis,X k ij are free parameters and non-diagonal matrices. Hereafter, we work in a basis in which the fermion (leptons and quarks) mass matrices are real and diagonal, where After rotation, one finds the following coupling strengths of the scalars h, h and H with fermions (leptons and quarks), respectively: where m f are the SM fermion masses and y f are the diagonal elements of the rotatedX f which are independent of the Yukawa couplings (y h f = m f /v) of the SM Higgsfermions interactions. δ manifestly becomes the scalar mixing angle between the mass eigenstates (H, h ) and the gauge eigenstates (H 0 2 , H 0 3 ). For neutrinos, we define n Ri = (U ν R ) ij ν Rj and n Li = (U ν L ) ij ν Lj such that the matrices M ν (= vX ν ) and m can be diagonalized as follows.
One can then define the following matrices: The part of the Lagrangian describing neutrino masses and interactions is then given by Consequently the neutrino mass Lagrangian becomes (19) and the neutrino mass matrix is given by The neutrino mass matrix m ν can be diagonalised, up to O(m Di /m ν R i ), by the neutrino mixing matrix N which can be written, up to corrections of O(m 2 Di /m 2 where Θ i = m Di /m ν R i . The neutrino mass eigenstates (physical states) are given by For the normal order (m ν1 < m ν2 < m ν3 ), the two mass squared differences of the light neutrinos determined from the oscillation data are ∆m 2 21 = (7.05 − 8.14) × 10 −5 eV 2 and ∆m 2 31 = (2.41 − 2.60) × 10 −3 eV 2 [74]. We have chosen a benchmark point, see Table I, so that it satisfies these values. Finally, the part of the Lagrangian specifying neutrino interactions is given by where the coupling strengths of the scalars h , H for vertices connecting active and sterile neutrinos, respectively, are as follows.
Additionally, the coupling strengths of the scalars h , H for vertices connecting two sterile states, respectively, are   Finally, we stress that all the Yukawa couplings of the light scalars (h , H)-fermion interactions (y h ,H f , f = , q, ν , N ) are free and independent of the Yukawa couplings (y h f = m f /v) of the SM Higgs-fermion interactions.

IV. THE INTERACTION IN MB AND LSND
In the process shown in Fig. 1, the heavy sterile neutrino N 2 is produced via the upscattering of a muon neutrino (ν µ = U µi ν i ) present in the beam, both for MB and LSND. Once N 2 is produced, it decays promptly to another lighter sterile neutrino N 1 and a light scalar h . In our scenario, N 1 is a long-lived particle that either escapes the detector or decays to lighter dark particles but h decays promptly to a collimated e + e − pair and produces the visible light that comprises the signal.
As shown in Fig. 1 [19] for R γ > 10, and the angular distribution (right panel) of the light due to the electron-like final state, for R γ > 1. The shaded blue region in both panels is our fit, and other shaded regions are the backgrounds.
plays an important role in producing the correct angular distribution in MB, as we discuss later. In our model, H and h predominantly couple to the first generation of quarks (u and d) and have negligible or tiny couplings to other families. The effective coupling (F N ) of either scalar to a nucleon (N ) can be written as [75][76][77] Here M N is the nucleon mass and the values of (f p Tu , f p T d , f n Tu , f n T d ) = (0.020, 0.041, 0.0189, 0.0451). In our scenario, f q = y H,h q , (q = u, d).
We include both the incoherent and coherent contribution in the production of N 2 in MB. For LSND, however, we consider only incoherent scattering from neutrons. The total differential cross section, for the target in MB, i .e., CH 2 , is given by The entire carbon nucleus (C 12 ) contributes in coherent scattering, with, however, decreasing contributions as |q 2 | = |(k − k) 2 | increases. To implement this, we employ a form factor, exp(−2b|q 2 |) [78]. Here b is a numerical parameter, which in the case of C 12 takes the value 25 GeV −2 [78,79]. The number of events is given by with E h ∈ [E h , E h + ∆E h ] and Φ ν is the incoming muon neutrino flux. η contains all detector related information like efficiencies, POT etc. All calculations for LSND, MB and the value of the muon g − 2 are carried out using the benchmark values in Table I. Finally, for these values, the calculated lifetimes of N 2 and h in the rest frame are 10 −17 s and 1.8 × 10 −12 s, respectively. Our results are presented in the next section.

V. RESULTS AND DISCUSSION
In this section we present the results of our numerical calculations, using the cross section for the process and the model described in Section III.
A. Results and discussion for MB and LSND Fig. 2 (top panels) shows the MB data points, SM backgrounds and the prediction of our model (blue solid line) in each bin. Also shown (black dashed line) is the oscillation best fit. The latest data set for the neutrino mode, corresponding to 18.75 × 10 20 POT, as detailed in [20] has been used in our fit. The left panel shows the distribution of the measured visible energy, E vis , plotted against the events for neutrinos. In our model, E vis is the same as E h . The angular distributions for the emitted light are shown in the right panel. The fit corresponds to benchmark parameter values shown in Table I. We have used fluxes, efficiencies, POT exposures, and other relevant information from [80] and references therein to prepare these plots. We see that very good fits to the data are obtained for both the energy and the angular distributions. The data points show only statistical uncertainties. We have assumed a 15% systematic uncertainty for our calculations. These errors are represented by the blue bands in the figures.
As mentioned earlier, the LSND observations measure the visible energy from the Cerenkov and scintillation light of an assumed electron-like event, as well as the 2.2 MeV photon resulting from coincident neutron capture on hydrogen. In our model, this corresponds to the scattering diagrams in Fig. 1 where the target is a neutron in the Carbon nucleus. Unlike the case of MB above, where both coherent and incoherent processes contribute to the total cross section, the LSND cross section we have used includes only an incoherent contribution. All necessary information on fluxes, efficiencies, POT etc for LSND has been taken from [19] and references therein. Fig. 2 (bottom-left panel) shows our results in comparison to the LSND data for R γ > 10, where R γ is a parameter defined by the LSND Collaboration (see, for instance [19]) that represents a likelihood ratio that the observed photon signalling the presence of the neutron was correlated as opposed to being accidental. This plot shows the energy distribution and the excess events in the data, as well as those resulting from our model using the same benchmark parameters as were used to generate the MB results. We find a total of 28.7 events from our model, compared to the 32 events seen by LSND for this choice of R γ . Fig. 2 (bottom-right panel) shows the angular distribution of the light due to the electron-like final state, for R γ > 1 and visible energies in the range 3 36 MeV < E vis < 60 MeV. In both panels, the blue shaded region is the result of our model, shown along with backgrounds and data.
Several points are pertinent to understanding the results obtained. We discuss them below: • All LSND events in our scenario stem from the high energy part of their DIF flux, which is kinematically capable of producing the N 2 (m N2 130 MeV). This flux originates in π + 's created in proton collisions in the LSND target (the experiment used two different targets over the running period, i .e., water and a high-Z material). This leads to a beam of ν µ 's, which interacts in the detector via ν µ CH 2 → n N 2 X → n N 1 h X → N 1 γ e + e − X (see Fig. 1). In the final step the photon is the correlated γ with an energy of 2.2 MeV signifying the capture of the neutron by a nucleus. The decays of both h and N 2 are prompt, while N 1 is either longlived and escapes the detector or decays to lighter invisible states. • In our scenario, both H and h act as mediators and contribute to the total cross section. The contribution of h is much smaller (∼ 10%) than that of H, since sin δ 0.1. However, this plays an important role in producing the correct angular distribution in MB. In particular, h is responsible for a coherent contribution which helps sufficiently populate the first (i .e., most forward) bin in the top-right panel of Fig. 2. • As a consequence of the heavy particle production (N 2 ) necessary, our model would not give any signal in KARMEN, which has a narrow-band DIF flux that peaks at ∼ 30 MeV, hence making it compatible with their null result.
• The DIF flux, in the oscillation hypothesis, generates electron-like events in energy bins beyond 60 MeV. Indeed, LSND saw 10.5 ± 4.9 such events (without a correlated neutron) in the range 60 MeV < E vis < 200 MeV, attributable to an oscillation probability of (2.9 ± 1.4) × 10 −3 [81]. Our model predicts 34 such events, which is within their acceptable range of uncertainty.
• LSND saw about 6 events with a correlated neutron in the energy range 60 MeV < E vis < 200 MeV, and our calculations yield 5.6 such events, in agreement with their observations. • As mentioned earlier, only incoherent neutron scattering contributes to the event counts in LSND.
We have assumed 8 MeV as the minimum energy transferred to a neutron in order to knock it out and register an event. Additionally, the masses of N 2 and N 1 are important factors in obtaining both the correct number and the correct distributions in this detector. Lowering the mass of N 2 increases the total events significantly, since it provides access to lower energies in the DIF flux spectrum. Decreasing the mass of N 1 shifts the event peak towards higher visible energies, and leads to higher numbers of correlated neutron events with energies > 60 MeV, which would conflict with what LSND saw. On the other hand, in MB the effects of N 2 and N 1 masses do not play as significant a role as they do in LSND, although the MB energy distribution improves if the N 1 mass is decreased from our current benchmark value.
• Finally, we note that the criteria as to when an e + e − pair constitutes a signal that may be counted as an electron-like event in both detectors are different. MB is not able to distinguish an e + e − pair from a single electron [52,82] if the invariant mass of e + e − < 30 MeV, or if the angle between the pair is 5 • or less. In our scenario, the mass of h , and hence that of the pair, is 17 MeV.
In LSND, the visible energies are quite low compared to those in MB. Hence, the opening angle of the e + e − pair can be large for the lower end of the visible energy (∼ 20 − 30 MeV). However, LSND did not attempt to search for e + e − pairs or γγ pairs, and for this reason it is reasonable to assume that it would reconstruct most e + e − pairs as a single electron event. In particular, because timing was their most powerful particle identifying variable, the fit to a Cherenkov ring would select the most significant ring, even for large angles between the e + e − pair. Therefore, e + e − pairs with correlated neutrons would explain the LSND excess [83], especially since no known e + e − or γγ backgrounds were expected in LSND. A more accurate calculation than the simple one performed here would incorporate the effects of fitting only the most energetic ring out of two which have a large angle between them. One effect of this would be to slightly increase the events in the middle bin (36 − 44 MeV) at the expense of those at higher energies (including those with energy > 60 MeV). It is evident from Fig. 2 (bottom-left panel) that this would improve the fit shown.

B. Muon anomalous magnetic moment
The one-loop contribution of the light scalars h , H to the muon g − 2 is given by [84,85] (29) where y φ µ is the coupling strength of φ-µ + -µ − interaction, which is defined in Eq. (15).
First, we note that m h , m H are fixed to fit the LSND and MB measurements. Also, both h and H contribute to the muon g − 2, ∆a µ and their ratio ∆a h µ /∆a H µ is proportional to tan 2 δ. In general, y µ and the angle δ would correspond to free parameters and they can be fixed to fit the ∆a µ central value.
For suitably chosen y µ and δ (y µ = 1.6 × 10 −3 and sin δ = 0.1), our benchmark yields values which lie in the experimentally allowed 2σ region, with ∆a µ = 2.24×10 −9 . For these values, the H contribution to ∆a H µ is dominant, with the h contribution being 16.6% of this quantity.

VI. CONSTRAINTS ON THE MODEL
This section is devoted to a discussion of constraints that the proposed scenario must satisfy, given the cou-plings of the extended scalar sector to fermions. We note here that in general, the off-diagonal couplings of the additional scalars in our model to down/up-type quarks are free parameters and can be very tiny, which is a relevant point that helps us stay safe from several existing bounds, as brought out below. A second relevant point in the discussion below is that we assume that the predominant decay mode for the lightest state among the dark neutrinos N i , i .e. N 1 , is to lighter dark sector particles.
Constraints from CHARM II and MINERνA: As discussed in [51], these experiments [86,87] constrain models attempting to explain the MB LEE and LSND based on results of high-energy ν-e scattering. A dark photon (Z ) model such as that discussed in [40] is tightly constrained for its chosen benchmark values, as shown in [51]. We also see that it is possible to evade this constraint provided the value of |U µ4 | (the mixing between the muon and the up-scattered sterile neutrino in the proposed model of [40]) stays equal to or below 10 −4 . In order to check that our model is safe from these constraints, we calculate the cross section contribution from our process, (with H, h as mediators) for CHARM II and MINERνA and compare its value with that for the model in [40], with |U µ4 | reduced to the safe value of 10 −4 . We find that the coherent cross section for our interaction stays more than an order of magnitude below this safe value, comfortably evading this constraint. We note that this is generically true for other recent models with scalar mediators, as also pointed out in [46,47].
We also note that elastic NC scattering of electrons with H, h as mediators is not a concern, since the final state contains a N 2 which promptly decays to an h N 1 and subsequently a prompt e + e − . This does not observationally resemble SM ν-e scattering.
Constraints from T2K ND280: As discussed in [53], the T2K near-detector, ND280 [88], is in a position to provide bounds on new physics related to the MB LEE. Relevant to our work here, the specific decay h → e + e − could be observable in this detector. Pair production can occur in the Fine-Grained Detectors (FGD), in particular. We have calculated the number of events for our process and find 9 events in FGD1, using a momentum cutoff of 300 MeV and an overall efficiency of 30%. This is comfortably below their bounds. In principle, at T2K, such events may occur in the TPC also. In our model, however, this decay is prompt, hence for detection in the TPC, the argon gas in it must act as both target and detection medium. Since the target mass is only 16 kg, however, the number of events is unobservably small in our case. We also note that the threshold for detection in the TPC is around 200 MeV.
Contributions to NC ν-nucleon scattering at high energies: Since the H and h in our model couple to neutrinos and quarks, a possible constraint arises from NC Deep Inelastic Scattering (DIS) of neutrinos on nucleons, to which these scalars would contribute as mediators. At high energies, IceCube and DeepCore are a possible laboratory for new particles which are produced via such scattering [89,90]. In the process shown in Fig. 1, the decay time of N 2 (leading to an e + e − pair) is short enough, to escape detection at these detectors. In terms of distances travelled, this corresponds to ∼ 1 m in DeepCore, and ∼ a few hundred meters in IceCube, even at very high energies. These lengths are are much smaller than the detector resolution necessary to signal a double bang event for both these experiments. In addition, we have checked that the high energy NC cross section stays several orders of magnitude below the SM cross section. We also note that N 1 in our model is assumed to decay predominantly to invisible particles, again escaping detection in large detectors.
Kaon and B-meson decay constraints: Prior to discussing specific cases, as a general remark, we note that in any heavy meson decay that involves u, d quarks, one can radiate an h which would promptly decay to an e + e − pair via the diagonal couplings between it and the quarks. While off-diagonal flavour changing couplings in our model are arbitrarily small, the first generation diagonal quark couplings to the scalars in our model are fixed by the requirements of fitting the LSND and MB data, and are approximately O(10 −5 ). These are small enough to suppress such decays by a factor O(10 −10 ), rendering them safe from existing upper bounds.
1. The BR(K L → π 0 e + e − ) < 2.8 × 10 −10 at a 90% C.L. has been measured at KTeV [91]. Hence in principle, the width for K L → π 0 h would contribute to this, while K L → π 0 H will not contribute, being kinematically forbidden. However, KTeV applies an invariant mass cut of 140 MeV for the e + e − pair, making the bound inapplicable due to kinematics. We note also, as mentioned already, off-diagonal couplings of h to d, s quarks in our scenario are tiny. Also, the BR(h → γγ) is negligible. Therefore, the constraints from K decays, e.g., K L,S → π 0 γγ [92] is not applicable.
2. The E949 Collaboration [93] and NA62 Collaboration [94] have measured the process K + → π +ν ν, which could be mimicked by the K + → π + h decay in our scenario. Since h decays primarily to e + e − , this means that it must be long-lived and escape the detector for this bound to apply. From [68], we see that given h mass of 17 MeV in our model, as long as lifetime is less than approximately 10 −10 s, one is safe from the constraint from invisibles. (In our model, h has a lifetime 1.8×10 −12 s.) Moreover, as mentioned above, off-diagonal couplings of h to d, s quarks in our scenario are tiny.
3. A light scalar coupled to muons can be emitted in the decay K + → µ + νφ. This is constrained by the NA62 [95], as discussed in [64]. H will evade these constraints because of its large mass, while h lies outside the constrained range due to its short lifetime and small coupling to muons. In addition, we note that data collected by the NA48/2 experiment [96] can in principle provide constraints via observation of K + → µ + νe + e − , as noted in [64]. This analysis has recently been done [96], but with an invariant minimum mass cut of 140 MeV for the e + e − pair, which makes it inapplicable to h .
4. The CHARM experiment [97,98] has measured the displaced decay of neutral particles into γγ, e + e − and µ + µ − . The relevant decays are K L → π 0 h and K + → π + h . Thus in our model, h can in principle be constrained by this experiment, but as discussed in [68], for m h 17 MeV, the lifetime in our case is much shorter than 10 −10 , which is the upper value set by this bound. Additionally, it is possible that CHARM, being sensitive to heavy neutral leptons given its dimensions, could have sensitivity to visible decays of N 1 . As noted earlier, however, N 1 decays primarily to invisible states. 5. The K µ2 experiment [99] has measured the K + → π + φ process. For our benchmark point (Table I), BR(K + → π + h ) 4.2 × 10 −12 [100], which is very small compared with the upper limit ∼ 10 −8 .
6. Similarly, the decay B → K * e + e − has been measured at the LHCb [101], BR(B → K * e + e − ) = (4.2 ± 0.7) × 10 −7 . In our model, this would correspond to B → K * h /H, with the latter going to e + e − . Given the b → s transition involved, and that couplings of H, h to b, s quarks can be arbitrarily small, we evade this constraint. Also, due to 2m e < m H < 2m τ , the decays B → K ( * ) H → K ( * ) µ + µ − are subjected to strong constraints from B → K ( * ) µ + µ − at the LHCb [102]. However, in our model, we evade these constraints because of the smallness of H coupling to b, s quarks. Finally, in our model, it is worth mentioning that the branching ratios of B → K ( * ) γγ/νν are negligible.
Constraints from neutrino trident production: The neutrino trident process [103] provides a sensitive probe of BSM physics at neutrino detectors, and has been measured [104][105][106]. It is relevant to our model given the couplings of the H, h to muons, which are used for our explanation to the muon g − 2 anomaly. Using the SM cross section and simple scaling, we have checked that our model is safe from this constraint. Pion decay constraints: H, h couple to quarks, hence H, h can mediate π 0 decay to e + e − . In the SM, this decay is loop-suppressed and consequently small. However, given the small couplings of the two scalars to the u, d quarks (∼ 10 −5 ) and the electron (∼ 10 −4 ), we find that we are safely below this constraint.
Collider bounds: At hadron colliders, the process Z → 4 proceeds via qq → Z * →¯ , along with a γ * attached to one of the external legs (either the quarks or leptons) and with the creation of a lepton pair from the γ * . This process has been measured at the LHC [107]. The Z and γ in principle can be replaced by H, h .
The bound from LHC, however, applies an + − invariant mass cut of 4 GeV, and hence does not apply to our situation.
Since H couples to leptons and quarks, H can be radiated from these in external and internal legs in any process, and then H → h h is possible, leading to two collimated pairs of e + e − , which will look like two leptons. An LHC search was conducted [108] and no significant deviation or excess was found. As discussed in [68], given the fact that the H to lepton/quark coupling is small ( < ∼ 10 −4 ), and also that the H → h h decay width is small due to the smallness of its coupling, the contribution will stay within the 1% level.
If H, h couple to b and s quarks, then the decay B s → µ + µ − can be mediated by them. This decay has been measured by both the LHCb and CMS (see [109][110][111][112][113]). However, in our model the couplings to b, s quarks can be arbitrarily small, hence this constraint can be avoided.
Constraints on y e and m h from dark photon searches: A generic dark photon search looks for its decay to a pair of leptons. One may translate such bounds [114,115] to constraints on a light scalar with couplings to leptons. Specifically, translated constraints relevant to our scenario arise from KLOE [116] and BABAR [117]. Current values of y e and m h in our scenario are safe from these bounds, but they will be tested in the future by Belle-II [118].
Constraints on y e and m h from electron beam dump experiments: As discussed in [63,64], a light scalar with couplings to electrons could be detected in beam dump experiments if it decays to an e + e − pair or to photons. For the mass range relevant here, the experiments E137 [119], E141 [120], ORSAY [121] and NA64 [122] can potentially provide restrictive bounds. While our present values are outside the forbidden regions, they will be tested in the future by the HPS fixed target experiment [123] which will scatter electrons on tungsten.
Constraints on y e and m H from dark photon searches: KLOE [124] searched for e + e − → U γ, followed by U decays to π + π − /µ + µ − leading to the constraint on y e (∼ 2 × 10 −4 ) at m U = 750 MeV. In our scenario, replacing U by H we note that the production of π + π − /µ + µ − by it in KLOE will be very suppressed due to its tiny coupling to u, d quarks and its predominant semi-visible decay (H → e + e − + / E). Moreover, both visible and invisible final states searches by BABAR [117,125] put upper limits on y e at m H , which can be evaded by H due to its predominant semi-visible decay.
Constraints on y µ and m H from colliders: BABAR has provided constraints [64,118] on these parameters via their search for e + e − → µ + µ − φ, where φ is a generic light scalar. Our values, while currently in conformity with these bounds, will be tested in the future by Belle-II [118]. In addition, BABAR [126] constrains a dark leptophillic light scalar with couplings which are proportional to m f /v. This set of constraints does not apply in our case since our couplings for h and H do not have this proportionality.
Contribution from the new scalars to the electron g − 2 anomaly: The (positive) one loop contribution in our model allows us to explain the observed value of ∆a µ . A similar (positive) contribution is made to ∆a e by both h and H, which we have computed, summed and found to be ∆a e = 4 × 10 −14 . This is well within the present uncertainties in this quantity. We note that our model allows the possibility of having negative offdiagonal Yukawa couplings to ∆a e by both h and H. This affords flexibility in varying this contribution and keeping it within acceptable limits, as well as possibly explaining the current ∆a e discrepancy at the one-loop level, as discussed in [47].

VII. CONCLUSIONS
Evidence for anomalous signals at low-energy noncollider experiments in general, and short baseline neutrino experiments in particular, has been gradually increasing over time, and has firmed up significantly over the past decade or so. Specifically, with reference to the LSND and MB excesses, it has gradually become evident that one may choose several different approaches towards understanding their origin, and these choices can lead down divergent and non-overlapping paths 4 .
An important premise underlying our effort in this paper is that a common, non-oscillation, new physics explanation exists for both LSND and MB. Furthermore, our effort is guided by the belief that such an explanation could not only yield a long-sought extension to the SM, but also delineate the contours of the portal connecting the SM to the dark sector, as well as shed some light on other related but as yet unresolved questions.
Pursuant to this, in the scenario presented here, the extension to the SM that these experiments lead to comprises of the well-known 2HDM. Access to the dark sector is achieved via mass-mixing with a (dark) relatively light singlet scalar and via the presence of heavier dark neutrinos in allowed gauge invariant terms in the Lagrangian. Two of the three CP-even scalars in the model are relatively light (m h 17 MeV and m H 750 MeV) and participate in the interaction that generates the excesses in LSND and MB, as well as contribute to the value of the muon g−2. Similarly, two of the three dark neutrinos not only participate in important ways in the interaction 5 in LSND and MB, but also, along with the third neutrino, generate neutrino masses via a simple seesaw mechanism.
The sub-GeV scalars in our model can be searched for in a variety of experiments. The masses of h and H lie especially close to existing bounds from electron beamdump experiments like E141 [120] and BABAR [117] respectively. Thus H can be searched for in Belle-II [118] and h in HPS [123]. The dark fermions in our model are amenable to searches in several upcoming experiments, e.g. DUNE ND [127] (for a more detailed discussion and references, see [128,129]).
In the near future, the MicroBooNE experiment [130][131][132] will provide first indications of whether the lowenergy electron-like event excesses in MB and LSND are due to electrons or photons. In the scenario presented here, the h has a very short lifetime prior to decay to an e + e − pair. At the energies under consideration, it would travel about 5 − 12 mm in the detector. Since tracks with a gap of greater than 1 cm would be interpreted as photons, most events resulting from our scenario would look like an excess of electrons in MicroBooNE, while the high energy ones could be mistaken for photons with short gaps. A dE/dx analysis would be required to actually detect that the events are e + e − pairs rather than electrons, which should also be possible with more data.
With respect to the scalar search for the h , we mention two existing experimental hints which are interesting: a) a significant excess in the 10−20 MeV invariant mass-bin of electron-like FGD1-TPC pairs detected by the T2K ND280 detector, (see Fig. 11 in [88]), and b) the higher than expected central value for the width Γ(π 0 → e + e − ) observed by the KTeV experiment [133], signifying the possible existence of a scalar with mass 17 MeV. Additionally, we would like to point out that the Kaon DAR search, planned at the JSNS 2 experiment [134,135] is in a position to provide a test of the proposal presented in this work via its flux of high energy ν µ [136].
In conclusion, we are hopeful that the long-standing and statistically significant anomalous results of LSND and MB, along with the connection established between them via the simple model presented here will help motivate a more focused search for these particles in ongoing and future experiments.