Numerical calculation of axisymmetric electrostatic modes for cold finite-length non-neutral plasmas
A numerical calculation of mode frequencies for cold, non-neutral plasmas is reported. The numerical method can be applied to any axisymmetric plasma shape in a trap. Here, it is used to study axisymmetric electrostatic modes in a long conducting cylinder. These modes were previously studied by Prasad and O'Neil [Phys. Fluids 26, 665 (1983)] and by Dubin [Phys. Rev. Lett. 66, 2076 (1991)]. In contrast to Dubin's calculation, the effects of a nearby cylindrical wall, including its influence on the shape of the plasma equilibrium, are considered. It is found that for plasmas with aspect ratios (length divided by diameter) near unity the numerical results can be approximately obtained by judiciously combining Dubin's calculation, and the Trivelpiece-Gould dispersion relation for infinitely-long geometry. For aspect ratios larger than about three, the Trivelpiece-Gould dispersion relation can be used in a simple way to obtain the numerically-computed mode frequencies with an accuracy of 1%, or better. The potential use of this calculation as a plasma diagnostic is also discussed, and it is argued that at the present level of accuracy (1-2%) its usefulness is marginal, but that an improvement by an order of magnitude might make it more interesting.
Numerical calculation of axisymmetric electrostatic modes for cold finite=length non=neutral plasmas Johnny K. Jennings, Ross L. Spencer, and K.Copyright Hansen Department of Physics and Astronomy. Brigham Y(mng University. Provo, Utah 84602 (Received 28 November 1994; accepted 20 March 1995) A numerical calculation of mode frequencies for cold, non-neutral plasmas is reported. The numerical method can be applied to any axisymmetric plasma shape in a trap. Here, it is used to study axisymmetric electrostatic modes in a long conducting cylinder. These modes were previously studied by Prasad and O'Neil [Phys. Fluids 26,665 (1983)] and by Dubin [Phys. Rev. Lett. 66, 2076 (1991)]. In contrast to Dubin's calculation, the effects of a nearby cylindrical wall, including its influence on the shape of the plasma equilibrium, are considered. It is found that for plasmas with aspect ratios (iength divided by diameter) near unity the numerical results can be approximately obtained by judiciously combining Dubin's calculation, and the Trivelpiece-Gould dispersion relation for infinitely long geometrj. For aspect ratios larger than about three the Trivelpiece Gould dispersion relation can be used in a simple way to obtain the numerically-computed mode frequencies with an accuracy of I%, or better. The potential use of this calculation as a plasma diagnostic is also discussed, and it is argued that at the present level of accuracy (1-2%) its usefulness is marginal, but that an improvement by an order of magnitude might make it more interesting. 1995 American Institute of Physics. I. INTRODUCTION Plasmas consist of many charged particles whose interactions are dominated by collective effects. Non-neutral plasmas are also dominated by collective effects, but do not satisfy the usual requirement that plasmas be quasineutra!. i\n example of a non-neutral plasma is the pure electron plasma, which has been studied extensivelv bv Malmberg and coworkers. I Much of what is known ;bout non-neutral plasmas is discussed in Davidson.2 The modes of oscillation of non-neutral plasmas are interesting, both in their own right and also because of their '.' A. .-... potential as non-destructive dJagnOstlcs. " tney are only useful in this way, however, if it is possible accurately to calculate mode frequencies for confinement geometries and plasmas of interest. The simplest calculation of the mode frequencies of such plasmas is in infinitely-long geometry with constant plasma density inside the plasma radius rp , with uniform magnetic field, and with zero temperature, where the Trivelpiece-Gould dispersion relation is obtained.5 The effects of finite-length were approximately _ __ .. T . (... . 'O. treatea by prasad ana U New lOr lOng, cyundflcauy-snaped plasmas in cylindrical traps, but without including the effects of realistic equilibrium shapes. Dubin7 solved the electrostatic fluid mode equation for cold plasmas of spheroidal equilibrium shape (with the conducting walls infinitely far away). We solve numerically for axisymmetric electrostatic modes in cold plasmas with a uniform magnetic field, including the effects of a conducting cylinder and of the correct equilibrium shape. Although the discussion in this paper is limited to the effects of cylindrical conducting wails, (he numerical method presented here can be extended to arbitrary axisynlmetric geometries. The mode equation for cold non-neutral plasmas is a partial differential equation which is elliptic outside the plasma and hyperbolic inside. Such problems in which the boundary separating the two regions is an arbitrary curve are always difficult to treat.s In this paper we describe two different ways of handling this discontinuity. The first uses the cold fluid equations directly and attempts to match between the elliptic and hyperbolic regions by using a complex finitedifferencing scheme near the boundary. The second uses a warm fluid model in which the sharp boundary is replaced by a diffuse one, the resulting equations are solved, and the results are extrapolated to zero temperature. Two numerical methods have been used with both models: a form of matrix shooting9 and a singularity search method; both involve the direct solution of the set of equations that result from finite differencing the mode equation at every point on a grid. When the plasma dimensions are small enough that the cylindrical 'Nulls are effectively far a\vay, our calculations reproduce Dubin's results to within about 3%. As the plasma increases in size, two things happen. First, the plasma is no longer spheroidal in shape. It becomes more rectangular, which the calculations presented here indicate increases the mode frequencies slightly. (Note that this effect accounts for part of the 3% disagreement with Dubin's calculation; with a finite number of grid points it is difficuit to remove the effect of the conducting wall.) Second, image charges in the wall reduce the perturbed electric field, which tends to decrease the mode frequencies. Neither effect is very large unless the plasma is either close to the wall or has a length longer than the cylinder diameter. For instance, for a plasma whose radius is 50% of the wall radius and whose aspect ratio (plasma half-length ::'p divided by plasma radius rp) is less than about 3, the frequencies of the several modes studied here are only about 15% less than the mode frequencies predicted by Dubin's analysis. And since the modes in longer plasmas are described to about the same accuracy by the Trivelpiece-Gould dispersion relation, rough estimates of mode frequencies are fairly easy to make. However, using mode frequencies as diagnostics can require quite high pre- 2630 Phys. Plasmas 2 (7). July 1995 1070-664X195/2(7)/2630/1 0/$6.00 1995 American Institute of Physics II. MODE EQUATION AND NUMERICAL METHOD in conjunction with other diagnostic methods. the paper and the warm-fluid calculation, which is used to get zero-temperature results by extrapolation, is described in the Appendix. Prasad and O'Neil showed that the thermal-equilibrium density of a cold axisymmetric plasma is constant inside the plasma and zero outside the plasma, the two regions being separated by a boundary of zero thickness (sharp h1lnJl_ln...r...l.."....:...1...T'-"U1..JTJ\ . 10 r.yL'.h...IP. ..."..JU "u:.1..1.C..,O...A..Copyright.,holn.Lnv:,.1...P.....r...l.. Ut)o. .f.'.: :J:U.LYL\1-.h....JT '"OJ'r" ,..........\l.11-1.'."..U. o".Copyright.', .r....p....npJr."o"".'c\"p.In1...tL_ given choice of the mode frequency w the set of equations nuities is that as_ (;) is v3..t.t1.ed, the strllcture of changes, cussed here, so the z-derivative term in Eq. (1) is negative inside the plasma, maleJng the equation hyperbolic in (r,z) inside the plasma. Thus, the equation is of mixed type, with the type change occurring along a curved boundary, making it very. difficult to solve,S even-with numerical methods. - The numerical solution of this mode problem proceeds as follows. We first write the finite-difference approximation to Eq. (1) at every point on a grid in rz-space, producing a set of algebraic equations that are linear in the values of 1> at the grid points. The resulting set of coupled equations has three difficulties. The first is that the linear system to be solved is singular when w is near a mode frequency, and near <.'l1 "h nnlntco th ':H"l'l1r-at nlul'lpri .... -al conlntinn nf l1n -ar couC'tPITlC' ........."'1.-1.- yvl.-I-I.-...... ....1-1- ..... "'............I.\.4...'" "-1.- ...... 1-"-1-"'1.-.1."'\.4.1. ... V.l. ......... .1.'-".I..1. V.l. J...J..l.l. ..... \.4.I. uJu..."'.I.I-.I.U is difficultY The second problem is that Eq. (1) is of mixed type, so standard iterative methods for solving this linear system do not converge well. These two problems are dealt with by keeping the grids as coarse as possible, by using a non-iterative banded linear system solver with high precision (which requires a lot of computer memory), and by putting up with somewhat noisy eigenfunctions. The third problem is that this is an eigenvalue problem, so the linear system is homogeneous and w is unknown. This problem is handled in two different ways, which will now be discussed. A. Two-dimensional matrix shooting The first technique fixes the value of the wave potential at an arbitrarily chosen grid point, called the pivot point. This is done by replacing the finite-difference equation for that grid point with an equation that forces cP at the pivot point to be one. This makes the linear system inhomogeneous, and is equivalent to choosing the amplitude of the mode. For a causing =0 curves in the rz-plane to move across the pivot point, inverting the value of the finite-difference approximation to the mode equation there. Hence, care must be taken to distinguish between sign reversals that are merely discontinuities of R(w) and those that indicate actual solutions of the set of finite difference equations. Some of this can then be solved. But for general choices of w, the solution has a kink at the pivot, similar to what happens to tent fabric when a pole is pushed up against it from inside. The residue R( w) is defined to be the finite-difference approximation to Eq. (1) evaluated at the pivot point divided by the maximum of IcPl on the grid, and is a measure of the size of the kink. To find mode frequencies, w is varied until the kink disappears, as indicated by the vanishing of R(w). (Note that this method is essentially the method of matrix shooting in two dimensions.9) At such values of w the obtained by the inversion should be t..he wave potential of an actual mode, but in practice is often rather inaccurate because it is difficult to maintain accuracy when solving a nearly-singular system. This means that frequencies are known quite a bit more accurately than are the mode eigenfunctions. An example of R(w) is shown in Fig. 1. Note that the zeroes of this function cannot simply be identified by searching for sign reversals, because ihere are jump discontinuities that cross the horizontal axis. The reason for these disconti- 1 a (' acP) ar( w;(r,z ) acPl- (1) -- r- +- 1- --0 r ar, ar, az L , w2 , az , cision, as in the work of Weimer et al., where the aspect ratio was found hy tal<-ing differences among closely-spaced frequencies3 and in the work of Tinkle et al.4 For such work accuracies of 10-15% are not good enough and better numerical methods like those described here are probably required. The cold sharp-boundary model discussed here seems to give accuracies in the 1-4% range while the warm-fluid model does somewhat better at 1-2%. This level of accuracy is probably not sufficient to be truly useful as a diagnostic, but might be sufficient to give useful infonnation when used ing the shape of the boundary, which is typically not spheroidal, but conforms more to the cylindrical shape of the outer conductor. The numerical method reported here makes it possible to calculate mode frequencies for such plasmas, including the effects of nearby walls on the dynamics. We use it to find axisymmetric electrostatic drift-fluid modes in cylindrical geometry, but we are confident that it couid be modified to handle other plasma models and more general The plan of the paper is as follows. Section II discusses the mode equation and the numerical method that is used to solve it and Section III compares this method to previous calculations. Section IV discusses the effect of a conducting cylindrical wall on the mode frequencies of non-neutral plasmas, and features a table comparing the results of the rather complicated calculations of this paper with simpler calculations. It is found that for aspect ratios of 3 or higher, the Trivelpiece-Gould dispersion relation can be used in a rather simple way to find all of the numerically determined fre- ""1'll""Pnf_'l C tn..., !an -IU:.::lI_ ,....f' llr f"'U n.f 10/" n1" hpttp.r p,....t1nn\J I"'nnf'lnF"lpc J ...,.0. ,L, 'OJ ....,_ _ ..., 'f _..., '-"0 where cP is the wave potential, where wp is the plasma fre- ....non",u /0"2,.,,. '-n2M /.". ... \/ IIIf 1:Y1ha.ra. 1/1_1.". .... \ ;C' f-h.CII 'iu\.u.....J '-Up \' ,-4) - "1 ' .U\.' ,-4/' C(p:,.,L, YVJ.J.'-.d...... 1 .0\' , J J. UJ.\.; equilibrium particle density, where q and M are the charge and mass, respectively, of the constituent particles, and where rand z are radial and axial cylindrical coordinates. Since no is constant inside the plasma and zero in the vacuum region, so is w;. Thus, Eq. (1) is elliptic outside the plasma. We will see later that w2< w; for the modes disgeomeLries as \vell. Prasad and 0' Neil6 give the drift-fluid mode equation for cold non-neutral plasmas in a uniform magnetic field parallel to the z-axis. This model assumes that all frequencies are much less than the cyclotron frequency. For axisymmetric modes (m = 0) their equation reduces to Phys. Plasmas, Vol. 2, No.7, July 1995 Jennings, Spencer, and Hansen 2631 (2) o z/r..... 1.950 in Fig. 2. The frequency of this solution is 0.231wp ' FIG. 2. The wave potential of the center-of-mass mode for a numerically computed equiiibrium. We found this mode using a grid resolution of 38 radial and 160 axial cells in a computation region that is 4.0 cm in radius (axis to wall) and 20.0 em in length (midplane to end-wall). The pivot point for the computation is close to the axis near the end of the plasma, but just inside. The plasma is 1./ Col in radius and 5.3 em in half-length. The frequency of this mode is 0.230w". This mode potential does not look like either a Trivelpiece-Gould mode or a Dubin ode, yet the Dubin calculation obtains nearly the correct frequency (0.242w,,). two hundred, The first has the appearance of being a dear single mode, but the mode shown in Fig. 3 looks like it might be the superposition of more than one mode. In such cases, the residues associated with the different modes could be cancelling each other rather than being zeroed independently, so that neither is actually a solution of Eq. (I), or the problem might simply be a result of solving a large system of linear equations that is nearly singular. In any case the wave potential shown in Fig. 3 is probably not a good approximation to a singie mode of Eq. (1). librium, pivot point a..'1d grid resolution as the center-of mi1Ss mode sho\vn FIG. 3. A solution of the nnite-difference equations with a frequency similar to that of the center-of-mass mode. Judging by its appearance. it is highly doubtful that it is a mode. This solution was obtained using the same equi- 11 I I \" I j 01-- " - i ---I -,IL- -L, --,,'I 0.18 0.23 0.28 FIG. 1. Residue versus w. The residue is expected to vanish when w is the frequency of a mode. Besides the zeroes, there are also sign reversals caused by jump discontinuities. In this example, the ratio of plasma radius rp to waH radius r. is 0.27 and the aspect ratio pirp is 4.82. The grid has 38 radial cells and 160 axial cells in a cylindrical computation region whose radius r. is 4.0 COl and whose length, from midplane to end-wall, is 20 Col. can be done by careful programming, but more often it is necessary to search for zeroes on a plot of R( w). The major drawback of this method for this problem is that the spectrum of Eq. (1) is very awkward. The mode frequencies of Eq, (1) are countably infinite in number and are all confined to the interval between a and wp ' This can be seen by multiplying Eq. (l) by cP and integrating over the interior of the conducting cylinder. After integrating by parts. the equation can be solved for w giving of lov/-order modes may also encounter other solutions of Eq. (I) with nearby frequencies, corresponding to higherorder modes. As an example of this difficulty Fig. 2 shows a plot of the wave potential of the center-of-mass mode obtained by our numerical method for a plasma shape obtained from a numerical equilibrium calculation. [The center-of-mass mode is basically just the axial sloshing of the plasma in the electrostatic confining field. In Dubin's notation it is the (l,0) r'l1(\t"lp 1 Pinl1r 1: hr\n,c QnAthAT Atlltinn nf th finitoO rliff,or ......... 'J ....,.J b H 'J V ""' ...., ... -"'J "' ...., "'- .1 ... 2_ 2 f v(r'Jcf>Ir'Jz)2 r dr dz w -wi' J[(J<pIIJr)2+(o<pIJz)2]r dr dz where wp is the value of the plasma frequency inside the cold plasma and where the subscript V on the integral in the numerator indicates that it is to be taken only over the volume of the plasma. The imegral in the denominator is taken over the entire computing region. Modes with short radial \vavelength have large average values of (ac/>I 8r) 2, 'vvhiIe modes with short axial wavelength have large average values of ( r'J