Constraints on Decaying Sterile Neutrinos from Solar Antineutrinos

Solar neutrino experiments are highly sensitive to sources of $\nu\to\overline{\nu}$ conversions in the $^8$B neutrino flux. In this work we adapt these searches to non-minimal sterile neutrino models recently proposed to explain the LSND, MiniBooNE, and reactor anomalies. The production of such sterile neutrinos in the Sun, followed the decay chain $\nu_4 \to \nu \phi \to \nu \nu \overline{\nu}$ with a new scalar $\phi$ results in upper limits for the neutrino mixing $|U_{e4}|^2$ at the per mille level. We conclude that a simultaneous explanations of all anomalies is in tension with KamLAND, Super-Kamiokande, and Borexino constraints on the flux of solar antineutrinos. We then present other minimal models that violate parity or lepton number, and discuss the applicability of our constraints in each case. Future improvements can be expected from existing Borexino data as well as from future searches at Super-Kamiokande with added Gd.


I. INTRODUCTION
Beyond neutrino mixing, the study of solar neutrinos provides important input to the standard solar Model (SSM) [1,2] and has been used to search for several phenomena beyond the SM of particle physics. One example is neutrino decay, originally proposed as an alternative solution to the solar neutrino problem [3,4]. Indeed, after precision measurements of the solar neutrino oscillation parameters by KamLAND [5], strong constraints on the lifetimes of ν 2 and ν 3 have been obtained [6,7] as data is consistent with no additional neutrino disappearance.
Recently, non-minimal neutrino decay models have received interest in the literature, where exotic decays of a relatively heavy (m exotic m active ) and mostly-sterile neutrino are invoked to explain longstanding experimental anomalies at short baselines (SBL). One category of models concerns "visible" sterile neutrino decay, where new sterile states are produced and decay back to visible active neutrinos. Originally proposed in Ref. [8] as an explanation to the Liquid Scintillator Neutrino Detector (LSND) anomaly [9,10], this scenario has now been revisited [11,12] in light of recent data of short-baseline ν µ → ν e and ν µ → ν e appearance at MiniBooNE [13][14][15], as well as ν e disappearance at reactors [16,17]. Due to the small mixing angles required in this explanation, the effects of attenuation in solar neutrino fluxes is small. Yet, the total number of heavy neutrinos produced is large, and if these states undergo sufficiently distinctive decays or scattering inside a detector, they can be searched for. In this article, we point out that if antineutrinos are produced in the decay of these heavy neutrinos, then they are strongly constrained by existing searches for neutrino-antineutrino transitions in solar neutrino experiments. * mhostert@umn.edu The flux of antineutrinos from the Sun at the MeV energies is negligible [18], which remains an excellent approximation down to tens of keV in energy [19]. Combined with the fact that the detection cross section for ν e is much larger and easier to measure compared to that of ν e , this makes solar neutrino experiments sensitive to very small fluxes of antineutrinos [20][21][22]. The current sensitivity reaches fluxes as small as a few times 10 −5 of the 8 B neutrino flux [23][24][25][26][27][28]. These searches have been discussed in the context of new physics, such as large ν → ν oscillations. This Lepton number (LN) violating process is rather small in most theories, being suppressed by (m ν /E) 2 , but can be enhanced due to spin-flavor precession [29,30]. The latter arises from the coupling of a large neutrino magnetic moment to the solar magnetic field, which induces ν e → ν x conversions, followed by flavor transitions into ν e due to matter effects. Another possibility to generate such LN violating signatures is neutrino decay. For instance, neutrino mass models where LN is a spontaneously broken global symmetry predict the existence of a pseudo-goldstone boson J, the majoron [31,32]. In these models, solar antineutrinos may be produced from the decay ν 2 → ν 1 J, which is enhanced in dense matter [33]. This possibility of production from neutrino decay is, in fact, quite general and can be realized in any LN violating model with neutrinos that decay sufficiently fast, be they ν 2 , ν 3 , or the new mostly-sterile state ν 4 discussed here (see Ref. [34] for an early discussion in the context of a 17 keV sterile neutrino).
In this work, we explore a new possibility where lepton number can, in fact, be conserved but the decay of a new light boson leads to a large flux of antineutrinos. We derive limits on the electron flavor mixing with ν 4 , working only with the gauge-invariant and parity-conserving model of Ref. [12]. Focusing solely on Dirac neutrinos, we show that our bounds exclude virtually all of the parameter space preferred that can simultaneously explain LSND and MiniBooNE, as well as the region of interest for reactor anomalies. They also disfavor most but not all parameter space suggested as a solution to the MiniBooNE anomaly. For models with Majorana neutrinos, the constraints become even stronger due to ν 4 → νϕ → ννν decays.
The paper is organized as follows. In Section II we review the benchmark model for decaying sterile neutrinos, and in Section III we discuss generic aspects of solar antineutrino searches. The resulting constraints, future prospects, and alternative search methods are then discussed in Section IV. We dedicate Section V to a survey of minimal alternative models for decaying steriles, and conclude in Section VI.

II. DECAYING STERILE NEUTRINO
The most significant deviations from the threeneutrino paradigm at SBLs are the LSND excess of ν e events, with a statistical significance of 3.8 σ when interpreted under a ν µ oscillation hypothesis, and the Mini-BooNE excess of ν e -like events, with a significance of 4.8 σ when interpreted under a ν µ → ν e and ν µ → ν e oscillation hypothesis. Reactors at very short-baselines also have some evidence of ν e disappearance [16,17], but in that case the neutrino flux predictions are highly uncertain and harder to control [35,36]. Despite the large significance of these anomalies, they remain unsolved. Their standard interpretation under oscillations of a eV-scale sterile neutrino leads to strong tensions between different data. This is driven mainly by the absence of anomalous results in ν µ disappearance experiments [37,38], as appearance and disappearance channels are strongly correlated in the oscillation scenario [39,40]. In addition, such new sterile states with eV masses are in strong tension with cosmological observations, which has prompted several studies to resolve this by invoking secret interactions [41][42][43][44][45].
Visible sterile neutrino decays are, therefore, a natural "next-to-minimal" explanation to SBL anomalies to consider. The advantages of this scenario are that it does not necessarily lead to strong correlations between appearance and disappearance channels, the mass scale of the new sterile state is not fixed by the oscillation length of the experiments, and that it already contains a secret interaction mechanism, possibly alleviating tension with cosmology. In this work, we focus on the decay of steriles with eV to hundreds of keV masses to a new scalar ϕ, as discussed in Refs. [8,11,12,46,47]. In all such visible decay scenarios, heavy neutrinos decay to mostly-active neutrinos via ν 4 → νϕ , where more neutrinos can be produced from ϕ → νν decay if ϕ is massive as in Ref. [12].
Here visible refers to the detectability of the decay products, in contrast to models where neutrinos decay to the wrong-helicity states that do not feel the weak interactions (up to tiny helicity-flipping terms proportional to m 2 ν /E 2 ν ). Such visible decays can explain the anomalous ν e -like events at SBL experiments by means of a sub-dominant population of ν 4 states in neutrino beams, which often decays to ν e -like daughters 1 . One typical prediction is that the spectrum of daughter ν e and ν e neutrinos is softer than the initial flux of ν 4 parents and associated neutrinos, skewing the effective flavor conversion towards lower energies. While this brings a mild improvement over the oscillation fit to the low energy excess observed at MiniBooNE, it leads to less satisfactory energy spectra at LSND, which is compatible with a signal that grows in energy. In addition, the neutrino flux at LSND comes from both π + and µ + decay at rest, yielding a large and monochromatic ν µ flux, and a spectrum of ν µ and ν e . Since only the ν e component is detected via the IBD process, the presence of a neutrino-to-antineutrino transition in the decay chain can convert the large ν µ flux to signal, since ν 4 states can be produced in pion as well as muon decays. This is a crucial point in the study of Ref. [12], which found improved compatibility between LSND and MiniBooNE regions of preference when this conversion is significant.
fashion by means of mixing between the heaviest neutrino state, ν 4 , and the active flavors. The relevant Lagrangian is given by where the neutrino mass mechanism is left unspecified and assumed to not play a role in the low-energy phenomenology. In the mass basis, the neutrino mass eigenstates are given by ν i = α U * αi ν α , with α ∈ {e, µ, τ, s}, and U a unitary mixing matrix. Under the assumption of parity conservation in the sterile sector, U is identical to the extended Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, now 4×4. We return to this issue in Section V. In the decays of ν 4 and ϕ, only the three lightest mass states are produced, and so the it is useful to define the low-energy flavor stateν s = 3 i=1 U αi ν i . For most processes of interest, however, the non-unitarity corrections introduced by working withν s instead of the full flavor states ν s is small and appears only at order |U α4 | 4 . Unless stated otherwise, we refer toν s as simply ν from now on, as the mass eigenstates have decohered on their way from Sun. The new scalar does not couple directly to the SM, and loop-induced couplings will ultimately depend on the UV completion of the model and its neutrino mass mechanism (see, for instance, Refs. [31,49]).
Due to mixing, heavy neutrinos with masses below the MeV scale would be produced in the Sun via the same processes responsible for ν e production at a rate |U e4 | 2 times smaller. Once produced, the ν 4 mass eigenstates immediately decay to a light neutrinos and the scalar boson. The scalar then decays to a neutrino-antineutrino pair, giving rise to our signal. Overall, the process of interest is For most cases of interest, E ν4 m 4 , so if ν 4 (ν 4 ) is produced via weak interactions, it will be left-handed (right-handed) polarized to a very good approximation. We then assume all heavy neutrinos to be polarized with a definite helicity h 4 = −1 for neutrinos and h 4 = +1 for antineutrinos. Nevertheless, due to the assumption of parity conservation for the ϕ interactions with neutrinos, both helicity flipping (HF) and helicity conserving (HC) decay channels are allowed. Assuming all neutrinos to be ultra-relativistic, we find the squared-amplitudes for polarized ν h4=−1 where r ϕ = m ϕ /m 4 . Integrated over phase space, both channels contribute identically to a total decay rate of Our decay rate is in agreement with Refs [12,50]. Note that helicity conserving decays prefer larger E 1 values, while helicity flipping decays prefer smaller values of E 1 . Therefore, for our present application, helicity-flipping decays are important since the antineutrinos from the subsequent scalar decay tend to be more energetic. Also important is the limit r → 1, where the scalar particle has most of the ν 4 energy regardless of the helicity structure of the decay. This is the scenario with the most energetic antineutrinos in the final state, for which a simultaneous explanation of MiniBooNE and LSND is most successful. The scalar decay length in the lab frame to leading order in the small mixing elements is As expected, the scalar decays are doubly suppressed by small mixing elements, and so it tends to decay more slowly than ν 4 . Nevertheless, the decay of both particles can be considered prompt within astrophysical objects. Finally, note that only due to parity conserving nature of the scalar interaction, both left-and righthanded antineutrinos are produced. In this case, only the right-handed antineutrinos (ν + ) are relevant for detection through weak interactions.

III. SOLAR ANTINEUTRINOS
The flux of MeV antineutrinos from the Sun in the SSM is negligibly small. The largest antineutrino flux at MeV energies comes from small fractions of long-lived radioactive isotopes in the Sun, namely 232 Th, 238 U, and mainly 40 K. This give rise to an antineutrino flux on Earth of about 200 cm −2 s −1 with E ν 3 MeV [18]. This component, however, is still 6 orders of magnitude smaller than the geoneutrino flux at the surface of the Earth at these energies, and can be safely neglected. At larger energies, photo-fission reactions produce an even smaller flux of antineutrinos of about 10 −3 cm −2 s −1 [18]. It is only down at the much lower energies of tens of keV that antineutrinos start being produced in thermal reactions at a similar rate to neutrinos with fluxes as large as 10 9 cm −2 s −1 [19].
Existing limits on the flux of solar antineutrinos are usually quoted in terms of an energy-independent probability P νe→νe of conversion of 8 B neutrinos into antineutrinos. The most stringent limits were obtained by Kam-LAND in 2011 [24] which was recently improved in 2021 [28], and by Borexino in 2019 [25] P Borexino all at 90% C.L. In addition, SuperKamiokande (SK) has derived limits on extraterrestrial ν e sources during phases I, II and III [51], but the high energy thresholds of E ν > 17.3 MeV make them irrelevant for the study of 8 B neutrinos. For SK phase IV (SK-IV), improvements to the trigger system were implemented and the detection of neutron capture on Hydrogen was made possible, lowering thresholds to E ν > 13.3 MeV [52]. The constraint on solar antineutrino flux was found to be P SK−IV−2013 νe→νe < 4.6 × 10 −4 . Recently, further improvements to the neutron tagging algorithm lowered this value to E ν > 8.3 MeV [27], and using the data 2008 to 2018 the limit was improved to A previous preliminary result was shown in Ref. [26] and an even more recent update was presented in Ref. [53]. Loading of Gd in the SK water tank is expected to greatly improve the neutron tagging efficiency, and would allow for much more stringent limits. With projections on the signal selection efficiency and background reduction, Ref. [27] finds that a limit of P SK−IV−Gd νe→νe (E ν ≥ 8.3 MeV) 2.2 × 10 −5 could be achieved with 0.2% Gd loading [27]. Lowering the energy threshold of the trigger could further improve these projections.
In addition to these, SNO has also set limits at the level of P SNO νe→νe (E ν ∈ [4, 14.8] MeV) < 8.3 × 10 −3 [23] at 90% C.L. All limits quoted above assume a total 8 B flux of 5.88 × 10 6 cm −2 s −1 , except Borexino which assumes 5.46 × 10 6 cm −2 s −1 , and KamLAND which assumes 5.94 × 10 6 cm −2 s −1 . At the lowest energies, a bound can also be obtained by noting that the number of elastic ν − e scattering events in solar neutrino experiments decreases if too many ν 4 states are produced, both due to lower ν e − e cross sections and suppressed ν e flux. These effects, however, are insensitive to variations of the total ν e flux below the tens percent level. The predictions from the sterile neutrino decay model are compared with the 8 B flux in Fig. 1. The independent bounds quoted by KamLAND, Borexino, and SK-IV are shown in Fig. 2 as a function of E ν .
The strength of the limits above is mostly due to the large cross section for Inverse beta decay (IBD) on free protons at MeV electron-antineutrino energies. Beyond dominating over the neutrino-electron elastic scattering cross section by about two orders of magnitude (see Fig. 1), this channel has a distinct signature that drastically reduces backgrounds. After produced, the positron annihilates and the final state neutron is quickly captured by the free protons. This results in a double-bang signal with a positron kinetic energy T e E ν − 1.8 MeV, and a delayed emission of a ≈ 2.2 MeV gamma. The cross section for this process is well understood at high [54] and low [55] energies, and relatively simple formulae that are valid for all energy regimes have been derived by Ref. [56]. In this work we implement the latter calculation, which is provided as machine-friendly data files by Ref. [57].
a. Backgrounds For Borexino, reactor neutrinos represent the largest source of backgrounds, but are effectively constrained by DayaBay measurements. Atmospheric neutrino events with genuine IBD scattering or otherwise inherit large uncertainties from the atmospheric neutrino flux and cross sections, but represent only a small contribution (6.5 ± 3.2 events). Finally, the 238 U and 232 Th geoneutrino fluxes are energetic enough to contribute to the IBD sample, but are only significant up to 3.2 MeV. Borexino omits the contribution of geoneutrinos from the Earth's mantle in their estimation, which is conservative. This component is the most likely explanation for the ≈ 2σ excess seen in the lowest energy bin [25] (see also their latest geoneutrino analysis [58]).
The reactor neutrino flux at KamLAND is dominant below 8.3 MeV, but contributes only about 2.2 events above that value. Due to the smaller overburden at Kam-LAND and SK, they suffer from larger spallation backgrounds, coming mainly from radioactive decays of 9 Li. After muon tagging and fiducial volume cuts, these are reduced to less than 5 events at both locations. The large number of neutrino-electron scattering events presents a background for SK. For this reason, a cut is applied requiring small shower angles with respect to the direction of the Sun, cos θ < 0.9. This does not impact IBD events as the positron angle with respect to the incom- ing neutrino is significantly larger ( cos θ ≈ 0) than in the predominantly forward process of elastic scattering. The observed event spectra and background predictions by the respective collaborations are shown in Fig. 3.

A. IBD Rates from Decays
The largest observable flux from sterile neutrinos would come from ν 4 states produced via weak interactions in the decay of 8 B. The number of IBD events in a given experiment can be computed as where N stands for total exposure of the experiment, Pν s→νe is the flavour transition probability for Solar antineutrinos averaged over the radius of the Sun (see Appendix A), and Note that Eq. (11) is the analogue of Eq. (9) from Ref. [12], and is simpler since we work with very long baselines and under the assumption that the number of initial ν µ states is negligible.  Fig. 4, but for mϕ/m4 = 0.5. In this regime LSND and MiniBooNE are harder to combine due to softer νe spectrum obtained via electron mixing. We show only the combination of MiniBooNE with other datasets, excluding LSND, as reported in Ref. [12]. The solar antineutrino spectrum is also softer, but can effectively constrain the preferred combined region at 99% C.L.
trino flux is 5.46 × 10 6 cm −2 s −1 . For low-metallicity models, our constraints on the new physics coupling are weakened by about 20%.
With the predicted number of IBD events at each solar neutrino experiment, we implement our statistical test (described in detail in Appendix B) to place limits on the active-heavy mixing angles. Our χ 2 test statistic models solar neutrino flux and experimental backgrounds uncertainties through bin-uncorrelated nuisance parameters with Gaussian errors. Both flux and background uncertainties are fixed at 10%, except for the SK-IV, for which we inflate those to 20%. For the KamLAND 2021 dataset, we combine the data into 2-MeV-wide bins when performing our fit. We have also performed a total rate fit for each of the energy bins to obtain model-independent limits on the solar antineutrino flux. We find similar to those provided by the collaborations, within 50%. Our limits are always weaker, and therefore more conservative, when compared to the official ones. For SK-IV, no model-independent limit was shown in the final article [27], so we show our own result in Fig. 2.

IV. RESULTS
We plot our 90% C.L. limits in Fig. 4 as a function of |U e4 | 2 and |U µ4 | 2 for m ϕ /m 4 = 0.9. On the same axes, we show the preferred (grey-shaded) regions obtained in Ref. [12] to explain LSND and MiniBooNE individually, as well as the combined fit to MiniBooNE, LSND and global data (except reactors and cosmology) as "All w/o cosmo". Weaker constraints from the OPERA [60] and KARMEN [61] neutrino experiments, as well as beta decay kink searches are shown as dashed grey lines. We pick two particular cases with the shortest ν 4 and ϕ life- times to compare against our limits, corresponding to m 4 Γ 4 = 1 eV 2 and m 4 Γ 4 = 10 eV 2 . These lifetimes are achieved for couplings close to the perturbativity limit, namely g 2 ϕ = (1.5) 2 × 10 −2 /(|U e4 | 2 + |U µ4 | 2 ) and g 2 ϕ = (12) 2 × 10 −2 /(|U e4 | 2 + |U µ4 | 2 ), respectively. It is clear that an explanation of LSND is in large tension with solar antineutrino searches for all three experiments we consider. It should also be noted that the region with large |U µ4 | 2 which is not excluded by our curves is excluded by MiniBooNE itself. As expected, KamLAND leads to the strongest bounds despite its large neutrino energy thresholds 8.3 MeV. Borexino and SK-IV bounds are competitive, with the latter performing better for harder antineutrino spectra.
In Fig. 5, we show our constraints for the case m ϕ /m 4 = 0.5. A simultaneous explanation of Mini-BooNE and LSND is more challenging now, and only a global fit including MiniBooNE but not LSND is available ("All w/o LSND"). In this case, we constrain the region preferred by MiniBooNE significantly. Lower values of m ϕ /m 4 are even more challenging from the point of view of explaining the SBL results, as ϕ becomes longer lived and helicity-flipping decays of ν 4 , which lead to soft daughter spectra, become more important.
Changing the heavy neutrino mass but keeping the ratio R = m ϕ /m 4 fixed leaves our bounds unaltered for all the mass range of interest. Lowering R, on the other hand, weakens our bounds slightly due to the softer solar antineutrino spectrum, although the weakening saturates once below m ϕ /m 4 ≈ 0.5. We show constraints on the electron mixing angle in Fig. 6. The strongest constraints are obtained for vanishing muon and tau mixing, and for R → 1, as in that case ϕ carries most of the energy of its parent particle ν 4 . These limits are independent of the absolute scale of m ϕ and m 4 , provided these masses are below the Q-values of the 8 B decays and above the light neutrino masses.
Accounting for perturbativity bounds on g ϕ and the baselines of LSND and MiniBooNE, a lower bound on m 4 of ≈ 100 eV can be obtained for neutrino mixing of the order of |U e4 | 2 ≈ 10 −3 . Below this value, ν 4 is too long-lived and does not lead to interesting signatures at short-baseline experiments. The preferred regions for MiniBooNE and LSND shift to larger mixing angles when either ν 4 or ϕ are longer-lived, while our constraints remain unaffected. The parameters used in Figs. 4 and 5 are chosen so as to minimize the decay length of new particles. Bounds from kink searches in beta decay become the strongest above m 4 5 keV, and peak searches in meson decay preclude an explanation of LSND with the current model at the MeV scale.
We note that a massless scalar has also been discussed as an explanation of the MiniBooNE and LSND anomalies [11], although it is disfavored with respect to the best-fit point of Ref. [12] at 99% C.L. Interpreting that study in the present lepton-number and parity conserving model, one would obtain no constraint from solar antineutrino searches. In that case, however, light neutrinos will also decay. We comment more on this possibility and others in V.

A. Future opportunities
Borexino has collected an additional O(500) days of data on top of the 2771 days already analyzed in Ref. [25].
In addition, the improvements made by the collaboration in the latest geoneutrino analysis [58], such as enlarged fiducial volume and improved background rejection, can be implemented in the solar antineutrino search. Just in terms of exposure, this represents an improvement of 40% with respect to the values we use.
a. Synergy with DSNB Solar antineutrino searches will become even more stringent with upcoming efforts to detect ν e events from the DSNB [62,63]. The SK detector is expected to detect this neutrino flux with the addition of Gd to its detector volume [64]. The large neutron capture cross sections on Gd and the emission of 8 MeV gammas will help reduce backgrounds and lower the ν e detection threshold to neutrino energies as low as the IBD threshold, provided E e > 0.8 MeV [65]. The increased detection efficiencies at lower energies, and reduced accidental and mis-reconstructed backgrounds will improve on the limits we set, being limited by intrinsic reactor and atmospheric ν e backgrounds, but also by an exponentially rising spallation background. Future large liquid-scintillator detectors, such as the Jiangmen Underground Neutrino Observatory (JUNO) [66] and the Jinping neutrino experiment [67], can also improve on current constraints with an expected threshold of E ν 8.5 MeV [68]. In the far future, observatories capable of accumulating larger numbers of DSNB events, such as the proposed detector THEIA [69], would play an important role in searching for solar antineutrinos. We note that in the event of a detection of the DSNB, one could also constrain the models considered here by requiring small DSNB absorption by relic neutrinos [70,71].
b. Light neutrino decay Even in the parity conserving model discussed so far, one can avoid solar antineutrinos by resorting to a massless ϕ. In that case, however, the light mostly-active neutrinos will decay. For typical parameters relevant for the SBL anomalies, this decay will happen within 1 AU, both visibly and invisibly. For instance, consider ν 2 → ν 1 ϕ and ν 3 → ν 1 ϕ decays with normal ordering and m 1 ≈ 0. For a coupling of g ϕ = 1, we find In the convention adopted by the literature, these correspond to τ 0 2 /m 2 = 4.4 × 10 −5 s/eV and τ 0 3 /m 3 = 1.3 × 10 −6 s/eV. For inverted ordering, both ν 1 → ν 3 ϕ and ν 2 → ν 3 ϕ decays are of the faster kind in Eq. (14). On top of the cosmological issues with such short lifetimes (for recent discussions, see Refs. [72,73]), the largest g ϕ values relevant for the allowed regions are already excluded by laboratory experiments, such as SNO [74], which is consistent with no disappearance of solar neutrinos (see Refs. [7,75]). Other datasets have also been discussed to constrain the lifetime of light neutrinos, including measurements of the flavor ratios of cosmic neutrinos [76] and of the Glashow resonance [77]. It should be noted, however, that light-sterile mixing parameters governing light neutrino decay are related to those of SBL anomalies in a model-dependent fashion. In principle, but not without fine-tuning, the correlation between |U e4 |, |U µ4 |, |U * s4 U sj |, and |U * si U sj | for i, j < 4 may be relaxed. In models with LN violation or LN charged scalars (see below), provided several constraints are satisfied, solar antineutrinos may become relevant again for massless ϕ as the light-neutrino decays ν 2 → ν 1 ϕ are open.

V. ALTERNATIVE MODELS: VIOLATING PARITY AND LEPTON NUMBER
Various other possibilities for visible sterile neutrino decay exist, depending on the Dirac or Majorana nature of neutrinos, as well as on the parity structure of the sterile neutrino sector. While we only focused on the parity conserving model discussed above, we would like to dedicate this section to understanding if other minimal extensions of the SM by a singlet sterile neutrino and a scalar are subject to our constraints. For clarity, we focus on SM extensions with a single new sterile neutrino: ν s = ν L s + ν R s in the Dirac case and ν R in the Majorana case. Our findings for the minimal models are summarized in Table I. a. Dirac L(ϕ) = 0 We start with a generalization of Eq. (1) by writing where index summation is understood. For complex g ϕ , this is the most generic parametrization of the vertex. We implicitly diagonalized the Dirac mass matrix by means of two unitary matrices V L and V R , defined by are the (chiral) mass eigenstates. Note that the enlarged PMNS matrix is defined as U PMNS = V L when the charged lepton mass matrix is diagonal. Abandoning the assumption of parity conservation in the sterile sector that was made previously, V L = V R , one can have allow for V L = V R by choosing different Yukawa couplings for ν L s and ν R s . By breaking parity at the level of the Dirac mass matrix, it is possible to independently tune the couplings appearing in the operators ν i P R ν 4 and ν i P L ν 4 . In practice, this allows to tune the rate for visible and invisible decays of ν − 4 neutrinos. The same is true for the decay of the scalars, which can be either visible or invisible, depending on V L and V R . In these models, a connection to the SBL anomalies through visible decays always predicts visible solar antineutrinos provided ϕ is heavy enough to decay.
b. Dirac |L(ϕ)| = 2 One can also introduce scalars carrying LN. These type of scalars have been usually discussed in the context of Majoron models, but for our current purposes, we assume no particular connection to neutrino masses. We consider a model with a Dirac field ν s , and a complex scalar ϕ carrying lepton number

Minimal Models
Parametric Limit Polarized ν − 4 decay Scalar decays Spectrum of ν + Expected Signals Dirac L(ϕ) = 0  < 4), and ν − 4 for the heavier left-handed polarized neutrino. The first two columns show the minimal model considered and the judicious choices of its parameters to achieve a certain parity structure. The third and fourth columns show the decay channels allowed in that model, separated by all possible helicity final states. The penultimate column shows the kind of visible solar antineutrinos (ν + ) energy spectrum is predicted in the model. This depends on the HF or HC nature of the ν − 4 decay. Note, however, that the spectrum may change in the limit mϕ/m4 → 1. For the minimal Majorana neutrino model, one always obtains a prediction for visible solar antineutrinos. See the main text for definitions.
In all generality, we can write where again we implicitly diagonalized the Dirac mass matrix with V L and V R . In this case, even for parity conserving mass matrices, one can violate parity in the neutrino-ϕ interactions by tuning the arbitrary g L and g R couplings. In this case, ν − 4 states produced in the Sun will always decay to visible antineutrinos provided g R = 0, independently of the decay products of ϕ * . For this model, one may attempt to explain SBL anomalies with only the decay products of scalar produced in ν − 4 decays by setting g R → 0. In that case, no visible solar antineutrinos appear.
c. Majorana neutrinos A final possibility is to abandon LN and work with Majorana neutrinos. In this case, a minimal model can be built with only ν R and a scalar ϕ. LN is violated by the ν R Majorana mass term, and the most general interaction Lagrangian in this case is where now we implicitly diagonalized the Majorana mass matrix by means of a single unitary matrix V = U PMNS . In this case, all light neutrinos as well as antineutrinos are visible due to the reduced number of degrees of freedom. Both HC and HF decays of ν 4 are controlled by the same parameters, and cannot be disentangled as easily. Solar antineutrinos could appear in this case if all other constraints are satisfied. d. Simplified models Finally, we note that Refs. [8,11] work with simplified models, and do not specify the origin of the L ⊃ g e ν e ν 4 ϕ vertex. Although this operator may arise from a Lagrangian as simple as Eq. (1), it can be considered more generically as a by-product of nonrenormalizable operators such as (LH) 2 ϕ and (LH)ν s ϕ. In these effective models, the active-heavy mixing necessary for ν 4 production in most accelerator experiments, |U µ4 | 2 , is independent from g e , which controls the decay rate of ν 4 → ν e ϕ. In this case, the only way to generate ν e appearance at LSND is via muon decays, µ + → e + ν e ν 4 . It also follows that the mixing |U e4 | 2 may be parametrically small, turning off ν 4 production in the Sun via mixing. Four-body decays of the type 8 B → 8 Be e + ν 4 ϕ are negligible as kaon decays constrain g 2 e |U µ4 | 2 < O(10 −7 ). If a vector particle is introduced instead, the cosmological history is yet even more involved. We do not study this case here, although our solar antineutrino bounds would also apply to parity-conserving scenarios with small modifications. Note that our constraints are not relevant for fully invisible sterile neutrino decays, as invoked to relax the tension between SBL appearance and disappearance tension in Refs [38,78].

VI. CONCLUSIONS
Puzzling results from some of the short-baseline neutrino experiments will eventually find an explanation with more data coming from the SBN program at Fermilab [79][80][81] and the π + decay-at-rest experiment at J-PARC, JSNS 2 [82]. At the moment it is possible to speculate that some form of new physics in the neutrino sector is responsible for the deviation of LSND and Mini-BooNE results from theoretical expectations within the minimal three-generation neutrino model. Among such speculations is the class of models where the excess of antineutrinos at LSND and excess of low-energy electronlike events at MiniBooNE is due to a prompt production and decay of dark sector particles. This new sector is likely to comprise a heavier, mostly sterile neutrino, that can be produced via neutrino mixing in meson decays and nuclear reactions. Such heavier neutrino can generate a cascade decay to an unstable bosonic mediator and light neutrino, giving rise to the admixture of electron antineutrinos in the flux at the end of the decay chain.
We have shown that up to some model dependence one should expect that regular nuclear processes in the Sun create an antineutrino flux. Such flux is stringently constrained by most of the solar neutrino experiments, at a O(few × 10 −5 ) level owing to a larger cross sections for the IBD processes, and additional structure to the signal that has been exploited to cut on backgrounds. After application of these constraints, our results disfavor large part of the parameter space of the model in Ref. [12], and disfavor this mechanism as an explanation of the LSND excess, while significantly narrowing possible parameter space for the MiniBooNE excess. The simulations used to produce the results in this paper are publicly available on github 2 .
In general, our limits disfavor large ν → ν transitions that could improve the combined fit of LSND and MiniBooNE data. Such transitions could in principle be avoided if ϕ is lighter than the lightest neutrino state, in which case, mixing angles and CP phases have to be fine-tuned to avoid light neutrino decays. The alternative models with parity violation or apparent LN violation presented in Section V may avoid ν → ν transitions even for massive ϕ, but would require a case-by-case study of the SBL physics and additional constraints.
Our constraints add to the existing list of problems of the decaying sterile neutrino solutions to the SBL puzzle. Chiefly among them is cosmology and astrophysics. As is well known, new and relatively strongly interacting states can be populated by the thermal processes leading to the modifications of observed quantities, such as primordial nucleosynthesis yields and/or total amount of energy density carried by neutrinos at late times. In addition, these models are likely to cause strong modifications to the supernovae neutrino spectrum. One reason for such modification is the possibility of the neutrino numberchanging processes, such as νν → ννϕ → νννν. Given relatively strong couplings in the models of Refs. [11,12], the underlying cross sections are far greater than weak interaction cross section, meaning that the neutrinos can share energy and maintain their chemical equilibrium immediately after they leave the star. The main physical effect, the degrading of average energy for the SN neutrinos, can be constrained with the observed signal of SN1987A. This has been explored to constrain neutrino self-interactions inside supernovae by the requirement that neutrinos carry sufficient energy to the outer layers of the collapsing star [83]. We point out, however, that a more general statement can be made regarding neutrino energy loss outside the dense environment, which is independent of the explosion mechanism. Details of this effect will be addressed in a future publication.
with A CC = ±2 √ 2E ν G F N e (r) for neutrinos (antineutrinos) proportional to the electron density at the production region. Since the sterile component does not feel neutral-current interactions, the neutral-current potential also modifies the relation above. Terms proportional to A NC = ∓ √ 2E ν G F N n (r), where N n (r) is the neutron number density in the Sun, are small, however. This is because they are proportional to the new small mixing angles, as well as due to the smaller number of neutrons in the Sun, N n /N e 1/2. In our simulation, we follow Ref. [84] and keep all corrections in θ 14 , θ 24 , θ 34 , and θ 13 in the transition probabilities. As an approximation, we assume the neutral-current potential to follow the same radial dependence as the charged-current one, with A NC /A CC = −1/4.
The full transition probabilites using Eq. (A1) in the energy region of the 8 B flux are shown in Fig. 7 for both neutrinos and antineutrinos. In the limit |U µ4 |, |U τ 4 | → 0, the scalar decays produce mostly-ν e states, and Pν s →νe P νe→νe . In this case, the flavor evolution is similar to the standard MSW effect, where neutrino (antineutrinos) undergo resonant (nonresonant) adiabatic flavor conversion, in which production of ν m 2 at the center of the Sun is enhanced (suppressed). For |U µ4 |, |U τ 4 | = 0, the situation is more complex, but the high energy behavior can be understood by taking the limit θ 12 → π/2 (0) in the mixing elements |U si | 2 for neutrinos (antineutrinos). Note that ν m 3 ∼ ν 3 , as it should be since ∆m 2 21 ∆m 2 31 .

Appendix B: Statistical Method
When deriving upper limits on the mixing angles, we minimize the following log-likelihood function where θ stands for the vector of physics parameters (e.g., |U α | 2 ), β the vector of nuisance parameters with individual entries β j and associated Gaussian errors σ j . As an approximation, we assume L to follow a χ 2 distribution when estimating our confidence intervals.
The most important systematics for our study are the uncertainties on the total 8 B solar neutrino flux and total backgrounds numbers. To be conservative, we assign each energy bin two normalisation systematics, one exclusive to the new physics prediction, modelling uncertainties in the solar flux, and one exclusive to backgrounds. All normalization systematics are assumed to be uncorrelated, which is conservative, and are assigned 10% Gaussian errors.