PHYSICS OF PLASMAS VOLUME 5, NUMBER 10 OCTOBER 1998 Studies of global stability of field-reversed configuration plasmas using a rigid body model H. Ji, M. Yamada, R. Kulsrud, N. Pomphrey, and H. Himuraa) Princeton Plasma Physics Laboratory, Princeton University, P.O. Box 451, Princeton, New Jersey 08543 (Received 22 April 1998; accepted 9 July 1998) Global stability of field-reversed configuration (FRC) plasmas has been studied using a simple rigid body model in the parameter space of s (the ratio of the separatrix radius to the average ion gyro-radius) and plasma elongation E (the ratio of the separatrix length to the separatrix diameter). Tilt stability is predicted, independent of s, for FRC's with low E (oblate), while the tilt stability of FRC's with large E (prolate) depends on s/E. It is found that plasma rotation due to ion diamagnetic drift can stabilize the tilt mode when s/E <= 1.7. The so-called collisionless ion gyro-viscosity also is identified to stabilize tilt when s/E <= 2.2. Combining these two effects, the stability regime broadens to s/E <= 2.8, consistent with previously developed theories. A small additional rotation (e.g., a Mach number of 0.2) can improve tilt stability significantly at large E. A similar approach is taken to study the physics of the shift stability. It is found that radial shift is unstable when E < 1 while axial shift is unstable when E > 1. However, unlike tilt stability, gyro-viscosity has little effect on shift stability. © 1998 American Institute of Physics. [S1070-664X(98)03110-3] I. INTRODUCTION The field-reversed configuration (FRC) is a unique toroidal magnetic confinement scheme in that there is no appreciable toroidal field. The plasma is confined purely by a poloidal field, which is produced by a toroidal plasma current. Thus the current flows in the direction perpendicular to the local magnetic field, sustaining maximum possible plasma beta close to unity. On the other hand, due to the lack of a center conductor and a confining toroidal field, FRC's are predicted to be unstable to many global magnetohydrodynamic (MHD) modes. However, FRC plasmas formed in theta-pinch devices exhibit remarkable global stability with a few exceptions. Much theoretical effort has been made to reveal stabilizing mechanisms of the predicted instabilities (the tilt mode in particular), including effects from plasma rotation, two-fluid, ion finite Larmor radius (FLR), energetic ions, and current profile. Although agreement between theory and experiment has improved over the years, few concrete physical pictures of stabilizing mechanisms have been given. In this paper, a simple equation of motion for each global mode is formulated and analyzed using a rigid body model of the FRC plasma. The strategy taken here is to elucidate semiquantitatively the essential physics for stabilizing mechanisms by using the simplest possible equations. Although the deduced marginal stability condition may not be sufficient due to the limited degrees of freedom of rigid body motion, the analyses described below should shed new light in understanding the fundamental physics of FRC stability. After a brief description of FRC models in Sec. II, tilt stability is analyzed in detail in Sec. III, including effects from the j x B torque, plasma rotation due to ion diamagnetic drift, ion gyro-viscosity, and E x B rotation. In Sec. IV, axial and radial shift stability is analyzed, followed by discussions and conclusions. II. MODELS OF FRC PLASMAS A. Solovev model of FRC plasmas The global modes of a plasma are often destabilized by the j x B force, which is usually a strong function of plasma shape, e.g., plasma elongation defined by the ratio of the separatrix length to the separatrix diameter. Here j is the internal current density of the plasma and B is the vacuum field produced by external coil currents. To quantify this force, a FRC equilibrium solution with a known vacuum field is needed. The simplest analytic model of FRC equilibrium with arbitrary elongation is the Solovev's solution given by psi = psi_0 [ 1 - (R Z / (R_0 Z_0))^2 - (R^2 / R_0^2 - 1)^2 ], where psi is the poloidal flux function, R_0 is the radius of the magnetic axis, and Z_0 is defined in Fig. 1(a). As is also obvious from Fig. 1(a), the length and radius of the FRC separatrix are L = 2 sqrt(2) Z_0 and R_s = sqrt(2) R_0, respectively, resulting in an elongation E = L / (2 R_s) = Z_0 / R_0. The trapped flux 2 pi psi_0 is related to the magnetic field at the edge B_0 = B_Z(R = R_s, Z = 0) by 2 pi psi_0 = 2 pi B_0 R_s^2 / 4. When E = 1, Solovev's solution reduces to the well-known spherical Hill's vortex with an analytic external solution. The vacuum solution for arbitrary E is obtained numerically by placing coils around the plasma. The coil currents are calculated by matching flux values at the separatrix (see the Appendix for details). One such example is shown in Fig. 1(a) a)JSPS research fellow on leave from Osaka University. where the internal (external) flux is represented by solid (dotted) lines. The vacuum flux together with coil locations is separately shown in Fig. 1(b). FIG. 1. Solovev's solution. (a) The internal (external) flux is represented by solid (dotted) lines. The rectangular box indicates rigid body model. (b) The vacuum flux produced by coils. The dotted line represents the separatrix. FIG. 2. (a) Three rotating axes (X,Y,Z) of the rigid body model. (b) Rotation of the vacuum field in the plasma frame. B. Cylindrical rigid body model In order to elucidate the essential physics of FRC global stability, we use an even simpler cylindrical rigid body model [a rectangular box in Fig. 1(a)] to analyze the motion of each mode. The cylinder is filled with a plasma of uniform density n, radius R_s, and length 2Z_0 = sqrt(2) E R_s. As a result, the moments of inertia of the cylindrical plasma with respect to each axis (see Fig. 2) are given by I_X = I_Y = (M_total R_s^2 / 4) * [ 1 + (2/3) E^2 ], I_Z = M_total R_s^2 / 2, where M_total = sqrt(2) pi R_s^3 E n m_i is the total mass (m_i = ion mass). One question which could arise is whether the global modes in this model are internal or external. An apparent answer is that they are external. However, we note that qualitatively they can be "internal" if the modeled region is smaller than the separatrix, i.e., R < R_s. In this respect, the nature of the analyzed modes can include internal behavior, as will be discussed later. Below, the analysis of the tilt and shift global modes using this simple model is described. III. TILT STABILITY The tilt instability has been regarded as the most dangerous in FRC's, although it has not been observed consistently in the traditional theta-pinch formation scheme. Theoretically, it has been shown to be unstable in plasmas with large E due to the destabilizing j x B force, but it can be stabilized by many non-MHD effects. In this section, the simplest possible model is constructed to reveal the stabilizing physical mechanisms of this mode. The simplest model relating the j x B force to the decay index n_decay (= -(R / B_Z) d B_Z / d R) of the external field produced by the coils is the current ring model used in the study of spheromak tilt stability. However, this model does not provide a link between n_decay and plasma elongation E, which is important in FRC tilt stability and will be dealt with in this paper. Instead, we use the Solovev model to calculate directly the tilt stabilizing or destabilizing j x B force as a function of E. A. j x B torque In FRC's, the internal current flows in the theta [theta_Z in Fig. 2(a)] direction and is denoted by j_theta = (-j_theta sin theta, j_theta cos theta, 0). This j_theta interacts with the vacuum field B_V = [B_X, B_Y, B_Z] = [B_R(R,Z) cos theta, B_R(R,Z) sin theta, B_Z(R,Z)], resulting in a torque "density" n(r) = r x (j_theta x B_V), where r is the position vector (R cos theta, R sin theta, Z). Before tilting, however, the total torque N = int n dV integrated over the whole plasma volume is zero. After the plasma tilts a small angle, a responding torque N_1 arises either to accelerate the tilting or to restore the plasma to its original equilibrium, depending on the vacuum field provided for each plasma shape (elongation). Without losing generality, this responding torque N_1 = int n_1 dV is calculated as a first order perturbation in tilting angle theta_X with respect to the X axis. Instead of tilting the plasma an angle theta_X with respect to the stationary background vacuum field, one simple way to calculate the perturbed torque n_1 is to use the plasma frame of reference and to tilt the background vacuum field at an angle -theta_X. In the plasma frame, the vectors r and j_theta are unperturbed; therefore, n_1 can be evaluated simply from the perturbed vacuum field B_V1, i.e., n_1 = r x (j_theta x B_V1). To first order in theta_X, this n_1 is the same both in the plasma frame and the vacuum field frame, as calculated below when the vacuum field is tilted over an angle theta_X counter-clockwise with respect to the plasma frame. The vacuum field perturbation at r comes from two effects [see Fig. 2(b)]: (1) a direction change due to tilting, and (2) a magnitude change due to the fact that a different field originally located at r' = (X', Y', Z') = (R' cos theta, R' sin theta, Z') moves into the current location r as a result of tilting: B_V = [ B_R(R,Z) cos theta, B_R(R,Z) sin theta, B_Z(R,Z) ]^T -> [ B_R(R',Z') cos theta, B_R(R',Z') sin theta cos theta_X - B_Z(R',Z') sin theta_X, B_Z(R',Z') cos theta_X + B_R(R',Z') sin theta sin theta_X ]^T. Here R' and Z' are related to R and Z by R' approx R + theta_X Z sin theta and Z' approx Z - theta_X R sin theta. Therefore, B(R',Z') - B(R,Z) approx (dB/dR)(R' - R) + (dB/dZ)(Z' - Z) approx theta_X sin theta ( Z dB/dR - R dB/dZ ), where B = B_R or B_Z. Then the first order change in the vacuum field is given by B_V1 = theta_X [ sin theta cos theta ( Z dB_R/dR - R dB_R/dZ ), sin^2 theta ( Z dB_R/dR - R dB_R/dZ ) - B_Z, sin theta ( Z dB_Z/dR - R dB_Z/dZ ) + B_R sin theta ]^T. With the use of div B = (1/R) d(R B_R)/dR + dB_Z/dZ = 0, dB_R/dZ = dB_Z/dR, the perturbed torque can be calculated and simplified to n_1 = theta_X j_theta [ -R B_Z (1 - n_decay) sin^2 theta, R B_Z (1 - n_decay) sin theta cos theta, 0 ]^T, where n_decay is a generalized decay index including effects from the j x B force off the mid-plane (Z != 0), n_decay = - (R / B_Z) [ dB_Z/dR - (Z / R^2) d/dR(2 R B_R + Z B_Z) ]. Clearly, the total torque N_1 = int int int n_1 R dZ dR dtheta does not have Y nor Z components while the X component is given by N_1X = -pi theta_X int int j_theta R^2 B_Z (1 - n_decay) dZ dR. By using j_theta = B_0 R (4 + 1/E^2) / (mu_0 R_s^2) in the Solovev's model, N_1X can be written as N_1X = -pi theta_X R_s^3 (B_0^2 / mu_0) (4 + 1/E^2) chi_tilt, chi_tilt = int int (R^3 B_Z (1 - n_decay) / (R_s^5 B_0)) dZ dR, (1) where the nondimensional parameter chi_tilt can be calculated numerically as a function of E, as plotted in Fig. 3. We note that, in general, j_theta and B_Z have opposite signs and here j_theta > 0 and B_Z < 0 have been chosen. This chi_tilt can be fit to 0.02 E + 0.342 - 0.225/E + 0.0425/E^2 - 0.00329/E^3. When a small theta_Y is introduced, the responding torque N_1Y has the same expression as N_1X but with theta_X replaced by theta_Y. Figure 4 shows the normalized j x B torque as a function of E. It can be seen that if E <= 0.5, the j x B torque is negative. In other words, it restores the plasma toward its equilibrium position due to a strong mirror vacuum field. When E >= 0.5, the FRC is tilt unstable, consistent with previous MHD studies. We note that the force from the plasma pressure gradient should not contribute to the tilting torque since it is balanced by the unperturbed j x B_int during tilting, where B_int is the field produced by the internal current. FIG. 3. The dimensionless parameters chi_tilt and chi_shift as functions of plasma elongation E. The lines are fitting functions: chi_tilt = 0.02E + 0.342 - 0.225/E + 0.0425/E^2 - 0.00329/E^3 and chi_shift = -0.0132 - 0.168/E + 0.259/E^2 - 0.0917/E^3 + 0.0104/E^4. FIG. 4. The total j x B torque normalized by theta_X (B_0^2 / 2 mu_0) R_s^3 as a function of E. B. Stabilizing effect from plasma rotation It is well known that plasma rotation in the theta direction can help stabilize the tilt mode. In this subsection, the simplest possible equations using the rigid body model are used to study this effect. The three-axis rigid body rotation is governed by Euler's equations, I_X ddot{theta}_X - (I_Y - I_Z) dot{theta}_Y dot{theta}_Z = N_X, I_Y ddot{theta}_Y - (I_Z - I_X) dot{theta}_Z dot{theta}_X = N_Y, I_Z ddot{theta}_Z - (I_X - I_Y) dot{theta}_X dot{theta}_Y = N_Z, which can be simplified since I_X = I_Y = I, N_X / theta_X = N_Y / theta_Y = N, and dot{theta}_Z = const = Omega with no net driving torque in the Z direction, i.e., N_Z = 0. Then the reduced equations are I ddot{theta}_X - (I - I_Z) Omega dot{theta}_Y - N theta_X = 0, (2) I ddot{theta}_Y + (I - I_Z) Omega dot{theta}_X - N theta_Y = 0. (3) Taking the derivative of Eq. (3) and substituting dot{theta}_Y from Eq. (2), we have a fourth order differential equation for theta_X, I^2 theta_X^(4) + (b^2 - 2 I N) ddot{theta}_X + N^2 theta_X = 0, where b = (I - I_Z) Omega. Assuming that the solutions have the form theta_X = theta_X0 exp(-i omega t), a fourth order algebraic equation for omega is obtained, I^2 omega^4 - (b^2 - 2 I N) omega^2 + N^2 = 0, which yields a solution of omega^2 = (b^2 - 2 I N +- sqrt(b^4 - 4 I N b^2)) / (2 I^2). The necessary and sufficient condition for omega to be real is [(I - I_Z) Omega]^2 >= 4 I N, (4) which provides the minimum plasma rotation to stabilize the tilt, as plotted in Fig. 5. A stable region appears at large elongation. This can be understood due to the fact that rotational stabilization becomes more effective when the difference in the moments of inertia between the tilting axis and the Z axis becomes larger. In FRC plasmas, a large ion pressure gradient grad p_i exists due to a large plasma beta dominated by ion beta. An inherent plasma rotation arises from the ion diamagnetic drift since ions carry most of the plasma momentum. The magnitude of this naturally occurring rotation can be estimated as M_diamag = V_theta / v_thi = grad p_i / (e n B v_thi) approx (R_s / Delta_R) * (1 / s), (5) where Delta_R is the radial scale length of p_i, v_thi is the ion thermal velocity, and s is defined as the ratio of R_s to average ion gyro-radius. [This s is approximately equivalent to S* defined as R_s / (c / omega_pi) and roughly four times larger than bar{s} = int_R0^Rs R dR / (R_s rho_i).] The stability condition can be determined in the parameter space of s and E using this diamagnetic rotation, as shown in Fig. 6 where Delta_R = R_s - R_0 is assumed. A stable window appears at low s and large E (s/E <= 1.7), consistent with a previous study where an ion diamagnetic rotation was introduced to a two-fluid model. For a given E, the plasma rotates faster for smaller s due to ion diamagnetic drift, while for a given s the rotational stabilization becomes more effective with larger E as shown in Fig. 5. FIG. 5. The stability diagram in M (rotation Mach number) and elongation. FIG. 6. The stability diagram in s and elongation with a stabilizing effect from the ion diamagnetic drift. C. Stabilizing effect from ion gyro-viscosity The above analysis of FRC stability is based on a rigid body model, which is appropriate for external tilting. However, as pointed out in Sec. II B, this approach also models stability properties of internal tilting since the modeled region can be smaller than the whole plasma region, i.e., any inner part of the plasma. When only an inner part of the plasma tilts, the internal structure or profiles are deformed. If the typical ion gyro-radius is the same order as or larger than the spatial scale of this deformation, responding forces can arise from ion kinetic effects, in addition to the MHD force dealt with in Sec. II A. The stabilizing effects of ion kinetic motion have been considered by many authors in various schemes, but little physical insight into the underlying mechanisms has been given. Below, we approach this problem with minimum mathematical complications in an attempt to elucidate the physics but at the acknowledged expense of obtaining only semi-quantitative results. The so-called ion gyro-viscosity is one particular kinetic effect from ion gyro-motion. If there are spatial variations in force, such as an electric field force, ions tend to experience the variations over larger areas than electrons. If these variations are linear in space, ions experience a larger force during one half of their cyclotron period and a smaller force during the other half. As a result, the average force experienced can be approximated by the force at the guiding center, resulting in a null kinetic effect. However, when the spatial variation is more than linear, such as quadratic, the guiding center approximation fails due to incomplete cancellation between forces during one gyro-motion. In this case, in addition to the force at the guiding center, a correction proportional to the second derivative (or curvature) of the spatial variation is needed. When the force is perpendicular to the local magnetic field, the correction can be expressed in a form proportional to the curvature of the corresponding perpendicular flow, resulting in an effective viscosity often referred to as ion gyro-viscosity, although it arises without collisions. In FRC plasmas, the second radial derivative (or radial curvature) of the toroidal flow due to ion diamagnetic drift is not negligible especially in the case of a hollow current profile since ions carry a large portion of the plasma current perpendicular to the local field. Therefore, a correction force proportional to d^2 V_theta / d R^2 arises in the radial direction, pointing inward. (Confusion could arise here since the corresponding force for diamagnetic flow is the ion pressure gradient force which appears only in the fluid equations, in contrast with an electric field which can be felt by each gyrating particle. However, the complete Braginskii treatment does give rise to such viscosity terms proportional to pressure curvature regardless of the nature of the corresponding force.) An intuitive cartoon of the gyro-viscous force before tilting is shown in Fig. 7(a). They balance each other, resulting in a null tilting torque. However, when only an inner portion of the FRC plasma tilts, one part of the plasma is compressed while another is decompressed, resulting in changes in the gyro-viscous force, as indicated in Fig. 7(b). These perturbations form a restoring torque against tilt. This restoring torque can be divided into two parts: one from the sides of the FRC plasma and another from the ends, and they can be evaluated separately. A more quantitative expression for the viscosity tensor in the strong field limit with nonuniform V_theta has been given by Braginskii, pi_RR = - (2/3) (eta_0 / R) d V_theta / d theta - eta_3 d V_theta / d R; pi_R theta = - (2/3) (eta_3 / R) d V_theta / d theta; pi_RZ = - 2 eta_3 d V_theta / d Z, where eta_0 = 0.96 n T_i tau_i and eta_3 = n T_i / (2 Omega_i) (tau_i is ion collision time and Omega_i is the ion cyclotron frequency). Then the ion viscous force in the R direction is given by F_R = - [ d pi_RR / d R + d pi_R theta / (R d theta) + d pi_RZ / d Z ] = (2/3) (eta_0 / R) d^2 V_theta / (dR d theta) + eta_3 d^2 V_theta / d R^2 + (2/3) (eta_3 / R^2) d^2 V_theta / d theta^2 + 2 eta_3 d^2 V_theta / d Z^2, where the last term is small due to the nature of rigid body motion, i.e., V_theta changes only linearly in the Z direction. The first term is due to the effect of so-called magnetic pumping or parallel viscosity, but it contributes only a force parallel to the tilting axis. Thus no net torque exists. The second and third terms are forces due to the gyro-viscosity mentioned above. The relative strength of this force to the ion pressure gradient grad p_i is found to be (R_s^2 / (2 Delta_R^2)) * (1 / s^2). The increased importance for smaller s is consistent with physical intuition. Similar expressions for the gyro-viscous force acting at the ends of FRC plasmas can be found with Delta_Z denoting axial scales of V_theta near the ends. For a small tilting angle theta_X, Delta_R and Delta_Z change by delta_R = Z theta_X sin theta and delta_Z = R theta_X sin theta, respectively. Then the perturbed forces are F_1X = 3 delta_R eta_3 (V_theta / Delta_R^3) sin theta [ 1 + 2 Delta_R^2 / (9 R_s^2) ] = (3 theta_X n T_i / (2 s^2)) * (Z R_s^2 / Delta_R^4) sin^2 theta [ 1 + 2 Delta_R^2 / (9 R_s^2) ], F_1Z = 3 delta_Z eta_3 (V_theta / Delta_Z^3) sin theta [ 1 + 2 Delta_Z^2 / (9 R^2) ] = (3 theta_X n T_i / (2 s^2)) * (R R_s^2 / Delta_Z^4) sin^2 theta [ 1 + 2 Delta_Z^2 / (9 R^2) ], FIG. 7. Schematic views of the ion gyro-viscous forces exerted on the FRC plasma: (a) before tilt, and (b) after tilt. The size of arrows indicates the force strength. where Eq. (5) has been used. Assuming F_1X acts on the outer portion of the plasma with thickness Delta_R, and F_1Z acts on each end of the plasma with thickness Delta_Z, the restoring torque can be calculated as N_GV = 2 int_0^(Z0-Delta_Z) int_0^(2 pi) (F_1XZ) R_s dtheta dZ Delta_R [ 1 - Delta_R / (2 R_s) ] + 2 int_0^(Rs-Delta_R) int_0^(2 pi) (F_1ZR) R dtheta dR Delta_Z = (3 theta_X n T_i R_s^3 / s^2) int_0^(2 pi) sin^2 theta dtheta [ (1 - f_R/2) (1 + 2 f_R^2 / 9) int_0^(Z0-Delta_Z) (Z^2 / Delta_R^3) dZ + int_0^(Rs-Delta_R) (1 + 2 Delta_Z^2 / (9 R^2)) (R^3 / (R_s Delta_Z^3)) dR ] = (3 pi theta_X / 4) (n T_i R_s^3 / s^2) [ (sqrt(2) / 3) (E / f_R)^3 (1 - f_R/2) (1 + 2 f_R^2 / 9) (1 - sqrt(2) f_Z / E)^3 + (1 - f_R)^4 / f_Z^3 + (1 - f_R)^2 / (9 f_Z) ], (6) where f_R = Delta_R / R_s and f_Z = Delta_Z / R_s. If N_GV is larger than the j x B torque, N_JxB from Eq. (1), the FRC plasma is tilt stable. Figure 8 shows the stability diagram again in parameter space s and E, where f_R = (R_s - R_0)/R_s and f_Z = E f_R are used. A stable window appears at low s and large E (s/E <= 2.2), similar to the case of stabilization due to diamagnetic rotation. This trend is consistent with a previous analysis, which employed the more thorough energy principle but did not give a detailed physical picture. Also plotted in Fig. 8 are contributions from the restoring torque, from the end (dashed line) and from the side (dash-dotted line) of the FRC plasmas shown in Fig. 7. The increased stability at low s and large E is because V_theta becomes larger at small s [Eq. (5)] and because the restoring torque from the side becomes more effective (due to a larger arm) with larger E. FIG. 8. Stability diagram in s and elongation with stabilizing effects from ion gyro-viscosity. FIG. 9. Stability diagram in s and elongation with stabilizing effects from ion diamagnetic drift and ion gyro-viscosity. D. Tilt stability with plasma rotation and ion gyro-viscosity Now we can examine FRC tilt stability combining the effects of plasma rotation and ion gyro-viscosity. We define the following dimensionless parameters, each of which represents the contribution from a different effect: K_rot = (E beta_i / sqrt(2)) * (1 - 2 E^2 / 3)^2 / (1 + 2 E^2 / 3), K_GV = 3 beta_i f_R^2 [ (sqrt(2) / 3) (E / f_R)^3 (1 - f_R/2) (1 + 2 f_R^2 / 9) (1 - sqrt(2) f_Z / E)^3 + (1 - f_R)^4 / f_Z^3 + (1 - f_R)^2 / (9 f_Z) ], K_JxB = 8 E (4 + 1/E^2) chi_tilt, where beta_i = n T_i / (B_0^2 / 2 mu_0). Then the stability condition can be deduced from Eq. (4) with an N which includes contributions from both the j x B force [Eq. (1)] and the gyro-viscosity force [Eq. (6)], written as K_JxB (f_R s)^2 <= (1 + f_R s M)^2 K_rot + K_GV, or (K_JxB - M^2 K_rot) (f_R s)^2 - 2 M K_rot (f_R s) - K_rot - K_GV <= 0, where M is a rotational Mach number in addition to the rotation due to ion diamagnetic drift. If M = 0, the stability condition can be reduced to s^2 <= (K_rot + K_GV) / (f_R^2 K_JxB), which is plotted in Fig. 9, where the stability window expands to s/E <= 2.8. If M != 0, i.e., there is an additional rotation due to E x B drift; then generally the stability improves as shown in Fig. 10. A small additional rotation (such as M = 0.2) can stabilize tilt significantly at large E. IV. SHIFT STABILITY Similar approaches can be taken to study the other two types of rigid body motion of a FRC plasma: axial shift and radial shift. The calculations are much simpler since they are planar motions. Contributions from both j x B and ion gyro-viscosity are considered. When the plasma shifts in the axial (Z) direction by xi_Z, the perturbed j x B force in the Z direction is given by F_1Z = -xi_Z int j_theta (dB_Z / dR) dV. Then the equation of motion is M_total ddot{xi}_Z = - 2 pi xi_Z R_s (B_0^2 / mu_0) (4 + 1/E^2) chi_shift; chi_shift = int int (R^2 / (R_s^3 B_0)) (dB_Z / dR) dZ dR, where the nondimensional parameter chi_shift can be explicitly calculated as a function of E and fit to -0.0132 - 0.168/E + 0.259/E^2 - 0.0917/E^3 + 0.0104/E^4 (see Fig. 3). Figure 11(a) shows the normalized j x B force as a function of E. If E <= 1, the axial shift is stable due to a restoring j x B force. When E >= 1, the FRC is unstable to the axial shift mode. When only an inner part of the FRC plasma shifts axially, the ion gyro-viscosity will provide a restoring force in the Z direction as in the tilt stability, F_GV = - 3 pi xi_Z (n T_i R_s / s^2) * (1 - f_R)^2 / f_Z^3, where the perturbed force is assumed to be active only at both ends of the plasma with a thickness Delta_Z. Then the stability condition is obtained by setting the gyro-viscous force equal to the j x B force. Figure 11(b) shows a stability diagram in the parameter space of s and E. Unlike the case of tilt stability, gyro-viscosity has little effect on the axial shift stability except for in very small s regimes. A very similar analysis can be applied to a small radial shift xi_X in the X direction. The perturbed j x B force in the X direction is given by F_1X = int xi_X j_theta (dB_Z / dR) sin^2 theta dV = pi xi_X R_s (B_0^2 / mu_0) (4 + 1/E^2) chi_shift. As shown in Fig. 12(a), the radial shift stability window is reversed compared to axial shift. The stabilizing ion gyro-viscous force, F_GV = 3 pi xi_X (n T_i R_s / (sqrt(2) f_R^3 s^2)) E (1 - sqrt(2) f_Z / E) (1 + 2 f_R^2 / 9) (1 - f_R / 2), is found to have little effect in the radial shift case, as shown in Fig. 12(b). FIG. 10. Stability diagram with stabilizing effects from ion diamagnetic drift, ion gyro-viscosity, and additional rotation: (a) M < 0 and (b) M > 0. FIG. 11. (a) The J x B force normalized by xi_Z (B_0^2 / 2 mu_0) R_s as a function of E for a small axial shift. (b) The stability diagram for the axial shift mode in the parameter space of s and E with ion gyro-viscosity included. V. DISCUSSIONS AND CONCLUSIONS Despite the very simple models used in the present studies, much physical insight has been gained regarding global stability of FRC plasmas. Tilt stability is predicted, independent of s, for FRC's with low E (oblate), while tilt stability of FRC's with large E (prolate) depends on s/E. It is found that plasma rotation due to ion diamagnetic drift can stabilize the tilt mode when s/E <= 1.7. The so-called collisionless ion gyro-viscosity also is identified to stabilize tilt when s/E <= 2.2. Combining these two effects, the stability regime broadens to s/E <= 2.8, consistent with existing theories. A small additional rotation (such as a Mach number of 0.2) can improve tilt stability significantly at large E, although it is not self-consistently included in equilibrium as in Ref. 4. A similar approach is taken to study the physics of shift stability. It is found that radial shift is unstable when E < 1 while axial shift is unstable when E > 1. However, unlike the tilt stability, gyro-viscosity has little effect on the shift stability. In an attempt to compare with experiments, Figure 13 plots the s-E tilt stability diagram together with existing experimental observations in theta-pinch devices. Without E x B rotation, some experimental data are outside the stable regime. But with a certain level of E x B rotation, all the data can be included in the stable regime, although a consistent E x B rotation has not been experimentally established in these devices. Of course, one should bear in mind that the above results are only semi-quantitative. Also likely is that the stability windows shrink when non-rigid-body motions are taken into account. A recent study shows that parallel viscosity due to both collisions or collisionless pitch-angle scattering can smooth out variations along the field line thus making the motion more rigid-body-like. On the other hand, the omitted effects of plasma compressibility (including both thermal and magnetic) and sheared flows are likely to broaden the stability window. Therefore, a more precise stability diagram would require a substantial theoretical and numerical effort, which is beyond the scope of the present work. Nonetheless, the present work has elucidated two important stabilization mechanisms in detail for the FRC tilt stability: plasma rotation (due to both ion diamagnetic drift and E x B drift) and ion gyro-viscosity. The tilt stability observed in past FRC plasmas formed by theta-pinches can be well explained by these two effects due to their large E. Another important conclusion is that FRC plasmas are always likely to be unstable to at least one global instability with any combination of s and E. Passive or active stabilizers (coils or conducting shells) are required to stabilize these global modes completely. Theta-pinch FRC's are subject to the axial shift instability, but it seems to be stabilized by mirror coils. FRC plasmas made by merging spheromaks are likely to be subject to the radial shift instability due to their oblate shape. A conductive shell is probably needed to stabilize this mode. The global stability of FRC plasmas with both oblate and prolate shapes and various s can be explored further in the proposed project SPIRIT (Self-organized Plasma with Induction, Reconnection, and Injection Techniques) based on the counter-helicity spheromak merging. FIG. 12. (a) The J x B force normalized by xi_X (B_0^2 / 2 mu_0) R_s as a function of E for a small radial shift. (b) The stability diagram for radial shift in the parameter space of s and E with ion gyro-viscosity included. FIG. 13. A comparison of theoretical predictions of tilt stability diagram with experimental observations of stable FRC plasmas. ACKNOWLEDGMENTS The authors thank Dr. L. Steinhauer and Dr. S. Jardin for their valuable suggestions and comments. APPENDIX: NUMERICAL CALCULATION OF EXTERNAL FIELDS FOR SOLOVEV EQUILIBRIA A free-boundary solution for the Solovev FRC equilibrium provides external fields for calculation of the tilt and shift instabilities. To obtain the solution we place N_C poloidal field coils on a closed contour surrounding the plasma and seek a set of coil currents which match appropriate boundary conditions at the plasma separatrix. The total flux, psi, anywhere in space can be written as separate contributions from the plasma and coil current sources, psi = psi_p + psi_c. Axisymmetric Green's functions relate the current sources to the fluxes: psi_p(X_i, Z_i) = int int J_phi(X,Z) G(X,Z; X_i, Z_i) dX dZ, where J_phi = (1/X) Delta* psi = - X p'. Here, p' = p_0 / psi_0 is the derivative of the pressure with respect to the poloidal flux, and is constant for the Solovev model. Similarly, the "external flux" provided by the poloidal field coils is given by psi_c(X_i, Z_i) = sum_{j=1}^{N_C} G(X_j, Z_j; X_i, Z_i) I_j. Specifically, the X_i, Z_i are chosen to be N_B equally spaced points on the plasma separatrix, excluding X = 0, and the N_C coils are equally spaced on (X/X_0)^2 + (Z/Z_0)^2 = 2 a_c^2, a contour conformal with the separatrix. Since psi = 0 on the separatrix, matching fluxes across the plasma-vacuum interface gives sum_{j=1}^{N_C} G_ij I_j = R_i, (A1) where R_i = - psi_p(X_i, Z_i), G_ij = G(X_i, Z_i; X_j, Z_j), i = 1, 2, ..., N_B. The coil currents, I_j, are found by solving Eq. (A1) by the method of Least Squares: Min [ W = sum_{i=1}^{N_B} sigma_i ( sum_{j=1}^{N_C} G_ij I_j - R_i )^2 + alpha_reg sum_{j=1}^{N_C - 1} (I_{j+1} - I_j)^2 ]. The regularization term multiplying alpha_reg avoids coil-to-coil current oscillation. Typical numerical parameters are N_C = 20, N_B = 20, a_c = 2.0, alpha_reg = 0.01. With suitably chosen weights, sigma_i, maximum errors Max_i |(sum G_ij I_j) / R_i - 1| < 0.005 are typically obtained (see Figure 1). REFERENCES 1 M. Tuszewski, Nucl. Fusion 28, 2033 (1988). 2 M. Tuszewski, D. C. Barnes, R. E. Chrien, J. W. Cobb, D. J. Rej, R. E. Siemon, D. P. Taggart, and B. L. Wright, Phys. Rev. Lett. 66, 711 (1991). 3 A. Mohri, J. Appl. Phys. 19, L686 (1980). 4 R. A. Clemente and J. L. Milovich, Phys. Fluids 26, 1874 (1982). 5 A. Ishida, H. Momota, and L. C. Steinhauer, Phys. Fluids 31, 3024 (1988). 6 D. C. Barnes, J. L. Schwarzmeier, H. R. Lewis, and C. E. Seyler, Phys. Fluids 29, 2616 (1986). 7 L. C. Steinhauer and A. Ishida, Phys. Fluids B 2, 2422 (1990). 8 R. Horiuchi and T. Sato, Phys. Fluids B 2, 2652 (1990). 9 A. Ishida, R. Kanno, and L. C. Steinhauer, Phys. Fluids B 4, 1280 (1992). 10 K. Nishimura, R. Horiuchi, and T. Sato, Phys. Plasmas 4, 4035 (1997). 11 Y. Nomura, J. Phys. Soc. Jpn. 54, 1369 (1985). 12 D. C. Barnes and R. D. Milroy, Phys. Fluids B 3, 2609 (1991). 13 J. W. Cobb, T. Tajima, and D. C. Barnes, Phys. Fluids B 5, 3227 (1993). 14 L. C. Steinhauer, A. Ishida, and R. Kanno, Phys. Plasmas 1, 1523 (1994). 15 R. Kanno, A. Ishida, and L. C. Steinhauer, J. Phys. Soc. Jpn. 64, 463 (1995). 16 A. Solovev, in Review of Plasma Physics (Consultants Bureau, New York, 1966), Vol. 6, p. 257. 17 M. J. Hill, Philos. Trans. R. Soc. London, Ser. A C/XXXV, 213 (1894); M. Y. Wang and G. H. Miley, Nucl. Fusion 19, 39 (1979). 18 C. Munson, A. Janos, F. Wysocki, and M. Yamada, Phys. Fluids 28, 1525 (1985). 19 J. L. Schwarzmeier, D. C. Barnes, D. W. Hewett, C. E. Seyler, A. I. Shestakov, and R. L. Spencer, Phys. Fluids 26, 1295 (1983). 20 See, e.g., Classical Dynamics, by J. B. Marion and S. T. Thornton, 3rd ed. (Saunders College Publishing, 1988), p. 386. 21 M. N. Rosenbluth, N. A. Krall, and N. Rostoker, Nucl. Fusion Suppl. 1, 143 (1962), and references therein. 22 S. I. Braginskii, in Review of Plasma Physics (Consultants Bureau, New York, 1966), Vol. 1, p. 218. 23 A. L. Hoffman, J. T. Slough, L. C. Steinhauer, N. A. Krall, and S. Hamasaki, in Plasma Physics and Controlled Nuclear Fusion Research, Proceedings of the 11th International Conference (International Atomic Energy Agency, Vienna, 1987), Vol. II, p. 541. 24 R. E. Siemon, W. T. Armstrong, D. C. Barnes, R. R. Bartsch, R. E. Chrien, J. C. Cochrane, W. N. Hugrass, R. W. Kenish, Jr., P. L. Klingner, H. R. Lewis, R. K. Linford, K. F. McKenna, R. D. Milroy, D. J. Rej, J. L. Schwarzmeier, C. E. Seyler, E. G. Sherwood, R. L. Spencer, and M. Tuszewski, Fusion Technol. 9, 13 (1986). 25 S. Okada, Y. Kiso, S. Goto, and T. Ishimura, Phys. Fluids B 1, 2422 (1989). 26 J. T. Slough, E. A. Crawford, A. L. Hoffman, R. D. Milroy, R. Maqueda, and G. A. Wurden, in Plasma Physics and Controlled Nuclear Fusion Research, Proceedings of the 14th International Conference (International Atomic Energy Agency, Vienna, 1993), Vol. II, p. 627. 27 N. Iwasawa, A. Ishida, and L. C. Steinhauer, see AIP Document No. E-PAPS:E-PHPAEN-10-031810 for 6 pages of text of "Global mode stability of field-reversed configurations." Order by PAPS number and journal reference from American Institute of Physics, Physics Auxiliary Publication Service, Carolyn Gehlbach, 500 Sunnyside Boulevard, Woodbury, NY 11797-2999. 28 L. C. Steinhauer and A. Ishida, Phys. Rev. Lett. 79, 3423 (1997). 29 Y. Ono, A. Morita, T. Itagaki, and M. Katsurai, in Ref. 26, p. 619. 30 M. Yamada, H. Ji, and P. Heitzenroeder, Bull. Am. Phys. Soc. 34, 2071 (1997).