Simulations of the instability of the m = 1 self-shielding diocotron mode in finite-length non-neutral plasmas
The "self-shielding" m = 1 diocotron mode in Malmberg-Penning traps has been known for over a decade to be unstable for finite length non-neutral plasmas with hollow density profiles. Early theoretical efforts were unsuccessful in accounting for the exponential growth and/or the magnitude of the growth rate. Recent theoretical work has sought to resolve the discrepancy either as a consequence of the shape of the plasma ends or as a kinetic effect resulting from a modified distribution function as a consequence of the protocol used to form the hollow profiles in experiments. Both of these finite length mechanisms have been investigated in selected test cases using a three-dimensional particle-in-cell code that allows realistic treatment of shape and kinetic effects. A persistent discrepancy of a factor of 2-3 remains between simulation and experimental values of the growth rate. Simulations reported here are more in agreement with theoretical predictions and fail to explain the discrepancy.
Simulations of the instability of the m 1 self-shielding diocotron mode in finite-length non-neutral plasmas Grant W. Mason and Ross L. Spencer Department of Physics and Astronomy, Brigham Young University, Provo, Utah 84602 Received 26 March 2002; accepted 29 April 2002 The self-shielding m 1 diocotron mode in Malmberg Penning traps has been known for over a decade to be unstable for finite length non-neutral plasmas with hollow density profiles. Early theoretical efforts were unsuccessful in accounting for the exponential growth and/or the magnitude of the growth rate. Recent theoretical work has sought to resolve the discrepancy either as a consequence of the shape of the plasma ends or as a kinetic effect resulting from a modified distribution function as a consequence of the protocol used to form the hollow profiles in experiments. Both of these finite length mechanisms have been investigated in selected test cases using a three-dimensional particle-in-cell code that allows realistic treatment of shape and kinetic effects. A persistent discrepancy of a factor of 2 3 remains between simulation and experimental values of the growth rate. Simulations reported here are more in agreement with theoretical predictions and fail to explain the discrepancy. 2002 American Institute of Physics. DOI: 10.1063/1.1488600 I. INTRODUCTION Non-neutral plasmas, typically ions or electrons, can be confined for long periods of time in a cylindrical Malmberg Penning trap similar to that shown in Fig. 1. A stiff axial magnetic field confines the particles radially, and charged rings at the ends of the otherwise grounded cylinder provide electrostatic longitudinal confinement. Diocotron modes are azimuthal drift waves in the cylindrical plasma that vary spatially as exp(im ). The theory of diocotron modes in nonneutral plasmas has its origins in seminal papers by Levy,1 Briggs, Daugherty and Levy,2 and the comprehensive treatment of non-neutral plasmas by Davidson.3 Of these modes, the m 1 mode occurs in two manifestations for nonmonotonic density profiles. The simplest wall mode occurs when the plasma is offset radially from the center and revolves around the center of the cylinder as a result of E B drift from the electric field of the induced charge on the walls and the longitudinal magnetic field. For hollow density profiles the azimuthal flow of the plasma exhibits shear and a rotation frequency profile 0(r) that rises with increasing radius from the center, peaks, than decreases to the wall. The mode frequency of the first kind of m 1 diocotron mode is the value of this frequency profile at the wall. A second kind of m 1 mode self-shielded has a frequency near the peak of the frequency profile. In the infinite length approximation, both modes are predicted not to be exponentially unstable for all radial density profiles of the plasma column. In contrast, when the plasma column is of finite length, the first kind of m 1 mode remains stable for all radial density profiles, but the second has been experimentally shown to be exponentially unstable for hollow density profiles.4 6 Several theoretical attempts have been made to understand the origin of the instability. Smith and Rosenbluth,7 Smith,8 Rasband et al.,9 Rasband,10 and Finn11 have considered various mechanisms that might account for the instability but have failed to account for the exponential character of the stability and/or the size of the growth rate. In particular, Smith8 has drawn attention to finite length effects and Finn et al.11 have drawn attention to the importance of the shape of the ends of the plasma based on an analogy to vortex stretching from topography variations in shallow fluid dynamics for geophysical flows. The theory of Finn et al. when adapted to vortex dynamics in non-neutral plasmas, demonstrates that the radial variation of the equilibrium plasma length causes compression of the plasma parallel to the magnetic field while conserving the line integrated density. Their theory predicts the observed exponential growth of the instability, but predicts a growth rate that is still somewhat more than a factor of 2 less than a test case taken by Finn et al.11 from data of Driscoll.5 Coppa et al.12 have refined the theory of Finn et al. by using a more rigorous definition of the length of the plasma column, by introducing an effective electrostatic potential to calculate the E B drift on a string of variable axial density, and by calculating the perturbation of the plasma length induced by density variations using a Green s function. They are able to separate various effects from one another for comparison, but their complete model predicts a slightly lower value than that of Finn et al. for the ratio of growth rate to real frequency for the test case, thus failing to remove the discrepancy between theory and experiment. Hilsabeck and O Neil13 treat finite length diocotron modes by relating perturbed charge density to perturbed potential using a Green s function analysis. Like Finn et al., they find that finite column length makes exponential growth possible. They also observe that experimental procedures to produce plasmas with hollow profiles involve lowering the confining ring potentials and dumping preferentially the particles in the tail of the original Maxwellian velocity distribu- PHYSICS OF PLASMAS VOLUME 9, NUMBER 8 AUGUST 2002 1070-664X/2002/9(8)/3217/8/$19.00 3217 2002 American Institute of Physics Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp tion, thus effectively truncating the velocity distribution near the center of the plasma. For certain distributions of axial energies, the instability can be substantially affected. They conclude that the instability cannot be understood in terms of plasma shape alone, but kinetic effects associated with non- Maxwellian velocity distributions created in the formation of the hollow profile plasmas might be sufficient to close the gap with experiment. In this paper we report the results of three-dimensional particle-in-cell simulations that also attempt to account for the remaining quantitative disagreement. We simulate the test case considered by Finn et al.We also consider the effect of modifying the velocity distribution as suggested by Hilsabeck and O Neil by simulating the experimental method used to produce hollow profiles. The simulations predict growth rates in rough agreement with the theory of Finn et al. Kinetic corrections increase the growth rates, but fail to lift the growth rate sufficiently to agree with the experimental test case. In Sec. II we briefly review the theory of the diocotron modes. In Sec. III we describe the model and numerical methods used in our study. In Sec. IV we apply our method to test cases of interest to the theory of Finn et al. and to the theory of Hilsabeck and O Neil. Finally, in Sec. V, we draw conclusions from our work. II. THEORY We first consider electrons confined in cylindrical geometry by an axial magnetic field and an electrostatic potential. The equations describing the motion are isomorphic to those of two-dimensional fluid flow1,2 in the limit that the length of the plasma is much greater than its radius. In the case of the non-neutral plasma, the particles are considered to bounce longitudinally while drifting azimuthally. Assuming a longitudinal wave number kz 0, the fundamental equations of the Drift Poisson Model then become in the infinite length approximation , t u 0, 1 u z B , 2 2 0 . 3 The equations can be linearized assuming 0 1ei(m t), 4 0 1ei(m t), 5 u u0 u1ei(m t), 6 leading to the diocotron mode equation, m 0 r 1 r d dr r d dr 1 m2 1 r2 mq 0Br n0 r 1 0. 7 For m 1 this equation has the remarkably simple discrete mode solution, 1 r 0 r , 8 with corresponding eigenfrequency 0(rwall) obtained from the boundary condition, 1(rwall) 0. This is the wall mode. For hollow profiles a second mode self-shielded mode also exists with frequency 0(rmax) for which 1 r 0 r , r rmax, 9 1 0,r rmax , 10 where 0(r) is the equilibrium rotation frequency profile and rmax is the radius at which the profile peaks. In both modes, the eigenfrequencies are real and both modes are neutrally stable in this infinite length approximation. Finn et al. identify two instability mechanisms when finite length plasma columns are considered. The first occurs when the shape of the end of the plasma is such that there is a radial variation of the equilibrium plasma length. In this case during motion there can be a compression of the plasma by the confining potential that conserves the line integrated density parallel to the magnetic field. The second mechanism is a perturbation of the plasma length when particles interact with the confining potential at the ends. Finn et al. demonstrate that both mechanisms give instability with comparable growth rates. The mode equation becomes11 m 0 r 1 r d dr r d dr 1 m2 1 r2 m n0 r 1 mq 0B n0 r L0 r L0 r 1 q 0 m 0 r n0 L0 1 , 11 where L0(r) and L0 (r) are, respectively, the equilibrium radial profile of the plasma length and its radial derivative. The functional is the first order correction to the plasma length caused by perturbations in the potential.11 To make the analysis tractable, Finn et al. approximate the equilibrium length of the plasma by a quadratic function. However, perturbations in the plasma length ( ) were implemented ignoring curvature of the ends as a simple approximation. In Sec. IV we describe results from our own code driftk that solves Eq. 11 if 0. See Appendix. FIG. 1. A Malmberg Penning trap. The axial magnetic field B confines charged particles radially and voltages V on the rings confine the non-neutral plasma longitudinally in the cavity space between the rings. 3218 Phys. Plasmas, Vol. 9, No. 8, August 2002 G. W. Mason 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 The work of Hilsabeck and O Neil and of Coppa et al. include refinements to the theory and implementation of Finn et al. Both accept an arbitrary plasma shape, use realistic axial boundary conditions, and incorporate perturbations in the plasma length self-consistently using a Green s function. Solutions are by numerical methods. Both efforts find similar growth rates for the self-shielded mode that are several times smaller than the comparable test case taken from experiment. Calculated real frequencies of the mode are slightly smaller than the maximum of the rotation profile in contrast to experimental values that are reported to be as much as 25% lower than the maximum of the profile.6,13 Hilsabeck and O Neil conclude that quantitative agreement with the measured growth rates and frequencies requires the inclusion of a kinetic effect which arises from the experimental method used to load the hollow density profiles. The experimental protocol is assumed to truncate the high-velocity tails of the longitudinal velocity distribution in a radially dependent way. This can be important because the fast particles penetrate into a region in the ends of the plasma where their 0(r) is reduced compared to slower particles that do penetrate so deeply. If the longitudinal velocity distributions have radial dependence, the dynamics of the mode can be altered. Hilsabeck and O Neil introduce the kinetic correction by linearizing the Hamiltonian for the system and keeping a term of order D 2 in the Debye length. They conclude that plasma shape and kinetic effects together could explain the discrepancy with experiment. III. SIMULATIONS In this paper we take a numerical approach by doing particle-in-cell simulations. The method has the advantage of incorporating realistic boundary and end conditions in detail while also providing diagnostic information about the plasmas that are otherwise unknown in the experiments or in the methods of Finn et al. and of Hilsabeck and O Neil. We perform numerical experiments with the intent of helping to understand whether plasma shape and/or kinetic effects are adequate to predict experimentally measured growth rates for the unstable m 1 diocotron mode in finite-length plasmas. Azimuthally symmetric equilibria are computed separately using a two-dimensional (r z) nonlinear Poisson solver.16 The two-dimensional density array is passed from the equilibrium code to the simulation code and interpolated onto a three-dimensional Cartesian grid. The threedimensional density is then represented by particles-in-cells PIC . The azimuthal symmetry of the distribution of particles is broken by density corrections or by small initial displacements of each particle chosen to seed a particular azimuthal mode using the infinite-length theory as guide for the mode shapes. In the present work, the plasma is typically represented by about 106 computational particles that, in turn, each represent several thousand plasma electrons. The computation is done in three-dimensional Cartesian geometry into which is embedded the confining cylinder. The grid used was Nx Ny Nz 65 65 129 for plasmas that were typically 0.30 m in length with a Debye length of 0.003 m. Shortlegged differential operators for the Laplacian operator are used at the cylindrical boundary so that the cylindrical shape is treated realistically. Likewise, boundary conditions are implemented realistically, with a grounded cylinder sandwiched between confinement rings at each end held at sufficient potential to confine the plasma see Fig. 1 . Beyond the rings longitudinally and away from the plasma on each end is a short buffer zone of grounded cylinder at the end of which periodic boundary conditions ( / z 0) are maintained to complete the boundary conditions for the computation region. Poisson s equation is solved by distributing density to the computational grid and using a three-dimensional multigrid algorithm15 to solve Poisson s equation. Particles are moved in the (x,y)-plane assuming E B drift motion and using a predictor-corrector algorithm. In the longitudinal z-direction we use Newton s Second Law and a leapfrog algorithm. Densities are distributed to the grid, fields are computed from Poisson s equation, particles are moved in response to the fields, and new densities are computed to begin the cycle anew. The total kinetic energy of the particles is monitored to ensure that the algorithm conserves energy throughout the computation and the Courant condition ( p t 1) is monitored to ensure that the code is numerically stable. The code was found to be very stable and once a somewhat optimized set of convergence parameters for the multigrid algorithm were chosen, the basic stepping code performed without further attention. For the test cases we describe in this paper, a 3 ns time step was typically used to follow about 10 30 cycles of the m 1 mode of interest. The code was run on an SGI Origin 2000 or IBM SP computer depending on computer resources required. Each simulation to determine a growth rate typically took about a week. IV. RESULTS In this section we describe results of several simulation calculations. We first describe a simulation of a test case of Finn et al. We then consider the effect on growth rate of changing the depth of the hollowness in the density profile. Finally, we shift to a second family of equilibria and investigate the effect on growth rate of creating plasmas with non-Maxwellian distributions of longitudinal velocity as is actually done in the experiments. We have considered two families of equilibria. The first is based on a test case used by Finn et al.11 The electron plasma had a radius rp 0.02 m, confined within a cylinder of radius rw 0.038 m. The magnetic field was 375 G and the confining ring potentials were 50 V. The central cylinder had a length of 0.32 m. The rings had a width of 0.03 m and in this case the length of the buffer zone was zero. The radial density profile of the plasma in the Finn theory was given by the parametrization, n0 r n0 0 1 r/rp 2 2 1 2 r/rp 2 12 for r rp and zero elsewhere. The radial profile of the length of the plasma was parametrized by Finn et al. as L0 r L0 0 1 r/rw 2 , 13 Phys. Plasmas, Vol. 9, No. 8, August 2002 Simulations of the instability of the m 1 . . . 3219 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 rw is the radius of the cylinder. The hollowness of the profile was controlled by the parameter and the curvature of the ends was described by the parameter . In the test case computed by Finn et al., was chosen to be 3, resulting in a ratio nmax /n0 1.28 and a value for of 0.25. We prepared a simulation using the same radial density profile, plasma radius, cylinder radius, confining ring potentials, and magnetic field as the Finn et al. test case. We chose a plasma temperature of 1.2 eV and a plasma length of 0.30 m. Under these conditions the value of from our equilibrium code was approximately 0.3 and the value of nmax was 6.28 1012m 3 corresponding to 0(max) 1.44 106 s 1. Plasma length profiles for the series of simulations based on the Finn et al. test case are shown in Fig. 2. Observe that the vertical scale is truncated. The radial variation of L0 is relatively small compared to the overall length of the plasma. The profiles in Fig. 2 do not appear to be strictly parabolic, so our value of kappa 0.3 is only a rough approximation to the Finn value of 0.25. At each time step of the simulation a longitudinally lineintegrated density function was formed and then Fourier analyzed to find the amplitude and phase of the m 1 mode. The phase signal as a function of time was used to measure the real frequency of the mode and an exponential function was fitted to the amplitude signal as a function of time to obtain the growth rate. The amplitude signal rises with an apparent exponential growth, but the signal wobbles slightly with time relative to the exponential. We have indicated with error bars an estimate of the uncertainty that this gave to the growth rate measurements. The wobbles may arise from interference in the amplitude signal with the m 1 continuum which would be expected to be present if our initial perturbation seed is not quite right. Figure 3 shows a typical amplitude signal with a corresponding exponential growth for comparison as well as the corresponding phase signal. The results of the simulations are shown in Fig. 4. The ratio of growth rate to real mode frequency is about 0.007 compared to the Finn et al. result of 0.009, the Coppa et al. result of 0.008, and the experimental value of 0.025. Our real frequency was 1.40 106 s 1 which is consistent with the expectation of a real frequency near the maximum of the radial rotation profile of our equilibrium (1.44 106 s 1). In each case in Fig. 4, the simulation gives a real frequency of the mode to within 1% 3% of the maximum of the rotation profile 0 max in contrast to the 25% reduction reported by Kabantsev. We completed the first family of simulations by maintaining the value of n0(max) and deepening the profile using FIG. 2. Radial length profiles for the plasmas identified with asterisks in Fig. 4. A parabolic fit (L(x) L0(1 (r/rw)2)) to the uppermost curve yields roughly a value 0.3. As nmax /n0 increases, the profiles are less and less parabolic. FIG. 3. Upper The phase signal in radians from which the real frequency r of the mode is determined. The phase signal from the simulation is remarkably linear. Lower An exponential growth curve compared to the amplitude signal for the m 1 mode from the particle-in-cell simulations. Error flags on the simulation results plotted in Fig. 4 are estimates of the uncertainty in growth rate determinations based on figures like this one. FIG. 4. Comparisons of growth rates to the theory of Finn et al. and an experimental test case open circles . Asterisks * mark simulation growth rates using Finn s density formula with 3.0, 5.66, 8.19, and 15.07. The corresponding ratios nmax/n0 are 1.28, 1.64, 2.00, 3.00. Hyphens - mark simulation points obtained by depleting the tails of the longitudinal velocity distributions. Corresponding Maxwellian cases are marked by x for comparison. Also included for comparison are predictions from our drift-kinetic code calculation open squares and results of the complete model of Coppa et al. marked with open diamonds. 3220 Phys. Plasmas, Vol. 9, No. 8, August 2002 G. W. Mason 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 the Finn et al. radial density parameterization with 5.66, 8.19, 15.07. The corresponding ratios nmax /n0 are 1.64, 2.00, 3.00. These growth rates are also shown in Fig. 4 as the points marked with the asterisk * symbol. The point at nmax /n0(0) 2 is coincident with the prediction of our drift-kinetic code see Appendix . The second family of simulations was done to address the possibility suggested by Hilsabeck and O Neil that the persistent discrepancy between experiment and theory may be a kinetic effect when particles of differing energies penetrate the confining potential at the ends to differing degrees. This effect may be enhanced because, experimentally, the hollow profiles are created from a nonhollow profile by temporarily lowering the end potentials. This procedure allows more particles to escape near r 0 than at larger radii, and also depletes the Maxwellian distribution of velocities in a radially dependent way. The non-Maxwellian distributions created in this way will depend heavily on the protocol used to create the plasmas and experimental velocity distribution data are not available for a specific test case. For purposes of simulation, we began with a flat-topped density profile. We again used a magnetic field of 375 G and a temperature of 1.2 eV. However, the radius of the cylinder was 0.05 m and the plasma length about 0.35 m. The central density plateau was 5 1012 m 3. The confining ring potentials were 200 V. The central cylinder had a length of 0.44 m, rings a width of 0.03 m, and a buffer zone length of 0.05 m. The equilibrium was calculated and loaded into the simulation code as before. However, when the code began, the ring voltages boundary conditions were reduced linearly over a 1 s period of time to some fraction of the original confining value chosen to reduce the potential near r 0 to a value close to the central potential of the plasma. The rings were then held down at this destination potential for 4 s, then linearly raised again to the original value over a final microsecond. The resulting density profile was hollow and not unlike the profiles obtained from the Finn et al. formula. The degree of hollowness nmax /n0(0) was controlled by the fraction applied to the 200 V ring potentials when the confining potentials were lowered. Since the potential remained down for several bounce times, virtually all particles with longitudinal velocities below a certain critical value were removed. The velocity distributions resulting from this protocol had a radial dependence and since equilibration times were much longer than the time of our simulation, the resulting velocity distributions as a function of radius were and remained non-Maxwellian through the course of the simulation. The hollowed density profiles created in this way are shown in Fig. 5. Figure 6 shows the corresponding root-mean-square velocity profiles. The velocity distributions themselves at each radius were not Maxwellian, but appeared to be approximately a combination of a core Maxwellian of one temperature, with a tail of different temperature. In the case of the least hollow plasma, the tail was cooler than the core distribution i.e., truncated . However, for the most hollow profile, the tail is warmer for the smallest radii, then crosses over to become cooler for a radius about halfway up the slope of the density hole, and finally becomes identical to the core temperature at the edge of the density hole. Thus, the radial velocity distributions can depend sensitively on the protocol that produces them. Three plasmas were created using the protocol described above with nmax /n0 1.5, 3.2, 6.2. The corresponding density profiles are shown in Fig. 5. The growth rates are shown in Fig. 4 marked with the hyphen - symbol. The 6.2 point is labeled as non-Maxwellian. For comparison, we took the radial density profiles from these three non-Maxwellian simulations and created distributions differing only in that the longitudinal velocity distributions were Maxwellian at all radii with temperature 1.2 eV . These are shown in Fig. 4 marked with the symbol. The non-Maxwellian growth FIG. 5. Simulated hollow line-integrated density profiles obtained by temporarily lowering the potentials on the confinement rings according to the protocol described in the text. Absolute line integrated densities at r 0 are 1.10 1012 bottom , 5.12 1011 middle , and 2.52 1011 top particles per square meter. FIG. 6. Radial root-mean-square velocity (vz) profiles for the simulated plasmas of Fig. 5. The corresponding velocity distributions at each radius are identified as non-Maxwellian in Fig. 4. Phys. Plasmas, Vol. 9, No. 8, August 2002 Simulations of the instability of the m 1 . . . 3221 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp rates are enhanced over their Maxwellian counterparts by factors of 1.9, 1.3, and 1.5 for nmax /n0 equal to 1.5, 3.2, and 6.2, respectively. To complete this second family of numerical experiments, we considered an alternative protocol for forming the hollow profile. We began with a markedly peaked equilibrium taken from Fig. 16 in the paper by Hilsabeck and O Neil.13 We then lowered the ring potentials for 2.5 s to a potential about equal to the central potential of the equilibrium. The line-integrated-density profiles before and after hollowing are shown in Fig. 7. The root-mean-square longitudinal velocity profile is shown in Fig. 8. The relative growth rate for this non-Maxwellian case is shown at a value of nmax /n0 1.3 in Fig. 4. Again, kinetic effects, to the degree that the protocol used corresponds to the experimental protocol, fail to remove the discrepancy with experiment. The growth rate is shown in Fig. 4 marked with a hyphen - , but lacking a crosspiece on the error bars. As an approximate check of the results of our simulations and to provide an additional comparison to results by Finn et al. and Coppa et al., we have used a separate linear drift-kinetic eigenvalue code see Appendix to compute growth rates. The code, which was originally written for the infinitely-long plasma approximation, was easily modified to include the L0 /L0 term on the right-hand side of Eq. 11 but not the term. The eigenvalue code uses equilibria calculated separately by the same code used to calculate equilibria for the simulations using density profiles from Eq. 12 see Fig. 4 . This independent calculation reproduces the Finn et al. growth rate almost exactly at nmax /n0 1.28, coincides exactly with the simulation point at nmax /n0 1.64, and meshes smoothly with the results of Coppa et al. using their complete model. The results of our eigenvalue code are shown as open squares in Fig. 4. The results of Coppa et al. are shown as open diamonds. Finally, to further check the performance of the code and our results, we repeated three simulations of the Finn test case with a doubled magnetic field, then with a more refined computational grid in the xy-plane, and finally with both the stable wall and unstable self-shielding modes seeded simultaneously. Doubling the magnetic field gave a value / r 0.007, which is virtually identical to the value of Fig. 4. This suggests that the bounce averaging of the longitudinal motion at the lower magnetic field is being done adequately. Doubling the computational grid resolution in the xy-plane from 65 65 129 to 129 129 129 increased / r to 0.008 which is just beyond the upper flag limit in Fig. 4. This value does not qualitatively change any of our results but does put us closer to the Finn and Coppa calculations. However, it is presently impractical for us to do extensive investigations at this resolution. Finally, we seeded both m 1 modes simultaneously at about the 1% level. In this case / r 0.008, i.e., no significant change in the growth rate of the unstable mode. V. CONCLUSIONS We have used particle-in-cell simulations to compute growth rates for the hollow, finite-length m 1 self-shielded diocotron mode. We have investigated test cases where the persistent discrepancy between theory and experiment may be a consequence of the shape of the ends of the plasma or of kinetic effects arising from non-Maxwellian velocity distributions introduced in the experimental preparation of the plasmas. In none of the test families were we able to achieve growth rates as large as the experimental value. Real frequencies computed in the simulations for the self-shielding mode were typically within about 3% of the maxima of the respective rotation profiles, in contrast to reported experimental reductions of 25%. Figure 4 summarizes our results. Our simulation of the growth rates for the experimental test case UCSD is a fac- FIG. 7. Original and simulated hollow line-integrated density profile for an alternative formation protocol described in the text. The hollow profile was obtained from the original shown in the figure by lowering the confinement potentials to about the central potential of the plasma for about 2.5 s. The curves are normalized to a line-integrated density of 1.20 1012m 2 at r 0. This non-Maxwellian case is shown with a hyphen - symbol for nmax /n0 1.3 in Fig. 4. FIG. 8. Radial root-mean-square velocity (vz) profile for the alternative formation protocol of Fig. 7. 3222 Phys. Plasmas, Vol. 9, No. 8, August 2002 G. W. Mason 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 tor of 3 lower than the experiment, but slightly less than the theoretical analyses of Finn et al. and Coppa et al. Observe, however, that the comparisons are made in a region of very steep dependence on nmax /n0(0). The simulation growth rates increase by a factor of about 2 as nmax /n0(0) increases to more than 6, but never reach the experimental value. Thus, we agree with the conclusion of Hilsabeck and O Neil that the end shape effects we considered are alone insufficient to remove the discrepancy with experiment. We have also simulated a non-Maxwellian longitudinal velocity distribution effect suggested by Hilsabeck and O Neil. The effect increases growth rates by 30% 90% compared to Maxwellian control cases. Nevertheless, the simulation is still a factor of 2.0 too low to remove the discrepancy with experiment at nmax /n0(0) 1.28. However, the two methods which we used to create the non-Maxwellian distributions may not correspond to the actual experimental protocols, leaving some room yet for further study. ACKNOWLEDGMENTS We gratefully acknowledge helpful discussions with J. M. Finn, D. del Castillo-Negrete, T. M. O Neil, D. H. E. Dubin, T. J. Hilsabeck, and A. A. Kabantsev. S. N. Rasband provided helpful results from a recently-developed cold-fluid code against which we compared test runs from the code described in this paper. We also gratefully acknowledge the use of facilities of the Brigham Young University Supercomputing Center. APPENDIX: DRIFT-KINETIC CALCULATION We solve Eq. 11 with the modified density derivative term of Finn and Castillo by a matrix shooting code called driftk . The code and underlying theory will be briefly discussed here. The code solves for linear modes in an infinitely long cylinder using the drift-kinetic equation f t vD f v f z q mp z f v 0, A1 where mp is the particle mass, v is the velocity in the z direction, and the drift velocity is vD z B . A2 Linearizing this equation according to r,z, 0 r 1 r exp im ikz i t , A3 f r,z, ,v f 0 r,v f 1 r exp im ikz i t , A4 with f 0 r,v n0 r vth 2 exp v vb r 2/2vth 2 , A5 where vb(r) is the beam velocity profile of the plasma if any and vth is the thermal velocity, then results in f 1 r,v m rB f 0 r kq mp f 0 v 1 m 0 r kvb r . A6 Using this expression for f 1 in Poisson s equation, 2 1 q 0 f 1dv A7 yields the following electrostatic drift-kinetic mode equation: 1 r d dr r d 1 dr m2 r2 1 k2 1 q2n0 r 0mpvth 2 W m 0 r kvb r kvth 1 mq 0Br m 0 r kvb r d dr n0 r W m 0 r kvb r kvth 1 1 0, A8 where W(z) is the plasma dispersion function using the notation of Ichimaru,14 W z 1 2 xe x2/2 x z dx. A9 Equation A8 is solved numerically by finite differencing on a radial grid to obtain a set of linear equations for 1 at the radial grid points. At the ends of the radial interval we use the boundary conditions d 1 dr 0 0 if m 0; 1 0 0 if m 1; 1 rwall 0. A10 Since the set of linear equations is homogeneous we must avoid the trivial solution 1 0 and this is done by requiring 1 rc 1 A11 for some arbitrary choice of critical radius rc between the origin and the conducting wall. Equation A11 replaces the finite-difference equation at the grid point corresponding to rc and the replaced difference equation is saved for later use to determine the mode frequency. The general result of solving the linear system with the replacement equation A11 consists of a function with a discontinuous derivative at r rc . The mode frequency is determined by adjusting its value until the original finite difference equation at r rc is satisfied, resulting in a smooth eigenmode. This adjustment is made by refining an initial guess using Powell s method as implemented in Numerical Recipes.15 For the unstable mode studied in this paper we set k 0, and m 1. For k 0 the plasma dispersion function vanishes. To include at least part of the physics discussed by Finn and Castillo, we make the replacement in Eq. A8 , Phys. Plasmas, Vol. 9, No. 8, August 2002 Simulations of the instability of the m 1 . . . 3223 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp d dr n0 r 1 L0 r d dr n0 r L0 r A12 to obtain m 0 r 1 r d dr r d 1 dr m2 r2 1 mq 0Br 1 L0 r d dr n0 r L0 r 1 0. A13 Thus, although the code driftk generally solves a much more general problem, the drift-kinetic mode equation reduces to Eq. 11 when k 0 and m 1. For each case that we study we run a separate equilibrium code16 to obtain the plasma length as a function of radius L0(r).We then fit a polynomial in r to the numerically determined L0(r) and use driftk as described above to solve for and 1(r). 1R. H. Levy, Phys. Fluids 8, 1288 1965 ; 11, 920 1968 . 2R. J. Briggs, J. D. Daugherty, and R. H. Levy, Phys. Fluids 13, 421 1970 . 3R.Copyright Davidson, Theory of Non-neutral Plasmas Benjamin, Reading, 1974 . 4C. F. Driscoll and K. S. Fine, Phys. Fluids B 2, 1359 1990 . 5C. F. Driscoll, Phys. Rev. Lett. 64, 645 1990 . 6A. Kabantsev and C. Driscoll, in Non-Neutral Plasmas III, edited by J. Bollinger, R. Spencer, and R. Davidson American Institute of Physics, New York, 1999 , pp. 208 213. 7R. A. Smith and M. N. Rosenbluth, Phys. Rev. Lett. 64, 649 1990 . 8R. A. Smith, Phys. Fluids B 4, 287 1992 . 9S. N. Rasband, R. L. Spencer, and R. R. Vanfleet, Phys. Fluids B 5, 669 1993 . 10S. N. Rasband, Phys. Plasmas 3, 94 1996 . 11J. M. Finn, D. del Castillo-Negrete, and D.Copyright Barnes, Phys. Plasmas 6, 3744 1999 . 12G. G. M. Coppa, A. D Angola, G. L. Delzanno, and G. Lapenta, Phys. Plasmas 8, 1133 2001 . 13T. J. Hilsabeck and T. M. O Neil, Phys. Plasmas 8, 407 2001 . 14S. Ichimaru, Basic Principles of Plasma Physics Benjamin, Reading, 1973 , Chap. 4. 15W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran, 2nd ed. Cambridge University Press, New York, 1992 , Chap. 10. 16R. L. Spencer, S. N. Rasband, and R. R. Vanfleet, Phys. Fluids B 5, 1738 1993 . 3224 Phys. Plasmas, Vol. 9, No. 8, August 2002 G. W. Mason 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