Proceedings of ITC/ISHW2007

Theory of Stellarators and Tokamaks in Three Dimensions
Paul R. Garabedian
Courant Institute of Mathematical Sciences, New York University
New York, New York 10012 USA
The NSTAB computer code applies the MHD variational principle to calculate equilibrium and stability of toroidal plasmas in three dimensions. Differential equations are solvedin a conservation form that describes force balance correctly across islands that are treated as
discontinuities. The method has been applied to stellarators, including LHD, and tokamak configurations such as DIII-D and ITER. Sometimes the solution of the equations turns out not
to be unique, and there may exist bifurcated equilibria that are nonlinearly stable when theory
predicts linear instability. With appropriate pressure profiles, the LHD plasma is found to bestable at β = 0.045, which correlates with recent experiments. Hence reactor values of pressure
can be confined stably in a stellarator configuration with robust flux surfaces. A similar analysis shows that tokamak configurations like ITER do not remain axially symmetric at finite β,for they develop helical islands, which may manifest themselves experimentally as neoclassicaltearing modes (NTMs) and edge localized modes (ELMs). These results motivate a continuingsearch for improved stellarator configurations which (a) employ relatively simple coils to generate good flux surfaces even at low aspect ratios typical of a tokamak; (b) retain the favorable
high β equilibrium and stability characteristics demonstrated in LHD; and (c) promise reduced
transport by virtue of quasiaxial symmetry.

more stable configuration with the external
field generated by a system of complicatedcoils in three dimensions. Modern methodsof computational science enable physicists
to perform parameter studies necessary todesign practical configurations for experiments to establish the feasibility of this concept of magnetic fusion.
Three-dimensional computer codes havebeen written to study equilibrium, stability
and transport of plasmas in toroidal geometry for magnetic fusion. Advanced numerical methods have been employed to solve
systems of differential equations that model
these configurations. Despite the simplicity of some of the models, good correlation
has been obtained with observations fromexperiments. We shall discuss comparisonswith the Large Helical Device built and operated in Japan and the International Ther-

1. Computational Science of MagneticFusion
For a long time it has been recognizedthat fusion of hydrogen to form helium iswhat powers the sun and the stars. Now
a community of scientists believes that asimilar process releasing energetic neutrons
can be exploited to construct nuclear reactors that might eliminate global warming by fossil fuels. A very hot plasma of
deuterium and tritium is to be confined in
toroidal geometry by a strong magnetic fieldthat separates the ionized gas from material
walls. The most successful fusion experiments have been axially symmetric tokamaks with net toroidal current induced in
the plasma itself to produce a poloidal fieldrequired for equilibrium. The stellarator is a1

Proceedings of ITC/ISHW2007

monuclear Experimental Reactor planned forconstruction in France by a consortium ofcountries. Afterwards we shall propose analternate device, the Modular Helias-like Heliac 2, a stellarator that seems to performbetter according to the predictions of a mathematical theory which we have developed.
2. EquationsThe Maxwell stress tensor enables one
to put the differential equations describing
force balance in the conservation form
hReferences
1. R. Alexander and P. Garabedian, Choiceof coils for a fusion reactor, Proc. Natl.
Acad. Sci. USA 104 (2007) 1225012252.
2. F. Bauer, O. Betancourt and P. Garabedian, Magnetohydrodynamic Equilibrium
and Stability of Stellarators, SpringerVerlag, New York, 1984.

i∇ · BB − (B 2 /2 + p)I = 0, ∇ · B = 0.

3. F. Bauer, P. Garabedian, D. Korn andA. Jameson, Supercritical Wing SecComparable finite difference equations, summed
tions II, Springer-Verlag, New York,over a test volume, telescope down to a cor1975.
rect statement of force balance over the boundary. This way we capture discontinuities in
4. D. Brennan, R. La Haye, A. Turnsolutions
bull, M. Chu, T. Jensen, L. Lao, S.
Kruger and D. Schnack, A mechanism
for tearing onset near ideal stability
boundaries, Phys. Plasmas 10 (2003)
1643-1652.

B = ∇s × ∇θ = ∇φ + ζ∇sof the MHD variational principle
δZ Z Z(B 2 /2 − p(s)) dV = 05. P. Garabedian and L. Ku, Reactors
with stellarator stability and tokamaktransport, Fusion Sci. Technol. 47(2005) 400-405.

for a plasma with separatrix defined by theFourier seriesr + iz = eiu

X∆mn e−imu+inv .

6. P. Garabedian and M. Taylor, Tokamak transport driven by quasineutrality and helical asymmetry, Nucl. Fusion 32 (1992) 265-270.

Islands occur where the spectrum Bmn defined byX
1=Bmn (s) cos [mθ − (n − ιm)φ]B2

7. S. Hirshman and O. Betancourt, Preconditioned descent algorithm for rapidcalculations of magnetohydrodynamic
equilibria, J. Comp. Phys. 96 (1991)
99-109.

activates small denominators of the parallel
current
J · B X mBmn (s)=cos [mθ − (n − ιm)φ] .
p′ B 2n − ιm
As an example, consider the RFP model
problem(Ψ2x )x = ηΨxxx,

8. A. Komori and H. Yamada, Overview
of the Large Helical Device, Plasma
Phys. Control. Fusion 42 (2000) 11651177.

Ψ = 1 − |x| .

2

Proceedings of ITC/ISHW2007

9. J. Nuehrenberg and R. Zille, Quasihelically symmetric toroidal stellarators, Physics Letters 129A (1988) 113116.
10. M. Taylor, A high performance spectral code for nonlinear MHD stability,
J. Comp. Phys. 110 (1994) 407-418.

FiguresFigure 2: Cycles of a calculation of the energy confinement time τE in milliseconds for
an NBI shot of the LHD experiment using aquasineutrality algorithm to adjust the electric potential Φ. Oscillations of Φ along
the magnetic lines model anomalous transport, so there is good agreement with theobserved value. Results from the LHD experiment have served to validate the numerical simulations of equilibrium, stability andtransport that are provided by runs of theNSTAB and TRAN computer codes. TheTRAN code models thermal transport byperforming a random walk among complicated orbits of the ions or electrons, which
are found from Runge-Kutta solutions of asystem of ordinary differential equations forguiding centers. The results depend primarily on the magnetic spectrum of the plasma,
which is obtained from runs of the NSTAB
code.

Figure 1: Poincaré sections of the flux
surfaces of a bifurcated, nonlinearly stable
LHD equilibrium at β = 0.048. The ripple in the flux surfaces on the right suggeststhat ballooning modes appear in the solution, but a reliable prediction of the β limit
cannot be made without more information
about the experiment. The existence of several solutions of the ideal MHD equilibrium
equations is considered to be evidence of linear instability.

3

Proceedings of ITC/ISHW2007

Figure 4: Poincaré section of the flux surfaces of a bifurcated ITER equilibrium at
β = 0.027 with net current bringing the rotational transform into the interval 0.93 >ι > 0.37. There are helical islands inthis three-dimensional solution of an axially symmetric MHD problem, which may
model ELMs. A forcing term was used earlyin the run to trigger a mode that appears
as a discontinuity in the three-dimensional
solution. This double precision calculation
converges to the level of round off error,
showing that the discrete problem has beensolved. Islands are captured numerically byfinite difference equations in a conservation
form that works despite the nested surface
hypothesis present in our mathematical formulation of the problem. The results areplausible on the long time scale of a magnetic fusion reactor.

Figure 3: Four Poincaré sections of the magnetic surfaces over one out of three artificial field periods of a bifurcated, nonlinearly stable ITER equilibrium at β = 0.03
with 0.9 > ι > 0.4 . The ripple in the flux
surfaces suggests that there may be NTMs
and ELMs in this three-dimensional solution of an axially symmetric problem. Extensive studies of numerical examples produce ample evidence that bifurcated equilibria can be expected to appear in mosttokamak problems, so provision should bemade for that. Convergent runs of theNSTAB code can capture islands whosewidths are smaller than the radial mesh sizeand compute bifurcated equilibria in examples where there is marginal linear instability, but global nonlinear stability, of theideal MHD model.

4Proceedings of ITC/ISHW2007

Figure 5: Cycles of a calculation of the energy confinement time τE in milliseconds for
the ITER tokamak using a quasineutrality
algorithm to adjust the electric potential Φ.
The three-dimensional effect of ripple associated with a system of just twelve toroidal
coils has been introduced to drive the radial electric field and cause the plasma tospin. It is not clear that two-dimensionalmodels describe transport in tokamaks adequately. The energy confinement time is calculated from an empirical relationship withthe particle confinement time. Without
three-dimensional terms we have not beenable to reconcile discrepancies between theion and the electron confinement times thatare computed by the Monte Carlo method
in tokamaks.

Figure 6: In magnetic fusion, hot deuterium
and tritium ions are combined to form helium and release neutrons intended to provide a commercial source of energy. The
color map of the hydrogen plasma shownin the figure displays a symmetry property
that enhances confinement. Twelve only
moderately twisted coils generate a magnetic field designed to keep the plasma instable equilibrium separated from material
walls. The side view of this compact stellarator shows that there is ample space between the coils for NBI heating. Our hope isthat a judiciously optimized QAS stellarator may overcome the troubles with poortransport and low ion temperature in conventional stellarators, and with poor stability and ELMs crashes in tokamaks.

5Proceedings of ITC/ISHW2007

Figure 7: Four out of twelve modular coils of
the MHH2 stellarator in a vacuum magnetic
field given by the Biot-Savart law. The coils
at the sides of the diagram are located atcorners over the full torus, so the distances
between all the coils can be estimated fromthe figure. Judicious filtering of the Fourier
series used to calculate filaments specifying
the geometry of the configuration defines
shapes that are not excessively twisted. Parameters have been adjusted to provide ample space around each coil, and the aspectratio of the plasma is 2.5. The spacing isadequate for superconducting coils in a reactor with major radius 7 m or 8 m, butthere is difficulty fitting the coils together
inside a small experiment at high magnetic
field. Through rigid motions, four copies of
the quadrant of coils shown in the plot canbe combined to give an accurate picture overthe full torus.

Figure 8: Cross section of a line tracing
calculation for a QAS stellarator displaying curves that define the control surface
for the coils and the shape of the plasma,
together with magnetic lines computed at
β = 0. The flux surfaces and islands outside
the separatrix show that the magnetic field
is well organized for a divertor. This works
because smooth coils were found to provide
the external field confining the plasma. It
is essential to choose the coils both so thatthe stellarator can be constructed without
too much difficulty and so that the flux surfaces remain robust when realistic changes
are made in physical parameters. The formulas we apply to represent the plasma surface and the shape of the coils are motivatedby a knowledge of conformal mapping andRunge’s theorem in the theory of analytic
functions of a complex variable.

6