The effect of these CP violating phases is subleading, but does not always yield conservative results on the mixings and thus the results should be taken with this caveat <cit.>.

The Super-K analysis found no significant evidence of a light sterile neutrino and reported constraints on $|U_{\tau 4}|$ and $|U_{\mu 4}|$ for $\Delta m^2_{41} > \SI{0.1}\eV^2$; concretely they limit $|U_{\mu 4}|$ to less than 0.041 and $|U_{\tau 4}|$ to be less than 0.18 at 90% C.L.
These results are shown in Fig. <ref> and have been superseded by constraints from ANTARES and IceCube; compare to Fig. <ref>.

Constraint on $|U_{\tau 4}|$ and $|U_{\mu 4}|$ for $\Delta m^2_{41} > \SI{0.1}\eV^2$ obtained by Super-K. Figure from Ref. <cit.>.

The ANTARES neutrino telescope <cit.>, a sub-gigaton-scale neutrino telescope, had been designed and optimised for the exploration of the high-energy universe by using neutrinos as cosmic probes. However, its energy threshold of about 20 GeV was sufficiently low to be sensitive to the first atmospheric oscillation minimum. The majority of the neutrino events have been recorded with energies between 100 GeV and few TeV, an energy range rich in signatures of eV-scale sterile neutrinos.

Survival probability of vertically up-going $\nu_\mu$ at ANTARES as a function of neutrino energy for different values of mixing angles $\theta_{24}, \theta_{34}$ and $\delta_{24}$ with $\Delta m^2_{41}=0.5$ eV$^2$, $\Delta m^2_{31}=2.5\cdot 10^{-3}$ eV$^2$ and $\sin^22\theta_{23}=1$.

This is illustrated in <ref> which shows the $\nu_{\mu}$ survival probability for maximal mixing of $\theta_{23}$ and different combinations of the mixing parameters $\theta_{24},\theta_{34}$ and $\delta_{24}$ with

\begin{eqnarray}
U_{\mu 4}  &=& e^{-i\delta_{24}}\sin\theta_{24},\\
U_{\tau 4} &=& \sin\theta_{34}\cos\theta_{24}.
\label{Umu4}
\end{eqnarray}

The fast oscillations due to $\Delta m^2_{41} \gtrsim 0.5$ eV$^2$ are unobservable due to the limited energy resolution of the detector, making $\Delta m^2_{41}$ not measurable.

The ANTARES neutrino telescope was located in the Mediterranean Sea, 40 km off the coast of Toulon, France, at a mooring depth of about 2475 m. The detector was completed in 2008 and took data until February 2022. ANTARES was composed of 12 detection lines, instrumenting a water volume of about 15 Mtons. ANTARES data collected from 2007 to 2016 with an total detector lifetime of 2830 days have been used to constrain $U_{\mu 4}$ and $U_{\tau 4}$ <cit.>. A total of 7710 low-energetic atmospheric neutrino candidate events have been selected, largely dominated by muon-(anti)-neutrino charge-current events, identified thanks to a long-range up-going muon track. Particular attention was paid to a consistent handling of the complex phase $\delta_{24}$ in conjunction with the neutrino mass ordering. The limits from the ANTARES analysis are shown in <ref>. As expected from <ref> ANTARES is particularly sensitive if both $U_{\mu 4}$ and $U_{\tau 4}$ are non-zero and improves existing  limits from <cit.> substantially.

90% (left) and 99% (right) CL limits for the 3+1 neutrino model in the parameter plane of $|U_{\mu 4}|^2=\sin^2\theta_{24}$ and $|U_{\tau 4}|^2=\sin^2\theta_{34}\cos^2\theta_{24}$ obtained by ANTARES (black lines), and compared to the ones published by IceCube/DeepCore <cit.> (red) and Super-Kamiokande <cit.> (blue). The dashed lines are obtained for NH and $\delta_{24}=0^\circ$ while the solid lines are for an unconstrained $\delta_{24}$ (ANTARES) or for IH and $\delta_{24}=0^\circ$ (IceCube/Deepcore) respectively. The colored markers indicate the best-fit values for each experiment.
The 1D projections after profiling over the other variable are also shown for the result of this work.

After profiling over the other variable, the following limits on the two matrix elements can be derived:

\begin{eqnarray}
|U_{\mu 4}|^2  &<& 0.007~(0.13)~ \mathrm{at}~ 90\%~(99\%)~\mathrm{CL},\\
|U_{\tau 4}|^2 &<& 0.40~(0.68)~ \mathrm{at}~ 90\%~(99\%)~\mathrm{CL}.
\end{eqnarray}

§.§.§ Radioactive Source Experiments

BEST The Baksan Experiment on Sterile Transitions (BEST) <cit.> was proposed to probe the possibility of the short-baseline electron neutrino disappearance using the same process studied by GALLEX and SAGE, Eq. (<ref>). The disappearance of neutrinos was suggested as an explanation of the deficit observed events measured by the previous radioactive source experiments (see Sec. <ref>).

Although the main experimental idea is the same, the BEST configuration grants the measurements to be taken at two different baselines. In this case, a spherical vase (inner target) is located inside a cylindrical container (outer target), both filled with liquid gallium (the detector). The radioactive source ($^{51}\rm{Cr}$) is placed inside the sphere, as shown in Fig. <ref>.

BEST detector configuration <cit.>.

Here again, the emitted neutrinos interact with the detector through reaction Eq. (<ref>), and the produced $^{71}\rm{Ge}$ atoms are extracted and counted for each vessel separately <cit.>. The ratio of the measured rate of $^{71}\rm{Ge}$ production at each distance to the expected one considering the cross section and experimental efficiencies are presented in Tab. <ref> <cit.>.

Inner                      Outer
$0.79\pm 0.05$     $0.77\pm 0.05$

Ratio of predicted and observed $^{71}\rm{Ge}$ event rates as measured by BEST using $^{51}\rm{Cr}$ with its inner and outer targets.

Recall that the cross section of the reaction in Eq. (<ref>) is used in the calculation of these ratios, and BEST used the Bahcall results <cit.>, including conservative uncertainties <cit.>. The ratios are $4.2\sigma$ and $4.8\sigma$ less than unity, respectively, supporting the gallium anomaly observed by other experiments.

Allowed regions obtained in the $\sin^22\theta$–$\Delta{m}^2_{41}$ parameter space from the analysis of the BEST results only (left) and from the BEST results combined with results from GALLEX and SAGE (right) <cit.>. Note that here $\sin^22\theta  = 4\left|U_{e4}\right|^2\left(1 - \left|U_{e4}\right|^2\right)$.

BEST performed an analysis of these results in the framework of short-baseline electron neutrinos disappearance due to neutrino oscillation governed by the survival probability in Eq. (<ref>), using $\sin^22\theta  = 4\left|U_{e4}\right|^2\left(1 - \left|U_{e4}\right|^2\right)$. The study leads to the allowed regions shown in Fig. <ref>, corresponding to $1\sigma, \,2\sigma$, and $3\sigma$ confidence levels (left plot), with the best-fit at ($\sin^22\theta $,  $\Delta m^2) = (0.42, 3.3$ eV$^2$). When combined with the results from GALLEX and SAGE (Tab. <ref>), the allowed regions are the ones illustrated in the right plot of Fig. <ref>, where correlated cross section uncertainties were considered. In this case, the best fit is located at $(\sin^22\theta, \Delta m^2) = (0.34, 1.25$ eV$^2)$ <cit.>.

Some of the calculations of the cross section of the process Eq. (<ref>) appear to decrease the significance of the gallium anomaly related to the GALLEX and SAGE experiments (Sec. <ref>); nonetheless, the results are still consistent with a possible short-baseline electron neutrino disappearance produced by active to sterile neutrino oscillations. On the other hand, the measurements by BEST are in agreement with the other source experiments and indicate gallium anomaly with larger significance.

§.§.§ Direct Neutrino Mass Measurements

High precision $\beta$ spectroscopy enables searches for eV- to MeV-scale sterile neutrinos.
The electron flavor neutrino emitted in $\beta$ decay does not have a well-defined mass, but is rather an admixture of the neutrino mass eigenstates.
The existence of a hypothetical sterile neutrino(s) implies the electron flavor neutrino may also contain a small admixture of the new fourth neutrino mass eigenstate, $m_4$.
The $\beta$-decay spectrum, $R_\beta(E)$ will be altered:

\begin{equation}
R_\beta(E) = \cos^2\Theta_s \, R_{\beta}(E, m_{\nu}^2) + \sin^2\Theta_s \, R_\beta(E,  m_4^2) \,
\label{eq:sterile}
\end{equation}

composed of both the spectrum corresponding to the effective light neutrino mass $m_\nu^2$ and a spectrum associated with $m_4$.
The maximal energy of the spectrum is $E_0 - m_n$, with the largest neutrino mass resulting in the lowest endpoint energy.
The amplitudes of the two decay branches are given by $\cos^2\Theta_s$ and $\sin^2\Theta_s$, respectively, where $\Theta_s$ is the active-to-sterile mixing angle. The resulting signature of a sterile neutrino is thus a kink-like distortion of the measured spectrum at an energy of $E_0 - m_4$.
For eV-scale sterile neutrinos, this signature naturally appears in the region of interest for direct neutrino mass measurements studying the beta endpoint.

The Karlsruhe Tritium Neutrino experiment (KATRIN) <cit.> provides high-precision electron spectrum measurement of tritium $\beta$-decay, $^3$H $\rightarrow$ $^3$He$^+$ + e$^-$ + $\bar{\nu}_\mathrm{e}$
(endpoint $E_0$ = 18.57 keV, half-life $t_{1/2}$ = 12.32 yr). KATRIN is designed to improve the sensitivity on the effective neutrino mass, $m_\nu$, to 0.2-0.3 eV (90$\%$ CL). Based on its first four-week science run in spring 2019 <cit.>, KATRIN can constrain the mass and mixture of a sterile neutrino that would manifest itself as a distortion of the $\beta-$electron spectrum. The signature would be a kink-like feature, as shown in a simulation presented in Fig. <ref>. First light sterile neutrino results is based on data from KATRIN's first high-purity tritium campaign, which ran from April 10 to May 13, 2019, at an average source activity of $2.45\cdot10^{10}$ Bq <cit.>.

a) Electron spectrum of KATRIN data $R(\langle qU \rangle)$ over the interval $\lbrack E_0 -40$ eV, $E_0 +50$ eV$\rbrack$ from all 274 tritium scans and the three-neutrino mixing best-fit model $R_\mathrm{calc} (\langle qU \rangle)$ (line). The integral $\beta$-decay spectrum extends to $E_0$ on top of a flat background $R_\mathrm{bg}$. 1-$\sigma$ statistical errors are enlarged by a factor 50. b) Simulation of an arbitrary sterile neutrino imprint on the electron spectrum. The ratio of the simulated data without fluctuation, including a fourth neutrino of mass $m_4=10$ eV and mixing $|U_{e4}|^2=0.01$, to the three-neutrino mixing model is shown (red solid line). c) Integral measurement time distribution.

The integral $\beta-$electron spectrum is scanned repeatedly in the range of $\lbrack E_0 -$90 eV, $E_0 +$50 eV$\rbrack$ by applying non-equidistant HV settings to the spectrometer electrode system.
Each scan lasted 2 h. At each HV set point, the transmitted electrons are counted.
Fig. <ref> shows the measurement time distribution. A total of 274 stable scans are selected with an overall scanning time of 521.7 h. Then they gather the 117 most similar pixels and combine them into a single effective pixel.
The resulting stacked integral spectrum, R($\langle qU \rangle$), includes $2.03\cdot10^6$ events, with $1.48\cdot10^6$ $\beta-$electrons expected below $E_0$ and a flat background ensemble of $0.55\cdot10^6$ events in the whole scan interval.
The search has a high signal-to-background ratio, rapidly increasing from 1 at $\langle qU \rangle=E_0-$12 eV to $>$70 at $\langle qU \rangle=E_0-$40 eV.
The modeled experimental spectrum $R_{\mathrm{calc}}(\langle qU \rangle)$ is the convolution of the differential $\beta$-spectrum $R_\beta(E)$  with the response function $f(E - \langle qU \rangle)$, and an energy-independent background rate $R_\mathrm{bg}$:

\begin{equation}
R_\mathrm{calc} (\langle qU \rangle) = A_\mathrm{s}\,\cdot\,N_\mathrm{T}\int R_{\mathrm{\beta}}(E) \cdot  f(E - \langle qU \rangle)~ dE + R_\mathrm{bg},
\end{equation}

where $A_{\mathrm{s}}$ is the tritium signal amplitude.
$N_\mathrm{T}$ denotes the number of tritium atoms in the source multiplied with the accepted solid angle of the setup $\Delta \Omega/4\pi = (1- \cos{ \theta_\mathrm{max}})/2$, with $\theta_{\mathrm{max}} = 50.4^{\circ}$, and the detector efficiency (0.95).
The function $f(E - \langle qU \rangle)$ describes the transmission probability of an electron as a function of its surplus energy $E-\langle qU \rangle$.
The Collaboration extends the experimental modeling and statistical analysis to constrain both the sterile neutrino mass squared $m^2_{4}$ and its mixing amplitude $|U_{e4}|^2$, following the same strategy as for our $m_\nu$ analysis <cit.>.
In the 3+1 active-sterile neutrino model extension the electron spectrum, $R_\beta$, is replaced by $R_\beta(E,m_\nu,m_4) = (1~-~|U_{e4}|^2) R_\beta(E,m^2_\nu) + |U_{e4}|^2 R_\beta(E,m^2_{4})$, where $U$ is the extended $4\times4$ unitary matrix, $R_\beta(E,m^2_\nu)$ is the differential electron spectrum associated with decays the include active neutrinos in the final state, and $R_\beta(E,m^2_{4})$ describes the additional spectrum associated to decays involving a sterile neutrino of mass $m_{4}$.

95% confidence level exclusion curves in the ($\sin^2(2\theta_{ee}),\Delta{m}^2_{41}$) plane obtained from the analysis of KATRIN data with fixed $m_\nu=0$. The green contour delimits the 3+1 neutrino oscillations allowed at 95% confidence level by the reactor and gallium anomalies <cit.>. KATRIN data improve the exclusion of the high $\Delta{m}^2_{41}$ values with respect to DANSS, PROSPECT, and STEREO reactor measurements <cit.>. Mainz <cit.> and Troitsk <cit.> exclusion curves <cit.> are also displayed for comparison. An estimation of KATRIN's final sensitivity  is represented by the dotted line.

The observable integral spectrum $R_\mathrm{calc}$ is henceforth modeled with six free parameters: the four original parameters ($A_\mathrm{s}$, $E_0$, $R_\mathrm{bg}$, $m_\nu^2$) <cit.>, $m_4^2$ and $|U_{e4}|^2$.
This extended model $R_{\mathrm{calc}}(\langle qU \rangle)$ is then fitted to the experimental data $R(\langle qU \rangle)$.

In order to mitigate bias, the full analysis is first conducted on a Monte Carlo (MC) data set before turning to the actual data. The fit of $R(\langle qU \rangle)$ with $R_\mathrm{calc} (\langle qU \rangle)$ is performed by minimizing the standard $\chi^2$-estimator.
In a “shape-only” fit, both $E_0$ and $A_{\mathrm{s}}$ are left unconstrained. To propagate systematic uncertainties, a covariance matrix is computed after performing ${\cal O}(10^4)$ simulations.
The sum of all matrices encodes the total uncertainties of $R_\mathrm{calc} (\langle qU \rangle)$, including HV set-point-dependent correlations.
The $\chi^2$-estimator is then minimized to determine the best-fit parameters, and the shape of the $\chi^2$-function is used to infer the uncertainties.
The fit range $\lbrack E_0 -$40 eV, $E_0 +$50 eV$\rbrack$ is chosen such that statistical uncertainties on $ |U_{e4}|^2 $ dominate over systematic uncertainties in the whole range of $m_{4}^2$ considered <cit.>.
The experimental result agrees well with the sensitivity estimates and are displayed in Fig. <ref>.
They are showing no evidence for a sterile neutrino signal and are compared with short-baseline neutrino oscillation experiments measuring the electron (anti-)neutrino survival probability $P(\Delta{m}^2_{41},\sin^2(2\theta_{ee}))$.
An estimation of KATRIN's five-year sensitivity is also presented.

The KATRIN experiment continues accumulating tritium data and improving the operating conditions.  Using additional data with lower background, preliminary results significantly strengthen the sensitivity, particularly at higher $\Delta m^2_{41}$ <cit.>.

§.§ Dark Sectors in Scattering and in the Beam

The MiniBooNE low energy excess has traditionally been interpreted
as a potential hint of neutrino oscillations at short baseline driven
by a new sterile state. However, simple interpretations (e.g. $3+1$ models)
are in tension with other oscillation measurements (see the discussion in the sections above), in
particular searches for muon flavor disappearance driven by the same
$\Delta m^2$. Given this tension, it is natural to consider other
possible new physics explanations for the LEE. These new physics
explanations must have detector signatures and production mechanisms
which are consistent with the existing MiniBooNE measurements.

The measured angular distribution of the events in the LEE constrains
the allowed interaction channels in the detector. For concreteness,
let’s suppose that the LEE is due to some new particle $X$. Simple
signatures involving elastic scattering on electrons (e.g. $X e^-
\rightarrow X e^-$) or fully visible decays (e.g. $X \rightarrow e^+
e^-$) produce events which are very forward, inconsistent with the
measurement which features a broad angular distribution. While
there is a highly significant forward component to the excess,
especially apparent in the 2021 MiniBooNE LEE result, the excess
remains significant across a broad range of angles, and one cannot
simply presume that the forward excess is due to new physics while the
rest of the excess is due to background underestimation. Considering
semi-visible decays (e.g. $X \rightarrow X' + \gamma$) lessens the
disagreement somewhat, but such models are still strongly disfavored,
with either too many forward events or too many backward events
depending on the $X$ mass. In light of these constraints, the only
class of models which can adequately reproduce the LEE angular
distribution are those involving inelastic scattering, similar to the
neutrino-nucleus scattering in the standard sterile neutrino

§.§.§ Beam Dump Searches

Neutrino experiments as well as dedicated electron or proton beam dump experiments can search for dark sector particles by looking for the decays or scatterings of states produced at the target station. This technique has been extensively explored in the literature in the search for light dark matter and mediators <cit.>. Measurements at LSND and MiniBooNE provide competitive limits in several light dark sectors.

Of particular relevance are the results of the MiniBooNE-DM run. MiniBooNE ran in a beam dump configuration where the proton beam was
aimed directly at the downstream beam dump instead of onto the
neutrino production target. In this mode, charged meson decay in
flight was suppressed, reducing the “background” neutrino flux and
enhancing sensitivity to new physics. No excess of events was observed
in beam dump mode, implying that any new physics production modes that
have a simple scaling with the number of protons on target like
neutral meson decays or continuum processes (i.e. bremsstrahlung) are
ruled out because they should produce a signature in beam dump mode.
This leaves production from charged meson decay in flight as the only
viable production mode for new particles that can explain the LEE.

New physics explanations of the MiniBooNE LEE based on inelastic scattering signatures in the detector are favored due to the
compatibility with the measured angular distribution of the excess. The MiniBooNE-DM results imply that such processes should be initiated by either neutrinos or new particles produced in charged meson decays at the target <cit.>.

§.§.§ Neutrino-Electron Scattering Measurement

Models that can explain the LSND or MiniBooNE anomaly through the production of new particles that decay to electromagnetic showers can be constrained by measurements of neutrino-electron ($\nu-e$) scattering. Since single photons and collimated $e^+e^-$ pairs appear as single showers, they can be searched for by looking for in the photon-like sidebands of this measurement. In Ref. <cit.>, the authors propose a technique to constrain $e^+e^-$ explanations of LSND and MiniBooNE using data from MINER$\nu$A and CHARM-II. These experiments were located in the NuMI and CERN SPS beams, respectively, and therefore cover a much broader and higher-energy neutrino flux than LSND and MiniBooNE (see left panel of <ref>).

Left: the neutrino flux at MiniBooNE, MINER$\nu$A, and CHARM-II are shown as shaded regions. In lines the upscattering cross section for producing heavy neutrinos through coherent and incoherent $Z^\prime$ scattering is shown. Right: The single shower event spectra as a function of $dE/dx$ at MINER$\nu$A (top) and of $E\theta^2$ at CHARM-II (bottom). Dashed lines indicate analysis cuts. Figures taken from Ref. <cit.>.

In the dark photon model discussed in Sec. <ref>, the HNL decays to overlapping $e^+e^-$ predict new signals in the large-$dE/dx$ sideband of $\nu-e$ scattering, with moderate values of $E \theta^2$. The energy deposition $dE/dx$, defined as the deposited energy in the first centimeters of a electromagnetic shower, helps discriminate between electron-like showers from photon-like showers, which have typically twice the energy deposited per cm than the former.
The variable $E \theta^2$, defined as the shower energy times the square of the shower angle with respect to the neutrino beam, is also used to reduce backgrounds since in the boosted electron in $\nu-e$ scattering obeys $E \theta^2 < 2m_e$, while backgrounds, mostly from $\pi^0$ decay and $\nu_e$ CCQE, can have much broader angular distributions. For the dark sector signatures, the resulting overlapping $e^+e^-$ pairs can be quite forward, especially when the dark photon is light. Examples of the dark sector predictions at MINER$\nu$A as a function of $dE/dx$ and at CHARM-II as a function of $E\theta^2$ are shown in <ref>.

Since the backgrounds have large uncertainties and the search is not tailored to the dark sector signals, the sensitivity is not enough to rule out all explanations of the MiniBooNE excess. However, due to the higher-energies at MINER$\nu$A and CHARM-II, explanations with large HNL masses, often preferred due to the broader angular distributions at MiniBooNE, can be robustly excluded. In <ref>, the limits are shown in the mixing of the HNL with muon neutrinos, $|U_{\mu 4}|^2$, and their mass $m_4$. Other parameters like kinetic mixing  $\epsilon$ and the dark sector coupling $\alpha_D$, as well as the dark photon mass $m_{Z^\prime}$ are fixed. Two curves are shown for each experiment: a solid curve for the nominal choice of background uncertainty, and a dashed curve corresponding to the case where uncertainties were inflated by a factor of a few (see caption). For MINER$\nu$A LE (ME), this corresponds to $30\%$ ($40\%$) background normalization uncertainty in the nominal case and $100\%$ uncertainty in the conservative case. For CHARM-II these correspond to $3\%$ and $10\%$, respectively. The uncertainties at CHARM-II are constrained by the sideband at large $E \theta^2$.

Future measurements, including antineutrino-electron scattering at MINER$\nu$A can improve on the limits discussed above. The limits can also be recasted onto models with photon final states, like those discussed in <ref>. Finally, HNLs produced via scalar mediators are less constrained by this technique since the upscattering cross section decreases at higher energies.

The MINER$\nu$A and CHARM-II constraints in the parameter space of the HNL mixing with muons as a function of its mass, at $90$% C.L. Solid lines show the nominal constraints, while dashed ones show constraints with inflated systematic errors on the background. The arrows on the vertical lines indicate where more than a given percentage of the total excess events at MiniBooNE are predicted in the forward-most angular bin. Figure taken from Ref. <cit.>.

§.§.§ Searches for Long-Lived Particles

The possibility that meson decays produce long-lived HNLs in neutrino beams is strongly constrained by direct searches at T2K <cit.>, MicroBooNE <cit.>, and ArgoNeuT <cit.>, as well as by other past-generation experiments such as PS-191 <cit.> (recently re-evaluated in Ref. <cit.> and  <cit.>), CHARM <cit.>, CHARM-II <cit.>, and NuTeV <cit.>. The constraints are placed on minimal models with a single HNL $N$ that interacts with the weak bosons only, with interaction strength suppressed by small mixing angles $U_{\alpha N}$. At first sight, this model differs from the light sterile neutrino only for the mass of the additional neutrino, which is supposed to be in the $10 - 500$ MeV range.
However, this difference results in a radically different phenomenology.
Because these new particles are so much heavier than the standard neutrinos, no oscillation is possible: the HNL mixes with the other neutrinos but loses coherence immediately with the rest of the wave packet.
As a result, in the minimal scenario, in which no other particle or interaction is present aside from one or more HNLs, the HNL is produced in the beam through mixing and decays in the detector.
The available decay channels depend on the mass and include $e^+e^- \nu$, $e \mu \nu$, $e \pi$, $\mu^+ \mu^- \nu$, $\mu \pi$.
However, this minimal model is ruled out by a combination of Big Bang Nucleosynthesis, which provides a lower bound on the mixing parameters <cit.>, and experiments with neutrino beams, which lead to upper bounds, with no available parameter space in between these bounds (left plot in Fig. <ref>).

In particular, the T2K Near Detector ND280, a modular detector able to resolve details of interactions track and identify individual final state particles, sets the best limits on these models.
For this analysis, the three TPCs filled with gaseous argon provide a low-density decay volume for the HNLs, with zero background from neutrino interactions, as illustrated in the left plot in Fig. <ref>.
The limits have been derived through a search with little expected background and zero observed events in every analysis channel <cit.>, later extended and combined with BBN constraints <cit.>.

Left: Exclusion plot for the minimal scenario, showing the mixing with the muon flavour versus the HNL mass.
Right: Exclusion plot for the non-minimal scenario where HNLs possess a magnetic moment.
It shows the mixing with the muon flavour versus the magnetic dipole moment for an HNL mass of 250 MeV.
The region of interest to explain the MiniBooNE excess is shown in green <cit.>.
The dark grey region is taken from or extrapolated by <cit.>.
Both plots are taken from <cit.>.

This model becomes more interesting if new interactions are present.
These new interactions allow new decay modes, hence shortening the lifetime and relaxing BBN limits.

For example, if the HNL possesses a magnetic moment, it could decay electromagnetically into $\nu \gamma$.
For magnetic moments of the order of PeV, this model could explain the MiniBooNE anomaly <cit.>.
However, this explanation is constrained by short-baseline experiments, like MicroBooNE and ND280 (right plot in Fig. <ref>).
Thanks to the high density of the liquid argon, MicroBooNE identifies single photons through conversion to $e^+e^-$ pairs.
ND280 is instead sensitive to the branching ratio into off-shell photons, which results in a genuine $e^+e^-$ pair.
In this case, the rate is lower by a factor of  100 but benefits from a zero background search.

Upscattering Recasts

As discussed in Secs. <ref>, HNLs can be short-lived if they interact via additional forces via a dark photon or scalar, for instance.
In these models, the lifetime of the HNL can range from tens of meters to sub-mm, depending on the choice of parameters.
In this regime, standard neutrinos scattering on nuclei can produce a HNL, as illustrated in the right panel of Fig. <ref>.
This particle propagates inside the detector and subsequently decays into an $e^+e^-$ pair plus another neutrino. Searches for $e^+e^-$ final states at neutrino detectors can be used to constraint this possibility. In particular, the near detector of T2K, ND280, and the MicroBooNE detector have both been used to search for the decay of long-lived particles into $e^+e^-$

Left: Schematic representations of an HNL decaying in flight in one of the ND280 GArTPCs.
Right: The HNL is produced in the dense layers of lead through upscattering and subsequently decays in one of the GArTPCs.
Figures from <cit.>.

In the scenario of a dark photon heavier than the HNL, lifetimes range between a few cm to several meters.
In this case, a dark neutrino would be produced in the Pi0 detector (POD), a very dense detector composed of high-Z materials such as lead, which lies upstream with respect to the TPCs.
The production is particularly copious as the process is coherent, thus scaling with $Z^2$ and benefiting the high Z materials.
In the scenario of a dark photon lighter than the HNL, the decay proceeds through an on-shell dark photon, resulting in a larger decay width and shorter lifetimes, sub-millimetre.
In this case, the HNL is produced and decays at the same point in the detector.
The Fine-Grained Detectors (FGD) of ND280 come into play for this study.
This plastic scintillator is dense enough to make the rate for production through upscattering significant while allowing precise tracking and identification of the $e^+e^-$ pair.
Because of the larger density, backgrounds from the beam or photons that convert inside the FGD are present.
On the other hand, the invariant mass of the $e^+e^-$ pair is an excellent quantity to discriminate this background as it peaks around the dark photon mass because of the on-shell decay - making this analysis a peak search.
Explanations of MiniBooNE under this model are not entirely ruled out yet <cit.>, but the next generation of analyses, together with an upgrade to the detector <cit.> and a much larger dataset <cit.> that will be collected by T2K, are expected to probe the entire parameter space of interest.

§.§ Conventional Explanations and Other Searches

Partly motivated by the lack of a single interpretation that can simultaneously explain all observed experimental anomalies, the possibility that the anomalies represent a collection of conventional explanations (or that in combination with new physics) has also been discussed; see, e.g. Ref. <cit.>, in which the possibility that the MiniBooNE anomaly is a combination of underestimated backgrounds is explored. In this case it was shown that no combination of varying the  known backgrounds can completely explain the MiniBooNE LEE, it can potentially reduce the significance to $3\sigma$ if several less well-known backgrounds saturate their prescribed errors.

Key among the anomalies is the MiniBooNE low-energy excess, where it has been shown repeatedly that the excess shape at low energy is incompatible with 3+$N$ light sterile neutrino oscillations, raising the question whether this excess is due to another effect. Furthermore, as this energy region is dominated by non-$\nu_e$ backgrounds, conventional interpretations for the excess that do not rely on new physics or introduce electrons into the MiniBooNE detector have been raised. The majority of these rest on the known fact that MiniBooNE—as a Cherenkov detector—had no ability to distinguish a single photon from a single electron, and in particular target the various photon backgrounds irreducibly contributing to and dominating the MiniBooNE observed $\nu_e$ CCQE rate at low energy.

In the reactor sector, the experimental search for sterile neutrinos were primarily motivated by the anomalously low rates of fluxes measured by the past reactor neutrino experiments.
Overestimated  rates based on mismodeled reactor models, against which the fluxes were compared, could also explain the observed discrepancy without invoking new BSM particle physics.

§.§.§ Constraints on Single-Photon Production

The majority of conventional interpretations to the MiniBooNE excess propose that the excess is due to mis-identified photons from various neutrino interactions that produce photons and no charged lepton in the final state. In MiniBooNE, these mis-identified backgrounds are primarily contributed by $\pi^0$ production, due to reconstruction inefficiencies (where one photon from the $\pi^0$ decay may be missed) or more rare processes such as NC single-photon production through $\Delta(1232)$ resonance production and subsequent radiative decay. Many of these processes can and have been constrained in situ at MiniBooNE. This was the case for the rate and momentum-dependence of misidentified NC $\pi^0$ decays, which were constrained by a high-statistics measurement of events with two reconstructed electromagnetic shower Cerenkov rings, representing the two photons from NC $\pi^0$ decay. Similarly, the “dirt” background component in MiniBooNE, which was NC $\pi^0$-dominated, was directly constrained by a high-statistics measurements of events close to the detector boundary.  These processes were further studied and disfavored as the source of the MiniBooNE anomaly by studying both the radial distribution of the excess as well as the timing of the events relative to the known beam bunch timing <cit.>. On the other hand, photons from the radiative decay of the $\Delta(1232)$ baryon produced in neutrino NC interactions (NC $\Delta\rightarrow N\gamma$, where $N=p,n$) were an irreducible background in MiniBooNE. This background was not constrained in situ, but rather the rate was constrained indirectly through its correlation with the NC $\pi^0$ measured rate (which proceeds predominantly via $\Delta\rightarrow N\pi^0$ decay).

The MicroBooNE collaboration has recently performed a direct search for single-photon events coming from neutrino-induced NC production of the $\Delta$(1232) baryon resonance with subsequent $\Delta$ radiative decay <cit.>. As discussed in Sec. <ref>, $\Delta$ decay is expected to be the dominant source of single-photon events in neutrino-argon interactions below 1 GeV. Although $\Delta$ radiative decay is predicted in the SM, and measurements of photoproduction <cit.> and virtual Compton scattering <cit.> are well described by theory, this process has never been directly observed in neutrino scattering. In a fit to the radial distribution of the MiniBooNE data with statistical errors only, an enhancement of NC $\Delta\rightarrow N \gamma$ by a normalization factor of $3.18$ (quoted with no uncertainty) was found to provide the best fit for the observed excess, and in good spectral agreement with the observed excess across a number of reconstructed kinematic variables <cit.>. The MicroBooNE collaboration searched directly for an excess of this magnitude, interpreted as an overall enhancement to the theoretically-predicted NC $\Delta\rightarrow N\gamma$ rate.

MicroBooNE utilized the strength of LArTPC neutrino detectors to search for this excess of single-photons from the $\Delta$ decay, with and without a proton track present in the interaction (referred to as $1\gamma1p$ and $1\gamma0p$ final states). The presence of the reconstructed track in the $1\gamma1p$ selection allowed for the tagging of the neutrino interaction vertex, and subsequent reconstruction of the photon conversion distance in argon, which led to improved background rejection compared to the $1\gamma0p$ selection, where the lack of associated hadronic activity prevented neutrino vertex tagging.

The analysis used approximately half of the total collected MicroBooNE data to date ($6.9\times 10^{20}$ POT).  It selected and simultaneously fitted the $1\gamma1p$ and $1\gamma0p$ samples together with two additional, mutually exclusive but highly correlated samples with high NC $\pi^0$ purity. These high-statistics NC $\pi^0$ samples effectively constrained the rate and systematic uncertainty of NC $\pi^0$ production in argon, which was the dominant background to the $1\gamma1p$ and $1\gamma0p$ selections.

MicroBooNE observed no evidence for an excess of NC $\Delta$ radiative decay, as shown in Fig. <ref>. The measurement ruled out the normalization enhancement factor of 3.18 as an explanation to the MiniBooNE low-energy excess at 94.8% CL (1.9$\sigma$), in favor of the nominal prediction for NC $\Delta\rightarrow N\gamma$. Note that this was a model-dependent test of the MiniBooNE excess, and therefore does not necessarily apply universally to any photon-like interpretation. Those include, e.g. coherent single-photon production, which is expected to be a rarer process in MicroBooNE than NC $\Delta\rightarrow N\gamma$, or BSM processes that manifest as single-photon events, such as co-linear $e^+e^-$ from the decay of new particles. Those will be the target of dedicated follow-up MicroBooNE analyses, as well as model-independent single-photon searches.

Summary of MicroBooNE's NC $\Delta \rightarrow N \gamma$ single-photon searches for both with (left) and without (right) an associated proton. Figures from <cit.>. No evidence of an enhanced NC $\Delta\rightarrow N\gamma$ rate is observed.

§.§.§ Reactor Flux Models

As discussed in Sec. <ref>, modeling the reactor  flux is quite challenging.
State-of-the art models used to compare against measured IBD yields employed the conversion approach and relied on $\beta$-decay measurements performed at ILL in the 1980s.
The presence of mistakes in these beta-decay measurements, as well as incorrect assumptions present in the conversion process, could lead to mis-modeled reactor  fluxes.
Recent and $\beta$-decay measurements have played a crucial role in further exploring  flux prediction issues as a possible source of the RAA.

IBD yields of  and  from Daya Bay (blue), RENO (red), and a combined fit to Daya Bay and RENO. Horizontal and vertical black lines show the predicted IBD yields for  and  respectively based on the Huber-Mueller (HM) model. While  yields agree with the model within $1\sigma$, a $>3 \sigma$ discrepancy can be noticed in  yields. Figure from Ref. <cit.>.

Daya Bay <cit.> and RENO <cit.> experiments – leveraging their IBD yield measurements as a function of evolving fission fractions in the reactors – were able to perform simultaneous IBD yield measurements of  and .
While continuing to observe the same time-integrated IBD yield deficit that defines of the RAA,  these results also showed that while  yields are in good agreement with the models,  is at 0.925 $\pm$ 0.015 ($>3\sigma$) of the modeled yield as shown in Fig. <ref>.
While these results are in slight tension ($\sim 1 \sigma$) with the average  IBD yields from the pure  measurements done using HEU reactors <cit.>, they are in approximate agreement ($\sim 0.5 \sigma$)  with the modern  yield measured by the STEREO experiment <cit.>.
It is also worth noting that summation-predicted IBD yields are in good agreement with the Daya Bay's flux evolution measurements <cit.>.

To test the possibility of mis-modeled IBD yields in the conversion model, Kopeikin et al. <cit.> performed simultaneous $\beta$-decay spectrum measurements of  and  at the Kurchatov Institute, generating a conversion prediction, referred to as the KI Model that uses completely different primary inputs than those used by the Huber model.
While the measurement was designed to achieve an extremely high degree of correlation between  and  measurements, it exhibited lower statistical precision than the ILL measurements of the 1980s.
The KI model found a consistently higher ratio of  to  $\beta$-decay rates over the full energy range compared to the ILL measurements as shown in Fig. <ref>.
Potential bias in $\beta$ decay measurements have also been pointed out in a separate study <cit.>.

Ratio of the cumulative $\beta$ spectra for  and as measured at ILL <cit.> (blue) and at KI <cit.> (red). KI sees consistently lower ratios than ILL. Figure from  <cit.>.

These new conversion prediction results, together with the flux evolution measurements from Daya Bay and RENO and the summation-conversion flux evolution mis-matches mentioned in Sec. <ref> collectively suggest mis-modeled  IBD yields as a major contributor to the RAA.
This collective picture is well-illustrated in Fig. <ref>, which shows predictions and -based measurements of IBD yields for the dominant fission isotopes  and , given in terms of a ratio with respect to the Huber-Mueller conversion prediction.
Conversion predictions based on the ILL beta measurements (Huber-Mueller and the HKSS model discussed in Sec. <ref>) appear to deviate from all other measurements and predictions.
Huber-Mueller models show substantial deficits with respect to direct IBD yield measurements, whereas discrepancies between EF and KI models with respect to the data are not large enough to claim the existence of any anomaly whatsoever.

It is important to stress that while these flux-based indications are fairly clear in their suggestion of flux prediction issues in the Huber-Mueller model and cast doubt on a pure sterile neutrino interpretation of the RAA, this global flux picture still leaves substantial room for short-baseline sterile oscillation phenomena at reactors.
As previously mentioned in Sec. <ref>, it is likely that the error bars assigned to most, if not all, predictions in Fig. <ref> are under-estimated, leaving ample room for sterile oscillation amplitudes around the 20% level or lower.
Moreover, multiple studies show that hybrid models containing both incorrect flux predictions and sterile neutrinos also provide a good fit to global IBD yield datasets <cit.>; some of these scenarios produce best-fit oscillation-induced deficits well in excess of 10%.
In light of this degeneracy between incorrect flux predictions and constant oscillation-induced deficits, short-baseline reactor measurements capable of directly probing the deficit's $L/E$ character are the better bet for cleanly elucidating the the role played by short-baseline sterile oscillations in the reactor sector.

Best-fit points and 95 % CL, and 99 % CL contours of  and  IBD yields obtained using integrated rates (red), evolution measurements by Daya Bay and RENO (purple), and combined integrated rates and evolution measurements (gray). The axes correspond to $r_{i}$ ($\overline{R}_{\text{HKSS}}$ as described above) for  and  along x and y direction respectively. Predictions from HM, EF, and HKSS models are shown in cyan, orange, and blue contours respectively. Black dashed line represents $r_{235} = r_{239}$. It is clear from the plot that while the measured IBD yields agree with the predictions from the EF (ab inito) model for both isotopes, only  agrees with the HM and HKSS models–both of which are conversion models relying on ILL $\beta$ spectra measurements–pointing towards an issue with the conversion approach-based  predictions. Figure from Ref. <cit.>.

§.§.§ MicroBooNE CC Measurement

MicroBooNE has additionally performed a search that explicitly tests the nature of the MiniBooNE excess in a physics-model-agnostic way. Specifically, MicroBooNE has carried out three independent analyses to investigate whether the MiniBooNE observed low-energy excess can be described by an effective enhancement of $\nu_e$ CC scattering at low energy given by unfolding the MiniBooNE observed excess distribution. This unfolding predicts a factor of 5–7 enhancement of $\nu_e$ CC interactions below 500 MeV in true neutrino energy.

The data used in this search
correspond to approximately half the data set collected by MicroBooNE during its entire operational run time in the Fermilab BNB. The search was carried out by three separate analyses, each targeting different exclusive and inclusive $\nu_e$ CC final states, and using separate reconstruction paradigms and signal selections. The $\nu_e$ final states explicitly targeted by MicroBooNE include a pion-less final state topology with one electron and with 0 or $N>1$ protons as part of the visible hadronic final state; a CCQE-like final state topology with one electron and with only 1 proton as part of the visible hadronic final state; and a CC-inclusive final state topology with one electron and with any number of charged pions or protons as part of the visible hadronic final state. While probing different event topologies with distinct event reconstruction and selection methods, the three independent analyses still share several common aspects, including the  signal model, a Geant4-based simulation of the neutrino beam, the detector response model, and a tuned variation of the GENIE v3 event generator incorporating the most up to date knowledge of neutrino scattering in the $<1$ GeV energy range. The results from each of the three analyses are presented in detail in <cit.>, and are summarized in <cit.>. The final observed and predicted distributions from each search are reproduced in Fig. <ref>, with the observed data to Monte Carlo ratios shown in Fig. <ref>.

While statistics-limited, the three mutually-compatible analyses collectively reported no excess of low-energy $\nu_e$ candidates, and were found to be either consistent with or modestly lower than the predictions for all $\nu_e$ event classes, including inclusive and exclusive hadronic final-states, and across all energies. With the exception of the pion-less, zero-proton
which was the least sensitive to a simple model of the MiniBooNE low-energy excess, MicroBooNE rejected the hypothesis that an enhancement of $\nu_e$ CC interactions at low energy is fully responsible for
the MiniBooNE low-energy excess at $>97$% CL for both exclusive and inclusive event classes. Additionally, MicroBooNE disfavored generic $\nu_e$ interactions as the primary contributor to the excess, with a $1\sigma$ ($2\sigma$) upper limit on the inclusive $\nu_e$ CC contribution to the excess of 22% (51%).

MicroBooNE's four targeted $\nu_e$ CC spectra. Figures taken from <cit.>.

Summary of MicroBooNE $\nu_e$ CC measurements, from <cit.>.

It should further be noted that while these searches do not explicitly test the possibility of light sterile neutrino oscillations, recent phenomenological studies have further analyzed the observed data within that context, including $\nu_e$ appearance and/or $\nu_e$ disappearance. The results of those are summarized in Sec. <ref>.

§ ASTROPHYSICAL AND COSMOLOGICAL INDIRECT PROBES

It has been known for a long time that astrophysical observations can provide powerful constraints on BSM physics.  Stellar evolution arguments have long been invoked to constrain couplings between light BSM particles and SM particles <cit.>.  In particular, the success of a core-collapse supernova and the synthesis of heavy elements via the $r$-process therein maybe turned into an argument for or against sterile neutrino states (see, e.g., Ref. <cit.>).  In the past two decades, however, the strongest astrophysical statements on light sterile neutrinos have come from cosmological measurements, i.e., those observations that probe the universe on the largest length scales.  We discuss in this section the relevant theoretical arguments and the most recent observational constraints.

§.§ Cosmology

The standard hot big bang model predicts that the three generations of SM neutrinos are held in a state of thermodynamic equilibrium with other SM particles via the weak interaction in the first second post-big bang.  At temperatures below about 20 MeV, the dominant equilibrating interactions are

\begin{equation}
\begin{aligned}
\nu e & \to \nu e, \\
\nu \bar{\nu} & \to e^+ e^-, \\
\nu \nu & \to \nu \nu.
\end{aligned}
\end{equation}

As the universe expands and cools, these interactions become less frequent.  When the universe cools to a temperature of ${\cal O}(1)$ MeV, the interaction rate per neutrino $\Gamma$ drops below the Hubble expansion rate $H$.  From this point onwards, the neutrinos are said to be “decoupled” from the thermal plasma.  The typical energy of the neutrino ensemble at decoupling is $E \sim 3 T \gg m_\nu$, i.e., the ensemble is largely ultra-relativistic at decoupling. Because of this, the neutrinos retain to a high degree of accuracy their relativistic Fermi-Dirac phase space distribution parametrized by a temperature and possibly a nonzero chemical potential.

Shortly after neutrino decoupling, the $e^+e^-$ plasma becomes non-relativistic at $T \sim 0.5$ MeV.  Here, kinematics favor the annihilation of $e^+e^-$ pairs.  The energy released in this process “reheats” the photon population.  The neutrino population, however, does not feel this reheating, because they have already decoupled: in other words, the weak interaction processes are no longer efficient at transporting the energy released from the annihilation to the neutrino sector.  The net effect is that the neutrinos emerge from the annihilation event a little colder than the photons. Assuming (i) ideal gases, (ii) instantaneous neutrino decoupling, and (iii) ultra-relativistic electrons/positrons at neutrino decoupling, one can show using entropy conservation arguments

\begin{equation}
\label{eq:nutemp}
T_\nu = \left(\frac{4}{11} \right)^{1/3} T_\gamma,
\end{equation}

where $T_\nu$ and $T_\gamma$ are the neutrino and the photon temperatures, respectively, after $e^+e^-$ annihilation.

Eq. <ref> is generally taken to define the standard neutrino temperature
after $e^+e^-$ annihilation.  In reality, however, dropping any one of the three aforementioned assumptions can result in percent-level corrections to the relation.  These corrections are commonly absorbed into the definition of the effective number of neutrinos $N_{\rm eff}$ via

\begin{equation}
\label{eq:neff}
\rho_\nu = N_{\rm eff} \, \frac{7 \pi^2}{120} \, T_\nu^4 = N_{\rm eff} \, \frac{7}{8} \left(\frac{4}{11} \right)^{4/3} \rho_\gamma.
\end{equation}

Here, $\rho_\nu$ is the actual total neutrino energy density after $e^+e^-$ annihilation (but still deep in the radiation domination era when the neutrinos are ultra-relativistic), the quantity $(7 \pi^2/120)T_\nu^4$ denotes the energy density in one family of thermally-distributed relativistic neutrinos with a temperature $T_\nu$ defined in Eq. <ref>, and $\rho_\gamma$ is the actual energy density in the photon population, i.e., what eventually becomes the cosmic microwave background (CMB) radiation.  Precision calculations of the SM prediction of $N_{\rm eff}$ find $N_{\rm eff}^{\rm SM} = 3.0440 \pm 0.0002$ <cit.>, including the effects summarized in Tab. <ref>.[Note that because the parameter $N_{\rm eff}$ is defined in relation to the neutrino energy density in an epoch when the neutrinos are ultra-relativistic, there is no ambiguity in the definition (Eq. <ref>) even if the neutrinos
should become non-relativistic at a later time because of their nonzero masses, provided the masses do not exceed the eV scale.]

Standard-model corrections to $N_{\rm eff}^{\rm SM}$     Leading-digit contribution

$m_e/T_d$ correction    $+0.04$

$\mathcal{O}(e^2)$ FTQED correction to the QED EoS    $+0.01$

Non-instantaneous decoupling+spectral distortion     $-0.005$

$\mathcal{O}(e^3)$ FTQED correction to the QED EoS    $-0.001$

Neutrino flavor oscillations     $+0.0005$

Type (a) FTQED corrections to the weak rates     $\lesssim 10^{-4}$

Leading-digit contributions from various SM corrections, in order of importance, that make up the final $N_{\rm eff}^{\rm SM}-3$ (adapted from Ref. <cit.>). The largest, $m_e/T_d$ correction results from dropping the assumption of an ultra-relativistic electron-positron population; finite-temperature quantum electrodynamics corrections (FTQED) to the QED equation of state (EoS) enter at ${\cal O}(e^2)$ and ${\cal O}(e^3)$, where $e$ is the elementary electric charge; the non-instantaneous decoupling+spectral distortion correction is defined relative to an estimate of $N_{\rm eff}^{\rm SM}$ in the limit of instantaneous decoupling assuming $T_d = 1.3453$ MeV <cit.>; and Type (a) FTQED corrections to the weak rates refer to neutrino-electron scattering rates corrected with thermal masses.

The effective number of neutrinos $N_{\rm eff}$ is of interest to cosmology primarily because energy density in relativistic particles affect directly the Hubble expansion expansion rate during the radiation domination era. In the epoch after $e^+e^-$ annihilation, the expansion rate is given by

\begin{equation}
\label{eq:hubble}
H^2(t)  = \frac{8 \pi G}{3} \left(\rho_\gamma + \rho_\nu \right),
\end{equation}

where $G$ is the gravitational constant.  With the photon energy density $\rho_\gamma$ having been measured to better than $0.1$% accuracy by the FIRAS instrument on board COBE <cit.>,
constraints on $H(t)$ in the early universe can be interpreted as bounds on the ratio $\rho_\nu/\rho_\gamma$ and hence $N_{\rm eff}$.

From a particle physics standpoint, a thermal
population of light sterile neutrinos is one possible cause of a $N_{\rm eff}$ that differs from the standard $N_{\rm eff}^{\rm SM}$ value.   However, it is
important to emphasize that, as far as the Hubble expansion rate $H(t)$ is concerned, any thermal background or non-thermal population (e.g., from decays) of non-photon light particles such as axions,  majorons, or even gravitons will contribute to $N_{\rm eff}$. Such scenarios have been considered by many authors, including Refs. <cit.>.
Likewise, any process that alters the thermal abundance
of neutrinos (e.g., a low reheating temperature) or affects directly the expansion rate itself (e.g., a time-dependent gravitational constant $G$) can mimic a non-standard $N_{\rm eff}$ value. Yet another way to change $N_{\rm eff}$ is to tinker with the photon energy density itself (via, e.g., interaction with millicharged particles <cit.> or late kinetic decoupling of the dark matter <cit.>),
while preserving the neutrino energy density.

In the case of a non-standard $N_{\rm eff}$ due to a BSM light particle, the mass of the new particle can also impact on the evolution of cosmological density perturbations via its role as a “hot dark matter”; the mathematical description of this effect goes beyond the Hubble expansion rate (Eq. <ref>).

In the following, we describe first how a thermal population of light sterile neutrinos can arise in the early universe through a combination of neutrino flavor oscillations and scattering with other SM particles.
We then discuss the various signatures of light sterile neutrinos in cosmological observables such as the light elemental abundances from big bang nucleosynthesis (BBN), the CMB anisotropies and the large-scale structure (LSS) distribution, constraints from current observations, as well as various proposals on how to get around them.

§.§.§ Light Sterile Neutrino Thermalization

If a sterile neutrino state mixes sufficiently strongly with any of the active  neutrino states, a thermal population of light sterile neutrinos that adds to $N_{\rm eff}$ can be produced prior to neutrino decoupling via a combination of active-sterile neutrino oscillations and collisions (i.e., elastic and inelastic scattering) with the primordial plasma of SM particles.  Roughly speaking, as the universe cools, an initial population comprising only active neutrinos can begin to oscillate into sterile neutrinos once the oscillation frequency, given by $\Delta m^2/(2 E)$, becomes larger than the Hubble expansion rate $H(t)$.  The role of collisions is then to force a neutrino into a flavor eigenstate and hence “measure” the flavour content of the ensemble.  Since the probability of measuring a sterile flavor is nonzero, signifying that an active neutrino has turned into a sterile state, collisions also play the role of refilling any gap in the active neutrino distribution vacated by the oscillation process.
This effect and the region of parameter
space leading to thermalization of the sterile neutrino was found by the early works of, e.g., <cit.>
If sterile neutrinos do become thermalized, then we expect them to have the same temperature as the
active neutrinos.

Nowadays, light sterile neutrino thermalization in the early universe can be computed precisely, using a generalized Boltzmann formalism developed in <cit.> which tracks the flavor evolution of a neutrino ensemble under the influence of neutrino oscillations and scattering.  Schematically, the generalized Boltzmann equation for the one-particle reduced density matrix of the neutrino ensemble, $\varrho(t,p)$, in the Friedmann-Lemaître-Robertson-Walker (FLRW) universe of standard hot big bang cosmology is given by

\begin{equation}
\partial_t \varrho - p H \partial_p \varrho = -{\rm i} [\mathbb{H},\varrho] + {\cal I}[\varrho],
\end{equation}

where $\partial_t$ and $\partial_p$ are partial derivatives with respect to the cosmic time $t$ and physical momentum $p$ respectively, $H$ is the Hubble expansion rate, $[\mathbb{H},\varrho] \equiv \mathbb{H} \varrho-\varrho \mathbb{H}$ denotes a commutator between the flavor oscillations Hamiltonian $\mathbb{H}$ and $\varrho$, and the collision integrals ${\cal I}[\varrho]$ encapsulate all non-unitary (scattering) effects on $\varrho$.  In the fully CP symmetric case, one set of density matrix $\varrho(t,p)$ suffices to describe the evolution of the whole neutrino ensemble including antineutrinos. If however the system is CP asymmetric, we will need to introduce a
separate one-particle reduced density matrix $\bar{\varrho}(t,p)$ for the antineutrino ensemble.  Here, we follow the convention of Ref. <cit.>, and define the density matrices using the “transposed notation, e.g., for a 3 active+1 sterile system one would have

\begin{equation}
\varrho(t,p) = \left(\begin{array}{cccc}
\varrho_{ee} & \varrho_{e \mu} & \varrho_{e \tau} & \varrho_{es} \\
\varrho_{\mu e} & \varrho_{\mu \mu} & \varrho_{\mu \tau} & \varrho_{\mu s}\\
\varrho_{\tau e} & \varrho_{\tau \mu} & \varrho_{\tau \tau} & \varrho_{\tau s} \\
\varrho_{s e} & \varrho_{s \mu} & \varrho_{s \tau} & \varrho_{s s} \\
\end{array}
\right), \qquad
\bar{\varrho}(t,p) = \left(\begin{array}{cccc}
\bar{\varrho}_{ee} & \bar{\varrho}_{\mu e} & \bar{\varrho}_{\tau e} & \bar{\varrho}_{se} \\
\bar{\varrho}_{e\mu} & \bar{\varrho}_{\mu \mu} & \bar{\varrho}_{\tau\mu} & \bar{\varrho}_{s \mu}\\
\bar{\varrho}_{e \tau} & \bar{\varrho}_{\mu \tau} & \bar{\varrho}_{\tau \tau} & \bar{\varrho}_{s \tau} \\
\bar{\varrho}_{e s} & \bar{\varrho}_{\mu s} & \bar{\varrho}_{\tau s} & \bar{\varrho}_{s s} \\
\end{array}
\right) .
\end{equation}

This convention enables the equations of motion to be expressed in a more compact form.

Working in the flavor basis, the oscillations Hamiltonian under this convention is

\begin{equation}
\begin{aligned}
\mathbb{H} (p,t,T(t)) =\, &   \pm \frac{U\mathbb{M}U^\dagger}{2p} + \sqrt{2} G_F \left[\mathbb{N}_\ell(T(t)) - \mathbb{N}_{\bar{\ell}}(T(t))+ \mathbb{N}_\nu(t)-\mathbb{N}_{\bar\nu}(t) \right] \\
& \mp 2\sqrt{2}G_F p
\left[
\frac{{\mathbb{E}}_\ell(T(t))+{\mathbb{E}}_{\bar\ell}(T(t))+{\mathbb{P}}_\ell(T(t))+{\mathbb{P}}_{\bar\ell}(T(t))}{m_W^2}
\frac{4}{3} \frac{{\mathbb{E}}_\nu(t)+{\mathbb{E}}_{\bar\nu}(t)}{m_Z^2}
\right],
\label{eq:hamiltonian}
\end{aligned}
\end{equation}

which contains a vacuum and an in-medium part.
The vacuum part consists of the neutrino squared-mass difference matrix $\mathbb{M}$ and the vacuum mixing matrix $U$, with the “$+$” sign pertaining to neutrinos and the “$-$” sign to the antineutrinos.
The in-medium part proportional to the Fermi constant $G_F$, i.e., the “matter potential”, contains a CP symmetric and a CP asymmetric correction to the neutrino dispersion relation.  The CP asymmetric part is similar to the usual matter potential found in, e.g., the Sun, with the matrix $\mathbb{N}_\ell-\mathbb{N}_{\bar{\ell}}$ denoting the asymmetry in the number density of charged leptons; in a 3+1 system $\mathbb{N}_\ell$ would take the form

\begin{equation}
\mathbb{N}_\ell \equiv \frac{1}{2 \pi^2} \int {\rm d} p \, p^2 \, {\rm diag}\Big(f_e,f_\mu, f_\tau,0  \Big) = {\rm diag}(n_e,n_\mu,n_\tau,0),
\end{equation}

where $f_\ell(p,T)$ is the (equilibrium) occupation number of the charged lepton $\ell$.  A similar expression exists also for $\mathbb{N}_{\bar{\ell}}$.  Typically, the charged-lepton asymmetries are of order $10^{-10}$, so this $\mathbb{N}_\ell-\mathbb{N}_{\bar{\ell}}$ term is not strictly a necessary ingredient.
Observe however the additional term proportional to ${\mathbb{N}}_\nu-{\mathbb{N}}_{\bar\nu}$, where

\begin{equation}
\mathbb{N}_\nu\equiv\frac{1}{2 \pi^2} \int \mathrm{d}p\,p^3 \, S_a \varrho S_a,
\end{equation}

and similarly for $\mathbb{N}_{\bar\nu}$.
This CP asymmetric term describes neutrino self-interaction in the presence of a large ($>10^{-5}$) neutrino asymmetry.  Current cosmological data constrain neutrino asymmetries only to  ${\cal O}(10^{-2})$ <cit.>, so the presence of a sizeable $\mathbb{N}_\nu-\mathbb{N}_{\bar{\nu}}$   term is an interesting possibility.  In standard calculations of light sterile neutrino thermalization, however, this asymmetry is set to zero.
The quantity $S_a$ is a diagonal matrix that projects out only the active neutrino states, since sterile states are by definition interaction-less.  In a 3+1 system it takes the form $S_a = {\rm diag}(1,1,1,0)$.

In the CP symmetric portion of the in-medium terms,
$m_W$ and $m_Z$ denote respectively the $W$ and $Z$ boson mass, while the terms $\mathbb{E}_\ell$ and $\mathbb{P}_\ell$ are momentum-integrals of some combinations of the charged-lepton energy $E_\ell(p)\equiv (p^2+ m_\ell^2)^{1/2}$ and occupation number $f_\ell(p,T)$, and $\mathbb{E}_{\bar\ell}$, $\mathbb{P}_{\bar\ell}$ are their antiparticle counterparts.  For a 3+1 system, we have

\begin{eqnarray}
\mathbb{E}_\ell & \equiv& \frac{1}{2 \pi^2} \int {\rm d} p \, p^2 \, {\rm diag}\Big(E_e f_e,E_\mu f_\mu, E_\tau f_\tau,0  \Big)
{\rm diag}(\rho_e,\rho_\mu,\rho_\tau,0),\label{eq:eell}\\
\mathbb{P}_\ell &\equiv& \frac{1}{6 \pi^2} \int {\rm d} p \, p^2 \, {\rm diag}\Bigg(\frac{p^2}{E_e} f_e,\frac{p^2}{E_\mu} f_\mu,\frac{p^2}{E_\tau} f_\tau,0\Bigg)
={\rm diag}(P_e,P_\mu,P_\tau, 0), \label{eq:pell}
\end{eqnarray}

where $\rho_\ell$ and $P_\ell$ are, respectively, the energy density and pressure of the charged lepton $\ell$.  Lastly, $\mathbb{E}_\nu$ is the equivalent of Eq. <ref> for an ultra-relativistic neutrino gas,

\begin{equation}
\label{eq:nurhosum}
\mathbb{E}_\nu\equiv\frac{1}{2 \pi^2} \int \mathrm{d}p\,p^3 \, S_a\varrho S_a,
\end{equation}

where $S_a$ is again a diagonal matrix that projects out only the active neutrino states.

Note that the $\mathbb{E}_\ell +\mathbb{P}_\ell$ term in Eq. <ref> differs from its usual presentation found in, e.g., equation (2.2) of <cit.>, which has $\mathbb{E}_\ell +\mathbb{P}_\ell$ replaced with $(4/3)\, \mathbb{E}_\ell$.
First reported in <cit.>, the former is in fact the more general result, while $(4/3)\, \mathbb{E}_\ell$ applies strictly only when the charged leptons are ultra-relativistic.

The collision integral ${\cal I}[\varrho]$ incorporates in principle all weak scattering processes wherein at least one neutrino appears in either the initial or final state.
All published calculations to date, however, account only for $2 \to 2$ processes involving
(i) two neutrinos and two charged leptons anyway distributed in the initial and final states, and
(ii) neutrino-neutrino scattering.
Then, schematically, ${\cal I}[\varrho]$ comprises 9D momentum-integrals that, at tree level, can be systematically reduced to 2D integrals of the form

\begin{equation}
{\cal I}[\varrho(p,t)]
\propto G_F^2 \int {\rm d}p_2\, {\rm d}p_3\, \Pi(p,p_2,p_3;t) \, F(\varrho), \label{eq:collision1}
\end{equation}

where $\Pi$ is a scalar function representing the scattering kernel, and $F$ is a phase space matrix including quantum statistics.

The vast majority of existing works on light sterile neutrino thermalization solve the above set of equations of motion for a 1 active+1 sterile system. Because of the complexity of the collision integrals ${\cal I}[\varrho]$, a variety of approximations have been introduced to simplify the integrals
and hence speed up the computation <cit.>, although of course it is also possible to solve the collision integrals exactly if percent-level precision is required <cit.>.
The publicly available code LASAGNA <cit.> solves a 1+1 system and also allows for the possibility of a CP asymmetry.  On the other hand,
multi-flavor (i.e., 3+1) effects have been  considered in Ref. <cit.> using the momentum-averaged approximation and in Ref. <cit.> using a multi-momentum approach.  Most recently, a publicly available multi-flavor code FortEPiaNO <cit.> has also become available, which is capable of solving a system with up to 3 active+3 sterile species, although it does not provide the CP asymmetric case.

The change in the effective number of neutrinos, $\Delta N_{\rm eff}  \equiv N_{\rm eff} - N_{\rm eff}^{\rm SM}$, at a temperature $T = 0.1$ MeV from active-sterile neutrino oscillations in a 1+1 system, where $\delta m^2$ is the squared mass difference and $\theta$ is the effective mixing angle between the active and the sterile state.  Figure taken from Ref. <cit.>.

The final effective number of neutrinos from active-sterile neutrino oscillations in a 3+1 scheme, under the constraint $\sum \log_{10} |U_{i4}|^2 = -13$.  Figure taken from Ref. <cit.>.

Irrespective of how exactly the system is solved (3+1, 1+1, approximate or exact collision integrals, etc.), however, the general conclusion is that the active-sterile neutrino squared mass difference  and mixing required to explain the short-baseline anomalies will inevitably lead to a fully thermalized light sterile neutrino population with a temperature similar to that of the active neutrinos.  In other words, in terms of the $N_{\rm eff}$ parameter, we expect the canonical light sterile neutrino solution to the short-baseline anomalies to lead to $N_{\rm eff} \sim 4$. Fig. <ref> and  <ref> show the prediction for $N_{\rm eff}$ as a function of the mixing parameters in a 1+1 and 3+1 scheme respectively.  We discuss in the next section the observable consequences of this light sterile neutrino population.

§.§.§ Observable Consequences

Three standard cosmological probes are sensitive to the presence of a light sterile neutrino population primarily through its contribution to increasing $N_{\rm eff}$ and hence the Hubble expansion rate, and/or through its nonzero mass and hence its role as a “hot dark matter”.  We describe these probes below.

Big Bang Nucleosynthesis

The discussion of light element formation from protons and neutrons via the process of big bang nucleosynthesis (BBN) (see Ref. <cit.> for a recent review) usually begins at a temperature around $T \sim 0.7$ MeV, when the weak processes

\begin{equation}
\begin{aligned}
\label{eq:betaprocess}
\nu_e + n & \leftrightarrow e^- + p,\\
\bar{\nu}_e + p &\leftrightarrow e^+ + n
\end{aligned}
\end{equation}

become inefficient compared with the Hubble expansion rate and equilibrium can no longer be maintained.  When these process go out of equilibrium, the ratio of neutron-to-proton number density also freezes out, to a value given by

\begin{equation}
\label{eq:np}
\left. \frac{n_n}{n_p} \right|_{\rm fr}   = \exp \left(- \frac{m_n-m_p}{T_{\rm fr}} \right),
\end{equation}

where $T_{\rm fr}$ is the freeze-out temperature, and $m_{n,p}$ are the neutron and proton masses.  For $T_{\rm fr} \simeq 0.7$ MeV, this ratio evaluates approximately to $1/6$. Free neutron decay over a lifetime of about 880 s, however, will reduce it to a smaller number by the end of BBN.

In standard BBN, the formation of elements commences at a temperature controlled by the baryon-to-photon ratio $\eta$, when the energy in the photon bath per baryon has become sufficiently low such that newly formed nuclei are no longer immediately broken apart.  For $\eta \sim 6 \times 10^{-10}$, this temperature is around $T \sim 0.1$ MeV.
The first element to be formed is Deuterium, followed by the production of heavier nuclei. Of particular note is Helium-4, whose mass fraction is defined as

\begin{equation}
\label{eq:yhe}
Y_p \equiv \frac{4 n_{\rm He 4}}{n_n+n_p},
\end{equation}

where $n_{\rm He4}$ is the number density of $^4$He.  Because $^4$He has the largest binding energy amongst the light elements, the bulk of all initially available neutrons will eventually end up bound in Helium-4 nuclei,
i.e., we expect $n_{\rm He 4} \simeq n_n/2$.  Then, to estimate $Y_p$ from Eq. <ref> we simply need to note that the neutron-to-proton ratio typically drops to about $1/7$ at the end of BBN via neutron decay.  From this we find $Y_p \simeq 0.25$.

Besides Helium-4 mass fraction, small amounts of Deuterium and $^3$He (D/H $\sim ^3$He/H $\sim O(10^{-5})$),
as well as traces of $^6$Li and $^7$Li, are expected to remain.  Unlike for $^4$He, however, there are no simple ways to estimate their abundances, and we must rely on solving a set of Boltzmann equations to track their number densities.  Several publicly available codes can perform this task, including AlterBBN <cit.>, PArthENoPE <cit.>, and PRIMAT <cit.>.  In standard BBN, barring experimental uncertainties in the nuclear reaction rates and the free neutron lifetime, the baryon-to-photon ratio $\eta$ alone enters these Boltzmann equations and is hence the sole free parameter in the determination of the elemental abundances.

A non-standard neutrino sector can alter this picture in two different ways.  Firstly, as can be seen in Eq. <ref>, electron neutrinos participate directly in the CC weak interactions that determine the neutron-to-proton ratio.  If because of non-standard physics these neutrinos should end up at $T \lesssim 1$ MeV with an energy spectrum that departs strongly from an equilibrium relativistic Fermi-Dirac distribution with zero chemical potential, then the equilibrium of the processes (Eq. <ref>) could shift to a different point and in so doing alter the neutron-to-proton ratio.  A particularly well studied example in this regard is the case of a nonzero electron neutrino chemical potential $\mu_e$, which shifts the neutron-to-proton ratio at weak freeze-out in a manner well described by

\begin{equation}
\left. \frac{n_n}{n_p} \right|_{\rm fr}   = \exp \left(- \frac{m_n-m_p}{T_{\rm fr}}-\frac{\mu_e}{T_{\rm fr}} \right).
\end{equation}

The effects of more general distortions to the $\nu_e$ and/or $\bar{\nu}_e$ energy spectra on $n_n/n_p$ need to be computed numerically  using a Boltzmann code.

Deuterium, Helium-4 and Lithium-7 abundances, as computed using PArthENoPE <cit.>, as a function of cosmological parameters: $\omega_b$ is the baryon density, $N_s$ is our $\Delta N_{\rm eff}  \equiv N_{\rm eff} - N_{\rm eff}^{\rm SM}$, and $\xi \equiv \mu_e/T$ is the electron neutrino degeneracy parameter.
Figure taken from Ref. <cit.>.

Secondly, neutrinos of all flavors influence the expansion rate of the universe prior to and during BBN through their energy densities via Eq. <ref>.  Therefore, if because of new physics the total neutrino energy should be larger than the standard $N_{\rm eff}^{\rm SM}$, the freeze-out of the processes (Eq. <ref>) would occur at a higher temperature.  This in turn pushes up the neutron-to-proton ratio via Eq. <ref> and hence the Helium-4 mass fraction as well via Eq. <ref>.
Fig. <ref> shows the Deuterium, Helium-4 and Lithium-7 abundances computed using PArthENoPE <cit.>, as a function of the baryon density $\omega_b$, the excess number of relativistic degrees of freedom $N_s \equiv \Delta N_{\rm eff}  \equiv N_{\rm eff} - N_{\rm eff}^{\rm SM}$, and the electron neutrino degeneracy parameter  $\xi \equiv \mu_e/T$.

For the specific problem of a light sterile neutrino with  mixing parameters compatible with hints from terrestrial experiments, only the second effect is relevant. This is because the relatively large mass-squared difference and mixing between the active and the sterile neutrino states essentially guarantee full thermalization of the sterile species prior to the decoupling of the active neutrinos (see Fig. <ref> and  <ref>).   In other words, $\Delta N_{\rm eff} \simeq 1$ and the phase space distribution of the sterile states follows closely the equilibrium relativistic Fermi-Dirac distribution of the active neutrinos. This also implies equipartition amongst the four neutrino flavors, such that any further active-sterile flavor oscillations after neutrino decoupling will not cause the $\nu_e$ phase space distribution to deviate from a thermal distribution.  Thus, to constrain such light sterile neutrino scenarios, we only need to extend standard BBN with one extra free parameter, namely, $N_{\rm eff}$.  We emphasize, however, that more general cases of active-sterile neutrino mixing would require that we determine the $\nu_e$ energy spectrum as well, in order to determine the full effect of light sterile states on the light elemental abundances.

Top: 68.27% and 95.45% contours in the $(\Omega_b h^2, N_{\rm eff})$-plane obtained from various data combinations.  Bottom: 1D marginalized posterior for $N_{\rm eff}$.
Figure taken from Ref. <cit.>.

Current BBN constraints Current astrophysical observations put the Helium-4 mass fraction at
$Y_p = 0.2449 \pm 0.0040$ <cit.>, the Deuterium abundance at ${\rm D}/{\rm H} = (2.527 \pm 0.030) \times 10^{-5}$ <cit.>, and the Lithium-7 abundance at $^7{\rm Li}/{\rm H} = (1.58 \pm 0.3) \times 10^{-10}$ <cit.>, while  the Helium-3 abundance is constrained to $^3{\rm He}/{\rm H}  < (1.1 \pm 0.2) \times 10^{-5}$ <cit.>.  Of these, only the D/H and $Y_p$ measurements have sufficient precision to probe $N_{\rm eff}$ during the BBN epoch.  How D/H and $Y_p$ probe $N_{\rm eff}$ can be seen in Fig. <ref>, which shows in the top panel the 68.27% and 95.45% contours in the $(\Omega_b h^2, N_{\rm eff})$-plane, where $\Omega_b h^2 \equiv \omega_b$ is the baryon density, obtained from
various data combinations, and in the bottom panel the corresponding 1D marginalized posterior for $N_{\rm eff}$.

As can be seen, the Helium-4 mass fraction alone is already quite sensitive to $N_{\rm eff}$, although this measurement is not particularly useful for pinning down the baryon-to-photon ratio $\eta$ and hence the baryon density $\Omega_b h^2$.  This is because $\eta$ has no strong influence on the rate at which $^4$He is formed, other than setting the initial time of BBN.
In contrast, the Deuterium abundance is strongly sensitive to both $N_{\rm eff}$ and $\eta$, which do directly affect the formation rate. However, because these two measurements have opposite degeneracy directions on the $(\Omega_b h^2, N_{\rm eff})$-plane, in combination they provide a good measurement of both $\Omega_b h^2$ and $N_{\rm eff}$.  For the latter, Ref. <cit.> finds

\begin{equation}
N_{\rm eff} = 2.88 \pm 0.27, \quad   (68\%~{\rm CL})
\end{equation}

using the code PRIMAT <cit.>. Using PArthENoPE <cit.> instead would have yielded a central value about 2% smaller <cit.>.
Thus, the conclusion here is that current measurements of primordial elemental abundances are completely consistent with the $N_{\rm eff}^{\rm SM}=3.0440$, and shows no evidence of any extra relativistic degrees of freedom.

Cosmic Microwave Background and Large-Scale Structure

Unlike BBN, probes of the universe’s inhomogeneities such as the CMB temperature and polarization anisotropies
and the large-scale matter distribution are not sensitive to the flavor content of the neutrino
sector, only to its contribution to the stress-energy tensor. If neutrinos are massless, then the $N_{\rm eff}$
parameter as defined in Eq. <ref> alone characterizes their effects on the universe’s evolution.
If neutrinos are massive, then, in addition to $N_{\rm eff}$ and the neutrino masses, in principle it is also necessary to know the exact form of the neutrino
momentum  distribution in order to solve the evolution equations for the inhomogeneities exactly.
However, unless the deviations from a thermal relativistic Fermi-Dirac spectrum is of order unity, it suffices to specify only the temperature of the distribution, as late-time cosmological probes are currently not very sensitive to spectral distortions in the neutrino sector.  For this reason, most existing analyses characterize the neutrino sector, including light sterile neutrinos produced as described in Sec. <ref>, only in terms of the neutrino mass spectrum and the $N_{\rm eff}$ parameter.  The light sterile states are assumed to share the same temperature as the standard active neutrinos.

One interesting variation to the above is the case in which an excess $N_{\rm eff} > N_{\rm eff}^{\rm SM}$ is due to a thermalized particle species that has temperature different from the standard neutrino temperature and/or has a different spin statistics, e.g., a thermalized bosonic particle species such as an axion that follows the Bose-Einstein distribution.  If this new particle is massless, then again it suffices to describe its phenomenology in terms of its contribution to $N_{\rm eff}$ alone.  If however the particle species is massive, then in addition to its mass, its temperature also plays a role in determining the “hotness” of the resulting hot dark matter, where the relation between the hot dark matter's temperature and abundance is fixed by the particle's spin statistics.

Cosmic microwave background primary anisotropies

The CMB primary anisotropies refer to the temperature and polarization fluctuations imprinted on the last scattering surface.  These are sensitive to the physics of the early universe up to the time of photon decoupling ($T \sim 0.1$ eV or redshift $z \sim 1000$), and differ from the secondary anisotropies which are additional spatial fluctuations gathered by the CMB photons as they free-stream from the last scattering surface to the observer and hence sensitive to late-time/low-redshift physics.

The effects of a non-standard $N_{\rm eff}$ value on the CMB primary anisotropies and its associated parameter degeneracies with the present-day Hubble expansion rate $H_0$ and physical matter density $\omega_{\rm m}$ in the context of $\Lambda$CDM cosmology
have been discussed extensively in, e.g., Ref. <cit.>.  Broadly speaking, if the particles that make up $N_{\rm eff}$ are ultra-relativistic at the time of CMB formation, then a non-standard $N_{\rm eff}$ can manifest itself in the following ways:

* The redshift of matter-radiation equality $z_{\rm eq}$ controls the ratio of radiation to matter at the time of photon decoupling and hence the evolution of the potential wells.  This in turn affects the peak height ratios of the CMB temperature anisotropy spectrum.  With seven acoustic peaks measured by the Planck mission <cit.>, the equality redshift $z_{\rm eq}$ has been measured to percent level precision in $\Lambda$CDM-type cosmologies.  As a probe of $N_{\rm eff}$, however, we note that

\begin{equation}
\label{eq:zeq}
z_{\rm eq} =  \frac{\omega_{\rm m}}{\omega_{\gamma}} \frac{1}{+0.227 N_{\rm eff}}-1,
\end{equation}

where $\omega_\gamma$ is the present-day photon energy density.  In other words, $N_{\rm eff}$ is exactly degenerate with the physical matter density $\omega_{\rm m}$, and measuring $z_{\rm eq}$ alone does not determine $N_{\rm eff}$.

* The angular sound horizon $\theta_s$, which determines the CMB acoustic peak positions, is another quantity sensitive to $N_{\rm eff}$.  Defined as $\theta_s \equiv r_s/D_A$, where $r_s$ is the sound horizon at photon decoupling and $D_A$ is the angular diameter distance to the last scattering surface, the parameter dependence of $\theta_s$ in $\Lambda$CDM cosmologies is as follows,

\begin{equation}
\label{eq:thetas}
\theta_s \propto \frac{\Omega_{\rm m}^{-1/2}}{\int^1_{a^*} \frac{{\rm d} a}{a^2 \sqrt{\Omega_{\rm m}a^{-3} + (1-\Omega_{\rm m})}}},
\end{equation}

where $\Omega_{\rm m} \equiv \omega_{\rm m}/h^2$, $h$ is the reduced Hubble expansion rate defined via $H_0 = 100 h\, {\rm km}\, {\rm s}^{-1}\, {\rm Mpc}^{-1}$, $a^*$ is the scale factor at photon decoupling, and we have held $z_{\rm eq}$ and the baryon density $\omega_{\rm b}$ fixed.  This relation implies that while $\theta_s$ constrains the parameter combination $\omega_{\rm m}/h^2$, it does not constrain $\omega_{\rm m}$ and $h$ individually.  Since there already exists an exact degeneracy between $N_{\rm eff}$ and $\omega_{\rm m}$ through $z_{\rm eq}$ (see Eq. <ref>), this additional $(\omega_{\rm m},h)$-degeneracy through $\theta_s$ immediately sets up a three-way degeneracy between $N_{\rm eff}$, $\omega_{\rm m}$, and $h$, which needs to be broken by some other means.  In more complex models, degeneracies between $N_{\rm eff}$ and a nonzero spatial curvature $\Omega_k$ or a non-canonical dark energy equation of state are also possible.

* The angular diffusion scale $\theta_d \equiv r_d/D_A$, where $r_d$ is the diffusion scale at photon decoupling, characterizes the scale at which the CMB temperature anisotropy power spectrum becomes suppressed due to diffusion damping (or Silk damping).  The phenomenon of diffusion damping occurs at the CMB damping tail, i.e., at multipoles $\ell \gtrsim 1000$, and was first measured by the Atacama Cosmology Telescope <cit.> and the South Pole Telescope <cit.>, and now by the Planck CMB mission <cit.>.

For fixed $z_{\rm eq}$, $\omega_{\rm b}$ and $a^*$, the angular diffusion scale has a parameter dependence

\begin{equation}
\theta_d \propto (\Omega_{\rm m} H_0^2)^{1/4} \, \theta_s,
\end{equation}

where $\theta_s$ is the angular sound horizon of Eq. <ref>.  Thus, a simultaneous measurement of $\theta_d$, $\theta_s$, and $z_{\rm eq}$ by a CMB mission such as Planck immediately constitutes a measurement of $\omega_{\rm m} \propto \Omega_{\rm m} H_0^2$ and hence $N_{\rm eff}$.  Fig. <ref> shows the signature of $N_{\rm eff}$ in the CMB damping tail.

Signature of $N_{\rm eff}$ in the damping tail of the CMB TT power spectrum.  Top: Here, the cosmological parameters have been adjusted such that the baryon density $\Omega_b h^2$, the redshift of matter-radiation equality $z_{\rm eq}$, the angular sound horizon $\theta_s$, and the normalisation at $\ell=200$ are fixed for all $N_{\rm eff}$ cases shown.  Middle: Like the top panel, but with the spectrum normalisation fixed at $\ell=400$ instead.  Bottom: Like the middle panel, but here the Helium-4 fraction $Y_p$ is also allowed to vary such that all $N_{\rm eff}$ cases end up with roughly the same angular diffusion scale $\theta_d$.
Figure taken from Ref. <cit.>.

Large-scale matter distribution   There are many different ways to quantify and probe the large-scale matter distribution.  The most basic quantity, however, is the present-day matter power spectrum $P(k)$.

The broad-band shape of the large-scale matter power spectrum in $\Lambda$CDM-type cosmologies by two quantities: the comoving wavenumber at matter-radiation equality,

\begin{equation}
\begin{aligned}
k_{\rm eq} & \equiv a_{\rm eq} H(a_{\rm eq}) \\
& \simeq 4.7 \times 10^{-4} \sqrt{\Omega_{\rm m} (1+ z_{\rm eq})}\, h \, {\rm Mpc}^{-1},
\end{aligned}
\end{equation}

which fixes the location of the “turning point” of $P(k)$, and the baryon-to-matter density fraction

\begin{equation}
f_b \equiv \frac{\omega_{\rm b}}{\omega_{\rm m}},
\end{equation}

which determines the suppression in power at $k > k_{\rm eq}$ due to baryon acoustic oscillations (BAO).  In addition, the matter power spectrum has small-amplitude oscillatory features which are the manifestation of the BAO themselves.  These oscillatory features originate in the same early universe physics as the acoustic peaks in the CMB anisotropy power spectra, and when analyzed together with CMB data, can act as a powerful standard ruler for distance measurements.

If some of the neutrinos are massive and become non-relativistic at late times, then they can constitute a fraction of the present-day dark matter content.  However, this neutrino dark matter is “hot”, in the sense that the neutrinos, although non-relativistic, come with a significant thermal velocity dispersion, which tends to hinder its gravitational  clustering on small scales.  In terms of the present-day large-scale matter power spectrum, keeping the total matter density fixed but replacing some of the cold dark matter with neutrino hot dark matter suppresses $P(k)$ at $k$-values larger than the free-streaming wave number by an amount dependent on  the neutrino fraction,

\begin{equation}
f_\nu \equiv \frac{\omega_\nu}{\omega_{\rm m}} = \frac{\sum m_{\nu}/94\, {\rm eV}}{\omega_{\rm m}}.
\end{equation}

Thus, neutrino masses also influence the overall shape of the matter power spectrum.

Since it is already possible to pin down $z_{\rm eq}$, $\Omega_{\rm m}$, and $\omega_{\rm m}$ using the CMB primary anisotropies, measurements of the matter power spectrum $P(k)$ generally do not improve the constraint on $N_{\rm eff}$ in the simplest $\Lambda$CDM+$N_{\rm eff}$ fit.  However, it must be noted that after the formation of the CMB primary anisotropies, these fluctuations are gravitationally lensed by the intervening matter distribution as the CMB photons propagate from the last scattering surface to the observer, contributing to the so-called CMB secondary anisotropies at multipoles $\ell \gtrsim 500$.  In other words, any CMB anisotropy signal at $\ell \gtrsim 500$ will always include some information about $P(k)$, and this is particularly useful for the purpose of constraining the neutrino mass sum $\sum m_\nu$.

Current CMB and LSS constraints  In a standard $\Lambda$CDM parameter inference, estimating cosmological parameter values from the CMB and related observations involve varying six free parameters related to cosmology: the baryon density $\omega_b$, the cold dark matter density $\omega_c$, the Hubble parameter $h$, the spectral index $n_s$ and amplitude $A_s$ of the primordial curvature perturbation, and the optical depth to reionization $\tau$.  Analyses of the Planck CMB data also require that we vary of order 20 nuisance parameters to model the foregrounds and instrumental systematics.  These are later marginalized.

To constrain radiation excess, at minimum we need to add $N_{\rm eff}$ as a free parameter to this list.  Doing so the Planck collaboration finds<cit.>

\begin{equation}
N_{\rm eff} = 3.00^{+0.57}_{-0.53} \quad (95\%~{\rm CL},~{\rm Planck~TT+lowE}),
\end{equation}

\begin{equation}
N_{\rm eff} = 2.92^{+0.36}_{-0.37} \quad (95\%~{\rm CL},~{\rm Planck~TTTEEE+lowE}),
\end{equation}

\begin{equation}
N_{\rm eff} = 3.11^{+0.44}_{-0.43} \quad (95\%~{\rm CL},~{\rm Planck~TT+lowE+lensing+BAO}),
\end{equation}

\begin{equation}
\label{eq:neffconstraintcmb}
N_{\rm eff} = 2.99^{+0.34}_{-0.33} \quad (95\%~{\rm CL},~{\rm Planck~TTTEEE+lowE+lensing+BAO}),
\end{equation}

using various combinations of the Planck temperature and $E$-polarization measurements (TTTEEE and lowE), the lensing potential extracted from the Planck temperature maps, as well as the BAO measurements from 6dFGS <cit.>, SDSS-MGS <cit.>, and BOSS DR12 <cit.>.
As can be seen, in all cases, the inference returns an estimate of $N_{\rm eff}$ that is remarkably consistent with the SM prediction of $N_{\rm eff}=3.0440$.

Since the CMB anisotropies are also sensitive to the large-scale matter distribution at low redshifts because of the weak gravitational lensing signal inherent in all CMB power spectrum, one can also derive a constraint on the neutrino mass sum $\sum m_\nu$ from the Planck CMB data at the same time as we constrain $N_{\rm eff}$.
In a 8-parameter $\Lambda$CDM+$N_{\rm eff}$+$\sum m_\nu$ fit, the Planck collaboration finds <cit.>

\begin{equation}
\begin{aligned}
\label{eq:combinedneffmnu}
N_{\rm eff}&  = 2.96^{+0.34}_{-0.33} \quad (95\%~{\rm CL},~{\rm Planck~TTTEEE+lowE+lensing+BAO}), \\
\sum m_\nu & < 0.12~{\rm eV}.
\end{aligned}
\end{equation}

Note that the fit assumes three degenerate neutrino mass eigenstates of equal abundances.  The role of $N_{\rm eff}$ is merely to dial up or down the abundances, which is why $N_{\rm eff}$ can go below the standard value of $N_{\rm eff}^{\rm SM} =3.0440$.  This combined fit is to be compared with the $N_{\rm eff}$ constraint quoted in Eq. <ref> in a 7-parameter fit of the same data combination, which has the same error bars (about 11%) and a central value off only by 2%.  It is also interesting to compare it with the constraint obtained on the neutrino mass sum $\sum m_\nu$ from 7-parameter $\Lambda$CDM+$\sum m_\nu$ fit <cit.>,

\begin{equation}
\sum m_\nu < 0.13~{\rm eV} \quad (95\%~{\rm CL},~{\rm Planck~TT+lowE+lensing+BAO}),
\end{equation}

\begin{equation}
\sum m_\nu < 0.12~{\rm eV} \quad (95\%~{\rm CL},~{\rm Planck~TTTEEE+lowE+lensing+BAO}),
\end{equation}

which is identical to the upper limit obtained from the 8-parameter fit (Eq. <ref>).
Thus, one can conclude from this comparison that  there is no strong degeneracy between the $N_{\rm eff}$ parameter and the neutrino mass sum $\sum m_\nu$ in the current generation of precision cosmological data.

What about constraints on the sterile neutrino mass $m_s$?  The Planck collaboration <cit.> also reports a “massive sterile neutrino” fit to their data in a scenario in which the three active neutrinos are assumed to have a fixed minimum mass sum of $\sum m_\nu=0.06$ eV and the effective sterile neutrino mass is defined as $m_{\nu, {\rm sterile}}^{\rm eff} \equiv \Omega_{\nu, {\rm sterile}} h^2 (94.1~{\rm eV})$. This effective mass is related to the physical sterile neutrino mass via

\begin{equation}
\label{eq:effsterile}
m_s = (\Delta N_{\rm eff})^{-1} m_{\nu, {\rm sterile}}^{\rm eff},
\end{equation}

assuming the sterile states have the same temperature as the SM neutrinos.
Imposing the priors  $\Delta N_{\rm eff} \geq 0$ and  $m_s < 10$ eV, they find the constraints

\begin{equation}
\begin{aligned}
\label{eq:effsterileconstraint}
N_{\rm eff}&  < 3.29 \quad (95\%~{\rm CL},~{\rm Planck~TTTEEE+lowE+lensing+BAO}), \\
m_{\nu, {\rm sterile}}^{\rm eff} & < 0.65~{\rm eV}.
\end{aligned}
\end{equation}

Taking $\Delta N_{\rm eff}$ to be a maximum allowed $0.29$, the mass bound  corresponds to an upper limit of $m_s < 2.24$ eV on the physical sterile neutrino mass.
Thus, while cosmological measurements do constrain $m_s$ in any interesting way, it cannot completely rule out a 1 eV light sterile provided thermalization is kept at below the 30% level.  In other words, the extent of thermalization as quantified by the $N_{\rm eff}$ parameter remains the limiting factor for the light sterile neutrino scenario in standard cosmology.

Lastly, we note that while the cosmological bound on $N_{\rm eff}$ given in Eq. <ref> is already strongly indicative that the short-baseline light sterile neutrino is in serious tension with precision cosmological measurements, it is nonetheless possible to analyse cosmological and oscillation data together in a consistent way.  This has been done most recently in Ref. <cit.>, which considers a 3+1 scenario, computes the corresponding light sterile neutrino thermalization using the thermalization code FortEPiaNO <cit.> and the associated $N_{\rm eff}$ with mixing parameters consistent with laboratory measurements, and feeds the output into a CMB analysis. Fig. <ref>
shows the constraints on the $(\Delta m_{41}^2, |U_{e4}|^2)$- and $(\Delta m_{41}^2, |U_{\mu4}|^2)$-planes.

Left: 2D marginalized 68% and 95% constraints on the mass splitting $\Delta m^2_{41}$ and mixing matrix element $|U_{e4}|^2$ from cosmology (blue), from the tritium $\beta$-decay end-point measurements by KATRIN (green), and from neutrinoless double-$\beta$-decay experiments (red). The preferred regions of reactor experiments are also indicated. Right: Cosmological 68% and 95% marginalized constraints on the $|U_{\mu 4}|^2$ (blue) versus constraints from the $\nu_\mu$ disappearance measurements of IceCube and MINOS+ (grey).
Figure taken from Ref. <cit.>.

§.§.§ Can We Evade Cosmological Constraints?

We have seen that the canonical 1 eV-mass light sterile neutrino motivated by the short baseline anomalies is in strong tension with cosmological measurements primarily because of the non-detection of a non-standard $N_{\rm eff} \sim 4$.
However, a criticism often levelled at cosmological constraints is that all parameter estimation via statistical inference are inherently dependent on the cosmological framework assumed in the inference exercise.  Given the large number of unknowns in cosmology, e.g., the nature of dark energy, inflation, etc., critics argue that there may exist a  corner of this vast unknown parameter space in which a completely thermalized light sterile neutrino state with a mass close to 1 eV might be permitted to live.

To investigate this possibility, the Planck collaboration has provided a large number of analyses of expanded cosmological parameter spaces, often in combination with external, non-CMB data sets.  These are available at the Planck Legacy Archive (https://pla.esac.esa.int/home).   Of particular interest to light sterile neutrinos  is the 10-parameter $\Lambda$CDM+$N_{\rm eff}$+$\sum m_\nu$+$w_0$+$n_{\rm run}$ fit, where $w_0$ is the equation of state parameter of the dark energy ($w_0=1$ for a cosmological constant), and $n_{\rm run}$ is the running of the scalar spectral index, a parameter related to the initial conditions of the universe ($n_{\rm run}=0$ in single-field inflation). Using the usual Planck TTTEEE+lowE+lensing+BAO data combination together with the Hubble parameter measurement of Ref. <cit.> and Supernova Ia data from the Pantheon sample <cit.>, the constraints on $N_{\rm eff}$ and $\sum m_\nu$ in this extended parameter fit are

\begin{equation}\label{eq:effsterileconstraint2}
\begin{aligned}
N_{\rm eff}&  = 3.11^{+0.37}_{-0.36} \quad (95\%~{\rm CL}), \\
\sum m_\nu & < 0.16~{\rm eV}.
\end{aligned}
\end{equation}

Clearly, while the bound on the neutrino mass sum $\sum m_\nu$ has relaxed somewhat and the central $N_{\rm eff}$ has shifted a little up, relative to the more limited bounds in Eq. <ref> and  <ref> the error bars on $N_{\rm eff}$ have not weaken significantly, and
the canonical light sterile neutrino is still in tension with precision cosmological measurements in this expanded parameter space.  The upward shift in $N_{\rm eff}$ can be attributed to the discrepancy between the Planck inference of the Hubble expansion rate and the local measurement of Ref. <cit.>.  Because of the degeneracy between $N_{\rm eff}$ and $H_0$ in the CMB primary anisotropies, combining Planck data with local measurements—the latter of which prefer a higher value of $H_0$—tends to drag up the inferred $N_{\rm eff}$ as well.  The neutrino mass sum $\sum m_\nu$, on the other hand, has long been known to be somewhat degenerate with the dark energy equation of state parameter $w_0$.  However,  the combination of BAO and Supernova Ia data can lift this degeneracy very effectively.

Thus, expanding the cosmological parameter space no longer appears to do much for the light sterile neutrino case (in the sense of allowing a larger $N_{\rm eff}$) the way it once did <cit.>.  In order to get around cosmological constraints, we need to introduce new physics that directly affects the cosmological phenomenology of the light sterile states.  Since the main problem of the canonical light sterile neutrino is that its thermalization in the early universe raises $N_{\rm eff}$ to an unacceptably large level for BBN and CMB/LSS, all known new physics solutions so far involve tampering with the  thermalization process, in order to maintain $N_{\rm eff}$ at as close to the SM value as possible.  A number of ideas have been proposed and explored throughout the years (though not all are guaranteed to work as desired), including

* Large chemical potentials or, equivalently, number density asymmetries for the active neutrinos <cit.>,

* Secret interactions of the sterile neutrinos <cit.>, and

* Low reheating temperature of the universe <cit.>.

Large chemical potentials for the active neutrinos As discussed earlier in Sec. <ref>, in a standard calculation of light sterile neutrino thermalization, the active neutrino asymmetries, defined as
$L_\alpha =(n_{\nu_\alpha}- n_{\bar{\nu}_\alpha})/n_\gamma$, are assumed to be zero.  However, if for some reason some of these asymmetries are large—usually taken to mean $L_\alpha > 10^{-5}$—then the CP asymmetric term $\sqrt{2} G_F(\mathbb{N}_\nu-\mathbb{N}_{\bar\nu})$ in the Hamiltonian (Eq. <ref>) can act as a large matter effect to suppress oscillations between the active and sterile states.  If this suppression is effective before neutrino decoupling, then it is possible to maintain $N_{\rm eff}$ at close to the SM prediction <cit.>.
For active-sterile neutrino mass splittings in the range $\Delta m^2 \simeq 0.2 \to 10~{\rm eV}^2$, the minimum neutrino asymmetries required to effect some degree of suppression are $L > 10^{-4} \to 5 \times 10^{-3}$ <cit.>. To significantly suppress thermalization, however, asymmetries as large as $L \sim 10^{-2}$ are required <cit.>.

Unfortunately, aside from the difficulty in explaining how such large neutrino asymmetries could have arisen in the first place, this solution also suffers from other undesirable effects, namely, significant distortion to the $\nu_e$ and $\bar{\nu}_e$ energy spectra.  While it is possible to suppress active-sterile oscillations with the choice of $L \sim 10^{-2}$ before neutrino decoupling, beyond this critical point vacuum oscillations will inevitably take over and distort the active neutrino energy spectra as a result <cit.>.  Ref. <cit.> has computed this distortion and its effect on the light elemental abundances.  They find that the larger the neutrino asymmetry employed to suppress light sterile neutrino thermalization, the larger the spectral distortion and the resulting Helium-4 mass fraction $Y_p$.  Thus, while large neutrino asymmetries do improve the outcome for $N_{\rm eff}$, at the same time they also affect at least one important observable in an undesirable way.  The solution is therefore far from fool-proofed.

Self-interaction or non-standard interaction for the sterile neutrino These solutions also work on the principle of suppressing sterile neutrino thermalization through the introduction of a non-standard matter potential for the sterile state in the oscillation Hamiltonian (Eq. <ref>).  They differ primarily in their coupling structures.

Ref. <cit.> consider an interaction of the form

\begin{equation}
\label{eq:vectorinteraction}
{\cal L}_{\rm int} = g_X \bar{\nu}_s \gamma^\mu P_L \nu_s X_\mu,
\end{equation}

where the sterile neutrino self-interaction is mediated by a MeV-mass vector boson $X$.  This leads to the addition of a matter potential to the Hamiltonian (Eq. <ref>) of the form

\begin{equation}
\mp 2 \sqrt{2} G_X p \left[\frac{4}{3}\frac{\mathbb{E}_{s}+\mathbb{E}_{\bar{s}}}{m_X^2}\right],
\end{equation}

where $G_X \equiv (\sqrt{2}/8) g_X^2/m_X^2$, $m_X$ is the mass of the $X$ boson, and
$\mathbb{E}_s$ and $\mathbb{E}_{\bar{s}}$ are defined like $\mathbb{E}_\nu$ and $\mathbb{E}_{\bar{\nu}}$ in Eq. <ref>, but with a projection matrix $S_a={\rm diag}(0,0,0,1)$ that singles out the sterile state for coupling.  The collision integrals ${\cal I}[\varrho(p,t)]$ also need to be modified appropriately.
For $g_X \gtrsim 10^{-2}$ and $g_X \lesssim 10$ MeV, the new interaction (Eq. <ref>) is able to suppress sterile neutrino thermalization and hence preserve $N_{\rm eff}$ at close to the SM value without altering BBN predictions.  However, at times after neutrino decoupling, the same interaction also leads to equipartition amongst the active and sterile states.  That is, if at neutrino decoupling the neutrino number densities are $(n_{e}, n_{\mu}, n_{\tau}, n_{s}) = (1,1,1,0)$, the secret interaction will redistribute it to $(3/4, 3/4, 3/4,3/4)$.  Thus, the mass of sterile state $m_s$ will nonetheless contribute to the hot dark matter energy density probed by the CMB anisotropies and the large-scale matter distribution, and be subjected
a substantially tighter constraint, about $m_s < 0.2$ eV, than is implied by Eq. (<ref>) or Eq. (<ref>).
This solution is therefore also not fool-proofed.

On the other hand, Ref. <cit.> consider a self-interaction of the mass eigenstate $\nu_4$ mediated by a massless pseudoscalar $\phi$:

\begin{equation}
\label{eq:pseudoscalarinteraction}
{\cal L}_{\rm int} = g_\phi \phi \bar{\nu}_4 \gamma_5 \nu_4.
\end{equation}

As with the massive vector boson case above, the secret interaction engenders a matter potential, which in turn suppresses the production of sterile states.  The transition from no to full thermalization happens in the range of coupling values  $10^{-6} < g_\phi < 10^{-5}$ <cit.>.  The authors further argue that because the secret interaction happens exclusively for the mass eigenstate $\nu_4$, the interaction cannot equilibrate the active and sterile states and whatever is the $\Delta N_{\rm eff}$ produced at neutrino decoupling is also the only component of the neutrino population that carries a mass of $m_s \simeq 1$ eV.   Thus, the scenario can easily evade both limits on $N_{\rm eff}$ and $m_s$.

Lastly, the solution of Ref. <cit.> invokes a coupling of the sterile state to an ultra-light real scalar field  of mass $m_\phi < 5 \times 10^{-17}$ eV that also contributes to the cold dark matter.  At early times, the coupling induces an effective mass for the sterile state, which suppresses active-sterile oscillations and hence thermalization of the sterile state in much the same way as the two scenarios discussed above.  After neutrino decoupling, the $\phi$ field starts to oscillate coherently.  Since unlike the MeV vector boson case the model does not lead to the equilibration of the active and sterile states, both limits on $N_{\rm eff}$ and $m_s$ can be easily avoided.

Low reheating temperature Low reheating temperature scenarios <cit.> refer to those cases in which the universe transitions to radiation domination at temperatures below $T \sim 10$ MeV.  This might happen because of a very low inflation energy scale, or because some  non-standard physics causes the universe to enter a period of matter domination immediately prior to the most recent phase of radiation domination, and the transition back to radiation domination takes place at $T \sim {\cal O}(1)$ MeV.

A low reheating temperature appears to be a viable way to evade cosmological constraints on light sterile neutrino states.  If reheating occurs at $T\sim {\cal O}(1)$ MeV, even the SM active neutrinos have barely enough time to interact  before neutrino decoupling happens.  Depending on how exactly reheating happens, some of the active neutrino species may not even reach equilibrium number or energy densities.  If light sterile neutrino thermalization was to happen at the same time, the shortage of active neutrinos in the plasma would also slow down the production rate.  For these reasons, it is possible to engineer a scenario in which the final $N_{\rm eff}$ is close to the SM value, while the ratio of sterile to active states remains smaller than 1 to 3.  Naively, this makes it possible to satisfy $N_{\rm eff}$ as well as $m_s$ bounds from CMB and/or BBN.

In practice, however, whether or not the solution works depends on the details of the reheating model.  Reference <cit.>, for example, finds that if the parent particle responsible for reheating
decays exclusively into electromagnetically interacting radiation, then a low temperature reheating can indeed render light sterile neutrinos consistent with measurements of the primordial elemental abundances.
If however the parent particle decays mainly into hadrons, then together with the presence of active-sterile neutrino mixing, the primordial synthesis of light elements can proceed in a way incompatible with observations for a wide range of the mass and the hadronic branching ratio of the parent particle.

In addition to low reheating scenarios, cosmological scenarios where entropy is conserved and the expansion rate is modified can significantly affect light sterile neutrino constraints. This includes, for example, scalar-tensor theories, and is discussed in detail in Refs. <cit.>, and also in <cit.> (for resonance sterile neutrino production).

§ FUTURE EXPERIMENTAL PROSPECTS

§.§ Decay-at-Rest Accelerator Experiments

§.§.§ JSNS$^2$ and JSNS$^2$-II

The JSNS$^2$ (J-PARC Sterile Neutrino Search at the J-PARC Spallation Neutron
Source) <cit.> and its second phase
JSNS$^2$-II <cit.>, aim to search
for neutrino oscillations with $\Delta m^2$ near
1 eV$^2$ at the J-PARC Materials and Life Science Experimental
Facility (MLF). Fig. <ref> shows the experimental setup and
search sensitivities.

The experimental setups of the JSNS$^2$ and JSNS$^2$-II
and their sensitivities.

An intense neutrino beam from muon decay at rest is produced by a
spallation neutron target with the 1 MW beam of 3 GeV protons
created by a Rapid Cycling Synchrotron (RCS).
Neutrinos come predominantly from $\mu^+$ decay :
$\mu^{+} \to e^{+} + \bar{\nu}_{\mu} + \nu_{e}$.
An oscillation of $\bar{\nu}_{\mu}$ to $\bar{\nu}_{e}$ through the
fourth mass eigen-state is searched for by detecting the Inverse-Beta-Decay (IBD)
interaction $\bar{\nu}_{e} + p \to e^{+} + n$, followed by
gammas from neutron capture on Gd.
The JSNS$^2$ detector, as the near detector in the JSNS$^2$-II setup,
contains 17 tonnes of Gd-loaded liquid scintillator and is located
24 meters away from the mercury target. The new far detector
of JSNS$^2$-II, currently under construction, is located outside
the MLF building with a baseline of 48 meters.
The far detector contains 32 tonnes of Gd loaded liquid scintillator
as a neutrino target.
Both JSNS$^2$ and JSNS$^2$-II employ a Hydrogen target for the neutrino source ($\mu$ decay-at-rest), and
a neutrino detection channel (IBD) identical to LSND.
With improvements made by the short pulsed beam and the neutron
capture signal, JSNS$^2$ and JSNS$^2$-II will provide
clean and direct tests of the LSND anomaly.
JSNS$^2$ started data taking in 2020, and accumulated
1.45$\times$10$^{22}$ Proton-On-Target (POT) by 2021, 13$\%$ of the
approved POT by J-PARC, which corresponds to 1 MW beam power for 3 years.
An extensive analysis is ongoing (e.g. <cit.>).
The search sensitivity of JSNS$^2$ with the full design POT
is shown in the bottom-middle plot of Fig. <ref>.
The construction of the far detector of JSNS$^2$-II started
in September 2021 and the aim is to start data taking in 2023.
It will provide additional sensitivity, especially in the low $\Delta m^2$ region with
1 MW beam power for 5 years.

Figure <ref> shows a picture of the JSNS$^2$ detector
during the installation in June 2020 and
the construction status of the new far detector of the JSNS$^2$-II.

The JSNS$^2$ detector during the installation (left) and
the construction status of the new far detector of the JSNS$^2$-II

§.§.§ COHERENT at the SNS

The COHERENT collaboration can perform a powerful test of oscillations of sterile neutrinos by considering NC disappearance.  Parameter space favored by a global fit of oscillation data to a 3+1 scenario is accessible to COHERENT in the near future with later data giving a much stronger constraint.

COHERENT measures coherent, elastic neutrino-nucleus scattering (CEvNS) and other low-energy neutrino scattering processes at the Spallation Neutron Source (SNS) at Oak Ridge National Lab (ORNL).  CEvNS is a neutral current process whose only signature is a low-energy nuclear recoil which was first measured by COHERENT on CsI in 2017 <cit.>.  The cross section is very large compared to other neutrino scattering cross sections below 50 MeV and is precisely predicted.  The SNS is an intense source of $\pi$ decay-at-rest neutrinos with energies 0-53 MeV, ideal for measuring CEvNS.  The width of the SNS beam, 360 ns FWHM, is small compared to the muon lifetime so that the neutrino flux separates in time to a prompt $\nu_\mu$ flux from $\pi^+\rightarrow\mu^+\nu_\mu$ and a delayed flux of $\bar{\nu}_\mu/\nu_e$ from $\mu^+\rightarrow e^+\bar{\nu}_\mu\nu_e$.  The $\nu_\mu$ flux is monoenergetic with $E_\nu=29.8$ MeV.  COHERENT builds and commissions several CEvNS detectors for operation at the SNS at baselines of 19.3 to 28 m.  These baselines place COHERENT detectors at the first oscillation maximum for the $\nu_\mu$ flux assuming the global best fit of $\Delta m^2_{41}$.  Since the neutrino flux at the SNS includes both $\nu_\mu/\bar{\nu}_\mu$ and $\nu_e$, CEvNS searches can simultaneously search for $\nu_\mu\rightarrow\nu_s$ and $\nu_e\rightarrow\nu_s$ disappearance with favorable sensitivity to both $\theta_{14}$ and $\theta_{24}$ with the same experiment.  Additionally, the largest systematic uncertainty, the neutrino flux normalization, is correlated between all detectors which mitigates its effect on a joint fit.

There are plans for three future detectors suitable for searching for sterile neutrinos through CEvNS disappearance in the near future.  The first is an upgrade of the CENNS10 detector which made the first CEvNS measurement on argon <cit.>.  This will be a liquid argon calorimeter with 610 kg of fiducial mass with a baseline of 28 m.  Detector performance is well understood from experience with CENNS10 operations and data.  A 10-kg CsI scintillation detector at a 19.3 m baseline is also planned.  This detector will be undoped and cooled to 77 K which can dramatically increase light yield while reducing background scintillation within the crystal <cit.>.  This ensures a low threshold, allowing tests of CEvNS disappearance with low-energy recoils.  Finally, a 50 kg germanium PPC detector at 22 m is planned as an upgrade to the 17 kg array currently being commissioned at the SNS.  The sensitivity of a joint fit using three years of data from all three detectors to search for a sterile neutrino through NC disappearance with CEvNS is shown in Fig. <ref>.  This would test the parameter space preferred by a global fit at 90$\%$ confidence <cit.>.

Sensitivity of COHERENT CEvNS detectors to constrain sterile neutrino parameter space assuming a 3+1 model compared to the LSND and MiniBooNE allowed regions.  A global fit to all short-baseline oscillation data is also shown.

Into the next decade, ORNL is investing in the SNS, doubling its power and constructing a new second target station (STS) with one in every four beam spills delivered to the STS to supplement work at the first target station (FTS).  Though the upgrade will not be completed until the 2030s, it will facilitate a strong search for sterile neutrinos.  The two targets would only be 140 m apart, meaning a large flux of neutrinos from both sources would pass through each CEvNS detector at the SNS.  Similar to two-detector oscillation experiments, this mitigates systematic uncertainties from neutrino interaction modeling and detector response by observing neutrinos from a near and far flux source.  The dominant remaining uncertainty comes from $\pi^+$ production in each beam target which is small.  With this control of uncertainties, a test of NC disappearance is possible at the 1$\%$-level.  Sensitivity to sterile neutrino oscillations in a 3+1 framework for a 10 t fiducial argon calorimeter running for five years, when placed 20 m from the STS and 120 m from the FTS, is shown in Fig. <ref>.  Exploiting flux from both targets, this large detector could test $\sin^22\theta_{\mu e}$ values of $10^{-5}$ at the global best fit $\Delta m^2_{41}$ and could test the LSND and MiniBooNE preferred regions for $\Delta m^2_{41}>0.04$ eV$^2$ at high confidence.

§.§.§ Coherent CAPTAIN-Mills

The physics program of the Coherent Captain Mills (CCM) Experiment comprises searches for new particles in the weak sector, including Dark Photons, Axion-like Particles (ALPs), and heavy neutral leptons in the keV to MeV mass range, extending the coverage of open parameter space for these searches at the order of magnitude level.
Many of these particles are invoked as alternative or additional explanations to oscillations involving sterile neutrinos as the source of MiniBooNE anomaly.
Thus, the results of CCM from the ongoing run at Los Alamos National Laboratory (LANL) have direct bearing on phenomenology presented in this white-paper.
Here, we describe the CCM detector, present a relevant CCM search for production of new bosons by charged meson decays <cit.> as an example of the impact of the results, and summarize other searches that can be performed.

The CCM experiment is relatively new to the scene of experiments to understand the phenomenology of short-baseline anomalies.
The experiment was conceived in 2017 and prototyped using “CCM120,” which tested 120 PMTs for the SBND liquid argon (LAr) experiment.     First physics results from CCM120 were recently published <cit.>.
In 2019, the LANL LDRD office and DOE Dark Matter New Initiative program recognized the relevance of the CCM rare-particle searches to dark matter studies and provided funding for an upgrade to 200 PMTs for “CCM200.”
This 5 t fiducial-volume (10 t total) LAr detector with 50% photomultiplier tube (PMT) coverage, seen in Fig. <ref>, was completed in Autumn 2021.
The detector is unusual for accelerator-based liquid argon experiments, in that it utilizes only light collection—no time projection chamber.
The detector is being commissioned now, and data totaling $2.25 \times   10^{22}$ POT will be collected in three runs between 2022-24, at the Lujan spallation neutron center.
This facility targets 100 microamps on tungsten of  800 MeV  protons with 275 ns spills at 20 Hz.
This is a prolific source of neutrinos
from stopped pion and muon decay, and, potentially, a source for production of new particles, such as ALPs <cit.>, that can be observed in CCM200 through interactions or decays.
CCM200 is located 90$^\circ$ off-axis and 20 m from the target.

The interior of the CCM200 detector.
Of the 200 PMTs, 80% are coated with wavelength shifter (TBP) to shift 128 nm scintillation light to the visible, leaving 20% uncoated (darker, more reflective PMTs in image), aiding discrimination of Cherenkov light.
TPB foils cover the walls.
The light-tight interior is surrounded by a veto region instrumented with PMTs.
CCM is now running.

The CCM200 design has a combination of features related to light collection that makes it powerful and unique.
The first is its PMT coverage (8" Hammamatsu R5912-MOD), which is orders of magnitude higher per unit volume than any LAr TPC experiment.
Furthermore, the large charge dynamic range of the PMT's and electronics enables the detector to have reconstructed energy sensitivity from  10 keV to over 200 MeV.
The second is the rate of PMT readout, which provides information at 500 MHz and is synced to the accelerator providing 2-ns absolute timing relative to the 275 ns beam pulse.
This is key for separating out early speed of light particles from the prolific beam related neutrons.
Third, as can be seen in the Fig. <ref>, 80% of the PMTs are coated with 1,1,4,4-Tetraphenyl-1,3-butadiene (TPB), while 20% are uncoated.
The TPB allows observation of the scintillation light from LAr, which emits 40,000 photons/MeV in zero electric field–$\times 4$ brighter than typical oil-based liquid scintillator— at 128 nm wavelength, by shifting to the visible to penetrate the PMT glass.
The uncoated tubes, which is unique to the CCM design,  allows clean observation of Cherenkov light.
An R&D goal of CCM is to make the first use of observed Cherenkov light on an event-by-event basis in an analysis to reject backgrounds, since, for a signal, the direction of the Cherenkov ring is opposite the Lujan target.

Left: example diagrams for three-body meson decay producing scalar or pseudoscalar(axion-like) particles that interact in the detector.
Right:  Allowed regions for scalers (green) and pseudoscalers (blue) that account for the MiniBooNE low energy excess, presented as a function of coupling versus new particle mass <cit.>.
Lines: predictions of 10,000, 1000, and 100 signal events in solid, dashed and dotted, respectively.

These design features make CCM particularly ideal for searches for electromagnetic signatures of new physics produced in the target,  which is a signature of popular explanations for the short-baseline anomalies.
The most popular new-physics explanation has been oscillations involving sterile neutrinos.
CCM can explore the recent large-mixing angle result from the BEST experiment <cit.> using $\nu_e$ disappearance for the pion decay-at-rest neutrino beam, since the threshold for $\nu_e$-argon scatters is 1.5 MeV.
In the longer term, an upgraded CCM complex can be modified to perform a two-detector search for $\nu_\mu \rightarrow \nu_e$ in the LSND range that may be motivated by JSNS$^2$ results <cit.>.
However, recently, community interest has turned to new particles to explain the observed anomalies.
Motivated by this, for this short review, we are featuring an example of CCM's new-particle-discovery capability.

As an example of an interesting new model that CCM can address, consider the proposed explanation of the MiniBooNE Low Energy Excess (LEE) <cit.> from three-body meson decay <cit.>.
The diagrams for production of a new scalar or pseudoscalar particles that will interact in the detector to produce a single photon  exchanging a light vector boson ($Z'$) with the nucleus are shown in Fig. <ref>, left.
The allowed region for the LEE is shown in Fig. <ref>, right, and this model (model 1) finds a good fit to both the angular and energy dependence of the LEE <cit.>.
Another version of this model (model 2) which also can fit the data involves the emission a light vector boson $Z'$ from the pion decay (just like the scalar/pseudoscalar).
$V$ will then  produce a photon by exchanging  a scalar with the nucleus at the detector.
This model explains the LEE. All these new mediators emerging from the charged pion decays, so far we have discussed, can be coupled to only quarks.

CCM will be able to probe both models.
The main production channel will be $\pi^0$ decay into  gamma and $Z'$.
In model 1, $Z'$ can be assumed to dominantly decay into a pair of scalars, which will subsequently produce a photon from the scattering at the detector as needed to resolve LEE.
$Z'$ also can decay into a photon and the scalar.
Here, our assumption is that $Z'$ does not decay into a pair of visible particles promptly.
Under this assumption, Fig. <ref> shows the  predicted number of signal events at CCM  in the allowed  parameter space  to explain the LEE.
Therefore, CCM200 has the capability for discovery, if this new physics is the source of the LEE.
Model 2 probes the MiniBooNE anomaly more directly since the same $Z'$ that emerges from the charged pion decay to address the anomaly also can be produced from the $\pi^0$ decay; therefore no assumption is needed to correlate the LEE and a possible signal at the CCM.
This possibility is under investigation at present.

Outside of models explaining the anomalies, CCM engages in a broad range of new-physics searches.
Limits on leptophobic dark photons from CCM120 <cit.>, the prototype run, will be extended by two orders of magnitude in CCM200.
Searches for the QCD axion can close the last remaining open-window at masses $>0.1$ MeV <cit.>.   Searches for “neutrissimos”–not so heavy neutral heavy leptons–that have focused on $>100$ MeV masses to address the LEE are being extended to lower masses in CCM.
Although these searches are not directly tied to explanations of the anomalies, a discovery would inevitably demand investigation on whether the observed new-physics is related.

In summary,  CCM is a small, fast-timescale experiment that is already taking data at LANL.
Its results have the potential to change our thinking about the anomalies.

§.§.§ PIP2-BD: GeV Proton Beam Dump at Fermilab’s PIP-II Linac

The completion of the PIP-II superconducting LINAC at Fermilab as a proton driver for DUNE/LBNF in the late 2020s creates an attractive opportunity to build a GeV proton beam dump facility at Fermilab dedicated to and designed from the ground up for HEP with excellent sensitivity to eV-scale sterile neutrinos via neutral current disappearance using the CEvNS reaction. Thus, relative to spallation neutron facilities tailored to neutron physics and optimized for neutron production operating at a similar proton beam power, a HEP-dedicated beam dump facility would be designed to suppress rather than maximize neutron production and implement a beam dump made from a lighter target such as carbon, which can have a pion-to-proton production ratio up to $\sim$2 times larger than heavier Hg or W targets. The facility could also accommodate multiple, 100-ton-scale high energy physics experiments located at different distances from the beam dump and at different angles with respect to the incoming proton beam.  This flexibility further improves the sensitivity of dark sector and sterile neutrino searches, by allowing relative measurements at different distances and angles to constrain uncertainties in expected signal and background rates.

The continuous wave capable PIP-II LINAC at Fermilab can simultaneously provide sufficient protons to drive megawatt-class $\mathcal{O}$(GeV) proton beams as well as the multi-megawatt LBNF/DUNE beamline. By coupling the PIP-II LINAC to a new Booster-sized, permanent magnet or DC-powered accumulator ring, the protons can be compressed into pulses suitable for a proton beam dump facility with a rich physics program. The accumulator ring could be located in a new or existing beam enclosure and be designed to operate at 800 MeV but with an upgrade path allowing for future operation in the GeV range. The accumulator ring would initially provide 100 kW of beam power, limited by stripping foil heating, and have a $\mathcal{O}(10^{-4})$ duty factor. One variant of this accumulator ring would be a $\sim$100 m circumference ring operating at 1.2 GeV with a pulse width of 20 ns and a duty factor of $\mathcal{O}(10^{-6})$, which would greatly reduce steady-state backgrounds.  Another is an accumulator ring coupled to a new rapid cycling synchrotron replacing the Fermilab Booster with an increased proton energy of 2 GeV and an increased beam power of 1.3 MW <cit.>.

PIP2-BD 90% confidence limits on active-to-sterile neutrino mixing compared to existing $\nu_\mu$ disappearance limits from IceCube <cit.> and a recent global fit <cit.>, assuming a 5 year run (left). Also shown are the 90% confidence limits for $\nu_\mu$ disappearance (left), $\nu_e$ disappearance (middle), and $\nu_e$ appearance (right), assuming the $\bar\nu_\mu$ and $\nu_e$ can be detected with similar assumptions as for the $\nu_\mu$.

Decay-at-rest neutrinos from a stopped pion beam dump provide an excellent source of $\nu_\mu$, $\bar\nu_\mu$, and $\nu_e$ with a time structure that can separate $\nu_\mu$ from $\bar\nu_\mu$ and $\nu_e$.

CEvNS provides a unique tool to definitively establish the existence of sterile neutrinos through active-to-sterile neutrino oscillations <cit.>.  Using CEvNS, we can explore both mono-energetic $\nu_\mu$ disappearance with $\text{E}_\nu$ = 30 MeV and the summed disappearance of $\nu_\mu$, $\bar\nu_\mu$, and $\nu_e$ to $\nu_S$, which can also put constraints on $\nu_\mu \rightarrow \nu_e$ oscillation parameters in a 3+1 sterile neutrino model. We consider here a setup consisting of identical 100-ton LAr scintillation detectors, located 15 m and 30 m away from a carbon proton beam dump with a 20 keV recoil energy threshold and an efficiency of 70%.  The 100-ton scale scintillation-only detector assumes a cylindrical volume with a 5 m height and 2.5 m radius, and we perform a full Geant4-based <cit.> scintillation photon simulation with wavelength shifting and propagation to photomultiplier tubes along the endcaps and side walls of the detector. Based on simulation studies, we assume the detectors have a position resolution given by $\sigma = \frac{40 ~\textnormal{cm}}{\sqrt{{T}/20~\textnormal{keV}}}$ in each spatial dimension, where T is the nuclear recoil produced by the neutrino interaction within the detector. This information allows the possibility for a “rate+shape” fit using five bins in the reconstructed neutrino propagation distance with a bin width of 1 m matched to the expected resolution of the reconstructed neutrino propagation distance. If the prompt $\nu_\mu$ can be separated from the delayed $\bar\nu_\mu$ and $\nu_e$, one can exploit the mono-energetic feature of the $\nu_\mu$ flux and perform a joint rate + shape disappearance fit of CEvNS events in the near and far detectors as a function of reconstructed position.

We consider here a setup consisting of identical 100-ton LAr scintillation detectors, located 15 m and 30 m away from a carbon proton beam dump with a 20 keV recoil energy threshold and an efficiency of 70%.  We assume the neutron background in this dedicated facility could be suppressed to a negligible level for this experiment and that the signal-to-noise ratio for the remaining steady-state backgrounds is 1:1. In Fig. <ref>, we calculate the 90% confidence limits on the $\nu_\mu \rightarrow \nu_S$ mixing parameter $\sin^22\theta_{\mu S}$ for a 5-year run of an upgraded 1.2 GeV proton accumulator ring operating with a pulse width of 20 ns, a duty factor of $\mathcal{O}(10^{-6})$, and a 75% uptime, assuming a 9% normalization systematic uncertainty correlated between the two detectors and a 36 cm path length smearing. Also shown are the 90% confidence limits for $\nu_\mu$ disappearance, $\nu_e$ disappearance, and $\nu_e$ appearance, assuming the $\bar\nu_\mu$ and $\nu_e$ can be detected with similar assumptions as for the $\nu_\mu$.

§.§.§ KPIPE at Fermilab

The KPIPE experimental concept, outlined in Ref. <cit.>, calls for a very long (120 m) and thin (1.5 m radius) cylindrical detector close to and oriented radially outward from an intense beam-dump source of monoenergetic 236 Mev $\nu_\mu$ from charged-kaon decay-at-rest ($K^+ \rightarrow \mu^+ \nu_\mu$, with branching ratio of 64%) to achieve sensitivity to short-baseline muon-neutrino disappearance. The idea is to search for an $L/E$-dependent oscillation wave using fixed-$E$ neutrinos with minimal background and only modest detector requirements.

The KPIPE detector, relying on liquid scintillator and silicon photomultipliers (or PMTs), is designed to look for 236 MeV $\nu_\mu n \rightarrow \mu^{-} p$ interactions, which provide a unique double-flash coincidence due to the muon decay following the initial prompt event. Mapping these interactions as a function of distance along the detector pipe, with a nominal, no-oscillation expectation of a $1/r^2$ rate dependence, provides sensitivity to muon-flavor disappearance. Given a beam dump, decay-at-rest neutrino source, the beam-based $\nu_\mu$ background (from decay-in-flight mesons) to these signal events is expected to be completely sub-dominant, at the 1-2% level. While cosmics can be considered a concern for such a surface or near-surface detector, this background can be mitigated by typical accelerator duty factors of $\sim10^{-6}-10^{-5}$ combined with the short charged kaon lifetime (13 ns). The monoenergetic neutrino source, combined with low decay-in-flight background and small beam duty factor, means that the signal-to-background ratio is expected to be well over 50:1 in the scenarios considered. This large ratio means that the detector requirements, in particular the photocoverage, can be quite modest. In fact, a preliminary estimate at Ref. <cit.> predicts that the entire KPIPE detector would cost $5M.