PENTrack---a simulation tool for ultracold neutrons, protons, and electrons in complex electromagnetic fields and geometries

Modern precision experiments trapping low-energy particles require detailed simulations of particle trajectories and spin precession to determine systematic measurement limitations and apparatus deficiencies. We developed PENTrack, a tool that allows to simulate trajectories of ultracold neutrons and their decay products---protons and electrons---and the precession of their spins in complex geometries and electromagnetic fields. The interaction of ultracold neutrons with matter is implemented with the Fermi-potential formalism and diffuse scattering using Lambert and microroughness models. The results of several benchmark simulations agree with STARucn v1.2, uncovered several flaws in Geant4 v10.2.2, and agree with experimental data. Experiment geometry and electromagnetic fields can be imported from commercial computer-aided-design and finite-element software. All simulation parameters are defined in simple text files allowing quick changes. The simulation code is written in C++ and is freely available at github.com/wschreyer/PENTrack.git.


Motivation
Precision experiments with particles at low energies require an excellent understanding of particle trajectories. Apparatus effects can influence the measurements and lead to false results.
Measurements of the neutron lifetime are a prime example. Neutron-lifetime experiments storing ultracold neutrons (UCNs) in material bottles recently have suffered from poorly understood apparatus effects [1,2], e.g. unaccounted losses of UCNs at the bottle walls, and their results often deviate beyond the quoted uncertainties [3]. They also deviate from results of beam experiments, which determine the neutron-decay rate from cold-neutron beams [4].
To improve this situation, next-generation experiments like PENeLOPE [5], UCNτ [6], and HOPE [7] plan to trap UCNs in complicated magnetic-field configurations and plan to detect the decay products-protons and electrons. To study apparatus effects of this type of trap, simulation tools are needed to track neutrons, protons, and electrons in inhomogeneous, time-dependent electromagnetic fields.
The search for an electric dipole moment of the neutron (nEDM) using trapped UCNs may give key constraints to CP-violating mechanisms necessary to explain the matterantimatter asymmetry of the universe [8]. Several new nEDM experiments are currently under construction. To reach the aspired sensitivity of 10 −27 e cm, simulations are needed to study apparatus effects, e.g. geometric phases, and optimize the performance of these experiments [9][10][11].
Existing simulation codes, e.g. STARucn and MCUCN [10,12], allow to simulate interactions of UCNs with matter and spin precession in weak magnetic fields. However, they cannot calculate UCN trajectories in strong magnetic fields and require a description of experiment geometries based on combinations of basic volumetric shapes, making implementation of complex geometries difficult. Geant4 [13], based on [14], allows the most comprehensive simulations of UCNs, neutron spins, protons, and electrons in electromagnetic fields.
For the PENeLOPE project, we developed the simulation tool PENTrack. It allows simulations of complete neutron-lifetime and nEDM experiments. The implemented physics processes cover UCN transport, UCN storage in material bottles and magnetic traps, spin precession of neutrons and co-magnetometer atoms, and tracking of protons and electrons in electromagnetic fields. It provides a flexible configuration interface and allows to load complex electromagnetic fields and experiment geometries directly from finite-element (FEM) and computer-aided-design (CAD) software. In this paper, we describe the underlying physics and algorithms, compare our results to experiments and other simulation tools, and provide examples for the optimization of experiments and the estimation of false results due to apparatus effects.

Equation of motion
Simulations of particles in electromagnetic fields require a numerical integration of their equation of motion. PENTrack performs a error-controlled fifth-order Runge-Kutta integration [15] of the relativistic equation of motion, 1ẍ of a particle with mass m, charge q, magnetic moment µ, and relativistic Lorentz factor γ in the reference frame of the experiment. The force includes i) gravitational acceleration |g| = 9.806 65 m s −2 in negative z direction, ii) Lorentz force of a magnetic field B and an electric field E, and iii) the force of a magnetic gradient ∇|B| on the magnetic moment with a polarization p of ±1.

Interaction with matter
Ultracold neutrons strongly interact with matter; interactions of protons and electrons with matter are not implemented yet. The latter particles are considered lost as soon as they hit a surface.
All materials in the simulation are described by a complex optical potential U = V − iW [16]. Its real part, depends on the number densities n i and bound coherent scattering lengths b i of each nucleus species i. The imaginary part, depends on the loss cross sections σ l,i for a given velocity v n . This cross section is the sum of absorption and inelastic-scattering cross sections, since inelastic scattering increases the energy of a UCN so far above the storage potential that it can be considered lost. Scattering lengths and absorption cross sections are tabulated in [17]. Inelastic-scattering cross sections, however, are often unknown.
Interaction with the surface of a material can lead to reflection, absorption, or transmission through the surface.
The reflection probability R for a UCN with kinetic energy E hitting a surface at an incident angle θ depends 1ẋ = dx/dt andẍ = d 2 x/dt 2 . on the energy component perpendicular to the surface, E ⊥ = E cos 2 θ: If the UCN is not reflected but transmitted through the surface and into the material, the velocity of the UCN undergoes refraction, changing its kinetic energy to E = E − V . The wave number k = 2m (E − V + iW )/h becomes complex and leads to an exponential decay of the amplitude exp [−Im(k )x]. The loss probability after a path length d in the material is then A UCN impinging on a surface can be scattered specularly or diffusely. PENTrack calculates the scattering distribution of the latter process using either a simple Lambert model or the microroughness model, as introduced in [18] and validated in [19]. To calculate the total microroughness-scattering probability one has to integrate the distribution over all scattering angles. PENTrack uses a fast Gauss-Kronrod integration [20], which only slightly impacts the performance.

Spin motion
Every spin-1 2 particle has a magnetic moment of size µ parallel or antiparallel to its spin vector S. Its motion in the reference frame of the experiment follows the Bargmann-Michel-Telegdi (BMT) equation [21], which describes a precession around the sum of the magneticfield vector in the rest frame of the particle, and the Thomas-precession axis, After integrating a step of the particle trajectory, PEN-Track separately integrates the BMT equation along this step. The separate integration of trajectory and spin precession allows each numerical integrator to choose the optimal internal step length for each process and improves performance.
Since the BMT equation is only valid in small magnetic fields, the user can define a threshold field only below which the BMT equation is integrated. Once a particle enters a field above this threshold, its spin collapses into one of its fully polarized eigenstates. The probability P to find the polarization p being parallel or antiparallel to the magnetic field is given by the projection of the spin onto the magnetic field: The sign of the polarization can also be flipped during reflection on surfaces. This can be accounted for by assigning a fixed spin-flip probability to each material.

Configuration
PENTrack allows to load maps of magnetic and electric fields calculated with commercial FEM tools like OPERA [22] on regularly spaced grids. Both three-dimensional maps of arbitrary fields and two-dimensional maps of rotationally symmetric fields are supported. Field values between grid points are calculated with tri-and bicubic interpolation [20,23]. We also implemented analytic magnetic fields, e.g. nearly homogeneous fields with small gradients and fields of straight conductors. Arbitrary time dependence of the fields can be described by a user-defined function.
Experiment geometries can be imported from virtually all CAD software as StL files [24]. StL files approximate surfaces with triangle meshes. The intersection of a particle trajectory with such a surface is detected using the CGAL library [25]. Each part of the geometry can be deactivated in user-defined time intervals, making it possible to simulate variable properties of valves and other moving parts.
All simulation parameters-material properties, geometric model files, field maps, particle sources, and particle spectra-are stored in simple text files, allowing quick changes and optimization. Random samples from initial distributions and random processes are generated by a Mersenne Twister random-number generator, provided by the Boost libraries [26].
Several variables of each tracked particle can be recorded: its position, velocity, and polarization at beginning and end of the simulation, at user-defined times, and when hitting a surface; its complete trajectory; and the trajectory of its spin vector.
PENTrack is written in C++ and its object-oriented structure simplifies the implementation of new electromagnetic fields, particles, and physics processes. It is based on open-source libraries and freely available at github.com/ wschreyer/PENTrack.git.

Validation
We performed two benchmark simulations to compare PENTrack with Geant4 v10.2.2 and STARucn v1.2. In addition, we modeled two UCN experiments and were able to replicate their results using PENTrack.

Comparison with other simulation tools
The first benchmark simulates transport of UCN through a vertical, U-shaped UCN guide with a square cross section of 5 × 5 cm 2 shown in fig. 1. We selected stainess steel with an optical potential of (184 − 0.0184i) neV as material of the guide and used either microroughness reflection 2 or a 5-% probability of diffuse Lambertian reflection. Both ends of the guide are covered with perfect absorbers. The simulations uncovered several flaws in the microroughness reflection implemented in Geant4 v10.2.2, resulting in incorrect scattering distributions and probabilities. Once these flaws were corrected, the transmission of UCNs through the guide agreed very well among all three programs, with a slightly lower transmission in Geant4  (fig. 2).
The second benchmark simulates the trajectories of UCNs in a strong magnetic field and in matter. The geometry consists of a guide tube with a length of 6 m and a diameter of 85 mm, coated with diamond-like carbon ( fig. 3). A cylindrically shaped source generates UCNs at one end. Four meters downstream of the source, a superconducting polarizer magnet generates an inhomogeneous magnetic field penetrating the guide tube. We placed an aluminium foil with a thickness of 0.1 mm in the center of the field and a polyethylene absorber at the far end of the guide.
UCNs with one polarization state, so-called low-field seekers, are repelled by the strong magnetic field and cannot penetrate the magnetic barrier. They are mainly absorbed on the source side of the guide. UCNs with the other polarization state, called high-field seekers, are attracted to the strong magnetic field and accelerated towards the foil. They are mostly absorbed by the foil or the absorber at the far end.
The simulated results showed that Geant4 did not correctly account for refraction when a UCN entered a material, resulting in too little absorption of UCNs in the foil. Once this flaw was corrected, the results of Geant4 and PENTrack agreed very well ( fig. 4).
All corrections to Geant4 will be included in version 10.3.

Comparison with experiments
The first experiment we simulated measured transmission of UCN through thin foils of pure aluminium [27]. We imported their time-of-flight geometry and reproduced the quoted time-of-flight spectrum of UCN with an initially cosine-distributed angle between their velocity and the guide axis. We assigned a diffuse Lambert-scattering probability of 10 % and the quoted loss cross section of 229 b at 6.2 m s −1 to the aluminium foils, resulting in an optical potential of (54.1 − 0.002 81i) neV. The simulated transmission rate of UCN through foils with different thicknesses agrees very well with the experimental data ( fig. 5). The 2 We assumed a roughness amplitude of 2.6 nm and a correlation length of 20 nm [19]. mean free path of (0.748 ± 0.016) mm determined from an exponential fit to the simulated data matches the experimental value of (0.725 ± 0.009) mm.
In a second simulation, we imported the geometry of the UCN source at the Research Center for Nuclear Physics (RCNP), Osaka University, as described in [28] and compared the resulting storage time of UCN with experimental data. The source is filled with superfluid helium at 0.8 K and enclosed by walls coated with NiP. Since these walls are cold, we assumed that no inelastic scattering takes place and assigned an optical potential of (213 − 0.0224i) neV to the walls. Following [29], we assumed a UCN-loss rate in the superfluid helium of v n n He σ l,He = 0.002 75 s −1 , resulting in an optical potential of (18.5 − 9.05 · 10 −10 i) neV. Guides from the source to a UCN valve and a detector are made of stainless steel with an optical potential of (183 − 0.0852i) neV. At the beginning of the simulation, the source generates UCNs with an energy spectrum proportional to √ E between 0 and 350 neV. After a certain storage time, the valve opens and UCNs are extracted into the detector. The UCN lifetime in the source of (79.6 ± 1.2) s, determined from the measurements with storage times of 50 s or more, matches the experimental result of (80.9 ± 0.4) s very well ( fig. 6).

Comparison with analytical calculations
To validate the simulation of relativistic particles, we simulated beta-decay electrons with energies of up to 780 keV and isotropic velocity distribution in orthogonal homogeneous electric and magnetic fields with strengths of 10 kV cm −1 and 0.1 T. The resulting E × B drift can be described analytically: In an inertial frame moving with velocity u = E × B/B 2 with respect to the fixed laboratory frame, the electrons follow a helical path along a reduced magnetic field B 1 − u 2 /c 2 [30, ch. 12.3]. The difference between the analytical and simulated paths linearly increases with the path length and 95 % of the simulated particles deviate by less than one nanometer after a flight path of one meter (fig. 7).
To validate the simulation of spins, we simulated the spin of a neutron during typical nEDM-experiment cycles. nEDM measurements are based on Ramsey's method of separated oscillating fields [32]. A short, oscillating magnetic-field pulse rotates the spins of polarized UCNs by an angle of π/2. The spins are then left to precess freely in homogeneous magnetic and electric fields. After a certain time, another π/2 pulse-in phase with the first pulse-again rotates the spin. We simulated nEDM-experiment cycles using π/2 pulses with an amplitude of 10 nT and varying frequency. After 50 s of free precession in a homogeneous 1-µT magnetic field, the simulation generates the analytically calculated Ramsey pattern [31] (fig. 8). Table 1 summarizes the single-threaded processing speed of the three different tools during the benchmark simulations. For simple simulations of UCN transmission, STARucn offers the highest speed. In comparison, PEN-Track is slower by a factor of five due to its flexible but computationally intensive geometry description. Geant4 is three to ten times slower than PENTrack.

Performance
If material interactions are modeled with microroughness reflection, PENTrack is slowed down by a factor of three compared to simulations with purely Lambertian reflection. Geant4 uses look-up tables for the microroughness distributions, which do not reduce processing speed but require an initialization time of 460 s. This time is not included in the calculations for table 1.
To speed up simulations of large numbers of particles, PENTrack was specifically designed to be run in several parallel instances on multi-core processors and computing clusters. One can assign a job number to each instance via command line, which is prepended to the corresponding output files. All output files can be merged with dedicated tools provided with PENTrack.

PENeLOPE
Using PENTrack, we were able to simulate a complete measurement cycle of the PENeLOPE experiment. PENe-LOPE uses a large, superconducting magnet ( fig. 9) to trap UCNs and measure the neutron lifetime.
At the beginning of each measurement cycle, UCNs are filled into the storage volume for 200 s, achieving 95 % of saturation. A valve then closes the storage volume and UCNs are trapped by its walls made from stainless steel. For another 200 s, absorbers in the storage volume remove UCNs with energies high enough to overcome the minimum trapping potential during magnetic storage of 115 neV. After this cleaning stage, the superconducting magnet is ramped up within 100 s and low-field seekers become trapped by the magneto-gravitational potential During magnetic storage, a proton detector at the top of the storage volume can directly observe the decay rate of the trapped UCNs, from which we can determine their lifetime in the trap. After a predefined storage time, the magnet is ramped down again and the remaining UCNs are counted by a UCN detector, providing a second measurement of their lifetime in the trap. In our simulations, ultracold neutrons were created at the converter surface of a UCN source and transported through UCN guides to the experiment. The results allowed us to optimize the vertical position of the experiment, the guide geometry, and the filling time with respect to the number of stored UCNs.
Simulations of the cleaning stage allowed us to optimize the geometry of the UCN inlet to shorten the time required to remove UCNs with energies high enough to overcome the magneto-gravitational trapping potential ( fig. 10).
While the superconducting magnet is ramped up, the slowly increasing magnetic potential increases the total energy of low-field seekers, E + U m , which could push their total energy above the trapping potential ( fig. 11). At the same time, high-field seekers are accelerated towards the walls and undergo many wall collisions ( fig. 12). Both effects would introduce losses of UCNs and their lifetime in the trap would become shorter than the inherent neutron lifetime. From our simulations, we estimated the energy increase of low-field seekers during ramping. We also determined the loss rate of high-field seekers during magnetic storage and tested strategies to remove them from the storage volume.
Simulations of the trajectories of protons and electrons from decays of trapped low-field seekers during magnetic storage showed that a voltage of at least −25 kV should be applied to the proton detector to efficiently extract and detect the decay protons with energies below 0.75 keV ( fig.  13). The decay electrons can hardly be influenced due to their much higher energy of up to 782 keV.

Geometric phases in nEDM experiments
In nEDM experiments, due to the electric field, a hypothetical electric dipole moment would slightly shift the precession frequency of the stored neutrons' spins. If the electric field is inverted, this shift is also inverted and a small phase difference between spins in opposite electric fields is accumulated over the free-precession time. A major uncertainty in nEDM experiments is caused by geometric phases, which can mimic the effect of an electric dipole moment. As shown in [9], such geometric phases can arise due to small gradients in the magnetic field combined with the relativisticẋ×E term in the BMT equation (8). To make this tiny effect visible, we simulated pairs of UCNs with identical trajectories subjected to electric fields with opposite direction. The UCNs were stored in a cylindrical nEDM chamber with a radius of 20 cm and a height of 10 cm. It was placed in a magnetic field with a strength of 1 µT and a rotationally symmetric vertical gradient of 10 pT cm −1 , and an electric field of ±10 kV cm −1 . Since the magnetic field has to obey the Maxwell equations, a vertical gradient leads to additional radial field components [9]. As shown in fig. 14, the spins of the stored UCNs accumulate a net phase difference, mimicking an electric dipole moment depending on the average velocity-although the UCNs have a random spatial distribution, have an isotropic velocity distribution, and undergo diffuse reflection on the chamber walls. This nicely replicates the calculations and simulations performed by [9].

Conclusions
PENTrack is a tool that allows comprehensive simulations of neutron-lifetime and nEDM experiments-including UCN transport, UCN storage in material bottles and magnetic traps, spin precession of neutrons and co-magnetometer atoms, and tracking of protons and electrons in electromagnetic fields. It provides a flexible configuration interface and allows to load complex electromagnetic fields and geometries from FEM and CAD software.
Detailed comparisons of results obtained with PEN-Track, STARucn v1.2, and Geant4 v10.2.2 showed very good agreement with STARucn and uncovered several flaws in Geant4. STARucn offers higher speeds than PEN-Track, but lacks support for microroughness reflection and magnetic fields. The speed of PENTrack is limited by its flexible geometry import, which, however, makes PEN-Track much more suitable for implementing complicated experiment geometries. The very general particle-simulation framework Geant4 offers similar functionality but has limited performance.