Since the dark photon interacts with the SM electric charge in an analogous way to the photon, the upscattering cross section can be calculated in an analogous way to the transition magnetic moment case, replacing the leptonic vertex and the propagator by a massive one. When the dark photon is light, coherent scattering with the nucleus is dominant and the rate is enhanced by the number of protons in the nucleus, $Z$, with respect to scattering on individual protons. For $m_{Z^\prime} \gtrsim 1$ GeV, the incoherent contributions start to dominate and the upscattering can kick a proton out of the Carbon nucleus in MiniBooNE, for example. In models where the HNL only mixes with muon neutrinos, the upscattering cross section can be shown to scale as $\sigma_{\nu_\mu A \to N A} \propto |U_{\mu 4} U_{D4}|^2$ by virtue of the unitarity of $U$. Therefore, $\sigma_{\nu_\mu A \to N A}$ can be readily constrained by existing limits on the active-heavy mixing $|U_{\mu 4}|^2$ for every choice of $|U_{D4}|^2$, which usually takes values close to unity.

The decay of the HNL is prompt when $m_{Z^\prime} < m_4$, where the dark photon is produced on shell and decays to $e^+e^-$. Decays to light neutrinos are suppressed by the mixing angle combination $|U_{Di} U_{Dj}|^2$, which is small for $i<4$. This regime produces very forward-going HNLs, which subsequently decay to boosted $e^+e^-$ pairs. While this produces more $e^+e^-$ events that mimic single photon or single electron showers, it also leads to very forward-peaked angular distributions, in contrast with MiniBooNE's observation <cit.>. The angular spectrum is less forward for $m_N \gtrsim 300$ MeV, however, at those masses constraints from high-energy experiments are severe <cit.>.

Left: the spectrum of dark neutrino events in neutrino and antineutrino energy spectrum, as well as in the angular distribution at MiniBooNE. Right: fit to the neutrino-energy distribution at MiniBooNE in a dark neutrino model with a light dark photon ($m_{Z^\prime} < m_4$). Figure from <cit.>.

For heavy dark photons, $m_{Z^\prime} > m_4$, the decay is a three-body one and, therefore, the HNL is much longer-lived. Ref. <cit.> proposed a model where $|U_{\tau 4}|^2 \gg |U_{\mu 4}|^2$, so that the decay process is effectively $N\to \nu_\tau e^+e^-$. A similar proposal was made in Ref. <cit.> where a model with two HNLs was used. In that case, the heaviest HNL decays to the intermediary state with a lifetime that is not suppressed by active-heavy mixing. Due to the heavy mediator, the upscattering happens with a larger $Q^2$, and the angular distribution can be less forward. The best agreement with the angular distribution is found when the dark photon interferes with the SM $Z$.

Left: Fit to the neutrino-energy distribution (left) and angular distribution (right) at MiniBooNE in a dark neutrino model with a heavy dark photon ($m_{Z^\prime} < m_4$). Figure from <cit.>.

Dark scalars

HNLs can also interact with additional scalars that play the role of the dark photon discussed above. In this case, the upscattering cross section lacks the t-channel singularity and does not asymptote to a constant, but rather falls as $1/E_{\nu}^2$ at large energies. Therefore, if the dark sector particles have masses of $\mathcal{O}(100)$ MeV, the upscattering process is largest at the energies of LSND and MiniBooNE and may avoid constraints from high-energy experiments altogether. Models of this type have been discussed in Ref. <cit.>.

As pointed out in Ref. <cit.>, the scalar mediator can also lead to the production of photons pairs via the decay chain $N\to \nu (S \to \gamma \gamma)$. At MiniBooNE, if $S$ is lighter than the pion, it will be more boosted and its decays to overlapping photons would mimic the excess signal. It is also possible that the branching ratios for $S \to e^+e^-$ and $S\to \gamma \gamma$ are both sizeable, in which case the excess would display a non-trivial shape in energy and angle.

Explaining LSND

Some of the models proposed in the literature above can, in principle, also explain the LSND anomaly. This sets a more stringent requirement on the theory since it requires the upscattering process to produce a neutron in the detector. The new signal should mimic the inverse-beta-decay signature with at LSND, with a prompt electromagnetic signal followed by delayed neutron capture. This already eliminates the dark photon as a potential solution, since the couplings with neutrons are much more suppressed than those with protons. In addition, depending on the mass of the heavy neutrino, the upscattering process may have too large of a threshold to be initiated by $\pi$ or $\mu$ decays at rest. Instead, for masses of $\mathcal{O}(100)$ MeV, one may take advantage of the number of pions that decay in flight. This requires a large cross section, since the at these energies is much smaller than the decay-at-rest one. Models of this type were proposed in Refs. <cit.>, where the heavy neutrinos are produced via the exchange of a new scalar boson with large couplings to neutrons.

Broad requirements to explain MiniBooNE and LSND

Such proposals must conform to strong  demands from  a) cross section requirements in order to yield a sufficient number of total events in both LSND and MB, b) the measured energy and angular distributions in both experiments, and finally, c) compatibility of the new physics model and its particle content with  bounds from an extensive swathe of particle physics experiments <cit.>.

We compare how scalar and vector mediators perform in helping to achieve a simultaneous  understanding of both anomalies. Our treatment is necessarily brief, and for details  and  a full set of references on all that follows the reader is referred to <cit.>.

We start by breaking up the interaction into two parts. We first consider the tree-level process leading to the up-scattering of an incoming muon neutrino, $\numu$, to a heavy neutral lepton ($N_2$) in the neutrino detector as shown in <ref>(a), with the underlying assumption that it subsequently decays promptly in the detector into either of the final states shown in  <ref>(b) and  <ref>(c). The mediator for the up-scattering process could be either a light neutral vector boson $Z'$ or a light CP-even scalar $H$. The relevant interaction Lagrangian for the up-scattering process in each case  is given by

\begin{eqnarray}
\mathcal{L}^{Z'}_{\rm{int}} \!\!&=&\!\! (C^{Z'}_{\nu} \bar{\nu}_{iL} \gamma^{\mu} N_{Lj}Z'_{\mu} +h.c.) +  C^{Z'}_{n} \bar{U}_{n} \gamma^{\mu}U_n Z'_{\mu},\\
\mathcal{L}^H_{\rm{int}} \!\!&=&\!\! (C^H_{\nu} \bar{\nu}_{iL} N_{Rj}H +h.c.) + C^H_{n} \bar{U}_{n} U_n H,
\end{eqnarray}

where $U_{n}$ is the nucleon field and $i, j=1,2,3$. For simplicity, we assume that for both mediators ($Z'$ and $H$), the coupling to a  proton is  the same as its  coupling to a neutron. The up-scattering cross section depends on the overall product of the coupling constants $C^{Z'}_{\nu} C^{Z'}_n$ ($C^{H}_{\nu} C^H_n$) for the vector (scalar) mediator.

(a)     (b)     (c)

Feynman diagrams for the production and the subsequent decay of the $N_2$ in LSND and MB.

$m_{Z'\!/\!H}$     $C^{Z'}_{\nu} C^{Z'}_n$     $C^{H}_{\nu} C^{H}_n$

50 MeV     1.04$\times 10^{-8}$     9.3$\times 10^{-8}$

1 GeV      8.5$\times 10^{-7}$     2.14$\times 10^{-6}$

300 MeV     $-$     5.35$\times 10^{-7}$

The overall coupling values for the vector and scalar mediators for different values of mediator masses ($m_{Z'\!/\!H}$) to produce 560 $N_2$ in the MB final state. The mass of $N_2$ is 100 MeV.

Based on the Lagrangian above and the benchmark parameter values in <ref>, we compute the coherent and incoherent cross sections (shown in <ref>) for both mediators,  required for generating the central value of the number of excess events in MB, which is $560$ <cit.>. From <ref>, we note that
the scalar and vector-mediated cross sections behave distinctly, and our representative calculations bring out the following qualitative points:

* For all cases the cross section initially rises as the energy is increased from its lowest values, however, it subsequently drops for a scalar mediator whereas it remains approximately flat with increasing energy for a vector mediator. This is true for both the coherent and incoherent parts. As the neutrino energy rises, it is this relatively rapid drop in the cross section  for $H$  that allows solutions with a scalar mediator to comfortably
skirt constraints <cit.> coming from CHARM II <cit.> and  MINER$\nu$A <cit.>, compared to the $Z'$ mediated process.

* It can also be seen  from <ref> that the coherent contribution dominates over the incoherent part for lighter mediator masses, whereas the reverse is true for the higher mass choice for mediators.

* For LSND, contributions to events come from the incoherent part of the cross section only, due to  the presence of a neutron in the final state. In the region in the left panel of  <ref>, we note that the  energy drops from MB ($\sim 800$ MeV) to LSND DAR flux values ($\sim 150-200$ MeV). Note that  for $m_{Z'\!/\!H}=1$ GeV (solid curves), while the incoherent cross section drops  for both mediators,  the vector cross section has lower values to begin with compared to the scalar. It also drops  more rapidly.  For example, it can be seen that the cross section for the $Z'$ drops an order of magnitude over this energy range for $m_{Z'}= 1$ GeV.

The incoherent (coherent) cross section per CH$_2$ molecule (C atom) as a function of incoming neutrino energy. The overall constants for different kinds of mediator masses are taken from <ref>.

For $m_{Z'\!/\!H}=50$ MeV (dashed curves), over this energy range, the incoherent scalar cross section is significantly higher than  that for the vector. In fact,  it increases as the energy is lowered, unlike its vector counterpart. This reduction  in the incoherent vector cross section at   energies ($< 800$ MeV) makes it more difficult for models with a vector to give a sufficient number of electron-like excess events at LSND, even though by using a high enough $m_{Z'}$ one may  successfully evade the CHARM II and  MINER$\nu$A bounds.  However, on the other hand, too low $m_{H}$ gives many more events than those observed in LSND,
both in the $20-60$ MeV visible energy range which recorded data, and beyond $60$ MeV, where only a limited number of events were seen.

* Finally, for  scalar mediators, especially those with low masses $m_H\simeq 100$ MeV, the cross section tends to rise  at low values of $E_\nu$. However, in such models, if  the primary decay modes of $N_2$ are to invisible daughters, as in <cit.>, the incoherent  interaction would mimic the neutral current interaction $\nu N \rightarrow\nu N$. This has  been measured at MB <cit.> at these energies and found to be in agreement with the SM, providing an important restriction on such models.

Overall, the cross section and mediator mass considerations for a common solution thus appear to favour scalar mediators over vectors.  Secondly, our representative calculations also point to a preference for  lighter (but not ultra-light)  mediators if both excesses are to have a simultaneous solution.

An examination of the angular distribution of MB is also useful from the point of view of  imposing requirements on proposed solutions. The  excess in MB  is distributed over all directions but is  moderately forward.   The  cross section responsible for the production of $N_2$  as a function of the cosine of the angle between the momentum direction of $N_2$ and the beam direction has been studied in a  bin-wise manner in <cit.> for both the coherent and incoherent contributions.  It was found that when $m_{Z'\!/\!H}=50$ MeV, almost all the produced $N_2$ are in the  most forward bin for both mediators. However, as the mediator mass is increased to $1$ GeV,  there was a shift in the distribution,  and the other bins were also populated for both types of mediators, even though there were qualitative differences between the two.
From this, at first it appears that  using a single scalar mediator and adjusting its mass to an intermediate value, as well as the mass of $N_2$ will allow us to find a common solution to the two anomalies as well as match the angular distribution in MB. However, further examination based on considerations related to the energy distributions in LSND and MB (for details, see <cit.>) reveals that this is not the case if good fits to both anomalies are desired.

Overall, as detailed in <cit.>, energy distributions in LSND and MB, the angular distribution in MB, when combined with the stringent constraints on light singlet scalars, all suggest the use of a scalar doublet, with one light and one moderately heavy partner.  We find this leads to a degree of angular isotropy while allowing a large number of events in the forward direction, consistent with observations. A combination of a moderately heavy and a light mediator complement each other  well when a common solution to the two anomalies is sought.   An example solution  to both  anomalies that incorporates  all the features that have been obtained in our study  has been provided in <cit.>.

Our insistence on a solution that addresses both anomalies simultaneously is, of course, a choice. It restricts  proposed solutions in ways that attempt to explain the anomalies individually do not.  However, it is noteworthy that  once we demand this, and adhere to the  dictates of the cross section,  the observed energy and angular distributions in both experiments  as well as the many constraints from various experiments <cit.>, then we are  led to a simple extension of the SM that i) resolves both anomalies, ii) provides a portal to the dark sector, iii) accounts for the experimentally observed value of the muon $g-2$ and iv) addresses the issue of neutrino mass via a Type I seesaw, in conformity with the  global data on the observed values of neutrino mass-squared differences in oscillation experiments, as shown in <cit.>.

§.§.§ Long-Lived Heavy Neutrinos

Heavy neutrino decays to single photons In Ref. <cit.>, the authors propose a solution to the MiniBooNE excess with a heavy neutrino that interacts with muon-neutrinos through mixing as well as through a transition-magnetic moment. The model is the same as the one presented in <ref>, but contrary to the solution in <ref>, the mixing angle $|U_{\mu N}|$ and transition magnetic moment $d_{\mu}$ are small, such that the heavy neutrino decays only in macroscopical distances. In particular, if the lifetime is larger than the distance between the target and the MiniBooNE detector, $c\tau^0 \gtrsim \mathcal{O}(500)$ m, then the heavy neutrino can be produced in pion and kaon decays via mixing, $K \to \mu N$, travel to the detector, and decay inside the active volume due to the transition magnetic moment, $N \to \nu \gamma$.

While the model was shown to successfully reproduce the energy spectrum of the MiniBooNE excess, the decay process tends to produce very forward signatures. A better agreement with the angular spectrum is achieved for larger masses of the $N$. However, if the heavy neutrino is too massive, its arrival in the MiniBooNE detector is delayed with respect to the SM neutrinos. This delay is proportional to the $(m_N/E_N)^2$, and should not exceed values much larger than several nanoseconds, as the excess events have been observed to be in the same time window as the beam neutrinos <cit.>. This sets a strong constraint on the mass of the heavy neutrino. For masses below $150$ MeV, where agreement with the timing requirements can be achieved, production in pion decays should also be considered.

Heavy neutrino decays to axion-like particles
In order to account for the MiniBooNE excess, Ref. <cit.> proposes an extension of the SM with Dirac HNL that couples to a leptophilic axion-like particle ($\ell$ALP).
The HNL, denoted by $N_{D}$ to emphasize its Dirac nature, mixes with the three SM neutrino flavors, $\nu_{\beta}$. The flavor and mass eigenstates, $\nu_{jL}$ and $N_{D}$, can be transformed into each other by means of a unitary matrix $U$ <cit.>:

\begin{align}
\nu_{\beta} = \sum^{3}_{j=1} U_{\beta j} \nu_{jL} + U_{\beta 4} N_{D}.
\end{align}

The relevant Lagrangian includes the interaction of $\ell$ALPs with sterile neutrinos and electrons <cit.>:

\begin{align}
\mathcal{L}_{a\ell} = -\frac{\partial_{\mu}a}{2f_{a}} \left(c_{N}\overline{\nu}_{D}\gamma^{\mu}\gamma^{5}\nu_{D} + c_{e}\overline{e}\gamma^{\mu}\gamma^{5}e\right),
\end{align}

where $f_{a}$ is the $\ell$ALP decay constant, $c_{N}$ and $c_{e}$ are the dimensionless parameters for the $\ell$ALP-sterile neutrino and $\ell$ALP-electron couplings, respectively.

$\ell$ALP scenario for the explanation of the MiniBooNE excess. $L$ is the travel distance of the sterile neutrino, $D$ is the diameter of the MiniNooBE detector, and $\theta_a$ is the scattering angle of the $\ell$ALP. Figure from <cit.>.

Our $\ell$ALP scenario is sketched in Fig. <ref>. In this framework, the sterile neutrino with a mass $m_N$ of $\mathcal{O} \left(100\right)$ MeV, the $\ell$ALP with a mass $m_{a}$ of $\mathcal{O} \left(10\right)$ MeV and an inverse decay constant $c_{e}/f_{a} \simeq \mathcal{O} \left(10^{-2}\right)~\mathrm{GeV}^{-1}$ are considered. The Dirac-type sterile neutrino $N_{D}$, produced from charged kaon decays through its mixing with the muon neutrino, travels a distance of 500m and decays into a $\ell$ALP and a muon neutrino inside the MiniBooNE detector.

The electron-positron pairs produced from the $\ell$ALP decays can be interpreted as electron-like events provided that their opening angle is sufficiently small. We verify that with appropriate choices of the parameters, the sterile neutrino and $\ell$ALP have the proper mean decay lengths consistent with the setup of the MiniBooNE experiment. We also make sure that the values of the parameters adopted in our model are allowed by the astrophysical and experimental constraints.

The comparison of the numerical results of the angular and visible energy spectra in the $\ell$ALP scenario and the data of the MiniBooNE experiment in the neutrino mode. Figure from <cit.>.

Following the approach applied in <cit.>, the angular and visible energy spectra in the neutrino mode are computed and compared with the results of the MiniBooNE experiment, as shown at the top of Fig. <ref>. At the bottom panel, the predicted total event numbers are obtained after summing over the spectra and shown as the contours consistent with the MiniBooNE excess at the $1\sigma$ to $3\sigma$ levels on the $\left(m_{N},\left\vert U_{\mu 4}\right\vert^{2}\right)$ and $\left(c_{e}/f_{a},\left\vert U_{\mu 4}\right\vert^{2}\right)$ planes, with the constraints obtained from other experiments.
We find that the scenario with the sterile neutrino mass in the range 150 MeV $\lesssim m_{N} \lesssim$ 380 MeV and the neutrino mixing parameter between $10^{-10} \lesssim \left\vert U_{\mu 4} \right\vert^{2} \lesssim 10^{-8}$ can explain the MiniBooNE excess.

Best-fit regions to the MiniBooNE excess in the HNL and $\ell$ALP parameter space. Limits and best-fit regions are shown on the $\left(m_{N},\left\vert U_{\mu 4}\right\vert^{2}\right)$ (left) and $\left(c_{e}/f_{a},\left\vert U_{\mu 4}\right\vert^{2}\right)$ (right) planes. Figure from <cit.>.

§.§.§ Dark Matter Particles

Since the neutrinos at MiniBooNE are produced primarily from charged meson decays and the decays of daughter muons of those charged mesons, neutrino-based solutions can accommodate the absence of any excess in the dump mode, in which the charged mesons are no longer focused by magnetic horns, unlike the neutrino and antineutrino modes. Essentially, the neutrino-based explanations work well because a key feature of the excess seems to be correlated to the focusing or suppression of charged mesons.

(a) \(\pi^+, K^+\);
[right=1.6cm of a] (b);
[above right=1.1cm of b] (f1) \(\nu_\mu\);
[below right=1.0cm of b] (c);
[above right=1.1cm of c] (f2) \(\phi, a\);
[below right=1.1cm of c] (f3) \(\mu^+\);
(a) – [scalar] (b) – [fermion] (f1),
(b) – [anti fermion] (c),
(c) – [scalar] (f2),
(c) – [anti fermion] (f3),

(a) \(\pi^+, K^+\);
[right=1.6cm of a] (b);
[above right=1.1cm of b] (f1) \(\nu_\mu\);
[below right=1.0cm of b] (c);
[above right=1.1cm of c] (f2) \(A^\prime\);
[below right=1.1cm of c] (f3) \(\mu^+\);
(a) – [scalar] (b) – [fermion] (f1),
(b) – [anti fermion] (c),
(c) – [boson] (f2),
(c) – [anti fermion] (f3),

3-body charged meson decay into a scalar, pseudoscalar, or vector. Analogous processes exist for $\pi^-$ and $K^-$ decay.

This poses a challenge to dark sector interpretations of the excess (e.g., using $\pi^0$ or dark bremsstrahlung production channels <cit.>), which have been more constrained and less holistic in their explanation of the excess with respect to their counterparts in neutrino BSM physics thus far. One solution proposed in Ref. <cit.> opens up the possibility for dark sector explanations by means of connecting the dark sector to the physics of charged meson decays, something that had previously been overlooked. The authors in Ref. <cit.> considered various DM scenarios involving couplings to muons, namely those in Eq. (<ref>),( <ref>), and (<ref>):

\begin{eqnarray}
\label{eq:s}
\mathcal{L}_{S} &\supset&  g_\mu \phi \Bar{\mu} \mu + g_n Z^\prime_\alpha \Bar{u} \gamma^\alpha u + \frac{\lambda}{4} \phi F^\prime_{\mu\nu} F^{\mu\nu} + \textrm{h.c.}, \\
\label{eq:p}
\mathcal{L}_{P} &\supset& i g_\mu a \Bar{\mu} \gamma^5 \mu + g_n Z^\prime_\alpha \Bar{u} \gamma^\alpha u + \frac{\lambda}{4} a F^\prime_{\mu\nu} \Tilde{F}^{\mu\nu} + \textrm{h.c.}, \\
\label{eq:v}
\mathcal{L}_{V} &\supset& e (\epsilon_1 V_{1,\mu}+\epsilon_2 V_{2,\mu}) J_{\rm EM}^\mu + (g_1 V_{1,\mu}+g_2 V_{2,\mu}) J_D^\mu+(g^\prime_1 V_{1,\mu}+g^\prime_2 V_{2,\mu}) J^{\prime\mu}_D.
\end{eqnarray}

Here, three massive bosons have been introduced; a long-lived scalar $\phi$ and pseudoscalar $a$, and a short-lived vector $A^\prime_\alpha$ decaying to DM fermions $\chi$, $\chi^\prime$ with $F^\prime_{\mu\nu} \equiv \partial_\mu A^\prime_\nu - \partial_\nu A^\prime_\mu$. The muonic couplings allow for the 3-body decays of the form $M \to \mu \nu X$ ($M=\pi^\pm, K^\pm$) became possible, as shown in Fig. <ref>. The 3-body nature of this decay mechanism is not phase-space suppressed in the contraction of $\nu^\dagger \mu_\uparrow$, as opposed to the ordinary 2-body decay which selects out only the combination $\nu_\downarrow \mu_\uparrow$. The branching ratios for scalar, pseudoscalar, and vector DM production are shown in Fig. <ref>. In the $\mathcal{L}_V$ scenario, the detector signature would then take place through DM upscattering, $\chi N \to \chi^\prime N (\chi^\prime \to \chi e^+ e^-$), while scenarios $\mathcal{L}_S$ and $\mathcal{L}_P$ consider long-lived $\phi$/$a$ scattering in the detector through a Primakoff-like process $\phi N \to \gamma N$ via a heavy mediator $Z^\prime$. The parameter space which fits the MiniBooNE excess for $\mathcal{L}_S$, $\mathcal{L}_P$, and $\mathcal{L}_V$ is shown in Fig. <ref>.

Branching ratios for the 3-body production of scalars, pseudoscalars, and vectors via a charged meson with mass $M = m_{\pi^+}, m_{K^+}$.

However, the muonic portal is not the only possibility, and DM production in the 3-body $\pi^\pm$ and $K^\pm$ decays from DM-quark couplings can also be treated. In this sense, the scope of the 3-body decay solution can cover several coupling schemes, and should be testable at other experiments with similar meson production capabilities.

The credible regions for fits to the MiniBooNE cosine spectrum with $\mathcal{L}_V$ (left) and $\mathcal{L}_{S,P}$ (right) are shown at 68% (dark-shaded) and 95% (light-shaded).

This has broad implications for accelerator facilities, such as LBNF (DUNE), that also have magnetic focusing horns that should be sensitive to forward-produced DM from the meson decays. Constraints on the parameter space from accelerator-based searches at CHARM, MINER$\nu$A, and T2K can also be considered, but their smaller POT and exposures do not give them sensitivity to the MiniBooNE excess. Other neutrino experiments such as CCM, JSNS, and COHERENT that produce stopped mesons and lack magnetic focusing horns can also probe the parameter space relevant for the excess, since while the DM signal from meson decays will be isotropic, their detectors are situated much closer to the beam targets to be sensitive to the DM flux.

§.§ Conventional Explanations

While the majority of the explanations explored above rely on some new or BSM physics to give rise to the various anomalies, the possibility that the origin lies in more conventional explanations, such as an underestimated background, mis-modelling in simulation or over-constrained cross-section uncertainties, must still be considered. While these explanations are generally difficult to test directly without access to (often) collaboration-internal experimental tools and data sets, there have been several attempts in recent years to test individual anomalies in this direction. In one such example <cit.>, it was shown that allowing a combination of theoretical uncertainties in different background channels to fluctuate in unison is not sufficient to resolve the MiniBooNE anomaly; however, it can reduce the significance of the MiniBooNE excess. In this section, we discuss several such possible “conventional” interpretations for the anomalies.

§.§.§ Single-Photon Production

The solution to the MiniBooNE puzzle may have important implications for our understanding of neutrinos and their interactions. In addition to the interpretation in terms of new physics, the MiniBooNE anomaly could be a manifestation of new forces of nature, while unaccounted or poorly modeled SM backgrounds cannot be entirely discarded. Once Cherenkov detectors like MiniBooNE misidentify single photon tracks as electrons, the excess of events could be due to their production through both SM and BSM mechanisms.

Feynman diagrams for NC single photon emission considered in the literature. The first two diagrams stand for direct and crossed baryon pole terms with nucleons and baryon resonances in the intermediate state. The third diagram represents $t$-channel meson exchange contributions.

In the SM, single photons can be emitted in NC interactions, NC$1\gamma$, on nucleons, $\nu (\bar{\nu})\, N \rightarrow \nu (\bar{\nu})\, \gamma \, N$, or on heavy nuclei, via incoherent or coherent (where the nucleus remains in its ground state) scattering. Theoretical models for the elementary NC$1\gamma$ <cit.> take into account $s$- and $u$-channel amplitudes with nucleons and $\Delta(1232)$, but also heavier baryon resonances, in the intermediate state, Fig. <ref>. The structure of nucleon pole terms at threshold is determined by symmetries. The extension towards higher energy and momentum transfers, required to predict cross sections at MiniBooNE, is performed by the introduction of phenomenologically parametrized  weak and electromagnetic form factors. The same strategy has been adopted for resonance terms. The $\Delta(1232)$ excitation followed by radiative decay is the dominant mechanism, as correctly assumed by MiniBooNE, but the contribution of non-resonant terms is also sizable. The uncertainty in the elementary NC$1\gamma$ cross section is dominated by the leading $N-\Delta$ axial transition coupling, $C_5^A(q^2=0)$, which is related to the $\Delta N\pi$ coupling (known from $\pi N$ scattering) by a Goldberger-Treiman relation, but has also been found to be $C^A_5(0) = 1.18 \pm 0.07$ <cit.> from data on $\pi$  production induced by neutrino scattering on hydrogen and deuterium. It is worth stressing that the $\Delta N \gamma$ couplings responsible for the resonance radiative decay are directly related to helicity amplitudes $A_{1/2,3/2}$ known with few-percent accuracy from photo-nucleon interactions <cit.>. Furthermore, owing to isospin symmetry, these quantities also constrain the vector part of the weak $\Delta$ production. This implies that large uncertainties in the $\Delta$ radiative decay couplings are at odds with hadron phenomenology but would also have an observable impact in weak pion production.  Among the non-resonant contributions, $t$-channel $\omega$-meson exchange, was proposed as a solution for the MiniBooNE anomaly <cit.> because of the rather large (although uncertain) couplings and the $\omega$ isoscalar nature, which enhances its impact on the coherent NC$1\gamma$ reaction. However, actual calculations found this contribution small compared to $\Delta(1232)$ excitation <cit.>. Nuclear effects, in particular the broadening of the $\Delta$ resonance in the nucleus, are important for single photon emission <cit.>. On $^{12}$C they cause a reduction of about $30$% in the cross section (see Fig. 9 of Ref. <cit.>).

SM prediction for NC1$\gamma$ events at  MicroBooNE  for $6.6 \times 10^{20}$ POT in $\nu$-mode.

With the ingredients outlined above, the SM single-$\gamma$ contribution to  the number events in the MiniBooNE detector and their distributions have been calculated <cit.> using the available information about the detector mass and composition (CH$_2$), neutrino flux and, quite significantly, photon detection efficiency. Results are to a large extent consistent with the (data based but relying on an improvable reaction model) MiniBooNE estimate (Figs. 4, 6 and 8 of Ref. <cit.>) and, therefore, insufficient to explain the excess. The impact of two-nucleon meson-exchange reaction mechanisms has been recently investigated <cit.> and found to be small (a factor of around 9 at $E_\nu=500$ MeV compared to single-nucleon mechanisms). In Ref. <cit.> it was estimated that the NC$1\gamma$ background should be enhanced by a factor between 1.52 and 1.62 over the MiniBooNE estimate, depending on the energy range and mode. Such an enhancement, shrinks the excess and significantly reduces the appearance-disappearance tension in global fits but is at odds with the earlier described theoretical calculations. An upper limit for the NC1$\gamma$ cross section on liquid argon has been  recently obtained by the MicroBooNE experiment <cit.>, as discussed in the following section. It disfavors that the excess could be solely attributed to this reaction channel but new results with higher statistics are required for a firm conclusion. Assuming $6.6 \times 10^{20}$ POT from the experiment's run plan, the distributions of the NC $1\gamma$ events calculated with the model of Ref. <cit.> are given in Fig. <ref> (adapted from Ref. <cit.>.). Comparison to future data shall offer valuable information about this process.

§.§.§ Reactor Flux Modeling

One conventional explanation for the Reactor Antineutrino Anomaly (RAA) has been gaining momentum thanks to the significant experimental and theoretical progress made in the past decade.
Reactor neutrino experiments that can measure oscillations without any reliance on reactor neutrino models (discussed further in Sec. <ref>) have been cutting away at RAA suggested sterile neutrino oscillation parameter space.
Additionally, neutrino, nuclear physics, and nuclear theory evidence has recently emerged suggesting that RAA may, at least in part, be caused by problems with reactor flux models.
While experimental developments are overviewed in Sec. <ref>, we will discuss here recent theoretical advancements in the flux prediction landscape that lend support for non-BSM origin to the RAA.

Model     $\bm{\sigma_{235}}$     $\bm{\sigma_{238}}$     $\bm{\sigma_{239}}$     $\bm{\sigma_{241}}$

$6.74 \pm 0.17$

$10.19 \pm 0.83$

$4.40 \pm 0.13$

$6.10 \pm 0.16$

$6.29 \pm 0.31$

$10.16 \pm 1.02$

$4.42 \pm 0.22$

$6.23 \pm 0.31$

$6.82 \pm 0.18$

$10.28 \pm 0.84$

$4.45 \pm 0.13$

$6.17 \pm 0.16$

$6.41 \pm 0.14$

$ 9.53 \pm 0.48$

$4.40 \pm 0.13$

$6.10 \pm 0.16$

Theoretical IBD yields of the four fissionable isotopes in units of $10^{-43} \text{cm}^{2}/\text{fission}$ predicted by different models <cit.>.

Model     2cRates     2cEvolution     2cRates + Evolution

$\bm{\overline{R}_{\text{mod}}}$     RAA     $\bm{\overline{R}_{\text{mod}}}$     RAA     $\bm{\overline{R}_{\text{mod}}}$     RAA

$0.936 {}^{ + 0.024 }_{ - 0.023 }$

$0.933 {}^{ + 0.025 }_{ - 0.024 }$

$0.930 {}^{ + 0.024 }_{ - 0.023 }$

$0.960 {}^{ + 0.033 }_{ - 0.031 }$

$0.975 {}^{ + 0.032 }_{ - 0.030 }$

$0.975 {}^{ + 0.032 }_{ - 0.030 }$

$0.925 {}^{ + 0.025 }_{ - 0.023 }$

$0.925 {}^{ + 0.026 }_{ - 0.024 }$

$0.922 {}^{ + 0.024 }_{ - 0.023 }$

$0.975 {}^{ + 0.022 }_{ - 0.021 }$

$0.973 {}^{ + 0.023 }_{ - 0.022 }$

$0.970 \pm 0.021$

Average ratio $\overline{R}_{\text{mod}}$ obtained in Ref. <cit.> from the least-squares analysis of the reactor rates in Tab. <ref> and of the Daya Bay <cit.> and RENO <cit.> evolution data for the IBD yields of the models in Tab. <ref>. The RAA columns give the corresponding statistical significance of the reactor antineutrino anomaly. The descriptions of all the models are given in text.

The statistical significance of the RAA depends not only on the magnitude of offset between reactor data and the Huber-Mueller prediction, but also on the size of the error bands applied to those predictions.
In the years following the inception of the RAA in 2011, a variety of reactor modeling studies have argued that the 2-3% error budget assigned to this prediction is likely under-estimated.
In particular, the role of forbidden beta transitions in altering the flux and spectrum reported by conversion predictions is not considered in the formulation of Huber-Mueller model error bands <cit.>.
When naive treatments of forbidden decay contributions are included in the prediction,  variations in flux of 4% or more are observed <cit.>.
This sentiment of under-estimated conversion prediction normalization uncertainties is fairly broadly shared across the community, as reflected, for example, in Ref. <cit.>.
On the other hand, a more recent conversion calculation (termed the HKSS model) that attempts to account for forbidden transition contribution using nuclear shell model-based calculations shows strong deviations from the Huber-Mueller model in spectral shape, but finds no major discrepancy in reported IBD yields <cit.>.

The past ten years have also brought about substantial development of state-of-the-art summation calculations, thanks to improved nuclear data evaluations and new nuclear structure measurements for a range of high-Q, high-yield fission  daughters <cit.>.
The improved summations can be compared to conversion calculations to provide an assessment of the latter's robustness.
As early as 2012, improved summation calculations were shown to generate reduced flux predictions with respect to earlier iterations <cit.>, indicating modest over-prediction of  fluxes in the Huber-Mueller prediction.
When comparing to conversion predictions, modern summation calculations were shown to predict different fuel-dependent flux variations <cit.>.
Most recently in 2019, Estienne, Fallot et al <cit.> used comprehensive improvements and updates in nuclear database information to generate a summation model (referred to as the EF Model) with a total predicted flux a few percent lower than the measured global IBD yield average, but with a  prediction $\sim$6% smaller than that predicted by the HM conversion model.
These comparisons are suggestive of a possible issue with conversion-predicted fluxes for individual isotopes.

Contours of the $2\sigma$ allowed regions in the ($\sin^2\!2\vartheta_{ee},\Delta{m}^2_{41}$) plane obtained from the combined neutrino oscillation fit of the reactor rates in Tab. <ref> and the Daya Bay <cit.> and RENO <cit.> evolution data. The blue, red, green, and magenta curves correspond, respectively, to the HM, EF, HKSS, and KI models in Tab. <ref>. Also shown are the contour of the $2\sigma$ allowed regions of the Gallium anomaly obtained in Ref. <cit.> from the combined analysis of the GALLEX, SAGE and BEST data (orange curve), and the $2\sigma$ bound obtained from the analysis of solar neutrino data in Ref. <cit.> (dark red vertical line).

The IBD yields predicted by these various models can be compared to the global IBD yields datasets listed in Tab. <ref> and to specialized Daya Bay <cit.> and RENO <cit.> dataset reporting changes in yield with variations in the content of the reactor cores; the latter datasets are described further in Sec. <ref>.
In particular, if comparisons are made under the hypothesis that the RAA is generated by sterile neutrino oscillations, the different reactor flux models generate substantially differing interpretations.
To illustrate, Fig. <ref> shows the contours of the $2\sigma$ allowed regions in the ($\sin^2\!2\vartheta_{ee},\Delta{m}^2_{41}$) plane of 3+1 active-sterile neutrino mixing parameters <cit.>.
One can see that an indication in favor of neutrino oscillations only for the HM and HKSS conversion models, which exhibit a significant reactor rate anomaly above $2\sigma$ (see Tab. <ref>).
Considering the EF model, for which the reactor rate anomaly is less than a few percent, the $2\sigma$ exclusion curves in Fig. <ref> allow only small values of $\sin^2 2\vartheta_{ee}$, including $\sin^2 2\vartheta_{ee}=0$, which corresponds to a lack of any statistically significant indication of sterile neutrino oscillations.
It should be noted here that error estimates for summation calculations are ill-defined, but are generally expected to be similar in magnitude to those provided by conversion predictions; efforts to provide more robust error envelopes are underway <cit.>.

For all the reactor flux models, it is also worth noting that upper bounds exist for the value of the mixing parameter $\sin^2 2\vartheta_{ee}$, with exact limits dependent on the value of $\Delta m_{41}^2$.
For $\Delta m_{41}^2 \gtrsim 2 \, \text{eV}^2$ the upper bounds for $\sin^2 2\vartheta_{ee}$ are between $0.14$ and $0.25$.
Figure <ref> shows that these bounds and the solar bound <cit.> are in agreement, but in tension with the large mixing <cit.> required to explain the anomaly of the GALLEX <cit.>, SAGE <cit.>, and BEST <cit.> Gallium experiments with short baseline neutrino oscillations.
This is a puzzling recent development in the phenomenology of short-baseline neutrino oscillations that may require an explanation that may extend beyond both conventional explanations and the simplest possible model of 3+1 active-sterile neutrino mixing.

§.§.§ The Gallium Anomaly and Interaction Cross-Section Uncertainties

The first analysis establishing the existence of the Gallium Anomaly <cit.> did not consider the uncertainties on the cross section of the detection process Eq. (<ref>) and, as mentioned in previous sections, this has been a source of subsequent investigation as a possible avenue for resolving this anomaly.

It is now clear that this quantity is of paramount relevance for the possible interpretation of the gallium anomaly. Different calculations of the cross section have been published after the seminal work by Bahcall <cit.>, and the values are summarized in Tab. <ref>.

$\sigma$ ($10^{-46}$ cm$^2$)     Ref.
$58.1  ^{+0.21}_{-0.16}$     <cit.>

$59.3  \pm 0.14$      <cit.>

$59.10 \pm 0.114$     <cit.>

$56.7  \pm 0.06$      <cit.>

$59.38 \pm 0.116$     <cit.>

Summary of the cross section values recently published.

Also, it has been pointed <cit.> that the rather large uncertainties come from the fact that only the cross section of the transition from the ground state of $^{71}\rm{Ga}$ to the ground state of $^{71}\rm{Ge}$ is known with precision from the measured rate of electron capture decay of $^{71}\rm{Ge}$ to $^{71}\rm{Ga}$. In fact, recent improvements on measurements of this transition <cit.> indicate that the gallium anomaly persists. However, electron neutrinos produced by processes in Eq. (<ref>), can also be absorbed through transitions from the ground state of $^{71}\rm{Ga}$ to two excited states of $^{71}\rm{Ge}$ at 175 keV and 500 keV.

Nuclear levels for the $^{71}\rm{Ga}$ transitions to $^{71}\rm{Ge}$.

When the aforementioned uncertainties are taken into account, different results for the total cross section for the radioactive sources are obtained, resulting in changes to the measured and expected $^{71}\rm{Ge}$ event rates. On the other hand, complete calculations of the cross sections of the interaction process, Eq. (<ref>), for neutrinos produced by $^{51}\rm{Cr}$ and $^{37}\rm{Ar}$
sources are given in Ref. <cit.>

\begin{equation}
\sigma = \sigma_{\rm{gs}} \left( 1 + \xi_{175}\frac{\rm{BGT}_{175}}{\rm{BGT}_{\rm{gs}}} + \xi_{500}\frac{\rm{BGT}_{500}}{\rm{BGT}_{\rm{gs}}} \right)
\,.
\label{eq:CrossSec_BGT}
\end{equation}

Here, $\sigma_{\rm{gs}}$ is the cross section of the transitions from the ground state of $^{71}\rm{Ga}$ to the ground state of $^{71}\rm{Ge}$, $\rm{BGT}_{\rm{gs}}$ is the corresponding Gamow-Teller strength, and $\rm{BGT}_{175}$ and $\rm{BGT}_{500}$ are the Gamow-Teller strengths of the transitions from the ground state of ${}^{71}\rm{Ga}$ to the two excited states of ${}^{71}\rm{Ge}$ at about 175 keV and 500 keV as shown in Fig. <ref>. The coefficients of $\rm{BGT}_{175}/\rm{BGT}_{\rm{gs}}$ and $\rm{BGT}_{500}/\rm{BGT}_{\rm{gs}}$ are determined by phase space:
$\xi_{175}({}^{51}\rm{Cr}) = 0.669$,
$\xi_{500}({}^{51}\rm{Cr}) = 0.220$,
$\xi_{175}({}^{37}\rm{Ar}) = 0.695$,
$\xi_{500}({}^{37}\rm{Ar}) = 0.263$

Table <ref> shows the different values of the rate when three different approaches (other than the one by Baxton, $R_{\rm{B}}$) are used to compute the cross sections <cit.>: $R_{\rm{HK}}$ uses information about the Gamow-Teller strengths from Haxton <cit.> and Krofcheck et al. <cit.>; for $R_{\rm{FF}}$, the corresponding numbers are taken from Frekers et al. <cit.>; $R_{\rm{HF}}$ uses $\rm{BGT}_{175}$ from Haxton and $\rm{BGT}_{500}$ from Frekers et al.; and $R_{\rm{JUN45}}$ uses calculations using nuclear shell-model wave functions obtained by exploiting recently developed two- nucleon interactions <cit.> (see also <cit.> for additional details).

GALLEX$_1$     GALLEX$_1$     SAGE$_{\rm{Cr}}$     SAGE$_{\rm{Ar}}$     Avg.
$R_{\rm{B}}$     $0.95 \pm 0.11$     $0.81 {}^{+0.10}_{-0.11}$     $0.95 \pm 0.12$     $0.79 \pm 0.08$     $0.86 \pm 0.05$

$R_{\rm{HK}}$	    $0.85 \pm 0.12$     $0.71 \pm 0.11$     $0.84{}^{+0.13}_{-0.12}$     $0.71 \pm 0.09$     $0.77 \pm 0.08$

$R_{\rm{FF}}$	    $0.93 \pm 0.11$     $0.79{}^{+0.10}_{-0.11}$     $0.93{}^{+0.11}_{-0.12}$     $0.77{}^{+0.09}_{-0.07}$     $0.84 \pm 0.05$

$R_{\rm{HF}}$	    $0.83{}^{+0.13}_{-0.11}	$     $0.71 \pm 0.11$     $0.83{}^{+0.13}_{-0.12}$     $0.69{}^{+0.10}_{-0.09}$     $0.75{}^{+0.09}_{-0.07}$

$R_{\rm{JUN45}}$	    $0.97 \pm 0.11$     $0.83 \pm 0.11$     $0.97 \pm 0.12$     $0.81 \pm 0.08$     $0.88 \pm 0.05$

Ratios of measured and expected $^{71}\rm{Ge}$ event rates in the GALLEX and SAGE experiments. The last column corresponds to the weighted average <cit.>. Information in the last row is from Ref. <cit.>.

The 90% allowed regions in the $|U_{e4}|^2$–$\Delta{m}^2_{41}$ plane
obtained from the analysis of the measured and expected $^{71}\rm{Ge}$ event rates <cit.>.

The last column of Tab. <ref> shows the corresponding weighted average for each case showing a gallium anomaly with a statistical significance of $3\sigma$, $2.9\sigma$, $3.1\sigma$, and $2.3\sigma$, respectively for the last four cases. This confirms the gallium anomaly and retains the indication in favor of a short-baseline disappearance of electron neutrinos, possibly due to neutrino oscillations. Figure <ref> shows the 90% contours in the $|U_{e4}|^2$–$\Delta{m}^2_{41}$ plane obtained from the analysis of the measured and expected $^{71}\rm{Ge}$ event rates, considering the neutrino survival probability Eq. (<ref>), where one can see that the squared-mass difference $\Delta{m}^2_{41}$ allowed values are $\sim1$ eV$^2$ or larger, for the different approaches to compute the $\nu_e$–${}^{71}\rm{Ga}$ cross section.

§.§ Summary of Interpretations

A summary of the interpretations detailed in this section is provided in Tab. <ref>.
The columns, from left to right, are the following: broad classes of models; specific models that fall in each class; the experimental signature of each model; which anomalies each model can address (LSND, MiniBooNE, Reactor and Gallium anomalies, respectively), and the corresponding references.
For convenience, in the first column, we also indicate the sections of this document that are related to each class of models.

Finally, in Tab. <ref> we summarize which experimental efforts can probe which signatures of new physics that can address the anomalies.
For concreteness, we focus on experiments that are either recent, are expected to be upgraded, or are still under proposal.
The leftmost column shows broad classes of experiments grouped arbitrarily by their source and type of experiment.
The other columns present which specific experiments, in each of these classes, can probe a given experimental signature resulting from different interpretations of the anomalies.
For clarity, the “Decays in flight” column title refers to particles produced by the decay of mesons, muons, or taus in flight.

2*Category     2*Model    2*Signature     4c|Anomalies     2*References
4-7             LSND     MiniBooNE     Reactors     Sources

[c]3*[-3em] [t]3 cm
Flavor transitions Secs. <ref>-<ref>, <ref>
(3+1) oscillations     oscillations                     Reviews and global fits <cit.>

(3+1) w/ invisible sterile decay     oscillations w/ $\nu_4$ invisible decay                       <cit.>

(3+1) w/ sterile decay     $\nu_4 \to \phi \nu_e$                       <cit.>

[c]2*[-2em] [t]3 cm
Matter effects  Secs. <ref>, <ref>
(3+1) w/ anomalous matter effects      $\nu_\mu \to \nu_e$ via matter effects                     <cit.>

(3+1) w/ quasi-sterile neutrinos      $\nu_\mu \to \nu_e$ w/ resonant $\nu_s$ matter effects                     <cit.>

[c]2*[-2em] [t]3 cm
Flavor violation  Sec. <ref>
Lepton-flavor-violating $\mu$ decays     $\mu^+\to e^+ \nu_{\alpha}\overline{\nu_e}$                     <cit.>

neutrino-flavor-changing bremsstrahlung     $\nu_\mu A \to e \phi A$                     <cit.>

[c]2*[-2em] [t]3 cm
Decays in flight Sec. <ref>
Transition magnetic mom., heavy $\nu$ decay     $N\to \nu \gamma$                      <cit.>

Dark sector heavy neutrino decay     $N\to \nu (X \to e^+e^-)$ or $N\to \nu (X\to \gamma \gamma)$                     <cit.>

[c]2*[-2em][t]3 cm
Neutrino Scattering Secs. <ref>, <ref>
neutrino-induced upscattering     $\nu A \to N A$,  $N\to \nu e^+e^-$ or $N \to \nu \gamma \gamma$                     <cit.>

neutrino dipole upscattering     $\nu A \to N A$,  $N\to \nu \gamma $                     <cit.>

[c]2*[-2em] [t]3 cm
Dark Matter Scattering Sec. <ref>
dark particle-induced upscattering      $\gamma$ or $e^+e^-$                     <cit.>

dark particle-induced inverse Primakoff     $\gamma$                     <cit.>

New physics explanations of the short-baseline anomalies categorized by their signature. Notation: – the model can naturally explain the anomaly,  – the model can partially explain the anomaly,   – the model cannot explain the anomaly.

3+1 Oscillations
Anomalous matter effects
Lepton flavor violation
Decays in flight
Neutrino-induced upscattering
Dark-particle-induced upscattering

DANSS upgrade, JUNO-TAO, NEOS II, Neutrino-4 upgrade, PROSPECT-II

Radioactive Source
BEST-2, IsoDAR, THEIA, Jinping

2|>p7.2 cm|IceCube upgrade, KM3NET, ORCA and ARCA, DUNE, Hyper-K, THEIA              2|>p4.8 cm|IceCube upgrade, KM3NET, ORCA and ARCA, DUNE, Hyper-K, THEIA

Pion/Kaon decay-at-rest
JSNS$^2$, COHERENT, CAPTAIN-Mills, IsoDAR, KPIPE         JSNS$^2$, COHERENT, CAPTAIN-Mills, IsoDAR, KPIPE, PIP2-BD             COHERENT, CAPTAIN-Mills, KPIPE, PIP2-BD

Beam Short Baseline
SBN             3|>p7.2 cm|SBN

Beam Long Baseline
3|>p9.6 cm|DUNE, Hyper-K, ESSnuSB      3|>p7.2 cm|DUNE, Hyper-K, ESSnuSB, FASER$\nu$, FLArE

Muon decay-in-flight
3|>p9.6 cm|$\nu$STORM         $\nu$STORM

Beta Decay and Electron Capture      KATRIN/TRISTAN, Project-8, HUNTER, BeEST, DUNE-$^{39}$Ar, PTOLEMY, $2\nu\beta\beta$

Summary of future experimental prospects to probe new physics explanations of the anomalies. We emphasize that all experiments can constrain the new physics models discussed in this paper in one way or another, but we list those that can provide a direct test of the respective model.

Summary of new physics scenarios that may explain all or a subset of the anomalies. Gallium anomaly column is missing. We did not refer to a cosmology anomaly, so we need to explain in the caption what that column means or otherwise remove it.

Add extension for Alves model (see commented out info in latex...)

§ BROADER EXPERIMENTAL LANDSCAPE

In this section, we review the broader landscape of existing experimental results with relevant sensitivity to interpretations of the LSND, MiniBooNE, Reactor, and Gallium anomalies.

§.§ Flavor Conversion

Other short-baseline, as well as long-baseline, neutrino experiments can look for anomalous flavor conversions. In the following sections, we discuss direct tests of MiniBooNE and LSND in $\nu_\mu \to \nu_e$ and $\overline{\nu}_\mu \to \overline{\nu}_e$ searches, as well as in $\nu_\mu$ and $\nu_e$ disappearance. We also discuss the $\nu_e$ and $\overline{\nu}_e$ disappearance in the context of the Gallium and reactor anomalies, where both null results and hints are observed.

§.§.§ Pion Decay-at-Rest Accelerator Experiments

The KARMEN (KArlsruhe Rutherford Medium Energy Neutrino) experiment was located at the highly pulsed
spallation neutron source ISIS of the Rutherford Laboratory (UK).  ISIS protons had an energy of 800 MeV and were delivered to the water-cooled Ta-D$_2$O target with a repetition rate of 50 Hz.  The time
structure of the ISIS protons (double pulses with a width of 100 ns separated by 325 ns) allowed a clear separation of $\nu_\mu$ events due to $\pi^+$ decay from $\bar\nu_\mu$ and $\nu_e$ events due to $\mu^+$ decay.

The KARMEN detector <cit.> was a segmented liquid scintillator calorimeter, located 17.7 m from the ISIS target at an angle 100 degrees relative to the proton beam. The active target consisted of 65 $\mbox{m}^3$ of liquid scintillator segmented into 608 modules with gadolinium-coated paper placed between modules for efficient detection of thermal neutrons.  KARMEN performed a search for $\bar\nu_\mu \rightarrow \bar\nu_e$ oscillations, analogous to the LSND search, using $p(\bar\nu_e,e^+)n$ and found measured rates agreed with background expectations <cit.>. The sterile neutrino 90% Confidence Interval (C.I.) obtained by the KARMEN measurement is shown in Fig. <ref> in relation to the LSND allowed regions. Although KARMEN did not see the LSND-like signal, it did not exclude the entirety of the LSND 99% allowed regions although it did strongly disfavour the larger $\Delta m^2 > 10 \text{ eV}^2$ solutions.

The KARMEN sterile neutrino 90% C.I. compared to other then contemporary experiments. Figure from <cit.>

§.§.§ Pion Decay-in-Flight Accelerator Experiments

Short-Baseline Experiments

While MiniBooNE's $\nu_e$ and $\bar{\nu}_e$ appearance searches in the Fermilab BNB have resulted in the observation of anomalous excesses, MiniBooNE BNB $\nu_\mu$ and $\bar{\nu}_\mu$ CC measurements have been relatively well understood, including a disagreement between data and the Monte Carlo prediction that has been attributed to cross-section effects and uncertainties <cit.>. As a result, MiniBooNE has been able to perform searches for $\nu_\mu$ and $\bar{\nu}_\mu$ disappearance, both exclusively, and inclusively, resulting in limits (using approximately the first half of data collected by MiniBooNE) as shown in Fig. <ref> <cit.>. The results complement those from prior short-baseline $\nu_\mu$ and $\bar{\nu}_\mu$ disappearance searches, from the CCFR <cit.> and CDHS <cit.> experiments.

MiniBooNE's 90% CL sensitivity (dashed line) and limit (solid line) for $\nu_{\mu}$ (top) and $\bar\nu_{\mu}$ (bottom) disappearance under a 3+1 scenario. Previous limits by CCFR (dark grey) and CDHS (light grey) are also shown. The figure is from <cit.>.

It should be noted that, in these searches, due to the lack of a near detector (ND), large flux and cross-section uncertainties limited MiniBooNE's sensitivity particularly in the high $\Delta m^2$ range, where (fast) oscillations are expected to lead to an overall normalization deficit, and are thus masked by systematic uncertainties on the overall $\nu_\mu$ and $\bar{\nu}_\mu$ CC rate normalization. More powerful searches were performed by MiniBooNE in combination with measurements by the SciBooNE experiment, as described next.

MiniBooNE/SciBooNE The SciBooNE experiment was located 100 m downstream from the BNB target and ran simultaneously with MiniBooNE from 2007-2008.  A simultaneous $\nu_\mu$ disappearance search under a (3+1) scenario was performed in the SciBooNE and MiniBooNE detectors with the BNB operating in forward horn current mode <cit.>.  A separate $\bar\nu_\mu$ disappearance search was performed in the SciBooNE and MiniBooNE detectors with the BNB operating in reverse horn current mode <cit.>. Exclusion (at 90% CL) regions for these two searches are shown in Fig. <ref>, and are consistent with exclusion limits from other accelerator-based neutrino experiments, including the short-baseline CCFR and CDHS experiments, and the long-baseline MINOS experiment <cit.>, as well as a MiniBooNE-only search limit.

MiniBooNE-SciBooNE 90% CL limits from a joint $\nu_\mu$ disappearance search (left) <cit.> and a joint $\bar\nu_\mu$ disappearance search (right) <cit.>.

MiniBooNE-NuMI In addition to neutrinos from the BNB, MiniBooNE has also been able to study neutrinos from the Fermilab-based NuMI beam, viewed by MiniBooNE at an off-axis angle of 6.3$^\degree$ <cit.>. MiniBooNE measured both $\nu_e$ and $\nu_{\mu}$ CCQE events from the NuMI off-axis beam (with no wrong-sign discrimination), as shown in Fig. <ref> (top). While an appearance or disappearance search with these data sets was not performed by MiniBooNE[The MiniBooNE-NuMI $\nu_e$ data set has been analyzed under a 3+1 appearance scenario in Ref. <cit.>, where it was shown that, despite the small observed data excess, the data showed no significant preference for oscillations over the null hypothesis.], the data was found to be in agreement with expectation. The NuMI off-axis beam at the MiniBooNE location is much higher in $\nu_e$ content than the BNB, as well as in neutrinos contributed from kaon decays in the beamline. Because of a higher intrinsic $\nu_e$ background, this data set was particularly limited in sensitivity to light sterile neutrino oscillations.

Top: Reconstructed energy distributions of $\nu_{\mu}$ (left) and $\nu_e$ (right) events in MiniBooNE from the NuMI off-axis beam. The figure is
from <cit.>. Bottom: Limit extracted from a fit to the $\nu_e$ distribution under a 3+1 hypothesis, from Ref. <cit.>.

Oscillation searches due to light sterile neutrinos are currently in development by the MicroBooNE collaboration.
Meanwhile, already-published MicroBooNE $\nu_e$ CC data sets <cit.> have been studied by the phenomenology community <cit.> and in combination with the MiniBooNE data by the MiniBooNE collaboration <cit.>.
Results from those studies, which involve a number of assumptions and approximations, are described in the following paragraphs.

Electron Neutrino Disappearance: It is possible

to search for anomalous $\nu_e \to \nu_e$ disappearance in MicroBooNE samples rich with intrinsic $\nu_e$ contributions <cit.>.
A preliminary analysis <cit.> of reported MicroBooNE results <cit.> finds a preference for $\nu_e$ disappearance in the various results (which are not statistically independent) at the ${\sim}2\sigma$ level.
This disappearance requires large mixing, $\sin^2\left(2\theta_{14}\right)\sim0.2-0.3$ with a mass-squared splitting of ${\sim}1.4$ eV$^2$, overlapping with the region preferred by the Gallium anomaly <cit.>.
However, a more recent analysis <cit.> incorporating systematic uncertainties from the collaboration's data release and accounting for neutrino energy reconstruction smearing finds this preference for $\nu_e$ disappearance to be statistically insignificant.

MicroBooNE recently reported their first electron neutrino data.
While much attention has been given to testing $\nu_e$ appearance in a primarily $\nu_\mu$ beam, it is also possible to look for $\nu_e$ disappearance due to the intrinsic $\nu_e$ contribution in the flux <cit.>.
The same intrinsic $\nu_e$ in the source was present in MiniBooNE as well, but given other backgrounds that can now be rejected by MicroBooNE, the signal to noise is significantly better with MicroBooNE than with MiniBooNE making MicroBooNE a great place to search for such a signal. An analysis of the various different MicroBooNE data sets <cit.> found that that several analyses, including the most sensitive one, show a slight dip that could be explained $\nu_e$ disappearance with a new frequency $\Delta m^2_{41}\sim1.4$ eV$^2$ and a large mixing angle of $\sin^2(2\theta_{14})\sim0.2-0.3$ <cit.>.
While the significance of this hint is low, $1-2\sigma$ depending on the analysis details, it interestingly is in the same region of parameter space hinted at by the Gallium anomaly which is now claimed at high significance <cit.> as well as another hint from the T2K near detector at low significance <cit.>.
These preferred parameters are in tension from solar neutrino data <cit.> and reactor antineutrino data <cit.> and are in strong conflict with data from the early universe <cit.>.
It is possible that the early universe constraint could be alleviated in light of the Hubble tension <cit.> which could be partially accommodated by a new light sterile neutrino <cit.>.
Future measurements by the rest of the short-baseline accelerator neutrino program can robustly test this scenario due to the variety of baselines.

MicroBooNE 3+1 Oscillations:  Ref. <cit.> has also provided the first analysis of MicroBooNE's $\nu_e$ results in the context of a light sterile neutrino. This includes both an analysis in the “simplified two-flavor” model, in which only $\nu_\mu \to \nu_e$ oscillations are taken into account, without significant $\nu_\mu \to \nu_\mu$ or $\nu_e \to \nu_e$ effects, consistent with MiniBooNE analyses prior to 2022 as well as the complete 3+1 model. It is found that, while MicroBooNE prefers the null hypothesis, i.e. no new sterile neutrino, the allowed regions from MiniBooNE at ${\sim}3\sigma$ are still allowed at the same confidence level by the MicroBooNE results, see left panel of Fig. <ref>.
Given its large dataset, the results of MicroBooNE's inclusive analysis <cit.> are found to be the most statistically powerful, while the CCQE sample provides a considerably weaker constraint.
Besides that, the MiniBooNE collaboration has performed a fit of both MiniBooNE and MicroBooNE CCQE data <cit.> to a 3+1 sterile neutrino model <cit.>, properly accounting for oscillations in the backgrounds.
As expected, the impact of adding MicroBooNE's CCQE sample to the fit has a marginal effect on the preferred regions for the 3+1 model, see right panel of Fig. <ref>.
While the two experiments' results are combined into a joint likelihood, the systematic uncertainties that the two experiments have in common, particularly those associated with their common neutrino flux from the Booster Neutrino Beam and the neutrino-nucleus cross section model, have not been combined.

The Combined Fit of MiniBooNE and MicroBooNE – The MiniBooNE collaboration has performed a joint fit of both MiniBooNE and MicroBooNE data to a 3+1 sterile neutrino model <cit.>. The fit is performed using both neutrino and antineutrino mode data from MiniBooNE's $\barparenb{\nu}_e$ and  $\barparenb{\nu}_\mu$ charged-current quasi-elastic (CCQE) analyses along with data from the recent MicroBooNE $\nu_e$ two-body CCQE analysis <cit.>. In MiniBooNE, this work considers the impact of $\barparenb{\nu}_e$ appearance and disappearance for signal events in the $\barparenb{\nu}_e$ CCQE channel, as well as $\barparenb{\nu}_\mu$ disappearance for signal events in the $\barparenb{\nu}_\mu$ CCQE channel and background events in the $\barparenb{\nu}_e$ CCQE channel. In MicroBooNE, this work considers the impact of $\nu_e$ appearance and disappearance for signal events in the $\nu_e$ CCQE channel. The $\nu_\mu$ background rate is very low in MicroBooNE's $\nu_e$ CCQE analysis, thus the impact of $\nu_\mu$ disappearance is negligible here. While the two experiments' results are combined into a joint likelihood, the systematic uncertainties that the two have in common (particularly those associated with their common neutrino flux from the BNB) have not been combined.

The joint fit maximizes a total likelihood function containing individual MiniBooNE and MicroBooNE likelihoods which account for both statistical and correlated systematic uncertainties in each experiment. This procedure finds a best fit point in sterile neutrino oscillation parameter space at $\sin^2(2\theta_{\mu e}) = 0.0316$ and $\Delta m^2 = 0.209\;{\rm eV^2}$, with $\Delta \chi^2 / {\rm dof} = 24.7/3$ with respect to the no-oscillation hypothesis. If one assumes Wilks' theorem for two degrees of freedom, closed contours exist in $\{\sin^2(2\theta_{\mu e}),\Delta m^2\}$ parameter space at the $4\sigma$ confidence level. Compared to a MiniBooNE-only fit, the addition of the MicroBooNE two-body $\nu_e$ CCQE result does not appreciably change the preferred regions in oscillation parameter space, as shown in Fig. <ref>.

Left: Preferred regions in the simplified two flavor model for MiniBooNE (gray), MicroBooNE inclusive analysis (blue), MicroBooNE CCQE analysis (red), KARMEN (green), and OPERA (brown) as indicated. For reference, the predicted sensitivity of SBN is also shown. Figure taken from Ref. <cit.>. Right: Preferred regions in the 3+1 model parameter space for a MiniBooNE-only fit and a joint fit including the MicroBooNE $\nu_e$ CCQE result.
Figure adapted from Ref. <cit.>.

NOMAD The NOMAD experiment <cit.>, which ran at CERN using protons from the 450 GeV SPS accelerator, employed a conventional neutrino beamline to create a wideband 2.5 to 40 GeV neutrino energy
source. These neutrinos were created with a carbon-based, low-mass, tracking detector located 600 m
downstream of the target. This detector had fine spatial resolution and could search for muon-to-electron and muon-to-tau oscillations. No signal was observed in either mode, and this experiment
set a limit to $\nu_\mu\rightarrow\nu_e$ appearance that excluded the majority of the 99% CL-allowed LSND region above 10 eV$^2$ at 90% CL, but had a significantly worse limit than KARMEN below 10 eV$^2$.

CCFR The CCFR experiment was carried out at Fermilab in 1984 <cit.>. The experiment made use of a narrow band beamline, with meson energies set to 100, 140, 165, 200, and 250 GeV, yielding $\nu_\mu$ and $\bar\nu_\mu$ beams that ranged from 40 to
230 GeV in energy. A two-detector experimental setup carried out a disappearance search, with the ND at 715 m and
the far detector (FD) at 1116 m from the center of the 352 m long decay pipe. The calorimetric detectors
were constructed of segmented iron with scintillator and spark chambers, and each had a downstream
toroid to measure the muon momentum. The data showed no evidence for a distance-dependent modification of the neutrino flux and ruled out oscillations of $\nu_\mu$ into any other single type of neutrino for $30<\Delta m^2<1000$ eV$^2$ and $\sin^2(2\theta)>0.02-0.20$.

CDHS The CDHS experiment <cit.> at CERN searched for $\nu_\mu$ disappearance with a two-detector design
of segmented calorimeters with iron and scintillator. The experiment used 19.2 GeV protons on a
beryllium target to produce mesons that were subsequently focused into a 52 m decay channel. The
detectors were located 130 m and 885 m downstream of the target. The experiment set a limit at 95% CL and set constraints that are comparable to the MiniBooNE $\nu_\mu$ disappearance limit described above, but extending to slightly lower $\Delta m^2$.

Long-Baseline Experiments
The phenomenology of sterile neutrino-driven oscillations at long baselines
consists of interference phenomena arising from at least two
distinct oscillation frequencies and several scale-determining mixing angles and phases.
In the context of a 3+1 model with the addition of one new flavor state $\nu_s$ and a new mass state $\nu_4$ extending the three-flavor PMNS matrix to a $4\times4$ unitary matrix, and assuming the parametrization $U=R_{34}R_{24}R_{14}R_{23}R_{13}R_{12}$, the extended matrix is written as:

\begin{equation}
U =
\begin{bmatrix}
U_{e1}    &  U_{e2}   &  e^{-i\delta_{13}}s_{13}c_{14}                                                                                & e^{-\delta_{14}}s_{14} \\
U_{\mu1}  & U_{\mu2}  & -e^{-i(\delta_{13}-\delta_{14}+\delta_{24})}s_{13}s_{14}s_{24} + c_{13}s_{23}c_{24}                                        & e^{-i\delta_{24}}c_{14}s_{24} \\
U_{\tau1} & U_{\tau2} & -e^{i\delta_{24}}c_{13}s_{23}s_{24}s_{34} + c_{13}c_{23}c_{34} - e^{-i(\delta_{13}-\delta_{14})}s_{13}s_{14}c_{24}s_{34}   & c_{14}c_{24}s_{34} \\
U_{s1}    & U_{s2}    & -e^{i\delta_{24}}c_{13}s_{23}s_{24}c_{34} - c_{13}c_{23}s_{34} - e^{-i(\delta_{13}-\delta_{14})}s_{13}s_{14}c_{24}c_{34}   & c_{14}c_{24}c_{34} \\
\end{bmatrix},
\end{equation}

which introduces three new mixing angles, $\theta_{14}$, $\theta_{24}$, and $\theta_{34}$, two new CP violating phases, $\delta_{14}$ and $\delta_{24}$, and one new linearly independent mass-splitting that is chosen to be $\Delta m^2_{41}$ in this case.
Neglecting $\Delta_{21}$ and using unitarity to rewrite any terms containing $U_{\alpha 1}$ or $U_{\alpha 2}$, the probability of $\nu_{\mu}$ survival in a 3+1 model can be written as

\begin{equation}\label{eq:SterileDisFull}
\begin{split}
P(\nu_{\mu} \rightarrow \nu_{\mu}) \approx 1 & - 4\left|U_{\mu 3}\right|^2\left(1 - \left|U_{\mu 3}\right|^2 - \left|U_{\mu 4}\right|^2\right)\sin^2\Delta_{31} \\
& - 4\left|U_{\mu 3}\right|^2\left|U_{\mu 4}\right|^2\sin^2\Delta_{43} -4\left|U_{\mu 4}\right|^2\left(1 - \left|U_{\mu 3}\right|^2 - \left|U_{\mu 4}\right|^2\right)\sin^2\Delta_{41}.
\end{split}
\end{equation}

Of particular interest are NC neutrino interactions, since as the three active flavors participate in the NC interaction at the same rate, the NC sample is insensitive to three-flavor oscillations between ND and FD. However, if a sterile neutrino exists, sterile-neutrino appearance would cause a depletion in the NC channel. The NC survival probability is thus defined as $1 - P(\nu_{\mu} \rightarrow \nu_{s})$. Similar to Eq. (<ref>), this can be written as

\begin{equation}\label{eq:SterileAppFull}
\begin{split}
1 - P(\nu_{\mu} \rightarrow \nu_{s}) \approx 1 & - 4\left|U_{\mu 3}\right|^2 \left|U_{s 3}\right|^2 \sin^2\Delta_{31} \\
& - 4\left|U_{\mu 4}\right|^2 \left|U_{s 4}\right|^2 \sin^2\Delta_{41} \\
& - 4\operatorname{Re}\left(Z\right)\left( \sin^2\Delta_{31} - \sin^2\Delta_{43} + \sin^2\Delta_{41} \right) \\
& - 2\operatorname{Im}\left(Z\right)\left( \sin2\Delta_{31} + \sin2\Delta_{43} - \sin2\Delta_{41} \right),
\end{split}
\end{equation}

where $Z = U^{*}_{\mu 4} U_{s 4} U_{\mu 3} U^{*}_{s 3}$.
The phenomenology of $\nu_{\mu}$-CC and NC disappearance driven by sterile neutrino oscillations is complicated at long baselines due to the interference of three-flavor oscillations and sterile oscillations, which does not occur at short baselines. To perform two-detector analyses typical of long-baseline accelerator neutrino experiments, effects at both short and long baselines must thus be understood.

For instance, in the case of the short-baseline oscillation, only probed by the ND, the oscillation probability for NC disappearance is approximately given by

\begin{equation}
1 - P(\nu_{\mu} \rightarrow \nu_s) \approx 1 - \cos^4\theta_{14}\cos^2\theta_{34}\sin^{2}2\theta_{24}\sin^2\Delta_{41},
\label{eqn:sblNC}
\end{equation}

where $\Delta_{ji} = \frac{\Delta m^2_{ji}L}{4E}$. For the typical beam neutrino energies and ND baselines, when $\Delta m^2_{41} < 0.05$, oscillations are not visible in the ND. Starting at $\Delta m^2_{41}$ $\sim$ 0.5 eV$^2$, oscillations begin to be visible at low energies in the ND, and as $\Delta m^2_{41}$ increases, the first oscillation maximum moves to higher energies. At sufficiently high $\Delta m^2_{41}$ values, the entire ND sees rapid oscillations that can no longer be resolved and are seen as a constant normalization shift described by

\begin{equation}
1 - P(\nu_{\mu} \rightarrow \nu_s) \approx 1 - \frac{1}{2}\cos^4\theta_{14}\cos^2\theta_{34}\sin^{2}2\theta_{24}.
\end{equation}

For $\nu_{\mu}$-CC at the ND, the oscillation probability can be approximated as

\begin{equation}
P(\nu_{\mu} \rightarrow \nu_{\mu}) \approx 1 - \sin^22\theta_{24}\sin^2\Delta_{41},
\label{eqn:sblCC}
\end{equation}

which behaves similarly to NC disappearance except it depends only on $\theta_{24}$, and in the rapid oscillation case the normalization shift is given by $(1/2)\sin^22\theta_{24}$.

However, when considering long baselines, terms oscillating at the atmospheric frequency cannot be neglected. Approximating the NC disappearance probability to first order in small mixing angles gives

\begin{equation}
\begin{aligned}
1 - P(\nu_{\mu} \rightarrow \nu_s) & \approx 1 - \cos^4\theta_{14}\cos^2\theta_{34}\sin^{2}2\theta_{24}\sin^2\Delta_{41} \\
& - \sin^2\theta_{34}\sin^22\theta_{23}\sin^2\Delta_{31} \\
& + \frac{1}{2}\sin\delta_{24}\sin\theta_{24}\sin2\theta_{23}\sin\Delta_{31}.
\end{aligned}
\end{equation}

In this expression, the first term is identical to the short-baseline approximation. The second and third terms both oscillate at the atmospheric frequency. If $\theta_{34}$ $>$ 0, the second term is non-zero, and if $\sin\delta_{24}$ and $\theta_{24}$ are non-zero, the third term will not be zero. In either case, this creates an oscillation dip visible at the FD regardless of the value of $\Delta m^2_{41}$. It is notable that this will happen even though NC disappearance cannot occur in the standard three-flavor paradigm. It should also be noted that the third term is CP-odd, since NC disappearance is effectively sterile neutrino appearance, so $\bar{\nu}$ data will also add to sterile sensitivity of long-baseline experiments.

For $\nu_{\mu}$-CC disappearance, expanding to second order in small mixing angles one finds

\begin{equation}
\begin{aligned}
P(\nu_{\mu} \rightarrow \nu_{\mu}) &\approx 1 - \sin^22\theta_{23}\cos2\theta_{24}\sin^2\Delta_{31}\\
&- \sin^22\theta_{24}\sin^2\Delta_{41},
\end{aligned}
\end{equation}

which can be rewritten as

\begin{equation}
\begin{aligned}
P(\nu_{\mu} \rightarrow \nu_{\mu}) &\approx 1 - \sin^22\theta_{23}\sin^2\Delta_{31} \\
& + 2\sin^22\theta_{23}\sin^2\theta_{24}\sin^2\Delta_{31} \\
& - \sin^22\theta_{24}\sin^2\Delta_{41}.
\label{eq:NuMuDisFull}
\end{aligned}
\end{equation}

The first term is the standard approximation for three-flavor $\nu_{\mu}$-CC disappearance. The second term also oscillates as a function of $\Delta m^2_{31}$, the atmospheric frequency, but it is driven by sterile mixing. Even at large $\Delta m^2_{41}$ values, this term does not enter into rapid oscillations. Thus, even for large mass splittings where the ND is well inside the rapid oscillation regime, the FD will still show shape variations in addition to normalization changes due to the terms oscillating at the atmospheric frequency with a magnitude that scales with $\sin^2\theta_{24}$. Furthermore, when considering $\nu_{\mu}$ disappearance alone, one could define an effective atmospheric mixing angle to account for both sterile and standard oscillations by combining the first two terms in Eq. (<ref>):

\begin{equation}
\sin^22\theta_{23}^{\mathrm{eff}} = \sin^22\theta_{23}\cos2\theta_{24}.
\end{equation}

While it may seem this would imply that the depth of the atmospheric dip would be insensitive to the large $\Delta m^2_{41}$ regime, it actually provides a constraint due to $\theta_{23}^{\mathrm{eff}}$ having been measured to be close to maximal. Due to the $\cos2\theta_{24}$ factor, a non-zero $\theta_{24}$ can only drive $\theta_{23}^{\mathrm{eff}}$ away from maximal.

In addition, experiments that can make precise measurements of  CC interactions can look for $\overset{(-)}\nu\!\!_e$-CC disappearance at the far detector following

\begin{equation}
\begin{aligned}
P(\overset{(-)}\nu\!\!_e \rightarrow \overset{(-)}\nu\!\!_e) &\approx 1 - \sin^22\theta_{13}\sin^2\Delta_{31} \\
& - \sin^22\theta_{14}\sin^2\Delta_{41},
\label{eq:NueDisFull}
\end{aligned}
\end{equation}

valid only if one assumes sterile-driven muon neutrino disappearance and electron neutrino appearance do not occur.

Finally, it is worth noting the near detectors of long-baseline experiments can be used to conduct short-baseline searches following the same methodologies described in Sec. <ref>

MINOS/MINOS+ The Main Injector Neutrino Oscillation Search (MINOS) experiment was a long-baseline neutrino oscillation experiment using the Neutrinos at the Main Injector (NuMI) neutrino beam and two detectors placed within a $735\,$km baseline. The NuMI beam is produced by collisions with a graphite target of protons accelerated to 120 GeV at Fermilab's Main Injector. The secondary products of these collisions, pions and kaons, are focused by two parabolic magnetic horns and eventually decay into muons and neutrinos inside a 675 m long decay pipe filled with helium. The muons are absorbed in the rock and neutrinos continue towards the 1 kton ND, 1 km downstream of the target, and beyond, towards the 5.4 kton FD. The relative positions of target and horns can be changed to tune the beam spectrum to lower or higher neutrino energies. The MINOS run concluded in 2012, with a total exposure of over $15\times10^{20}\,$protons-on-target (POT) in neutrino and antineutrino mode since the start of data taking in 2005. The MINOS+ experiment operated the MINOS detectors using the NuMI beam upgraded from 320 kW to 700 kW of beam power, part of the NOvA experimental setup. Instead of the low-energy configuration used for the MINOS run, MINOS+ used the NOvA medium-energy configuration, which for MINOS+ corresponds to a neutrino energy spectrum peaked around 7 GeV, as shown in Fig. <ref>.

The NuMI neutrino energy spectrum for the MINOS+ medium-energy tune, shown as the red solid line. The  spectrum shown in blue is obtained with the same tune at a 14 mrad offset from the beam axis (MINOS+ is on-axis). For comparison, the spectrum corresponding to the NuMI low-energy tune used by MINOS is shown as the gold histogram.

Exposure of the MINOS+ detectors to a beam peaked above the three-flavor oscillation maximum provided excellent sensitivity to new physics through precise measurements of muon neutrino disappearance between the ND and FD. MINOS+ operated from 2014 through 2016 and accumulated only neutrino mode data. A combined analysis of the $10.6\times10^{20}$ POT of MINOS neutrino data and $5.8\times10^{20}$ POT of MINOS+ neutrino data using a two-detector fitting technique placed stringent limits on sterile driven muon neutrino disappearance within a 3+1 model, as shown in Fig. <ref> <cit.>.

The MINOS and MINOS+ 90% Feldman-Cousins exclusion limit compared to the previous MINOS result <cit.> and results from other experiments. The Gariazzo et al. region is the result of a global fit to neutrino oscillation data <cit.>.

The null results from MINOS/MINOS+ are one of the primary drivers, along with IceCube results, of the large tension between appearance and disappearance data when attempting to explain observations purely through sterile neutrino mixing. This is further evidenced in combinations with reactor data looking for electron (anti)neutrino disappearance, as described in Sec. <ref>

NOvA The NuMI Off-Axis $\nu_e$ Appearance (NOvA) experiment is a long-baseline accelerator neutrino experiment based at Fermilab and the Far Detector Laboratory at Ash River, Minnesota. NOvA has as its primary goal to measure three-neutrino mixing parameters, including determination of the neutrino mass ordering, by looking for appearance of electron neutrinos or antineutrinos, and disappearance of muon neutrinos or antineutrinos, using the NuMI neutrino beam produced at Fermilab. This is accomplished by using two detectors separated by 810 km, placed 14 mrad off the NuMI beam axis. Due to the off-axis placement, the detectors sample a narrow range of neutrino energies between 1 and 4 GeV, peaked at 2 GeV as shown in Fig. <ref>. This configuration is chosen to drastically reduce feed-down of NC interactions of higher-energy neutrinos, which typically represent the dominant background to the measurement of  CC interactions in on-axis experiments.

The plot on the left displays the predicted neutrino energy spectra of the NuMI beam  during the run, including the on-axis spectrum sampled by MINOS+ as the dotted line, as well as the spectrum at different off-axis positions. The 14 mrad off-axis position of the NOvA detectors is identified by the red line. The plot on the right shows neutrino energy as a function of the parent pion energy and the angle between the neutrino and the decaying pion.

The 0.33 kton ND is located underground next to the MINOS ND hall at Fermilab, while the 14 kton FD is positioned at the surface in Ash River, Minnesota. Both detectors are composed of extruded 32-cell PVC modules filled with liquid scintillator. The cells are read out by 32-pixel avalanche photodiodes (APDs).
NOvA began collecting data in 2014 and has so far accumulated large samples in both neutrino-dominated and antineutrino-dominated modes.

NOvA placed constraints on sterile neutrinos via searches for differences in the rate of NC neutrino interactions between the Near and Far detectors. The analysis was based on  $6.05\times10^{20}$ protons-on-target taken in neutrino-dominated mode, and 95 NC candidates were selected at the Far Detector compared with $83.5\pm9.7(\rm{stat.})\pm9.4(\rm{syst.})$ events predicted assuming mixing only occurs between active neutrino species. Therefore, NOvA found no evidence of active-sterile neutrino mixing. Interpreting these results within a 3+1 model results in constraints on the sterile mixing angles of $\theta_{24}$ $<$ $20.8^{\circ}$ and $\theta_{34}$ $<$ $31.2^{\circ}$ at the 90% C.L. for $0.05$ eV$^2\leq \Delta m^2_{41}\leq 0.5$ eV$^2$, the range of mass splittings for which no significant oscillations over the ND baseline are expected <cit.>. The  energy spectrum of the NC selected candidates at the Far detector is shown along with a comparison of the allowed regions for the sterile matrix elements obtained by NOvA with similar constraints by SuperK and IceCube-DeepCore in Fig. <ref>.

Left: Calorimetric energy spectrum for NC neutrino candidates in the NOvA FD data compared with three-flavor prediction; Right: The NOvA 68% (dashed) and 90% (solid) Feldman-Cousins non-excluded regions (shaded) in terms of $|U_{\mu4}|^{2}$ and $|U_{\tau4}|^{2} $, where it is assumed that $\cos^2\theta_{14}=1$ in both cases, compared to SuperK and IceCube-DeepCore constraints.

NOvA has also reported results on the first search for sterile antineutrino mixing in an antineutrino beam, using an exposure of $12.51\times10^{20}$ protons-on-target from the NuMI beam at Fermilab running in antineutrino-dominated mode. NOvA observed 121 NC antineutrino candidates at the FD, compared to a prediction of $122\pm11\rm{(stat.)}\pm15\rm{(syst.)}$ assuming mixing between only three active flavors. Therefore, no evidence for $\bar{\nu}_{\mu}\rightarrow \bar{\nu}_s$ oscillations is observed.  In this case the 3+1 model constraints on the mixing angles are found to be $\theta_{24} < 25^{\circ}$ and $\theta_{34} < 32^{\circ}$ at the 90% C.L. for $0.05~\rm{eV}^2 \le \Delta m_{41}^2 \le 0.5~\rm{eV}^2$ <cit.>. The antineutrino energy spectrum at the FD and obtained allowed regions are shown in Fig. <ref>.

Left: Calorimetric energy spectrum for NC antineutrino candidates in the NOvA FD data compared with three-flavor prediction; Right: The NOvA 68% (dashed) and 90% (solid) Feldman-Cousins non-excluded regions (shaded) in terms of $|U_{\mu4}|^{2}$ and $|U_{\tau4}|^{2} $, where it is assumed that $\cos^2\theta_{14}=1$ in both cases, compared to SuperK and IceCube-DeepCore constraints.

T2K The Tokai-to-Kamioka experiment <cit.> is a neutrino oscillation experiment based in Japan. T2K's primary goal is to measure three flavor neutrino oscillation parameters. An intense source of (anti)neutrinos is produced at J-PARC and is directed toward a series of detectors placed 280 from the target, and a massive detector 295 away (Super-Kamiokande, SK <cit.>). The neutrino beam composition is predominantly muon neutrino flavor, with a small (0.5%) fraction of electron neutrinos <cit.>; the beamline elements can be configured to produce a predominantly antineutrino beam as well. The “near detectors” include the INGRID detector <cit.>, WAGASCI detectors <cit.>, and ND280 detector suite. These detectors observe interactions from slightly different neutrino energy spectrums, with INGRID's peak neutrino energy approximately at 1 (on-axis), ND280 at 0.6 (2.5 degrees off-axis), and WAGASCI at 0.8 (1.5 degrees off-axis).
T2K has operated since 2009 and has produced a series of results on light sterile neutrinos using the ND280 detectors and SK with subsets of the data taken. T2K plans to install improvements to the ND280 detector <cit.>, and, with improvements to the beamline, T2K will have further opportunities for updated or expanded analyses between detectors.

CC electron neutrino interactions were selected in ND280 to test for 3+1 disappearance, motivated by radiochemical experiments <cit.> and discrepancies in reactor experiments energy spectra <cit.>. Isolating CC candidates was challenging, with significant backgrounds to the selection from  photons (from, for example, the inactive material surrounding the active scintillator target). These backgrounds were controlled by dedicated photon selections and with additional systematic uncertainty estimation on production of $\pi^0$ from neutrino interactions on non-scintillator materials. The analysis also assumes no disappearance, and systematic uncertainty from the flux and cross section models were reduced using dedicated selections as is done for T2K three flavor oscillation analyses. T2K excludes regions of interest at 95% CL with 5.9e20 POT between approximately $\sinstt{ee}>0.3$ and $\dm{\rm eff} > \SI{7}{\square\electronvolt}$ <cit.>. An updated analysis of large mixing regions would be useful in light of recent MicroBooNE results <cit.>.

A search for light sterile neutrinos using SK has been performed in 2019 <cit.>, using 14.7e20 (7.6e20) POT in neutrino (antineutrino) mode. It is focusing on a 3+1 model where a single sterile neutrino is added and assumed to be mixing with $\nu_2$ and $\nu_3$ mass eigenstates through new parameters $\theta_{24}$ and $\theta_{34}$ (T2K has limited sensitivity to other parameters, that are thus neglected $\theta_{14}=\delta_{14}=\delta_{24}=\SI{0}{\degree}$).

The analysis consists of a simultaneous fit of the CC muon, electron and NC neutrino samples. While CC channels are mainly sensitive to the new mass splitting $\Delta m^2_{41}$ and to $\theta_{24}$, the NC channel measures the active neutrino survival probability and is also sensitive to $\theta_{34}$. The five CC analysis samples are the same as the one used in standard oscillation analyses <cit.>, but it has been the first time NC$\pi^0$ (single $\pi^0$ production where the pion decays and produces two Cherenkov rings) and NC $\gamma$-deexcitation samples are used in oscillation measurements.

Most of the systematic parameters are constrained in the same way as in the 3 flavor analysis (using e.g. ND constraints), but additional 30% normalization uncertainties are added for the new NC samples. The joint maximum-likelihood fits allows to draw exclusion limits in the $(\sinst{24}, \dm{41})$ plane: the most stringent limits on $\theta_{24}$ are obtained for $\dm{41} < \SI{3e-3}{\square\electronvolt}$. Similarly, limits in the $(\sinst{24}, \cos^2 \theta_{24} \sinst{34})$ plane have been obtained although more statistics in the NC samples and additional systematic studies are needed to further improve the measurement.
A comparison of the exclusion regions obtained with other existing constraints is shown in Fig. <ref>.

One limitation of work to date on T2K is the completeness of the assessment of interaction model uncertainties as applied to short-baseline analyses. T2K analyses so far assume no disappearance, however the interaction model systematic uncertainties are assessed based on external and ND280 measurements. Those measurements are placed close to production and therefore could be sensitive to a disappearance signal, potentially biasing a dedicated disappearance search.  T2K studied the possible impact of a subset of interaction model uncertainties on a ND280 disappearance result <cit.> and found it to be robust, but this does not consider a full re-assessment of where external data is used to inform the model. Current efforts in T2K cross-section measurements and the implementation of ab initio computations in the context of three-flavor analysis would greatly benefit such studies as well. Certain event selections, like coherent neutrino interactions or low $\nu$ selections <cit.> may have better theoretical understanding for a single detector analysis, but can be challenging to use due to statistics or acceptance for T2K.

In the future, joint analyses, including disappearance, appearance/disappearance may also be performed either between several near detectors (with different technologies and/or neutrino energy spectra) or benefiting from the complementary coverage of near versus far facilities. The upgraded near detector, for example, will have improved detection thresholds and sensitivity to neutrons, possibly enabling the selection of NC events at the near site.

The T2K 90% CL. exclusion limits on $\sinst{24}$ as a function of 41 compared to other experiments. Areas on the right are excluded.

OPERA and ICARUS
The Oscillation Project with Emulsion-tRacking Apparatus (OPERA) was a long-baseline accelerator neutrino experiment that sampled the CERN Neutrinos to Gran Sasso (CNGS) beam with a detector placed at the Gran Sasso laboratory (LNGS) 730 km from the neutrino production source, having as primary physics goal to observe $\numu\rightarrow\nutau$ appearance. This goal required the use of emulsion detection technology, such that the detector was a hybrid apparatus made of a nuclear emulsion/lead target complemented by electronic detectors.
The detector “target” region had a total mass of about 1.25 kt and was composed of two identical sections, each with 31 walls of emulsion cloud chamber bricks, interleaved by planes of horizontal and vertical scintillator strips used to select bricks in which a neutrino interaction had occurred. The exquisite spatial resolution afforded by the emulsion layers enabled identification of the characteristic “kink” due to decay of the final-state $\tau$ particle from a $\nutau$ interaction. The target region was complemented by a magnetized muon spectrometer.

Reconstruction of neutrinos interacting within the emulsion layers was conducted by automated scanning robots. Due to the high-energy tau production threshold ($E_{\nu_\tau} \gtrsim 3.5$ GeV), the CNGS neutrino flux was distributed at higher energies than other LBL experiments, as shown in Fig. <ref> <cit.>. OPERA collected CNGS beam data from 2008 to 2012, with an integrated exposure of $17.97\times10^{19}$ protons-on-target. A total of 19505 neutrino interaction events in the detector target were recorded by the electronic detectors, of which 5603 were fully reconstructed in the emulsion layers.

Fluxes of the different CNGS beam neutrino components at LNGS.

OPERA used its full data sample to conduct a search for sterile neutrino mixing within the context of a 4-flavor (3 active + 1 sterile) model based on both the $\nutau$ and $\nue$ appearance channels, with the $\numu$ CC disappearance channel not being considered given the low sensitivity to sterile mixing for that channel due to ambiguities with potential NC disappearance <cit.>.
Defining $\sin^{2}2\theta_{\mu\tau} = 4|U_{\tau4}|^2|U_{\mu4}|^2$, $\sin^{2}2\theta_{\mu e} = 4|U_{e4}|^2|U_{\mu4}|^2$, exclusion regions of $\Delta m^{2}_{41}$ versus $\sin^22\theta_{\mu\tau}$ and $\sin^22\theta_{\mu e}$ were computed and are shown in Fig. <ref>. The result is restricted to positive $\Delta m^{2}_{41}$ values since negative values are disfavored by results on the sum of neutrino masses from cosmological surveys <cit.>.
The results are consistent with no active-sterile neutrino mixing. For $\Delta m^{2}_{41} > 0.1 \textrm{ eV}^{2}$, the upper limits on $\sin^{2}2\theta_{\mu\tau}$ and $\sin^{2}2\theta_{\mu e}$ are set to 0.10 and 0.019 both for the case of Normal Ordering and Inverted Ordering. The values of the oscillation parameters $\left(\Delta m^{2}_{41} = 0.041 \textrm{ eV}^{2},\, \sin^{2}2\theta_{\mu e} = 0.92\right)$ corresponding to the MiniBooNE combined neutrino and antineutrino best-fit <cit.> are excluded with a $p$-value of $8.9\times10^{-4}$, corresponding to a significance of 3.3 $\sigma$.

Left: OPERA 90% C.L. exclusion region in the $\Delta m^{2}_{41}$ and $\sin^{2} 2\theta_{\mu\tau}$ parameter space for the normal (NH, solid line) and inverted (IH, dashed line) hierarchy of the three standard neutrino masses. The exclusion regions by NOMAD <cit.> and CHORUS <cit.> are also shown.
Right: OPERA 90% C.L. exclusion region in the $\Delta m_{41}^{2}$ and $\sin^{2}2\theta_{\mu e}$ plane is shown for the normal (NH, solid line) and inverted (IH, dashed line) hierarchy of the three standard neutrino masses. The plot also reports the 90% C.L. allowed region obtained by LSND <cit.> (cyan) and MiniBooNE combining $\nu$ and $\bar{\nu}$ mode <cit.> (yellow). The blue and red lines represent the 90% C.L. exclusion regions obtained in appearance mode by NOMAD <cit.> and KARMEN2 <cit.>, respectively. The 90% C.L. exclusion region obtained in disappearance mode by the MINOS and DayaBay/Bugey-3 joint analysis <cit.> is shown as green line. The black star ($\star$) corresponds to the MiniBooNE 2018 best-fit values for the combined analysis of $\nu$ and $\bar{\nu}$ data.

The ICARUS (Imaging Cosmic And Rare Underground Signals) experiment operated the T600 liquid argon time projection chamber detector, with a fiducial mass of 447 tons, exposed to the same CNGS beam as OPERA, and located at the same 730 km from the neutrino production target in the LNGS. The baseline and CNGS beam typical energy determine that ICARUS is primarily sensitive to sterile mixing in the region where $L/E \sim 36.5$ km/GeV. ICARUS conducted a search for sterile neutrino mixing by looking for excess $\nue$ appearance in the data sample collected in 2010 and 2011. This search assumed a simplified two-flavor model (one active, one sterile) and found no evidence of sterile neutrino oscillations, as detailed in Ref. <cit.>. However, it was pointed out soon after publication of the ICARUS results that the limits obtained on $\sin^{2}2\theta_{\mu e}$ by using a two-flavor approximation (which works well for a short-baseline experiment) in a long-baseline setup neglects sizable four-flavor effects, induced by the interference of the new large squared-mass splitting $\Delta m^2_{41}$ with the atmospheric one. The analysis also neglected contributions to the four-flavor oscillation probabilities arising from the intrinsic $\nue$ component of the CNGS beam. It is estimated that these four-flavor effects weaken the reported ICARUS constraints by up to a factor of 3 <cit.>.
The ICARUS T600 detector has since been moved to Fermilab, where it is being operated as the “Far” detector for Fermilab's Short-Baseline Neutrino program.

§.§.§ Reactor Neutrino Experiments

A conclusive way to test whether RAA is due to mixing between active and sterile neutrinos is by searching for baseline-dependent sterile neutrino induced spectral variations.
In this section, we summarize the results from the experiments that performed relative spectral search for 3+1  oscillations induced by active to sterile mixing.
Relative oscillation searches curtail the dependence on the reactor  spectrum model and detection efficiency by comparing the ratio of energy spectrum at different baselines to the corresponding predicted ratio under oscillation hypothesis.
Baseline-dependent spectral measurement is done either by using a segmented detector, a single-volume movable detector, or by use of multiple detectors placed at different baselines.
At reactor  energies, since eV-scale neutrinos induce oscillations at meter scales, purpose-built experiments typically place their detectors as close to the reactor as possible, ideally $<$10 m.
Furthermore, experiments built to search for $\theta_{13}$ are also sensitive to  sterile neutrino oscillations but at lower oscillation frequencies, typically for $\dm{13}<0.1 \mathrm{eV}^2$.

All the experiments discussed in this section use scintillator detectors and the inverse beta decay interaction process, where a characteristic signature is provided by the timing and spatial coincidence between prompt positron annihilation signal and delayed neutron capture signal producing two separate flashes ($\Delta t \sim \mu$s).
Refer to Sec. <ref> for more details on the IBD mechanism.
Experiments also typically use an external neutron capture agent (other than hydrogen), either Gadolinium or $^6$Li, to increase IBD efficiency by increasing the neutron capture cross-section.
$^6$Li has an added benefit of constraining the spatial extent of the delayed signal since the capture products have high $dE/dx$ in scintillators.

In addition to the delayed coincidence, experiments also use various combinations of active and passive shielding, detector segmentation, and pulse shape discrimination (PSD) for background reduction.
PSD capable scintillators generate pulse shapes that are particle-dependent and could in principle be used for both prompt and delayed selection.
Whenever possible, experiments also measure reactor-off data to constrain reactor-uncorrelated IBD-mimicking backgrounds.
Following is a discussion of the individual reactor experiments that searched for sterile neutrino induced oscillations.
A list of purpose-built experiments and their features pertinent to eV-scale oscillation search can be found in the Tab. <ref>.

Experiment     Baseline (m)     Reactor     Reactor     Detector     Target     Sterile $\nu$

Type     Power (MWth)     Size         Search Strategy

DANSS <cit.>      11–13     LEU     3000     1 m$^3$     Segmented PS with Gd coating     Multi-Site
NEOS <cit.>       24     LEU     2800      1 m$^3$     Single-volume GdLS + PSD     Single-Site
Neutrino-4 <cit.>     6–12     HEU     90      2 m$^3$     Segmented GdLS     Multi-Site/Zone
PROSPECT <cit.>      7–9     HEU     85     4 m$^3$      Segmented $^6$LiLS + PSD     Multi-Zone
STEREO <cit.>      9–11     HEU     57     2 m$^3$     Segmented GdLS + PSD     Multi-Zone

Details of the short-baseline experiments that have been specifically built to search for eV-scale sterile neutrinos. Short baselines are preferable for oscillation searches involving eV-scale sterile neutrinos. HEU reactors are typically $\sim 10x$ smaller than commercial LEU power reactors and are preferred for eV-scale oscillation searches since they have smaller source positions oscillation washout. Detector segmentation, capture agent, and pulse shape discrimination capabilities are highly beneficial for background reduction.

Short-Baseline Experiments

DANSS The DANSS (Detector of AntiNeutrinos based on Solid State Scintillator) experiment <cit.> consists of 1 m3 highly segmented, plastic scintillator detector with each scintillator bar wrapped in Gd-loaded reflector.
Light and signal readout is performed using wavelength-shifting fibers and a combination of pixel photon detectors and PMTs operated at room temperature.
The experiment samples from a 3 GW$_{th}$ low enriched uranium (LEU) reactor at Kalinin Nuclear Power Plant in Russia.
As opposed to most other reactor experiments, the DANSS detector is located below the reactor core, benefiting from its 50 meter water equivalent (m.w.e) overburden.
The detector is placed on a lifting platform which enables measurements at the baselines of 10.9 m, 11.9 m, and 12.9 m.
Due to its close proximity to the reactor core, the detector measures a high flux of around 4000 events per day.

The DANSS detector was commissioned in  2016 and started taking data in October 2016. The experiment has collected around 4 million signal events over three years at three different positions.
The biggest source of background that constitutes 13.8% of the IBD signal events is the accidental coincidence background, along with non-negligible cosmogenic backgrounds and IBD backgrounds from neighboring reactors. The data excludes a large area in sterile neutrino parameter space ($\Delta m^{2}$, $\sin^{2}{2\theta_{ee}}$) and most interestingly excludes the best fit point of RAA and Ga experiments at more than $5\sigma$ level as shown in Fig. <ref>. Additionally, the IBD positron spectrum was measured and compared to the simulated Huber-Mueller spectrum.  Although the measured spectrum disagrees with DANSS's Huber-Mueller-derived prediction, final oscillation measurements are largely independent of choice of models due to the ratios taken at different baselines <cit.>.

Left panel: Exclusion curves at 90% C.L. (filled area) and 90% C.L. sensitivity contours (dashed line) are shown. Expected regions from
RAA and GA are also shown.
Right panel: Exclusion curves for 3+1 neutrino oscillations in the  $\sin^{2}{2\theta_{14}}$, $\Delta m^{2}_{14}$ parameter space.

NEOS The Neutrino Experiment for Oscillation at Short baseline (NEOS) aims to search for light sterile neutrinos by detecting electron antineutrinos ($\bar{\nu_{e}}$) from a reactor at a very short baseline. The NEOS detector consists of a liquid scintillator (LS) target ($\sim$1000 liter with $\sim$0.5% Gd loading by weight), two buffer zones filled with mineral oil where 19 8 inch PMTs per a zone are attached, and a muon veto system. NEOS had been installed in the tendon gallery of the 5$^{\mathrm{th}}$ reactor of the Hanbit Nuclear Power Plant in Korea for two periods, from Aug. 2015 to May 2016 (NEOS-I) and from Sept. 2018 to Oct. 2020 (NEOS-II). For the NEOS-II, target GdLS was newly produced and there was a minor modification in the muon veto system. With a 24 m baseline from the core of the reactor (2.8 GW$_{th}$ in 100%) about 2 000 inverse beta decay (IBD), $\bar{\nu_{e}} + p \rightarrow e^+ + n$, events per day are observed. Signal to background ratio is about 22 thanks to relatively good overburden ($\sim$20 m.w.e.). Background contribution from the nearest neighboring reactor ($d = 256$ m) is found to be less than 1% of the total $\bar{\nu_e}$ flux from the 5$^{\mathrm{th}}$ reactor. Calibration data using radioactive sources, $^{22}$Na, $^{137}$Cs, $^{60}$Co, PoBe, and $^{252}$Cf, had been taken regularly.

NEOS-I <cit.> using 180 (46) live days of reactor-ON (-OFF) data excluded the RAA best fit at 90% C.L. by comparing the prompt energy spectrum of Daya Bay where the sterile neutrino oscillation effect is averaged out. The best-fit was found at (sin$^2 2\theta, \Delta m^2)$ = (0.05, 1.70 eV$^2$) with $\chi^2$/NDF for 3$\nu$ and 4$\nu$ are 64.0/61, and 57.5/59, respectively. The corresponding p-value is estimated to be 22%. The well-known “5 MeV excess” was clearly observed as well in NEOS-I. Recently a joint analysis <cit.> between NEOS-I and RENO was performed, yielding a slightly improved result beyond its previous result using early Daya Bay data.  More details on the NEOS-I and RENO joint analysis is discussed in a following sub-section (Sec. <ref>).

Neutrino-4 Neutrino-4 is the only reactor neutrino experiment performing relative spectral comparison that has reported evidence of oscillations.
The detector is a 1.8 m$^3$ Gd-doped liquid scintillator detector divided into 10 rows each row consisting of 5 sections each of size 0.225 m $\times$ 0.225 m $\times$ 0.85 m.
The detector samples from the 57 MW$_{th}$ SM-3 reactor, and HEU in Dimitrovgrad, Russia.
A movable platform enables the detector to sample baselines from 6 m–12 m.
The short reactor on (off) cycles of 8–10 (2-5) days enables the experiment to perform rapid signal and background measurements.

Left: 1$\sigma$ (blue), 2$\sigma$ (green), and 3$\sigma$ (yellow) suggested regions from the Neutrino-4 sterile neutrino oscillation search.  Right: $L/E$ distribution of background-subtracted IBD rates reported by the Neutrino-4 Experiment.  Plots from <cit.>

The detector collected data for five years from 2016–2021 with $\sim$300 events/day.
The latest analyzed dataset consists of 720 (860) calendar days of reactor on (off) data with a signal-to-background ratio of $\sim$0.5.
The experiment performed a sterile neutrino search in the L/E space and observed evidence for 2.9$\sigma$ neutrino oscillation effect with the best-fit at 14=7.3$\pm$0.13(stat) $\pm$1.16(syst) and $\sin^2\theta_{14}$=0.36$\pm$0.12(stat) (Fig. <ref>).
It is important to point out that these results are widely debated by several groups <cit.>.
Additionally, PROSPECT and STEREO experiments <cit.> also individually disfavor the best-fit point and a significant portion of the Neutrino-4 suggested parameter space.

PROSPECT The Precision Reactor Oscillation and SPECTrum (PROSPECT) experiment <cit.> is a U.S.-based reactor neutrino experiment installed at short baselines of 6.7–9.2 m from the 85 MW High Flux Isotope Reactor (HFIR) at Oak Ridge National Laboratory.
The 4-ton 6Li-loaded PSD-capable liquid scintillator <cit.> detector is composed of a two-dimensional grid of
optically isolated segments <cit.> as shown in Fig. <ref>.
Low mass, highly reflective, rigid separators <cit.> were used to achieve segmentation and photomultiplier tubes (PMTs) enclosed in mineral oil-filled acrylic housings were installed at either ends of the segments for signal readout.

The detector shielding was optimized based on the neutron and $\gamma$ ray background measurements performed at HFIR <cit.> and consisted of top-heavy hydrogenous shielding to reduce cosmogenic backgrounds and a fixed local lead shield to mitigate reactor-specific backgrounds.
PROSPECT was the first on-surface reactor neutrino experiment to achieve a signal-to-background (S:B) > 1, thanks to the high background suppression enabled by the detector design.

Left: Schematic of the PROSPECT detector, shielding, and electronics in the HFIR Experiment Room. The center of the detector is located at a distance of 7.84 m from the center of the reactor core (red). Right: Cutout of the PROSPECT active segmented detector enclosed in the containment and shielding. Also shown is the vertical local lead shield adjoining the reactor pool wall that was installed to mitigate the reactor-related backgrounds.

PROSPECT was commissioned in February 2018 and started collecting data in March.
The first oscillation search result was published in 2018 <cit.> with a relatively small dataset of 33 (28) reactor on (off) live days which was followed by a result with a longer dataset consisting of 96 (73) reactor on (off) live days composed of $>$ 50k signal events in 2020.
The oscillation search was done by performing relative comparison of baseline-dependent spectra minimizing the dependence on the reactor neutrino model <cit.>.

A combination of the compact HFIR reactor core (a cylinder of diameter 0.435 m and height of 0.508 m) with the fine detector segmentation (14.5 cm $\times$ 14.5 cm cross-section) enabled high sensitivity to oscillation frequencies of $\Delta m^2>1$ eV$^2$ .
PROSPECT observed no statistically significant indication of to sterile neutrino oscillations.
Using Feldman-Cousins technique <cit.>, PROSPECT excluded a significant portion of the RAA suggested parameter space at 95% CL and the best-fit point at 2.5 $\sigma$ as shown in Fig. <ref>.
A complementary technique using the Gaussian CLs method <cit.> also provides similar exclusion.

Results from PROSPECT's search for $\anue$ to sterile neutrino oscillations. Exclusion contours were drawn using Feldman-Cousins (black) and Gaussian CLs (red) methods. Green and yellow bands show the 1 $\sigma$ and 2 $\sigma$ PROSPECT sensitivities to the sterile neutrino oscillations. Also shown for comparison are the RAA suggested parameter space and the best-fit point. Figure from <cit.>

The PROSPECT detector was decommissioned in 2020 after an unexpected HFIR down time.
During the course of data taking, it was also observed that a number of segments have slowly lost functionality due to failures of PMT electronics induced by the ingress of the liquid scintillator into the PMT housings.
Due to these reasons, the results shown above–corresponding to the full PROSPECT dataset–are primarily dominated by statistics.

However the collaboration is carrying out two major analysis modifications by leveraging distinctive detector features.

a) Splitting the data taking period into multiple time-frames and b) allowing the use of segments with single PMTs.
The former would allow for increase in statistics by allowing segments with non-functioning PMTs to be used for part of the dataset and the latter increases the S:B by improving particle identification and consequently reduces the backgrounds.
Each of the analysis improvements are individually expected to increase the effective statistics by $\sim$50% and enable considerable improvement in the statistical power of the PROSPECT's oscillation search.

STEREO The STEREO experiment <cit.> is a Gd-doped liquid scintillator detector located at the ILL in Grenoble, France. The detector uses a high-flux 58 MW$_{th}$ research reactor that consists of highly enriched Uranium. The detector is located 9 m away from the reactor core and is segmented detector where position of different segments of detector serves as a different baselines varied between 9-11 m as shown in Fig. <ref> (Left).

The IBD signal events are selected using cut-based approach where selection cuts on energy and time variables are optimized by trading off between detection efficiency and background rejection. Furthermore, antineutrino signal rates separated from the remaining background using Pulse Shape Discrimination (PSD) where pulses generated from neutrons have longer tails compared to that of gamma. Hence, a ratio of the pulse tail to total charge is used to mitigate neutron-related background.

The test of sterile neutrino oscillations is performed using ratios of energy spectrum at six different segments to that of first segment and therefore making the measurement independent of absolute normalization and of the prediction of the reactor spectrum. With 119 (211) days of reactor-on (-off) data, STEREO found no evidence for sterile neutrino oscillations and the results are compatible with null oscillation hypothesis. The data excluded RAA best fit point at $>4\sigma$ level and exclusion curve with 3+1 neutrino oscillations scenario is shown in Fig. <ref> (Right).

(Left) STEREO setup. 1–6: target cells (baselines from core: 9.4–11.1 m). (Right) Exclusion contour (red) and exclusion sensitivity (blue) at 95% C.L. Overlaid are the allowed regions of the RAA (grey) and its best fit point (star).

Medium-Baseline Experiments

Daya Bay Daya Bay's unique configuration makes it an excellent experiment to search for sterile-active neutrino mixing <cit.>. In this experiment, electron antineutrinos emitted from six $2.9$ GW$_\mathrm{th}$ nuclear reactors are detected in eight identically-designed antineutrino detectors (ADs) placed underground in two near experimental halls (EHs) and one far hall.
The two near halls, EH1 and EH2, are located $\sim$350-600 m away from the reactors, whereas the far hall, EH3, is located $\sim$1500-1950 m from the reactors.

Daya Bay's latest constraints in the $(\sin^2 2\theta_{14},\Delta m^2_{41})$ parameter space, obtained with a 1230-day data set, are shown in Fig. <ref> <cit.>. Two complementary analysis methods are used to set the exclusion contours, one relying on the Feldman-Cousins (FC) frequentist approach and the other on the CL$_s$ approach. Daya Bay is most sensitive to $\sin^2 2\theta_{14}$ in the $10^{-3}$ eV$^2$ $\lesssim |\Delta m^2_{41}| \lesssim 0.3$ eV$^2$ region, where a distortion from the standard three-neutrino oscillation framework would be visible through a relative comparison of the rate and energy spectrum of reactor antineutrinos in the different EHs. For $|\Delta m^2_{41}| \gtrsim 0.3$ eV$^2$, the oscillations are too fast to be resolved, and the sensitivity arises primarily from comparing the measured rate with the expectation, resulting in less stringent limits.

Daya Bay ceased operations in December 2020 after collecting data for over 3000 days. Throughout this time it accumulated the largest sample of reactor antineutrinos to date, consisting of more than six million events. This sample is still being analyzed and final results are expected to be released by early 2023. The sizable increase in statistics, combined with potential reductions in systematic uncertainties, implies that significant improvements over the existing limits are expected. The new constraints will likely remain the best in the world for the foreseeable future in the $|\Delta m^2_{41}| \lesssim 0.3$ eV$^2$ region, which no experiment in the horizon is expected to cover at the time of writing.

Left: Feldman-Cousins (FC) 90% C.L. and 90% CL$_s$ exclusion regions from an oscillation analysis of 1230 days of Daya Bay data. The dashed red line shows the 90% C.L. median sensitivity along $1\sigma$ and $2\sigma$ bands. The excluded region for the original Bugey-3 limit <cit.> is shown in green, while the resulting CL$_s$ contour from Daya Bay and its combination with the reproduced Bugey-3 results are shown in grey and black, respectively. From <cit.>.  Right: RENO's 95 % C.L. exclusion contour for the oscillation parameters $\sin^2 2\theta_{14}$ and $|\Delta m^2_{41}|$. The black solid contour represents an excluded region obtained from spectral distortion between near and far detectors. The green shaded band represents expected 1$\sigma$ exclusion contours due to a statistical fluctuation. The blue dotted contour represents its median.

Double Chooz Double Chooz <cit.> consists of two nearly identical gadolinium-doped liquid scintillator detectors located close to the nuclear power plant comprising two 4.25 GW nuclear reactors.
The near (far) detectors are located underground at an overburden of 120 m (300 m) at a distance of 469 m and 355 m (1115 m and 998 m) from the two reactor cores.
The detector-reactor locations are such that the relative contributions to both the detectors from the reactors are very similar which helps reduce reactor-related uncertainties.

This analysis includes three datasets amounting to a total 5 year long dataset.
The first (FD-I) dataset consists of 455.21 days of livetime collected with the far detector before the commissioning of the near detector.
The second (FD-II) and third (ND) datasets are collected during the same period of time and consists of 362.97 days and 257.96 days of livetime respectively.
The livetime for ND is lower than FD-II because of the higher muon rate causing larger deadtime in the near detector.
In order to obtain a measurement independent of absolute flux predictions, the experiment directly compares the event rates measured in the two identical detectors which helps in canceling most of the reactor flux and detection efficiency related uncertainties.
The experiment does not see any indications of sterile neutrino oscillations and set exclusion limits in similar regions of 14 as the other $\theta_{13}$ experiments.

RENO The RENO collaboration has reported a search result for light sterile neutrino oscillations.
The search is performed using six 2.8 GW$_{\text{th}}$ reactors
and two identical detectors located at 294 m (near) and 1383  m (far), respectively, from the center of six reactor cores at the Hanbit Nuclear Power Plant Complex in Yonggwang.
The reactor flux-weighted baseline is 410.6 m for the near detector and 1445.7 m for the far detector, respectively.
The near (far) underground detector has 120 m (450 m) of water equivalent overburden.

A spectral comparison between near and far detectors was performed to search for reactor $3+1$ light sterile neutrino oscillations

The RENO sterile analysis uses roughly 2200 live days of data taken in the period between August 2011 and February 2018 amounting to 850 666 (103 212) $\overline{\nu}_e$ candidate events

in the near (far) detector.

The details of pull terms and systematic uncertainties are described in Ref. <cit.>.

Exclusion regions at 95 % confidence level are set for $\Delta \chi^2 > 5.99$ and are shown in Fig. <ref>.
Exclusion contours obtained using the Gaussian CL$_s$ method <cit.> show negligible difference with the pictured $\Delta\chi^2$ method.
Fig. <ref> shows the 95 % C.L. exclusion contour and median sensitivity including 1$\sigma$ band due to statistical fluctuations.

Comparison of the exclusion limits.
The right side of each contour indicates excluded region.
The black solid line represents the 95 % C.L. exclusion contour using spectral distortion between near and far spectra.
For the comparison, Daya Bay's <cit.> 90 % C.L. (green), Bugey's <cit.> 90 % C.L.(blue), KARMEN + LSND <cit.> 95 % C.L.(magenta)  and NEOS's <cit.> 90 % C.L.(brown) limits on $\overline{\nu}_e$ disappearance are also shown.
The blue shaded area represents a 95 % C.L. allowed region from Planck data <cit.>.

The limit of $\sin^2 2\theta_{14}$ is mostly determined by a statistical uncertainty, while the systematic uncertainties dominate in the $|\Delta m^2_{41}| \lesssim 0.06$  eV$^2$.
The uncertainty of background is a dominant systematic source in the $0.003 \lesssim |\Delta m^2_{41}| \lesssim 0.06$  eV$^2$ region, and the energy-scale uncertainty  is a major limiting factor in the $|\Delta m^2_{41}| \lesssim 0.008$  eV$^2$ region.
The RENO result provides the most stringent limits on sterile neutrino mixing at $|\Delta m^2_{41}| \lesssim 0.002$  eV$^2$ using the $\overline{\nu}_e$ disappearance channel.
Adding data taken since 2018 and reducing the above systematic uncertainties will improve the results significantly.

Joint Fits of Reactor Neutrino Experiments

As discussed above, the conclusive way to test whether RAA is due to active-sterile mixing is by searching for sterile neutrino induced spectral variations as a function of baseline.
With the exception of the Neutrino-4 experiment, no experiment has claimed to observe statistically significant hints of oscillations.
Nonetheless, modest hints of oscillation the other experiments have been reported, and it is worth considering global fits of these datasets to determine build a broader context from their individual results.
A joint analysis <cit.> using Monte Carlo simulations with data from a combination of the short-baseline reactor neutrino experiments DANSS, NEOS, Neutrino-4, PROSPECT, and STEREO shows that the combination of these datasets is statistically compatible with the three neutrino model.
The exclusion curve from this joint fit is shown in Fig. <ref>.

Left: 2$\sigma$ CL Feldman-Cousins (FC) exclusion curve of the combined reactor neutrino data. Also overlayed are the 2$\sigma$ exclusion FC curves from solar data, $\nu_{e}-^{12}$C scattering data from LSND and KARMEN, 95% C.L. exclusion from the KATRIN experiment, and 2$\sigma$ gallium suggested region. Right: 1$\sigma$ , 2$\sigma$, and 3$\sigma$ CL FC suggested regions for a combination of gallium anomaly and relative reactor measurements. The suggested regions are at higher  mainly driven by the BEST and Neutrino-4 experiments. Figures from Ref.<cit.>.

Although a combination of relative reactor neutrino datasets is incompatible with sterile neutrino induced oscillations, they are in good agreement with the Gallium Anomaly.
This is primarily driven by the data from BEST and Neutrino-4 experiments which prefer oscillations with 14$> 5 \textrm{ eV}^2$–a region where other reactor experiments have minimal sensitivity.
Experiments with sensitivity in the 14$ > 5 \textrm{ eV}^2$ are needed to fully address these suggested regions.
This could be achieved by a combination of upcoming reactor experiments (see Sec. <ref>) and $\beta$-decay experiments (see Sec. <ref>, <ref>).

Joint Analysis with Accelerator Experiments
It is impossible for a single experiment to cover all the parameter space of interest to experimentalists and phenomenologists, which motivates the undertaking of joint analyses carried out by the members of the experimental collaborations that properly treat systematic uncertainties and their correlations. Of particular interest is the combination of Daya Bay's data with those of other reactor experiments operating at shorter baselines to cover a wide range of $|\Delta m^2_{41}|$ values. A case in point is the joint fit of the Daya Bay and Bugey-3 data presented in Ref. <cit.> that results in the black contour of Fig. <ref>. Furthermore, powerful constraints on sterile-driven neutrino oscillations can also be extracted from combining data from reactor experiments with data from long-baseline accelerator experiments. In a 3+1 scenario, reactor experiments primarily measure $|U_{e4}|^2 = \,\,\sin^2\theta_{14}$ through electron antineutrino disappearance, while long-baseline accelerator experiments are typically most sensitive to $|U_{\mu4}|^2 = \,\,\sin^2\theta_{24}\cos^2\theta_{14}$ through measurements of muon (anti)neutrino disappearance. The product of the two matrix elements provides the amplitude of short-baseline electron neutrino appearance in a primarily muon neutrino source:

\begin{equation}
P_{\overset{(-)}\nu\!\!_\mu\rightarrow\overset{(-)}\nu\!\!_e}^{SBL} = \,4|U_{e4}|^2|U_{\mu 4}|^2\sin^2\left(\frac{\Delta m^2_{41}L}{4E}\right), \label{eq:PmueU}
\end{equation}

where $4|U_{e4}|^2|U_{\mu4}|^2 = \sin^22\theta_{14}\sin^2\theta_{24} \equiv \sin^22\theta_{\mu e}$.
Possibly the most representative examples of this type of combination are found in Refs. <cit.>, where constraints on $\theta_{14}$ from the Daya Bay and Bugey-3 experiments are combined with constraints on $\theta_{24}$ from the MINOS/MINOS+ experiments to constrain the effective mixing parameter $\sin^{2} 2\theta_{\mu e}$. This work culminated in the most stringent constraints to date from disappearance searches on active to sterile neutrino oscillations, and probed the parameter space allowed by the LSND and MiniBooNE anomalies, as shown in Fig. <ref>.

Comparison of the MINOS, MINOS+, Daya Bay, and Bugey-3 combined 90% $\mathrm{CL_s}$ (left) and 99% $\mathrm{CL_s}$ (right) limits on ${\sin}^22{\theta}_{\mu e}$ to the LSND and MiniBooNE 90% and 99% C.L. allowed regions, respectively. The limit also excludes the 90% and 99% C.L. regions allowed by a fit to global data by Gariazzo et al. where MINOS, MINOS+, Daya Bay, and Bugey-3 are not included <cit.>, and the 90% and 99% C.L. regions allowed by a fit to all available appearance data by Dentler et al. <cit.> updated with the 2018 MiniBooNE appearance results <cit.>.

As highlighted in Ref. <cit.>, there are attractive opportunities in combining Daya Bay's data with other current and future experiments, such as PROSPECT, STEREO, NEOS, and JUNO-TAO. It is worth noting that the Daya Bay collaboration plans to publicly release its full data set once all final results have been released <cit.>, allowing such combinations to occur even well after the collaboration has dissolved. Similarly, following end of data taking in 2016, the MINOS/MINOS+ $\mathrm{CL}_s$ surfaces remain available for use in future combinations.

§.§.§ Atmospheric Neutrino Experiments

IceCube Neutrino telescopes, such as the IceCube Neutrino Observatory in the South Pole, play a unique role in searches for new physics associated with the short-baseline anomalies <cit.>.
For atmospheric neutrino oscillation baselines, the $L/E$ range of the short-baseline anomalies corresponds to neutrino energies of order 100 GeV-1 TeV.
Because this closely matches the peak of the detected neutrino flux at IceCube, it is natural to expect any Lorentz invariant (LI, implying scaling as $L/E$) phenomenon connected with the short-baseline anomalies <cit.> may be meaningfully testable at IceCube.

Sensitivity to certain models is even further enhanced by fortuitous features of this energy range, notably: 1)  For sterile neutrinos with 0.1-10 eV$^2$ mass splittings, a resonant matter effect leads to dramatic enhancement of oscillations at $\sim$1 TeV, leading to sensitivity far exceeding that for vacuum-like oscillations <cit.>; 2) Any new physics model that invokes non-LI effects scaling positively with $E$ will be enhanced at IceCube relative to all other experiments at lower $L/E$ - these include anomalous decoherence <cit.> and Lorentz violation <cit.> models; 3) This energy range offers unique access to the $\nu_\tau$ appearance sector <cit.>.
Other notable features of the IceCube event sample are its “broadband” nature, covering five decades of energy of atmospheric neutrinos (10 GeV-1 PeV) and 3.5 decades of baselines (20-12,750 km);  high statistical precision, owing to the high total exposure of a billion ton detector operating stably for ten years collecting 70,000 atmospheric neutrino events per year, and well controlled cross section uncertainties due to the predominance of deep-inelastic scattering interactions.
These features have enabled world-leading sensitivity to the parameters governing three-flavor oscillations as well as non-standard oscillation models including sterile neutrinos <cit.>, tests of low-energy manifestations of quantum gravity <cit.>, neutrino decay <cit.>, and non-standard interactions <cit.>.

Present generation light sterile neutrino searches at IceCube
IceCube has made powerful sterile neutrino searches in both high ($\geq 400$ GeV) and  low ($\leq$ 60 GeV) energy ranges.
The former targets the matter resonance <cit.> expected for $\Delta m^2\sim$ ${\cal O}$(1 eV$^2$) splittings, and is one of the world's most sensitive in the $\nu_\mu$ disappearance channel at eV$^2$-scale mass splittings.
The latest generation analysis <cit.> uses a sample of 305735 reconstructed $\nu_\mu$ events and excludes mixing angles down to $\sin^2 2\theta_{24} \sim$ 0.02 at $\Delta m^2\sim 0.2$ at 99% CL.
At 90% CL the analysis yields a closed contour that may be interpreted as a statically weak hint of a signal, with best-fit point at $\sin^2 2\theta_{24}$ = 0.10 and $\Delta m^2_{41}$ = 4.5 eV$^2$.
This result is consistent with the no sterile neutrino hypothesis with a p-value of 8.0%.  The 90%  CL contour is shown in Fig. <ref>, left.

At low energy, sterile neutrino mixing within an extended neutral lepton mixing matrix enhances the standard oscillation probability proportionally to the matter column density traversed <cit.>.
The effect is approximately independent of $\Delta m^2$, as oscillation cycles are irresolvable within detector energy resolution.
IceCube has tested for this effect using a multi-flavor sample over an energy range of 10-60 GeV <cit.>, with the strongest effect expected at an energy of $\sim$20 GeV for upgoing muons.
The analysis yielded no evidence of anomalous oscillations, setting a limit on the extended PMNS matrix elements $|U_{\mu 4}|^2=\sin^2 \theta_{24}$ and $|U_{\tau 4}|^2=\sin^2 \theta_{34}\cos^2 \theta_{34}$, marked IC2017(NO) in Fig. <ref>, center.

Next generation sterile neutrino searches at IceCube

The sterile neutrino sensitivity at IceCube has yet to be exhausted, with near-term improvements expected from event samples already in hand.
At low energies, a sterile neutrino search using the full ten year dataset of $\geq$300,000 events with $E_{\nu}\leq 100$ GeV spanning all flavors is underway <cit.>, promising unique sensitivity in the  $|U_{\mu 4}|^2$, $|U_{\tau 4}|^2$ plane (Fig. <ref>, center).
At higher energies ($\geq$ 400 GeV), attention to date has been focused on searches for $\nu_\mu$ disappearance through non-zero $\theta_{24}$, motivated by the necessity of its finite value if sterile neutrino osculations were responsible for the short-baseline ($\nu_\mu\rightarrow\nu_e) / (\bar{\nu}_\mu\rightarrow\bar{\nu}_e $) anomalies ($\theta_{14}= \theta_{34}=0$ leads to the most conservative limits on $\theta_{24}$ <cit.>).
Efforts are now underway to incorporate the high energy cascade event sample into these analyses, which includes topologies associated with $\nu_e$ and $\nu_\tau$ CC and all flavor NC interactions.
This extension promises sensitivity to both $\nu_e$ and $\nu_\tau$ appearance signatures associated with non-zero $\theta_{14}$ and $\theta_{34}$, respectively.
Preliminary studies <cit.> of the sensitivity in this channel suggest that  $\nu_\tau$ appearance signatures are discoverable for values of  $\theta_{34}$ consistent with  world data and IceCube's existing $\theta_{24}$ limits; and that values of $\theta_{14}$ consistent with reactor <cit.> and gallium/BEST <cit.> anomalies may yield observable $\nu_e$ appearance signatures.
The expected sensitivity of the combined high energy $\nu_\mu$ disappearance and cascade appearance signatures is shown in Fig. <ref>, right.
Augmentations of the sterile neutrino searches using machine learning techniques, starting-event topologies, and improved reconstruction methods are also underway, expected to provide continuing improvements over the coming Snowmass period.

Left: Results from IceCube's high-energy muon-neutrino disappearance search, Middle: results from IceCube low- and high-energy analyses; Right: expected sensitivity with cascades.

Super-K Super-Kamiokande (Super-K, SK) has performed a search for light sterile neutrinos using approximately 4,000 live-days of atmospheric neutrino data <cit.>.
The SK analysis focuses on light sterile neutrinos with mass-squared differences greater than $\SI{0.1}\eV^2$.
In the energy range of the Super-K analysis, predominantly below 10 GeV, such a large mass-square difference implies that the oscillation effects due to a light sterile neutrinos are averaged out.
Thus, Super-K analysis is insensitive to the mass-square difference, but only to the mixing elements.
In the Super-K analysis they choose to constraint the mixing elements: $|U_{\tau 4}|$ and $|U_{\mu 4}|$.
In principle, Super-K has also sensitivity to $|U_{e 4}|$, but this is a subleading effect.
Additionally, the Super-K also is sensitive to the additional CP violating phases that appear in the presence of a sterile neutrino.