Vibrational modes of thin oblate clouds of charge
A numerical method is presented for finding the eigenfunctions (normal modes) and mode frequencies of azimuthally symmetric non-neutral plasmas confined in a Penning trap whose axial thickness is much smaller than their radial size. The plasma may be approximated as a charged disk in this limit; the normal modes and frequencies can be found if the surface charge density profile sigma(r) of the disk and the trap bounce frequency profile wz(r) are known. The dependence of the eigenfunctions and equilibrium plasma shapes on nonideal components of the confining Penning trap fields is discussed. The results of the calculation are compared with the experimental data of Weimer et al. [Phys. Rev. A 49, 3842 (1994)] and it is shown that the plasma in this experiment was probably hollow and had mode displacement functions that were concentrated near the center of the plasma.
Vibrational modes of thin oblate clouds of charge Thomas G. Jenkinsa) and Ross L. Spencer Department of Physics and Astronomy, Brigham Young University, Provo, Utah 84602 Received 15 February 2002; accepted 9 April 2002 A numerical method is presented for finding the eigenfunctions normal modes and mode frequencies of azimuthally symmetric non-neutral plasmas confined in a Penning trap whose axial thickness is much smaller than their radial size. The plasma may be approximated as a charged disk in this limit; the normal modes and frequencies can be found if the surface charge density profile (r) of the disk and the trap bounce frequency profile z(r) are known. The dependence of the eigenfunctions and equilibrium plasma shapes on nonideal components of the confining Penning trap fields is discussed. The results of the calculation are compared with the experimental data of Weimer et al. Phys. Rev. A 49, 3842 1994 and it is shown that the plasma in this experiment was probably hollow and had mode displacement functions that were concentrated near the center of the plasma. 2002 American Institute of Physics. DOI: 10.1063/1.1482765 I. INTRODUCTION A Penning trap see Fig. 1 uses a combination of electric and magnetic fields to confine particles of a single sign of charge, and this versatile device has been used in a variety of applications.1 12 These have included early laser cooling experiments,5 as well as related experiments in single-ion spectroscopy6 and metrology.7 These traps have also been used to study the collective behavior of single-species non-neutral plasmas, including phenomena such as rotational equilibria8 and phase transitions in strongly coupled ion plasmas.9 12 In the ideal version of this trap, described using cylindrical coordinates, a strong magnetic field B in the z-direction confines the particles radially and hyperbolic electrodes produce an electrostatic potential, r,z z 2mp 2q z2 r2 2 0 , 1 which confines particles axially. Here mp and q are, respectively, the mass and charge including the sign of the confined particles, 0 is the potential at the center of the trap hereafter assumed to be zero , and z is the bounce frequency of particles oscillating in the z-direction. In an ideal trap the bounce frequency is uniform in r. An analytic theory describing the dynamics of fluid modes in a non-neutral, zero-temperature plasma confined in a Penning trap has been developed by Dubin.13 This theory exploits the simple spheroidal shape of these plasmas14,15 to carry out a novel boundary-value problem in spheroidal coordinates. The mode frequencies of such a plasma are expressed as functions of the confining fields of the trap and the plasma aspect ratio defined for a plasma with central axial half-width zp and radius rp as zp /rp . Experimental results8,16,17 have shown good agreement with these theoretical predictions. This theory was employed by Weimer et al.,18 who used the observed mode frequencies in a plasma which had been thinned by radial transport to deduce the plasma aspect ratio. InWeimer s experiment, a pure electron plasma consisting of approximately 43,000 particles was confined in a Penning trap and held at a temperature of about 4 K. After a time sufficient for the plasma to come to global thermal equilibrium, modes were excited and observed in the plasma. For each identifiable mode, the measured mode frequency was used to calculate the plasma aspect ratio. Although the agreement was fairly good, calculated values of for different modes using Dubin s theory were found to disagree with each other by as much as 20%. It was shown through numerical simulations by Mason et al.19 that thermal effects, as well as nonideal confining fields in the Penning trap, might account for the discrepancy between Dubin s theory and this experiment. Paulson and Spencer20 extended this work, and were able to calculate, without the use of a simulation, the effects of finite plasma temperatures and nonideal confinement fields on the plasma equilibria. The basis for their calculation involved an approximation valid for small-aspect ratio plasmas. In this limit, it becomes meaningful to think of the plasma as a charged disk with a surface charge density profile (r). In this paper, the same surface-charge density approximation is employed to derive an eigenvalue equation which has, as its solutions, the eigenfrequencies and incompressible eigenmodes of an azimuthally symmetric plasma of small aspect ratio. It should be noted that the effects of surface charge induced on the conducting walls of the trap are ignored in this paper, as they are in Dubin s theory and in the equilibrium calculations of Ref. 20. The mathematical form of this eigenvalue equation is a singular integral equation of a type first studied analytically by Carleman;21,22 in this paper it is solved numerically. In Sec. II of this paper we derive a Present address: Plasma Physics Laboratory, Princeton, New Jersey 08543; electronic mail: tjenkins@pppl.gov PHYSICS OF PLASMAS VOLUME 9, NUMBER 7 JULY 2002 1070-664X/2002/9(7)/2896/13/$19.00 2896 2002 American Institute of Physics Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp this equation. In Sec. III we present the numerical method used to solve the equation, and demonstrate that this method reproduces the ideal trap results of Dubin.13 In Sec. IV we briefly review and generalize the equilibrium calculation of Paulson et al. for nonideal nonquadrupole confining fields, and use this generalization to find (r) in global thermal equilibrium for the trap used by Weimer.18 In Sec. V we discuss the qualitative behavior of the normal modes of these plasmas in ideal and nonideal confinement fields. The validity of these results at finite temperatures is discussed, and commentary relating to Weimer s experiment is also given. In Sec. VI we conclude the paper. II. DERIVATION OF THE MODE EQUATION Consider a plasma in cylindrical coordinates whose radius is so much larger than its thickness in z that it can be described by an equilibrium surface charge density profile (r). This is the natural end state of a cloud of charge in a Penning trap because transport increases the radius, and as it does so the axial electric field of the trap compresses the cloud in z. We seek a mathematical theory of the drumheadlike modes of such a plasma, restricting our attention in this paper to modes with azimuthal symmetry. There are two important effects that determine the frequencies and eigenfunctions of these modes. The first is that when the cloud is no longer confined near the trap center, the higher-order multipole moments of the trap produce a nonuniform axial bounce frequency profile z(r), so that particles at different radii vibrate at different frequencies. The second is that the vibrational motion of the cloud is also affected by mutual repulsion between different parts of the (r) profile. Ignoring modes with internal variation in z, we let the mode displacement function be purely in the z-direction, i.e., z(r) z(r)z . Such a mode is indicated in Fig. 2, where it can be seen that the displacement in z from equilibrium of different parts of the profile is denoted by the function z(r). To find the normal mode shapes and the mode frequencies of this plasma, we need to be able to calculate the axial electric field E1z r, z(r) which the plasma creates at each ring of charge composing (r). This field can be obtained from the cylindrical Green s function H(r,r ,z z ) ignoring the effects of induced charge in the trap walls which gives the electrostatic potential at (r,z), due to a ring of surface charge (r ) with width dr located at (r ,z ), as d r,z 1 0 H r,r ,z z r r dr . 2 Hence, we may find the z-component of the electric field at radius r and infinitesimal axial displacement z z(r) due to the rest of the displaced charge in the mode as E1z r, z r 1 0 0 H r,r , z r z r z r r dr . 3 In this equation, 0 is the permittivity of free space and the function H(r,r ,z z ) is given by H r,r ,z,z 1 2 m rr K m , 4 where m 4rr r r 2 z z 2 , 5 and where K(m) is the complete elliptic integral of the first kind.20 Carrying out the differentiation in Eq. 3 yields FIG. 1. A Penning trap, which is used to confine charged particles. A voltage difference applied between the endcaps left and right and the ring middle produces an electrostatic field approximating Eq. 1 . This potential confines the particles in z, while a uniform magnetic field provides radial confinement. FIG. 2. The surface shows the deformation of a thin plasma in an axisymmetric vibrational mode. The origin is shifted upward along the z-axis to better show r, r , and the axial displacement functions z(r) and z(r ). Phys. Plasmas, Vol. 9, No. 7, July 2002 Vibrational modes of thin oblate clouds of charge 2897 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp H z 1 2 m rr E m z z r r 2 z z 2 , 6 where E(m) is the complete elliptic integral of the second kind. Expanding Eq. 6 in small z z to just first-order and substituting the result in Eq. 3 yields, for the z-component of the electric field at z(r) due to all of the other rings of charge at z(r ), E1z r 1 0 P 0 1 1 r r E m r r 2 r z r z r r dr , 7 where P denotes a principal value integral. Note that the apparent quadratic singularity in the integral is actually a first-order singularity because the factor in the numerator containing z(r) vanishes when r r. With a way of calculating the perturbed electric field created by the plasma, we may now write down the z-component of Newton s second law at each radial position r in the surface charge distribution: mp z r qEext(z) r,z qE1z r,z . 8 Requiring a normal mode so that z(r) is proportional to e i t, using the known profile z(r) of the axial bounce frequency, and using Eq. 7 then leads to 2 z r z 2 r z r q 0mp P 0 1 r r E m r r 2 r z r z r r dr . 9 This is a singular integral equation of a kind first studied by Carleman, who had a simple kernel to deal with and was able to proceed analytically. The presence in our problem of the complete elliptic integral of the second kind with a complicated argument almost certainly means that numerical methods are required to solve it. Note, however, that the small-aspect-ratio analytic results of Dubin must be reproduced by this equation, so perhaps an analysis is possible. One simple result can be verified in Eq. 9 , and it is that the center of mass mode in an ideal trap, for which z(r) and z(r) are both constant, has frequency z , as expected. To obtain further results requires that we proceed numerically. III. NUMERICAL METHOD Beginning with our mode equation, Eq. 9 , we assume that for the plasma and confining fields under consideration, the surface charge density profile (r) and the bounce frequency profile z(r) are known. Beyond the plasma radius rp , the vanishing of (r) makes the integrand in Eq. 9 equal to zero, so we replace the upper limit of integration with rp . The integrand has a principal-value singularity at r r. However, the singularity can be removed by subtracting a carefully chosen term from the integrand to make it nonsingular at r r and then adding this term back in as a separate integral. Doing so, we obtain 2 z r z 2 r z r q 0mp 0 rp r r r E m r r 2 z r z r z r r r r dr q z r 0mp P 0 rp r E m r r 2 r2 dr , 10 noting that z (r) is the derivative of z(r) with respect to r. The second integrand in this expression also has a principal-value singularity at r r , so we again add and subtract appropriately to make a nonsingular integral and a simple principal-value integral, P 0 rp dr r r ln rp r r , 11 which yields an integral equation that is equivalent to Eq. 9 but which is in proper form to be approximated numerically: 2 z r z 2 r z r q 0mp 0 rp r r r E m r r 2 z r z r z r r r r dr q z r 0mp 0 rp r E m r r 2 r2 r 2 r r dr q z r r 2 0mp ln rp r r . 12 We now discretize Eq. 12 by letting the variables (r,r ) correspond to the discrete variables (rm ,rn) and by converting the integrals to discrete sums using the midpoint method. For the discrete variables rm and rn we use a cellcentered grid containing N grid points. On this grid, the grid spacing is r rp /N, and the position of the kth grid point is given by rk (k 1/2) r; k 1,2,...,N. Note that this choice of grid eliminates the r 0 and r rp singularities in the logarithmic term. We also assume hereafter that all sums run from 1 to N, except as explicitly stated. We obtain 2 z rm z 2 rm z rm q r 0mp n rn rm rn Emn rn rm 2 z rn z rm z rm rn rm rn q z rm r 0mp n rn Emnrn rn 2 rm 2 rm 2 rn rm q z rm rm 2 0mp ln rp rm rm , 13 where 2898 Phys. Plasmas, Vol. 9, No. 7, July 2002 T. G. Jenkins and R. L. Spencer Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp Emn E 4rmrn rn rm 2 , 14 and where E(m) is the complete elliptic integral of the second kind. These summations are indeterminate at rn rm , so we separate out these terms explicitly using the Kronecker delta mn . Using L Hospital s rule and Emn mn E(1) mn mn since E(1) 1 we obtain 2 z rm 0 z 2 rm z rm 1 q r 0mp n m rn rn rm rn Emn rn rm 2 z rn 2 q r z rm 0mp n m rn rn rm rn Emn rn rm 2 3 q r z rm 0mp n m rn rn Emn rn 2 rm 2 4 q r 4 0mp rm z rm 5 q z rm r 0mp n m rn Emnrn rn 2 rm 2 rm 2 rn rm 6 q z rm r 0mp rm 2 rm 4rm 7 q z rm rm 2 0mp ln rp rm rm , 15 in which we have introduced the factors i 1,(i 0,...,7) simply as a means of keeping track of terms in what follows. The left-hand side of Eq. 15 may be interpreted as a matrix operator if we write z(rm) mn z(rn). Similarly converting the right-hand side of Eq. 15 to a matrix, as discussed in the Appendix, this equation may be written as Omn z rn 2 mn z rn , 16 where Omn i 0 7 Gmn (i) 17 see the Appendix for the definitions of the matrices Gmn (i) . The procedure outlined in the Appendix does not work, however, on the top and bottom rows of Omn because the derivatives in Eq. 15 are difficult to represent at the first and last grid points. To handle these two troublesome rows we choose to apply boundary conditions. The symmetry of the modes requires z (r) 0 at r 0, and this condition can be represented on our grid by the relation 2 z r1 3 z r2 z r3 0. 18 There is no corresponding natural boundary condition at r rp , but we can at least be neutral there, and avoid having to compute too close to the square-root singularity at r rp , by requiring that z(rN) be equal to the quadratic extrapolation of the three previous points, i.e., z rN 3 z rN 1 3 z rN 2 z rN 3 . 19 We experimented with several conditions like this and this one is adequate. These conditions are implemented in the matrix equation by replacing the top and bottom rows of the identity matrix mn on the right side of Eq. 16 with zeros, calling this new matrix Rmn . Additionally, we replace the top row of the Omn matrix with the row 2s,3s, s,0,0,0, . . . ,0 and replace its bottom row with 0, . . . ,0, s,3s, 3s,s . The quantity s is the largest element of the matrix Omn , and serves as a factor to scale the matrix elements in the top and bottom rows of Omn to be approximately the same order of magnitude as the elements in the other rows. This rescaling keeps the condition number of Omn sufficiently low that numerical solutions of Eq. 16 are possible, and is of no consequence otherwise since the right-hand side of the equation corresponding to these rows is zero, that is, the first and Nth rows of Rmn are full of zeros . We obtain Omn 2s 3s s 0 0 0 O21 O22 O23 O2,N 2 O2,N 1 O2,N O31 O32 O33 O3,N 2 O3,N 1 O3,N ] ] ] ON 1,1 ON 1,2 ON 1,3 ON 1,N 2 ON 1,N 1 ON 1,N 0 0 s 3s 3s s , 20 s max max Omn . 21 Phys. Plasmas, Vol. 9, No. 7, July 2002 Vibrational modes of thin oblate clouds of charge 2899 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp This leads to the generalized eigenvalue equation, Omn z rn 2Rmn z rn , 22 which we solve using Matlab 6. A good way to test the accuracy of our numerical method is to check that it agrees with Dubin s theory,13 which gives the eigenmodes and eigenfrequencies for a cold spheroid in a Penning trap with constant z(r) z0 . Assuming that the plasma is a spheroid of constant density, as required by Dubin s theory, we obtain r 2qncoldzp 1 r rp 2 2qncold rp 1 r rp 2 , 23 where ncold p 2 0mp q2 , 24 and p is the plasma frequency. Bollinger et al.16 showed that for cold spheroids of small aspect ratio, z0 and p are related to second order in by z0 2 1 2 2 2 p 2. 25 This means that Eqs. 23 25 can be combined to give r 2 z0 2 0mp rp q 1 r rp 2 , 26 to first order in . Putting this form for (r), together with its derivative, into Eq. 15 , we obtain a numerical method for a cold spheroid of small aspect ratio. The constancy of z(r) allows us to move the 0-term in Eq. 15 to the other side of the equation. Dividing by z0 2 then turns Eq. 22 into an eigenvalue equation of the form O mn z rn Rmn z rn , 27 where 1 2 z0 2 , 28 O mn Omn z0 2 Imn , 29 and Imn is the identity matrix. Note that O mn is proportional to and thus to , so the small aspect ratio cancels out of Eq. 27 . Dubin s calculation to first order in gives the first 4 eigenvalues of this equation as 1,0 0; 3,0 5 8 ; 5,0 161 128 ; 7,0 969 512 . 30 Table I shows the comparison between the eigenvalues computed with Matlab 6 and these first-order analytic ones. Since the calculation presented here is accurate to first order in we might hope to reproduce these results exactly. The limitations of our simple midpoint integration rule and the presence of a singular (r) profile at r rp seem to limit this calculation to a relative accuracy of about 10 4. Fortunately this is sufficiently accurate for this calculation to be of physical significance. IV. COLD PLASMA EQUILIBRIA IN NONIDEAL CONFINING FIELDS The methods of the preceding section may be used to numerically find the normal modes and mode frequencies of a cold plasma in a nonideal not strictly quadrupole confining field, provided that the plasma aspect ratio is small enough that a description in terms of (r) is valid. But to proceed we need a way of calculating the surface charge density profile. Such a calculation is presented in the Appendix of Paulson and Spencer.20 It should be noted that in general it is not possible to calculate (r) because an experimenter can load almost any desired charge distribution. But under the conditions of Weimer s experiment, where the plasma was slowly expanding due to collisions, we expect the system to be close to a state of global thermal equilibrium. So it makes sense to follow the thermal equilibrium calculation of Paulson and Spencer, which will be briefly reviewed here. A. Equilibrium conditions It is traditional in the Penning trap literature to describe the electrostatic field of the trap in terms of vacuum harmonics. This has the advantage that a single set of coefficients can describe both the midplane potential e(r,0) and the particle bounce frequency z(r). But the disadvantage of this description is that the harmonic expansion does not converge for spherical radii larger than the distance to the nearest singularity in the induced surface charge density on the electrodes see, for example, Jackson23 , and in Weimer s experiment this distance is 3.5 mm, the distance from the center of the trap to the conical end cap. Since the ring electrode had a radius of 5 mm, plasmas that approach this electrode can extend into the nonconvergent region, invalidating the standard expansion. To circumvent this difficulty, we numerically compute e(r,z) using a grid and Weimer s geometry, then make separate least-squares polynomial fits to the midplane potential e(r,0) and to the axial bounce frequency z(r) in the form TABLE I. The difference between the numerically obtained eigenvalues see Eqs. 27 28 for a thin spheroid and those predicted by Dubin s theory taken to first order in . The number of grid points in the numerical calculation is N, the numerical eigenvalue is and the first-order eigenvalue from Dubin s theory is *. The subscripts on the eigenvalues are Dubin s notation for the axisymmetric modes. The computation was carried out using the eigenvalue routine in Matlab 6. N 10 1*0 30 3*0 50 5*0 70 7*0 10 10 15 2.3 10 3 9.2 10 2 4.4 10 1 20 4 10 15 2.0 10 4 1.4 10 2 9.3 10 2 40 4 10 15 2.9 10 4 1.4 10 3 1.5 10 2 80 1.5 10 14 2.4 10 4 2.4 10 4 1.5 10 3 160 1.3 10 13 1.3 10 4 3.1 10 4 2.8 10 4 2900 Phys. Plasmas, Vol. 9, No. 7, July 2002 T. G. Jenkins and R. L. Spencer Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp e r,z 0 n 1 a2nr2n 31 and z 2 r n 1 b2nr2(n 1). 32 In passing we note that if the usual description of the potential in terms of vacuum harmonics with spherical radius R and polar angle , e R, n 1 C2nR2nP2n cos 33 does converge; then we have a2n C2nBn ; b2n 4q mp C2nBnn2, 34 where Bn P2n 0 1 n 2n ! n! 24n . 35 In Appendix B of Paulson and Spencer the surface charge density functions that produce polynomial potentials are found: r Vn r2n rp 2n is produced by r 4 0Vn rp sn r/rp , 36 where the function sn(x) is given by sn x m 0 n 1 1 m 1 BnBm 1 n m 1 1 x2 m 1/2, 37 where (m 1 n ) is the binomial coefficient. These relations and the condition of global thermal equilibrium then lead to the following two equations that determine (r), given the external field coefficients a2n and the number of particles N in the cloud: r 4 0 rp V s1 r/rp n 2 a2nrp 2nsn r/rp , 38 N 8rp 0 q 4 3V n 2 m 0 n 1 a2nrp 2n 1 m n m 1 BnBm 1 2m 3 . 39 Note that Eq. 39 determines the voltageV which is needed to make (r) in Eq. 38 determinate. In a trap where only C2 and C4 are important, the equilibrium plasma surface charge distribution is given by r 3Nq 2 rp 2 1 64C4rp 5 0 9Nq 2/5 r2/rp 2 1 r2/rp 2. 40 As discussed in Paulson and Spencer, there are limits on the values C4 may have if the surface charge distribution is to all of one sign, as it must be in a Penning trap. These limits are 15 64 C4rp 5 0 Nq 45 128 , 41 where at the upper limit (0) 0; for C4 beyond this limit the center of the plasma would be oppositely charged. At the lower limit the outside edge of the plasma makes the transition to being oppositely charged. Exceeding the upper limit probably means that the plasma becomes hollow, with a space near the center that has no charge. Dropping below the lower limit is rather problematic in the case of a cold plasma, but with finite temperature Paulson and Spencer show that the state of global thermal equilibrium in this case would be for the plasma to condense on the ring electrode, so perhaps this negative limit gives the radius beyond which a plasma with C4 and q of opposite sign can no longer be confined in thermal equilibrium. We also note see Ref. 20 that for plasmas cold enough that the Debye length is small compared to the axial plasma thickness, the surface charge density (r) is the same as the plasma thickness. B. Equilibria in Weimer s experiment We now apply this equilibrium theory to the plasmas in Weimer s experiment. But Eq. 38 is the surface charge density of a zero-temperature plasma in thermal equilibrium, so the question now arises as to whether it makes sense to use it to describe experiments on thin plasmas with finite temperature. Paulson and Spencer show that the cold formula for (r) is a good approximation to warm distributions provided that D 8 zprp D0 2 40, 42 where D0 2 kT mp z 2 43 and where rp and zp are the plasma radius and half thickness of a cold spheroid with density given by Eqs. 24 and 25 . Since the total number of particles in such a spheroid is given by N 4 3 rp 2zpncold , 44 we may convert the condition in Eq. 42 into a condition on the plasma radius: rp 3Nq2 1280kBT 0 . 45 Phys. Plasmas, Vol. 9, No. 7, July 2002 Vibrational modes of thin oblate clouds of charge 2901 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp In the experiment N 43000 and T 4 K, giving the condition rp 5.3 mm. 46 The ring electrode in the experiment had a radius of 5 mm, so we expect this calculation of to work for the plasmas in the experiment. But applying this equilibrium theory to Weimer s experiment requires more than just the extra C4 coefficient. The plasma radii in this experiment were so large that a full expansion is necessary, and, indeed, the expansion in vacuum harmonics does not even converge for plasma radii beyond about 3.5 mm. But least-squares fitting a polynomial to the midplane potential e(r,0) and to z 2(r) to obtain a2n and b2n does work, and was accomplished by representing the electrodes in the experiment on a grid and doing an electrostatic solve to calculate (r,z). The non-neutral plasma equilibrium code of Spencer et al.24 was used for this calculation. It has the ability to specify potentials on arbitrary line segments in the r,z plane and hence can handle complex electrode structures like those used in Weimer s experiment. The midplane potential e(r,0) and z(r) were then passed into Matlab 6, which was used to make the least-squares fit. With the guard ring voltage set to 2.9 V, as was done in the published experiments, the coefficients a2n and b2n in Eqs. 31 and 32 were found to be the values given in Table II. Using these a2n coefficients in Eq. 38 allows us to calculate the surface charge density profile in Weimer s experiment as a function of plasma radius. The results of this calculation are shown in Fig. 3 together with two profiles of (r) computed with a global thermal equilibrium code that uses a grid.19,24 The agreement between the two is quite good except at the outer edge of the plasma, where the lack of a thermal tail in the cold equilibria causes disagreement. Recall that this good agreement was predicted by Eqs. 42 45 . V. GENERAL BEHAVIOR OF THE MODES With a method for calculating plasma equilibria, we may now use the methods described in Sec. III to calculate the normal modes of these equilibria. To do so we need to know the functions (r) and z(r), whose calculation is discussed in Sec. IV. The calculations presented here will involve equilibria described by Eq. 38 , but the method will work on other equilibria as well. We report, however, that when the outer radius of the plasma is ill-defined, as it is for the profiles of with thermal tails calculated by Paulson and Spencer,20 the method of Sec. III gives poor results. But the method seems to work well as long as (r) goes to zero at a finite radius and should work equally well for analytic profiles and fits to experimental profiles. We now discuss briefly the mode labeling convention used by Dubin and employed by Weimer in discussing the results of his experiment.13,18 We then discuss the behavior of the normal modes in simple nonideal confining fields C2 and C4 only , and finally discuss the normal modes for Weimer s experiment using the equilibrium sequence shown in Fig. 3 and the values in Table II. A. Mode labeling convention In Dubin s theory, the normal modes of a cold, uniformdensity plasma spheroid are labeled by the two integers (l,m), satisfying l m because the associated Legendre functions Pl m(x) are involved . The integer m belongs to the familiar angular variation function eim , while l describes the radial variation of the eigenfunction along the surface of the plasma. The modes discussed here all have m 0 and our assumption of rigid displacement in z requires that l be odd.18 The number of radial nodes in z(r) is given by (l 1)/2. For the special case l 1, there are no radial nodes, so the (1,0) mode for constant z(r) is simply the axial center-of-mass mode of the plasma, i.e., a rigid shift along the z -direction. B. Behavior of the modes in nonideal fields We first study the behavior of the 1,0 mode fundamental mode of a zero-temperature small aspect ratio plasma in a nonideal confining field. For simplicity, we assume throughout this subsection that only the C4 coefficient is im- FIG. 3. A sequence of surface charge density profiles for increasing plasma radius is shown by the solid curves. They are calculated from the coefficients a2n as described in the paper using the geometry of Weimer et al. Ref. 18 . The open circles give surface charge density profiles corresponding to two of the cases in the paper of Mason et al. Ref. 19 . They were computed by running the grid equilibrium code of Spencer et al. Ref. 24 which was also used by Mason . TABLE II. Fitting coefficients for the midplane potential and the axial bounce frequency in Weimer s experiment with the guard ring voltage set to 2.9 V. n a2n b2n 1 2.04 105 1.43 1017 2 4.09 102 3.81 1020 3 1.38 107 3.50 1025 4 7.45 1012 1.06 1031 5 2.61 1017 6.91 1035 2902 Phys. Plasmas, Vol. 9, No. 7, July 2002 T. G. Jenkins and R. L. Spencer Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp portant. Equation 40 gives the equilibrium surface charge density of the plasma under these conditions, while Eq. 41 gives the allowable range of values of C4 . As previously noted, the fundamental mode of a lowaspect ratio plasma in the ideal quadrupole fields of a Penning trap is a rigid shift; that is, z(r) is a constant. The entire plasma oscillates axially as a rigid body in the trap, with frequency z . If the coefficient C4 is present, then if C2C4 0, z(r) decreases with radius, while if C2C4 0, then z(r) increases with radius. See Fig. 4 and Eqs. 32 and 34 . Exploring these two cases by numerical experimentation using the method described in Sec. III, we find that the fundamental mode frequency moves similarly up or down, having a value intermediate between the extremes of z(r) across the plasma see the dashed lines in the upper panel of Fig. 4 . In doing these calculations it became apparent that the eigenfunctions were behaving like evanescent wave functions. For instance, notice in the lower panel of Fig. 4 that the displacement functions z(r) tend to decrease as they extend into the radial interval where the mode frequency is greater than z(r). More dramatic examples of this effect appear later in the paper. This behavior is backwards from the usual evanescent behavior of waves in a forbidden region where frequencies below some threshold are cut off, but the same physics is actually involved as discussed below in Sec. VD . The reason that it works backwards is that the coupling between different parts of the medium in this problem is repulsive instead of attractive. Figures 5 and 6 show that attenuation is found for higher modes as well. In both figures the upper frame displays Dubin s eigenfunctions, which are given by the formulas z r (1,0) 1, 47 z r (3,0) 1 5 2 r2 rp 2 , 48 z r (5,0) 1 7 r2 rp 2 63 8 r4 rp 4 , 49 z r (7,0) 1 27 2 r2 rp 2 297 8 r4 rp 4 429 16 r6 rp 6 . 50 FIG. 4. The properties of the fundamental 1,0 mode are displayed for two different nonideal cases with only C2 and C4 nonzero. Both cases have 43,000 electrons: 75 grid points in the mode calculation, and rp 3 mm; their surface charge density profiles (r) are readily calculable from Eq. 40 . The upper frame shows the profiles of z(r) as solid curves and the fundamental mode frequencies as dashed lines for two different choices of C4 . Case a has C2 4.2511 105 and C4 5.63 108 in SI units which is halfway from C4 0 to the upper completely hollow limit in Eq. 41 . Case b has C2 4.2511 105 and C4 3.74 108 in SI units which is halfway from C4 0 to the lower limit in Eq. 41 . The lower frame shows the mode displacement function z(r) for both a and b . Note that it decays in the region where z(r). FIG. 5. The attenuation of the first four modes for the hollow plasma case a in Fig. 4 is shown. The upper frame shows Dubin s eigenfunctions distinguished from each other by their number of zeros while the lower frame shows the numerically computed eigenfunctions also distinguished by the number of zeros . Note that they are substantially attenuated toward the outside edge of the plasma compared to the ideal eigenfunctions. FIG. 6. The attenuation of the first four modes for the plasma with C4C2 0 case b in Fig. 4 is shown. The upper frame shows Dubin s eigenfunctions while the lower frame shows the numerically computed eigenfunctions. Note the substantial attenuation toward r 0. Phys. Plasmas, Vol. 9, No. 7, July 2002 Vibrational modes of thin oblate clouds of charge 2903 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp In the plots these functions have been normalized for convenient comparison with the corresponding nonideal eigenfunctions. Figure 7 shows that this attenuation behavior persists right up to the extreme limits of Eq. 41 . In these limits, it is interesting to note that the plasma surface-charge density (r) becomes very small in certain regions of the plasma at the plasma center or at the plasma edge, depending on the sign of C2C4 , and that the modes are rapidly attenuated outside the region of space where (r) is small. This concentration of displacement toward the thinned-out center or edges of the plasma probably affects the ability of these modes to be detected by measurements of induced charge on the Penning trap electrodes as in Weimer s experiment18 , and judging by these extreme-case plots of z(r) one might think that this would be an important effect. But the electronic signals depend not on the displacement itself, but on the amount of displaced charge. A measure of this quantity is the displaced charge per unit radius q (r) z(r) (r)r, and this quantity is displayed for both of the cases of Fig. 7 in Fig. 8. Notice that the concentration effect is greatly reduced, though still present. As seen in Fig. 3, the profiles in Weimer s experiment should have been hollow, so q (r) should have been concentrated near the center of the plasma. Weimer reported that when their trap was detuned so that C4 had the opposite sign, they did not observe modes. In this case the charge displacement should have been more peaked near the outside edge of the plasma, but whether this effect has anything to do with their inability to detect modes with this detuning requires further study. Copyright Radial attenuation of the modes The physical process which causes the radial attenuation of the modes can be better understood by careful examination of Eq. 8 , repeated here: 2 z r z 2 r z r q mp E1z r,0 . 51 In this equation, the last term on the right represents the repulsive coupling between all of the charged rings composing the plasma. The actual behavior of the coupling, as we have seen from Eq. 9 , is quite complicated, and the origin of the mode attenuation is not easily seen. However, if we consider a similar system where only nearest-neighbor rings are coupled by a force of form F k z to give the correct sign for repulsive coupling and separated by an infinitesimal distance r, Eq. 51 becomes 2 z r z 2 r z r k mp z r r z r k mp z r z r r . 52 Note: a system like this, consisting of hacksaw blades dynamically coupled by ring magnets, has recently been built and studied25 and will be discussed in Sec. VD. The second term on the right is the repulsive force of the ring at r r acting on the ring at r, while the third term is the repulsive force of the ring at r r acting on the ring at r. This equation can be simplified to obtain 2 z 2 r z r k r2 mp z r r 2 z r z r r r2 . 53 In this form, the last term in Eq. 53 is a numerical approximation to the second derivative of z with respect to r, and the equation can be put into the form FIG. 7. The profiles of surface charge density dashed and axial displacement solid of the fundamental mode are shown for the two extreme cases in Eq. 41 . The upper frame is for the upper limit and the lower frame is for the lower limit. FIG. 8. For the same extreme cases shown in Fig. 7, the displaced charge per unit radius q (r) z(r) (r)r is displayed as a solid curve while the profile of is displayed as a dashed curve. Note that the concentration effects at r 0 and r rp are greatly reduced. 2904 Phys. Plasmas, Vol. 9, No. 7, July 2002 T. G. Jenkins and R. L. Spencer Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp mp 2 z 2 r k z r r2 z r . 54 This equation is analogous except for a minus sign to the time-independent Schro dinger equation in one dimension for a particle in a potential U(r) z 2(r). The mode displacement function z(r) corresponds to the particle s wave function, and 2 corresponds to the particle s energy E. In the quantum mechanical system, the wave function of a particle exhibits oscillatory behavior in regions of space where E U(r). Likewise, in the classically forbidden region where E U(r), the particle s wave function in the quantum system dies exponentially in most simple cases. In the plasma, however, because of the extra minus sign, the mode displacement function z(r) exhibits oscillatory behavior when 2 z 2(r), and attenuates exponentially when 2 z 2(r), as seen in Figs. 4 6. The oscillatory behavior of the higherorder normal modes of the plasma is exactly analogous to the excited states of the quantum particle s wave function. Although the coupling in this example is somewhat contrived, the effect is analogous to the effect which arises from the complicated coupling term in Eq. 9 perturbations in the trap fields cause the mode displacement function to attenuate into regions where 2 z(r)2. D. Comparison with the results of Weimer et al. In Dubin s theory,13 taken to first-order accuracy in aspect ratio , the frequencies of the first four odd, azimuthally symmetric normal modes of cold plasma spheroids satisfy the following: 1,0 2 z 2, 55 3,0 2 z 2 1 5 8 , 56 5,0 2 z 2 1 161 128 , 57 7,0 2 z 2 1 969 512 . 58 In Weimer s experiment,18 these results were used to estimate the plasma aspect ratio after the frequencies corresponding to the various modes had been identified. With 1,0 as an estimate for z , as in Eq. 55 , Eqs. 56 58 can each be solved for using the measured values for 3,0 , 5,0 , and 7,0 . Weimer et al. found that the experimental values of obtained through this process always satisfied the inequality 3,0 5,0 7,0 59 as in Fig. 5 in their paper , with the difference between these estimates on the order of 20%: 3,0 7,0 3,0 3,0 0.2. 60 To compare our calculation with their experimental results we used the values of a2n in Table II to generate a sequence of (r) profiles for many different plasma radii in the range 1 mm 3.94 mm. At this largest radius the plasma became completely hollow, i.e., (0) 0. Notice that according to the first two values in Table II, C2 a2 ( 1/2) and C4 a4 (3/8) should have opposite signs, which does not make a hollow plasma according to the discussion of Eq. 41 . This simple rule involving the sign of C2C4 , which works if only C2 and C4 are important, does not work here because the other terms are more important at large plasma radius. These profiles were then used to calculate the frequencies 1,0 , 3,0 , 5,0 , and 7,0 , and these values were used in Eqs. 55 58 to calculate 3,0(rp), 5,0(rp), 7,0(rp). This theoretical simulation of the experiment is shown in Fig. 9. The first thing to notice is that the computed values are in the correct order of Eq. 59 and that the computed relative range of is about 13%, about the value obtained in Mason s simulations,19 and about a factor of 2 low compared to the experiment. In the experiment the range of values of computed s when all four modes could be measured was only about 0.001 0.004 . This range is shown in the inset in the upper frame of Fig. 9 and is displayed in expanded form in the lower frame. Notice that this calculation thus predicts that the experimental plasma radius only varied between about 3.6 mm and 3.9 mm during the 25 minutes when all four modes were measured in Fig. 5 of the experimental paper.18 Looking at this range of profiles in Fig. 4 of this paper shows that this calculation also predicts that the experimental plasma was quite hollow during the measurable period of time. This hollow behavior, when taken together with the result of Paulson and Spencer20 that the plasma thickness cannot be smaller than a distance on the order of D0 see Eq. 43 , means that the quantity measured inWeimer s experiment was not simply the ratio of the plasma halfthickness at r 0 to the outer plasma radius. FIG. 9. Computed values of 3,0 , 5,0 , and 7,0 as a function of plasma radius are shown for a finely-spaced sequence of equilibria following the more coarse pattern of Fig. 3. The upper curve in each frame is 3,0(r), the middle curve is 5,0(r), and the lower curve is 7,0(r). This is the same ordering as in the experiment, but the relative spread between the lower curve and the upper curve is about a factor of 2 too small. The lower frame is an expanded view of the curves in the inset box in the upper frame. Phys. Plasmas, Vol. 9, No. 7, July 2002 Vibrational modes of thin oblate clouds of charge 2905 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp Recent work by Spencer and Robertson25 indicates that in systems comprised of weakly coupled oscillators it is hard to find modes if the individual oscillators all have the same frequency as in an ideal trap , but that detuning makes the modes more well defined. They also show evidence that in the presence of perturbations to the system like error fields due to stray patches of charge on electrodes or machining errors in Penning traps , detuning helps the modes retain their integrity in spite of the perturbations. To test whether this idea has relevance for the modes in the experiment, we added random perturbations of various strengths to the z(r) profile in our mode calculations. We find that when errors are added at a relative level of about 0.25%, the curves in Fig. 9 are changed substantially, as shown in Fig. 10. We were surprised to find that this effect not only makes the curves of computed s be more like the jagged curves of the experiment which was expected , but also that went up. In fact, by choosing the perturbation level properly it is easy to get curves that look about like those of the experiment and that have about the right value for / 3,0 of about 20%. In the simulations of Mason et al., the factor of two discrepancy in the spread could be accounted for by increasing the temperature by a factor of 4, but this error field estimate provides an alternative explanation for the increased spread. It is difficult to pursue this idea further without some knowledge of the level of error fields in the experiment, but it does suggest that error fields may have played a role in what was observed in the experiment. Note: we do not study image charge effects here, but their magnitude was estimated in Mason et al., and it was found that their inclusion probably decreases the spread. This would require perturbation levels even higher than 0.25% to account for the spread. VI. CONCLUSION When a finite-temperature plasma is confined in a Penning trap, transport and radial expansion lead to a state of global thermal equilibrium in which the plasma is very thin. A mode equation for the azimuthally symmetric, incompressible fluid modes of such a plasma has been derived and solved numerically in the limit that the plasma can be described as a thin layer of surface charge density (r). When the plasma temperature is zero and the trap fields are ideal, near-perfect agreement with Dubin s zero-temperature theory13 is obtained to first order in . Our mode equation, however, can be solved for arbitrary confining fields and temperatures provided the functions (r) and z(r) surface charge density and axial bounce frequency profiles are known. Additionally, the computation time required by this method is significantly less than that required for a particle simulation. We have also examined the dependence of the plasma shape and the normal mode eigenfunctions on nonideal components of the trap fields. We find that the amplitudes of the normal modes tend to be large in regions of the plasma where (r) is small equivalent to regions where z(r) and that the amplitude is evanescent in regions where z(r) corresponding to larger values of (r) . The equilibrium and mode calculations have also been applied to the experiments of Weimer et al.18 and we find that their plasmas were probably hollow. We also reproduce the ordering of the -values they calculated, but we cannot reproduce the amount of spread in these values unless we add random perturbations to the equilibrium fields in an ad hoc manner. ACKNOWLEDGMENT The authors wish to thank John Bollinger, National Institute of Standards and Technology, Boulder, for his help and encouragement on this project. APPENDIX: THE MATRICES The representation of Eq. 15 as a sum of matrix operators is discussed in this appendix. Each of the 8 terms corresponding to i for i 0,1,..,7 in this equation can be written as a matrix, with the top and bottom rows full of zeros because boundary conditions will be applied in these rows, corresponding to r 0 and r rp . The matrices will be denoted by the symbols Gmn (i) and each will be discussed in turn. The matrix corresponding to 0 is, of course, simply Gmn (0) mn z 2 rm A1 except for the top and bottom rows, which are full of zeros . The 1 term has a simple interpretation as a matrix multiplication Gmn (1) z(rn) with Gmn (1) q r 0mp rn rm rn Emn rn rm 2 rn ; m n, A2 Gmn (1) 0; m n, i.e., a full matrix excepting the rows reserved for boundary conditions with zeros down the main diagonal. Using the substitution z(rm) mn z(rn), the 2 term can be interpreted as a diagonal matrix Gmn (2) multiplying z(rn), where FIG. 10. These curves are just like those of Fig. 10 except that z(r) has been randomly perturbed at the 0.25% level in an attempt to assess the effect of electrostatic field errors. This effect makes the curves look more like those of the experiment and also increases the relative spread, in this case to a level about the same as the experiment. 2906 Phys. Plasmas, Vol. 9, No. 7, July 2002 T. G. Jenkins and R. L. Spencer Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp gm (2) k Gmk (1) , Gmn (2) 0 0 0 0 0 g2 (2) 0 0 0 0 g3 (2) 0 0 0 0 g4 (2) ] ] ] ] 0 0 0 0 . A3 Similarly, the 3 term is written as Gmn (3) z(rn), where the substitution, z rm z rm 1 z rm 1 2 r , A4 leads to gm (3) k rk rm Gmk (1) and Gmn (3) 1 2 r 0 0 0 0 g2 (3) 0 g2 (3) 0 0 g3 (3) 0 g3 (3) 0 0 g4 (3) 0 ] ] ] ] 0 0 0 0 . A5 In the 4 term we use z rm z rm 1 2 z rm z rm 1 r2 , A6 to obtain gm (4) q r rm 4 0mp , A7 Gmn (4) 1 r2 0 0 0 0 g2 (4) 2g2 (4) g2 (4) 0 0 g3 (4) 2g3 (4) g3 (4) 0 0 g4 (4) 2g4 (4) ] ] ] ] 0 0 0 0 . Following similar procedures with the 5 , 6 , and 7 terms yields gm (5) q r 0mp k m rk Emkrk rk 2 rm 2 rm 2 rk rm , Gmn (5) 1 2 r 0 0 0 0 g2 (5) 0 g2 (5) 0 0 g3 (5) 0 g3 (5) 0 0 g4 (5) 0 ] ] ] ] 0 0 0 0 , A8 gm (6) q r 0mp rm 2 rm 4rm , A9 Gmn (6) 1 2 r 0 0 0 0 g2 (6) 0 g2 (6) 0 0 g3 (6) 0 g3 (6) 0 0 g4 (6) 0 ] ] ] ] 0 0 0 0 , gm (7) q rm 2 0mp ln rp rm rm , A10 Gmn (7) 1 2 r 0 0 0 0 g2 (7) 0 g2 (7) 0 0 g3 (7) 0 g3 (7) 0 0 g4 (7) 0 ] ] ] ] 0 0 0 0 . 1F. M. Penning, Physica Amsterdam 3, 873 1936 . 2J. R. Pierce, Theory and Design of Electron Beams Van Nostrand, Princeton, NJ, 1949 , p. 40. 3H. G. Dehmelt, Adv. At. Mol. Phys. 3, 53 1967 . 4R.Copyright Thompson, Adv. At., Mol., Opt. Phys. 31, 63 1993 . 5D. J. Wineland, R. E. Drullinger, and F. L. Walls, Phys. Rev. Lett. 40, 1639 1978 . 6D. J. Wineland and W. M. Itano, Phys. Lett. A 82, 75 1981 . 7J. J. Bollinger, J. D. Prestage, W. M. Itano, and D. J. Wineland, Phys. Rev. Lett. 54, 1000 1985 . 8D. J. Heinzen, J. J. Bollinger, F. L. Moore, W. M. Itano, and D. J. Wineland, Phys. Rev. Lett. 66, 2080 1991 . 9T. B. Mitchell, J. J. Bollinger, D. H. E. Dubin, X.-P. Huang, W. M. Itano, and R. H. Baughman, Science 282, 1290 1998 . 10J. N. Tan, J. J. Bollinger, B. Jelenkovic, and D. J. Wineland, Phys. Rev. Lett. 75, 4198 1995 . 11J. J. Bollinger and D. J. Wineland, Phys. Rev. Lett. 53, 348 1984 . 12S. L. Gilbert, J. J. Bollinger, and D. J. Wineland, Phys. Rev. Lett. 60, 2022 1988 . 13D. H. E. Dubin, Phys. Rev. Lett. 66, 2076 1991 . 14W. D. MacMillan, The Theory of the Potential Dover, New York, 1958 , pp. 17 and 45; O. D. Kellogg, Foundations of Potential Theory Ungar, New York, 1929 , p. 195. 15D. J. Wineland, J. J. Bollinger, W. M. Itano, and J. D. Prestage, J. Opt. Soc. Am. B 2, 1721 1985 . 16J. J. Bollinger, D. J. Heinzen, F. L. Moore, W. M. Itano, D. J. Wineland, and D. H. E. Dubin, Phys. Rev. A 48, 525 1993 . 17T. B. Mitchell, J. J. Bollinger, X.-P. Huang, and W. M. Itano, Opt. Express 2, 314 1998 . Phys. Plasmas, Vol. 9, No. 7, July 2002 Vibrational modes of thin oblate clouds of charge 2907 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp 18C. S. Weimer, J. J. Bollinger, F. L. Moore, and D. J. Wineland, Phys. Rev. A 49, 3842 1994 . 19G. W. Mason, R. L. Spencer, and J. A. Bennett, Phys. Plasmas 3, 1502 1996 . 20D. L. Paulson and R. L. Spencer, Phys. Plasmas 5, 345 1998 . 21T. Carleman, Ark. Mat., Astron. Fys. 16, 19 1922 . 22F. G. Tricomi, Integral Equations Interscience, New York, 1957 , p. 186. 23J. D. Jackson, Classical Electrodynamics, 3rd edition Wiley, New York, 1999 , p. 103. 24R. L. Spencer, S. N. Rasband, and R. R. Vanfleet, Phys. Fluids B 5, 1738 1993 . 25R. L. Spencer and R. D. Robertson, Am. J. Phys. 69, 1191 2001 . 2908 Phys. Plasmas, Vol. 9, No. 7, July 2002 T. G. Jenkins and R. L. Spencer Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp