Zeeman deceleration beyond periodic phase space stability

In Zeeman deceleration, time-varying spatially inhomogeneous magnetic fields are used to create packets of translationally cold, quantum-state-selected paramagnetic particles with a tuneable forward velocity, which are ideal for cold reaction dynamics studies. Here, the covariance matrix adaptation evolutionary strategy is adopted in order to optimise deceleration switching sequences for the operation of a Zeeman decelerator. Using the optimised sequences, a 40% increase in the number of decelerated particles is observed compared to standard sequences for the same final velocity, imposing the same experimental boundary conditions. Furthermore, we demonstrate that it is possible to remove up to 98% of the initial kinetic energy of particles in the incoming beam, compared to the removal of a maximum of 83% of kinetic energy with standard sequences. Three-dimensional particle trajectory simulations are employed to reproduce the experimental results and to investigate differences in the deceleration mechanism adopted by standard and optimised sequences. It is experimentally verified that the optimal solution uncovered by the evolutionary algorithm is not merely a local optimisation of the experimental parameters—it is a novel mode of operation that goes beyond the standard periodic phase stability approach typically adopted.


Introduction
The precise control of the velocity of beams of atoms and molecules using Stark and Zeeman decelerators has led to multiple benefits in a range of fields spanning atomic and molecular physics, including high resolution spectroscopy [1][2][3], quantum-state resolved collisional scattering [4][5][6], lifetime measurements [7,8] and the exploration and exploitation of properties of cold matter [9][10][11][12][13]. Travelling-wave Zeeman decelerators have seen the velocities of paramagnetic species manipulated by confining these species in a true three-dimensional travelling well [14,15], with Zeeman-decelerated species also successfully loaded into magnetic traps [16,17]. For all such applications, the optimisation of the density and number of particles in the decelerated beam is an important factor contributing to the experimental viability, in addition to attaining a narrow velocity distribution. In this paper we describe the use of an evolutionary strategy approach to optimise the timed sequence of magnetic fields applied in a Zeeman decelerator to achieve maximum throughput of velocitycontrolled decelerated atoms.
Zeeman deceleration is an experimental technique in which a supersonically expanded beam of paramagnetic particles, such as radicals or metastable electronically excited species, typically passes through a sequence of solenoid coils. The forward velocity of the beam is reduced by exposing it to a tailored sequence of magnetic field pulses. The presence of unpaired electrons in paramagnetic species creates the magnetic dipole moment which couples to the applied magnetic field, lifting the degeneracy of the magnetic sublevels. Those levels for which the energy increases (decreases) with the field are known as low-field seeking (LFS) (high-field seeking (HFS)) states. Zeeman deceleration takes advantage of the interactions that particles in these magnetic sublevels have with time-varying, spatially inhomogeneous magnetic fields: only paramagnetic species in selected LFS quantum states are decelerated [13,18].
Conventionally, the deceleration switching sequence-that is, the timing of the switching on and off of the magnetic field inside each solenoid coil composing the decelerator-is determined by considering the trajectory of an idealised 'synchronous' particle. The trajectory of this synchronous particle under the influence of the magnetic fields is calculated, and the switching time of each solenoid coil is determined by when the particle reaches a selected axial position within that coil relative to the coil centre. This position is normally expressed as a phase angle [9] and conventionally such devices are operated at contant phase angle (and constant maximum current) throughout the deceleration sequence. Changing the selected phase angle enables a range of different final velocities of the decelerated species to be attained, albeit with variable transmission efficiency. The switching sequences produced by this standard approach are capable of decelerating other particles close to the synchronous particle in phase space, but they are not necessarily the optimal sequences; there are limitations to using phase angle alone when seeking to control the properties of the resulting beam. The use of a varying phase angle (evolving through the sequence) in order to keep the deceleration per stage constant has been investigated, and it was demonstrated that this could enhance the output beam density [16,[19][20][21]. Given that the individual solenoid pulse timings can be varied arbitrarily and independently in a given apparatus, it is possible that more complex timing sequences could beneficially adapt the characteristics of the output beam for specific applications. In order to fully optimise the switching sequences to achieve the maximum number of decelerated particles, we must consider more than just the synchronous particle. A suitable approach to investigate the full optimisation of the sequence, subject to certain practical constraints, is to use an evolutionary strategy.
Evolutionary strategies are a class of numerical optimisation techniques, drawing inspiration from Darwinʼs evolutionary theory of natural (genetic) selection. Since genetic algorithms were popularised by Holland in 1975 [22] as a way of solving computationally intractable problems, genetic and evolutionary algorithms have been implemented in a wide range of disciplines spanning the mathematical, physical and medical sciences. An abstract adaptation of biological genetic processes-selection, recombination, evolution and mutation-is employed to modify a set of parameters that influence a measurable outcome, and to locate the optimal solution for that set of parameters required to achieve the closest match to the desired outcome. From an initial set of parameters, an evolutionary algorithm carries out a process of fitness-based ranking and selection, followed by recombination to produce a successor population of parameter sets-the next generation. During recombination, the best performing parent 'chromosomes' are preferentially selected and their 'genetic material' is recombined and passed on to the successor population. As this process is iterated, successive generations are formed and the average fitness of the potential solutions typically improves until convergence is achieved (or some stopping criterion is reached). In this way, an evolutionary algorithm 'evolves' the optimal solution to a given problem. A covariance matrix adaptation evolutionary strategy (CMA-ES) is the evolutionary algorithm adopted in this work [23].
Applying this approach to Zeeman deceleration, each trial solution will be a set of values for the parameters that define the switching sequence of the coils (i.e.pulse durations and timings), and we can define the fitness of a given sequence as the fraction of particles that exit the decelerator with a longitudinal velocity within, for example, 10 m s 1  of the target velocity. To achieve optimisation, we need to look at more than a single synchronous particle; we must consider the trajectories of the full set of particles representing the initial distribution in phase space through the decelerator apparatus. For a constant current of 243A applied to the solenoid coils, optimisation of the pulse durations represents a 12-dimensional problem for the apparatus in this work. While iterative steepest-decent methods, such as the Levenberg-Marquardt algorithm [24][25][26], are useful when the initial parameter set is close to optimal, stochastic methods are better suited at locating the global minimum of high-dimensional systems. Furthermore, evolutionary strategies do not require any information about the gradient of the target function at a given set of parameters, meaning that the function itself need not be continuous or differentiable. Evolutionary strategies are thus broadly applicable and typically cost-effective to implement. While evolutionary algorithms have already been shown to increase the number of decelerated particles in both Stark [27] and Zeeman [19] deceleration, their full potential was perhaps not fully exploited in these earlier studies. Wiederkehr et aladopted a CMA-ES approach to maximise the number of Zeeman-decelerated particles whilst minimising their spread in velocity and position. They reported an increase in the intensity of decelerated D atoms-attributed to improved phasespace acceptance and optimising the balance between transverse focusing and defocusing throughout the deceleration process-but concluded that the evolutionary strategy-optimised sequences were not a significant improvement over the conventional constant nominal phase angle mode of operation [19].
In this work, we build on previous studies and utilise the CMA-ES evolutionary algorithm to fully optimise the operation of a Zeeman decelerator for collision studies. We describe how evolutionary algorithms can generate deceleration sequences that not only achieve more efficient deceleration than a comparable standard sequence-yielding a higher number of decelerated particles within a given narrow velocity range-but also enable significantly lower final velocities to be achieved than would be possible with a standard sequence. We demonstrate how evolutionary algorithms can be employed to enhance our ability to manipulate the resulting decelerated beam, exercising far more control over the beam than is possible through phase angle manipulation alone. Finally, we investigate the deceleration mechanism employed by the standard and optimised sequences, identifying the key differences between the two approaches and how the optimised sequences yield such significantly improved results.

Methods
The Zeeman decelerator apparatus (figure 1), previously described in [20], consists of two differentially pumped high-vacuum chambers separated by a 2.0 mm diameter skimmer. In the first (source) chamber, hydrogen atoms are generated within a supersonic expansion. The second (detection) chamber contains the decelerator solenoid coils and time-of-flight (ToF) mass spectrometer for deceleration and detection of the H atoms, respectively.
Hydrogen atoms are generated in a pulsed supersonic beam by ArF excimer laser photolysis at 193nm of a beam of NH 3 seeded at a 1:9 ratio in Kr. The photolysis occurs within a quartz capillary fixed over the beam orifice on the face plate of a pulsed solenoid valve. The H atoms and other photofragments expand out of the capillary and travel towards the skimmer. The valve body is cooled to 238 K by a flow of cold nitrogen gas, reducing the peak in the velocity distribution of the gas packet to 500 m s −1 . A subset of the hydrogen atoms (in LFS states) are decelerated by applying a sequence of timed current pulses to the decelerator solenoid coils, synchronised relative to the photolysis laser pulse.
After passing through the decelerator, ground state H atoms are ionised by (2+1) resonance-enhanced multi-photon ionisation (REMPI) at 243nm via the 2s S 2 1 2 ( )state, with the resulting ions detected by a microchannel plate detector. The velocity distribution of the neutral H atoms is established by scanning the delay between the photolysis laser and the REMPI laser to collect traces that are referred to here as ToF profileswith the ToF being the flight time of the H atoms between the axial positions where the beam is intersected by the two laser pulses. These ToF profiles are recorded both with and without the magnetic fields applied in alternating shots, achieved by switching the current pulse sequence on and off. The spatial distribution of hydrogen atoms leaving the decelerator is recorded by translating the REMPI laser beam along the molecular beam axis by a micrometre-adjustable periscope. This is illustrated in figure 1, where the range of REMPI laser positions is explicitly indicated.

Zeeman deceleration
The Zeeman decelerator is composed of a series of 12 solenoid coils that are rapidly pulsed with high currents (243 A), creating a maximum on-axis magnetic field of up to 1.8 T within the coils. In the presence of an external magnetic field, the energy levels of atomic hydrogen are split. In ground state hydrogen 1s S 2 1 2 ( ), the interaction between electron spin and nuclear spin gives rise to hyperfine splitting, yielding two LFS states and two HFS states (figure 2(a)). As the beam of paramagnetic particles travels toward the centre of each coil, particles experience a positive magnetic field gradient and the potential energy of particles in LFS states increases at the expense of their kinetic energy ( figure 2(b)). These LFS particles can be permanently slowed down if the current is switched off before the particles reach the centre of the coil. In this way, each coil removes an amount of kinetic energy equal to the Zeeman energy of the state at the position when the coil is switched off. Thus a precisely timed switching sequence applied to the coils can slow a packet of quantum-state-selected H particles to an arbitrary final velocity.
Particle trajectory simulations Three-dimensional numerical particle trajectory simulations serve to both guide the selection of standard experimental switching sequences and to aid interpretation of the experimental ToF profiles. An arbitrary number of particles in specific quantum states can be simulated, with their initial positions and velocities randomly drawn from a normal distribution with parameters that match the experimental conditions. We assume that particles are equally distributed amongst the four Zeeman sub-levels of H 1 S 2 1 2 . The geometry of the simulated set-up-in addition to the shape of the current pulse profiles-is based on the experimental parameters. Simulations are conducted with and without the magnetic fields present within the coils; when a current is applied to the coils, the resulting magnetic field properties are derived from the analytical solution for a current loop [28]. A velocity Verlet algorithm integrates the equations of motion every 100ns, with acceleration forces calculated using the particles' quantum state and the magnetic field distribution. Particles are removed from the simulation when they strike a surface (such as a decelerator coil) or reach the detection region, defined as the volume of space illuminated by the detection laser. The detection laser is aligned with the x axis (perpendicular to the atomic beam, which propagates along the z axis) and is modelled as a Gaussian in the z and y dimensions, with a standard deviation of 0.8mm and 0.4mm respectively. The simulated signal is given by the number of trajectories passing through the laser detection volume in a given time interval.

Generation of standard pulse sequences
The conventional, or standard, current pulse sequence applied to the decelerator coils is intuitive: it consists of pulses of monotonically increasing duration applied to each coil in sequence (see top panel in figure 3(a)). This arises because standard switching sequences calculated from the trajectory simulations, as decribed above, consider the passage of a single synchronous hydrogen atom in the M F =1 state through the apparatus (where M F is the projection of the total angular momentum F onto the local magnetic field axis). The idealised synchronous particle travels along the axis of the decelerator. When the synchronous particle reaches a desired axial position within each coil the magnetic field within the coil is switched off. This is repeated for all twelve coils, and in this way a complete switching sequence is constructed. As the synchronous particle travels through each coil, it gets increasingly decelerated and therefore it takes it longer to reach the switch-off position in the next coil, thus giving rise to the sequential increase in successive pulse durations.
As can be seen in figure 3(a), the currents applied to each coil in our apparatus exhibit non-zero rise and fall times of 7 μs. This is because the magnetic field strength changes at a finite rate, dictated by the time constant of the inductor-resistor circuit formed by the coil and its switching electronics. The switch-on time of each coil is chosen so that the rising field of a coil is overlapped in time with the decaying field of the previous coil; an overlap time of 6 μs is adopted here. This temporal overlap of adjacent pulses maintains a non-zero quantisation magnetic field, thus avoiding Majorana spin-flip transitions to HFS states. Maximum deceleration is achieved at a phase angle of  f denotes the position of the particle when the coil is switched off (i.e.the position of the particle when the current reaches zero, termed the effective phase angle). synchronous particle lie in a phase-stable region and are also decelerated, forming a compact bunch that travels through the decelerator. The size of this phase-stable region depends on the choice of 0 f : the closer 0 f is to 90°, the greater the kinetic energy that is removed at each stage but, at the same time, the smaller the size of the stable phase-space volume, which means that fewer particles are decelerated. Non-standard modes of operation have improved the transverse phase stability of Zeeman deceleration, increasing the density of particles confined to the plane perpendicular to the molecular beam axis by periodically assigning deceleration and focusing coils [20]. As first demonstrated with a Stark decelerator [29], the interplay between longitudinal and transverse motion can reduce the phase space acceptance of the decelerator stages and thereby limit the efficiency of the deceleration process. Operation in, for example, the s=3 configuration (where only every third stage is used for deceleration, with the other stages focusing the beam) effectively decouples the longitudinal and transverse dynamics, increasing the number of decelerated particles that can be obtained. There are, however, limitations associated with using modes of operation such as s=3; with some coils facilitating transverse focusing in place of deceleration, far less kinetic energy can be removed compared with the standard mode of operation. Indeed, it has been previously proposed that the s=3 mode of operation is not feasible for Zeeman deceleration [21,30]. Recent work has seen the combination of a series of alternating magnetic hexapoles and solenoid coils in a multi-stage Zeeman decelerator, enabling phase stable operation over a wide range of velocities for applications in scattering experiments [31]. Here, we describe an alternative approach: optimising pulse sequences using evolutionary algorithms.
Evolutionary optimisation of pulse sequences Evolutionary algorithms seek to optimise a set of parameters with respect to a fitness function: the 'fitness' of each set of parameters (which give rise to a deceleration sequence) in a population is computed, after which a new generation of sequences is created based on the results from the previous generation. The parameters that are optimised here are the durations of the pulses of current applied to each of the twelve coils. A full switching sequence is generated from a set of durations by applying the necessary constraints. Each pulse must overlap with the previous one to maintain a non-zero quantisation magnetic field, as described above for the standard sequences (the same overlap time of 6 μs is adopted). Additionally, for the apparatus employed here, each pulse must not exceed 60 μs in duration. The upper limit of the pulse duration (60 μs) arises from the limited cooling efficiency of the solenoid coils, given the large amount of heat that needs to be dissipated when repeatedly pulsing such high currents (243A) at 10Hz. The lower limit of the pulse duration is 14 μs, given by the 7 μs rise and fall times of the current pulses. The total number of particles arriving within 10 m s 1  of the target velocity is optimised, with the fitness of each sequence evaluated by simulating the trajectories of 50 000 particles in the M 1 F = Zeeman sublevel and recording the number that reach the detection region with a velocity within the target range.
A CMA-ES approach is used to optimise the switching sequence so as to maximise the number of decelerated particles reaching the laser detection volume. A standard pulse sequence (the top trace in figure 3( = in the evolutionary algorithm approach. A population of λ candidate sequences is obtained, where the pulse durations that compose each sequence are sampled from twelve multivariate normal distributions, one for each coil. The mean and standard deviation of each of these distributions are computed from the preceding generation or, in the case of the first generation, from the durations of the standard pulse sequence that is being optimised with an arbitrarily selected standard deviation of 10 μs. The constraints on pulse duration and switching times are implemented by replacing any pulse duration that exceeds (is lower than) the maximum (minimum) pulse duration with the maximum (minimum) permitted pulse length of 60 μs (14 μs); this method was found to be more efficient than simply discarding any sequences that fall outside the required constraints, or imposing a penalty proportional to the distance the candidate sequence lies outside the permissible region. The sequences thus generated are evaluated and ranked in fitness order based on the total number of particles decelerated to within the target velocity range. The mean values and standard deviations of the multivariate distributions from which the durations are drawn are updated at every generation, and they are calculated from the weighted average of the μ best performing individual sequences. Typical settings involve 2 m l  , corresponding to half of the sequences being discarded and the other half being used to update the distributions and normalised such that the weight of an individual sequence and is proportional to its ranking. The updated mean values and standard deviations are then used with the covariance matrices to generate twelve new multivariate normal distributions from which λ sequences are created for the next generation, achieving the selection and recombination characteristics of an evolutionary algorithm.
The covariance matrix update is at the core of the CMA-ES method. An accurate estimate of the covariance matrix requires a large population of parameter sets. However, this comes at the expense of significantly increased computing time. A successful compromise-yielding faster convergence with high reliability-is achieved by using information from the best performing sequences in previous generations to update the covariance matrix. In this way, the best performing sequences of the first generation are used to create the second generation of sequences, with the worst performing sequences discarded. In essence (after the first generation), the distibution of pulse durations for each coil is established using the mean and standard deviation of the best performing sequences from the previous generation. This process of parameter optimisation means that the CMA method is nominally parameter-free, only requiring an initial choice of mean and standard deviationvalues which are quickly adapted by the optimisation process. Throughout this iterative process, the globally best performing sequence is stored and replaced each time a better sequence is found. The optimisation process is deemed to have converged when the best performing sequence has not been updated for a given number of generations, chosen to be 50 000 generations in this work. The final best performing sequence is then selected and utilised experimentally.
Optimising for sequentially lower target velocities (below v 200 z final = m s −1 , the second trace in figure 3(a)), the nearest optimised evolutionary algorithm sequence is taken as the starting sequence from which the new first generation is created. Thus the evolutionary algorithm-optimised sequence for v 200 m s = -(bottom trace in figure 3(a)); attempts to realise lower final velocities fail to yield any particles in the target velocity range during the first generation of sequences and so cannot be optimised.

Comparison of sequences
As figure 3(a) clearly illustrates, evolutionary algorithm optimisation yields markedly different sequences compared to those generated by the standard trajectory simulations (which consider a single synchronous particle). The standard sequence (top trace) features monotonically increasing current pulse lengths applied to each coil; as the synchronous particle slows down, it takes longer to travel through each successive coil. As has already been noted, the maximum experimentally feasible current pulse duration and magnetic field gradient across all coils imposes an upper limit on the final velocity that can be achieved with standard deceleration sequences; starting from an initial velocity of 500 m s −1 , H atoms in the ground state can be decelerated to a minimum final velocity of 205 m s −1 .
In contrast, sequences produced by evolutionary algorithms have pulse durations that vary in an oscillatory manner across consecutive coils. This results in some earlier coils being switched on for a longer time, enabling lower final velocities to be achieved compared to the standard sequence. The evolutionary algorithm sequences also see very short pulses applied to coils 4 or 5, with this mid-decelerator coil aiding longitudinal focusing (rather than playing a primarily decelerating role). Thus for the same constraints of minimum and maximum pulse length, and the same incoming beam, evolutionary algorithm sequences can decelerate H atoms down to a velocity of 75 m s −1 (see figure 3(a), bottom trace)-some 130 m s −1 lower than the minimum velocity achievable with the standard sequence. This is due to the fundamentally different approaches taken when generating the standard versus evolutionary algorithm sequences. In essence, the evolutionary strategy optimises sequences considering the properties of the entire beam and the passage of all particles through the decelerator, whereas the standard sequence considers the trajectory of only one synchronous particle.

Experimental results
The current profiles generated by the standard sequence and a range of evolutionary algorithm sequences are implemented experimentally, with the resulting ToF traces presented in figure 3(b). The experimental ToF traces (solid lines) are shown beside the corresponding pulse sequence in column (a). ToF traces recorded with no deceleration (i.e. with no current applied to the decelerator coils, referred to here as 'undecelerated' traces) are shown in red, with the experimental traces shown as a solid red line. The broad, featureless profiles of the undecelerated beams are the same across all panels, in excellent agreement with the simulated ToF profiles (shaded red traces in figure 3(b)), confirming that the properties of the initial beam and the decelerator geometry are well described in the simulations.
The 'decelerated' ToF traces (i.e.where current is applied to the coils, solid blue lines) show a characteristic multi-peak profile. The first and most intense peak (with the earliest arrival time) corresponds to the part of the beam that gets focused as it travels through the decelerator-without being significantly decelerated or accelerated. The final peak corresponds to the part of the beam that gets decelerated to the target velocity, with the arrival time of this final peak having a linear relationship with the inverse target velocity of the sequence. The features between the first and final peaks in the decelerated ToF traces arise from particles that are partially decelerated by the magnetic fields. The multiple peaks and features of the experimental decelerated ToF profiles are well described by the simulations (shaded blue profiles in figure 3(b)). Importantly, the experimental ToF traces provide clear evidence that the optimised sequences are feasible; decelerated particles are detected at the anticipated ToF in all cases, even for velocities down to 75 m s −1 . In all instances, the intensity and spread of the decelerated peak is reproduced by the simulations.
For each deceleration sequence, ToF traces are recorded at five different detection distances (as indicated by the series of REMPI detection laser positions in figure 1) to facilitate an independent verification of the velocity of H atoms contributing to the final peak. One such series of ToF profiles at a range of detection distances is provided in figure 4(a), corresponding to the evolutionary algorithm sequence with a target final velocity of 200 m s −1 . By plotting the change in detection distance against the change in ToF (schematically represented by the dashed green line across the panels in figure 4(a)), one can estimate the mean experimental velocity of the particles contributing to the decelerated peak. With this approach, the velocity is experimentally determined to be (210±10)m s −1 , in agreement (within error) with the target velocity of (200±10)m s −1 utilised in the optimisation procedure. The slight overestimate in the experimental velocity established with this simple approach is attributed to the neglect of dispersion, which causes the peak in density to move with a faster apparent velocity than the mean velocity of the packet. The mean longitudinal velocity and temperature of the undecelerated beam is estimated from the simulations to be 480 m s −1 and 2.0 K respectively. Any minor discrepancies between the simulated and experimental ToF profiles can be attributed to differences between the simulated ideal gas pulse and the experimental beam.
The simulation analysis enables one to disentangle the contribution that decelerated H atoms make to the final ToF distribution. This confirms that state selectivity is attained by the decelerated particles, and that the final desired velocity is achieved. The analysis is illustrated for the case of the evolutionary algorithm-generated sequence for a final velocity of (200±10)m s −1 . The panels on the left-hand side column of figure 4(b) plot the final velocity of the simulated particles against their ToF-shown for the undecelerated (first row) and decelerated (second row) beams, with the contribution from LFS (third row) and HFS (fourth row) particles explicitly shown. From the 'Dec (LFS)' plot, it can be seen that the majority of the particles in the decelerated peak have a ToF arrival time of approximately 860 μs, and that this packet of decelerated LFS particles travels with a final velocity centred on 200 m s −1 . A minor contribution to the decelerated peak arises from the small number of particles (in both LFS and HFS states) with an initial velocity around 300 m s −1 that are not decelerated by the magnetic fields. Thus it should be noted that-while all particles travelling at approximately 200 m s −1 are in LFS quantum states-not all particles forming part of the 'decelerated' peak in the ToF trace have been decelerated, as there is a minor contribution from initially slow-moving particles in the incoming beam. The panels on the right of figure 4(b) show plots of the final velocity of the simulated particles against their initial velocity. Here, it can be seen that all particles that are decelerated to a final velocity of approximately 200 m s −1 are LFS particles that started with an initial velocity between 400 and 480 m s −1 . Repeating the same analysis for all sequences confirms that the desired final velocity is reached as anticipated from the sequences in all cases.

Discussion
Simulation analysis Evolutionary algorithms can improve on the standard switching sequences to achieve greater deceleration than would otherwise be possible, given the experimental contraints. In order to understand how this is achieved, simulations of the standard pulse sequence yielding a final velocity v 205 m s  figure 3(a), respectively). Figure 5 shows the phase-space plots of the trajectories of particles contributing to the decelerated peak, both before entering the decelerator (left panel) and after all particles in the phase-stable region have exited the final coil (right panel). The standard sequence behaves as expected: as it is based on the trajectory of a single synchronous particle with an initial velocity of 500 m s −1 , particles with a similar velocity and position to the synchronous particle also fall within the phase-stable region. In contrast, the optimised evolutionary algorithm sequence considers an ensemble of particles with a distribution of velocities and positions. It preferentially addresses the slower moving particles in the initial distribution, meaning that less translational energy needs to be removed to achieve the same final velocity. The particles decelerated by the evolutionary algorithm sequence emerge from the decelerator with a similar spread in velocity but a larger spread in position than those decelerated by the standard sequence ( figure 5, right panel). This is further illustrated in figure 6, where it can be seen that the spatial distribution in the x and y dimensions is comparable but that there are differences in the spread along the z axis after particles exit the decelerator. The evolutionary algorithm sequence yields decelerated particles with a 14mm spread of positions along the z axis, whereas the standard sequence distribution along z is only 8mm. Despite the distribution in z extending an additional 6mm, simulations indicate that the average density of decelerated particles is only reduced by approximately 20% in the evolutionary algorithm sequence compared to the standard sequence. While considering the density of particles produced by the decelerator is important, optimising the number of particles produced by a deceleration sequence is far more important for our applications. Ultimately, we seek to combine the decelerator with an ion trap to undertake collision studies between paramagnetic species and trapped ions. Provided the distribution in the xy plane matches the acceptance of the ion trap-which could potentially be achieved through the inclusion of additional focusing components-the spread in z is less important to our intended application than the number of particles that reach the reaction region. To this end, the evolutionary algorithm sequence far exceeds the performance of the standard sequence: the number of Only the particles in the initial beam that are decelerated to velocities within the relevant ranges are shown. Particles decelerated to within±10 m s −1 of the target velocity are indicated by dark blue and dark red dots for the evolutionary algorithm and standard sequences, respectively. A wider contour of particles decelerated to the target velocity±70 m s −1 is also shown for reference, indicated by light blue and light red dots for the evolutionary algorithm and standard sequences, respectively. The dashed lines in the left panel represent the boundaries of the initial phase-space distribution of the undecelerated beam. Note that the dashed lines are not vertical due to the intial expansion of the incoming beam before it reaches the first decelerator coil (with faster particles travelling further than slower particles). decelerated particles is 40% higher for the optimised sequence for a comparable final velocity. The success of the evolutionary algorithm sequence can be attributed to its increased phase-space acceptance and to its ability to address particles that are already travelling slower and therefore require less deceleration. The latter is observed across all of the evolutionary algorithm sequences studied; the lower the target velocity to be optimised for, the lower the initial velocity of the particles that are decelerated. Qualitatively, this can be seen in the progressively delayed switch-on time of the current pulse applied to the first coil in figure 3(a). To achieve lower final velocities using a standard sequence, a synchronous particle with a lower initial velocity could be selected, provided a lower effective phase angle is also employed to ensure that the maximum pulse duration does not exceed 60 μs (for the decelerator used in this work). A synchronous particle with an initial velocity of, for example, 450 m s −1 -comparable to the modal velocity selected by the evolutionary algorithm-optimised sequence yielding the same final velocity-could be decelerated to 200 m s −1 using a standard sequence, with 53 0 f = .
However, our simulations demonstrate that this approach yields significantly fewer decelerated particles. While one might anticipate an increase in the number of decelerated particles transmitted through the decelerator when the phase angle is reduced, this effect is outweighed by the decrease in the number of particles in the initial beam that have a low enough velocity to be addressed by the sequence. Indeed, for the same final velocity, the standard sequence with = -) and 65% fewer particles than the evolutionary algorithm-optimised sequence. Thus employing standard sequences with a lower phase angle and lower initial velocity synchronous particle results in significantly fewer decelerated particles. Such sequences cannot be designed to reach lower final velocities than standard sequences with 90 0 f = .
The velocity profile as a function of time inside the decelerator also looks distinctly different for the standard and evolutionary algorithm sequences (see figure 7). For the standard sequence, the packet of decelerated atoms is kept together in phase space and experiences a periodic magnetic field pattern and thus a stepwise decrease in mean velocity at each stage. For the optimised sequence, on the other hand, the deceleration appears to occur in two stages: in the first stage (coils 1-5), particles that are initially moving faster are targeted, whilst the velocities of the slower particles are barely affected by the magnetic fields at all. In this first stage, longitudinal focusing (in the z direction) is accompanied by the defocusing of v z in addition to transverse defocusing (defocusing in the xy plane). This results in the spread of velocities that can be observed just before coil 6 is switched off (sixth vertical yellow line). In the second stage, all particles in the phase-stable bunch are decelerated by the remaining coils and v z gets re-focused, as shown by the increased density of particles around v 200 m s z 1 = towards the end of the deceleration process. In order to keep the entire packet together and refocus the velocity spread, the initially slower particles (which are now the faster particles in the packet, owing to the deceleration of the previously faster particles) are preferentially decelerated at a higher rate. The velocities of all particles trend downwards for each successive coil from coil 6 to the end of the decelerator. While this description of the two-stage deceleration is an oversimplification-explicitly considering only the extremes of the initial velocity distribution-it helps to illustrate the differences in how the evolutionary algorithm and standard sequences achieve deceleration of the phase-stable packet, and the way that the evolutionary algorithm is able to capture a broader range of the initial phase space distribution. The two-stage deceleration process adopted by the evolutionary algorithm sequence can also be observed in the distribution of phase angles exhibited by all particles in the decelerated packet at the switch-off time of each coil (see figure 8)-corresponding to the positions of particles in the coil at the time when the magnetic field is switched off. The range of phase angles in particles addressed by the standard sequence (left column) are fairly uniform and narrow for all coils, with the peak in the distribuion a little below 90 f =  for the earlier coils and slowly shifting to lower phase angles for progressively later coils. In the evolutionary algorithm sequence (right column) there is a clear discontinuity between coils 5 and 6, the point at which the first stage of deceleration (outlined above) ends and the second stage begins. In the first stage, corresponding to coils 1-5, the distributions peak before 90 f =  and the peak shifts to lower phase angles as seen in the standard sequence-although the decrease of the peak phase angle occurs at a higher rate. In contrast to the standard sequence, the distributions are much broader due to the slower particles in the bunch that are lagging behind the faster ones. In the second stage, from coil 6, an appreciable number of particles have travelled further than the middle of the coil ( 90 f = ) by the time the coil is switched off and, as a result, they get decelerated and then slightly accelerated again. As close inspection of figure 7 reveals, some of the initially fastest moving particles are periodically decelerated and then accelerated as they pass through coils 6-12-although there is still a net deceleration of all particles as they pass through each of the coils in the second stage. This unintuitive approach serves to bunch the velocity distribution, returning the desired narrow distribution of velocities at the exit of the decelerator.
The simulated temporal evolution of the decelerated bunch of particles is provided as supplementary material available online at stacks.iop.org/NJP/19/083016/mmedia. These supplementary material videos show the change in the phase-space distribution of particles addressed by the standard and evolutionary algorithm sequences as they move through the decelerator coils, giving a more complete picture of the different deceleration mechanisms.
Performance of the evolutionary algorithm sequences Optimised evolutionary algorithm sequences outperform standard sequences due to two key factors. Firstly, the optimised sequences can circumvent the issues faced by the standard sequences by targeting particles that are already travelling slower without needing pulses that are too long to be implemented experimentally. Given that 500 m s −1 is the most probable velocity of particles in the beam entering the decelerator, addressing particles with v 500 m s z initial 1 < means that fewer particles are available for deceleration. This trade-off between addressing already slower particles and maintaining the highest possible number of decelerated particles is automatically optimised by the algorithm. Secondly, the evolutionary sequences can target a larger spread of initial velocities since the optimisation is performed by considering a large number of particles with a range of initial positions and velocities instead of one synchronous particle, with focusing effects automatically accounted for. As a result, the evolutionary algorithm sequences have a larger phase-space acceptance and achieve significantly enhanced numbers of decelerated particles than standard sequences, for a comparable final velocity. The beauty of using evolutionary algorithms is that one can tweak the criteria against which the pulse sequence parameters are optimised. Should a narrower distribution in phase space be required, one could choose to optimise the sequences both in terms of final velocity and spread of position by adding an additional criterion to the fitness function. This is a key difference between the optimised evolutionary algorithm sequences and the standard sequences. Typically, deceleration sequences can be designed to rotate the phase space distribution to minimise the spread in z, or to compact the velocity distribution at the expense of the distribution in z; it is difficult to devise standard sequences that minimise the spread in both velocity and space. Through the use of evolutionary algorithms, we can select a larger part of the phase space distribution of the incoming beam in order to optimise both parameters simultaneously. The two-stage deceleration approach seen with the evolutionary algorithm sequence allows species to enter the phase-stable region as they move through the apparatus-increasing the density of the packet of decelerated particles throughout the second stage of deceleration such that the number of species decelerated to within the target velocity range is maximised at the end of the decelerator.
The inclusion of additional criteria was not explored in this implementation, for the dual reasons of wanting to keep the optimisation as simple as possible and because a narrow spatial distribution (at least along the z axis) is not necessary for the intended applications of the decelerator described in this work. Additional focusing elements positioned after the decelerator-which will shortly be incorporated into the apparatus-will provide one with the ability to focus the decelerated packet of particles for reaction studies or other applications. Indeed, there are no limits in principle to the criteria against which the parameters can be optimised using evolutionary algorithms; provided that at least one particle meets the criteria during the first generation of sequences, the evolutionary strategy approach will return the best sequence for a given set of restraints. Computational time could potentially become a limiting factor for the optimisation of very complicated systems. The optimisation procedure described in this work takes on the order of 1 hour per sequence, compared to seconds to generate a standard sequence, performed on a quad-core processor.

Conclusion
Optimised sequences of current pulses applied to a twelve-coil Zeeman decelerator have been generated using a CMA-ES algorithm, for a range of target final velocities. Experimental measurements confirm that sequences optimised using an evolutionary algorithm yield both higher numbers of decelerated particles (for the same final velocity) and enable lower final velocities to be reached than is possible with standard sequences. Specifically, optimised sequences achieve a 40% increase in the number of decelerated particles for a final velocity of 200 m s −1 , and they can produce particles with a modal kinetic energy that is 98% lower than that of the incoming beam-removing significantly more kinetic energy than can be achieved with standard sequences. These results are obtained using the same incoming beam and by imposing the same set of experimental restrictions for coil current (243A), repetition rate (10 Hz) and maximum pulse duration (60 μs). The experimental ToF traces are well reproduced by three-dimensional particle trajectory simulations, which reveal some surprising differences in how deceleration is achieved in the standard and optimised approaches.
Standard sequences consider the passage of a single synchronous particle travelling with (typically) the modal initial velocity of the incoming beam. In contrast, optimised sequences target particles that are travelling slower than the synchronous particle-and therefore require less deceleration to reach the same final velocitywhile at the same time compensating for the loss of intensity in the beam by addressing a wider range of intial velocities. By selecting progressively slower portions of the initial distribution, optimised sequences can decelerate particles to lower final velocities without exceeding the maximum pulse duration. This successful deceleration of particles with a range of initial velocities is achieved by a non-trivial deceleration mechanism that could not have been conceived using a standard optimisation method.
In conclusion, using evolutionary algorithms to optimise deceleration sequences leads to significant increases in both the number of decelerated particles and the amount of deceleration achievable with a 12-stage Zeeman decelerator. The approach described here is highly versatile, robust and easily implemented. Employing a CMA-ES optimisation method could greatly improve the efficiency and operation of similar devices, potentially circumventing the need for longer decelerators in order to access lower final velocities. The method could be used to optimise the sequence for a wide range of desired outputs by modifying the optimisation criteria (i.e. fitness function) of the algorithm. For example, the number of decelerated particles arriving at a specific position within a given velocity range could be maximised whilst the number of particles with velocities outside the target range at the desired position is simultaneously minimised-achieving both deceleration and spatial purification of the beam at a specific post-decelerator position. Optimisation using evolutionary algorithms is extremely powerful in generating superior deceleration sequences with improved abilities to decelerate more particles and reach lower velocities. The use of such evolutionary strategies could lead to significant advances in controlling the properties of cold beams formed by the manipulation of magnetic or electric fields.