Hamiltonian description for magnetic field lines in fusion plasmas: A tutorial

Creating and confining hot plasmas are the foundation of fusion studies.^1,2 As the increase in the temperature helps to create the plasma, magnetic fields are able to confine it in suitable containers. The spatial structure of magnetic field lines is an important ingredient in many theoretical analyses of magnetically confined plasmas in toroidal devices like tokamaks, stellarators, reversed field pinches, etc.³ In tokamaks, the magnetic field responsible for the confinement results from the superposition of the toroidal field, generated by external coils wound around the entire torus, and the poloidal field, due to plasma current itself.⁴ However, the equilibrium magnetic field can be modified by plasma oscillations or by external coils used to control instabilities.⁵ 

An interesting situation is where the magnetic field is time-independent, as the case in MHD equilibrium configurations.¹ Starting from a symmetric plasma equilibrium configuration with an ignorable coordinate (e.g., the toroidal angle in tokamaks), the magnetic field line equations can be cast in the form of canonical equations, if the ignorable coordinate plays the role usually assigned to physical time in classical mechanics.⁶ Furthermore, as the magnetic field is divergence free, we can describe the field lines using a two dimensional area-preserving map, with respect to a surface of section of the torus at a fixed toroidal angle.⁷ The resulting phase space of the field lines is identical to a Hamiltonian phase space, indicating that the field lines act, at least locally, as trajectories.⁷ Hence, the dynamics described by the corresponding Hamiltonian represents not a true motion but instead a magnetostatic structure parameterized by the time-like coordinate.⁸ The main advantage of making this analogy is to use the powerful toolbox of Hamiltonian theory to investigate the magnetic field line structure, in particular, if nonsymmetrical perturbations are considered.⁹ 

The analogy between magnetic field lines and a Hamiltonian system has been first pointed out by Kruskal, in 1952.^7,10,11 Kruskal proposed and iterated an area-preserving map, similar to the standard map, in order to describe the magnetic field lines of stellarators.^7,11 This connection between field lines and Hamiltonian formalism was also recognized simultaneously but independently in the United States by Donald Kerst (the inventor of betatron)¹² and in Soviet Union.¹³ Nevertheless, an explicit and generalized Hamiltonian description was only proposed later, by Whiteman¹⁴ and Boozer.¹⁵ 

Even though magnetic field lines can be described by low-dimensional Hamiltonian systems, the numerical integration of the motion equations can be computationally costly.¹⁶ For this reason, explicit area-preserving maps, derived from the Poincaré map of the magnetic field line, are often used. These Hamiltonian maps inform us about the global and fine scale structure of the edge magnetic topology in toroidal systems. They serve as an important tool for studying the kinetic and fluid transport process as plasma turbulence and MHD stability.⁵ In addition, these maps permit long-time examination of individual trajectories for the statistical analysis of the field lines and the investigation of transport with a reasonable computational time.² 

The theory of magnetic field lines in confined plasma devices could not guarantee the regularity of the field lines.¹¹ An example is the study of magnetostatic perturbations produced by coils placed outside the plasma: the resulting magnetic field lines may present unexpected and complicated behaviors like periodic, quasi-periodic, and even chaotic orbits.^7,17 The latter, in particular, represent a local destruction of the magnetic surfaces that confine the plasma.^18,19

The main motion of the plasma particles is along the field lines while slowly drifting across the equilibrium fields due to the lines curvature, the particle rotation around the lines, and electric drift.^20–22 Thus, as the fastest motion is along the field lines, the particle escape to the wall can be predicted by the field line configuration.²³ Observing magnetic field lines in a surface section of a tokamak, they can be closed lines within the trace of a toroidal magnetic surface or they can fill a two dimensional domain. For the first case, the field line is regular, while the second case indicate chaos.¹¹ Chaotic behavior is related to the topology of the magnetic field lines and the dynamics of the particles gyrating along these lines, as well as the turbulent transport, ray dynamics, and radio frequency heating.^11,15 Furthermore, the non-uniform particle transport at the tokamak plasma edge can be roughly estimated from the field lines escaping to the wall.^23,24 Broader reviews between Hamiltonian chaos and fusion plasmas can be found in Refs. 5, 11, and 16.

From a classic mechanics framework, a general description in curvilinear coordinate system was proposed by Whiteman¹⁴ and extended later by Boozer,¹⁵ Cary and Littlejohn,⁸ and Elsässer.²⁵ The general formalism described by Whiteman has been applied to a variety of coordinate systems: cylindrical,²⁶ helical,²⁷ spherical,²⁸ and pseudotoroidal.²⁹ Various applications of the magnetic field line Hamiltonian have been made by Freis et al.³⁰ and Hamzeh³¹ for a toroidal machine called Levitron and by Lichtenberg^32,33 in an investigation of the m = 1 island on sawtooth oscillations in tokamaks.

In addition to its applications in fusion plasmas, the Hamiltonian description of magnetic field lines provides a nice non-mechanical example of the usefulness of the Hamiltonian formalism to intermediate and advanced students. Moreover, the magnetic field line problem has the unique feature that the corresponding phase space actually coincides with the configuration space, which facilitates the visualization of complex dynamical concepts like Kolmogorov-Arnold-Moser (KAM) tori, homoclinic tangles, and so on. On the other hand, the basic material on the Hamiltonian description of magnetic field lines is often available only in publications targeted to the plasma physics experts, which creates an additional difficulty for an interested novice reader.

In order to overcome the latter problem, we wrote this tutorial as an aid to students and researchers interested to master the basic ideas of the Hamiltonian description for the magnetic field lines. Moreover, we present some representative applications of this formulation so as to illustrate its usefulness in plasma physics problems. We emphasize that this paper is not a review of this subject. While we focused on fusion plasmas, the methods can also be used in plasmas of astrophysical and geophysical interest, provided we have situations of MHD equilibrium with adequate stability properties.

This paper is organized as follows: in Sec. II, we show the derivation of a variational principle for magnetic field lines and the role played by gauge invariance. The Hamiltonian description, in general, curvilinear coordinates is presented in Sec. III. In Sec. IV, we describe in some detail an application of the description in cylindrical coordinates to a large aspect-ratio tokamak with an ergodic magnetic limiter (EML), using canonical perturbation theory to derive an analytical formula for the width of magnetic islands, which is an expression of practical interest for stability and transport theoretical studies of tokamak plasmas.¹ In Sec. V, we present an application of the general formulation for magnetic flux coordinates, which are widely in numerical codes for computer simulation of plasmas,¹⁵ displaying an application to a magnetic field line map (tokamap) proposed by Balescu and co-workers.³⁴ The last section is devoted to our Conclusions.

When using the above formulas, one must have in mind that quite often the contravariant components of the magnetic field have not the same dimensions as the field itself due to the metric coefficients. In the forthcoming section, our aim is to present representative applications of this description, using different coordinates systems, as the cylindrical and magnetic flux coordinates.

Let us consider a toroidal plasma with major radius R₀ and minor radius a. In the local (or pseudotoroidal) system of coordinates $( r , θ , ϕ )$⁠, θ and $ϕ$ are the poloidal and toroidal angles, respectively, and r is the radial distance to the magnetic axis, which is a circle of radius R₀ centered at the torus major axis.

The torus aspect ratio is $ε = R 0 / a$⁠. In the large aspect ratio approximation (⁠ $R 0 ≫ a$⁠), we can neglect the toroidal curvature and consider the torus as a periodic cylinder of radius a and length $2 π R 0$⁠. In this case, it is possible to use cylindrical coordinates $( r , θ , z )$⁠, where $z = R 0 ϕ$ is the rectified toroidal circumference. Due to the periodicity, we identify all points for which z is an integer multiple of $2 π R 0$⁠.

Divertors are devices used in tokamak experiments with the purpose to displace the interactions between the plasma particles and the tokamak wall, thereby avoiding direct contact between them and improving plasma confinement.^4,38 Initially, the divertors were designed to act directly over the toroidal or the poloidal field of the plasma and they required additional coil currents of the magnitude of plasma currents or even larger.³⁹ This led to experimental limitations and technological problems for the tokamak operation.³⁹ As an alternative, Karger and Lackner proposed the helical divertor, which requires smaller currents and possesses helical symmetry, which generates a magnetic field that resonates with the field at a surface in the plasma boundary, which is diverted.³⁹ 

The resonance created by the divertor can also lead to a chaotic motion in the plasma edge, a process called “ergodization.”^40,41 The term “ergodic,” however, has been later replaced by “chaotic,” which is a more adequate description of the area-filling orbit created when the invariant manifolds stemming from unstable periodic orbits intercept in a complicated way forming the homoclinic tangle.⁴² The chaotic field lines increase the diffusion coefficient at the boundary of the plasma, reducing the plasma contamination,⁴⁰ controlling the MHD oscillations,⁴¹ reducing thermal flux density,⁴³ and controlling plasma–wall interaction.⁴⁴ For more information about the theoretical and experimental development of helical divertors, confer Refs. 39, 40, 43, 45, and 46.

The external magnetic fields generated by the helical divertor create magnetic islands that can overlap and, consequently, form a stochastic layer in the plasma edge. Therefore, such a divertor is also called ergodic magnetic limiter (EML).⁴³ The EML is a filamentary current ring with length $ℓ$⁠, wound around the torus (with radius a). There are two types of current segments in a EML (Fig. 1): straight segments parallel to the magnetic axis and curved segments along the poloidal direction. There are m pairs of segments, such that two adjacent segments carry a current I_L in opposite directions,⁴⁷ which produces a resonant helical field.⁴¹ 

In order to verify whether or not these conditions are satisfied, we will take parameters from the TCABR (Tokamak à Chauffage Alfvén Brésilien at Instituto de Física, Universidade de São Paulo, Brazil) presented in Ref. 49, where $R 0 = 0.61$ m and a = 0.18 m. The plasma current is about $I p = 0.1 MA$⁠, and the toroidal field at the magnetic axis is $B 0 = 1.1 T$⁠. The safety factor is $q a = 2.95 ≈ 3$ at the plasma edge, and $q 0 = 1.5$ at the magnetic axis. An EML has been installed in TCABR with m = 3 pairs of wires with length $ℓ = 0.1 m$⁠, carrying a current about $I L = 2500 A$⁠.⁵⁰ These values imply that $ε = 0.025$ and $ξ = 0.163$ are small enough to justify treating the EML field as a Hamiltonian perturbation upon the equilibrium magnetic surfaces. Ergodic limiters have also been used in other tokamaks as in Textor,² Tore-Supra,⁴⁶ and Text.⁴³ 

The pendulum Hamiltonian describing the resonance is integrable because we have averaged out the non-resonant terms in (55). If we include these terms again, the system will become quasi-integrable and the pendulum separatrices will no longer join smoothly, but rather will present an infinite number of homoclinic and heteroclinic points. The dynamics near such points is known to be chaotic, and, as a result, instead of separatrices, the islands will have a thin stochastic layer of chaotic motion.⁵¹ As long as the intensity of perturbation is small enough, these locally chaotic layers do not connect themselves, preventing large-scale chaotic transport of field lines. If the limiter current, however, is larger than a critical value, the locally chaotic layers merge together forming a globally chaotic region, and allowing large-scale chaotic transport. Pendulum approximation have been used to estimate the islands' width and to apply Chirikov criterion to find the critical perturbation amplitudes to create a chaotic layer in the plasma.^52,53 We remark that Escande and Momo have derived a similar formula for the island half-width without using Fourier components like $A m ( J * )$⁠, but rather using the magnetic flux through a ribbon, the edges of which are lines passing by the elliptic and hyperbolic points in a given magnetic island.³⁶ 

The EML Hamiltonian (55) exhibits an explicit dependency on the ignorable variable $ϕ$⁠. For non-autonomous systems, the trajectories belong to an extended phase space, where the time is treated as a coordinate.⁵¹ For the quasi-integrable Hamiltonian system in (58), the solutions are in a three dimensional phase space and the flow is parameterized by the variable $ϕ$⁠. In this way, a state is determined by three variables: J, θ, and $ϕ$⁠.

Instead of studying the solution of the system in a three dimensional geometric space, we can reduce the dimensionality of the problem by the construction of a Poincaré surface. The Poincaré map is formed by the intersection of the solutions in the surface defined at a constant value of $ϕ$⁠. In this way, we have the values of J and θ for when the magnetic field lines cross the surface, i.e., for each complete turn in the toroidal direction.

Integrating Eqs. (73), for initial conditions [ $I ( ϕ = 0 ) , θ ( ϕ = 0 )$], using a symplectic Euler⁵⁴ method and defining the Poincaré section at $ϕ = 0 ( mod 2 π )$⁠, we construct the Poincaré sections in Figs. 3(a) and 3(b), for a scenario of small (⁠ $ε = 0.025$⁠) and large (⁠ $ε = 0.15$⁠) limiter currents, respectively.

For the Poincaré section in Fig. 3(a), the limiter current corresponds to 2.5% of the total plasma current I_L. In this scenario, we observe regular solutions in most part of the space, represented by the rotational circles and the three islands (oscillatory circles). We observe a thin chaotic layer acting as a “separatrix” of the island. Since ε is too small, the width of this chaotic layer is so tiny it resembles a separatrix curve, as the one presented in the pendulum phase space in Fig. 2. If the limiter current is increased until it corresponds to 15% of the total plasma current [six times the current in Fig. 3(a)], we have the Poincaré section shown in Fig. 3(b). The second Poincaré section also shows three islands, and the separatrix between them is replaced by a thick chaotic layer. The chaotic behavior emerges with the increase in the perturbation parameter ε, i.e., with the increase in the current in the ergodic limiter.

Finally, we would like to estimate until what value of ε Eq. (70) is a good approximation for the half-width of the islands in the phase space. We numerically solve the system (73), construct the Poincaré section, and compare the half-width of the islands for each ε with the value predicted by Eq. (70). The results are presented in Fig. 4, and we observe that the pendulum approximation is valid for the half-width of the islands in the phase space until $ε ≈ 0.1$⁠. For higher values of ε, the value of I_max does not increase at the same rate proportional to $ε 1 / 2$⁠. For $ε > 0.1$⁠, we observe small increases and decreases in I_max, which represents the enlargement and the following destruction of the island.

In order to put these equations into the form (116) and (117), it is necessary to solve them first for $ψ n + 1$⁠. Although this can be done analytically in some cases, it is always possible to use root-finding methods to do so numerically. Another important point, emphasized by Balescu and co-authors, is that the field line map must have two properties: (i) if $ψ 0 > 0$⁠, then $ψ n > 0$⁠, for all values of n; (ii) if $ψ 0 = 0$⁠, then $ψ n ≥ 0$ for all n.³⁴ The former property comes from the definition of the coordinate ψ, which must be a definite positive number, whereas the latter is related to the fact that ψ = 0 stands for the magnetic axis (which is a degenerate magnetic surface).

We observe that the tokamap exhibits mainly periodic and quasi-periodic solutions for K = 1, indicated by the existence of only islands and rotational circles in the Poincaré section of Fig. 5(a). Like in Sec. IV C, if the perturbation strength is too small, the size of a chaotic layer in the neighborhood of an island is so tiny that it can be revealed only by magnifications of the Poincaré section. From the winding number profile calculated in $θ = 0.5$ [Fig. 5(b)], we observe a defined ω_n for almost every value of ψ, and the profile monotonically decreases. The possible exceptions consist of tiny intervals for which the orbit is chaotic.

In Fig. 5(b), we choose three different plateaus and highlight them with the colored rectangles. The correspondent island of the plateau is shown with the same color in the Poincaré section of Fig. 5(a). From the values of ω_n, we identified a direct relation with the period of the islands. The winding number is related to the period τ of the island by $ω n = 1 / τ$⁠. For example, the red island of period 1 presents a winding number equal to $ω n = 1 / 1$⁠. The green and pink islands present winding numbers equal to $ω n = 0.5 = 1 / 2$ and $ω n = 0.33 … = 1 / 3$⁠, respectively. These periodic islands are on rational tori, since we can write their frequencies as a ratio between two integer numbers.

Increasing the perturbation parameter to K = 3.5, we have the Poincaré section and the winding number profile showed in Fig. 6, where we observe that some regular solutions are replaced by stochastic layers, represented by the chaotic seas around the green and the orange lines. The chaotic behavior is restricted, and the chaotic regions are not connected. If we increase K, these chaotic regions will eventually enlarge and merge together into a single area-filling chaotic orbit. Following the same methodology as for Fig. 5, we computed the winding number profile, at $θ = 0.5$⁠, and highlight the plateaus of constant ω_n values. Again, we observed the directed relation between the winding number value and the period τ of the islands. The red, green, and pink islands of Fig. 5 are also seen here, with $ω n = 1 / 1 = 1 , ω n = 1 / 2 = 0.5$⁠, and $ω n = 1 / 3 = 0.33 …$⁠, respectively. We also highlight two other chains of islands, orange and blue, with $ω n = 2 / 3 = 0.66 …$ and $ω n = 2 / 5 = 0.4$⁠, respectively. For these last two islands, we have that the trajectory always “jumps” one island during the time evolution, i.e., if we choose an initial condition in the blue island close to $θ = 0.5$ (the third island counting from left to right), the next point will be in the fifth island, the second iteration will be in the second island, and so on. The chaotic regions are represented by the “gaps” in the winding number profile, since the limit (135) does not converge.

The Tokamap has been used to interpret the particle escape to the wall in Textor Tokamak. In particular, the theoretically obtained fractal distribution of field lines at the plasma edge is similar to the one measured in this tokamak.^2,23