Investigation of two-frequency Paul traps for antihydrogen production

Radio-frequency (rf) Paul traps operated with multifrequency rf trapping potentials provide the ability to independently confine charged particle species with widely different charge-to-mass ratios. In particular, these traps may find use in the field of antihydrogen recombination, allowing antiproton and positron clouds to be trapped and confined in the same volume without the use of large superconducting magnets. We explore the stability regions of two-frequency Paul traps and perform numerical simulations of small, multispecies charged-particle mixtures that indicate the promise of these traps for antihydrogen recombination.


INTRODUCTION
The measurable properties of hydrogen (H) and antihydrogen (H) atoms are expected to be identical as postulated by the combined charge (C), parity (P), and time (T) reversal symmetry [1]. One of the most promising tests of this symmetry is the precise comparison of the optical and microwave spectra of hydrogen and antihydrogen. The spectrum of hydrogen has been extensively studied [2,3], but precise measurements for antihydrogen are complicated by the small quantities ofH available and the technical complexity of the experimental apparatus [4]. The efficient production and trapping of cold and neutral antimatter systems is therefore a topic of great interest.
Production of antihydrogen requires the ability to trap antiprotons and antielectrons (positrons) in the same volume. The state of the art is dominated by Penning traps, where a constant homogeneous magnetic field and inhomogeneous static electric field allow for confining particles of mass m and charge Q. The ALPHA experiment [5,6] and the ATRAP experiment [7,8] rely on a variation of a Penning trap for initial particle confinement. Penning traps have the advantage of robust trapping for a wide range of charge-to-mass ratios, while also facilitating a high charge density of positrons for efficient three-body recombination. A large trap volume and superconducting magnet creates a high magnetic trap depth ( 1 K) for the resulting neutral antiatoms. A limitation, however, is the inability to trap the oppositely charged particles in equilibrium in the same volume due to the use of a DC potential for confinement along the axial trap direction. Recombination is achieved by injecting antiprotons into the positron cloud [9]. The resulting antiatoms are typically created with energy above the magnetic trap depth, and most antiatoms are lost during recombination. Typical yields in the ALPHA apparatus are several trapped antiatoms per attempt every ≈15 minutes [10,11]. The ASACUSA experiment has pursued an alternative to spectroscopy on trapped atoms with a CUSP trap [12], which uses an anti-Helmholtz field to generate a beam of spin polarized antihydrogen for eventual microwave spectroscopy atoms [13,14].
A solution to the problem of equilibrium charge overlap was previously explored in a hybrid Penning-Paul trap [15]. In that work a magnetic field and DC potential of a Penning trap confined protons and the radiofrequency Paul trap potential compensated the axial DC potential for electrons. The method still relied on a strong magnetic field for radial confinement, and to our knowledge this technique has not been continued or extended to antimatter systems.
Two-dimensional confinement of electrons has been achieved in planar devices [16], but there are no reports on confinement of ions with high charge-to-mass ratio Q m := Q/m in a three-dimensional trap. In order to achieve such three-dimensional confinement, the stability parameters of electrons -or in our case positrons -would need to be worked out. For antihydrogen formation, we approach the problem of simultaneous three-dimensional particle confinement of antiprotons and positrons with the idea of a two-frequency Paul trap. This trap design is aimed to combine the stability parameters of both particles and would also allow for charge overlap inside the trap.
A Paul trap provides a dynamical trapping potential in all three space directions, and works for positive and negative charges equally well. The problem arises from the vastly different charge-to-mass ratio of antiprotons (p) and positrons (ē). The stability of a Paul trap is characterized by dimensionless stability parameters a and q [17], which are related to the static and dynamic amplitudes, (1). Particularly interesting are planar all-rf Paul traps indicated by the geometry in lower right. Such a geometry is suitable for miniaturization with modern atom chip technology [21,22]. Atom chip technology may then also support deep traps for the produced neutral antihydrogen.
respectively, of the confining potential. Both parameters scale linearly with the charge-to-mass ratio, Q m . A is stable for 0 < q < 0.9 in case of a ≈ 0, with optimal trapping achieved around q = 0.5. A trap optimized for trapping antiprotons will have an effective q ≈ 900 for positrons and is fully unstable.
A Paul trap optimized for positrons with large Q m is theoretically stable for antiprotons, but suffers from poor equilibrium charge overlap. Particle confinement is characterized by the pseudopotential U ∝ m(a + q 2 )Ω 2 r 2 , where m is the particle mass, Ω is the frequency of the trap potential, and r is the distance from the trap center [18]. If antiprotons and positrons confined in the same region thermalize due to the Coulomb interaction and a ≈ 0, the characteristic cloud radius of antiprotons will be a factor of mp/mē ≈ 45 larger due to the dependence of q on Q m . The larger cloud radius of antiprotons will also make them more susceptible to anharmonicities of the trapping potential. We note that the ASACUSA collaboration reported work for several years on a largevolume, superconducting resonant-cavity Paul trap for antihydrogen production [19]. More recent reports indicate the intention to use this trap for spectroscopy of antiprotonic helium,pHe + [20].
In this paper we discuss features of a two-frequency Paul trap with an infinite, perfect quadrupole potentialthat allows simultaneous confinement of antiprotons and positrons and allows the antiproton and positron cloud sizes to be matched. Trap frequencies are chosen such that positrons are confined by the high-frequency component of the trap potential and protons are primarily confined by the low-frequency component. This allows the pseudopotentials for antiprotons and positrons to be adjusted independently. Our work was partially inspired by the preliminary discussion of two-frequency Paul traps in Ref. [23].

TWO-FREQUENCY PAUL TRAP
The quadrupole potential of a two-frequency Paul trap takes the form where r 0 is a geometric scale for the trap. We will choose frequencies such that the fraction Ω 2 /Ω 1 is a number η ≥ 1. The potential can be created by a system of hyperbolic electrodes with cylindrical symmetry, or approximated by more practical geometries as indicated in Fig. 1.
In the initial discussion only motion along the xdirection is considered. From the symmetry of the potential these results will also hold for the y-direction, and may be extended to the z-direction by scaling stability parameters by a factor of −2. The equation of motion for a charged particle in the potential of Eq. (1) can be written (2) where τ = Ω 2 t/2 is a normalized time, are low-frequency (q 1 ), high-frequency (q 2 ), and DC (a) Mathieu parameters. The time derivative indicated bÿ x(τ ) is with respect to the time variable τ . Equation (2) is a specific example of a Hill differential equation: a second order, linear differential equation with periodic coefficients.

Qualitative discussion
Setting q 1 = 0 recovers the well known Mathieu equation for a single-frequency trap. A trapped particle undergoes high-frequency motion at multiples of the trap drive frequency, Ω 2 , in addition to a slow macromotion at a secular frequency of ω 2 = (1/4)(a + q 2 2 /2)Ω 2 2 . In an optimal trap a ≈ 0, and the secular frequency is ω ≈ q 2 Ω 2 /(2 √ 2) [17]. If we operate the trap at q 1 ≈ 0, q 2 ≈ 0.5 we can choose the ratio of trap drive frequencies, η, large enough so that the secular oscillation frequency ω Ω 1 . In this regime it is possible to treat a non-zero q 1 as a slowly varying DC term in addition to a.
We now consider Eq. (2) from the perspective of two charged particles, A and B, with opposite charges and masses m A < m B . To facilitate the discussion we introduce the notation q A,B 1,2 and a A,B to distinguish trap parameters for the light particle, A, and the heavy (2). The inset shows the calculated FFT used to extract secular frequencies for q1 = 0.004, q2 = 0.02. Vertical lines indicate the driving frequencies. The lowest-frequency peak is the secular frequency. We do not include Coulomb interaction here but do so in Eq. 18.
The same relationship holds for a, although this DC parameter will be set to zero for most of the manuscript.
If the high-frequency confinement is optimized for the lighter particle A and for η < m B /m A , the secular oscillation frequency of the heavy particle B due to q B 2 will be slower than Ω 1 . Therefore we must also consider the dynamic pseudo potential of q B 1 for B. We may write the effective potentials experienced by each particle as [17] where we typically choose η to be large enough that q A 1 is close enough to zero that its effect on particle A can be ignored. The validity of the pseudopotential approximation is discussed in detail in Refs. [24,25]. The pseudopotential for particle B in Eq. (5) can be arrived at by alternately considering the limiting cases where q B 1,2 go to zero. When q B 2 = 0, Eq. (2) can be rewritten with another time transformation τ → τ = Ω 1 t/2 to obtain the η 2 factor multiplying q B 1 . We confirmed the scaling of these pseudopotentials by numerical determination of the secular frequency, shown in Fig. 2. Equation (5) also illustrates the ability of a two-frequency trap to create overlapping yet independent potential wells for two Sketch of pseudopotentials for particles A and B assuming a frequency ratio of η = 200 and mass ratio, mB/mA ≈ 1836 (matching that of positrons and antiprotons). Parameters were q A 2 = 0.37, q A 1 = 0 (solid lines) and q A 2 = 0.37, q A 1 = 0.02 (dashed lines). The small perturbation of UA(x) by q A 1 oscillates with frequency Ω1 as indicated by the shaded band and will average to zero.
charged particles with an appropriate choice of trap parameters and frequency ratio, visually demonstrated in Fig. 3. This opens the possibility for combining and separating different species of charged particles with high efficiency.
In the remainder of the manuscript we discard the superscript notation. Where relevant, it is assumed that q 1,2 and a refer to trap parameters for the lighter charged particle.

Floquet Theory
To determine the stability of a two-frequency trap we define the vector u(τ ) = [x(τ ),ẋ(τ )]. Equation (2) may then be written in matrix forṁ where If η is a rational number it can be represented as an irreducible fraction m/n, where m and n are both integers and m ≥ n. In this case the matrix P(τ ) has periodicity T = mπ such that P(τ + T ) = P(τ ). General closed-form solutions of Eq. (7) do not exist, however it is possible to use Floquet theory to make statements about the existence of bound solutions for particular values of equation parameters. The existence of bound solutions implies a stable trap. A discussion of Floquet theory may be found in most differential-equations texts, such as Ref. [26], and application of Floquet theory to Paul traps in Refs. [24,25,27].
Here we simply state that the boundary of stability regions may be found by identifying parameters for which the solution x(τ ) has periodicity T or 2T . A general solution with period 2T contains all solutions with period T , so we choose where the sum over k extends from −∞ to +∞. Equation (9) and Eq. (8) lead to the identity The only way for Eq. (10) to hold for all τ is if each element of the sum satisfies this relation independently. Equation (10) may therefore be summarized as a matrix equation where the elements of the infinite matrix D are given by

Stability diagrams
Equation (11) is equivalent to the statement det(D) = 0.
Stable trap operating parameters can be identified by finding parameters that satisfy Eq. (13). Although D is an infinite matrix, a matrix of size (10m + 1) × (10m + 1) centered around D 00 was found to be a sufficient approximation for evaluating stability. We find that parameters where det (D) > 0 correspond to stable traps, and parameters where det (D) < 0 correspond to unstable traps. This means only det (D) needs to be evaluated, and Eq. (13) does not need to be solved exactly. Using larger matrices changes the stable area by less than 0.1%. The matrix evaluated can still be large, but most elements are zero and programs such as Mathematica or MATLAB have efficient tools for computations with sparse matrices. The stability diagrams in q 1 , q 2 space are shown in Fig. 4 for integer and rational frequency ratios. These diagrams show a structure of unstable resonances that increase in density with the frequency ratio. Near the q 2 axis these unstable features correspond to a parametric resonance condition between the secular oscillation frequency ω of the particle due to q 2 and the frequency The unstable resonances are finely resolved in the left diagram and extend all the way to the q2-axis, while in the diagram to the right the stability arms merge before reaching the q2axis. For further comparison between the two methods, the analytic calculation took 15 s and the numerical calculation took 9700 s, a difference of nearly three orders of magnitude. Regarding the scaling of the computation time, both methods take significantly longer for larger η. Ω 1 . These resonances become infinitely thin, but extend all the way to the q 2 axis. For large frequency ratios this structure indicates that a damping mechanism will be necessary for long-term stable operation of a twofrequency trap. The stability region for rational numbers shares general features with the closest integer ratio diagrams, with a significantly denser resonance structure as can be seen in Fig. 4.
It is important to note that stability for the light particle does not guarantee stability for the heavy particle. Stability calculations can easily be evaluated for both particles, however simultaneous stability can be obtained with a general guideline. Revisiting particles A and B: if m B /m A η, stable values of q 2 and q 1 for the light particle are also stable for the heavy particle if q 1 η 2 m A /m B < 0.9. This can be seen by setting q 2 = 0 in Eq. (2) for the heavy particle, and recovering the regular Mathieu equation with the transformation τ = τ η.

NUMERICAL SIMULATIONS
In the previous section we require the determinant of the matrix D in Eq. (11) to be zero. With this method we are forced to approximate an infinite matrix with a large but finite representation.
Another issue of the matrix determinant solution is that we are bound to rational frequency ratios η, which means that in practice an irrational frequency ratio can only be approximated by rounding to a nearby rational number. While the matrix method works for any rational number in principle, even short decimal numbers require the evaluation of impractically large matrices. For example, finding stable regions for η = 3.14159 would require the evaluation of a 314160 × 314160 matrix for comparable accuracy to the stability diagrams in Fig. 4. In this case direct numerical integration of the equation of motion is more efficient.
The method we chose for this purpose is numerical integration of Eq. (2) in Mathematica. The numerical method then not only makes solutions with an irrational η possible, but also allows for modifications of the equation of motion to incorporate effects such as damping, multiple interacting particles, real trap geometries, or electrical noise that heats the particles.
We evaluate the equation of motion for a time interval [τ 0 , τ 1 ], during which the amplitudes of a particle's oscillation A 0 at τ 0 and A 1 at τ 1 are either of the same order of magnitude or escalate to a difference of many orders of magnitude for, respectively, stable or unstable (q 1 , q 2 ) combinations. This time interval's length has to be chosen in respect to the available computation power. Simultaneously the precision with which we can determine stability increases with the length of this time interval due to the fact that a solution close to the border of stability takes a long time to diverge. We compromised between a short solution time and precision by choosing ∆τ = τ 1 − τ 0 < 1000 (with units 2/Ω 2 ). The computer had 64 GB RAM and 10-core processor, running simulations for a time interval of [0, 600].  c)). The stability region broadens with increased damping in such a way that the area of stability between two unstable resonances extends, widening the stable q2 region and the maximum of stable q1. (bottom row) Stability diagrams of a charged particle in a two-frequency Paul trap with η = 45 and an additional magnetic field B0 along the z-axis (see text), calculated with the matrix determinant method. Values of the magnetic field parameter correspond to p = 0.1 (d)), p = 0.3 (e)) and p = 0.7 (f)). The stability region grows larger with increased magnetic field coupling p so that at q2 = 0 the maximum value of q1 increases, while the maximum value of stability for q2 decreases.
We define the parameter s(τ 0 , τ 1 ) to determine whether a particular combination of (q 1 , q 2 ) is stable. This stability parameter is evaluated by integrating the square of the solution over [τ 0 , τ 1 ] and dividing it by the intervals length ∆τ : We choose s(τ 0 , τ 1 ) < 10000 to identify a combination of (q 1 , q 2 ) as stable. This condition catches oscillations that are close to the parametric resonances and have large amplitudes. For s(τ 0 , τ 1 ) > 10000 we expect the amplitude to continue growing for times beyond the integration limit τ 1 , but ultimately the value of this condition is arbitrary and chosen based on available computing resources. This stability parameter s(τ 0 , τ 1 ) is evaluated by numerically integrating the coupled equation s (τ ) = 2x(τ )/∆τ simultaneously with the equation of motion. This significantly reduces computational overhead. Using these rules to decide if a particular combination of (q 1 , q 2 ) is stable or not, we can modify the original equation of motion and consider extensions to the stability question.

Numeric integration vs. matrix determinant
Without damping and coupling the numerical and matrix determinant methods agree, but the ability to treat other cases -in general modifications to Eq. (2) -and the fact that we can use irrational frequency ratios η is what sets the numerical method apart from the matrix determinant method. A downside of these numerical calculations for differential equations is the larger amount of time it takes to do a calculation, and the lower precision around sharp features. While the analytic solutions take around five minutes to yield results, the numerical solutions take anywhere from three to over ten hours.
The stability diagrams from numerical integration and matrix determinant calculations are compared in Fig. 5. The thin resonances extending to the q 2 -axis are less visible in the numerically calculated diagrams due to the limited (q 1 , q 2 ) resolution and finite integration time. These resonances are infinitely thin near the q 2 -axis -which means hardly resolvable -and evolve over large time periods that would require excessive computation power to evaluate.
Stability diagrams for irrational values of η share general features with diagrams for nearby integer values of η, but contain complex structures of instability that appear to have a fractal nature of scale-invariant complexity as shown in Fig. 6.

Equation of damped motion
Adding a damping term to Eq. (2) yields a new differential equation that can now model effects such as laser cooling or coupling of the particles' mechanical motions to a cold resonant circuit [28,29], where b = β/mΩ 2 and β is the damping parameter. Equation (15) is numerically integrated and the stability assessed by the same threshold as in Eq. (14). The effect of damping terms is pictured in Fig. 7.

Equations of motion with magnetic field
For antihydrogen production the ion trap must be operated in the presence of a magnetic field to trap the resulting neutral particles. A magnetic trap uses an inhomogeneous magnetic field to create a potential well. As a first approximation of the effect of the magnetic field on the charged particles, we consider the effect of a uniform magnetic field B 0 along the z-axis of the trap. Due to the Lorentz force the x-and y-motion are no longer independent and we get a pair of coupled equations, y(τ ) + pẋ(τ ) + (a − 2q 1 cos 2η −1 τ − 2q 2 cos 2τ )y(τ ) = 0, (17) where p = (2B 0 Ze)/(Ω 2 m) is a dimensionless magnetic field parameter related to the cyclotron frequency. Stability can be evaluated using the matrix determinant method after making a coordinate transformation to a frame rotating with frequency p/2 [17]. In this rotating frame the magnetic field acts as an extra DC potential. The effects of different values of p are shown in the stability diagrams of Fig. 7. In general the magnetic field increases stability due to the radially directed Lorentz force. This effect has been observed in single-frequency Paul traps combined with Penning traps [15].

ANTIHYDROGEN PRODUCTION IN A TWO-FREQUENCY PAUL TRAP
As discussed previously, a benefit of a two-frequency trap for charged particles with vastly different chargeto-mass ratios is the independent control over the trapping potential for each species. In a single-frequency trap the heavy particle experiences a much weaker trapping potential, resulting in a larger volume of confinement and poor overlap with the light particle cloud. Adding a low-frequency field allows heavy charged particles to be compressed without affecting the light charged particles. This also opens the possibility of independent transport of charged species within the same volume. Using a nested electrode structure the two charged species may be initially trapped in different volumes and then merged. This is advantageous for antihydrogen production, where the current procedure is to trap antiprotons and positrons in independent potential wells and then inject the antiprotons into the positron cloud.
We ran numerical calculations of trapping that account for the inter-particle Coulomb interactions. For a collection of N positrons and N antiprotons we introduce the variables r i = (x i , y i , z i ) and R k = (X k , Y k , Z k ) to indicate the position vector for the i-th positron or k-th antiproton, respectively. The equations of motions that we solve are: where the ρ = m p /m e is the ratio of masses, and Γ is a dimensionless Coulomb constant given by This constant is defined in cgs units (e = 4.8 × 10 −10 statC), and the characteristic length scale l 0 is chosen to be 10 −4 cm. This formulation of the equations of motion is inspired by the work in Ref. [30], where antihydrogen production in a single-frequency trap was considered. In that work production of transient, classically bound antihydrogen states was observed in numerical simulations for a single-frequency Paul trap optimized for positron confinement. Reference [30] also claims a trapping mechanism related to the attraction between trapped positrons and antiprotons. While the Coulomb interaction certainly provides a significant attractive force between antiprotons and positrons in close proximity, we believe the effect to be overstated in [30] and they primarily observe the weak but still-significant  confining force of an infinite-range dynamic potential on antiprotons. This is demonstrated in Fig. 8, where a single-frequency Paul trap optimized for positrons provides weak confinement for antiprotons without positrons in the trap.
In the first simulation, we consider a mixture composed equally of positrons and antiprotons. The equations are numerically integrated for 2 to 10 particles in total.The root-mean-square (rms) radius, |r i | 2 , of each particle's orbit in a simulation is evaluated and the average of these values over all particles is calculated. The simulation for each particle number is repeated 15 times with a randomly chosen set of initial velocities corresponding to a 4 K temperature. However, the temperature of particles inside the trap during an experiment is not necessarily 4 K, as heating mechanisms such as due to the RF drive in the presence of nonlinearities and the Coulomb interaction in comparison to the resistive cooling time constant may be significant.The results of these simulations for positrons and antiprotons are shown in Fig. 8. The values of q 2 = 0.37 and q 1 = 0.024 for the positrons were chosen as stable operating regions for all three frequency ratios and along all three axes. The trap drive frequencies are chosen such that Ω 2 /2π = 600 MHz and η = 170. Results were calculated with and without the q 1 term, and clearly show that for large frequency ratios the antiproton cloud can be compressed considerably without having any significant effect on the positron cloud.
In the second simulation we consider a mixture of 10 positrons and two antiprotons using the same parameters as before. The simulation shows that recombination is likely to occur at an enhanced rate with a twofrequency trap, due primarily to the increased overlap of the antiproton and positron clouds. Figure 9 shows positron and antiproton orbits with and without the lowfrequency potential. The low-frequency potential does seem to moderately affect the positron orbits, possibly due to the increased rate of energy changing collisions between positrons and antiprotons. The average rms distance between antiprotons and positrons is plotted in Fig. 10.
To quantify relative recombination rates, we calculate the energy of each reduced-mass positron-antiproton pair, where v rel is the relative velocity of a given positronantiproton pair and r rel is the distance between them. A negative energy corresponds to a classically bound state. A list containing the energy of every positron-antiproton pair E ij (t) is calculated as a function of time, and the number of negative energies is recorded at every point in the simulation. Figure 10 shows the cumulative number of bound pairs produced during the simulation, with and without the low-frequency confining potential. A non-zero q 1 results in bound antiproton-positron pairs appearing 5× more often than for q 1 = 0.

Charge density
The Coulomb interaction between positrons and antiprotons is conservative, and to create a bound state from initially unbound particles a third party must remove energy from the system. The spontaneous emission of photons is one possible mechanism, but is a slow process compared to the characteristic close-interaction time for charged particles. The primary mechanism for antihydrogen production in the ALPHA experiment is a three-body scattering process that relies on a highdensity positron plasma [31]. More than 10 6 positrons are trapped with a density of > 10 7ē /cm 3 [10].
We estimate the achievable positron densities in a Paul trap by assuming a force balance between the Coulomb and pseudopotential forces for a positron at the edge of the positron cloud, where the equation is written in cgs units. The antiproton density is assumed negligibly low. This leads to an expected charge density ρē = (3/8π)mω 2 /(e 2 ). If we assume a trap drive frequency of Ω 2 /2π = 6 × 10 8 Hz, and q 2 = 0.37, the positron secular frequency is ω/2π ≈ 80 MHz and the maximum positron density is ρē ≈ 2 × 10 8ē /cm 3 , significantly higher than in the ALPHA experiment. This suggests that three-body recombination is also a viable option for an all rf trap.
Extrapolating the conclusions of our simulations from several charged particles to several million charged particles is not straightforward, for instance instability may arise due to excitation of collective particle oscillations by The appearance of bound pairs is five times more frequent with the non-zero q1. Note that the average rms amplitude does not increase significantly, suggesting that heating on the recombination timescale can be neglected. the dynamic potential. These problems may be avoided, however, by using more compact trap configurations, for instance planar traps fabricated with atom-chip technology [21,22]. These can localize charged particles more precisely and permit optical access without restrictions enforced by large magnets. Increased overlap and localization would facilitate studies of other recombination mechanisms, such as resonantly enhanced photoinduced recombination [32][33][34]. Recombination in a smaller volume may also simplify direct laser cooling of ground-state antihydrogen that may be produced [35,36]. Together these methods may reduce the number of positrons and antiprotons necessary, and increase the precision and rate of experiments with trappedH.

Heating, cooling and non-linear resonances in real Paul traps
Real Paul traps have effects that perturb the quadrupole potential we considered. Perfect hyperbolic electrodes are hardly realizable and often approximated by spherical surfaces. Those real trap geometries make for a potential different from the perfect quadrupole and add structures of instabilities to the stability diagrams [37]. Mathematically, the additional instabilities can be calculated as resonances of the secular frequency with terms of the series expansion of the potential. Space charge effects of the particles being confined in a small volume of the trap also contribute to deviations from the assumed harmonic potential. The effects of nonlinear resonances caused by deviations from the perfect quadrupole in Paul traps are well known for single frequencies, but have to be worked out for a two-frequency trap.
In our simulations using only up to ten particles we do not see any signatures of these instabilities, but for larger particle numbers inside the trap, the probability of collisions between particles is increased and can place the momentum of the particles out of phase with the potential. The particles then gain energy from the field instead of being forced on a stable orbit. This energy gain translates to a higher equilibrium temperature of the particle cloud and can have considerable influence on the cloud size. This rf heating effect is dependent on the particle number, T ∝ N 2/3 [38].
How big the influence of rf heating, trap imperfections and (sympathetic) cooling in our trapping constellation is, is a topic of future investigations.

SUMMARY
We have discussed the potential of two-frequency Paul traps for the simultaneous trapping of positrons and antiprotons for recombination to antihydrogen. Stable regions in the trap parameter space have been identified and confirmed using independent methods based on Floquet theory and direct numerical integration of the equations of motion. Floquet theory provides stability maps for any rational frequency ratio, while numerical integration provides stability maps of reduced precision for any possible frequency ratio. Additional effects such as those of damping and magnetic fields were also investigated. We have further confirmed that two-frequency potentials enable charged particles with very different charge-mass ratios to be trapped simultaneously in volumes of similar size, a significant improvement over single-frequency Paul traps. The influence of this control on the rate of antihydrogen production is a topic of continued investigation.
The feasibility of two-frequency Paul traps for antihydrogen recombination is a topic that merits further study; a number of important questions need to be answered. What effect does micromotion in an rf trap have on the energy spectrum of the produced antihydrogen? How will the trap electric fields contribute to ionization loss of Rydberg states? What trap depths can be achieved with real electrode geometries? Can atomic ion species be used for sympathetic cooling?
Investigation of recombination dynamics in twofrequency traps can be pursued initially with ions such as 9 Be + or 40 Ca + and electrons. These positive ions have convenient laser-cooling wavelengths, which simplifies many technical challenges of detection and cooling. If the techniques can eventually be proven with protons and electrons, the extension to antihydrogen is mainly a question of technical complexity and availability of antiprotons.
This work was supported in part by the DFG DIP project Ref. FO 703/2-1 1 SCHM 1049/7-1. DB, FSK, and NL acknowledge financial support by the Cluster of Excellence PRISMA at the Johannes-Gutenberg Universität Mainz. FSK acknowledges support from the DFG within the project BESCOOL. NL was supported by a Marie Curie International Incoming Fellowship within the 7th European Community Framework Programme.