A New Design for a Traveling-Wave Zeeman Decelerator: I. Theory

The concept of a novel traveling wave Zeeman deccelerator based on a double-helix wire geometry capable of decelerating paramagnetic species with high efficiency is presented. A moving magnetic trap is created by running time-dependent currents through the decelerator coils. Paramagnetic species in low-field-seeking Zeeman states are confined in the moving traps and transported to the end of the decelerator with programmable velocities. Here, we present the theoretical foundations underlying the working principle of the traveling-trap decelerator. Using trajectory simulations, we characterise the performance of the new device and explore the conditions for phase-space stability of the transported molecules.


Introduction
In recent years, significant efforts have been invested into developing methods for the production of molecules at cold (<1 K) and ultra-cold translational temperatures(<1 mK) [1,2,3]. Besides the development of direct laser cooling [4], different techniques have emerged for cooling molecules based on their interactions with electric or magnetic fields [5]. These developments have been motivated by prospects of studying molecular collisions and chemical reactions at low and precisely controllable collision energies [6,2,7,8], of precision spectroscopic measurements for testing fundamental physical concepts [9,10], and of new approaches to quantuminformation processing and quantum simulation [11,12,13]. Methods based on the deceleration of supersonic molecular beams are particularly well suited for collision experiments since the final longitudinal velocity of the sample can be tuned over a wide range with narrow velocity spreads [6,14,15,16,17,18,19,20]. In this context, the Zeeman deceleration method relies on the state-dependent interaction of neutral paramagnetic atoms or molecules with time-dependent inhomogeneous magnetic fields [21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37] and is thus suitable for open-shell systems such as molecular radicals or metastable atoms and molecules [38,6,5].
Here, we present a novel traveling-wave Zeeman decelerator recently developed in our laboratory. The new decelerator operates with traveling magnetic traps containing molecules which are adiabatically slowed down, which is conceptually similar to the traveling-wave Stark decelerator [39,40,41,42,43,44] but works for paramagnetic atoms and molecules. The present paper focuses on the theory and numerical analysis of the operational principle of the device. The experimental implementation and its demonstration are described in the accompanying article [45]. As areas of application of the new decelerator we envisage cold-collision experiments, trap loading for further cooling and studies of the chemistry of trapped particles at very low temperatures.

Longitudinal dynamics of the moving trap
The decelerator features a modular design [45] with each module consisting of 32 copper wires wound around a cylinder in two layers where the inner (outer) layer consists of 16 right-handed (left-handed) helices, see Figure 1. The wires are labelled as indicated in the figure. By supplying time-varying currents to the double-helix coil geometry, a traveling magnetic wave along the coil axis is formed. To illustrate the principle, we start by analyzing the magnetic fields generated by a single helix within a deceleration module. To a good approximation, the magnetic field on the axis of a coil generated by nn + Figure 1. Schematic of a decelerator module consisting of two stacked layers of 16 wires wound in right-and left-handed helices. The two layers are depicted separately in blue and yellow for clarity. n − (n + ) ∈ {0, .., 15} labels the individual wires in the outer, left-handed (inner, right-handed) helix layer.
z y x Figure 2. (a) Schematic of a single helical wire in a decelerator module. (b),(c) and (d) Magnetic field components B +x0 (z), B +y0 (z) and B +z0 (z) along the central coil axis (z, x = y = 0) generated from a single wire. Data points represent the results of numerical calculations based on the Biot-Savart law (see text). The lines are a fit of the data to Equation (1) in the text. a current of 1 A supplied to a right-handed wire helix can be expressed as B +x0 (z) = a + sin(kz), and to a left-handed helical coil as The x i = x, y, z components of the magnetic-field vector generated by a single helical wire, B +/−x i n , are shown in Figure 2. The subscript n = 0, ..., 15, denotes the wire index and + (−) stands for a right-(left-)hand orientation of the wires. k = 2π/λ is the wave number of the traveling wave, where λ = 14 mm is the periodicity of the helices in the present implementation. The parameters a +/− represent the amplitude of the magnetic field at a probe current of 1 A. Their values are derived in the supplementary material, where also the validity of Equation (1) and Equation (2) is demonstrated. It is straightforward to obtain equivalent expressions for the magnetic fields produced by the other wires with a rotational transformation by an angle n∆ where ∆ = 2π/16 due to 16 wires being distributed evenly around a cylinder. For the right-and left-handed helices, the components of the magnetic field are accordingly given by B +xn (z) = a + sin(kz + n∆), B +yn (z) = a + cos(kz + n∆), Time-dependent currents of the form are applied to the n-th wire, where φ +/− (t) is the time-dependent phase of the currents supplied to the right-(left-) handed layer. The additional phase n∆ serves as a compensation to the geometrical arrangement of the wires in order to minimise higher-order harmonics in the synthesised magnetic-field profiles. c +/− are the current amplitudes for left-(+) and right-(-) handed wires, respectively, and can be tuned individually. The magnetic-field components generated by the entire double-layer coil geometry are expressed as c + sin φ + (t) + n∆ a + cos(kz + n∆) By introducing new parameters A avg and ∆A which are the sum and the difference, respectively, of the geometry-weighted current amplitudes applied to the two layers, c + a + and c − a − can be expressed as The phases φ + (t), φ − (t) can be written as The ratio c + a + /c − a − depends on the radius and the current amplitude of each layer and can be tuned in a way such that ∆A ≈ 0. Inserting Equation (9) and Equation (10) into Equation (6) and Equation (7) and assuming ∆A ≈ 0 yields where N = 16. Thus, the magnitude of the magnetic field on the central axis of the full double-helix assembly reduces to the form of a traveling wave: As can be seen in Equation (11) and Equation (14)), the dynamics of the magnetic wave is controlled by the two time-dependent phases ∆φ(t) 2 and φ avg (t). The phase ∆φ(t) serves to control the position of the wave along the coil axis. By applying a time dependence of the form  the traveling magnetic wave can be programmed to decelerate, accelerate or propagate with constant velocity along the z coordinate. The parameters A, B and C are given by A determines the deceleration (acceleration) and B encodes the initial velocity v I of the traveling wave. The constant phase shift C can be set to zero for practical purposes. v F stands for the final velocity of the traveling wave and L is the length of the decelerator. The on-axis trap dynamics as a function of time-dependent phase ∆φ(t) 2 is illustrated in Figure 3 (a) where the propagation of the magnetic wave is illustrated by plots of the magnetic field for four values of the phase ∆φ 2 = π 16 , 5π 16 , 9π 16 and 13π 16 . So far, only the magnetic-wave dynamics on the central coil axis was discussed. The transverse, i.e., perpendicular to the central axis, components of the magnetic field were calculated numerically as explained in the supplementary material. The transverse dynamics of the magnetic wave is controlled by the time-dependent phase φ avg (t). This is illustrated in Figure 3 (b) where the calculated magnetic field |B(x, y, 0, t)| is plotted for four different values ∆φ avg = 0, π 4 , π 2 and 3π 2 . As can be seen in Figure 3 (b), the magnetic field in the transverse direction exhibits a deep minimum along one direction capable of confining magnetically low-field-seeking species. In the perpendicular direction, the minimum is very shallow so that particles can escape along this axis. Changing the phase ∆φ avg varies the orientation of the trap in the transverse direction. This effect can be understood as follows: assume a stationary on-axis magnetic wave, i.e., ∆φ 2 (t) = 0 , then the magnitude of the magnetic field in the transverse direction (at z = 0) generated by the n-th right-and left-hand-oriented helix is given by where B +/−n (x, y, 0) is the magnetic field in the xy plane (z = 0) generated by the probe current supplied to the n-th helix in right (left) layer. The right-and left-handed n-th helices are maximally contributing to the total magnetic field when ∆φ avg +n∆ = kπ 2 , k = 0, 1, 2, ... If the phase ∆φ avg (t) changes with time as ∆φ avg = ω avg t and an angular frequency ω avg , the time at which the n-th helices are maximally contributing is t n = kπ−n∆ 2ωavg . At the time t n+1 = t n + ∆ 2ωavg , the maximal contribution comes from the (n + 1)-th helices and the transverse field was rotated by an angle ∆. In this way, the transverse field is rotated with time at the angular frequency ω avg . The two phases ∆φ 2 (t) and φ avg (t) can be chosen independently from each other, leading to a decoupling of the longitudinal and transverse motions of the trap. This feature prevents trap losses due to the motional coupling which are characteristic of the conventional Stark and Zeeman decelerators [46,47]. Experimentally, both phases are independently controllable through the time-dependent phases of the supplied currents φ + (t) and φ − (t), Equation (10). As an example, numerically calculated magnetic fields in zx, zy and xy planes at t = 0 generated by time-varying currents as given by Equation (5) and a current amplitude I = 300 A are shown in Figure 4 (a)-(c). The magnetic field exhibits confining characteristics along the z and y coordinates, but only weak trapping along the x coordinate as discussed above (see also Figure 5).

Transverse stability
In this section, the transverse stability of particles confined in a rotating trap is discussed. The equations of motion along the x and y directions are given by the coupled differential equations Here, B 0 x(t), y(t), z(t), t is the magnetic field generated by a probe current of 1 A supplied to the coil assembly. The magnetic field generated by an arbitrary current is obtained by scaling the α is a dimensionless parameter given by α = dc, where d = µ eff M is the magnetic-dipole-moment-to-mass ratio, µ eff is the effective magnetic dipole moment and M is the mass of the particle (in atomic units). µ B is the Bohr magneton and m u is the atomic mass constant. To discuss the dynamics of the particles in the rotating trap, the magnetic field |B 0 | along the x and y coordinates at t = 0 are plotted in Figure 5. The calculated magnetic field is shown by the red lines. A polynomial fit using the function f fit = a 0 + a 2 x 2 i + a 4 x 4 i + a 6 x 6 i + a 8 x 8 i + a 10 x 10 i , x i = x, y is shown by the black dash-dotted lines. By visual inspection of Figure 5 and from the fitted parameters listed in Table 1, it becomes apparent that the magnetic field along both directions exhibits a high degree of anharmonicity with 4-th-and 6-th-order terms still being significant. As a result of the anharmonicity, analytical Table 1. Polynomial fit coefficients for the magnetic field |B 0 | shown in Figure 5 along the x and y coordinates.
x coordinate y coordinate a solutions to the coupled differential equations Equation (18) and Equation (19) do not exist. Therefore, the stability of the solutions as a function of the parameter α and the rotational frequency ω avg was examined numerically. Solutions were propagated in time up to t = 50 ms using an adaptive 4-th order Runge-Kutta algorithm. The positions of the particles at the end of the propagation were examined. If the position was found to lie within the trapping region, the trajectory was considered as stable, and otherwise as unstable. For given parameters α and ω avg , the equations were solved numerically and the parameter β = B max 0,α x(t),y(t),0,t B max 0 x(t),y(t),0,t was calculated. The pointed brackets denote a time average. β represents the ratio of the maximum of the time-averaged magnetic field B 0 for which the trajectory of a particle is stable and the time-averaged maximum of the magnetic field B 0 , or in other words, the ratio of the effective trapping potential and the maximum trapping potential in the limit ω avg → ∞ for a given particle. The resulting stability diagram of the trap is shown in Figure 6. Solutions of the equations of motion were explored for a frequency range ω avg = 2π × (0 − 10) kHz, which corresponds to the operational range of the experimental decelerator. The dependence of the parameter β on the parameters α and ω avg is illustrated by a color map in Figure 6. Regions of unstable trajectories appear in black. Seven distinct regions were found within the parameter space studied, four regions with stable and three with unstable trajectories. White dashed lines delineate the borders between stable and unstable regions, and were produced by fitting a harmonic function kω 2 avg to the data where k is a fitting parameter. The highest stability (maximum β) is achieved for low values of α where trajectories are stable across the majority of the range of frequencies. With increasing α, the threshold frequency indicated by the thick white dashed line above which there are no regions of unstable behaviour increases and is proportional to √ α. The existence of regions of stable trajectories both at low and high rotational frequencies allows, in principle, for a selective confinement of multiple molecular species in the deceleration process. This property is illustrated by the example of the simultaneous trapping of H atoms in the 1 2 S 1/2 (m s = 1/2) state (blue dash-dotted line) and OH molecules in the X 2 Π 3/2 (ν = 0, m J = 3/2) state (red dash-dotted line). At low rotational frequencies (e.g. ω avg = 2π × 0.

Numerical trajectory simulations
In order to model the particle dynamics within the entire decelerator and characterize its deceleration performance, a numerical trajectory simulation code in the Python programming language was developed. In the simulation, 10 5 OH molecules were initialised and their trajectories through the decelerator were calculated under different operating conditions. The forces F i acting on each particle i along its trajectory x i (t)) were obtained from its Zeeman energy W i according to F i (x i (t)) = −∇W i (x i (t)). The Zeeman energy level structure of OH [48] in its X 2 Π 3/2 , v = 0 ground electronic and vibrational state in a magnetic field was calculated according to where Λ and Σ are the quantum numbers of the projection of the electron orbital and spin angular momenta on the molecular axis, g s is the electron-g factor (g s ≈2), and M J and Ω ef f are the quantum numbers of the projection of the total angular momentum J (with quantum number J) along the external magnetic field and the molecular axis, respectively [48]. The Zeeman energy level structure for the lowest rotational level, J = 3/2, is shown in Figure 8. The low-field seeking quantum states with the largest Zeeman energy at a given magnetic field strength which can be decelerated are X 2 Π 3/2 , |Ω| = 3/2, e and f, M J = 3/2. In the simulations, only molecules in the M J = 3/2 state were considered. The molecules were initially generated within a cylinder with a radius of 2 mm and a length of = 11 mm, which mimics the initial spatial distribution of our experimental molecular beam [49], with a transversal velocity spread of 45 m/s and a longitudinal velocity spread of 90 m/s. The mean forward velocity was initially chosen to be 450 m/s in line with our experiments [45]. Thus, the particles fill an initial 6-dimensional phasespace volume of 0.5 × 10 8 mm 3 (m/s) 3 . Due to the time-dependence of the currents and complex structure of the wire geometry, the magnetic fields could not be fully calculated prior to a simulation run and were thus calculated on the fly. In order to reduce the computation time, the calculations were GPU-accelerated using the NumbaCUDA library in Python [50]. In the simulations, the decelerator was assumed to be 1.792 m long (32 deceleration modules [45]) followed by a 15 mm region of free flight after which the molecules were detected. Exemplary results of time-of-flight profiles of the molecules extracted from the trajectory simulations are displayed in Figure 9. Here, the molecules were decelerated to nine different final velocities in the range 449 m/s -50 m/s corresponding to decelerations of 0.25 km/s 2 -55.80 km/s 2 . Depending on To illustrate the phase-space stability of the new deceleration method, we show in Figure 10 (a)-(e) 2D longitudinal phase-space distributions of particles in a single trap at the end of the decelerator for decelerations ranging from 0.25 km/s 2 to 55.8 m/s 2 corresponding to final velocities in the range 449 m/s to 50 m/s. Panel (f) corresponds to a deceleration of -27.6 m/s 2 , i.e., an acceleration to a final velocity of 550 m/s. A normalized particle density is shown as a color map, the calculated separatrices delineating the phase-stable regions are represented by white solid lines. As evidenced by Figure 10, the phase-space volume transported through the decelerator decreases with increasing deceleration, as is also the case with previous implementations of Starkand Zeeman decelerators [46,47].
Additionally, we explored the 2D phase-space dynamics in both transverse directions for the case of a stationary transverse trap (ω avg = 0 kHz) and a trap rotating at ω avg = 10 kHz. Comparatively long simulation times of 50 ms were chosen to allow the particles to fully explore the phase space. Figure 11 illustrates transverse phase-space distributions extracted from the trajectory simulations. The density of the particles is represented by a color map. Figure 11 (a) and (b) shows the phase-space distributions of   in agreement with transverse magnetic field distribution outlined in Section 2.2, Figure 4 and Figure 5, where it was shown that stable trapping regions exist for both ω avg =0 kHz and 10 kHz along the y direction. Along the x direction at a rotational frequency ω avg =0 kHz, the phase-space distribution shows free-flight characteristics due to very small confinement. At a rotational frequency ω avg =10 kHz, there is no appreciable difference in phase-space structure between the x and y directions as the particles are trapped in the same time-averaged potential along both directions during the deceleration process.
In addition, the complete 6D phase-space acceptance was explored numerically for a range of decelerations. Results of these calculations are shown in Figure 12 (a). The 2D longitudinal phase-space acceptance is shown by the black line, while the transverse phase-space acceptances both for a stationary (ω avg = 0 kHz) and rotating (ω avg = 10 kHz) trap are indicated by the blue lines. The full 6D phase-space acceptance for a range of decelerations is shown in Figure 12 (b). An upper limit of the complete 6D phasespace volume is approximated by a product of the 2D volumes along each coordinate. The dashed traces in Figure 12 (b) show the thus calculated maximum 6D phase-space acceptance for ω avg = 10 kHz (blue dashed line) and ω avg = 0 kHz (red dashed line). Additionally, the full 6D phase-space volume at the end of the deceleration region was extracted from the numerical trajectory simulations (blue and red full lines in Figure 12 (b)). The difference between the phase-space acceptances estimated from the product of the 2D acceptances and the ones obtained from the trajectory simulations is attributed to transverse instabilities of trajectories far from the centre of the trap. The maximum phase-space acceptance volumes calculated for the highest deceleration (55.8 kms 2 corresponding to the final velocity of around 50 ms) are 8 · 10 4 mm 3 m/s 3 for ω avg = 10 kHz and 5 · 10 4 mm 3 m/s 3 for ω avg = 0 kHz. For comparison, the maximum phase-space acceptance volume of a typical Stark decelerator at a final velocity of 50 m/s is 8 · 10 3 mm 3 m/s 3 [51,52], and 2 · 10 7 mm 3 m/s 3 for the Zeeman decelerator reported in [28]. The phase-space acceptance of the current assembly could be further enhanced by increasing the currents supplied to the coils.

Conclusions
We have developed the theory describing the working principle of a novel traveling-trap Zeeman decelerator including analytical expressions for the magnetic fields generated from a double-helix coil geometry in one dimension. Two independent parameters governing the dynamics of the moving magnetic traps were introduced. The transverse stability of the decelerator was explored in detail. It was shown that trap stability in relation to the rotation frequency ω avg enables the selectivele trapping of molecules according to the magnetic-dipole-moment-to-mass ratio and the current amplitude supplied to the decelerator. The deceleration efficiency is limited by the deceleration applied to the molecules which is an intrinsic property of decelerators based on conservative forces.