Higher order theory of quasi-isodynamicity near the magnetic axis of stellarators

A stellarator^1–3 is an attractive magnetic confinement concept to achieve controlled thermonuclear fusion. Its flexibility and generality make it unique in its ability to confine plasmas avoiding common shortcomings of other confinement designs.^2,4 However, to achieve good confinement, the field must be optimized.⁵ Omnigeneous stellarators^3,6–10 are one such class of optimized stellarators, one that minimize the rapid loss of particles. Quasisymmetric stellarators^11–13 are a particular subclass.

In omnigeneous stellarators, in the absence of collisions, particles do not drift away from the device, on average, which avoids the fast loss of particles. This property is achieved only when the magnitude of the magnetic field satisfies certain conditions.¹⁰ Even with these constraints provided, the class of approximately omnigeneous stellarators has a significant freedom compared to the more restrictive case of a tokamak.^14,15 This makes the design of stellarators versatile but necessarily complicates their systematic study.

A possible way to approach the study of omnigeneous stellarators is to consider the behavior of these configurations near their magnetic axes. This simplifying near-axis expansion approach^16–21 has been a recurrent analytic tool in the study of stellarators. In this paper, we use it to present the consequences of omnigenity on the magnetic field magnitude near the magnetic axis. We specialize to the concept of quasi-isodynamicity, continuing the work initiated in Ref. 22.

The basic definition of omnigenity we build upon is presented in Sec. II. Deviating from Ref. 22, we use an alternative, physically transparent definition based on particle drift motion, which facilitates the work that follows. In Sec. III, we present the perturbative approach to assess quasi-isodynamicity near the axis to an order higher than what had been previously worked out. There, the property of quasi-isodynamicity is reduced to symmetry properties on the asymptotic form of $| B |$⁠. Section IV treats the issues of pseudosymmetry and the possibility of topological defects in $| B |$⁠. Section V closes with a brief summary and conclusions.

We define omnigenity as the condition on the magnetic field such that the bounce average radial drift of trapped particles vanishes. This notion goes back to the seminal work by Hall and McNamara.⁸ Here, we present one particularly physical form of the concept.

Consider a magnetic field with nested toroidal flux surfaces labeled by the toroidal flux $2 π ψ$⁠. On every surface, the magnetic field satisfies $B · ∇ ψ = 0$ and, thus, can be conveniently written in the Clebsch form²³ as $B = ∇ ψ × ∇ α$⁠, where α is the field line label $α = θ − ι φ$⁠. Here, $φ$ and θ are the toroidal and poloidal Boozer angles,²⁴ respectively, and ι is the rotational transform.³ The use of Boozer coordinates is possible as we consider B to be in equilibrium, satisfying $j × B = ∇ p$⁠, where $j = ∇ × B$ is the current density and p is the scalar pressure. It will be convenient to think of the magnitude of the magnetic field, B, as a separate function of the two angular coordinates at each magnetic flux surface.

With this setup, we are interested in studying the motion of trapped particles. Charged particles, to leading order in guiding center theory,²⁵ move along magnetic field lines as they gyrate about them. They do so by conserving energy, W, as well as the first adiabatic invariant, the magnetic moment μ. This conservative motion is what leads to the notion of trapped particles (also known as bouncing particles). From the definition of $μ = W ⊥ / B$⁠, where $W ⊥$ is the kinetic energy of the particle in the direction normal to the field, particles with a particular value of W and μ are excluded from regions for which $B > W / μ$⁠. As particles approach such regions, they are forced back. The result is particles that bounce back and forth along field lines between points at which $B = W / μ ≐ 1 / λ$⁠. Those points are called bouncing points and come in pairs dictated by $| B |$⁠. Formally, these are related through a bouncing function $η ( φ , θ , ψ )$ (see Fig. 1). We define this function to be the toroidal angle value of a bouncing point paired to one at $φ$ along a field line α, that is, the point with the same value of B on the other side of the well. This should not be confused with η in the previous work by Cary and Shasharina⁷ and Plunk et al.²² The parameter λ can be interpreted as defining different classes of bouncing particles.

We now study the property of omnigenity close to the magnetic axis of a toroidal magnetic configuration. This requires an asymptotic description of the magnetic field and an appropriate handling of the omnigenity property in Eq. (7). In that form, unlike in previous approaches,²² breaking omnigeneity has a direct interpretation in terms of particle drifts [see Eq. (4)], which also makes the connection to other important measures such as $ϵ eff$ straightforward,^28,29 since the factor given in Eq. (6) appears explicitly.

The expression in Eq. (8) presents two key features of the problem. The first is the coupling of the harmonics of the field and the powers of ϵ, necessary to enforce analyticity at the magnetic axis. In the second place, the magnetic field on axis is generally a function of the toroidal angle. We will assume B₀ to have such a non-trivial $φ$ dependence, which amounts to specializing to quasi-isodynamic (QI) configurations within the class of omnigeneous fields. These correspond to optimized omnigenous configurations with poloidally closed $| B |$ contours. This choice necessarily excludes the two other subclasses of omnigenous fields, those with toroidally and helically closed contours, as they can only be achieved near the magnetic axis if quasisymmetry is also satisfied.²² As in the treatment leading to Eq. (7), for simplicity, we shall assume there to be a single magnetic well, with its minimum located at $φ = 0$⁠.

Once we have the magnetic field magnitude in the form of Eq. (8), we are in a position to evaluate the radial drift measure Y. We must (i) first express Eq. (5) in Boozer coordinates and then (ii) consider its expansion in ϵ.

The second equation, Eq. (15b), describes the change in the bounce map as a result of the first order perturbation. If the change in B to the right and the left of the well is not the same, then the map function changes. This is generally the case, even when stellarator symmetry is invoked.

The conditions for omnigenity are obtained by bringing the expansion of the drift Y and the bounce map together.

These are precisely the same symmetry conditions obtained in Ref. 22 using the Cary–Shasharina construction⁷ (see also  Appendix C). The amplitude of the magnetic field correction B₁ must be, in a sense, odd about the minimum of B₀, while the function ν should be even.

Following the parity condition, d = 0 at the minimum of B₀. The same argument follows for the global maximum of B₀. In the scenario of multiple wells, the QI condition in Eqs. (21) does not appear to require the vanishing of d at maxima, which are not global. We must, however, be careful, as the asymptotic description fails in the neighborhood of the turning points. To enforce omnigeneity in that region, we must enforce pseudosymmetry to first order: that is, limit the radial drift of barely and deeply trapped particles to be $O ( ϵ α )$⁠, with $α > 1$⁠. As shown in  Appendix A, this requires B₁ to vanish at all turning points of B₀. At the minima, enforcing Eqs. (21) is enough. This result is consistent with Ref. 22, and it emphasizes the physical origin of the requirement, unlike in the cited work.

It is noteworthy that, as recognized in Ref. 22, it is impossible to satisfy the QI conditions, Eqs. (21), exactly. Doing so necessarily makes $| B |$ non-periodic, following the secular behavior of the field line label α (which at fixed θ gives $α → α − 2 π ι$ after one toroidal turn). In general, periodicity requires $ν ( 2 π ) = ν ( 0 ) + 2 π ( ι − N )$ for some integer N (related to the self-linking number of the axis^29,32), a difference that must be zero for omnigeneous fields, Eq. (21b). Thus, omnigeneity must be broken at the (global) top of the wells. Although this may appear to negate any higher order study of omnigeneity, there are practical strategies to address this. For instance, one is free to define intervals in $φ$ around minima of B₀ where Eq. (21b) is, if desired, satisfied exactly, and thus, all trapped particles residing in this interval behave in an omnigenous way to that accuracy. In addition, the pseudosymmetry requirement reduces the losses from breaking omnigeneity near the top, as it minimizes the magnitude of the radial drifts there. On a practical level, it is worth noting that approximately QI solutions found by numerical optimization can deviate significantly from first order solutions, and it is certainly useful that such deviations manage to preserve omnigenity. We, thus, continue with the procedure to higher order, assuming Eqs. (21) to be satisfied.

To cast these conditions in a form that is closer to the symmetry conditions at first order, we may define functions $B ̃ 2 s$ and $B ̃ 2 c$ so that

That is, omnigenity allows the second harmonics of the magnetic field to have a freely chosen part, so long as the choice satisfies the same symmetry as d at first order. We have defined $B ̃ 2 s$ and $B ̃ 2 c$ in this particular form as it separates $B 2 s$ and $B 2 c$ naturally into two components with “opposite” symmetry; the explicit, lower order terms in Eqs. (27) have the symmetry $f ( φ ) = − η 0 ′ ( φ ) f [ η 0 ( φ ) ]$⁠. In other words, $B 2 s$ and $B 2 c$ each have a fixed even-parity part necessary to outweight the non-omnigeneous influence of B₁ at second order, as well as a free odd-parity part.

We shall often be interested in learning about the reduced set of configurations that possess stellarator symmetry. Stellarator symmetry implies the invariance of $| B |$ under a set of discrete, parity-like maps. Here, we need only consider the symmetry around the minimum of the magnetic well, $( θ , φ ) → ( − θ , − φ )$⁠, which we take to be the symmetry point. This simplifies the QI conditions.

Getting to this point in the second order we were not very cognizant of the behavior of the expansion near points of $B 0 ′ = 0$⁠. We saw that to leading order, the behavior at these points required B₁ to vanish wherever $B 0 ′$ did. However, we gained nothing concerning the behavior of B in the neighborhood of these points. We know, though, that this is especially important for expressions like that of the perturbed map, η₁, Eq. (23), as there could be a region of asymptotic breakdown around $B 0 ′ = 0$ where η₁ is not well behaved. Misbehavior in this region could be seen to be a consequence of a change in $| B |$ contours, and potentially, the change of their topology. Whether this occurs or not depends on the behavior of B₀ and B₁ in the neighborhood of the turning points. The appearance of such “defects” (see Fig. 2) is in fact observed in many QI design efforts.^30,33,34

The magnetic field in the neighborhood of these points may be modeled as $B 1 ∼ φ v$ and $B 0 ′ ∼ φ u − 1$⁠. We refer to the indices v and u as the order of the zeros of B₁ and B₀, respectively. When the perturbation B₁ is flat enough, $v ≥ u$⁠, then the extrema of $| B |$ remain straight in the $( θ , φ )$ plane, and the topology of the $| B |$ contours is preserved (see Fig. 2). In this convenient case, the asymptotics in the neighborhood of the turning point are correct. However, this choice of order of the zeroes is generally not a necessary condition for quasi-isodynamicity. In particular, this choice retains the straightness of $B min$⁠, which we know is not implied by omnigenity. In fact, only the global maximum must be straight in QI configurations, but this conclusion only arises when periodicity is also considered.^7,10 To assess the behavior near the turning points and assess the physical requirement on the order of the zeroes, we must bring the notion of pseudosymmetry onto the scene.

For a magnetic field to be truly omnigeneous to the second order, the net drift of particles must be negligible to this order. Requiring the radial drift of particles to be of that order (or larger) at turning points of $| B |$ (along field lines) provides several conditions on the field, and the details of which may be found in  Appendix A. To be able to do this consistently, the formulation of the problem in terms of the geometric quantity Y is essential. The behavior near turning points depends on the values of the indices v and u. These may be organized in four different categories (see Fig. 4) according to the presence of topological defects and whether pseudosymmetry (zero radial drift at turning points) is achieved to order $ϵ 2$⁠. The condition that preserves the topology near extrema (⁠ $v ≥ u$⁠) is consistent with omnigeneity provided $∂ θ B 2 = 0$ at the turning points. To order $ϵ 2$⁠, defects on $| B |$ contours arising from $2 v > u > v > 1$ preserve omnigeneity provided that $∂ θ B 2 = 0$ at the turning points of B₀. Those defects remain asymptotically small (higher order than $ϵ 2$⁠).

Not all topological defects described by the near-axis expansion are equally disruptive, though. All other combinations of v and u, i.e., $v ≤ u / 2$⁠, are too disruptive, with the exception of one very special case, perhaps the simplest. For v = 1 and u = 2, the choice of $B 2 c$ consistent with QI, Eq. (32c), precisely cancels the leading disruptive contribution from the puddle. This marginal case is the only case in which the second order can directly interact with the first order to amend omnigeneity.

Interestingly, all the scenarios in which the minimum of $| B |$ is not straight (see white broken lines in Fig. 2) do also exhibit topological defects in the form of puddles. The puddle size (which one may think of in terms of the variation of B along the line of minima) scales like $ϵ 2$ in the special case and with a larger power in the case of the “allowable” puddles. In addition, only in the last special case does the $φ = 0$ line not correspond to a contour of constant $| B |$⁠. The two lower examples of Fig. 2, corresponding to this special case, demonstrate that the addition of the second order field B₂ can have the effect of “healing” the topological defects, in the sense that the large first-order (first harmonics in θ) islands are broken and replaced by higher-order structures. Although there is no hint that puddles can be completely eliminated at any finite order in ϵ, this phenomenon resolves the apparent contradiction between the asymptotic construction of QI fields and the exact concept of QI fields with a non-straight minimum, as the latter absolutely forbids such topological features. This leaves open the possibility of an exactly QI solution near the bottom of the well.

To make the above theoretical discussion more concrete, we include in Fig. 3 an illustration of the scenarios in Fig. 2 using global equilibria computed with VMEC. Details on these configurations are presented in  Appendix B. The main features identified theoretically are apparent in the contours with a clear hierarchy in the magnitude of topologic defects. In the special case, $u = 2 , v = 1$⁠, we see that the field variations beyond first order affect critically the behavior near the minimum of $| B |$⁠, changing the behavior of puddles and contours nearby.

In comparing the lower panels of these two figures, one may notice some qualitative differences in, for instance, the location of the puddles and the degree of shaping of the minimum contour. These features may be attributed to the fact that the fields in Fig. 3 are evaluated at a finite, sizeable radius, and thus, higher order contributions from the radial expansion are more significant compared to Fig. 2, which illustrates the asymptotic scenario (formally, a very small $r ≪ 1$⁠). This results in a more complicated θ-independent part of $| B |$⁠, and a more sizeable m = 2 contribution to $| B |$⁠, and hence the resulting differences. Note also that in the asymptotic scenario of Fig. 2, the m = 2 component was chosen to have a rather particular form to illustrate the possibility of healing of the puddles asymptotically, while we have not made any such choice in the global equilibrium.

In this paper, we derive the conditions of quasi-isodynamicity on the magnetic field magnitude $| B |$ near the magnetic axis. We do so by asymptotic expansion of the difference in radial drift at opposing portions of magnetic wells, providing a clear physical approach to the problem. This allows us to obtain QI conditions on the second order components of $| B |$ in the distance from the magnetic axis.

The approach and results in this paper set the ground for further exploration of quasi-isodynamic configurations and their properties using the near-axis framework, reformulating, and extending the original work in Ref. 22. This includes consideration of appropriate shape choices, MHD stability, etc.

Here, we have only considered the implications of omnigenity on the magnetic field magnitude. This is only part of the whole problem of constructing equilibrium fields, which includes also the full description of B, requiring the solution of the so-called MHD constraint equations.^18,21,35 An analysis of the consequences of the second order QI conditions in this context will be presented in the future.

The authors would like to thank Katia Camacho-Mata for providing VMEC equilibrium examples to illustrate the work on puddles. The authors would also like to acknowledge fruitful discussion with Per Helander and Rogerio Jorge.

The authors have no conflicts to disclose.

Eduardo Rodriguez: Conceptualization (lead); Formal analysis (lead); Investigation (equal); Writing – original draft (lead); Writing – review & editing (equal). Gabriel G. Plunk: Formal analysis (supporting); Investigation (equal); Writing – review & editing (equal).

The data that support the findings of this study are openly available in Zenodo, Ref. 36.

In this appendix, we consider the asymptotic considerations at and in the neighborhood of turning points of B₀. Such turning points are special as deeply trapped and barely trapped particles spend an infinite amount of time at them. For the deeply trapped this is obvious, as these particles are unable to exist anywhere else. For the barely trapped, it follows from $v ∥ ∼ 1 − ( 1 − B 0 ″ φ 2 ) ∼ φ$⁠, and the divergence of $∫ d φ / v ∥ ∼ log φ$ at the turning point. The consequence is that, to confine both classes of particles, we must make the radial drift exactly vanish at those points. Formally, $B × ∇ B · ∇ ψ = 0$ wherever $B · ∇ B = 0$⁠. This is known as the condition of pseudosymmetry.^30,31

A magnetic field that is pseudosymmetric over a given magnetic flux surface, will possess contours of constant $| B |$ all with the same topology.¹⁰ That is to say, the representation of $| B |$ as a function of $( θ , φ )$ should not present any contour that closes within a field period, i.e., features that resemble puddles, which can be regarded as topological defects. It is, thus, tempting to require such puddles not to be present at any of the asymptotic orders in which we have considered our QI construction. In the spirit of the asymptotic approach in this paper, though, it is only consistent with treat this pseudosymmetric condition asymptotically. This cannot be done through an asymptotic analysis of Eq. (7), both because the condition yields no information at the turning point of $| B |$⁠, and the expansion itself breaks down in the neighborhood.

Instead, to keep the behavior at the extrema accountable, order by order, we shall assess the location of the extrema of $| B |$ along field lines and evaluate the radial drift there. Assessing the magnitude of the radial drift we may then deem the field consistent or inconsistent with omnigeneity (and thus also pseudosymmetry) to the right order. In practice, this requires an asymptotic expansion of both $B · ∇ B$ and $B × ∇ B · ∇ ψ$⁠. Fortunately, we already have these in Eqs. (11) and (12), as we needed them to construct Y. Thus, all that remains, order-by-order, is (i) to find the turning points $φ t = φ t ( 0 ) + φ t ( 1 ) + ⋯$ and (ii) to evaluate the drift there.

In summary, then, to second order, the pseudosymmetry condition requires one of the following three:

1. For $v ≥ u$⁠: $∂ θ B 2 ( θ , φ t ( 0 ) ) = 0$⁠, which preserves the topology of the contours of $| B |$ to this order.

2. For $2 v > u > v > 1$⁠: $∂ θ B 2 ( θ , φ t ( 0 ) ) = 0$⁠, breaks the topology of $| B |$-contours with the appearance of puddles, but the derived break-down of omnigeneity is higher order.

3. For $u = 2 , v = 1$⁠: $∂ θ B 2 = ∂ θ ( ( B 1 ′ ) 2 / 2 B 0 ″ )$ at $( θ , φ t ( 0 ) )$⁠. This is satisfied by the QI conditions to second order around the bottoms of the wells. It also gives puddles.

The above consideration gives a sense of the importance of the turning points and the behavior of the various near-axis functions about them. We saw that in certain cases, the QI conditions derived in this paper do not apply close to the turning points. An example of that was the $2 v > u > v > 1$ case. Misbehavior near the minimum affects not only the deeply trapped population but also the remainder classes, which must physically traverse the region. To estimate by how much, consider that some region $Δ φ ∼ ϵ α$ is spoiled near the minimum. In that region, the drift to be order $ϵ β$⁠. Then, we expect the effect on the bounce averaged drift to be $O ( ϵ α + β )$⁠. In the $2 v > u > v > 1$ situation (with t large enough), $α ∼ 1 / ( u − v )$ and $β ∼ u / ( u − v )$⁠. Thus, the spoiling of QI will be order $( u + 1 ) / ( u − v ) ≥ 2 + 2 / ( v − 1 ) > 2$⁠, limit in which $u = 2 v − 1$⁠. (The lower limit of $u = v + 1$ gives an order of three or larger.)

In this appendix, we present the equilibrium configurations used for constructing the $| B |$ contours in Fig. 3. To do so, we follow the conventional approach to construct equilibria from near-axis constructions, which consists on evaluating the near-axis fields at a finite radius, constructing a flux surface boundary, and passing that to a global equilibrium solver, in this case VMEC³⁷ in vacuum. For the purpose of the figure, we consider the construction of these equilibria using different near-axis choices, informed by the insight gained in this paper, solving the near-axis at first order and evaluating the corresponding flux surface. This is the conventional approach, as it has been used in the construction of quasisymmetric configurations^19,38,39 and more recently also quasi-isodynamic ones.^22,33,34 The suite SIMSOPT⁴⁰ is then used to analyze the fields.

Following that approach, we specify in this section the near-axis constructions that were used to build those equilibria, as well as the radial distance r at which they were evaluated.

• QI configuration with $u = 2 , v = 3$⁠: this configuration was chosen to have a third order zero in the curvature of the axis, but a quadratic magnetic field profile B₀. To achieve the right third order zero for the curvature, we follow the prescription in Ref. 33 for how to choose the Fourier harmonics describing the magnetic axis. The axis is described, in cylindrical coordinates, by

  $R = 1 + 0.30899 cos ϕ − 0.31461 cos 2 ϕ − 9.8315 × 10 − 2 cos 3 ϕ + 3.3708 × 10 − 2 cos 4 ϕ + 1.4045 × 10 − 2 cos 5 ϕ ,$

  (B1a)

  $Z = − 1.1000 sin ϕ + 0.1000 sin 3 ϕ ,$

  (B1b)

  where $ϕ$ is the cylindrical angle.The magnetic field is chosen $B 0 = 1 + 0.15 cos φ$⁠, with the first order variation $d ¯ = d / κ = 0.73$⁠. The deviation from QI at the first order is introduced (to guarantee periodicity) following Ref. 33, with the choice of k = 2. The equilibrium surface was evaluated at r = 0.05, corresponding to an aspect ratio of A = 20 (or in VMEC, $A VMEC = 14.4$⁠). A rendition of the configuration boundary is shown in Fig. 5.

• QI configuration with $u = 6 , v = 3$⁠: this configuration was chosen to have the same magnetic axis as the previous third-order zero configuration, but with a magnetic field profile B₀ in this case with u = 6. This is chosen to illustrate the difference in the presence of topological effects arising when going from the “no-puddle” to “forbidden puddle” regime (see Fig. 4). The only difference is then the choice for $B 0 = 1 + 0.15 cos φ + 0.06 cos 2 φ + 0.01 cos 3 φ$⁠, which can straightforwardly be shown to have the desired flat minimum. The equilibrium surface was evaluated at r = 0.05, corresponding to an aspect ratio of A = 20 (or in VMEC, $A VMEC = 14.7$⁠). A rendition of the configuration boundary is shown in Fig. 6.

• QI configuration with $u = 2 , v = 1$⁠: this configuration was chosen to be of the “special” kind in which second order effects can most notably affect the leading order puddle and minimum structure. For this case, we chose a configuration presented in Ref. 33, which has the right order of zero and magnetic field. For completeness,

  $R = 1 + 0.058824 cos 4 φ ,$

  (B2a)

  $Z = 0.3921 sin 2 φ + 4.90 × 10 − 3 sin 4 φ .$

  (B2b)

  The B₀ and $d ¯$ profiles are the same as the previous quadratic minimum one. The omnigeneity control parameter was set to k = 2. The equilibrium surface was evaluated at r = 0.05, corresponding to an aspect ratio of A = 20 (or in VMEC, $A VMEC = 18.3$⁠). A rendition of the configuration boundary is shown in Fig. 7.