Nonlinear dynamics in an alternating gradient guide for neutral particles

Neutral particles can be guided and focussed using electric field gradients that focus in one transverse direction and defocus in the other, alternating between the two directions. Such a guide is suitable for transporting particles that are attracted to strong electric fields, which cannot be guided using static fields. Particles are only transmitted if their initial positions and transverse speeds lie within the guide's phase space acceptance. Nonlinear forces are always present in the guide and can severely reduce this acceptance. We consider the effects of the two most important nonlinear forces, a term in the force that is cubic in the off-axis displacement, and a nonlinear term which couples together the two transverse motions. We use approximate analytical techniques, along with numerical methods, to calculate the influence of these nonlinear forces on the particle trajectories and on the phase space acceptance. The cubic term alters the focussing and defocussing powers, leading either to an increase or a decrease of the acceptance depending on its sign. We find an approximate analytical result for the phase space acceptance including this cubic term. Using a perturbation method we show how the coupling term leads to slow changes in the amplitudes of the transverse oscillations. This term reduces the acceptance when it reduces the focussing power, but has little influence when it increases that power. It is not possible to eliminate both nonlinear terms, but one can be made small at the expense of the other. We show how to choose the guide parameters so that the acceptance is optimized.


Introduction
A neutral particle in an electric field gradient feels a force. The field polarizes the particle and, if it has a gradient, pulls on the induced dipole moment. In its ground state, a neutral particle is always attracted to regions of strong electric field. In other quantum states, the attraction may either be to strong fields, or to weak fields. Here, we are concerned with the guiding and focussing of strong-field-seeking particles. Static electric fields cannot do the job and so a dynamic focussing scheme needs to be used [1,2]. An alternating gradient (AG) focusser is a series of lenses, each focussing in one transverse direction and defocussing in the other. The focus and defocus directions alternate. Particles have stable trajectories through the focusser because they are close to the axis when they travel through the defocussing lenses, and further away when passing through the focussing ones.
The purpose of an alternating gradient guide is to transmit a large number of particles over arbitrarily long distances with as little loss as possible. Such a guide might employ an electrode structure that is uniform along the beamline, with the alternation achieved by switching the applied voltages -alternation in time. Alternatively, the voltages can be fixed with the alternating gradients set up by the electrode structure itself -alternation in space. Early experiments with alternating gradient structures used two-rod and four-rod geometries to focus polar molecules, e.g. [3,4]. These same ideas were later applied to cold atoms, demonstrating that an atomic beam or fountain could be focussed this way [5,6]. Of considerable interest at present is the application of AG focussing to the Stark deceleration of strong-field seeking molecules. Prototype alternating gradient decelerators have been demonstrated experimentally [7,8,9], and their focussing properties investigated in detail both theoretically and experimentally [10]. For heavy molecules to be decelerated to rest, these decelerators will need to employ very many stages of AG focussing, typically 100 lenses or more. Recently, a long AG guide with a double bend has been used to guide slow molecules from an effusive source [11]. As buffer gas cooling technology advances, it is likely that such guides will be more widely used as efficient tools for extracting molecules from the cold buffer gas [12,13]. An AG guide has also been used to separate conformers of C 6 H 7 NO by exploiting the selectivity of the guide to the ratio of dipole moment to mass [14]. The theoretical analysis we present here is suitable for describing all of these linear guiding structures. It can also be applied directly to other alternating gradient structures such as a storage ring for strong-field seeking molecules [15] or linear ac traps of the kind demonstrated for both molecules [16] and atoms [17]. The physics we present also carries over directly to ac traps with cylindrical symmetry because the same terms appear in the equations of motion [18].
The anharmonic forces present in AG guides and traps for neutral particles tend to be detrimental to the phase space acceptance, often reducing this by a factor of 10 or more. Similarly, these nonlinearities reduce the quality of the images in an AG focussing system. This problem has long been recognized. In reference [4] for example, molecules were guided in a four rod geometry with the rod positions carefully chosen to optimize the linearity of the forces, resulting in a large improvement in transmission over previous two rod geometries. The optimization of electrode geometry with respect to minimizing nonlinear forces has been discussed in [19] and [10]. Numerical simulations have been used to predict the performance of several types of AG guides [11], decelerators [10] and traps [18,16] all confirming the hugely deleterious effects of the nonlinearities.
In this paper we use a mixture of analytical and numerical techniques to calculate the effect of the nonlinear forces on the particle trajectories and the phase space acceptance of an AG guide. We start by reviewing the linear theory since this is the foundation for the analysis of nonlinearities. It is well known that the motion consists of an oscillation driven at the frequency of alternation, called the 'micromotion', superposed on a slower oscillation of larger amplitude called the 'macromotion'. A good approximate description of the macromotion is obtained using an effective potential. The anharmonic terms modify this effective potential and so alter the phase space acceptance. The leading anharmonic terms in this potential are a quartic term along with a term that couples together motion in the two transverse directions. The quartic term in the potential leads to terms in the force that are cubic in the transverse coordinates. We treat this quartic term analytically, deriving an expression for the phase-space acceptance in the presence of this term. To understand the effect of the transverse coupling term we apply a perturbation method which reveals the nature of the coupled trajectories and shows that this term tends to decrease the phase space acceptance even when it increases the depth of the effective potential. We use numerical techniques throughout to support the theory, test the approximations and obtain results when the perturbation approach ceases to be valid. We apply our analysis to some example AG guides for neutrals.

Equations of motion
Our starting point is equation (8) of reference [10], which gives a suitable form for the electrostatic potential inside an AG guide, Φ(x, y), as a function of the transverse coordinates x and y: Φ(x, y) = Φ 0 a 1 x r 0 + a 3 x 3 − 3xy 2 3r 3 0 + a 5 x 5 − 10x 3 y 2 + 5xy 4 5r 5 0 . ( This equation is derived by writing down a multipole expansion of the potential, and then eliminating terms that are unsuitable for a guide. The scale factors, r 0 and Φ 0 characterize the size of the electrode structure and the applied voltages, while a 1 , a 3 and a 5 are dimensionless constants representing the sizes of the dipole, hexapole and decapole terms in the expansion. A neutral particle inside the guide will be polarized by the electric field to a degree that depends on the particle polarizability and on the magnitude of the field. The interaction of this induced polarization with the electric field changes the potential energy of the particle; this is the Stark shift, W . If the electric field is spatially inhomogeneous, so too is W , and there is an associated force acting on the particle, F = −∇W . For our present purpose, we do not need to know the detailed dependence of the Stark shift on |E|, the magnitude of the electric field. We need only know that, for sufficiently small electric fields, the Stark shift is quadratic in |E|, W = − 1 2 α|E| 2 , while for sufficiently large electric fields it becomes linear, W = −µ|E|. The constants of proportionality in the two cases are known as the polarizabilty, α, and the dipole moment, µ. Whether the electric field is considered small or large depends on the polarizability of the particle. For example, the Stark shift of a CaF molecule is in the linear regime at a field of 100 kV/cm, whereas the Stark shift of a Li atom is quadratic in the same field.
We now impose the additional constraint that a 5 ≪ a 3 ≪ a 1 , meaning that the electric field in the guide is dominated by a constant term, and that the decapole term is much smaller than the hexapole term. We then obtain the following approximate expression for the n'th power of the electric field magnitude, valid throughout the region where x, y < r 0 : In this equation, we have set the redundant scale parameter a 1 to unity and have introduced scaled transverse coordinates X = x/r 0 , Y = y/r 0 . E 0 is the electric field at the centre of the guide. The neglected terms are of order a 3 a 5 , a 3 3 and higher. Since n can be any number, we see from (2) that the potential in which the particles move always has the same functional dependence on the transverse coordinates, irrespective of whether the Stark shift is linear, quadratic, or intermediate. Although we have treated the case of a long guide, similar terms appear in the treatment of cylindrically symmetric traps [18], and so our analysis could also be readily applied to that case.
In writing the equations of motion, we may either use the axial coordinate, z, or the time, t, as the independent variable. The former is most appropriate for long guiding structures, while the latter tends to be used for trapping geometries. We will use z throughout, writing the equations of motion as where the primes denote differentiation with respect to the scaled axial coordinate Z = z/r 0 , and η, χ and ζ are dimensionless parameters that depend on the character of the Stark shift, but not on the spatial coordinates. For a linear Stark shift, χ = 2|a 3 | − 6a 5 /|a 3 |, ζ = 2a 5 /|a 3 |, and For a quadratic Stark shift, χ = |a 3 | − 6a 5 /|a 3 |, ζ = |a 3 | + 2a 5 /|a 3 |, and We have introduced the quantity R, which characterizes the Stark shift on the axis of the guide in terms of the forward kinetic energy, and will be useful later in the paper. We have written the equations of motion for a lens that focusses in the y-direction when a 3 is positive. For focussing in the x-direction, simply interchange X and Y . We note from (3) that the dominant term in the force is linear in the transverse coordinate, focussing in one plane and defocussing in the orthogonal plane. In addition to the desirable linear term, there is a term cubic in the transverse displacement, and also a nonlinear transverse coupling term. By a suitable choice of the parameters a 5 and a 3 , it is possible to eliminate one or other of the nonlinear terms, but simultaneous elimination of both nonlinear terms is impossible without also setting η to zero. For example, when the Stark shift is linear, the cubic term in the force can (in principle) be eliminated by choosing a 5 = 0; alternatively, the coupling term can be removed by choosing a 5 = a 2 3 /3; the elimination of both terms is clearly not possible. Thus, the lenses in an alternating gradient guide for neutral particles are necessarily aberrant. More generally, the values of ζ and χ are linearly related via χ + 3ζ = 2n|a 3 |. This relationship shows that it is also impossible for both χ and ζ to be negative, implying, as we will show later, that no matter how the electrodes are arranged at least one of the terms is detrimental to the transmission of the guide.
In (3), we are considering the case where a 5 ≪ a 3 ≪ 1, and so the coefficient of the linear term is much larger than the coefficients of the nonlinear terms. It is therefore tempting to suppose that the effect of the nonlinearity on the trajectories is small. However, in going from a focussing lens to a defocussing lens, the linear term changes sign (as it must) but the nonlinear terms do not. The alternating gradient guide can support a stable beam because the envelope of the beam is always smaller in the defocussing lenses than in the focussing lenses. For this reason, the net effect of the linear forces is to confine the molecules, but this trajectory-averaged confining force is considerably weaker than the linear forces of the individual lenses. By contrast, the nonlinear forces are identical in the two lenses and simply add up. It follows that, even though the nonlinear terms in the individual lenses are small, their role is crucial in understanding the dynamics of the guided particles.

Motion in the linear approximation
Throughout this paper we consider a linear AG guide consisting of alternate converging and diverging lenses each of length L, separated by gaps of length S. To understand the effect of the nonlinear forces, we first analyze the motion in the absence of the nonlinearities. This idealized motion has been discussed in some detail in reference [10] and is based on the original work of Courant and Snyder [20] who discussed the same problem in the context of charged particle focussing. We give a brief summary of this theory here, and extend the discussion of reference [10].

Trajectories
Neglecting the nonlinear terms, we write the equation of motion in one dimension as Here, K(Z) = η 2 G(Z), where G(Z) is a function whose value is +1 in a focussing lens, −1 in a defocussing lens, and 0 in a drift space. We know the solution in each region of constant G. For example, the solution inside a focussing lens is given by The solution inside a defocussing lens is identical to (7), with the replacement of η by iη, and the transfer matrix denoted by D(η, Z). For a region of free space, let η ⇒ 0 in (7) and call the matrix O(Z).
Knowing the solutions in each of the separate regions, we can construct the full trajectories piecewise, but the Courant-Snyder formalism offers a far better approach. This involves looking for a general solution of the form where β is a Z-dependent amplitude function that has the same periodicity as the AG array, ψ is a Z-dependent phase, and ǫ, δ, A 1 and A 2 are defined by the initial conditions. By substitution we find that (8) is a valid solution to (6) provided and It appears that all we have achieved is the replacement of the differential equation for X(Z) by a more complicated differential equation for β(Z). The advantage is that we will never have to solve this new differential equation. We already know the solution to (6), but piecewise, and our aim is to express this in a convenient form. From (8) we have where Using (8) and (11) we find the relationship between the coordinates X, X ′ at position Z and those at position Z + L cell , L cell being the periodicity of the AG array. We make use of the periodicity constraint on β, β(Z + L cell ) = β(Z) (and hence also on α) and thus obtain: where γ = (1 + α 2 )/β and is the phase advance per unit cell. Since the integral is taken over a full period, µ is independent of Z. The matrix appearing in (13) is known as the Courant-Snyder matrix, M. The procedure to obtain β(Z) is now straightforward. We illustrate this procedure for a simple guide with no gaps. Consider an alternating sequence of focussing and defocussing lenses, each of length L, with Z = 0 defined to be at the entrance of a focussing lens. We first obtain an explicit expression for the transfer matrix between a point Z inside a focussing lens and the equivalent point one lattice unit further downstream. This matrix is M = F (η, L − Z).D(η, L).F (η, Z). We then calculate 1 2 Tr(M ), which, from the Courant-Snyder matrix is equal to cos µ. The explicit result is cos µ = cos(ηL) cosh(ηL). Finally, we calculate β by equating the upper right hand element of M to that of M, obtaining for 0 ≤ Z ≤ L. We can apply the same procedure to a point Z inside a defocussing lens (i.e. L < Z < 2L). The result for β(Z) is identical to (15) with the replacement of η by iη and the replacement of Z by Z − L. Since β has period 2L, we now know it everywhere, and since ψ(Z) can be calculated directly from β(Z) we have the complete general solution to the equation of motion in the form of (8). This same procedure can equally well be applied to more complicated guiding structures.

Phase-space acceptance
Returning to the case of a general guide, particles enter with a range of transverse positions and angles and we wish to know which of them are transmitted. Consider a particle whose initial conditions are given by ǫ and δ, as in (8). Using the first line of (8) to form the quantity X 2 + (αX + βX ′ ) 2 , and recalling that β ′ = −2α, ψ ′ = 1/β and 1 + α 2 = γβ, we find the invariant Suppose we record the transverse position and angle of the particle at the longitudinal positions Z + N L cell , for a large number of integer values of N . Then a plot of these points on a phase-space diagram, whose axes are X and X ′ , traces out the ellipse defined by (16). If we repeat this procedure for a different value of Z, the shape of the ellipse evolves periodically with Z according to (16), but its area remains the same, πǫ. A distribution of particles having ǫ i ≤ ǫ fills the ellipse. If the guide has a constant transverse aperture with walls at X = ±1, then it is clear from (8) that all particles having ǫ i < ǫ will be transmitted indefinitely provided that √ β max ǫ < 1, where β max is the maximum value of β(Z). If this condition is satisfied, the transverse extent of the beam never exceeds the transverse aperture of the guide. Particles that have ǫ i > 1/β max may also be transmitted by the guide provided that the turning points of the macromotion (the cosine factor in (8)) never coincide with the turning points of the micromotion (the amplitude function in (8)). That condition may occur for a large number of particles if µ/π is a rational number, µ being the phase advance per unit cell defined by (14). However, these "resonances" in the phase-space acceptance become very narrowly centred on these special values of µ once the guide is long compared to the wavelength of the macromotion, and we can usually ignore them. Then, the transverse phase-space acceptance is simply π/β max . As discussed in [10], the maximum value of β always occurs in the centre of a focussing lens. Returning to our specific example of a gapless guide, this is explicitly seen in (15) at Z = L/2. It follows that, in this case, the transverse phase space-space acceptance in one dimension is Since we are using Z as the independent variable, the acceptance has the units of XX ′ , i.e. r 0 · radians. There are, of course, two transverse directions, and the 2D transverse acceptance is simply the square of the above formula. So far, we have supposed that the transverse aperture formed by the electrodes of the guide is uniform along its entire length. For many electrode geometries, this is a poor approximation. An interesting case is that of alternating two-rod lenses, discussed in [10] and used in all experimental work on alternating gradient deceleration reported to date. Here, the electrodes lie in the defocussing plane and there is no aperture at all in the focussing plane. The maximum beam envelope still occurs at the centre of the focussing lens, but there is no aperture at this point for the beam to crash into. Instead, the appropriate limiting value of β to insert into the acceptance equation is its value at the entrance to the defocussing lens. As a result, the phasespace acceptance increases. For the case of a gapless guide with an aperture at the defocussing lenses only, the second term in the denominator of (17) is replaced by cos(ηL) sinh(ηL). The acceptance in one transverse direction of a gapless guide is plotted in figure 1, for both kinds of transverse aperture. We see that when the aperture is defined only by the defocussing lenses (dashed line) the 1D acceptance peaks at a slightly higher value of ηL where it is almost twice as large. Later in this paper, we will make a "small micromotion approximation" which amounts to neglecting the modulation of β, replacing it by the constant value L cell /µ. In this approximation, the phase space acceptance in one dimension becomes πµ/L cell . For a gapless guide, this approximate acceptance is shown by the dotted line in figure 1. We see that the approximate result agrees well with the exact results when ηL ≪ 1, because the micromotion amplitude is indeed small in this regime. As ηL increases, the approximate result overestimates the acceptance. At ηL = 1 for example, the approximate acceptance is 31% larger than the exact value with uniform aperture, and 2% larger for an aperture at the defocussing lenses only.
We turn now to some specific examples. Consider first a beam of ground state CaF molecules travelling at 350 m/s into a gapless guide having E 0 = 100 kV/cm and a 3 = 0.2. In this case, η = 0.053, and with r 0 = 1 mm the acceptance reaches its maximum value of 39 mm mrad when the lenses are 24 mm in length. Thus, the acceptance of the guide, which can transport the molecules over an arbitrarily long distance, is the same of that of a 1 mm 2 square cross-section pipe of length 26 mm. As a second example, consider Li atoms launched with a speed of 10 m/s from a magnetooptical trap (MOT) into the same guide. Here, η = 0.31, the optimal lens length is 4 mm and the transverse acceptance is 230 mm.mrad, equivalent to 2.3 mm.m/s in position-velocity space. This acceptance is larger than the phase-space area occupied by the atoms in a typical Li MOT, indicating that all the atoms can be successfully transported.

Fill factor
The previous section showed how to calculate the transverse phase-space acceptance, but often the source is unable to fill this acceptance area completely. The final output of particles depends on the overlap between the acceptance and the phase-space area occupied by the particles at the guide's entrance, which we call the fill factor. As an illustration, consider a beam formed in a supersonic expansion, passing through a small aperture (a skimmer) and then into the guide. For the purposes of calculation it is convenient to proceed in the opposite direction. We calculate the phase-space acceptance at the entrance of the guide, project this back through the skimmer and onto the source to find out what fraction is filled by the source. In practice, the source is usually large enough and divergent enough to completely fill the part of the acceptance that fits through the skimmer, and we assume this to be the case in the calculation that follows. Let us write the phase-space ellipse accepted by the guide as Here, A 1 = β max γ(0), B 1 = β max α(0) and C 1 = β max β(0), where Z = 0 is defined to be at the entrance of the guide. We project this ellipse through a distance −D and so obtain a new ellipse with coefficients A 2 , B 2 , C 2 . The relationship between the This is easily derived using the drift matrix O(Z) (see (7)). At this point, we can calculate the area of that part of the phase-space ellipse that fits through the skimmer aperture. The aperture provides a constraint in the X direction, but not in the X ′ direction. Integrating the ellipse between X = ±R, we obtain the area ‡, Repeating this for the other transverse direction, and multiplying together the two results yields the 2D transverse phase-space acceptance of the entire guide and aperture setup. This result, divided by the 2D acceptance of the guide alone, (π/β max ) 2 , gives the fill factor. Figure 2 shows an example of how the fill factor decreases as D increases, stressing the importance of minimizing the distance between the guide and any aperture placed upstream.

An effective potential approach to the nonlinear dynamics
In order to solve the nonlinear problem, whose solution is X(Z), let us define new coordinates, w and s, in terms of the function β(z) which appears in the solution to the linear problem, equation (8): The phase advance s increases by 1 each time the particle advances through one unit cell. The quantity w represents the macromotion of the particle. Applying the transformation to (6) we obtain (22) where the primes still denote differentiation with respect to Z. Note that we have made no assumptions about the form of X(Z) in this transformation. With the help of (10) and (12) we see that 1/2 β β ′′ − α 2 + Kβ 2 = 1, and so our equation of motion in the transformed coordinates reduces to that of a harmonic oscillator of angular frequency µ,ẅ where the dots denote differentiation with respect to s. In the transformed coordinates the micromotion has disappeared and the macromotion appears as motion in a harmonic potential whose curvature is proportional to µ 2 .
Re-introducing the nonlinear terms, we write (3) in the form where we use the function G(Z) introduced earlier so as to handle an arbitrary sequence of converging lenses, diverging lenses and drift spaces. Applying the same coordinate transformation as above, this becomes In (25a), w 1 , β 1 and s 1 are identical to the w, β and s introduced earlier, and the dot indicates differentiation with respect to s 1 . The quantities w 2 , β 2 and s 2 in (25b) are the corresponding quantities for the y motion, and differentiation is with respect to s 2 . Naturally, β 1 and β 2 are related: So far, our analysis has been exact. In order to make progress, let us assume that the guide is operated well within the range of stability. In this regime the amplitude of the micromotion is small compared to the macromotion amplitude and the envelope function, β(Z), is approximately constant. We set β 1 = β 2 = β 0 . It follows from (14) and (21) that β 0 = L cell /µ and s 1 = s 2 = Z/L cell . In addition, the wavelength of the macromotion is large compared to the length of a unit cell in this regime, and so G 2 oscillates rapidly relative to the macromotion wavelength. We will replace it by its mean value G 2 . In this way, the same formalism will apply to guides with or without gaps. For many practical guides these approximations will turn out to be useful and not very restrictive.
Making these replacements in (25a) and (25b), and then transforming back to the more familiar variables, X and Y , we obtain where differentiation is with respect to s, interchange of X and Y yields the second of the two equations of motion, and Patience will reward the reader's curiosity about the seemingly superfluous little epsilon.
In the small micromotion approximation, we can think of the particles moving in an effective potential, V eff . From (26), and the equation obtained by interchange of X and Y , we find this effective potential to be given by We also note that the transverse phase space acceptance obtained by applying the small micromotion approximation in the harmonic case is simply πµ. Converting from our present coordinates, (X, dX ds ), back to (x, dx dz ), this acceptance becomes πµ η r 0 /(ηL cell ). This result was already discussed below (17) in the context of figure 1.

The cubic term
In this section we derive analytical expressions for the trajectories of the particles and for the phase space acceptance of the guide, with the cubic term included. We set b 2 = 0 in (26), multiply by X ′ , and then integrate with respect to s, obtaining where h is a constant of the motion (the total energy), and we identify the first term as the kinetic energy and the second and third terms as the potential energy. Rearranging for X ′ and integrating gives where p = µ 2 /2h and q = ǫb 1 /4h. To complete the solution, we factorize the expression under the square root to the form Then, the integral takes the standard form for an elliptic integral of the first kind, F , and the solution is § Finally, this expression can be inverted to give X(s) in terms of a Jacobi elliptic function: This equation provides an analytical form for the trajectory of a particle through the guide, in the presence of the cubic nonlinearity, but averaged over the micromotion. Figure 3(a) shows some examples of how the trajectories depend on the size of the cubic term, calculated using (33). The figure shows that as ǫb 1 increases, the wavelength increases and the trajectories evolve from a sinusoidal shape towards a more square shape. These changes reflect the way the effective potential changes with ǫb 1 . Part (b) of the figure shows the effective potential along the X-axis for the same parameters used to obtain the trajectories plotted in (a). As the value of ǫb 1 increases, the steepness of the effective potential well decreases, and so the wavelength of the macromotion increases. In the extreme case where ǫb 1 = µ 2 , the particle will stop at the maximum of the potential and remain there forever. From (31) and the discussion that follows it, we find the wavelength of the macromotion to be X max being the turning point of the motion where X ′ = 0, found to be X max = 1/ √ a + . The wavelength evaluates to   where K is the complete elliptic integral of the first kind. It is illuminating to write down the series expansion for Λ, or its inverse, the spatial frequency ν, in powers of the small quantity ǫb 1 : For positive ǫb 1 , the spatial frequency decreases with increasing ǫb 1 and with increasing oscillation amplitude. To lowest order, the change depends on the square of the oscillation amplitude. Particles that have a large amplitude will oscillate more slowly through the guide than those that remain close the guide's axis. Figure 3(c) and (d) illustrate the exact spatial frequency versus X max and ǫb 1 . These plots shows the quadratic dependence on the oscillation amplitude and the linear dependence on ǫb 1 . We now calculate the phase-space acceptance of the guide in the presence of the cubic term. The two transverse directions are uncoupled at this point in our discussion, so it is convenient to calculate the acceptance in 1D, squaring the result to get the total 2D acceptance. When ǫb 1 is positive, the effective potential in which the particles move has a turning point at Particles with enough energy to reach this turning point cannot be confined in the guide. A particle that has just enough energy to reach the top of the effective potential well has an energy h c = 1/2µ 2 X 2 c − 1/4ǫb 1 X 4 c = µ 4 /(4ǫb 1 ). To determine the phasespace acceptance we need to distinguish between two cases. In the first case, the turning point is 'inside' the physical aperture and so the acceptance is limited by the turning point. In the second case, the turning point is 'outside' the physical aperture, or there is no turning point at all (because ǫb 1 < 0), and the acceptance is limited by the physical aperture. In the first case, the equation defining the area in phase space where particles stay inside the potential is The enclosed area is the phase space acceptance and is easily calculated by rearranging for X ′ and integrating under the curve from X = 0 to X = X c where the curve crosses the X-axis. The result is where we have used equation (37). The phase-space acceptance scales as the cube of the phase-advance per unit cell, a measure of the trapping strength, and is inversely proportional to ǫb 1 , the coefficient of the nonlinearity.
In the second case, where the acceptance is limited by the physical aperture at X = X A , the equation defining the accepted phase-space area is By integrating the expression for X ′ with respect to X from X = 0 to X = X A , we find the area enclosed by this curve to be where κ = X 2 A /(2X 2 c − X 2 A ), and K and E are the complete elliptic integrals of the first and second kinds. We note that, for X A = 1, this equation reduces to A = πµ in the limit where ζ → 0. This is the result we already found when applying the constant beta approximation in the harmonic limit. We find it useful to factor out the πµ, writing our result as the product of the acceptance in the harmonic limit, A 0 , and a factor, F ζ , that accounts for the cubic nonlinearity: where when X c ≤ X 0 , and when X c > X 0 . By writing the result in this way, we can choose to use the exact result for A 0 , derived using the procedure detailed in section 3.2, rather than πµ which we already know from figure 1 is not very accurate unless ηL is small. In this way, we effectively use the small micromotion approximation to handle only the nonlinear part of the problem. This approach, which we apply throughout, is justified by the high accuracy of the results it gives when compared to numerical simulations. Together, (42b) and (42c) give the phase-space acceptance for all values of ǫb 1 . The accuracy of this analytical approach depends on the validity of the approximation we have made, namely that β is approximately constant. To explore this, we can calculate the acceptance numerically for some specific cases, and then see how well the numerical and analytical results compare. To calculate the acceptance exactly, we set χ = 0 in (24) and then solve this equation numerically for a large number of particles, taking care to ensure that the initial phase-space distribution is large enough to completely fill the acceptance area. Particles whose trajectories exceed the transverse boundaries of the guide (at X = ±1), are removed from the calculation so that we are left with the set of particles transmitted by the guide. We calculated for a guide consisting of 200 lenses. This was long enough to show that the phasespace acceptance is independent of the guide's length, provided this is longer than the wavelength of the macromotion. Figure 4 shows how the acceptance of a gapless guide depends on the value of ζ. It compares the exact acceptance, calculated numerically, with the analytical result, F ζ , given by (42b) and (42c). Part (a) shows the results obtained with ηL = 0.5 for which the phase advance per cell is µ = 0.14. The numerical results have been normalized to the exact acceptance obtained when ζ = 0, which in this case is 0.426ηr 0 , while the analytical result is already normalized. When ζ is positive, the acceptance decreases rapidly with increasing ζ, as the turning point of the effective potential moves inwards. For negative values of ζ, the aberration increases the curvature of the effective potential and so increases the acceptance. The analytical and numerical results show excellent agreement for all values of ζ calculated. For positive ζ the results agree to within the error bars on the numerical values, which are approximately the size of the points. For negative values of ζ, the analytical result is consistently 2% larger than the numerical results. Part (b) of the figure shows the results obtained with ηL = 1 for which µ = 0.58. Again, the numerical results have been normalized to the aberration-free acceptance, which is 0.701ηr 0 . Here, our small micromotion approximation begins to break down, and the true acceptance is less well represented by the analytical result. Nevertheless, the analytical result differs by 20% at most from the true value. The analytical approximation is very useful since it can be evaluated instantly, unlike the numerical simulation. The full 2D acceptance is simply the square of the 1D acceptance that we've calculated.
We now examine the effect of the cubic term on some example guides. At present, all the experimental work on alternating gradient deceleration of polar molecules has used two-rod lenses. In these decelerators, two rods, 6 mm in diameter, with centres separated by 8 mm, are aligned with their axes parallel to the molecular beamline, and the experiments are performed in the linear Stark shift regime. In this case a 5 /a 3 = a 3 = 1/7 and so ζ = 2/7. Usually, the length, S, of the drift spaces is approximately equal to the length of the lenses, and so we take ηS = ηL. We use our analytical formulae to calculate how the cubic term reduces the 2D phase space acceptance in this geometry, normalizing the results to the values obtained when ζ = 0. When ηL = 0.5, the normalized acceptance, F ζ , is a miserable 3%. This grows to 18% when ηL = 0.75 and to 60% when ηL = 1. We see that the nonlinearity is hugely detrimental if the guide is operated at low values of ηL where the phase advance per cell, and hence the depth of the effective potential, is small. As ηL is increased towards the edge of stability, the nonlinearity becomes less important. This geometry is even worse for the electric guiding of cold atoms whose Stark shift is quadratic, because the quadratic Stark shift gives ζ = 3|a 3 | instead of 2|a 3 |.
In the four rod guide considered in [10], rods of radius 2.3r 0 were placed with centres on the corners of a square of side 6.6r 0 . The rods on one side of the square were at positive high voltage, and those on the other side at negative high voltage, so the arrangement has two-pole symmetry. In this arrangement the cubic term is small when the Stark shift is linear, ζ = 0.028, and the normalized acceptance, F ζ , is 20% when ηL = 0.5 and 79% when ηL = 1. The guide is not so suitable when the Stark shift is quadratic, for then ζ = 0.17. The structure used in [11] to guide slow molecules from an effusive source was built from four rods of radius r 0 with centres on the corners of a square of side 3r 0 . The rods on two diagonally opposite corners were at high voltage, ±V , with the other two grounded. Fitting the electrostatic potential obtained in this geometry to the multipole expansion, (1), we obtain a 3 = 0.59 and a 5 /a 3 = −0.055. So this geometry has ζ = −0.11 for polar molecules with linear Stark shifts, meaning that the cubic term enhances the acceptance over the perfect linear case. However, the large value of a 3 means that our series expansion of the electric field, equation (2), converges very slowly. The expansion does not provide an accurate representation of the field unless terms of higher order are included and so our conclusion about the effect of the cubic term is not so helpful in this case. We will see later that the lowest order transverse coupling term is severely detrimental for this geometry.

Multiple scales analysis
In the preceding section we found analytical solutions for the motion in the presence of the cubic nonlinearity. We can no longer do so when b 2 = 0 in equation (26). Nevertheless, provided the nonlinear terms are small, we can elucidate the motion of the particles using a perturbation technique known as multiple scales analysis. We are particularly concerned with stable motion in the guide, i.e. with those particles that remain inside the guide indefinitely. We recognize that the coupling term will lead to an exchange of energy between the x and y oscillations, and so the oscillation amplitudes will be modulated. If ǫb 2 is small enough, this modulation will occur on a much longer scale than the wavelength of the macromotion. To handle this, we look for a solution in the form of a perturbation expansion using two separate length scales. We follow references [21,22,23], and refer the reader to these books for more details on the multiple scales method.
The independent variable in (26) is s = Z/L cell . We introduce a second length scale, σ = ǫs, treat X and Y as functions of both s and σ, and then expand X and Y as a power series in ǫ: X(s, σ) = X 0 (s, σ) + ǫX 1 (s, σ) + . . .. In order that the initial conditions on X and Y be satisfied uniquely in this expansion, we choose to set X n (0) = Y n (0) = X ′ n (0) = Y ′ n (0) = 0 for all n = 0. The derivatives are Substituting into (26) we obtain a differential equation whose solution should not depend on the arbitrary quantity ǫ. We therefore have a separate equation for each power of ǫ, these being and similarly for Y 0 , Y 1 , by simple interchange of the symbols X and Y . The solutions for X 0 and Y 0 are where c.c stands for the complex conjugate and r 1 , r 2 , θ 1 , θ 2 are real functions of σ, but independent of s. These solutions are now substituted into (44), giving, after a little algebra, Consider for a moment the solution to the above differential equation with the right hand side replaced by zero (the corresponding homogeneous equation). Its solution has terms in e iµs and e −iµs , terms which also appear in the right hand side of the inhomogeneous equation, (46). Terms in an inhomogeneous equation which are themselves solutions to the associated homogeneous equation lead to secular terms in the solution. These secular terms grow with s more rapidly than the corresponding solution of the homogeneous equation by at least a factor of s [23], i.e. as s cos(µs) in this case. If our solution is to represent stable motion, corresponding to particles trapped in the guide indefinitely, we must have to eliminate the secular terms by setting to zero the coefficient, c 1 , of e iµs in (46). We are free to do this, since the r's and θ's are otherwise undetermined functions of σ. Setting c 1 = 0 and then separating out the real and imaginary parts, we obtain the following differential equations describing the slow variation of r 1 , r 2 , θ 1 , θ 2 : By solving this set of differential equations we obtain the slow evolution of the amplitudes and phases, which we can then substitute into equations (45) to get the full solution to lowest order in ǫ. As we will see later, this procedure is remarkably accurate. A great deal of insight into the motion of the particles can be acquired from the above set of equations. By adding together (47a) and (47b) we see that r 2 1 + r 2 2 = E 0 , a constant, proportional to the total energy at this order of approximation. Subtracting (47a) from (47b) and (47c) from (47d) yields where ξ is the fractional energy associated with oscillations in the x direction, ξ = r 2 1 /E 0 , and γ = θ 2 − θ 1 is the phase difference between the two transverse oscillations. The solutions of these equations tell us how the phase difference and the partitioning of the energy evolve with the axial coordinate. Even more instructive is to eliminate the independent variable by dividing (49) by (48), and then integrating to give where C is a constant. This equation is the central result of this analysis. It tells us that there is an invariant quantity involving only ξ and γ. For a given particle, the value of C is fixed by the the initial conditions. The particle is then confined to travel on the contour in (ξ, γ) space defined by the above equation. This tells us how the amplitudes of the two transverse oscillations will change, allowing us to determine whether a particle will remain in the guide or crash into its boundaries.

Small coupling coefficient
In this section, we use the results obtained above to understand the effect of the nonlinear coupling on the motion of particles in the guide. To isolate the effect of the coupling, we set b 1 = 0. To demonstrate that the approximation developed above is a useful one, we begin by looking at an example trajectory. The solid line in Figure 5 shows the exact trajectory of a particle travelling through 120 unit cells of a gapless guide that has ηL = 1, ζ = 0 and χ = 0.03. The trajectory shows the usual superposition of micromotion and macromotion. The transverse coupling is responsible for the slow modulation of amplitude and relative phase evident in the figure. The dashed line in the figure shows the result of the multiple scales approximation. The parameters given above correspond to µ = 0.585 and ǫb 2 = (ηL cell ) 2 χ = 0.12, and the dashed curve is obtained from (45) using the solutions to (47a) -(47d). We see that the approximate solution does a good job of describing the true motion of the particle, including the modulation of the amplitude and phase. In fact, in order that the approximate and exact trajectories could be distinguished in the figure, we had to choose a case where our approximations begin to break down -the micromotion amplitude is not so small relative to the macromotion, and ǫb 2 is not very small compared to µ 2 . In situations where the approximations hold more accurately, the two trajectories are virtually indistinguishable. Now that we have illustrated the usefulness of this approximate approach, we continue with our analysis, returning to the key result given in (50). When b 1 = 0, this equation reduces to representing a set of contours of a form that is independent of b 2 (except as an overall scaling factor). Figure 6 is a plot of these contours along with a colour map The two maxima correspond to motion along the diagonals of the guide. The oscillations in x and y have equal amplitudes and are either exactly in phase (γ = 0), or exactly in anti-phase (γ = π). Even though the coupling is nonlinear, we have found the normal modes familiar from the theory of linearly coupled harmonic oscillators. Indeed, applying the multiple scales method to a pair of linearly coupled harmonic oscillators results in the relationship (1 + cos(2γ))ξ(1 − ξ) = constant, which is remarkably similar to (51). The frequencies of the normal modes can be found by differentiating their phase, φ(s), with respect to s. The phase of the X oscillation at position s is, according to (45), φ(s) = µs + θ 1 (σ), and so the oscillation frequency is ω = dφ/ds = µ + ǫ dθ 1 /dσ. With the help of (47c) the normal mode frequencies are found to be We see that the frequency is lowered when χ (and hence ǫb 2 ) is positive because the coupling term reduces the trap depth along the diagonals of the guide. The opposite is true when χ is negative. Note that the normal modes with γ = 0 and γ = π have equal frequencies, unlike the linear coupling case where the coupling removes the frequency degeneracy. There is another important difference too: the frequency of the normal mode differs from µ by an amount that depends on the square of the oscillation amplitude (the total energy). This is similar to the frequency shift caused by the cubic term discussed earlier, and quite different from a linear coupling whose normal mode frequencies are independent of amplitude. The two maxima in figure 6 are encircled by contours with 1/4 < C/b 2 < 3/4. In this range, the two transverse oscillations exchange energy completely. Looking at the contours near γ = 0, the phase difference is zero when the amplitude of one mode has its maximum value, zero again when it has its minimum value, and reaches a maximum when the amplitudes of the two modes are equal. The phase difference never exceeds π/2. Again, this type of coupled motion is similar to that of linearly coupled oscillators. For a given value of C/b 2 in this range, the maximum value of ξ is found at γ = 0 and is The saddle point in figure 6 is a coupled stationary mode corresponding to circular motion in the xy plane -we call it a circular mode. Its stationary nature is unstable in the sense that an infinitesimally small deviation from this point results in oscillations with time-varying amplitudes and phases. We will show later that the circular modes can become stable when b 1 = 0.
When C/b 2 < 1/4 the contours are below the saddle point and so no longer encircle the maxima. In this range, although the amplitude of one mode does grow at the expense of the other, the amplitudes do not exchange completely. The contours are open and the phase difference increases continually. Of particular interest is the fact that motion with ξ = 0 or ξ = 1 is of this type, and does not yield any exchange of amplitude at all. The motions uncouple when all the energy is in one direction, in stark contrast to the classic motion of linearly coupled oscillators which exchange energy fully even when one oscillator begins at rest. For a given value of C/b 2 < 1/4, the large ξ contour has maxima at γ = 0 and π, whose value, ξ m1 , is given by (53). The small ξ contour has a maximum at γ = π/2 with value Due to the exchange of energy between the transverse modes, there will be some particles whose maximal amplitudes of oscillation take them outside the boundaries of the guide, even though their initial amplitudes are within these boundaries. These particles are lost from the guide in the presence of coupling, and would not have been in its absence. We note again that the value of ǫb 2 dictates the rate at which the amplitudes exchange, but not the degree of exchange. If the guide is infinitely long, and ǫb 2 is small enough that our approximations hold, the loss of acceptance is independent of the size of ǫb 2 . A larger value only results in the particles being lost more quickly, but those same particles will eventually be lost even when the coupling term is infinitesimally small. Similarly, and somewhat counter-intuitively, the same particles are lost whether the nonlinear term increases or decreases the depth of the guide (along the diagonals).
At the present level of approximation, it is straightforward to determine whether a particle with given initial conditions will be transmitted through the guide. We first relate the invariant quantity, C/b 2 , to the initial conditions. To lowest order in ǫ, the initial values of a 1 , θ 1 are related to the initial conditions via and similarly for a 2 , θ 2 . With these, we can form ξ(0) and γ(0) and hence C/b 2 from (51). An equivalent approach is to realize that, to lowest order in ǫ Knowing C/b 2 we use (53) and (54) to find the maximum value of ξ for the particle, and hence the maximum values of r 1 and r 2 . These maxima will eventually be reached (unless ǫb 2 = 0 exactly), and so the particle will be guided indefinitely if these are smaller than the boundaries of the guide, and will be lost with certainty if they are larger. We applied this algorithm to a large number of particles to find the set that are transmitted by a guide in the presence of a very small transverse coupling. We find that the coupling reduces the acceptance of an infinitely long guide to 72% of its acceptance in the absence of nonlinear forces. This is a general result that depends only on the coupling term being small enough that our lowest order approximation is valid. In the absence of coupling, the acceptance ellipses in the (X, X ′ ) and (Y, Y ′ ) phase spaces are uniformly dense -all particles that lie inside both ellipses are transmitted. This ceases to be true when the coupling is turned on. In fact, the sizes of the ellipses that enclose the stable particles do not change, but these particles are no longer distributed uniformly within them. For example, a particle that is close to the boundary of the ellipse in (X, X ′ ) may be driven out of the guide if it also happens to lie close to the boundary in (Y, Y ′ ), but it will be successfully transmitted if it happens to lie close to the origin in (Y, Y ′ ). Thus, the phase space acceptance in one direction has a dense region around the origin, becoming less dense towards the boundaries.
By eliminating the independent variable, σ, from (48) and (49), we have obtained a great deal of insight into the coupled motion of the particles. However, we do not yet know the wavelength of the slow oscillations of amplitude and phase. To find this, we need to study the solutions of (48) and (49) with b 1 = 0. There appears to be no straightforward analytical solution, so we solved the equations numerically for a set of initial conditions corresponding to the full range of the contours seen in figure 6, and for several values of the single parameter appearing in these equations, ǫb 2 E 0 /(4µ). In this way, we found the wavelength for the slow energy-exchanging oscillations (in units of s) to be Here, g is a number that depends only on the value of C/b 2 , i.e. on the particular (ξ, γ) contour that the particle travels along. For most of the contours, the value of g is between 5 and 10. The shortest wavelengths (g ≃ 5) are obtained for contours near the two maxima in figure 6, and near the edges of the plot where ξ is 0 or 1. The longest wavelengths are found for values of C/b 2 near 1/4 where the crossover from the closed contours to the open ones occurs. Taking an example where µ = 0.5, ǫb 2 = 0.1, E 0 = r 2 1 + r 2 2 = 1 and g = 5 we find λ s to be 100 unit cells of the guide.

Larger coupling coefficient
Having understood how the acceptance is reduced by a very small nonlinear coupling, we now consider larger values of the coefficient χ. Applying the multiple scales method to higher orders in ǫ becomes quite awkward mathematically, so we turn to numerical methods. Then we have a choice: we could either solve the exact differential equations given by (24), or we could use the small-micromotion approximation and solve (26) instead. In the former case, the step size needs to be chosen small relative to the micromotion wavelength, while for the latter it need only be small relative to the macromotion wavelength. For this reason, the approximate method is computationally faster by the ratio of the two wavelengths. Figure 7 shows the results obtained by both methods for a gapless guide containing 100 lenses each having ηL = 1. The two methods give similar results, the greatest difference being for large positive values of χ where neglect of the micromotion underestimates the true acceptance.
In line with our discussion in section 7.1, the acceptance drops on either side of the χ = 0 position, due to the exchange of energy between transverse modes. This loss mechanism is fully established once χ is large enough that a cycle of the slow amplitude and phase oscillations is completed within the length of the guide. The width of the sharp peak depends inversely on this length, a longer guide resulting in a narrower central peak. Away from the central peak there is a striking asymmetry between positive and negative values of χ. This happens because there are two effects that cause the acceptance to depend on χ, one which reduces the acceptance for both positive and negative χ, and another which reduces the acceptance on the positive side but increases it on the negative side. The first effect is the exchange of energy between the two transverse oscillations, already discussed above. Importantly, the wavelength of the amplitude and phase oscillations, λ s , is a function of the particle's transverse energy, E 0 , and its value of C/b 2 . Those particles with C/b 2 near 1/4 will have longer energy-exchanging wavelengths than others, and so this loss mechanism continues to be active as |χ| increases. As a result the central peak in the figure has broad wings where there continues to be more particle loss with increasing |χ|, reducing the acceptance on both sides of the central peak. The second effect is a direct result of the change in shape of the effective potential with the value of χ. As χ increases, the nonlinear term makes the effective potential shallower along the diagonals of the guide, allowing an increasingly large fraction of the particles to escape along these directions. For negative values of χ, the guide becomes deeper along the diagonals, allowing extra particles to be confined. It so happens that, for negative χ, the two effects nearly cancel so that the acceptance is nearly a constant on the negative side. It is also important to note that the coupling term leaves the shape of the effective potential unaltered along the principal axes. For all negative χ, the potential depth is shallowest along these axes. In the presence of the coupling, particles tend to explore many different regions of the guide, and so for many particles, increasing the potential depth along the diagonals does not improve the confinement since these particles can still escape along the principal axes. Figure 8 shows an example of the motion in the XY -plane when ηL = 1, comparing χ = −0.05 to χ = 0. In the absence of coupling the particle moves in a closed ellipse. The coupling results in a far more complicated trajectory, with the particle visiting far more points in the plane.
Finally in this section, we consider the effect of the transverse coupling on the example guides that we examined earlier. In the case of the two-rod guide used in many previous experiments on focussing, guiding and decelerating, χ is always negative and therefore has little effect on the acceptance. This, along with its simplicity, is a great advantage of the two-rod geometry. The four rod guide with two-pole symmetry discussed in reference [10] has χ = 0.19 when the Stark shift is linear, and this reduces the acceptance to about 40% of the value it takes in the absence of the coupling (for ηL = 1 and ηS = 0). When the Stark shift is quadratic χ = 0.055 in this geometry and so the normalized acceptance is about 70%. The four rod bent guide used for guiding molecules in reference [11] has χ = 1.5 when the Stark shift is linear. This very large value for χ is a great disadvantage of this geometry since it reduces the acceptance enormously. It is, of course, necessary to consider the higher order terms to get a complete picture in this case.

Small nonlinear coefficients
We return now to the general case where both b 1 and b 2 are nonzero. Figure 9 shows contours of (50) for several different values of the ratio b 1 /b 2 , along with colour maps to indicate the values, C/b 2 , of the contours. As before, the fixed points in these plots correspond to the coupled stationary modes. The ones at (ξ, γ) = ( 1 2 , 0), ( 1 2 , π) are the diagonal in-phase and anti-phase normal modes, while the fixed point at ( 1 2 , π 2 ) represents the clockwise and anticlockwise circular modes. When b 1 /b 2 is small, as in (a), the plot closely resembles the b 1 = 0 case shown in figure 6. The main difference is that a greater fraction of the contours encircle the two maxima that correspond to the normal modes. Part (b) of the figure shows the contours for b 1 /b 2 = 0.5. The saddle point has been replaced by a minimum meaning that the circular modes, as well as the diagonal modes, are stable. All motion involves a complete exchange of energy between one transverse direction and the other. As b 1 /b 2 is increased further, the area of the plot containing contours centred on ( 1 2 , π 2 ) grows. Eventually, once b 1 /b 2 = 1 as in part (c) of the figure, this area includes all the contours. With a ratio greater than 1, e.g. b 1 /b 2 = 2 as shown in (d), the diagonal modes become unstable, and a larger fraction of the contours are associated with these modes.
It is clear from figure 9 that the coupled stationary modes can be stable or unstable, depending on the ratio of b 1 to b 2 . To find the stability condition, we form the Hessian matrix of the left-hand-side of (50). A fixed point is a stable one if the determinant of this matrix is positive. Using this procedure it is straightforward to show that the normal modes are stable provided b 1 /b 2 < 1, and the circular modes are stable when b 1 /b 2 > 1/3.
We can also easily find the spatial frequencies of the coupled stationary modes, following the same procedure as before. The frequency in the X direction is ω = µ + ǫ dθ 1 /dσ. With the help of (47c) we find the normal mode frequency to be and the frequency of the circular mode to be Comparing these with equations (36) and (52), we see that, to this order in ǫ, the change in frequency resulting from the inclusion of both cubic and coupling terms together is just the sum of the changes introduced by each term alone.

Larger nonlinear coefficients
To find the 2D phase space acceptance over a wide range of ζ and χ we solved the equations of motion for a large number of particles using the effective potential given by (29). The results for a 100-lens gapless guide with ηL = 1 are shown by the points in figure 10, for −0.025 ≤ ζ ≤ 0.15 and −0.3 ≤ χ ≤ 0.3. For all values of ζ, the results are rather similar to the ζ = 0 result shown in figure 7 -there is a narrow spike centred at χ = 0, a decline in acceptance for positive χ and a plateau for negative χ. This similarity suggests a hypothesis, which we could write symbolically as A(ζ, χ) = A 0 F ζ (ζ)F χ (χ). In words, the hypothesis is that the total acceptance is the product of the one obtained from the linear theory, A 0 , the factor F ζ given by (42b) and (42c), which is independent of χ and accounts for the cubic term, and a third factor, F χ , that accounts for the coupling term independently of ζ and is plotted in figure 7. The lines in figure 10 are plots of this product form, showing that for a large region of the parameter space this simple method agrees with the full numerical calculation to within a few percent. The largest discrepancies occur when both χ and ζ are large and positive. Here, the numerical calculations show that the acceptance is considerably larger than the product form suggests, e.g. 70% larger when ζ = 0.15, χ = 0.3. Away from this region, the next largest discrepancies occur when χ is small and negative. The numerical calculations show that in this region the acceptance tends to be larger than the product form predicts, by up to 20%. We are, at last, ready to put everything together to find the best way of building alternating gradient guides, decelerators and traps. We first found that, in the linear theory, the 2D phase-space acceptance is proportional to η 2 , and so to the value of |a 3 |. However, making a 3 as large as possible may not be the best strategy, because of the detrimental effects of positive ζ and χ which are related to one another through a 3 : χ + 3ζ = 2n|a 3 |. If we make a 3 too large, the increased acceptance suggested by the linear theory will not be realizable due to the resulting large value of χ. Decreasing χ by increasing ζ does not necessarily help, since large positive ζ also reduces acceptance.
In the effective potential picture, increasing a 3 increases the depth of the potential, but increasing either ζ or χ decreases the depth. We are only free to choose two of the three parameters, a 3 , ζ, χ. What values should we choose to optimize the acceptance?
We will look at some examples for particular values of the operating parameters, ηL and ηS.
We start with the two-rod geometry. This is the easiest case to consider since there is only one free parameter, a 3 , determined by the the ratio of the rod radius, R, to the gap between the rods, 2r 0 : a 3 = (r 0 /R)/(2 + r 0 /R) and a 5 = a 2 3 . The two nonlinear coupling coefficients are therefore ζ = 2|a 3 | and χ = −4|a 3 | when the Stark shift is linear, and ζ = 3|a 3 |, χ = −5|a 3 | when the Stark shift is quadratic. The calculation is hugely simplified by χ being negative because we know that the acceptance is insensitive to χ when it is negative, and we have an approximate analytical expression to describe the effect of the cubic nonlinearity. Our question can therefore be answered analytically with good accuracy. Figure 11(a) shows the transverse acceptance versus a 3 for a two-rod gapless guide with ηL = 1, in the case of a linear Stark shift. The unit of the acceptance is Rr 2 0 and is independent of a 3 . See equation (4) for the definition of R. The solid line in the figure is the analytical result obtained in the absence of nonlinearities, (17). The dashed line in the figure is also a completely analytical result -we take the result given by the solid line, multiply by the correction factor given by (42b) and (42c) to account for the cubic nonlinearity, and finally multiply by 0.9 to account for the negative value of χ, as suggested by figure 7. As a 3 is increased from zero the acceptance first increases linearly, just as expected from the theory when the nonlinearities are neglected. However, as a 3 increases further, the growing cubic nonlinearity halts and then reverses the ascent of the acceptance, which reaches its maximum value near a 3 = 0.03. We note that the rod radius needs to be 16r 0 to obtain this value of a 3 . The maximum acceptance is approximately 0.008Rr 2 0 . We will examine this value in the context of some typical experiments at the end of the paper. The round points in the same figure show the results of numerical calculations of the transverse acceptance using the effective potential approach. Our analytical result agrees very well with these numerical calculations. The square points show the results obtained from a full numerical calculation, including the micromotion. As in previous examples, the inclusion of the micromotion increases the acceptance, but the overall shape of the plot remains unchanged, and the maximum of the acceptance still occurs close to the same value of a 3 . We have shown the result for a gapless guide operated at ηL = 1 which is close to the optimal lens length (see figure 1). Results for other cases are just as easily obtained using the same procedure. In general, the optimal value of a 3 increases as ηL is increased towards the edge of stability. Figure 11(b) provides the same information for a guide with gap lengths equal to the lens lengths. This is the situation typically found in AG stark decelerators. Again, we chose the parameters close to the optimal ones in the linear theory, specifically ηL = ηS = 0.7. The figure shows exactly the same trends as for the gapless case, but the acceptance is now maximized when a 3 ≃ 0.07 − 0.08, requiring rods of radii 6-7r 0 . Again, the numerical simulations based on the effective potential agree very well with the analytical result, while the full numerical simulations show that the true acceptance is a little larger.
Returning to the more general case, we will take a 3 and ζ to be the two free parameters, χ being fixed once these two are chosen. We will concentrate on the example of a gapless guide operated with ηL = 1. Then, there are no new calculations to do since all the information we need is already given in figure 10. For each value of a 3 and ζ, we calculate the acceptance in the linear case, which is just linear in a 3 , and then read off the multiplication factor from figure 10.
We present the results, for both linear and quadratic Stark shifts, in figure 12. Let us concentrate first on the linear case, part (a) of the figure. We know from figure 4 that negative values of ζ increase the acceptance above the value given by the linear theory, and it is tempting to conclude that we should design the guide with negative ζ. However, negative ζ implies large positive values for χ which greatly reduces the acceptance. As a result, the optimal acceptance is small and occurs for small values of a 3 , as illustrated by the results for ζ = −0.025 shown in figure 12(a). The figure shows that increasing the value of ζ increases the maximum acceptance, and also the optimal value of a 3 . This is because the beneficial effect of a smaller χ outweighs the detrimental effect of a larger ζ. The acceptance reaches its maximum value of 0.025Rr 2 0 when ζ ≃ 0.075 and a 3 ≃ 0.11. At this point, χ is very close to zero. For values of a 3 below the optimum, χ is negative and so has very little influence on the acceptance, which consequently increases rather linearly with a 3 . As soon as χ becomes positive further increases in a 3 only reduce the acceptance because of the associated increases in χ. When ζ > 0.075, the same trends are seen -the acceptance increases linearly with a 3 up to the point where χ is zero, and then decreases. Because increasing ζ reduces the acceptance, the slope of this linear rise is smaller than before, the acceptance is maximized at higher values of a 3 , and this maximum value is reduced. For comparison, the result obtained in the absence of nonlinearities is shown by the solid line in the figure.
The results for a quadratic Stark shift are shown in figure 12(b). The figure shows exactly the same characteristics as for the linear Stark shift, but with all features occurring at lower values of a 3 . The acceptance obtains its maximum value when ζ ≃ 0.075 and a 3 ≃ 0.06. As before, χ ≃ 0 at this point, and this is also true of the turning points for higher values of ζ. The maximum acceptance is 0.011Rr 2 0 , about half the value obtained for the linear Stark shift.
The recipe for designing an AG guide should now be clear -find an electrode structure that produces a small positive value of ζ (close to 0.075), and a χ very close to zero, and choose lens lengths such that ηL is close to 1. If there is a need for gaps between the lenses, such as in an AG decelerator, it is best to keep the gaps short, in which case the above conclusions are not much altered.

Conclusions
In this paper we have shown how nonlinear forces alter the trajectories of neutral particles in an AG guide, and have shown how to calculate the phase space acceptance of guides with these nonlinear forces included. We have found how to design the guide so that the acceptance is optimized. To conclude, we apply our results to some example experiments, revisiting those discussed in Sec. 3.2.
Consider a beam of CaF molecules entering an AG guide with a forward speed of 350 m/s. The guide has r 0 = 1 mm, E 0 = 100 kV/cm, ηL = 1 and ηS = 0, and the electrode geometry is designed so that a 3 and ζ have optimal values, as in figure 12. The Stark shift will be linear in this case, and the ratio of the Stark shift at the centre of the guide to the forward kinetic energy is R = 1/88. The transverse acceptance will be approximately 0.025Rr 2 0 = (17 mm mrad) 2 . The equivalent trap depth of this guide is approximately 30 mK. The available acceptance is well matched to the phase space area occupied by the beam in a typical supersonic beam experiment where the molecules pass through a small skimmer before entering the guide. We can therefore expect the guide to transmit a large fraction of the total beam. For a second example we consider Li atoms launched with a speed of 10 m/s from a MOT, into an AG guide. The guide has the same parameters as above, except that a 3 is halved so that the acceptance is optimized for the quadratic Stark shift. In this case R = 0.47, and the transverse acceptance is approximately 0.01Rr 2 0 = (70 mm mrad) 2 . The equivalent trap depth of the guide is approximately 50 µK, similar to the temperature of the atoms in the MOT. Again, we can expect the guide to transport a large fraction of the atoms. Unless the length of the guide is shorter than the macromotion wavelength, in which case the transmission will be larger than the above estimates, the acceptance does not depend on the length. The guide can be made arbitrarily long, delivering atoms and molecules to an experiment far removed from the source. Our examples demonstrate that a carefully designed AG guide has an acceptance that is adequate for many common experimental situations, making it a useful tool for cold atom and molecule physics.