Thermal equilibrium of warm clouds of charge with small aspect ratio
Global thermal equilibrium computations are presented for non-neutral plasmas whose radial size is much larger than their axial thickness. Axial and radial density profiles are computed for both ideal and nonideal Penning trap fields. Simple results are obtained in the limits of both low and high central density. Comparison is made to the grid calculations of Mason et al. [Phys. Plasmas 3 (5), 1502 (1996)].
Thermal equilibrium of warm clouds of charge with small aspect ratio Deborah L. Paulson and Ross L. Spencera) Department of Physics and Astronomy, Brigham Young University, Provo, Utah 84602 Received 12 August 1997; accepted 22 October 1997 Global thermal equilibrium computations are presented for non-neutral plasmas whose radial size is much larger than their axial thickness. Axial and radial density profiles are computed for both ideal and nonideal Penning trap fields. Simple results are obtained in the limits of both low and high central density. Comparison is made to the grid calculations of Mason et al. Phys. Plasmas 3 5 , 1502 1996 . 1998 American Institute of Physics. S1070-664X 98 00402-9 I. INTRODUCTION When charged particles are loaded into a Penning trap, the initial density profile can be almost anything, depending on the details of the loading process. Slow transport, however, gradually causes the radius of the cloud to increase and drives the plasma toward global thermal equilibrium. If the plasma is cold enough that the Debye length is small compared to the cloud, the approach to global thermal equilibrium also causes the density profile to become uniform and the plasma shape to become a thin oblate spheroid.1 A warm plasma, however, behaves differently. Because of the kinetic energy of the particles, the plasma reaches a limit at which its thickness is on the order of the Debye length in the axial confining field. At this point, the plasma thickness can no longer decrease. As the radius grows, particle conservation requires the density to drop. This effect was observed in numerical simulations by Mason2 et al. and is the subject of this paper. Such pancake-like plasmas were studied in an experiment performed at the National Institute of Standards and Technology NIST in Boulder, Colorado by Weimer1 et al. In the experiment consisting of about 43 000 electrons that were reported to be held at a temperature of approximately 4 K the tenuous nature of the plasma made it necessary to find a nondestructive method of measuring its physical properties. The method used by Weimer et al. involved measuring the mode frequencies of the plasma and comparing them with the cold fluid theory of Dubin3 to find the plasma aspect ratio , defined as the ratio of the axial half-width zp to the radius rp of the plasma. This method worked reasonably well, but resulted in a systematic error of approximately 20% in the determination of using different modes. The Dubin theory is based on a zero-temperature, constant-density, non-neutral plasma confined axially by an electrostatic quadrupole potential. This externally imposed potential causes the cold plasma to assume a spheroidal shape. Given this uniform-density spheroidal plasma, Dubin used cold fluid theory to calculate the oscillation frequencies of the plasma. These frequencies are functions of the electrostatic confining field and the plasma aspect ratio. Since the plasmas in Weimer s experiment were too warm for thermal effects to be negligible, a possible source of the disagreement between the experimental results and Dubin s theory is the nonzero temperature. This possibility was studied in numerical simulations by Mason2 et al. They found that the experimental plasmas should have had significant variations in density due to thermal effects. Mason et al. were able to study these thermal effects for some of Weimer s plasmas by using numerical grid calculations, but both computer memory and running time limitations made it difficult to handle the very thin plasmas of the experiment. In this paper a calculation that does not require a grid, and hence uses only small amounts of computer memory, will be described. The key idea of this calculation is to exploit the small aspect ratio of these thin plasmas to separate the calculation of the density profile into axial and radial parts. Because the calculation is somewhat analytic, it also provides interesting physical insights as well as analytic results in special limits. Although Weimer s experiment served as the motivation for this study, there is one difficulty in comparing the results for a thin spheroid to Weimer s data. It was assumed in the experiment that the plasma remained spheroidal as the radius increased. Mason et al. demonstrated, however, that the nonideal components of the Penning trap field caused a distortion in the plasma shape, particularly at a small aspect ratio. For the sake of simplicity, most of the calculation presented here uses ideal confining fields. The more complicated nonideal fields encountered in experiments are discussed in the Appendices of this paper. In Sect. II we discuss global thermal equilibrium for a thin plasma confined by an external quadrupole field. Both axial and radial density profile calculations are presented. In Secs. III and IV we discuss the special cases of low and high central density, and in Sec. V, we conclude the paper. The same calculations for nonideal fields are discussed in Appendices A and B. II. GLOBAL THERMAL EQUILIBRIUM FOR THIN PLASMAS A plasma in global thermal equilibrium obeys the following form of Poisson s equation,4,5 where the potential e p represents the sum of the external potential e and the potential due to the plasma p , 2 q 0 n e q /kT Cr2, 1 a Electronic mail: spencerr@maxwell.byu.edu PHYSICS OF PLASMAS VOLUME 5, NUMBER 2 FEBRUARY 1998 1070-664X/98/5(2)/345/12/$15.00 345 1998 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 where k is Boltzmann s constant, q is the particle charge, and n is the central density. The constant C is related to the canonical angular momentum and is given by Copyright (m/2kT) ( Copyright, where is the rigid-rotor frequency of the global thermal equilibrium and where c is the cyclotron frequency including the sign of the charge.4,5 Since we are not interested here in , it is convenient to work only with C. It is assumed that the plasma is isolated in infinite space in the presence of an ideal external electrostatic quadrupole field. The boundary condition on p at infinity is that it has the constant value that makes p(0,0) 0. The external potential is the quadrupole field of a Penning trap, which can be written as e r,z m 2q z 2 z2 1 2 r2 , 2 where z is the Penning frequency, the axial bounce frequency of a single particle of mass m, and charge q in the trap. This frequency is related to the density ncold of a cold spheroid of aspect ratio by6 z F p , 3 where p 2 q2ncold 0m . 4 For low aspect ratio oblate spheroids pancakes 6 F2 1 2 2 2. 5 Cold spheroids will often be mentioned in the remainder of this paper. In fact, the point of most of this paper is to study what happens to a cold spheroid of radius rp , halflength zp , and central density ncold when its temperature is raised. This problem is quite difficult to solve for spheroids of an arbitrary aspect ratio, requiring the use of grid codes.5 Note, however, that in the limit that the Debye length is small compared to the plasma thickness, Dubin has computed thermal corrections to cold equilibria.7 However, when the aspect ratio is small, it is possible to find an approximate separable form for the axial and radial density profiles for arbitrary values of the Debye length. These separate profile calculations will now be presented. A. Axial density profile: Thin approximation When the axisymmetric plasma is thin, the axial variation of the plasma potential in the neighborhood of the plasma is greater than the variation in the radial direction. This can be illustrated by solving the simple electrostatic problem of the distribution of potential across a uniformly charged infinitely thin disk with surface charge density and radius a. The potential difference between the center and the outer edge is found to be a 0 1 2 . 6 To apply this result to a thin plasma of density n and thickness zp , we use qnzp , so that 1 r r r p r 1 2 qnzp 0a , 7 where the radial derivatives have simply been estimated by dividing by a2. The ratio of this part of 2 p to the full 2 p qn/ 0 is thus on the order of zp /a , which is small. In other words, 2 p z2 1 r r r p r . 8 It is easy to show that this result holds to an accuracy of order for cold spheroids, and it was also checked for small aspect ratio thermal plasmas by using a grid code.5 This inequality is the basis of the thin approximation used throughout this paper and is synonymous with a low aspect ratio. Using this thin approximation, the electrostatic problem can be separated into axial and radial parts. The axial equation is obtained simply by neglecting the radial part of 2 in Eq. 1 . To obtain an equation for the axial density profile, it is useful to describe the radial density profile in the plasma midplane with the dimensionless function (r): n r,0 r n0 , 9 where the scale density n0 is given by n0 z 2 0m q2 . 10 This density is nearly the same as ncold n0 F2( )ncold but makes the resulting formulas simpler. It also simplifies comparisons with experiments because it depends only on the easily measured frequency z . It is also useful to define the dimensionless potential g(z) at each radius: g z q r,z kT q r,0 kT , 11 where is the total potential. Although g(z) is actually a function of both r and z, it is used to define the axial density profile at a given radius, so in order to simplify the notation it will be referred to simply as g(z). With this definition for g(z), the density profile can be written as n r,z n0 r e g z . 12 In a further effort to make the equations dimensionless, we use a natural length scale D0 . This is a Debye length that depends on n0 and on the temperature T of the plasma: D0 0kT q2n0 kT m z 2. 13 The reason that this unusual definition for the Debye length is natural is that the axial confining electric field of the Penning trap is equivalent to the field of a uniform slab of constant charge density n0 with a sign opposite to that of the confined particles. This is similar to the case of global thermal equilibrium in infinite cylinder geometry,4,5 where the confining magnetic field is equivalent to a uniform charge 346 Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson 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 density. The length D0 is useful in the thin spheroid calculation discussed here because the Debye length defined in terms of the plasma central density varies both with plasma size and temperature, whereas D0 only depends on temperature and the Penning trap field. This length also has physical meaning here, for we shall see later that in the limit of low plasma density the plasma thickness is of order D0 . With this definition for D0 , we change to the dimensionless variable z/ D0 and obtain an approximation to Eq. 1 at each radius in the form of an axial ordinary differential equation. This differential equation is derived using 2 ext 0 and 2 p 2 p / z2. When Eq. 11 is used to evaluate 2 p / z2, the following differential equation for g is obtained: d2g d 2 1 r e g. 14 Because both the external and plasma potentials are chosen to vanish at the origin, it is clear that g(0) 0; and since both of these potentials are even functions of z, we have g (0) 0. These are the initial conditions for Eq. 14 . Note that e g is the dimensionless axial density profile at a given radius. This second-order ordinary differential equation can be solved numerically to obtain axial density profiles for various choices of (r). Figure 1 is a set of axial profiles at r 0 obtained numerically for several values of (0) ranging from 0.1 up to 0.9. Note that in the limit of (r) 1, the density profile becomes n r,z n0 r e 1/2 2, 15 so that in this limit the axial density profile is a Gaussian with a scale length of D0 . Hence, it is clear that the plasma can never be any thinner than about D0 . This low-density limit is discussed further in Sec. III. Notice also that the solutions to Eq. 14 have the property that until is very close to unity, the width of the density profile e g is characterized by D0. As approaches unity, Eq. 14 implies a uniform axial density near the center of the plasma. This prediction that the plasma density becomes uniform as approaches 1 is in error, for the correct limiting value of obtained from Eq. 3 is 1/F2( ) 1 /2. This error, however, is only of order , as expected from using the thin approximation. B. Radial density profile In order to obtain the form for the midplane radial density profile n(r,0), we must first determine the midplane plasma potential, p(r,0) see Eq. 1 . Because of the thin approximation we have already made Eq. 8 , it is difficult to use Eq. 1 in differential form to obtain p (r,0). Instead, we will use the electrostatic Green s function appropriate for p(0,0) 0, p r,0 0 2 G r,0,r ,z G 0,0,r ,z qn0 r e g z r dr dz , 16 where the cylindrical Green s function is given by G r,z,r ,z 1 42 0 m rr K m . 17 Here K(m) is the complete elliptic integral of the first kind and m 4rr r r 2 z z 2 . 18 In the limit that the Debye length D0 is small compared to the plasma radius, we again use the thinness of the plasma to treat it approximately as a disk with surface charge density (r) q n(r,z)dz. In this limit, e g is highly localized so that it may be treated as an unnormalized delta function. The axial integration may then be performed to obtain p r,0 0 r 0 H r,r r dr , 19 where H r,r 2 0 G r,0,r ,0 G 0,0,r ,0 . 20 It is possible to obtain a fairly simple formula for (r) by using the axial density profile calculation of Sec. II A as follows: r qn0 r e g dz, 21 2qn0 D0 r 0 e g d . 22 This integral can be done more easily in the variable g. To obtain the integral in this form, multiply Eq. 14 by g and integrate to obtain FIG. 1. A series of axial density profiles is shown for the following parameters: a 43 000 electron plasma at T 4 K, confined in an ideal trap with z 3.867 108 s 1. The shape of the axial profile depends on the value of 0 (0). Profiles shown are for a scale density n0 of 4.7 1013 m 3. They correspond to 0 values of 0.9, 0.75, 0.5, and 0.1. Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson and R. L. Spencer 347 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp g & g r e g 1 1/2, 23 so that the charge density becomes r &qn0 D0 r 0 e g dg g r 1 e g . 24 For simplicity, we define the function S x & 0 e g dg g x 1 e g . 25 With the dimensionless radial variable r/rp , where rp is the maximum radial extent of the cold plasma, the surface charge density may be written as qn0 D0 S . 26 Before proceeding further, it is useful to have available an approximate form for the function S(x) valid for x 0.9. An approximate form is given by the following: S x 2 1 p1x 1 p2x 2 1 p3 1 x 1 p4 1 x ln 1 x , 27 where the four parameters are best given by for accuracy to about 0.01% p1 0.823 28, p2 0.266 73, p3 1.927 65, p4 1.404 25. Defining the dimensionless midplane plasma potential h( ) by h q p ,0 kT , 28 the electrostatic equation in the plasma midplane becomes h rp D0 0 S H , d , 29 where depends on h( ) through the form for the density given in Eq. 1 : r 0 exp Cr2 1 4 D0 2 r2 h r , 30 and 0 (0). This equation for (r) can be simplified by defining D Crp 2 rp 2 4 D0 2 , 31 so that 0eD 2 h . 32 For fixed 0 and D, Eqs. 29 and 32 can be solved simultaneously to obtain the radial density profile . In the calculation described here this is done by using successive substitution with under-relaxation on Eq. 29 :5 hn 1 1 hn rp D0 0 n S n H , dp , 33 where hn is the nth iterate and where is a small underrelaxation parameter on the order of D0 /rp . For specified values of 0 and D, radial density profiles can be obtained from Eq. 33 , but it is desirable to be able to vary the temperature and have the computed plasma equilibria maintain the same total number of particles N and the same canonical angular momentum p as a reference cold spheroid. Imposing these constraints determines 0 and D, and their calculation will be discussed in the next section. Copyright Constraints It is interesting to compare this model to a cold spheroidal plasma to answer the question of what happens to the cold spheroid as its temperature is raised. If the temperature is raised without adding or removing particles or applying any torques, then both N and p must be conserved. In the limit of a small Larmor radius (rL rp), p r2 . This limit is assumed in the following derivation. For a cold spheroid, N and r2 are given by N 4 3 rp 3 ncold and r2 2 5 rp 2. 34 Note that knowledge of these two quantities and the Penning frequency of the trap z is sufficient to determine rp , zp , and n0 through Eq. 34 together with Eqs. 3 and 4 . In the remainder of this paper when zp appears, it is assumed that it was determined from N, r2 , and z in this way. For a thermal density profile, integration over changes Eqs. 34 to the following constraint equations. Insisting that particle number be conserved and using F( ) 1 yields 0 S dp 2zp 3 D0 , 35 where zp is the cold axial half-width such that zp rp . Insisting on r2 conservation yields 0 S 3 d 0 S d 2 5 . 36 These two constraint equations determine the values of 0 and D. These constants are found by simultaneously solving both the electrostatic equations Eqs. 29 and 32 and the constraint equations Eqs. 35 and 36 . There are probably many ways to do this, but one successful algorithm will be described here. 1 Begin by iterating on Eq. 33 a few times. 2 Solve Eq. 36 for D using a nonlinear zero finder, and then go back to 1 , repeating 1 and 2 together until an adequate level of convergence is achieved. 3 This procedure is followed by an adjustment of 0 based on writing Eq. 35 in the form 0,n 1 2 3 0,n zp D0 1 0 n S n d , 37 348 Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson 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 after which the algorithm returns again to 1 . This entire procedure is continued until both of the constraints and the electrostatic equation are well satisfied. In difficult cases 0 approaching unity it is necessary to under-relax in step 3 as well. To test the validity of this calculation, the results from it have been compared to the results of memory-intensive grid calculations.2,5 A global thermal equilibrium was computed on a grid with 100 radial grid points and 1200 axial grid points for a pancake with the same N and r2 as a cold spheroid with an aspect ratio 0.01. The temperature corresponded to D0zp 1. A graphical comparison between the two methods for this case is shown in Fig. 2. The central density of the grid calculation was about 1% higher than that of the calculation described here and the radial and axial profiles agreed to the same precision. The error caused by using the thin approximation described in this paper is expected to be of order , which is consistent with the error in this and other comparisons. Since the thin approximation becomes more accurate as the aspect ratio of the plasma is decreased, the radial density profile calculation described in this section should be good to better than 1% for aspect ratios less than 0.01. III. LOW CENTRAL DENSITY LIMIT It is possible to obtain the radial density profile in the special limit of low central density ( 0 1) using analytical methods. As 0 approaches zero, Eq. 29 indicates that h( ) also approaches zero so that the profile see Eq. 32 becomes a Gaussian, 0eD 2. 38 Using this form in the constraint equations, Eqs. 35 and 36 , and noting that S(0) 2, 0 and D are determined to be 0 10zp 3 2 D0 and D 5 2 . 39 This result combined with the low-density limit for the axial profile see Eq. 15 , yields a density profile that is double Gaussian, n , 0n0e 1/2 2 5/2 2. 40 This limit of small 0 applies potentially to plasmas in an experiment such as Weimer s. As the plasma sits in the trap, the root-mean-square rms plasma radius grows due to transport.1 A cold plasma under these conditions would maintain its density and flatten (zp 0) to compensate for the increasing radius. When the plasma is warm, however, its thickness, or better, its thinness, is limited by the Debye length D0 . Because the plasma cannot become thinner than about D0 as the radius grows, the central density is forced to drop in order to conserve the particle number. In time, the plasma would reach a state where the small 0 limit describes its density profile. However, this limit is not reached until h D 2. 41 If 0 is proportional to zp / D0 and h( ) is proportional to 0rp / D0 see Eqs. 39 and 29 , then Eq. 41 implies that this limit is not reached until rpzp D0 2 1. 42 Numerical experiments with the equilibrium calculation described in Sec. II show that when the factor rpzp / D0 2 0.1, the density profile is within a few percent of the approximate double Gaussian in Eq. 40 . The approach to this limit can be seen in Fig. 3. For instance, for case d in this figure we have rp 2.5 r2 0.0632 m, ncold n0 4.70 1013 m 3, zp 5.46 10 8 m, and D0 6.37 10 5 m. Hence, for this case, rpzp / D0 2 0.85, which is not small. Still, the profile shown in Fig. 3 d is Gaussian to about 10%, D 2.8 about 10% different from 2.5 , and 0 from Eq. 39 is 0.001 14, about 10% higher than the actual 0 . IV. LOW-TEMPERATURE LIMIT In the opposite limit of the calculation, where 0 1/F2( ) 1, lies the cold plasma. At T 0, global thermal equilibrium gives a constant density spheroid.6,8,9 Following Bollinger6 et al. we may therefore write for p , p r,0 m 4q z 2 p 2 r2, 43 which in the thin limit where 1 becomes using Eq. 3 p 8 mzp qrp p 2r2, 44 FIG. 2. A comparison with the memory intensive grid calculation of Mason et al. is shown. The case is a plasma consisting of 41 888 electrons at T 2.1 K in an ideal trap with z 5.5905 107 s 1. The cold spheroid density ncold is 1 1012 m 3. The cold aspect ratio is 0.01 and the cold plasma radius is 0.01 m. The top graph is a plot of the radial density profile. The open circles correspond to the data from the grid calculation, while the solid line represents the results from the calculation described in this paper. The lower graph is a plot of the axial density profile taken at the radius marked on the top graph with a solid square. Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson and R. L. Spencer 349 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp and the surface charge density of the thin cold spheroid of radius rp , (r) q n(r,z)dz, can be written as r 2qncoldzp 1 r2 rp 2. 45 Now imagine the temperature rising slightly. As long as D0 zp the plasma approximately maintains its cold spheroid shape. But, as the temperature continues to rise so that D0 approaches zp , the axial density profile becomes rounded, as shown in Fig. 1. As this evolution occurs the constraint on r2 must be maintained. Because the radial density profile is proportional to exp D 2 h( ) and both D and h( ) are very large for D0 zp rp see Eqs. 29 and 31 , to maintain r2 we must have h( ) D 2. To find the surface charge density that could give rise to such a potential, we need only refer to the cold spheroid solution in Eqs. 44 and 45 . These equations show that the charge density that would make a quadratic potential across a disk is r 1 r2 rp 2. 46 In fact, the constraint on N, which implies the same total charge for both thermal and cold spheroids, is all that is needed to see that the charge density for a thermal plasma with D0 zp must be about the same as that for a cold spheroid, i.e. warm cold . 47 Such a condition may seem to be obvious, since this limit relies on D0 being very small, i.e., the spheroid being fairly cold, but it is, in fact, an interesting result. If h( ) D 2, then the charge density is nearly that of the cold spheroid. But in the low-temperature limit both h and D are very large, so even though D 2 h( ), this does not imply that is flat or even nearly flat, since the difference between these two is the argument of an exponential. This means that even though may stay approximately the same as it is for a cold spheroid, the density profile can change quite drastically. In fact, as the axial density profile is altered with an increase in temperature, the change in the radial density profile is the means by which is able to remain the same. This condition that ( )warm ( )cold requires that the potential p(r) also be the same as in the case of the cold spheroid, giving h 8 rpzp D0 2 2. 48 Finally, because h( ) D 2, we have D 8 rpzp D0 2 , 49 valid when D is large. This approximation for D will be accurate as long as the surface charge density is close to that of the cold spheroid. We have already assumed that rpzp / D0 2 1, so that D must be very large. The question is now this: how large must D be for this approximation i.e., the equality of the surface charge densities to remain valid? One way to test the useful range of this approximation is to calculate by iteration using as in Eq. 45 . Here 0 is first found by iteration using 0 (0) see Eqs. 25 26 , 0,n 1 2zp D0S 0,n . 50 This is then used to obtain , n 1 0S 0 S n 1 2. 51 Numerical comparisons between density profiles calculated in this way and density profiles from the calculation described in Sec. II show that this approximation is good to within a few percent if D 40. This is when the profiles are very close to the cold spheroid profiles. Figure 4 shows a comparison between surface charge densities for a cold spheroid and a thermal spheroid with D 58.1. The approach to this limit can be seen in Fig. 3. For instance, for case b in this figure we have rp 2.5 r2 0.003 79 m, ncold n0 4.73 1013 m 3, zp 1.51 10 5 m, and D0 2.01 10 5 m. Hence, for this case, rpzp / D0 2 142. The approximate D from Eq. 49 is 55.6, about 4% below its true value, and the value of 0 obtained from the iteration indicated in Eq. 50 is 0.468, about 1% below the actual value of 0.473. This same limit is discussed for nonquadrupole confining fields in Appendices A and B. Some final comments on the approach to T 0 are in order here. As the temperature drops toward zero, 0 approaches 1 and S( 0) becomes logarithmically singular see FIG. 3. A series of radial density profiles is shown for the following parameters: a 43 000 electron plasma is confined in an ideal trap with z 3.867 108 s 1. In this series, the root-mean-square rms plasma radius is varied from 1.6 to 40 mm. The outer profile a corresponds to a rms plasma radius of 1.6 mm with 0 0.777 at a temperature of T 4 K. At this small radius, the profile is approaching a uniform density step. The next profile b is for a rms plasma radius of 2.4 mm with 0 0.473 also at T 4 K. The third profile Copyright is for a rms plasma radius of 10 mm and has a low density: 0 0.0364, with T 4 K. The fourth profile (d) is a special case to demonstrate the Gaussian limit. The rms plasma radius is 40 mm, the temperature is 40 K, and 0 0.001 04. 350 Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson 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 Eq. 27 . The presence of this singularity makes the iteration scheme indicated in Eqs. 50 and 51 not converge. However, if these two equations are solved by a more powerful algorithm, such as Newton s method, it is found that both 0 and approach unity, except right at the edge where 1 2 0. The correct limit is 1/F2( ) 1 /2, so once again using the thin approximation gives an answer correct to order . This scheme is important computationally because the more complex algorithm described in Secs. II B and II C converges very poorly for 0 0.8. Therefore, the only way to explore the approach to cold equilibrium is to use the approximate equations of this section. V. CONCLUSIONS When a cloud of charge in a Penning trap expands radially due to transport, it becomes thin axially. An important transition occurs when the plasma thickness approaches the length D0 kT/m z 2. When this thickness is reached, radial expansion cannot further reduce the thickness and thermal effects begin to dominate the equilibrium properties of the plasma. This can be an important effect, even in cryogenic systems. In the limit that the plasma thickness is much smaller than the plasma radius, an approximate theory that gives the equilibrium density as the product of an axial profile and a radial profile gives good results, both for ideal and nonideal electrostatic confining fields see the Appendices . In this thin limit numerical calculations are required, but they are much faster than full two-dimensional calculations. In the limit that the Debye length is sufficiently small ( D0 rpzp), it is found that even though the plasma density profile is quite different from the zero-temperature step-like profile, the line density profile n dz is nearly the same as at zero temperature. In this case a simple nonlinear algebraic equation see Eqs. 50 51 may be solved to find density profiles, both for ideal and nonideal confining fields. ACKNOWLEDGMENTS The authors acknowledge helpful discussions with John Bollinger and Grant Mason, as well as the constructive comments of the reviewer of this paper. APPENDIX A: NONIDEAL CONFINING FIELDS A more general approach to this calculation is to consider a thermal plasma in a nonideal nonquadrupole confining field. To do so requires a knowledge of the external field in the midplane of the plasma. This field e may be expressed using Legendre Polynomials: e R, n 1 C2nR2nP2n cos , A1 where (R, ) are the radius and polar angle in spherical coordinates. For a thin plasma, high-order axial terms such as z4,z6,..., may be neglected. Newton s Second Law for a charged particle of mass m and charge q then implies that in this field, C2 m 2q z 2 0 , A2 so the external potential may be expressed as e m 2q z 2 r z2 m 4q z 2 0 r2 1 e r , A3 where z(r) is no longer a constant, but a function of r in the midplane of the plasma, and e(r) is a dimensionless parameter describing the higher-order corrections to the quadrupole field. For an ideal quadrupole, z(r) z(0) and e(r) 0. For the nonideal field, z(r) and e(r) are calculated from the trap parameters: z 2 r z 2 0 1 n 2 2n2 C2n C2 P2n 0 r2n 2 z 2 0 r , A4 e r n 2 2C2n C2 P2n 0 r2n 2, A5 where P2n 0 1 n 2n ! n! 24n . A6 Note that (r) in Eq. A4 is a dimensionless radial profile defined in order to simplify the notation in the equations that follow r 1 n 2 2n2 C2n C2 P2n 0 r2n 2. A7 Once z(r) and e(r) are obtained from the Cn coefficients, the calculation proceeds analogously to the perfect quadrupole calculation. FIG. 4. A comparison between surface charge densities of a cold spheroidal plasma and a thermal plasma with D 58.1, corresponding to the radial density profile labeled b in Fig. 3, is shown. The cold corresponding to this case cuts off sharply at 1. Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson and R. L. Spencer 351 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp 1. Modified equations Because there is no longer a cold spheroid to compare with, the plasma is specified by the quantities N, r2 , and T. The main changes in the equations arise from the fact that z and D0 are no longer constants, but functions of radius, z(r) and D0(r). For convenience, the scale density n0 is given by n0 z 2 0 0m q2 , A8 and the Debye length D0(r) is given by D0 r kT m z 2 r . A9 Defining the dimensionless coordinates, and , z D0 0 and r rrms , A10 where rrms is the root-mean-square rms radial extent of the plasma, the rest of the calculation is relatively straightforward. With the new external field, the axial equation see Eq. 14 becomes d2g d 2 r r e g. A11 The surface charge density (r) can then be written as r qn0 D0 r r S r / r . A12 For the radial equation see Eq. 30 , we obtain r 0 exp Cr2 r2 4 D0 2 0 1 e r h r , A13 where e(r) is calculated from the Cn coefficients as indicated in Eq. A5 . Defining D see Eq. 31 as D Crrms 2 rrms 2 4 D0 2 0 , A14 and (r) as r r2 4 D0 2 0 e r , A15 we have see Eq. 32 0 exp D 2 h . A16 The electrostatic equation see Eq. 29 then becomes h rrms D0 2 0 0 S / D0 H , d , A17 and can be solved as before. The new constraint equations see Eqs. 35 and 36 are N 2n0rrms 2 0 D0 S / d A18 and 0 D0 S / 3 d 0 D0 S / d 1. A19 2. Low-density limit for nonideal confining fields As before, we consider the special limit of low central density. Once again, we see that the dimensionless plasma potential h( ) is small, and we can approximate as 0 exp D 2 . A20 Using this form for in the constraint equation for the number of particles yields 0 N 2 3/2n0rrms 2 0 D0 exp D 2 d . A21 Unfortunately, this equation is more complicated than its analog in the ideal quadrupole field. It is, however, possible to solve it numerically after using the angular momentum constraint and a nonlinear zero finder to solve for D. Already, however, one important feature of nonuniform external fields can be seen in Eq. A20 . If the nonideal terms in are positive, then at large radius they will overcome the D 2 term, giving a global thermal equilibrium state that is unconfined. This is discussed further in Sec. 4 of this appendix. 3. Low-temperature limit for nonideal confining fields The same argument made in Sec. IV holds for the case of the nonideal confining field, so that once again we find warm cold , A22 for sufficiently small D0( ). As before, the limit of low temperature implies that the argument of the exponent in Eq. A16 is approximately zero, D 2 h 0, A23 or h D 2 . A24 This means that the plasma potential is no longer proportional to 2, but has contributions from all even powers of . The charge density ( )cold that would produce such a potential is calculated in Appendix B. By combining Eq. A22 with the constraints on N and r2 , it is possible to determine D numerically and to perform the following iteration in analogy with Eq. 51 to determine : n 1 cold qn0 D0 S n / , A25 352 Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson 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 valid for low temperatures ( 0 0.5). Using this iteration at different values of , it is possible to obtain the profile in this limit without performing the difficult electrostatic integration in Eq. A17 . Figure 5 shows a comparison between the radial density profiles calculated using Eq. A25 and the full electrostatic calculation. The case shown is for a temperature of T 1 K, with 31 592 electrons in a nonideal confining field. The radius rrms was 1.9 10 3 m, and 0 was 0.637. The coefficients of the external field as determined by analyzing the vacuum field from the code of Ref. 5 were C2 ,C4 , 4.1220 105, 8.6641 107, 6.5565 1013, 3.2567 1017, 3.7686 1022, 5.4069 1027... . A26 Agreement is good to better than 1% for most of the profile. Beyond the radius of the equivalent cold plasma, however, the approximate iteration of course fails to give the correct nonzero value for in the thermal tail of the density distribution. 4. Possible effects of nonideal fields on radial confinement It is interesting to note that Weimer et al. discovered that a large C4 actually seemed to slow the radial expansion of the plasmas, making it possible to measure the mode frequencies for a longer time.1 Equations 1 , A3 , and A5 indicate that C4 can significantly affect the radial density profile. For instance, if C4 and C2 are of the same sign, the non-ideal part of the external field can help drive the density toward zero at large radius. For C4 and C2 of opposite sign, however, this same non-ideal part of the field tends to make the density profile become infinite at large radius. This effect is most easily seen in the low-density limit see Eq. A20 . At radii for which the nonideal corrections in the function are not small, it will dominate the behavior of . Note: strictly speaking, always dominates at large radius, but real systems are of finite size, so only radii less than the system radius should be considered. If is negative then it will cause the density profile to cut off sharply, making the global thermal equilibrium state toward which transport is driving the system better confined radially. When this effect is important, the term D 2 in Eq. A20 must tend to become positive to keep the particle count constant. This effect can give rise to hollow density profiles, as discussed in Appendix B and in Mason2 et al. This occurs typically when C4 has the same sign as C2 see Eq. A5 , although the higher-order terms often matter too. If, however, is positive C4 and C2 of opposite sign, typically , the density profile will tend to blow up at large radius, possibly giving an unconfined global thermal equilibrium state if this effect is important at the radius of the confining conductors. In practice, this would probably mean that evolution toward global thermal equilibrium would tend to make the plasma concentrated near the outer ring of the trap in this case. This influence of the external field has been seen in numerical grid experiments we have performed in which a real Penning trap geometry was used to solve for the vacuum external field for various trap tunings different guard ring voltages; see Weimer1 et al. . By solving for the potential and subtracting away the ideal components of the field, we found the higher-order contributions to the field. Then, by plotting exp q e /kT , where e represents only the higherorder contributions, it is easily determined whether the trap will provide radial confinement or radial deconfinement for the plasma in global thermal equilibrium. Using the same geometry as in the experiment by Weimer et al., we found that with the guard ring voltage at 0.0 V, the higher-order terms in the external potential were such as to cause deconfinement, but by raising the voltage on the ring, making C4 larger and of the same sign as C2 , radial confinement was enhanced. This effect was clear in the numerical experiments, but it was not possible to duplicate the numerical values of C4 given in Weimer et al. This may be due to a difference between the experimental apparatus and the numerical version of it used in the grid calculation, or perhaps due to different ways of calculating C4 . It would be nice to perform another experiment to test these ideas. APPENDIX B: COLD EQUILIBRIA IN NONIDEAL FIELDS To fully explore the problem of what thin thermal plasma equilibria look like when the confining field is not ideal, it is necessary to solve the equilibrium problem in such fields at zero temperature. As T approaches zero, Eq. 1 requires that p r,0 e r,0 CkT q r2, B1 FIG. 5. A comparison between the radial density profile from the low-temperature limit iteration Eq. A25 and from the full electrostatic calculation is shown. The iterated profile cuts off at a of approximately 1.40, while the profile from the full calculation cuts off at a of 1.44. The profile corresponds to a plasma with 31 592 electrons held at T 1 K having a rrms of 1.9 10 3 m. See Fig. 6 in Appendix B for a comparison of ( )warm to ( )cold for this case. Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson and R. L. Spencer 353 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp which means that the form of the potential produced by the plasma is known. The constant C must be determined selfconsistently by the radial size of the plasma, but the r dependence of the term in which it appears is known. Because the plasma equilibria are assumed to be thin, the cold equilibrium problem can be stated this way: given a known distribution of potential p(r) across a disk of radius a, what surface charge density distribution (r) produces it? This problem can be approached by assuming that the desired p is represented by a power series in r2 and by finding the corresponding for each power. Hence, the problem of cold thin plasma in global thermal equilibrium can be solved if the following electrostatic problem can be solved: given a charged disk of radius a, what surface charge distribution (r) will produce a potential (r) V0(r/a)2n on the surface of the disk? The solution of this problem is the main subject of this Appendix, but physical consequences are discussed at the end of the Appendix. Remember that throughout the paper it is assumed that (0) 0, which means that the constant potential at the center of the disk is always subtracted away. 1. Surface charge distributions for cold equilibria To solve for the surface charge distribution, the axisymmetric Green s function used in Sec. II B is inconvenient because it involves an elliptic integral. It is easier to obtain an analytic connection between (r) and (r) by performing a two-dimensional integral over the surface of the disk using a cylindrical coordinate system , whose origin is at x r. Using this coordinate system, the distance r from the center of the disk is given by r r2 2 r cos 2 B2 and the equation of the edge of the disk is r cos a2 r2 sin2 . B3 The potential at a distance r from the center of the disk is then given by r 1 4 0 0 2 0 r d . B4 We already know the answer in the case of (r) r2 from Eqs. 44 45 in Sec. IV: r 0 1 r2/a2 r 0a 16 0 r a 2 . B5 In the process of obtaining this result from Eq. B4 an integral of a special type is encountered, which leads to the solution of the desired problem for all powers of r2. This integral is Xn 1/2 d 2n 2 ! n 1 ! 2 4k n 1 k 2C B X Copyright r 0 n r! r 1 ! 2r 2 ! 4kX r d X , B6 where X A B Copyright 2; k 4C 4AC B2 . B7 It will be shown that this integral can be applied to the problem given by Eq. B4 by choosing surface charge distributions proportional to (a2 r2)n 1/2 note that other choices may be possible , and choosing X a2 r 2 a2 r2 2r cos 2. B8 This means that X vanishes on the edge of the disk, i.e., X ( ) 0 see Eq. B3 , and that k 1/(a2 r2 sin2 ). The integral in Eq. B6 can now be used to perform the integration in Eq. B4 . In this integral the following simplifications occur: 1 the complicated term containing the sum on r vanishes at the upper limit because X vanishes there, and 2 the lower limit gives a function that integrates to 0 when the integral in Eq. B4 is performed. The remaining term can be integrated: 0 d X sin 1 2 2r cos 4 a2 r2 sin2 0 . B9 The upper limit is simply /2 while the lower limit gives a function of that is annihilated by the integral in Eq. B4 . Hence, the following useful integration formula related to Eq. B4 is obtained: 0 2 d 0 Xn 1/2 d 2n 2 ! n 1 ! 2 4 n 1 2 0 2 a2 r2 sin2 n 1 d . B10 The remaining integral can be performed by using the binomial expansion: 0 2 a2 r2 sin2 n 1 d 2a2 n 1 m 0 n 1 Bm n 1 m r a 2m , B11 where (m n 1) is the binomial coefficient and Bm 1 m 2m ! m! 24m P2m 0 , B12 and Pm(x) is the Legendre polynomial. Using this result, we may now write down the potential (r) produced by the family of surface charge densities that are proportional to (a2 r2)n 1/2: r n 1 r2/a2 n 1/2 r 4 na 0 Bn 1 m 0 n 1 Bm n 1 m r a 2m . B13 Now we are quite close to solving the problem of finding what surface charge distribution makes the potential Vo(r/a)2n across its surface, because we have found a family 354 Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson 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 of surface charge densities that make potentials that are polynomials in (r/a)2. Note from Eq. B13 that any given value of n produces a polynomial (r) whose highest power is (r/a)2(n 1). Hence, to produce (r) V0(r/a)2n we may choose n 1 to match it in Eq. B13 , then find the proper values of m for 0 m n 1 to eliminate the lower powers of (r/a)m in the expression for (r) in Eq. B13 . This linear algebra problem may be solved to obtain the following: r V0 r/a 2n is produced by r 4 0V0 a sn r/a , B14 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. B15 Equations B14 and B15 may now be used to solve the problem of what (r) corresponds to thin global thermal equilibrium as T approaches zero. It will be assumed that the external confining field is known through the coefficients C2n in the expansion e R, n 1 C2nR2nP2n cos , B16 where (R, ) are the radius and polar angle in spherical coordinates. It is also assumed that the total number of particles N in the equilibrium is known, as well as the rms plasma radius rrms . As in the ideal case, the r2 term in e and the angular momentum term involving C in Eq. B1 combine together to give the following form for p(r) at z 0: p r V r2 a2 n 2 C2na2nP2n 0 r a 2n B17 where V corresponds to the constant D defined in Sec. II B and where a is the as-yet-unknown outer radius of the cold plasma. These two constants are to be determined by the particle number and the rms plasma radius. To make the connection between (V , a) and (N,rrms), and to find (r), we use Eqs. B14 , B15 , and B17 : r 4 0 a V s1 r/a n 2 C2na2nP2n 0 sn r/a . B18 The particle number and rrms may now be determined by taking radial moments of (r) as follows: N 2 q 0 a r r dr; P 2 q 0 a r r3 dr; rrms P/N. B19 Since Eq. B18 gives an explicit form for (r), these integrals may be performed to obtain the following formulas for N and P: N 8a 0 q V 1 3B1 2 n 2 m 0 n 1 C2na2n 1 m n n m 1 Bm 1 2m 3 , B20 P 16a3 0 q V 1 15B1 2 n 2 m 0 n 1 C2na2n 1 m n n m 1 Bm 1 2m 3 2m 5 . B21 Unfortunately, this is as far as analysis can take us. Even if only one nonideal term (C4) is included, a seventh-order polynomial equation must be solved to determine a. Therefore, to determineV and a it is necessary to solve Eqs. B19 numerically. Note, however, that if a is known, then the equation for P may be ignored andV may easily be obtained from Eq. B20 . Figure 6 shows how this calculation compares with a plasma equilibrium calculation done with the finitedifference code described in Ref. 5. The calculation was carried out in the geometry of Weimer s experiment,1 as used in the calculations of Mason2 et al. The voltage on the guard ring was 8.7 V, z 3.81 108 s 1 and the coefficients of the external field are as given in Appendix A see Eq. A26 . The plasma consists of 31 592 electrons at a temperature of 1 K, with rrms 1.9 10 3 m. The calculation was performed on a grid with 300 radial points and 1000 axial points. The two calculations should not agree exactly because of the nonzero temperature in the grid calculation, but they are quite similar. FIG. 6. The computed surface charge densities from the T 0 calculation of Appendix B b and from a low-temperature (T 1 K) grid calculation a are compared. The case is a plasma consisting of 31 592 electrons in a nonideal trap. Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson and R. L. Spencer 355 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp 2. Cold plasma shapes in presence of nonzero C4 One interesting application of this result is found in the calculation of the plasma shape at T 0. If the plasma is thin enough that the nonideal analog of the approximation F( ) 1 is valid see Eq. 5 , then at each radius a T 0 plasma would condense to the cold density, ncold r z 2 r 0m q2 , B22 where z 2(r) is given in Eq. A4 . However, if both n(r) and (r) are known, then the axial half-width zp(r) is easily determined to be zp r r 2qncold r . B23 This makes it possible to calculate plasma shapes in nonideal traps after the plasma has expanded into a pancake shape. In the following discussion it is assumed that only C4 is important to simplify things. In general, the higher-order terms are also important. For example, in a trap with nonzero C4 , the plasma shape at T 0 would be given by zp r zideal r 1 3C4r2/C2 1 32a5C4 0 2 3 r/a 2 4 15 3Nq , B24 where zideal(r) is zp(r) for the ideal quadrupole: zideal r 3Nq 1 r/a 2 8a2C2 0 . B25 Careful examination of Eq. B24 provides an interesting physical picture of what happens to the plasma as C4 is changed. If C4 is of the same sign as C2 and therefore the same as q , then increasing the magnitude of C4 causes the center of the plasma (r 0) to cave in. The critical point where (0) 0 occurs at C4 45N q 128 0a5 , B26 or, defining an average axial half-width, zp 0 a ncold r zp r r dr 0 a ncoldr dr , B27 the critical point may be expressed as C4a2 C2 45 32 zp a 45 32 . B28 For a2C4 /C2 any larger than this, Eq. B24 is no longer valid, as it implies at r 0 negative charge density, or charge of the opposite sign, which would not be confined axially. This probably means that the plasma would form into a ring for C4 beyond this limit. Because the plasmas being considered are thin, is very small ( 1), and Eq. B28 provides a stringent limit on C4 . The reader may have noticed the possibility of a zero in the denominator for C4r2/C2 1/3. This implies a loss of axial confinement at the outer edge of the profile if a2C4 /C2 gets as big as 1/3, or a2C4 C2/3. However, it is evident from Eq. B28 that the critical point at which the plasma is no longer a disk is reached long before this occurs. If C4 and C2 have opposite signs, increasing C4 has a curious effect. If C4 is large enough of the same order as in Eq. B28 , then Eq. B24 implies negative charge density beyond a certain radius. This suggests that as the plasma expands radially and approaches thermal equilibrium in the presence of some nonzero C4 , it may form into rings around a central disk, or collect on the walls of the trap, as discussed in Sec. A 4. Some of these physical effects may be seen in the radial density profile of a warm (T 0) plasma as C4 same sign as C2 is changed. As C4 approaches the critical value in Eq. B28 , both the warm and cold plasmas exhibit the same general behavior: the charge density (0) approaches zero. The cold theory predicts that the plasma will maintain uniform density, thereby being forced to cave in near r 0 until it eventually becomes a ring in order to satisfy the condition that (0) 0. The warm plasma, unable to get any thinner than about D0 , must also satisfy this condition, and it does so by dropping the central density until it, too, has formed a ring. For the T 0 case shown in Fig. 5, C4 is at about one-twentieth of its critical value, suggesting that if C4 were increased by about a factor of 20, the central density would drop to zero. We have verified this prediction in numerical experiments. 1C. S. Weimer, J. J. Bollinger, F. L. Moore, and D. J. Wineland, Phys. Rev. A 49, 3842 1994 . 2G. W. Mason, R. L. Spencer, and J. A. Bennett, Phys. Plasmas 3, 1502 1996 . 3D. H. E. Dubin, Phys. Rev. Lett. 66, 2076 1991 . 4S. A. Prasad and T. M. O Neil, Phys. Fluids 22, 278 1979 . 5R. L. Spencer, S. N. Rasband, and R. R. VanFleet, Phys. Fluids B 5, 4267 1993 . 6J. J. Bollinger, D. J. Heinzen, F. L. Moore, Wayne M. Itano, D. J. Wineland, and D. H. E. Dubin, Phys. Rev. A 48, 525 1993 . 7D. H. E. Dubin, Phys. Rev. E 53, 5268 1996 . 8D. J. Wineland, J. J. Bollinger, W. M. Itano, and J. D. Prestage, J. Opt. Soc. Am. B 2, 1721 1985 . 9J. B. Jeffries, S. E. Barlow, and G. H. Dunn, Int. J. Mass Spectrom. Ion Proc. 54, 169 1983 . 356 Phys. Plasmas, Vol. 5, No. 2, February 1998 D. L. Paulson 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