Loading a linear Paul trap to saturation from a magneto-optical trap
J. E. Wells, R. Blümel, J. M. Kwolek, D. S. Goodman, W. W. Smith
LLoading a linear Paul trap to saturation from a magneto-optical trap
J. E. Wells,
1, 2
R. Bl¨umel, J. M. Kwolek, D. S. Goodman,
2, 4 and W. W. Smith W. M. Keck Science Department of Claremont McKenna,Pitzer, and Scripps Colleges, Claremont, California 91711 Department of Physics, University of Connecticut, Storrs, Connecticut 06269 Department of Physics, Wesleyan University, Middletown, Connecticut 06459 Department of Sciences, Wentworth Institute of Technology, Boston, Massachusetts 02115 (Dated: October 6, 2018)We present experimental measurements of the steady-state ion number in a linear Paul trap(LPT) as a function of the ion-loading rate. These measurements, taken with (a) constant Paul trapstability parameter q , (b) constant radio-frequency (rf) amplitude, or (c) constant rf frequency, shownonlinear behavior. At the loading rates achieved in this experiment, a plot of the steady-state ionnumber as a function of loading rate has two regions: a monotonic rise (region I) followed by a plateau(region II). Also described are simulations and analytical theory which match the experimentalresults. Region I is caused by rf heating and is fundamentally due to the time dependence of therf Paul-trap forces. We show that the time-independent pseudopotential, frequently used in theanalytical investigation of trapping experiments, cannot explain region I, but explains the plateauin region II and can be used to predict the steady-state ion number in that region. An importantfeature of our experimental LPT is the existence of a radial cut-off ˆ R cut that limits the ion capacityof our LPT and features prominently in the analytical and numerical analysis of our LPT-loadingresults. We explain the dynamical origin of ˆ R cut and relate it to the chaos border of the fractal ofnon-escaping trajectories in our LPT. We also present an improved model of LPT ion-loading as afunction of time. PACS numbers: 37.10.Ty, 52.27.Jt, 52.50.Qt
I. INTRODUCTION
The loading dynamics of the linear Paul trap (LPT) areof interest in relation to recent work measuring the totalcharge exchange and elastic collision rate of atomic ionswith their parent atoms in a hybrid trap [1, 2]. In thesehybrid-trap measurements of collisions between atomsand dark ions, those without optically accessible tran-sitions, the fluorescence of the atoms in the MOT wasmonitored in the presence of trapped ions. Knowledge ofboth the number of trapped ions and the size of the ioncloud is necessary for finding the total interaction rate.Both of these papers attempted to create a model forthe loading of ions into a Paul trap. However, neithergroup was able to create a satisfactory model of Paultrap loading. This paper aims to understand and modelthe loading of a linear Paul trap at the conditions usedin current hybrid trap experiments.Hybrid traps consist of a neutral atom trap coincidentwith an ion trap. A variety of neutral atom traps havebeen used in hybrid trap experiments: magnetic traps [3–5], magneto-optical traps (MOT) [1, 2, 6–15], and opticaldipole traps [16–18]. In contrast, every hybrid trap ex-periment listed above used a Paul trap to hold the ionicspecies, except [10], which used an octupole trap.In hybrid trap experiments for dark ions and their par-ent atoms, the ion-atom total elastic and inelastic inter-action rate per atom γ ia is given by [2] γ ia = k ia N I CV ia . (1)In this equation k ia is the total elastic and charge- exchange collision rate constant, N I is the average num-ber of trapped ions overlapping the atom cloud, C is afunction that describes the concentricity of the atom andion clouds, and V ia is the effective overlap volume of thetwo clouds [2].In these experiments the ion trap is saturated, whichhas two benefits. The first is that this maximizes the rateand thus the experimental resolution of the rate constant, k ia , because the ion number and volume will be as largeas possible.The second benefit to a saturated ion trap is that theion number, the concentricity, and the overlap volumewill all be time independent on average. In the case ofdark ions, the ion number can only be measured destruc-tively; this measurement would be extremely difficult ifthe ion number were time dependent. For bright ions,those with optical transitions, these quantities could bemeasured in a time dependent way using the fluorescence.Typically, however, bright ion experiments have focusedon the charge-exchange rate constant and have not mea-sured the total rate. If the technique described in [1, 2]were used in conjunction with the methods for measuringthe charge exchange rate constant, then the elastic rateconstant could be found as well.Neither of the dark ion experiments developed a satis-factory model for predicting the steady-state populationof the saturated Paul trap as a function of the ionizationrate. In the model used in Lee, et al. [1], an ad hoc termwas required to prevent the ion number from becominginfinite at large photoionization intensities and γ ia wasproportional to the Paul trap loss rate (cid:96) . In Goodman, et al. [2], we found experimentally that γ ia was not a r X i v : . [ phy s i c s . a t o m - ph ] J a n proportional to (cid:96) . We also found, by considering how thenumber of atoms in the MOT depended on the photoion-ization intensity, that the equation for the ion number asa function of photoionization intensity is naturally finiteas the photoionization intensity tends to infinity.Ion loading in a Paul trap as a function of time istypically fit to the solution of dN ( τ ) dτ = Λ − (cid:96) N ( τ ) − (cid:96) N ( τ ) , (2)where N is the number of trapped ions, τ is the time, Λis the loading rate, (cid:96) is the one-body loss rate, and (cid:96) is the two-body loss rate. All quantities in (2) are in SIunits. The model (2) is used, e.g., in [1, 2, 19, 20], thoughsometimes (cid:96) is set to zero without significantly affectingthe fit. The loading rate Λ is constant and set by theion source. The one-body and two-body loss rates (cid:96) and (cid:96) are usually assumed to be independent of the loadingrate. However, in Goodman, et al. [2] we found thatnot to be the case. Additionally, while this model canfit well, as in [1, 2], in other cases it overshoots the rise,undershoots the knee of the curve, and overshoots theplateau, as can be seen in Fig. 1. This trend is presenteven in experiments that do not have a hybrid trap andtherefore do not load from a MOT (see, e.g., Ref. [19]).Neglecting the two-body loss mechanism in (2), i.e., for (cid:96) = 0, the traditional model from the solution of (2) pre-dicts that the steady-state ion number is directly propor-tional to the loading rate. In our previous work, Bl¨umel, et al. [21], simulations showed that the steady-state ionnumber depends on the loading rate non-monotonically.There are four regions in a plot of the steady-state ionnumber as a function of loading rate. At small loadingrates, the steady-state ion number increases monotoni-cally, but nonlinearly, in region I. In region II, the steady-state ion number plateaus. Region III is characterized bya dip, where the steady-state ion number decreases withincreasing loading rate. Finally, in region IV, the steady-state ion number once again increases monotonically withthe loading rate, but much faster than in region I, andwith a power that does not agree with the predictionsof the model defined in (2). This behavior was presentfor simulations of both the linear and three-dimensionalPaul trap geometries.Despite being in use for over 50 years [22] in manydisciplines, including mass spectrometry, biology, chem-istry, and physics, the fundamental loading behavior ofthe Paul trap is still unknown. This work examines theloading behavior of the Paul trap experimentally, analyt-ically, and with simulations in regions I and II, where theloading rates are achievable experimentally in every hy-brid trap system and many mass spectrometry systems.This paper is organized as follows: Section II describesthe experimental apparatus and technique. Section IIIdescribes the experimental results and Sections IV andV compare those results with the simulations and theanalytical models. We conclude in Section VI. Figure 1. (Color online) Fits of ion loading as a function oftime based on (2) overshoot the rise, undershoot the knee, andovershoot the plateau. This is the case regardless of whetherthe ion number is measured via fluorescence at low loadingrates (top, reproduced from Ref. [19] with permission fromThe Royal Society of Chemistry) or measured using a channelelectron multiplier at high loading rates from a hybrid trap byour group (bottom). Where the error bars (showing statisticalerrors) are not seen, they are smaller than the correspondingplot symbols.
II. APPARATUS AND METHOD
The apparatus used in this experiment has been de-scribed elsewhere [2], so we will only briefly describe ithere. A diagram of our system can be seen in Fig. 2. Ourhybrid trap uses a linear radio frequency (r.f.) Paul trapwith a segmented design as the ion trap and a sodiumMOT as the neutral trap. The MOT is a vapor cell de-sign where the background sodium vapor is generated bya getter source. Even with constant loading from the get-ters, the background pressure is below ≈ . × − Torr.There are two MOT transitions in sodium involving dif-ferent hyperfine levels [23]. We use the type-II MOTtransition to take advantage of the higher trapped atomnumber achieved using that transition. The MOT ismade by retro-reflecting the three trapping beams and yz x
405 nm PI Laser 589 nm MOT LaserLPTMOT Ion Cloud
Figure 2. (Color online) A diagram of the hybrid apparatusused in these experiments. The anti-Helmholtz coils, withcurrent directions shown by arrows, are located outside of thevacuum chamber (not shown). the repumper is obtained from a sideband put on the 589-nm beam using an electro-optic modulator prior to split-ting them into three beams. The anti-Helmholtz coils arelocated outside of the vacuum chamber.The MOT can be characterized in two ways: using aphotomultiplier tube (PMT) or using a CMOS camera.Both the PMT and the CMOS camera can be used tomeasure the trapped atom number, but the camera canalso be used to measure the MOT radius. The MOT has a1/ e -density radius of ≈ .
75 mm and holds on the order of10 atoms. The PMT yields temporal information aboutthe loading of the MOT.The ions are created by a two-step process: atomsare resonantly excited to the 3P state by the trappinglaser beams of the MOT, then they are ionized by a sec-ond laser beam at 405 nm. The first resonant step as-sures that the sample of loaded ions is pure, from thesame species as the MOT neutrals.The beam size is fixedat a 1 /e -intensity radius of 1 . . z = 48 . r = 19 mm and the vertical distance between the surface of the rodsis 8 . r e = 8 . r e /r = 0 .
9, slightly smaller than the ideal ratio of1.1468 according to Ref. [24].Our trap also deviates from the ideal Paul trappingpotentials because the central segment is longer than theoptimal length for a quadratic axial confining potential.Instead, the static axial confinement is a superposition ofa quadratic and a quartic potential. As shown in our ex-periments and by the results of our simulations reportedin Appendix B, trapping is possible even for electrodegeometries that result in trapping potentials that have astrong quartic component and thus deviate considerablyfrom the ideal quadrupole trapping potential. The price,as shown in Appendix B, is the emergence of determin-istic chaos, even on the single-ion level, which reduces,and ultimately defines, the trapping volume of the LPT.The number of ions in the Paul trap can be de-structively measured using a channel electron multiplier(CEM). By putting DC voltages on the endcaps in adipole configuration, the ions can be directed out of thetrap along the axial direction toward the CEM. The pro-grammed CEM extraction sequence captures a constantfraction of the number of ions trapped immediately be-fore extraction. This extraction efficiency is built intoour calibration [2], which thus allows us to determine thenumber of ions immediately before extraction. We use acustom LabVIEW virtual instrument (VI) to control theexperimental timing and to record the measurements.The basic experimental scheme is to allow the MOTto load to steady state from the background vapor in thepresence of the photoionizing (PI) beam, then the iontrapping potentials are turned on. This way the loadingrate from the MOT is constant during the entire time theion trap is being loaded. The ions are extracted after aset delay time of 0 .
001 s after the loading ends and theion number is measured by the CEM during a 1 s interro-gation time. Since we can only measure the ion numberdestructively, we take many measurements at differentloading times to create a time series of the ion numberas a function of loading time.It is impossible to completely remove the ion loss mech-anisms of the LPT, so the loading rate can never be mea-sured independently of the loss rate. Whatever the exactform of the differential equation for the loading of theLPT, the loss rate depends on the number of trappedions and the loading rate is independent of the numberof trapped ions. Therefore, to measure the loading rate,the ion number is measured for very short loading times( ≈ .
05 s or shorter), where the loading can be describedby [see (2)] dNdτ ≈ Λ . (3)An example of this is shown in Fig. 3. We see that in thelimit of small loading times the number of ions trappedis a linear function of loading time. I on S i gna l ( V ) Time (s)
Figure 3. (Color online) Ion signal as a function of time forshort loading times, which minimizes the effects from losses.The ion signal is proportional to the number of trapped ions;the slope is the loading rate from the MOT. Where the errorbars (showing statistical errors) are not seen, they are smallerthan the corresponding plot symbols.
To measure the steady-state ion number, the loadingtime is set to be long enough to ensure that the ion pop-ulation in the trap has reached equilibrium. At the low-est loading rates this time could be up to 3 minutes; atthe highest loading rates the trap is saturated in mil-liseconds. Several consecutive measurements of the ionnumber are taken at a single loading time, then averagedto find the steady-state ion number. The loading rateis controlled by changing the 405 nm PI beam intensityand the process is repeated until the accessible part ofthe steady-state ion number as a function of loading rateis mapped out. With no PI beam present, no ion signalis measured.The CEM measures a voltage that is proportional tothe ion number, but the exact proportionality depends onthe high-voltage gain applied to it. The most straightfor-ward way to calibrate the CEM voltage would be to usean optical method to count the ion population and recordthe corresponding CEM voltage. However, because thereare no optically accessible transitions in Na + , the num-ber of ions could not be measured directly. Furthermore,our CEM is designed to have a large bias current, whichallows it to detect large ion signals. Consequently, theCEM can only operate in an analog detection mode andcannot be calibrated using pulse counting methods likethose used in Ref. [25].Instead, we used the same indirect method that wehave used previously [2], which we will briefly describehere. The ion loading rate is measured two ways andthen compared. First, the method described in Fig. 3is used, with the assumption that the fraction of ionsmeasured by the CEM, whatever it may be, does not Figure 4. (Color online) Two-way measurement of the ionloading rate as a function of the intensity of the photo-ionizinglaser for the purpose of finding the calibration factor betweenthe ion signal and the number of ions. Measurements usingthe CEM (black squares) are compared to measurements ofthe atom loss rate from the MOT in the presence of the pho-toionization laser (solid red line) multiplied by a fitted scalingfactor, which is the reciprocal of the CEM calibration factor.Where the error bars (showing statistical errors) are not seen,they are smaller than the corresponding plot symbols. change. The second method compares the loading of theMOT as a function of time with and without the PI beampresent. The increase of the loss rate in the presence ofthe 405 nm beam is equivalent to the loading rate of thePaul trap, under the assumption that every ion createdfrom the MOT is trapped. This assumption is good whenthe MOT is smaller than the trapping volume of the Paultrap. To convert the PMT voltage to units of number-of-atoms, we use the two-level atom model to determine theexcited-state population of the type-II MOT. One wouldexpect this approximation to be especially poor in thecase of the type-II MOT, where multiple hyperfine levelsin the excited state play a role in the pumping transition.We find however that the two-level model used with amodified saturation intensity of 37 . , comparedto the theoretical saturation intensity of 13 . ,fits quite well.The two ion-loading rates are plotted and comparedusing a one-parameter fit, where the fitting parameter isthe number of volts per ion. As seen in Fig. 4, the fitshows good agreement with the data for all CEM biasvoltages. This justifies the assumptions required for thecalibration. III. EXPERIMENTAL RESULTS
As discussed in Sec. II, the electric potential in ourLPT is not a pure, ideal quadrupole potential, but hasa significant admixture of a quartic component. As dis-cussed further in Appendix B, the quartic componentin our potential leads to single-ion chaos, which ren-ders our trap unstable in the radial direction from about r = 5 mm on, where r = (cid:112) x + y is the radial dis-tance from the axis of the LPT. This might suggest thata proper description of our LPT is possible only in termsof a sum of quadratic and quartic terms. However, sinceup to r = 5 mm the quartic terms are small comparedto the quadratic component, a simpler description of ourLPT potential as a quadratic potential with a cut-off at r = ˆ R cut ≈ r < ˆ R cut , the electric potential of our LPT, written in SIunits, is a pure quadrupole potential, given by [11, 26]: φ ( (cid:126)r, τ ) = V rf cos(Ω τ ) (cid:18) x − y r (cid:19) + ηV end z ( z − x − y ) , (4)where (cid:126)r = ( x, y, z ) is the position vector of an ion inthe trap, τ is the time, V rf is the rf voltage applied tothe electrodes of the trap, Ω = 2 πf is the angular fre-quency of the applied rf voltage ( f is the lab frequency), r and z are defined in Sec. II, η = 0 . V end is the voltage applied tothe end-segments of the trap. For r > ˆ R cut the form ofthe potential is not needed, since ions are rapidly ejectedfrom the trap (see Appendix B) as soon as they crossthe chaos border at r = ˆ R cut . Measured with respectto ˆ R cut , and in pseudo-potential approximation (see Ap-pendix A), the depth of the LPT potential is given by D = (cid:20) e V m Ω r − eηV end z (cid:21) ˆ R . (5)Experimentally, the trap depth D can be changed bychanging the rf amplitude V rf , the angular rf frequency Ω(via its dependence on the lab frequency f = Ω / (2 π )), orthe end-cap potential V end . In the experiments reportedin this paper, V end = 30 V is kept constant and only V rf and f are varied. Equation (5) was derived for a sin-gle trapped ion; if a second ion is trapped, the Coulombrepulsion changes the effective trap depth each ion ex-periences. The anti-trapping space-charge effect frommultiple trapped ions makes the trap depth no longerexpressible analytically. However, as supported by ouranalysis of ion numbers in the saturation region II be-low, the functional dependence of the trap depth on V rf and f appears to be the same.It is convenient to express the trap depth in terms ofthe single-particle stability parameter q = 4 eV rf m Ω r (6) S t ead y S t a t e I on N u m be r ( ) Loading Rate (Ions/rf Cycle)
Figure 5. (Color online) Loading curves taken at constant q = 0 . V rf = 13 V and f = 450 kHz for the squares (black), V rf = 16V and f = 500 kHz for the circles (red), and V rf = 19 . f = 550 kHz for the triangles (blue). Where the errorbars (showing statistical errors) are not seen, they are smallerthan the corresponding plot symbols. so that D = (cid:20) qeV rf r − eηV end z (cid:21) ˆ R . (7)We also define the dimensionless loading rate λ = 2 π Λ / Ω , (8)i.e., the number of ions loaded per rf cycle.In Fig. 5, loading curves of the steady-state ion num-ber as a function of loading rate λ are shown for threesets of rf parameters that each result in q = 0 .
3. In eachcase, both the monotonic rise in region I and the plateauof region II are visible, predicted and previously observedin [21]. Since q is kept constant in all three cases shownin Fig. 5, the depth of the LPT potential in these threecases is most conveniently evaluated according to (7). Inaddition, (7) shows that, for constant q , D is independentof the frequency and depends only on V rf and geometricconstants of the LPT. Because of the quadratic form (4)of the LPT potential, the trapped ion cloud is an ellip-soid with semi-major axes equal to ˆ R cut in the x and y directions. In the z direction we have D = mω z ˆ Z ,where ˆ Z cut is the extent of the ion cloud in the z directionand ω z is the pseudo-oscillator frequency in the z direc-tion. Since ω z is determined by the static potential dueto the end-caps of the LPT, ω z is a constant and, there-fore, ˆ Z cut ∼ √ D . Thus, the volume of the trapped ioncloud is V = (4 π/
3) ˆ R ˆ Z cut ∼ ˆ R √ D . According toPoisson’s equation of electrostatics, the density ρ of thetrapped ions is proportional to the Laplacian of the trap-ping potential, i.e., according to (4), ρ ∼ ∇ φ ∼ V / Ω . S t ead y S t a t e I on N u m be r ( ) Loading Rate (Ions/rf Cycle)
Figure 6. (Color online) Loading curves taken at constant V rf = 16 V and plotted on a lin-log scale. The trap settingswere q = 0 .
22 and f = 580 kHz for the squares (black), q =0 .
26 and f = 535 kHz for the circles (red), q = 0 .
30 and f = 500 kHz for the upright triangles (blue), and q = 0 . f = 450 kHz for the inverted triangles (turquoise). Theplots are in order of their trap depths. Where the error bars(showing statistical errors) are not seen, they are smaller thanthe corresponding plot symbols. This is all we need to predict, up to a proportionalityconstant, the number, N , of stored particles in the LPT.On the basis of the above discussion we have N = ρV ∼ ˆ R V D / / Ω . (9)More details on the derivation of (9) can be found inSec. V B.We can now use (9) for a consistency check of our ex-perimental results in Fig. 5. Denoting by N , N , and N . , the particle numbers corresponding to the cases V rf = 13 V, V rf = 16 V, and V rf = 19 . N : N : N . , which we can predict on the basis of(9), even without knowledge of the proportionality con-stant in (9). Indeed, assuming that ˆ R cut depends onlyweakly on V rf and Ω, the cut-off radius ˆ R cut cancels whentaking ion-number ratios on the basis of (9). Thus, thecomputation of ratios involves only V rf , Ω, and known ge-ometric constants. This way we obtain the explicit, ana-lytical prediction N : N : N . = 1 : 1 .
63 : 2 .
48. Thisprediction may be compared with the actual ion num-bers in the plateau regimes read off from Fig. 5. With N ≈ , N ≈ , N . ≈ , N : N : N . = 1 : 1 .
69 : 2 .
44. Thisis in excellent agreement with the theoretical prediction.As mentioned above, taking ratios has the advantage ofeliminating ˆ R cut , which is difficult to obtain directly inour experiments, since the trapped Na + ions are dark.Region I and region II can likewise be seen in the S t ead y S t a t e I on N u m be r ( ) Loading Rate (Ions/rf Cycle)
Figure 7. (Color online) Loading curves taken at constant f = 500 kHz and plotted on a lin-log scale. The trap settingswere q = 0 .
22 and V rf = 12 V for the squares (black), q = 0 . V rf = 14 V for the circles (red), q = 0 .
30 and V rf = 16 Vfor the upright triangles (blue), and q = 0 .
37 and V rf = 20 Vfor the inverted triangles (turquoise). The plots are in orderof their trap depths. Where the error bars (showing statisticalerrors) are not seen, they are smaller than the correspondingplot symbols. curves shown in Fig. 6, which are taken at a constantrf amplitude, V rf = 16 V. Shown in Fig. 6 are four load-ing curves with f = 580 kHz, 535 kHz, 500 kHz, and450 kHz. Using the frequencies as labels, we read off N ≈ , N ≈ , N ≈ , N ≈ , N : N : N : N = 1 : 1 .
54 : 2 .
92 : 3 .
31. The theo-retical prediction for these ratios, according to (9), is N : N : N : N = 1 : 1 .
57 : 2 .
14 : 3 .
29. Justlike for the results shown in Fig. 5, and except for thecurve with f = 500 kHz, the experimental ratios matchthe theoretical predictions very well. At present it is notclear why the third curve, at f = 500 kHz, is an outlier inthis sequence. This is the more puzzling that this curveis the same as the corresponding curve shown in Fig. 5,where it fits the sequence in Fig. 5 very well. That thiscurve does not fit well in Fig. 6 is also immediately obvi-ous from the visual context in Fig. 6. This curve producesa gap in region II of Fig. 6, whereas a more even spacingwas expected, mirroring the small decrements in frequen-cies corresponding to the four curves shown in Fig. 6.Several loading curves were also taken at a constantrf frequency; they are shown in Fig. 7. Akin to Figs. 5and 6, the three lower curves all show both region I andregion II. Only the loading curve taken at the largest rfvoltage ( V rf = 20 V) does not look like it has reached sat-uration (region II) yet, an impression confirmed by ourratio test to be conducted next. In the case of Fig. 7there is a wrinkle in our theoretical analysis of the case V rf = 12 V in that for the experimental parameters thedensity (7) comes out negative, which means that wedo not obtain a real square root in (9) and therefore N cannot be predicted. This is not a disaster. It sim-ply means that for V rf = 12 V and f = 500 kHz, ourtrap is operated so closely to the global instability bor-der of the LPT that the pseudo-potential analysis (seeAppendix A) is not accurate enough in this borderlinecase to make accurate predictions. For our ratio test weopted to ignore this case and normalize to the secondcurve in Fig. 7, i.e., the case V rf = 14 V, f = 500 kHz,for which we obtain a positive trap depth of substan-tial magnitude for which our pseudo-potential analysis isvalid. Using voltages as labels, like we did in the caseof Fig. 5, we predict N : N : N = 1 : 1 .
95 : 4 . N ≈ , N ≈ , N ≈ , N : N : N = 1 : 2 .
00 : 3 .
47. Similarly tothe cases discussed in connection with Figs. 5 and 6, theion-number ratio of the two cases of Fig. 7, which are con-fidently in region II, is very close to its predicted value. Incontrast, the predicted ratio for N /N is much largerthan the experimentally observed ratio, confirming oursuspicion that for V rf = 20 V and at λ = 3 ions/rf-cyclethe loading curve in this case is still climbing (still in re-gion I), and its saturated ion number, expected to occurat higher values of λ than those shown in Fig. 7, willeventually be larger than N = 660 , D , perhapsby neglecting the term proportional to V end , in order toturn (9) into a more concise formula. Alas, as the case V rf = 12 V in Fig. 7 vividly illustrates, this is not possi-ble. The two terms, i.e., the terms involving V rf and V end in (5) [(7), respectively], are of similar magnitude, whichmakes it impossible to neglect one with respect to theother. Thus, the square-root behavior in (9) is essential.Our theoretical estimates and predictions above arebased on a single-ion picture (the pseudo-potential anal-ysis presented in Appendix A), and do not include anyspace-charge- or many-body effects. The overall excel-lent agreement of our predictions for the ion-number ra-tios in region II of the cases shown in Figs. 5, 6, and 7leads us to conclude that in our LPT experiments theseeffects are either negligible or lead to a simple renormal-ization of our expression for N that results in an overallconstant that cancels upon taking ratios. In Sec. V B,based on the single-particle pseudo-potential picture, wemake predictions of the absolute magnitude of N in re-gion II, which agree very well with the experimentallyobserved values. This might argue for the space-charge-and many-body effects to be negligible. However, sinceˆ R cut enters these formulas multiplicatively, we cannot besure whether the renormalization constant is not simplyabsorbed in our effective ˆ R cut , used in Sec. V B. Onlydirect experimental observation of ˆ R cut can resolve thisissue. This, however, due to the optical darkness of theNa + ions used in our experiments, is currently beyond our experimental capabilities. Nevertheless, the excel-lent agreement of the experimental ion-number ratios inregion II with our theoretical predictions supports thevalidity of our experimental LPT ion-loading curves. IV. SIMULATIONS
Model simulations have already been done [21] thatconfirm the existence of the four different dynamicalregimes. In addition it was shown in [21] that the phe-nomenon is robust with respect to the statistical dis-tribution of time between loading events, temperature,loading mechanisms, and the geometry of the absorbingboundary. So here the emphasis is not so much on prov-ing the existence of the four dynamical regions, or theirrobustness, but to see whether the simulations can qual-itatively, and to some extent quantitatively, describe theexperimentally observed characteristics of regions I andII.Since independence of the loading statistics has alreadybeen demonstrated in [21], we focus in this paper on thecase of uniform loading statistics. Qualitatively, our re-sults stay valid, which we checked explicitly, if differentloading statistics, such as Gaussian statistics, are used.The Newtonian equations of motion of a singly chargedion of mass m in the trap is m d (cid:126)rdτ = − e(cid:126) ∇ φ ( (cid:126)r, τ ) , (10)which, written out in components, results in m d dτ xyz = − V rf e cos(Ω τ ) xr + ηeV end z x V rf e cos(Ω τ ) yr + ηeV end z y − ηeV end z z . (11)Defining the dimensionless time t = (cid:18) Ω2 (cid:19) τ, (12)the set of equations may be written as ¨ x ¨ y ¨ z = − q cos(2 t ) x + bx q cos(2 t ) y + by − bz , (13)where the dots indicate differentiation with respect to di-mensionless time t , the dimensionless control parameter q is defined in (6), and b = 4 eηV end m Ω z . (14)If more than one ion of charge e and mass m are storedin the trap, the ions interact via the Coulomb force, re-sulting in the following set of coupled equations: ¨ x i + 2 q cos(2 t ) x i − bx i ¨ y i − q cos(2 t ) y i − by i ¨ z i + 2 bz i = N (cid:88) j =1 j (cid:54) = i (cid:126)r i − (cid:126)r j | (cid:126)r i − (cid:126)r j | , (15)where i = 1 , . . . , N counts the number of particles in thetrap at time t , and (cid:126)r is measured in units of l = (cid:18) e π(cid:15) m Ω (cid:19) / , (16)where (cid:15) is the electric permittivity of the vacuum.It is the introduction of the unit of length (16) thatallows us to normalize the coefficient in front of theCoulomb force on the right-hand side of (15) to 1, andthus arrive at the set of equations (15) that depends onlyon the two dimensionless parameters q and b . q is an ad-justable control parameter that, depending on the trapvoltage V rf and frequency f , may be set to a value inthe interval 0 < q (cid:46) .
9, where q ≈ . q , in order to achieve trapping, the parameter b hasto satisfy 0 < b < q / N of trappedions, the coupled set of equations (15) is well conditioned.Therefore, a standard 4th order Runge-Kutta method[28] is enough to reliably integrate (15).While our numerical model (15) captures the essenceof the experimental LPT, and its parameters q and b areadjusted to their experimental values, our model is never-theless an idealization. The electrodes in our experimentare not hyperbolic surfaces, as required if (15) is expectedto be exact, and while (15) assumes a quadratic poten-tial in z direction, the potential in our LPT also has aquartic component [11, 26]. Still, the proportions of thenumerical trap, as expressed in (15), are correct and weexpect that our model captures the essential parts of thephysics in our experimental LPT.A final comment concerns the number of particles weare able to simulate compared with the number of par-ticles in our experimental trap. In order to accumulatesufficient statistics, and given our computer resources, wefound that 2,000 simultaneously stored ions are a practi-cal upper limit for our numerical simulations. Althoughorders of magnitude smaller than the experimental num-ber of particles in our trap, 2,000 particles is not a smallnumber, and using scaling relations, to be discussed be-low, we are able to compare our simulations not just qual-itatively, but also quantitatively with our experimentalresults.Our simulations proceed in the following way. Fora given loading rate λ we generate a time sequence t j , j = 1 , ..., M of loading events, which have a Poisso-nian distribution whose average corresponds to the spec-ified loading rate λ . For each individual parameter set-ting we check that M is large enough so that we aredeeply in the saturated regime where we are able to ex-tract the steady-state number of ions, N s , with excellentstatistics. Denoting by (cid:104) . . . (cid:105) t the time average in thesaturated regime, we also compute the statistical spread∆ N s = (cid:112) (cid:104) N s (cid:105) t − (cid:104) N s (cid:105) t , which characterizes the ion-number fluctuations in the saturated regime due to the Poissonian loading process. At each loading event t j anew ion with zero initial velocity is created in the trapat a random location inside of a spherical loading zoneof radius ˆ R load (in SI units), representing the creation ofions from the MOT via the photoionizing 405 nm laser.Between loading events, i.e. for t in t j < t < t j +1 , weintegrate the ion trajectories in the trap, including thenewly created ion, according to the system (15). Once t j +1 is reached, we eliminate all ions from the trap whosepositions at t j +1 lie beyond a pre-specified absorbingboundary B. In [21] we already showed that the quali-tative shape of the N s ( λ ) curves does not depend on thegeometry of B. Therefore, making use of this freedom,we chose in this paper a cylindrical absorbing boundaryB with radius ˆ R cut in the x - y direction and total length2 z in the z direction, i.e., ions are absorbed in the z direction if they exceed | z | = ˆ Z cut = z = 24 . q , all we need for our simulations arethe parameters b , R load = ˆ R load /l , R cut = ˆ R cut /l , and Z cut = ˆ Z cut /l . The radius of the loading zone ˆ R load isdefined by the size of the type-II MOT, which, accordingto Ref. [2], is r a = 0 .
75 mm. The electrodes are posi-tioned at r = 9 . R cut . According to the discussion above, we cannot sim-ulate the full-sized experimental LPT, since it typicallyholds more ions than we are able to realistically simu-late. Accordingly, we simulate a scaled-down version ofour LPT in which all linear dimensions are scaled by afactor 0 < σ <
1, resulting in R load = σ ˆ R load /l = 88 . σ [ f (MHz)] / ,R cut = σ ˆ R cut /l = 1119 . σ [ f (MHz)] / ,Z cut = σ ˆ Z cut /l = 2827 . σ [ f (MHz)] / . (17)As discussed in Sec. III, in a harmonic trap, i.e., a trapwith time-independent quadratic trapping potentials inall three directions, and at zero temperature, the chargedensity ρ el in the trap is constant, which follows imme-diately from ρ el ∼ ∇ φ . In this case the particle numberin the trap scales with the volume of the trap, i.e., N ∼ σ . (18)Although our experimental trap is not exactly harmonic,and the temperature is finite, we nevertheless expect that(18) holds to a good approximation. Therefore, in orderto compare with our experimental results, we scale N s obtained from our simulations according to N scaled s = N simulated s /σ , (19)which allows a direct comparison between the results ofour simulations with experimental results of the satu-rated number of ions in the trap.We are now ready for the numerical simulations.Since they are expensive, we focused on simulating theconstant- q experiments, performed with q = 0 .
3, and de-scribed in Sec. III. Figure 8 shows the results of our sim-ulations. I on N u m be r Loading Rate (Ions/rf cycle)
Figure 8. (Color online) Simulation results for the modelLPT performed at constant q = 0 . f = 450 kHz, V = 13 V:open, blue circles; f = 500 kHz, V = 16 V: filled, green cir-cles; f = 550 kHz, V = 19 . σ = 1 /
40. These simulation results may be compared with thethree corresponding LPT experiments shown in Fig. 5. Theheavy solid line (purple) is the analytical result for N s ( λ ) inregion IV [21]. The lengths of the error bars, equal to 2∆ N s ,characterize the statistical fluctuations of the ion number inthe saturated regime. Compared with the experimental results shown inFig. 5 we see that in our simulations the maximum ofregion II occurs at a loading rate which is about a factor10 lower than in the experiments. However, the locationof the region-II maximum depends on the scaling factor σ and shifts to higher loading rates as σ approaches 1where the simulated ion trap size becomes identical tothe experimental Paul trap. This effect can be seen inFig. 9, where simulations are plotted at several differentvalues of σ , along with several experimental curves. Ad-ditionally, the simulations in Fig. 9 have been scaled upin ion number [see (19)]. The simulated ion curves matchthe experiments in the ion number, but the regions occurat different ion loading rates. The reason for this discrep-ancy is not clear. It is possible that effects not included inthe idealizations made for the simulations, such as strayfields causing excess micromotion or the presence of thegrounded vacuum chamber, are the cause. V. ANALYTICAL THEORY
In this section we present an analytical theory for re-gions I and II. In Sec. V.A we show that in region I thesaturated ion number N s ( λ ) follows a power law in λ ,which we confirm experimentally. We also compute theapproximate exponent of the power law, which agrees I on N u m be r ( ) Loading Rate (Ions/rf cycle)
Figure 9. (Color online) Simulations of loading curves at threescaling factors: σ = 1 /
10 small triangles (blue), σ = 1 / σ = 1 / q = 0 . V rf = 13 V, and f = 450 kHz. The scatter of the experimen-tal data points gives an idea of the statistical and systematicvariations in our experimental loading data. For these sim-ulations, ˆ R cut = 4 . well with our experiments. In Sec. V.B we present atheory for region II. This theory explains the plateau be-havior of N s ( λ ) in region II and also fits the temporalbehavior of N s ( λ, t ) better than all other theories so fardescribed in the literature. A. Region I
In this subsection we present a simple, analyticallysolvable model for the steady-state ion population N s ( λ )as a function of loading rate λ . Our model predictsmonotonic power law behavior in region I. Under cer-tain reasonable assumptions, based on heating ratesobtained from molecular-dynamics simulations of non-neutral plasmas published in the literature [29], the ex-ponent derived from our analytical model is close to theexponent observed in our experiments. Since our ana-lytical calculations assume spherical symmetry, our ana-lytical results for region I are primarily applicable to thethree-dimensional quadrupole Paul trap (3DPT). Com-paring with experiment, we found that these results alsohold well for the LPT.In region I we need to consider the stationary state inwhich for each particle loaded one particle escapes. Inthe stationary state the spatial probability distribution,0 ρ ( (cid:126)r ), of the ions in the trap is approximately Gaussianwith a width that is proportional to √ T , where T is thetemperature. Since, according to E ∼ kT , where k is theBoltzmann constant, the average energy E of a stored ionis proportional to T , the width of the spatial Gaussian isproportional to √ E .In order to obtain analytical results in closed form, wereplace the Gaussian distribution by a flat distributionwith a sharp cut-off, i.e. we represent the spatial density ρ ( (cid:126)r ) by a homogeneous sphere according to ρ ( (cid:126)r ) = N s πw , | (cid:126)r | ≤ w, , | (cid:126)r | > w, (20)where N s is the number of ions in steady state, w = α √ E (21)is the radial width of the density distribution, and α isa constant. Since a Gaussian is a steeply descendingfunction in the wings, the approximation (20) is benignand does not change the scaling of N s ( λ ) in λ .We start with the situation in which the probabilitysphere, due to rf heating [29–31], has expanded beyondthe location R cut of the absorbing boundary B, just farenough for the total excess probability beyond R cut tointegrate to 1 particle. Denote by κ the heating rate perparticle. Then, the energy E per particle is E = E + κ ∆ τ, (22)where ∆ τ is the time that has passed since the last load-ing (ion creation) event and E is the energy per particleimmediately after the last loading event. Since we arein the steady state, exactly ∆ τ = 1 / Λ has passed on av-erage between ion creation and ion loss, where Λ is theloading rate. Therefore, from (22), E = E + κ Λ , (23)and, according to (21), w = α (cid:16) E + κ Λ (cid:17) / . (24)The excess width is δw = w − R cut = α (cid:16) E + κ Λ (cid:17) / − R cut . (25)In steady state, the excess width δw corresponds to ex-actly 1 particle. Therefore, denoting by δV the volumeof the shell of width δw :1 = ρδV = ρ πR δw = (cid:18) N s πR (cid:19) πR (cid:20) α (cid:16) E + κ Λ (cid:17) / − R cut (cid:21) = 3 N s R cut (cid:20) α (cid:16) E + κ Λ (cid:17) / − R cut (cid:21) . (26) In region I, i.e., for small Λ, the dominant term in (26) is κ/ Λ. Therefore, in region I, we may approximately write:1 = 3 N s αR cut (cid:16) κ Λ (cid:17) / . (27)In order to compute the dependence of N s on Λ, we needto know how the heating rate κ depends on N s . In orderto answer this question we turn to Fig. 2 of [29]. Thisfigure shows the heating rate H of N -ion clouds in steadystate as a function of cloud size ˆ s . This figure is relevantsince, up to the choice of units, κ and H are identical.Therefore, the N -scaling of κ is the same as the N -scalingof H . Although Fig. 2 of Ref. [29] was computed fora 3DPT, it nevertheless gives us a first idea on the N s scaling of κ for the LPT that we focus on in this paper.Since our trap has an effective radius R cut , we need toextract heating rates H from Fig. 2 of [29] as a functionof N at constant cloud size ˆ s . The most striking featureof the heating data shown in Fig. 2 of [29] is that theheating rate curves for different N are parallel to eachother and have about the same spacing when doublingthe number of particles N . Therefore, since Fig. 2 of [29]shows H on a log scale, both features combined showthat, at given ˆ s , independently of ˆ s , H follows a powerlaw in N . On the basis of the data displayed in Fig. 2 of[29], we find, at constant ˆ s : H ∼ N / . (28)Therefore, because of κ ∼ H , we obtain: κ = βN / s , (29)where β is a constant. Using this result in (27), we obtain1 = 3 αβ / N / s R cut Λ / . (30)Using (8), we solve (30) for N s in terms of λ : N s = (cid:20) R cut α (2 πβ/ Ω) / (cid:21) / λ / . (31)Thus, this simple model predicts N s ∼ λ . , whichmay be compared with the experimental region-I result N s ∼ λ . . Since the numerical value of the exponentpredicted by our model depends on the scaling of theheating rate κ in N , which was not separately determinedfor our LPT, the value of the exponent predicted by (31)is less important than the prediction that N s follows apower law. Therefore, we treat the value of the exponentas a fit parameter. Fits of the data from Fig. 5 using thismodel can be seen in Fig 10. The difference in the valueof the exponent is likely due to the difference in heatingrate between the ideal 3DPT of [29] and the heating ratein the experimental LPT.Since rf heating is one of the central ingredients in thismodel, both the prediction of a power law in itself and1 S t ead y S t a t e I on N u m be r ( ) Loading Rate (Ions/rf cycle)
Figure 10. (Color online) The data points from Fig. 5 shownwith fit using N s = Aλ . . We see that the power-law formof the fit function, as predicted by our analytical model, fitsthe data in region I very well. Where the error bars (show-ing statistical errors) are not seen, they are smaller than thecorresponding plot symbols. the approximate agreement of the power law exponentwith our experimental results indicate that rf heatingis the factor that governs the behavior of N s ( λ ) in re-gion I. This is corroborated by Fig. 11, which shows acomparison between N s obtained as a result of solvingthe fully time-dependent set of equations of motion (15),i.e., the equations of motion including rf heating (aster-isks in Fig. 11) and N s obtained as a result of solving thetime-independent pseudo-oscillator equations of motion(50), which, due to the lack of explicit time-dependence,are not capable of simulating rf heating (filled circles inFig. 11). Clearly, while the pseudo-oscillator model is ca-pable of reproducing regions II, III, and IV, it completelyfails to reproduce region I, which can only be attributedto a lack of rf heating in the pseudo-oscillator model,since otherwise this model contains all many-body forcesexactly as in the full set of equations (15). Thus, weproved conclusively that it is rf heating that determinesboth the power law and the power law’s exponent in re-gion I. B. Region II
In region II, the steady-state ion number plateau isdetermined by the depth of the trap. The pseudopoten-tial approximation, described in Appendix A, is valid inregion II. The pseudopotential well is completely filledwith ions in this region. The expression for the depth ofthe pseudopotential predicts a steady-state ion number I on N u m be r Loading Rate (Ions/rf cycle)
Figure 11. (Color online) Simulation of the LPT: Com-parison between the time-dependent model with rf switchedon (asterisks) and the pseudo-potential model with only thepseudo-potential present (filled circles). We simulated thecase V rf = 16 V, f = 450 kHz, with scaling factor σ = 1 / N s , characterizethe statistical fluctuations of the ion number in the saturatedregime. Where error bars are not seen, they are smaller thanthe plot symbols. in region II of N IIs = 4 π (cid:15) m × e f [Mhz] × ˆ R cut [mm] q (cid:114) q b − . (32)There are no adjustable parameters in this expression.However, ˆ R cut is not easily determined for dark ions. Thetraditional derivations of the pseudopotential depth as-sume that ˆ R cut = r and that the ions are only lostwhen they collide with or move beyond the trap elec-trodes [32, 33]. This has been shown not to be the casein recent hybrid trap work [2, 9]. Finding ˆ R cut directlyallows the trap depth to be calculated without makingany approximations. In principle, ˆ R cut can be found ina variety of ways.The first is to follow the method used in [2, 9], wherethe idealized single-particle trap depth is equated to theenergy of a simple harmonic oscillator with spring con-stant k = mω , where ω is the secular frequency of theionic motion. This method has the advantage of be-ing simple, but the drawback is that it relies on thesingle-particle trap depth, which is a credible approxi-mation, but may not be accurate enough. Indeed, thissimple model disagreed with subsequent hybrid trap ex-periments, described in [2], by 25%.A second method, mentioned in the caption of Fig. 9,is to find the ˆ R cut value which makes the simulationsbest match the experimental results. This requires some2computationally demanding simulations of the trappedions.A third method is to find an ˆ R cut that makes (32) fitbest at one trap setting. For sodium ions, (32) becomes N IIs = 548344 f [Mhz] ˆ R cut [mm] q (cid:114) q b − . (33)By selecting a single data set and finding the ˆ R cut thatmakes both sides of (33) approximately equal, we are ableto fit nearly all of our data at least as well as the othermethods, and in some cases much better, using a muchsimpler procedure (see Table I). The data set with thesmallest ion number, which was taken at q = 0 . V rf =16 V, and f = 500 kHz, is the one exception; using theseparameters in this model returns an imaginary numberof ions. These settings also have the smallest numberof ions at steady state in region II, i.e., approximately80000, which is still a large number for most Paul trapexperiments. It is also unclear what differentiates thesettings where this value of ˆ R cut fits well and the onesthat do not.There has been no previous study of why ˆ R cut (cid:54) = r ,but it has been hypothesized that it could be due to con-tributions from higher-order multipoles [34]. In this pa-per we confirm that it is the admixture of the quarticcomponent of the confining potential in z direction thatis responsible for the reduction of ˆ R cut to ˆ R cut < r . Infact, while a single ion in an ideal trap, i.e., a trap inwhich only quadrupole potentials are present, is neverchaotic, even when driven by the rf trap fields, a singleion in a potential with a quartic admixture shows a tran-sition in space from regular, confined motion, to chaotic,unconfined motion (see Appendix B). This brings us to afourth method for determining ˆ R cut , described in detailin Appendix B. According to this method, ˆ R cut is identi-cal with the single-ion chaos border. This means that for r < ˆ R cut the ion’s motion is perfectly regular and the ion,in the absence of noise, is perfectly trapped. However, assoon as r exceeds ˆ R cut , the ion’s motion becomes chaotic.As a consequence, the ion is free to explore spatial regionswith r > ˆ R cut , which quickly leads to an encounter withthe electrodes at which the ion is absorbed. Thus, ˆ R cut is determined by a method that is based on the dynam-ics of a single ion and therefore allows a very quick andefficient determination of ˆ R cut for various trap settings.This method does not only have technical advantages forthe determination of the value of ˆ R cut . It also solves thepuzzle of the very existence of ˆ R cut , identifying its originas a fundamental, purely dynamical effect, a chaos transi-tion, whose exact location is determined by the strengthof the admixture of higher multipole fields to the LPT’squadrupole trapping field.The loading of the Paul trap as a function of time wasalso examined. In our previous work [21], the loading asa function of time in region IV was determined from the Table I. The experimental ion number compared with the ionnumber predicted by (33) with ˆ R cut = 3 .
74 mm. The cellswith N/A indicate an imaginary number predicted. q f[MHz] b PredictedIon Number MeasuredIon Number PercentDifference0.30 0.550 0.022 574273 566162 1.40.30 0.500 0.026 387308 381897 1.40.30 0.450 0.032 231851 229396 1.10.22 0.500 0.026 N/A 81654 N/A0.26 0.500 0.026 185254 198434 6.60.37 0.500 0.026 882666 672169 31.30.37 0.450 0.032 594951 455767 30.50.26 0.535 0.023 271704 205591 32.20.22 0.580 0.019 163397 132166 23.6 equation t = − R λ / (cid:26) ln ( α − ˜ N / ) −
12 ln ( ˜ N + α ˜ N / + α )+ √ (cid:32) N / + αα √ (cid:33) − √ (cid:18) √ (cid:19)(cid:41) , (34)where α = ˜ R / ˜ λ / , and the tilde indicates that thequantity corresponds to the loading zone, which in thiscase is the MOT volume. Therefore, for example, ˜ R isthe radius of the loading zone, which in this case is theMOT radius.Equation (34) was derived for region IV and one shouldnot expect it to fit in regions I and II. Indeed, applying itto plots at extremely small loading rates results in poorfits and the prediction of imaginary ion numbers. How-ever, in region II, (34) fits the loading data at least aswell, if not better, than the traditional model (2). Thisincludes versions of the traditional model where (cid:96) = 0or (cid:96) = 0, as can be seen in Fig. 12. All forms of thetraditional model demonstrate the overshoot of the rise,undershoot of the knee, and overshoot of the plateau seenalso in Fig. 1. While (34) overshoots the rise as well, itfits the knee and the plateau better than the traditionalmodel. Further, unlike the traditional model, it is de-rived from the underlying physical process and dependson quantities that can be experimentally measured inde-pendently of the model itself. VI. CONCLUSION
In this work, two of the four regions in the load-ing curve of a Paul trap, predicted in [21], have beenconfirmed experimentally. Additionally, both simula-tions and analytical models match reasonably well to thesteady-state ion number and shape of the curves. Thetwo regions are relevant and accessible to the majorityof Paul trap experiments. Also proposed are new, exper-imental and computational methods for finding ˆ R cut , a3 I on S i g . ( V ) Time (s)
Figure 12. (Color Online) Ion number vs. time fitted tothe model (34) (blue), the traditional model (2) (green), thatmodel with (cid:96) = 0 (purple), and that model with (cid:96) = 0(red). The traditional model displays the familiar overshoot-undershoot-overshoot pattern, but the model (34) fits boththe knee and the plateau. Where the error bars (showingstatistical errors) are not seen, they are smaller than the cor-responding plot symbols. quantity necessary for finding the actual trap depth (adifficult experimental prospect), the total collision rateconstant for dark ions [1, 2], and for describing the be-havior of the Paul-trap loading in region III [21]. Furtherwork remains to be done to fully understand the loadingof the Paul trap. An upcoming paper will examine thebehavior of regions III and IV through experiment, sim-ulations, and analytical theory. VII. ACKNOWLEDGEMENTS
W.W.S. would like to acknowledge NSF support (inpart) from grant 1307874.
APPENDIX A: LPT PSEUDOPOTENTIAL
In this appendix we derive the single-particle pseu-dopotential for the LPT used in the analysis of our ex-periments. The pseudopotential is used in our numericalsimulations to prove (a) that the power-law behavior ofthe loading curves in region I is a consequence of rf heat-ing, which we prove in reverse by demonstrating thatthe pseudo-potential equations of motion, lacking an rfterm, cannot explain region I, and (b) that the existenceof region II does not depend on the time-dependence ofthe rf drive of the trap, i.e., as demonstrated in Sec. V B,the time-independent pseudopotential alone is capable ofexplaining the plateau in region II. We start with the equations of motion (13), where x , y ,and z are in units of l [see (16)], t is in units of 2 / Ω [see(12)], and q , b are the dimensionless control parametersdefined in (6) and (14), respectively. Focusing on the x component of (13), and following the procedure outlinedin [35], we split the x coordinate into a large-amplitude,slowly varying component X ( t ), called the macromotion,and a small-amplitude, fast-oscillating component ξ ( t ),called the micromotion according to x ( t ) = X ( t ) + ξ ( t ) . (35)Defining the cycle average (cid:104) f ( t ) (cid:105) = 1 π (cid:90) π/ − π/ f ( t + t (cid:48) ) dt (cid:48) , (36)and in line with the physical meanings of X and ξ , weassume (cid:104) X ( t ) (cid:105) = X ( t ) , (cid:104) ξ ( t ) (cid:105) = 0 , (cid:104) ¨ X ( t ) (cid:105) = ¨ X ( t ) , (cid:104) ¨ ξ ( t ) (cid:105) = 0 . (37)Focusing first on the time-dependent (rf) part of (13) andusing the decomposition (35), we have¨ X ( t ) + ¨ ξ ( t ) = − q cos(2 t )[ X ( t ) + ξ ( t )] . (38)Because ¨ ξ ( t ) dominates the left-hand side of (38) and X ( t ) dominates the right-hand side, we may write ap-proximately ¨ ξ ( t ) = − q cos(2 t ) X ( t ) . (39)Since, according to (37), X ( t ) is assumed to be constantover one rf cycle, we may integrate (39) immediately,resulting in ξ ( t ) = q t ) X ( t ) , (40)where we set the integration constants to zero. This isnecessary for consistency, since these constants lead tonon-oscillating, slow terms that are assumed to be con-tained in X ( t ). To compute X ( t ), we take the cycle av-erage of (38). Assuming that X ( t ) and ξ ( t ) are uncor-related, i.e., (cid:104) X ( t ) ξ ( t ) (cid:105) = 0, and (cid:104) cos(2 t ) X ( t ) (cid:105) = 0, wearrive at¨ X ( t ) = − q (cid:104) cos(2 t ) ξ ( t ) (cid:105) = − q (cid:104) q (2 t ) X ( t ) (cid:105) = − q X ( t ) , (41)where we used (cid:104) cos (2 t ) (cid:105) = 1 /
2. This equation of motionfor X ( t ) may be derived from the potential U eff ( X ) = q X (42)via ¨ X ( t ) = − ∂U eff ( X ) ∂X . (43)4To obtain U eff ( X ) in SI units, we multiply (42) with theunit of energy E = ml Ω . (44)The force proportional to b in (13) may be derived fromthe potential U stat ( X ) = − b X . (45)Notice that this potential is deconfining. Combining (42)and (45) results in the total pseudopotential U pp ( X ) = U eff ( X ) + U stat ( X ) = (cid:18) q − b (cid:19) X (46)acting on the macromotion coordinate X of an ion in theLPT. Since the y equation of (13) is formally identicalwith the x equation, we obtain immediately: U pp ( Y ) = (cid:18) q − b (cid:19) Y , (47)where Y is the macromotion coordinate of a trapped ionin y direction. To obtain the pseudopotential for the z coordinate of a trapped ion, all we need to do is to set q to zero and replace b → − b in the above derivations toobtain U pp ( Z ) = bZ , (48)where Z is the macromotion coordinate of an ion in z direction. Clearly, in order to achieve trapping in the x and y directions, we need the pseudo-oscillator potentialsin x and y directions to be confining, which requires thecoefficients in front of the X and Y terms in (A12) and(A13) to be positive, which, in turn, requires b < q / z direction, we need b in (A14) to be positive. Combining these two conditions,we obtain the condition0 < b < q / ¨ x i + q x i − bx i ¨ y i + q y i − by i ¨ z i + 2 bz i = N (cid:88) j =1 j (cid:54) = i (cid:126)r i − (cid:126)r j | (cid:126)r i − (cid:126)r j | . (50)Notice the change of sign in the y -equation part of (50)with respect to (15), which is consistent, since the rf field,on average, produces a confining force in y direction,which, in the pseudopotential equations (50), requiresa “+” sign in front of the q y i term. APPENDIX B: THE DYNAMICAL ORIGIN OF ˆ R cut In this appendix we show that ˆ R cut has a purely dy-namical origin. It is explained as a chaos border due tothe quartic admixture in the z potential P ( z ) of the trap.A fit of P ( z ) on the axis of the trap yields P ( z ) = 3 . × − z − . × − z +3 . × − z + 4 . × − z + 0 . , (51)where z is in mm and P ( z ) is in volts. Assuming cylindri-cal symmetry and neglecting the small terms asymmetricin z proportional to z and to z , we extend P ( z ) into the x and y directions, i.e., P ( z ) → P ( x, y, z ), by requiring ∇ P ( x, y, z ) = 0. We obtain P ( x, y, z ) = 0 . V + 1 .
986 V z (cid:18) z − r (cid:19) + 11 .
082 V z (cid:18) z + 38 r − z r (cid:19) , (52)where z = 24 . P ( z ) → P ( x, y, z )is unique, once cylindrical symmetry is assumed. We areaware of the fact that cylindrical symmetry can be trueonly close to the LPT’s axis, since closer to the rods, wehave a four-fold symmetry, which breaks rotational in-variance around the LPT’s z axis. However, close to theLPT’s axis, (52) is an acceptable analytical approxima-tion, which, for r (cid:46) P ( x, y, z ) by less than 30%.On the basis of (52) we obtain the following single-ionequations of motion ¨ x + 2 q cos(2 t ) x − b x + b (cid:2) ( x + xy ) − xz (cid:3) ¨ y − q cos(2 t ) y − b y + b (cid:2) ( y + yx ) − yz (cid:3) ¨ z + 2 b z + b (cid:2) z − z ( x + y ) (cid:3) = 0 , (53)where b = 1 .
986 eV mπ z f , b = 11 .
082 eV l mπ z f (54)and l is defined in (16). Integrating the system ofequations (53) for many initial conditions, we foundthat the single-ion dynamics governed by (53) exhibitstrapped and escaping trajectories. We illustrate thisin the following way. For q = 0 . f = 450 kHz (thecase shown in Fig. 9), we determined the lifetimes L (in rf cycles) of 72,000 trajectories with initial condi-tions x n = − n × . n = 1 , . . . , y = 0,and z m = −
20 mm + m × . m = 1 , . . . , x direction, while trajecto-ries that visit any of the colored regions quickly escape.Therefore, from Fig. 13, we conclude that ˆ R cut ≈ Figure 13. (Color online) Color-coded fractal of escape times.Denoting by L the number of rf cycles it takes for an initialcondition ( x, z ) to reach | x | ≥ r = 9 . L < ≤ L < ≤ L < ≤ L < ,
000 (blue), and
L > ,
000 (black). The black areacorresponds to initial conditions that lead to ion trajectoriesthat never escape. The black area is bounded in x directionby ˆ R cut ≈ x = 0, protruding fromthe fractal, are also shown in black, since they correspond toon-axis equilibrium points that, too, never escape. This is consistent with the value ˆ R cut ≈ . R cut is identified asa chaos border. Therefore, far from caused by any non-controllable effects, such as patch fields, stray fields, ornoise (although these effects certainly may modify ˆ R cut ),the reduced trapping capacity of our LPT, characterizedby ˆ R cut , is a purely deterministic, dynamic effect, whichis fundamentally related to the shape of the trapping po-tential of our LPT in z direction. While the investigationof the properties of the escape fractal shown in Fig. 13 isan interesting project in itself, it is beyond the scope ofthis paper and not necessary for the purpose of explain-ing the dynamical origin of ˆ R cut . We will report moreresults on the escape fractal, including chaos and orderin our LPT, elsewhere. [1] S. Lee, K. Ravi, and S. A. Rangwala, Phys. Rev. A ,052701 (2013).[2] D. S. Goodman, J. E. Wells, J. M. Kwolek, R. Bl¨umel,F. A. Narducci, and W. W. Smith, Phys. Rev. A ,012709 (2015), 10.1103/PhysRevA.91.012709.[3] C. Zipkes, S. Palzer, C. Sias, and M. K¨ohl, Nature ,388 (2010).[4] C. Zipkes, S. Palzer, L. Ratschbacher, C. Sias,and M. K¨ohl, Phys. Rev. Lett. , 133201 (2010),10.1103/PhysRevLett.105.133201.[5] S. Schmid, A. H¨arter, and J. H. Denschlag,Phys. Rev. Lett. , 133202 (2010), 10.1103/Phys-RevLett.105.133202.[6] A. Grier, M. Cetina, F. Oruˇcevi´c, and V. Vuleti´c,Phys. Rev. Lett. , 223201 (2009), 10.1103/Phys-RevLett.102.223201.[7] W. G. Rellergert, S. T. Sullivan, S. Kotochigova,A. Petrov, K. Chen, S. J. Schowalter, and E. R. Hud-son, Phys. Rev. Lett. , 243201 (2011), 10.1103/Phys-RevLett.107.243201.[8] F. H. J. Hall, M. Aymar, N. Bouloufa-Maafa, O. Dulieu,and S. Willitsch, Phys. Rev. Lett. , 243202 (2011).[9] K. Ravi, S. Lee, A. Sharma, G. Werth, and S. Rangwala,Nature Communications , 1126 (2012).[10] J. Deiglmayr, A. G¨oritz, T. Best, M. Weidem¨uller, andR. Wester, Phys. Rev. A , 043438 (2012).[11] I. Sivarajah, D. S. Goodman, J. E. Wells, F. A. Narducci,and W. W. Smith, Phys. Rev. A , 063419 (2012).[12] F. H. Hall, M. Aymar, M. Raoult, O. Dulieu, and S. Willitsch, Molecular Physics , 1683 (2013).[13] S. Haze, S. Hata, M. Fujinaga, and T. Mukaiyama, Phys.Rev. A , 052715 (2013), 10.1103/PhysRevA.87.052715.[14] W. W. Smith, D. S. Goodman, I. Sivarajah, J. E. Wells,S. Banerjee, R. C¨ot´e, H. H. Michels, J. A. Mongtomery,and F. A. Narducci, Applied Physics B , 75 (2014).[15] S. Haze, R. Saito, M. Fujinaga, and T. Mukaiyama,Phys. Rev. A , 032709 (2015), 10.1103/Phys-RevA.91.032709.[16] A. H¨arter, A. Kr¨ukow, A. Brunner, W. Schnitzler,S. Schmid, and J. H. Denschlag, Phys. Rev. Lett. ,123201 (2012), 10.1103/PhysRevLett.109.123201.[17] A. Kr¨ukow, A. Mohammadi, A. H¨arter, and J. H.Denschlag, arXiv:1602.01381 [physics] (2016), arXiv:1602.01381.[18] Z. Meir, T. Sikorsky, R. Ben-shlomi, N. Akerman,Y. Dallal, and R. Ozeri, arXiv:1603.01810 [cond-mat, physics:physics, physics:quant-ph] (2016), arXiv:1603.01810.[19] M. Kitaoka, K. Jung, Y. Yamamoto, T. Yoshida, andS. Hasegawa, Journal of Analytical Atomic Spectrometry , 1648 (2013).[20] S. Hasegawa, Proceedings of the International School ofPhysics “Enrico Fermi” , 65 (2015).[21] R. Bl¨umel, J. E. Wells, D. S. Goodman, J. M. Kwolek,and W. W. Smith, Phys. Rev. A , 063402 (2015),10.1103/PhysRevA.92.063402.[22] W. Paul, Reviews of Modern Physics , 531 (1990).[23] E. L. Raab, M. Prentiss, A. Cable, S. Chu, and D. E. Pritchard, Phys. Rev. Lett. , 2631 (1987).[24] D. R. Denison, Journal of Vacuum Science and Technol-ogy , 266 (1971).[25] K. Ravi, S. Lee, A. Sharma, G. Werth, and S. A. Rang-wala, Applied Physics B , 971 (2012).[26] D. S. Goodman, I. Sivarajah, J. E. Wells, F. A. Narducci,and W. W. Smith, Phys. Rev. A , 033408 (2012).[27] M. Abramowitz and I. A. Stegun, Handbook of mathe-matical functions with formulas, graphs, and mathemati-cal tables , National Bureau of Standards Applied Mathe-matics Series, Vol. 55 (For sale by the Superintendent ofDocuments, U.S. Government Printing Office, Washing-ton, D.C., 1964).[28] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.Flannery,
FORTRAN Numerical Recipes , 2nd ed. (Cam-bridge University Press, Cambridge [England] ; NewYork, 1992).[29] J. D. Tarnas, Y. S. Nam, and R. Bl¨umel, Phys. Rev. A , 041401 (2013), 10.1103/PhysRevA.88.041401.[30] Y. S. Nam, E. B. Jones, and R. Bl¨umel, Phys. Rev. A , 013402 (2014), 10.1103/PhysRevA.90.013402.[31] R. Bl¨umel, J. M. Chen, E. Peik, W. Quint, W. Schleich,Y. R. Shen, and H. Walther, Nature , 309 (1988).[32] P. K. Ghosh, Ion Traps , International Series of Mono-graphs on Physics No. 90 (Oxford University Press, NewYork, 1995).[33] F. G. Major, V. N. Gheorghe, and G. Werth,
ChargedParticle Traps: Physics and Techniques of Charged Parti-cle Field Confinement , Springer series on atomic, optical,and plasma physics No. 37 (Springer, New York, 2005).[34] R. Alheit, S. Kleineidam, F. Vedel, M. Vedel, andG. Werth, International Journal of Mass Spectrometryand Ion Processes , 155 (1996).[35] L. D. Landau and E. M. Lifshitz,
Mechanics , 3rd ed.,Course of Theoretical Physics (Elsevier, ButterworthHeinemann, Amsterdam, 2007).[36] B. B. Mandelbrot,
The fractal geometry of nature (W.H.Freeman, San Francisco, 1982).[37] Y.-C. Lai and T. T´el,