Time dependent simulation of the flow reduction of D$_2$ and T$_2$ in the KATRIN experiment

The KArlsruhe TRItium Neutrino experiment (KATRIN) aims to measure the effective electron anti-neutrino mass with an unprecedented sensitivity of 0.2 eV/c$^2$, using beta-electrons from tritium decay. Super-conducting magnets will guide the electrons through a vacuum beam-line from the windowless gaseous tritium source through differential and cryogenic pumping sections to a high resolution spectrometer. At the same time tritium gas has to be prevented from entering the spectrometer. Therefore, the pumping sections have to reduce the tritium flow by at least 14 orders of magnitude. This paper describes various simulation methods in the molecular flow regime used to determine the expected gas flow reduction in the pumping sections for deuterium (commissioning runs) and for radioactive tritium. Simulations with MolFlow+ and with an analytical model are compared with each other, and with the stringent requirements of the KATRIN experiment.


Introduction
The KArlsruhe TRItium Neutrino experiment (KATRIN) has been designed to determine the effective mass of electron anti-neutrinos with an unprecedented sensitivity of 200 meV/c 2 at 90% CL, using electrons from tritium β-decay [1,2]. The analysis is focused on the last few eV below the 18.6 keV endpoint of the β-spectrum.
The experiment is located at the Karlsruhe Institute of Technology (KIT), Campus North, near Karlsruhe, Germany. The approximately 70-m long beamline is shown in Fig. 1. The experiment can be divided into two main sections. The Source and Transport Section is responsible for the production and adiabatic transport of tritium β-particles to the Spectrometer and Detector Section. Their energy is determined in the integrating, electrostatic Main Spectrometer, which can provide high energy resolution with a wide open solid angle acceptance for β-electrons, emitted isotropically in the windowless gaseous tritium source (WGTS). The expected signal rate in the last eV of the β-spectrum is about 10 −2 counts per second; thus, a necessary precondition of reaching the neutrino mass sensitivity goal is a similarly low background rate [1]. This stringent requirement entails a thorough understanding and mitigation of background sources along the entire beamline of the experiment. Of particular importance are background electrons produced in the main spectrometer. One such background source are β-decays of tritium passing from the WGTS through the connecting beam-line into the spectrometer.
The WGTS has been designed to produce more than 10 11 β-particles per second [3]. A constant flow of 95% pure T 2 gas is injected at the center of the 10-m long beam-tube with a pressure of 3.4 × 10 −3 mbar at a temperature of 30 K [4]. Superconducting magnets adiabatically guide half the emitted electrons through the differential [5] and cryogenic [6,7] pumping sections towards the spectrometer section, while reducing the tritium flow by at least 14 orders of magnitude [1].
Over 99% of the gas is already pumped out by the turbo-molecular pumps (TMP) of the first stages of the differential pumping systems, which are integrated at both ends of the source cryostat (DPS1-R and DPS1-F). The detailed rarefied gas flow simulation has been described by Kuckert et al. [4]. In the remaining text, WGTS refers to the entire source cryostat, including the DPS1-R and DPS1-F subsections. In this paper, we describe the simulation in the molecular flow regime of the flow reduction factors for tritium in the second stage of the differential pumping section (DPS2-F), and in the cryogenic pumping section (CPS). Section 2 describes the general design and the vacuum system of both pumping systems. The Test Particle Monte Carlo (TPMC) simulation with MolFlow+ [8] is described in Sec. 3. The large flow reduction along the beam tube forced us to subdivide the simulation of a long beam-tube into several independent steps and concatenate the results in the subsequent analysis. In the analysis of the CPS simulation the time dependence of the reduction factor was introduced, assuming a slow migration of the adsorbed and redesorbed tritium towards the spectrometer section. An alternative simulation method for the CPS is introduced in Sec. 4. The geometry is hard-coded in C++, optimizing the speed of the simulation. It also allows the simulation of timedependent properties that are difficult to characterize with MolFlow+ in a single pass. In Sec. 5, the results are presented and discussed in Sec. 6.

The Differential Pumping Section
The Differential Pumping Section (DPS2-F) has to fulfill three different tasks. The first task is to reduce the tritium gas flow between the WGTS and the CPS by 5 orders of magnitude. The second task is to guide the β-electrons adiabatically from the WGTS downstream to the Cryogenic Pumping Section (CPS) along the magnetic flux tube. Each tritium decay in the WGTS ionizes on average about 10 tritium molecules, which are also guided by the magnetic flux tube through the beamline. The third task is to prevent these ions from reaching the spectrometer section, where they would increase the background rate.

Geometry
The DPS2-F has five pump ports (PP0, PP1, ..., PP4), which are aligned perpendicular to the interconnecting beamlines (BT1-5). Pump port 0 was added at a late stage during the design of the pumping section, which is why in previous publications the DPS2-F is mostly described with only four pump ports [9,10]. Incoming and outgoing beamlines of PP1-4 change direction at the pump ports by 20 • (see Fig. 2). This geometry prevents a direct line of sight and increases the number of collisions with the walls, and thus, the pumping probability for the neutral molecules. Electrostatic dipoles (halfpipe-shaped stainless-steel electrodes) and ring-shaped electrodes remove tritium ions by either drifting them towards the walls where they are neutralized, or reflecting them back towards the source.
The beamline of the DPS-2F, between the WGTS and the CPS cryostat, is 7.3-m long. Each beam tube consists of a central tube, connected via bellows to a flange on each side. The central tubes have an inner diameter of 100 mm. The smallest diameter of the beamline is 85 mm, defined by the beamline instrumentation. The entire DPS2-F has a weight of about 10 4 kg and is fixed to the floor with earthquake protections. A computeraided design (CAD) drawing of the DPS2-F is shown in Fig. 2.

Vacuum system
The required gas flow reduction factor is achieved by six TMPs (Leybold MAG-W 2800). Two TMPs are connected to PP0. The other four TMPs are located at the lower ends of PP1-4. Each TMP can be separated Figure 2: Shown is the modular setup of the DPS2-F as a CAD drawing with a half-cut along the horizontal plane. Each beam tube (silver) is located in a warm bore of a superconducting solenoid (blue). The beam tubes are connected via the pump ports (green) to one beamline. The valves (orange) are located between the pump ports and the entrance of the TMPs (yellow). from the beamline by a DN 250 mm gate valve (VAT series 10). UHV vacuum gauges (MKS 421 inverted magnetron) are mounted on each pump port. The pressure along the DPS2-F beamline from PP0 to PP4 drops from approximately 10 −6 mbar to 10 −8 mbar. This absolute pressure reduction is not to be confused with the reduction of the tritium flow and the associated reduction of the partial pressure simulated here.
The superconducting magnets surrounding the beamline provide a guiding field of up to 5.5 T [11]. A passive magnetic shield encases each TMP to prevent eddy currents from heating up, and possible crashing, of the fast-moving rotors [12].

The Cryogenic Pumping Section
The last part of the transport and pumping section is the Cryogenic Pumping Section (CPS), which has to reduce the residual gas flow by more than seven orders of magnitude. For this purpose, a cold argon frost layer prepared on the inner beamline surface, with an area of about 2 m 2 that is maintained at 3 K, adsorbs the incoming tritium molecules.

Geometry
In Fig. 3 a CAD drawing of the CPS is shown. The 12-ton CPS cryostat built by ASG Superconductors S.p.A. is about 6.5 m long and 4 m high. The beam tube elements of the CPS are subdivided into seven sections with a total length of roughly 7 m (inner diameter: ≥75 mm) and two pump ports. Each section is surrounded by a superconducting solenoid that produces the 5.6 T magnetic field [11] to guide β-electrons adiabatically through the CPS. The second and fourth beam Figure 3: CAD drawing of the CPS cryostat in 3/4 section. The goldplated beam tube is surrounded by seven superconducting magnets (in red). The LHe vessel (4.5 K) provides a reservoir of 4.5 K cold helium, which is used for the cooling of the magnets and beam tube. The cold trap can be seen in blue between pump port 1 and the cold gate valve (CGV). The CGV is a safety valve operated inside the CPS at a temperature of 4.5 K. The differential pumping section (DPS2-F) is connected on the left side, the Pre-Spectrometer on the right side. tube are rotated by 15 • from the longitudinal spectrometer axis, so that neutral tritium molecules would hit the beam tube wall, where they are adsorbed. Each beam tube is a stainless steel tube with gold plated on the inner surface. Additionally, there are 90 circular fins (length of 2 mm) inside each of the beam tube sections 2-4 enlarging the inner surface [6].

Vacuum system
The main part of the vacuum system of the CPS is the 3-meter long cold trap (sections 2-5), covered by an argon frost layer at 3 K. The gold plated inner surface provides a clean surface for argon frost crystallization, as well as reducing the diffusion of hydrogen isotopes into the stainless steel. In order to reach 3 K in the LHe vessel (3 K) the 4.5 K helium is pumped down to a pressure of 0.16 bar and circulated through the cooling loop of the cold trap. For safety reasons the argon layer will be regenerated after 60 days of measurement time, which corresponds to an accumulated tritium activity below 1 Ci stored on the cold trap.
The other beam tube elements are cooled with liquid nitrogen and therefore operated at about 77 K. At PP1 and PP2 there are cold cathodes (MKS 421 inverted magnetron) to monitor the pressure in front of and behind the cold trap. TMPs are installed to both pump ports, but are turned off during standard KATRIN operation.

Cold trap temperature profile
The tritium molecules adsorbed on the frost layer can redesorb either via beta decay or by thermal desorption. Therefore the temperature of the argon frost layer is an important parameter for the pumping efficiency of the CPS. The correlation of the mean sojourn time τ des and the temperature T is [13] where τ 0 is the adsorbed particle's period of oscillation perpendicular to the surface, E des is the desorption energy for one mole of adsorbed gas and R = 8.314 J mol K −1 the molar gas constant. For monitoring the temperature three rhodium-iron sensors with 50 mK accuracy are attached to each beam tube section of the cold trap. During the first activation of the 3 K-cooling system, the measured temperatures on the beam tube did not reach the expected 3 K [14], but varied between 3.4 K and 6.2 K.
In order to investigate the origin of the temperature discrepancy, the heat transfer module of the commercial simulation program COMSOL Multiphysics R was used with a finite-element-method simulation. A CAD model of the cold trap geometry is imported; the model includes the magnetic coils 2-5, the inner radiation shields, part of the 3 K-cooling loop connected to the beam tube, and the beam tube. The simulation is initialized with a fixed temperature of 4.5 K for the magnets and 3 K for the pipes of the cooling loops around the beamline. To reduce the calculation time, radiation is only allowed between opposite surfaces, e.g. between magnetic coils and inner radiation shield, while the radiation inside the beam tube is neglected. The heat load from elements, which are not explicitly simulated in the geometry model, is taken into account by assuming a uniform thermal black body radiation with a specific temperature. In COMSOL Multiphysics R this parameter is called ambient temperature. In order to minimize the differences between the simulated and the measured temperatures, the ambient temperatures for the different beam tube sections vary between 70 and 90 K. In this way the non-negligible radiation of the pump ports is included.
In Fig. 4, the deviations to the measurement results are shown. The errors correspond to the temperature gradient along the 30.0 mm×14.8 mm×19.5 mm copper sensor housing connected to the beamline. The largest discrepancies are located in the regions near to the beam tube bends and in beam tube section 5. The first one can be explained by the higher radiative heat loads on these regions due to the gaps between the inner radiation shields and the magnets. In beam tube section 5, a 180 mm stretch of the cooling loop is not brazed to the beam tube due to a repair in this area, which leads to a large area on the right side that has a temperature above 5 K. Regions near the cooling loop reach the expected 3 K while the temperature increases further away. Due to the narrower windings of the cooling loop towards the end of the beam tube sections, there are areas which are homogeneous at 3 K. Hot spots arise as a result of bolts connecting the warmer inner radiation shield to the beam tube sections 2-4. Except for these hot spots, most of the beamline areas, in particular those with the fins, are in a temperature range between 3 and 4 K. This simulated temperature profile is used in the next sections to calculate the reduction factor of the cold trap.

The TPMC Models
The TPMC models have been simulated with the software package MolFlow+ (version 2.5.6) [15]. The software is designed for particle tracking in the molecular flow regime. Its basic concepts are described in the following section.

Simulation with MolFlow+
MolFlow+ tracks test particles through the model of a vacuum chamber build up by a polygon mesh, of socalled facets. The particles only interact with the walls.
When they hit a facet, they can either be adsorbed, reflected or transmitted. The properties of a facet are defined by various parameters, such as the temperature, the sticking factor defining the adsorption probability, the type of reflection, and the opacity. A facet can also be defined as a desorbing source of gas, where new particle tracks originate. Each track is simulated through a series of diffuse reflections and transmissions at facets until it is finally adsorbed. Each facet has three counters that are assigned for desorptions, hits and adsorptions; the counters are incremented accordingly when hit by a particle. A pump, such as a TMP, is represented by a circular facet the size of its entrance flange and a sticking factor α ∈ [0, 1] equaling its pumping probability for the simulated gas type. If a particle is reflected, we assume diffuse reflection, following the cosine law. Apart from the opaque facets representing the walls and internal structures of the vacuum chamber, the user can also define virtual facets, which are transparent. They do not affect the path of a particle, but count its transmission as a hit. A one-sided facet counts only particles impinging from one direction, while a two-sided one counts all particles. Virtual facets can be placed anywhere in the chamber, allowing us to monitor the pressure in the volume or to determine the transmission of particles from one part of the model to another.
After the simulation of an appropriate number of particle tracks, the results, represented by the three counters of each facet, are in general used to determine conductances, effective pumping speeds, and partial pressures. Details of the analysis methods for the DPS2-F and the CPS are described in the following sections.

The DPS2-F Model
For DPS2-F particle tracking simulations, the MolFlow+ model shown in Fig. 5 is used [7]. It comprises about 19000 facets. The models of the WGTS [4] and the DPS2-F join at the entrance of PP0. The TMP sticking factors α = 0.252 correspond to the pumping speed of TURBOVAC MAG W 2800 pumps for a gas mass of m = 6 g mol −1 , which has been chosen as the particle mass for the simulation of tritium molecules. It should be mentioned that the main systematic uncertainty of the simulations is the pumping probability α; the uncertainty is estimated using an empirical model for the pumping probability by Malyshev [10] to be 20%. The geometry of the electrodes of the ion detection and removal system inside the beam tube is included.
The particles desorb from the entrance facet of PP0 and are tracked through the complete geometry until they get pumped by one of the six TMPs or reach the entrance valve V2 towards the CPS. A sticking factor of α = 1.0 has been assigned to the V2 facet. The entrance facet of PP0 is transparent. On the upstream side the geometry is terminated by the last beam tube element of DPS1-F, which ends in a facet with sticking factor α = 1.0, representing the TMPs of the second stage of the DPS1-F.

Concatenation of subsequent sections
The simulation of a flow reduction factor R = Φ in /Φ out of 5 orders of magnitude and more is very time consuming. Therefore, the geometry of the DPS2-F was subdivided into four different parts (see Fig. 5). Each part was simulated independently, and the individual results were concatenated subsequently to derive the transmission probability, hits (for pressure) and adsorptions [7]. The reduction factor can be calculated by dividing the number of started particle tracks N des,E 0 at facet E 0 of PP0 by the number of particles N ads,V 2 adsorbed (leaving DPS2-F) at the gate valve V 2 between the DPS2-F and the CPS: Here N hit,G 2−4 are the number of particles passing the virtual facets in the downstream direction between the parts indicated in Fig. 5. Apart from the first part, the particle tracks were started at the entrances G 1−3 of the preceding beam tube section. Thus, the resulting solid-angle distribution of the velocities of the incoming tracks at facets G 2−4 was closer to the one expected for a single-pass simulation of the entire geometry. All simulations were done with the full model, changing only the desorbing facet where the tracks started. The simulations ran until the statistical uncertainty of the hits in the respective concatenating facets was better than 1%.

The CPS Model
The MolFlow+ model of the CPS is shown in Fig. 6. It consists of about 58000 facets. The geometry of the CPS starts at valve V2 and ends at valve V4. Seven elements build up the complete beam tube. The TMPs at the two pump ports are turned off during standard operation, so the only pumping mechanism is cryosorption. For particles reflected back into the DPS2-F the last beam tube element is included up to the last DPS2-F pump port (PP4). Particles reaching this pump port are considered to be pumped out with high probability (α = 1). This scheme ensures that the concatenation of the DPS2-F and CPS models is simulated with matching boundary conditions. A particle leaving the DPS2-F in the previous simulation through valve V2 is adsorbed (α = 1), not taking into account that a fraction of the particles is actually reflected back into the DPS2-F. This is done in the subsequent CPS simulation, making sure that the back-reflection is not taken into account twice.
In order to reduce the computing time, the CPS model is subdivided into four parts, which are simulated separately and concatenated afterwards, similar to the DPS2-F (see Sec. 3.2.1). Since the sticking factor of the argon frost layer is not precisely known and also depends on the initial coverage, simulations were performed for α = 0.0 to 0.7 in steps of 0.1. The upper value of 0.7 is the expected sticking factor for a well prepared argon layer at 3 K [16]. For facets not belonging to the cold trap, the sticking factor is set to α = 0.0. At both ends of the model (DPS2F-PP4 and V4) it is set to α = 1. In this case, the particles are either pumped out at DPS2F-PP4 or enter the Pre-Spectrometer volume, which has a diameter approximately twenty times larger than the CPS beam tube and a more than two orders of magnitude higher effective pumping speed for tritium than the conductance between the Pre-Spectrometer and the cold section of the CPS.
The tracking of gas particles starts at valve V2. Before a particle reaches the cold trap it is only diffusely reflected when hitting a wall. Once a particle is adsorbed on a facet of the cold trap it sticks on it forever and the tracking path in MolFlow+ ends. According to Eq. 1 in Sec. 2.2.3 this model is not appropriate since particles can leave the facet by thermal desorption after some characteristic sojourn time τ des . The redesorbed gas slowly migrates towards the end of the CPS, adding to the gas leaving through valve V4 towards the Pre-Spectrometer. Since the incoming gas flow from the DPS2-F is assumed to be constant, the reduction factor R CPS = Φ in /Φ out of the CPS decreases over time.
Therefore a model has been developed, which combines the results from a multitude of MolFlow+ simula- tions and, in a second step, adds the effect of a finite sojourn time τ des on the adsorbed tritium molecules. The basic idea is to subdivide the cold trap of the CPS beam tube into n = 102 smaller segments, and consider each segment as an individual cryo pump, where particles can adsorb and redesorb.

Segmentation of the cold trap
The number of adsorbed particles A i (t) sitting on the surface of segment i can either originate from the incoming gas flow Φ in through the DPS-2F or from the desorption off other segments j. The change in the number of particles adsorbed on segment i can be described by a system of coupled differential equations The first term describes the adsorption rate on segment i from Φ in . The parameter U ads,i des,V2 is the adsorption probability on segment i for particles entering the CPS through valve V2. The second term is the sum of adsorptions on segment i of particles desorbing from all segments j of the cold trap. The desorption rate A j (t)/τ des is proportional to the number of adsorbed particles A j (t) and the adsorption probability U ads,i des, j . The third term subtracts the rate of desorbing particles from segment i. All adsorption probabilities can be combined in a matrix U determined by MolFlow+ simulations.
The segments in the cylindrical part of the beam tubes were chosen to be 3.2 cm long, which corresponds to the distance between four consecutive fins (see Fig. 6). The cones at the ends of each beam tube element were simulated as longer segments. In addition, four simulations were necessary for simulating the gas inlet at valve V2, and the concatenations of the four parts, similar to the DPS2-F, for calculating the direct transmission through the CPS from V2 to V4. For the direct gas flow through V2, the results from these concatenating simulations were also used to calculate the adsorption probabilities U ads,i des,V2 for segments beyond H2 (see Fig. 6). For each segment i the counters of the corresponding facets were added up to the number of desorptions N des,i , adsorptions N ads,i and hits N hit,i . With these numbers the probability matrices for desorptions (U) and pressure (hit matrix V) can be calculated: N des, j is the probability that a particle desorbing from segment j is adsorbed on segment i.
N des, j is the probability that a particle desorbing from segment j hits segment i. In order to attain better statistics in the region, where most of the gas is adsorbed, the concatenations of the CPS geometry were also applied for calculating the hit V des, j ads,i and adsorption U des, j ads,i probabilities for desorptions from the segments of the first cold trap section ( This case-by-case analysis was constructed in such a way that the solid angle under which a particle enters the next part of the CPS is comparable to a single-pass simulation. For elements j ∈ [28, 102], the concatenation has not been applied since the first test simulations indicated a coverage reduced by several orders of magnitude compared to the first CPS simulation part.

Time-dependent gas flow
In order to describe time-dependent processes the coupled differential equations 3 of the amounts of adsorbed gas on each segment are numerically integrated over time with discrete steps of ∆t. The gas inlet into the CPS starts at t 0 = 0 s. The inlet rate Φ in from the DPS2-F stays constant over the whole time. In the simulation, this gas inlet is represented by the desorption of particles from facet V2. The time difference between t 0 and any other time t n is subdivided into n intervals with length ∆t. At the time t 0 = 0 s, it is assumed that there are no particles in the system at all (A i (0) = 0), which defines the boundary conditions for the numerical integration. After the first iteration at t 1 = ∆t the amount of gas adsorbed on a specific segment i, is described by For any time t n = n · ∆t > t 1 the number of adsorbed particles on segment i is defined as for i ∈ [1, 102]. The first term represents the number of adsorbed particles at the time t n−1 . The second term takes the gas inlet between t n−1 and t n into account. The factor A j (t n−1 )∆t/τ des in the third term equals the amount of desorbed gas from segment j. By multiplying this with U ads,i des, j , one gets the number of particles desorbed from segment j and adsorbed on segment i during the time interval ∆t. For the gas Φ out leaving the CPS towards the spectrometer section through valve V4 (i = 103), we neglect any gas coming back to the CPS due to the large pumping speed of the Pre-Spectrometer. With we can describe the time-dependent flow reduction as The pressure at segment i in the TPMC simulation can be determined from the number of hits N hit,i , the area of the segment F i , the mean thermal velocityc, the number of desorbed particles N des, j from segment j, and the actual gas flow (or outgassing rate) into the chamber Q j (t) [17]: Multiplying the particle flow Φ in with the Boltzmann constant k B and the temperature T , the pressure after the first iteration at t 1 = ∆t is and for t n = n · ∆t: The pressures at CPS-PP1 and CPS-PP2 are particularly relevant since these are measurable quantities. From simulations, the pressure ratio is This ratio is correlated with the flow reduction: By simulating both the flow reduction and the pressure ratio, the ad hoc factor k(t n ) can thus be determined. This result is essential for interpreting measured pressure values at the pump ports.

Including the beamline temperature profile
The COMSOL Multiphysics R simulation of the CPS cold trap temperature profile in Sec. 2.2.3 revealed inhomogeneities of several Kelvin. These inhomogeneities have a non-negligible influence on the pumping efficiency of the cold trap; in particular, the mean sojourn time τ des is affected. This can be included in the analysis of the MolFlow+ simulations by calculating the weighted mean of allτ des in the system: The weight A i is the area of one of the corresponding beamline surface elements with temperature T i of the COMSOL Multiphysics R mesh. With a fixed desorption energy E des , the flow reduction can now be calculated on an absolute time scale. Since the magnitude of the desorption energy E des is not known, a range from 1200 J mol −1 to 1600 J mol −1 is investigated. The lower boundary is taken from the estimation given in [18], the upper boundary can be estimated from the first retention measurements which will not be discussed in detail within this publication.

Tritium decay
For radioactive adsorbates, the sojourn time can no longer be described solely by the desorption time τ des given in Eq. 1. One has to take the influence of radioactive decays inside of the adsorbens into account. In addition to the released decay products, a tritium decay can also induce the desorption of other atoms in its vicinity, including both tritium and argon. The amount of desorbed tritium η(s) from a single β-decay inside the argon frost layer depends on the surface coverage s. It can be described with the following formula investigated by Malyshev [19]: where η max denotes the upper limit for the desorption yield, and s m the kink between the linear rise and the plateau where η(s) reaches saturation. The values of these parameters are estimated with η max = 10 3 T 2 /decay and s m ≈ 4 × 10 14 T 2 cm −2 [19]. Given this relation, the effective desorption time can be derived from the differential equation describing the rate of desorbing particles by both thermal desorption and radioactive desorption. The latter is described by the decay constant λ T . The variable σ is used to describe the gas composition and can take values from 0 to 2 depending on the amount of tritium atoms of the adsorbed isotopologues (H 2 : σ = 0, HT: σ = 1, T 2 : σ = 2, ...). The probability density distribution ρ(t) for the sticking time is then given by (16) with the effective time constant With the additional decay-induced desorptions, the effective sojourn time will no longer increase with falling temperatures but converges towards a limit as is shown in Fig. 7.

Semi-Analytical Tracking Model of the CPS
Since tritium decay plays a major role in decreasing the sojourn time of a tritiated molecule inside the cold pump, this effect needs to be taken into account when simulating the tritium reduction factor of the CPS. The currently available simulation programs for TPMC do not meet the requirements of simulating large reduction factors in combination with radioactive adsorbates and strongly inhomogeneous temperature profiles. This was the reason for developing a custom C++-based Semi-Analytical Tracking Model. By disassembling the geometry given in Fig. 6 into its basic geometric primitives (namely cylinders, cones and cuboids), the motion and interaction points with the surface can be calculated analytically. For the desorption process, a multistage Monte Carlo sampling is needed to determine both the sojourn time and the direction in which a particle desorbs. The former is done by sampling from the distribution given in Eq. 16, the latter with a cosine law sam-pling [20]. The advantage of the semi-analytical tracking is that the model offers the possibility to integrate three-dimensional models not only for the temperature, but also for the surface coverage along the beamline, which will be discussed in more detail in Sec. 5.2.2. The integration of a temperature profile will provide a more precise simulation result than just assuming an average beam tube temperature, since there are regions that are more likely to be hit because of the rotation of the beam tube sections. Therefore a very detailed temperature model from the simulations in Sec. 2.2.3 with more than 12 000 equally distributed temperatures along beam tube sections 2 to 5 of the CPS is used for the simulations. Compared to the TPMC simulations in Sec. 3.3, the Semi-Analytical Tracking Model requires a new simulation for each desorption energy by changing τ des and cannot be computed from one set of simulations.

Reduction factor calculation
With the Semi-Analytical Tracking Model, the migration time of a single particle along the entire CPS is calculated. The simulation of a particle track can have four different outcomes: i) The particle leaves the CPS through valve V4 into the Pre-Spectrometer. Only these events contribute to the determination of the reduction factor. ii) The particle is reflected back into the DPS2-F, where it reaches PP4 and is pumped out by the TMP. iii) The simulation is aborted because the particle takes longer to leave the CPS than the initially set maximal migration time. iv) The particle decays while still in the CPS, most likely being adsorbed on the argon frost layer.
Only 1% of the tritium decays during the nominal operation time of 60 days between subsequent regenerations of the argon frost layer, which is why the loss of tritium due to its decay is neglected in this simulation. The results can also be used for stable isotopes, such as hydrogen or deuterium, if the mean sojourn time τ des is much longer than the time of flight between two hits of the walls.
To extract the reduction factor of the cold trap, a significant amount of particle tracks and the corresponding migration times need to be calculated. Storing the information into a histogram, normalized to the total number of simulated particle tracks, gives a probability density distribution m(t) for the time t which a molecule needs to migrate along the geometry. Since m(t) depends on the desorption energy, tritium purity, sticking probability and the temperature, these simulations have to be repeated for each parameter setting. Once the probability density has been simulated, it can be transformed into a time-dependent tritium flux Φ out (t) from the downstream end of the CPS to the Pre-Spectrometer. This is done by integrating over the migration probability m(t) and multiplying with the expected constant tritium flux Φ in = 10 12 molecules s −1 [7] from the DPS2-F into the CPS The reduction factor R(t) after time t of continuous operation is defined as the ratio of the outgoing and incoming flux: The transmission probability density and its integral (Eq. 18) are shown in Fig. 8. Although the distributions were simulated with up to 2.5 × 10 11 events each, the simulations provide virtually no events in the region between 0 and 60 days for reduction factors R(60 d) > 10 12 . In this case a linear interpolation between the first non zero bin and the origin (where the tritium flow is expected to be zero) is applied.

Results
In the following, the simulation results of the flow reduction of both DPS2-F and CPS are presented. In Secs. 3 and 4, two different approaches were introduced. The DPS2-F results are solely based on MolFlow+ simulations, while the CPS results were complemented by the Semi-Analytical Tracking Model.

DPS2-F
The DPS2-F simulation has been split into four parts as described in Sec. 3.2. Table 1 gives an overview of the results. Concatenating all four simulations to an overall reduction factor yields R DPS2−F = (1.577 ± 0.008 stat. ) × 10 5 (20) with the statistical uncertainty calculated by using binomial statistics. To estimate the systematic uncertainty of the TMP pumping probability (see Sec. 3.2) a dedicated simulation with α = 0.202 was performed. The result R lower DPS2−F = (8.99 ± 0.05 stat. ) × 10 4 (21) gives a lower limit for the reduction factor with a 20% reduction of pumping probability.

CPS 5.2.1. MolFlow+
Two different scenarios were simulated with MolFlow+. The first one is the standard neutrino mass measurement; the parameter of interest is the reduction factor. The other scenario is the commissioning measurement with D 2 ; the parameters of interest are the pressures at both pump ports. In the commissioning simulation, the inlet valve V2 and the outlet valve V4 are assumed to be closed, while they are opened during the neutrino mass measurement simulation.
Combining the resulting values R(t) for the neutrino mass measurement simulation with p PP1 (t) and p PP2 (t) for the commissioning measurement simulation, the factor k(t) of Eq. 12 was determined.
The simulation results for the three relevant parameters are displayed in Figs. 9 to 11. Since the mean sojourn time τ des is unknown, the time axes are normalized to τ des . The time interval for the iterative integration was set to ∆t = 0.01 · τ des . As expected, the reduction factor and the pressure ratio show a similar time-dependent behavior. Over a period of 2τ des both parameters decrease by 2 to 3 orders of magnitude, except for α = 0.0 where they stay constant as there is no cryosorption at all. Lower sticking coefficients result in lower reduction factors and pressure ratios. The results for the ad hoc factor k(t) are important for interpreting the D 2 commissioning measurements. The simulated values lie between 8.5 and 21.5, and stay more or less constant.
For the nominal KATRIN operation sufficient tritium suppression for the whole 60-day run period  Figure 11: Factor k = R/(p PP1 /p PP2 ) for different sticking coefficients from α Ar = 0.0 to 0.7, simulated with MolFlow+. Here R is the reduction factor for standard tritium operation and p PP1 /p PP2 is the pressure ratio for the D 2 commissioning measurement scenario. For α Ar = 0.0, the value is constant at k ≈ 5.9. is of paramount importance.
To understand the long-term suppression, the x-axis in Fig. 9 has to be multiplied with a constant τ des . Therefore, the desorption energy E des is fixed, and the inhomogeneous beamline temperature profile is included according to Sec. 3.3.3. This has been done for α = 0.7 and three different desorption energies E des = 1200 J mol −1 ,1400 J mol −1 , and 1600 J mol −1 . The results are shown in Fig. 12.

Semi-Analytical Tracking Model
Because of the dependence of τ eff on η(s), a detailed knowledge of the surface coverage s along the segments of the CPS is required to simulate the impact of radioactive decays on the tritium reduction factor. With the cold trap temperature profile implemented in the simulation code, the surface coverage turns from a smooth distribution into the density map shown in Fig. 13. The cor-relation between the lower temperatures and a high surface coverage can be easily explained by the mean sojourn time of the molecules adsorbed on the argon frost, which depends strongly on the local temperature (see Eq. 16). The inhomogeneous temperature profile leads to an enhanced migration from regions of higher temperatures to regions of lower temperatures. As shown in Fig. 14, the mean surface coveragē decreases almost exponentially from about 10 15 T 2 cm −2 at the upstream entrance to 10 9 T 2 cm −2 or less at the downstream end of the CPS. Here s i denotes the surface coverage of bin i in Fig. 13 and n i the corresponding amount of molecules. Calculating the average weighted by particles within a given bin instead of its area is required to correctly simulate β-desorptions as described in Sec. 3.3.4. To cover the expected range of the desorption energy (see Sec. 3.3.3), a set of simulations reaching from E des = 1200 J mol −1 to E des = 1600 J mol −1 has been performed. For non-radioactive gases with desorption energies above 1200 J mol −1 , this simulation method is not suitable since the probability function used to calculate the reduction (see Eq. 19) has no entries close to the region from t = 0 days to t = 60 days. Therefore the uncertainty of any extrapolation would span several orders of magnitude.
But as soon as radioactive desorption is considered, the sojourn time and the reduction factor are drastically reduced. In this case, the Semi-Analytical Tracking Model is valid even for higher desorption energies and lower temperatures where τ eff converges towards a constant value, as shown in Fig. 7. This can be seen in Fig. 15 where the tritium reduction factor rises slower with larger desorption energies for T 2 than it does for isotopologues with only one tritium atom (HT, DT).
In order to reduce the complexity of the simulation code, a static mean surface coverages after 60 days of each beam tube section is used. These values were obtained from preceding simulations for each binding energy without taking radioactive desorptions into account (see Fig. 14). Usings is justified, since the influence of radioactive decays is only significant in regions with low temperature (see Fig. 7), where the majority of the molecules is adsorbed. Because no time dependency is implemented, the results for the reduction factor have to be seen as a conservative lower limit. Here the influence of the desorption energy on the reduction factor is only in the range of two orders of magnitude. Even for the very conservative assumption of a static surface coverage and a desorption energy of 1200 J mol −1 , the simulation yields a reduction factor of at least 2.6 × 10 10 , which exceeds the requirements by three orders of magnitude.

Discussion
With the specified reduction factor of 10 5 for the DPS2-F, the results obtained with the MolFlow+ simulations are right on target. The lower limit for the reduction factor takes into account the conservative uncertainty of the TMP pumping probabilities for T 2 .
On the other hand, the MolFlow+ results for the CPS show a reduction factor that exceeds the specified goal of 10 7 by several orders of magnitude. For the safe operation of KATRIN, it is important to understand the accuracy of these results.
The impact of the cold trap capacity on the sticking probability has not been included in the model. After 60 days of KATRIN operation, approximately 1.7% of the total cold trap capacity is reached [7]. Further, it is assumed that the first 1.7% of the trap is covered to 100% (α = 0), while the rest is free of adsorbates. Under the assumption, the reduction factor drops exponentially along the cold trap; the final value of 2.7×10 11 (see Fig. 12) for the lower limit of E des would be reduced by 36%.
A second source of systematic uncertainties is the concatenation algorithm of the four independently simulated parts of the CPS MolFlow+ geometry. By comparing simulations with a small sticking coefficient α = 0.1 for a single-pass simulation and with concatenation, the total error is estimated to be less than a factor of two.
The analysis of the time dependence with MolFlow+ simulations provides results for a time scale normalized to the sojourn time τ des . In order to extract results on an absolute time scale, an appropriate range for the unknown desorption energy E des has to be estimated. Reasonable values lie between 1200 J mol −1 and 1600 J mol −1 . In addition the temperature inhomogeneity of the cold trap was included by calculating an Figure 13: Two-dimensional representation of the tritium coverage after a pumping time of 60 days, simulated with the Semi-Analytical Tracking Model. The z axis follows the central axis of the beam tube. The azimuthal angle Θ covers the full circumference of the beam tube. The apparent inhomogeneities are due to the inhomogeneous temperature profile along the cryogenic pumping section.  Figure 15: Time-dependent tritium reduction when assuming no radioactivity for a desorption energy of 1200 J mol −1 (solid line) and for the case of a tritiated adsorbate for even higher desorption energies (dashed/dotted lines), simulated with the Semi-Analytical Tracking Model. The composition of the adsorbate is described with the variable σ which is 1 for HT/DT and 2 for pure T 2 . effective mean sojourn timeτ des (see Sec. 3.3.3). This is justified by the large number of adsorptions on different segments when particles migrate through the cold trap.
Considering all of the relevant systematic uncertainties, the results' accuracy is assumed to be on the orderof-magnitude level. The simulations are very important for characterization measurements of the CPS cold trap with D 2 because the reduction factor cannot be measured directly. The measurements might also help to reduce the large uncertainties of the input parameters. For the standard KATRIN operation, additional β-induced desorptions from tritium decays have to be taken into account, which would reduce the assumed sojourn time even further. Despite all these uncertainties in the MolFlow+ simulations, the expected reduction factor still exceeds the nominal value by several orders of magnitude.
Compared to the D 2 simulations with MolFlow+, the Semi-Analytical Tracking Model for tritiated isotopologues has a significantly lower reduction factor, as shown in Fig. 15. The reduction factor for tritium still exceeds the requirements by far. A necessary simplification had to be made in order to attain the reduction factor for radioactive adsorbates. This simplification includes the assumption of a time independent, mean surface coverage per beam tube section, as described in Sec. 5.2.2. This approach overestimates the surface coverage and results in a conservative limit.
Another systematic uncertainty arises from the linear interpolation of the cumulative probability density function of the migration time (see Fig. 8) for t < 60 days. If no event is produced in this region, a conservative lower limit of 2.5 × 10 11 for the reduction factor can be inferred.
In Fig. 16, the results of both simulation methods are compared for D 2 simulations for an assumed desorption energy of 1200 J mol −1 . The two very different approaches show similar results. After 60 days, the results R = 2.7 × 10 11 for MolFlow+ and R = 4.0 × 10 11 for the Semi-Analytical Tracking Model agree to within a factor of two. Taking into account the complexity and different approximations of both methods, the results can still be considered to be in good agreement.

Conclusions
In the KATRIN experiment, it is mandatory to reduce the tritium gas flow between the WGTS and the Pre-Spectrometer by at least 14 orders of magnitude. Tritium decays inside the spectrometer section would otherwise increase the background rate of the experiment, limiting the ultimate sensitivity for the neutrino mass. For this reason, the differential pumping section DPS2-F and the cryogenic pumping section CPS are located between the WGTS and the spectrometer section to reduce the tritium flow accordingly.
An initial simulation of the temperature profile of the CPS cold trap with COMSOL Multiphysics R revealed inhomogeneities that were taken into account in the subsequent vacuum simulations. For these simulations in the molecular flow regime, two different models were used: a TPMC simulation with MolFlow+ and a Semi-Analytical Tracking Model developed in C++.
The simulations help to infer actual flow reduction factors from the measured pressure ratios close to the inlet and outlet flanges of the pumping sections. The initial simulations of the cold trap of the CPS with deuterium included time-dependent migration along the beamline due to thermal re-desorption. The effect of βdecay-induced desorption was added by decreasing the desorption time τ des accordingly. Even with the large systematic uncertainties due to the wide spread of the possible input parameters in the simulations, both methods reach the same conclusion that the combined reduction factor of the pumping sections exceeds the design value by several orders of magnitude.
The preliminary results of the ongoing measurements with deuterium support the findings of these simulations. Further measurements and detailed comparisons with the models will help us to reduce the uncertainties of the input parameters, which would ultimately lead to more accurate predictions.

Acknowledgements
We acknowledge the support of the Helmholtz Association (HGF), the German Ministry for Education and Research BMBF (05A17VK2), the Helmholtz Alliance