Modes and quasi-modes for m 1,2 in a gyrokinetic model for a non-neutral plasma S. Neil Rasband and Ross L. Spencer Department of Physics and Astronomy, Brigham Young University, Provo, Utah 84602 Received 27 October 1998; accepted 16 February 1999 Modes and quasi-modes for m 1,2 are studied in a gyro-kinetic model for a pure-electron plasma. Only z-independent perturbations are considered. Numerical methods are used to solve the relevant differential equations for smooth, analytic density profiles. Different temperatures and representative profiles are considered and comparison is made with the familiar cold fluid model from which the results depart but little, except at higher temperatures. A continuum component to the spectrum, present in the cold-fluid model, remains in the gyro-kinetic model to the order considered. 1999 American Institute of Physics. S1070-664X 99 04005-7 I. INTRODUCTION Some recent investigations have demonstrated renewed interest in damped quasi-modes for z-independent perturbations of a non-neutral plasma studied nearly 30 years ago by Briggs, Daugherty and Levy.1 Corngold2 published an analytic analysis for some special density profiles and the present authors published a study using numerical methods for more general profiles.3 This recent theoretical interest has been stimulated, at least in part, by experimental observations of such modes by Pillai and Gould4 and investigations of similar behavior at University of California, San Diego UCSD .5 We have been partly motivated by a desire to understand how some of the special features of the cold fluid CF model, the continuum modes and quasi-modes, are modified when additional physical effects are included. The plan of this article is to discuss the physical models in Sec. II, the numerical codes and tools used in our study in Sec. III, and in Sec. IV we consider m 1 and m 2 modes in representative hollow and monotonic equilibrium density profiles. Section V contains our conclusions. II. THE COLD FLUID AND GYROKINETIC MODELS The cold fluid CF model has been documented and discussed in the references given earlier. We use the definition of the Laplace transform given by f 0 f t ei tdt, 1 where the inversion integral is along a contour in the upper half of the complex -plane lying above all singularities. In the CF model for the Laplace transform (1)(r, ) of the perturbed potential (1)(x,t) we use the differential equation in the form 1 r d dr r d 1 dr m2 r2 1 1 r2 1 mr p 2 m 0 n 0 n 0 4 iqn 1 r,0 m 0 , 2 where m is the axial mode number from the assumed dependence of the form exp(im ). The unperturbed density is denoted as n(0)(r), its derivative dn(0)/dr as n(0) , and the corresponding plasma frequency p 2(r) 4 q2n(0)(r)/M, with q the charge and M the mass. The rotation profile 0(r) is given by v(0) r 0(r) , where v(0) is the equilibrium drift velocity, v(0) (q/M )z (0), and qB/Mc is the signed gyrofrequency. Angular frequencies with a tilde over them denote that they have been scaled by the gyrofrequency, e.g., p 2 p 2/ 2. The scaled transform frequency is / . The quantity n(1)(r,0 ) denotes the initial density perturbation. For the gyrokinetic GK model the reader is referred to the Appendix where a sketch and summary of an earlier calculation by one of us is given.6 In the GK model the differential equation for the Laplace transform (1)(r, ) takes the form b0 1 r d dr r d 1 dr m2 r2 1 b1 r2 1 b2 r d 1 dr b3r d dr r d dr n 0 r 1 I2 b4r d dr n 0 r 1 I1 4 iqn r,0 I1 , 3 where b0 r 1 2 p 2 2 p 2 2 0 p 2 m 2 4r n 0 n 0 I2 , 4 PHYSICS OF PLASMAS VOLUME 6, NUMBER 5 MAY 1999 1070-664X/99/6(5)/1435/7/$15.00 1435 1999 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 b1 r mr p 2 n 0 n 0 I1 1 p 2 2 2 0 0 p 2 2 2 m r 2 0 p 2 2 n 0 n 0 m 0 m2 2 r2 I2 , 5 b2 r rn 0 n 0 p 2 1 2 0 p 2 m 0I1 2 2 p 2 2 0 p 2 , 6 b3 r 2m 2 p 2 n 0 r2 , 7 b4 r 2m 0 2 p 2 n 0 . 8 The right-hand side of Eq. 3 represents a combination of initial value terms. Since we are not interested in following the time development of specific particular initial states, this right-hand side is arbitrary to an extent and we choose to write it in a form similar to the right-hand side of Eq. 2 . We recall that is simply an ordering parameter, a marker, to keep track of ordering in the gyro-expansion. In calculations is set equal to 1. The coefficient functions b0(r), . . . ,b4(r) of the differential Eq. 3 , as given in Eqs. 4 8 , are the same as given in Ref. 6 with a few minor corrections in second order terms. Other symbols occurring in these equations are 2 vT 2 / 2, where vT 2 2kT/M; T is temperature and k is the usual Boltzmann constant. The quantities I1 and I2 are phase space integrals coming from the integration of the perturbed distribution function to obtain the perturbed density. If we define g 1 0 2 0 p 2 9 and r,W 0g 2W2 p 2n 0 4r n 0 , 10 then I1 2 vT 2 0 W exp W2/vT 2 m r,W dW I/G, 11 I2 2 vT 4 0 W3 exp W2/vT 2 m r,W dW 1 zI /G, 12 where G 2m 2 p 2 0 n 0 4rn00 13 and z m 0g /G, 14 with n00 and p 2(0) denoting values at r 0 of the unperturbed density and scaled plasma frequency, respectively. The quantity I(z) denotes the integral I z 0 e z d ezE1 z , 15 where E1(z) is the exponential integral in the notation of Abramowitz and Stegun,7 E1 z z e t t dt ln z n 1 1 nzn nn! arg z , 16 where 0.57721 is Euler s constant. In the limit that 0 Eq. 3 becomes 1 r d dr r d 1 dr m2 r2 1 m r p 2 n 0 n 0 I1 1 4 iqn r,0 m 0 , 17 where in this limit I1 1/( m 0), showing that Eq. 17 and therefore Eq. 3 are entirely consistent with Eq. 2 and the CF model. A. Persistence of a continuous spectrum in the GK model Setting the right-hand side of Eq. 3 equal to zero gives a mode equation for (1): b0 1 r d dr r d 1 dr m2 r2 1 b1 r2 1 b2 r d 1 dr b3r d dr r d dr n 0 r 1 I2 b4r d dr n 0 r 1 I1 0. 18 Let rs denote a value of the radial coordinate where a given is such that m 0(rs),0 rs rwall . We look at solutions to Eq. 18 in the neighborhood of rs and let x r rs . With an expansion of all derivatives Eq. 18 can be written in the standard form A(r, )d2 (1)/dr2 B(r, )d (1)/dr C(r, ) (1) 0 where the coefficient functions are lengthy algebraic expressions. Singularities in this mode equation occur through the phase space integrals I1 , I2 , and their derivatives. For our purposes we note that A(r, ) is given by the expression for b0 given in Eq. 4 but with the 1/4 replaced by 5/4. Thus in the limit that x 0 we have z 0 and since I2 z ln z, A(r, ) limits to just a number. Furthermore, it has no zeros for 0 r rwall . The singularities are in the coefficient C(r, ). Since as x 0 I1 ln z, the most singular term comes from dI1 /dr dI1 /dz 1/z 1/( m 0g) and is thus similar to the singularity in the CF mode equation Eq. 2 with the right-hand side set equal to zero . 1436 Phys. Plasmas, Vol. 6, No. 5, May 1999 S. N. Rasband 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 To elucidate the continuity conditions on (1) and d (1)/dr at singularities we consider specifically an integration of the CF mode equation across the singular layer. The singular term is proportional to 1/x and because this singularity is nonintegrable, no continuity conditions on d (1)/dr are obtained. Consequently, the jump in d (1)/dr across the singular layer can be arbitrarily large. As a consequence, to the right of x 0 the constants multiplying the two linearly independent solutions to the mode equation in the CF model can be chosen to provide continuity for (1) at r rs and for (1) 0 at r rwall . This holds for any choice of in the range 0(r), 0 r rwall , and hence the continuous part of the spectrum. In the GK model this picture does not change substantially. If we take the first order limit of Eq. 18 1, 2 0 , the singularity is still proportional to 1/x. Keeping terms to second order, I(z) is proportional to ln x and the lowest order term in the mode equation from the GK expansion is integrable. However, second order terms, which now enter into the mode equation, have nonintegrable singularities of the form 1/x precluding convergence of the GK expansion at x 0. This leads again to an absence of a continuity condition on d (1)/dr at rs . As a specific example, the term in Eq. 18 with the coefficient b4 has a term with dI1 /dr which expands near the singular layer into a number of terms, one of which is proportional to 1/x. Thus for the same reason as in the CF model there is a continuous range of s for which a solution exists. The usual Frobenius analysis near the singular layer of the CF limit for Eq. 18 shows the solution in the CF model to be of the form 1 r 1 rs 1 c1x ln x c2x . 19 For the GK model, discarding the terms with coefficient functions b2 ,b3 ,b4 and keeping only the lowest order term in in b1 , we find 1 r 1 rs 1 c1x2 ln x c2x . 20 Thus we see that in the complex -plane the line of points along the real axis where m 0(r)g(r) for 0 r rwall represents a line of branch points, and this is true for either the CF or the GK model. Thus the rationale for deforming the r-contour off the real axis into the complex plane to uncover quasi-modes, as discussed in Refs. 1 3, is as valid for the differential equation 18 as it is for differential equation 2 , with right-hand side equal to 0. There is one additional complexity in solving Eq. 3 as compared to Eq. 2 . The function I(z) occurring in the coefficients of Eq. 3 has a branch cut from the origin to , normally taken along the negative real axis. In solving Eq. 3 from r 0 to r rwall it is necessary to avoid crossing the branch cut in the complex z-plane. Recall that z is a function of r through Eq. 14 , which can be a complicated path if the integration path is deformed into the complex r-plane. This becomes an issue particularly when considering hollow density profiles with deformed r-contours. In practice we monitor the path of z(r) in the complex plane as defined by Eqs. 13 and 14 and choose the branch cut for I(z) so that the path z(r) does not cross it. We restrict consideration of profiles or complex deformations to those for which such a choice is possible. III. THE NUMERICAL TOOLS We have used a number of codes to explore the behavior of solutions to Eq. 3 and the homogeneous eigenvalue equation obtained from Eq. 3 by setting the right-hand side equal to zero. To study the eigensolutions to the homogeneous equation we have used a code C1 based on finite elements and a Galerkin approximation to the differential equation. The potential perturbation (1) is expanded in set of cubic B-splines and then the homogeneous set of equations obtained from the Galerkin approximation is solved by matrix shooting.8 The purpose of this code is to find the quasi-modes by deforming the interval over which the differential equation is solved into the complex plane following the method suggested in Ref. 1. The interval 0 r rwall is analytically continued into the complex plane by choosing r s rwall s ih s , h 0 0; h 1 0; h s 0; 0 s 1. 21 The effect of the substitution represented by Eq. 21 is to push, i.e., analytically continue, the curve m (r) into the lower half of the complex -plane. With a sufficiently large deformation the complex frequencies of any existing quasimodes are left exposed above it. Some additional detail on the numerical method for finding the quasi-modes is given in Ref. 3. In order to follow the changes in for modes and quasimodes as profile parameters are changed, we have coupled the code described above with a continuation algorithm given by Allgower and Georg.9 This code C2 allows us to explore mode dependencies on profile shapes and identify modes and quasi-modes with those obtained analytically for sharp-boundary profiles. The third code C3 we use solves Eq. 3 for a given n (r,0 ) and a sequence of s approximating the Bromwich contour which surrounds the line of branch points m 0(r) for 0 r rwall . Again (1) is decomposed into a set of cubic B-splines with the appropriate boundary conditions and then a Galerkin approximation to the differential equation is made. The norm of this solution, as a function of the s, shows a peak at values of with a real part near the real part of a quasi mode frequency. The Laplace inversion is then carried out to obtain the electric field at the wall. The decay in time of this field can give a good estimate of the imaginary part of the quasi mode frequency when only a single quasi mode is present. This procedure is not sensitive to the choice made for n (r,0 ). Phys. Plasmas, Vol. 6, No. 5, May 1999 S. N. Rasband and R. L. Spencer 1437 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp IV. REPRESENTATIVE DENSITY PROFILES A. Hollow density profile For our numerical studies we consider only electrons and the first representative profile n0(r) we consider is depicted in Fig. 1 along with the corresponding rotation profile 0(r). This profile is given by the formula n 0 r n00 1 r rp 2 exp r rp , 22 where controls the hollowness of the profile, rp controls the position of the peak, and controls the steepness of the exponential decline in the density. The values 1.35, rp 1.79, and 4.5 were chosen for the profile depicted in Fig. 1. The numerical values are selected to give profiles with some semblance of experimentally measured profiles.10 The value of the density at the center is taken as n00 3.6 1012m 3. The value of the constant axial magnetic field is taken as 375.0 G. The peak value of the rotation profile for these values is 1.099 106 sec 1. 1. Resonances for m 1 Figure 2 depicting results obtained from C3 shows that we can expect two modes: one the usual diocotron mode with a frequency near 0(rwall) and the other near to the peak value of 0(r). The actual resonance frequencies for the GK model are Doppler shifted by the quantity g of Eq. 8 and are slightly different from the corresponding frequencies of the CF model, regardless of the temperature. Also in the CF model the mode with the frequency resonant at the peak of the 0 profile is the special mode with (1)(r) Ar( 0) that has been included in many studies of CF resonances,6,11 13 and has no imaginary part to the resonance frequency. In the GK model this mode is weakly unstable. Table I compares the frequencies for the CF and the GK models. For the GK model results for two values of the temperature are given; larger temperatures give, of course, a greater difference between the values of the CF and GK models. 2. Resonances for m 2 Using the same density profile as given in Fig. 1 but now with m 2, code C3 gives frequencies 0.974 max and 1.94 max as estimates for the real part of resonance frequencies. Figure 3 displays the results of this calculation in a way different from that of Fig. 2. The norm of (1) is plotted in arbitrary units along the normal direction away from the point on the Bromwich contour where Eq. 3 is solved. The modes and quasi-mode are also indicated in this figure. Code C1 gives the first as a quasi-mode with (0.97753 i(0.01098)) max and the second as an unstable mode with (1.9341 i(0.01592)) max , for a temperature of 1.0 eV. For the two temperatures of 1.0 eV and 100.0 eV we have used code C2 to compute curves for the variation of the mode frequencies as the hollowness of the density profile is changed. To obtain the curves the depth of the central depression in the density is changed by varying the parameter in Eq. 22 from 1.35 to 0.0. This corresponds to varying the ratio of the peak density to the density at the center (r 0) from 1.36 to 1.0. Figures 4 and 5 show the variation in the mode and quasi-mode frequencies for both the CF and the GK T 1.0 and 100.0 eV models. FIG. 1. An analytic radial density profile and the corresponding rotation profile. These profiles are scaled by their maximum values. The values at the center are n00 3.6 1012 m 3 and 0(0) 0.8686 106 sec 1 and at the peak max 1.099 106 sec 1. FIG. 2. For m 1 and the density profile of Fig. 1 the norm of the Laplace transform (1)(r, ) in arbitrary units is plotted as a function of normalized frequency / max as varies around the Bromwich contour depicted in Fig. 3. TABLE I. Mode frequencies (m 1) for hollow profile of Fig. 1 where max 1.099 106 sec 1. Model Diocotron/ max Peak mode/ max CF 0.2673 0.9998 GK 1.0 eV 0.2676 0.9990 i 0.0 009 492 GK 100.0 eV 0.2766 0.9895 i 0.008 591 1438 Phys. Plasmas, Vol. 6, No. 5, May 1999 S. N. Rasband 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 B. Double-hump profile As a second representative type, we consider a monotonically decreasing profile but with two shoulders as it decreases to zero. Figure 6 shows a typical example with the corresponding 0(r) profile. This profile is given by the analytic form: n 0 r n00 2 1 na r nb r , 23 where na(r) and nb(r) are functions of the form n r e 2kr2 1 e 2k rwall 2 r2 1 / 2 1 e 2krwall 2 e 2kr2 / 2 . 24 The parameter is given by the formula 2e 2krp 2 1 2e 2k rwall 2 rp 2 e 2krwall 2 1 e 2krp 2 e 2krwall 2 e 2k rwall 2 rp 2 , 25 where k determines the sharpness of the step and rp the position of the step; at r rp we require n(rp) 1/2. In the functions na(r) and nb(r) we choose respectively 0.5rwall and 0.75rwall for rp where rwall 3.81 cm with k 1 for the profile of Fig. 6. 1. Resonances for m 1 Figure 7 displays results from C3 which shows clearly the existence of the familiar diocotron wall mode ( real 0.35 106 sec 1) and suggests the existence of a heavily damped mode or quasi-mode with real 0.64 106 sec 1. Code C1 finds the wall mode at (0.35306 i(0.0)) 106 sec 1 (0.40649) max , where max 0(r 0), but finds no mode or quasi-mode associated with the bump cen- FIG. 3. The complex -plane showing the Bromwich contour dotted curve on which the differential equation 3 is solved, the mode eigenfrequencies asterisks and quasi-mode star found with code C1, the line of branch points along the real axis heavy bold , and finally in arbitrary units the log( (1) ) is plotted along the normal direction away from the Bromwich contour. FIG. 4. Frequency variation for the quasi-mode as the ratio of the peak density to the central density varies from 1.36, as depicted in Fig. 1, to 1.0 where there is no central depression in the density. Different line styles distinguish the different curves for the cold fluid CF and the gyrokinetic GK models and the GK curves are labeled with the chosen temperature. FIG. 5. Frequency variation for the unstable mode as the ratio of the peak density to the central density varies from 1.36, as depicted in Fig. 1, to 1.0 where there is no central depression in the density. Different line styles distinguish the different curves for the cold fluid CF and the gyrokinetic GK models and the GK curves are labeled with the chosen temperature. FIG. 6. An analytic radial density profile with two shoulders and correspondingly two quasi-modes. The values at the center (r 0) are given by n00 3.6 1012 m 3 and max 0(0) 0.8686 10 6 sec 1. Phys. Plasmas, Vol. 6, No. 5, May 1999 S. N. Rasband and R. L. Spencer 1439 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp tered on 0.64 106 sec 1. Results from running C3 on a corresponding CF model show a similar bump. However, for m 1 no other mode than the wall mode exists. The second step in the density profile seems to have an enhanced effect on the dynamics without giving rise to a mode or quasimode. 2. Resonances for m 2 Figure 8 from code C3 suggests two quasi-modes with the real part of the frequencies approximately 0.90 106 sec 1 and 1.37 106 sec 1. Table II gives the results from code C1 in terms of the central frequency for three different temperatures. We see little change from the CF model, 1% even for 100.0 eV. V. CONCLUSIONS In qualitative terms there is little to be gained by using the GK model instead of the CF model, certainly for temperatures in the range of those seen experimentally. With the exception of the weakly unstable m 1 mode for hollow profiles, the modes and quasi-modes have frequencies, both real and imaginary parts, that differ by at most a few percent in the two models. The instability for the m 1 modes in hollow profiles has been noted before6 and gives a growth rate only about 10% of that seen experimentally.10 The GK expansion sheds no further light on the continuum modes extant in the CF model. The procedure of inverting the Laplace transform by integrating around the Bromwich contour, as embodied in code C3, proved to be very useful in obtaining initial estimates for the frequencies of modes and quasi-modes. APPENDIX: RESULTS FOR THE GK MODEL In this appendix we summarize results from Ref. 6 necessary to obtain the differential equation 3 for the gyrokinetic model of a nonneutral plasma. We reference equations from this article by giving equation numbers preceded by an R, e.g., Eq. R12 . We assume that the unperturbed, time independent distribution function is a simple Maxwellian: F 0 r,U,W n 0 r M3vT 3 3/2 exp U2 W2 /vT 2 , A1 and that perturbed and unperturbed fields and distributions are independent of the coordinate z along the direction of the imposed magnetic field. From the Laplace transform of the Vlasov equation for the perturbed distribution function, Eq. R14 , we obtain for the Laplace transform,F (1)(x, ;U,W), of the perturbed distribution function FIG. 8. The complex -plane showing the Bromwich contour dotted curve on which the differential equation 3 is solved, the quasi-mode frequencies stars found with code C1, the line of branch points along the real axis heavy bold , and finally in arbitrary units the log( (1) ) is plotted along the normal direction away from the Bromwich contour. FIG. 7. For m 1 and the density profile of Fig. 6 the norm of the Laplace transform (1)(r, ) in arbitrary units is plotted as a function of normalized frequency / max as varies around the Bromwich contour depicted in Fig. 8. TABLE II. Quasi-mode frequencies (m 2) for the double-hump profile of Fig. 6 where max 0.8686 106 sec 1. Model 1 / max 2 / max Temperature eV 2/rwall 2 CF 1.036 i 0.002501 1.577 i 0.02445 NA NA GK 1.036 i 0.002507 1.577 i 0.02444 1.0 5.56 10 6 GK 1.037 i 0.002510 1.579 i 0.02445 10.0 5.56 10 5 GK 1.048 i 0.002433 1.592 i 0.02419 100.0 5.56 10 4 1440 Phys. Plasmas, Vol. 6, No. 5, May 1999 S. N. Rasband 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 F 1 x, ;U,W q M eim m F 0 r m r 1 m 1 r m r p 2 1 2 2 m r W2 4 D2 1 m 0 2 1 r p 2 m 1 r m r p 2 1 i m F 1 0 F 0 r q M eim 1 p 2 1 r 0 . A2 The quantity is defined in Eq. 8 and D2 1 1 r r r 1 r m2 r2 1 . A3 From Eq. R21 in the z-independent case we obtain for the Laplace transform of the perturbed density n 1 x, I1 F 1 2 z 0 I1 F 1 2 z 1 I1 F 0 2 2 I1 F 1 z 0 2 ve 0 ve 0 I1 F 0 2 z 0 z 1 ve 0 ve 1 ve 1 ve 0 Dve Dt 0 I1 F 1 Dv e Dt 1 I1 F 0 2I3 F 1 , A4 where the velocity integrals I1( ) and I3( ) are given by the formula Ik 2 M3 dU 0 WkdW. A5 We also recall that ve (q/M )z , ve , and D/Dt / t ve . Substituting F (1) from Eq. A2 into Eq. A4 and using the definitions given, we obtain for the Laplace transform of the perturbed density n 1 x, q M eim n 0 m r n 0 n 0 1 I1 m r p 2 n 0 n 0 1 I1 n 0 n 0 1 r 2D2 1 2 m r 1 2m r 2 0 p 2 n 0 n 0 m 0 m2 r2 2I2 2 0 0 p 2 I1 1 r n 0 n 0 2 0 p 2 m 0 2 I1 2 r 2 0 p 2 D2 1 2 p 2 0 m 2 4r n 0 n 0 I2 rm 0 2 n 0 d dr n 0 r 1 I1 m 2 rn 0 d dr r d dr n 0 r 1 I2 initial value terms. A6 This result for n (1) is then substituted into Poisson s equation for the Laplace transform (1). This results in Eq. 3 . 1R. J. Briggs, J. D. Daugherty, and R. H. Levy, Phys. Fluids 13, 421 1970 . 2N. R. Corngold, Phys. Plasmas 2, 620 1995 . 3R. L. Spencer and S. N. Rasband, Phys. Plasmas 4, 53 1997 . 4N. S. Pillai and R. W. Gould, Phys. Rev. Lett. 73, 2849 1994 . 5A.Copyright Cass, Experiments on Vortex Symmetrization in Magnetized Electron Plasma Columns, thesis, 1998. Available in: Reprints in Pure Electron Plasmas as 2D Fluids, from Department of Physics, University of California San Diego. 6S. N. Rasband, Phys. Plasmas 3, 98 1996 . 7M. Abramowitz and I. Stegun, Handbook of Mathematical Functions Dover, New York, 1972 , Chap. 5. 8J. P. Freidberg and D. W. Hewett, J. Plasma Phys. 26, 177 1981 . 9E. L. Allgower and K. Georg, Numerical Continuation Methods Springer- Verlag, New York, 1990 , pp. 296 311. 10C. F. Driscoll, Phys. Rev. Lett. 64, 645 1990 . 11R. A. Smith and M. N. Rosenbluth, Phys. Rev. Lett. 64, 649 1990 . 12R. A. Smith, Phys. Fluids B 4, 287 1992 . 13R. H. Levy, Phys. Fluids 11, 920 1968 . Phys. Plasmas, Vol. 6, No. 5, May 1999 S. N. Rasband and R. L. Spencer 1441 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp