Analysis methods for the first KATRIN neutrino-mass measurement

We report on the data set, data handling, and detailed analysis techniques of the first neutrino-mass measurement by the Karlsruhe Tritium Neutrino (KATRIN) experiment, which probes the absolute neutrino-mass scale via the $\beta$-decay kinematics of molecular tritium. The source is highly pure, cryogenic T$_2$ gas. The $\beta$ electrons are guided along magnetic field lines toward a high-resolution, integrating spectrometer for energy analysis. A silicon detector counts $\beta$ electrons above the energy threshold of the spectrometer, so that a scan of the thresholds produces a precise measurement of the high-energy spectral tail. After detailed theoretical studies, simulations, and commissioning measurements, extending from the molecular final-state distribution to inelastic scattering in the source to subtleties of the electromagnetic fields, our independent, blind analyses allow us to set an upper limit of 1.1 eV on the neutrino-mass scale at a 90\% confidence level. This first result, based on a few weeks of running at a reduced source intensity and dominated by statistical uncertainty, improves on prior limits by nearly a factor of two. This result establishes an analysis framework for future KATRIN measurements, and provides important input to both particle theory and cosmology.

the energy threshold of the spectrometer, so that a scan of the thresholds produces a precise measurement of the high-energy spectral tail. After detailed theoretical studies, simulations, and commissioning measurements, extending from the molecular final-state distribution to inelastic scattering in the source to subtleties of the electromagnetic fields, our independent, blind analyses allow us to set an upper limit of 1.1 eV on the neutrino-mass scale at a 90% confidence level. This first result, based on a few weeks of running at a reduced source intensity and dominated by statistical uncertainty, improves on prior limits by nearly a factor of two. This result establishes an analysis framework for future KATRIN measurements, and provides important input to both particle theory and cosmology. The absolute mass scale of the neutrino remains a key open question in contemporary physics, with farreaching implications from cosmology to elementary particle physics. Despite numerous efforts along three complementary lines of approach (observational cosmology, the search for neutrinoless double-β decay, and direct searches using the kinematics of weak-interaction processes such as single β decay or electron capture), only upper bounds on the neutrino mass have been found so far (see, e.g., [1][2][3] for reviews on these subjects). Meanwhile, neutrino flavor-oscillation experiments (e.g., [4,5]) have firmly established the existence of non-zero neutrino masses.

CONTENTS
With the advent of precision cosmology, corresponding bounds on neutrino masses have been dramatically improved, and now form the tightest constraints available. Yet, cosmological bounds on m i (the sum of the distinct neutrino-mass eigenvalues m i ) are derived using the paradigm of the cosmological standard model (ΛCDM), and the values obtained vary with the selection of data sets included in the analysis. The Planck collaboration has inferred robust bounds from cosmicmicrowave-background power spectra alone: m i < 0.26 eV (95 % confidence level, CL), which can be further improved to m i < 0.12 eV (95 % CL) by including lensing and baryon-acoustic-oscillation data [6]. Meanwhile, laboratory searches for neutrinoless double-β decay are sensitive to the neutrino-mass scale, under the assumption that neutrinos are Majorana particles that make the dominant contribution to the decay mechanism. Here, the observable is the coherent sum of weighted neutrino mass values m ββ = | U 2 ei m i |, where U ei denotes the electron-flavor element coupled to the i th neutrino-mass state in the neutrino mixing matrix. Presently, the most sensitive limits on m ββ are set by searches in 76 Ge (GERDA, 0.07 -0.16 eV) [7] and in 136 Xe (KamLAND-Zen, 0.06 -0.17 eV) [8]. The ranges of these 90% confidence limits arise from uncertainties in nuclear-matrix elements.
Direct laboratory-based measurements are an indispensable model-independent probe of the neutrino-mass scale, resting solely on the determination of kinematic parameters. Two weak processes particularly suitable for this quest are the electron capture of 163 Ho [9,10] and the β decay of tritium: The kinematics of these decays provide access to the effective neutrino-mass square value, an incoherent sum over the weighted squares of the mass values m i (i = 1, 2, 3): Historically, the Mainz and Troitsk experiments used tritium to set the previous most stringent direct upper limit at m ν < 2 eV (95 % CL) [11,12] with a highaccuracy shape measurement of the β-decay spectrum in the vicinity of its kinematic endpoint (E 0 = 18.57 keV for molecular tritium, T 2 ). Meanwhile, the mass splittings measured in oscillation experiments impose a lower limit on this observable. Depending on the ordering of the pattern of neutrino-mass eigenstates ν i , this floor is either approximately 8 meV (normal ordering) or 50 meV (inverted ordering) -see, e.g., Ref. [13].
The Karlsruhe Tritium Neutrino (KATRIN) experiment [14,15] is further improving this approach to target a neutrino-mass sensitivity of 0.2 eV (90 % CL) after five years of measurement time; note the change to 90% confidence level. This goal requires an improvement of about two orders of magnitude in the m 2 ν observable. To accomplish this challenging measurement, KATRIN relies on the proven technology of the MAC-E filter (Magnetic Adiabatic Collimation with an Electrostatic filter, developed for neutrino-mass measurements by the Mainz and Troitsk groups [16,17]) and a large β-decay luminosity provided by a gaseous molecular tritium source (following pioneering work at the Los Alamos experiment [18]). After commissioning and characterizing the complex 70 mlong electron beamline, initially with monoenergetic calibration sources [19] and subsequently with first-tritium β electrons [20], the KATRIN collaboration has recently reported an improved upper limit on the neutrino mass of m ν < 1.1 eV (90 % CL) based on an initial four-week science run [21]. This result yields an improvement of about a factor of two with respect to the best previous direct bound.
In this work, we present a detailed account of the data set acquired, data-handling and analysis techniques applied, and statistical inference methods employed to derive this result. In the following we will use the term "KATRIN Neutrino Mass run 1" (KNM1) to label the inaugural four-week science campaign that marks the first operation of KATRIN at high tritium purity, at about a quarter of the nominal tritium source strength. During KNM1, an integrated β spectrum was acquired over a "full" energy interval stretching from about 90 eV below to about 50 eV above the endpoint E 0 . The actual neutrino-mass analysis was performed in a narrower interval, [E 0 -37 eV, E 0 +49 eV], in which the measurement is statistics-dominated. Within this 86 eV analysis interval, the data set comprises a total ensemble of 2.03 × 10 6 events after data-quality selection cuts. The ensemble was collected over a measurement time of 521.7 h and is composed of 1.48 × 10 6 β decay electrons below E 0 and 0.55 × 10 6 events in a flat background over the entire analysis interval.
We begin this paper with an overview of the experimental setup (Sec. II) and the configuration in which the KATRIN beamline was operated, including data handling and measurement strategy (Sec. III). (For reference, Table I lists abbreviations frequently used in this paper.) Two key ingredients of the analysis, the β-spectrum model and the instrument response function, are presented in Secs. IV and V. Section VI summarizes relevant sources of background and their characteristics. General underlying principles of the analysis, in which data from the individual detector pixels and β-spectrum scans are combined into a single spectrum for fitting, are given in Sec. VII.
Section VIII presents a detailed assay of individual systematic uncertainties. Section IX documents the strategy employed for blind analysis, describes two complementary methods employed to propagate the systematic uncertainties into the neutrino-mass fit, and shows the resulting spectral fit and uncertainty breakdown. Section X details the construction of the confidence belt and the derivation of the neutrino-mass upper limit via the Feldman-Cousins [22] and the Lokhov-Tkachov [23] approaches. Our Lokhov-Tkachov result of m ν < 1.1 eV (90 % CL), presented in Ref. [21], was obtained using Frequentist methods. In this work we also present a derivation of the upper limit based on Bayesian methods, yielding a limit of m ν < 0.9 eV (90% C.I.) (Sec. XI). This method uses a different approach to deal with the unphysical region of negative neutrino-mass squared.
In Sec. XII, as a consistency check of KATRIN's absolute energy scale, we show that the effective endpoint value E 0 obtained from the fit to the β spectrum agrees with independent measurements of the Q-value through the 3 He-T mass difference.
We conclude by summarizing our findings (Sec. XIII) and discussing them in the wider context of contemporary neutrino-mass probes (Sec. XIV). Figure 1 gives an overview of the KATRIN apparatus. Briefly, in order to ensure sufficient statistics, a bright tritium source produces some 2.45 × 10 10 decays each second in the KNM1 configuration. In order to perform a fine-grained energy analysis near the tritium endpoint, the energies of the resulting β electrons are analyzed by a pair of MAC-E-filter spectrometers [16,17]. These basic functions require the support of extensive systems for handling the tritium gas, maintaining vacuum conditions, ensuring adiabatic electron transport, mitigating or eliminating backgrounds, detecting β electrons, and calibrating and monitoring the apparatus as a whole. The resulting 70 m beamline is described in detail in Ref. [24]; here, we offer a brief summary. T 2 gas from a temperature-and pressure-controlled buffer vessel at 313 K is cooled to 30 K and continuously injected via a capillary into the center of the source sys-tem. The resulting Windowless Gaseous Tritium Source (WGTS) freely streams to both ends of the system, where it is continuously pumped away with turbomolecular pumps. This results in a stable pressure distribution inside the source beam tube [25]. Once the T 2 gas is pumped away, it flows over a PdAg membrane filter that is permeable only to hydrogen isotopes. A constant fraction of the circulating gas is also removed at this stage for later purification, and is replaced with highly pure T 2 directly after the filter. The purified gas is fed back to the temperature-and pressure-controlled buffer vessel, forming a closed loop. The loop system is integrated with the infrastructure of the Tritium Laboratory Karlsruhe, which provides tritium purification of exhaust gas, tritium storage, and fresh tritium supply for KATRIN [26][27][28].

II. KATRIN EXPERIMENTAL SETUP
Within the 10 m-long, 90 mm-diameter source beam tube [29], tritium decays produce β electrons that are guided along magnetic field lines [30] through the rest of the experimental beamline. At the upstream end, the WGTS terminates in a gold-plated rear wall, which can be held at a fixed potential and/or illuminated with ultraviolet light to liberate photoelectrons. At the downstream end, the windowless nature of the source is essential to avoid catastrophic energy loss, but necessitates other means for the confinement of tritium. The β electrons are first guided around magnetic chicanes through two pumping stages, namely a differential pumping system and a cryogenic pumping system, which collectively reduce the partial pressure of tritium by more than 14 orders of magnitude [31]. Specially designed electrodes within the differential stage [32] prevent the transmission of tritium ions.
β electrons must then pass through a pair of MAC-Efilter spectrometers, operated in tandem. Each MAC-E filter is characterized by strong magnetic fields at the entrance and exit, with a region of weak magnetic field in the center. Since the magnetic moment is conserved in the adiabatic transport of the electrons through the beamline, the electron momenta rotate to become approximately parallel to the magnetic field lines, producing a broad, roughly collimated beam. A longitudinal retarding potential therefore analyzes the total kinetic energy of the electrons at the central "analyzing plane," at which the magnetic field is the weakest. Electrons below the resulting energy threshold are reflected upstream, toward the source; electrons above the energy threshold are transmitted downstream, toward the spectrometer exit. The transmission function of the spectrometers was extensively calibrated prior to the measurement (Sec. V).
The first MAC-E filter in the tandem pair, the prespectrometer [33], has a fixed energy threshold at 10 keV and removes the bulk of the low-energy electrons. Immediately downstream, the main spectrometer is the high-resolution, adjustable-threshold filter that analyzes the integral β spectrum. Each data-taking "scan" (Sec. III G) consists of a sequence of main-spectrometer retarding-potential settings, with a new threshold of in-Main spectrometer

Rear system
Source system Differential pumping system

Detector system
Monitor spectrometer

Cryogenic pumping system
Pre-spectrometer FIG. 1. Overview of the 70 m KATRIN beamline. Moving downstream, from left to right, the major components are: the rear system, the source system, the differential pumping system, the cryogenic pumping system, the pre-spectrometer, the main spectrometer, and the detector system. The monitor spectrometer monitors the retarding potential of the main spectrometer.
tegration at each setting. The electropolished interior stainless-steel surface of the main spectrometer is lined with two layers of inner, wire electrodes, providing fine shaping of the electric fields and, when operated at a negative potential offset from the main-spectrometer vessel, electrostatic rejection of low-energy secondary electrons from the main-spectrometer surface [34]. The vessel potential is supplied by a commercial system, with additional regulation and post-regulation designed and built by the collaboration to suppress 50 Hz mains noise and other sources of interference [35]. Air-cooled magnetic coils, mounted on a framework surrounding the main spectrometer, compensate for the Earth's magnetic field, fringe fields of the solenoids, and residual magnetization [36]. The ultra-high vacuum in the spectrometer is maintained by non-evaporable getter strips and turbomolecular pumps [37]. Liquid-nitrogen-cooled copper baffles are positioned across the pump ports to suppress background electrons due to radon decay in the main volume [38,39]. To mitigate backgrounds from the Penning trap between the two MAC-E filters, a conductive electron catcher is inserted into the inter-spectrometer region at each change in the set voltage of the main spectrometer [40]. This device removes trapped electrons that would produce secondary ions and electrons.
Electrons that pass through the main spectrometer undergo additional acceleration via the post-acceleration electrode, improving rejection of non-spectrometer backgrounds. When they reach the detector system, they are counted in the focal-plane detector (FPD) [41], a monolithic silicon p-i-n diode segmented into 148 equalarea pixels. The FPD and its readout electronics are elevated to the post-acceleration potential, and preamplified signals are transmitted to the data-acquisition (DAQ) system via optical fiber. Each FPD pulse is digitized in a 12-bit analog-to-digital converter (ADC), and its amplitude and timing are reconstructed online by the sequential application of two trapezoidal fil-ters [41,42]. These values are then recorded using the Object-oriented Real-time Control and Acquisition (ORCA) framework [43], which can also communicate directly with the main-spectrometer high-voltage system using a web-based database tool [44]. Pulse amplitudes are translated into energies in near-time processing (Sec. III F), based on the results of regular calibration runs with an 241 Am photon source.
Multiple calibration and monitoring systems provide essential information during both neutrino-mass scans and dedicated runs [45]. In the tritium loops feeding the source, a laser-Raman spectroscopy system (LARA) [46][47][48] monitors the relative concentrations of hydrogen isotopologs, particularly T 2 , DT, and HT, within the source gas. In the rear system upstream of the source, an electron gun (e-gun), following the design of a similar e-gun used for testing the main spectrometer [49], serves as an angle-and energy-selective calibration source. This egun delivers electrons through an aperture in the rear wall at the upstream end of the source. Observed in the FPD, these electrons test the response function of the experiment as a whole. Two radioactive, in-vacuum calibration sources are also available: gaseous 83m Kr that can be circulated within the source when its temperature is elevated to about 100 K [50], and a condensed 83m Kr source that can be inserted into the cryogenic pumping system [51].
Upstream of the rear wall, the β-induced x-ray spectroscopy system continuously monitors the source activity: silicon drift detectors view x-rays produced by β electrons scattering in the rear wall [52]. Further downstream, within the cryogenic pumping system, a forward beam monitor provides complementary activity monitoring [53]. This monitor includes two silicon p-i-n diodes for electron rate and spectrum measurements, a Hall sensor, and a temperature gauge. A vacuum manipulator allows these sensors to be positioned radially within the beam; normally, the forward beam monitor is positioned at the outer edge of the β electron flux.
The main-spectrometer retarding potential, which defines the energy analysis, is continuously monitored both by a voltage divider with demonstrated part-per-million (ppm) precision [54][55][56][57] and by the refurbished MAC-E filter from the historical Mainz experiment [11]. Now relocated to KATRIN, this monitor spectrometer references the main-spectrometer retarding potential to an atomic standard via synchronous scans of a 83m Kr conversion line [58].
Prior to the KNM1 neutrino-mass run, the full KA-TRIN beamline was commissioned with photoelectrons, ions, and 83m Kr conversion electrons in 2016-2017 [19], and with small amounts of tritium in D 2 carrier gas in 2018 [20]. Subsequently, in another campaign with D 2 , the electron gun was commissioned and gas properties of the source were investigated [59]. KNM1 marked the first time that the inner surfaces of the injection capillary and source system were exposed to large amounts of tritium. Radiochemical reactions between T 2 and these metal surfaces produced both CO and tritiated methane, which condensed on the cold metal surface of the capillary and partially obstructed tritium flow over time.
To improve stability during this burn-in period, KA-TRIN operated at a reduced column density of ρd exp = 1.11 × 10 17 molecules/cm 2 .

III. THE KNM1 MEASUREMENT CAMPAIGN
In this section we describe the operating conditions of the KATRIN experiment during its first high-purity tritium campaign (KNM1), which took place from 10 th April to 13 th May 2019. In particular, we characterize the system performance in terms of the source-gas isotopic purity (Sec. III A) and column density (Sec. III B), as well as the reproducibility, homogeneity, and stability of the electron starting potential in the source (Sec. III C) and the retarding potential in the analyzing plane (Sec. III D). We also discuss the detection of β electrons and the definition of a region of interest (Sec. III E) as well as the processing and analysis pipeline for the data (Sec. III F). The requirements for system stability arise from the method adopted to measure the tritium β spectrum by repeatedly scanning the retarding potential in alternating up and down sweeps (Sec. III G), and from the fact that KNM1 data from all pixels and all scans are combined into a single spectrum for fitting. In the final analysis, then, experimental parameters are essentially averaged over both space (across the detector) and time (across like scan steps throughout the KNM1 data-taking period). Later on, Sec. VII explores the justification for this analysis method in the statistics-dominated KNM1 data set.
For the KNM1 campaign, the sequence of scan steps, each consisting of a retarding-potential set point distributed in the interval [E 0 −91 eV, E 0 +49 eV], resulted in a typical scan duration of 2.5 h. Therein, each scan step corresponds to a measurement time varying from 17 s for high-rate points deeper in the spectrum to 576 s near the endpoint region, as will be shown later in Fig. 7.

A. Tritium source parameters
The average source activity during KNM1 neutrinomass data-taking was about 2.45 × 10 10 Bq, maintained by a column density of 1.11 × 10 17 molecules/cm 2 . This was achieved by a cumulative tritium throughput of 4.9 g/d.
The gas injected into the source consists mainly of molecular T 2 . Due to initial impurities and exchange reactions with the stainless-steel piping and vessel, the other hydrogen isotopologs (H 2 , HD, HT, D 2 , and DT) are also present in minor fractions. A PdAg membrane (permeator) in the tritium loop [60] continuously filters the circulated tritium gas to prevent the recirculation of built-up impurities. The relative fractions c x of the six hydrogen isotopologs are continuously monitored by LARA, downstream of the permeator. The relative molecular isotopolog fractions c x and the atomic tritium purity ε T are defined as: where N x is the number of molecules of isotopolog x in the source, and the sums are over all six isotopologs. The tritium purity is monitored with better than 10 −3 statistical precision [48]. The time evolution of the relative fractions of the three tritiated isotopologs injected into the source during KNM1 is shown in Fig. 2. On average, the concentrations of the tritiated species throughout the campaign were c T2 = 0.953, c HT = 0.035, and c DT = 0.011; these values are used in the final neutrino-mass analysis. The resulting tritium purity is ε T = 0.9758(13) [48]. The prominence of HT as a secondary species is due to exchange reactions with H atoms that are naturally present in stainless-steel piping [59], and the residual presence of DT is due to the isotope-separation process used to purify the tritium [61]. The inactive species (H 2 , HD, and D 2 ) are only present in trace amounts, as they are strongly suppressed by shifts of the chemical equilibrium in the presence of high-surplus T 2 .

B. Column density
The column density ρd determines the number of tritium atoms N T in the source where A is the cross-sectional area of the WGTS, and the factor of 2 is necessary because ρd is defined in terms of the number of T 2 molecules. The column density further defines the s-fold scattering probabilities P s of electrons, traveling parallel to magnetic field lines through the entire tritium source, with the gas molecules: The product ρdσ, where σ is the cross section for inelastic scattering of electrons from molecular tritium (Sec. V B), gives the expected number of scatterings. It must be known with high accuracy for the analysis [25]. The precise absolute value of ρdσ is obtained from measurements with the narrow-angle, quasi-monoenergetic egun located in the rear system. This e-gun produces a high-intensity beam of electrons via the photoelectric effect according to the principle described in Ref. [49]. On their path towards the detector, the electrons traverse the source, where they can undergo inelastic scattering and in the process lose energy. Only those electrons with sufficient remaining energy to surpass the spectrometer potential are counted in the detector. By measuring the electron rate at different retarding potentials and fitting a model response function (Sec. V) to these data, we may make a precise determination of ρdσ.
E-gun electrons differ from β-electrons with regard to their starting positions and their energy and angular dis-tributions. For this reason a modified response function, including a precise description of the e-gun beam characteristics, is used in the column-density determination. The e-gun electron rate is measured at retarding potentials where the impact of the column density is the strongest. The mean energy of the e-gun electrons is set to 18.78 keV, allowing a clean separation from β electrons that could bias the column-density determination. During KNM1, ρdσ was determined with the e-gun on a weekly basis, achieving relative uncertainties of less than 0.9 %.
As described in Sec. II, the first exposure of the inner loop to T 2 resulted in the production of gas species which condensed on the surface of the injection capillary. This obstruction caused the tritium injection flow and column density to drift over time at constant tritium injection pressure. By lowering the column density to be a factor of approximately 5 smaller than the nominal column density ρd nom , and by increasing the tritium injection pressure several times during KNM1, these drifts were kept lower than 3 %. To ensure precise monitoring of the column density during the whole measurement period, the e-gun measurements were combined with continuous ρd fluctuation data from a mass-flow meter with 200 sccm full-scale range [63], applied to the tritium injection flow. The reproducibility of the flow meter during KNM1 is conservatively estimated to be 1.5 µbar · l/s. Based on simulations that show a linear relation between ρdσ and the tritium injection flow for a narrow throughput range [64], a linear calibration function is suitable to relate the measured throughput to ρdσ.
With this strategy, we determine the column density with high precision for all tritium data-taking. The time evolution and distribution of the column-density values are shown in Fig. 3; the average value of ρdσ is 0.404 at the molecular tritium endpoint. Using the cross-section value from Eq. 17 further below, this value translates to an average column density of ρd = 1.11 × 10 17 molecules/cm 2 .

C. Electron starting potential
The starting potential of the β electrons is provided by a cold and strongly magnetized plasma in the WGTS. The magnitude of the potential depends on the boundary conditions at the rear wall and the grounded beam tube. By optimization of the rear-wall set voltage, a homogeneous, stable plasma potential can be created. This is important because both spatial inhomogeneities and temporal fluctuations of the plasma potential distort our spectrum in a manner analogous to the neutrino mass. Indeed, the shift in neutrino-mass squared due to an error ∆σ 2 in the Gaussian variance of a continuous variable (such as the starting potential of the β electrons) is given at leading order by [65]:  3. Evolution of the column density during KNM1; the uncertainty is dominated by systematics arising from the relationship of e-gun data to the measured throughput, and from fluctuations in the latter quantity. The visible decrease of the column density over time is caused by conductance changes of the tritium injection capillary. By increasing the tritium injection pressure several times, the column density was stabilized.
Since we combine all pixels and all scans for our KNM1 fits (Sec. VII), our analysis does not account for inhomogeneities or temporal fluctuations, and the full variance of the electron starting potential therefore contributes via Eq. 7. The source plasma is generated by the weakly selfionizing tritium gas. According to simulation, each β electron creates on average 36 secondary electrons, and thus 36 positive ions, through scattering interactions. Throughout the central part of the WGTS, the ions have a mean free path of less than 0.5 m for momentum transfer with the neutral gas. Consequently, the flow of neutral tritium gas drives the ions toward both ends of the source. The low-energy, secondary electrons follow the ion motion in order to maintain quasi-neutrality, facilitated by their much higher mobility along the magnetic field lines. While the ions quickly become fully thermalized to the meV scale, the energy spectrum of secondary and β electrons ranges from meV to keV.
The electric potential inside the plasma depends on the surface potentials at its boundaries. These are determined in turn by their intrinsic work functions φ, which are expected to differ by several 100 mV [66], and by the applied bias voltages. As the beam tube is grounded (U bt = 0 V), only the rear-wall bias voltage U RW remains to compensate the work-function differences. At an optimal U RW , the radial and longitudinal inhomogeneities of the plasma potential both vanish, as expected from simulations with the assumption of negligible work-function inhomogeneities [67].
The optimal rear-wall bias voltage was determined by measuring the β-rate at various U RW settings. Comparing these rates to reference spectra, we extracted the dependence of the spectral endpoint E 0 on the FPD ring number -which correlates to radius in the source.
For U RW = −150 mV, a flat radial E 0 distribution was found. Also, the measurement of the plasma-induced cur-rent on the rear wall showed no drifts and less noise than at other bias voltages. U RW was therefore set to −150 mV for the measurement campaign.
The systematic effect of remaining spatial inhomogeneities and fluctuations of the plasma potential can be constrained by studying the line widths and positions of quasi-monoenergetic conversion electrons from gaseous 83m Kr co-circulating in the T 2 gas [68]. The L 3 -32 line at 30 472.2(5) eV is particularly interesting for this study. First, it is located above E 0 . Second, the 37.8(5) % branching ratio into this final state leads to a high signal-to-noise ratio [69]. Third, it possesses a small intrinsic line width of Γ ≈ 1 eV. In a previous campaign using gaseous 83m Kr in the absence of tritium gas [19], the KATRIN experiment measured an L 3 -32 line position of E L3-32 = 30472.604 ± 0.003 stat ± 0.025 sys eV and a Lorentzian line-width of Γ L3-32 = 1.152 ± 0.007 stat ± 0.013 sys eV [70]. This effective line position includes a shift arising from the absolute work-function difference between the source and the main spectrometer.
After the KNM1 neutrino-mass campaign ended, plasma studies were performed for two days with co-circulating 83m Kr and T 2 .
It should be noted that the column density during neutrino-mass measurements was only 22 % of the nominal value of 5.0 × 10 17 molecules/cm 2 , while during the plasma study it was about 30 % of the nominal value.
The krypton admixture did not affect general plasma properties, such as charged-particle density or electric potentials, because the partial pressure and activity (≈ 3 MBq) of krypton were several orders of magnitude below those of tritium (≈ 33 GBq). However, the plasma was affected by the beam-tube temperature of 100 K, elevated from the nominal 30 K. This higher temperature was necessary to prevent the krypton from freezing, but also increased the temperature of the dominant low-energy part of the electron energy distribution [71]. The plasma temperature is known to strongly influence the rate of electron-ion recombination at the meV scale. As the recombination rate is much stronger at 30 K, we expect plasma effects at elevated source temperature to be more prominent. We thus use results obtained during the krypton measurement at 100 K to set an upper limit of the scale of possible plasma effects. The β-decay electrons and non-thermalized electrons make only minor contributions to the number density, but their dominant role in the energy density of charged particles requires a detailed investigation.
The intrinsic Lorentzian line width was measured with gaseous 83m Kr in the absence of tritium, with the experimental conditions as similar as possible to the L 3 -32 measurements with co-circulating T 2 / 83m Kr (described above). By comparing these two measurements and assuming an energy-independent background, we find that the presence of T 2 results in a Gaussian line broadening of < 80 mV for rear-wall settings in the range −350 mV< U RW < 350 mV. The collaboration is currently investigating the impact of a possible radial-dependent back-ground, which could arise due to detector effects.
The impact of these findings on the neutrino-mass measurement is discussed in Sec. VIII C.

D. Analyzing-plane potentials
The threshold energy for electrons to pass through the MAC-E filter is determined by the value of the retarding potential U at the analyzing plane. Any unknown instabilities in the retarding potential directly affect the energy scale of the tritium spectrum and can introduce systematic effects on m 2 ν . To first order, significant, unaccounted-for continuous inhomogeneities and fluctuations effectively broaden the spectrum as seen in Eq. 7. Our KNM1 analysis does not account for inhomogeneities or fluctuations in U , so that the full variance is seen in the broadening. For the target sensitivity of KATRIN, the energy scale must be stable to within 60 meV or 3 ppm on a baseline retarding potential of −18.6 kV.
To achieve this, we have constructed a dedicated measurement chain, including precision high-voltage dividers with proven long-term stability on the ppm level over one year [56]. A custom-built post-regulation system [35] ensures stability at higher frequencies, up to 1 MHz.
In order to stack multiple scans for the KNM1 analysis (Sec. VII B), it is not only necessary to have 3 ppm monitoring, but also to achieve comparable precision in both the stability at each scan step and the reproducibility of the retarding potential from scan to scan. Figure 4 shows the achieved high-voltage stability while acquiring data at individual scan steps over the full measurement interval. This stability is on average below 15 mV, significantly exceeding requirements. The observed increase in standard deviation as a function of scan-step duration is described well by a simple statistical model that combines a random-walk diffusion process with a feedback loop. The reproducibility of retarding potentials from scan to scan follows a Gaussian distribution with a width of σ =34(1) mV. This limitation of the reproducibility is directly related to the digital-to-analog converter inside the post-regulation setup; for measurement phases after KNM1, finer-grained regulation is in place.
The retarding potential is continuously monitored during the measurements. Therefore, at each scan step, the time evolution of the retarding potential is known with ppm precision. Neglecting this in the analysis introduces an additional broadening of the energy scale, leading to a neutrino-mass shift of ∆m 2 ν = −3 × 10 −3 eV 2 . This shift is less than half the allotment for the high-voltage-related systematic uncertainty in the KATRIN uncertainty budget for full five-year statistics [15], and can be neglected in the KNM1 analysis.

E. Electron counting and region of interest
The FPD records a low-resolution, differential spectrum of electrons that have passed the high-resolution energy threshold set by the main spectrometer. Measuring the integrated tritium β spectrum for KNM1, and thereby extracting the neutrino mass, requires an accurate count of electrons that arrive at the FPD within an energy region of interest (ROI) during each scan step. The ROI cut allows rejection of backgrounds and noise events generated near or in the FPD.
When electrons strike the FPD, its pixels are triggered individually, with thresholds set just above the noise floor at around 5 keV. As described in detail in earlier work [41], the energy and timing for each pulse are reconstructed online using a double trapezoidal filter and then recorded; FPD waveforms are not saved during normal operations. The shaping length of the trapezoidal-filter pair is set to 1.6 µs, optimizing the energy resolution at around 1.8 keV (full width at half maximum, FWHM). During β scans, rates are too low for significant pileup, but severe pileup during high-rate e-gun measurements can result in deadtime when multiple coincident events drive the baseline out of the ADC dynamic range. This effect is mitigated by individually adjusting the gain of each channel to approximately 5 ADC counts per keV, preserving good energy resolution while defining a dynamic range (up to 400 keV) sufficient to accommodate pileup. These settings were implemented in the DAQ firmware prior to the KNM1 measurement. Simulations of the readout chain show that the fraction of time during which the baseline is shifted out of the ADC input FIG. 5. FPD pixel selection for KNM1. All 117 selected pixels, colored in solid green, are used for the analysis. The two pixels filled with horizontal blue lines are excluded due to shadowing by the forward beam monitor. The six pixels filled with vertical orange lines are excluded due to intrinsic noisy behaviour. All pixels filled with gray circles are excluded since they are partially shadowed by beamline components. range is less than 0.05 % for 50 kcps of 28.6 keV electrons, a 100-fold improvement compared to previous settings.
Out of the 148 pixels, we define a list of 117 selected detector pixels, distributed as shown in Fig. 5. The excluded pixels are either noisy, or shadowed by beamline instrumentation in the β-electron path along the magnetic flux tube.
Electrons that transit the spectrometer (Sec. V) receive an additional 10 keV of kinetic energy from the postacceleration electrode, and 120 eV from the bias voltage applied to the FPD. For a retarding potential around 18.6 kV, this results in a broad peak in the FPD energy spectrum at around 28 keV (Fig. 6). Background electrons and β electrons share this characteristic energy spectrum in the FPD, since the primary background during KNM1 arises from low-energy electrons that are created inside the main spectrometer and then accelerated by the retarding potential (Sec. VI). The FPD energy scale is calibrated with a 241 Am gamma source every two weeks.
Our ROI is defined as [14 keV, 32 keV], as measured by the FPD (Fig. 6). The upper bound of the ROI is determined simply from the peak position and the peak width; the lower bound is determined for stability and robustness. In contrast to earlier studies that considered backgrounds originating near the detector [41], the choice of a low-energy KNM1 ROI lower bound does not reduce the signal-to-background ratio, since an energy cut cannot differentiate between β electrons and mainspectrometer background. A cut far away from the peak, where the spectrum shape derivative is small, improves stability against fluctuations of energy scale and resolution. Consequently, corrections for peak-position depen- dence on retarding potential are negligible.
The specific lower bound of the ROI, 14 keV, was chosen so as to cancel two effects that arise from charge sharing, in which energy from a single incident electron is divided between two neighboring pixels. If a pixel loses more than half the event charge, its loss from the ROI decreases the effective rate; if a pixel receives more than half the event charge, its inclusion in the ROI increases the effective rate. With the FPD threshold set at half the peak energy, these two effects exactly compensate each other.

F. Data pipeline
Following each pixel trigger (Sec. III E), the DAQ records the trigger timestamp from a 20 MHz clock and the energy information as raw ADC counts integrated over the shaping time of the trapezoidal filter. A scan is divided into scan steps. Each scan step is defined by its HV set point, and its duration is determined according to the measurement-time distribution of the scan (Sec. III G). Prior to acquisition start at each scan step, handshakes between the DAQ and the HV control system ensure that the HV read-back value has reached the setpoint value within a defined accuracy of 50 mV, as measured by a four-point moving average over the last 8 s. The inter-spectrometer electron catcher is inserted and removed during this change of scan steps, so that it does not obstruct the beamline during data-taking. A series of pulse-per-second (PPS) pulses from a precision clock synchronized to the Global Positioning System (GPS) defines both the start and stop times of scan steps, providing boundary time accuracy better than 1 ns. The 50 ns digitization timestamps are also phase-locked to 10 MHz pulses from the same precision clock. The readout system is capable of handling a pixel rate of 100 kcps and a total rate of 3 Mcps. Therefore, no deadtime is expected for the actual tritium scan, which has a maximum count rate of 7 kcps. A typical two-hour scan produces roughly 120 MB of data.
Immediately after completion of a scan, data files are processed automatically. This processing includes the transfer to storage computers, time-wise event sorting, conversion to offline data formats, and indexing into a run database, followed by automated user-side analysis including the reduction of data in user-specified data files. Except for the handshakes between the DAQ and HV systems, slow-control channels are independent from the tritium scans. Each slow-control sensor has a defined recording interval, typically between 2 and 10 s. This is a heterogeneous system for which timestamps are taken from computer timestamps synchronized to the Network Time Protocol (NTP). In the offline analysis, special care is taken for synchronization among different slow-control channels, as well as between the DAQ and slow controls.
An intermediate data layer, consisting of user-side shared data storage with version management, splits the data analysis chain. The first half of the chain covers analysis at the event and time-series levels, and the second half provides higher-level analysis including model fitting. For each scan, results of the first-level analysis are summarized in digest files that contain analyzed FPD counts with efficiency corrections, individual scan steps, calibrated slow-control values (including LARA isotopolog concentration and column density, and analyzed rates extracted from β-induced x-ray spectroscopy and the forward beam monitor), and data-quality flags. Some experimental parameters, such as beamline alignment information and magnetic-and electric-field values determined by measurements and simulations, are shared across all scans in a given measurement period; each such period is summarized in a digest file containing the values of these parameters.
During data-taking, acquisition occasionally began before the HV readback values achieved stability due to minor synchronization errors. The first two seconds of every scan step were removed from the data to address these issues. Count-rate, livetime, efficiency, and stability calculations are performed after these data-quality cuts.
G. Acquisition of the integral β decay spectrum KATRIN measures the integral tritium β decay spectrum by sequentially applying different retarding energies qU , or equivalently HV settings, to the main spectrometer and counting the rate of transmitted β electrons, R(qU ), with the FPD. Our choice of the scan steps -that is, the HV set points and the measurement time at each set point -maximizes the sensitivity for m 2 ν by focusing on a narrow region where the impact of the neutrino mass on the spectrum is most pronounced. The location of this region depends on the experimental conditions; in the KNM1 campaign, it lies at E 0 −14 eV [72]. Figure 7 shows the measurement-time distribution used during this campaign, developed using a nominal value of E 0 = 18 574 eV. The spectrum is scanned repeatedly over the range [E 0 −91 eV, E 0 +49 eV] by sequentially applying the non-equidistant HV values (each constituting one scan step) to the main spectrometer. A complete set of measurements at all 39 scan steps is defined as a scan. Each scan over this energy range takes approximately 2.5 h and is performed in alternating upward and downward directions. This mitigates the effects of any time-dependent drifts of the slow-control parameters. As explained in Sec. VIII I, the analysis interval is limited to an energy range of [E 0 −37 eV, E 0 +49 eV], consisting of 27 scan steps. A brief, additional scan step at E 0 −201 eV is used for rate-stability monitoring.
For each tritium scan, we apply quality cuts to relevant slow-control parameters to select a data set with stable run conditions. As Sec. VII describes in detail, data from all active detector pixels are summed, effectively converting the detector wafer into a single, uniform pixel for analysis. Furthermore, all 274 scans are combined by summing counts from like scan steps, forming a single spectrum for fitting.
The 27 scan steps within the analysis interval cover a total measurement time of 521.7 h, corresponding to 2.03 × 10 6 events. Table II summarizes key operational parameters and figures for events and scans, covering both the full interval and the analysis interval. The evolution of the integrated β-decay luminosity over the course of KNM1 is displayed in Fig. 8.

IV. TRITIUM-SPECTRUM MODELING
The KNM1 analysis relies on a model of the measured spectrum, which convolves the theoretical β spectrum (outlined in this section) with the experimental response function (details in Sec. V). We first describe the general theory of β-decay in Sec. IV A, along with some straightforward corrections. To account for the physics of KA-TRIN's molecular source (T 2 with some HT and DT), we then address the molecular final-state distribution (FSD) in detail in Sec. IV B. Since an error in the FSD variance across our measurement interval will (to first order) shift the extracted, squared neutrino-mass value according to Eq. 7 in the previous section, we have invested substantial effort in checking and extending our treatment of the FSD.
A. Theoretical β-spectrum of molecular tritium In KATRIN's molecular source, the β decay parent in Eq. 1 becomes T 2 , with a molecular decay product 3 HeT + . To model the resulting differential β spectrum, we begin with a point-like Fermi interaction, which causes the weak decay, and then apply the sudden approximation, in which the Coulomb interaction of the β electron with the remaining molecular system 3 HeT + is neglected. The validity of this approximation was demonstrated in Refs. [73,74].
Choosing the center-of-mass coordinate frame to align with the momentum of the neutrino and integrating over the experimentally unresolved neutrino and electron directions and neutrino energy, the decay rate into the nuclear and molecular configuration f of the daughter 3 HeT + at a given electron kinetic energy E reads [73] in natural units with c = = 1. m e and m ν are the electron and neutrino masses, respectively; ε f (E) has the form of the neutrino energy after energy conservation has been enforced by the Heaviside function Θ(ε f (E) − m ν ). |T f | 2 is the transition matrix element to the nuclear and molecular state f . Since the derivation of the decay rate is performed in the center-of-mass frame, which almost perfectly coincides with that centered on the decaying molecule, there is no need to integrate over the recoil momentum of the molecule; the recoil kinetic energy is naturally added as a constant energy loss.
|T f | 2 may be factorized in the sudden approximation as where T weak f 2 is independent of the electron energy for the superallowed tritium β decay. Similarly, the leptonic part T lep f 2 is independent of the electron energy in the sudden approximation. As is customary, however, the Fermi function F (E, Z = 2) (as given in Ref. [75]) is included in this factor. This allows a partial incorporation of the influence of the Coulomb interaction during the decay by accounting for the charge of an isolated 3 He daughter nucleus, leading to an effectively Coulombdistorted sudden approximation. Meanwhile, T mol f 2 is equal to the probability ζ f that 3 HeT + populates the unresolved set of molecular electronic, vibrational, and rotational states with energy V f . Since the motion of the center of mass of 3 HeT + must balance the neutrino and electron momenta, T mol f 2 theoretically depends on the electron energy after the integrations are performed. The KNM1 analysis interval is narrow enough to neglect this dependence.
After evaluating |T f | 2 according to Eq. 9, summing over the possible final nuclear states, and explicitly summing over the included range of molecular states, we obtain The prefactors include the energy-independent quantities G F (the Fermi constant), Θ C (the Cabibbo angle), and |M nucl | 2 (the nuclear matrix element). Meanwhile, where the reduced endpoint E 0 represents the total maximum electron kinetic energy in the case of a massless neutrino. While E 0 is retrieved from the fit during the neutrino-mass analysis (Sec. IX), the internal molecular excitation energies V f and the corresponding population probabilities ζ f come from computation (see Sec. IV B). The values of all constants are as in Ref. [72].
Beyond the molecular effects discussed in detail in Sec. IV B, theoretical corrections to the tritium β decay spectrum arise at the particle, nuclear, and atomic levels (see Ref. [76] for details). Of these, we include only the radiative corrections [77] in this work; these have by far the largest effect on the high-energy tail of the β electron spectrum.
Finally, the electron spectrum R β is Dopplerbroadened due to the finite motion of tritium molecules in the source. To account for this effect, we replace each discrete final state with a Gaussian centered at the finalstate energy V f , normalized to ζ f and with a standard deviation of 94 meV according to the Doppler broadening at 30 K. Effects due to the bulk gas flow are negligible.
For effects that give rise to continuous modifications of the spectrum, such as the molecular final-state distribution and Doppler broadening, a mistake in the modeled variance will introduce a bias on the extracted neutrinomass squared according to Eq. 7.

B. Final-state distribution (FSD)
Within the sudden approximation, the β decay effectively corresponds to a sudden change of the nu-clear charge of one of the tritium nuclei. This induces electronic and vibrational excitations of the daughter molecular ion 3 HeT + , possibly including its dissociation and/or ionization. Furthermore, the departing β electron and neutrino induce external (translational) and internal (rotational, vibrational, and -to a smaller extent, neglected here -electronic) excitations.
Since only the energies of the β electrons are analyzed by KATRIN, the undetected energy associated with the remaining molecular system must be computed ab initio by first solving the Schrödinger equation for the initial and final molecular systems, and then computing the transition probabilities ζ f = T mol f 2 to the molecular daughter states f thus found. Earlier calculations either focused on lower temperatures than KATRIN's 30 K [78], thus artificially constraining the population of initial molecular states, or did not include all the tritiumcontaining isotopologs [79]. In the following, we provide only a minimal description of the new computations carried out for the initial gas states relevant to KATRIN; a detailed publication is in preparation [80]. The theoretical prediction of the dissociation probability of the daughter 3 HeT + ion, following β decay, has recently been experimentally verified [81].

Solutions to the molecular Schrödinger equation
As in previous works, these computations adopt two fundamental approximations. First, the Coulombdistorted version of the sudden approximation neglects the interaction of the β electron with all but the daughter nucleus 3 He + in the β decay. Second, the Born-Oppenheimer approximation allows a separate treatment of the electronic and nuclear motions that define the full, internal molecular Schrödinger equation.
Our solution of the Schrödinger equation describing the nuclear motion uses the isotopolog-independent Born-Oppenheimer electronic potentials generated according to Ref. [82] and presented explicitly in Ref. [83]. Mass-dependent corrections are applied for the electronic ground states of specific isotopologs -T 2 , DT, HT, 3 HeT + , 3 HeD + , and 3 HeH + -and the potential curves are extended up to an internuclear separation of 20 a 0 , with a 0 the Bohr radius. Because of the rotational symmetry of the corresponding Schrödinger equation, the solutions for nuclear motion are expanded as products of spherical harmonics and radial functions. They are then augmented by the rotational barrier for non-zero initial angular momenta J i .
The electronic ground state of the daughter molecule supports about 300 rotational/vibrational bound states and a large number of predissociative resonances in the dissociation continuum. We have therefore adopted a new approach for solving the nuclear motion in these electronic potentials. Expanding the radial part in B-spline functions and adopting vanishing boundary conditions at the end of the radial grid, the solution of the Schrödinger equation is turned into a generalized matrix eigenvalue problem and requires only the diagonalization of a very sparse matrix. The spectral density and energy range of the resulting discretized spectrum may be controlled by the size of the adopted spherical box and the number of B-splines.

Energy-resolved FSD
With the newly obtained nuclear-motion solutions, and the isotopolog-independent Born-Oppenheimer electronic overlaps defined in Refs. [83] (final electronic ground state n = 1) and [82] (final electronic states n ∈ 2; 6 ), the transition probabilities between the initial and final states of interest in the KNM1 analysis interval can be obtained by integrating the matrix elements over the internuclear separation vector. The transition operator, which can be expanded into spherical Bessel functions, depends on this vector.
Compared to earlier work, our new calculation extends the results of Ref. [79] from the first 6 to the first 13 bound electronic states, and employs more accurate molecular masses than Refs. [78,79]. These more accurate masses are used in the Hamiltonian, in the fraction of the recoil momentum imparted onto the spectator nucleus -which selects the population of the states due to the molecular β decay via the transition operator -, and in the recoil energy of the whole molecular system. For the electronic excited final states n ∈ 1; 6 , we have been able to reproduce the results of Refs. [78,79] for the published initial states of T 2 (J i ∈ 0; 3 ), DT (J i ∈ 0; 1 ), and HT (J i = 0), when using the old kinematic inputs. Figure 9 shows a comparison of the current distribution with Ref. [79] for transitions from the most populated T 2 initial state at T = 30 K. The new distribution of transitions to the electronic ground state is ∼3 meV lower on average than that in Ref. [79]; this difference mostly originates from the updated recoil momentum as a consequence of the more accurately determined endpoint. Electronic final states with n > 6, combined with the electronic continuum, contribute negligibly -at the 10 −4 level -to the KNM1 analysis interval, with its lower bound at E 0 -37 eV. In our new calculation, these have been adapted for energy-scale changes from the calculations in Ref. [78]. The n > 6 bound states were omitted in Ref. [79], explaining the slightly higher transition probabilities of the new distribution around 40 eV.
For the KATRIN analysis, we consider all J i ∈ 0; 3 for all three decaying isotopologs and weight their respective contributions based on the source temperature. The Boltzmann distributions are calculated at 30 K. However, for the homonuclear T 2 molecule, the resulting J i probability must be multiplied by nuclear-spin probabilities characteristic of 700 K. The molecules in the tritium loop dissociate when they arrive at the permeator (Sec. II), which is operated at 700 K. After diffu- sion through the permeator, the atoms recombine into molecules with an ortho-para ratio of 0.75, characteristic of that temperature. The time for natural conversion to a lower-temperature ortho-para ratio is many orders of magnitude longer than the O(1 s) passage time of the molecules through the 30 K region of the injection capillary and source tube, so the T 2 gas retains an ortho-para ratio of 0.75. Weighting based on the relative concentrations of T 2 , DT and HT, as measured during KNM1, is performed at a subsequent stage of the analysis.

V. RESPONSE FUNCTION MODELING
The observed KNM1 tritium integral spectrum R(qU ) is the convolution of the differential β electron spectrum R β (E) from Eq. 10 with the instrumental response function f (E − qU ), with an added energy-independent background rate R bg : Here, N T,eff denotes the effective number of tritium atoms in the source, as adjusted by the detector efficiency and by the solid-angle acceptance of the setup ∆Ω/4π = (1 − cos θ max )/2, where θ max ≈ 50.5°as discussed below. A s is the signal amplitude.
As shown in Fig. 10, the response function f (E − qU ) [72] describes the probability of transmission of an electron with initial energy E through the beamline as a function of its surplus energy E − qU relative to the retarding potential U . Below, we discuss its calculation in detail. First, Sec. V A defines the response function and describes the effects of the beamline electromagnetic fields on the β electrons. We then treat the inelastic scattering cross section for β electrons (Sec. V B) and develop a model of energy loss experienced in flight through the KATRIN apparatus (Sec. V C).

A. Response and transmission functions
The transmission condition for any electromagnetic configuration of the KATRIN MAC-E filter determines whether an electron with starting energy E and starting angle θ is transmitted through a retarding potential U : Here, θ = ∠( p, B) is defined as the initial pitch angle of the electron, the polar angle of its momentum relative to the magnetic field: p 2 ⊥ = E sin 2 θ · (γ + 1) · m e . The Lorentz factor γ arises from its relativistic motion and has a maximum value of about 1.036 at E 0 . Meanwhile, B min = 0.63 mT is the magnetic field in the analyzing plane, B max = 4.23 T the maximum field of the beam line, and B S = 2.52 T the source magnetic field.
Only electrons with sufficient surplus energy satisfy the transmission condition and are included in the measured integral spectrum. The KATRIN main spectrometer achieves a magnetic-field ratio B min /B max ≈ 1/6700 ≈ ∆E/E, corresponding to a filter width (energy resolution) of ∆E = 2.8 eV at 18.6 keV. The maximum acceptance angle θ max = arcsin B S /B max ≈ 50.5°limits the range of pitch angles contributing to the integral spectrum. The magnetic fields and the retarding potential are provided by detailed field calculations using the Kassiopeia software [85]. To compute the precise electromagnetic fields across the analyzing plane, we use an as-built geometry of the beamline magnets with a detailed three-dimensional model of the main spectrometer. The resulting transmission conditions can be included in the model individually for each active pixel.
The detailed response function of the KATRIN apparatus is calculated from Eq. 13, as modified by energy losses between source and analyzing plane [72]: For an ensemble of electrons, f (E − qU ) depends on the acceptance angle θ max and the amount of neutral gas the electrons pass in the WGTS, which is described by the scattering probability P s (θ) and the inelasticscattering energy-loss function f s ( ) for a given number of scatters s. As Sec. V C will discuss in detail, we measure f (E − qU ) using monoenergetic electrons with small angular spread, and thus deduce f s ( ). Briefly, these electrons are produced in the e-gun with surplus energies E − qU spanning a 50 eV range. They follow the magnetic-field lines and pass through the integral column density ρd of the source. This allows us to observe single (s = 1) and multiple (s > 1) electron scatterings in the source. The scattering probability P s (θ = 0 • ) (Eq. 6) follows a Poisson distribution with the expected number of scatterings given by the product of the effective column density ρd and the inelastic-scattering cross section σ (Sec. V B).
In an isotropic source like the WGTS, electrons are emitted with an angular distribution ω(θ)dθ = sin θdθ, and we can define an integrated transmission function T (E, U ): Although analysis of non-isotropic e-gun data requires the full expression in Eq. 14, the neutrino-mass analysis in this work exploits the isotropic nature of the tritium β-source and uses the simplified response function In principle, the response function is slightly modified due to the dependence of the path length, and therefore the effective column density, on the pitch angle of the β-electrons [72]. The resulting effect on the measured endpoint is small compared to the overall uncertainties of the electric potential of the source, and this effect is not taken into account in the current analysis. Synchrotron energy losses of β-electrons in the high magnetic field in the source and transport systems are included as an analytical correction to the transmission function [72].

B. Inelastic-scattering cross section
The theoretical total inelastic-scattering cross section of electrons with T 2 molecules in the high-energy Born approximation can be written as [86][87][88]: where R H = 13.606 eV is the Rydberg energy, a 2 0 = 28.003 × 10 −18 cm 2 the Bohr radius squared, and E nr denotes the non-relativistic kinetic energy of the electron: E nr = 0.5 m e β 2 , with β 2 = 1 − m 2 e /(m e + E) 2 and E the relativistic kinetic energy of the electron. At the spectral endpoint for molecular tritium β decay, we take E = E 0 = 18.575 keV and E nr = 17.608 keV.
The dominant parameter M 2 tot can be calculated reliably and with high accuracy, since it is a special electron expectation value for the ground-state hydrogenmolecule wave function. For the three isotopologs, we have [89,90] [88]: c tot = 1.18. With these numbers, we obtain σ inel [T 2 ](E 0 ) = 3.64 × 10 −18 cm 2 , with an estimated uncertainty of 0.5 %. It must be noted that this theoretical cross section differs from the measured value, 3.40(7) × 10 −18 cm 2 [91], by 7 % (3.5 σ). However, it is ρdσ, directly measured by the e-gun as described in Sec. III B, which is used in the neutrino-mass analysisnot σ as a separate input.

C. Energy-loss function
Electrons traversing the WGTS can scatter elastically or inelastically from tritium molecules before being analyzed in the main spectrometer. (Here, "elastic" scattering refers to interactions that do not change the electronic state of the molecule.) While elastic scattering only causes a small broadening of the measured response function (∼0.03 eV), inelastic scattering can result in energy losses from ∼11 eV up to E/2, where the lower bound is associated with the lowest electronic excitations in T 2 .
Small inelastic energy losses, in particular, can move electrons emitted at energies close to the endpoint (the sensitive region for m 2 ν ) into a region still within the analysis interval extending 37 eV below the endpoint. Precise knowledge of the energy loss spectrum is, therefore, a crucial input for the KATRIN response function. During planning, its uncertainty was estimated to be one of the dominant systematics of the experiment [15]. A detailed paper on the energy-loss determination is in preparation [92].
Various electronic excitations, in combination with rotational and vibrational states of the T 2 molecule, result in a rich spectrum up to the ionization threshold at 15.486 eV [93]. Prior to this work, there were no calculations of the energy-loss spectrum with the required accuracy. We therefore measured the energy-loss function with the e-gun installed in the rear system of the KA-TRIN beamline. In contrast to β electrons originating within the source, these calibration electrons start with an adjustable kinetic energy chosen close to the endpoint of the tritium β spectrum and traverse the full length of the source. The dependence of the energy-loss function on the kinetic energy of the electrons can be neglected within the small fit window around the endpoint at ∼18.6 keV.
The e-gun uses a pulsed ultraviolet laser to create photoelectrons from a gold layer deposited onto the front face of an optical fiber. These electrons are then accelerated in an electric field with an adjustable angle to the local magnetic field lines. The electron energy is continuously scanned, in alternating directions, between 5 eV below and 55 eV above the main-spectrometer energy threshold qU .
The e-gun was operated in two different modes: a fast mode with a 100 kHz laser repetition rate to obtain a quasi-continuous electron beam used to record integral spectra as shown in Fig. 11 (top panel) and a slower mode with a 20 kHz repetition rate, in which the electron start times were synchronized with the DAQ to record time-offlight (TOF) spectra as shown in Fig. 11 (center panel).
This TOF information allows us to record a differential energy spectrum by applying a TOF cut on individual events [94]. Electrons with energies close to qU take significantly longer to reach the detector since they are decelerated to almost zero kinetic energy near the analyzing plane. Selecting electrons with flight times between 35 µs and 50 µs, as illustrated in Fig. 11 (bottom panel), effectively turns the main spectrometer from a high-pass filter into a narrow band-pass filter with a width of ∼0.02 eV. Apart from effects of multiple scattering and finite energy resolution, this method provides direct access to the electron energy-loss spectrum.
The energy-loss function is parametrized by a semiempirical model using three Gaussians to describe the three groups of lines created by excitations of the (2pσ 1 Σ + u ), (2pπ 1 Π u ) and (3pπ 1 Π u ) molecular states  around 12.6 eV [95] and the binary-encounter-dipole (BED) model [96] to describe the continuous ionization tail at energy losses above 15.5 eV (Fig. 11 center). The model has nine parameters given by the mean, width, and strength of each Gaussian. The normalization of the tail is chosen such that one obtains a smooth continuation of the Gaussian part of the model at the ionization energy.
To fit the measured TOF spectra, the model function is first convolved several times with itself, to account for multiple inelastic scatterings in the source, and then with the measured spectrum of unscattered electrons (peak at 0 eV in Fig. 11 center). This spectrum of electrons which have not undergone inelastic scattering naturally includes the effects of elastic scattering and the filter width of the main spectrometer. The resulting curves for single and multiple scattering are then weighted with the Poisson-distributed scattering probabilities and summed. The expectation value of this Poisson distribution is a nuisance parameter in the fit. A combined fit of TOF spectra taken at different column densities must also account for differences in the e-gun laser intensity between the individual measurements, leading to changes in the count rate. Additional normalization factors are therefore included as nuisance parameters in the fit. Finally, additional background components are included in the fit. Background electrons produced by the impact of positive ions onto the photocathode of the e-gun, for example, do not exhibit a TOF structure and appear in the differential spectrum as a small additional component with the shape of an integral energy-loss spectrum. The scaling factors of this background are additional nuisance parameters.
We performed a combined fit to four TOF datasets measured at different column densities. Each dataset contains about 12 hours of data, resulting in ∼6 × 10 5 events surviving the TOF cut. The nine model parameters of interest are shared between all datasets, whereas each dataset has its own nuisance parameters as described above.
The resulting best-fit parametrization is shown in Fig. 11 top and center for the integral and differential data, respectively. The same energy-loss function describes all four datasets well and the fit has a reduced χ 2 close to one.
Uncertainties used in this work are of a statistical nature only. However, more advanced combined fits that also take into account the integral energy-loss measurements yield the same parameter values within their statistical uncertainties.
Systematic uncertainties in the energy-loss determination are largely canceled by alternating up-and downward scans. A study of systematic effects on the parameter uncertainties has been undertaken using a Monte Carlo (MC) approach and taking into account disturbances like column-density drifts, background events, detector pileup and the binning of the continuous voltage ramp. These systematic uncertainties are negligible for the KNM1 analysis.
An improved parametrization of the energy-loss function and its uncertainties is under investigation for future, more sensitive neutrino-mass campaigns.

VI. BACKGROUND
The rate of background events during KNM1 was dominated by the two steady-state mechanisms described in Sec. VI A. In Sec. VI B, we also consider a background dependent on the duration of the corresponding scan step.

A. Steady-State Background
The steady-state background originates from excited or unstable neutral atoms which can propagate freely in the ultra-high-vacuum environment of the main spectrometer. It has two primary causes.
First, a significant part of the steady-state background arises from hydrogen Rydberg atoms sputtered from the inner spectrometer surfaces by 206 Pb recoil ions following α decays of 210 Po. These processes follow the decay chain of the long-lived 222 Rn progeny 210 Pb, which was surface-implanted from ambient air (activity ∼1 Bq/m 2 ) during the construction phase. A small fraction of these Rydberg atoms is ionized by black-body radiation when propagating through the magnetic flux tube. The resulting sub-eV scale electrons are accelerated to qU by the MAC-E-filter, adding a Poisson component to R bg .
The second significant steady-state background mechanism originates with α decays of single 219 Rn atoms (t 1/2 = 3.96 s) emanating from the non-evaporable-getter pumps. Each decay releases a large number of electrons up to the keV scale. If the decay occurs in the magnetic flux tube, these electrons are stored due to their significant transverse momenta. They subsequently produce secondary electrons by scattering on the residual gas until they have cooled to energies of a few eV, when they can escape; both primary and secondary electrons contribute to R bg at qU [97]. Since several background electrons may originate from each 219 Rn decay in the magnetic flux tube, this background source is not purely Poissonian. Liquid-nitrogen-cooled copper baffles at the ports to the getter pumps mitigate this effect by preventing 219 Rn from diffusing into the sensitive volume [39,98]. Due to the formation of a thin layer of H 2 O covering the baffle surface, the retention of 219 Rn was hampered such that R bg retains an observable non-Poissonian component during KNM1.
In KNM1, the overall steady-state background rate, R bg , is continuously measured through the energyindependent part of the spectrum R( qU ). The whole spectrum is fitted, leading to a value over the 117 selected pixels of R bg = 0.293(1) cps that is largely constrained by the 5 scan steps above the expected E 0 . This value is consistent with data from independent background runs. Full fit results are given in Sec. IX D.
The background is not distributed uniformly across the detector, as shown in Fig. 12. The decrease of R bg towards smaller radii can be explained by radiative deexcitation of the Rydberg atoms as they propagate inside the main spectrometer. Further from the spectrometer wall, fewer Rydberg atoms are therefore available for ionization by the thermal radiation.
The steady-state background was monitored for each β-scan with the five dedicated background-region scan steps. Figure 13 shows the time evolution of these background measurements during KNM1. A linear fit was applied to the data in order to test the long-term stability of the background. The slope of −0.01(8) mcps/d is compatible with a background that is stable over long time scales.
The non-Poissonian component of R bg causes a broadening of the event distribution of the five background-  region scan steps, amounting to 6.4 % compared to the prediction from pure Poisson statistics (Fig. 14). Our model predicts a background that is independent of qU near E 0 . To test this expectation, we performed a dedicated background-only measurement, without an active tritium source, in June 2018. As shown in Fig. 15, qU was scanned in 26 steps over an interval of 16.975 keV to 18.615 keV. We then fit a line with a free slope parameter to these data. The resulting best-fit slope, −2.2(43) mcps/keV, is compatible with zero, and we take its uncertainty as an overall uncertainty on our assumption of a qU -independent background (Sec. VIII G).  Background measurement without tritium source in a region near E0. The slope of a linear fit to the data is compatible with zero, supporting our assumption that the background is independent of qU . Fits to the immediate E0 region, where scan steps are more evenly spaced, also find no significant trend in qU .

B. Background Dependence on Scan-Step Duration
With both the pre-spectrometer and main spectrometer held at negative retarding potentials, a Penning trap inevitably forms in the strong magnetic field of the grounded inter-spectrometer region. Electrons trapped in this region slowly lose energy by ionizing residual gas molecules. The resulting ions may escape into the main spectrometer, where they can create background electrons when their own collisions with the residual gas or the vessel wall release ionization electrons, Rydberg atoms, or photons. The intense WGTS feeds the Penning trap when β electrons produce positive ions on their way into the pre-spectrometer; these ions sputter Rydberg atoms from the pre-spectrometer walls, and the Rydberg atoms in turn produce low-energy ionization electrons that fill the trap [40]. This mechanism may also play a role in main-spectrometer backgrounds, when β electrons scatter further downstream and the resulting ions strike the main-spectrometer walls.
During each transition to a new scan step, an electron catcher is briefly inserted into the beamline to remove stored electrons from the Penning trap. At higher pre-spectrometer potential, this has been shown to provide a statistically significant reduction in the baseline background [40]. However, since the electron catcher is inserted only at the beginning of a scan step, the Penning trap continues to fill until a new electron-catcher actuation at the beginning of the next scan step. The corresponding rise of the background rate is strongly influenced by surface conditions and by the achieved pressure between the spectrometers. In principle, however, this mechanism can produce a background that effectively increases in rate for longer-duration scan steps (see measurement-time distribution in Fig. 7). This effect was observed in a subsequent KATRIN scientific run, but for KNM1 -the initial science run, with pristine surfaces and lower column density -no statistically significant dependence on scan-step duration was observed. Section VIII G will address the impact on the neutrino-mass measurement.

VII. ASSEMBLING SPECTRAL DATA FOR KNM1
Data are acquired in a sequence of O(2 h) scans and the integral spectrum (Eq. 12) is recorded with the FPD. In the final analysis (Sec. IX), the spectral fit uses four free parameters: the signal amplitude A s , the effective β-decay endpoint E 0 , the background rate R bg , and the squared neutrino mass m 2 ν . In this analysis we leave E 0 and A s unconstrained, which is equivalent to a "shapeonly" fit. The 4-parameter fit procedure over the averaged scan steps qU compares the experimental spectrum R( qU ) to the model R model ( qU ).
Spectra from all of the scans and pixels have to be combined in the final analysis without loss of information. In the following we describe the strategy applied to combine all these data prior to the final spectral fit to extract the effective neutrino mass.

A. Pixel combination
During KNM1, the electric potential and magnetic field in the analyzing plane of the main spectrometer were not perfectly homogeneous, but varied radially by about 140 mV and 2 µT, respectively, and to a much smaller extent azimuthally. The pixelation of the detector allows us to account for these spatial dependencies. Each pixel has a specific transmission function and records a statistically independent tritium β-electron spectrum. In this analysis, we combine these pixel-wise spectra into a single effective pixel by adding all counts and assuming an average transmission function for the entire detector. The averaging of fields leads to a negligible broadening of the spectrum which does not affect the filter width, and carries a negligible bias of O(10 −3 eV) on m 2 ν . Combining all 274 scans that passed data-quality cuts, single-pixel fits were performed resulting in an endpoint E fit 0 for each pixel, as shown in Fig. 16. We find no systematic spatial (i.e. pixel) dependence of E fit 0 . The standard deviation from the mean endpoint is 0.16 eV, which is consistent with statistical fluctuations. This indicates a good description of the electric potential and magnetic field in the analyzing plane, and the absence of a significantly spatially dependent electron starting potential. We therefore merge the data of all 117 selected pixels used in the analysis (Fig. 5).

B. Scan combination (stacking)
Combining all pixels in a uniform fit, we can now consider the stability of the fit parameters with respect to possible temporal variations. We investigate all four free parameters in the fit. For single scans of 2 hours, the accumulated statistics are not sufficient to significantly constrain the neutrino mass. Therefore, the neutrino mass is fixed to zero. The 274 fit values show excellent stability over the course of a month (Fig. 17). The standard deviation from the mean endpoint is 0.25 eV, which is again consistent with statistical fluctuations. In order to constrain the neutrino mass, the statistics of all 274 scans must be combined. Based on our stability results, we achieve this by merging the data of all 274 scans into a single stacked, integral spectrum. In the underlying process, the events at like scan steps are summed and the corresponding retarding-potential values are averaged over all scans. This procedure yields one high-statistics integral spectrum with the same number of scan steps as a single scan. Since this method does not correct for scan-to-scan variations of slow-control parameters, it relies on good time stability and excellent reproducibility of the individual HV settings from scan to scan. The Gaussian spread of these HV settings is on average σ = 34(1) mV (better than 2 ppm) (Sec. III D). The scan stacking results in a minor systematic effect, which is included in the analysis.

VIII. SYSTEMATIC UNCERTAINTIES
Systematic uncertainties generally arise from parameter uncertainties that enter into the calculation of the integral spectrum, and from instabilities of experimental parameters. The KNM1 analysis heavily relies on a precise description of the spectral shape, including all relevant systematic effects and a robust treatment of their  uncertainties. Any erroneously neglected effect or uncertainty can lead to a systematic shift of the deduced neutrino mass [99]. The individual systematics are described in detail below. A summary of these systematic uncertainties is given in Table III, while their ultimate impacts on the m 2 ν uncertainty budget are collated in Table IV.

A. Tritium concentration
The concentration of the tritium isotopologs in the source affects the model in two different ways.
First, the total activity is directly correlated to the tritium purity described in Eq. (4). The absolute number does not impact the neutrino-mass measurement, as the signal normalization is a free fit parameter. Changes during a given scan, however, could introduce a slight spectral distortion which would bias the measurement. As described in Sec. III A, the tritium purity was measured continuously by the LARA laser-Raman spectroscopic system. The precision was determined from the shot noise √ N x of the Raman signal and then propagated to ε T and c x ; the resulting precision of better than 2 × 10 −3 for each scan was reported in Ref. [21]. Scan-to-scan fluctuations of the tritium purity amount to 0.39 × 10 −3 after accounting for anti-correlations between the isotopologs.
Second, each of the three tritium isotopologs also has a slightly different FSD. Systematic uncertainties on their relative fractions, mainly determined by the trueness of the LARA calibration, thus propagate into the spectral shape. The impact on m 2 ν from this effect is less than 2 × 10 −4 eV 2 and is thus negligible for KNM1.

B. Column density and expected number of scatterings
The determination of the expected number of scatterings, ρdσ, is described in Sec. III B. The total uncertainty on ρdσ arises from three separate contributions: the limited precision of single column-density measurements made with the e-gun; uncertainty on the throughput measurement, arising from fluctuations of the gas throughput and imperfect reproducibility of the flow meter; and the scaling of the inelastic-scattering cross section to a lower electron energy via Eq. 17. This last operation is necessary because the e-gun is operated at an energy of 18.78 keV, well above E 0 , for measurements of the column density -but the β electrons, at lower energies, have a slightly different scattering cross section. We take 18.575 keV as a representative value for our observed β electrons; the variation of the inelastic-scattering cross section within the analysis interval is negligible.
Taking these three contributions into account leads to a total systematic uncertainty on ρdσ of less than 0.85 % for all scan steps.

C. Electron starting potential
Spatial inhomogeneities and temporal fluctuations of the starting potentials of the β electrons would lead to a shift of the neutrino mass according to Eq. 7. As discussed in Sec. III C, the intrinsic width of the 83m Kr L 3 -32 line is a diagnostic tool to investigate these effects, probing the plasma-potential distribution.
In the KNM1 analysis, we treat the fitted Gaussian line broadening in the presence of a T 2 plasma as a conservative upper limit for the inhomogeneity of the plasma potential, yielding a negative m 2 ν shift with magnitude less than 0.013 eV 2 .
Electrons undergo inelastic scattering as described in Sec. V. The s-fold scattering probabilities for each β electron depend on the longitudinal position of its creation. As a result, the populations of β electrons with different scattering multiplicities also have different distributions of starting positions, and therefore different distributions of starting potentials if the plasma potential is inhomogeneous. Analysis of the positions of the krypton L 3 -32 lines of unscattered and singly scattered electrons shows that a plasma-induced mutual shift of these positions cannot be larger than 70 mV. The corresponding additional m 2 ν -shift can be neglected for KNM1. We thus conclude that the effective L 3 -32 broadening parameters given above serve as a very conservative upper limit of plasma effects in the neutrino-mass analysis.
In addition to the 83m Kr spectroscopy method, radial plasma inhomogeneities can be inferred directly from the neutrino-mass data by radial evaluation of E 0 . The spectral fit from twelve separate detector rings (see Fig. 5 for detector structure) revealed a slope of −2(5) mV/ring, consistent with a slope of zero.
A full propagation of the plasma model and its uncertainty was not included in the KNM1 analysis, primarily due to the immaturity of the plasma model as applied to the low KNM1 column density. Adding this O(−0.01 eV 2 ) uncertainty in quadrature to the total systematic uncertainty does not yield significant leverage on the total budget.
The neutral-gas density strongly affects the charge densities from secondary electrons and ions, as well as other plasma parameters. For this reason, we are currently investigating the effect of different column densities, gas temperatures, source magnetic-field strengths, and changing boundary conditions on plasma parameters. This will inform the consideration of plasma effects in the data analysis for upcoming campaigns, in which the gas throughput will be higher by a factor of up to four.

D. Detector efficiency
Although numerous physical and detector effects can reduce the detector efficiency, any effects which do not depend on the retarding potential U will not affect the KNM1 fit results due to the overall, free scaling parameter for each spectrum and the uniform, all-pixel fit.
The overall FPD detection efficiency within the ROI has been estimated by both simulation and commissioning analysis to be approximately 95 %, with an uncertainty of a few percent, and per-pixel variations of about the same size. For KNM1, the ROI is fixed regardless of U (Sec. III E). However, the shape of the FPD energy spectrum changes with U , primarily due to the β-electron energy threshold at qU . Additional distortions are due to energy-or rate-dependent detector effects: energy loss in the dead layer, charge sharing among pixels, pileup, and back-scattering of electrons and their subsequent reflection back toward the FPD by local electric and magnetic fields.
We have studied the effects of these spectral shape changes using a reference spectrum for each pixel, acquired at U 0 = −18 375 V. For each scan step at U i = U 0 + ∆U , the reference spectrum is shifted by the corresponding q∆U and a count correction is calculated. As |U | decreases, the corrections become larger, with a maximum size of about 0.05 %. We estimate the error relative to these correction factors at less than 0.05 %, determined by comparing spectral shapes at nearby U values. In the KNM1 analysis, we apply these corrections to FPD counts while neglecting the corresponding uncertainty.
Pileup events also result in event loss, since the energy is erroneously reconstructed above the upper bound of the ROI. We assume that pileup events arise from random coincidences; each coincidence produces a total energy deposit that is an integer multiple of 28.6 keV, within the shaping time of the trapezoidal filter. We calculate and apply the corresponding correction to the event rate for each pixel and scan step, up to a maximum correction factor of 0.02 % at low |U | and, correspondingly, high rate. Our conservative estimate of the relative error on these correction factors is less than 18 %, based on the shape of the measured FPD energy spectrum and a simulation of the trapezoidal filter. This error is negligible.
Our final consideration is electron backscattering from the FPD. The majority of backscattered electrons are reflected back to the FPD, either by magnetic fields in the detector system, or by the electric potentials of the postacceleration electrode or main spectrometer. Even with multiple backscatters, the electron returns to the same pixel each time, always arriving well within the shaping time of the trapezoidal filter, so that the detector does not register the event as separate hits. Our spectralshape calculations include the resulting reconstructedenergy shifts, due to multiple transits of the detector dead layer and hits distributed within the shaping time. However, an additional correction is in principle needed for those few backscattered electrons which have enough energy to surmount the qU threshold and escape towards the source. Simulations show that the resulting event loss is less than 0.01 % for the KNM1 analysis window. This effect is therefore neglected in this analysis.

E. Final-state distribution
The uncertainty estimation on the FSD is based on differences between the theoretical ab initio calculations from Saenz et al. [78] and Fackler et al [100]. The differ-ence between the calculations for the ground-state variance is found to be small, of O(1 %) [101]. However, the descriptions of the electronic excited states and the electronic continuum exhibit larger discrepancies.
We conservatively estimate the uncertainty on the variance of the ground state (excited states and continuum) to be 1 % (4 %). The uncertainty on the normalization of the ground to excited-state populations is taken as 1 %.
Our narrow analysis interval, extending 37 eV below E 0 , is dominated by electrons from the ground-state distribution. Consequently, the uncertainty on the FSD only contributes on the order of O(10 −2 ) eV 2 to the total systematics budget on m 2 ν within our analysis interval.

F. Response function
Response-function-related systematic uncertainties are connected with the electromagnetic fields that define the transmission function (Eq. 16) and with the energy-loss function. The electromagnetic fields are computed from a simulation of the beamline magnets and the mainspectrometer vessel.
a. Magnetic fields Systematic uncertainties on the magnetic field at the analyzing plane arise from residual magnetic fields in the spectrometer hall, e.g. due to magnetized materials, and from model imperfections. A sensor network was used to compare measured fields at the spectrometer vessel to simulation results. Our assessment of the maximum deviation yields a conservative systematic uncertainty of ∆B min /B min = 1 %.
The maximum magnetic field, located at the exit of the main spectrometer, was measured in 2015 at the center of the magnet bore [30] and compared to simulations. We include a conservative systematic uncertainty of ∆B max /B max = 0.2 %.
The source magnetic field was measured in 2009 by the manufacturer with Hall probes on the central axis and compared to simulations. We include a conservative systematic uncertainty of ∆B S /B S = 2.5 %.
b. Electric potentials Since any offset of the simulated retarding potential at the analyzing plane is compensated by the free endpoint parameter, no additional systematic uncertainty is assigned for the spectral fit.
c. Energy-loss function The uncertainty of the energy-loss parametrization is obtained from fits to the measurements described in Sec. V C. For each of the 9 parameters describing the energy-loss function, an individual fit uncertainty is determined. As stated in Sec. V C, the contribution of systematic effects is about one order of magnitude lower than the uncertainties related to the current statistics of the e-gun measurements. As a result, only statistical fit uncertainties are considered for this analysis. Correlations between the energy-loss parameters are taken into account, reducing the overall uncertainty of the energy-loss function with respect to the uncorrelated case.
The systematic effect on m 2 ν due to the uncertainties of the energy-loss function is determined to be below 0.01 eV 2 .

G. Background
The steady-state background enters the uncertainty budget in two independent ways: rate and shape.
The background rate distribution, as shown in Fig. 14 shows an over-dispersion of 6.4 % compared to the Poisson expectation. This enters the analysis as an additional uncorrelated uncertainty on the background rate, effectively increasing the statistical error in the region with E > E 0 − 15 eV.
As described in Sec. VI A, we expect the background to be flat with respect to the retarding potential. In this analysis we assess the slope uncertainty via a slope parameter, which makes a first-order correction to the constant expectation. Based on the dedicated measurement in June 2018 (Fig. 15), the slope parameter is consistent with zero, within an uncertainty of 5 mcps/keV. In the final spectral fit (Sec. IX), we use a central value of 0 mcps/keV.
A Penning-induced background (Sec. VI B) may increase over the course of each scan step, effectively introducing a higher background for scan steps with longer duration. Since longer scan steps are concentrated near E 0 − 14 eV (Sec. III G), the net effect is a shape distortion of the background shape. An analysis of KNM1 scan steps yields a best-fit linear time slope of (−3.8 ± 4.4) µcps/s, which would result in a systematic uncertainty of 0.15 eV 2 on the squared neutrino mass. This systematic was not taken into account in the spectral fit (Sec. IX), but would not alter the statisticsdominated final uncertainty.

H. Stacking
The averaging of the scan steps within the stacking techniques introduces a small bias on m 2 ν and E 0 . In order to quantify these biases, we construct an Asimov dataset [102] by simulating 274 statistically unfluctuated "MC twin" spectra, incorporating the actual variation of slow-control parameters (including measured highvoltage values, isotopic compositions, and column densities) between scans. Later on, the MC spectra are combined into a single integral spectrum through the stacking procedure, as described in Sec. VII. As a last step, we fit this stacked MC spectrum. Comparing this fit result to the MC truth yielded a 1σ stacking uncertainty of 14 × 10 −2 eV 2 in one analysis approach (Sec. IX B), and 5 × 10 −2 eV 2 in the other (Sec. IX C), as shown in Table IV further below. The discrepancy between the two approaches arises from different treatments of the individual contributions to this subdominant uncertainty; the stacking method and error treatment will be optimized in the analysis of future neutrino-mass campaigns, in which scan-to-scan fluctuations are also expected to be smaller.

I. Neutrino-mass fit range
The full spectrum was recorded over a large energy range down to E 0 − 91 eV. Several systematic uncertainties, like those related to inelastic scattering and the FSD, increase further away from the endpoint, while the statistical uncertainty decreases. The optimization of the neutrino-mass fit range is performed using MC twin simulations of KNM1 (Sec. VIII H), assuming a zero neutrino mass and using the set of systematics presented earlier in the section (Table III). The lower bound of the fit interval is then varied between E 0 − 91 eV and E 0 − 30 eV, and two fits are performed in turn. The first fit considers statistical uncertainty only, while the second fit uses both statistical and systematic errors. For each pair of fits, the systematic uncertainty is deduced by subtracting the statistical uncertainty in quadrature from the total error. As a result, both statistical and systematic uncertainties become equal for the fit range starting at about E 0 − 70 eV, and systematic uncertainties become dominant when including data below E 0 − 70 eV. Moreover, the overall sensitivity only marginally improves by including data at energies below E 0 − 40 eV.
This study addresses only the dependence of the measurement precision on the fit range. It does not address the accuracy of the determination of the neutrino mass, since the same model is used for the fit and for the MC twins. Indeed, further than about E 0 − 40 eV, the electronic continuum -with less well-validated modelingdominates the FSD (Sec. IV B). Therefore, before unblinding the data (Sec. IX A, below), we fixed the analysis interval to cover the region of E 0 − 37 eV (22 scan steps) and E 0 + 49 eV (5 scan steps).

IX. SPECTRAL FIT
In this section we discuss our blinding method (Sec. IX A) and present two approaches for inferring the value of the neutrino mass squared m 2 ν and the endpoint E 0 simultaneously, based on fitting the integrated β spectrum (Eq. 12) assembled as described in Sec. VII. In both approaches, the spectrum is fitted using a shapeonly analysis with four free parameters. In addition to m 2 ν and E 0 , these are the signal amplitude A s and the background rate R bg .
The first approach (Sec. IX B) uses a standard χ 2 estimator and covariance matrices to encode all uncertainties. The second approach (Sec. IX C), Monte-Carlo propagation, repeats the final fits many times, for each fit choosing randomized input values for the systematic nuisance parameters.
Three analyses were performed, each with its own spectrum calculation and analysis software: two us-ing the covariance-matrix approach, and one using the MC-propagation approach. The analyses were performed blind and give consistent results, as described in Sec. IX D. The resulting breakdown of systematic uncertainties is given in Table IV, below. Section X uses these spectral results to derive frequentist bounds on the neutrino mass, while Sec. XI uses the same data to derive Bayesian bounds.

A. Blinding strategy
For the KNM1 analysis we enforced blind analysis procedures to fix data selection, analysis cuts, and model composition before the model was fitted to the data. This standard technique is designed to avoid observer's bias.
For this first KATRIN m 2 ν limit, we employed model blinding rather than data blinding. The fit results are highly dependent on the molecular FSD (Sec. IV B); in particular, the value of m 2 ν depends on the width of the distribution of transitions to the electronic ground state of the daughter molecule 3 HeT + . Using an FSD with too large a width pushes m 2 ν towards higher values, while too narrow a width pushes it towards lower values. Indeed, historically, inaccurate FSD models were likely responsible for artificially negative m 2 ν results from the Los Alamos [18] and Livermore [103] experiments, a problem which is resolved by using the more modern theory described in Sec. IV B [101].
If we fit the data with a model using an FSD groundstate width that has been picked randomly within a suitable interval, the true value of m 2 ν cannot be retrieved. That is, the analysis is blind to its parameter of interest, while the remaining three parameters are left essentially unaffected [64]. The range of possible ground-state widths was chosen so that the sensitivity of the KATRIN blind analysis could not improve upon the results of previous direct m 2 ν measurements [11,12]. In addition, because the endpoint fit parameter only depends -to a good approximation -on the mean of the FSD, leaving that mean value untouched ensured that the endpoint could still be used during a blind analysis, e.g. for comparison with other independent measurements (Sec. XII).
In practice, the theoretical electronic ground-state manifold of the FSD was swapped with a Gaussian distribution function, constructed with the true mean and a randomly chosen width. To prevent accidental unblinding, the adjusted FSD was provided as an independent software module synchronized with the main fitting software.
The second measure to mitigate biasing is to perform the full analysis, including parameter fitting, using MCbased data sets first, before turning to the experimental data. For each experimental scan i we generate an MC twin (Sec. VIII H) from its averaged slow-control parameters to calculate the expected rate R β (E) i with the corresponding response function f (E − qU ) i and background rate R bg,i . Analyzing the MC twins allows us to verify the accuracy of our parameter inference by recovering the correct input MC values for m 2 ν . This MC dataset is used to assess statistical (σ stat ) and systematic (σ syst ) uncertainties and to compute our expected sensitivity. It is also used to benchmark the independent analysis codes. At this stage, all model inputs and systematic uncertainties are frozen.
Before the unblinding via incorporation of the unmodified FSD, a final benchmark was successfully performed on the data with the blinded FSD to verify that the independent analysis codes eventually lead to very consistent results. After this final test, the "true" FSD was revealed to the collaboration for the final neutrino-mass analysis of the data. The first, overnight fits -using the independent analysis codes -already yielded preliminary, consistent results the very next morning.

B. Covariance-matrix approach
Here, we report on our results using the covariancematrix approach to include and propagate systematic uncertainties in the neutrino-mass fit. The spectrum calculation code and methods used for this analysis are described in detail in Ref. [104].
The free fit parameters in our analysis, θ, are inferred from the data points {R i } by minimizing the negative logarithm of the ratio of the Poisson likelihood function to the saturated model where the summation is over scan steps i. The model points, denoted by R model i , depend on both the model parameters θ and the systematic nuisance parameters η (including column density and tritium isotopolog concentrations). In the fit the nuisance terms η are fixed according to our best knowledge of operational parameters averaged over KNM1.
Since the β spectrum measured in this first KATRIN science run comprises a large number of observed events in each scan-step bin, the negative Poisson likelihood function (Eq. 18) is replaced by the standard χ 2 estimator The covariance matrix C describes the correlated and uncorrelated model uncertainties, including both statistical and systematic uncertainties. This fit procedure has been extensively tested by injecting fake neutrino-mass signals in simulated pseudo-experiments. It was verified that the fit results provide an unbiased estimation of the injected parameters.
Systematic uncertainties on the nuisance parameters η are propagated using covariance matrices. For this purpose the values of η are randomized according to their associated probability density functions. Correlations between parameters are taken into account. Subsequently, O(10 4 ) sample spectra {R sample } are simulated [20,105,106]. For each sample-spectrum calculation, a different η is drawn from the set {η sample }.
The signal normalization A s , being a free fit parameter, is not considered in the uncertainty propagation. Therefore, all fluctuations in {R sample } that translate solely into an overall signal normalization uncertainty must be eliminated. The transformation of {R sample } into shapeonly sample spectra is achieved by normalizing the statistics of each sample spectrum to the statistics of the average sample spectrum.
Finally, the shape-only covariance matrix is estimated from {R sample } using the sample covariance as an estimator. For any set of uncorrelated systematic effects, the associated covariance matrices can be calculated independently of one another. The sum of all matrices encodes the total uncertainties on the model points R model and their scan-step-dependent correlations.
In the fit, χ 2 (θ) is minimized to determine the bestfit parametersθ, whereas the profile of the χ 2 function is used to infer the uncertainties onθ. Once the covariance matrices are pre-calculated, the spectral fit and major diagnoses can be performed within a few hours on a standard personal computer.
The data and results of this fit are displayed in Fig. 19. Of the four free parameters, the signal amplitude A s is unconstrained for the shape-only analysis. The effective β-decay endpoint E 0 can be related to the Q-value after final corrections of the energy scale (Sec. XII). The background rate R bg is primarily constrained by the 5 HV scan steps above E 0 . The squared neutrino mass m 2 ν can be varied freely and therefore can take any positive or negative value.
We find a best-fit value of m 2 ν = (−0.98 + 0.89 − 1.06 ) eV 2 with a goodness of fit of χ 2 = 21.4 for 23 degrees of freedom (d.o.f.). This corresponds to a p-value of 0.56, meaning that there is a probability of 56 % to retrieve a χ 2 -value at least as large as the one obtained.
The total uncertainty budget of m 2 ν is first calculated on an Asimov data set assuming the null hypothesis. Based on the final fit applied to these simulated data, we derive m 2 ν = 0.00 +0.78 −0.94 eV 2 . The relative impact of each systematic effect is assessed by performing a series of fits, each one including solely the selected effect in addition to statistical uncertainties (stat+1 test). The statistical uncertainty is then subtracted in quadrature. The same breakdown is then calculated using the unblinded data, and is in excellent agreement with our MC expectations. This data-driven uncertainty breakdown is shown in Table IV. As expected, the total uncertainty is largely dominated by σ stat (0.94 eV 2 ) as compared to σ syst (0.30 eV 2 ).

C. Monte-Carlo-propagation approach
Here we report the fit results using the MCpropagation approach to propagate systematic uncertainties. The spectrum-calculation code used is described in Ref. [107] while the method is adapted from Refs. [108,109].
In the MC-propagation method, we repeat the fitting process ∼ 10 4 times, each time with newly randomized input values for the systematic nuisance parameters η that are held fixed during that fit. Compared to the wellknown approach of free nuisance parameters constrained with pull terms, this method has two key advantages for the KATRIN analysis. Foremost, the computationally expensive response function does not have to be recomputed with varying η during the fit. In addition, the minimization is technically simplified due to the reduced number of free parameters.
To retrieve an initial estimate of the best-fit valuesθ data of our four fit parameters θ (that is, m 2 ν , E 0 , A s , R bg ), we fit the original data with the additional parameters η fixed to our best knowledge from the experiment. Next, we generate MC spectra assuming the valuesθ data for our model and a Poisson distribution of the counts. We then fit each of these statistically random- ized MC spectra, retrieving one sample of valuesθ stat sample for our free parameters. The resulting distribution of {θ stat sample } can be used to infer the statistical uncertainty of θ.
Our next step is to assess the systematic uncertainties, beginning by varying the values of η according to their uncertainties. The model is initialized with the random values η sample . We then fit the randomized model to our reference spectrum, which assumes the best estimate for η andθ data . In principle, the resulting distribution of {θ syst sample } reflects the systematic uncertainty, taking into account only the external information on η. However, the data may also contain information to constrain η. To account for this, we also fit the randomized model to the data to retrieve the likelihood value L(θ syst sample ). This likelihood value is used to weight the corresponding sampleθ syst sample . The resulting weighted distribution {θ syst sample } weight is then used to retrieve the systematic uncertainty on θ as proposed in Ref. [110]. At this point we would like to note that this systematics-only distribution is solely used to calculate a breakdown of the uncertainties and does not enter into the final confidence interval.
In the final step, we combine the statistics-and systematics-only steps described above.
As in the systematics-only approach, we initialize our model with randomized values for the nuisance parameters η sample . Instead of fitting it to the unfluctuated best estimate, we now fit this model to statistically randomized spectra to retrieve the valuesθ tot sample of our parameters of interest. This model is then also fit to the unmodified data spectrum to retrieve the likelihood L(θ tot sample ). We infer the combined statistical and systematical uncertainty from the distribution of {θ tot sample } weight , which is weighted by these likelihood values.
Initially, we apply this method to the MC twin data described in Sec. VIII H). From the statistics-only fit, we derive m 2 ν = 0.00 +0.75 −0.90 eV 2 . Including the systematic uncertainties described in Sec. VIII, the best-fit value becomes m 2 ν = 0.00 +0.76 −0.96 eV 2 . This is only a slight change with respect to the statistics-only analysis.
After freezing the method and inputs on MC spectra, we repeat the analysis on the data. Here the statistics-only fit to the data gives a best-fit value of Using the MC propagation of uncertainty, it is possible to analyze the impact of individual systematic effects on the parameters of interest. Table IV, further above, shows the uncertainty budget on m 2 ν for KNM1.

D. Fit results
The results of the two independent methods of Secs. IX B and IX C agree to within a few percent of the total uncertainty. As a best-fit value for the squared neutrino mass, we quote m 2 ν = −1.0 +0.9 −1.1 eV 2 . This bestfit result corresponds to a 1 σ statistical fluctuation to negative values of m 2 ν . Assuming the true neutrino mass is zero, the probability to retrieve a best-fit value as negative as ours is 16 % and is thus fully compatible with statistical expectations. The total uncertainty budget of m 2 ν is largely dominated by σ stat (0.97 eV 2 ) as compared to σ syst (0.32 eV 2 ). The dominant contributions to σ syst are found to be the non-Poissonian background from radon and the uncertainty on the background slope. Uncertainties on the column density, energy-loss function, FSD, and magnetic fields play a minor role in the budget of σ syst . Likewise, the uncertainties induced by fluctuations of ε T and HV parameters during a scan are negligibly small compared to σ stat .
For the effective β-decay endpoint we find a best fit value of 18 573.7(1) eV. Figure 21 shows the interplay between m 2 ν and E 0 . The large correlation (0.97) between the two parameters is in line with expectation [3,99].
For completeness, we report here that our best-fit background rate is R bg = 293(1) mcps. The signalnormalization parameter A s absorbs the rate effects of our systematic uncertainties, and does not have a straightforward interpretation.

X. FREQUENTIST BOUNDS ON THE NEUTRINO MASS
The result of a neutrino-mass experiment is commonly presented in form of a confidence interval for the neutrino mass, or an upper limit if the lower boundary of the confidence interval is zero. These values are used by the community for constraining phenomenological models, developing theoretical predictions, and comparing the results of different experiments, and as input parameters to both terrestrial experiments and cosmological observations.
There are several methods of constructing the confi-dence intervals with additional information on the estimated parameter. To account for the physical bound of m 2 ν ≥ 0, despite the fact that m 2 ν is unconstrained in the fit, we perform full Neyman constructions using the methods of Lokhov and Tkachov and of Feldman and Cousins, for completeness. Both of these methods avoid empty confidence intervals for negative best-fit estimates m 2 ν . In each case, we apply both of our spectral analysis approaches (described in Sec. IX B and Sec. IX C) to incorporate statistical and systematic uncertainties into the calculated Monte Carlo quantities. This results in two calculations of each type of confidence interval, which agree with each other in both cases. We briefly compare the Feldman-Cousins and Lokhov-Tkachov methods below.
In the Feldman-Cousins method [22], the likelihood ratio determines the order in which the estimates m 2 ν are added to the acceptance region for an assumed value of m 2 ν , thereby constructing the confidence interval. This ordering principle avoids empty intervals, but at the same time results in more stringent limits for negative best-fit estimates that are further from zero, as in Fig. 22a. This yields an excessively strict upper limit in the case of statistical fluctuations in one direction, or in the presence of an unknown systematic bias as seen in most neutrinomass experiments of the early 1990s (see Fig. 26). While our best-fit result is statistically compatible with zero, we decided after unblinding to pursue an alternative approach to ensure a conservative handling of fluctuations.
Following the prescription of Lokhov and Tkachov [23], a new estimator m 2 ν can be defined such that The estimator is by definition as close as possible to the unknown true non-negative value of the m 2 ν , which is the fundamental aim of the statistical estimation. The confidence interval for the new estimator m 2 ν is then constructed according to the Neyman procedure, which guarantees the correct coverage. The non-physical values of the best-fit estimate m 2 ν are indistinguishable and give the same confidence interval from zero to the experimental sensitivity (Fig. 22b). Therefore more negative values of m 2 ν , obtained due to a statistical fluctuation or an improperly treated systematic contribution, do not yield better upper limits. This makes it possible to compare the upper limits of different measurements directly without the need to know the best-fit estimate, as long as m 2 ν is not significantly positive.
In order to allow the squared-neutrino-mass estimator to become negative in either analysis, the differential spectrum shape must be extended into the unphysical region of m 2 ν < 0. In previous experiments [11,12] the extension was made by modifying the differential spectrum shape so that the χ 2 function became symmetric around m 2 ν = 0. Such a modification depends on the particular shape of the χ 2 function and consequently on the experimental setup. In the present analysis we take the differential spectrum shape in Eq. (8) without any modification for m 2 ν < 0. This leads to a χ 2 function with an asymmetric shape, as shown in Fig. 20. The Lokhov-Tkachov method yields the same upper limit for all m 2 ν < 0. Therefore, by construction, the upper limit does not depend on a particular choice of the extension.
Using the Lokhov-Tkachov construction we derive an upper limit of m ν < 1.1 eV (90 % CL) as the central result of this work. For comparison, the Feldman-Cousins method yields the upper limit m ν < 0.8 eV (90 % CL). We have also derived upper limits at 95 % CL for comparison to the Mainz [11] and Troitsk [12] Feldman-Cousins results. In the Lokhov-Tkachov method, this becomes m ν < 1.2 eV (95 % CL); using Feldman-Cousins, as was done by Mainz and Troitsk, we find m ν < 0.9 eV (95 % CL).

XI. BAYESIAN BOUND ON THE NEUTRINO MASS
Bayesian analysis methods provide an alternative means of handling the unphysical, m 2 ν < 0 region. We used the MC-propagation model and data framework, described in Sec. IX C, to set a first limit using Bayesian techniques. Posterior probability distributions were constructed according to Bayes' theorem, using Markovchain Monte Carlo methods within the Bayesian Analysis Toolkit (BAT) [111]. We use uniform priors, flat in probability, for A s , E 0 , and R bg ; this choice is most straightforward for analysis of stacked spectra. An informative prior, restricting the result to only physically allowed m 2 ν values (equal to or larger than zero), is used to ultimately obtain an upper credibility limit on the neutrino mass in a Bayesian interpretation. In the allowed region, this prior is flat in m 2 ν space. Future work will investigate alternate choices of prior, including a prior flat in m ν .
First, we extract statistical uncertainties and compare with other analysis methods using the basic model, including the four-parameter set θ with flat prior probabilities. The global mode (maximum value) of the 4dimensional posterior for m 2 ν is found at −1.0 eV 2 . The two-sided 1σ interval, with equal probability on either side, is obtained from the posterior distribution marginalized for m 2 ν as [−2.1, −0.3] eV 2 . Four of the leading systematic uncertainties are included in this analysis, and are incorporated into the fit in various ways. A background slope is included as a fifth free parameter with a Gaussian prior probability centered around zero and a width given by its uncertainty. Non-Poissonian background counts are included by widening the underlying likelihood distribution in each scan step according to background measurements (Sec. VI). Variations of the response due to uncertainties in the magnetic field or the column density were too computationally ex- pensive at the time of the analysis. Instead of including these as free parameters in the model, multiple independent fits were parallelized on a computing cluster. Each fit was started with the input systematic fixed at a different value, following a Gaussian distribution with a width given by the parameter uncertainty. The median values of the output posterior distributions were used to obtain parameter estimates with systematic uncertainties. The same results are obtained by combining the Markov chains of the individual fits into a single chain, and subsequently performing the same parameter-estimation procedure. Additional systematics will be analyzed in future work.
The present dataset is strongly dominated by statistical uncertainties, and individual systematic effects are largely masked below 0.1 eV 2 by numerical uncertainties. These uncertainties come from the finite number of Markov-chain Monte-Carlo samples and are on the order of 0.006 eV 2 in the 1σ posterior width. Hence, the systematic budget was investigated with Asimov data, artificially increasing the amount of data and thus enhancing each included systematic effect with respect to statistical uncertainties.
Taking these four explicitly included systematic uncertainties into account, the most probable m 2 ν value was found at −1.0 eV 2 and the two-sided, 1σ, probabilitysymmetric interval at [−2.2, −0.3] eV 2 . Using Table IV to estimate upper bounds on the primary excluded systematics -scan fluctuations and the FSD -we find that they affect the total uncertainty on this most probable value by about 1%.
To determine the limit on the neutrino-mass, we then perform the same fits with a flat prior in m 2 ν ≥ 0. The m 2 ν marginalized posterior distribution is shown in Fig. 23. The best-fit value is found at m 2 ν = 0. The 90 % quantile of the marginalized posterior distribution is at 0.78 eV 2 . The Bayesian upper limit is thus m ν < 0.9 eV (90% C.I.). The constant prior probability in m 2 ν -space gives equal probability for statistical fluctuations in the data. In our case, the Bayesian 90% credibility limit is numerically closer than the Feldman-Cousins 90% confidence limit to the sensitivity of the experiment and to the Lokhov-Tkachov limit, as is often observed in the presence of larger statistical fluctuations.
As an additional test, the positive flat prior was slightly modified by knowledge from oscillation experiments, allowing only m ν > 8 meV (normal ordering) or m ν > 50 meV (inverted ordering) [13]. The posterior quantiles show no numerical difference, as is expected with the current data.

XII. Q-VALUE MEASUREMENT
A consistency check of the energy scale of KATRIN can be performed by extracting the experimental Q-value for molecular tritium from KATRIN data, and comparing it to Q-values based on Penning-trap measurements of the 3 He-T atomic mass difference. The Q-value represents the amount of kinetic energy released in β decay for zero neutrino mass; Fig. 24 shows its relationship to the mass difference and the binding energies of the atomic and molecular states involved in T 2 β decay. In equation form, we have:  The KATRIN result for the Q-value in molecular tritium β decay is derived from the best-fit value of E 0 with corrections for the center-of-mass molecular recoil of the 3 HeT + daughter ion, as well as the relative offset of the electron starting potential in the source to the work function of the inner electrode of the main spectrometer.
For the effective endpoint, our two fitting methods both obtain a best-fit value of E 0 = 18 573.7(1) eV (Sec. IX D). The recoil of the 3 HeT + molecule is given by E rec = E 2 0 + 2E 0 m e m HeT + = 1.720 eV.
E-gun data were used to investigate the work function of the inner electrode system of the main spectrometer. First, the work function of this electron source was measured with the Fowler method [115] to be Φ egun = 4.44(5) eV. Next, a transmission function was measured with photoelectrons from the e-gun traveling through an evacuated source. Knowing the energy of the transmission edge and the work function of the e-gun, we can estimate the work function of the inner-electrode system as Φ IE = 4.1(2) eV.
The β-electron starting potential inside the tritium source is defined by the cold and strongly magnetized plasma within its boundary conditions at the rear wall and the grounded beam tube (Sec. III C). By assuming that the magnitude of the plasma potential is small, as indicated by the 83m Kr measurement campaign, we treat the electron starting potential as mainly defined by the bias voltage and work function of the gold-plated rear wall, especially at small radii.
The work function of the rear wall was measured with the Fowler method prior to KNM1. Due to the illumination conditions, only the inner two-thirds of its area could be used for the measurement. The resulting raw, mean value from this measurement is Φ vac RW = 4.29 eV. However, this measurement was performed with an evacuated source. Previous measurements with deuterium gas indicate that the work function changes by about −100 meV when the rear wall is exposed to hydrogen isotopes in the source, as is the case during tritium operation. This estimate of the in situ work function of the rear wall has a large uncertainty, which we estimate at about ±200 meV. Further, during KNM1 the rear wall was set to a voltage of U RW = −150 mV, which is numerically equivalent to an increase of the work function by 150 meV. These considerations lead us to estimate an actual rear-wall work function of Φ RW = 4.34 (20) eV during KNM1.
We assume an additional uncertainty of ±100 mV for the sum of all involved voltages. The main contribution to this is the uncertainty of the absolute voltage of the main spectrometer, ∆U abs = ±94 mV [56]. The dominant uncertainty for the Q-value determination is the possibility of a plasma potential in the source that differs from the rear-wall potential. We assume an uncertainty of U plasma = ±400 mV because we cannot directly probe the plasma potential under KNM1 operational con-  25. Comparison of the Q-value of molecular tritium found in this work to values derived from Penning-trap measurements of the atomic-mass difference. In chronological order, the values of the Penning-trap measurements are those reported in Refs. [116], [117], [118] and [112].
ditions. Our final result is then: Q(T 2 ) KNM1 = E 0 +E rec −0.2 eV±0.5 eV = 18 575.2(5) eV. (25) Q(T 2 ) KNM1 and Q(T 2 ) ∆M agree within uncertainties. Figure 25 shows a comparison of the obtained Q-value in KATRIN with values derived from Penning-trap measurements. The consistency of the Q-values underlines the robustness of the energy scale in our scanning measurement of the T 2 β spectrum.

XIII. RESULTS AND DISCUSSION
In this work we have presented the first neutrino-mass measurement campaign of the KATRIN experiment. The acquired high-precision T 2 β decay spectrum, containing a total of 2 million electrons in an energy range of [E 0 −37 eV, E 0 +49 eV], was compared against a model of the theoretical spectrum, incorporating relevant experimental effects such as electromagnetic fields, backgrounds, and scattering. The experiment was operated at a reduced column density. Taking into account both the reduced activity and the reduced scattering probabilities, the β electrons recorded in the ROI during our four-week KNM1 campaign correspond to just 9 days of measurement time at the full, design source strength.
The full analysis was carried out applying a multi-stage blinding scheme. All analysis inputs were fixed on MC twin copies of the data; the spectrum model was blinded with a modified molecular final-state distribution; and finally the full analysis was performed using two independent analysis techniques (covariance matrix and MC propagation) which revealed a high degree of consistency.
We find excellent agreement of the calculated spectrum with the data. The covariance-matrix fit method obtains a goodness of fit of of χ 2 = 21. The effective spectral endpoint E 0 , which is inferred from the spectral fit alongside m 2 ν , can be related to the nuclear Q-value using the molecular recoil and the offset between the source potential and spectrometer work function. Our analysis gives a Q-value of 18 575.2(5) eV, which is in excellent agreement with measurements based on the 3 He-3 H atomic mass difference [112]. While the neutrino-mass result does not depend on the absolute energy scale of the spectrum, this consistency check is still of major importance to our understanding of the obtained spectra.
The best fit of the squared neutrino mass was found at m 2 ν = −1.0 +0.9 −1.1 eV 2 . The uncertainty is largely dominated by the statistical error of σ stat (m 2 ν ) = 0.97 eV 2 . If one were to assume the true neutrino mass to be equal to zero, the probability of obtaining this fit result given our total error budget is 16 %. The best-fit results of the covariance-matrix and MC-propagation techniques agree within 2 %.
We have applied three methodologies to derive an upper limit on the neutrino mass, based on the best-fit result. The Lokhov-Tkachov limit construction was developed in particular for direct neutrino-mass experiments [23]. By construction, in the case of a negative best-fit value of m 2 ν it yields the experimental sensitivity as an upper limit. Based on this technique we find m(ν e ) < 1.1 eV (90 % CL). The standard Feldman-Cousins technique for confidence-belt construction [22] yields an upper limit of 0.8 eV (90 % CL). Finally, we also apply Bayesian inference methods to the neutrinomass search, excluding negative values of m 2 ν through a flat, positive prior. The Bayesian result is presented in this work for the first time, yielding a 90 % credibility interval of 0 to 0.9 eV.
The newly obtained upper limit on the neutrino mass improves the previous best direct bounds by a factor of nearly two (Fig. 26, top). The effective 9 days of measurement time of this first neutrino-mass campaign (out of a total planned measurement time of 1000 days) led to an improvement of the statistical uncertainty on m 2 ν by a factor of two compared to the final results of the Troitsk and Mainz experiments [11,12] (Fig. 26, bottom), while the systematic uncertainties are reduced by a factor of six (Fig. 26, center).
The systematic error budget is expected to improve with future measurement campaigns. Most notably, new means to further suppress the background rate are now in place. These will increase the signal-to-background ratio and at the same time reduce the dominant systematic uncertainties related to the dependence of the background on time and retarding potential. Furthermore, in this first measurement the activity stability suffered from a burn-in phase, in which the structural material was exposed to tritiated gas for the first time. Subsequent to this first campaign, significant improvements of the activity stability have been demonstrated at an increased intensity about four times the KNM1 source strength. Finally, sub-dominant systematic effects, such as uncer- The total uncertainty is reduced by a factor of three. The historical measurements plotted here are: Los Alamos (1991) [18], Tokyo (1991) [119], Zurich (1992) [120], Mainz (1993) [121], Beijing (1993) [122], Livermore (1995) [103], Troitsk (1995) [123], Mainz (1999) [124], Troitsk (1999) [125], Mainz (2005) [11], Troitsk (2011) [12].
tainties in the final-state distribution, have been conservatively estimated for this analysis. Our knowledge of these systematics is expected to improve significantly in our future commissioning and measurement phases.

XIV. CONCLUSION
The new upper limit m ν < 1.1 eV (90 % CL) from KATRIN's first science run improves upon previous work [11,12] by almost a factor of two, based on a measuring period of only four weeks while operating at reduced column density -equivalent to just 9 days at nominal source strength.
In the coming years, KATRIN will soon reach the first sub-eV sensitivity, and finally tackle its ultimate design sensitivity of 0.2 eV (90 % CL). In addition, the precise measurement of the tritium spectrum allows searches for physics beyond the Standard Model, including righthanded weak currents [126] and sterile-neutrino admixtures with masses from the eV [127,128] to the keV scale [76]. KATRIN's model-independent probe of the neutrino mass is of paramount importance for both particle physics and cosmology. In particle physics, this measurement narrows the allowed range of quasi-degenerate neutrino-mass models.
In cosmology, it provides laboratory-based input for studies of structure evolution in ΛCDM and other cosmological models. In the absence of a definitive observation of dark matter, the neutrinomass scale is unique as a ΛCDM parameter that is directly observable in the laboratory.
Upcoming cosmological probes are expected to achieve a determination of the sum of the neutrino masses over the next 5 to 10 years, making this laboratory measurement particularly important for obtaining a consistent picture of the neutrino as both particle and dark-matter constituent in the universe. This first KATRIN result serves as a milestone towards this goal.