3D toroidal physics: Testing the boundaries of symmetry breaking^a)

There is increasing interest in fusion devices that incorporate 3D, symmetry-breaking effects in their equilibrium state and in the development of new simulation methods that can address these effects. The motivations for introducing the 3D magnetic fields include plasma confinement and stability optimization (stellarators), plasma control (tokamak edge-localized and resistive wall instabilities), self-organized states (reversed field pinch single-helicity states), and engineering/economic constraints (limited number of toroidal field coils, asymmetrical particle/energy sources and non-uniform ferrous steel structure near the plasma).

Symmetry is a compelling and attractive feature for magnetic confinement systems. It assures the existence of magnetic flux surfaces, good particle confinement from conserved invariants of the motion, and simplifies analysis, for example, allowing many issues to be dealt with in some depth using analytic theory. However, deviations from symmetry have become integral to the operation and design of toroidal devices and constitute one of the major remaining areas for innovation that can lead to improved plasma performance. These 3D and symmetry-breaking effects complicate analysis, both with respect to theoretical modeling and in the interpretation of experimental diagnostics. The goal of this tutorial paper is to describe some of the new methods and the issues that are important to understanding the physics of 3D systems.

Symmetry-breaking effects have been present in magnetic fusion confinement systems for many years. Simple magnetic mirror systems had, in principle, the highest inherent degree of symmetry (axial) of any devices, with deviations only present from coil positioning errors and construction tolerances. However, as these evolved to minimum-B and tandem mirrors, 3D effects were of interest^1–4 as configurations were changed to improve interchange stability and reduce end losses. The early tokamaks tested various levels of symmetry-breaking effects from different arrangements of toroidal field coils and non-symmetrical breaks in the conducting vacuum vessel. In some cases⁵ (ORMAK) a large number, 56 toroidal field (TF) coils, were used in attempt to make ripple strength negligible; for other designs, such as the original Alcator device, non-axisymmetric groupings of coils led to relatively high ripple strengths of several percent even near the plasma center.⁶ Discrepancies between ion temperatures derived from neutron measurements and those from charge exchange were attributed⁶ to ripple losses. Existing and planned tokamak designs now typically use in the range of 18–22 TF coils. Tokamaks have recently introduced 3D coils for control of edge localized mode (ELM) instabilities; in ITER, toroidally non-uniform distributions of ferritic steel from test blanket modules and structures related to the neutral beams will introduce low levels of long toroidal wavelength ripple.⁷ Stellarators have been 3D since their invention,⁸ first in order to create vacuum rotational transform, and more recently to improve confinement, stability, bootstrap current consistency, and turbulent transport. Some of the earliest stellarators had issues with poor flux surfaces and confinement⁹ caused by error fields and symmetry deviations; however, subsequent generations have improved flux surface integrity either by minimizing low order resonant (ɨ = 1/q = n/m) field errors, designing in more shear in the rotational transform, or, in the case of low shear configurations, avoiding low order rational surfaces altogether. Confinement has also been significantly improved by new optimization techniques, such as quasi-symmetry,¹⁰ and quasi-omnigeneity,^11,12 which will be described in more detail later. Finally, reversed field pinches, while nominally toroidally symmetric in their earlier forms, have recently demonstrated spontaneous helical 3D states,^13,14 which have led to longer discharge durations and the generation of internal transport barriers. These equilibria have internal helical structures that are similar¹⁵ to the “snake” structures that have been observed¹⁶ and modeled¹⁷ in the JET tokamak.

This evolving acceptance of 3D effects in magnetic fusion configurations is motivated by the fact that not all symmetry-breaking mechanisms count in the same way and plasma performance is tolerant of such effects up to some level. For example, long wavelength/resonant perturbations are generally more dangerous to the magnetic configuration than shorter wavelengths and phenomena such as collisions, turbulence, finite orbit width effects, ambipolar electric fields, and flow suppression of islands provide some shielding of 3D effects. In future fusion reactors, while the toroidal field coils will likely be further from the plasma (resulting in lower TF ripple), the effects of ports and other non-uniformity in the surrounding steel structures will introduce novel forms of ripple for which there is no current experience. 3D coils in reactors also are one of the few available tools for plasma control and optimization that do not involve large levels of recirculating power¹⁸ (assuming such coils can be superconducting). However, the lower plasma collisionalities and larger energetic particle populations that will characterize reactors imply levels of symmetry-breaking that are tolerable in current experiments may need to be reexamined in light of reactor parameters. This motivates continued advances in 3D simulation tools and diagnostics.

Fortunately, as a result of the historical development of stellarators, and more recent interest in 3D effects in tokamaks and RFPs, a good basis has been established for these efforts. The purpose of this paper is to provide a tutorial introduction to current and historical work on 3D toroidal physics in tokamaks, stellarators, and reversed field pinches. A basic toolkit of methods that can be applied in a unified way to all 3D systems will be described and current research trends discussed that are expected to lead to improvements in the modeling and simulation of these systems. The topics of 3D symmetries, design of 3D configurations, equilibrium, confinement and stability will be covered, with various examples given in each area. This is not intended to be a comprehensive review. Excellent reviews^18–21 for those seeking more depth in this area have recently been published.

This tutorial paper is organized as follows. In Sec. II, the extended forms of symmetry (field period and stellarator symmetry) that are relevant to 3D configurations are described. In Sec. III, the characterization of 3D systems and the design methods used to synthesize them are discussed. Section IV presents a basic formulation of the 3D plasma equilibrium problem and describes computational approaches that are being used to obtain equilibria. Section V first discusses methods for analyzing single particle confinement and then outlines their application to configuration optimization, energetic particle transport, neoclassical transport, and anisotropic transport in chaotic field line regions. Finally, Sec. VI describes the stability analysis of 3D systems, including energetic particle instabilities, micro-turbulence, stability of the 3D tokamak edge, and MHD stability of stellarators, with discussions included of possibilities for stability optimization.

The most general definition of symmetry is that an object observed from different viewpoints looks the same. Magnetic confinement systems such as simple mirrors and idealized axisymmetric tokamaks possess such symmetries about a single axis of rotation. This is denoted as continuous symmetry since such systems look the same for arbitrary angular rotations about a symmetry axis. The next level of symmetry is field period symmetry; this characterizes idealized stellarators and rippled tokamaks. This is a discrete form of axial symmetry and implies that the device will only look the same if it is rotated through discrete angular increments. For example, the LHD stellarator,²² with 10 field periods, would appear the same if rotated through a toroidal angle of 36°, a tokamak with 20 planar toroidal field coils, an angle of 18°, etc. The lowest form of field period symmetry would be one field period. This only implies periodicity over a 2π change in toroidal angle, which all toroidal devices must conform to. Recent examples of devices with one field period are ITER with the effects of ferrous materials in the test blanket modules (TBM) taken into account⁷ (isolated to a 120° sector of the tokamak) and DIII-D with TBM simulation coils²³ (localized to one toroidal location). Also, if random field errors are taken into account, one field period would be the default degree of symmetry for a toroidal device. At the other extreme, a perfectly axisymmetric tokamak can be thought of as having an infinite number of field periods. In Figure 1, examples of (a) continuous and (b) discrete symmetry are given. From a computational perspective, devices with higher field periods require less data to achieve the same degree of resolution, since equilibrium quantities only need be represented over one field period; alternately, in Fourier representations, this implies that expansions in toroidal mode number go in increments of the field period. Low field period configurations are then the most challenging to model, even in the case when the 3D effects are weak. Examples of these include compact stellarators, tokamaks with ELM coils, and tokamaks with TBM effects.

The second form of symmetry that characterizes toroidal systems is known as “stellarator symmetry.” This can be thought of as a generalization of reflection symmetry to helical systems. For stellarators, this implies an observer will see the same view looking one way around the torus as looking the opposite direction and standing on one's head [i.e., quantities at a point (R, ζ, and z) are the same as at (R, −ζ, and −z), where R is the major radius, ζ is the toroidal angle, and z is the elevation]. For tokamaks, this type of symmetry implies up-down symmetry for each cross-sectional slice at a fixed toroidal angle. Stellarator symmetry for 3D systems is not known to have any direct physics advantages,²⁴ but simplifies the analysis of such systems in the sense that only half of the volume needs to be taken into account in representations of equilibrium magnetic fields, coil geometries, etc. For example, in Fourier space, the magnetic field magnitude can be expanded as a series in cos(mθ − nζ), without the need to include sin(mθ − nζ) terms. Since most stellarators have tried to adhere to this type of symmetry (i.e., at least in the idealized limit, not including field errors), many of the modeling tools originally developed for stellarator physics assume stellarator symmetry. Tokamaks, routinely deviate from stellarator symmetry by tolerating up-down asymmetries, as for example, are the case with single-null divertors. Increasingly, stellarator analysis tools are being generalized to be applicable diverted tokamaks with 3D effects that violate stellarator symmetry.

The design techniques for 3D configurations have evolved from largely intuitive/semi-analytical approaches that begin with multipolar helical or window-pane coil concepts to computationally intensive optimization methods in high dimensional design spaces. The former is characteristic of the early stellarator designs and more recent tokamak 3D control coil designs; the computational approach has been used first for the W7-AS,²⁵ W7-X,²⁶ and HSX²⁷ devices and for recent modular stellarator design efforts. One of the strengths of 3D shaping is the very large design space that is available. The available 3D plasma shaping parameters are limited by the fact that field structures that can be produced by collections of magnetic dipoles outside the plasma drop off strongly as one moves into the plasma. Since approximately a 1-m separation is needed in a reactor between the coils and the plasma (for the blanket and shield) and engineering constraints will set some minimum radius of curvature for coils, this precludes producing short wavelength variations in equilibrium magnetic fields beyond some level. Estimates have been made²⁸ of the maximum number of parameters that can effectively be used by performing SVD (singular value decomposition) of coil-to-plasma shaping transfer functions; these studies have concluded that ∼30–40 shape parameters should be available. However, this is still a very large design space compared to axisymmetric tokamaks, where the 2D shape is typically determined by 3–4 parameters (aspect ratio, elongation, triangularity, indentation, etc.). For example, if one was to perform a thought experiment by arbitrarily quantizing each shape parameter into 10 levels, 2D systems (i.e., axisymmetric tokamaks) would have about 10⁴ possibilities, while 3D systems (optimized stellarators) would have 10⁴⁰ possibilities. While many of these may be of no interest or redundant, even if only a small fraction are attractive, the size of the 3D design space remains quite large.

Currently, the most advanced methods for 3D design use numerical optimization methods, such as Levenberg-Marquardt,²⁹ genetic evolution,³⁰ or differential evolution.³¹ Since 3D equilibria are uniquely determined by the outer flux surface shape in combination with the pressure profile plus either the rotational transform or plasma current profiles, one approach for designing 3D systems has been to vary the shape of the outer magnetic flux surface to minimize a range of physics and engineering target functions.³² The profile functions can either be fixed or varied. The target functions typically would include particle confinement, transport, MHD stability, micro-turbulence, aspect ratio, etc. Once a sufficiently optimized configuration is obtained, a set of modular coils can then be derived³³ on a specified coil-winding surface surrounding the plasma. These are first obtained as continuous sheet currents, and then steps are taken to discretize the current distribution, resulting ultimately in 3D current paths that can be practically realized as coils. Due to the fact that the discretized coils may not precisely reproduce the original optimized surface, tradeoffs must be made and the process may go through a succession of iterations. A more recent approach has been directly optimizing from the coil-set,³⁴ with its geometry described either by discrete points on the coils or by mathematical functions, such as splines or Fourier series. In this method, the external magnetic fields produced by the coils are used to compute a 3D plasma free-boundary equilibrium, physics/engineering target functions are calculated, the coil geometry and currents are varied, the equilibrium and targets are recomputed, and the cycle continues until a desired configuration is found.

There are a variety of ways to categorize 3D systems. For this section, the focus will be on methods used for producing rotational transform, since this is the first basic issue to be considered in a design effort. The earliest stellarators were designed with the goal of producing all of the needed rotational transform from the vacuum magnetic fields, allowing steady-state confinement without the need for externally driven plasma currents. The two primary methods for achieving this (other than internal toroidal plasma current) have been recognized^8,20,35 since the earliest development of stellarators: (a) using a helical magnetic axis, which introduces torsion, or (b) keeping the magnetic axis planar, but helically rotating the outer magnetic surfaces. In Figure 2, these approaches are shown, along with the case of a 3D tokamak, which uses only internal current to produce transform; here, TJ-II³⁶ given as an example of (a) and LHD²² as an example of (b). The rotational transform profiles (ɨ = 1/q, where q is the tokamak safety factor) associated with the configurations of Figure 2 are given in Figure 3. More recent optimized stellarators may use some mixture of these two approaches. Even devices such as LHD with perfectly planar magnetic axes in the vacuum state will tend to acquire helical deformations in the axis as the plasma beta is increased. As can be seen, the two stellarators produce rotational transform profiles that increase in going outward from the magnetic axis, while standard tokamak transform profiles decrease in going toward the plasma edge (i.e., q increases going outward from the axis). An intuitive way to understand how 3D shaping can produce vacuum rotational transform has been introduced,¹⁸ noting the analogy between magnetic fields (satisfying ∇•B = 0) and incompressible fluid flows (satisfying ∇•v = 0). An incompressible fluid flowing through these shapes would develop helical streamlines in a similar way as the magnetic field line trajectories follow helical paths. The two methods shown here for producing vacuum transform introduce separate tradeoffs. The helical axis approach is more efficient, leading to larger levels of transform (good for confinement of passing particles) at moderate aspect ratios and lower field periods. However, the planar axis approach produces more shear in the transform profile, which minimizes resonant magnetic island widths and can be stabilizing for some instabilities. As will be discussed further in Sec. V, both approaches have moderate to high levels of magnetic field ripple, which leads to higher levels of neoclassical transport and energetic particle loss. This has been improved by performing optimizations that target these physics characteristics. This approach has led to devices such as HSX,³⁷ W7-X,³⁸ QPS,⁴¹ and NCSX.⁴⁰ Another issue is that within each approach, if one wants to increase the rotational transform, while preserving similar shaping within each field period, there are two equilibrium scaling laws: ɨ/N_fp= constant and ⟨R⟩/⟨a⟩/N_fp = constant that are relevant. These imply that in order to increase transform one must increase the number of field periods and aspect ratio simultaneously.

An alternate design method that was identified for avoiding the trend to higher aspect ratio and/or stronger helical axis is the so-called hybrid approach, which accepts that some fraction of the transform can be produced by plasma current (but allowing a lower current level than in the equivalent tokamak) to augment the vacuum rotational transform present from 3D shaping. This plasma current could be provided either through the self-generated bootstrap current or driven by external sources (e.g., Ohmic, RF, and beams). This concept is used in the CTH experiment³⁹ and was incorporated to varying degrees in the NCSX⁴⁰ and QPS⁴¹ compact stellarator projects. Their characteristic flux surface shapes and cross sections are shown in Fig. 4 with the associated transform profiles given in Figure 5. As indicated, relatively low aspect ratios were achieved with moderate levels of rotational transform. The examples shown have stellarator-like shear in the iota profiles, but configurations were also found with tokamak-like shear. The QPS design incorporated a racetrack shaped 2-field period helical axis, while NCSX had closer to a planar axis in a 3-field period design.

Another approach is to use tilted planar toroidal field coils.^42,43 In the case of two field periods (the CNT device at Columbia University),⁴⁴ this produces vacuum rotational transform. As more tilted coils are added (as in the Proto-Circus device, also at Columbia),⁴⁵ a configuration results that can be considered a “rotational transform amplifier.” i.e., while it may not have vacuum rotational transform, as plasma current is added it produces higher levels of rotational transform (lower q) than a similar tokamak at the same aspect ratio. Calculations have shown that these configurations also can offer lower magnetic field ripple⁴⁶ than the equivalent tokamak with the same number of TF coils.

The 3D reversed field pinch relies on inductively-driven currents and spontaneous MHD self-organization to produce a single-helicity configuration^13,14 with a helical magnetic axis having low levels of q (very high levels of rotational transform) that can pass through zero and reverse near the plasma edge. These configurations derive most of their rotational transform from plasma currents, with only small components from the 3D effects; however, they have demonstrated improved plasma sustainment and internal transport barriers with the spontaneous 3D equilibrium states. Typically, the internal flux surfaces of single-helicity RFPs are strongly 3-dimensional near the axes, but then become more axisymmetric near the plasma outermost surface.

In the case of tokamaks, all rotational transform is produced from plasma currents. 3D perturbations are used for other purposes and can be either introduced intentionally (in the case of ELM coils) or of necessity (TF ripple, TBMs, ferritic materials, etc.). In the case of ELM coils, the so-called “window-pane” coils, as shown in Figure 6, are added; these surround the plasma and produce field components that are approximately normal to the plasma outer surface. The purpose of these coils is to locally break-up the outer flux surfaces, forming an edge stochastic layer, as illustrated in the right-hand image in Figure 6; the goal is to limit the strong plasma pressure gradients in the pedestal region. This can suppress ELM instabilities, which, in larger systems such as ITER and reactors, will need to be avoided since they can lead to severe transient heat loads on plasma-facing components. Figure 7 shows the types of perturbations in the outer magnetic flux surface of a tokamak that (a) ELM coils can produce and (b) TF ripple + TBMs (test blanket modules) can produce. Here, the colors and contour lines indicate magnetic field strength and the non-axisymmetric shape changes have been enhanced by a factor of 50.

Accurate 3D toroidal magnetic equilibria, self-consistently including the effects of the plasma pressure and currents, form the basis for all subsequent plasma physics analysis and modeling. While, as will be explained, 3D equilibria offer a number of unique challenges, there are now efficient computational tools available that provide equilibria over a hierarchy of approximations and function in nearly as routine a manner as axisymmetric tokamak equilibrium models.

The 3D equilibrium problem was first analyzed mathematically by Grad⁴⁸ and Lortz.⁴⁹ They determined that the existence of solutions of Eqs. (1)–(3) could only be guaranteed for the special cases of closed field line (bumpy torus) and axisymmetric (tokamak) systems. This issue continues to be pursued, with recent work by Weitzner⁵⁰ showing that equilibria could be obtained in a simple topological torus model through expansions in a small parameter measuring the deviations of pressure surfaces from a simple nested surface structure. Basic differences between 2D and 3D equilibria arise from the characteristics of the governing partial differential equations (PDE). Axisymmetric equilibria are described by the Grad-Shafranov equation, which is a nonlinear elliptic 2D PDE, and there are a number of accurate, efficient solution methods available, along with many useful analytic approximations. 3D equilibria, however, lead to a mixed nonlinear hyperbolic/elliptic system. As discussed in Ref. 50, equation systems of this type only appear in a few other applications (transonic flow, lower hybrid wave propagation) and general solution methods are relatively undeveloped. 3D magnetic equilibria also can have singular current sheets at rational surfaces, which lead to magnetic islands, and chaotic field line regions.

The existence of singular current sheets and magnetic islands in both stellarators and 3D tokamaks has stimulated development of a number of models beyond VMEC. These aspects of 3D equilibria become especially important as the plasma β is increased. The PIES (Princeton Iterative Equilibrium Solver)⁵⁸ is a direct solution approach based on the Picard iteration algorithm. PIES iterates between computing currents and magnetic fields, requiring some degree of averaging over old and new solutions to maintain numerical stability. The HINT-2 model⁵⁹ is also a direct solver, iterating between maintaining pressure constant along field lines, calculating magnetic fields, and currents. HINT-2 uses rotating helical coordinates, which are an efficient representation for heliotrons. The SPEC model⁶⁰ uses a piecewise continuous approximation in the radial flux variable and avoids current sheets by setting the pressure gradient to zero on rational surfaces. SIESTA⁶¹ is an extension of the VMEC variational approach to include islands. This is achieved by generalizing the VMEC independent variables to include a third magnetic field component normal to magnetic flux surfaces and a 3D pressure function. Small levels of resistivity are introduced to allow changes in magnetic field topology required for island formation. SIESTA utilizes a direct solver, based on coordinates determined by VMEC. A stellarator equilibrium and associated Poincaré plots indicating several island chains based on the SIESTA model is shown in Fig. 8.

For the case of tokamaks with weak 3D effects, various models have recently been developed. The simplest approach is to simply superimpose vacuum 3D magnetic fields on a conventional 2D equilibrium. This model has often been used for fast ion loss calculations in tokamaks with 3D perturbations. However, this approach misses self-consistent 3D modifications of flux surface shapes and plasma shielding/amplification effects. Beyond this, several models have used MHD equations to calculate the response of initially 2D tokamak equilibria to 3D fields that are imposed as inhomogeneous driving terms. Such models include the M3D-C1,⁶² MARS-F,⁶³ and IPEC⁶⁴ solvers. MARS-F and IPEC are linear, single fluid, ideal codes, while M3D-C1 is resistive, can be run in either linear or nonlinear modes and can use either single or two fluid models with plasma rotation. Plasma rotation/flows can have important effects on limiting the size of magnetic islands.^65,66 Such models so far have only been applicable to tokamaks since they are based on the assumption of small perturbations away from axisymmetric equilibria. A code comparison/benchmark activity⁶⁷ has been initiated to compare predictions of a number of these models for 3D tokamak equilibria. A more recent comprehensive review⁶⁸ has compared results from IPEC,⁶⁴ MARS-F,⁶³ M3D-C1,⁶² VMEC,⁵¹ NSTAB,⁶⁹ and HINT-2.⁵⁹ Also, high resolution magnetic probe arrays have been recently installed⁷⁰ on DIII-D to provide experimental validation data for this activity. Calculation of 3D effects is especially important in the edge/pedestal and divertor scrape-off regions of tokamaks. Understanding transport in the edge/pedestal region and nearness to peeling-ballooning stability boundaries is critical to sustainment of enhanced confinement regimes and mitigating edge-localized instabilities. The goal of 3D perturbations from external ELM control coils is to regulate the pedestal pressure gradients so that ELMs are suppressed without destroying confinement. In future, larger tokamaks ELMs can transfer large transient heat loads to plasma-facing components and will need to be suppressed. The observations that as 3D fields are applied some experiments see ELMs suppressed,^71,72 while others see ELM triggering,⁷³ make the pedestal region both of great interest, but also remains a significant challenge for modeling. In the outer divertor scrape-off region, 3D fields lead to strike point splitting,⁷⁴ novel structures known as homoclinic tangles^75–78 and possible loss of detachment at the divertor plates. Due to the extremely high scrape-off layer heat fluxes⁷⁹ anticipated in future tokamaks, good modeling of this region is also quite important. In addition to the substantial recent interest in modeling the edge and scrape-off regions, improved diagnostics⁸⁰ are now available, based on soft x-ray emissivity measurements. In order to resolve detailed features, energy filtering, reconstruction, and image analysis techniques are employed; also, by alternating the sign of the 3D magnetic field perturbations, differencing methods can be used to improve the contrast of the images. A recent example of a comparison⁸¹ between the predicted effects of 3D perturbations, using the M3D-C1 model, and the emissivity measurements is shown in Fig. 9.

With the breaking of toroidal symmetry, conserved invariants of particle motion, such as the canonical toroidal angular momentum (P_ϕ), are no longer available. In its place, the longitudinal invariant J = m∫v_||dl plays an analogous role as P_ϕ; bounce-averaged drift velocities can be derived from the invariant function J by: v_d = (1/eτ_bB²) B × ▿J, where τ_b = bounce time = ∫dl/v_||. The longitudinal invariant has been extensively used in the analysis of neoclassical transport, plasma stability, and as an optimization target function. If the energy and pitch angle are fixed, then contours of constant J represent the puncture points orbit trajectories onto a 2D constant toroidal angle plane analogous to the use of constant P_ϕ contours to project tokamak orbits onto a toroidal slice. This property was used in one of the first applications⁸² of the STELLOPT stellarator optimization code, and is known as quasi-omnigeneity. For systems with low levels of rotational transform per field period, J can be simplified⁸³ to J*, which replaces the integration along field lines with an integration along the toroidal angle at fixed values of poloidal angle and magnetic flux. The use of J for particle orbits in 3D systems is somewhat more complex than (P_ϕ) for several reasons. First, the evaluation of P_ϕ is more straightforward (without the need to perform field line integrals), while J must be computed by performing integrals along field lines for a sequence of pitch angles (also for a range of energies, if an ambipolar electric field is included). Second, for configurations with multiple magnetic wells (trapping regions) along field lines, multiple branches of J must be considered and transition probabilities between them taken into account.⁸⁴ For these reasons, most recent studies of confinement in 3D configurations have been based on direct integrations of the particle equations of motion. The guiding center trajectory equations have been derived in both canonical and non-canonical formulations. These form the basis for many areas of simulation for 3D systems, including neoclassical transport, energetic particle confinement/stability, hybrid fluid/particle models, and gyrokinetics.

Due to the form chosen for the Hamiltonian, the poloidal/toroidal angular coordinates (θ, ζ) have been chosen, so that the magnetic field lines in each magnetic flux surface are straight lines, and the particle trajectory equations only involve B (magnitude of the magnetic field) and its derivatives. Such coordinates are normally referred to as Boozer or magnetic coordinates,⁸⁵ and result in the Jacobian = 1/ $g=B2/(G+ɨI)(dχdρ)$⁠, where ρ is the magnetic flux surface label. This orbit model can also be extended⁸⁶ to include perturbation field components with the form $δB=∇×(αB)$⁠, as might be used to model magnetic islands or the presence of MHD or Alfvénic instabilities. The advantage of this choice is that less data are required to represent the magnetic field (i.e., only B and its derivatives are required, rather that 3 components of B with 3 derivatives of each component). Since this also means that fewer field evaluations must be made at each orbit step and memory requirements are smaller, it can increase efficiency of numerical trajectory calculations. Also, as will be described later, this form for the orbits allowed the recognition of the magnetic field structures known quasi-symmetries, in which the net radial drift over a bounce orbit could be brought to zero, leading to confinement characteristics similar to those of an axisymmetric tokamak. The disadvantage of the orbit model of Eqs. (7)–(10) is the need to transform to a particular choice of coordinates that are best suited to the case of equilibria with nested, closed flux surfaces. Methods for extending such coordinates outside of the last closed flux surface are discussed in Ref. 87. Also, Boozer coordinates are generally not optimal for representation of the geometries of plasma-facing structures, such as metallic walls, divertor plates, etc. at the plasma edge.

The above two orbit models are based on the guiding center approximation. For energetic particle populations and regions of the plasma where parameters vary rapidly (e.g., edge and divertors), a number of models^89,102 are beginning to incorporate full Lorentz model orbits (i.e., including the rapid gyromotion). A related approach⁹⁰ is to use a hybrid of guiding center and Lorentz, with guiding center orbits followed in the central region of the plasma and then switching over to Lorentz orbits near the edge.

The orbit equations result in four ordinary differential equations (ODEs) per particle (six in the case of Lorentz), can be solved with a variety of methods (e.g., Runge-Kutta, variable time step, symplectic) and are readily adapted to parallel computers. For applications to transport and energetic particle confinement, the time integration of the orbit trajectories is periodically paused and discrete Langévin (Monte Carlo) collisions⁹¹ are applied to the velocity space variables in order to simulate drag and diffusion from Coulomb collision processes. 3D configurations introduce a larger number of orbit classes than axisymmetric systems. Several examples are shown in Figure 10 based on integrating the Hamiltonian orbit equations. Figures 10(a) and 10(b) display several orbit classes for the LHD stellarator, including passing, trapped-passing transitional (trapped in some field periods, but not others), and deeply trapped orbits that are localized to the magnetic field minima and precess around the torus following these minima. For these plots, an ambipolar electric field was present; otherwise the latter two classes of orbits would have been rapidly lost. Typically, the transitional and deeply trapped orbits have the highest loss rates due to uncompensated radial drifts and lead to the ripple neoclassical transport scaling (∼1/ν, where ν = collision frequency), which are unique to 3D systems. In Figures 10(c)–10(e), several orbit classes that are important for 3D tokamaks are displayed. This is an ITER case that has both TF coil ripple and ripple from iron in adjacent TBM. The red colored trajectories alternate between being locally ripple trapped and trapped over several ripple wells due to the interaction between the TF and TBM ripple effects. The purple orbit trajectory is a more normal trapped orbit, bouncing, and precessing within the 1/R variation of |B|. However, due to the ripple, there is a slight jitter in the location of its banana tip locations. The banana tips do not remain attached to the original flux surface, as would be the case in a perfect axisymmetric system, and this leads to a form of ripple transport⁹² that can be important for energetic particle species.

The simplified dependency on B of the orbit equations in magnetic (Boozer) coordinates coupled with flexibility of 3D shaping has led to significant progress in confinement improvement in 3D systems. In order for optimization methods to target confinement, it was necessary to develop proxy functions that both were accurate reflections of the confinement physics and efficient to evaluate. An early approach to this problem was to target the extremal values (along the field lines) of the magnetic field strength. By recording the minimum and maximum values of B within a field period along each field line passing through a fixed toroidal angle cross section, 2D contours⁸³ of these B_min and B_max values could be plotted and compared with magnetic flux surfaces. By aligning the B_min contours with adjacent flux surfaces, confinement of deeply trapped particles could be improved, while aligning B_max contours with flux surfaces lessen the number of transitional, partially trapped/partially passing orbits. A subsequent concept was then to align contours of the longitudinal invariant, J, with flux surfaces. Such configurations are known as omnigenous⁹³ and offer the possibility that, by using evaluations of J at several different values of the pitch angle parameter λ = μ/ε, confinement could be improved over a larger region of velocity space than possible with the B_min and B_max optimizations. If exactly obtainable, such an optimization strategy could reduce bounce averaged cross-field drifts to zero (they would still be locally non-zero) and remove⁹⁴ the ripple transport (D ∼ 1/ν) regime scaling that is discussed in Sec. V D. The J optimization approach has also been shown to be equivalent to maintaining an equidistant separation⁹⁵ between contours of constant B within each flux surface. An example of applying the omnigenous optimization target is given in Ref. 32, where the related function J*, which replaces the integration along field lines with an integration along toroidal angle, was used. Unfortunately, configurations that exactly satisfy the omnigenous criterion have not been found, but this optimization target has proven useful³¹ in guiding the optimizer towards shapes that offer confinement improvement.

The next level of confinement improvement¹⁰ is the quasi-symmetry approach; these require that the magnetic field has the form B = B(ρ, mθ − nN_fpζ), where m and n are integers, N_fp is the number of field periods, ρ is a flux surface label, and θ, ζ are the poloidal, toroidal angles in Boozer coordinates. Coordinate-free versions of this criterion have also been derived.^96,97 This magnetic field structure implies that an ignorable coordinate is present, leading to a conserved canonical momentum, similar to an axisymmetric tokamak. If this could be achieved, all of tokamak transport theory could be applied directly, through a procedure known as an isomorphic transformation.⁹⁸ Hypothetically, three separate quasi-symmetries can be considered: toroidal (n = 0), poloidal (m = 0), and helical (m, n both non-zero). However, it has not been possible to exactly achieve any of these symmetries in realistic devices and this fact has been supported by theoretical analyses.⁹⁹ Poloidal quasi-symmetry cannot be attained exactly in a toroidal system due to its incompatibility with toroidal curvature. However, racetrack shaped two field period configurations, such as QPS,⁴¹ and those with long straight sections joined by short curved regions can approach¹⁰⁰ this type of variation in the straight sections and reduce the losses of trapped particles. Quasi-toroidal symmetry is likely to some extent incompatible with vacuum transform since the only two known methods of creating vacuum transform (discussed in Sec. III) involve deviations in real space from toroidal symmetry. However, hybrid configurations, such as NCSX,⁴⁰ have achieved a good degree of quasi-toroidal symmetry by allowing a portion of the rotational transform to come from plasma current. A high degree of quasi-helical symmetry has been achieved in the design of the HSX device, which uses n = m = 1 and N_fp = 4. In principle, other forms of quasi-helical configurations are possible by choosing different ratios of n/m than 1. As was the case with the omnigenous approach to optimization, targeting quasi-symmetry has yielded significant confinement improvement; even though these criteria have not been achieved in any exact sense.

Confinement studies of the energetic particle (EP) component (i.e., beams, RF tails, alpha particles, runaway electrons) are an important application of the above orbit models in 3D systems due to concerns about loss of heating efficiency and the possible transfer of substantial energy fluxes onto plasma facing components. Since EPs are relatively collisionless, they are more sensitive to symmetry-breaking effects than thermal plasma components. For the case of tokamaks with 3D effects, the impact of magnetic field perturbations from localized test blanket modules have been studied both for ITER and in DIII-D, where a toroidally localized set of coils were used to simulate the effect of test blanket module ripple. For ITER, the field perturbations arise from the fact that test blanket modules will contain significant masses of ferritic steel. Ferritic steel alloys are likely to remain present in future fusion reactors due to their superior properties¹⁰¹ in high temperature/high radiation environments. A variety of orbit Monte Carlo slowing-down models have been developed and applied to both ITER with TBMs and to DIII-D with TBM simulation coils. These follow a large number of alpha particle (ITER) or beam ion (DIII-D) orbits in the presence of energy and pitch angle scattering. Heat loads are accumulated as orbit trajectories leave the plasma and intersect the wall. An example of modeling the DIII-D TBM coil-induced losses is shown in Figure 11 using the Spiral full Lorentz orbit model.¹⁰² A variety of models have been applied to this case, and show generally good agreement with infrared camera measurements of localized heating in the graphite tiles that are present near the TBM mock-up coil location. Based on this validation, many of the same models^7,90,163 have been applied to the ITER TBM case. An example¹⁶³ of localized alpha particle losses for the ITER TBM case is shown in Figure 12. Under quiescent plasma conditions (no instabilities taken into account), the predicted TBM induced heat fluxes of escaping alpha particles incident on the first wall have remained under the prescribed limits (for ITER, this was set at 50 kW/m²). Efforts in this type of modeling are continuing as the effect of various plasma instabilities are taken into account and also losses are followed down into the divertor chamber.

In the case of stellarators and RFPs, EP confinement has been extensively evaluated with respect to its impact on heating efficiency, first wall heat loads, and for comparison with fast ion loss diagnostics. For heliotron/torsatron configurations, such as LHD, ATF, and CHS, significant improvements in EP confinement are possible by shifting the plasma inward. This lines up the minimum values of B as one follows a field line, lowers the ripple strength and results in better-centered fast ion orbits; such effects were originally noted in Ref. 103 using analytical field models and are known as sigma-optimization. This effect was analyzed for the LHD stellarator with the GNET global particle simulation model¹⁰⁴ and indicated about a factor of ∼3 improvement in trapped fast ion densities as the magnetic axis is shifted in from 3.75 m to 3.53 m. Similar improvements were seen in experimental charge exchange measurements of the fast ion component. In parallel, this form of orbit optimization also applies to the neoclassical transport properties of the thermal plasma. GNET also includes capabilities for simulating RF heating and has been applied to both to ICRF ion tails and ECRF produced suprathermal electron populations. Fast ion and alpha particle losses have been evaluated in a variety of more recent modular coil optimized stellarators. Neutral beam losses were calculated¹⁰⁵ using 5D Monte Carlo methods for the compact NCSX stellarator, indicating 10%-15% energy losses. Alpha particle losses have been evaluated in a number of stellarator reactor designs,^106–108 indicating that with sufficient optimization, very low levels of alpha loss are attainable. Neutral beam losses were evaluated in W7-X¹⁰⁹ and studies made of their sensitivity to plasma β and adjustments of in the main and auxiliary coil currents, indicating that improved beam ion confinement results as the plasma β is raised. Studies of fast ion confinement for LHD in the presence of MHD¹¹⁰ and Alfvén instabilities¹¹¹ have been done using both experimental scintillator-based lost-fast ion probe data (SLIP) and modeling with a guiding center orbit trajectory code (DELTA5D¹⁶³) that includes fluctuating TAE mode structures with down-chirping frequencies. In both the measurements and modeling, two loss regimes were found: at low TAE amplitudes, the fast ion flux to the wall scaled as δB_θ,TAE (convective), while at higher TAE amplitudes the losses scaled as (δB_θ,TAE)² (diffusive). In the case of RFPs, neutral beam injection has been used on the MST device and studies made of fast ion confinement in both the conventional axisymmetric¹¹² and the spontaneous single helicity (SHAx) modes¹¹³ of operation. Higher losses of beam ions were observed in the single helicity mode, especially for the hydrogen ions, due to the symmetry-breaking (ripple) effects in the central region of the plasma.

3D configurations introduce unique neoclassical transport physics characteristics not present in axisymmetric tokamaks. These are a consequence of the uncompensated radial drifts that characterize deeply trapped particles and the large orbit excursions that are present in the trapped-passing region of velocity space. Also, the breaking of axisymmetry requires that ambipolarity (Γ_ion = Γ_electron) must be explicitly imposed, leading to a nonlinear equation for the radial electric field; in the tokamak ambipolarity is automatically satisfied and the electric field is only explicitly determined if higher order gyro-viscous effects are introduced or if plasma rotation is included.¹¹⁴ 

An example of the different transport regimes in a 3D configuration is shown in Figure 13. This is based on a calculation of the monoenergetic density diffusion coefficient done with the DKES (Drift Kinetic Equation Solver)¹¹⁵ for the LHD stellarator with a finite level of ambipolar electric field included. The Pfirsch-Schlüter (PS) and plateau regimes are present in both stellarators and tokamaks. The ripple transport (or 1/ν regime) is unique to 3D systems; this scaling with collision frequency becomes dominant when collisions are sufficiently infrequent that the uncompensated radial drifts of the deeply trapped orbits and excursions of transitional orbits are active at transporting them away from flux surface in between collisions. In this regime, increasing collisionality actually improves confinement by limiting the movement of these classes of orbits away from their initial flux surface locations. An analogy can be made to the way in which a large crowd of people in a room can slow down an individual's progress towards an exit door. Next, as collisionality is lowered even further, such orbits experience weak E × B and B × ∇B drift components that can compensate the radial drifts, causing their orbits to close poloidally, and leading to the superbanana regime. In this regime, transport scales with collisionality in a more normal way, which can be related to a collision frequency times an effective step size. In between the ripple and superbanana regimes isolated regimes with different collisionality scaling enter in as D transitions smoothly from a D ∝ 1/ν to a D ∝ ν scaling. In Figure 14(a), a similar plot is given for the monoenergetic density diffusion coefficient in LHD, showing the dependencies on both the collisionality and ambipolar electric field parameters (ν/v, E_r/v). At E_r = 0 or small finite E_r, there is no low ν superbanana regime for the collisionalities examined here, only the 1/ν ripple regime. As E_r is increased, a superbanana regime enters in, and as it is further increased, the 1/ν ripple regime eventually disappears. As indicated, for this case, and for most stellarators, there is a very strong variation with E_r at low collisionalities. In Figure 14(b), the same calculation is done for a tokamak (DIII-D in this case) with ELM coil 3D perturbed fields included. As can be seen, the stellarator-like dependence of D on electric field only appears at low collisionalities, well below the plateau regime. Also, the general level of the diffusion coefficient is lower at equivalent collisionalities than for the stellarator case.

A simple and convenient way to characterize low collisionality neoclassical transport across different 3D systems is the effective ripple¹¹⁶ parameter ε_eff. This provides a single quantity that contains the configuration dependent part of the 1/ν asymptotic E_r = 0 ripple monoenergetic transport coefficient D/D_plateau = (4/3π)²(2ε_eff)^3/2/ν^*, where D_plateau is the plateau value of the diffusion coefficient¹²⁵ and ν^*= R₀ν/(ɨv), and thus is useful for comparing different configurations. In Figure 15, ε_eff^3/2 is plotted for a number of stellarators, several tokamaks with 3D effects, and an RFP single helicity (SHAx) state. As indicated, the earlier stellarators had higher values (10⁻³< ε_eff^3/2< 10⁻¹) than the more recent optimized stellarators (10⁻⁴< ε_eff^3/2< 10⁻³); tokamaks with only TF ripple are quite low, while those with ELM coils are about a factor of 10 below the best optimized stellarators. The RFP single helicity case is in the range of the un-optimized stellarators near the center of the plasma, but drops to the level of 3D tokamaks near the edge, reflecting the localization of the helical core to the center region, accompanied by a nearly axisymmetric outer region.

As mentioned above, neoclassical transport in 3D configurations is not automatically ambipolar (Γ_i, ion particle flux, = Γ_e, electron particle flux); the balance of Γ_i with Γ_e in a 3D system will only occur for particular values of electric field for which ∑e_aΓ_a = 0. Figure 16 shows an example of this phenomenon, taken from a transport analysis¹¹⁷ of the HSX stellarator. For this case, there are three points where Γ_i = Γ_e. The usual naming convention has been to identify the root with inward pointing electric field (E_r < 0) as the ion root since it can be related to the ion flux rising to meet the electron flux, while the farthest right hand root (E_r > 0) is the electron root since the electron flux is dropping to meet the ion flux. In between there is often a third root that is not stable due to the form of the electric field equation: ∂E_r/∂t ∝ Γ_e − Γ_i (this is derived from the radial component of Ampere's law). The instability of this third root is related to the fact that if E_r is perturbed to a larger (smaller) value, then Γ_i decreases (increases) and Γ_e increases (decreases); this further increases (decreases) E_r, resulting in linear instability. However, for the points identified as ion and electron roots, the responses to E_r perturbations are in a restoring/stabilizing direction; increasing (decreasing) perturbations in E_r lead to decreasing (increasing) changes in ∂E_r/∂t. Due to the local diffusive approximation of the neoclassical transport model described above, the ambipolarity condition must be solved on each flux surface. There has been interest¹¹⁸ in scenarios in which the ion root is present over part of the plasma and the electron root in other regions since this would create strong localized electric field shearing that could play a role in transport barrier formation. Methods for driving the ambipolar balance to the electron root have been of interest since the electric field is larger causing the particle and energy fluxes of both species tend to be lower at the electron root than the ion root; also, impurity accumulation has been predicted^119,120 to be less of an issue for the electron root.

The above discussion of neoclassical transport is focused for simplicity on the monoenergetic density transport coefficient D₁₁ that goes with the density gradient. The full neoclassical transport description involves at least a 3 × 3 matrix of coefficients and includes coefficients that go with temperature gradient and electric field terms for both density and heat transport as well as bootstrap current. A number of different models have been developed to obtain these transport coefficients in 3D systems, along with earlier analytic theory;¹²¹ research in this area is ongoing and will be described only in a general way here. As referred to above, the DKES model¹¹⁵ provides this matrix of coefficients for arbitrary values of collisionality and electric field. This model and other early approaches utilized various approximations to reduce the dimensionality of the kinetic equation. DKES used a radially localized approach, incorporated a pitch-angle scattering collision operator, assumed a distribution function close to Maxwellian and made an incompressible approximation to the E × B velocity: E × B/B²≈ E × B/². This allowed a variational formulation to be derived in which the energy and radial variable were parameters, providing upper/lower bounds on transport coefficients. The kinetic equation was solved using Fourier expansions in poloidal/toroidal angle and Legendre polynomials in pitch angle. The particle and energy fluxes and the bootstrap current could then be obtained by performing energy integrations over the monoenergetic coefficients. Another approach that was contemporaneous with DKES was the GSRAKE model,¹²² which used a ripple-averaged kinetic equation to reduce the dimensionality. GSRAKE required a specific model for the magnetic field, the multi-helicity model¹²² to perform the averaging. In addition, particle based Monte Carlo models^123,124 have been developed to obtain local diffusion coefficients, or, in some cases, to test the diffusive transport approximation.¹²³ These different approaches have been benchmarked for a range of different 3D configurations and covered in detail in Ref. 125.

As mentioned above, the early transport calculations were based on including only a Lorentz pitch-angle scattering operator. This prevented the conservation of momentum at a microscopic level and precluded accurate prediction of plasma flows and bootstrap currents. While inclusion of a full momentum-conserving collision operator would have made the early methods computationally unwieldy, moments methods were developed^126–129 for correcting this issue at a macroscopic level. This allowed the earlier transport coefficients based on the DKES model to be used and momentum conservation to be achieved at least at a fluid moment level. These approaches were extended to include impurities,¹²⁸ higher order energy moment expansions,¹²⁹ and incorporated into computational models, such as PENTA.¹³⁰ As optimized quasi-symmetric stellarators were developed, interest increased in predictions of plasma flows¹³¹ and bootstrap currents.¹³² Plasma flow directions and magnitudes were a distinguishing feature of different optimization approaches and a reliable prediction of bootstrap current is important for stellarator hybrids, which derive significant fractions of their rotational transform from plasma current. Also, stellarators, such as W7-X, with island divertors that require well-defined levels of rotational transform at the plasma edge need accurate evaluations of bootstrap current. More recently, further extensions to stellarator neoclassical theory have been made^133,134 (SFINCS model) that allow 4D models (toroidal/poloidal angles, pitch angle, energy) to be solved directly with momentum-conserving collision operators and full drift trajectories. This has been compared¹³³ with the DKES model, indicating both regions of agreement and disagreement, depending on the size of the electric field. Another development has been inclusion¹³⁵ of poloidal electric field effects, which can be important, especially in the case of impurity transport. Finally, neoclassical transport theory^136,137 specific to the regimes of 3D tokamaks (δB/B < ρ_ion/a) has been derived and used to obtain the neoclassical toroidal viscosity; effects on toroidal viscosity have been of significant interest in tokamak experiments since toroidal flows are readily measured and may influence MHD stability by increasing the probability of locked modes.

As discussed in Sec. IV, there is no guarantee that 3D equilibria will have nested closed flux surfaces. In particular, there are areas near the edge of the plasma and in the divertor scrape-off region that are of significant interest for modeling, but in which the magnetic field lines are chaotic and large magnetic islands are present. When this is coupled with the highly anisotropic heat conductivities (χ_||/χ_⊥ may exceed 10¹⁰) that characterizes fusion systems, an entirely new transport regime must be addressed. In these regimes, the normal cross-field transport that is addressed by neoclassical transport becomes largely irrelevant and transport is dominated by parallel conduction along the wandering field lines. Since these fields are non-integrable, it is impossible to construct discrete coordinate meshes that are well-aligned with the magnetic field structures. Also, attempts to solve such highly anisotropic transport equations using finite differences or finite element methods experience numerical pollution and eventually are limited by an accumulation of errors. These facts generally force such approaches to use lower values of χ_||/χ_⊥ than are realistic. A new paradigm, the Lagrangian Green's (LG) function method,¹³⁸ has been developed to address these challenges. This involves transforming the parallel heat equation into a non-local closure using a Green's function solution. In this form, the standard diffusive transport model can be solved efficiently and accurately even for fields with complex filamentation and braiding;¹³⁹ the technique can also be readily generalized to other fractional diffusion models. Applications have been made to heat pulse propagation¹⁴⁰ in chaotic fields, and non-diffusive transport characteristics analyzed. It has been shown that partial transport barriers are possible even with chaotic magnetic fields.

3D configurations share many of the same stability issues (pressure/current driven MHD, microturbulence, energetic particle instabilities) as axisymmetric systems. However, they also provide unique phenomena and challenges of their own and offer new opportunities for control and optimization of stability. In this section, some of the more recent topics of interest are discussed, along with the tools/methods used for stability analysis.

Fast ion driven instabilities were first observed¹⁴¹ in the W7-AS stellarator, driven by neutral beam injection. This was also the first use of ECE (Electron cyclotron emission) Imaging to diagnose the 2D mode structure of these instabilities. They were dominated by (m, n) = (3, 1), (5, 2), (2, 1), and (5, 3) modes that tracked variations in the ɨ profile and followed an Alfvénic frequency scaling [ω ∝ 1/(m_ion n_ion)^1/2] as the plasma ion mass density was varied. Since then there have been many observations of both fast ion and suprathermal electron driven Alfvén modes in different stellarators, including CHS,¹⁴² LHD,¹⁴³ TJ-II,^144,154 H-1N,¹⁴⁵ Heliotron-J,¹⁴⁶ and HSX.¹⁴⁷ In addition, Alfvénic modes have been observed in reversed field pinches,¹⁴⁸ both in the normal operational mode and in single helicity 3D states. As with tokamaks, interest in these instabilities is motivated by concerns about damage to plasma-facing components in reactors, loss of heating efficiency, and diagnostic use (e.g., inferring properties of the ɨ profile and dynamical changes in the ion density).

In similarity with Alfvén instabilities in axisymmetric systems, the analysis of these instabilities begins with calculations of the Alfvén gap structure,¹⁴⁹ identification of eigenmodes with a global structure^150,151 within open gaps, perturbative stability analyses, and, more recently, development of hybrid fluid-kinetic and gyrokinetic models. The modes can be classified into several categories: modes in continuum gaps produced by poloidal couplings (TAE, EAE), modes in continuum gaps produced by helical and toroidal mode couplings (HAE, MAE), extremal modes that occur above or below a continuum minimum/maximum (GAE, RSAE), and beta-induced and acoustic coupled modes (BAE, BAAE). The naming conventions arise from the type of mode interaction responsible for the Alfvén spectral gap within which the mode resides. The “AE” stands for Alfvén Eigenmode and the prefixes T = toroidal, E = elliptical, H = helical, M = mirror, B = beta, BA = beta acoustic. The modes created by helical/toroidal mode couplings (HAE, MAE) are unique to 3D configurations and are present due to the toroidal variations in the equilibrium magnetic field.

This leads into the topic of mode families,^152,157 which is basic to all aspects of macroscopic stability analysis of 3D systems. For the axisymmetric tokamak, the toroidal mode number may be considered a good quantum number since the equilibrium fields include only n = 0 and provide no coupling between instabilities having different n values. Linear instabilities at arbitrary values of n may then be considered independently; they will only couple with each other nonlinearly. In the case of 3D systems, the equilibrium has a field-period periodicity (as described in Sec. II), i.e., it will include not only n = 0, but also n = ±N_fp, ±2N_fp, ±3N_fp,… This fact implies that the toroidal mode number is no longer a good quantum number, but groups of n's known as mode families must be considered together. The coupling rules are that for a mode family centered on n, a sequence of other toroidal modes, denoted by n′ also need to be included, where n′ ± n = kN_fp, with k = 0, 1, 2,… As a result of this coupling, a sequence of n's will need to be included even for linear stability analysis of 3D systems (up to some truncation point); however, there will only be a finite set of mode families to consider since they repeat themselves beyond a certain point. For even values of N_fp, there are 1 + N_fp/2 mode families and for odd values of N_fp, there are (N_fp – 1)/2 + 1 mode families. In Figure 17, the toroidal mode numbers contained in the n = ±1 mode family, and the number of different mode families are plotted for devices with different numbers of field periods. Typically, the toroidal mode couplings are less dominant for devices with high field periods, due to the large jumps in n, than for low field period systems. In modeling lower field period configurations, the analysis also becomes more demanding as toroidal mode couplings become more closely spaced and more n's must be included.

Several examples of Alfvén modes in stellarators are given in Figure 18. Figures 18(a)–18(c) present the 3D coupled continuum and RSAE mode structure from an LHD case that had a weakly reversed shear region¹⁵³ in the ɨ profile near the axis; such a profile was created using off-axis neutral beam current drive. This was an n = 1 dominated mode that had very weak coupling to adjacent n's, as is expected for a high field period device (N_fp = 10). Figures 18(d)–18(f) present the 3D coupled continuum from the lower field period (N_fp = 4) TJ-II device¹⁵⁴ and an edge-localized HAE mode structure. For this mode and the other gaps indicated, more significant toroidal coupling was present due to the lower field periods and lower aspect ratio of TJ-II.

The stability analysis of energetic particle instabilities in 3D equilibria has progressed from semi-analytical models,¹⁵⁵ to perturbative approaches, hybrid particle-MHD, and gyrokinetic models. The first perturbative approach was the CAS3D-K model,¹⁵⁶ which utilized stable Alfvén gap eigenfunctions from the CAS3D MHD δW code¹⁵⁷ and solved the linearized fast, thermal ion and electron kinetic equations for growth and damping rates, keeping the mode structure and real frequency fixed. A next step in the development of perturbative models (i.e., fixed MHD mode structure and real frequency) was to include finite orbit width effects using δf methods that evaluate wave-particle energy transfer rates. By following a large number of fast ion orbits in the presence of the AE mode of interest, and accumulating averages of energy transfer rates, growth rates could be inferred. This approach has been incorporated in codes such as AE3D-K,¹⁵⁸ Venus,¹⁵⁹ and EUTERPE.¹⁶⁰ The hybrid approach is to couple a particle δf kinetic equation solution to an MHD model, and this is used in the MEGA¹⁶¹ code, which has been applied to LHD. Finally, several global gyrokinetic PIC codes are under development for 3D systems (GTC¹⁶² and EUTERPE¹⁶⁰).

The study of energetic particle driven instabilities in 3D tokamaks is at an earlier stage of development than for stellarators. High-resolution continuum studies¹⁶³ have demonstrated the fine scale TF ripple–produced gap structures that are present as high n AE modes couple with ripple with similar wavelengths. Experiments on NSTX¹⁶⁴ have recently studied both enhanced losses of energetic particle populations with 3D ELM control fields present and regimes where fast ion confinement is not strongly influenced, but AE stability may be suppressed by enhanced continuum damping from the toroidal mode couplings introduced by the 3D coils. Figure 19 displays the Alfvén continuum changes that have been predicted with and without 3D effects from the STELLGAP code¹⁴⁹ for the n = 2 Alfvén mode. Suppression of the n = 2 mode was observed on NSTX with activation of the ELM coils in a regime, where fast ion confinement was not strongly affected.¹⁶⁵ As indicated in Figure 19(b), the 3 field period 3D equilibrium couplings of the n = 2, 5, and 8 modes near the edge create a localized modification in the continuum structure, which could increase continuum damping in the frequency range (red line) where the n = 2 instability was observed.

In similarity with tokamaks, 3D configurations experience plasma transport from microturbulence and this generally dominates the transport physics. Confinement scaling databases have been accumulated for stellarators and scaling laws inferred^166,167 that are similar to tokamak gyro-Bohm scalings with an H-mode enhancement factor; a configuration scaling function, based on neoclassical effects, has also been identified.¹⁶⁶ Both linear¹⁶⁸ and nonlinear simulations¹⁶⁹ of core microturbulence candidate modes, such as the ITG (Ion Temperature Gradient) and TEM (Trapped Electron Mode) instabilities have been made for a variety of stellarators. Recently, the possibility of optimizing stellarators to reduce turbulent transport has been explored. A variety of proxy target functions have been investigated,¹⁷⁰ including linear growth rate estimates, quasilinear/mixing length estimates of the heat flux, and direct linear calculations of growth rates from the GENE model. Applications have been made to NCSX, W7-X, and HSX. A recent optimization, based on using the flux-surface version of GENE has resulted in a turbulence optimized stellarator configuration¹⁷¹ that was evolved from W7-X as a starting point. The differential evolution option¹⁷² of STELLOPT was used. The optimization history, demonstrating a significant improvement in the ITG target function, is shown in Figure 20(a) and the resulting improvement in the post-evaluated ion heat flux is shown in Figure 20(b) between the initial W7-X shape and the optimized MPX configuration. The possibility of optimizing stellarators for reduced levels of microturbulence transport is an exciting new direction for configuration design.

The tokamak edge pedestal and scrape-off layer are important, but challenging regions for 3D modeling. Axisymmetric models, such as the ELITE¹⁷³ model, have been developed that successfully predict characteristics of plasma gradients in the pedestal region, and the ELM instabilities they drive by balancing off the marginal stability boundaries for peeling and ballooning instabilities. Here, peeling instabilities are global external current gradient-driven instabilities, while ballooning instabilities are localized pressure gradient driven modes. The edge is also the region that is most likely to be modified by the 3D magnetic field perturbations from ELM coils. These will, at low levels, introduce corrugations in flux surfaces [recall Figure 7(a)], and at higher levels, produce islands and chaotic field lines. To the extent, that corrugations are resonant with the field line helicity, they can also facilitate unstable external kink modes. Work on this area of 3D physics is still at an early phase, due to the difficulties associated with resolving regions with high shear, many coupled modes/resonances, and the low field period characteristics of the 3D equilibria. Two recent approaches are the direct application¹⁷⁴ of global δW MHD stability analysis and the development of coupled high-n (ballooning) and intermediate-n (peeling) models that include 3D fields.¹⁷⁵ In Ref. 174, the CAS3D model was applied to low n (≤4) modes for several cases: a tokamak with a helical core; RMP coils in ASDEX upgrade; and ITER with TBMs. Both stabilizing and destabilizing effects of the ELM coils and TBM effects were seen. In Ref. 175, an analogue of the ELITE peeling-ballooning model was developed for 3D systems by using multiple-scale techniques and energy minimization methods.

MHD instability analysis for stellarators has predominantly focused on pressure-driven modes, due to the small levels of plasma current that characterize most stellarators. Exceptions to this are the hybrid systems, which augment the vacuum transform with that provided by limited amounts of plasma current; these systems must also take into account current driven modes.¹⁷⁶ The simplest criteria for pressure-driven modes are those obtained in the interchange limit (k_||/k_⊥ ≪ 1) and include the magnetic well condition (⁠ $V″<0 where V′=specific volume=∮dl/B$⁠) and the Mercier criterion.¹⁷⁷ These lead to conditions local to each flux surface that are readily evaluated from plasma equilibrium quantities. Ballooning stability models, such as COBRA,¹⁷⁸ have also been developed and provide efficient evaluations that can be used in optimizations. Stellarator optimizations have typically targeted at least the magnetic well, Mercier, and ballooning stability conditions. For global linearized MHD, models (CAS3D¹⁵⁷ and Terpsichore¹⁷⁹) based on the δW variational principle are used; these have been applied both to pressure and current driven instabilities. Finally, several initial value MHD models, such as NORM,¹⁸⁰ FAR3D,¹⁸¹ and M3D,¹⁷⁶ have been applied both linearly and nonlinearly. A topic of long-standing interest for MHD models is the “soft β-limit” characteristic of stellarators and the fact that experiments seem to routinely exceed stability limits predicted by theory. For example, LHD has achieved¹⁸² peak β's of 5.1% and W7-AS has achieved¹⁸³ ∼7%. This has recently been analyzed¹⁸⁰ for the case of LHD using the NORM nonlinear MHD code and a multi-scale approach that allows the pressure profile to slowly evolve in a self-organized manner uncoupled from the MHD fields. This results in many locally flattened regions in the pressure, but avoids disruptive behavior and achieves higher β's, albeit lower than those seen in the experiment. For the both experiments, some degree of magnetic well self-stabilization may also have contributed. Magnetic well stabilization and second stability access was a key feature of the ATF stellarator design, and evidence¹⁸⁴ of this was observed as heating power was increased. Second stability for ballooning stability in compact quasi-omnigenous stellarator configurations was reported, leading to hybrid systems with very high β limits¹⁸⁵—volume averaged ballooning β limits up to 23%. Second stable regimes were also found for QPS configurations.¹⁸⁶ This ability for stellarators to achieve high β's experimentally without strong disruptions and in optimized designs offers attractive opportunities for reactor designs, given the strong dependence of reactor cost on plasma β.

Symmetry-breaking and 3D physics effects exist to some degree in all toroidal plasma confinement devices and are increasingly used to achieve control and improved plasma performance objectives. Exploration of the 3D design space is one of the major remaining areas of opportunity for the optimization of fusion confinement devices. Specific directions that are currently of interest include: (1) the tokamak 3D edge—optimization for ELM suppression, divertor structures, control of scrape-off layer, and detachment; (2) improved RFP confinement and sustainment through spontaneous helical states; and (3) new directions in stellarator optimization—transport, microturbulence, MHD stability, and energetic particle physics. The increasing availability of computational resources and development of new theoretical methods is well aligned with the needs for improved modeling of 3D physics issues. The significant interest that is now directed toward 3D effects in tokamaks, stellarators, and RFPs also provides excellent opportunities for testing and validating new simulation models that have been designed with application to all of these systems taken into account. This tutorial paper has presented an overview of a range of 3D physics topics that apply to all toroidal confinement devices. Those interested in greater depth on specific topics are directed to the extended list of references provided and ongoing research that is actively reported on in these areas.

This material based on work was supported both by the U.S. Department of Energy, Office of Science, under Contract No. DE-AC05-00OR22725 with UT-Battelle, LLC and under the U.S. DOE SciDAC GSEP Center. This tutorial has benefitted from discussions with many colleagues, including: Jeffrey Harris, Andreas Wingen, Alessandro Bortolon, Steve Hirshman, Diego del-Castillo-Negrete, Jeremy Lore, John Canik, Yasushi Todo, Axel Könies, Lee Berry, John Koliner, John Sarff, Harold Weitzner, Allen Boozer, Allan Reiman, Nate Ferraro, Gerrit Kramer, Harry Mynick, Pavlos Xanthopoulos, Matt Landremann, Per Helander, Anthony Cooper, Kazuo Toi, Sam Lazerson, Mike Zarnstorff, and Kouji Shinohara.

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).