Analytical study of electrostatic ion beam traps

The use of electrostatic ion beam traps require to set many potentials on the electrodes (ten in our case), making the tuning much more difficult than with quadrupole traps. In order to obtain the best trapping conditions, an analytical formula giving the electrostatic potential inside the trap is required. In this paper, we present a general method to calculate the analytical expression of the electrostatic potential in any axisymmetric set of electrodes. We use conformal mapping to simplify the geometry of the boundary. The calculation is then performed in a space of simple geometry. We show that this method, providing excellent accuracy, allows to obtain the potential on the axis as an analytic function of the potentials applied to the electrodes, thus leading to fast, accurate and efficient calculations. We conclude by presenting stability maps depending on the potentials that enabled us to find the good trapping conditions for oxygen 4+ at much higher energies than what has been achieved until now.


I. INTRODUCTION
In the last few years a variety of electrostatic devices for storing and handling low energy ion beams have been designed and operated. Electrostatic storage rings [1-3], cone traps [4], Orbitraps [5], and electrostatic ion beam traps (EIBT) [6,7] are now used to study atomic and molecular metastable states or molecular fragmentation, photodissociation or mass spectrometry (see, e.g., Ref. [8] for a review). The design and study of these instruments relies nowadays mainly on computer simulations. An EIBT, as designed by D. Zajfman and collaborators [6,9], and independently by W.H. Benner [7], is a purely electrostatic trap composed of two electrostatic mirror-Einzel lens combinations, as represented on Fig. 1. This trap has many interesting features [10] : on one hand, it offers trapping of energetic particles (keV) in a well defined direction and on the other hand it is small, relatively inexpensive and has a field-free region where ions move freely and where measurements can easily be performed. It can also be used as a moderate-resolution mass spectrometer [11].
In this paper, our aim is to provide a method enabling the determination of analytical solution to the electrostatic potential in any axially symmetric configuration using the elegant method of conformal mapping. We will present the method on the EIBT, but it can be used on other sets of electrodes.
The main reason for quadrupole traps' extensive use in precision experiments, lies in the fact that their fields can easily be described by an analytical formula. It enables a deeper understanding of many subtle phenomena like frequency shifts due to space charge. The inventors of the ion trap resonator used either a matrix approach [11] or numerical simulations [12]. The former are useful for acquiring a qualitative understanding, but cannot provide detailed insight in the operating conditions, while the later may not provide enough numerical accuracy to follow the particles during several tens of thousands of oscillations inside the trap, or be too time and resources consuming to explore many different potential configurations. The tracking of particles in the EIBT is rather difficult because of the combination of a long free-flight zone between the mirrors and of two areas in which the particles slow down, stop and reverse their course on very short distances, while being subjected to strong, rapidly varying, electrostatic fields. The simulation time becomes rapidly a limitation when many potential configurations must be studied while looking for new operating conditions. Instead of two tuning parameters as in a Paul trap, we have to fix the potentials of five electrodes in a symmetric configuration (ten when each side of the trap is set differently). The space to explore is therefore too large for numerical simulations in which the field is determined by usual finite element methods, where one has to make a different calculation for each set of parameters. Moreover numerical errors accumulate and perturb the trajectory of the particle over long trapping times. In the sequel, we will show how to find a formula, depending only on these parameters, which is able to give the electrostatic potential with good accuracy. We will also show some practical applications.
Before we present our method, we will just review two ways of calculating electrostatic potentials and explain why they are difficult to use in our case.
Green's functions often yield analytical results because they allow reduction of the solution of the Dirichlet problem to the calculation of the following integral : where M 0 is the point where the potential is evaluated, ǫ 0 is the vacuum permeability, n parametrizes the direction orthogonal to the surface, U S (M ) is a given potential distribution over the surface S and G is a Green's function. However, Green's functions are only known for simple geometries, which limits the analytical approach. A very interesting method called quasi-Greens' functions has been developed in [13]. This method is based on the division of a complicated geometry into different simple shapes. The main drawback is that the final expression is given as an infinite sum whose coefficients have no closed form. In practice, the given expression, although analytical, is much more complicated than the one we will present in this article.
Another technique, is the charge ring method [14][15][16][17]. The integral form of the Poisson equation is applied to N rings representing the geometry. If N is large enough, we get a set of linear equations Φ = AQ , where Φ is the potentials applied to the rings, Q is the charge induced on each ring and A is a matrix depending only on the geometry. Once the inverse of A has been determined, the charge of each ring is know and the potential at a point r that is not on the boundary is given by : where s i is the area of ring i. Hundreds of rings a usually provide a accuracy of the on-axis potential of the order of 10 −4 (one order of magnitude better than the method proposed here, see Sec. III A). However, it is at least one order of magnitude slower to compute : even though the inverse of A is calculated only once for a given geometry, one has to evaluate numerically N integrals, which is much longer than to evaluate common mathematical functions. We have implemented both methods on the same geometry and the conformal method was 64 times faster than the ring method. The choice between those two techniques will rely on the need to improve accuracy or speed. This article is organized as follows : in Sec. II, we present an approximate method to obtain the analytical potential of an axisymmetric set of electrodes, having the same radius. In Sec. III, we explain how to use Schwarz-Christoffel method to alter the metric and obtain a set of electrodes with the same radius and in Sec. III A we solve the problem in this new space using the method of section II. Section III B contains a summary of the key steps of the whole method as well as a discussion on the improvement of the accuracy. Finally, in Section IV we show that the dynamics of the ions in the EIBT is governed by an Hill's equation and we present a stability map showing what experimental parameters lead to an efficient trapping.

II. SEPARATION OF VARIABLES AND BERTRAM'S METHOD
We start from the Laplace equation in cylindrical coordinates : where V = V (r, z) is the potential at radius r from the axis and at a distance z from the center of the trap. Using the method of separation of variables, i.e., assuming V (r, z) = R(r)Z(z), we obtain [18] : where k is a real constant.
(2) is a particular case of the general Bessel equation [19], whose solution is a Bessel function of first kind : J 0 (kr). We can use the general solution given by a Fourier-Bessel series of the form : enabling us to take into account the boundary conditions. Thus, if the potential at some radius R is known as a function of z, then a(k) can be found by means of the Fourier transform : Given that any solution of the Laplace equation in a cylindrical symmetry is also of the general form : it is sufficient to determine V (0, z) along the axis. Following Bertram [20], we assume that if we know the potential V (R, ζ) at a distance R from the axis then, the potential on the axis is well approximated by the formula : where the constant ω = 4A 0 = 1.3152, and A 0 is the first coefficient of the following development [20] : In our case, if all the electrodes had the same radius, we could use directly this method, taking V (R, ζ) as a piecewise linear function of the set of potentials However, since the radius of the first four electrodes is R = 8mm, different from the Einzel electrodes where R z = 13mm, this method does not work in the area where the radius changes. We tried to introduce a smooth function R(z) in the integral, but the difference between those results and a finite element solution always shows a large discrepancy as illustrated on  The previous method is rather simple and efficient to find the potential. Its only limitation is the need to have identical radii for all electrodes. In the next section, we will show how to use conformal mapping to place ourselves in the space where the borders have a constant radius.

III. CONFORMAL MAPPING
Conformal mapping is widely used in applied physics and chemistry. One can cite the design of airfoils : the Joukowsky transformation [21] reduces the study of the laminar flow on a complicated profile to the much easier study of a cylinder in the transformed flow, or the study of diffusive flow at micro-ring electrodes in analytical chemistry [22,23]. Applications in electrostatic potential determination are also known (see, e.g., [24][25][26][27]). However, all these works use conformal mapping in a geometry where one dimension can be considered infinite. A section perpendicular to this infinite dimension is then mapped. In this paper, we show that conformal mapping can be used to find the electrostatic potential in a space limited axisymmetric geometry by mapping the plane parallel to the axis and rotating it.
In order to use the results of II, the first step is to find the holomorphic function that maps section A, described on Fig. 3, on a rectangle. Holomorphic functions are of great interest because angles are conserved under those transformations : equipotential lines stay orthogonal to field lines.
We first use the Schwarz-Christoffel transformation [28] to map domain A, parametrized by z = x + iy onto the upper half plane where t = r + is, as represented on where K 1 and K 2 are constants to be determined. With an appropriate choice of the origin (f(0)=0), we have K 2 = 0. For K 1 , we use the boundary conditions : going from point A to point E in the Z-plane (see Fig.4) corresponds to a large semicircle of radius ρ → +∞ and θ from 0 to π in the T -plane : When taking the limit ρ → +∞, this reduces to : and so K 1 = R z /π. The second boundary condition, going from B to B ′ in the Z-plane, is expressed by integrating around BB ′ with ρ → 0 and θ going from 0 to π : Choosing a = 1 only fixes the origin in the T -plane, and it implies √ b = Rz R . Finally, we make the substitution : and the integration gives : The mapping from the W-plane to the T-plane is much simpler : and we can finally link the Z-plane to the W-plane by the following transformation : = ζ(w), The transformation z = ζ(w) gives a one-to-one mapping between points of the Z-plane and points of the W -plane and it is a conformal transform as the composition of two conformal transforms. Since ζ is not invertible, we have found a good approximating function defined on three intervals and inverse of which is : where z 0 = 0.007mm and P (z) is a 6th-order polynomial connecting continuously the two asymptotic expansions. The points at ±z 0 represent the limit of validity of the two asymptotic expansions. There is one last thing to do before one can solve the problem in the W -plane. Fig. 5 shows the image of the rectangle defined as w = u + iv with −20 < u < 8 and −π < v < 0. On the vertical bold lines, u is constant and we will call these lines iso-u. We see that the isou are distorted near 0. Therefore, the distance between the points w i = ζ −1 (z i ) is not the same as the distance between the points z i : the metric is not conserved on the border. Consequently, the width of the electrodes near the point C in the W -plane is not the physical width, which leads to errors in the potential. In order not to modify the metric on the border, one should first apply the function to the points z i so that their image have the same distance between them in the W -plane and in the Z-plane (physical plane). This function uses the fact that, on the axis, the metric is not modified from one plane to the other. It is the same idea as in [26] where a "space dependent diffusion coefficient" is used to account for the fact that real space is compressed/expanded unevenly to fit the W -plane. In the sequel, we shall use the following notation :z i = β(z i ).

A. Solving in the W-plane
Using ζ(w) we can now apply Bertram's method in the W-plane where the radius is constant. Following [20], we assume that the potential varies linearly between the electrodes. We shall come back to this approximation in the last part of this section. We now use (6) in the Wplane where the radius is constant R = π. On the border, we have : with u i = ζ −1 (z i ), wherez i = β(z i ). z i are the positions of the points on Fig. 3, V i is the potential at z i . Π(u i → u i+1 ) is a function equal to zero everywhere except between u i and u i+1 where it is equal to one. Re-placing in (6), we obtain : where : . Now that we have V W , the potential along the axis in the Z-plane is given by : The error between the numerical solution and our result is around 1%, which is enough for many applications. There are two main sources of errors. First, the calculation does not take into account the field leaking at the front and at the rear of the set of electrodes. It can be seen on Fig. 3 that the potential is not exactly zero on the axis after the neutral electrodes preceding V z and following V 1 . This explains the two dark negative zones on Fig. 7. Second, on the same figure, we see that the hypotheses we have made, concerning the variation of the potential between the electrodes, induce an error of about 100V on the border. Yet, since this error is oscillating along each border, it compensates and the error on the axis is only about 20V (Fig. 6). An attempt to enhance the potential at the borders is described in [29]. Finally, there is an error coming from the cylindrical term of Eq. (1), which is not invariant under conformal mapping. However, it is always at least one order of magnitude smaller than the two previous error types.
The best way to increase accuracy is to use a domain that takes into account the shape of the electrodes. Instead of the mapping of Fig. 5, we could have used the one of Fig. 8. Besides showing that this method can be applied to complicated geometries, Fig. 8 enables us to emphasize a crucial point : the mapping might not be analytic for a given geometry. However, since the geometry does not change, it is enough to calculate the inverse map numerically one time, find the polynomial that fits this numerical solution on the axis and then, using Bertram's formula we have an analytic expression whose parameters are the electrodes' potentials. With this mapping, the accuracy is approximately 0.1%. Solving numerically a Schwartz-Christoffel map will not be treated here because we chose a complete analytical case to show all the details of the method. One could refer to [28,30] for an extensive review on all the numerical techniques involved. Note that all the results presented in the sequel, use the simple case of Fig. 8.
B. Summary of the method 1. Write the the Schwartz-Christoffel transformation adapted to the geometry to obtain z = ζ(w). For more details see [28].
2. If the inverse w = ζ −1 (z) is not straightforward, it can be well approximated by a polynomial.
3. Let z i = x i + iR(x i ) be the points of the Z-plane defining the position of each side of the electrodes, one has to apply the function β(z) = ζ(ζ −1 (z)−iπ) to obtain the set of pointsz i = ζ(ζ −1 (z i ) − iπ) 4. Transpose the problem from the Z-plane to the W -plane where w i = ζ(z i ) = u i + iv i . 5. Use the formula (10) to obtain V W (u, −P i).
6. The potential of one point z = x + 0 × i on the axis in the Z-plane is V Z (z, 0) = V W (ζ −1 (z)).

IV. APPLICATION : STABILITY MAP
The main difficulty to tune an EIBT is due to the large number of parameters implied in its manipulation : five potential values on each side, the energy of the ions, their temperature and their charge-to-mass ratio. The designers of the trap used an optical model consisting of mirrors and lenses [10], yet, since the focal length is not linked to the values of the potential, the behavior of the EIBT, given a set of potentials, is unpredictable. Until now, experimentalists had to use simulation soft-  Figure 9: (Color online) These curves show 1 − |∆| as a function of Vz (potential applied to the lens) in two distinct cases : (right curve, +) is the trapping efficiency of Ar + at 4.2keV with the same conditions as in [32]. The set of points is an reproduction of the experimental data in their Fig.3 a. (left  curve, o) represents the trapping efficiency of Ar + at 1.2keV with the same conditions as in [10]. The set of points is an reproduction of the experimental data in their Fig. 8. There is very good agreement between theory and experiment.
ware like SIMION R [31] to determine optimal trapping conditions [10]. Since a finite element calculation has to be achieved each time a parameter is changed, this method turns out to be fastidious. Worse, simulating the trajectory of ions going back and forth is much more difficult than simulating a beam because the errors on the position accumulate and become comparable to the size of the trap. Even with the smallest step and a recursive method, SIMION R [31] was not able to enforce energy conservation after a few hundreds oscillations (the error was about 100eV). From Eq. (3), we can limit ourselves to the second order. Higher order terms can be neglected as long as r < 8mm, which is the aperture of our trap. The potential in the trap is then given by where V (z) is the potential along the axis. It follows that the trajectory of one ion in the trap is determined by the following set of equations : The second term of the right side of equation (12) is small compared to the first (for r < 8 mm) and can thus be neglected. We obtain the longitudinal motion of the ion z(t). Substituting z(t) in (13), we obtain an Hill's equation [33] : the term in parentheses being a periodic function of period T . This equation arose in the study of the moon's dynamic and the usual method to discuss the stability of its solution consists of calculating infinite determinants [34]. The principal matrix of this equation is : where ψ 1 (t, t 0 ) is the solution of (14) with initial conditions ψ 1 (t 0 , t 0 ) = 1 andψ 1 (t 0 , t 0 ) = 0 and ψ 2 (t, t 0 ) with ψ 2 (t 0 , t 0 ) = 0 andψ 2 (t 0 , t 0 ) = 1. Liouville's formula [35] shows that : det M (t, t 0 ) = 1 (16) and therefore the characteristic equation of the monodromy matrix M (t 0 + T ) is given by [35] : where : [36], we know that if ∆ 2 > 1, one of the two solutions is unbound but if ∆ 2 < 1 there are two solutions : where p ± (t + T ) = p ± (t), γ = Im( 1 T log(x + )) and In conclusion, Hill's equation is stable, and thus trapping can be observed, when |∆| < 1. We now compare these theoretical results with experiment. Figure 9 shows a comparison with the results published by other groups using the same kind of trap. These curves show 1 − |∆| as a function of V z (potential applied to the lens).They show a perfect agreement with two independent groups [32], [10]. We also tried to reproduce the data of [37] using H + 2 at 1.0keV and a negative potential on the Einzel lens. We also predict three stability intervals, approximately at the same position, because we have not taken into account the differences in the geometry of their trap. However we notice that the method works also with negative potentials. Our method enables to plot these curves in less than a second whereas the cited authors had to make a SIMION R [31] simulation for each point. Figure 10 shows the stability map depending on two parameters V 1 and V z (the respective potentials of the rear electrode and of the Einzel lens). This is the equivalent to the famous Ince-Strutt diagram used to tune quadrupolar traps [38]. The white dots indicate settings where trapping was experimentally observed. We fixed the value of V 1 and scanned V z . Trapped ions go through a ring at the center of the trap and induce a current. We then analyze this amplified current with a spectrum analyzer : if we see a peak corresponding to the oscillating movement of the ions, we mark this position with a white dot. The radius of each dot is proportional to the trapping efficiency. There is a shift between theory and experiment, which was first thought to be caused by a technical problem : it is difficult to monitor with a good accuracy our power supplies up to 8kV, especially because they raise in a   few nanoseconds. However, a recent improvement of the model, taking space charge effects into account, seems to explain this shift. These results will be presented somewhere else. Nonetheless, we could not explain the absence of trapping on the upper part of the stability zone II. We tested to see if higher order terms of Eq. (3) could account for this, without success. No trapping can be seen in zone III, but as shown in Fig. 11, this region corresponds to very peculiar closed trajectories, which seem more theoretical than observable.

V. CONCLUSION
In this article, we have given a method to calculate analytical solution to the Laplace equation in axially symmetric devices. This method is very general and can be used for various parts constituting a beam line. We successfully applied it to the EIBT and showed that the electrostatic potential is given by a formula depending on the five electrodes potentials. The formula has been compared with many finite element calculation where the potentials have been changed on the whole achievable range (from 0V to 8000V for each of the five potentials) and, in the region where the ions can move, the error is never greater than 1%. This analytical expression of the potential inside this trap provides users with a much more powerful tool to study and optimize this novel kind of trap. As an example, we study the stability of the trap and show that they agree with experiments giving a fast and easy way to predict trapping parameters.