Direct neutrino-mass measurement with sub-electronvolt sensitivity

Since the discovery of neutrino oscillations, we know that neutrinos have non-zero mass. However, the absolute neutrino-mass scale remains unknown. Here we report the upper limits on effective electron anti-neutrino mass, mν, from the second physics run of the Karlsruhe Tritium Neutrino experiment. In this experiment, mν is probed via a high-precision measurement of the tritium β-decay spectrum close to its endpoint. This method is independent of any cosmological model and does not rely on assumptions whether the neutrino is a Dirac or Majorana particle. By increasing the source activity and reducing the background with respect to the first physics campaign, we reached a sensitivity on mν of 0.7 eV c–2 at a 90% confidence level (CL). The best fit to the spectral data yields mν2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mbox{}}}{m}_{\nu }^{2}{{\mbox{}}}$$\end{document} = (0.26 ± 0.34) eV2 c–4, resulting in an upper limit of mν < 0.9 eV c–2 at 90% CL. By combining this result with the first neutrino-mass campaign, we find an upper limit of mν < 0.8 eV c–2 at 90% CL. In its second measurement campaign, the Karlsruhe Tritium Neutrino experiment achieved a sub-electronvolt sensitivity on the effective electron anti-neutrino mass.

T he discovery of neutrino flavour oscillations 1,2 proves that neutrinos must have a mass, unlike originally assumed in the standard model of particle physics. Neutrino oscillation experiments have shown that the weakly interacting neutrino flavour eigenstates ν f , where f ∈ {e, μ, τ} for electron, muon and tau-neutrino, are admixtures of the three neutrino-mass eigenstates ν i with mass eigenvalues m i i ∈ {1, 2, 3}. Although neutrino oscillation experiments can probe the differences of squared neutrino-mass eigenvalues Δm 2 ij , the absolute neutrino-mass scale remains one of the most pressing open questions in the fields of nuclear, particle and astroparticle physics today. In this paper, we report a measurement of the effective electron anti-neutrino mass defined as m 2 ν = ∑ i |U ei | 2 m 2 i , where U ei are elements of the Pontecorvo-Maki-Nakagawa-Sakata matrix that describes the mixing of neutrino states.
The neutrino masses are at least five orders of magnitude smaller than the mass of any other fermion of the standard model, which may point to a different underlying mass-creation mechanism 3 . The determination of the neutrino mass would, thus, shed light on the fundamental open question of the origin of particle masses. Despite the smallness of their masses, neutrinos play a crucial role in the evolution of large-scale structures of our cosmos due to their high abundance in the Universe 4,5 . A direct measurement of the neutrino mass could, hence, provide a key input to cosmological structure formation models. In this respect, cosmological observations themselves provide a stringent limit on the sum of neutrino masses of ∑m i < 0.12 eV (95% confidence level (CL)) 6,7 (here we use the convention c = 1 for the speed of light). However, these limits strongly rely on the underlying cosmological assumptions 8,9 . An independent measurement of neutrino mass could help in breaking the parameter degeneracies of cosmological models. A powerful way to probe this neutrino property in the laboratory is via a search for neutrinoless double-beta (double-β) decay. In contrast to m ν , the effective mass in double-β decay is given by m ββ = i U 2 ei m i . This neutrino-mass interpretation is only valid under the assumption that neutrinos are their own anti-particle (Majorana particle) and that light neutrinos mediate the decay. The current most stringent limits derived from different isotopes are m ββ < 79−180 meV ( 76 Ge) (ref. 10 ), m ββ < 75−350 meV ( 130 Te) (ref. 11 ) and m ββ < 61−165 meV ( 136 Xe) (ref. 12 ), and the spread is related to uncertainties in the model-dependent nuclear matrix element calculation.
The most direct way to assess the neutrino mass is via the kinematics of single-β decays or electron capture processes. This method is independent of any cosmological model and of the mass nature of the neutrino, that is, it may be a lepton of the Majorana or Dirac type. The neutrino masses m i lead to a reduction in the maximum observed energy of the decay and a small spectral-shape distortion close to the kinematic endpoint of the β-decay spectrum. In the quasi-degenerate mass regime, where m i > 0.2 eV, the mass splittings are negligible with respect to masses m i , and the observable value can be approximated as m 2 ν = ∑ i |U ei | 2 m 2 i (refs. 13,14 ). The Karlsruhe Tritium Neutrino (KATRIN) experiment 15,16 exploits the single-β decay of molecular tritium as T 2 → 3 HeT + + e − +νe (1) and currently provides the best neutrino-mass sensitivity in the field of direct neutrino-mass measurements with its first published limit of m ν < 1.1 eV (90% CL) 17,18 . KATRIN is designed to determine the neutrino mass with a sensitivity of close to 0.2 eV (90% CL) in a total measurement time of 1,000 days (ref. 15 ). Another class of experiments is based on the electron capture of 163 Ho, where the decay energy is measured with micro-calorimeters 19,20 . Note that electron capture experiments based on 163 Ho measure the mass of the neutrino ν and β − experiments based on tritium measure that of the anti-neutrino ν. New ideas exist to extend the sensitivity of tritium-based neutrino-mass experiments beyond the KATRIN design sensitivity by new technologies, such as cyclotron radiation emission spectroscopy and the development of atomic tritium sources 21,22 .
In this work, we present the second neutrino-mass result of KATRIN, reaching an unprecedented sub-electronvolt sensitivity and limit in m ν from a direct measurement.

The KATRIN experiment
The design requirements to detect the small signature of a neutrino mass in the last few electron-volts of the β-decay spectrum are a high tritium activity (1 × 10 11 Bq), a low background rate (≲0.1 counts per second (cps)), an electron-volts-scale energy resolution and an accurate (0.1% level) theoretical prediction of the integral spectrum.
The KATRIN experiment tackles these challenges by combining a high-activity molecular tritium source with a high-resolution spectrometer of the magnetic adiabatic collimation and electrostatic (MAC-E)-filter type 16 . The experiment is hosted by the Tritium Laboratory Karlsruhe, allowing the safe supply of tritium at the 10 g scale with continuous tritium reprocessing 23,24 . Figure 1 shows the 70-m-long KATRIN beamline. The isotope tritium has a short half-life of 12.3 years and a low endpoint of 18.6 keV, both of which are favourable properties to achieve high count rates near the endpoint. Moreover, the theoretical calculation of the super-allowed decay of tritium is well understood [25][26][27] . The strong gaseous tritium source (nominal activity, 1 × 10 11 Bq) is established by the throughput of 40 g per day of molecular tritium through the 10-m-long source tube, which is cooled to 30 K to reduce the thermal motion of tritium molecules.
A system of 24 superconducting magnets 28 guides the electrons out of the source towards the spectrometer section and to the focal-plane detector. Between the source and main spectrometer, differential 29 and cryogenic 30 pumping systems reduce the flow of tritium by 14 orders of magnitude.
High-precision electron spectroscopy is achieved with the spectrometer of MAC-E-filter type 31,32 . The spectrometer acts as a sharp electrostatic high-pass filter, transmitting only electrons (with charge q = −e) above an adjustable energy threshold qU, where U is the retarding potential applied at the spectrometer. A reduction in the magnetic-field strength by about four orders of magnitude from the entrance (exit) of the spectrometer B src = 2.5 T (B max = 4.2 T) to its centre, the analysing plane (B ana = 6.3 × 10 −4 T), collimates the electron momenta. This configuration creates a narrow filter width of ΔE = 18.6 keV × (B ana /B max ) = 2.8 eV and allows for a large angular acceptance, with a maximum pitch angle for the β-decay electrons of θmax = arcsin √ (Bsrc/Bmax) = 50.4° in the source. The pitch angle refers to the angle between the electron's momentum and the direction of magnetic field at the position of the electron. A 12.6-m-diameter magnetic coil system surrounding the spectrometer finely shapes the magnetic field and compensates the Earth's magnetic field 33,34 .
The transmitted electrons are detected by a 148 pixel silicon-PIN-diode focal-plane detector installed at the exit of the spectrometer 35 . By measuring the count rate of transmitted electrons for a set of qU values, the integral β-decay spectrum is obtained. The main spectrometer is preceded by a smaller pre-spectrometer, which operates on the same principle and transmits only electrons above 10 keV, reducing the flux of electrons into the main spectrometer. The upstream end of the beamline is terminated with a gold-plated rear wall, which absorbs the non-transmitted β-electrons and defines the reference electric potential of the source. The rear wall is biased to a voltage O (100 mV) to minimize the difference in the surface potential to that of the beam tube, which minimizes the inhomogeneity of the source electric potential.
The rear section is equipped with an angular-and energy-selective electron gun 36 , which is used to precisely determine the scattering probability of electrons with the source gas, governed by the product of column density (number of molecules per square centimetre along the length of the source) and scattering cross section. Furthermore, we use the electron gun to measure the distribution of energy losses for 18.6 keV electrons scattering off the molecular tritium gas, providing one of the most precise energy-loss functions for this process to date 37 . Another key calibration source is gaseous krypton, which can be co-circulated with the tritium

Fig. 1 | Illustration of the 70-m-long KATRIN beamline.
The main components are labelled. The transport of β-electrons and magnetic adiabatic collimation of their momenta p are illustrated. a-f, View into the tritium source depicts three systematic effects: molecular excitations during β-decay (a), scattering of electrons off the gas molecules (b) and spatial distribution of the electric potential in the source U src (r, z) (c). The view into the spectrometer illustrates the main background processes arising from radon decays inside the volume of the spectrometer 47-49 (d), highly excited rydberg atoms sputtered off from the structural material via α-decays of 210 Po (refs. [44][45][46] ) (e) and positive ions created in a Penning trap between the two spectrometers 57 (f). Low-energy electrons, created in the volume as a consequence of radon decays or rydberg-atom ionizations, can be accelerated by qU ana towards the focal-plane detector, making them indistinguishable from signal electrons 66 .
gas 38 . Mono-energetic conversion electrons from the decay of the metastable state 83m Kr are used to determine spatial and temporal variations in the electric potential of the tritium source. These variations are caused by a weak cold-magnetized plasma, which arises from the high magnetic field (2.5 T) and a large number of ions and low-energy electrons (~1 × 10 12 m −3 ) in the tritium source. The methods of calibration are described in more detail in Methods. The beamline is equipped with multiple monitoring devices. The throughput of tritium gas within the tritium source tube and tritium circulation loop is measured by a flow meter. A laser Raman system continuously monitors the gas composition, providing a measurement at the ≤0.05%-precision level each minute. A silicon drift detector system installed in the transport section and a β-induced X-ray system at the rear section 39 continuously monitor the tritium activity, yielding a result at the 0.03%-precision level each minute. The high voltage of the main spectrometer is continuously measured at the parts-per-million level with a high-precision voltage divider system 40,41 and an additional monitoring spectrometer 42 . The magnetic fields are determined with a high-precision magnetic-field sensor system 43 .
The current background level of KATRIN of about 220 mcps mainly originates from the spectrometer section. The dominant source of background arises from α-decays of 210 Po (refs. [44][45][46] in the structural material of the spectrometer vessel. The recoiling 206 Pb creates highly excited Rydberg states at the inner spectrometer surfaces, which can be ionized during propagation in the inner volume by thermal radiation. The second source is 219 Rn and 220 Rn decays in the spectrometer volume, creating primary electrons that are magnetically trapped and slowly cool down by ionizing the residual gas in the spectrometer [47][48][49] . The resulting low-energy electrons of both sources are accelerated by retarding energy qU ana towards the focal-plane detector, making them indistinguishable from signal electrons using the energy information only. After the successful commissioning of the complete KATRIN beamline in the summer of 2017 (ref. 50 ), the first tritium operation was demonstrated with a small tritium activity of 5 × 10 8 Bq in mid-2018 (ref. 51 ). During the first KATRIN neutrino-mass (KNM1) campaign in 2019 (ref. 17 ), the source was operated in a 'burn-in' configuration at a reduced activity of 2.5 × 10 10 Bq, which is required when structural materials are exposed to high amounts of tritium for the first time. Major technical achievements of the second KATRIN neutrino-mass (KNM2) campaign are the operation of the tritium source at its nominal activity of 9.5 × 10 10 Bq and improved vacuum conditions by baking of the spectrometer 52 to 200 °C for approximately two weeks that led to a reduction in the background by 25% to 220 mcps. The thermal conditioning of the surfaces reduces the coverage of weakly bound atoms (which contribute to the background) and removes water molecules from the cryogenic copper baffles (which improves the capture efficiency for radon emanating from the main getter pumps located behind the baffles 53 ).
We, thus, increased the β-electron-to-background ratio by a factor of 2.7 with respect to the first campaign. In the last 40 eV of the integral spectrum, we collected a total number of 3.7 × 10 6 β-electrons. Figure 2a compares the spectra of both neutrino-mass campaigns. A direct comparison of the experimental parameters is given in Table 1.

Measurement of tritium β-decay spectrum
The integral β-decay spectrum is obtained by repeatedly measuring the count rate R data (qU i ) for a set of 39 non-equidistant high-voltage settings U i , creating retarding energy qU i for electrons with charge q. The retarding energy is adjusted in the range of qU i ∈ (E 0 -300 eV, E 0 + 135 eV), where E 0 = 18,574 eV is the approximate spectral endpoint. Note that for the spectral fit, only 28 out of those points in the range of qU i ∈ [E 0 − 40 eV, E 0 + 135 eV] are used. Data points further below the endpoint are used to monitor the activity stability, complementing the other monitoring devices mentioned above. The time spent at each high-voltage set point (called the scan step) varies between 17 and 576 s, which corresponds to a total measurement The graph illustrates a reduced background rate, higher signal strength and overall higher statistics (smaller error bars). The fit description is given in equation (2). b, Normalized residuals for the fit (blue line) to the data from KNM2. The shaded areas indicate statistical and total uncertainties. The contribution of systematic uncertainties is derived with the covariance matrix method, explained in the main text. c, Data points with statistical error (multiplied by a factor of 10) for the 12 Fig. 2e, the measurement-time distribution compensates the steeply falling count rate towards the endpoint and peaks at approximately 10 eV below the endpoint, where a neutrino-mass signal would be at the maximum. Based on predefined experimental conditions and data-quality criteria, 361 golden scans were chosen for the presented analysis. Since the high-voltage values can be set with a reproducibility at the sub-parts-per-million level, these scans are later combined into a single spectrum by adding the counts at a given set point.
The focal-plane detector ( Fig. 1, inset) is segmented into 148 individual pixels to account for spatial variations in the electromagnetic fields inside the source and spectrometer, as well as of the background. Removing malfunctioning pixels and those in the outer rings shadowed by parts of the beamline upstream, 117 pixels have been selected. For the presented analysis, these golden pixels are grouped into 12 concentric rings, resulting in 336 data points R data (qU i , r j ), where i ∈ {1,…,n qU = 28} and j ∈ {1,…,n rings = 12} per scan. Figure 2c displays the spectra for each of the 12 detector rings, where all the scans have been combined.

Data analysis strategy
To infer the neutrino mass, we fit the spectral data with a spectrum prediction, given by an analytical description of the β-decay spectrum and the experimental response function, described in the following sections.
Spectrum calculation. The prediction of detection rate R(qU i , r j ) is given by a convolution of the differential β-decay spectrum R β (E) with experimental response function f(E, qU i , r j ) and background rate R bg (qU i , r j ): Here N T is the signal normalization calculated from the number of tritium atoms in the source, maximum acceptance angle and detection efficiency. Also, A s ≈ 1 is an additional normalization factor, which is used as a free parameter in the fit. The differential β-decay spectrum R β (E) is given by Fermi's theory. In the analysis, we include radiative corrections and the molecular final-state distribution 54,55 in the differential β-decay spectrum. Doppler broadening due to the finite thermal motion of tritium molecules in the source, as well as energy broadenings due to spatial and temporal variations in the spectrometer and source electric potential, is emulated by Gaussian broadenings of the final-state distribution.
The response function f(E, qU i , r j ) includes energy losses due to scattering and synchrotron radiation in the high-magnetic-field regions, as well as the transmission of electrons through the main spectrometer. Compared with previous KATRIN analyses, we now use a modified transmission function to account for the non-isotropy of β-electrons leaving the source. The angular distribution of electrons leaving the source is slightly non-isotropic due to the effective path length of electrons through the source, which depends on the pitch angle. A detailed description of the spectrum calculation can be found in ref. 56 and Methods.
Unbiased parameter inference. We infer m 2 ν by fitting the experimental spectrum R data (qU i , r j ) with prediction R(qU i , r j ) by minimizing the standard χ 2 = R data C -1 R function, where C contains the variance and covariance of the data points. In addition to the neutrino mass squared (m 2 ν ), the parameters A s (r j ), R bg (r j ) and effective endpoint E 0 (r j ) are treated as independent parameters for the 12 detector rings, leading to a total number of (1 + (3 × 12)) = 37 free parameters in the fit. The introduction of ring-dependent parameters was chosen to allow for possible unaccounted radial effects. In particular, the effective endpoint E 0 (r j ) could account for a possible radial dependence of the electric potential in the source. However, the final fit revealed a negligible (<100 meV) radial variation in the endpoint. Another advantage of ring-dependent parameters is to avoid the averaging of transmission function over all the rings, which would introduce energy broadening and hence reduce the resolution.
The following analysis procedure was implemented to minimize the potential for human-induced biases. As the first step, the full analysis is performed on a Monte Carlo copy of each scan, simulated based on the actual experimental parameters (such as magnetic fields and column density), obtained by external calibration measurements. In the next step, the fit is performed on the experimental dataset, but with a randomly broadened molecular final-state distribution, which imposes an unknown bias on the observable m 2 ν . In the final step, the analysis of experimental data with unmodified final-state distribution is executed. To prevent human-induced errors, each analysis step was performed by four independent analysis teams, using different software and strategies. β-electron-to-background ratio 3.7 9.9 KNM1 refers to the first KATrIN campaign results 17 and KNM2 refers to this work. The total number of β-electrons is counted in the last 40 eV interval of the integral spectrum, which is used for the spectral fit. The β-electron-to-background ratio is given by the ratio of this number and the background counts in the same energy range, that is, E 0 − 40 eV to E 0 . Consistent results were obtained at each analysis step before proceeding to the next stage.
Systematic effects. The analytical description R(qU i , r j ) of the integral β-decay spectrum contains various experimental and theoretical parameters, such as magnetic fields, column density and probability for given molecular excitations, which are known with a certain accuracy. Different techniques (based on covariance matrices, Monte Carlo propagation, nuisance parameters and Bayesian priors) are applied to propagate these systematic uncertainties in the final result. A detailed description of these methods can be found in Methods. The neutrino-mass result presented here is dominated by statistical uncertainty. The largest systematic uncertainties are related to background properties and source electric potential. First, radon decays in the main spectrometer lead to a non-Poissonian background-rate overdispersion 49 of about 11%, effectively increasing the statistical uncertainty. Second, mechanisms for background generation may show a retarding-potential dependence of the background, parametrized by a slope (m qU bg = ( 0 ± 4.74) mcps keV -1 ). Third, a removal of stored electrons from a known Penning trap between the spectrometers 57 after each scan step can lead to a slowly increasing background rate during each scan (m tscan bg = (3 ± 3) µcps s -1 ) and thus to a scanstep-duration-dependent background contribution. Finally, spatial and temporal variations in the source electric potential modify the spectral shape and therefore lead to a relevant systematic uncertainty for the neutrino-mass measurement. The impact of all the systematic effects on the neutrino mass is listed in Table 2 and described in detail in Methods.

Best-fit result and upper limit
The χ 2 minimization reveals an excellent goodness of fit with a χ 2 per degree of freedom of 0.9, corresponding to a p value of 0.8. For the best fit of the neutrino mass squared, we find m 2 ν = 0.26 +0.34 −0.34 eV 2 based on the Monte Carlo propagation technique (Fig. 3). The independent analysis methods agree within about 5% of the total uncertainty. The total uncertainty on the fit is dominated by the statistical error followed by uncertainties in the background parameters and source electric potential. The full breakdown of uncertainties can be found in Table 2.
Based on the best-fit result, we obtain an upper limit of m ν < 0.9 eV at 90% CL using the Lokhov-Tkachov method 58 . The Feldman-Cousins technique 59 yields the same limit for the obtained best fit. This result is slightly higher than the sensitivity of 0.7 eV due to the positive fit value, which is consistent with ~0.8σ statistical fluctuation assuming a true neutrino mass of zero. We also perform a Bayesian analysis of the dataset, with a positive flat prior on m 2 ν . The resulting Bayesian limit at 90% credibility is m ν < 0.85 eV.
The ring-averaged fitted effective endpoint is E 0 = (18,573.69 ± 0.03) eV. The Q value defines the energy release in a nuclear reaction. Taking into account the centre-of-mass molecular recoil of T 2 (1.72 eV), as well as the absolute electric source potential ϕ src (σ(ϕ src ) = 0.6 V) and work function of the main spectrometer ϕ MS (σ(ϕ MS ) = 0.06 eV), we find a Q value of (18,575.20 ± 0.60) eV, which is consistent with the previous KATRIN neutrino-mass campaign, namely, (18,575.20 ± 0.50) eV (ref. 51 ), and the calculated Q value from 3 He− 3 H atomic-mass difference of (18,575.72 ± 0.07) eV (ref. 60 ). The good agreement underlines the stability and accuracy of the absolute energy scale of the apparatus.
We combined the neutrino-mass results from this work with the previously published KATRIN (KNM1) results 17,18 . A simultaneous fit of both datasets yields m 2 ν = (0.1 ± 0.3) eV 2 and a corresponding upper limit of m ν < 0.8 eV at 90% CL, based on the Lokhov-Tkachov or Feldman-Cousins technique. The same result is obtained when   67 , Tokyo (1991) 68 , Zürich (1992) 69 , Mainz (1993) 70 , beijing (1993) 71 , Livermore (1995) 72 , Troitsk (1995) 73 , Mainz (1999) 13 , Troitsk (1999) 74  multiplying the m 2 ν distributions from Monte Carlo propagation or adding the χ 2 profile of the individual fits. As both datasets are dominated by statistics, correlated systematic uncertainties between both campaigns are negligible. Furthermore, we investigated a Bayesian combination, where the KNM1 posterior distribution of m 2 ν is used as the prior for the KNM2 analysis, yielding consistent results. More details on the combined analyses can be found in the Supplementary Information.

Conclusion and outlook
The second neutrino-mass measurement campaign of KATRIN, presented here, reached sub-electronvolt sensitivity (0.7 eV at 90% CL). Combined with the first campaign, we set an improved upper limit of m ν < 0.8 eV (90% CL). We, therefore, have narrowed down the allowed range of quasi-degenerate neutrino-mass models and we have provided model-independent information about the neutrino mass, which allows the testing of non-standard cosmological models 61,62 . Figure 4 displays the evolution of the best-fit m 2 ν results from historical neutrino-mass measurements up to the present day.
Compared with its previous measurement campaign, the KATRIN experiment has decreased the statistical and systematic uncertainties by about a factor of three and two, respectively. With the total planned measurement time of 1,000 days, the total statistics of KATRIN will be increased by another factor of 50. A reduction in the background rate by a factor of two will be achieved by an optimized electromagnetic-field design of the spectrometer section 45,63 . Moreover, by eliminating the radon-and Penning-trap-induced background, the background-related systematic effects are expected to be substancially reduced. Together with a high-rate 83m Kr calibration scheme and a method to improve the determination of magnetic fields, we will minimize the dominant systematic uncertainties to reach the target sensitivity of m ν in the vicinity of 0.2 eV.
High-precision β-spectroscopy with the KATRIN experiment has proven to be a powerful means to probe the neutrino mass with unprecedented sensitivity and to explore physics beyond the standard model, such as sterile neutrinos 64 . Together with cosmological probes and searches for neutrinoless double-β decay 65 , the upcoming KATRIN data will play a key role in measuring the neutrino-mass parameters.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41567-021-01463-1.

Methods
Here we describe the data analysis chain starting from data processing to high-level fit and limit setting. Moreover, we provide details on one of the key calibration campaigns, concerning the source electric potential. A more extensive description of the KATRIN analysis procedure is provided elsewhere 18 .
Data processing, selection and combination. The first step of the analysis chain is data preparation. Raw data are combined into integral spectral data points, which are then fitted with an analytical spectrum prediction including the response of the experiment.
Rate determination. Electrons transmitted through the main spectrometer are further accelerated by a post-acceleration electrode with potential U PAE = 10 keV before they are detected by the focal-plane detector 35 . The latter provides a high detection efficiency (>95%) and moderate energy resolution (2.8 keV full-width at half-maximum at 28 keV). The total rate per pixel at a given retarding potential qU i is determined by integrating the rate in a wide and asymmetric region of interest (ROI) of 14 ≤ E ≤ 32 keV. The asymmetric ROI is chosen to account for energy losses in the dead layer of the Si PIN diodes 78 , for partial energy deposition due to the backscattering of electrons off the detector surface as well as charge sharing between pixels. The detector efficiency slightly depends on qU i . Three effects are considered.
(1) The differential energy spectrum of the transmitted electrons at the focal-plane detector is shifted (and slightly scaled) according to qU i . Since the same ROI is used for each qU i , some electrons are no longer covered by the fixed ROI when lowering the potential, which effectively changes the detection efficiency. Based on reference measurements, the relative reduction in detection efficiency at E 0 − 300 eV, with respect to the efficiency at E 0 , is determined to be δ ROI = 0.00048, with an uncertainty of 0.16%. The data are corrected according to the efficiency at the given qU i value.
(2) The count rate at the focal-plane detector as well as the probability of pile-up depends on qU i . The energy of most pile-up events is added such that they are not covered by the ROI, thereby effectively changing the detector efficiency with qU i . Based on a random-coincidence model assuming a Poisson-distributed signal, the relative reduction in efficiency at E 0 − 300 eV is estimated to be δ PU = 0.00017, with an uncertainty of 18%. The data are corrected accordingly. (3) Finally, a fraction of about 20% of all the electrons impinging on the detector surface are backscattered. The majority of these electrons are back-reflected to the same detector pixel (within the integration time of the energy filter) by the pinch magnetic field or retarding potential. However, for low retarding potentials and small energy depositions in the focal-plane detector, the backscattered electrons have a chance of remaining undetected by overcoming the retarding potential a second time in the direction towards the tritium source. The lower the qU i , the higher is the probability for lost electrons, effectively changing the detection efficiency. Based on Monte Carlo simulations, we estimate the efficiency reduction to be δ BS < 0.001 at E 0 − 300 eV. We neglect this effect in this analysis. Systematic uncertainties related to efficiency corrections are negligibly small for the presented analysis, which only considers the last 40 eV of the tritium spectrum.
Selection of golden pixels and scans. The focal-plane detector is segmented into 148 pixels of equal area. For the analysis, 117 pixels have been chosen. The rejected pixels show undesirable characteristics such as broadened energy resolution or higher noise levels (a total of 6 pixels), partial shadowing by the silicon-drift-detector-based beam monitor (a single pixel) and decreased rate due to misalignment of the beamline with respect to the magnetic-flux tube (a total of 24 pixels) that stems from tiny deviations in the orientation of the 24 superconducting magnets from the design 28 .
Out of the 397 scans recorded, 361 passed a strict quality assessment and were included in the neutrino-mass analysis. The other runs were rejected for reasons of failed set points of spectrometer electrodes (9 runs), downtime of the laser Raman system (23 runs) and other operational conditions (4 runs). The golden pixel and run selection is performed before the unblinding of the data.
Data combination. The 117 golden pixels are grouped into 12 detector rings and the detector counts recorded within each ring are summed to obtain 12 independent spectra. All the golden scans are combined by adding the counts recorded at the same retarding energies. This method leads to a total number of data points of n data = n qU × n rings = 336.

Spectrum calculation.
To infer the neutrino mass, the integral spectrum is fitted with the analytical spectrum calculation including the experimental response. As shown in equation (2), the theoretical spectrum is composed of the differential tritium spectrum R β (E) and experimental response function f(E, qU i ). The differential spectrum is given by where G F denotes the Fermi constant, cos 2 Θ C is the Cabibbo angle, |M nucl | 2 is the energy-independent nuclear matrix element and F(E, Z′ = 2) is the Fermi function.
where E 0 denotes the maximum kinetic energy of the electron in the case of zero neutrino mass and V f describes the molecular excitation energies populated with probabilities ζ f . Also, E and m e denote the kinetic energy and mass of the β-electron, respectively. Beyond the molecular effects, further theoretical corrections arise on the atomic and nuclear level 79 . Only radiative corrections to the differential spectrum are relevant for this analysis, which are included in the analytical description, but not shown here 56 . Furthermore, we emulate the effect of Doppler broadening, as well as spatial and temporal source and spectrometer electric potential variations, by broadening the final-state distribution with a Gaussian distribution (discussed in the main text).
The experimental response function depicts the probability of an electron with starting energy E to reach the focal-plane detector. It combines the transmission function T (E − ϵ, θ, U) of the main spectrometer and electron's energy losses ϵ due to inelastic scattering with the tritium molecules in the source. The scattering energy losses are described by the product of the s-fold scattering probabilities P s (θ), which depend on the path length through the source and hence on the pitch angle θ (angle between the energy momentum and magnetic field at the position of the electron), and the energy-loss function f s (ϵ) for a given number of scatterings s. The energy-loss function is experimentally determined with the electron gun installed at the rear of the source. A typical response function and the corresponding energy-loss function is shown in Extended Data Fig. 1. The integrated transmission function for an isotropic source of electrons is given by It is governed by the magnetic fields at the starting position B src of the electron, the maximum field B max in the beamline and the magnetic field in the spectrometer's analysing plane B ana . The acceptance angle of the MAC-E-type filter is given by θmax = arcsin √ (Bsrc/Bmax) = 50.4°. Synchrotron energy losses 56 of β-electrons in the high magnetic field in the source and transport systems are included as an analytical correction to the transmission function (not shown here).
Systematic uncertainties. The analytical description R(qU i , r j ) of the integral β-decay spectrum (equation (2)) contains various signal-and background-related parameters, which are known with a certain accuracy. In the following, we describe these parameters, their uncertainties and their treatment in the neutrino-mass analysis.
Signal-related systematic effects. Spectrum prediction includes uncertainties in the nine parameters of the empirical energy-loss function (the individual relative uncertainties of parameters σ eloss,k are between 0.016% and 3.800%); a relative uncertainty of the product of column density and scattering cross section (σ ρd×σ = 0.2500%); relative uncertainties of the theoretical description of the molecular final-state distribution (σ FSD = O (1.000%)); and relative uncertainties of the magnetic field in the source (σ Bsrc = 1.700%), in the analysing plane (σ Bana = 1.000%) and the maximal field (σ Bmax = 0.100%). The variation in β-decay activity during a scan was determined from the product of column density (as determined from the tritium throughput) and tritium purity or alternatively from the beam monitor (silicon drift detector system in the transport section). The obtained standard deviation of the relative activity is negligibly small (σ scan = 0.035%).
The spatial and temporal variations in the source electric potential were not included in the first neutrino-mass campaign. With the increase in source column density from 1.11 to 4.23 × 10 17 cm -2 in this campaign, the creation rates and densities of the charge carriers, as well as the fraction of scattered electrons, increased accordingly, which makes plasma effects more relevant 15,80,81 . A detailed description of the plasma calibration is given in the main text.
Background-related systematic effects. Electrons created during radon decay in the main spectrometer volume can initially be magnetically trapped. Each of these trapped electrons can create a cluster of 10-100 secondary electrons by scattering off the residual gas in the main spectrometer 49,66 , which is operated at a pressure of 10 −11 mbar. These secondary electrons arrive at the focal-plane detector within a time window of about 1,000 s, hence leading to a non-Poissonian rate distribution 49 . The observed background rate can be modelled by a Gaussian distribution, with a width that is 11.2% wider than expected from a purely Poisson distribution. This overdispersion is treated as an increased statistical uncertainty in the analysis.
As the transmission conditions for the background electrons slightly depend on the retarding-potential setting, a small retarding-potential dependence of the background can occur. In the analysis, we allow for a linear dependence of the background on the retarding potential and use dedicated test measurements to constrain the possible slope to m qU bg = (0 ± 4.7) mcps keV -1 . Finally, the pre-and main spectrometers, both being operated at a high voltage, create a Penning trap between them. Stored electrons and the subsequently produced positive ions, which can escape the trap into the main spectrometer, are a source of background, as illustrated in Fig. 1. To mitigate this background, the trap is emptied with an electron-catcher system 57 after each scan step. However, a potentially small increase in the background rate within a scan step cannot be excluded, which could lead to a background dependence on the duration of the scan step. By fitting a linear increase to the rate evolution within all the scan steps, we find a slope of m tscan bg = (3 ± 3) µcps s -1 , which is included in the β-decay spectrum fit. This effect was first observed in later physics runs, in which the scan-step duration was significantly increased. It was, therefore, only included after the data were already unblinded to a subgroup of the collaboration. As the evaluation of the size of this uncertainty was provided by an independent task group of the collaboration, the reported result remains free of bias.
Source-potential calibration. The absolute electric potential of the source does not affect the spectral shape of the measured spectrum. An unknown absolute source potential is largely absorbed by the effective endpoint, which is a free parameter in the fit. A change in the effective endpoint mostly leads to a shift in the spectrum and has a negligible effect on the spectral shape. Accordingly, an unknown radial variation in the electric potential is in good approximation absorbed by the ringwise endpoint parameters. Consequently, these effects have a minor impact on the neutrino-mass analysis. However, any temporal variation, like high-frequency plasma instabilities, slow drifts of the absolute potential or longitudinal variations of the source potential can lead to spectral distortions, which are parameterized by Gaussian broadening σ P and parameter Δ P that quantifies the longitudinal rearto-front asymmetry of the electric potential of the source. This asymmetry of the potential results in a shift in the energy spectrum associated with the scattered electrons (predominantly originating from the rear of the source tube) compared with the spectrum of the unscattered electrons (predominantly originating from the front of the source tube). It enters into the model via equation (4), where the energy-loss function f s=1 (ϵ) of singly scattered (s = 1) electrons is shifted by Δ P relative to the energy-loss function f s=0 (ϵ) of unscattered (s = 0) electrons.
Both parameters are assessed with the help of co-circulating 83m Kr gas, assuming that the possible plasma instabilities or longitudinal plasma profile are not affected by its presence in a minute concentration. The spectroscopy of its mono-energetic conversion electron lines reveals information about the broadening σ P of the lineshape, from which an upper limit of Δ P is derived 82 .
The calibration was performed at an elevated source temperature of T = 80 K to prevent the condensation of Kr gas (compared with the T = 30 K set point used during the neutrino-mass measurements). The tritium circulation loop of the tritium source 83 can be operated in two modes: (1) in a mode with the direct recycling of the krypton-tritium mixture, which is limited to a maximum column density of 40% (2.08 × 10 17 cm −2 ), but delivers a high 83m Kr rate; (2) in a mode of fractional direct recycling, which can be operated up to 75% (3.75 × 10 17 cm −2 ) of the nominal column density, but with only about 0.5% of the maximum 83m Kr activity compared with the mode described in (1).
The plasma parameters are inferred by combining the measurements of internal conversion lines N 2,3 -32 and L 3 -32. The energy of the conversion electrons, emitted from a particular subshell of the 83m Kr atom, is 30,472.6 eV for the L 3 -32 singlet and 32,137.1 eV (32,137.8 eV) for the N 2 -32 (N 3 -32) doublet 84 . The L 3 -32 line, on one hand, has a high intensity, but its natural linewidth is not precisely known 85 . The N 2,3 -32 lines, on the other hand, have a low intensity, but their natural linewidth is negligible compared with spectral broadening caused by variations in the electric potential. Thus, any broadening of the N 2,3 -32 lines beyond the known spectrometer resolution and thermal Doppler broadening can be assigned to variations in the electric potential within the source. Based on the N 2,3 -32-line doublet measurement at 40% of the nominal column density, we find σ 2 P (40%) = (1.4 ± 0.3) × 10 -3 eV 2 (Extended Data Fig. 2). Combined with L 3 -32 line measurements at both 40% and 75% of the nominal column density, we assess the relative change in broadening and find σ 2 P (75%) = (8.0 ± 8.2) × 10 -3 eV 2 . Finally, we conservatively extrapolate by exponentially scaling this value to 84% of the nominal column density, leading to σ 2 P (84%) = (12.4 ± 16.1) × 10 -3 eV 2 . The plasma parameters and uncertainties deduced from the calibration measurements at 80 K cover the conditions during the neutrino-mass measurements at 30 K. The effect of different temperatures is smaller than the assumed uncertainty. For future physics runs, the tritium source will be operated at 80 K.
The contribution to broadening from slow drifts and shifts as determined from the run-wise analysis of the tritium spectra during KNM2 only yields a value of 10 −3 eV 2 . This value is significantly smaller and can be determined with much higher precision than broadening σ 2 P obtained from the 83m Kr measurement. In the analysis, we limit broadening σ P to positive values. We construct the probability density function of the asymmetry parameter Δ P based on phenomenological considerations on the relation of both plasma parameters given by |Δ P | ≤ σ P /1.3 (ref. 82 ). For many samples of σ 2 P , we draw a value for Δ P from a uniform distribution in the range of −σ P /1.3 ≤ Δ P ≤ σ P /1.3. The resulting distribution can be approximated by a Gaussian distribution centred around 0 meV with a width of 61 meV.
We expect to reduce the systematic uncertainties in future campaigns by operating the tritium source at the same temperature during the neutrino-mass and krypton measurements and by using an ultrahigh intensity 83m Kr source (with about 10 GBq of the 83 Rb parent radionuclide), which will make the N 2,3 -32 line measurement possible at a nominal column density (mode b).
Parameter inference. We infer the parameters of interest via a minimization of the χ 2 function where R data (qU, r) gives the measured count rates at a retarding potential qU i for the detector ring r j , R(qU, r) gives the predictions of these rates and C is the covariance matrix that includes the statistical uncertainties and can also be used to describe systematic uncertainties. The usage of a χ 2 minimization is justified, as the numbers of electrons per scan step qU i and detector ring r j are sufficiently large (>700) to be described by a Gaussian distribution instead of a Poisson distribution. The fit has (1 + (3 × 12)) = 37 free parameters Θ, including a single parameter for the neutrino mass squared m 2 ν , and 12 ring-wise values for each of the three spectral parameters, namely, normalization factor A s , background rate R bg and effective endpoint E 0 . In addition, the spectrum depends on systematic parameters η (Extended Data Table 1), such as the column density, tritium isotopologue concentrations, magnetic fields and so on. These parameters are known with a certain accuracy, and their uncertainty needs to be propagated to the final neutrino-mass result. Four different methods are used for the KATRIN analysis, as discussed below.
Pull method. In the 'pull method' , systematic parameters η i are treated as free parameters in the fit and introduced as nuisance terms in the χ 2 function: The nuisance terms allow the parameter to vary around its best estimation η i according to its uncertainty ση i as determined from external measurements. This method is computationally intensive due to the complexity in calculating the tritium spectrum and minimization with respect to multiple free parameters. For example, it is not practical to treat the uncertainties of the molecular final-state distribution, which is given as a discrete list of excitation energies and corresponding probabilities, with this method. The advantage of this method is that we make the maximum use of the data. If the spectral data contain information about the systematic parameters η, it is automatically taken into account.
Covariance matrix method. As shown in equation (7), the standard χ 2 estimator includes covariance matrix C, which can describe both statistical and systematic model uncertainties. The diagonal entries describe uncertainties that are uncorrelated for each R(qU i , r j ), whereas the off-diagonal terms describe the correlated uncertainties between R(qU i , r j ). Covariance matrix C is computed by simulating O(10 4 ) β-decay spectra, with systematic parameters η i varied according to their probability density functions in each spectrum. From the resulting set of spectra, the variance and covariance of the spectral points R(qU i , r j ) are determined.
As the covariance matrices for individual or combined systematic effects are computed before fitting, this method is efficient with respect to computational costs. The dimension of the covariance matrix is given by the number of data points. Therefore, the efforts for matrix calculation and inversion can be diminished by reducing the number of data points.
Monte Carlo propagation method. Generally, in the Monte Carlo propagation technique, the uncertainties of parameters η i are propagated by repeating the fit about 10 5 times, with the systematic parameters varied according to their probability density functions each time. The method then returns the distributions of the fit parameter Θ, which reflects the uncertainty of the systematic parameters.
More precisely, when assessing the systematic effects alone, a simulated spectrum without statistical fluctuations is created (based on the best-fit parameters), which is then fitted 10 5 times with a model whose systematic parameters of interest are varied each time. Contrarily, when evaluating the statistical uncertainty alone, 10 5 statistically randomized Monte Carlo spectra are created and fitted with a constant model. Accordingly, for obtaining the total uncertainty, both steps are combined.
To extract information on parameters η from the data, each entry in the resulting histogram of the best-fit parameters is weighted with the likelihood of the corresponding fit. The final distributions are then used to estimate the best-fit values (mode of distribution) and uncertainty (integration of distribution up to 16% from both sides).
The advantage of this method is that the number of free parameters is kept at a minimum, which facilitates the minimization procedure. The larger number of fits, however, is time consuming and requires the usage of large computing clusters.
Bayesian inference. In Bayesian inference, one computes the posterior probability for the parameters of interest Θ from a prior probability and a likelihood function according to Bayes' theorem. For the KATRIN analysis, the Bayesian approach has the advantage that the prior knowledge of neutrino mass can naturally be applied via a corresponding informative prior. For example, the prior of m 2 ν can be chosen to be flat and positive, which restricts the posterior distribution to the physically allowed regime and accordingly changes the credible interval.
Ideally, all systematic effects η i would be included as free parameters constrained with our prior knowledge on the parameter. However, due to the computationally expensive spectrum calculation and the fact that the Bayesian inference requires a large number of samples in the Markov chain for sampling the posterior distribution, only the qU-dependent background is currently treated in this way. All the other systematic uncertainties are included with a covariance matrix in the likelihood or by model variation. For the model variation, a large set of Markov chains is started with randomized but fixed model variations. The randomizations are drawn from the systematic uncertainty distributions. The resulting set of posterior distributions is averaged.
Limit setting. We present two frequentist methods and one Bayesian method for the construction of an upper limit of the neutrino mass. For the former, we use the classical Feldman-Cousins 59 and Lokhov-Tkachov 58 belt constructions (Extended Data Fig. 3). In the Feldman-Cousins technique, the acceptance region for m 2 ν is determined by ordering this estimator according to the likelihood ratio L( m 2 ν | max(0, m 2 ν )) for a given best-fit neutrino mass squared m 2 ν . This method leads to more stringent upper limits for increasingly negative best-fit values. The Lokhov-Tkachov method avoids this feature by using the standard Neyman belt construction for positive values of m 2 ν and defining the experimental sensitivity as the upper limit in the non-physical regime of m 2 ν < 0. The Bayesian 90% credible interval is obtained by integrating the posterior distribution of m 2 ν from zero to m 2 ν limit , such that the total probability is 90% (Extended Data Fig. 4). Note that the interpretation of the limit obtained in this way is different from the frequentist confidence limits, and hence, the numerical values may not coincide.
Results of individual strategies. The frequentist analyses are performed by three independent teams, which use differing implementations of spectral calculation and different strategies to propagate systematic uncertainties. A Bayesian analysis is performed, which interfaces to one of the spectrum calculation software. The analysis by independent teams is a powerful means to cross-check individual analyses. The following four data analysis strategies were applied to the second dataset of KATRIN.
• Strategy 1 is based on a C++ framework (C++14 / Gnu compiler) 86 using the Minuit minimizer. It mainly employs the pull method for handling systematics. We note that for the presented analysis, the input value for the 'scan-step-duration-dependent background' was corrected after the official unblinding of the data. • Strategy 2 is implemented in a MATLAB framework (R2020a) and exclusively uses the covariance matrix approach to propagate systematic uncertainties 87 . We note that for the presented analysis, the 'scan-step-duration-dependent background' systematic was implemented only after the official unblinding of the data. • Strategy 3 is based on a C++ framework (C++17 / Gnu compiler) using a custom-developed minimizer 88,89 . In this strategy, the Monte Carlo propagation of uncertainties is mostly applied. • Strategy 4 performs a Bayesian interpretation of the data. For this approach, the spectrum calculation software of 'Strategy 3' is interfaced with the Bayesian analysis toolkit 90 . Here most of the systematic uncertainties are treated via the model variation technique.
• An overview of the strategies is found in Extended Data Table 2. The resulting best fit and systematic uncertainty breakdown are listed in Extended Data Table 3.

Data availability
Source data are provided with this paper. Further, we provide the χ from one of the analysis strategies only for the KNM2 and combined KNM1/KNM2 datasets.

Code availability
The code used in this study is available from the corresponding authors upon reasonable request.