Princeton Plasma Physics Laboratory PPPL-4075 Numerical Study of Field-reversed Configurations: The Formation and Ion Spin-up E.V. Belova, R.C. Davidson, H. Ji, M. Yamada, C.D. Cothran, M.R. Brown, and M.J. Schaffer June 2005 Prepared for the U.S. Department of Energy under Contract DE-AC02-76CH03073. PPPL Report Disclaimers Full Legal Disclaimer This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, nor any of their contractors, subcontractors or their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or any third party's use or the results of such use of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof or its contractors or subcontractors. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. Trademark Disclaimer Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof or its contractors or subcontractors. PPPL Report Availability This report is posted on the U.S. Department of Energy's Princeton Plasma Physics Laboratory Publications and Reports web site in Fiscal Year 2005. The home page for PPPL Reports and Publications is: http://www.pppl.gov/pub_report/ Office of Scientific and Technical Information (OSTI): Available electronically at: http://www.osti.gov/bridge. Available for a processing fee to U.S. Department of Energy and its contractors, in paper from: U.S. Department of Energy Office of Scientific and Technical Information P.O. Box 62 Oak Ridge, TN 37831-0062 Telephone: (865) 576-8401 Fax: (865) 576-5728 E-mail: [email protected] National Technical Information Service (NTIS): This report is available for sale to the general public from: U.S. Department of Commerce National Technical Information Service 5285 Port Royal Road Springfield, VA 22161 Telephone: (800) 553-6847 Fax: (703) 605-6900 Email: [email protected] Online ordering: http://www.ntis.gov/ordering.htm Numerical Study of Field-Reversed Configurations: The Formation and Ion Spin-up E. V. Belova 1), R. C. Davidson 1), H. Ji 1), M. Yamada 1), C. D. Cothran 2), M. R. Brown 2), M. J. Schaffer 3) 1) Princeton Plasma Physics Laboratory, Princeton NJ, USA, 2) Swarthmore College, Swarthmore PA, USA, 3) General Atomics, San Diego CA, USA E-mail: [email protected] Short title: Numerical Study of FRCs PACS numbers: 52.55.Lf, 52.65.Kj, 52.65.Rr, 52.35.Vd Abstract Results of three-dimensional numerical simulations of field-reversed configurations (FRCs) are presented. Emphasis of this work is on the nonlinear evolution of magnetohydrodynamic (MHD) instabilities in kinetic FRCs, and the new FRC formation method by counter-helicity spheromak merging. Kinetic simulations show nonlinear saturation of the n = 1 tilt mode, where n is the toroidal mode number. The n = 2 and n = 3 rotational modes are observed to grow during the nonlinear phase of the tilt instability due to the ion spin-up in the toroidal direction. The ion toroidal spin-up is shown to be related to the resistive decay of the internal flux, and the resulting loss of particle confinement. Three-dimensional MHD simulations of counter-helicity spheromak merging and FRC formation show good qualitative agreement with results from the SSX-FRC experiment. The simulations show formation of an FRC in about 20 - 30 Alfvén times for typical experimental parameters. The growth rate of the n = 1 tilt mode is shown to be significantly reduced compared to the MHD growth rate due to the large plasma viscosity and field-line-tying effects. 1. Introduction The field-reversed configuration (FRC) is a compact toroid with little or no toroidal field. It offers a unique fusion reactor potential because of its compact and simple geometry, translation properties, and high plasma beta [1]. At present, the most important issues are FRC stability with respect to low-n (toroidal mode number) MHD modes, and the development of new FRC formation and current-drive methods. The traditional theta-pinch formation method usually produces highly kinetic FRCs with relatively low flux, small S* (the FRC kinetic parameter S* is the ratio of the separatrix radius to the ion skin depth) and large elongation E. A theoretical understanding of the observed FRC stability properties has proven to be elusive, due to the complicated interplay of several non-ideal MHD effects [2-4], including finite-Larmor-radius (FLR) effects, the Hall term, and plasma flow effects. Advanced numerical simulations are required to describe the self-consistent stability properties of kinetic FRCs [3,5,6]. Results of a set of such simulations are presented in this paper. A slow FRC formation technique, based on counter-helicity spheromak merging has demonstrated the advantage of this approach compared with traditional theta-pinch formation methods [7]. The counter-helicity spheromak merging method allows formation of the configuration with large S*, thus permitting experimental studies of large-S* FRC stability properties. The SSX-FRC experiment [8] is designed to study FRC formation by the counter-helicity spheromak merging method, and to examine the general issue of FRC stability properties at large S*. Three-dimensional MHD simulations have been performed in support of the SSX-FRC experiment, and show good agreement with the experimental results. 2. Study of FRC nonlinear stability properties Numerical studies of the nonlinear evolution of magnetohydrodynamic (MHD) instabilities in kinetic, prolate FRCs (theta-pinch-formed FRCs) have been performed using the 3D nonlinear hybrid and MHD simulation code HYM [5]. The stability properties of MHD modes with toroidal mode numbers n ≥ 1 have been investigated, including finite-ion-Larmor radius (FLR) and rotational effects. It has been demonstrated that due to the strong FLR stabilization of the higher-n modes, the n = 1 tilt mode is the most unstable mode for nearly all experimentally-relevant non-rotating FRC equilibria [3]. An empirical FLR scaling of the tilt mode linear growth rate has been obtained: γ = C V_A / (R_s E) exp(-3 E ρ_i / R_s), where γ_mhd = C V_A / (R_s E) is the MHD growth rate, R_s is the separatrix radius, E is the separatrix elongation, ρ_i is ion thermal Larmor radius, and C ≈ 2 is a constant. Nonlinear kinetic simulations performed for a set of FRC equilibria with E = 4 - 6 and S* = 10 - 80 show that the n = 1 tilt mode saturates nonlinearly without destroying the configuration, provided the FRC kinetic parameter is sufficiently small, S* ≲ 20. In addition to the saturation of the tilt mode, nonlinear hybrid simulations show that the ions spin-up toroidally in the ion diamagnetic direction, and the n = 2 mode grows in the nonlinear phase of the simulation. Initial conditions for the simulations have been set at t = 0 so that all of the equilibrium toroidal current is carried by the electrons, and the ions have a non-rotating Maxwellian distribution. These initial conditions are consistent with experimental observations just after FRC formation. However, as the simulation proceeds, the ions gradually begin to rotate, and near the end of the simulation run the ion toroidal flow velocity becomes comparable to the ion diamagnetic velocity. Therefore, the saturation of the tilt instability occurs in the presence of a significant ion toroidal rotation, and it is accompanied by the growth of the n = 2 rotational mode, which is often seen in experiments [1]. The saturation of the n = 1 tilt mode and the growth of the n = 2 rotational mode can be seen in Fig. 4 of Ref. [3], where the results of nonlinear hybrid simulations are shown for a configuration with E = 6.25 and S* ≈ 20 (for comparison, the parameters for typical theta-pinch FRC experiments [1] are E = 5 - 8 and S* ≲ 20). In the nonlinear phase, the ion rotation rate is comparable to the linear growth rate of the tilt mode. Therefore, the ion toroidal spin-up, in addition to driving the rotational instability, is likely to contribute significantly to the saturation of the tilt instability. 2.1. Ion spin-up A set of 2D (axisymmetric) nonlinear hybrid simulations has been performed in order to study the resistive evolution of the kinetic FRC, and investigate the mechanism of the ion toroidal spin-up. The initial ion distribution function is taken to be f = f(ε) = A exp(-ε/T_0), where ε = m_i v^2 / 2 + eφ is the ion energy, T_0 is the (uniform) ion temperature, and φ is the electrostatic potential. Since the ion pressure is scalar, a kinetic equilibrium for the hybrid runs has been constructed similar to the ideal MHD equilibrium, i.e., by solving the Grad-Shafranov equation for a chosen pressure profile. The cold-fluid description is used for the electrons, and quasi-neutrality is assumed. The numerical simulations have been performed for an FRC with E = 4, S* = 20, and an elliptical separatrix shape. The value of the normalized resistivity at the field null is η_o = 1/S = 10^-4, where S = V_A R_c / η is the Lundquist number, and a resistivity profile has been used with η inversely proportional to the plasma density. The simulations show that there is a significant particle loss associated with the resistive decay of the poloidal flux. Figure 1 shows the time evolution of the trapped poloidal flux ψ_0(t), and the number of ions inside the separatrix region ψ < 0. It can be seen that the decay of the poloidal flux results in the loss of a significant fraction of the particles at t ≳ 50 t_A. Here, the Alfvén time is defined as t_A = R_c / V_A, where R_c is the radius of the flux conserving shell, and V_A is the characteristic Alfvén velocity. Analysis of the particle phase-space shows that the particles with initial values of p_φ such that 0 < p_φ ≲ Δψ are lost from the closed-field-line region when the absolute value of the trapped poloidal flux ψ_0 is reduced by Δψ. Here p_φ = m_i/e R v_φ - ψ is the canonical toroidal angular momentum of the ion, and the signs are chosen such that p_φ > 0 is the approximate confinement condition. Since most of the particles with small p_φ have negative toroidal velocity, the particle loss results in a net flux of the negative momentum away from the separatrix region, and therefore there is a net positive ion rotation inside the separatrix (where the positive direction corresponds to the current direction). Therefore, the ion toroidal spin-up in the simulations is related to the resistive decay of the internal flux, and the resulting loss of weakly-confined particles. Figures 2a shows the time evolution of the toroidal angular momentum of all ions, L = ∫ R ∫ v_φ f d^3v d^3x, and the angular momentum of the ions inside the separatrix region obtained in the simulations shown in Fig. 1. It can be seen that the net angular momentum is approximately conserved due to the imposed periodic boundary conditions in z-direction. In contrast, the angular momentum of the part of the plasma confined inside the separatrix is positive, and increases in time as the configuration decays. (Due to the imposed periodic boundary condition, particles which leave the closed-field-line region remain on the open field lines, so that the plasma outside the separatrix spins-up in the negative direction.) The maximum value of the ion flow velocity is plotted in Fig. 2b. The peak value of V_φ is about 0.2-0.3 V_A, which is comparable to the ion diamagnetic velocity for S* ≈ 20. In the final state, the ions carry a significant fraction of the total current. Similar values of the ion toroidal flow velocity have also been obtained in nonlinear three-dimensional (3D) simulations [3]. Both 2D and 3D simulations with zero initial ion rotation demonstrate the formation of an approximate rigid-rotor profile inside the separatrix in about 40-60 Alfvén times, depending on the plasma resistivity. Poloidal contour plots and radial profiles of the ion toroidal flow velocity are shown in Fig. 3. Figure 3 shows that the ion flow velocity changes sign outside the separatrix, and that there is a significant velocity gradient close to the separatrix at the plasma edge. The n = 2 and n = 3 rotational instabilities are found to be localized in the vicinity of maximum velocity shear near the edge, and have a similar structure to the external modes. The growth rates of the rotational modes are found to be larger in smaller-S*, more kinetic configurations. The details of the ion toroidal spin-up determine the nonlinear evolution of these instabilities. Toroidal ion spin-up has always been observed during the quasi-steady-state decay phase of FRC experiments [1]. The measured ion rotation frequency is comparable to the ion diamagnetic frequency α = Ω_i / Ω_di ~ 1 by the time when about one-half of the plasma has been lost due to decay, and the n = 2 rotational instability begins to grow. In addition, recent observations in the TS-3 and TS-4 experiments [9] suggest that the ion toroidal spin-up plays an important role during the observed nonlinear stabilization of the n = 1 tilt mode. These observations are in a good agreement with the results of 2D and 3D hybrid simulations described above and in Ref. [3]. A significant number of theoretical and experimental studies have been performed to investigate the plasma spin-up in FRCs [1]. Two major physical mechanisms which have been proposed to explain the observed ion rotation are: particle loss [10], and the end-shortening of the radial electric field [1,11]. The simulation results presented here provide good agreement with the observed ion rotation, even though the end-shortening mechanism is not included in our simulations due to the imposed periodic boundary conditions. The experimental estimates of the particle loss time and the spin-up time are also consistent with the assumption of spin-up resulting from the particle loss, at least for smaller FRC devices [12]. The conclusion of Reference [11] that the particle loss will result in spin-up in the wrong direction when electric-field effects are taken into account is incorrect, because (1) it neglects the magnetization current of the lost ions; and (2) it considers rotation of the lost ions relative to the laboratory frame, rather than their rotation relative to the bulk ion population. The linear velocity profile (V_φ ~ R) obtained in 2D simulations (Fig. 3) suggests that the distribution function f_i of the ions inside the separatrix evolves towards an exponential rigid-rotor distribution function, which corresponds to a shifted local Maxwellian distribution. The evolution of f_i towards a Maxwellian distribution in a collisionless model and in the absence of 3D instabilities may indicate that the stochasticity of the ion orbits plays a significant role in FRC relaxation. Earlier one-dimensional (neglecting axial variation) hybrid simulations of the ion toroidal spin-up have modeled the particle loss by removing the simulation particles whose gyrocenter was close to the separatrix [13]. The resulting plasma rotation was localized near the separatrix, in contrast to the nearly-rigid-rotation profiles found in the 2D simulations described above (Fig. 3). This difference in the rotation profiles can be explained by the complexity of the ion orbits in the two-dimensional poloidal field, and (perhaps) the model chosen for the particle loss in Ref. [13]. 3. Counter-helicity spheromak merging simulations FRC formation by the counter-helicity spheromak merging method has been developed in the TS-3 experiments in Japan [7,9]. These experiments have shown that an FRC is formed after the opposing toroidal magnetic fields of two merging spheromaks are annihilated, and the plasma is heated to form an FRC-like pressure profile. This method allows the slow formation of an FRC with large poloidal flux, as well as studies of magnetic reconnection and MHD relaxation in high-beta plasmas [14,15]. The SSX-FRC experiment is designed to study the new FRC formation method by counter-helicity spheromak merging, and the FRC stability properties for large values of S*. In addition, the effects of the residual (axially-antisymmetric) toroidal field on macroscopic stability properties are being studied. Experimental results demonstrate the formation of an FRC-like configuration with S* > 35, and indicate the presence of the global n = 1 instability, consistent with the tilt-mode instability [8]. However, the observed growth rate of this instability is smaller by a factor of 6-8 than that of the ideal MHD growth rate. In addition, the experimental formation studies consistently show the presence of an axially-antisymmetric toroidal field, which does not completely annihilate during the reconnection. The underlying reasons for this are not yet understood. Numerical simulations using the HYM code have been performed to investigate this issues, and to study 3D spheromak merging for experimentally-relevant parameters. Two- and three-dimensional MHD simulations of counter-helicity spheromak merging have been performed for the SSX-FRC geometry and parameter range. The simulations have been performed using the 3D nonlinear resistive MHD version of the HYM code [5] with high resolution up to 127 × 513 grid points in the poloidal (R, Z) plane. The boundary conditions are taken to correspond to a cylindrical flux conserver which is perfectly conducting on the perturbation (fast) time-scale, but allowing for equilibrium poloidal magnetic field penetration, thus taking into account the field-line-tying effects present in the experiment. The initial spheromak formation by plasma guns have not been simulated. Instead, the initial conditions for the numerical studies have been chosen to correspond to the experimental conditions at the beginning of the spheromak merging process. Thus the initial conditions in the simulations correspond to two, low-beta, nearly force-free spheromaks with opposite toroidal magnetic fields, placed close to the device midplane. The SSX-FRC experimental parameters are as follows: flux conserver radius and length are R_c = 20 cm and L = 60 cm, respectively, the edge magnetic field is B_0 = 1 kG, the plasma density is 10^15 cm^-3, the plasma temperature (after the FRC formation) is T ≈ 30 eV, and the FRC kinetic parameter is S* = R_c/λ_i = 28. The Alfvén time is defined as t_A = R_c / V_A, which can be estimated as t_A = 2.8 µs. It has been observed that the poloidal flux at the midplane rises to its maximum value in t ~ 10 - 15 t_A (after the spheromaks are ejected at t = 20 µs), and that a high-beta FRC-like configuration is formed at the midplane at t ≳ 20 - 25 t_A. Growth of the n = 1 instability is seen at t ≳ 18 t_A, when the perturbation amplitude becomes larger than the level of the n = 1 turbulence, i.e., at about 10% of poloidal magnetic energy density. The characteristic growth time of the n = 1 mode is estimated to be 10 - 14 t_A, and the n = 1 component of the magnetic field energy becomes the dominant component at t ≳ 30 t_A. The experimental parameters and results are described in greater detail elsewhere [8]. 3.1. Axisymmetric simulations A set of axisymmetric simulations of counter-helicity spheromak merging have been performed in order to study the dependence of the reconnection rate and the toroidal field annihilation on values of the plasma resistivity and viscosity. For simplicity, the resistivity and viscosity profiles have been assumed to be uniform in the simulations. The Lundquist number S = V_A R_c / η based on classical resistivity can be estimated as S ≈ 10^3 for T_e ~ T_i ≈ 12 eV. For these temperatures, the ions are collisional with ω_ci τ_i ~ 1, and there are several estimates for the normalized plasma viscosity coefficient ν = 1/Re, where Re is the Reynolds number. The estimates of Braginskii's unmagnetized and weakly-magnetized ion viscosity are ν_0 ≈ 8 × 10^-3 and ν_1 ≈ 2 × 10^-3, respectively, and the gyroviscosity is ν_gyro ≈ 5 × 10^-3. Due to the strong dependence on the ion temperature, there is a large uncertainty in the above values. Nevertheless, it can be seen that the SSX-FRC plasma is in the viscosity-dominated regime with ν > η. Numerical simulations for η = 10^-3 and µ = 10^-3 - 4 × 10^-3 show formation of an FRC in about 20-30 Alfvén times, where µ = n̂ ν is the dynamic viscosity, and n̂ is the normalized plasma density. Evidently, large toroidal and poloidal flows with flow velocity up to ~ 0.5 - 1 V_A (based on the edge field) are generated during the reconnection phase (Fig. 4b). The plasma pressure is significantly increased by Ohmic and viscous heating, and the FRC-like pressure profile is formed due to convective transport by poloidal flows (Fig. 4a). The flow velocity reduces by an order-of-magnitude after the FRC formation is complete. These results agree with previous numerical studies of axisymmetric spheromak merging [16,17]. Figure 5 shows the results of five MHD simulation runs performed for values of η = 2.5 × 10^-4 - 4 × 10^-3, and µ = 10^-3. It can be seen that the time for a complete reconnection depends weakly on the resistivity (this time can be approximately determined by the strong reduction in the kinetic energy in Fig. 5a). Thus the reconnection is complete and the FRC forms in about 20 t_A for larger values of η, and in about 25-30 t_A for smaller η. This is consistent with previous theoretical studies of driven magnetic reconnection [16,18]. On the other hand, the details of the time evolution of the plasma kinetic energy (Fig. 5a) and magnetic energy (Fig. 5b) vary strongly with η. For larger values of η, the kinetic energy peaks at t ~ 5 t_A, and the flow velocity (mostly toroidal) increases with η. This initial increase in the kinetic energy is caused by the force imbalance present in the initial conditions. For values η ≥ 0.004, the configuration decays faster than the FRC forms. In the smaller resistivity cases (η ≤ 0.002), there is a maximum in the flow energy which occurs at t ~ 20 - 25 t_A. This peak in the kinetic energy correlates with the large reconnection rate and the fast reduction in the magnetic energy seen in Fig. 5b. The toroidal flow velocity in these cases is comparable to the characteristic Alfvén velocity, V_φ ~ V_A. Figure 6 shows the results of a set of MHD simulation runs performed to investigate effects of plasma viscosity for η = 5 × 10^-4. As expected, larger viscosity results in a strong reduction of the plasma flow (Fig. 6a). The reduced flow velocity, on the other hand, is related to a slower reconnection of the toroidal field (Fig. 6b). Figure 7 shows contour plots of the toroidal field from the simulations with η = 0.001, and different values of viscosity µ = 0.001, and µ = 0.004. It can be seen that larger values of plasma viscosity result in incomplete reconnection and significant residual toroidal fields, which are present at t ~ 20 - 30 t_A. Figure 7 also shows the reversal of the initial toroidal field on the outer flux surfaces at t ≳ 20 t_A associated with the so-called "sling-shot effect" [14,15]. For the case with η = µ = 0.001 (Fig. 7a), the magnitude of the reversed field is about 10% of the maximum toroidal field at t = 25 t_A. This effect is weak in the large viscosity case (Fig. 7b) due to the reduction in the toroidal flow velocity. The width of the reconnection layer estimated from the axial profile of J_R is δ ≈ 0.8 cm for the low-viscosity case with µ = 0.001, but it increases to δ ≈ 1.6 cm for µ = 0.004. For comparison, the experimentally measured reconnection layer is about 2 cm to 3 cm wide [19]. Therefore, the numerical results indicate that the relatively large values of plasma viscosity in the experiments (µ ≳ 0.002) may be responsible for the incomplete reconnection of the toroidal magnetic field which is observed in SSX-FRC. The simulations predict a reduced "sling-shot" effect in the SSX-FRC experiment, and a reconnection layer width of about 2 cm, which is comparable with the experimentally measured one. 3.2. Three-dimensional simulations Three-dimensional MHD simulations have been performed to study the effects of plasma viscosity, self-generated flows and magnetic field-line-tying effects on the unstable global modes. A random initial perturbation is applied at t = 0. The n = 1 tilt mode is found to be a dominant mode in all cases, and the configuration is tilted at the end of the simulations. Figure 8a shows the time evolution of the n = 1 mode obtained in four simulation runs with η = 10^-3 and µ = 5 × 10^-4 - 4 × 10^-3. In the linear phase, t ≲ 20 t_A, the growth rate is largest for small viscosity, and the growth rate reduces as µ increases. The simulations demonstrate that large values of viscosity have a stabilizing effect on the n = 1 tilt mode, as well as the higher-n MHD modes. The nonlinear slow-down of the tilt instability, seen at t ≈ 20 - 25 t_A in low-viscosity runs, occurs when the amplitude of the n = 1 mode becomes comparable to that of the n = 0 mode. The slowing-down of the tilt motion in other cases (µ = 0.002 and µ = 0.004) observed at t > 20 t_A is due to the evolution of the background quasi-equilibrium, and the resistive decay of the magnetic field, which reduces the instability drive. Figure 8b shows the time evolution of the total kinetic energy obtained for the same set of simulations as shown in Fig. 8a. The growth of the n = 1 tilt mode is not seen until t > 20 t_A, when the amplitude of this mode becomes comparable to the n = 0 component of the kinetic energy. Note that by that time the growth rate of the n = 1 mode is reduced compared to its linear value. The calculated growth rates obtained from the simulation shown in Fig. 8a for µ = 0.004 are: γ ≈ 0.33 · V_A / R_s for t < 20 t_A, and γ ≲ 0.15 · V_A / R_s for t > 20 t_A. For comparison, an estimate for the ideal MHD growth rate is γ_0 ~ (0.7 - 1.3) V_A / R_s. Several factors may contribute to the reduction of the instability growth at later times (t ≳ 20 t_A), including nonlinear mode interactions and an increase of the separatrix elongation. Additional time-evolution plots are shown in Fig. 9 where the peak values of the poloidal and toroidal fields, the radial current density, the axial flow velocity, and the toroidal flow velocity are shown for 3D simulations with η = 0.001 and (a) µ = 0.001 and (b) µ = 0.004. The growth of the tilt instability can be seen in these plots at t > 25 t_A for µ = 0.001, and at t > 35 t_A for the µ = 0.004 simulation. The toroidal magnetic field reconnection rate is proportional to the radial current J_R, which is significantly smaller in case (b) for t ≲ 22 t_A. Figure 9 shows that the maximum toroidal flow velocity is reduced approximately by a factor of two when µ is increased from 0.001 to 0.004. The possible stabilizing effect of the toroidal flows, generated during the reconnection process, appears to be less significant than the stabilizing effect of viscosity, because the growth rates for the larger V_φ cases (i.e., smaller µ) are larger than the growth rates for smaller V_φ (larger µ). The finite residual toroidal field, on the other hand, may contribute to the reduction of the growth rate of the n = 1 mode in the high-viscosity cases (µ > 10^-3). The simulations shown in Figs. 8 and 9 have been performed for realistic boundary conditions, including the effects of magnetic field line-tying. Another set of simulations have been performed with different boundary conditions, i.e., by neglecting these effects. A comparison between the runs with different boundary conditions show that line-tying effects increase the reconnection time by 5 - 10 t_A, and reduce the n = 1 mode growth rate. It has also been found that the peak toroidal flow velocity is reduced by a factor of two due to line-tying effects. In summary, the simulations show that both the large plasma viscosity and the field-line-tying boundary conditions reduce the growth rate of the tilt mode. In addition, there is a nonlinear reduction of the instability drive at the time when the mode amplitude becomes experimentally observable. These results may provide an explanation for the experimentally-measured growth rates that are about 6 - 8 times smaller than the ideal MHD growth rate. 4. Conclusions The hybrid simulations presented here show that, while ion FLR effects determine the linear stability properties of non-rotating FRCs, the inclusion of nonlinear and ion-toroidal-flow effects is necessary for a satisfactory description of plasma behavior in low-S* FRC experiments. In particular, it has been shown that the ion toroidal spin-up plays an important role in FRC nonlinear evolution, including that of the n = 1 tilt mode. The 3D hybrid simulations have been able to reproduce all major experimentally observed stability properties of kinetic (theta-pinch-formed) FRCs. Namely, the scaling of the linear growth rate of the n = 1 tilt instability with S*/E parameter has been demonstrated for a class of elongated elliptic FRCs [3]; and ion toroidal spin-up, the nonlinear saturation of the tilt mode, and the growth of the n = 2 rotational mode have been demonstrated. It has been shown that the loss of ions with a preferential sign of toroidal velocity due to the resistive decay of the poloidal flux results in the ion toroidal spin-up, which reproduces very well the experimentally-observed ion rotation. The time scale of the ion spin-up is determined by the flux decay time. The MHD version of the HYM code has also been used to study FRC formation by the counter-helicity spheromak merging in support of the SSX-FRC experiment [8], and contributed to interpretation of several puzzling experimental observations. A parameter scan, which is not easily accessible experimentally, has been carried out using numerical simulations, and it has provided a better qualitative understanding of the several of the experimental results. In particular, the persistence of the residual toroidal field, and the slower-than-MHD growth of the tilt instability have been shown to be related to the large plasma viscosity and line-tying effects in the SSX-FRC experiments. Acknowledgments This research was supported by DOE contract DE-AC02-76CH03073. Calculations were performed at the U.S. National Energy Research Supercomputing Center. References [1] Tuszewski M. 1988 Nucl. Fusion 28 2033 [2] Ishida A., Momota H. and Steinhauer L. C. 1988 Phys. Fluids 31 3024 [3] Belova E. V., Davidson R. C., Ji H. and Yamada M. 2004 Phys. Plasmas 11 2523 [4] Barnes D. C. 2002 Phys. Plasmas 9 560 [5] Belova E.V. et al 2000 Phys. Plasmas 7 4996 [6] Ohtani H., Horiuchi R. and Sato T. 2003 Phys. Plasmas 10 145 [7] Ono Y. et al 1997 Phys. Plasmas 4 1953 [8] Cothran C. D. et al 2003 Phys. Plasmas 10 1748 [9] Ono Y. et al 2003 Nucl. Fusion 43 649 [10] Eberhagen A. and Grossmann W. 1971 Z. Phys. 248 139 [11] Steinhauer L. C. 2002 Phys. Plasmas 9 3851 [12] Tuszewski M. et al 1982 Phys. Fluids 25 1696; 1988 Phys. Fluids 31 946 [13] Harned D. S. and Hewett D. W. 1984 Nucl. Fusion 24 201 [14] Yamada M. et al 1990 Phys. Rev. Lett. 65 721 [15] Ono Y. et al 1996 Phys. Rev. Lett. 76 3328 [16] Watanabe T.-H., Sato T. and Hayashi T. 1997 Phys. Plasmas 4 1297; Sato T., Oda Y. and Otsuka S. 1983 Phys. Fluids 26 3602 [17] Lukin V. S. et al 2001 Phys. Plasmas 8 1600 [18] Biskamp D. and Welter H. 1980 Phys. Rev. Lett. 44 1069 [19] Kornack T. W., Sollins P. K. and Brown M. R. 1998 Phys. Rev. E 58 R36 Figure captions Fig.1. Time evolution of (a) the normalized value of trapped poloidal flux, and (b) the normalized number of ions inside the separatrix obtained from 2D hybrid simulations with S* = 20 and E = 4. Fig.2. (a) Time evolution of the normalized angular momentum of all ions, and the ions inside the separatrix; and (b) the maximum value of the ion toroidal flow velocity obtained from same simulations as in Fig. 1. Fig.3. (a) Contour plots of the ion toroidal velocity in the r - z plane at t = 40 t_A and t = 80 t_A obtained from 2D hybrid simulations with S* = 20 and E = 4; and (b) Radial profiles of the ion toroidal flow velocity at the FRC midplane at t = 20, 40, and 80 t_A. The separatrix radius is R_s/R_c ≈ 0.6. Fig.4. Contour plots of (a) the pressure, and (b) the toroidal velocity at t/t_A = 10, 14, 20, 25 obtained from 2D MHD simulations of counter-helicity spheromak merging for µ = η = 0.001. The velocity is normalized to the characteristic Alfvén velocity; the contour values correspond to V_max = 0.4 and V_min = -0.2 for t = 10 - 20, and V_max = 0.2 and V_min = -0.1 for t = 25. Fig.5. Time evolution of (a) the normalized kinetic energy and (b) the magnetic energy obtained from 2D MHD simulations of counter-helicity spheromak merging for µ = 0.001 and several values of resistivity. Fig.6. Time evolution of (a) the normalized kinetic energy and (b) the magnetic energy (total and toroidal) obtained from 2D MHD simulations for η = 5 × 10^-4 and different values of viscosity: µ = 0.0005 (dotted line), µ = 0.001 (solid), and µ = 0.002 (dashed). Fig.7. Contour plots of the toroidal magnetic field at the poloidal plane obtained from 2D MHD simulations with η = 0.001 and (a) µ = 0.001 (t/t_A = 10, 20, 30); and (b) µ = 0.004 (t/t_A = 20, 30, 40). The maximum toroidal magnetic field value (normalized to the initial edge field) is shown for each plot. Fig.8. Time evolution of (a) the n = 1 mode energy, and (b) the total kinetic energy obtained from four 3D MHD simulation runs of counter-helicity spheromak merging for µ = 5 × 10^-4 - 4 × 10^-3. Fig.9. Time evolution of the maximum values of the poloidal and toroidal magnetic fields, the radial current density, the axial flow velocity, and the toroidal flow velocity obtained in 3D simulations with η = 0.001 and (a) µ = 0.001, and (b) µ = 0.004. External Distribution (05/16/05) Plasma Research Laboratory, Australian National University, Australia Professor I.R. Jones, Flinders University, Australia Professor João Canalle, Instituto de Fisica DEQ/IF - UERJ, Brazil Mr. Gerson O. Ludwig, Instituto Nacional de Pesquisas, Brazil Dr. P.H. Sakanaka, Instituto Fisica, Brazil The Librarian, Culham Science Center, England Mrs. S.A. Hutchinson, JET Library, England Professor M.N. Bussac, Ecole Polytechnique, France Librarian, Max-Planck-Institut für Plasmaphysik, Germany Jolan Moldvai, Reports Library, Hungarian Academy of Sciences, Central Research Institute for Physics, Hungary Dr. P. Kaw, Institute for Plasma Research, India Ms. P.J. Pathak, Librarian, Institute for Plasma Research, India Dr. Pandji Triadyaksa, Fakultas MIPA Universitas Diponegoro, Indonesia Professor Sami Cuperman, Plasma Physics Group, Tel Aviv University, Israel Ms. Clelia De Palo, Associazione EURATOM-ENEA, Italy Dr. G. Grosso, Instituto di Fisica del Plasma, Italy Librarian, Naka Fusion Research Establishment, JAERI, Japan Library, Laboratory for Complex Energy Processes, Institute for Advanced Study, Kyoto University, Japan Research Information Center, National Institute for Fusion Science, Japan Professor Toshitaka Idehara, Director, Research Center for Development of Far-Infrared Region, Fukui University, Japan Dr. O. Mitarai, Kyushu Tokai University, Japan Mr. Adefila Olumide, Ilorin, Kwara State, Nigeria Dr. Jiangang Li, Institute of Plasma Physics, Chinese Academy of Sciences, People's Republic of China Professor Yuping Huo, School of Physical Science and Technology, People's Republic of China Library, Academia Sinica, Institute of Plasma Physics, People's Republic of China Librarian, Institute of Physics, Chinese Academy of Sciences, People's Republic of China Dr. S. Mirnov, TRINITI, Troitsk, Russian Federation, Russia Dr. V.S. Strelkov, Kurchatov Institute, Russian Federation, Russia Kazi Firoz, UPJS, Kosice, Slovakia Professor Peter Lukac, Katedra Fyziky Plazmy MFF UK, Mlynska dolina F-2, Komenskeho Univerzita, SK-842 15 Bratislava, Slovakia Dr. G.S. Lee, Korea Basic Science Institute, South Korea Dr. Rasulkhozha S. Sharafiddinov, Theoretical Physics Division, Insitute of Nuclear Physics, Uzbekistan Institute for Plasma Research, University of Maryland, USA Librarian, Fusion Energy Division, Oak Ridge National Laboratory, USA Librarian, Institute of Fusion Studies, University of Texas, USA Librarian, Magnetic Fusion Program, Lawrence Livermore National Laboratory, USA Library, General Atomics, USA Plasma Physics Group, Fusion Energy Research Program, University of California at San Diego, USA Plasma Physics Library, Columbia University, USA Alkesh Punjabi, Center for Fusion Research and Training, Hampton University, USA Dr. W.M. Stacey, Fusion Research Center, Georgia Institute of Technology, USA Director, Research Division, OFES, Washington, D.C. 20585-1290 The Princeton Plasma Physics Laboratory is operated by Princeton University under contract with the U.S. Department of Energy. Information Services Princeton Plasma Physics Laboratory P.O. Box 451 Princeton, NJ 08543 Phone: 609-243-2750 Fax: 609-243-2751 e-mail: [email protected] Internet Address: http://www.pppl.gov