Ultracold Neutron Storage Simulation Using the Kassiopeia Software Package

The Kassiopeia software package was originally developed to simulate electromagnetic fields and charged particle trajectories for neutrino mass measurement experiments. Recent additions to Kassiopeia also allow it to simulate neutral particle trajectories in magnetic fields based on their magnetic moments. Two different methods were implemented: an exact method that can work for arbitrary fields and an adiabatic method that is limited to slowly-varying fields but is much faster for large precession frequencies. Additional interactions to simulate reflection of ultracold neutrons from material walls and to allow spin-flip pulses were also added. These tools were used to simulate neutron precession in the Paul Scherrer Institute's neutron electric dipole moment experiment and predict the values of the longitudinal and transverse relaxation times as well as the trapping lifetime. All three parameters are found to closely match the experimentally determined values when simulated with both the exact and adiabatic methods, confirming that Kassiopeia is able to accurately simulate neutral particles. This opens the door for future uses of Kassiopeia to prototype the next generation of atomic traps and ultracold neutron experiments.


Introduction
At sufficiently low energies, neutrons can be reflected from material walls by the coherent strong interaction [1][2][3]. Such neutrons can be stored in material traps, and are termed ultracold neutrons (UCNs).
The ability of UCNs to be stored for prolonged periods makes them of interest in several experimental areas. They can be used to measure the neutron lifetime [4,5], which continues to be of interest due to the remaining discrepancy between different experimental measurements [6]. UCNs can also be used to constrain the neutron electric dipole moment, which is related to the strong CP problem and can constrain beyond the standard model physics [7][8][9]. Currently, the best upper bound on the neutron electric dipole moment comes from measurements on UCNs, described in [10].
Several software packages have previously been used to simulate UCNs, including Geant4 [11,12], PENTrack [13], STARucn [14,15] and MCUCN [16]. The topic of this paper, the Kassiopeia toolkit [17], provides particular set of advantages and additional options for the users. In this paper, we demonstrate the ability of Kassiopeia, a software package originally developed for neutrino experiments, to accurately simulate UCN storage. In section 2, we describe some general features of Kassiopeia, not specific to UCNs. We then describe neutral particle tracking features that have been added to Kassiopeia in order to simulate UCNs in section 3. Next, we describe a specific experiment that we simulated and the parameters of our simulation in section 4, and describe our results in section 5. In section 5, we have also highlighted the advantages of our toolkit and compared the processes used in vetting various other toolkits. Finally, in section 6 we consider possible future applications of Kassiopeia's neutral particle tracking capabilities.

Spin features in Kassiopeia
At each integration step of a trajectory calculation, two parameters of the particle's evolution must be calculated: the time derivative of the spin, and the force on the particle as a function of the spin and surrounding fields. We therefore next describe how these are calculated for the exact spin trajectory.
The relativistic generalization of a classical spin is the four-vector S = (S 0 , S) given by, for a particle with classical spin s and velocity v = βc and Lorentz factor Γ (to avoid confusion with the gyromagnetic ratio γ) [25], S 0 = Γβ · s, ( 1 a ) The equations of motion for a relativistic spin with four-velocity U α in fields F μν , using the particle's proper time τ , are given by the BMT equation: The spin-dependent contribution to the particle's motion is given by the magnetic dipole force term, which is implemented nonrelativistically: This exact spin tracking method is completely general, except for the assumption that relativistic corrections to the magnetic dipole force term are negligible. However, the step size for numerically integrating these equations of motion is limited by the precession rate of the spins, γ|B|. This is typically fast relative to other timescales of a particle's motion, which can make the exact spin tracking method impractically slow for experiments with long timescales. Given the high computational cost of simulating tracks with the exact spin trajectory, an adiabatic spin trajectory was also implemented. The adiabatic spin trajectory records particle spins using two variables instead of four: an aligned spin m that gives the component of the spin along the magnetic field at the particle's position, and a spin angle φ that gives the orientation of the spin around that field. This spin angle is defined with respect to two unit vectors orthogonal to the local magnetic field, e 1 and e 2 , defined, for a magnetic field B = (B x , B y , B z ) in laboratory coordinates, as: Given these axes, φ gives the angle of the projection of the spin into the plane perpendicular to B away from e 1 , in the direction of e 2 : The relationship between these vectors is illustrated in figure 1. These two parameters, m and φ, along with the position x, then give the adiabatic spin equations of motion (excluding contributions independent of spin) [26]: where x is position, γ is the particle's gyromagnetic ratio, M is its mass, s is the magnitude of its total spin, and we have defined: for convenience. This adiabatic spin trajectory is accurate when particles are non-relativistic and spin precession is fast compared to the rate at which the field at the particle's location changes. Assuming the spin precession rate is dominated by the first term γ|B|, this requires γ|B| |ẋ · ∇B|/|B|, which should hold for all applications discussed in this work. During a single integration time step, the exact trajectory requires the spin to rotate by much less than a radian, so the left-hand side of that inequality approximately sets the maximum step size for exact spin integration. The adiabatic trajectory can handle large rotations per time step but assumes a constant field during the step, so the right-hand side limits the step size for adiabatic spin integration. In this limit, therefore, the adiabatic spin trajectory can be much faster. These were run in a version of the neutron precession chamber described later in this work, but with the magnetic field scaled up by a factor of 1000, resulting in a precession frequency of approximately 30 MHz. Note that the plot necessarily suffers from aliasing, since the adiabatic time step was chosen to be larger than the precession period. However, sampling was done simultaneously so that points at the same time can be compared.
A comparison of the outputs from the adiabatic and exact trajectories is shown in figure 2. Note that the exact trajectory becomes increasingly inaccurate as the integration time approaches the precession period 33 μs, even for 10 ms tracks. There may, however, be other limitations on the step size that prevent this from being the case. In particular, for the UCN tracking described in this work, weak magnetic fields mean that the limiting factor on the time step is error in position, not spin, so adiabatic and exact spin tracking require similar computation time.
In this article, we present our results using both methods. The exact method could be especially applicable in non-adiabatic cases, e.g. tracking the spin-precession of co-magnetometer 199 Hg atoms in room temperature neutron EDM experiments [27], where typical magnetic fields of order ∼ 1 μT are used. In this case, γ Hg |B|τ c < 1, since the 199 Hg atoms are relatively hot, where τ c = 2R/v RMS , R is the transverse size of the storage chamber, v RMS is the root mean square velocity of the species under consideration ( 199 Hg atoms in this case), and γ is the gyromagnetic ratio of the atoms. On the other hand, in order to track the UCNs in a similar magnetic field environment, one can use the adiabatic method, since γ n |B|τ c > 1, owing to the low velocity of the UCNs. This ability to choose the tracking method, allows greater flexibility in the applicability of this toolkit.

UCN storage
Having established that the exact and adiabatic spin tracking methods agree under appropriate conditions, we next consider how well Kassiopeia can reproduce experimental data from UCN precession chamber in references [10,[27][28][29]. The neutron precession chamber under consideration is a cylinder of radius 23.5 cm and height 12 cm. The interior rounded surface is primarily deuterated polystyrene (dPS), with small windows of deuterated polyethylene, while the two flat ends of the cylinder (the electrodes) are covered in diamond-like carbon (DLC).
UCN reflections from surfaces are described by three features: the reflection probability, the depolarization probability, and the distribution of outgoing angles for reflected neutrons. The first two are described in [30], while the third is described in [31]. We will summarize these and discuss our estimates of the associated parameters.
The UCN reflection probability at an energy E and angle to the surface normal θ is given by [32]: where E ⊥ = E cos 2 θ is the component of energy from normal momentum, V f is the real part of the surface's optical potential, and η is another constant of the surface. For dPS, these values are known only approximately, with η ∼ (1-3) × 10 −4 [30] and V f ∼ 161 neV [33], while a range of values for DLC is η ∼ (0.7-1.7) × 10 −4 [34] and V f ∼ 286 neV [35]. The UCN depolarization probability is the simplest to describe, and is given by a constant β. For dPS, this is known very roughly as β ∼ (1.3-8.9) × 10 −6 [30], and again a range of values for DLC is β ∼ (0.7-2.2) × 10 −6 [34]. Detailed theoretical treatments of the distribution of outgoing angles for reflected UCNs can be found in [31,32,36]. Since the behavior of UCNs in the precession chamber should not depend significantly on the exact angular distribution, we derive an approximate form. The angular distribution of stored neutrons is averaged out and does not play an effect in rotationally symmetric storage chambers such as the cylinder considered in this work [16].
The roughness of surfaces leads to non-specular reflection. In our simulation, which involves storage of UCNs in a cylindrical chamber, the role of roughness can be played by adjusting the loss (η) and depolarization (β) parameters [32]. However, non-specular reflection can be much more significant in other contexts, especially in neutron waveguides. We have implemented the micro-roughness model in our toolkit to allow for such simulations.
We start from the micro-roughness model result which, up to a normalization constant, can be modeled as: Here, θ i and θ f are the incident and outgoing angles to the normal, φ f is the change in direction around the normal, w is the surface roughness, k is the neutron wave vector, and, where k c = 2mV f / is the wave number corresponding to the optical potential and we consider only neutrons with k < k c , since higher-energy neutrons quickly escape the precession chamber. Then, where we defined Δθ = θ f − θ i . This gives a simple approximate distribution of outgoing angles that can be efficiently sampled. Note that this derivation assumes that the change in incident and outgoing angles is small. This is a rough approximation in our case, but we found that our simulations were insensitive to the exact parameters of the reflected distribution, so this should not significantly affect our results. For our simulations, we used β = 6 × 10 −6 , η = 1 × 10 −4 , V f = 150 neV, and w = 30 nm, which were consistent with the available estimated ranges of each value (see [37,38] for measurements of w). Note that the first two parameters were chosen so as to be most consistent with data in references [39][40][41][42]; we emphasize that we did not attempt simulations with multiple values of these parameters in order to find those most consistent with the measured depolarization timescales and storage times. In the future, it might be possible to use our simulations and the known timescales as a means of actually determining these parameters, but this was beyond the scope of this work. Extraction of these input parameters from the published values of T 1 , T 2 , and storage time constant in references [39][40][41][42], respectively, will be a topic of a future article.
For each test simulation, the precession chamber's interior had an essentially constant magnetic field of approximately 1 μT along the cylinder radius, with O(nT) deviations along each coordinate axis. In an experiment dedicated to measuring the electric dipole moment, the apparatus would also contain an electric field in either axial direction, but this was not included in our simulation since electric effects are below the sensitivity of the experiment [28]. The axial magnetic fields were different (and in particular were in opposite directions) for the two electric field directions, and we used a set of magnetic field measurements of each field configuration [44] along with a combination of natural neighbors interpolation (to reach a rectilinear lattice) and cubic interpolation (within lattice cells) to estimate the magnetic field everywhere.
During a particular run of our simulation, a number of neutrons would be added to just above the bottom of the precession chamber. While neutrons may not actually all start near the bottom surface, their motion around the trap is sufficiently fast relative to the simulation times that the choice of the origin location of the neutrons should not affect the results significantly. The initial direction of the neutrons were picked isotropically. For the sake of this simulation, neutron kinetic energies were generated from a Gaussian distribution with a mean of 150 neV and a standard deviation of 50 neV. This is an energy distribution that is harder than the energy distribution estimated in reference [45]. Softening of the energy spectra during storage is chiefly due to the loss of higher energy neutrons. Our simulation includes this loss channel, so the choice of initial energy spectra has a negligible effect on the relative parameters we use to simulate in this paper. It may be possible to match the measured energy distribution more closely by varying the parameters of β, η, V f and w, described above, but this was beyond the scope of this work. During a given run, all neutrons would also begin with the same spin, which varied by experiment as described below. Three example tracks inside of the simulated precession chamber are shown in figure 3.

Results and discussion
We chose to vet the updated Kassiopeia toolkit by simulating UCN storage, since there was experimental data available to benchmark its performance and accuracy. We compared our simulation's results with three previously measured parameters associated with the precession chamber under consideration: the longitudinal relaxation time T 1 [39,42], the transverse relaxation time T 2 [40,42], and the effective trap lifetime τ including neutron decay and wall losses [41,42]. Using a similar simulation but with the converse of the process we have employed here, the parameters of η, V f and w associated with the trap could be obtained just by using the data of neutron storage loss as a function of storage time alone, as was done in reference [20,43]. Previous works neglected the effects of depolarization. Particularly, the improvement we present in this work also takes into consideration the effects of depolarization per bounce parameter of β, allowing us to simulate the longitudinal and transverse depolarization rates in addition to the trapping lifetime.
A total of 17 280 neutrons were simulated up to maximum times of 500 s. This was done using four parallel c4.8× large instances of amazon elastic compute cloud (EC2), or 144 total cores. Kassiopeia does not directly support parallelization of particle tracking, so we instead ran 144 separate simulations of 120 neutrons each and then combined the outputs. This took approximately 2500 CPU-hours, though further optimization is likely possible.

The longitudinal relaxation time T 1
The longitudinal relaxation time, T 1 , is the time constant of the decay of average polarization for neutrons with initial spins that are either aligned or anti-aligned relative to the local magnetic field. It should be noted that, even for UCNs, the μT-order magnetic fields in the precession chamber do not significantly affect neutron energies since μ N |B| is far less than the kinetic energy. On a microscopic level, longitudinal relaxation is associated with inhomogeneities of the magnetic field, and with depolarization due to reflections [27].
Since the adiabatic equation of motion (6b) for fully aligned or anti-aligned spins always gives a zero derivative of the aligned spin, our simulations instead started spins at 1 • or 179 • to the magnetic field, though we found that any angle up to a few degrees led to indistinguishable results, as the dominant contribution to T 1 came from reflection depolarization.
We simulated 4320 neutrons for 500 s using both the adiabatic and exact trajectories with half initially aligned and half initially anti-aligned, and compared our average polarization as a function of time to that obtained in references [39,42]. Our results can be seen in figure 4, where a relative polarization of 1  [39,42] with the results of our adiabatic and exact simulations, including 1σ confidence intervals for the simulations and 1σ error bars for the experimental data. Here, a 0.95 relative polarization corresponds to an average spin along the z-axis equal to 0.95 times the initial average spin along the z-axis. Each data set averages over the two field configurations and over initially aligned and anti-aligned neutrons. We see good agreement over the entire period; the deviation of the exact data at long times is not statistically significant. Table 1. The results of χ 2 testing of our results against the experimental data in references [39][40][41][42], including left-sided p-values. Since remaining atom counts and polarizations are highly correlated (an upward deviation at one time will tend to lead to upward deviations at later times), ratios of counts and polarizations at successive times were compared instead. The generally small χ 2 values suggest that we are overestimating our errors, though it is unclear how this could be the case. Conversely, the T 1 adiabatic lifetime data's disagreement with the experimental data may indicate that our assumed value of η was not correct. indicates an average polarization equal to the initial value (±1 in simulations, but smaller magnitudes in references [39,42]). As figure 4 shows, both the adiabatic and exact simulations were able to accurately reproduce the measured longitudinal depolarization rate. Both simulations are within 1σ of the experimental data (using the combined error) until around 350 s, and the adiabatic simulations continue to closely match the experiment for the entire run. The apparent discrepancy between the experimental data and the exact trajectory at high times is not statistically significant (see table 1).

The transverse relaxation time T 2
The transverse relaxation time, T 2 , is the time constant of the decay of average polarization for neutrons with initial spins that are aligned along an axis perpendicular to the local magnetic field. The dominant contribution to transverse relaxation of neutrons in the precession chamber is the variation in the magnetic field over the trap volume [27]. Note that spin-spin interactions, which may be the dominant contribution to T 2 in other contexts, are not significant for neutrons in this context owing to their low density [46]. However, spin-spin interactions between two neutrons cannot be implemented in Kassiopeia at this time, as it tracks particles one at a time.
In the T2 measurements reported in references [40,42], neutrons began with polarizations aligned or anti-aligned with the magnetic field and then T 2 was measured between two π/2 pulses that rotated them into and out of the perpendicular plane defined by the direction of the main ∼ 1 μT holding magnetic field. While ideal spin rotating pulses are implemented in Kassiopeia, they were not used for this simulation since we were interested only in the behavior occurring within the field-perpendicular plane.  [40,42] with the results of our adiabatic and exact simulations, including 1σ confidence intervals for the simulations and 1σ error bars for the experimental data. Here, a 0.95 relative polarization corresponds to an average spin along the axis of average precession equal to 0.95 times the initial average spin along the initial polarization axis. Each data set averages over the two field configurations and over two initial spin directions. We see good agreement over the entire period; the deviation of the exact data at long times is not statistically significant.
We simulated 4320 neutrons for 500 s using both trajectory types and two different initial spin directions perpendicular to the magnetic field. We then compared our average polarization as a function of time to that obtained in references [40,42]. Our results can be seen in figure 5, with relative polarization defined as for T 1 , but along the axis corresponding to average precession over that time instead of along the z-axis.
As figure 5 shows, both simulations were also able to accurately reproduce the measured transverse depolarization rate. The adiabatic simulation matches all of the experimental values, however the exact simulation results appear to deviate from the data at long times. However, this discrepancy is again not statistically significant (see table 1), as our simulation errors are highly correlated: neutrons that deviate significantly from the average polarization (due to depolarization on reflection, for example) are likely to remain in the trap for some time, lowering the average polarization at those later times as well.

The effective trap lifetime τ
Since neutron reflection from material surfaces is not guaranteed, especially at higher energies (see equation (8)), the effective lifetime of neutrons in the precession chamber is significantly lower than their decay lifetime. Neutrons may also be lost due to holes or other imperfections in a physical trap. Storage lifetimes may be extracted independently from both the simulations dedicated to T 1 or T 2 , but the storage lifetime should be nearly independent of the neutrons' spins in a ∼ 1 μT field. The results from both sets of simulations, as well as from references [41,42], are shown in figure 6. Figure 6 once again shows good agreement between each of the four simulations and the experimental data. The adiabatic T simulation deviates significantly from the experimental data (see table 1), but the lifetime is strongly dependent on the chosen value of η and the initial energy distribution-both of which are known with limited precision for the chamber under consideration-so this discrepancy is not practically significant.

Summary of results
A summary of our simulated data's agreement to the experimental data in references [39][40][41][42] is presented in table 1, showing generally good agreement. The apparent overestimation of errors is expected to be the result of correlations between our errors at different times. This demonstrates that Kassiopeia can reproduce the effective lifetime as well as both decoherence times for neutrons within a precession chamber, confirming that it can accurately simulate neutral particles.

Comparison of the simulation toolkits
As other authors have noted before us [16], it is critical to have multiple simulation toolkits to allow for comparisons and thereby build confidence in the physics results. Kassiopeia toolkit offers a versatile new resource to augment the various existing simulation toolkits. In this sub-section, we compare the features and validation strategies of other toolkits while highlighting the advantages of the updated Kassiopeia toolkit.  [41,42] with the results of our adiabatic and exact simulations, including 1σ confidence intervals for the simulations and 1σ error bars for the experimental data. Averages for each data set are taken as in figures 4 and 5. Good agreement is seen along the entire data set from references [41,42]. Note that the exact data uses a shifted y-axis for visual clarity. The experimental values are arbitrarily normalized as the initial neutron count is not known.
All of the other toolkits, i.e. Geant4 [11,12], PENTrack [13], STARucn [14,15] and MCUCN [16], are limited to slowly varying fields to accommodate the adiabatic case. Additionally, all but the PENTrack toolkit are restricted to weak fields, and are unable to simultaneously simulate co-magnetometer atoms. Kassiopeia allows for non-adiabatic systems to be simulated. While PENTrack does allow certain co-magnetometer atoms to be simulated, our new Kassiopeia package is built to handle neutral particles, atoms, and molecules. Given that most co-magnetometer species are in fact non-adiabatic systems, Kassiopeia is the only toolkit currently available that is able to realistically simulate co-magnetometer atoms. Furthermore, while MCUCN and PENTrack have implemented reflection over rough surfaces, it is left up to the user to extraneously add roughness to materials in Geant4UCN and StarUCN, in order to simulate non-specular reflection. Also, MCUCN and StarUCN require magnetic fields defined as a combination of basic volumetric shapes. In Kassiopeia (and PENTrack), generic finite element field maps can be directly used.
Previous neutron simulation toolkits have taken a variety of approaches to verifying their simulation results. In order to validate Geant4UCN, reference [11] compares their simulation results-collision frequencies and transmission in a neutron guide-with prior numerical code. Similarly, StarUCN [14] has compared various scenarios with MCUCN and Geant4UCN. PENTrack [13] has presented comparisons of transmission in neutron guide and neutron absorption in foils and absorbers with Geant4UCN and StarUCN. MCUCN [16], on the other hand, presents comparisons with analytic solutions for the evolution of the energy spectra and center-of-mass of the ensemble of stored UCNs.
No previous toolkits has presented checks involving spin precession tracking, as we have done here by vetting the depolarization time constants using the Kassiopeia toolkit. In this article, we have presented comparisons between Kassiopeia results and real world data, demonstrating the capabilities of this new toolkit.

Other applications
Having established that Kassiopeia can accurately simulate neutral particles' motion in at least one context, we next consider a few potential future applications of these features.
The application most directly stemming from the simulations in this paper is simulated prototyping of future UCN experiments [47][48][49][50][51]. In this work, we simulated neutron precession using a previously measured magnetic field, but Kassiopeia also includes extensive field calculation methods, which can be used to calculate the magnetic fields inside of a trap based on its design [18,19]. Kassiopeia could therefore be used to estimate the effective lifetime and relaxation times of future UCN experiments.
Such simulation toolkits are vital in the analysis of many experimental results, and provide inputs that are an integral part of the final result. For example, searches for neutron to mirror-neutron oscillations depend on the mean time-of-flight of the stored neutrons inside the chamber. The mean time-of-flight, in turn, depends intimately on the time evolution of the energy spectra of the stored neutrons. A simple model involving a single loss-coefficient for the entire chamber was simulated using MCUCN to extract the energy evolution and thereby the mean time-of-flight of the stored neutrons in references [20,43]. Kassiopeia can further improve such simulations by performing realistic simulations of co-magnetometer atoms simultaneously with UCNs.
Neutron simulations have, in fact, already had an impact on measurements of the neutron lifetime. Several Monte Carlo (MC) simulations have been developed to model the early UCN experiments that used material storage bottles [21,22]. These MC simulations have adjusted the corresponding lifetime by a few seconds, reducing the tension between the beam and the bottle lifetime measurements [23,24]. There are no independent simulations verifying the claims of these MCs, however. The Kassiopeia toolkit may be able to provide a new, comprehensive check of these results.
Another potential application is simulating tritium storage methods in the Project 8 neutrino mass experiment [52]. Currently, atomic tritium is stored using a purely magnetic trap. We have previously simulated versions of this trap in order to obtain estimates of necessary currents and effective lifetimes, but development of the final trap configuration is ongoing and further simulation work could be done in this area. Furthermore, one current proposal is to use a combined magneto-gravitational trap instead, which has not been simulated as of writing.
Neutral particle simulation can also be applied to atomic experiments. Experiments such as ACME [53] and others such as [54] use molecular beams to measure the electric dipole moment. Molecules are prepared in particular spin states and then precess through constant electric and magnetic fields; a dependence of precession on the electric field can indicate an electric dipole moment. Kassiopeia may be able to simulate molecular evolution within these fields to aid with estimating systematic errors.