1 Introduction

Right-handed neutrinos are a minimal and well-motivated extension of the Standard Model of Particle Physics (SM) [1]. Right-handed neutrinos, as opposed to the known left-handed neutrinos, would not interact in any SM interaction and are therefore called sterile neutrinos. The introduction of right-handed partners to the left-handed neutrinos provides a natural way to create neutrino masses [2]. No gauge symmetry of the SM forbids the introduction of a Majorana mass term of arbitrary scale for the right-handed neutrino. As a consequence, new neutrino-mass eigenstates arise, which are mostly sterile, but can have an admixture of the active SM neutrinos [3]. The size of the admixture is typically given by \(\sin ^2\theta \), where \(\theta \) refers to the active-to-sterile mixing angle. In the following the new mass eigenstates are referred to as sterile neutrinos.

Very light sterile neutrinos in the eV-mass range are motivated by long-standing anomalies in short-baseline-oscillation experiments [2, 4, 5]. While the so-called “Reactor Anomaly” [6] can be explained with renewed calculations of the reactor neutrino spectra and fluxes [7], a new anomaly was reported in a recent result from the BEST experiment [8]. Sterile neutrinos in the keV scale are viable candidates for dark matter [1, 9,10,11]. For very large masses (>GeV), sterile neutrinos could solve the puzzle of the lightness of active neutrinos via the see-saw mechanism and may shed light on the matter/anti-matter asymmetry of the universe [12,13,14].

A notable feature of sterile-neutrino dark matter is that it can act as effectively cold or warm dark matter depending on its production mechanism in the early universe [15]. This property can help mitigate tensions between predictions of purely cold dark-matter scenarios and observations of small-scale structures in the universe. The existence of sterile-neutrino dark matter is strongly bound by indirect searches and cosmological observations, which limit their mixing amplitude with active neutrinos to \(\sin ^2\theta < 10^{-6}{-}10^{-10}\) in a mass range of \((1{-}50) \, \textrm{keV}\) [1, 11, 16,17,18,19]. These limits can be model-dependent and could potentially be circumvented [20]. Current laboratory-based limits are orders of magnitude weaker [21,22,23,24,25].

2 Sterile neutrino searches with KATRIN

Sterile neutrinos can leave a signature in a \(\upbeta \)-decay spectrum [26,27,28]. If the sterile neutrino mass \(m_4\) is smaller than the endpoint energy \(E_0\) of the decay, the emission of a sterile neutrino along with the \(\upbeta \) electron is kinematically allowed.

Generally, the \(\upbeta \)-decay spectrum \(R_{\beta }(E)\) is a superposition of spectra corresponding to the different neutrino-mass eigenstates with masses \(m_i\) that contribute to the electron neutrino flavor eigenstate. Due to the tiny mass differences the three known light neutrino-mass eigenstates cannot be resolved with current experiments and are instead approximated by the squared effective electron neutrino mass \(m_\upnu ^2 = \sum _{i=1}^3 |U_{\upnu i}|^2 m_i^2\), where \(U_{\upnu i}\) denote elements of the Pontecorvo–Maki–Nakagawa–Sakata matrix.

In contrast, a heavy mass state \(m_\textrm{4}\), would lead to a distinct decay branch, with a maximal energy of \(E_\textrm{max}= E_0-m_\textrm{4}\) and an amplitude given by the mixing of the fourth mass eigenstate with the electron neutrino flavour, here parameterized as \(\sin ^2\theta \). As a result, the complete tritium \(\upbeta \)-decay spectrum

$$\begin{aligned} R_\upbeta (E) = \cos ^2\theta \, R_\upbeta (E, m_\upnu ^2) + \sin ^2\theta \, R_\upbeta (E, m_4^2) \end{aligned}$$
(1)

is a superposition of the effective active decay branch \(R_{\beta }(E, m_\upnu ^2)\) and the sterile branch \(R_{\beta }(E, m^2_4)\). The sterile neutrino signature thus appears as a kink-like feature and spectral distortion, as illustrated in Fig. 1.

Fig. 1
figure 1

Illustration of a keV-scale sterile-neutrino signature in the tritium \(\upbeta \)-decay spectrum \(R_\upbeta (E)\). The position of the kink-like signal is determined by the mass of the sterile neutrino \(m_4\) (here \(m_4=10\,\textrm{keV}\)) and the amplitude is governed by the mixing amplitude \(\sin ^2\theta \) (here \(\sin ^2\theta = 0.1\)). This value is unrealistically large, and was chosen for illustrative purpose only

Full size image

The Karlsruhe Tritium Neutrino experiment (KATRIN) [29] has one of the strongest tritium sources used for scientific research. The primary goal of the experiment is to probe the effective electron anti-neutrino mass with a sensitivity of \(<0.3\,\textrm{eV}\) at \(90\%\) confidence level after 1000 days of data taking. This is achieved by analyzing the shape of the tritium \(\upbeta \)-decay spectrum near the endpoint at \(E_0 = 18.6\,\hbox {keV}\), where the impact of the neutrino mass is maximal. Recently, KATRIN published its first sub-eV limit on the effective electron anti-neutrino mass of \(0.8\,\textrm{eV}\) (\(90\%\) CL) [30,31,32], based on the first two high-tritium-activity data-taking campaigns.

Several studies have shown that a KATRIN-like measurement also provides a promising sensitivity to eV- and keV-scale sterile neutrinos [27, 28, 33]. Based on the first two KATRIN measurement campaigns, improved limits could be set on eV-scale sterile neutrinos [34, 35].

While an eV-scale sterile neutrino leaves a signature within the standard measurement interval of KATRIN, which extends to about 100 eV below \(E_0\), the signature of a keV-scale sterile neutrino lies further away from the endpoint, outside of this interval. Consequently, a search for keV-scale sterile neutrinos requires an extension of the measurement interval which bears several challenges. One of them is the fact that the count rates deep in the spectrum exceed the level that can be resolved by the KATRIN focal-plane detector system [36, 37]. Nevertheless, at the cost of reduced statistics, it is possible to extend the measurement interval by reducing the source activity [33].

In 2018, the KATRIN beamline was operated for the first time with a small amount of tritium gas [38]. For safety reasons, the isotopic abundance of tritium in the deuterium carrier gas was set to only \(0.5\%\) in this commissioning campaign. The reduced tritium activity provided a unique opportunity to explore the spectrum in a wide energy range down to \(1.6\, \textrm{keV}\) below the endpoint. In this work we present the search for sterile neutrinos in the \(0.01{-}1.6\,\hbox {keV}\) mass range based on this 12-day-long series of commissioning measurements.

Fig. 2
figure 2

The experimental setup of the 70-m-long KATRIN beamline. Gaseous molecular tritium is inserted through capillaries at the center of the Windowless Gaseous Tritium Source (WGTS) (b). \(\upbeta \)-electrons created in the 10-m long WGTS are guided with a system of superconducting solenoids through the transport section (c) towards the spectrometer section. The pre-spectrometer (d) can pre-filter electrons and the main spectrometer (e) transmits only electrons above a sharp adjustable transmission edge. The 148-pixel focal plane detector (f) counts the transmitted electrons as a function of the main spectrometer’s transmission edge. Non-transmitted electrons are eventually absorbed in the rear wall of the rear-section (a) of the beam line

Full size image

3 The KATRIN working principle

The KATRIN experiment consists of a 70-m-long beamline (Fig. 2), combining a high-activity (up to \(10^{11}\) decays per second) gaseous molecular tritium source with a high-resolution (\(\nicefrac {\varDelta E}{E}=\mathscr {O}\)(1 eV)) spectrometer to obtain an integral \(\upbeta \)-decay spectrum [29].

The windowless gaseous tritium source (WGTS) consists of a \(10\,\textrm{m}\) long stainless-steel tube with a diameter of \(90\,\textrm{mm}\). Highly purified tritium gas is injected continuously at the center of the WGTS and diffuses to the up- and downstream end of the source tube where it is pumped out and fed back to the tritium loop system that is integrated in the infrastructure of the Tritium Laboratory Karlsruhe [39].

The source and spectrometer sections of the KATRIN beamline are connected by the so-called transport section. Here, differential [40] and cryogenic [41] pumping systems reduce the tritium flow by more than 14 orders of magnitude, while the electrons are guided adiabatically to the spectrometers by a system of superconducting magnets.

The high-resolution main spectrometer selects the electrons according to their energy, by applying an Electrostatic filter (E) combined with a Magnetic Adiabatic Collimation (MAC). The MAC-E filter only transmits electrons with a longitudinal kinetic energy (kinetic energy component associated with the motion parallel to the magnetic field lines) larger than the retarding energy qU, where U is the precisely adjustable voltage of the spectrometer [42] and q refers to the electron charge. A magnetic field, which decreases by approximately 4 orders of magnitude from the ends to the center of the spectrometer, transforms the total kinetic energy of the electrons into longitudinal energy. The MAC-E filter technology thus combines a large angle acceptance of \(51^{\circ }\) with an sharp filter width of 2.8 eV. The main spectrometer is preceded by a smaller pre-spectrometer, which also works according to the MAC-E filter principle, and transmits only electrons above \(10\,\hbox {keV}\), to reduce the flux of electrons into the main spectrometer.

Electrons that overcome the retarding potential in the main spectrometer are counted at the focal-plane detector (FPD) [36, 37]. The FPD is a monolithic silicon array, radially and azimuthally segmented in 148 pixels. By measuring the count rate at different retarding energies, the integral \(\upbeta \)-decay spectrum is obtained. In order to increase the signal-to-background ratio, the transmitted electrons are accelerated by a post-acceleration electrode (PAE) with an electrostatic potential of \(U_\textrm{PAE}=10 \, \textrm{kV}\) before impinging on the detector surface.

4 The first tritium campaign

The First Tritium (FT) campaign, which inaugurated the KATRIN experiment, was a commissioning campaign to demonstrate the stable operation of the integral system and test different analysis strategies. A technical description of the measurement campaign and the results with respect to stability and analysis techniques can be found in [38].

4.1 Tritium source operation

During the FT campaign, the WGTS was operated at a column density (gas density \(\rho \) integrated over the length d of the source) of \(\rho d = 4.46\times 10^{17} \, \textrm{molecules}/\textrm{cm}^2\), with a reduced tritium activity of \(500\,\textrm{MBq}\), which corresponds to 0.5% of the activity used for neutrino-mass measurements. This activity limitation was achieved by mixing traces of tritium (in the form of DT) with pure deuterium (\(\text {D}_2\))  [43, 44]. This gas mixture was circulated through the WGTS via the main tritium loop [45]. At all times, the gas composition was monitored by a Laser Raman spectroscopy (LARA) system [46, 47] and by the Forward Beam Monitor (FBM) [48]. In the FT experimental configuration, the downstream end of the KATRIN beam line was terminated by a stainless-steel gate valve, rather than the gold-plated rear wall, which was later in place for the neutrino-mass measurements.

Fig. 3
figure 3

The top panel shows simulated data points corresponding to a \(\upbeta \)-decay spectrum including a fourth, sterile mass eigenstate with \(m_4=400\,\textrm{eV}\) and a mixing amplitude of \(\sin ^2 \theta = 0.01\). The dashed orange line displays the spectrum prediction with these sterile neutrino parameters, and the solid green line shows a spectrum without sterile neutrino. In order to illustrate the sterile neutrino signature, the middle panel displays the ratio of the spectrum with sterile neutrino to a spectrum without sterile neutrino (scaled by \(\sin ^2 \theta = 0.99\), to account for the difference in total normalization). At an energy of 400 eV below \(E_{0}\) the additional sterile-neutrino branch kicks in and distorts the overall spectrum for energies \(E < E_{0}-400\,\textrm{eV}\). The bottom part of the figure shows the accumulated measurement time distribution of all analyzed tritium scans

Full size image

4.2 Spectrometer operation

KATRIN obtains the integral \(\upbeta \)-decay spectrum in so-called scans, i.e. by sequentially applying different retarding energies \(qU_i\) to the main spectrometer and counting the number of transmitted \(\upbeta \)-electrons \(N(qU_i)\) with the focal plane detector. During the FT campaign, the spectrum was measured at 26 different retarding potentials in the range of \(E_0 - 1600\,\hbox {eV} \le qU_i \le E_0 + 30\,\hbox {eV}\). Figure 3 shows the measurement-time distribution during FT data taking. The sequence of applied retarding potentials is either increasing (up scans) or decreasing (down scans). Applying up scans and down scans in an alternating fashion optimizes the averaging of possible drifts of slow-control parameters and also minimizes the time for setting the retarding potentials.

The FT measurement entails 122 scans, each of which with a duration of 1–3 h, leading to a total measurement time of 168 h. The \(\upbeta \)-decay spectrum obtained in each individual scan, was analyzed separately to test the stability of the system. The obtained effective endpoint of each spectrum shows an excellent stability, consistent with purely statistical fluctuations [38]. Moreover, we demonstrate that the inferred endpoint does not depend on the scanning direction, the column density, and the analysis window [38].

5 Spectrum calculation

The expected integral \(\upbeta \)-decay spectrum is composed of two main parts: (1) the theoretical differential \(\upbeta \)-electron spectrum \(R_\upbeta (E)\) and (2) the experimental response function \(f_\textrm{calc}(E, qU_i)\). The total calculated rate \(R_\textrm{calc}(qU_i)\) at a given retarding energy \(qU_i\) is given by

$$\begin{aligned} R_\textrm{calc}(qU_i) = A_{\textrm{s}} N_{\textrm{T}} \int _{qU_i}^{E_0}R_\upbeta (E) f_\textrm{calc}(E, qU_i) \ \textrm{d}E + R_{\textrm{bg}}, \end{aligned}$$
(2)

where \(A_{\textrm{s}}\cdot N_{\textrm{T}}\) is the signal normalization, which includes the number of tritium atoms in the source, the maximum acceptance angle of the MAC-E filter and the detection efficiency. \(R_{\textrm{bg}}\) denotes the retarding-potential-independent background rate [49]. Both \(A_{\textrm{s}}\) and \(R_{\textrm{bg}}\) are treated as free nuisance parameters of the fit.

5.1 Differential \(\upbeta \)-decay spectrum

As described above, the \(\upbeta \)-decay spectrum \(R_\upbeta (E)\) is given as a superposition of the active \(R_\upbeta (E, m_\upnu ^2)\) and sterile \(R_\upbeta (E, m_4^2)\) decay branches. Each one of them is described by Fermi’s theory

$$\begin{aligned} R_\upbeta (E, m_\upnu )= & {} C \cdot F(Z',E) \cdot (E+m_\textrm{e}) \cdot p \cdot \sum _i P_i\nonumber \\{} & {} \cdot (E_0-E-E_i)^2 \cdot \sqrt{1-\biggl (\frac{m_\upnu }{E_0-E-E_i}\biggr )^2}, \nonumber \\ \end{aligned}$$
(3)

where \(C= \frac{G_F^2}{2\pi ^3}\cos ^2\varTheta _C |M_{\text {nucl}}|^2\) with \(G_{\textrm{F}}\) denoting the Fermi constant, \(\varTheta _{\textrm{C}}\) the Cabibbo angle, and \(M_{\text {nucl}}\) the energy-independent nuclear matrix element. The \(F(E, Z')\) represents the Fermi function with \(Z'=2\) for the atomic number of helium, the daughter nucleus in this decay. E, p, and \(m_{\text {e}}\) denote the kinetic energy, momentum, and mass of the \(\upbeta \)-electron, respectively.

After the \(\upbeta \)-decay of tritium in a DT molecule, the daughter molecule \(\mathrm {^3HeD^+}\) can end up in an electronic ground state or excited state, each of which is broadened by rotational and vibrational excitations of the molecule [50]. As a consequence, this excitation energy \(E_i\) reduces the available kinetic energy for the electron. Thus the differential \(\upbeta \)-electron spectrum is a superposition of spectra, corresponding to all possible final states, weighted by the probability \(P_i\) for decaying into a certain final state i. For this analysis, we use the latest calculation of Saenz et al. for the isotopologue DT [51].

The molecular final-state distribution depends slightly on the \(\upbeta \)-decay energy. Mainly, the mean and width of the ground-state distribution depends on the recoil energy of the daughter molecule, which in turn depends on the \(\upbeta \)-decay energy [50]. By taking into account this energy dependence in the theoretical calculation of the integral \(\upbeta \)-decay spectrum, we find that \(R_\textrm{calc}(qU_i)\) is altered by less than \(0.007\%\) for all retarding energies. Hence, we neglect the energy dependence of the final-state distribution in this analysis. Doppler broadening due to the thermal motion of tritium molecules in the source, which is operated at 30 K, is emulated as a broadening of the molecular final-state distribution [52].

5.2 Response function

The experimental response function

$$\begin{aligned} f_\textrm{calc}(E, qU_i)= & {} \int _0^E T(E-\varepsilon , qU_i) \left( P_0\,\delta (\varepsilon ) + P_1\,f(\varepsilon ) \right. \nonumber \\{} & {} \left. + P_2\,(f\otimes f)(\varepsilon ) +\cdots \right) \, \textrm{d}\varepsilon , \end{aligned}$$
(4)

is the probability of an electron with a starting energy E to reach the detector. It combines the transmission function T of the MAC-E filter and the electron’s energy losses \(\varepsilon \) in the source. The transmission function T reflects the resolution of the main spectrometer and is governed by the magnetic fields at the starting position of the electron, the maximum field in the beamline, and the magnetic field in the spectrometer’s analyzing plane. Energy losses due to inelastic scattering with the deuterium molecules in the source are described by the product of the s-fold scattering probabilities \(P_s\) and the energy-loss function \(f(\varepsilon )\) convolved \((s-1)\) times with itself. We consider an energy-dependent cross-section, but treat the energy-loss function \(f(\varepsilon )\) as energy independent. Here we use an energy-loss function measured in situ for deuterium [53]. Synchrotron energy losses of \(\upbeta \)-electrons in the high magnetic field in the source and transport section are included as a correction to the transmission function. Furthermore, the response function is slightly modified due to the dependence of the path length (and therefore effective column density) on the pitch angle (angle between the electron’s momentum vector and the magnetic field line) of the \(\upbeta \)-electrons. A detailed description of the theoretical spectrum calculation can be found in  [52] and  [31].

5.3 Wide-interval corrections

Beyond the tritium spectrum calculation, described above, we investigate specific effects relevant at energies further away from the endpoint, outside the nominal KATRIN analysis window.

5.3.1 Detection efficiency

The total detection efficiency is of minor relevance as it only affects the normalization of the measured spectrum and not its shape. In contrast, a retarding-potential-dependent detection efficiency alters the shape of the integral spectrum. Figure 4a displays the retarding-potential dependence of the detection efficiency. It includes the following effects:

Region-of-interest coverage

In order to count the events at a given retarding potential, the measured rate at the focal-plane detector is integrated in a wide and asymmetric region of interest (ROI) of \(14\,\hbox {keV} \le E + qU_{\text {PAE}} \le 32\,\hbox {keV}\), where E is the \(\upbeta \)-electron energy and \(U_{\text {PAE}}=10\,\hbox {keV}\) is the post-acceleration voltage. This ROI is chosen to account for the moderate energy resolution of about \(3\,\hbox {keV}\) (full-width-half-maximum) and the low-energy tail of the spectrum due to the energy loss of electrons in the dead layer and backscattering from the detector surface [37]. The same ROI is used for each retarding-potential setting. As the mean of the electron peak shifts with the retarding potential, some electrons move out of the fixed ROI, which effectively changes the detection efficiency. This change of detection efficiency is experimentally determined based on reference measurements, and is corrected accordingly [54]. For this effect, we interpret the variation of the correction between all detector pixels used in the analysis as the uncertainty. Assuming a detection efficiency of \(\varepsilon _{\text {roi}}=1\) at \(E_0\), we find a relative detection efficiency at \(1\,\textrm{keV}\) below \(E_0\) of \(\varepsilon _{\text {roi}}=0.99911 \pm 0.00036\).

Pile-up

As the counting rate at the focal-plane detector depends on the retarding potential, so does the probability of pile-up. Most pile-up events occur outside the ROI, thereby effectively changing the detector efficiency [55]. We estimate the detection efficiency \(\varepsilon _{\text {pu}}\) with a two-fold random coincident model, according to

$$\begin{aligned} \varepsilon _{\text {pu}}(R)= \left( 1-\frac{\alpha }{2}\right) \cdot e^{-2WR}+\frac{\alpha }{2}, \end{aligned}$$
(5)

where R is the Poissonian-distributed signal rate, \(1-\frac{\alpha }{2} = 0.79\pm 0.02\) denotes the pile-up event rejection ratio, and \(\hbox {W} = 1.826(0.026)\,\upmu \hbox {s}\) denotes the effective window length of the trapezoidal energy filter used to determine the energy of each event [37]. The uncertainty of this correction is determined by the uncertainty of these model parameters, listed in Table 1. At \(1\,\textrm{keV}\) below \(E_0\), pile-up reduces the detector efficiency to \(\varepsilon _{\text {pu}} = 0.99952 \pm 0.00001\).

Fig. 4
figure 4

a Multiplicative spectral corrections due to retarding-potential-dependent detection efficiency (green), magnetic trapping of \(\upbeta \)-electrons in the WGTS (blue), and backscattering of \(\upbeta \)-electrons on the rear wall of the WGTS (red). b Relative Poisson statistical uncertainty of the spectral data points (gray), relative statistical uncertainty arising from the deuterium-tritium (DT) source activity fluctuations (light blue), relative spectral uncertainties arising from the three corrections displayed in a (blue, red, green). Note that the latter three uncertainties are correlated between the spectral data points

Full size image
Table 1 Summary of systematic uncertainties. We list the uncertainties of the input parameters used to construct the covariance matrices. Note, that the correlations of the parameters describing the energy loss function is not shown in the table. In most cases the uncertainty corresponds to a 1-\(\sigma \) Gaussian uncertainty. An exception is the case of the detection efficiency correction due to the region-of-interest-coverage (see Sect. 5.3), here we estimate the uncertainty from the variation of the detection efficiency between different pixels. Also, for the case of the rear-wall backscattering the uncertainty does not correspond to a 1-\(\sigma \)-Gaussian uncertainty, but instead we consider the difference of the correction between two GEANT-4 physics libraries as a measure of the uncertainty, see main text for details
Full size table

Backscattering

A significant fraction of about \(20\%\) of all electrons impinging on the detector surface are backscattered [56]. For low retarding potentials and small energy depositions in the detector, these backscattered electrons have a chance of getting lost by overcoming the retarding potential a second time. The lower the retarding potential, the higher is the probability to lose electrons this way, effectively changing the detection efficiency [55]. We estimate this effect by Monte Carlo simulations with the KATRIN-specific simulation packages Kassiopeia [57] and Kess [56] which is optimized for the tracking of keV-scale electrons in electro-magnetic fields and in silicon, respectively.

We estimate the uncertainty of this correction by varying the input parameters of the simulation according to a Gaussian distribution, which 1-\(\sigma \) width is given by the uncertainties of the input parameters. These are the magnetic field at the position of the detector \(B_\textrm{det}\) and at the pinch magnet \(B_\textrm{pch}\), which is situated at the exit of the main spectrometer. The uncertainties of the magnetic fields are estimated via comparisons of measurements and simulations [38] and are quoted in Table 1. Moreover, as the Si-crystal lattice orientation relative to the electron’s incident angle is not precisely known, we allow for an uncertainty of the amplitude of the elastic backscattering peak. We conservatively vary the amplitude obtained by Monte-Carlo simulations by \(+50\%\), emulating the two extreme cases of anomalous transmission and absorption [58]. At \(1\,\textrm{keV}\) below \(E_0\), backscattering reduces the detector efficiency to \(\varepsilon _{\text {bs}} = 0.99893 \pm 0.00027\).

5.3.2 Rear-wall backscattering

Another effect which is negligible in the case of an endpoint analysis is the detection of \(\upbeta \)-electrons which are backscattered at the rear wall of the beamline and still reach the focal-plane detector. During the FT measurement campaign a stainless-steel gate valve terminated the beamline. In the backscattering process, the electrons lose some amount of energy, which typically forbids them to be transmitted through the main spectrometer. However, for low retarding potentials, there is a non-negligible probability for this transmission to occur [33]. The backscattering of tritium \(\upbeta \)-decay electrons from the stainless-steal plate was simulated with GEANT4, providing the backscattering probability as well as the energy and angle distribution of backscattered electrons. The corresponding correction to the integral \(\upbeta \)-decay spectrum is depicted in Fig. 4a.

As for the case of the detector backscattering, we estimate the uncertainty of this correction by varying in the simulation the magnetic fields at the rear wall \(B_\textrm{rw}\) and in the source section \(B_\textrm{s}\) according to Gaussian distributions, which 1-\(\sigma \) widths are given by the magnetic field uncertainties listed in Table 1. In addition, we estimate a theoretical uncertainty arising from the GEANT4 simulation, by computing the correction with different physics packages (i.e. the emlivermore and emstandardSS packages) and interpreting the difference as a measure of the uncertainty. This uncertainty is propagated by randomly drawing from a Bernoulli distribution. At \(1\,\textrm{keV}\) below \(E_0\), we find a multiplicative correction to the observed rate of \(\varepsilon _{\text {rw}} = 1.00097 \pm 0.00096\).

5.3.3 Magnetic trapping

The source beam line exhibits small local magnetic field minima, arising from the small gaps between adjacent superconducting coil units. Electrons starting with a pitch angle larger than a certain threshold in such local magnetic field minima can be magnetically trapped. Frequent elastic and inelastic scattering change their angle and they eventually escape from the trap with reduced energy. If the retarding potential of the spectrometer is low enough, these electrons have a chance to reach the detector [33]. Based on a Monte Carlo simulation with Kassiopeia, we calculate the corresponding correction to the integral \(\upbeta \)-decay spectrum, as displayed in Fig. 4a.

As described before, we obtain the uncertainty on the correction by varying the relevant simulation input parameters, namely the source \(B_\textrm{s}\) and pinch \(B_\textrm{pch}\) magnetic field, the gas density in the source \(\rho \), and the parameters of the energy loss function, according to a Gaussian distribution. The corresponding uncertainties listed in Table 1. At \(1\,\textrm{keV}\) below \(E_0\), we find a multiplicative correction to the observed rate of \(\varepsilon _{\text {mt}} = 1.00510 \pm 0.00017\).

5.3.4 Non-adiabaticity

At low retarding potentials of the MAC-E filter, some electrons have a comparatively high surplus energy. This is of concern, since the magnetic guiding field drops from about \(5\,\hbox {T}\), at the entrance to about \(6\times 10^{-4}\) T in the center of the spectrometer. If an electron experiences an excessive change of the magnetic field within one cyclotron circle, it exhibits non-adiabatic motion. The non-adiabatic motion causes a chaotic change of the pitch angle and hence a possible magnetic reflection at the exit of the spectrometer. Eventually this can lead to a reduction of the number of transmitted electrons [33]. A full Monte Carlo simulation with Kassiopeia shows that in the realistic magnetic field settings of the FT campaign, non-adiabatic effects can indeed occur at more than \(1\,\textrm{keV}\) below the endpoint. However, averaged over all radii in the spectrometer, this effect leads only to a small reduction of the rate of less than \(0.01\%\) for all retarding potentials used in this measurement and can thus be neglected.

6 Data selection and combination

The full FT data set is composed of a large number of individual \(\upbeta \)-decay spectra: (1) As mentioned above, the integral tritium spectrum is recorded in 122 scans to accommodate temporal changes of slow-control parameters, such as the source activity. (2) Each of the 148 pixels of the focal-plane detector measures a statistically independent tritium \(\upbeta \)-decay spectrum, to take into account radial and azimuthal variations of the electric and magnetic fields in the analyzing plane.

For this analysis, we combine a selection of 82 “golden” scans, excluding scans that were performed at different experimental settings, such as at a different column density or with different HV set points. The combination is performed by adding the counts recorded at each retarding potential set point, called scan step, to construct a high-statistics single spectrum with \(n_{\text {scan-step}}=26\) data points.

Equivalently, for all “golden” scans, we combine 119 “golden” pixels, excluding pixels which do not record the full flux of electrons due to misalignment. The pixels are combined in a single effective pixel, by adding all counts and assuming an average response function for the entire detector. Simulations have shown that these assumptions lead to a negligible error on the fitted parameters [38]. A full description of the data quality criteria can be found in [38].

7 Method of the fit and confidence level setting

The calculated model spectrum \(\vec {R}_\textrm{calc}\) is fit to the data \(\vec {R}_\textrm{data}\) by minimizing

$$\begin{aligned} \chi ^2(\alpha _i) = (\vec {R}_\textrm{calc}(\alpha _i)-\vec {R}_\textrm{data})^{\textsf{T}} C^{-1} (\vec {R}_\textrm{calc}(\alpha _i)-\vec {R}_\textrm{data}), \end{aligned}$$
(6)

with respect to \(\alpha _i\), which includes the parameters of interest and the nuisance parameters. The unconstrained nuisance parameters in this analysis are the signal normalization, the effective endpoint of the spectrum, and an overall background rate.

C is the covariance matrix, which contains both statistical and systematic uncertainties. The statistical uncertainty appears on the diagonal of the matrix and is determined by the total counts at each retarding potential set point. The covariance matrix describing systematic uncertainties is constructed by computing the spectrum prediction about \(10^4\) times while varying the systematic parameters at each execution. Most systematic parameters are varied according to a Gaussian distribution, which width corresponds to the 1-\(\sigma \) uncertainty of the respective parameter. One exception is the error propagation of the rear-wall backscattering (see Sect. 5.3), where we vary the prediction of two physics libraries of GEANT-4, according to a Bernoulli distribution. A second exception is the detector efficiency correction due to the ROI cut (see Sect. 5.3), where we vary the correction according to the variation of the detection efficiency of the detector pixels. As a result of this procedure, we obtain the covariance matrix C, i.e. both the variance and the covariance of all spectral data points \(\vec {R}_\textrm{calc}\), arising from a specific systematic effect. The matrix is then included in the \(\chi ^2\)-function, as can be seen in equation (6).

In order to search for the signal of a sterile neutrino, we follow the standard Neyman procedure for constructing the frequentist confidence intervals [59]. To this end, the fit is repeated on a fine grid of fixed tuples of the sterile neutrino mass \(m_4\) and mixing amplitude \(\sin ^2\theta \). At each grid point the value of \(\chi ^2\) is computed, by performing the fit with \(m_4\) and \(\sin ^2\theta \) fixed to the value of the grid point. According to Wilks’ theorem [60], the \(95\%\) confidence level (C.L.) exclusion limit is constructed by determining the \(\varDelta \chi ^2=\chi ^{2}(m_4, \sin ^2\theta )-\chi ^2_{\textrm{min}} < 5.99\) contour, where \(\chi ^2_{\textrm{min}}\) corresponds to the global best fit. The applicability of Wilks’ theorem for a sterile neutrino search with KATRIN was tested with Monte-Carlo simulations for the null hypothesis and positive signals [61].

8 Statistical sensitivity and impact of systematic uncertainties

As a first step of the analysis we perform a sensitivity study based on a Monte-Carlo copy of the FT data, where we assume no sterile neutrino signature. This allows us to assess the statistical sensitivity and the impact of individual systematic effects. We consider the standard KATRIN systematic uncertainties, described in detail in [38], and uncertainties arising from the wide-range corrections, described in Sect. 5.3. All systematic uncertainties are listed in Table 1.

Figure 5 displays the statistical sensitivity at the \(95\%\) C.L. and the sensitivities when including individual and all systematic uncertainties. The \(95\%\) C.L. statistical sensitivity reaches down to a value of \(\sin ^2\theta < 5\times 10^{-4}\) at \(m_4 = 1000\,\hbox {eV}\). Including all systematic uncertainties the best sensitivity is reduced to \(\sin ^2\theta < 2\times 10^{-3}\).

The study shows that the most limiting uncertainty arises from the uncertainty of the tritium activity in each scan step. This uncertainty corresponds to an additional statistical error and thus can mimic a sterile neutrino signature, which explains the large impact on the sensitivity. With a relative magnitude of \(5\times 10^{-4}\) it dominates over the Poisson error (arising from the counting statistics) for retarding energies of \(qU < E_0 - 400\, \textrm{eV}\), see Fig. 4b. This uncertainty is dominated by the statistical uncertainty of the LARA and FBM systems. Accordingly, it will be reduced when operating at higher activity and with longer measurement time.

The relative uncertainty arising from electrons that scatter off the rear wall and then reach the detector amounts to \(1\times 10^{-3}\) at \(qU < E_0-1000\, \textrm{eV}\), see Fig. 4b. Even though the magnitude of this uncertainty is larger than the one from activity fluctuations, its impact on the sterile-neutrino sensitivity is smaller. This is due to the fact that these uncertainties are strongly correlated between the different scan steps, thus preventing this correction from mimicking a kink-like sterile-neutrino signature [27]. Similarly, the total uncertainty arising from the detector efficiency is of similar size to the activity fluctuations. However, this uncertainty is strongly correlated between the data points, and thus its impact on the sterile-neutrino sensitivity is mitigated.

The next largest effect arises from scattering of \(\upbeta \) electrons in the source section. Relevant parameters to describe energy losses due to scattering are the column density, the cross section, and the parameterized energy-loss function [53], which uncertainties are given in Table 1. The impact on the sterile-neutrino sensitivity is relatively large in this analysis, since the uncertainty on the column density during the FT campaign was rather high (3\(\%\)). This is due to the fact that a calibration of the absolute column density was not available at that time. In subsequent KATRIN campaigns, the product of column density and cross-section is directly measured with an electron gun, reaching a precision of \(<1\%\) [31].

The uncertainties of the magnetic fields along the KATRIN beamline lead to a minor but still visible impact on the sterile-neutrino sensitivity. With an uncertainty of \(2.5\%\), the source magnetic field has the largest impact on the sterile neutrino sensitivity. With the help of new calibration methods, the magnetic-field uncertainties were reduced by up to one order of magnitude in later KATRIN measurement campaigns. An impact of similar size can be attributed to the uncertainties of the molecular final state distribution. Here, we assume an uncertainty on the order of \(1\%\) on the probability to decay into the electronic ground-state and the broadening due to rotational and vibrational states, see Table 1. This uncertainty is larger than what is stated by experts on the theoretical calculations [62] and it was chosen as a conservative estimation [31].

Interestingly, an uncertainty on the retarding-potential dependence of the background does not lead to a visible effect on the sterile-neutrino sensitivity. The retarding-potential dependence of the background is modeled by a so-called background slope, which is constrained by dedicated background measurements to less than 5 mcps/keV. Similarly, the uncertainty arising from initially trapped electrons in local magnetic field minima the source section leads to a negligible effect on the sterile-neutrino sensitivity in this analysis.

Fig. 5
figure 5

\(95\%\) confidence level (C.L.) sensitivity to sterile neutrinos based on a Monte-Carlo copy of the first tritium data set. The statistical-only sensitivity is displayed by the solid blue line. The sensitivity including both statistical and systematic uncertainties is shown in black. The dashed lines show the impact of the statistical and individual systematic uncertainties: DT activity fluctuations (dark-blue-loosely-dashed), rear-wall backscattering (dark-red-dotted), all detector efficiency corrections (dark-red-dashed), scattering in the source (green-dash-dotted), all magnetic fields (red-densely-dash-dotted), molecular final state distribution (blue-loosely-dotted), background slope (yellow-densely-dash-dot-dotted), magnetic traps (magenta-dash-dot-dotted)

Full size image

9 Best fit and exclusion limit

Fig. 6
figure 6

The best fit of all 82 spectra combined (by adding the counts at each retarding energy) with the normalized residuals expressed in standard deviation \(\sigma \). The best fit is found for \(m_4 = 71.2\,\hbox {eV}\) and \(\sin ^2\theta = 0.017\) with a goodness-of-fit of \(\chi ^2/\textrm{ndof}=14.79/21\), and a corresponding p-value of 0.83

Full size image

The statistics of the full data set amount to \(1.2 \times 10^{9}\) \(\upbeta \)-electrons. The corresponding spectrum with the best fit including a sterile neutrino and all systematic uncertainties, shows an excellent agreement of the model with the data, as shown in Fig. 6. We find \(\chi ^2/\textrm{ndof}=14.79/21\), and a corresponding p-value of 0.83.

Following the procedure outlined in Sect. 7, we scan the parameter space (\(m_4\), \(\sin ^2\theta \)) and determine the minimal \(\chi ^2\)-value at each grid point. The best fit is found for \(m_4 = 71.2\,\hbox {eV}\) and \(\sin ^2\theta = 0.017\). This best fit is compatible with the null hypothesis, as the \(\chi ^2\)-difference to the null hypothesis is \(\varDelta \chi ^2 = 5.13\), corresponding to a significance of \(92.3\%\) (for two degrees of freedom). Based on this result, we determine the \(95\%\) C.L. exclusion limit, as shown in Fig. 7. For a mass of \(m_4 = 300\,\hbox {eV}\) we find the strongest exclusion limit of \(\sin ^2\theta < 5\times 10^{-4}\) at \(95\%\) CL.

In Fig. 7 we also display the result of the eV-scale sterile neutrino search based on the first two high-tritium-activity data taking campaigns of KATRIN [34, 35], illustrating the complementary of the eV- and keV-scale sterile neutrino searches.

Fig. 7
figure 7

\(95\%\) C.L. exclusion limit obtained based on the first tritium data set of KATRIN (blue). The sensitivity obtained via MC simulation is shown by the dashed black line. We improve the current laboratory limits [21,22,23,24,25] (colored shaded areas) on the active-to-sterile mixing amplitude in a mass range of \(0.1\,\textrm{keV}< m_4 < 1.0\,\textrm{keV}\) by up to an order of magnitude. The orange line shows the \(95\%\) C.L. exclusion limit obtained from the search for eV-scale sterile neutrinos based on the first two KATRIN data taking campaigns [34, 35]

Full size image

Finally, we compare our achieved exclusion limit with previous laboratory-based sterile-neutrino searches [21,22,23,24,25]. The Troitsk nu-mass experiment provides the leading limit for sterile-neutrino masses of \(m_4 < 0.1\,\textrm{keV}\), based on a re-analysis of their neutrino-mass data [24]. An upgrade of the experiment [63] allowed the extension of the measurement interval, setting a new limit for sterile-neutrino masses in the range of \(0.1 \, \textrm{keV}< m_4 < 2\,\textrm{keV}\) [64]. With the analysis presented in this work, we can improve this limit in a mass range of \(0.1\,\textrm{keV}< m_4 < 1.0\,\textrm{keV}\).

10 Conclusion and outlook

In this work we have performed a search for keV-scale sterile neutrinos with a mass of up to \(1.6\,\textrm{keV}\), based on the first commissioning run of the KATRIN experiment. The analysis includes a careful study of possible systematic uncertainties that occur when extending the nominal KATRIN measurement interval, which is restricted to a region close to the tritium endpoint.

As a result we exclude an active-sterile mixing amplitude of \(\sin ^2\theta < 5\times 10^{-4}\) for a sterile neutrino mass of \(m_4 = 300\, \textrm{eV}\). With this work, we improve currently leading laboratory-based bounds in a mass range of \(0.1\,\textrm{keV}< m_4 < 1.0\,\textrm{keV}\). This result establishes a major milestone for the keV-scale sterile-neutrino program of KATRIN and sets the groundwork for future high-statistics measurements.

Currently, a new detector system for KATRIN, the TRISTAN detector, is being developed, which is designed to allow KATRIN to extend the measurement interval to several keV below the endpoint and further improve the laboratory-based sensitivity to keV-scale sterile neutrinos [65]. This technique will exploit a combination of differential and integral spectral measurements to exclude large classes of systematic effects [65].