Monte Carlo simulations of the electron-gas interactions in the KATRIN experiment

At the KATRIN experiment, the electron antineutrino mass is inferred from the shape of the $\beta$-decay spectrum of tritium. Important systematic effects in the Windowless Gaseous Tritium Source (WGTS) of the experiment include the energy loss by electron scattering, and the extended starting potential. In the WGTS, primary high-energy electrons from $\beta$-decay produce an extended secondary spectrum of electrons through various atomic and molecular processes including ionization, recombination, cluster formation and scattering. In addition to providing data essential to the simulation of energy loss processes, the electron spectrum also provides information important in the simulation of plasma processes. These simulations will then provide an insight on the starting potential. Here, a Monte Carlo approach is used to model the electron spectrum in the source for a given magnetic and electric field configuration. The spectrum is evaluated at different positions within the WGTS, which allows for a direct analysis of the spectrum close to the rear wall and detector end of the experiment. Alongside electrons, also ions are tracked by the simulation, resulting in a full description of the currents in the source.


Introduction
The KArlsruhe TRItium Neutrino (KATRIN) experiment is designed to measure the electron antineutrino mass to unprecedented sensitivity down to 0.2 eV [1]. At the experiment, the mass is inferred from the shape of the β-decay spectrum of tritium. The contribution of the mass on the shape is most dominant in the endpoint region. Thus, the KATRIN experiment is designed to provide a high energy resolution, as well as a high event rate in this region of interest. The energy resolution is reached through a spectrometer with an energy resolution of 2.8 eV at 18.6 keV kinetic energy of the electron [2]. The measurement of the spectrum is performed down to 40 eV below the endpoint [2]. This part of the spectrum is well described by the convolution of the response function of the experiment and the theoretical beta spectrum, see [2] for more information. This paper however focuses on the evaluation of charged particle properties far below the endpoint, in the high luminous tritium source, called the Windowless Gaseous Tritium Source (WGTS). Gaseous tritium is circulating at an ambient temperature of 80 K, see section 2.1. A strong magnetic field of 2.5 T is directed along the WGTS. The source tube of the WGTS is closed off at one side by the so-called rear wall, a gold-plated disc. At the other side, the source tube of the WGTS is continued in the differential pumping section (DPS), which houses additional pumps for tritium removal and electrodes for ion flux reduction.
When tritium decays, the resulting primary electrons follow the magnetic field lines and are either terminated directly at the rear wall or reach the DPS section of the experiment. All electrons except for those with high energy are repelled here at the dipole electrodes, which are at negative potential ( DP,min = −85 V [3]). The high energy component of the electrons reaches the spectrometer. However, only a small portion of these electrons (relative factor of approximately 10 −10 [3,4]) reach the detector, and contribute to the neutrino mass measurement. The majority is reflected back to the source. So in total, most of the electrons will eventually be terminated at the rear wall.
While the description of the electron movement seems straightforward, the tritium gas is also a target for the primary electrons from tritium beta decay, producing secondary electrons and altering the primary electron spectrum [5]. This leads to a further complication; previous studies [6,7] have already shown that the density of electrons is sufficient to assume a plasma inside the WGTS. The plasma potential changes the relative energy of primary electrons to the detector. A precise knowledge of the electron spectrum and in turn also the plasma itself is necessary to understand the signal at the detector, see section 2.2.
Studies of the total electron spectrum [6] in the source exist for more than a decade, but have proven to be not sufficient to answer some questions: the electron spectrum at arbitrary positions (which is required for the simulation of kinetic plasmas) and the current of all charged particles (not only electrons) to the various walls of the experiment. Only some particle currents can be measured by the experiment, namely the total current towards the rear wall and the ion current towards the DPS [3]. Thus, the non-measurable currents must be determined through simulation.
The basic simulation concept of using a Monte Carlo approach has proven to be successful by Nastoyashchii et al. [6], but since then more precise interaction cross sections are known. Due to recombination processes and undetermined ion currents, which influence the development of the plasma, it is self-evident that the simulation by Nastoyashchii et al. of electrons alone does not provide a self-consistent picture of the processes in the WGTS.
In this paper, we will outline which physical processes are implemented in the newly developed code named KARL (and which cross sections are used to model these processes) and how this is realized in the numerical experiment. After a description of validation methods, the code is then used to provide electron spectra, particle densities and particle currents at multiple positions inside the source at current measurement conditions. The code base for KARL was initiated in the Masters's Thesis of C. Reiling [8], which provides the groundwork for the simulation of the relevant physical processes inside the source. Overview of the beam line. Tritium gas is injected in the center of the windowless gaseous tritium source (WGTS) cryostat and pumped off at both sides. The neutral gas flow is further reduced in the differential pumping section (DPS) and the cryogenic pumping section (CPS). The simulations presented below will cover the 16 m long WGTS cryostat, where the creation and interaction probability of electrons and ions is the highest. Figure adapted from [11].

The experiment
At the KATRIN experiment, the neutrino mass is inferred from the measurement and subsequent spectral shape analysis of tritium β-decay electrons [2]. The experiment has taken data since 2018. The data of the first two measurement campaigns combined yielded the world's best neutrino mass limit from β-spectroscopy with 0.8 eV (90% CL) [9]. The target sensitivity will be reached after five years of operation [1], reducing the statistical uncertainty by a factor of approximately ten to the previous result by Kraus et al. [10].

Overview and principle of work
The KATRIN experiment uses single β-spectroscopy of molecular tritium T 2 to determine the electron neutrino mass. The experimental setup is shown in figure 1. It consists of seven main components, which will be described in the following.
In the Windowless Gaseous Tritium Source (WGTS) β-electrons are created by β-decay of T 2 , with an endpoint energy of 0 ≈ 18.6 keV [12]. Molecular tritium gas is continuously injected in the center of the source and differentially pumped-off at both ends using turbo molecular pumps (TMP) [13]. The gas decays with an activity of 10 11 Bq [1]. Strong magnetic fields guide the electrons along the beam tube [14]. "At the upstream end, the WGTS terminates in a gold-plated rear wall" [2]. At the downstream end, the WGTS is connected to the differential pumping section [15] and the cryogenic pumping section [16]. These two sections further reduce the flux of remaining neutral tritium by a factor of at least 10 15 [3]. This value exceeds the design requirement of 10 14 , posted by consideration of residual tritium background in the spectrometers [1]. The high energy electrons, which leave the pumping sections, are then filtered at the spectrometers. The spectrometers use the MAC-E filter principle [3]. The transversal kinetic energy of the electrons is converted to longitudinal kinetic energy. The movement is opposed by an electric potential of the spectrometers. Only the small fraction of electrons with an energy greater than the electric potential can pass and are measured at the detector. The neutrino mass is then inferred from the shape of the measured spectrum.
In total, most of the electrons leaving the WGTS in the direction of the DPS are either reflected by the large potential of the spectrometers or are reflected by the relatively small potential of the Figure 2. "Cross-sectional view of the WGTS. The gas density profile inside the beam tube is shown in green" [17], with flow direction indicated by the arrows. The neutral gas flow, which will be obtained through simulations by Kuckert et al. [13], will be used as an input for our simulations presented below. Figure from [17]. dipole electrodes of the DPS, which are installed in the DPS to filter out the ions through the use of the × drift [3].
The WGTS consists of a cylindrical source tube, with a length of 16 m and a radius of 4.5 cm, see figure 2. Tritium is injected in the center, with a stable injection pressure of up to 10 −3 mbar [3]. The gas then streams to both sides of the inlet and is mostly pumped out at the four attached pumping ports through turbomolecular pumps at approximately 5 m and 6.5 m distance from the inlet. In total, the tritium density drops from the center of the WGTS towards its sides by a factor of approximately 1000 [13]. The pumped gas is purified and ultimately fed back into the source. Overall, the WGTS can produce a column density of 5 × 10 21 m −2 for standard KATRIN operation [3] with a stability of 1 part in 1000 [9].
The WGTS is cooled down to reduce thermal Doppler blurring of the β-spectrum "and to ensure a high tritium density at low flow rate. On the other hand, a low tritium pressure and small clustering of tritium molecules is required" [18]. The measurement campaigns KNM1 and KNM2 were therefore performed with a 30 K setting [4,9]. All following campaigns are performed at 80 K "to allow for equal source conditions in tritium and krypton-83m commissioning measurements" [18]. Thus, the results presented here were calculated for the 80 K setting.
The WGTS cryostat also houses seven superconducting magnets, which provide the constant guiding magnetic field of 2.5 T [17] for the electrons.

Plasma effects in the WGTS
The final neutrino mass sensitivity of 0.2 eV can only be reached through a precise examination of all systematic effects [2]. One of these effects is the plasma inside the source. It is generated through β-decay and subsequent impact ionization of neutral tritium gas (all relevant processes are listed in section 4). The resulting charged particle density is sufficient to assume a plasma inside the WGTS [7]. The effects of this plasma can influence the motion of the electrons.
The plasma, in conjunction with boundaries, exerts a potential throughout the source, also called starting potential, which can show a non-homogeneous structure. The final energy of the β-electrons selected by the spectrometers can therefore be dependent on the original position of β-decay inside the source. Hence, the plasma potential can distort the integral spectrum measured at the detector, and has to be accounted for in the analysis by a quantitative understanding of the plasma. Only the inhomogeneous part of the potential is of special interest, as the homogeneous part shifts the entire β-spectrum and can be easily taken into account.
In the past, the diffusion ansatz was used by L. Kuckert [7] to provide a quantitative model of the plasma. Recent studies indicate, that this ansatz may not include all significant plasma effects, such as the high energy tail of the spectrum, the transition between a collisional and noncollisional plasma, as well as all time dependent effects. Therefore, a new kinetic simulation is being developed within the KATRIN collaboration. A precondition for any theoretical plasma study is the knowledge of electron and ion distributions (and velocities) at any point within the tritium source. In the following sections, we will describe how these density and velocity distributions can be calculated.

Simulation approach
Electrons within the WGTS are generated through β-decay and through subsequent secondary processes. For each β-decay, this forms a Markov chain of events, with further branching event chains for secondary particles. This scenario is best simulated using the Monte Carlo technique.
The Monte Carlo algorithm was implemented in a newly developed code, named KARL. The code base for KARL was initiated by C. Reiling [8], but significant changes were made since then and erroneous calculations corrected. The algorithm of the new code version is described in section 3.1. The physics is outlined in section 4 and the first results presented in section 7.

Monte Carlo scheme
Each particle chain starts at the β-decay, where the position of the new electron and corresponding ion is sampled from the known neutral particle density distribution [13]. The energy is sampled from the Fermi spectrum of the β-decay [19]. The β-decay particles are then added to the particle queue and the simulation starts its iterative cycle, see figure 3. In each iteration step, a particle is selected at random from the particle queue. Then the mean free path, see equation (3.2), is determined from all possible interacting partner densities at the position of the particle and the energy of the particle. The particle is then moved according to the mean free path and the electromagnetic fields to a new position, see section 3.3. If the particle hits the simulation boundaries, it will be deleted. Otherwise, an interaction is selected randomly from the available interactions, see section 3.2. The interaction products are then added to the particle queue and the iteration starts from the beginning. The particle loop is repeated until there are no more particles in the particle queue. Then, new β-decay particles are created and processed until the specified target decay number is reached. The particles of the simulation can be used to determine particle densities during the simulation, see section 3.4. These densities can then be used to simulate electron-ion interactions.

Two particle scattering
In the simulation, particles interact with other particles. The probability of this interaction is given by the total cross section (TCS). In all cases presented here, it is dependent on the energy of the  . Monte Carlo scheme of the KARL code. β-electrons and ions are added to the particle queue. A particle is picked from the queue and processed. The iteration step either ends in a particle interaction or through termination at the simulation boundaries. New β-decay particles are added when the particle queue is empty. The tritium density is specified before the simulation. The particle densities of electrons and ions are determined during the simulation.
interacting partners. In the following, the TCS is always given in the inertial frame of reference, where the target particle is at rest. The probability for an interaction of a particle after passing a distance assuming a constant density of target particles T is given by Eq. (3.1) can be used to determine the average distance between two interactions, called mean free path. It calculates as This equation needs to be adapted if more than one interaction channel exists. For this reason, a new quantity, named the total macroscopic cross section (TMCS), is defined through where Σ ( ) = i ( ) · T,i is named the macroscopic cross section (MCS) of the interaction channel . T,i is the corresponding target particle density. The mean free path then calculates to rendering the TMCS a direct indicator for the mean free path. After interaction, the final state of a specific property of the particles can be important for the simulation. The statistical distribution of the final state is given by the differential cross section (DCS) d d . The probability for observing a particle within the parameter space between [ , ] is given by In the context of this work, the relevant properties that define the final state are the energy of the particles after collision and their scattering angle.

Reaction rates
The cross section of some reactions is not available from the literature. This occurs mostly for rare and low energetic interactions. In these cases, reaction rate coefficients can be used to provide a rough estimate of the cross section. For an unimolecular elementary step, the reaction rate is calculated by where is the initial density of particles. From a microscopic point of view, this rate can be identified with the collision rate = =v , wherev is the mean velocity of the interacting particle and the mean free path. The velocity can be calculated from the kinetic energy of the interacting particle. The mean free path is calculated using equation (3.2). In total, the cross section can be written as where is the mass of the interacting particle. Thus, the reaction rates coefficients can be used as an approximate estimate of the energy dependency of the cross section.

Charged particle motion
The particle motion is bound by the mean free path and the electromagnetic background fields in the source. The actual traveled distance is determined in the simulation through a pseudo random number ∈ (0, 1) and the mean free path taking into account the randomness of the real motion. The magnetic field is assumed to be temporally and spatially homogeneous. Therefore, the movement can be performed by the drift approximation [20]. In the drift approximation algorithm, the movement of the particle is described by the superposition of the movement of the guiding center and the gyro motion of the particle around the guiding center. The movement of the particle is subdivided into small spatial steps Δ to account for changes in the electromagnetic field during movement. The step size is chosen by the step size of the electromagnetic field data, to avoid movement across multiple electric field cells. The discretization of the spatial dimension can be translated into a discretization of the time domain to where |ì | is the absolute velocity of the particle before movement. The movement algorithm is performed in three steps. First, the position in longitudinal direction is calculated by where the index denotes the timestep at which the quantity is given. is the velocity of the particle parallel to the magnetic field. Second, the electric field during the motion is approximated by the mean of the field value at the original position and at the new position. This field is then used to calculate the movement perpendicular to the magnetic field and used to determine the new velocity +1 . The perpendicular particle position changes to where ⊥ is given by the so called E-cross-B drift [20] ì ⊥ = ì × ì 2 . (3.14) Third, the parallel velocity changes corresponding to the Lorentz force to where is the electric field strength parallel to the magnetic field, is the mass of the particle and its charge. The movement of a particle is stopped, either when the particle reaches an absorbing surface or if it crosses a density bin border, see section 3.4. If the particle movement is stopped at a density bin border, the movement step is performed again from the beginning. This method ensures a more precise description of particle movement through a cloud of particles with varying density. The absorbing surfaces for particle deletion depend on the particle type, see section 2.1. Ions, on the one hand, are deleted each time they hit the beam tube walls, the rear wall, or leave the source at the DPS side. The latter represents the filtering of the ions by the dipole electrodes in the DPS. Electrons, on the other hand, are deleted, when they hit the beam tube walls or the rear wall. They are reflected back to the simulation domain, when they cross the boundary towards the DPS. This procedure represents the reflection of the electrons both at the DPS dipole electrodes or the spectrometer potential. Backscattering of electrons and ions at the metal surfaces is not considered in the current version of the code. . Schematic representation of the source segmentation for density fields and particle current calculation. The source is segmented into evenly spaced segments in longitudinal, radial and azimuthal direction predefined by the user. The segment boundaries are used as predefined planes for the determination of particle currents, see section 3.5.

Secondary particles as target field
Electrons and ions mostly interact with the neutral gas in the source due to the large neutral particle density. These interactions depend solely on the primary particle and the neutral gas properties (density, velocity distribution). The gas density is high enough that it can be assumed that the neutral gas flow is not influenced by charged particle interactions [3]. Thus, each Monte Carlo event would be independent of the other, which is common for most Monte Carlo codes. Nevertheless, electrons and ions can interact with each other in the source, see section 4. Therefore, each Monte Carlo event must be dependent on the other events. This dependence is mimicked in the simulation by so-called density fields. These density fields determine the number density of each particle specie during the simulation. These fields are then used to determine the TMCS, see equation (3.3).
In the simulation, the source is segmented into evenly spaced segments in longitudinal, radial and azimuthal direction, predefined by the user, see figure 4 for a schematic representation. In each iteration step of the particle queue, it is recorded how long a particle with index stays in one of the predefined segments, labeled by the time . The number density of a species is then proportional to the total time of all particles of the corresponding species spent in the segment and the volume of the segment. In case of recombination, an additional term has to be added corresponding to the number of recombination of the recombination partner in the segment, resulting in where sim is called the simulated time (formula adapted from [8]). It describes the simulated physical time that is past after simulating a fixed number of beta decays . It can therefore be directly determined by the activity of the source through sim = . (3.17) The activity of the source is not dependent on the interaction of the charged particles in the source, and can therefore be evaluated from the predefined amount of neutral gas in the source. The number density can also be determined depending on the energy of the particle. These number densities are then binned to a predefined energy scale, resulting in the spectrum of the particles. Thus, the density fields can be used in the calculation of the TMCSs, see equation (3.3), and also used for sampling a new particle with the simulated spectrum. The sampled particles can then be used in the interaction step, see section 4, to evaluate the outcome of a particle interaction of two different types of particles. Therefore, these new sampled particles are called interaction partners.
The sampling of an interaction partner is performed in three different steps: sample particle energy, roll velocity direction, copy position of primary particle. First, the energy of the interaction partner is sampled from the energy distribution of the corresponding density field at the position of the primary particle. Secondly, a velocity direction is picked at random. The velocity of the new particle is then scaled according to the energy and mass of the particle. Finally, the new particle is generated with this velocity and the position of the primary particle.

Particle tracking and diagnostics
The Monte Carlo approach of the simulation needs high enough particle statistics. The statistic is determined by the number of generated β-decay events , specified by the user. Each of these events produces multiple secondary particles. Recording of each position and particle would result in too much data. Thus, it was decided on three different output methods: density fields, see previous section 3.4, absolute particle numbers of conversion and termination, and particle currents through virtual barriers.
The density fields with number densities of the previous section are recorded each time after a specified number of generated β-decay events. The last output is performed after decays. This output is used for the final determination of the particle number densities of the simulation, which is one of the main objectives of the code. The other outputs are only used for diagnostics of the code.
Special counters are implemented to track each time a particle is terminated at the simulation boundary, converted to another particle species or recombines. The termination counter is selective on the boundary type: rear wall, tube wall or DPS. These counters can be used for diagnostics of the code or for tests of the code.
The second main objective of the code is to determine the particle currents in the source. This task is performed through so-called virtual barriers: predefined planes within the source, where each particle is registered, when it crosses from one side to another. The total current through the segment faces is then given by The current through the faces can be determined depending on the energy and type of species crossing the faces. This is performed through a selective binning to a predefined energy scale during the simulation. The absolute current is then given by the sum over all energy bins.
The accuracy of the simulation depends on the number of decays . This number is inserted directly in the calculation of the current and the density, but it also contributes indirectly through the production of secondary particles. The specific contributions are difficult to calculate directly because of the multitude of different interactions. Hence, the accuracy of the simulation was determined through multiple simulations with the same number of decays but with different random number generator seeds. The mean of the normalized standard deviation of each energy bin was calculated. For a counting experiment, it is expected that the standard deviation scales with 1/ √︁ , which was found to be in good agreement with our test simulation results. To reduce the statistical uncertainty to 0.1 %, it was decided to use 1 × 10 6 β-decays in the following simulations.

Physical processes
In this section, the most relevant physical processes for charged particles in the source are described. The focus lies here on the implementation in the code. For the interactions, the cross sections were calculated by analytical formulas or through numerical tables. The corresponding parameters and values were obtained from the literature. Initially, this collection was compiled by C. Reiling [8].
Due to the nature of the Monte Carlo approach, uncertainties of the cross sections can only be considered in the code through an evaluation of multiple simulations with varied cross sections, which will be done in the future. The interested reader is referred to the original publications for the evaluation of the uncertainties for the parameters.

Particle injection
The injection of particles mimics the β-decay of the tritium molecules in the source into an electron eand an ion molecule T He + . In the code, the ion is described as a T + 2 ion. Isotopic differences are assumed to be small [5]. Additionally, due to charge transfer processes from T He + to T T + and dissociation, the lifetime of T He + is negligible in this context [21]. The position of the new β-electron and ion is sampled from the known density distribution of the neutral gas through the acceptance-rejection method [22]. The energy of the β-electron is sampled through the Fermi distribution of β-decay [19] in conjunction with the acceptance-rejection method. The velocity direction is assumed to be isotropically distributed. The energy of the ion is sampled from a Maxwell-Boltzmann distribution with an additive drift velocity of the neutral gas. The recoil of the ion is neglected due to the large mass difference.

Interactions of electrons with tritium molecules
β-electrons are created in the source with energies of up to 18.6 keV. In contrast, the molecular tritium gas has a temperature of 80 K, which corresponds to energies of 6.9 meV. Thus, e --T 2 interactions need to be described on a large energy range. Density-wise, the WGTS is dominated by molecular tritium. Hence, interactions with molecular tritium are very probable, and the possible interactions need to be modeled precisely.
A collection of the implemented models for e --T 2 interactions is shown in the following sections, see also figure 5. There is almost no data available for interactions with tritium. Isotopic differences of the interaction probability are assumed to be small [5]. Therefore, models of e --H 2 interactions are adopted in the code.   [23]. Data for ionization and electronic excitation taken from [24]. Data of elastic scattering from [23], simulation using ELSEPA [25] and extrapolation, see below. Figure based on [8].

e --T 2 elastic scattering
Elastic scattering with molecular tritium is the dominant interaction channel for low energy electrons ( < 60 eV), see also figure 5. For energies above 60 eV, other interactions, such as ionization and electronic excitation, become dominant. Nevertheless, elastic scattering has still an effect on the particle motion. There was no model available for the cross section of elastic scattering on the large energy range. Therefore, the TCS was derived from three different sources, see also figure 6. The majority of the energy range (0.02 eV < < 100 eV) was covered by the experimental studies from [23] and implemented as an analytical function. The cross section for higher energies > 100 eV were derived from the ELSEPA simulation code [25] and implemented as numerical tables. No experimental or theoretical cross sections were found for energies below < 0.02 eV. The shape of the cross section of [23] indicates a continuous trend towards lower energies. Thus, the cross sections of Yoon was extrapolated towards lower energies as a first approximation, see also figure 7. The analytical expression of this extrapolation (from [8]) reads as The final state of the elastically scattered electrons depends on the two scattering angles and . The orientation of the tritium molecules is assumed to be isotropically distributed. Thus, Cross section (cm 2 ) extrapolated experiment calculated Figure 6. Visualization of the composition of the cross section of e --T 2 elastic scattering. The cross section is derived from three models. Cross sections for energies > 100 eV were calculated using ELSEPA [25]. Cross sections for energies in the range 0.02 eV < < 100 eV were adapted from experimental data from [23]. Cross sections for energies < 0.02 eV are extrapolated values from experimental data from [23].   the polar angle is sampled from the interval [0, 2 ). The distribution of the azimuthal angle is determined through the use of the ELSEPA code [25]. The simulations showed, that the azimuthal angle can be assumed to be isotropically distributed for energies below < 10 eV. For higher energies, the azimuthal angle is sampled from the tabulated data of the ELSEPA calculations.
The position of the electron does not change in the scattering process. Nevertheless, the guiding center of the gyro motion changes its positions. Thus, this change of the guiding center is reflected in the simulation. The new position of the guiding center is shifted by the length Δ with  where is the projected scattering angle to the magnetic field direction. The electron performs many gyrations between interactions. Therefore, the polar angle of the gyration can be considered as randomly distributed. Thus, the shift is applied randomly in the x and y directions.

Electron impact ionization of T 2
Most of the electrons and ions in the source are not created through β-decay, but through electron impact ionization. There exist multiple different ionization channels, which mostly depend on the state of the tritium molecule before, during and after ionization [24]. Because the tritium gas is cooled down to 80 K, it is assumed here that the molecule is in the ground state, prior to the interaction. Similarly, it is assumed that the maximum reachable state of the molecule is the first excited state. In this case, the study of Janev et al. [24] provide analytical formula for three different ionization channels: non-dissociative ionization, and dissociative ionization over the ground and first excited state of the ion, see also figure 8. The most probable channel is non-dissociative ionization e -+ T 2 e -+ T + 2 + e -, with a threshold energy of 15.42 eV. Each electron creates an additional electron and a tritium ion molecule. In dissociative ionization, first an electron is removed from the molecule, which then resides in an excited state. This state then decays into a tritium ion and a tritium atom e -+ T 2 e -+ T +* 2 + e -T + + T + 2 e -, with threshold energies of 18.15 eV and 30.6 eV for the ground and the first excited state respectively. In the simulation, the energy of the ions is assumed to follow the distribution of the gas flow. Therefore, the energy of the ion is sampled from the Maxwell-Boltzmann distribution of the neutral gas. The energy of the two final electrons is correlated. In the code, the energy of the first electron is sampled from the analytical formula of the DCS from [26]. The energy of the second electron is then calculated correspondingly as the difference between the initial energy, the ionization energy and the final energy of the first electron. The distribution of the azimuthal and polar angle of both electrons is adopted from [27].

Electronic excitation of T 2 by e -
Tritium can be electrically excited by collision with electrons. The energy difference of the levels is on the order of O (10 eV). Thus, high energy electrons > 100 eV can lose a significant amount of their energy due to excitation. Similar to electron impact ionization, the cross section of the interaction depends on the initial and final state of the molecule. Again, it is assumed, that the molecule is in the ground state prior to the interaction and that the maximal state of the molecule is the first excited state. In this case, the study of Janev et al. [24] provides three main categories of electronic excitation: • Excitation to dipole allowed singlet state • Excitation to dipole forbidden singlet state

• Excitation to triplet states
Each of the categories can be subdivided by additional quantum numbers of the molecule. In the code, the cross section of the most dominant interactions from [24] were implemented as an analytical formula. This way, 15 different excitation channels were resolved.
In the simulation, it is assumed that the tritium molecule does not gain any additional kinetic energy due to the large mass difference. Thus, the electron only loses the energy of the excitation process. The direction of motion is therefore also left unchanged.

Rotational and vibrational excitation of T 2 by e -
The energies of the rotational and vibrational excitation levels of T 2 and H 2 are assumed to be different, due to the significant mass difference of the isotopes. Studies of Yoon et al. [23] and Janev et al. [24] show, that rotational excitation and vibrational excitation have a non-negligible effect on the total cross section. Therefore, in the code, the cross sections of H 2 are used as a first approximation for rotational and vibrational excitation of T 2 . The cross sections are implemented in the code as analytical formula, where the data of rotational excitation is taken from [23] and the data of vibrational excitation is taken from [24].
Similar to the electronic excitation, it is assumed that the tritium molecule does not gain any additional kinetic energy due to the large mass difference. Thus, the electron only looses the energy of the excitation process and the direction of motion is left unchanged.

Recombination of electrons with tritium ions
The direct interaction between ions and electrons is improbable, due to the expected low densities of electrons and ions ( ≈ ≈ 10 6 cm −3 ) [7]. Therefore, only recombination is considered in the code, because this reaction changes the total number of electrons and ions.
In the code, three different ion species are distinguished, see section 4.3. Therefore, only recombination processes with these species are discussed here. There was no data found on the recombination cross sections of tritium. Therefore, the cross section of recombination of tritium is approximated by the cross section of hydrogen.
Atomic ions can directly recombine with electrons through radiative recombination, which produces a neutral particle and a photon. The cross section of this interaction is low [28] in comparison to the other recombination channels, see also figure 9, and it is therefore not considered further in the code.
Tritium molecule ions can undergo two-body recombination and three-body recombination with additional electrons. Matveyev et al. [30] derived a mean rate coefficient 3,rec = 4.7 × 10 −31 m 6 s −1 , (4.5) for three-body recombination at 80 K. The mean free path of the reaction, assuming even a thousand times higher electron density of ≈ 1 × 10 9 cm −3 , is calculated to be on the order of O (10 3 m), which is significantly larger than the length of the source ( WGTS = 16 m [3]). Therefore, three-body recombination is not considered further. Janev et al. [24] show that most of the tritium ions undergo dissociative recombination through In the code, the cross section of recombination is implemented directly from the analytic formulas of Janev et al. [24]. The cross section data is depicted in figure 9.
In the code, the recombination process is implemented as a deletion of the impact particle. The reduction of the density of the corresponding interaction partner is handled by the increase of the additional term in equation (3.16). Neutral particles are not tracked in the simulation, and are therefore not considered further.

Interactions of tritium ions
The elastic scattering probability of ions with neutral gas is much larger than for electrons [24]. Thus, they scatter more often and will most likely adopt the velocity distribution of the neutral gas. Therefore, the energy distribution of ions is not tracked in the simulation, but approximated as a Maxwell-Boltzmann distribution. Only the absolute density is recorded. In the same manner, only ion-gas interactions are implemented, which are relevant at thermal energies.
In the code, three different cluster sizes are distinguished: T + , T + 2 and T + 3 . Higher order clusters are neglected here, due to missing cross section data or reaction rates. Similar to e --T 2 interactions, there is no sufficient data on interactions of tritium ions. Models of H ions are adopted throughout this section.

Elastic Scattering with T 2
Tritium ions T + and T + 3 can scatter elastically from the neutral gas. At low energies, elastic scattering is the dominant interaction channel. Thus, it is assumed that the ions in the source reach thermal equilibrium quickly. In the code, the elastic cross sections from [31] are implemented as analytical formulas.
T + 2 does not scatter elastically with the neutral gas. Nevertheless, these ions perform charge exchange processes with the neutral gas. This way, the T + 2 also adopts the velocity distribution of the neutral gas. In the code, the cross section of charge exchange from [24] was implemented as an analytical formula.
The kinematic of elastic scattering is treated similar to the one of scattering of electrons with the neutral gas. Also, the guiding center position of the ions is changed in the interaction. The scattering angle is assumed to be isotropically distributed due to the relatively small collision energy.
In the charge exchange process, a new ion is created with a velocity sampled from the velocity distribution of the neutral gas. The position of the initial ion is adopted for the new ion, with the additional shift of the guiding center. The direction of motion is picked at random.

Ion cluster formation
Data from [32] shows that tritium cluster ions exist at 80 K. The formation of these clusters is also represented in the code. As already mentioned, three different cluster sizes are distinguished in the code: T + , T + 2 and T + 3 . T + can perform radiative clustering or ternary association clustering [33]. In the radiative clustering process, the binding energy is released by emitting a photon. In the ternary association clustering process, the additional energy is transferred to a second tritium molecule. There was no data found on the cross section of the clustering processes. Thus, data from the rate coefficient from [33] and [32] was used to calculate a rough estimate of the cross section through equation (3.8).
T + 2 can form T + 3 through collsion with T 2 either by transfering one atom of the T 2 molecule to the T + 2 ion or by transfering a T + ion from the T + 2 ion towards the T 2 molecule (directly adapted from proton and atom transfer of hydrogen interactions) [24]. The sum of the cross sections of both processes is provided by [24] and implemented as an analytical formula in the code.
In the code, the new ion is created at the same position as the original ion. There was no data available on the differential cross section. Therefore, as a first approximation, the new ions were created with the same direction of motion and the same kinetic energy as the initial ion.

Validation
The code was tested thoroughly, including a comparison of specific simulations with expectations. These expectations were formulated for scenarios where an analytical expression of the simulation output can be derived. For these simulations, each interaction was investigated detached from the others. The corresponding tests are summarized in the following.
Elastic scattering Elastic scattering was tested through an injection of monoenergetic particles into a neutral gas reservoir with a Maxwell-Boltzmann distributed velocity. If elastic scattering is implemented correctly, then these particles will adopt the velocity distribution of the neutral gas independent of the initial energy. The tests show that the particles indeed adopt the velocity distribution of the neutral gas. This result was found to be independent on the initial energy of the particles.
Ionization Ionization was tested through an injection of mononenergetic electrons. The energy of the electrons was chosen in such a way that only one ionization could occur due to the energy loss in the ionization process. The energy distribution of the ionization products were compared to the expected energy distribution. An agreement was found within the limits of the statistical uncertainty. The ionization process can occur either through a dissociative or non-dissociative process. Thus, the number of T + ions in relation to the T + 2 ions should match the relation between the cross sections of the reactions. The numbers of the particles were compared and found to be in agreement with the expectations.
Excitation Excitation was tested through an injection of monoenergetic electrons. In the excitation processes the electrons lose energy depending on the excitation state while keeping all other properties, see section 4.2 and 4.2. Therefore, a comparison of the initial and final energy will result in the energy loss through excitation.
Recombination was tested through an injection of monoenergetic ions and electrons at one end of the simulation with the energies and respectively. All other interactions were disabled and a magnetic field in the -direction was applied. Thus, there is no dependence on the x and y position. The expected density of the ions ( , ) and electrons ( , ) can be calculated from the coupled continuity equation of electrons and ions with a particle sink. The differential equations in this case then read as In the limit of large times it is expected that the density reaches equilibrium. So this coupled differential equation can be solved by the linear equation where is an integration constant. A linear relation between the electron and ion density is expected.
The slope of the test simulation can now be compared to the expectations from the theoretical cross sections. It was found that both values are in good agreement. Additionally, it was recorded how many electrons and ions recombine in the simulation. In the case where electrons and ions have the same velocity, they induce the same density in a simulation without recombination. Thus, in a simulation with recombination it is expected that the number of simulated electrons, which recombine is similar to the number of the simulated ions, which recombine. In a test simulation, it could be shown that the two numbers agree well. In total, it could be shown that the treatment of electron-ion interaction through density fields works as expected.

Physical environment
Neutral gas density The density distribution inside the source is generated by the injection of neutral gas in the center and subsequent pumping of the gas at the pump ports. Thus, the density distribution is shaped by the underlying geometry of the setup, the injection pressure and the pumping rate. A corresponding model was developed by Kuckert et al. [13]. The calculations of Kuckert et al. were performed at the design conditions of the source ( = 30 K and in = 0.337 Pa [1]). The current measurements of the KATRIN experiments are now performed at higher temperatures and higher injection pressure ( = 80 K and in = 0.7 Pa). This results in a lower injection density of in ≈ 6.34 × 10 14 cm −3 , using the ideal gas law. It is assumed that the density profile does not change significantly, and the density can be scaled accordingly to the inlet pressure. The resulting density distribution can be found in figure 10. The calculation of Kuckert et al. were only performed in the direction of the detector. Due to the symmetry of the setup, it can be assumed that the density is distributed similarly in the direction of the rear wall.
A simplified density profile was used for the Monte Carlo simulation. Here, the density is described by six linear functions, see figure 10.

Neutral gas velocity
The velocity distribution of the neutral particles is influenced by the inlet pressure and temperature, the pumping speed, and the temperature of the beam tube walls. Each contribution is weighted differently according to the density of the neutral particles. In the center, for example, there are many molecule-molecule interactions. Here, the velocity distribution is dominated by the particle flow. In the rear wall chamber, the density is reduced drastically, such that molecule-molecule interactions become negligible and collisions with the tube walls are dominant.
In the viscous flow regime, the velocity distribution can be described by a Maxwell-Boltzmann distribution with an additional bulk velocity in the longitudinal direction [7]. The temperature of this distribution is given by the temperature of the beam tube walls. The bulk velocity describes the streaming process from the inlet in the center of the WGTS to both sides. Simulations of Kuckert et al. [13] show that the bulk velocity is dependent on the longitudinal as well as the radial position. The bulk velocity increases from the inlet to both sides. Radially, the bulk velocity is reduced at the tube walls, corresponding to wall interactions which slow down the gas movement. At lower densities, the velocity can no longer be described through a drifting Boltzmann distribution, because interparticle collisions are reduced. In first order, the particles follow a Maxwell-Boltzmann distribution, with a temperature of the beam tube walls.

Magnetic field
The magnetic field inside the source is generated by nine super conducting magnets. These magnets surround each tube element of the source. These magnets provide a strong field and are set to 2.5 T with a low drift < 0.03 % per month [14] and fluctuations below 0.2 % per data acquisition run [3]. For the simulation, the magnetic field is assumed to be homogeneous throughout the source.

Results
The results presented below were generated through a simulation with a maximal tritium density of 6.3 × 10 14 cm −3 at a temperature of 80 K. The source tube was approximated by a cylinder with a length of 16 m and a radius of 4.5 cm, simplifying the segmentation of the source, see section 3.4. The injection of tritium gas is assumed to be at the center of the source tube; 8 m from the rear wall. No electric background field was used, and the magnetic field was set to 2.5 T. 10 6 β-decays were simulated, resulting in a statistical uncertainty of 0.1 %, see section 3.5. This uncertainty will be used for all simulation results presented below. The source was segmented into 32 z bins and 15 radial bins. No segmentation was used in the azimuthal direction due to the expected azimuthal symmetry of the result. Density and currents were recorded at 100 energy bins with a logarithmic subdivision.
Electron spectrum The electron spectrum was recorded at 32 different positions in the source. Three distinct positions are compared here: a segment, which is directly connected to the tritium inlet; a segment, which is next to the rear wall; and a segment, which is connected to the DPS. The corresponding spectra can be found in figure 11. Three different energy regions can be identified (distinguished by vertical lines in the figure): a so-called thermal region ( < 6 × 10 −2 eV), where the electrons reached thermal equilibrium with the neutral gas, a so called β-decay region ( > 10 2 eV), where the spectrum is governed by the Fermi distribution of β-decay and a so-called secondary region in between the other two regions. The spectrum in the β-decay region shows little deviation between the three different positions. Therefore, the β-electrons either travel through the source unobstructed, or lose the main portion of their energy through ionization and subsequent scattering. There is a slight difference at the lower end of the β-region between the rear wall and the other two positions. Electrons which reach the rear wall segment can have traveled at maximum two times through the source through reflection at the spectrometer, while electrons at the DPS traversed the source only once. Therefore, more β-electrons have scattered before reaching the rear wall segment, which reduces the density in this energy range.
The spectral shape in the thermal region is similar for each segment. These electrons have scattered multiple times with the neutral gas and thus adopted the temperature of the gas. The shape can be described by the Maxwell-Boltzmann distribution [34], also depicted in figure 11, scaled by the factor √ taking into account the logarithmic binning and the speed of particles leaving the simulation domain. Differences at the low energy side arise from the reduced cross section of elastic scattering at lower energies, see figure 5. Particles with lower energies scatter less, increasing their time spent in the simulation domain, which in turn manifests in a higher density, compare equation (3.16). This behavior was verified in targeted test simulations with fixed cross section, which resulted in a good agreement of the Maxwell-Boltzmann distribution with the simulated data in the lower energy regime down to 10 −6 eV.
The similarity to the Maxwell-Boltzmann distribution also shows that recombination is only a secondary effect. Otherwise, there would be much fewer particles at low energies, where the cross section of recombination is increased significantly, see figure 9.
The absolute density is the largest in the center, while it decreases towards both sides. In the center, there are many collisions. Thus, the electron density is increased here. The density at the DPS is higher than at the rear wall. This is caused by the reflection of electrons at the DPS, see section 2.
The secondary region shows two distinct features. A continuous decrease from 5 × 10 −2 eV to 10 eV and an edge at 10 eV. The energy of the edge is consistent with the threshold energy for ionization and electronic excitation. Electrons with a lower energy can no longer take part in these interactions. The continuous decrease arises as a combination between elastic scattering and rotational and vibrational excitation. The rotational excitation has a threshold energy of approximately 5 × 10 −2 eV, which is consistent with the end of the secondary region, see figure 5.
Particle densities The density of each particle species can be derived from the sum over the density spectrum. The corresponding distribution along the source tube can be found in figure 12. It can be seen that the density profile of each specie is different from each other. The absolute charge density is non-zero in this simulation, despite the fact that electrons and ions are created in pairs. This is not surprising because electrons scatter less and are moving faster, and therefore leave the source tube faster than the ions. Therefore, their density is reduced in comparison to the ions. In the experiment, it can be assumed to first order that there exists quasi neutrality in the source. Thus, the density differences between electrons and ions will induce an electric field. This field will be calculated in a different dedicated simulation.
The electron density is dominated by the thermal part of the spectrum. β and secondary electrons make up only a maximum of 4 × 10 2 cm −3 . The electron distribution shows a maximum in the center of the source and falls off to both sides. The probability of scattering is the highest in the center of the source. Hence, the electron density is the largest here. The scattering probability is decreased to both sides as the neutral gas density decreases. The electron density is asymmetric towards the inlet, despite the symmetric shape of the density. This is caused by the reflection of electrons at the DPS, which increases the density in this region. Nevertheless, the density before the DPS is not increased drastically, which suggests that electrons can pass the high density region of the inlet.
The ion distribution is different for the three different ion species. The most abundant ions are T + 3 ions. This is expected due to the high conversion cross section of T + and T + 2 towards T + 3 . The density of the T + 3 ions is governed by ionization through electron impact and by elastic scattering with the neutral gas. Ionization is most probable in the center, where the neutral gas density is the highest. The resulting density of ionization is superposed by the streaming process, caused by the elastic scattering. The motion of ions is aligned with the bulk velocity of the gas. Thus, the ions are transported to both sides of the source, increasing the density in these regions. The density drops 0 2 4 6 8 10 12 14 16 Distance to rear wall z (m) in the region in front of the rear wall and in front of the DPS. This is caused by the reduction of the neutral gas density and the drop out of the neutral gas drift. The density of T + is distributed similar to the density of T + 3 . The shape of the spectrum is created by the same effects as for T + 3 , namely elastic scattering with neutral gas and particle transport. The density of T + is lower than the density of T + 3 . In this simulation, T + can only be created through dissociative ionization. The cross section of this reaction is significantly lower than for non-dissociative ionization. Therefore, less T + is created in comparison to T + 3 , which is created through the efficient conversion of T + 2 . Additionally, T + can be converted to T + 3 through ternary association, which lowers the density of T + in comparison to the density of T + 3 . The density of T + 2 is significantly different from the other two ion species. The shape is almost constant, and the density is significantly lower. This behavior can be explained by the combination of two effects. First, the cross section of cluster formation is higher than the cross section for charge transfer [24]. Therefore, T + 2 ions exist only for a short time between creation and their first interaction. Thus, the density of T + 2 is significantly lower. Secondly, the probability for cluster formation is dependent on the tritium density, while the creation mechanism through beta decay and ionization is also dependent on the tritium density. In total, both effects cancel out and the T + 2 density therefore shows no large longitudinal variation. Similar to T + 3 and T + , there is a slight decrease of the density in front of the rear wall and the DPS. Again, the mean free path of the particles has increased enough that some particles can leave the source, which reduces the density.
Particle currents The particle current was determined through 33 virtual planes evenly spaced longitudinally. The planes were aligned in such a way that the first and last longitudinal plane coincide with the rear wall and the passage to the DPS. The corresponding longitudinal current profile can be found in figure 13. The shape of the electron current coincides with expectations. Electrons are reflected at the DPS by the potential of the dipole electrodes and the spectrometer, see section 2. The corresponding current is therefore zero. Electrons are most likely to be terminated at the rear wall. Thus, the electron current is the highest here. The electron current decreases continuously from the rear wall towards the DPS. Hence, the electrons can stream from the upstream side towards the downstream side and are not hindered from passing by the neutral gas flow.
The ions are created in the center and stream to both sides. There is no driving force of the ions towards the other side of the inlet. Thus, the ion current at the inlet (8 m in the simplified geometry) is zero. The current leaving towards the DPS and towards the rear wall is almost identical. The differences can be explained by the increase of β-electrons through reflection at the DPS side.
The current of both electrons and ions reaches a plateau between 0 m to 2 m, and 14 m to 16 m. This behavior is directly linked to the reduced neutral gas density in these regions, see figure 10. A negligible amount of ions and electrons are created here through beta decay and ionization. Thus, the current in these regions is a result of the creation of electrons and ions in the center of the source and the subsequent magnetic guiding towards the plateau regions. Hence, the current does not increase further in these regions.
The sum of the ion current leaving the source at the DPS and at the rear wall is below the electron current leaving at the rear wall. This can only be explained through an ion current, which is directed towards the beam tube. The corresponding current, displayed in figure 14, was detected 0 2 4 6 8 10 12 14 16 Distance to rear wall z (m) in the simulation through radial virtual planes located at the beam tube walls. This current is caused by the position change of the guiding center after elastic scattering. The position change is dependent on the Larmor radius of the particle, see equation (4.2), which in turn scales with the square root of the mass. Additionally, the cross section of ion scattering is much higher than for electron scattering. Therefore, more collisions occur, which can shift the guiding center, increasing the maximal position change in radial direction [20]. Because of the larger position change and larger scattering probability, the radial ion current is expected to be significantly higher than the radial electron current, which is also reflected in the simulation results.
The absolute current to the tube is also influenced by the density profile of electrons and ions. It can be seen that the shape of the current towards the beam tube is similar to the shape of the neutral particle density. This is caused by the increased elastic scattering probability in the center due to the increased neutral gas density. The T + current is below the T + 3 current, which is caused by the density difference and the different Larmor radius of the two ion species. The current of T + is of the same order of magnitude as that for the electrons, even though the masses of both species is significantly different. This behavior is caused by the significant difference in the density.
There is a contribution of T +

Discussion
The latest results of the electron spectrum in the KATRIN source was published by Nastoyashchii et al. [6] in 2005. They derived an electron spectrum of the KATRIN source through a Monte Carlo simulation. The spectrum was derived through an evaluation of the total time particles spent in the entirety of the source. Position dependence of the spectrum was not included, and therefore the spectrum can only be construed in first order as a mean spectrum of the electrons in the source. In this paper, we derived electron spectra at different positions in the source, see figure 11. A comparison of our results to Nastoyashchii et al. is depicted in figure 15. Nastoyashchii et al. only derived a relative spectrum and did not specify the absolute scale. Here, we assume, that the high energy region of their and our simulation coincide. Thus, their spectrum is scaled accordingly, such that the data point of Nastoyashchii et al. at 6886 eV has the same value as the energy bin from 5447 eV to 7060 eV of our calculation from the central position. It can be seen that the distribution of Nastoyashchii et al. shows similar features as our simulation: a thermal region, a β-region and a secondary region. However, there exist significant differences, which will be discussed in the following.
The results of Nastoyashchii et al. only provide a mean particle density of the source. It can be seen in our simulations that there exist significant differences in the spectrum at the different positions in the source. These differences might change the description of the plasma. Nastoyashchii et al., for example, assume implicitly in their affiliated calculation of the plasma potential that the high energy electrons of the spectrum can be neglected and it can solely be described by a thermal distribution (equation (1) in their publication). However, our simulations show that this assumption is not true for every position in the source, especially when considering that the maximum energy density in the beta region ( ∼ 10 eVcm −3 ) is higher than the energy density of the thermal region in front of the rear wall ( ∼ 1 eV cm −3 ). Thus, a different plasma description is necessary, which must also include high energy electrons.
There exists a spike (one data point) at 10 eV in the results of Nastoyashchii et al., which does not exist in our simulation, see figure 15. The existence of this peak is important, because there might be plasma instabilities linked to it (Penrose criterion [35]). Nastoyashchii et al. do not specify the origin of this peak and provide no information on the statistics of the simulation. We did not find any physical explanation of the peak. To test whether statistical fluctuations could produce such a peak, we performed our simulations with different random number seeds. No peak was found.
Although the thermal region of Nastoyashchii et al. seems to have the same shape and position as our calculations, there is a difference between both. For comparison, a least squares fit of a Maxwell-Boltzmann distribution [34], scaled by the factor √ taking into account the logarithmic binning and the speed of particles leaving the simulation domain, was performed on both datasets in the energy range of 1.3 × 10 −4 eV to 2.0 × 10 −2 eV, resulting in Nast. = 37.2 K and Karl = 27.7 K. The corresponding distributions are also depicted in figure 15. The expected temperature is 30 K, the same as the temperature of the tritium gas. The slightly lower temperature and increased density in the lower energy regime of our calculation can be retraced to the energy dependency of the elastic cross section, confirmed by test simulations, compare to section 7. The result of Nastoyashchii et al. shows both a higher temperature and lower densities in the lower energy regime than expected. It is unclear, which mechanism could cause this behavior, and Nastoyashchii et al. did not specify, where the discrepancy comes from.
The simulation of Nastoyashchii et al. were performed at fixed source conditions ( = 30 K and 100% 50% 10% Figure 16. Comparison of the electron spectrum at different inlet densities. 100 % correspond to an inlet density of 6.3 × 10 14 cm −3 . The source temperature was set to = 80 K and the magnetic field strength to = 2.5 T. tube walls, which is caused by the gyro center movement through elastic scattering, see figure 14. This is an improvement over the calculation of Nastoyashchii et al., which did not calculate particle currents in their simulation, especially not in radial direction. Our simulations can also be used to investigate the behavior of this current, especially under a change of the source magnetic field strength, see figure 17. It can be seen that the calculated current towards the beam tube is dependent on the source magnetic field. This behavior is expected, because the gyro radius of the particles is increased at lower magnetic field strength, and thus increasing the position change per collision, see equation (4.2). The simulation with varied magnetic field also shows significant differences in the radial particle densities and particle currents, see figure 17. These densities and currents can now be used to investigate the radial profile of the potential in the source through dedicated field solvers. This theoretical investigation can be supplemented by measurements with the radially resolved detector [3]. Thus, an experimental investigation of the simulation results will be possible.

Conclusion and outlook
In this paper, we presented the first results of a new Monte Carlo code, which simulates the density and energy distributions of charged particles in the KATRIN source. The distributions were obtained through tracking single particles in the source until they either recombine or are absorbed by the metal surfaces. The single particles can interact indirectly with each other through density fields, 2.5 T 2 T 1 T Figure 17. Comparison of the radially resolved ion density and current density at different magnetic field strength at 5 m distance to the rear wall. The source temperature was set to = 80 K and the inlet density to 6.3 × 10 14 cm −3 . see section 3.4. Therefore, this approach allows for the simulation of electron-ion interactions alongside electron-neutral and ion-neutral interactions. The simulated electron spectra showed three regions: a β-region, a thermal region and a secondary region. The relative density of these regions was found to be dependent on the position in the source, see figure 11. We also determined the total density of electrons and ions depending on the position in the source, see figure 12. Significant differences were found between the densities of the different species. Additionally, we determined the particle current depending on the position in the source through virtual barriers, see figures 13 and 14. A non-negligible ion current towards the beam tube was found.
Our simulations are a significant improvement to the previous simulation results of Nastoyashchii et al. [6]. Notably, the electron spectrum can be determined depending on the position in the source. Additionally, we can now investigate the particle behavior at various source conditions, which can be directly compared to measurement data. The source conditions include among others the magnetic field strength, the gas temperature and the neutral particle density. Investigations of the influence of the interactions and their cross section can also be performed.
The simulations presented in this paper were performed without an electric field in the source. This field arises in the experiment from surface potential differences and through the plasma generated by electrons and ions. The results of this paper can provide a starting point of the electric field calculations. The interplay of electric field and particle density is very intricate due to the nature of plasmas and demands special treatment in the future.