1 The need for and the origin of the rotational transform

Magnetostatic equilibria are characterized by the following equations:

(1.1a-c)\begin{equation} \boldsymbol{\nabla} p = {\boldsymbol{J}} \times {\boldsymbol{B}}, \quad \boldsymbol{\nabla} \times {\boldsymbol{B}} = {\boldsymbol{J}}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{B}} = 0, \end{equation}

where ${\boldsymbol {B}}$ is the magnetic field, ${\boldsymbol {J}}$ is the current density and $p$ is the scalar pressure. With the assumption that $p$ is not constant in a small region and the surface lies in a bounded volume of space, the boundary must be topologially toroidal (Kruskal & Kulsrud Reference Kruskal and Kulsrud1958). The phase portrait of ${\boldsymbol {B}}$, where the magnetic field lines are treated like integral curves of a Hamiltonian dynamical system, is characterized by the topology of the level sets of $p$. A magnetic field line action can be defined (Cary & Littlejohn Reference Cary and Littlejohn1983), which serves as a starting point for the Lagrangian integration carried out in this paper.

For magnetic confinement of plasmas in toroidal geometries, that the magnetic field lines rotate poloidally (the short way) as they rotate toroidally (the long way) around the torus is essential for cancelling charged particle drifts, which would otherwise lead to loss (Spitzer Reference Spitzer1958). The number of poloidal rotations that a field line achieves per toroidal period is called the rotational transform, $\iota$ (Spitzer Reference Spitzer1958).

Kruskal & Kulsrud consider the case where ‘$p$ is reasonably smooth and not constant in any region’, and this entails that the magnetic field lines lie on nested flux surfaces. In the following, we do not place such a constraint on the pressure and we do not require that the magnetic field has nested flux surfaces. The existence of nested magnetic surfaces, their link with continuous symmetries and the preservation of magnetic surfaces under small symmetry-destroying perturbation are widely discussed in Kolmogorov (Reference Kolmogorov1954), Arnol'd (Reference Arnol'd1963) and Moser (Reference Moser1962) and also in texts and reviews: Moser (Reference Moser1973), Arnold (Reference Arnold1978), Lichtenberg & Lieberman (Reference Lichtenberg and Lieberman1992) and Meiss (Reference Meiss1992). Herein, we only require that the magnetic field is ‘toroidal’, which we describe as follows. We consider a magnetic configuration in the domain $\varOmega$, which is assumed to be a solid torus, with ${\boldsymbol {B}}\boldsymbol {\cdot } \boldsymbol {\nabla } \phi >0$ everywhere where $\phi$ is the toroidal angle, and with ${\boldsymbol {B}}\boldsymbol {\cdot }{\boldsymbol {n}}=0$ where ${\boldsymbol {n}}$ is normal to the boundary. For any poloidal section $\varSigma$ of $\varOmega$, the Poincaré return map is defined by the intersection of field lines and $\varSigma$ after one toroidal period. Brouwer's fixed-point theorem (Brouwer Reference Brouwer1910) ensures that there will be at least one fixed point of the Poincaré first return map.

Taking such a point as the origin of polar-like coordinates, e.g. $(r,\theta )$, we may consider the effect of iterating the mapping on nearby points, where the point $(r_n,\theta _n)$ gets mapped to $(r_{n+1},\theta _{n+1})$. If, upon iterating the return map, the points rotate around the fixed point, then the fixed point is called elliptic. (A more rigorous definition of elliptic fixed points is given in Meiss (Reference Meiss1992).) There may be more than one fixed point, and not all fixed points are elliptic. Even though the methods of Lagrangian integration as described below can be applied more generally, for clarity of exposition we restrict our attention to the case of an isolated elliptic fixed point. Stellarators are typically designed to have one easily identifiable elliptic fixed point about which most field lines rotate, and this is called the magnetic axis.

This sequence of angular displacements, $\Delta \theta = \theta _{n+1}-\theta _n$, enables one to define the rotational transform of a field line:

(1.2)\begin{equation} \iota := \lim_{N\uparrow \infty}\frac{1}{2{\rm \pi} N}\sum_n^N\Delta\theta_n. \end{equation}

Even though the angle coordinate becomes degenerate at the origin in polar-like coordinates, the common approach is to define the on-axis rotational transform by taking the limit as the starting point gets closer to the axis.

Note that the this definition of the rotational transform is purely geometric, in that the rotational transform effectively measures how many times a given magnetic field line links the magnetic axis. A discussion of the rotational transform as an asymptotic linking of neighbouring field lines can be found in Arnold (Reference Arnold2014).

Mercier (see Mercier Reference Mercier1964; Helander Reference Helander2014) expressed the on-axis rotational transform as an integral along the axis:

(1.3)\begin{equation} \iota = N + \frac{1}{2{\rm \pi} } \oint \frac{{\rm d}\ell}{\cosh{\eta}} \left( \frac{J_0}{2B_0}+\delta' - \tau\right) , \end{equation}

with $\ell$ denoting the arc length. The on-axis current density and magnetic field are denoted by $J_0$ and $B_0$, respectively, and the torsion of the axis by $\tau$. The eccentricity of flux surfaces is described by $\eta$ and $\delta$ is a parameter describing their rotation around the axis, with $\delta ' := {\rm d}\delta /{\rm d}\ell$. Here $N$ is an integer coming from the phase of the rotation term $\delta$ (see Pfefferlé et al. Reference Pfefferlé, Gunderson, Hudson and Noakes2018). Note that Spitzer (Reference Spitzer1958) had identified independently that a way to generate some rotational transform is to give torsion to an axis. This expression provides invaluable insight, showing how rotational transform can be produced by plasma currents, as is used by tokamaks, or by geometrical shaping, as is used by stellarators.

One may also encounter different definitions of the rotational transform. Assuming that the magnetic configuration possesses magnetic surfaces, ensuring that the toroidal and poloidal magnetic fluxes $\psi$ and $\chi$ respectively, can be properly defined, then as shown by Mercier et al. (Reference Mercier1974), $\iota$ can also be expressed as the ratio of differential fluxes:

(1.4)\begin{equation} \iota = \frac{{\rm d}\chi}{{\rm d}\psi}. \end{equation}

Mercier's expression (1.3) for instance was derived starting from the flux definition (1.4). The fluxes were expanded in a power series of the distance to the axis by constructing custom polar coordinates, so-called Mercier coordinates, which are described in Appendix C.

In this paper, we present a derivation of Mercier's formula using Lagrangian integration. This approach expresses the rotational transform as a Floquet exponent of the field lines. The relationship between the rotational transform and the Floquet exponent was described by Greene (Reference Greene1979), and more recently by Duignan & Meiss (Reference Duignan and Meiss2021). In § 2, the magnetic field-line action is defined. Stationary curves of the action are shown to be field lines, enabling one to identify a magnetic axis. Assumption is made in this paper that the axis is elliptical, but the formalism can be applied to the hyperbolic case. In § 3, Mercier's formula (1.3) is derived from the second variation of the action. The result is obtained through a near-axis expansion of the null eigenspace of the second variation operator. The periodicity enables the use of a Floquet description of the solutions. This same result is derived through the theory of the Hill infinite determinant, analogous to Schrödinger's equation in a periodic potential. In both cases, the rotational transform is shown to be a Floquet exponent of the null eigenspace of the second variation operator. In § 4, a discrete formalism is introduced so that the Floquet exponents can be solved for numerically.

2 Description of the magnetic field-line action

Let us start by defining the magnetic field line action, as introduced by Cary & Littlejohn (Reference Cary and Littlejohn1983). This is defined as a line integral that, for a given magnetic field, depends only the integration contour. Herein, the integration contours are assumed to be curves that close after one toroidal period. Let us consider a closed, differentiable curve $\mathcal {C} \subset \mathbb {R}^3$ with total length $L$, closing after one toroidal transit. The latter assumption that the curve closes after one toroidal period is a necessary condition for a magnetic axis. Curve ${\mathcal {C}}$ is parametrized by the $C^1$, $L$-periodic vector-valued function

(2.1)\begin{equation} \left. \begin{gathered} {\boldsymbol{x}}:[0,L]\rightarrow \mathbb{R}^3\\ \ell \mapsto {\boldsymbol{x}}(\ell)\in {\mathcal{C}}, \end{gathered} \right\} \end{equation}

$\ell$ being the arc length and we note ${\boldsymbol {x}}' := \textrm {d}{\boldsymbol {x}}/\textrm {d}\ell$. The action is defined as the circulation along $\mathcal {C}$ of the magnetic vector potential ${\boldsymbol {A}}$, with ${\boldsymbol {B}} = \boldsymbol {\nabla }\times {\boldsymbol {A}}$ in $\varOmega$:

(2.2)\begin{equation} \mathcal{S} := \oint_{\mathcal{C}}\,{\rm d}\ell\, {\boldsymbol{A}} \boldsymbol{\cdot} {\boldsymbol{x}}'. \end{equation}

Properties of the magnetic field are accessible through calculus of variations, performed on the action. The variations are performed with respect to changes in the geometry of the curve. The first variation with respect to a variation $\delta {\boldsymbol {x}}$ is

(2.3)\begin{equation} \delta \mathcal{S}[\delta {\boldsymbol{x}}] = \oint_{\mathcal{C}} \,{\rm d}\ell\, {\boldsymbol{x}}^\prime \times {{\boldsymbol{B}}} \boldsymbol{\cdot} \delta {\boldsymbol{x}}, \end{equation}

which shows that stationary curves are tangential to the magnetic field and hence field lines. We used that $\delta {\boldsymbol {A}} [\delta {\boldsymbol {x}}]=\delta {\boldsymbol {x}} \boldsymbol {\cdot }\boldsymbol {\nabla } {\boldsymbol {A}}$. From (2.3), a magnetic axis can be identified. We focus in particular on elliptical axes, but one can apply the following formalism to hyperbolic axes as well. The difference would be seen in the flux-surface functions around the axis, and in the nature of the Floquet exponents. Let us denote an elliptical axis by ${\mathcal {C}}_a$. Additional properties of the field appear at higher orders of variations of $\mathcal {S}$. For the rotational transform in particular, the second-order variation $\delta ^2 {\mathcal {S}}$ is of interest.

3 Verification of Mercier's formula for the Floquet exponent

3.1 The second variation as an operator

In order to express the on-axis rotational transform, the second-order variation of the action applied to the axis needs to be derived. Assume that an elliptical axis ${\mathcal {C}}_a$ has been found as a stationary curve of the action. The formalism is identical for hyperbolic axes, but in this paper, we focus on the elliptic case. The second variation of the action performed from the axis ${\mathcal {C}}_a$ is

(3.1)\begin{equation} \delta^2 \mathcal{S}[\delta {\boldsymbol{x}}] = \oint_{\mathcal{C}_a} \,{\rm d}\ell\,\delta({\boldsymbol{x}}^\prime \times {{\boldsymbol{B}}} \boldsymbol{\cdot} \delta{{\boldsymbol{x}}}) = \oint_{\mathcal{C}_a} \,{\rm d}\ell\,\delta {\boldsymbol{x}} \boldsymbol{\cdot} (\delta {\boldsymbol{x}}^\prime \times {\boldsymbol{B}} + {\boldsymbol{x}}' \times \delta {\boldsymbol{B}}), \end{equation}

where $\boldsymbol {f}':=\textrm {d}\boldsymbol {f}/\textrm {d}\ell$ for any $\boldsymbol {f}$. Using $\delta {\boldsymbol {B}} = \delta {\boldsymbol {x}} \boldsymbol {\cdot } \boldsymbol {\nabla } {\boldsymbol {B}}$ and the Einstein summation convention, we write the second variation as an operator:

(3.2)\begin{equation} \delta^2 \mathcal{S} = \oint_{{\mathcal{C}}_a} \,{\rm d}\ell\, \delta {\boldsymbol{x}}^i \frac{\delta^2 S}{\delta {\boldsymbol{x}}^i \delta {\boldsymbol{x}}^j} \delta {\boldsymbol{x}}^j , \end{equation}

where

(3.3)\begin{equation} \frac{\delta^2 \mathcal{S}}{\delta {\boldsymbol{x}}^i \delta {\boldsymbol{x}}^j} = \epsilon_{ijk}{\boldsymbol{B}}^k\frac{{\rm d}}{{\rm d}\ell} + \epsilon_{imk}{\boldsymbol{x}}'^m \partial_j {\boldsymbol{B}}^k, \end{equation}

and $i,j,k,m \in \{1,2,3\}$, which in matrix form reads

(3.4)\begin{equation} {\boldsymbol{\mathsf{M}}} := \frac{\delta^2 \mathcal{S}}{\delta {\boldsymbol{x}} \delta {\boldsymbol{x}}}={-}\left(\boldsymbol{\mathsf{I}}\times {\boldsymbol{B}} \right) \frac{{\rm d}}{{\rm d}\ell}+ {\boldsymbol{x}}' \times (\boldsymbol{\nabla} {\boldsymbol{B}})^\mathsf{T}. \end{equation}

The identity tensor is denoted by $\boldsymbol{\mathsf{I}}$. Note that the covariant tensor $\boldsymbol{\mathsf{M}}$ can be easily symmetrized (see Hudson & Dewar Reference Hudson and Dewar2009). However, the direction in which the derivative $\textrm {d}/\textrm {d}\ell$ is taken has to be chosen carefully. We continue with the non-symmetric form (3.4).

We show that the null eigenspace of $\boldsymbol{\mathsf{M}}$ is of particular interest. Indeed, by periodicity of the field lines considered to compute the action, the eigenspaces of the operator $\boldsymbol{\mathsf{M}}$ are related to geometric properties of the neighbouring field lines, including the rotational transform. Let ${\boldsymbol {v}}$ be a null eigenfunction of $\boldsymbol{\mathsf{M}}$ such that

(3.5)\begin{equation} \boldsymbol{\mathsf{M}}{\boldsymbol{v}} =\boldsymbol{0}. \end{equation}

For ${\boldsymbol {v}}$ to be non-trivial, (3.5) rewrites

(3.6)\begin{equation} \det\left( \boldsymbol{\mathsf{M}} \right) =0. \end{equation}

We will demonstrate that the condition (3.5) is satisfied by solutions $\boldsymbol {v}$, such that their Floquet exponent is the on-axis rotational transform, and (3.6) will serve in the discrete formalism introduced in § 4. To solve (3.5), the near-axis formalism developed by Mercier, Solov'ev and Shafranov (Solov'ev & Shafranov Reference Solov'ev and Shafranov1970; Mercier et al. Reference Mercier1974) is used. Near-axis formalism in the inverse coordinate approach, where the flux surface $\psi$ is used as a coordinate, has proved to be very successful in understanding quasi-symmetry (see Garren & Boozer Reference Garren and Boozer1991; Landreman & Sengupta Reference Landreman and Sengupta2018, Reference Landreman and Sengupta2019; Rodríguez, Sengupta & Bhattacharjee Reference Rodríguez, Sengupta and Bhattacharjee2023). However, for this work, the near-axis formalism in direct or Mercier–Solov'ev–Shafranov coordinates (Jorge, Sengupta & Landreman Reference Jorge, Sengupta and Landreman2020b,Reference Jorge, Sengupta and Landremana; Sengupta et al. Reference Sengupta, Rodriguez, Jorge, Bhattacharjee and Landreman2024) is more relevant.

3.2 Derivation of the Floquet exponent from the second variation

The operator is expanded in the Solov'ev–Shafranov coordinates, introduced in Appendix C.2. The latter set of coordinates is closely related to the Mercier coordinates, presented in Appendix C.1.

In what follows, it is assumed that the magnetic axis is a closed curve ${\mathcal {C}}_a\subset \mathbb {R}^3$ parametrized by the vector-valued function ${\boldsymbol {r}}_0$, with the arc length $\ell$ as parameter, and a total length $L$. The basis of expansion is taken to be $\{\boldsymbol {e}_1, \boldsymbol {e}_2, \boldsymbol {e}_3\}$, whose expressions in terms of $\boldsymbol {\mathcal {N}}$, $\boldsymbol {\mathcal {B}}$ and ${\boldsymbol {t}}$ are given in Appendix C.2, but the results will always be expressed in terms of the Solov'ev–Shafranov vectors $\{\boldsymbol {\mathcal {N}}, \boldsymbol {\mathcal {B}}, {\boldsymbol {t}}\}$. Additionally, the following notation is adopted: $f'=\textrm {d}f/\textrm {d}\ell$ for any $f$.

3.2.1 Magnetic field expansion near an elliptic magnetic axis

From (3.4), the magnetic field ${\boldsymbol {B}}$ needs to be described in the vicinity of the axis ${\boldsymbol {r}}_0$ for $\boldsymbol{\mathsf{M}}$ to be expanded in the near-axis formalism. The expansion of ${\boldsymbol {B}}$ in Solov'ev–Shafranov coordinates is carried in the limit $x,y\ll 1$, or in other words, is limited to linear terms only. The magnetic field can be written in contravariant form as

(3.7)\begin{equation} \sqrt{|g|}{\boldsymbol{B}}= \sqrt{|g|}B^1 \boldsymbol{e}_1+\sqrt{|g|}B^2 \boldsymbol{e}_2+\sqrt{|g|}B^3 \boldsymbol{e}_3 , \end{equation}

with

(3.8a-c)\begin{equation} \sqrt{|g|}B^1= a_1 x+a_2 y , \quad \sqrt{|g|}B^2= b_1 x+b_2 y, \quad \sqrt{|g|}B^3= B_0+c_1 x+c_2 y, \end{equation}

where $a_1,a_2,b_1,b_2, c_1,c_2$ are periodic functions of $\ell$, $B_0$ the zeroth-order magnitude of ${\boldsymbol {B}}$ and $\sqrt {|g|} = h$ as defined in Appendix C.2. Therefore, in the limit $x,y\ll 1$, $h=1$ and the magnetic field reads

(3.9)\begin{gather} {\boldsymbol{B}} = \left(B_0 +c_1 x+c_2 y\right){\boldsymbol{t}} + \left[a_1 x+(a_2 +u' B_0 )y\right]\boldsymbol{\mathcal{N}}\nonumber\\ +\left[(b_1-u' B_0) x+b_2 y\right]\boldsymbol{\mathcal{B}}+ \boldsymbol{\mathcal{O}}(x^2+y^2). \end{gather}

A direct calculation shows that

(3.10)\begin{equation} \left. \begin{aligned} \boldsymbol{\nabla} {\boldsymbol{B}} & =\kappa B_0 {\boldsymbol{t}} {\boldsymbol{n}} +a_1 \boldsymbol{\mathcal{N}}\boldsymbol{\mathcal{N}}+b_2 \boldsymbol{\mathcal{B}}\boldsymbol{\mathcal{B}}+B_0' {\boldsymbol{t}}{\boldsymbol{t}} + (a_2+u' B_0) \boldsymbol{\mathcal{B}}\boldsymbol{\mathcal{N}}\\ & \quad +(b_1-u' B_0)\boldsymbol{\mathcal{N}}\boldsymbol{\mathcal{B}} + (c_1 \boldsymbol{\mathcal{N}}+c_2 \boldsymbol{\mathcal{B}})\boldsymbol{t},\\ \boldsymbol{x'}\times (\boldsymbol{\nabla} \boldsymbol{B})^\mathsf{T} & =\kappa B_0 {\boldsymbol{b}} {\boldsymbol{t}} +\boldsymbol{\mathcal{B}} \left[a_1 \boldsymbol{\mathcal{N}}+(a_2 +u' B_0) \boldsymbol{\mathcal{B}} \right]-\boldsymbol{\mathcal{N}} \left[(b_1-u' B_0) \boldsymbol{\mathcal{N}}+b_2 \boldsymbol{\mathcal{B}} \right], \end{aligned} \right\} \end{equation}

which yields, using dyadic algebra

(3.11)\begin{equation} \left. \begin{gathered} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{B}} = \boldsymbol{\mathsf{I}}\boldsymbol{:}\boldsymbol{\nabla} {\boldsymbol{B}} = B_0' + a_1+b_2,\\ \boldsymbol{\nabla}\times {\boldsymbol{B}} ={-}\boldsymbol{\mathsf{I}}\boldsymbol{\overset{\times}{.}}\boldsymbol{\nabla} {\boldsymbol{B}} = (b_1-a_2-2 u' B_0) {\boldsymbol{t}}+ (c_2 -\kappa B_0 \sin\delta )\boldsymbol{\mathcal{N}}\\ - (c_1-\kappa B_0 \cos \delta)\boldsymbol{\mathcal{B}}. \end{gathered} \right\} \end{equation}

Here $\kappa$ denotes the local curvature of the axis ${\boldsymbol {r}}_0$. Additionally, (C10a,b) has been used to rewrite ${\boldsymbol {b}}$ in terms of $\boldsymbol {\mathcal {N}}$ and $\boldsymbol {\mathcal {B}}$. Requiring that $\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {B}}=0$ and $\boldsymbol {\nabla } \times {\boldsymbol {B}}=J_0 {\boldsymbol {t}}$, one finds

(3.12a-d)\begin{equation} a_1+b_2={-}B_0', \quad \frac{b_1-a_2}{2B_0}=u'+\frac{J_0}{2B_0}, \quad c_1 =\kappa B_0 \cos \delta, \quad c_2 =\kappa B_0 \sin \delta. \end{equation}

So far, only two constraints have been derived for the four expansion functions $a_1,a_2,b_1,b_2$. Turning to the Mercier representation (D2a,b) and (D7a–c) of ${\boldsymbol {B}}$ brings in additional conditions. Together with the definitions of $\boldsymbol {\mathcal {N}}$ and $\boldsymbol {\mathcal {B}}$,

(3.13)\begin{align} \frac{{\boldsymbol{B}}}{B_0} & =\quad {\boldsymbol{t}}\left[1+\kappa \sqrt{x^2+y^2} \cos{(u-\delta)}\right] \nonumber\\ & \quad+ \boldsymbol{\mathcal{N}} \left[ -\frac{1}{2}\left(\frac{B_0'}{B_0}+\eta'\right)x + \left(\frac{-J_0}{2 B_0}+ \frac{\dfrac{J_0}{2B_0}- \tau +\delta'}{ \cosh{\eta}} \sinh{\eta}\right) y\right]\nonumber\\ & \quad+ \boldsymbol{\mathcal{B}} \left[ \left(\frac{J_0}{2 B_0}+ \frac{\dfrac{J_0}{2B_0}- \tau +\delta'}{ \cosh{\eta}} \sinh{\eta}\right) x-\frac{1}{2}\left(\frac{B_0'}{B_0}-\eta'\right)y \right]\nonumber\\ & \quad + \boldsymbol{\mathcal{O}}(x^2+y^2), \end{align}

where $\delta$, as explained in Appendix C, describes the rotation of elliptical flux surfaces around the expansion axis ${\boldsymbol {r}}_0$. As for $\eta$, it describes the eccentricity of the flux surfaces around ${\boldsymbol {r}}_0$. Defining

(3.14)\begin{equation} \varOmega_0(\ell)=\frac{\dfrac{J_0}{2B_0}- \tau +\delta'}{ \cosh{\eta}}, \end{equation}

the comparison between the Mercier representation (3.13) and the Solov'ev–Shafranov representation (3.9) yields for the expansion coefficients

(3.15)\begin{equation} \left. \begin{aligned} \frac{a_1}{B_0} & ={-}\frac{1}{2}\left(\frac{B_0'}{B_0}+\eta'\right), \quad \frac{a_2}{B_0}={-}\left(u'+\frac{J_0/2}{ B_0}\right)+\varOmega_0\sinh{\eta},\\ \frac{b_2}{B_0} & ={-}\frac{1}{2}\left(\frac{B_0'}{B_0}-\eta'\right), \quad \frac{b_1}{B_0}=\left(u'+\frac{J_0/2}{ B_0}\right)+\varOmega_0\sinh{\eta}, \end{aligned} \right\} \end{equation}

which satisfy the relations (3.12a–d). From the definition of $\varOmega _0$, we can further simplify $a_2$ and $b_1$ as

(3.16a,b)\begin{equation} \frac{a_2}{B_0}={-}\varOmega_0 {\rm e}^{-\eta}, \quad \frac{b_1}{B_0}=\varOmega_0 {\rm e}^{\eta}. \end{equation}

As a check of correctness for the above expressions, one can show that ${\boldsymbol {B}}\boldsymbol {\cdot } \boldsymbol {\nabla } \psi =0$, where

(3.17)\begin{align} \sqrt{g} {\boldsymbol{B}}\boldsymbol{\cdot} \boldsymbol{\nabla} & = \sqrt{g}B^1\partial_x +\sqrt{g}B^2\partial_y+\sqrt{g}B^3\partial_\ell\nonumber\\ & = B_0\partial_\ell + (a_1 x + a_2 y)\partial_x + (b_1 x + b_2 y)\partial_y. \end{align}

The flux surface function $\psi$ for the elliptic case is given in the Mercier representation by (D10) or equivalently in the Solov'ev–Shafranov representation by

(3.18)\begin{equation} \psi(x,y)=B_0({\rm e}^{\eta}x^2 + {\rm e}^{-\eta}y^2). \end{equation}

Note that in the case of a hyperbolic axis, the plus sign from the flux function (3.18) would be replaced by a minus sign.

3.2.2 The second variation tensor and its null eigenspace

Now that the magnetic field has been properly expressed and expanded in Solov'ev–Shafranov coordinates, the obtained ${\boldsymbol {B}}$ and $\boldsymbol {\nabla } {\boldsymbol {B}}$ can be substituted in (3.4) to expand the second variation tensor. From (3.10), we find that it is given by

(3.19)\begin{equation} \left. \begin{aligned} \boldsymbol{\mathsf{M}} & := \boldsymbol{\mathsf{M}}_{{1}} \frac{{\rm d}}{{\rm d}\ell}+\boldsymbol{\mathsf{M}}_{{2}},\\ \boldsymbol{\mathsf{M}}_{{1}} & = (\boldsymbol{\mathcal{N}} \boldsymbol{\mathcal{B}}-\boldsymbol{\mathcal{B}} \boldsymbol{\mathcal{N}})B_0,\\ \boldsymbol{\mathsf{M}}_{{2}} & = \kappa B_0 \left[\boldsymbol{\mathcal{B}}\cos{\delta} - \boldsymbol{\mathcal{N}}\sin{\delta} \right]{\boldsymbol{t}} +\boldsymbol{\mathcal{B}} \left[a_1 \boldsymbol{\mathcal{N}}+(a_2+u' B_0) \boldsymbol{\mathcal{B}} \right]\\ & \quad -\boldsymbol{\mathcal{N}} \left[(b_1-u' B_0) \boldsymbol{\mathcal{N}}+b_2 \boldsymbol{\mathcal{B}} \right]. \end{aligned} \right\} \end{equation}

Observing the dyadic form of $\boldsymbol{\mathsf{M}}$, it can be shown that tangential components $v^t$ of the eigenvectors ${\boldsymbol {v}}=(v^t, v^{{\mathcal {N}}}, v^{{\mathcal {B}}})^{\mathsf {T}}$ will not contribute to $\boldsymbol{\mathsf{M}}{\boldsymbol {v}} =\boldsymbol {0}$. Alternatively, the underlying argument is that the arc length parametrization of ${\mathcal {C}}$ constrains the admissible variations $\delta {\boldsymbol {x}} = {\boldsymbol {v}}$, removing the dependence in the tangential component. Therefore, $v^t$ can be absorbed by redefining $v^{{\mathcal {N}}}$ and $v^{{\mathcal {B}}}$. Thus, the null eigenvectors are chosen to be of the form ${\boldsymbol {v}}= v^{{\mathcal {N}}} \boldsymbol {\mathcal {N}}+ v^{{\mathcal {B}}} \boldsymbol {\mathcal {B}}$. We note that

(3.20)\begin{equation} \frac{{\rm d}}{{\rm d}\ell}{\boldsymbol{v}} = \left(v^{{\mathcal{N}}\prime}+u' v^{{\mathcal{B}}}\right)\boldsymbol{\mathcal{N}} + \left(v^{{\mathcal{B}}\prime}-u' v^{{\mathcal{N}}}\right)\boldsymbol{\mathcal{B}} +\kappa\left[v^{{\mathcal{N}}} \cos{(u-\delta)} + v^{{\mathcal{B}}} \sin{(u-\delta)}\right]{\boldsymbol{t}}. \end{equation}

The tangential terms $\propto {\boldsymbol {t}}$ are not relevant since they do not contribute to $\boldsymbol{\mathsf{M}} \boldsymbol {v} =\boldsymbol {0}$. The equations for $v^{{\mathcal {N}}},v^{{\mathcal {B}}}$ then read

(3.21)\begin{equation} B_0 \begin{pmatrix} 0 & +1 \\-1 & 0 \end{pmatrix} \dfrac{{\rm d}}{{\rm d} \ell} \begin{pmatrix} v^{{\mathcal{N}}}\\ v^{{\mathcal{B}}} \end{pmatrix} + \begin{pmatrix} -b_1 & -b_2\\ a_1 & a_2 \end{pmatrix} \begin{pmatrix} v^{{\mathcal{N}}}\\ v^{{\mathcal{B}}} \end{pmatrix} =0. \end{equation}

From (3.21), by periodicity of $a_i/B_0$ and $b_i/B_0$, $i=1,2$, the system can immediately be rewritten in the form of a periodic system:

(3.22a,b)\begin{equation} \frac{{\rm d}{\boldsymbol{v}}}{{\rm d}\ell}=\boldsymbol{\mathsf{A}}(\ell) {\boldsymbol{v}}, \quad {\boldsymbol{v}} = \left(\begin{array}{c}v^{{\mathcal{N}}} \\v^{{\mathcal{B}}}\end{array}\right). \end{equation}

Here, $\boldsymbol{\mathsf{A}}(\ell )$ is periodic in $\ell$. The Floquet theorem, (A6), can be applied to conclude that the solution must be of the form

(3.23)\begin{equation} {\boldsymbol{v}} = \boldsymbol{\mathsf{U}}(\ell){\rm e}^{\boldsymbol{\mathsf{C}} \ell/L}, \end{equation}

where $\boldsymbol{\mathsf{U}}$ is a symplectic periodic matrix (with period $L$) and $\boldsymbol{\mathsf{C}}$ is a constant Hamiltonian matrix (Duignan & Meiss Reference Duignan and Meiss2021). The eigenvalues of $\boldsymbol{\mathsf{C}}$ are the Floquet exponents, which must be purely imaginary near an elliptic axis and real in the hyperbolic case.

Although the system is already in a form that allows one to solve for the exponents $\nu$ by means of (3.23), identifying the matrices $\boldsymbol{\mathsf{C}}$ and $\boldsymbol{\mathsf{U}}$ can be somewhat troublesome. For that reason, (3.21) can be alternatively rewritten in the following ordinary differential equation (ODE) form:

(3.24)\begin{equation} \frac{{\rm d}v^{{\mathcal{N}}}}{a_1 v^{{\mathcal{N}}} + a_2 v^{{\mathcal{B}}}}= \frac{{\rm d}v^{{\mathcal{B}}}}{b_1 v^{{\mathcal{N}}} + b_2 v^{{\mathcal{B}}}} = \frac{{\rm d} \ell}{B_0}, \end{equation}

which shows that the eigenvector ${\boldsymbol {v}}$ satisfies the general characteristic equation along the magnetic field:

(3.25)\begin{equation} \frac{{\rm d}{\mathcal{X}}}{\sqrt{g}B^1({\mathcal{X}},{\mathcal{Y}},\ell)} = \frac{{\rm d}{\mathcal{Y}}}{\sqrt{g}B^2({\mathcal{X}},{\mathcal{Y}},\ell)}= \frac{{\rm d}\ell}{\sqrt{g}B^3({\mathcal{X}},{\mathcal{Y}},\ell)}. \end{equation}

The eigenvector ${\boldsymbol {v}}$ components $v^{{\mathcal {N}}},v^{{\mathcal {B}}}$ can be identified with ${\mathcal {X}},{\mathcal {Y}}$, the displacements of the magnetic field line from the closed field line ${\boldsymbol {r}}_0$ along the rotated normal and binormal directions. The expansion coordinates $(x,y)$ ought not to be confused with $({\mathcal {X}},{\mathcal {Y}})$, which are the solutions of the characteristic ODEs (3.25). The equations for $({\mathcal {X}},{\mathcal {Y}})$ are

(3.26)\begin{equation} \left. \begin{aligned} & {\mathcal{X}}'+\frac{1}{2}\left(\frac{B_0'}{B_0}+\eta'\right){\mathcal{X}} + \varOmega_0 {\rm e}^{-\eta}{\mathcal{Y}}=0,\\ & {\mathcal{Y}}'+\frac{1}{2}\left(\frac{B_0'}{B_0}-\eta'\right){\mathcal{Y}} - \varOmega_0 {\rm e}^{+\eta}{\mathcal{X}}=0. \end{aligned} \right\} \end{equation}

Introducing the variables $X,Y$ defined through

(3.27a,b)\begin{equation} {\mathcal{X}}= \frac{1}{\sqrt{B_0}}{\rm e}^{-\eta/2}X(\ell), \quad {\mathcal{Y}}= \frac{1}{\sqrt{B_0}}{\rm e}^{+\eta/2}Y(\ell), \end{equation}

which is possible since the periodicity is conserved, the system reduces to the simple harmonic oscillator:

(3.28a-c)\begin{equation} X' + \varOmega_0 Y =0, \quad Y'-\varOmega_0 X =0, \quad \varOmega_0 = \frac{\dfrac{J_0/2}{B_0}- \tau +\delta'}{ \cosh{\eta}}, \end{equation}

with ‘time-dependent’ frequency $\varOmega _0(\ell )$, $\ell$ being the time-like parameter. Using the complex variable $Z=X+\textrm {i}\, Y$, the system reshapes as a single complex ODE:

(3.29)\begin{equation} Z'-{\rm i}\, \varOmega_0 Z =0 \quad \Rightarrow \quad Z(\ell)= Z_0\exp{\int_0^\ell {\rm i}\,\varOmega_0(s) \,{\rm d}s}. \end{equation}

Separating the periodic and non-periodic parts of the exponential, we obtain

(3.30)\begin{equation} \left. \begin{gathered} Z(\ell) = Z_p(\ell) {\rm e}^{{\rm i}\, \nu \ell/L}, \quad Z_p(\ell)= Z_0 \exp{\int_0^\ell {\rm i}\,\tilde{\varOmega}_0(s) \,{\rm d}s},\\ \bar{\varOmega} := \frac{1}{L}\int_0^L \varOmega_0(s) \,{\rm d}s, \quad \tilde{\varOmega}_0= \varOmega_0 - \bar{\varOmega}. \end{gathered} \right\} \end{equation}

Comparing (3.30) with the Floquet form (3.23), we find that $\nu$ is given by

(3.31)\begin{equation} \nu = \int_0^L \varOmega_0(s) \,{\rm d}s \pmod{2{\rm \pi}} = \oint \frac{\dfrac{J_0(s)/2}{B_0(s)}- \tau(s) +\delta'(s)}{ \cosh{\eta(s)}} \,{\rm d}s \pmod{2{\rm \pi}} = 2{\rm \pi} \iota.\end{equation}

This matches the expression for the rotational transform (1.3) up to a factor $2{\rm \pi}$, as $\nu$ represents an angle and $\iota$ the number of turns that this angle constitutes. It follows that solving for the null eigenspace of the second variation tensor expanded in the Solov'ev–Shafranov near-axis formalism yields the correct on-axis rotational transform. The exact same approach can be followed in the hyperbolic case.

3.2.3 Derivation of the Floquet exponent from the Hill's determinant equation

Here we pursue an alternative approach to obtain the Floquet exponent, using the theory of the Hill's infinite determinant. This way of solving for the Floquet exponents of a periodic system has been known for a long time (Magnus Reference Magnus1953); however, to the authors’ knowledge, it is the first time that such an approach is used from the magnetic action. We start with the system (3.26) written as

(3.32a,b)\begin{equation} {\mathcal{X}}'=\frac{1}{B_0}(a_1 {\mathcal{X}} +a_2 {\mathcal{Y}}), \quad {\mathcal{Y}}'=\frac{1}{B_0}(b_1 {\mathcal{X}} +b_2 {\mathcal{Y}}). \end{equation}

Eliminating ${\mathcal {Y}}$ from (3.32a,b), we obtain the following second-order ODE for ${\mathcal {X}}$:

(3.33)\begin{equation} {\mathcal{X}}'' +2 C_1 {\mathcal{X}}' + C_2{\mathcal{X}} =0, \end{equation}

where

(3.34a-c)\begin{equation} 2C_1 ={-}\frac{a_2'}{a_2}+2\frac{B_0'}{B_0}, \quad C_2 = {\mathcal{D}} -\frac{a_2}{B_0}\left(\frac{a_1}{a_2}\right)', \quad {\mathcal{D}} = \frac{1}{B_0^2}(a_1 b_2 -a_2 b_1). \end{equation}

We can eliminate the first derivative term from (3.33) by the change of variables

(3.35)\begin{equation} {\mathcal{X}} = \exp{\left(-\int C_1 \,{\rm d}\ell\right)}\varPsi = \frac{\sqrt{a_2}}{B_0}X, \end{equation}

which leads to

(3.36a,b)\begin{equation} \varPsi'' + \omega^2 \varPsi=0, \quad \omega^2 \equiv C_2 -C_1'-C_1^2. \end{equation}

Equation (3.36a,b) is in the form of Hill's equation or a Schrödinger equation with a periodic potential $\omega ^2$. We note that the linear transformation (3.35) implies that both ${\mathcal {X}}$ and $\varPsi$ have the same Floquet exponent. This is because the multiplication factor $\sqrt {a_2}/B_0$ is periodic in nature and therefore does not change the Floquet exponent.

Leveraging the periodicity of $\omega ^2$, $\omega ^2(\ell +L)=\omega ^2(\ell )$, we can Fourier expand

(3.37)\begin{equation} \omega^2(\ell)=\sum_{k\in \mathbb{Z}}\varOmega_k\,{\rm e}^{{\rm i}\,\ell({2{\rm \pi}}/{L})k}, \end{equation}

where $\{\varOmega _k\}_{k\in \mathbb {Z}}$ denote the Fourier coefficients of $\omega ^2$. Let the fundamental solutions of (3.36a,b) be given by $\varPsi _{\pm }(\ell )$ such that

(3.38a-d)\begin{equation} \varPsi_+(0)=1, \quad \varPsi_{-}(0)=0, \quad \varPsi'_+(0)=0, \quad \varPsi'_{-}(0)=1, \end{equation}

where these conditions have been chosen so that the basis functions $\varPsi _{\pm }$ are linearly independent and with unit Wronskian. The Floquet solutions are given by

(3.39)\begin{equation} \left. \begin{aligned} \varPsi_+(\ell) & = {\rm e}^{+{\rm i}\, \nu \ell}\sigma_+(\ell), \quad \sigma_+(\ell+L)=\sigma_+(\ell), \quad \sigma_+(0)=1,\\ \varPsi_{-}(\ell) & = {\rm e}^{-{\rm i}\, \nu \ell}\sigma_{-}(\ell), \quad \sigma_{-}(\ell+L)=\sigma_{-}(\ell), \quad \sigma_-(0)=0, \end{aligned} \right\} \end{equation}

where $\sigma _+$ and $\sigma _{-}$ satisfy the periodicity condition and the linear independence, and $\nu$ denotes the Floquet exponents for which we seek. Nothing more needs to be known about the functions $\sigma _{\pm }$ since their coefficients do not appear in the final expression for the Floquet exponents. Let us now Fourier expand the Floquet solution as

(3.40)\begin{equation} \varPsi_+{=} {\rm e}^{{\rm i}\, \nu \ell}\sum_{n\in \mathbb{Z}}b_n\, {\rm e}^{{\rm i}\,({2{\rm \pi} n }/{L})\ell} = \sum_{n\in \mathbb{Z}} b_n \exp\left({{\rm i}\, \ell \left(\frac{2{\rm \pi} n }{L}+ \nu\right)}\right), \end{equation}

where $\{b_n\}_{n\in \mathbb {Z}}$ is the set of Fourier coefficients of $\sigma _+$. The Fourier expansion of (3.36a,b) reads

(3.41)\begin{align} & \varPsi_+^{''}(\ell)+\omega^2(\ell)\varPsi_+(\ell)={-}\sum_{n\in \mathbb{Z}}b_n \left( \nu + \frac{2{\rm \pi} n}{L}\right)^2 \exp\left({{\rm i}\, \ell \left(\nu + \frac{2n{\rm \pi}}{L}\right)}\right)\nonumber\\ & \qquad + \left(\sum_{k\in \mathbb{Z}}\varOmega_k\,{\rm e}^{{\rm i}\,\ell({2{\rm \pi}}/{L})k} \right)\left(\sum_{n\in \mathbb{Z}} b_n \exp\left({{\rm i}\, \ell \left(\frac{2{\rm \pi} n }{L}+ \nu\right)}\right)\right)\nonumber\\ & \quad =\sum_{n\in \mathbb{Z}}\left(\sum_{k}\varOmega_kb_{n-k}-\left(\frac{2{\rm \pi} n }{L}+ \nu\right)^2b_n\right) \exp\left({{\rm i}\,\ell\left(\frac{2{\rm \pi} n }{L}+ \nu\right)}\right)\nonumber\\ & \quad =0. \end{align}

In order for (3.41) to be verified, the following condition has to hold:

(3.42)\begin{equation} \sum_{k\in \mathbb{Z}}\varOmega_kb_{n-k}-\left(\frac{2{\rm \pi} n }{L}+ \nu\right)^2b_n=0,\quad n \in \mathbb{Z}, \end{equation}

which can be rewritten in the following matrix form, after dividing (3.42) by $\varOmega _0-(2{\rm \pi} n/L+\nu )^2$:

(3.43)\begin{equation} \sum_{m\in \mathbb{Z}}B_{nm}b_m=0, \quad n \in \mathbb{Z}, \end{equation}

where the matrix $B$ is defined as

(3.44)\begin{equation} B_{nn}=1, \quad B_{nm}=\frac{\varOmega_{n-m}}{\varOmega_0-\left(\dfrac{2{\rm \pi} n}{L}+\nu\right)^2}. \end{equation}

This yields the following determinant equation:

(3.45)\begin{equation} \det{B_{nm}} = 0. \end{equation}

The above determinant can be considered a function of the Floquet exponents $\nu$:

(3.46)\begin{equation} \det{B_{nm}(\nu)}:=\Delta(\nu)\equiv \det{\left(\delta_{nm} + \frac{\varOmega_{n-m}}{\varOmega_0-\left(\dfrac{2{\rm \pi} n}{L}+\nu\right)^2}\right)}, \quad n,m \in \mathbb{Z}. \end{equation}

For (3.45) to be verified, we seek for the Floquet exponents $\nu$ such that $\Delta (\nu )=0$. Following Wang, Guo & Xia (Reference Wang, Guo and Xia1989), we now show that the exponents can be deduced from the very simple expression (3.50). One notes that $\Delta (\nu )$ is a $2{\rm \pi}$-periodic function, with poles in $\nu = \pm \sqrt {\varOmega _0} - 2{\rm \pi} n/L$, $n \in \mathbb {Z}$. It can be shown that the determinant is absolutely convergent in the whole $\nu$-plane, except in these poles. Therefore, $\Delta$ is meromorphic. Moreover, as the imaginary part $\mathrm {Im}{(\nu )}\rightarrow \pm \infty$, $\Delta (\nu )\rightarrow 1$. Let us now define the complex function $f$ as

(3.47)\begin{equation} f(\nu):= \cot{\frac{L}{2}(\nu - \sqrt{\varOmega_0})} - \cot{\frac{L}{2}(\nu + \sqrt{\varOmega_0})}. \end{equation}

It is useful to introduce $f$ since it has the same poles and periodicity as $\Delta$. Moreover, it is bounded as $\mathrm {Im}{(\nu )}\rightarrow \pm \infty$. This way, there has to exist a constant $K\in \mathbb {C}$ such that

(3.48)\begin{equation} D(\nu) \equiv \Delta(\nu) + Kf(\nu) \end{equation}

has no singularity in the whole $\nu$ plane. Together with the fact that it is bounded as $|\nu |\rightarrow \infty$, according to Liouville's theorem, $D$ is a constant function. In the limit $|\nu |\rightarrow \infty$, we see that $D=1$. To determine $K$, take $\nu =0$:

(3.49)\begin{equation} \nu =0 \implies K = \frac{1-\Delta(0)}{f(0)} = \frac{1-\Delta(0)}{2\cot{\dfrac{\sqrt{\varOmega_0}L}{2}}}.\end{equation}

Using the value of $K$ from (3.49), the Floquet exponent equation reduces to

(3.50)\begin{equation} \sin^2{\nu\frac{L}{2}} = \Delta(0)\sin^2{\frac{\sqrt{\varOmega_0}L}{2}}. \end{equation}

Although (3.50) is very simple, one obstacle remains to compute the exponents: one has to evaluate the infinite determinant, $\Delta (0)$. We refer to Wang et al. (Reference Wang, Guo and Xia1989) for some approximations of $\Delta (0)$. For instance, when $\varOmega _n$ are sufficiently small, $\Delta (0)$ can be approximated by the order-3 determinant with $B_{00}$ as central element (we recall that the determinant involves summation over all $\mathbb {Z}$), providing

(3.51)\begin{equation} \Delta(0)\simeq 1+\frac{2 \varOmega_1^2}{\varOmega_0\left( 4 - \varOmega_0\right)^2}+\frac{2 \varOmega_1^2 \varOmega_2}{\varOmega_0\left( 4 - \varOmega_0\right)^2}-\frac{\varOmega_2^2}{\left( 4 - \varOmega_0\right)^2}. \end{equation}

However, this approximation breaks down when the coefficients become too large as $|n|$ increases. As of the exponents computed, they might not all be suitable for $\iota$. One has to discard the results that are not relevant. Moreover, we emphasize that the exponents may be shifted by $2k {\rm \pi}/L$, with $k$ an integer, without changing the mathematics of the system, by periodicity. However, one has to carefully choose the adapted value for $\iota$, by setting the appropriate phase shift – see (1.3).

In a nutshell, the second variation of the magnetic action, at an elliptical axis, has been expanded in the Solov'ev–Shafranov near-axis formalism. It has been shown that the null eigenspace of this operator yields the correct on-axis rotational transform, upon applying Floquet theory to solve for the latter. On the other hand, solving for the null eigenspace has led to a system that could be rewritten in the form of a Hill equation. The Floquet exponents appear as solutions of an infinite determinant equation, that can only be approximated analytically. The derivation of the on-axis rotational transform as a Floquet exponent of the null eigenvectors of $\delta ^2{\mathcal {S}}$ has also been carried on in Mercier coordinates, to support the previous conclusions – see Appendix D.

Magnetic confinement device design needs fast and accurate computation of $\iota$. Currently, the most widely used method to compute $\iota$ is field-line tracing (Todoroki Reference Todoroki2003). Field-line tracing methods compute $\iota$ as an infinite time limit of following a field line. In practice, it can be done by following a field line sufficiently long to achieve convergence. As an alternative, one can also compute Greene's residue of the axis, which is related to $\iota$ (see Greene Reference Greene1979; Hanson & Cary Reference Hanson and Cary1984; Hudson Reference Hudson2004). Even though it still involves ODE integration, we can limit ourselves to just one circuit around the torus.

In the following section, we introduce a new method to compute $\iota$ as a solution to a discrete problem involving the magnetic action. The action is discretized in a way similar to that in Mackay & Meiss (Reference Mackay and Meiss1983) and Hudson & Suzuki (Reference Hudson and Suzuki2014).

4 Discrete formalism: piecewise action

We have seen that the problem described in the previous section, involving the second-order variation of the magnetic action to determine the on-axis Floquet exponent, leads to the operator equation (3.5). In another approach, the magnetic axis can be discretized, and the rotational transform can be determined from the multipliers of the latter curve. The multipliers have been shown to be linked to the residue of the curve (see Greene Reference Greene1979; Mackay & Meiss Reference Mackay and Meiss1983). This discrete approach is explored in what follows. We consider

(4.1)\begin{align} \mathcal{S} & = \sum_{i=1}^{n-1} \int_{\mathcal{C}_i} {\boldsymbol{A}} \boldsymbol{\cdot} {\rm d}\boldsymbol{\ell}\nonumber\\ & = \sum_{i=1}^{n-1} S({\boldsymbol{x}}_{i},{\boldsymbol{x}}_{i+1}), \end{align}

meaning that the curve has been discretized with $n\in \mathbb {N}$ points. So far, the way that ${\mathcal {C}}_i$ are defined is not important. Only the endpoints matter. The particular case where ${\mathcal {C}}_i$ are segments is given in Appendix E. The following notation is adopted for the derivatives:

(4.2)\begin{equation} \left. \begin{aligned} \boldsymbol{S}_1^{[i,i+1]} & := \boldsymbol{\nabla}_{{\boldsymbol{x}}_i} S^{[i,i+1]}=\boldsymbol{\nabla}_{{\boldsymbol{x}}_i} S({\boldsymbol{x}}_i, {\boldsymbol{x}}_{i+1}),\\ \boldsymbol{S}_2^{[i,i+1]} & := \boldsymbol{\nabla}_{{\boldsymbol{x}}_{i+1}} S^{[i,i+1]}=\boldsymbol{\nabla}_{{\boldsymbol{x}}_{i+1}} S({\boldsymbol{x}}_i, {\boldsymbol{x}}_{i+1}), \end{aligned} \right\} \end{equation}

as well as for the second-order derivatives:

(4.3)\begin{equation} \boldsymbol{\mathsf{S}}_{21}^{[i,i+1]} := \boldsymbol{\nabla}_{{\boldsymbol{x}}_{i+1}}\boldsymbol{S}_1^{[i,i+1]}, \end{equation}

and similarly for $\boldsymbol{\mathsf{S}}_{12}$, $\boldsymbol{\mathsf{S}}_{11}$ and $\boldsymbol{\mathsf{S}}_{22}$. Let us also define generalized periodic orbits of type $(q)$, as orbits with

(4.4)\begin{equation} {\boldsymbol{x}}_{i+q} = {\boldsymbol{x}}_i, \end{equation}

for some $q\in \mathbb {N}$. Therefore, the magnetic axis as discretized above is a general periodic orbit of type $(n)$. The terminology ‘orbit’ for field lines is justified by the Hamiltonian behaviour of the magnetic field. Moreover, the discretized action (4.1) satisfies the periodicity condition

(4.5)\begin{equation} \sum_{i=1}^{n-1} S({\boldsymbol{x}}_{i+n},{\boldsymbol{x}}_{i+n+1}) = \sum_{i=1}^{n-1} S({\boldsymbol{x}}_{i},{\boldsymbol{x}}_{i+1}) + C, \end{equation}

where the constant $C=0$. In the one-dimensional case, as dealt with in Mackay & Meiss (Reference Mackay and Meiss1983), the existence of extrema of the action is ensured by an additional convexity condition on the Lagrangian, $-L_{12}>0$. In our three-dimensional case, this can be generalized in that the second-order derivative tensors $\boldsymbol{\mathsf{S}}_{12}^{[i,i+1]}$ ought to be negative definite for all $i$:

(4.6)\begin{equation} {\boldsymbol{x}}^{\mathsf{T}}{\boldsymbol{\mathsf{S}}_{12}}^{[i,i+1]}{\boldsymbol{x}} <0, \quad {\boldsymbol{x}} \in \mathbb{R}^{*3}, \quad 1\leq i \leq n.\end{equation}

The fact that $C=0$ together with the condition that $\boldsymbol{\mathsf{S}}_{12}^{[i,i+1]}$ guarantees that the action of periodic orbits of type $(n)$ is bounded below and the Poincaré–Birkhoff theorem (Birkhoff Reference Birkhoff1913) ensures the existence of at least two stationary trajectories among the space of all generalized periodic paths of type $(n)$, one that minimizes the action and is elliptic, and one that is a saddle called minimax and is usually hyperbolic but can be alternating hyperbolic, where hyperbolic and elliptic describe the behaviour of nearby trajectories.

For an $(n)$-periodic curve that extremizes the action, the latter has to be stationary with respect to an arbitrary variation in its geometry $\delta {\boldsymbol {x}}_i$, and the stationarity condition can be expressed in terms of the previously defined derivatives:

(4.7)\begin{equation} \left. \begin{gathered} \delta S [\delta {\boldsymbol{x}}_i] = \left[\boldsymbol{\nabla}_{{\boldsymbol{x}}_i} S^{[i-1,i]}+ \boldsymbol{\nabla}_{{\boldsymbol{x}}_i} S^{[i,i+1]} \right]\cdot \delta {\boldsymbol{x}}_i=0,\\ \Leftrightarrow\quad \boldsymbol{S}_2^{[i-1,i]} + \boldsymbol{S}_1^{[i,i+1]} = \boldsymbol{0}. \end{gathered} \right\} \end{equation}

Note that (4.7) holds for any $1\leq i \leq n$, and the stationarity condition is expressed for each point in terms of the two nearest neighbours. This way, a magnetic axis can be identified. Similarly to the continuous case dealt with in § 3, the neighbouring field lines have to satisfy

(4.8)\begin{equation} \boldsymbol{\nabla}_{{\boldsymbol{x}}_{i-1}} \boldsymbol{S}_2^{[i-1,i]} \delta {\boldsymbol{x}}_{i-1} + \boldsymbol{\nabla}_{{\boldsymbol{x}}_{i+1}}\boldsymbol{S}_1^{[i,i+1]} \delta {\boldsymbol{x}}_{i+1} + \boldsymbol{\nabla}_{{\boldsymbol{x}}_{i}}\left(\boldsymbol{S}_2^{[i-1,i]} + \boldsymbol{S}_1^{[i,i+1]}\right) \delta {\boldsymbol{x}}_{i} = \boldsymbol{0}, \end{equation}

which we rewrite in terms of the second-order derivatives of the action as

(4.9)\begin{equation} \boldsymbol{\mathsf{S}}_{12}^{[i-1,i]} \delta {\boldsymbol{x}}_{i-1} + \boldsymbol{\mathsf{S}}_{21}^{[i,i+1]} \delta {\boldsymbol{x}}_{i+1} + \left( \boldsymbol{\mathsf{S}}_{22}^{[i-1,i]} + \boldsymbol{\mathsf{S}}_{11}^{[i,i+1]}\right) \delta {\boldsymbol{x}}_{i} = \boldsymbol{0}, \quad 1\leq i\leq n. \end{equation}

From Mackay & Meiss (Reference Mackay and Meiss1983), we know that the multipliers $\lambda$ of an ($n$)-periodic orbit are defined by the existence of a tangent orbit satisfying

(4.10)\begin{equation} \delta {\boldsymbol{x}}_{i+n} = \lambda \delta {\boldsymbol{x}}_i, \end{equation}

so the following holds:

(4.11a,b)\begin{equation} \delta {\boldsymbol{x}}_{n+1} = \lambda \delta {\boldsymbol{x}}_{1}, \quad \delta {\boldsymbol{x}}_0 = \lambda^{{-}1}\delta {\boldsymbol{x}}_n. \end{equation}

They can be written in their Floquet form (Greene Reference Greene1979):

(4.12)\begin{equation} \lambda = {\rm e}^{{\rm i}\,\nu}, \end{equation}

where the Floquet exponent $\nu$ describes the average rotation angle per period of the ($n$)-orbit, therefore, the rotational transform $\iota$. Since (4.9) is valid for any $1\leq i\leq n$, together with (4.11a,b), it can be rewritten in tensor form:

(4.13)\begin{equation} \left(\begin{array}{cccc}\left( \boldsymbol{\mathsf{S}}_{22}^{[01]} + \boldsymbol{\mathsf{S}}_{11}^{[12]} \right) & \boldsymbol{\mathsf{S}}_{12}^{[12]} & & \lambda^{{-}1}\boldsymbol{\mathsf{S}}_{21}^{[01]}\\ \boldsymbol{\mathsf{S}}_{21}^{[12]} & \boldsymbol{\mathsf{S}}_{12}^{[23]} & \ddots & \\ & \ddots & \ddots & \boldsymbol{\mathsf{S}}_{12}^{[n-1,n]} \\ \lambda\boldsymbol{\mathsf{S}}_{12}^{[n,n+1]} & & \boldsymbol{\mathsf{S}}_{21}^{[n-1,n]} & \left( \boldsymbol{\mathsf{S}}_{22}^{[n-1,n]} + \boldsymbol{\mathsf{S}}_{11}^{[n,n+1]} \right)\end{array}\right) \left(\begin{array}{c} \delta {\boldsymbol{x}}_{1}\\ \vdots \\ \\ \vdots \\ \delta {\boldsymbol{x}}_{n} \end{array}\right) = \boldsymbol{0}, \end{equation}

where the block-tridiagonal form with corners arises naturally. Denoting by $\boldsymbol{\mathsf{M}}$ the matrix of second derivatives, we get an equation for the multipliers:

(4.14)\begin{equation} \boldsymbol{\mathsf{M}}(\lambda)\delta {\boldsymbol{x}} = \boldsymbol{0}. \end{equation}

The blank spaces in $\boldsymbol{\mathsf{M}}$ are blocks of $0$. Note that each block in $\boldsymbol{\mathsf{M}}$ is of size $3\times 3$ since $\boldsymbol{\mathsf{S}}_{lm}$ contain the second derivatives of a line integral embedded in three-dimensional space. Thus, for an $(n)$-periodic curve, $\boldsymbol{\mathsf{M}}\in \mathcal {M}_{3n\times 3n}$, and $\delta {\boldsymbol {x}} \in \mathbb {R}^{3n}$. For (4.14) to hold, the determinant of $\boldsymbol{\mathsf{M}}$ ought to be zero to avoid the trivial solution $\delta {\boldsymbol {x}} = \boldsymbol {0}$.

For a block-tridiagonal matrix defined as in (4.13) with $\lambda \in \mathbb {C}$, the determinant can be expressed analytically (Molinari Reference Molinari2008):

(4.15)\begin{equation} \left. \begin{aligned} \det \boldsymbol{\mathsf{M}}(\lambda) & = \frac{({-}1)^{3n}}{(-\lambda)^3}\det\left(\boldsymbol{\mathsf{T}}_S-\lambda \boldsymbol{\mathsf{I}}_{6} \right)\det\left( \prod_{i=1}^{n}\boldsymbol{\mathsf{S}}_{12}[i,i+1]\right),\\ \boldsymbol{\mathsf{T}}_S & = \prod_{i=1}^{n}\left(\begin{array}{cc} -{\boldsymbol{\mathsf{S}}^{{-}1}_{12}}^{[i,i+1]}(\boldsymbol{\mathsf{S}}_{22}^{[i-1,i]}+\boldsymbol{\mathsf{S}}_{11}^{[i,i+1]}) & -{\boldsymbol{\mathsf{S}}_{12}^{{-}1}}^{[i,i+1]}\boldsymbol{\mathsf{S}}_{12}^{[i-1,i]} \\ \boldsymbol{\mathsf{I}}_{3} & \boldsymbol{\mathsf{0}}\end{array}\right), \end{aligned} \right\} \end{equation}

which requires one to define the so-called transfer matrix $\boldsymbol{\mathsf{T}}_S$. Our generalized convexity condition (4.6) that ensures the existence of minimizing orbits guarantees that the determinant of the transfer matrix can be computed as $\boldsymbol{\mathsf{S}}_{12}^{[i,i+1]}$ is invertible $\forall i$. The multipliers $\lambda$ are then given by the solutions of

(4.16)\begin{equation} \lambda^{{-}3} \det\left(\boldsymbol{\mathsf{T}}_S-\lambda \boldsymbol{\mathsf{I}} \right) = 0, \end{equation}

so they are the non-zero eigenvalues of $\boldsymbol{\mathsf{T}}_S$. Once the multipliers have been computed, the Floquet exponent can be easily determined from (3.31).

Alternatively, note that (4.13) gives a recursive relation for $\delta {\boldsymbol {x}}_i$:

(4.17)\begin{equation} \left. \begin{gathered} \left(\boldsymbol{\mathsf{S}}_{22}^{[01]} + \boldsymbol{\mathsf{S}}_{11}^{[12]} \right) \delta {\boldsymbol{x}}_1 + \boldsymbol{\mathsf{S}}_{12}^{[12]}\delta{\boldsymbol{x}}_2 + \lambda^{{-}1}\boldsymbol{\mathsf{S}}_{21}^{[01]}\delta{\boldsymbol{x}}_n =0,\\ \lambda\boldsymbol{\mathsf{S}}_{12}^{[n,n+1]}\delta{\boldsymbol{x}}_1 + \boldsymbol{\mathsf{S}}_{21}^{[n-1,n]}\delta{\boldsymbol{x}}_{n-1} + \left(\boldsymbol{\mathsf{S}}_{22}^{[n-1,n]} + \boldsymbol{\mathsf{S}}_{11}^{[n,n+1]} \right) \delta {\boldsymbol{x}}_n = 0,\\ \boldsymbol{\mathsf{S}}_{21}^{[k-1,k]}\delta{\boldsymbol{x}}_{k-1} + \left(\boldsymbol{\mathsf{S}}_{22}^{[k-1,k]} + \boldsymbol{\mathsf{S}}_{11}^{[k,k+1]} \right) \delta {\boldsymbol{x}}_k + \boldsymbol{\mathsf{S}}_{12}^{[k, k+1]}{\boldsymbol{x}}_{k+1} =0; \quad 2\leq k\leq n-1. \end{gathered} \right\} \end{equation}

Equation (4.17) can be rewritten as

(4.18)\begin{equation} \left(\begin{array}{c} \delta {\boldsymbol{x}}_{n+1} \\ \delta {\boldsymbol{x}}_{n} \end{array}\right) = \boldsymbol{\mathsf{J}} \left(\begin{array}{c} \delta {\boldsymbol{x}}_{1} \\ \delta {\boldsymbol{x}}_{0}\end{array}\right), \end{equation}

with

(4.19)\begin{equation} \boldsymbol{\mathsf{J}} := \prod_{k=1}^{n}\left(\begin{array}{cc} -{\boldsymbol{\mathsf{S}}^{{-}1}_{12}}^{[k,k+1]}(\boldsymbol{\mathsf{S}}_{22}^{[k-1,k]}+\boldsymbol{\mathsf{S}}_{11}^{[k,k+1]}) & -{\boldsymbol{\mathsf{S}}_{12}^{{-}1}}^{[k,k+1]}\boldsymbol{\mathsf{S}}_{12}^{[k-1,k]} \\ \boldsymbol{\mathsf{I}} & \boldsymbol{\mathsf{0}}\end{array}\right). \end{equation}

The multipliers of a periodic orbit are the eigenvalues of the derivatives of the return map around the orbit, $\boldsymbol{\mathsf{J}}$ (Mackay & Meiss Reference Mackay and Meiss1983), so upon comparison with (4.15), $\boldsymbol{\mathsf{J}}=\boldsymbol{\mathsf{T}}_S$ confirming the result from (4.16).

The strength of this method resides in the fact that the problem is reduced to finding the eigenvalues of a $6\times 6$ matrix. In fact, the size of the operator matrix $\boldsymbol{\mathsf{M}}$ increases linearly with the number of discretization points, but of interest are solely the eigenvalues of the matrix $\boldsymbol{\mathsf{T}}_S$, whose size is $6\times 6$, no matter how many discretization points have been used. Although the number of matrices in the product $\boldsymbol{\mathsf{J}}$ also increases linearly with the number of discretization points, the most commonly used algorithms to determine the eigenvalues of a matrix have a cost that is generally cubic (Francis Reference Francis1961, Reference Francis1962; Kublanovskaya Reference Kublanovskaya1962). Since an $m\times m$ matrix has at most $m$ distinct eigenvalues, reducing its size is crucial for computational efficiency. The derivation of the aforementioned results in the case where the curve of interest is discretized by piecewise linears is given in Appendix E.