Computing the m=1 diocotron frequency via an equilibrium calculation in non-neutral plasmas Ross L. Spencer Department of Physics and Astronomy, Brigham Young University, Provo, Utah 84602 (Received 19 May 2004; accepted 12 August 2004; published online 22 October 2004) The m=1 diocotron mode in non-neutral plasmas has long been thought of as a shifted equilibrium, and its frequency has been approximately calculated in this way by Fine and Driscoll [Phys. Plasmas 5, 601 (1998)]. This article shows that this idea can be coupled with a standard axisymmetric equilibrium calculation on a grid to calculate the frequency of this mode to very high precision including both finite-length and thermal effects, provided that the Debye length is small enough. As the Debye length begins to approach the plasma size not only does the shifted equilibrium calculation fail to predict correctly the frequency of the mode, but the idea that the mode is a simple shift of the original equilibrium also becomes invalid. 2004 American Institute of Physics. [DOI: 10.1063/1.1803840] A conceptually straightforward way to calculate the diocotron mode frequency in non-neutral plasmas that has been investigated approximately by Fine and Driscoll,1 is to treat the mode as a shifted equilibrium (m=1 perturbation). In this scheme the plasma is shifted sideways in the x direction by an amount D and then both the perturbed density n1 r, z cos and perturbed electrostatic potential 1 r, z cos are calculated. The plasma and external/ induced-charge potentials are then separated in the way given below: 0 r,z = 0p r,z + 0e r,z , 1 r,z = 1p r,z + 1e r,z , 1 where the subscript e indicates the external potential (due both to equilibrium confining fields and induced-charge fields) and the subscript p indicates the potential of the plasma alone. The external components of the field are then allowed to act on the plasma and the first-order contributions to the electric force are integrated over the volume of the plasma to find the net force on the shifted plasma: F1 = q n1 0e + n0 1e dV. 2 In this integral only the component of the force in the direction of the shift will survive (the x component in this case) so the diocotron mode frequency can be computed from the net E B drift of the shifted plasma from D = F1x NqBD , 3 where N is the number of particles in the plasma, q is the charge on each particle, and B is the magnitude of the uniform confining magnetic field. (This equation is essentially the same as Eq. (5) in Fine and Driscoll.1) This conceptually simple method is not so straightforward in computational practice because it is difficult to separate the plasma and external components of the field when the field is computed on a grid. (Fine and Driscoll did this analytically by means of a Green s function technique.) It would be better if there was a method for finding F1x due to external sources that used the full field on a computational grid, for then the m=1 diocotron frequency could be computed accurately as part of a standard non-neutral plasma equilibrium calculation.2 A method for carrying out this calculation is the subject of this paper. To find this method we simply write down the obvious full-field extension of Eq. (2), then separate the external and plasma parts, F1 = q n1 0 + n0 1 dV, 4 or F1 = q n1 0e + n0 1e dV q n1 0p + n0 1p dV. 5 If the second integral were to vanish then we would have F1 = F1 and the diocotron frequency could be calculated from the full field on the grid. That the second integral is zero is obvious on physical grounds. It is the net self-force on the plasma after it has been shifted and allowed to come into equilibrium, and the plasma cannot exert a net force on itself, as pointed out by Fine and Driscoll.1 Of course, having the integral vanish on physical grounds and having it vanish on a computational grid are not the same thing, but it will be seen later in this paper that ignoring it leads to accurate results. The computation of the diocotron mode frequency from an equilibrium calculation on a grid proceeds as follows: (a) First compute an axisymmetric equilibrium.2 (b) Mathematically shift the equilibrium density in the x direction and compute the perturbed density n1 r, z cos and perturbed potential 1 r, z cos for this shifted equilibrium using a perturbed equilibrium calculation on an r-z grid. Copyright Finally, PHYSICS OF PLASMAS VOLUME 11, NUMBER 11 NOVEMBER 2004 1070-664X/2004/11(11)/5356/4/$22.00 5356 2004 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 compute the gradients of 0 and 1 and use them in Eqs. (3) and (4) to find the diocotron frequency. Each of these steps will now be discussed in detail. (a) Solve the nonlinear Poisson equation that describes a non-neutral plasma equilibrium.2 (b) As a first approximation to the solution of the perturbed equilibrium problem it is possible simply to shift the plasma in the x direction, resulting in a perturbed density given by the formula n1 = Dx n0 = D n0 r cos . 6 Given this density, the m=1 Poisson equation 1 r r r 1 r + 2 1 z2 1 r2 = qn1 0 7 can then be solved. [Note that the factor cos has been suppressed for both n1 and 1 in Eq. (7).] Numerical experiments show that this method actually comes reasonably close to the correct answer, but because it does not allow the perturbed plasma to relax to axial thermal equilibrium, it is not yet correct. To obtain a more accurate result recall that the general form for the density in equilibrium is2 n r,z = nmid r exp q r,z r,0 /kT , 8 where nmid r is the radial density profile in the plasma midplane. In the shifted equilibrium the perturbed density, resulting from small changes in nmid and is given by n1 r,z, = n1mid r exp q 0 r,z 0 r,0 /kT q 1 r,z 1 r,0 kT n0 r,z , 9 which may then be used in Poisson s equation (7). The unknown profile n1mid r can be determined by noting that when the plasma is shifted and then relaxed axially, the relaxation takes place along the magnetic field lines, which preserves ndz. So n1mid r can be determined by requiring that it produce a thermally relaxed perturbed density with the same profile of n1dz as that of a rigid shift n1dz = D n0 r dz 10 at each radius. With this specification of n1, Eqs. (7), (9), and (10) form a nonlinear system of equations for 1 (through the constraint on n1mid) which can be solved by an iteration process similar to that used for unshifted equilibrium problems.2 Note, however, that the difficulties which arise in the general equilibrium problem because q /kT appears as the argument of the exponential function are avoided here because we are linearizing in small 1. This makes Eq. (7) almost linear, except for the relatively minor adjustment required on n1mid r . The iteration technique used for determining n1mid r is simply to solve Eqs. (9) and (10) for n1mid * at each step of the electrostatic relaxation and substitute it back into Eq. (7) for the next step. It is found that underrelaxing on this iteration is usually necessary: n1mid m+1 = n1mid * + 1 n1mid m , 11 where m denotes iteration level and 0.1. Copyright With the perturbed density and potential in hand the diocotron frequency can be computed from Eqs. (3) and (4) as follows. First calculate the x derivatives of the equilibrium and perturbed potentials 0 x = cos 0 r 12 and 1 x = cos2 1 r + sin2 1 r . 13 Then substitute these results into Eq. (4) to obtain F1x = q n0 r,z n1 r,z 0 r + 1 r + 1 r rdrdz, 14 where the integration in has already been performed in the volume integral. Numerical experiments using a three-dimensional (3D) drift-kinetic simulation4 show that this method works very well (better than 1% accuracy) as long as the Debye length is small, as shown in Fig. 1. In the simulations an axisymmetric equilibrium was loaded into the simulation program after which gradually, over the first few hundred time steps, the plasma was shifted off axis by a total of about 3% of the FIG. 1. The ratios of computed diocotron frequencies to results from a 3D time-dependent simulation are displayed. The open circles are for the numerically computed perturbed equilibria described in this paper. The crosses are from the approximate analytic calculation of Fine and Driscoll (Ref. 1). Phys. Plasmas, Vol. 11, No. 11, November 2004 Computing the m=1 diocotron frequency 5357 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp conducting wall radius. (The time step in the simulation is restricted by a Courant condition on the plasma frequency, so typically p/2.) After the shift was finished the plasma was allowed to evolve freely and the position of the centerof- mass was monitored. When the run was finished the center-of-mass data were analyzed to find the instantaneous angular frequency, which is the diocotron frequency. Examples of the time evolution of this frequency are displayed in Fig. 2 and the simulation parameters used to produce Fig. 1 are displayed in Table I. In all cases the equilibrium profiles were roughly flattops with thermal edges (widths on the order of the Debye length) in both r and z. The relevant parameter involving the Debye length that determines whether the diocotron frequency can be accurately calculated by the shifted-equilibrium method seems to be the parameter of Peurrung and Fajans:3 D 2 RpLp , 15 where D is the Debye length, Rp is the plasma radius, and Lp is the full plasma length. Figure 1 shows a comparison between the perturbed equilibrium calculation, the approximate calculation of Fine and Driscoll, and results from the 3D simulation for 22 different equilibria with various radii, lengths, and profiles. As long as is below about 0.01 the perturbed equilibrium calculation is accurate to about 1% and the approximate method of Fine and Driscoll is accurate to about 3% (about the level of accuracy quoted in their paper when they checked their calculation against experimentally measured frequencies). But for larger values of neither calculation works well. The reason for the inaccuracy of the shifted-equilibrium method when is not small is, as pointed out by Peurrung and Fajans, that particles of different axial kinetic energy have different rotation frequencies due to the way they turn around axially in the end of the thermal plasma. They saw significant smearing effects in end-on pictures of the plasma shape as approached unity, but this same mechanism apparently affects the shifted-equilibrium calculation at much smaller values of . Figure 2 shows simulation traces of the instantaneous angular frequency of the plasma center-ofmass for two plasmas in the data set of Fig. 1, one with =0.0024 and the other with =0.081. Notice in the case of the larger that a rigid shift seems not to be the correct perturbation to produce a mode initially, but that things settle down later on. For small , however, an excellent mode is obtained. Figure 3 shows the result of dividing the particles in the simulation into four different axial energy groups and monitoring the instantaneous x position of the center-of-mass of each group. Notice that for small each group behaves in the same way, indicating that a simple shift produces a mode. But for the case of larger the four groups are disorganized at first. Only later on do they synchronize, indicating that a diocotron mode has been established. TABLE I. The plasma parameters used in the simulations displayed in Fig. 1 are listed in order of increasing . The central density is n0, the temperature is T, the midplane plasma radius (half-density point) is Rp, the plasma half length (half-density point) is Zp, the magnetic field is B, the central rotation frequency is 0, and the simulation diocotron frequency is D. In all cases the plasma was contained by a grounded cylinder of radius 4 cm with end rings of length 3 cm set at a confining voltage of 100 V. The density profile was nearly flattopped with thermal edges in r and z a few Debye lengths wide. n0 1012/m3 T (eV) Rp (cm) Lp (cm) B (T) 0 106/s D 106/s 3.00 0.10 2.00 16.7 0.01 0.000 555 2.71 0.768 1.00 0.10 2.00 23.3 0.01 0.001 19 0.905 0.251 1.00 0.10 2.00 20.6 0.01 0.001 34 0.905 0.254 3.00 0.10 2.00 6.55 0.01 0.001 41 2.71 0.931 3.00 0.10 2.00 4.50 0.01 0.002 05 2.71 1.030 0.10 0.01 2.00 12.9 0.0025 0.002 14 0.362 0.109 1.00 0.10 2.00 11.7 0.01 0.002 36 0.905 0.276 1.00 0.10 2.00 10.8 0.01 0.002 55 0.905 0.280 1.00 0.20 2.00 10.8 0.01 0.005 10 0.905 0.285 1.00 1.00 2.00 23.5 0.01 0.0118 0.905 0.305 1.00 0.40 1.62 10.8 0.01 0.0126 0.905 0.245 1.00 1.00 1.91 11.4 0.01 0.0254 0.905 0.298 1.00 1.00 1.77 12.1 0.01 0.0259 0.905 0.279 1.00 2.37 1.77 24.2 0.01 0.0306 0.905 0.270 1.00 4.33 1.77 25.3 0.01 0.0534 0.905 0.295 1.00 2.37 1.77 13.1 0.01 0.0567 0.905 0.316 1.00 6.87 1.77 26.4 0.01 0.0813 0.905 0.318 1.00 4.00 1.91 13.3 0.01 0.0874 0.905 0.377 1.00 4.33 1.77 14.2 0.01 0.0954 0.905 0.352 1.00 4.00 1.62 13.2 0.01 0.103 0.905 0.353 1.00 10.0 1.77 27.4 0.01 0.114 0.905 0.340 1.00 10.0 1.77 16.2 0.01 0.193 0.905 0.416 5358 Phys. Plasmas, Vol. 11, No. 11, November 2004 Ross 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 Note. One of the workers suggested that the rigidity parameter (axial particle bounce frequency divided by the central rotation frequency) might also be relevant to this study. Unfortunately, nearly all of the simulations had the same values for the magnetic field and the density, so rigidity and increased and decreased together as the temperature was varied. To explore this issue at least somewhat, several more simulations were run at both small and large in which everything was kept the same except the magnetic field, which varies the rigidity while keeping constant. It was found that for below about 0.01 decreasing the rigidity by a factor of 10 (down to values near 1, or below) caused only small changes in the ratio of the simulated diocotron frequency to that obtained from the equilibrium calculation (about 1%, which is the level of accuracy for the cases studied here). But with near 0.1, decreasing the rigidity down to values near 1 causes a larger shift in this ratio. (In the case with =0.0813 when B was changed from 0.01 to 0.00033 T the simulation frequency changed from 3.18 105 to 1.02 107 / s, a ratio shift of 6%.) So the rigidity seems to matter somewhat when the equilibrium prediction does not work well, but is of little significance when it does. In conclusion, it is found that a relatively simple perturbed equilibrium calculation can be added to the standard axisymmetric equilibrium calculation used by many researchers in non-neutral plasma physics to compute the m =1 diocotron mode frequency, provided that the Debye length parameter = D 2 / RpLp introduced by Peurrung and Fajans is on the order of, or smaller, than 0.01. Such equilibrium calculations are already used to help diagnose experimentally produced plasmas, and this extension provides an additional tool for better understanding of these plasmas. 1K. S. Fine and C. F. Driscoll, Phys. Plasmas 5, 601 (1998). 2R. L. Spencer,S. N. Rasband, and R. R. Vanfleet, Phys. Fluids B 5, 4267 (1993). 3A. J. Peurrung and J. Fajans, Phys. Fluids B 5, 4295 (1993). 4G. W. Mason, Phys. Plasmas 10, 1231 (2003). FIG. 2. Simulation results for the angular velocity of the plasma center-ofmass are displayed. In the upper panel the Debye length is small, with =0.0024, while in the lower panel it is larger, with =0.081. Both simulations were seeded with a rigid sideways displacement, and it is clear that in the lower panel this procedure does not give a clean mode. FIG. 3. The x position of the center-of-mass for four different equallypopulated energy groups for two different simulations are displayed (lowest; solid line, low; dashed, medium: dotted, high; dash-dot). When is small the four groups are synchronized for a rigid shift perturbation, and the perturbed equilibrium calculation accurately predicts the diocotron frequency. When is larger they are not synchronized at first, but gradually become synchronized after the perturbed density evolves. In this case the diocotron frequency after synchronization is not well-predicted by the rigid-shift/ perturbed equilibrium calculation (the equilibrium prediction is high by 20%). Phys. Plasmas, Vol. 11, No. 11, November 2004 Computing the m=1 diocotron frequency 5359 Downloaded 04 Mar 2009 to 128.187.0.164. Redistribution subject to AIP license or Copyright; see http://pop.aip.org/pop/copyright.jsp