Deconvolution of the energy loss function of the KATRIN experiment
Volker Hannen Irina Heese Anna Sejersen Riis Kathrin Valerius Christian Weinheimer






 hannen@uni-muenster.de


Institut für Kernphysik, Westfälische Wilhelms-Universität Münster, Wilhelm-Klemm-Str. 9, 48149 Münster, Germany


Department of Physics and Astronomy, University of Aarhus, Ny Munkegade, DK-8000 Aarhus C, Denmark


Institut für Kernphysik, Karlsruher Institut für Technologie, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany


The KATRIN experiment aims at a direct and model independent determination of the neutrino mass with 0.2 eV/c^2 sensitivity (at 90% C.L.) via a measurement of the endpoint region of the tritium beta-decay spectrum. The main components of the experiment are a windowless gaseous tritium source (WGTS), differential and cryogenic pumping sections and a tandem of a pre- and a main-spectrometer, applying the concept of magnetic adiabatic collimation with an electrostatic retardation potential to analyze the energy of beta decay electrons and to guide electrons passing the filter onto a segmented silicon PIN detector.

One of the important systematic uncertainties of such an experiment are due to energy losses of -decay electrons by elastic and inelastic scattering off tritium molecules within the source volume which alter the shape of the measured spectrum. To correct for these effects an independent measurement of the corresponding energy loss function is required. In this work we describe a deconvolution method to extract the energy loss function from measurements of the response function of the experiment at different column densities of the WGTS using a monoenergetic electron source.



Neutrino mass, electron scattering, deconvolution









§ INTRODUCTION


The KArlsruhe TRItium Neutrino (KATRIN) experiment aims at determining the neutrino mass in a model independent way from the kinematics of tritium -decay. The observable in this case is an ”average electron anti-neutrino mass” given by the incoherent sum of neutrino mass eigenstates weighted by the squared elements of the mixing matrix. 
The experiment combines a Windowless Gaseous Tritium Source (WGTS) and a high resolution electrostatic retarding spectrometer (MAC-E filter) to measure the spectral shape of -decay electrons close to the endpoint energy at 18.6 keV with an unprecedented precision. KATRIN's sensitivity to the neutrino mass will be 0.2 eV/c^2 (at 90% C.L.) after 3 years worth of data taking <cit.>. 
An observed mass signal of 0.35 eV/c^2 will have a 5σ significance at the expected level of statistic and systematic uncertainties. 
In order to reach the desired sensitivity, all systematic effects of the measurement must be well under control with the major systematic uncertainties being allowed to contribute no more than Δ m^2 = 0.0075  eV^2/ c^4 to the systematic error budget.
An overview of the KATRIN experiment is shown in figure <ref>. The experiment starts with the WGTS where 



molecular  gas is injected at the center of the source and removed at both ends by turbo-molecular pumps. The  gas is kept at a constant temperature of 30 K within the source that is operated at a column density of 5· 10^17 cm^-2. The operational parameters of the source cryostat are monitored by a complex sensor network and a dedicated calibration and monitoring section at the rear of the source system <cit.>. 
About 10^10 -decay electrons are emitted per second into the accepted forward solid angle with pitch angles less than θ_max = 51^∘ and are guided magnetically through the transport section to the spectrometer tandem consisting of pre- and main-spectrometers. The task of the transport section made up of a differential pumping section and a cryo-pumping section is to suppress the flow of  molecules into the direction of the spectrometers by at least a factor of 10^14 in order to reduce experimental background from tritium decays within the spectrometers. A first energy discrimination is performed by the pre-spectrometer which rejects the low energy part of the  spectrum (up to 300 eV below the endpoint) and thereby reduces the rate of electrons going into the main spectrometer to approximately 10^3  s^-1. Like the pre-spectrometer the main spectrometer operates as a so-called MAC-E filter <cit.> and has the task to perform a precise energy analysis of the decay electrons.
In a MAC-E filter electrons are guided magnetically against an electrostatic retardation potential that can only be surpassed by electrons with  sufficiently high longitudinal energy with respect to the electric field. Here the longitudinal energy is given by E_∥ = E_ kin·cos^2θ with E_ kin being the kinetic energy of the electron and θ being the angle between electron momentum and  magnetic field direction. The transverse energy is accordingly given by E_⊥ =  E_ kin·sin^2θ.

The spectrometer acts as a high pass filter with a transmission function describing the observed electron rate as a function of the electron surplus energy (see section <ref>). To reduce the amount of transversal energy of the electrons that is not analyzed by the spectrometer, the technique of magnetic adiabatic collimation is used. The idea is that the magnetic guiding field drops by several orders of magnitude from the entrance of the spectrometer to the analyzing plane, where the electric potential reaches its maximum. If the gradient of the magnetic field is small enough, such that the field is approximately constant along one cyclotron loop of the electron movement, the magnetic moment of the cyclotron motion μ = E_⊥ / B (non-relativistic) is constant, and as B drops the transversal energy of the electrons is converted into longitudinal energy E_∥ that can be analyzed by the spectrometer. 
By varying the electric potential of the spectrometer it is then possible to scan the relevant region around the endpoint energy of tritium -decay and accumulate a spectrum. 

Electrons with sufficient energy to pass the spectrometer are finally detected by a 148 pixel silicon PIN detector <cit.> at the end of the setup.
Among the main systematical uncertainties of the experiment are energy losses from inelastic scattering of electrons in the source, fluctuations of the source density, fluctuations of the spectrometer analyzing potential, uncertainties in the transmission function and uncertainties in the final state distribution of the daughter molecules left after the decay reaction. A sophisticated calibration and monitoring system is being set up to keep the aforementioned systematic effects under control.

While there is some information on the energy loss of 18.6 keV electrons in gaseous tritium or quench condensed deuterium from the former neutrino mass experiments in Troitsk and Mainz <cit.>, precise experimental information on energy losses of electrons with energies near the endpoint of the tritium  spectrum are only available for molecular hydrogen as target gas <cit.>.
A measurement of the energy differential scattering cross section of 18.6 keV electrons off molecular tritium is therefore highly desirable. Such a measurement can be performed using a monoenergetic source of electrons mounted upstream of the WGTS to determine the response function of the overall experiment at different column densities of the source. A deconvolution method suitable to extract the energy loss function from the data will be presented in the following sections.
Once the energy loss function for tritium is known with sufficient accuracy, the same measurement setup can be used for an independent check of the column density of the WGTS during intervals between the regular measurement cycles of the KATRIN experiment <cit.>. 




§ ENERGY LOSS FUNCTION



The processes contributing to the energy loss of electrons traversing the molecular tritium gas within the WGTS are excitation of rotational and vibrational states of the T_2 molecules, excitation of electronic molecular states, dissociation and ionization of the molecules.

Aseev et al. <cit.> report on measurements of energy losses of electrons in gaseous tritium and in quench condensed deuterium films. Because of the limited energy resolution of a few eV the shape of the energy loss spectrum was not directly extracted from the data in their analysis, but approximated by a Gaussian representing electronic excitations and dissociation and a one-sided Lorentzian curve representing the continuum caused 
by ionization of the molecules. The parameters of the two functions were then adapted to fit the observed integral energy spectra obtained with an 18.6 keV mono-energetic electron source for gaseous tritium or from 17.8 keV mono-energetic conversion electrons from a ^83 mKr film covered by various thicknesses of D_2 absorbers. In both cases, energy losses caused by rotational and vibrational excitations of the molecules without electronic excitation could not be resolved and were neglected.

More detailed information is available for the scattering of 25 keV electrons from molecular hydrogen gas <cit.> where direct measurements of the energy loss function with resolutions down to 40 meV have been performed. 









This information about the scattering of electrons from molecular hydrogen has been implemented into a computer code by F. Glück <cit.> that can be used in simulations to generate energy losses Δ E and scattering angles Δφ in individual scattering events. The spectral shape produced with this routine is shown in figure <ref>. 



It is used in a toy Monte Carlo simulation of the WGTS to evaluate the deconvolution methods described in the following sections.

The probability for an electron of kinetic energy E to lose a specific amount of energy Δ E in a single scattering event is described by the differential energy loss function dσ/ dΔ E. For our purpose, we normalize the function by the total inelastic scattering cross section σ_ tot, obtaining[The integral over the energy losses runs up to E/2 since the incoming electron and the secondary electron in an ionisation process (assuming E is larger than twice the ionisation energy) are identical quantum particles.]

    f(Δ E) = 1/σ_ tot· dσ/ dΔ E     with    ∫^E/2_0 f(Δ E)   dΔ E = 1   .

The total inelastic scattering cross section for 18.6 keV electrons off gaseous tritium is given by σ_ tot( T_2) = (3.40 ± 0.07)· 10^-18  cm^2 <cit.>. 
The above mentioned code by F. Glück <cit.> for scattering of 18.6 keV electrons off hydrogen gives a total inelastic cross section of σ_ tot( H_2) = 3.7 · 10^-18  cm^2.




§ DECONVOLUTION METHOD


In the following sections we describe suitable mathematical methods to extract the energy loss function of 18.6 keV electrons in gaseous tritium from a series of measurements of the overall response function of the experiment at different column densities of the WGTS.



 §.§ Response function







Figure <ref> displays simulated response functions of the KATRIN experiment at different column densities ρ d = 0 and 
ρ_1 d < ρ_2 d < ρ_3 d assuming a mono-energetic electron source with narrow angular emission characteristics, i.e. pitch angles w.r.t. the magnetic field lines θ≤ O(1^∘). 
Shown is the transmission probability as a function of the nominal surplus energy E_s = E - qU of the electrons, i.e., the difference between the setpoint energy E of the electron source and the retardation potential qU of the main spectrometer. 
The response function at non-zero column density of the tritium source is given by a summation over contributions corresponding to n-fold (i.e. no, single, double, etc.) scattering of the electrons within the tritium source, weighted by the probabilities P_n for n-fold scattering


    R(E_s) = P_0 · T_e(E_s) + P_1 · T_e(E_s) ⊗ f(Δ E) + 
                P_2 · T_e(E_s) ⊗  f(Δ E) ⊗ f(Δ E) + …  .


Here T_e(E_s) is the experimental transmission function of the main spectrometer for the given electron source and f(Δ E) the sought after energy loss function neglecting small scattering angles. 
Given that we only consider a small energy interval of up to 50 eV below the endpoint energy of 18.6 keV of the  spectrum, the dependence of the energy loss function f(Δ E) on the kinetic energy of the electrons can be neglected.
The scattering probabilities are normalized such that ∑_n=0^∞ P_n = 1.

The experimental transmission function T_e(E_s) is determined from a measurement of the response function R(E_s) with an empty tritium source and hence without scattering. It is then equal to the analytical transmission function of the spectrometer T(E_s) convolved with a function S_e describing the energy spread and angular distribution of the electron source


    T_e(E_s) = . R(E_s)|_ρ d=0 = T(E_s) ⊗ S_e   .


In the simulations, the energy spread of the source is described by a Gaussian smearing of the energy setpoint with a width of σ_e = 0.2 eV. It is assumed that the electron source has a small angular divergence with starting angles θ_e ≤ 0.5^∘ := θ_e,max relative to the magnetic field direction at its location at the rear end of the WGTS. Within this narrow cone, the emission angles of the electrons are assumed to be isotropically distributed. 
Suitable electron sources that emit single electrons at adjustable total energy and adjustable emission angle have been developed and tested within the KATRIN collaboration <cit.>.
The above mentioned numerical values for the energy and angular spread are compatible with the characteristics of the photo-electron source used during the commissioning of the KATRIN main spectrometer <cit.>.

The analytical transmission function of the spectrometer T(E_s) for such a source is given by the following relation <cit.>


    T(E_s) = {[                                                       0                                            for   E < qU; 1 - √(1 - E-qU/EB_e/B_A)/1 - √(1 - E_⊥,A, max/EB_e/B_A)                          for   qU ≤ E ≤ qU + E_⊥,A, max;                                                       1                              for   qU + E_⊥,A, max < E, ].


where B_e is the magnetic field at the electron source, B_A the magnetic field at the analyzing plane of the spectrometer and E_⊥,A, max =  Esin^2(θ_e, max)B_A/B_e is the maximum remaining transversal energy component in the analyzing plane.


Defining

    ϵ_0(E_s)     =     T_e(E_s) 
    ϵ_1(E_s)     =     T_e(E_s) ⊗ f(Δ E) 
    ϵ_2(E_s)     =     T_e(E_s) ⊗ f(Δ E) ⊗ f(Δ E) 
       …

as the n-fold scattering functions we can rewrite equation <ref> to obtain


    R(E_s)= P_0 ·ϵ_0(E_s) +  P_1 ·ϵ_1(E_s) +  P_2 ·ϵ_2(E_s) + …  .


If we manage to determine the single scattering function ϵ_1(E_s) from measured response functions we can, with the knowledge of T_e(E_s), extract the energy loss function f(Δ E) using suitable deconvolution methods. 




 §.§ Scattering probabilities


The mean free path of the electrons within the tritium gas inside the WGTS can be expressed in terms of a mean free column density (ρ d)_ free which the electrons pass before an interaction and which is calculated taking the inverse of the total scattering cross section (ρ d)_ free = 1/σ_ tot.
The actual column density seen by an electron traversing the WGTS at an angle θ relative to the symmetry axis, i.e. the magnetic field axis, is given by ρ d / cosθ. Neglecting possible scattering angles Δφ in the scattering processes for the moment, the mean number of expected scatterings is

    μ(θ) = ρ d/(ρ d)_ freecosθ =  ρ d  σ_ tot/cosθ 
                 = μ_0/cosθ  .

The probability for an n-fold scattering is given by a Poissonian distribution:

    P_n(μ(θ)) = μ^n(θ)/n!exp(-μ(θ))    with   n = 0,  1,  2, …  .

We have to take into account that electrons generated by the electron source follow an angular distribution which is, for our purpose, assumed to be isotropic within a narrow interval between 0^∘≤θ_e ≤θ_e, max. 
If we further take into account that the magnetic field at the location of the electron source will be lower than the field within the WGTS, the starting angles have to be transformed according to 

    θ = arcsin( sinθ_e ·√(B_ WGTS/B_e)) ≈θ_e ·√(B_ WGTS/B_e)  ,

resulting in angles θ of electron momenta relative to the magnetic field direction within the WGTS.

To obtain average scattering probabilities we weigh the values from equation <ref> with g(θ) = sinθ, corresponding to the probability function of an isotropic distribution, integrate over the given range of angles and normalize 

    P_n(μ_0)     =    ∫^θ_ max_0 g(θ) P_n(μ(θ))   dθ / ∫^θ_ max_0 g(θ)   dθ
        =    ∫^θ_ max_0 sinθ·(μ_0/cosθ)^n/n!exp(-μ_0/cosθ)   dθ /  
           (1 - cosθ_ max)   .

Equation <ref> delivers approximate values for the scattering probabilities, as it does not take into account changes in the direction of the electrons during scatterings. Secondly this equation also assumes a homogeneous distribution of tritium molecules in transverse direction within the WGTS. 
Compared to scattering probabilities extracted from the simulations accounting for scattering angles in elastic and inelastic scattering described in section <ref> the deviations to the results calculated with <ref> were found to be on the <10^-3 level, however. The small difference is due to the fact, that the scattering angles for inelastic scattering in our energy loss range of interest (Δ E < 50 eV) are strongly forward peaked with a mean of Δφ=0.5^∘ and elastic scattering is a subdominant process (see figure <ref>, right). 


To obtain more precise values for the scattering probabilities, a detailed simulation using a 3-dimensional description of the column density within the WGTS is required, which is beyond the scope of this paper.




 §.§ Extraction of the single scattering function



In order to determine the single scattering function ϵ_1(E_s) we have to perform measurements of the response function R(E_s) at different column densities. Neglecting multiple scattering events with more than three interactions of the electrons with the tritium gas inside the WGTS[The probability of a 4-fold scattering process with the maximum energy loss under consideration Δ E < 50 eV is below 8% at the maximum column density of 5· 10^17  cm^-2 used in the simulations.], we can set up a system of linear equations for measurements at three column densities labeled a, b and c:


    R^a(E_s) - P_0^a · T_e(E_s)     =     P_1^a ·ϵ_1(E_s) + P_2^a ·ϵ_2(E_s) + P_3^a ·ϵ_3(E_s)   ,
    
      R^b(E_s) - P_0^b · T_e(E_s)     =     P_1^b ·ϵ_1(E_s) + P_2^b ·ϵ_2(E_s) + P_3^b ·ϵ_3(E_s)    ,
    
      R^c(E_s) - P_0^c · T_e(E_s)     =     P_1^c ·ϵ_1(E_s) + P_2^c ·ϵ_2(E_s) + P_3^c ·ϵ_3(E_s)   ,


which we can write as a matrix equation:


    R⃗ - P⃗_0 · T_e(E_s) = 𝐏·ϵ⃗   with  𝐏 = ([ P_1^a P_2^a P_3^a; P_1^b P_2^b P_3^b; P_1^c P_2^c P_3^c ])   .


Taking higher scattering orders into account would require additional measurements at further non-zero column densities and would increase the dimension of the system of linear equations to be solved. Whether the inclusion of only three scattering orders provides sufficiently accurate results will be evaluated in section <ref>.
Multiplying with the inverse of 𝐏, which is calculated using the Gauss-Jordan algorithm from the ROOT software package <cit.>, we obtain 

    ϵ⃗ = 𝐏^-1·(R⃗ - P⃗_0 · T_e(E_s))

from which we can calculate the single scattering function ϵ_1(E_s).




 §.§ Deconvolution of the energy loss function



As described in section <ref>, the single scattering function is the result of the convolution of the experimental transmission function of the spectrometer with the energy loss function. This convolution is calculated taking the integral

    ϵ_1 (E_s) = T_e(E_s) ⊗ f(Δ E) = ∫^E/2_0 T_e(E_s-Δ E)f(Δ E)   dΔ E   .

In our case, where the values of the functions in question are only known at N equally distributed discrete measurement points defined by the applied retardation voltage U_i, the integral is replaced by a sum

    ϵ_1 (E-qU_i) = ∑^N-1_j=0 T_e(E-qU_i-Δ E_j)f(Δ E_j)   .

The latter equation can be rewritten in N × N matrix form 

    ϵ⃗_⃗1⃗ = 𝐓_𝐞·f⃗  ,

where the 𝐓_𝐞 matrix is constructed from the discrete transmission function T_e(E_s,i = E-qU_i) as[The zeroes in the right upper corner of 𝐓_𝐞 are caused by the transmission condition E-qU_i-Δ E_j ≥ 0.]

    𝐓_𝐞 = ( [   T_e(E_s,0)            0                                      ⋯            0;   T_e(E_s,1)   T_e(E_s,0)            0                         ⋯            0;   T_e(E_s,2)   T_e(E_s,1)   T_e(E_s,0)            0            ⋯            0;            ⋮            ⋮            ⋮            ⋮                         ⋮; T_e(E_s,N-1) T_e(E_s,N-2)            ⋯                             T_e(E_s,0) ])   .


One could now try to solve equation <ref> by multiplying with the inverse of the 𝐓_𝐞 matrix. The latter, however, is close to being singular and cannot easily be inverted numerically. We therefore have to apply more sophisticated methods to deconvolve the energy loss function from the matrix equation. In the following two methods are applied to the problem: the so-called Singular Value Decomposition (SVD) <cit.> and the iterative Stabilized Biconjugate Gradient method <cit.>.




  §.§.§ Singular Value Decomposition



The Singular Value Decomposition (SVD) is a method to deal with systems of linear equations given by a matrix equation 𝐀·x⃗ = b⃗ that are either singular or numerically very close to singular and is able to provide useful, although not necessarily unambiguous, solutions to the given problem. 
It is based on the theorem that any M× N matrix 𝐀 with M≥ N can be written as a product of an M× N column-orthogonal matrix 𝐔, an N× N diagonal matrix 𝐖 whose elements are the so-called singular values w_i ≥ 0, and the transpose of an N× N orthogonal matrix 𝐕 <cit.>:

    𝐀 = 𝐔·𝐖·𝐕^𝐓 = 
     𝐔·( [ w_1   0;     w_2;       …;   0 w_N ]) ·𝐕^𝐓  .

As 𝐔 and 𝐕 are orthogonal the inverse of equation <ref> can be written as

    𝐀^-1 = 𝐕·𝐖^-1·𝐔^𝐓 = 𝐕· [diag(1/w_i)] ·𝐔^𝐓  .

Problems arise when some of the singular values w_i are either zero or so small that their values are dominated by numerical rounding errors. Using the SVD method it is still possible to construct an approximate solution vector x⃗ that will minimize the residual r given by 

    r ≡ |𝐀·x⃗ - b⃗|  .

For that purpose all diagonal elements of 𝐖^-1 where the singular values w_i are below a chosen threshold value w_ thr are set to zero, thereby removing infinite or problematically large matrix elements. The matrix constructed using the modified 
diagonal matrix 𝐖̃^-1 is the so-called pseudoinverse matrix Ã^-1, and the solution vector x⃗ is then given by

    x⃗≈Ã^-1·b⃗ = 𝐕·𝐖̃^-1·𝐔^𝐓·b⃗  ,

which, translated to our original problem of deconvoluting the energy loss function from the measured single scattering function (see equation <ref>), becomes

    f⃗≈𝐓̃_𝐞^-1·ϵ⃗_⃗1⃗ = 𝐕·𝐖̃^-1·𝐔^𝐓·ϵ⃗_⃗1⃗  .

What remains to be settled is the optimal threshold value w_ thres for suppression of the problematic singular values. 
This can only be determined by investigating the influence of the deconvolved energy loss function on the extracted neutrino mass values in simulated neutrino mass runs of the KATRIN experiment. Such a study, applying a toy Monte Carlo simulation of the experiment, is presented in section <ref>.




  §.§.§ Stabilized Biconjugate Gradient method



As an alternative to the SVD method we tested the so called Stabilized Biconjugate Gradient (Bi-CGSTAB) method described by Sleijpen and Fokkema <cit.>. 
The bi-conjugate gradient method iteratively solves linear sets of equations 𝐀·x⃗ = b⃗ where 𝐀 is an N× N matrix. In each iteration step, labeled k, the approximate solution x⃗_k is modified by some search correction that depends on the true residual r⃗_k = b⃗ - 𝐀·x⃗_k and some “shadow residual” r⃗̃⃗_⃗k⃗
calculated using the transpose 𝐀^T.
The residuals are forced to converge by making r⃗_k orthogonal to the shadow residuals r⃗̃⃗_⃗j⃗ for j < k.


For an in-depth description of the algorithm we refer to reference <cit.>. In our simulations we use the implementation of the Bi-CGSTAB algorithm provided as part of the Meep software package for finite-difference time-domain simulations developed at the Massachusetts Institute of Technology (MIT) <cit.>. As for the SVD method, the energy loss function resulting from a deconvolution using the Bi-CGSTAB algorithm is evaluated in a toy Monte Carlo simulation of the experiment, as presented in section <ref>.




§ EVALUATION USING TOY MONTE CARLO SIMULATION



In order to test the deconvolution methods described in section <ref>, we performed Monte Carlo simulations of response function measurements at different column densities. In these simulations we assumed a perfectly homogeneous gas distribution within the WGTS and used the model of F. Glück <cit.> to generate energy losses in scattering events of 18.6 keV electrons with molecular hydrogen. Table <ref> provides an overview of the input parameters of the simulation.



The simulated response functions are displayed in figure <ref> (left). From these response functions we can 



extract the n-fold scattering functions as described in section <ref> with the results shown in figure <ref> (right). The single scattering function ϵ_1(E_s) is then the input for either one of the two deconvolution methods described in section <ref>.




 §.§ SVD results


To optimize the result obtained for the energy loss function by the SVD method, we have to determine the optimum value of the threshold w_ thr below which the corresponding matrix elements of the inverse diagonal matrix 𝐖^-1 are discarded.
For that purpose we scan over a range of possible values for w_ thr and with each value calculate the deconvolved energy loss function f⃗_ SVD and the difference |f⃗_ mod - f⃗_ SVD| to the input energy loss model f⃗_ mod. 
Figure <ref> displays on the left the result of this scan. 



The threshold values are given in percent of the maximum singular value present in the diagonal matrix 𝐖. 
Deconvolved energy loss functions obtained at threshold values of 0.2%, 0.3% and 0.6% are shown together with the input model. The suppression of low singular values in the SVD acts by damping numerical fluctuations in the deconvolved energy loss function. While the result obtained with w_ thr = 0.2% exhibits a significant amount of high frequency noise, a value of w_ thr = 0.6% produces a smoother result at the expense of washing out structures observed in the input model.
Besides this smoothening of the deconvolved function we also note an oscillatory response to large spikes in the input model, that is best visible in the lower right plot of figure <ref>.

An optimum resemblence of the deconvolved function to the input model is found at a threshold value of w_ thr = 0.3%.

It should be noted that the fact, that the SVD method cannot yield the fine structure of the energy loss function (see figure <ref>) in the interval 11 eV≤Δ E ≤ 16 eV at that threshold value, does not matter, since KATRIN`s transmission function for an isotropic electron source as the tritium source will have a width of 0.93 eV. We will investigate the influence on the neutrino mass measurement in section 4.3.



 §.§ Bi-CGSTAB results


Figure <ref> displays the result obtained for the energy loss function using the iterative Bi-CGSTAB algorithm described in section <ref>.



While the method does seem to reproduce some of the fine-grained structure related to inelastic electron excitations starting at Δ E ≥ 12 eV, it is obviously much noisier than the results obtained with the SVD method. To counteract the noise, we tried to filter the deconvolved function with a second-order low pass Butterworth filter with a cut off frequency of 1 eV^-1, with the result shown in green in figure <ref>. The filtered function does exhibit lower noise and still follows the general features of the input model. 




 §.§ Influence on the measured neutrino mass



To assess the influence of the deconvolution process on the determination of the neutrino mass from the experimental beta spectrum measured by the KATRIN experiment, a large number of such measurements have been simulated and analyzed using the different energy loss functions displayed in figures <ref> and <ref>. From the mean of the distribution of the extracted ^2 values we can then estimate the systematic uncertainty caused by the deconvolution process.


The energy spectrum of nuclear  decay can be calculated starting from Fermi's golden rule and has the following form <cit.> (with units ħ = c = 1):


    dΓ/ dE    =    G_ F^2 cos^2θ_ C/2π^3  |M|^2   F(E,Z+1)   p   (E + m_e)  
          ·∑_j 𝒫_j   E_,j √(E_,j^2-^2) Θ(E_,j - )   ,


where G_ F is the Fermi coupling constant, θ_ C the Cabibbo angle and M the nuclear matrix element of the transition. The Fermi function F(E, Z+1) takes into account the final state interaction of the emitted electron with the daughter nucleus of charge Z+1 and p(E+m_e) is the phase space factor of the outgoing electron. The product of the neutrino momentum and its energy given by E_ν,j = E_0-E^*_j- E is the phase space of the emitted anti-neutrino, which shapes the  spectrum near its endpoint. The neutrino phase space factor has to be summed up over all final states E^*_j of the daughter molecule that are populated with probabilities 𝒫_j. The inclusion of the Θ function in eq. <ref> ensures that E_,j -  > 0.
The observable ^2 that can be extracted from the spectral shape near the endpoint is defined by an incoherent sum over the neutrino mass eigenstates m_i weighted by the matrix elements of the U_PMNS mixing matrix <cit.> known from oscillation experiments:


    ^2 = ∑^3_i=1 | U_ei |^2 m_i^2   .


The simulation of the observed integral spectrum measured by KATRIN uses the electron spectrum described by equation <ref> and takes into account a number of experimental effects:


 
 
  * The final state distribution of the ^3HeT^+ daughter molecules calculated by Saenz et al. <cit.> as shown in figure <ref> (left).
 
  * Up to four-fold scattering of the electrons within the WGTS according to the same energy loss model as used in the simulations of the deconvolution process.
 
  * The nominal transmission function of the main spectrometer.
 
  * An experimental background rate of 10^-2 counts per second in the energy region of interest at the FPD.
 
  * Statistical fluctuations according to an effective three years worth of data taken at the design value for the column density of the source of 5× 10^17/ cm^2 and with a measurement time distribution as described in the KATRIN Design Report <cit.>. 

For the purpose of this study, a vanishing neutrino mass = 0 is assumed.




With these inputs, 1000 hypothetical KATRIN measurements were simulated and subsequently fitted taking into account the same physical effects as in the simulations, but using one of the deconvolved functions to describe the energy losses and using as free parameters 


  * the squared neutrino mass ^2,

  * the spectral endpoint energy E_0,

  * the experimental background rate, and

  * the overall amplitude of the spectrum.

    
Figure <ref> (left) shows the resulting distribution of extracted ^2 values obtained with the 



energy loss function from the SVD deconvolution at a threshold value of w_ thr = 0.3 %. The mean of the distribution is at μ = (0.0053 ± 0.0005)  eV^2 yielding the systematic uncertainty of the deconvolution method, compatible with the KATRIN error budget for this systematic effect of 0.0075  eV^2 <cit.>. The systematic uncertainties corresponding to usage of the energy loss functions obtained with the different threshold values of the SVD algorithm and with the iterative Bi-CGSTAB algorithm (with and without filter), respectively, are compared in figure <ref> (right).
The energy loss functions deconvolved using the SVD algorithm result in significantly lower systematic errors than with the iterative Bi-CGSTAB algorithm and are, for all three threshold values tested, within the error budget defined for this uncertainty in the KATRIN design report <cit.>.
As a cross check, the fits were also run using the input model itself for the energy loss function and the resulting ^2 distribution was indeed centered around zero with μ = (0.0003 ± 0.0005)  eV^2.




§ SUMMARY AND OUTLOOK


A method to deconvolve the energy loss function of -decay electrons in the tritium source of the KATRIN experiment from measurements of the response function at different column densities with a mono-energetic electron source has been developed and tested in simplified Monte Carlo simulations of the experiment. Two different algorithms to deconvolve the energy loss function from the single scattering function extracted from the measurements have been tested: Firstly, the so-called Singular Value Decomposition (SVD), and secondly, the Stabilized Biconjugate Gradient (Bi-CGSTAB) method as an example of an iterative algorithm.
Both methods allow to obtain an approximation of the real energy loss function with varying levels of detail and of numerical fluctuations in the deconvolved result. The numerical noise can be dampened in the case of the SVD algorithm by setting an appropriate threshold value for the suppression of small singular values and in the case of the Bi-CGSTAB algorithm by filtering the deconvolved function. When applying the deconvolved energy loss function to the analysis of simulated integral  spectra, we find that the SVD method delivers usable results that induce a systematic error below the error budget of 0.0075  eV^2 allocated for this contribution in the KATRIN design report. The Bi-CGSTAB method, however, significantly exceeds the limit, both for the direct and for the filtered result of the deconvolution.

An option to improve the result from SVD deconvolution that has been looked at during simulations is to split the energy loss function into several energy intervals and perform a sectionwise deconvolution of the data. In this way it is possible to choose lower threshold values in the regions of rotational and vibrational states and of inelastic electron excitations where more structure is expected and to dampen numerical fluctuations by higher threshold values in between these regions and in the ionization tail, where we know that the curve is smooth.
With such an approach it was indeed possible to get a closer approximation of the fine structure observed in the input model of the energy loss function. However, the energy loss function deconvolved in this manner led to larger systematic errors in the fits of the neutrino mass spectra and was therefore discarded.

Further improvements of the deconvolved results are possible by increasing the statistics of the response function measurements. 
In the simulations we have assumed to collect 10 million events at each voltage step of the measurements. Assuming an event rate of 25 kHz this would result in about 10 days measurement time plus overhead, e.g., for switching between voltage settings during the measurements and for changing the column density of the WGTS in between runs. 

Moreover, one might consider including four-fold scattering in the extraction of the single scattering function (see section <ref>) to increase the precision of the subsequently deconvolved energy loss function. However, also this action comes at the cost of increased measurement time, as this would require to measure at yet an additional, higher column density setting of the source.

Finally it could be investigated if going beyond the present precision might be possible by applying Bayesian methods as used to reconstruct spectral functions in lattice QCD <cit.>.




§ ACKNOWLEDGMENTS
 
This work was supported by the German Federal Ministry of Education and Research (BMBF) under grant no. 05A11PM2 and 05A14PMA. 
KV wishes to acknowledge support by the Helmholtz Association (Young Investigators Group VH-NG-1055).




§ REFERENCES


99

KAT04 J. Angrik et al., (KATRIN Collaboration), 
   KATRIN Design Report 2004, 
   http://bibliothek.fzk.de/zb/berichte/FZKA7090.pdfWissenschaftliche Berichte, FZ Karlsruhe 7090 

WGTS M. Babutzka et al., 
  “Monitoring of the properties of the KATRIN Windowless Gaseous Tritium Source”,
  http://dx.doi.org/doi:10.1088/1367-2630/14/10/103046New J. Phys. 14 (2012) 103046

Pic92 A. Picard et al.,
  A solenoid retarding spectrometer with high resolution and transmission for keV electrons,
  http://dx.doi.org/10.1016/0168-583X(92)95119-CNucl. Inst. and Meth. B 63 (1992) 345-358

Ams15 J.F. Amsbaugh et al.,
  Focal-plane detector system for the KATRIN experiment,
  http://dx.doi.org/10.1016/j.nima.2014.12.116Nucl. Inst. and Meth. A 778 (2015) 40-60

Ase00 V.N. Aseev et al.,
  Energy loss of 18 keV electrons in gaseous T_2 and quench condensed D_2 films,
  http://dx.doi.org/10.1007/s100530050525Eur. Phys. J. D 10 (2000) 39-52

Gei64 J. Geiger,
  Streuung von 25 keV-Elektronen an Gasen,
  http://dx.doi.org/10.1007/BF01380873Zeitschrift für Physik 181 (1964) 413-425

Uls72 R.C. Ulsh et al.,
  Bethe surface, elastic and inelastic differential cross sections, Compton profile, and binding effects for H_2 obtained by electron scattering with 25 keV incident electrons,
  http://dx.doi.org/10.1063/1.1680755J. Chem. Phys. 60 (1974) 103

Glueck F. Glück, computer code scatter, private communication, 2008

Val11 K. Valerius et al.,
  Prototype of an angular-selective photoelectron calibration source for the KATRIN experiment,
  http://dx.doi.org/10.1088/1748-0221/6/01/P01002JINST 6 (2011) P01002

Bec14 M. Beck et al.,
  An angular-selective electron source for the KATRIN experiment,
  http://dx.doi.org/10.1088/1748-0221/9/11/P11020JINST 9 (2014) P11020

Beh16 J. Behrens et al.,
  A pulsed, mono-energetic and angular-selective UV photo-electron source for the commissioning of the KATRIN experiment,
  to be submitted to Eur. Phys. J. C

root R. Brun and F. Rademakers,
  ROOT - An object oriented data analysis framework,
  http://dx.doi.org/10.1016/S0168-9002(97)00048-XNucl. Inst. and Meth. A 389 (1997) 81-86

numrec W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
  Numerical Recipes. The Art of Scientific Computing,
  http://numerical.recipes/3rd ed., New York: Cambridge University Press (2007)

Sle93 G.L.G. Sleijpen and D.R. Fokkema, 
  BiCGstab() for linear equations involving unsymmetric matrices with complex spectrum,
  https://eudml.org/doc/118652Electronic Transactions on Numerical Analysis 1 (1993) 11-32

meep A.F. Oskooi et al.,
  Meep: A flexible free-software package for electromagnetic simulations by the FDTD method,
  http://dx.doi.org/10.1016/j.cpc.2009.11.008Computer Physics Communications 181 (2010) 687-702

Ott08 E.W. Otten and C. Weinheimer, 
  Neutrino mass limit from tritium  decay,
  http://dx.doi.org/10.1088/0034-4885/71/8/086201Rep. Prog. Phys. 71 (2008) 086201

PDG14 K.A. Olive et al. (Particle Data Group),
  Review of Particle Physics,
  http://dx.doi.org/10.1088/1674-1137/38/9/090001Chin. Phys. C 38 (2014) 090001 

Sae00 A. Saenz, S. Jonsell and P. Froelich, 
  Improved Molecular Final-State Distribution of HeT^+ for the -Decay Process of T_2,
  http://dx.doi.org/10.1103/PhysRevLett.84.242Phys. Rev. Lett. 84 (2000) 242 

burnier13 Y. Burnier and A. Rothkopf,
    Bayesian Approach to Spectral Function Reconstruction for Euclidian Quantum Field Theories,
    http://dx.doi.org/10.1103/PhysRevLett.111.182003Phys. Rev. Lett. 111 (2013) 182003 


