3.1. The coordinate system

Our coordinate system $\{s,\theta ,\zeta \}=\{q^1,q^2,q^3\}$ is defined (inversely) by the position vector $\boldsymbol{r}=\boldsymbol{r}(s,\theta ,\zeta )$ . In the curvilinear coordinates, the tangent basis vectors are $\displaystyle {\boldsymbol{e}_i={\partial \boldsymbol{r}}/{\partial q^i}}$ and the reciprocal basis vectors are $\boldsymbol{e}^i=\boldsymbol{\nabla } {q^i}$ . For any vector $\boldsymbol{w}$ , superscripts denote the components in the expansion over the tangent basis, $\boldsymbol{w}=w^i\boldsymbol{e}_i$ with $w^i=\boldsymbol{w}\boldsymbol{\cdot }\boldsymbol{e}^i$ (contravariant components), and subscripts denote the components in the expansion over the reciprocal basis, $\boldsymbol{w}=w_i\boldsymbol{e}^i$ with $w_i=\boldsymbol{w}\boldsymbol{\cdot }\boldsymbol{e}_i$ (covariant components). The components of the metric tensor are $g_{ij}=\boldsymbol{e}_i\boldsymbol{\cdot }\boldsymbol{e}_j$ , and its inverse has components $g^{ij}=\boldsymbol{e}^i\boldsymbol{\cdot }\boldsymbol{e}^j$ . The Jacobian is therefore $\sqrt {g}=\boldsymbol{e}_1\boldsymbol{\cdot }(\boldsymbol{e}_2\times \boldsymbol{e}_3)=\boldsymbol{e}^1\boldsymbol{\cdot }(\boldsymbol{e}^2\times \boldsymbol{e}^3)$ . In this paper, we assume $\zeta$ denotes the direction of the symmetry; i.e. all the variable fields are independent of $\zeta$ (the 2-D assumption). Using a general coordinate system, we encode the geometry in the metric tensor. We shall assume that the $\{s,\theta ,\zeta \}$ coordinates, as well as the position vector $\boldsymbol{r}$ , are dimensionless; otherwise, the subscripted and superscripted components might not carry the physical dimension of their corresponding vector quantities (Truesdell Reference Truesdell1953). In this paper we assume $s\in [0,1]$ , $\theta \in [0,2\pi )$ and $\zeta \in [0,2\pi )$ . We note that these coordinates are not necessarily straight field line coordinates.

To find the position vector (and therefore the coordinate transformations), we first define the position of the two enclosing boundaries of the domain, labelled by $s=0$ and $ s=1$ , as $\boldsymbol{\mathcal{R}}_-(s,\theta )$ and $\boldsymbol{\mathcal{R}}_+(s,\theta )$ . The position vector is then extended from these boundaries to the interior of the domain by linear interpolation; i.e. by defining $\boldsymbol{r}=(1-s)\boldsymbol{\mathcal{R}}_- +s\boldsymbol{\mathcal{R}}_+$ .

The magnetic field lines are defined as the solution to the differential equation $({\mathrm{d} q^i}/{\mathrm{d} \tau _0})=B^i$ , where $\tau _0$ is a time-like ordering parameter. The Poincaré map can be obtained by solving the field line equations and collecting the points at which orbits intersect the $\zeta =0$ plane. However, assuming that $B^\zeta$ is never zero, it is easier to integrate the alternative equations

(3.1) \begin{align} \frac {\mathrm{d} s}{\mathrm{d} \zeta }&=\frac {B^s}{B^\zeta },\\[-10pt]\nonumber \end{align}

(3.2) \begin{align} \frac {\mathrm{d} \theta }{\mathrm{d} \zeta }&=\frac {B^\theta }{B^\zeta }, \end{align}

to obtain the Poincaré map on $\{s,\theta \}$ plane. The 2-D assumption $\partial _\zeta = 0$ guarantees that the field lines described by these equations are integrable. The 2-D assumption also restricts all boundary perturbations to have toroidal mode number $n = 0$ . Consequently, the only possible magnetic resonance occurs at $\iota = 0$ . Higher- $m$ harmonics modify the structure of this resonant layer but do not generate additional island chains. We note that because $\{s,\theta ,\zeta \}$ are not straight-field-line coordinates, the rotational transform can be a function of both $s$ and $\theta$ . In the following numerical calculations, we mainly report the $\theta$ -averaged rotational transform.

3.2. Generalised Beltrami equation

Equation (2.8d ) is the generalised Beltrami equation. To solve this equation, we start by representing the divergence-free magnetic field as $\boldsymbol{B}=\boldsymbol{\nabla }\times {\boldsymbol{A}}$ where $\boldsymbol{A}$ is the vector potential. Using the gauge freedom, we can assume $\boldsymbol{A}=\bar {\boldsymbol{A}}+\boldsymbol{\nabla } {G}$ so that

(3.3a) \begin{align} A_s(s,\theta ,\zeta )=\bar {A_s}(s,\theta ,\zeta )+\frac {\partial G}{\partial s}\;(s,\theta ,\zeta )&=0,\\[-10pt]\nonumber \end{align}

(3.3b) \begin{align} A_\theta (1,\theta ,\zeta )=\bar {A_\theta }(1,\theta ,\zeta )+\frac {\partial G}{\partial \theta }\,(1,\theta ,\zeta )&=\frac {\psi _\zeta }{2\pi },\\[-10pt]\nonumber \end{align}

(3.3c) \begin{align} A_\zeta (1,0,\zeta )=\bar {A_\zeta }(1,0,\zeta )+\frac {\partial G}{\partial \zeta }\,(1,0,\zeta )&=-\frac {\psi _\theta }{2\pi }, \end{align}

where $1$ and $0$ are the $s$ coordinates of the two enclosing boundaries, respectively. The quantities $\psi _\theta$ and $\psi _\zeta$ are two constants that are provided as inputs. Equations (3.3b ) and (3.3c ) provide boundary conditions for $A_\theta$ and $A_\zeta$ at $s=1$ boundary. In addition, in two dimensions ( $({\partial}/({\partial \zeta })) = 0$ ) the $\boldsymbol{B}\boldsymbol{\cdot }\hat {\boldsymbol{n}}=\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla } {s}$ reads

(3.4) \begin{align} \frac {\partial\kern-1.1pt A_\zeta }{\partial \theta } (s=0)=\frac {\partial\kern-1.1pt A_\zeta }{\partial \theta } (s=1)=0. \end{align}

Therefore, $A_\zeta$ must be constant on both $s=0$ and $s=1$ boundaries. From (3.3c ) for the $s=1$ boundary this constant is $-\psi _\theta /(2\pi )$ . For the $s=0$ boundary, we use the boundary conditions

(3.5a) \begin{align} \int _0^{2\pi } A_\theta (0,\theta )\mathrm{d}\theta &=0,\\[-10pt]\nonumber \end{align}

(3.5b) \begin{align} A_\zeta (0,\theta )&=0. \end{align}

Using the boundary values of $\boldsymbol{A}$ , one can show that $\psi _\theta$ and $\psi _\zeta$ are the total enclosed poloidal and toroidal fluxes, respectively; i.e.

(3.6a) \begin{align} \int _{0}^{1}\int _0^{2\pi }B^\theta (s,0,\zeta )\sqrt {g}\mathrm{d} s \mathrm{d} \zeta &=\int _0^{2\pi } (-A_\zeta (s=1,0)+A_\zeta (s=0,0))\mathrm{d}\zeta =\psi _\theta,\\[-10pt]\nonumber \end{align}

(3.6b) \begin{align} \int _{0}^{1}\int _0^{2\pi }B^\zeta (s,\theta ,0)\sqrt {g}\mathrm{d} s \mathrm{d} \theta &=\int _0^{2\pi } (A_\theta (s=1,\theta )-A_\theta (s=0,\theta ))\mathrm{d} \theta =\psi _\zeta. \end{align}

Using the gauge condition (3.3a ), the $s$ and $\theta$ components of (2.8d ) read

(3.7a) \begin{align} -\mu _B \frac {\partial\kern-1.1pt A_\zeta }{\partial \theta }-\nu \frac {\partial u_\zeta }{\partial \theta }+\frac {\partial\kern-1.4pt B_\zeta }{\partial \theta }&=0,\\[-10pt]\nonumber \end{align}

(3.7b) \begin{align} \mu _B \frac {\partial\kern-1.1pt A_\zeta }{\partial s}+\nu \frac {\partial u_\zeta }{\partial s}-\frac {\partial\kern-1.4pt B_\zeta }{\partial s}&=0, \end{align}

where $B_\zeta =g_{31}B^s+g_{32}B^\theta +g_{33}B^\zeta$ . These two equations can be combined in

(3.8) \begin{align} \mu _BA_\zeta +\nu u_\zeta -B_\zeta =\mathcal{C}, \end{align}

where $\mathcal{C}$ is a constant determined by the boundary conditions. Using (2.8b ), (3.8) can be written as

(3.9) \begin{align} \mu _BA_\zeta +\nu v_\zeta +\big[\big(M^{Rx}_\parallel \big)^2-1\big] B_\zeta =\mathcal{C}, \end{align}

where $M^{Rx}_\parallel$ is the parallel Alfvén Mach number defined by

(3.10) \begin{align} \big(M^{Rx}_\parallel \big)^2\equiv \frac {(u^{Rx})^2}{c_A^2}=\frac {\nu ^2}{\rho }, \end{align}

and $c_A\equiv ({B^2/\rho })^{1/2}$ is the Alfvén speed. Equation (3.8) should be solved in conjunction with the $\zeta$ component of (2.8d ):

(3.11) \begin{align} -\mu _B B^\zeta -\frac {\nu }{\sqrt {g}} \left (\frac {\partial u_\theta }{\partial s}-\frac {\partial u_s}{\partial \theta }\right )+\frac {1}{\sqrt {g}} \left (\frac {\partial\kern-1.4pt B_\theta }{\partial s}-\frac {\partial\kern-1.4pt B_s}{\partial \theta }\right )=0. \end{align}

3.3. The equation of motion

Neoclassical considerations motivate the study of flow that is approximately tangent to contours of $B$ , i.e. $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }B \approx 0$ (Shaing & Callen Reference Shaing and Callen1982; Hinton & Wong Reference Hinton and Wong1985; Helander Reference Helander2007). In 2-D equilibria where $B$ varies predominantly across $(s,\theta )$ , this suggests that the dominant component of the flow lies along the symmetry direction ( $\zeta$ ). To achieve this, we prescribe $\boldsymbol{v}$ to be

(3.12) \begin{align} \boldsymbol{v}=v^\zeta \boldsymbol{e}_{\zeta }, \end{align}

so that the only non-zero component of $\boldsymbol{v}$ is $v^\zeta$ . In two dimensions, this assumption also implies $\boldsymbol{\nabla }\boldsymbol{\cdot }{\rho \boldsymbol{u}}=({{\partial (\sqrt {g}\rho v^\zeta )}/{\partial \zeta }})=0$ ; i.e. the continuity (2.8e ) is satisfied. Using (2.8b ), $\boldsymbol{u}$ can be decomposed into the poloidal and toroidal components

(3.13) \begin{align} \boldsymbol{u}=\frac {\nu }{\rho }(B^s\boldsymbol{e}_s+B^\theta \boldsymbol{e}_\theta )+\left(\frac {\nu B^\zeta }{\rho }+v^\zeta \right)\boldsymbol{e}_\zeta. \end{align}

Therefore, like the magnetic field, the total flow $\boldsymbol{u}$ is generally helical. To be consistent with the expectations of neoclassical transport theory, one can assume $\nu \ll 1$ and $\rho \sim O(1)$ (hence $M^{Rx}_\parallel \ll 1$ ) so that the poloidal flow as defined by (3.13) is much smaller than $v^\zeta$ .

Another useful decomposition of $\boldsymbol{u}$ is into its field-aligned and cross-field components

(3.14) \begin{align} \boldsymbol{u}=\boldsymbol{u}_\parallel +\boldsymbol{u}_\perp = \frac {\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{B}}{B^2}\boldsymbol{B}+ \frac {1}{B^2}\boldsymbol{B}\times ( \boldsymbol{u}\times \boldsymbol{B}). \end{align}

Using (2.8b ) and (3.12),

(3.15a) \begin{align} \boldsymbol{u}_\parallel &=\left(\frac {\nu }{\rho }+\frac {B_\zeta v^\zeta }{B^2}\right)\boldsymbol{B},\\[-10pt]\nonumber \end{align}

(3.15b) \begin{align} \boldsymbol{u}_\perp &=v^\zeta \boldsymbol{e}_\zeta -\frac {B_\zeta v^\zeta }{B^2}\boldsymbol{B}, \end{align}

and the parallel and perpendicular components of the flow are

(3.16a) \begin{align} u_\parallel &=\frac {\nu B^2+\rho v^\zeta B_\zeta }{\rho B},\\[-10pt]\nonumber \end{align}

(3.16b) \begin{align} u_\perp &=\frac {|v^\zeta |}{B}\sqrt {g_{33}B^2-B_\zeta ^2}=\frac {|v^\zeta |}{B}\sqrt {g_{33}B_p^2-g_{31}B^sB_\zeta -g_{32}B^\theta B_\zeta }, \end{align}

where $B_p\equiv \sqrt {B_sB^s+B_\theta B^\theta }$ is the norm of poloidal field. We define the ‘cross-field flow ratio’ as $({u_\perp }/{u_\parallel })$ .

Because the O point of a magnetic island is a fixed point of the Poincaré map, (3.1) and (3.2) imply that $B^s = B^\theta = 0$ at the O point. Consequently, (3.16b ) predicts that the cross-field component of the flow vanishes locally and the flow becomes field-aligned at the island centre. This prediction is consistent with experimental observations in tokamaks and stellarators, which show that cross-field flow vanishes or is strongly suppressed at the O point of low-level magnetic islands (Estrada et al. Reference Estrada, Ascasíbar, Blanco, Cappa, Hidalgo, Ida, López-Fraguas and van Milligen2016; Jiang et al. Reference Jiang2018). Importantly, this behaviour is not imposed in the present RxMHD formulation but arises self-consistently from the variational structure of the model together with the constrained-flow ansatz (3.12). This provides additional justification for the physical consistency of the ansatz adopted here.

Using (3.12) one can show that $(\boldsymbol{\nabla }\times {\boldsymbol{u}})\times \boldsymbol{v}=-v^\zeta \boldsymbol{\nabla } {u_\zeta }$ and therefore (2.8c ) reads

(3.17) \begin{align} v^\zeta \boldsymbol{\nabla } {u_\zeta }=\boldsymbol{\nabla } {h_p}+\frac {\boldsymbol{\nabla } {(u^{\zeta }u_{\zeta })}}{2}, \end{align}

where $h_p\equiv ({(u^s u_s+u_\theta u^\theta)}/{2})+\tau \ln ({{\rho }/{\rho _\varOmega }})$ is the poloidal part of the $h$ . Taking the curl of (3.17), and using (2.8b ) we get

(3.18) \begin{align} \boldsymbol{\nabla } {v^\zeta }\times \boldsymbol{\nabla } {u_\zeta ^{Rx}}+v^\zeta \boldsymbol{\nabla } {v^\zeta }\times \boldsymbol{\nabla } {g_{33}}=0. \end{align}

Equation (3.18) constitutes a solvability condition for the steady-state semirelaxed equilibrium equations. In the absence of time dependence, the system (2.8e ) is generically underdetermined, and not all prescribed constrained flows $\boldsymbol{v}$ lead to admissible equilibria. The solvability condition expresses a necessary compatibility between the constrained flow $v^\zeta$ , the relaxed field-aligned flow $u^{Rx}_\zeta$ , and the geometry through the metric coefficient $g_{33}$ . This result plays a central role in what follows: it dictates which forms of $v^\zeta$ are admissible in a given geometry and explains why flow parameters affect equilibrium structure differently in slab, cylindrical and toroidal configurations.

An alternative form of (3.18) can be calculated by using (2.8b ) and (3.8) to obtain

(3.19) \begin{align} \boldsymbol{\nabla } {u_\zeta ^{Rx}}=-\frac {1}{\nu }\boldsymbol{\nabla } {(\mu _BA_\zeta -B_\zeta +\nu v_\zeta )}. \end{align}

Using this relation, (3.18) reads

(3.20) \begin{align} \boldsymbol{\nabla } {(-\mu _BA_\zeta +B_\zeta )}\times \boldsymbol{\nabla } {v^\zeta }=0. \end{align}

To satisfy this equation, it is sufficient that $v^\zeta$ is an arbitrary function of $-\mu _BA_\zeta +B_\zeta$ . Using (3.8) and (3.12), one can reduce (3.11) to

(3.21) \begin{align} & -\mu _Bg^{31}B_s-\mu _Bg^{32}B_\theta -\mu _Bg^{33}\left (\frac {\mathcal{C}-\mu _BA_\zeta -\nu v^\zeta g_{33}}{\nu ^2/\rho -1}\right )+\frac {1}{\sqrt {g}}\frac {\partial }{\partial s}\left [\frac {(\nu ^2/\rho -1)}{\sqrt {g}}\frac {\partial\kern-1.1pt A_\zeta }{\partial s}\right ]\nonumber \\[7pt]&\quad +\frac {1}{\sqrt {g}}\frac {\partial }{\partial \theta }\left [\frac {(\nu ^2/\rho -1)}{\sqrt {g}}\frac {\partial\kern-1.1pt A_\zeta }{\partial \theta }\right ]-\frac {\nu }{\sqrt {g}}\left (\frac {\partial (g_{32}v^\zeta )}{\partial s}-\frac {\partial (g_{31}v^\zeta )}{\partial \theta }\right )=0. \end{align}

3.4. The numerical methods and convergence plots

In each geometry, (3.8), (3.11) and (3.17) are solved to obtain an equilibrium solution. As shown later, (3.18) allows (3.17) to be reduced to an algebraic equation prior to the numerical analysis. Numerical solutions are obtained using Dedalus v3 (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020), a Python-based framework that enables symbolic specification of partial differential equations. The underlying numerical solver employs tau-spectral methods, which are widely used in magnetic confinement equilibrium studies.

For the 2-D studies presented here, we use Chebyshev expansions in the $s$ direction and Fourier expansions in the $\theta$ direction. To assess the numerical accuracy of our computations, we perform a series of convergence tests. We denote by $N_s$ the number of Chebyshev polynomials used in the $s$ direction and by $N_\theta$ the number of Fourier modes used in the $\theta$ direction. Figure 1 shows the residue of the ideal MHD force-balance equation (2.9) as a function of $N_s$ and $N_\theta$ separately: in figure 1(a) $N_s$ is varied with $N_\theta =64$ fixed, and in figure 1(b) $N_\theta$ is varied with $N_s=64$ fixed. In each case the resolution that is held fixed is large enough that it does not limit the residual, so that each panel isolates the convergence in a single direction. In both directions, the residual initially decreases exponentially over the first four refinement steps, as expected for Chebyshev–Fourier (pseudo-)spectral methods applied to smooth functions (Boyd 2001). Beyond this, the decrease plateaus at a floor of order $10^{-12}$ , where the finite solver tolerance and machine round-off error become the dominant sources of numerical error. The convergence test here is reported only in the toroidal geometry corresponding to the example solution of § 5. We have, however, observed similar convergence in all numerical solutions reported in this study. Although we did not need to use this capability in this study, Dedalus v3 is also a parallel framework designed for high-performance computing. However, it has some limitations that may make it less suitable for large-scale magnetic confinement projects. In particular, it does not support custom coordinate systems of the type used in this work. As a result, we were required to express the governing equations explicitly in terms of scalar partial derivatives and to handle metric components manually. Another limitation of Dedalus v3 is that it relies on direct matrix inversion to solve the linear systems arising from spectral discretisation and Newton iterations, rather than more modern iterative solvers such as GMRES (Saad Reference Saad2003). Consequently, it does not provide any solution for an underdetermined system of equations, as such systems will lead to a rank-deficient matrix that cannot be inverted.

Figure 1.

Numerical convergence of the Dedalus solver in the toroidal geometry. (a) Force-balance residual versus the Chebyshev resolution $N_s$ , with $N_\theta =64$ fixed. (b) Force-balance residual versus the Fourier resolution $N_\theta$ , with $N_s=64$ fixed.

Figure 2.

(a) The HKT slab geometry. (b) Cylindrical geometry. The perturbed boundaries are shown in blue.

4.1. The HKT slab geometry

In the Cartesian coordinate system, the position vector is $\boldsymbol{r}=x\hat {x}+y\hat {y}+z\hat {z}$ , where $\hat {x}$ , $\hat {y}$ and $\hat {z}$ are the usual unit vectors of Cartesian coordinates. The domain of our HKT slab is shown in figure 2. The left-hand boundary and right-hand boundaries, are positioned at $\boldsymbol{\mathcal{R}}_-=\delta \cos {3\theta }\hat {x}+\theta \hat {y}+\zeta \hat {z}$ and $\boldsymbol{\mathcal{R}}_+=(1-\delta \cos {3\theta })\hat {x}+\theta \hat {y}+\zeta \hat {z}$ , respectively. Using linear interpolation, the position vector in the whole domain is $\boldsymbol{r}= [(1-s)\delta \cos {3\theta }+s(1-\delta \cos {3\theta })]\hat {x}+\theta \hat {y}+\zeta \hat {z}$ . Therefore, transformation from the $\{x,y, z\}$ coordinates to the curvilinear coordinates $\{s,\theta ,\zeta \}$ is

(4.1) \begin{align} x&=(1-s)\delta \cos (3\theta )+s(1-\delta \cos {3\theta }), \\[-10pt]\nonumber\end{align}

(4.2) \begin{align} y&=\theta ,\\[-10pt]\nonumber \end{align}

(4.3) \begin{align} z&=\zeta . \end{align}

In the HKT slab geometry, the metric components are $g_{ij}=({\partial x}/{\partial q^i})({\partial x}/{\partial q^j}) + ({\partial y}/{\partial q^i})({\partial y}/{\partial q^j}) + ({\partial z}/{\partial q^i})({\partial z}/{\partial q^j})$ which gives

(4.4) \begin{align} g_{ij}=\left ( \begin{array}{c@{\quad}c@{\quad}c} g_{11} & g_{12} & 0 \\[3pt] g_{21}& g_{22} & 0 \\[3pt] 0 & 0 & 1 \\ \end{array} \right ), \end{align}

where $g_{11}=(1-2 \delta \cos (3 \theta ))^2$ , $g_{12}=g_{21}=3 \delta (2 s-1) (\sin (3 \theta )-\delta \sin (6 \theta ))$ and $g_{22}=9 \delta ^2 (1-2 s)^2 \sin ^2(3 \theta )+1$ . Therefore, we have

(4.5) \begin{align} g_{33}=1, \end{align}

and (3.18) reads

(4.6) \begin{align} \boldsymbol{\nabla } {u_\zeta ^{Rx}}\times \boldsymbol{\nabla } {v^\zeta }=0. \end{align}

In the neighbourhood of points where $({\partial u^{Rx}_\zeta }/{\partial s})\neq 0$ , one can invert the function $u_\zeta =u_\zeta ^{Rx}(s,\theta )$ to $s=s(u^{Rx}_\zeta ,\theta )$ and write $v^\zeta = v^\zeta (u^{Rx}_\zeta ,\theta )$ , which upon using (4.6), implies $v^\zeta$ is a function of only $u_\zeta ^{Rx}$ ; i.e.

(4.7) \begin{align} v^\zeta =\mathcal{F}\big(\nu u^{Rx}_\zeta \big)=\mathcal{F}\big(\big(M^{Rx}_\parallel \big)^2B_\zeta \big) \end{align}

for an arbitrary function $\mathcal{F}$ . Note that (4.7) alone is sufficient for (4.6) to hold. Also, from (2.8b ), (4.7) is equivalent to $u^\zeta =\mathcal{F}(u^{Rx}_\zeta )$ . Although $\mathcal{F}$ is an arbitrary function, one can argue that in its simplest form, it should be a linear function. Assuming the function $\mathcal{F}$ is analytical in a neighbourhood of $(M^{Rx}_\parallel )^2B_\zeta =0$ , it can be expanded as

(4.8) \begin{align} \mathcal{F}\big(\big(M^{Rx}_\parallel \big)^2B_\zeta \big)=\sum _{n=0}^\infty c_nB_\zeta ^n\big(M^{Rx}_\parallel \big)^{2n}. \end{align}

If $M_\parallel ^{Rx}\ll 1$ , in (4.8) one can ignore the $(M_\parallel ^{Rx})^m$ term for $m\geqslant 4$ ; i.e.

(4.9) \begin{align} v^\zeta =\alpha u^{Rx}_\zeta + \omega , \end{align}

where $c_0 \nu =\omega$ , and $c_1 \nu =\alpha$ . Equation (4.9) also means

(4.10) \begin{align} u^\zeta =(\alpha +1) u^{Rx}_\zeta +\omega . \end{align}

Physically, $\omega$ is like the angular frequency of plasma rotating around the axis of symmetry, while $(\alpha +1) u^{Rx}_\zeta$ is a small perturbation which is taken as a function of $u^{Rx}_\zeta$ to satisfy (4.6). Using this equation, (3.8) reads

(4.11) \begin{align} \mu _BA_\zeta +\big[(\alpha +1) \big(M^{Rx}_\parallel \big)^2-1\big] B_\zeta =\mathcal{C}, \end{align}

where $\omega$ is absorbed in the $\mathcal{C}$ . This equation has an obvious singularity in $(M^{Rx}_\parallel )^2=1/(\alpha +1)$ , or equivalently $(\alpha +1)\nu ^2=\rho$ . Because we are assuming $(M^{Rx}_\parallel )^2\ll 1$ , this singularity is avoided as long as $\alpha$ is not too large. Nevertheless, we require a high numerical resolution to resolve the equilibrium points near this singularity accurately.

Using (4.5) and (4.9), (3.17) reads

(4.12) \begin{align} \frac {u^{Rx,s} u^{Rx}_s+u_\theta ^{Rx} u^{Rx,\theta }}{2}+\tau \ln {\frac {\rho }{\rho _\varOmega }}+(\alpha +1)\frac {\big(u^{Rx}_\zeta \big)^2}{2}=0, \end{align}

where any constant on the right-hand side is absorbed in the arbitrary constant $\rho _\varOmega$ .

We solve (4.11) to (4.12) for $A_\theta$ , $A_\zeta$ and $\rho$ , subject to the boundary conditions of (3.3b )–(3.5b ). Other quantities in (2.8) can then be obtained from these three quantities. Although this reduced system of equations does not depend on $\omega$ , the flow will depend on it. Using (3.16a ) and (3.16b ) the cross-field flow ratio is calculated as

(4.13) \begin{align} \frac {u_\perp }{u_\parallel }=\frac {|\nu \alpha B_\zeta B_p+\rho \omega B_p |}{\nu B^2+\nu \alpha B_\zeta ^2+\rho \omega B_\zeta }. \end{align}

Figure 3(a–d) shows a solution in the HKT slab geometry using the parameters listed in table 1. In these parameters the $\nu$ is chosen as a relatively small number so that, with $\rho$ of order unity in our dimensionless units, $M^{Rx}_\parallel \ll 1$ as we discussed in § 3.3. On the other hand, $\tau$ is chosen as a relatively large dimensionless value, corresponding to a hot, low-collisionality plasma in our units, while keeping all quantities of order unity. The magnetic fluxes $\psi _\theta$ and $\psi _\zeta$ are tuned so that the rotational transform profile passes through $\iota =0$ and an island is formed in the domain. This, therefore, allows us to study the interaction of magnetic islands with flow. As we see in figure 3(a), at about $s=0.3$ the rotational transform passes through $\iota =0$ and therefore the resonant islands appear in the Poincaré plot. In this paper, the $\lt \gt$ sign always denotes the $\theta$ -averaged values. The $\theta$ -averaged norm of the magnetic field and flow is maximum at the island position. The cross-field flow vanishes at $s\approx 0.3$ , so that the flow becomes field-aligned there, which corresponds to the position of the O points of the island. As mentioned, this is expected from (3.16b ) regardless of the geometry.

Table 1.

Parameters used in numerical experiments.

Table 1. Long description

A table with two columns and seven rows. The first column is labeled ‘Parameter’ and the second column is labeled ‘Value’. The parameters and their corresponding values are as follows: Row 1: Parameter mu sub B, Value 1; Row 2: Parameter nu, Value 0.1; Row 3: Parameter tau, Value 8; Row 4: Parameter psi sub theta, Value 0.2; Row 5: Parameter psi sub zeta, Value 1; Row 6: Parameter rho sub capital omega, Value 1.

Figure 3.

Example solutions in an HKT slab geometry. (a) Poincaré plot of the magnetic field, (b) magnetic field norm, (c) flow norm and (d) cross-field flow ratio. Here $\delta =0.02$ and $\alpha =10$ used along with the parameters of table 1.

Figure 3. Long description

Panel A: A Poincaré plot of the magnetic field. The plot shows nested magnetic flux surfaces with varying geometries. The x-axis is labeled ‘s’ ranging from 0.0 to 1.0, and the y-axis is labeled ‘θ’ ranging from 0 to 6. An inset graph shows a linear relationship between two variables. Panel B: A heatmap of the magnetic field norm. The x-axis is labeled ‘s’ ranging from 0.0 to 1.0, and the y-axis is labeled ‘θ’ ranging from 0 to 6. The color scale on the right indicates the magnetic field norm values, with a color gradient from purple to yellow. An inset graph shows a curve with a peak. Panel C: A heatmap of the flow norm. The x-axis is labeled ‘s’ ranging from 0.0 to 1.0, and the y-axis is labeled ‘θ’ ranging from 0 to 6. The color scale on the right indicates the flow norm values, with a color gradient from purple to yellow. An inset graph shows a curve with a peak. Panel D: A heatmap of the cross-field flow ratio. The x-axis is labeled ‘s’ ranging from 0.0 to 1.0, and the y-axis is labeled ‘θ’ ranging from 0 to 6. The color scale on the right indicates the cross-field flow ratio values, with a color gradient from purple to yellow. An inset graph shows a curve with a trough.

We emphasise that all quantities in this work are treated as dimensionless, consistent with the dimensionless coordinates introduced in § 3 and the choice $\mu _0=1$ . Equivalently, all variables and parameters are expressed in a single self-consistent system of units in which a reference length, magnetic field and mass density are set to unity. The Lagrange multipliers $\mu _B$ , $\nu$ and $\tau$ are therefore likewise dimensionless in these units. The values in table 1 are not intended to model a specific device; rather, they are chosen so that the equilibrium lies in the regime of interest for this study, namely a small parallel Alfvén Mach number, $(M^{Rx}_\parallel )^2=\nu ^2/\rho \ll 1$ (with $\rho$ of order unity in these units), together with magnetic fluxes that place a resonant surface within the domain so that an $\iota =0$ island forms. The purpose of the numerical examples is to demonstrate that the RxMHD steady-state problem is well-defined and to expose the coupling between cross-field flow and geometry, rather than to reproduce the parameters of a particular experiment. We note that, because all quantities are of order unity in these units, the corresponding local plasma beta $\beta = 2p/B^2$ is also of order unity in the numerical examples. This is a consequence of the order-unity parameter choice adopted for numerical convenience and is not a prediction of the model: a low- $\beta$ regime more representative of tokamaks and stellarators can be accessed by choosing $p \ll B^2$ , for example by choosing a smaller $\tau$ or larger magnetic fluxes.

4.2. The cylindrical geometry

We next consider a hollow cylinder geometry shown in figure 2. In the cylindrical geometry, the position vector is $\boldsymbol{r}=R\hat {r}(\phi )+Z\hat {z}$ , where $\hat {r}$ and $\hat {z}$ are the unit vector basis of the cylindrical coordinate. The inner boundary and outer boundaries, are positioned at $\boldsymbol{\mathcal{R}}_-=r_-\hat {r}+\zeta \hat {z}$ and $\boldsymbol{\mathcal{R}}_+=(r_+-\delta \cos {3\theta })\hat {r}+\zeta \hat {z}$ , respectively. Using these relations, the transformation from the $\{R,\phi , Z\}$ coordinates to the curvilinear coordinates $\{s,\theta ,\zeta \}$ is

(4.14) \begin{align} R&=(1-s)r_-+s[r_+-\delta \cos (3\theta )], \\[-10pt]\nonumber\end{align}

(4.15) \begin{align} \phi &=\theta , \\[-10pt]\nonumber\end{align}

(4.16) \begin{align} z&=\zeta , \end{align}

where $r_-$ and $r_+$ are the inner and outer radii of the hollow cylinder. In the cylindrical geometry, the metric components are $g_{ij}=({\partial R}/{\partial q^i})({\partial R}/{\partial q^j})+({\partial \hat {r}}/{\partial q^i}\boldsymbol{\cdot })({\partial \hat {r}}/{\partial q^j})R^2+({\partial z}/{\partial q^i})({\partial z}/{\partial q^j})$ , which gives

(4.17) \begin{align} g_{ij}=\left ( \begin{array}{c@{\quad}c@{\quad}c} g_{11} & g_{12} & 0 \\[6pt] g_{21} & g_{22} & 0 \\[6pt] 0 & 0 & 1 \end{array} \right ), \end{align}

where $g_{11}=(\delta \cos (3\theta )+r_- - r_+)^{2}$ $g_{12}=g_{21}= -3\,\delta \, s\, \sin (\theta )(2\cos (2\theta )+1) (\delta \cos$ $ (3\theta )+r_- - r_+)$ and $g_{22}= (r_- (s-1)-r_+ s)(r_- (s-1)-r_+ s+2\delta s \cos (3\theta ))+\delta ^{2} s^{2}$ $(5-4\cos (6\theta ))$ . Because in the cylindrical geometry we have $g_{33}=1$ , the condition of (3.18) is satisfied in the same way that it was in the HKT slab; i.e. by defining $v_\zeta$ in the form of (4.9). Therefore, the system of equations is the same system that was solved for the HKT slab.

Figure 4(a–d) shows a solution of the model in the cylindrical geometry, using the same parameters as the HKT slab. In comparison with the slab geometry, we see in figure 4(a) that the rotational transform has a convex profile, likely because of the reduction in toroidal current cross-section due to the added curvature in the poloidal direction. As a result, the island is moved towards the interior boundary to $ s \approx 0.17$ . Also, because the added curvature compressed the field lines close to the interior boundary, the norm of the magnetic field is slightly elevated in that region (figure 4 b). This leads to an increase in the flow norm close to the $s=0$ boundary, because the relaxed flow is proportional to the magnetic field (figure 4 c). The increase in the relaxed flow increases the parallel flow, and as a result, the cross-field flow ratio decreases close to the $s=0$ boundary, as seen in figure 4(d). In the cylindrical geometry, we can again see that the cross-field flow vanishes at the O point of the island.

Figure 4.

Example solutions in cylindrical geometry. (a) Poincaré plot of the magnetic field, (b) magnetic field norm, (c) flow norm and (d) cross-field flow ratio. Here $\delta =0.02$ , $\alpha =10$ , $r_-=0.5$ and $r_+=1$ used along with the parameters of table 1.

4.3. Effect of the $\alpha$ on the cross-field flow ratio, equilibrium profiles and magnetic islands

In the slab and cylindrical geometries, the system of equations does not depend on $\omega$ and therefore the value of $\omega$ does not affect the magnetic field or pressure profiles. Although $\omega$ contributes to the toroidal flow, this contribution is evident from (4.10). The effect of the particular choice of $\alpha$ , however, is not immediately obvious, because it explicitly appears in the system of equations. For this reason, in this section we set $\omega = 0$ and vary $\alpha$ to investigate how the equilibrium profiles and the flow characteristics are affected by this parameter. The choice $\omega = 0$ is also physically motivated. From (4.10), if $\omega$ is interpreted as the dominant contribution to the toroidal flow and $\alpha u_\zeta ^{Rx}$ is treated as a small perturbation, then $\omega = 0$ is expected in slab and cylindrical geometries. This is because, in these geometries, the toroidal direction is bounded by two walls which prevent substantial ion flow. Additionally, in this section we restrict our investigation to the cylindrical geometry. Since the system of equations is identical for the slab and cylindrical geometries, the effects observed here apply equally to both geometries.

Figure 5(a) illustrates the cross-field flow ratio as a function of $\alpha$ in the cylindrical geometry. The cross-field flow ratio in this figure is reported as the total average of (4.13) in the whole domain. As expected from (4.13), when $\alpha =0$ the cross-field flow vanishes and the cross-field flow ratio is zero, i.e. the flow becomes completely field-aligned. Therefore, this case is similar to the MRxMHD equilibrium with flow Qu et al. (Reference Qu, Dewar, Ebrahimi, Anderson, Hudson and Hole2020). In figures 5(b) and 5(d) we see that in this case, the flow and pressure profiles are almost flat (though some small changes across the domain still exist). As we increase $\alpha$ from zero, the cross-field flow increases and the cross-field flow ratio grows. As this ratio increases, the change of pressure, flow and magnetic field profiles across the domain increases. As $\alpha \rightarrow \infty$ , $({u_\perp }/{u_\parallel })\rightarrow (|{B_p}/{B_\zeta }|)$ . The rotational transform also changes as $\alpha$ (and therefore, the cross-field flow) increases. However, from figure 5(c), one can see that this change is relatively small. Rather, the rotational transform profile is crucially affected by the poloidal and toroidal magnetic fluxes. In our numerical experiments, we also observed that the width of the islands is not significantly affected by flow parameters such as $\alpha$ . It is likely because the islands are primarily determined by the domain geometry, which, because of $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}=0$ boundary condition, determines the geometry of magnetic field lines. As we will see, this conclusion does not hold in the geometries with a finite curvature along the $\zeta$ axis, where the flow and the geometry become strongly coupled.

Figure 5.

Effect of $\alpha$ on the (a) cross-field flow ratio, (b) flow norm, (c) magnetic field and (d) pressure, in cylindrical geometry.

Figure 5. Long description

Panel A: A line graph shows the cross-field flow ratio as a function of alpha. The x-axis is labeled alpha, and the y-axis is labeled (u_r / u_total). The graph shows an increasing trend as alpha increases from 0 to 10. Panel B: A line graph displays the flow norm as a function of s for different values of alpha (0, 5, 10). The x-axis is labeled s, and the y-axis is labeled (u). The graph shows three lines, each representing a different alpha value, with the flow norm decreasing as s increases from 0 to 1. Panel C: A line graph illustrates the magnetic field as a function of s for different values of alpha (0, 5, 10). The x-axis is labeled s, and the y-axis is labeled (B). The graph shows three lines, each representing a different alpha value, with the magnetic field decreasing as s increases from 0 to 1. An inset graph shows the derivative of the magnetic field with respect to s. Panel D: A line graph presents the pressure as a function of s for different values of alpha (0, 5, 10). The x-axis is labeled s, and the y-axis is labeled (p). The graph shows three lines, each representing a different alpha value, with the pressure increasing as s increases from 0 to 1.

It is interesting to note that, in contrast to ideal axisymmetric equilibrium codes such as FLOW (Guazzotto et al. Reference Guazzotto, Betti, Manickam and Kaye2004), the toroidal flow in our model can vanish, and flow can become exclusively poloidal. In the slab and cylindrical geometries, this can be achieved by setting $\alpha =-1$ and $\omega =0$ in (4.10). In such cases, (3.17) reduces to $\boldsymbol{\nabla } {h_p}=0$ or $ (({u^s u_s+u_\theta u^\theta })/{2})+\tau \ln ({{\rho }/{\rho _\varOmega }})=0$ , which should be solved alongside (3.8) and (3.11). This possibility is likely a consequence of not enforcing the ideal constraint $\boldsymbol{\nabla }\times {(\boldsymbol{u}\times \boldsymbol{B})}=0$ , as done in Guazzotto et al. (Reference Guazzotto, Betti, Manickam and Kaye2004), which inherently couples the poloidal direction with the toroidal direction.

5.1. Effect of the angular frequency $\omega$ on the cross-field flow ratio and the geometry of magnetic island

In (5.6) and (5.9), one can see that, unlike the slab and cylindrical geometry (where $g_{33}=1$ ), in geometries where $g_{33}$ is not constant – such as in helical and toroidal geometries – the $\omega$ is not absorbed in the arbitrary constants and therefore appears in the equations. More importantly, in (3.21), (5.6) and (5.9), the $\omega$ is coupled with the geometry; i.e. it is coupled with the $g_{33}$ and $g_{32}$ components of the metric. Because of this, one expects that adjusting the cross-field flow through $\omega$ changes geometrical effects such as islands. This point will be clarified in this section. Using (3.16a ) and (3.16b ) one can find the norm of the parallel and cross-field flows as

(5.10a) \begin{align} u_\parallel &=\frac {\nu B^2+\rho \omega B_\zeta }{\rho B} ,\\[-8pt]\nonumber\end{align}

(5.10b) \begin{align} u_\perp &=\frac {|\omega |\sqrt {R^2B^2-B_\zeta ^2}}{B} \end{align}

and the cross-field flow ratio as

(5.11) \begin{align} \frac {u_\perp }{u_\parallel }=\frac {|\omega | \rho \sqrt {R^2B^2-B_\zeta ^2}}{\nu B^2+\rho \omega B_\zeta } .\end{align}

Figure 8.

(a–g) The evolution of primary and secondary magnetic islands versus $\omega$ . (h) The primary island’s width, measured across the O point in the $s$ direction, cross-field flow ratio averaged in the island region from $s=0.1$ to $s=0.6$ , and the absolute value of the first and second Fourier harmonics of $B^s$ .

To demonstrate the interaction of the flow and the magnetic island, we change the angular frequency $\omega$ and solve the model in the toroidal geometry for the parameters also used in figure 7(a–d). Figure 8(a–g) shows the evolution of the Poincaré plots of the magnetic field lines as $\omega$ changes from –3 to 3. Corresponding to each value of $\omega$ , in figure 8(h) we have shown the width of the magnetic island, the cross-field flow ratio and the amplitude of the first and second Fourier harmonics of $B^s$ . The width of the primary island is measured across the O point from separatrix to separatrix in the $s$ direction, using the Poincaré section. We note that for all values of $\omega$ in these plots, the toroidal magnetic field $B^\zeta$ remains positive, and therefore, the negative sign of $\omega$ points to the counter-toroidal-field direction. In figure 8(a), we see a large island around $s\approx 0.35$ for $\omega =-3$ . This island is initially located at $\theta \approx \pi$ , located at the inboard of the torus in figure 6. In figure 8(b), we see that as $\omega$ increases to $-1.25$ , the width of this island shrinks. In figure 8(c), we see that as we further increase $\omega$ , the primary island slightly shrinks while a secondary island appears inside it. As we keep increasing $\omega$ , the width of the primary island across the O point decreases, and for $\omega =-0.65$ , the primary island in figure 8(d) disappears and only the secondary islands remain in the plot. Upon increasing $\omega$ to $-0.5$ , the secondary islands move in opposite directions towards each other, and they merge at around $\theta =0$ (equivalent to $\theta =2\pi$ ) and recreate a primary island around that point (figure 8 e). From there, by increasing the $\omega$ , the secondary island disappears, and the width of the primary island grows. This behaviour of the primary and secondary islands can be explained based on the first and second harmonics of $B^s$ in figure 8(h). In this figure, we see that as $\omega$ increases from $-3$ , the width of the primary island and the value of $b^s_{m=1}$ decrease. The value of the second Fourier harmonics $b^s_{m=2}$ , however, remains almost constant. Between $\omega \approx -1$ and $\omega \approx -0.4$ , the second harmonic is dominant. The Poincaré plots of figures 8(d) and 8(e) belong to this range of $\omega$ , in which the primary island either disappears or is small. For $\omega \gtrsim -0.4$ , the first Fourier harmonic and the island width grow. Figure 8(h) also shows the correlation between the primary island width and the cross-field flow ratio ( ${u_\perp }/{u_\parallel }$ ). The cross-field flow ratio decreases with a proportional rate to the primary island’s width for $-0.65\geqslant \omega \geqslant -3$ and reaches its minimum at $\omega =-0.65$ . For $\omega \gt -0.65$ , the cross-field flow ratio grows quickly as the island width also grows. These results show that the flow, and especially the cross-field flow ratio, can significantly affect the geometry of magnetic islands and even break a primary island into a secondary island.

The qualitative dependence of the island size can be captured from (3.21), after some approximations. In the toroidal geometry (3.21) reads

(5.12) \begin{align} & -\frac {\mu _B}{R^2}\left (\frac {\mathcal{C}-\mu _BA_\zeta -\nu \omega R^2}{\nu ^2/\rho -1}\right )+\frac {1}{\sqrt {g}}\frac {\partial }{\partial s}\left [\frac {(\nu ^2/\rho -1)}{\sqrt {g}}\frac {\partial\kern-1.1pt A_\zeta }{\partial s}\right ]\nonumber \\[5pt]& \qquad +\frac {1}{\sqrt {g}}\frac {\partial }{\partial \theta }\left [\frac {(\nu ^2/\rho -1)}{\sqrt {g}}\frac {\partial\kern-1.1pt A_\zeta }{\partial \theta }\right ]=0, \end{align}

with $\sqrt {g}=w(sw+r_-)(R_0+(sw+r_-)\cos \theta )$ and $w\equiv r_+-r_-$ . We note that the poloidal flux function, $\varPsi (s,\theta )\equiv \int _{0}^{s}\int _0^{2\pi }B^\theta (s,\theta ,\zeta )\sqrt {g}\mathrm{d} s \mathrm{d}\zeta =-2\pi A_\zeta$ , and therefore, (5.12) can be seen as a modified Grad–Shafranov equation written in general curvilinear coordinates. As we mentioned, we are interested in the $(M^{Rx}_\parallel )^2=({\nu ^2}/{\rho })\ll 1$ regime, and therefore (5.12) reads

(5.13) \begin{align} \frac {\mu _B}{R^2}(\mathcal{C}-\mu _BA_\zeta -\nu \omega R^2)-\frac {1}{\sqrt {g}}\frac {\partial }{\partial s}\left (\frac {1}{\sqrt {g}}\frac {\partial\kern-1.1pt A_\zeta }{\partial s}\right ) -\frac {1}{\sqrt {g}}\frac {\partial }{\partial \theta }\left (\frac {1}{\sqrt {g}}\frac {\partial\kern-1.1pt A_\zeta }{\partial \theta }\right )=0. \end{align}

This equation is now linearised and decoupled from (5.7). To further approximate this equation, we use our numerical observation that close to the island $A_\zeta$ will always show a parabolic shape with an extremum at the O point. The term $(\partial/\partial s)((1/\sqrt {g})(\partial\kern-1.1pt A_\zeta /\partial s))$ is therefore approximated as $\mathcal{C}'/\sqrt {g}$ , with $\mathcal{C}'$ being a constant found from the boundary condition (3.3c ) ( $\mathcal{C}'= -({10\psi _\theta }/{3\pi })$ ). Equation (5.13) therefore reads

(5.14) \begin{gather} \frac {\mu _B}{R^2}(\mathcal{C}-\mu _BA_\zeta -\nu \omega R^2)-\frac {\mathcal{C}'}{(\sqrt {g})^2} -\frac {1}{\sqrt {g}}\frac {\partial }{\partial \theta }\left (\frac {1}{\sqrt {g}}\frac {\partial\kern-1.1pt A_\zeta }{\partial \theta }\right )=0. \end{gather}

Furthermore, we assume the aspect ratio is large ( $({(r_-+sw)}/{R_0})\ll 1$ ) and use expansions $1/R\approx (1/R_0)(1-\cos \theta )(1-(r_-+sw)\cos \theta /R_0)$ and $1/R^2\approx (1/R_0^2)(1-2(r_-+sw)\cos \theta /R_0)$ . Around the primary island, we expect the first Fourier harmonics to be dominant and therefore, we truncate the Fourier series at the first harmonics and write $A_\zeta \approx a_{\zeta 0}(s)+a_{\zeta 1}(s)e^{i\theta }$ . Substituting this relation in (5.14) and demanding that the zeroth and first Fourier harmonics of the left-hand side of this equation vanish, we get

(5.15a) \begin{align} \frac {\mathcal{C}\mu _B}{R_0^2} -\frac {\mathcal{C}'\left (2 R_0^2+\varrho ^2\right )}{2 R_0^4 w^2 \varrho ^2} -\frac {\mu _B^2}{R_0^2}\,a_{\zeta 0} +\left (\frac {3 r_{-}-s w}{2 R_0^3 w^2 \varrho ^2} +\frac {\mu _B^2 \varrho }{R_0^3}\right )a_{\zeta 1} &=\mu _B\nu \omega ,\\[-10pt]\nonumber \end{align}

(5.15b) \begin{align} \frac {\dfrac {\mathcal{C}'}{w^2}-\mathcal{C}\mu _B\varrho ^2}{R_0^3\varrho } +\frac {\mu _B^2 \varrho }{R_0^3}\,a_{\zeta 0} +\left (\frac {2 R_0^2\!\left (\varrho ^{-2}-\mu _B^2 w^2\right )+1}{2 R_0^4 w^2}\right )a_{\zeta 1} &=0, \end{align}

where $\varrho \equiv r_- + s w$ . In these equations, $\omega$ only appears on the right-hand side of (5.15a ). Because of this, the effect of $\omega$ on the first Fourier harmonics (and therefore, the primary island’s width) is due to the coupling between $a_{\zeta 0}$ and $a_{\zeta 1}$ in (5.15a ). This coupling vanishes in the limit $({(r_-+sw)}/{R_0})\rightarrow 0$ , i.e. when there is no curvature in $\zeta$ direction, and the torus reduces to a cylinder. We conclude that the effect of the flow on the $a_{\zeta 1}$ and therefore the primary island’s width is due to the toroidal geometry, namely due to the curvature of the torus in the $\zeta$ direction. Solving (5.15a ) and (5.15b ), we can find the $a_{\zeta 0}(s)$ and $a_{\zeta 1}(s)$ coefficients and therefore $A_\zeta$ . Finally, noting that $B^s={(\partial\kern-1.1pt A_\zeta }/{\partial \theta })/\sqrt {g}$ , we find the amplitude of the first Fourier harmonic of $B^s$ as

(5.16) \begin{align} b_{m=1}^s=\left | \frac {(r_-+s w)^3 \left (\mathcal{C}'+2 \mu _B \nu R_0^4 w^2 \omega \right )}{4 R_0^5 w \left (\mu _B ^2 w^2 (r_-+s w)^2-1\right )-4 \mu _B ^2 R_0^3 w^3 (r_-+s w)^4}\right |. \end{align}

The dependence of the $b_{m=1}^s$ on $\omega$ is in the form of a shifted absolute value function, which is qualitatively consistent with the full numerical results in figure 8(h). Moreover, in (5.16) the shift of the absolute value is proportional to the $\mathcal{C}'=-({10\psi _\theta }/{3\pi })$ . We validated this prediction with the full numerical solutions. In numerical solutions, we observed that for positive values of $\psi _\theta$ the absolute value form of $b^s_{m=1}$ is shifted towards negative $\omega$ , while for negative $\psi _\theta$ it is shifted towards the positive $\omega$ . For $\psi _\theta =0$ , however, $b^s_{m=1}$ at $\omega =0$ does not vanish, and a relatively small island will still exist in contradiction with (5.16). This discrepancy is likely rooted in our large-aspect-ratio assumption and the first-order expansions that we used to capture the qualitative behaviour of $b_{m=1}^s$ . Another effect that needs analysis at least to the second Fourier harmonic is the secondary islands and their dominance around the $-1\lesssim \omega \lesssim -0.4$ regime.