# Steady plasma flows in a periodic non-symmetric domain

Harold Weitzner 1  Wrick Sengupta 2,3<EMAIL_ADDRESS>1Courant Institute of
Mathematical Sciences, New York University, New York, New York 10012, USA
2Department of Astrophysical Sciences, Princeton University, Princeton, NJ,
08543 3Princeton Plasma Physics Laboratory, Princeton, NJ, 08540

###### Abstract

Steady plasma flows have been studied almost exclusively in systems with
continuous symmetry or in open domains. In the absence of continuous symmetry,
the lack of a conserved quantity makes the study of flows intrinsically
challenging. In a toroidal domain, the requirement of double-periodicity for
physical quantities adds to the complications. In particular, the
magnetohydrodynamics (MHD) model of plasma steady-state with the flow in a
non-symmetric toroidal domain allows the development of singularities when the
rotational transform of the magnetic field is rational, much like the
equilibrium MHD model. In this work, we show that steady flows can still be
maintained provided the rotational transform is close to rational and the
magnetic shear is weak. We extend the techniques developed in carrying out
perturbation methods to all orders for static MHD equilibrium by Weitzner
(Physics of Plasmas 21, 022515 (2014)) to MHD equilibrium with flows. We
construct perturbative MHD equilibrium in a doubly-periodic domain with nearly
parallel flows by systematically eliminating magnetic resonances order by
order. We then utilize an additional symmetry of the flow problem, first
discussed by E. Hameiri in (J. Math. Phys. 22, 2080 (1981) Sec. III), to
obtain a generalized Grad-Shafranov equation for a class of non-symmetric
three-dimensional MHD equilibrium with flows both parallel and perpendicular
to the magnetic field. For this class of flows, we can obtain non-symmetric
generalizations of integrals of motion, such as Bernoulli’s function and
angular momentum. Finally, we obtain the generalized Hamada conditions, which
are necessary to suppress singular currents in such a system when the magnetic
field lines are closed. We do not attempt to address the question of
neoclassical damping of flows.

##  1 Introduction

Flows in magnetic confinement devices are known to have a strong influence in
suppressing turbulence and affecting neoclassical (Hinton & Wong, 1985;
Helander, 1998) and turbulent transport (Burrell, 1997), as well as stability
(Chu et al., 1995). In particular, sheared flows like zonal flows can reduce
turbulent transport by shearing the turbulent eddies (Lin et al., 1998;
Fujisawa, 2008; Terry, 2000). Study of the interaction of plasma flows and
turbulence is critical (Stroth et al., 2011; Groebner et al., 1990) in
understanding the edge region of tokamaks and the L-H transition. It has been
experimentally verified that the radial electric field plays a key role in the
transition (Burrell, 1999). Flows can also have applications in healing
magnetic islands in stellarators (Hegna, 2012). Flows are also an essential
object of study in space and astrophysical plasmas (Goedbloed et al., 2004;
Beskin, 2009; Istomin & Pariev, 1996; Davis & Tchekhovskoy, 2020).

If the flow speed is assumed on the order of the sound speed, ideal MHD is a
good description of plasmas since gyroradius corrections do not enter (Morozov
& Solov’ev, 1980; Freidberg, 2014; Goedbloed et al., 2004). The theory of
steady axisymmetric plasma flows is well developed (Hameiri, 1983; Hassam,
1996; Guazzotto & Betti, 2005; McClements & Hole, 2010; Tasso &
Throumoulopoulos, 1998; Abel et al., 2013; Beskin, 1997). Toroidal angular
momentum conservation, which follows from Noether’s theorem, facilitates the
reduction of axisymmetric flows to the study of a modified Grad-Shafranov (GS)
equation which includes the effect of steady plasma flow (Hameiri, 1983;
Hassam, 1996; Tasso & Throumoulopoulos, 1998; Goedbloed & Lifschitz, 1997;
Beskin, 1997). Similarly, for plasma systems with helical symmetry, a helical
Grad-Shafranov with large flows can be obtained (Andreussi et al., 2012;
Villata & Tsinganos, 1993). Global solutions of such flow modified Grad-
Shafranov equations provide useful models of plasma jets in space
(Bogoyavlenskij, 2000a; Villata & Ferrari, 1994; Bogoyavlenskij, 2000b;
Beskin, 1997).

While the standard GS equation for static equilibrium is always elliptic, the
flow-modified GS equation can be elliptic, parabolic, or hyperbolic, depending
on the flow’s Mach number. Therefore, even with symmetry, complicated flows
such as transonic flows (Morawetz, 1985) that modify the characteristics in a
given domain can appear (Lifschitz & Goedbloed, 1997).

There is considerable interest in understanding the effects of symmetry
breaking three-dimensional perturbations and stability of three-dimensional
MHD flows in astrophysical plasmas (Igumenshchev et al., 2003; McKinney &
Blandford, 2009; Istomin & Pariev, 1996; Igumenshchev, 2008). In particular,
2D and 3D simulations of magnetically arrested disks show very different
behaviors (White et al., 2019; Igumenshchev, 2008). However, high-resolution
3D simulations are required to carefully capture the effects of symmetry-
breaking (White et al., 2019).

In the absence of symmetry, reduction to a Grad-Shafranov like equation fails
in general. Non-symmetric toroidal MHD equilibrium with flows share the same
mathematical difficulty as the static non-symmetric MHD equilibrium, namely,
the appearance of singularities on rational surfaces.

Although large steady shear-flows are desirable in various confinement
geometries, steady flows show a strong alignment with a symmetry direction
when present (Helander, 2007; Helander & Simakov, 2008; Helander, 2014; Spong,
2005). In axisymmetric tokamaks, the persistent flow is in the toroidal
symmetry direction (Throumoulopoulos et al., 2006; Tasso & Throumoulopoulos,
1998). While the toroidal angular momentum in a tokamak is typically conserved
over several inverse collision frequency (Helander et al., 2012), poloidal
flows are damped due to neoclassical effects. Axisymmetry makes tokamaks
intrinsically ambipolar, i.e., the radial fluxes of ions and electrons are
equal and independent of the radial electric field in the leading order
gyroradius expansion, leading to flow in the symmetry direction being
unconstrained.

In the absence of intrinsic amibipolarity, the rotation is damped to the value
required for ambipolar radial transport and typically leads to flows of the
order of diamagnetic flows (Helander & Simakov, 2008; Helander, 2014). Plasma
rotation in tokamaks and stellarators, therefore, differ significantly because
stellarators do not, in general, have a symmetry direction and are not
intrinsically ambipolar (Helander et al., 2012). The only class of
stellarators that have been shown to support large flows are quasisymmetric
(Tessarotto et al., 1996; Canik et al., 2007; Spong, 2005; Helander, 2007;
Helander & Simakov, 2008). In a quasisymmetric stellarator, continuous
symmetry comes from the symmetry of the magnitude of the magnetic field. As a
result of the symmetry, such stellarators are also approximately intrinsically
ambipolar (Helander & Simakov, 2008). However, it is difficult, if not
impossible, to obtain exact global quasisymmetry in a volume (Garren & Boozer,
1991; Landreman & Sengupta, 2019). The breakdown of quasisymmetry leads to
flow damping (Spong, 2005; Simakov & Helander, 2011; Helander & Simakov,
2008). In the literature, it has been further argued in the literature (Sugama
et al., 2011; Simakov & Helander, 2011; Tessarotto et al., 1996) that
quasisymmetry is not enough to accommodate the large centrifugal forces from
sonic flows, and one needs axial/helical symmetries.

As a result of the difficulties imposed by the three-dimensional geometry, a
lack of symmetry that leads to magnetic resonances on rational surfaces and
stringent restrictions on plasma flows from neoclassical physics, a complete
mathematical description of the characteristics of a generic 3D ideal MHD
equilibrium, let alone stability of such systems is not available in the
standard literature. The analytic description of nonsymmetric 3D flows has
employed various approximations such as incompressibility (Kamchatnov, 1982)
or large aspect ratio (Kovrizhnykh & Shchepetov, 1989, 1980) have been
invoked.

Several variational formalisms for 3D plasma flow problem exist (Hameiri,
1998; Greene & Karlson, 1969; Ilgisonis & Pastukhov, 2000; Vladimirov &
Moffatt, 1995). A nontrivial result from the variational formalism is the
proof that under certain conditions ideal MHD allows more than one kind of
cross-helicity (Hameiri, 1981; Ilgisonis & Pastukhov, 2000; Vladimirov &
Moffatt, 1995). We shall call the additional symmetry vector associated with
the conserved quantity Hameiri’s vector because it was first discovered in
(Hameiri, 1981) and later derived independently in (Ilgisonis & Pastukhov,
2000) and (Vladimirov & Moffatt, 1995) (see discussion in (Hameiri, 1998) and
(Ilgisonis & Pastukhov, 2000)). A major disadvantage of the flow-modified
variational formalism is that, unlike the static limit, the energy functional
for the flow problem does not possess a minimum unless the equilibrium flow is
small enough or is parallel to the magnetic field (Hameiri, 1998).

In this work, instead of addressing the flow damping problem in a 3D system,
we take a step back and analyze the ideal MHD equilibrium with steady flows in
a non-symmetric geometry. Our goal here is to construct a class of
perturbative solutions to non-symmetric ideal MHD with steady flows and obtain
consistency conditions for steady flows when the magnetic fields close on
themselves. The flows that we consider are larger than diamagnetic flows but
are not sonic. We show that the same procedure that allowed us to obtain non-
symmetric ideal MHD static equilibrium in a low shear stellarator can be
generalized to MHD with such flows. In particular, we highlight the crucial
role of closed field lines Grad (1971); Weitzner (2014, 2016) in avoiding
magnetic resonances and going to arbitrarily high order in the perturbative
analysis. We also show that for a special class of flows that possess the
additional symmetry, the Hameiri’s vector, one can obtain a generalized Grad
Shafranov equation (GGS) similar to the quasisymmetric generalized Grad
Shafranov equation first derived in (Burby et al., 2020a) (see also related
works by Burby et al. (2020b); Constantin et al. (2020)). The GGS reduces to
the standard flow-modified Grad-Shafranov equation in symmetric geometries. We
explicitly obtain the necessary conditions that need to be satisfied on
rational surfaces. If these conditions are not satisfied current singularities
can develop on rational surfaces (Solov’ev & Shafranov, 1970; Grad, 1967;
Boozer, 1981).

The outline of the paper is as follows. First, in section 2, we study the
ideal MHD equilibrium with flows and identify the characteristic surfaces.
Then, in section 3 we study the perturbation of a class of exact solutions and
obtain the dispersion relation. Then, we outline in section 4 how the
perturbation theory can be carried out to all orders for nearly parallel
flows. Next, we utilize Hameiri’s vector to derive the GGS and the necessary
conditions on rational surfaces in section 5. Finally, we discuss the
implications of our results in section 6.

##  2 Preliminaries

We employ the standard model of ideal magnetohtdrodynamics to represent the
plasma in steady flow state but we add the simplifying assumption that the
entropy is a constant in the entire plasma. The system is

$\displaystyle\bm{\nabla}\cdot(\rho\bm{u})$ $\displaystyle=0$ (2.1a)
$\displaystyle\rho\bm{u}\cdot\bm{\nabla}\bm{u}+\bm{\nabla}p$
$\displaystyle=\bm{J}\times\bm{B},\quad\bm{J}=\bm{\nabla}\times\bm{B}$ (2.1b)
$\displaystyle\bm{\nabla}\cdot\bm{B}$ $\displaystyle=0$ (2.1c)
$\displaystyle-\bm{E}=\bm{\nabla}\Phi$ $\displaystyle=\bm{u}\times\bm{B}.$
(2.1d)

It is often convenient to replace (2.1d) by the equivalent expression

$\displaystyle\bm{\nabla}\times\left(\bm{u}\times\bm{B}\right)=0$ (2.1d*)

While (2.1d) and (2.1d*) are equivalent, the form (2.1d) makes explicit the
need to select specific EMFs, equivalently periods associated with the
potential $\Phi$. As noted we add the assumption

$\displaystyle p=p(\rho)$ (2.2)

to close the system. For concreteness, we shall assume an adiabatic equation
of state of the form $p\sim\rho^{\gamma}$, where $\gamma$ is the ratio of
specific heats. If we introduce

$\displaystyle\Pi\equiv p+\frac{1}{2}\bm{B}\cdot\bm{B}$ (2.3)

then an equivalent form for (2.1b) is

$\displaystyle\rho\bm{u}\cdot\bm{\nabla}\bm{u}+\bm{\nabla}\Pi=\bm{B}\cdot\bm{\nabla}\bm{B}$
(2.1b*)

We start with a formal mathematical exploration of the system (2.1a)-(2.2). It
is clear there are eight scalar first-order differential equations for the
eight scalar unknowns $\rho,\Phi,\bm{u},\bm{B}$. If we were to use (2.1d*)
instead of (2.1d), there would be one fewer unknown but also a hidden
solvability condition for (2.1d*) that adds an unknown and leads back to
(2.1d). We now wish to determine the “type” of the system, elliptic,
hyperbolic, parabolic, mixed, or whatever may be. A closely related issue is
the description of the standing waves in the linearized system. We linearize
about a constant state for $\rho,\bm{u}$ and $\bm{B}$ and look for waves
propagating as $\exp{(i\bm{k}\cdot\bm{x})}$. Such solutions are special cases
of simple waves in fluid dynamics (Courant & Friedrichs, 1999). With
$c^{2}\equiv dp/d\rho$, we obtain the following linear system

$\displaystyle(\bm{k}\cdot\bm{u})\delta\rho+\rho(\bm{k}\cdot\delta\bm{u})$
$\displaystyle=0$ (2.4a)
$\displaystyle\rho(\bm{k}\cdot\bm{u})\delta\bm{u}+c^{2}\bm{k}\delta\rho$
$\displaystyle=\delta\bm{B}(\bm{k}\cdot\bm{B})-\bm{k}(\delta\bm{B}\cdot\bm{B})$
(2.4b)
$\displaystyle\bm{k}\times(\delta\bm{u}\times\bm{B})+\bm{k}\times(\bm{u}\times\delta\bm{B})$
$\displaystyle=0$ (2.4c) $\displaystyle\bm{k}\cdot\delta\bm{B}$
$\displaystyle=0.$ (2.4d)

It follows easily that given the first order perturbations $\partial\rho$,
$\delta\bm{u}$ and $\delta\bm{B}$ , they satisfy the system (2.7a)-(2.7d)

From (2.4c) we find

$\displaystyle\delta\bm{u}(\bm{k}\cdot\bm{B})-\bm{B}(\bm{k}\cdot\delta\bm{u})-\delta\bm{B}(\bm{k}\cdot\bm{u})=0$
(2.4c*)

while (2.4b) together with (2.4a) yields

$\displaystyle\rho(\bm{u}\cdot\bm{k})^{2}\delta\bm{u}-\rho\bm{k}c^{2}(\bm{k}\cdot\delta\bm{u})=\left(\bm{k}\times\left(\delta\bm{u}(\bm{k}\cdot\bm{B})-\bm{B}(\bm{k}\cdot\delta\bm{u})\right)\right)\times\bm{B}$
(2.5)

Taking components of (2.5) in $\bm{k}$ and $\bm{B}$ directions, we obtain

$\displaystyle\rho\left((\bm{k}\cdot\bm{u})^{2}-c^{2}k^{2}\right)\bm{k}\cdot\delta\bm{u}=k^{2}B^{2}(\bm{k}\cdot\delta\bm{u})-k^{2}(\bm{k}\cdot\bm{B})(\bm{B}\cdot\delta\bm{u})$
(2.6a)
$\displaystyle\rho(\bm{k}\cdot\bm{u})^{2}(\bm{B}\cdot\delta\bm{u})-\rho(\bm{k}\cdot\bm{B})c^{2}(\bm{k}\cdot\delta\bm{u})=0.$
(2.6b)

A simple calculation shows that the condition for a solution is the
“dispersion relation” :

$\displaystyle\left(\rho(\bm{k}\cdot\bm{u})^{2}\right)^{2}-\rho(\bm{k}\cdot\bm{u})^{2}\left(\rho
c^{2}+B^{2}\right)k^{2}+k^{2}(\bm{k}\cdot\bm{B})^{2}\rho c^{2}=0$ (2.7)

We observe that this condition is satisfied identically for all vectors
$\bm{k}$ perpendicular to $\bm{u}$ and $\bm{B}$. The expression also shows
connections to the fast-slow wave combination found in ideal
magnetohydrodynamic wave propagation studies. Note, however, that no Alfvén
wave appears. Moreover, the magnitude of $\bm{k}$ does not appear in the
system; only its direction occurs. We study the dispersion relation (2.7) and
its solutions in much more detail in Appendix A. In particular, we show that
real physical solutions of (2.7) exist.

An alternate interpretation of the relation is also significant. We might pose
the question: What are the characteristic surfaces for the steady flow system?
Characteristic surfaces are defined by the following property. Suppose one
specifies $\rho,\bm{u}$ and $\bm{B}$ on a given surface. Can one then
determine the normal derivative of these quantities on that surface? If so,
the surface is not characteristic. On a characteristic surface, the flow
variables are not arbitrary; they must be related. These surfaces define
important properties of any solution. Suppose $\chi=0$ is such a surface and
define $\bm{k}$ as the gradient of $\chi$. Hence $\bm{k}$ is normal to the
surface. If one further defines $\delta\rho,\delta\bm{u}$ and $\delta\bm{B}$
as the gradients of these variables dotted into the normal to the surface,
then these quantities satisfy the system (2.7), except that there would be
inhomogeneous terms added to the equations, terms given solely in terms of the
variables on the given surface. With this interpretation of the system, the
relation is the condition that one cannot give the unknown functions
arbitrarily on the surface in question. Such peculiar surfaces, the
characteristic surface, are thus defined by (2.7).

A physical problem without such real surfaces is typically elliptic. When they
exist, there are usually hyperbolic properties of the system as well. Their
appearance indicates the complexity of steady magnetohydrodynamic flow. Flows
with both elliptic and hyperbolic characteristics can occur in ideal
magnetohydrodynamic equilibria. The mixed nature (neither fully elliptic nor
hyperbolic) of MHD leads to difficult unresolved mathematical issues.

##  3 A simple linearized problem

In order to gain some insight into the nature of steady flow states we examine
a simple linearized problem. We start from an exact solution of the system
(2.1a)-(2.1d)

$\displaystyle\bm{B}$ $\displaystyle=(0,f(x),g(x))$ (3.1a)
$\displaystyle\bm{u}$ $\displaystyle=(0,v(x),w(x))$ (3.1b) $\displaystyle p$
$\displaystyle=p(x)$ (3.1c) $\displaystyle\Pi$
$\displaystyle=p(x)+\frac{1}{2}(f^{2}+g^{2})$ (3.1d) $\displaystyle-\bm{E}$
$\displaystyle=\Phi_{,x}=v(x)g(x)-w(x)f(x).$ (3.1e)

We linearize about the state and introduce the perturbed variables with the
structure

$\displaystyle\bm{B}^{(1)}$
$\displaystyle=(B^{(1)}_{x}(x)\sin{(my+nz)},B^{(1)}_{y}(x)\cos{(my+nz)},B^{(1)}_{z}(x)\cos{(my+nz)})$
(3.2a) $\displaystyle\bm{\bm{u}}^{(1)}$
$\displaystyle=(u^{(1)}_{x}(x)\sin{(my+nz)},u^{(1)}_{y}(x)\cos{(my+nz)},u^{(1)}_{z}(x)\cos{(my+nz)})$
(3.2b) $\displaystyle(p^{(1)},\Pi^{(1)},\Phi^{(1)})$
$\displaystyle=(p^{(1)}(x),\Pi^{(1)}(x),\Phi^{(1)}(x))\cos{(my+nz)}$ (3.2c)

We next give the linearized system of equations for the many amplitudes in
(3.2). The two divergent relations are

$\displaystyle B^{(1)}_{x,x}-(mB^{(1)}_{y}+nB^{(1)}_{z})=0$ (3.3a)
$\displaystyle-\rho^{(1)}(mv(x)+nw(x))+(\rho(x)u^{(1)}(x))_{,x}-\rho(mv^{(1)}+nw^{(1)})=0$
(3.3b)

while $\Pi^{(1)}$ is specified as

$\displaystyle\Pi^{(1)}=c^{2}\rho^{(1)}+f(x)B^{(1)}_{y}+g(x)B^{(1)}_{z}$ (3.4)

where $c^{2}$ is the sound speed $dp/d\rho$ at the corresponding value of
$\rho(x)$. Ohm’s law reads

$\displaystyle\Phi^{(1)}_{,x}$
$\displaystyle=v^{(1)}g+vB^{(1)}_{z}-w^{(1)}f-wB^{(1)}_{y}$ (3.5a)
$\displaystyle-m\Phi^{(1)}$ $\displaystyle=w(x)B^{(1)}_{x}-u^{(1)}g(x)$ (3.5b)
$\displaystyle-n\Phi^{(1)}$ $\displaystyle=-v(x)B^{(1)}_{x}+u^{(1)}f(x)$
(3.5c)

Finally the three momentum balance equations are

$\displaystyle u^{(1)}(mv+nw)\rho(x)+\Pi^{(1)}_{,x}$
$\displaystyle=B^{(1)}_{x}(mf+ng)$ (3.6a)
$\displaystyle\rho\left(u^{(1)}v^{\prime}-(mv+nw)v^{(1)}\right)-m\Pi^{(1)}$
$\displaystyle=B^{(1)}_{x}f^{\prime}(x)-(mf+ng)B^{(1)}_{y}$ (3.6b)
$\displaystyle\rho\left(u^{(1)}w^{\prime}-(mv+nw)w^{(1)}\right)-n\Pi^{(1)}$
$\displaystyle=B^{(1)}_{x}g^{\prime}(x)-(mf+ng)B^{(1)}_{z}$ (3.6c)

We introduce the notation

$\displaystyle\bar{f}=mf+ng,\quad\bar{v}=mv+nw$ (3.7)

We see that $\rho^{(1)}$ occurs only in (3.3b) and (3.4). It is convenient to
eliminate $\rho^{(1)}$ from the problem and we replace (24) by

$\displaystyle\bar{v}\Pi^{(1)}=c^{2}\left(p(mv^{(1)}+nw^{(1)})-(pu^{(1)})_{,x}\right)+f(x)B^{(1)}_{y}+g(x)B^{(1)}_{z}$
(3.4*)

Provided the series retains the structure in the angles given in (3.2) in all
orders $m=n=0$ is never possible in the two equations: Ohm’s law (3.5b) and
(3.5c). Thus, we may replace (3.5b),(3.5c) by the following

$\displaystyle-\bar{v}\Phi^{(1)}$ $\displaystyle=u^{(1)}\left(vg-
wf\right)=u^{(1)}\Phi^{(0)}_{,x}$ (3.5b*) $\displaystyle-\bar{f}\Phi^{(1)}$
$\displaystyle=\left(uf-vg\right)B^{(1)}_{x}=-B^{(1)}_{x}\Phi^{(0)}_{,x}$
(3.5c*)

We may carry out similar reductions of pressure balance and from (3.5b) and
(3.5c) we obtain

$\displaystyle\rho
u^{(1)}\bar{v}^{\prime}-\rho\bar{v}(mv^{(1)}+nw^{(1)})-(m^{2}+n^{2})\Pi=B^{(1)}_{x}\bar{f}^{\prime}-\bar{f}(nB^{(1)}_{y}-mB^{(1)}_{z})$
(3.6b*) $\displaystyle\rho
u^{(1)}(nv^{\prime}-mw^{\prime})-\rho\bar{v}(nv^{(1)}-mw^{(1)})=B^{(1)}_{x}(nf-
mg)-\bar{f}(nB^{(1)}_{y}-mB^{(1)}_{z})$ (3.6c*)

The $x$ component of (3.6a) has the simpler form

$\displaystyle(vg-wf)\Pi^{(1)}_{,x}=(\bar{f}^{2}+\bar{g}^{2})\Phi^{(1)}$
(3.6a*)

At this point we have five equations (3.3a,3.5a,3.4*) and (3.6b*,3.6c*) in
which we may eliminate $u^{(1)}$ and $\bm{B}^{(1)}_{x}$ by means of
(3.5b*,3.5c*) and with the unknowns
$(\Phi^{(1)},\Phi^{(1)}_{,x},v^{(1)},w^{(1)},B^{(1)}_{y},B^{(1)}_{z})$ and
$\Pi^{(1)}$. Hence if we consider $\Phi^{(1)}$ and $\Pi^{(1)}$ as given then
we may solve for $\Phi^{(1)}_{,x}$. Thus, we obtain a first order homogeneous
linear differential equation for $\Phi^{(1)}$ of the form

$\Phi^{(1)}_{,x}=a(x)\Phi^{(1)}+b(x)\Pi^{(1)}.$

However, the construction fails with the vanishing of the determinant of the
system of five equations for
$(\Phi^{(1)}_{,x},v^{(1)},w^{(1)},B^{(1)}_{y},B^{(1)}_{z})$ in terms of
$\Phi^{(1)}$ and $\Pi^{(1)}$. The zeros of this determinant constitute the
resonant singularities of the system. We then obtain the determinant of this
system by replacing $B^{(1)}_{,x}$ and $u^{(1)}$ by $\Phi^{(1)}$ and ignoring
the terms proportional to $\Phi^{(1)}$ or $\Pi^{(1)}$. We find

$\displaystyle\Delta=\begin{vmatrix}\frac{\bar{f}}{\Phi^{(0)}_{,x}}\quad\quad
0\quad\quad\quad 0\quad\quad-m\quad\quad-n\\\ 1\quad\quad-g\quad\quad\quad
f\quad\quad\quad w\quad\quad-v\\\ 0\quad\quad\quad\rho\bar{v}\quad\quad\quad
0\quad\quad-\bar{f}\quad\quad\quad 0\\\ 0\quad\quad\quad
0\quad\quad\quad\rho\bar{v}\quad\quad\quad 0\quad\quad-\bar{f}\\\ -\frac{\rho
c^{2}\bar{v}}{\Phi^{(0)}_{,x}}\quad\rho c^{2}m\quad\>\>\rho c^{2}n\quad\quad
f\bar{v}\quad\quad g\bar{v}\end{vmatrix}$ (3.10)

and after straightforward reduction we obtain

$\displaystyle\Delta=(\rho\bar{v}^{2}-\bar{f}^{2})(\rho\bar{v}^{2}(c^{2}+(f^{2}+g^{2})/\rho)-c^{2}\bar{f}^{2})$
(3.11)

For the particular background state, resonances are associated with Alfvén
waves, the first factor, and a particular form of the fast-slow wave, the
second factor.

It is clear that for a general equilibrium there will be resonances when
$\bar{v}=\pm\bar{f}/\sqrt{\rho}$ or
$\pm\bar{f}c/\sqrt{c^{2}\rho+(f^{2}+g^{2})}$. For wide ranges of the ratio of
the mode indices, $m/n$ resonant surfaces will appear. Just as it was possible
to find formal expansions in mode amplitudes of ideal MHD equilibrium, it may
be of interest to explore the possibility of formal expansions of flows with
no resonant singularities.

In a more mathematically precise language we can state the main idea as
follows. One starts with a homogeneous linear system of five equations in
seven unknowns of the form $AX=0$, where $X$ is a 7-component vector and $A$
is a 5x7 rectangular matrix. Next one write $X=(x,y)$, where $x$ denotes the
first 5 components of X, and $y$ denotes the last two components, $\Phi^{(1)}$
and $\Pi^{(1)}$. Then the value of $y$ is fixed, which leads to the
inhomogeneous linear system $ax=b$, where $a$ is the restriction of $A$ to the
subspace $y=0$ and $b=A(0,y)$. When $\text{det}(a)$ is non-zero, one therefore
gets expressions for each component of $x$ that are linear in $y$. When
$\text{det}(a)=0$, there is no solution for $x$ in general unless $b$ is in
the image of $a$. The quantity $\Delta$ in eq. (3.10) is just that
$\text{det}(a)$.

Such a possibility depends on arranging the state so that it is in resonance
everywhere or nowhere. Note that when $\Delta=\text{det}(a)=0$, there is still
a solution for $x$ provided $y$ can be adjusted so that $A(0,y)$ is in the
image of $a$. Such a situation is possible for parallel flows, i.e., flows
where $\bm{B}$ and $\bm{u}$ are parallel, or

$\displaystyle\bm{u}=(0,v(x),w(x))=\lambda(x)(0,B_{y}(x),B_{z}(x)).$ (3.12)

In this case $\bar{v}=\lambda\bar{f}$, so that (3.11) becomes

$\displaystyle\Delta=\bar{f}^{2}\left(\rho(\lambda(x))^{2}-1\right)\left(\rho(\lambda(x))^{2}(c^{2}+(f^{2}+g^{2})/\rho)-c^{2}\right)=0$
(3.13)

If $\lambda(x)$ is chosen so that the second or the third term vanishes, then
$\Delta$ is identically zero for every mode. Such flows may well be of some
interest, but they seem extremely pathological and do not appear to be
realizable. However, the case $\bar{f}=0$ is of more interest and similar to
the equilibrium case. In particular we assume the second factor never vanishes
and $\bar{f}=0$ for $m=M,n=N$ where $M$ and $N$ are relatively prime. Thus, we
assume

$\displaystyle\bm{B}^{(0)}=(0,N,-M)\mu(x)$ (3.14)

where $\mu(x)$ is an arbitrary function, while $\bm{u}$ is given by
(3.12,3.14). We expand in amplitude about this state and describe expansions
to all order in the amplitude.

##  4 Nearly parallel flows

In the last section we showed that perturbations of steady flows are typically
subject to appearances of singularitiies on particular surfaces. However, for
parallel flows given by (3.12,3.14) singular surfaces need not appear. In this
section we examine the flows with these properties and we lay out the argument
that one may construct a formal series solution to all orders without the
appearance of any singularities when the lowest order system is a parallel
flow. We do not give the full details of the proof of expansions to all orders
but we consider the conclusion valid. The process closely parallels the
development in (Weitzner, 2014), where such a formalism is fully developed.

The expansion parameter of the series is the amplitude of the flow and field
components. We assume that the lowest order field and flow state is given by
(3.12,3.14). In first order we add fields and flows which is a sum of terms
given by (3.2) subject to the condition that none of them is resonant. We may
then continue to construct an expansion order by order. The series we
construct are assumed to have the same structure of angular dependence in all
orders. The linearized system used to determine the characteristics of the
system, closely parallels (2.4).

In this expansion at each order we must solve a system similar to
(3.3a),(3.4*),(3.5a), (3.5b*),(3.5c*),(3.6a),(3.6b*),(3.6c*) where sums of
products of lower order terms may appear added to the right hand side of the
equations. When the mode is not resonant so that $\Delta$ given by (3.13) is
not zero, one may solve the linear homogeneous system. Our task is to indicate
how one should be able to select series so that order by order there is no net
resonant term arising in the combination of sums of products of lower orders
which form the inhomogeneous terms of the generalization of the system.

Before we address the central issue in this work we must characterize the
structure of the resonance more fully. In resonance the quantities $\bar{v}$
and $\bar{f}$ are identically zero. Thus, from (3.5b*),(3.5c*), (3.6a*) and
(3.8) we find that

$\displaystyle u^{1}(x)$ $\displaystyle=B_{x}^{(1)}=0$ (4.1a)
$\displaystyle\Phi^{(1)}$ $\displaystyle=0\quad\quad(\text{from
\eqref{VbfbPhiOne}})$ (4.1b) $\displaystyle\Pi^{1}$
$\displaystyle=0\quad\quad(\text{from \eqref{26as}}).$ (4.1c)

Finally, the system is closed with the two conditions

$\displaystyle mB^{(1)}_{y}+nB^{(1)}_{z}$ $\displaystyle=0$ (4.2a)
$\displaystyle mv^{(1)}+nw^{(1)}$ $\displaystyle=0$ (4.2b)

The state given by (4.1), (4.2) is the resonant mode in any order, provided

$\displaystyle\frac{m}{n}=\frac{M}{N}.$ (4.3)

We finally turn to the inhomogeneous system where we expand the solution to
some order in the amplitude say order $P$. We obtain a system of the form
(3.3a),(3.4*), (3.5),(3.6) where the index one is replaced by $P$ and
inhomogeneous terms are added to the right hand sides of the equations. These
terms are sums of products of terms of order lower than $P$. We can solve the
system provided that there is no net contribution at a resonant term, for if
there were such a term we could not guarantee the construction of a periodic
solution of the system. We must show that the series can be arranged so that
no net inhomogeneous term is present.

We modify the structure of the solution in orders $(P-2)$ and $(P-1)$ so that
no net singular terms arise in order $P$. We assume that before modification
there are terms in order $P$ with $m$ and $n$ satisfying (4.3). We add a
resonant term in order $(P-2)$ with angular dependence $\exp{i(\mu y+\nu z)}$
and undetermined $x$ dependence. There are two such independent additions
$B_{y}^{(P-2)}$ and $v^{(P-2)}$, say. We assume that there are non-resonant
terms in the first order of the form (3.2). When these terms beat against each
other there will be terms in order $P-1$ with angular dependence
$\exp{i((\mu\pm m)y+(\nu\pm n)z)}$. Finally, in order $P$, terms from order
$(P-1)$ beat again against the first order terms with angular structure
$\exp{i(my+nz)}$. We choose the as yet undetermined $x$ dependence so as to
satisfy the two constraints. The quantities
$u^{(P)},B_{x}^{(P)},\Phi^{(P)},\text{ and }\Pi^{(P)}$ which were all zero in
the solution of the homogeneous data problem are all non-zero and are
determined.

##  5 Analysis of steady flows that admit an additional symmetry vector

We have so far outlined how one can perturbatively construct non-symmetric
steady plasma flows in ideal MHD. We considered flows that are nearly parallel
to avoid resonances. Compared to non-symmetric MHD equilibrium without flows,
the steady flow system (2.1) is far more complicated due to extra equations
and nonlinearities from the flow variables. In particular, it is not clear if
the resonances can be avoided. Fortunately, analytical progress can be made if
the steady flows possess a symmetry vector, the Hameiri vector (Hameiri,
1983), which is closely related to the quasisymmetry vector (Burby et al.,
2020a; Rodriguez et al., 2020). In the following, we employ a more formal
approach to investigate non-symmetric 3D flows that come equipped with the
Hameiri vector. We shall now state our goals and the main results obtained in
this section.

Our primary objective is to investigate the properties of a class of non-
symmetric 3D flow system that has flow components both parallel and
perpendicular to the magnetic field and admits Hameiri’s symmetry vector.
Under an additional assumption that the density does not change along this
vector, we show that non-symmetric generalizations of Bernoulli’s law, angular
momentum conservation, a generalized Grad-Shafranov (GGS) equation, and
associated generalized Hamada conditions can be obtained. Therefore, our main
results show that the description of these special non-symmetric flows
parallels symmetric flows.

The organization of this section is as follows. We begin with a discussion of
the Hameiri vector and point out its close connections with the quasisymmetry
vector in section 5.1. We then discuss the key assumptions that we make in
order to carry out the analysis in 5.2. We derive the generalizations of
Bernoulli’s law in section 5.3. We show that a generalized Grad-Shafranov
equation can be constructed following the approach of Burby et al. (2020a) in
section 5.4. We discuss the generalization of Hamada conditions in section 5.5
and summarize the results in section 5.6.

###  5.1  Hameiri’s vector C and weak quasisymmetry

We briefly summarize a key result due to Hameiri, which is essential to our
work. Assuming a steady flow with nested toroidal flux surfaces labelled by
$\psi$, Hameiri (1983) showed that ideal MHD allows an additional cross-
helicity of the form $\int d^{3}\bm{r}\>\bm{u}\cdot\bm{C}$, where
$d^{3}\bm{r}$ denotes volume integral over 3D space, if there exists a vector
$\bm{C}$ such that

$\displaystyle\bm{C}\cdot\bm{\nabla}\psi=0,\quad\bm{\nabla}\cdot\bm{C}=0,\quad\bm{\nabla}\times\left(\bm{B}\times\frac{1}{\rho}\bm{C}\right)=0,\quad\bm{\nabla}\times\left(\bm{u}\times\bm{C}\right)=0.$
(5.1)

Hameiri has shown that in axisymmetry, $\bm{C}$ is $R^{2}\bm{\nabla}\varphi$
in the usual $(R,Z,\varphi)$ cylindrical coordinates used in tokamak
literature. Furthermore, from (5.1) we can see that the vector $\bm{C}$ is
analogous to a magnetic field since it is frozen-in with the flow and lies on
flux surfaces.

We shall make an additional assumption that the density is constant in the
direction of the Hameiri vector i.e. $\bm{C}\cdot\bm{\nabla}\rho=0$. Writing
$\bm{C}=\rho\bm{Q}$, we get an equivalent set of conditions on $\bm{Q}$,

$\displaystyle\bm{Q}\cdot\bm{\nabla}\psi=0,\quad\bm{\nabla}\cdot(\rho\bm{Q})=0,\quad\bm{\nabla}\times\left(\bm{B}\times\bm{Q}\right)=0,\quad\bm{\nabla}\times\left(\rho\bm{u}\times\bm{Q}\right)=0.$
(5.2)

The conditions on $\bm{Q}$ from (5.2) reads

$\displaystyle\bm{Q}\cdot\bm{\nabla}\rho=0$ (5.3a)
$\displaystyle\bm{B}\times\bm{Q}=\bm{\nabla}\psi$ (5.3b)
$\displaystyle\bm{\nabla}\cdot\bm{Q}=0,$ (5.3c)

while the conditions on $\bm{u}$ from (5.2) and (2.1) are

$\displaystyle\bm{\nabla}\cdot(\rho\bm{u})=0$ (5.4a)
$\displaystyle\bm{u}\times\bm{B}=\bm{\nabla}\Phi(\psi)$ (5.4b)
$\displaystyle\bm{\nabla}\times\left(\rho\bm{u}\times\bm{Q}\right)=0.$ (5.4c)

We want to point out that the conditions given in (5.3) are closely related to
the “weak quasisymmetry” conditions obtained in (Rodriguez et al., 2020). If
instead of $\bm{Q}\cdot\bm{\nabla}\rho=0$, we assume
$\bm{Q}\cdot\bm{\nabla}B=0$ and $\rho=\rho(\psi,B)$ we get exactly the weak
quasisymmetry (QS) conditions. This close connection of ideal MHD flows with
quasisymmetry was pointed out earlier by Helander (Helander, 2007, 2014).

The analogy with “weak QS” becomes more evident when we look at the internal
consistency of (5.3). For given $\bm{B}$ and $\rho$, the system (5.3) is an
overdetermined system for $\bm{Q}$ as it represents four constraints for the
three components of $\bm{Q}$. From (5.3b) we get

$\displaystyle\bm{Q}=\frac{1}{B^{2}}\left(\left(\bm{B}\cdot\bm{Q}\right)\bm{B}-\bm{B}\times\bm{\nabla}\psi\right).$
(5.5)

The component $\left(\bm{B}\cdot\bm{Q}\right)$ can be chosen such that (5.3a)
is satisfied. The divergence-free condition (5.3c) then imposes the following
condition on $\bm{B}$

$\displaystyle\bm{J}\cdot\bm{\nabla}\psi=\bm{B}\cdot\bm{\nabla}\left(\bm{B}\cdot\bm{Q}\right)-\bm{Q}\cdot\bm{\nabla}B^{2},$
(5.6)

which must be satisfied in order for (5.3) to have a consistent solution. Such
a consistency condition indeed appears in QS systems (Rodriguez et al., 2020;
Burby et al., 2020a) without the last term since $\bm{Q}\cdot\bm{\nabla}B=0$.

If we assume that $\bm{Q}$ is given we can construct a $\bm{B}$ consistent
with $\bm{Q}$. Using (5.3b), we can write $\bm{B}$ in terms of $\bm{Q}$ as

$\displaystyle\bm{B}=\frac{1}{Q^{2}}\left(\left(\bm{B}\cdot\bm{Q}\right)\bm{Q}+\bm{Q}\times\bm{\nabla}\psi\right),$
(5.7)

which, is analogous to the symmetry flux coordinates in a tokamak (D’haeseleer
et al., 1991) and the form used by Burby et al. (2020a) to study the strong
form of QS. To ensure that $\bm{\nabla}\cdot\bm{B}=0$ is satisfied by (5.7),
we must have

$\displaystyle\bm{Q}\cdot\bm{\nabla}\left(\bm{B}\cdot\bm{Q}\right)+\bm{\nabla}\cdot\left(\bm{Q}\times\bm{\nabla}\psi\right)-\bm{B}\cdot\bm{\nabla}Q^{2}=0.$
(5.8)

In the following, we shall assume that there exists a Hameiri vector of the
form $\bm{C}=\rho\bm{Q}$, where $\bm{C}\cdot\bm{\nabla}\rho=0$ and $\bm{Q}$
satisfies (5.3) and (5.8) and is given.

###  5.2 Underlying assumptions

Before proceeding further, we shall discuss the various assumptions that we
make in our analysis. We begin with our fundamental assumption that there
exist non-symmetric plasma flows that possess nested flux surfaces. In our
formal exploration, we are not forced to make any specific assumptions about
the rotational transform of the magnetic field or the magnetic shear. The
nestedness of flux surfaces is not guaranteed but is often made in the
literature. In the following, we shall derive the necessary conditions (the
Hamada conditions) required to suppress current singularities on the rational
surfaces and preserve the nested surfaces. It will be apparent that the Hamada
conditions can not possibly be satisfied in a toroidal volume unless magnetic
shear is weak. With the analysis of previous sections in mind, we now make a
conscious choice of weak magnetic shear in the following.

Our second important assumption is the existence of Hameiri’s symmetry vector
$\bm{C}$. Theoretically, $\bm{C}$ can be obtained by applying Noether’s
theorem to the one-fluid ideal MHD Lagrangian (Ilgisonis & Pastukhov, 2000).
In particular, $\bm{C}$ is related to the relabelling symmetry of ideal MHD
(Ilgisonis & Pastukhov, 2000; Hameiri, 1998; Vladimirov & Moffatt, 1995).
However, the system of equations determining $\bm{C}$, (5.1), being a
nonlinear overdetermined system, it is difficult to prove the existence of
such a vector. Moreover, unlike QS, a systematic study of $\bm{C}$ has not
been carried out to the best of our knowledge.

As discussed in details in (Hameiri, 1998), an important necessary condition
that needs to be satisfied in order for $\bm{C}$ to exist in a toroidal domain
(Hameiri, 1983) is that on each closed magnetic field line,

$\displaystyle\oint\frac{d\ell}{B}\rho=m(\psi).$ (5.9)

Here, $m(\psi)$ is assumed to be a smooth function and $d\ell$ denotes the
differential element along the magnetic field. If condition (5.9) is not
satisfied then the flow must be nearly-parallel (Hameiri, 1998).

We then assumed that $\bm{C}\cdot\bm{\nabla}\rho=0$, which implies that the
density (and hence pressure from (2.2) ) possess a continuous symmetry along
$\bm{C}$. This choice has been made solely for simplicity and analytical
tractability. Let us briefly discuss some of the consequences of the
assumption of $\bm{Q}\cdot\bm{\nabla}\rho=0=\bm{Q}\cdot\bm{\nabla}p$. Firstly,
we note that (5.9) is satisfied identically. Secondly, if $\bm{Q}$ lines do
not close on themselves then density and pressure must not vary on a flux
surface i.e. $p=p(\psi),\rho=\rho(\psi)$. Therefore, such systems must be
subsonic (Tasso & Throumoulopoulos, 1998) since large centrifugal forces lead
to variations of density and pressure on flux surfaces. When $\bm{Q}$ lines
are closed, flows do not have to be subsonic, and pressure variations on flux
surfaces are allowed. Although we shall not make any assumption on the closure
of $\bm{Q}$ lines, it is interesting to note that in QS, the symmetry lines
close helically or toroidally depending on the type of QS.

The final assumption that $\bm{Q}$ satisfies (5.8), is due to the fact that
(5.7) is divergence-free only in axisymmetry (Burby et al., 2020a). The
condition (5.8) is not really an additional constraint on $\bm{Q}$ (or
$\bm{C}=\rho\bm{Q}$), since it follows directly from (5.3b). Any $\bm{Q}$
satisfying (5.3b) with a divergence free physical magnetic field necessarily
satisfies (5.8). We have included the condition (5.8) only because of its
usage in subsequent calculations.

###  5.3  Generalized Bernoulli’s law and angular momentum

It is well known (Hameiri, 1983; Throumoulopoulos et al., 2006) that five
scalar functions are needed to describe the steady axisymmetric ideal MHD
flows. Out of these, two functions are given by the electrostatic potential
$\Phi^{\prime}(\psi)$ and the entropy (assumed constant). The remaining three
functions, denoted by $(\Lambda(\psi),H(\psi),I(\psi))$, are associated with
the parallel component of the flow, the Bernoulli law, and angular momentum,
respectively. In the following, our goal is to derive these quantities
systematically for non-symmetric flows with the help of the Hameiri vector
$\bm{C}=\rho\bm{Q}$.

For a given magnetic field $\bm{B}$ and $\bm{Q}$ that satisfy (5.3) and (5.6),
it can be shown (Hameiri, 1998) that a steady state flow satisfying (5.4) is
given by

$\displaystyle\bm{u}=\frac{\Lambda(\psi)}{\rho}\bm{B}-\Phi^{\prime}(\psi)\bm{Q},$
(5.10)

which clearly resembles the form of the axisymmetric steady flow of a tokamak
(Hameiri, 1983). We also get the flux function $\Lambda(\psi)$ from (5.10).

From (5.10), we note that any two vectors from $\bm{u},\bm{B}$ and $\bm{Q}$
can be used as basis vectors on a flux surface. We shall choose $\bm{B}$ and
$\bm{Q}$ since the commutator of the two derivatives
$\bm{B}\cdot\bm{\nabla},\bm{Q}\cdot\bm{\nabla}$ vanish (as shown in Appendix
B). Choosing $\bm{B},\bm{Q},\bm{\nabla}\psi$ as basis vectors, we shall now
obtain the components of the ideal MHD force balance equation (2.1b).

We find it convenient to rewrite (2.1b) in the form of a vorticity equation

$\displaystyle\bm{\nabla}\left(\frac{1}{2}u^{2}+h\right)=\bm{u}\times\left(\bm{\nabla}\times\bm{u}\right)+\frac{1}{\rho}\bm{J}\times\bm{B}.$
(5.11)

Here, we have defined $h$ such that $\bm{\nabla}h=(1/\rho)\bm{\nabla}p(\rho)$,
and we have used the standard vector identity
$\bm{u}\cdot\bm{\nabla}u=\bm{\nabla}(u^{2}/2)-\bm{u}\times(\bm{\nabla}\times\bm{u})$.
Dotting with $\bm{B}$ and $\bm{Q}$ respectively, and using ((5.4b),(5.10)) we
get

$\displaystyle\bm{B}\cdot\bm{\nabla}\left(\frac{1}{2}u^{2}+h\right)+\bm{\nabla}\cdot\left(\Phi^{\prime}(\psi)\bm{u}\times\bm{\nabla}\psi\right)=0$
(5.12a)
$\displaystyle\rho\bm{Q}\cdot\bm{\nabla}\left(\frac{1}{2}u^{2}+h\right)-\bm{J}\cdot\bm{\nabla}\psi+\bm{\nabla}\cdot\left(\Lambda(\psi)\bm{u}\times\bm{\nabla}\psi\right)=0.$
(5.12b)

To evaluate the divergence terms that appear in (5.12) we use (5.4b) and
(B.2). We note that the assumption $\bm{Q}\cdot\bm{\nabla}\rho=0$ helps us to
simplify (5.12b).

Simplifying (5.12) using (5.3a), (5.6) and (B.2) we get

$\displaystyle\bm{B}\cdot\bm{\nabla}\left(\frac{1}{2}u^{2}+h+\Phi^{\prime}(\psi)\left(\bm{u}\cdot\bm{Q}\right)\right)-\bm{Q}\cdot\bm{\nabla}\left(\left(\bm{u}\cdot\bm{B}\right)\Phi^{\prime}(\psi)\right)=0$
(5.13a)
$\displaystyle\bm{B}\cdot\bm{\nabla}\left(\left(\bm{B}-\Lambda(\psi)\bm{u}\right)\cdot\bm{Q}\right)-\bm{Q}\cdot\bm{\nabla}\left(\frac{1}{2}\rho
u^{2}+B^{2}-\Lambda(\psi)\left(\bm{u}\cdot\bm{B}\right)\right)=0.$ (5.13b)

In axisymmetry, the $\bm{Q}\cdot\bm{\nabla}$ terms in (5.13) vanish
identically, and we get two homogeneous magnetic differential equations.
Solving the homogeneous magnetic differential equations, we get two flux
functions which denote Bernoulli’s law and momentum conservation in the
symmetry ($\bm{Q}$) direction (Hameiri, 1983; Tasso & Throumoulopoulos, 1998).
The situation is similar here except that the magnetic differential equations
are non-homogeneous; therefore, we need to check for consistency conditions.
We postpone the discussion on the consistency conditions until section 5.5.

Assuming that the solvability conditions are satisfied, we solve (5.13) and
obtain the generalized Bernoulli and the generalized angular momentum
equations

$\displaystyle\frac{1}{2}u^{2}+h+\Phi^{\prime}(\psi)\left(\bm{u}\cdot\bm{Q}\right)=H(\psi)+\mathcal{H}\quad\quad\quad\text{(Bernoulli)}$
(5.14a)
$\displaystyle\left(\bm{B}\cdot\bm{Q}\right)-\Lambda(\psi)\left(\bm{u}\cdot\bm{Q}\right)=I(\psi)+\mathcal{I},\quad\quad\text{(angular
momentum)}$ (5.14b)

where, $\mathcal{H},\mathcal{I}$ defined by

$\displaystyle\mathcal{H}=\bm{Q}\cdot\bm{\nabla}\int\frac{d\ell}{B}\left(\bm{u}\cdot\bm{B}\right)\Phi^{\prime}(\psi),\quad\mathcal{I}=\bm{Q}\cdot\bm{\nabla}\int\frac{d\ell}{B}\left(\frac{1}{2}\rho
u^{2}+B^{2}-\Lambda(\psi)\left(\bm{u}\cdot\bm{B}\right)\right),$ (5.15)

are single-valued and doubly-periodic functions if Hamada conditions (see
(5.24) below) are satisfied. They vanish identically in the axisymmetric case.
As is well-known (Beskin, 2009; Hameiri, 1998, 1983) availability of conserved
quantities are very helpful in characterizing flow patterns.

###  5.4  Generalized Grad-Shafranov equation for flows with the Hameiri
vector

We are now in a position to obtain the GGS from the $\bm{\nabla}\psi\cdot$
component of the vorticity equation (5.11). The step-by-step derivation is
provided in Appendix D. The steps to derive the GGS are identical to the ones
in the derivation of the flow-modified classical Grad-Shafranov equation in
axisymmetry or helical symmetry. However, there are important differences
which we shall now discuss. Firstly, in the absence of axial/helical symmetry,
the vector

$\displaystyle\bm{w}=\left(\bm{\nabla}\times\bm{Q}\right)\times\bm{Q}+\bm{\nabla}Q^{2},$
(5.16)

which denotes the deviation of $\bm{Q}$ from axial/helicalsymmetry vector
(Burby et al., 2020a) appears in the GGS. We can check that (5.8), the
condition required to ensure $\bm{\nabla}\cdot\bm{B}=0$ when $\bm{B}$ is
expressed in terms of $\bm{Q}$ (equation (5.7)), is identical to

$\displaystyle\bm{Q}\cdot\bm{\nabla}\left(\bm{B}\cdot\bm{Q}\right)-\bm{B}\cdot\bm{w}=0.$
(5.17)

Obviously, if $\bm{Q}$ is an axial/helical symmetry both terms individually
vanish and $\bm{B}$ is automatically divergence-free. Secondly, the geometric
quantity

$\displaystyle\omega_{Q}\equiv\bm{Q}\cdot\bm{\nabla}\times\bm{Q},$ (5.18)

which doesn’t appear in the axisymmetric GS, but appears in the helically
symmetric GS, also appear in the GGS. Finally, terms appear with the force-
like quantity

$\displaystyle\bm{F_{Q}}=\bm{J}\times\bm{Q}+\bm{\nabla}\left(\bm{B}\cdot\bm{Q}\right),$
(5.19)

which vanishes identically in static MHD (Burby et al., 2020a).

With the definitions ((5.16),(5.18) and (5.19)) in mind, we can write the
flow-modified GGS in the following form

$\displaystyle\bm{\nabla}\cdot\left(\left(1-\frac{\Lambda^{2}}{\rho}\right)\frac{1}{Q^{2}}\bm{\nabla}\psi\right)+\rho\frac{dH}{d\psi}+$
$\displaystyle\left(\bm{u}\cdot\bm{B}\right)\frac{d\Lambda}{d\psi}+\frac{\left(\bm{B}\cdot\bm{Q}\right)}{Q^{2}}\frac{dI}{d\psi}-\rho\left(\bm{u}\cdot\bm{Q}\right)\frac{d\Phi^{\prime}}{d\psi}+\mathcal{N}=0,$
(5.20)

where,

$\displaystyle\mathcal{N}=\frac{I+\mathcal{I}}{Q^{2}}\left(\frac{\omega_{Q}}{Q^{2}}+\frac{\bm{F_{Q}}\cdot\bm{\nabla}\psi}{|\bm{\nabla}\psi|^{2}}\right)$
$\displaystyle+\left(\rho\bm{\nabla}\mathcal{H}+\frac{\left(\bm{B}\cdot\bm{Q}\right)}{Q^{2}}\bm{\nabla}\mathcal{I}\right)\cdot\frac{\bm{\nabla}\psi}{|\bm{\nabla}\psi|^{2}}$
$\displaystyle+\frac{\bm{w}\cdot\bm{\nabla}\psi}{Q^{2}}\left(\frac{1}{Q^{2}}\left(1-\frac{\Lambda^{2}}{\rho}\right)-\frac{\rho\Phi^{\prime}\left(\bm{u}\cdot\bm{Q}\right)}{|\bm{\nabla}\psi|^{2}}\right).$
(5.21)

In axisymmetry (Hameiri, 1983), $\mathcal{N}=0$. For helical symmetry, an
additional term of the form $(I\omega_{Q})/Q^{4}$ appears. To study force-
balance other than ideal MHD force-balance, we can replace $\bm{F_{Q}}$ by
other forces (Rodriguez et al., 2021). For the case of flows only in the
symmetry direction $\Lambda(\psi)=0$ and for strictly parallel flows
$\Phi^{\prime}(\psi)=0$. One can easily extend the GGS to include other non-
relativistic forces.

Unlike in symmetric geometries, the coefficients and derivatives in the flow-
modified GGS depend on all three spatial variables (Constantin et al., 2020;
Burby et al., 2020a). Hence, we need to impose the additional condition

$\displaystyle\bm{Q}\cdot\bm{\nabla}\psi=0.$ (5.22)

Whether such solutions to the GGS can be constructed is a challenging open
problem.

###  5.5  The generalized Hamada conditions

We now return to the consistency conditions that must be satisfied in order
for (5.13) to have non-singular smooth solutions. When the rotational
transform is irrational, we can flux-surface average (5.13). Since both
$\bm{B}$ and $\bm{C}$ are tangential to the flux surfaces, there is no
inconsistency. On the other hand, for rational rotational transform, the
vanishing of the two equations under closed field line $\oint d\ell/B$
integral is not automatic and leads to nontrivial constraints (Newcomb,
1959)\- the so-called Hamada conditions (Hamada, 1962; Helander, 2014).

We recall that in ideal MHD equilibrium with scalar pressure $p(\psi)$, the
integral constraints that must be satisfied on rational surfaces are (Grad,
1967; Solov’ev & Shafranov, 1970; Grad, 1971)

$\displaystyle\oint\frac{d\ell}{B}=c_{1}(\psi),\quad\oint\frac{d\ell}{B}B^{2}=\oint\bm{B}\cdot\bm{d\ell}=c_{2}(\psi),\quad
c^{\prime}_{2}(\psi)+p^{\prime}(\psi)c_{1}(\psi)=0,$ (5.23)

where, $c_{1}(\psi),c_{2}(\psi)$ are single-valued continuous functions of
$\psi$ that reduce to constants in the vacuum limit. We give a straightforward
proof of these conditions in Appendix C.

Hamada conditions (5.23) are identically satisfied in symmetric geometries but
are not satisfied in general for non-symmetric geometry with arbitrary
pressure and rotational transform profiles (Boozer, 1981; Weitzner, 2016),
leading to singular currents on rational surfaces (Loizu et al., 2015b, a). In
the following, we shall obtain the generalized Hamada conditions for steady
flows with Hameiri’s vector $\bm{C}=\rho\bm{Q}$.

We note that (5.9) is already in the form of a Hamada condition. To find the
other conditions, we carry out the $\oint d\ell/B$ integrals of (5.13)
(denoted by $\langle\quad\rangle$ brackets ). Using the commutation of the two
derivatives to get two homogeneous equations of the form
$\bm{Q}\cdot\bm{\nabla}\langle A\rangle=0$, which implies that $A$ must be a
flux function. Thus, we obtain two consistency conditions

$\displaystyle\langle\bm{u}\cdot\bm{B}\rangle$
$\displaystyle=\oint\bm{u}\cdot\bm{d\ell}=C_{1}(\psi)$ (5.24a)
$\displaystyle\bigg{\langle}B^{2}+\frac{1}{2}\rho
u^{2}-\Lambda\left(\bm{u}\cdot\bm{B}\right)\bigg{\rangle}$
$\displaystyle=\oint\bm{B}\cdot\bm{d\ell}+\oint\frac{d\ell}{B}\left(\frac{1}{2}\rho
u^{2}\right)-\Lambda(\psi)C_{1}(\psi)=C_{2}(\psi),$ (5.24b)

which show that the circulations of $\bm{u}$ and $\bm{B}$ along the closed
field line are required to be flux functions in order that (5.13) be solvable.

The final Hamada condition is obtained by closed field line averaging of
(5.20) on a rational flux surface, which leads to the following equation for
$C_{2}^{\prime}(\psi)$ (see Apendix E)

$\displaystyle
C^{\prime}_{2}(\psi)+C_{1}(\psi)\frac{d\Lambda}{d\psi}=H\partial_{\psi}\langle\rho\rangle-\Phi^{\prime}\partial_{\psi}\langle\rho\left(\bm{u}\cdot\bm{Q}\right)\rangle-\langle\mathcal{N}\rangle+\partial_{\psi}\bigg{\langle}\rho\mathcal{H}+\frac{\mathcal{I}}{Q^{2}}\left(\bm{B}\cdot\bm{Q}\right)\bigg{\rangle}.$
(5.25)

Equations (5.24) and (5.25), together with (5.9) describe integral constraints
that need to be satisfied on each rational flux surface. These are the
generalized Hamada conditions that include the effect of steady flows. In
Appendix C we show that when these consistency conditions are satisfied
magnetic resonances on rational surfaces are avoided. Given the complicated
nature of the generalized Hamada conditions, it is doubtful that they will be
satisfied by a generic steady flow. Furthermore, since the Hamada conditions
need to be satisfied on each rational surface, the magnetic shear needs to be
sufficiently weak to avoid passage through multiple low-order rational
surfaces.

###  5.6 Summary of section 5

In section 5, we discussed a special class of non-symmetric steady flows
endowed with Hameiri’s symmetry vector $\bm{C}=\rho\bm{Q}$, where
$\bm{C}\cdot\bm{\nabla}\rho=0$ and $\bm{Q}$ satisfy the overdetermined system
(5.3). If the integral constraints (5.9), (5.24) and (5.25), are satisfied
then the description of such flows involve a GGS (5.20) and four flux-
functions similar to axisymmetric flows (Hameiri, 1983; Tasso &
Throumoulopoulos, 1998) namely,
$\Lambda(\psi),\Phi^{\prime}(\psi),H(\psi),I(\psi)$. Note that only the first
four appear in the GGS since we have assumed barotropic equation of state with
constant entropy. The flux function $m(\psi)$ is needed to ensure existence of
Hameiri’s symmetry vector. The integral constraints (5.24) are necessary to
avoid current singularities (see Appendix C). They furnish two more flux-
functions $C_{1}(\psi),C_{2}(\psi)$ which measure the circulations of $\bm{u}$
and $\bm{B}$ in a closed field line system. Physically, they amount to the
constraint that the flows should not lead to accumulation of charges inside
the closed magnetic flux tubes on rational surfaces. Finally, non-symmetric
geometry introduces several extra terms grouped together as $\mathcal{N}$ such
as $\mathcal{H},\mathcal{I},\bm{w}\cdot\bm{\nabla}\psi$ etc.

##  6 Discussion

We studied several classes of non-symmetric ideal MHD equilibrium with flows
larger than the diamagnetic flows. Such flows could occur in stellarators with
We showed that the techniques developed in carrying out perturbation methods
to all orders for static MHD equilibrium Weitzner (2014) can also be extended
to MHD equilibrium with flows. The basic idea is that if the field lines are
closed, i.e., the rotational transform is rational and the magnetic shear is
weak, one can systematically eliminate magnetic resonances at each order by
utilizing the “free-function” that one gets by solving the magnetic
differential equation from lower orders. Plasma flows introduce extra
resonances, which are absent in static MHD. It is, in general, very
complicated to eliminate all such resonances. However, for the class of
equilibrium with nearly parallel flows, one can do so. It is to be noted that
exact solutions with large parallel flows can be constructed (Kamchatnov,
1982) when the flow is considered incompressible. We have argued that
compressible, nearly parallel flows are also possible in the MHD model of
plasma.

Although the present work on the analysis of magnetic resonances has been
carried out on rectangular coordinates with periodic boundary conditions,
extensions to polar and toroidal coordinates should be fairly direct. The
principal issue is behavior near the magnetic axis. To deal with this
complication, one must require that components of the fields and flows vanish
sufficiently rapidly near the axis. The interactions through the nonlinear
terns typically do not destroy these properties. We leave the detailed near-
axis analysis for the future.

We have studied a particular class of perpendicular flows by taking a more
formal approach. Our approach is based on utilizing a symmetry vector found by
Hameiri (Hameiri, 1981, 1998). The generalizations of Bernoulli’s law and
conservation of angular momentum were also obtained thanks to Hameiri’s
symmetry vector. We obtain a flow-modified generalized Grad-Shafranov equation
for such flows and the constraints resulting from the closed magnetic field
lines. It is clear from the nature of the constraints that if such non-
parallel flows exist, they must be special. There exists a close connection
between Hameiri’s vector and the quasisymmetry vector that will be further
discussed elsewhere (Rodriguez et al., 2021).

We have not attempted to estimate the damping of the flows due to neoclassical
effects, which is an essential but difficult question (Simakov & Helander,
2009). We shall address this in detail in a forthcoming paper (Rodriguez et
al., 2021). We only make a few observations here. First, the neoclassical
effects are tied to the parallel electric field, which does not appear in
ideal MHD. If the parallel electric field can be shown to be self-consistently
small or higher-order in the Larmor radius, then the MHD description prevails,
and flow damping would be minimized. Therefore, our results are valid within
the framework of ideal MHD. Second, the techniques developed here can be
easily extended to more sophisticated MHD models that are relevant to
astrophysical plasmas.

Our results here do not necessarily contradict earlier results (Simakov &
Helander, 2009; Sugama et al., 2011) that argue that large flows can not be
supported even if the stellarator is quasisymmetric but not completely
axisymmetric. We show that it is challenging to avoid magnetic resonances
unless the flow is nearly parallel. Except for very special flows, the
generalized Hamada conditions on rational surfaces can not be satisfied,
severely limiting the possibility of steady non-symmetric flows in generic
stellarators. We expect that kinetic constraints will further restrict the
class of available solutions.

##  Acknowledgements

The authors would like to sincerely thank the anonymous reviewer for their
excellent, meticulous, and constructive suggestions that have significantly
improved the organization and presentation of the results in the paper. In
addition, WS would like to thank E. Rodriguez, A.Bhattacharjee, A.B.Hassam,
E.J.Paul, and J.Juno for helpful discussions. This research was partly funded
by the US DOE grant no. DEFG02-86ER53223 and Simons Foundation/SFARI (560651,
AB).

##  Appendix A Analysis of the dispersion relation (2.7)

Our goal here is to show that real solutions of the dispersion relation (2.7)
exist. We shall rewrite the dispersion relation in terms of dimensionless
physical quantities and analyze the parametric space of these parameters,
which support real solutions of (2.7). We shall also look at some simple,
special cases which will allow us to gain insight into the parametric space.

Let us consider a generic flow of the form

$\displaystyle\bm{u}=u_{\parallel}\hat{\bm{B}}+\bm{u_{\perp}},\quad\hat{\bm{B}}=\frac{\bm{B}}{B},\quad
u_{\parallel}=\bm{u}\cdot\hat{\bm{B}}.$ (A.1)

Associated with this flow, we can define the parallel and perpendicular Mach
numbers as follows

$\displaystyle
M_{\parallel}=\frac{u_{\parallel}}{c}=\frac{\bm{u}\cdot\bm{B}}{cB},\quad
M_{\perp}=\frac{u_{\perp}}{c}=\frac{|\bm{u}_{\perp}|}{c},$ (A.2)

$c$ being the sound speed.

As noted in the text, the dispersion relation is independent of the magnitude
of $\bm{k}$. Therefore, we can divide (2.7) by $(\rho c^{2}k^{2})^{2}$ to get

$\displaystyle(k_{u})^{4}-(k_{u})^{2}\left(1+\frac{1}{\beta}\right)+(k_{B})^{2}\frac{1}{\beta}=0,$
(A.3)

where,

$\displaystyle\frac{1}{\beta}=\frac{B^{2}}{\rho c^{2}},\quad
k_{u}=\bm{\hat{k}}\cdot\frac{\bm{u}}{c},\quad
k_{B}=\bm{\hat{k}}\cdot\hat{\bm{B}}.$ (A.4)

We will not consider any components of $\bm{\hat{k}}$ orthogonal to $\bm{u}$
and $\bm{B}$ since they trivially satisfy (A.3). Hence, we will assume that
the unit vector $\hat{k}$ is coplanar with $\bm{u},\bm{B}$ and therefore,
admits the following orthogonal decomposition

$\displaystyle\bm{\hat{k}}=\cos{\theta}\>\hat{\bm{B}}+\sin\theta\>\hat{\bm{u}}_{\perp},\quad\hat{\bm{u}}_{\perp}=\frac{\bm{u}_{\perp}}{u_{\perp}}.$
(A.5)

Furthermore, (A.4) together with (A.2) implies that

$\displaystyle k_{B}=\cos\theta,\quad
k_{u}=M_{\perp}\sin{\theta}+M_{\parallel}\cos{\theta}.$ (A.6)

To determine $\theta$ we substitute (A.6) into the dispersion relation (A.3).
which leads to the following quartic equation for $\xi=\tan\theta$

$\displaystyle a_{4}\xi^{4}+a_{3}\xi^{3}+a_{2}\xi^{2}+a_{1}\xi+a_{0}=0,$ (A.7)
$\displaystyle
a_{4}=M_{\perp}^{2}\left(M_{\perp}^{2}-\left(1+\frac{1}{\beta}\right)\right),\quad\quad
a_{0}=M_{\parallel}^{2}\left(M_{\parallel}^{2}-\left(1+\frac{1}{\beta}\right)\right)+\frac{1}{\beta},$
$\displaystyle
a_{3}=4M_{\parallel}M_{\perp}\left(M_{\perp}^{2}-\frac{1}{2}\left(1+\frac{1}{\beta}\right)\right),\quad
a_{1}=4M_{\parallel}M_{\perp}\left(M_{\parallel}^{2}-\frac{1}{2}\left(1+\frac{1}{\beta}\right)\right),$
$\displaystyle
a_{2}=6\left(M_{\parallel}M_{\perp}\right)^{2}-\left(1+\frac{1}{\beta}\right)\left(M_{\parallel}^{2}+M_{\perp}^{2}\right)+\frac{1}{\beta}.$

We note that there are only three physical parameters that come into the
problem: plasma beta, $\beta=\rho c^{2}/B^{2}$, parallel and perpendicular
Mach numbers $(M_{\parallel},M_{\perp})$ given by (A.2). In the following, we
shall assume that all three physical quantities are real and positive.

The theory of quartic equations is well-developed, and the necessary
conditions for a quartic equation to have real roots can be determined
straightforwardly (Rees, 1922). When the quartic discriminant is negative, the
equation has two real and two imaginary roots. We have plotted the associated
parametric region in Figure (1(a)) and (1(b)). The region where four real
roots can exist is shown in Figure (1(c)). These two regions are complimentary
and hence the equation always has a real root. The details are not
particularly illuminating and will be omitted here. Fortunately, a lot of
insight can be gained by studying some simple cases, which we shall do next.

(a) region with 2 real roots

(b) side-view

(c) region with 4 real roots

Figure 1: a) and b) The region of parameter space where (A.7) has 2 real and 2
imaginary roots as viewed from two sides. c) The region with four real roots.
Since the latter is complimentary to the former, the equation has at least one
real root. Along the lines $M^{2}_{\parallel}=1,1/\beta,1/(\beta+1)$ and
$\beta=1$ multiple roots co-exist (see text).

We begin with the special case of purely parallel flows i.e. $M_{\perp}=0$. In
this case, the equation simplifies considerably and we get

$\displaystyle\xi^{2}=\left(\frac{M_{\parallel}^{2}-\frac{1}{\beta}}{\frac{1}{\beta+1}-M_{\parallel}^{2}}\right)\frac{\left(1-M_{\parallel}^{2}\right)}{\left(1+\frac{1}{\beta}\right)}.$
(A.8)

As $M_{\parallel}^{2}\to 1/(\beta+1),\xi\to\infty$, and the angle
$\theta=\arctan{\xi}\to\pi/2$, which implies that $\bm{k}$ becomes orthogonal
to $\hat{B}$. On the other extreme, when $M_{\parallel}^{2}=1$ or
$M_{\parallel}^{2}=1/\beta$, $\xi=0$, i.e. $\bm{k}$ is parallel to $\hat{B}$.
The case $M_{\parallel}^{2}=1/\beta$ is discussed in (Kamchatnov, 1982). These
special solutions appear as lines of discontinuity in Figure (1(b)) and
(1(c)). The quartic determinant is zero along these curves leading to multiple
roots. The solutions corresponding to $M_{\parallel}^{2}=1/(\beta+1),1/\beta$,
for flows of the type (3.12) were described as pathological solutions in
(3.13) which lead to $\Delta=0$. Physically, the system can develop shocks
along these curves.

Omitting the three special cases discussed above, we find that real solutions
with $\xi>0$ exist for $M_{\parallel}$ less (greater) than one, provided
$M_{\parallel}^{2}$ lies inside (outside) of the interval
$(1/(\beta+1),1/\beta)$. We note that for large values of $\beta$ subsonic
parallel ($M_{\parallel}^{2}<1$) flows are rare compared to the supersonic
parallel flows, in accordance with Figure (1(a)).

Next, we consider the case of purely perpendicular flows, i.e.,
$M_{\parallel}=0$. If

$\displaystyle M_{\parallel}=0\quad\text{and}\quad
M_{\perp}^{2}=1+\frac{1}{\beta}\quad,\quad\xi^{2}=\frac{1}{1+\beta+\frac{1}{\beta}}.$
(A.9)

Since $\xi^{2}$ is positive definite, $\xi$ is real and the flow is
perpendicular and supersonic.

If $M_{\parallel}=0$ but $M_{\perp}$ does not satisfy (A.9), equation (A.7)
reduces to the following bi-quadratic equation

$\displaystyle\xi^{4}+\frac{(1+\beta)}{M_{\perp}^{2}\beta}\left(\frac{\frac{1}{\beta+1}-M_{\perp}^{2}}{M_{\perp}^{2}-\frac{1+\beta}{\beta}}\right)\xi^{2}+\frac{1}{M_{\perp}^{2}\beta}\frac{1}{M_{\perp}^{2}-\frac{1+\beta}{\beta}}=0.$
(A.10)

For the roots to be real, the discriminant of (A.10) must be positive, which
implies

$\displaystyle
M_{\perp}^{2}>\frac{\frac{\beta+1}{\beta-1}+2\sqrt{\beta}}{\beta-1}\>\>\text{or}\>\>M_{\perp}^{2}<\frac{\frac{\beta+1}{\beta-1}-2\sqrt{\beta}}{\beta-1}.$
(A.11)

We find that $\beta=1$ represents a discontinuity in the solution as depicted
in Figure (1(a)) and (1(b)). For values of $\beta$ close to zero we find that
perpendicular flows can exist only if $M_{\perp}\approx 1$. For large $\beta$
we find that the first condition is easily satisfied but not the second one.

Furthermore, for $\xi$ to have four real roots we need

$\displaystyle
M_{\perp}^{2}>\left(1+\frac{1}{\beta}\right)\>\>\&\;\>M_{\perp}^{2}>\frac{1}{\beta+1}.$
(A.12)

On the other hand, when $M_{\perp}^{2}<1+1/\beta$, $\xi$ has at least two real
roots provided (A.11) is satisfied. Therefore, both supersonic and subsonic
perpendicular flows can exist.

Figure 2: Depiction of the dispersion relation (A.3) as a curve on a unit
sphere with coplanar unit vectors
$(\bm{\hat{B}},\bm{\hat{k}},\bm{\hat{u}_{\perp}})$. As the unit vector
$\hat{u}_{\perp}$ moves on the red circle, $\bm{\hat{k}}$ traces out the blue
circle. $\xi=\tan{\theta}$ satisfies (A.7).

Finally, we note that $\bm{\hat{k}},\hat{\bm{B}},\bm{\hat{u}_{\perp}}$ can be
plotted on a unit sphere as shown in Figure 2. We choose $\hat{\bm{B}}$ to be
in the $z$ direction and, therefore, $\bm{\hat{u}_{\perp}}$ must lie in the
$x-y$ plane. The coplanar vector $\bm{\hat{k}}$ makes an angle $\theta$ with
$\bm{B}/B$, where $\theta=\arctan{\xi}$. Since the equation for $\xi$, (A.7),
depends on $\bm{u}$ only through the combinations
$M_{\parallel}\propto\bm{u}\cdot\bm{B}$ and
$M_{\perp}\propto|\bm{u}_{\perp}|$, it is independent of the direction of
$\hat{\bm{u}}_{\perp}$. As a result, the curve that $\bm{k}$ traces out is
symmetric around $\hat{B}$ and hence a circle on the unit sphere.

##  Appendix B Useful identities and expressions

We collect some useful identities and expressions in this appendix that have
been used in deriving some of the key results in the text.

A useful expression is that of the commutator
$[\bm{B}\cdot\bm{\nabla},\bm{Q}\cdot\bm{\nabla}]$. Employing standard Einstein
summation convention, $\bm{B}\cdot\bm{\nabla}=B_{i}\partial_{i}$,
$\bm{Q}\cdot\bm{\nabla}=Q_{j}\partial_{j}$, together with (5.4b) implies that
$\epsilon_{klm}B_{l}Q_{m}=\partial_{k}\psi$. It is straight-forward to show
that

$\displaystyle[\bm{B}\cdot\bm{\nabla},\bm{Q}\cdot\bm{\nabla}]f=\bm{B}\cdot\bm{\nabla}(\bm{Q}\cdot\bm{\nabla}f)-\bm{Q}\cdot\bm{\nabla}(\bm{B}\cdot\bm{\nabla}f)=\bm{\nabla}\cdot\left(\bm{\nabla}\psi\times\bm{\nabla}f\right)=0$
(B.1)

The following identity is useful to simplify the divergence of
$\bm{u}\times\bm{\nabla}\psi$

$\displaystyle\bm{u}\times\bm{\nabla}\psi=\bm{u}\times(\bm{B}\times\bm{Q})=\left(\bm{u}\cdot\bm{Q}\right)\bm{B}-\left(\bm{u}\cdot\bm{B}\right)\bm{Q}$
(B.2a)
$\displaystyle\bm{\nabla}\cdot\left(\bm{u}\times\bm{\nabla}\psi\right)=\bm{B}\cdot\bm{\nabla}\left(\bm{u}\cdot\bm{Q}\right)-\bm{Q}\cdot\bm{\nabla}\left(\bm{u}\cdot\bm{B}\right).$
(B.2b)

##  Appendix C Closed field line consistency conditions and current
singularities

The close relation between current singularities that can develop on rational
surfaces in a nonsymmetric toroidal domain in ideal MHD, and the necessary
conditions in the form of integrals over closed magnetic field lines such as
$\oint d\ell/B=f(\psi)$ are well known (Hamada, 1962; Grad, 1967; Solov’ev &
Shafranov, 1970; Boozer, 1981; Loizu et al., 2015a). We show here that to
avoid current singularities on rational surface the integral constraints of
the form (5.23) for ideal MHD equilibrium, and (5.24), (5.25) for steady
flows, are sufficient conditions.

We shall use a generic $(\psi,\alpha,\varphi)$ coordinate system, where $\psi$
is the flux-label, $\alpha$ is the field-line label, $\varphi$ denotes another
angle like coordinate. We will assume that the magnetic field lines are closed
and shear is weak such that both $\bm{\nabla}\alpha$ and $\bm{\nabla}\varphi$
are single-valued. A more general derivation for finite shear can be carried
out using generalized Boozer coordinates.

In the $(\psi,\alpha,\varphi)$ coordinates, we can represent $\bm{B}$ in the
following covariant and contravariant forms

$\displaystyle\bm{B}=B_{\psi}\bm{\nabla}\psi+B_{\alpha}\bm{\nabla}\alpha+B_{\varphi}\bm{\nabla}\varphi,\quad\bm{B}=\bm{\nabla}\psi\times\bm{\nabla}\alpha.$
(C.1)

The functions $B_{\psi},B_{\alpha},B_{\varphi}$ are all single-valued
functions since $\bm{B}$ is single-valued and so are the gradients of
$(\psi,\alpha,\varphi)$. The Jacobian
$\mathcal{J}^{-1}=\bm{\nabla}\psi\times\bm{\nabla}\alpha\cdot\bm{\nabla}\varphi=\bm{B}\cdot\bm{\nabla}\varphi$
connects $B_{\varphi}$ to $B^{2}$ through

$\displaystyle B_{\varphi}=\mathcal{J}B^{2}.$ (C.2)

The current obtained by taking the curl of the covariant form of $\bm{B}$
satisfies

$\displaystyle\bm{J}=\bm{\nabla}B_{\psi}\times\bm{\nabla}\psi+\bm{\nabla}B_{\alpha}\times\bm{\nabla}\alpha+\bm{\nabla}B_{\varphi}\times\bm{\nabla}\varphi$
(C.3)
$\displaystyle\bm{J}\times\bm{B}=\bm{\nabla}\psi\left(\bm{B}\cdot\bm{\nabla}B_{\psi}-\mathcal{J}^{-1}\partial_{\psi}B_{\varphi}\right)+\bm{\nabla}\alpha\left(\bm{B}\cdot\bm{\nabla}B_{\alpha}-\mathcal{J}^{-1}\partial_{\alpha}B_{\varphi}\right)$
(C.4)

We now consider a general force-balance equation of the form

$\displaystyle\bm{J}\times\bm{B}=\bm{F},\quad\bm{F}=F_{\psi}\bm{\nabla}\psi+F_{\alpha}\bm{\nabla}\alpha.$
(C.5)

Note that $\bm{F}$ does not have a $\bm{\nabla}\varphi$ component because of
parallel force balance condition $\bm{F}\cdot\bm{B}=0$.

If the current $\bm{J}$ simultaneously satisfies
$\bm{J}=\bm{\nabla}\times\bm{B}$ and $\bm{J}\times\bm{B}=\bm{F}$, we can
equate $\bm{J}\times\bm{B}$ from (C.5) to (C.4) to get the following magnetic
differential equations for $B_{\psi}$ and $B_{\alpha}$ that enforce force-
balance

$\displaystyle\bm{B}\cdot\bm{\nabla}B_{\psi}-\mathcal{J}^{-1}\partial_{\psi}B_{\varphi}=F_{\psi},\quad\bm{B}\cdot\bm{\nabla}B_{\alpha}-\mathcal{J}^{-1}\partial_{\alpha}B_{\varphi}=F_{\alpha}.$
(C.6)

Eliminating $B_{\varphi}$ we obtain another magnetic differential equation

$\displaystyle\bm{B}\cdot\bm{\nabla}\left(\partial_{\alpha}B_{\psi}-\partial_{\psi}B_{\alpha}\right)=\mathcal{J}^{-1}\left(\partial_{\alpha}(\mathcal{J}F_{\psi})-\partial_{\psi}(\mathcal{J}F_{\alpha})\right)$
(C.7)

Since $B_{\psi}$ and $B_{\alpha}$ are single-valued, the well-known necessary
and sufficient conditions for the solvability of magnetic differential
equations (Newcomb, 1959) leads to the following constraints

$\displaystyle\partial_{\psi}\oint d\varphi B_{\phi}+\oint
d\phi\mathcal{J}F_{\psi}=0,\quad\partial_{\alpha}\oint d\varphi B_{\phi}+\oint
d\phi\mathcal{J}F_{\alpha}=0$ (C.8) $\displaystyle\partial_{\alpha}\oint
d\phi\mathcal{J}F_{\psi}-\partial_{\psi}\oint d\phi\mathcal{J}F_{\alpha}=0$

Using (C.2) and $d\ell/B=\mathcal{J}d\varphi$ we see that these constraints
can be written as

$\displaystyle\partial_{\psi}\oint\bm{B}\cdot\bm{d\ell}+\oint\frac{d\ell}{B}F_{\psi}=0,\quad\partial_{\alpha}\oint\bm{B}\cdot\bm{d\ell}+\oint\frac{d\ell}{B}F_{\alpha}=0.$
(C.9)
$\displaystyle\partial_{\alpha}\oint\frac{d\ell}{B}F_{\psi}-\partial_{\psi}\oint\frac{d\ell}{B}F_{\alpha}=0$

We note that the last consistency condition is not independent of the first
two since it is obtained by eliminating $B_{\phi}$ from the first two
conditions. However, once the integral constraints (C.9) are satisfied, we can
solve the MDEs for $B_{\psi}$ and $B_{\alpha}$. Therefore, the constraints
(C.9) are necessary and sufficient for force-balance.

For ideal MHD equilibrium $\bm{F}=p^{\prime}(\psi)\bm{\nabla}\psi$ and we
recover (5.23) from (C.9). For steady flows both $F_{\psi}$ and $F_{\alpha}$
are nonzero and we recover (5.24b) and (5.25).

To show the sufficiency of the integral conditions, we start with the solution
of the force-balance condition and impose the divergence-free condition on the
current. Writing $\bm{J}=\bm{J}_{\perp}+\bm{B}(j_{||}/B)$, and taking the
divergence, we obtain a magnetic differential equation for the parallel
component of current,

$\displaystyle\bm{B}\cdot\bm{\nabla}\left(\frac{j_{||}}{B}\right)+\bm{\nabla}\cdot\bm{J}_{\perp}=0.$
(C.10)

Therefore, if the Newcomb condition

$\displaystyle\oint\frac{d\ell}{B}\bm{\nabla}\cdot\bm{J}_{\perp}=0$ (C.11)

is not satisfied, then current singularity follows (Loizu et al., 2015a).

Using (C.5) and (C.1) we can obtain an expression for the perpendicular
component of the current

$\displaystyle\bm{J}_{\perp}=\frac{\bm{B}\times\bm{F}}{B^{2}}=\bm{B}\left(\frac{B_{\psi}F_{\alpha}-F_{\psi}B_{\alpha}}{B^{2}}\right)+\mathcal{J}^{-1}B_{\varphi}\left(F_{\psi}\bm{e}_{\alpha}-F_{\alpha}\bm{e}_{\psi}\right),$
(C.12)

where, $\bm{e}_{\psi}=\mathcal{J}\bm{\nabla}\alpha\times\bm{\nabla}\varphi$
etc are the basis vectors. The divergence of $\bm{J}_{\perp}$ is

$\displaystyle\bm{\nabla}\cdot\bm{J}_{\perp}=\bm{B}\cdot\bm{\nabla}\left(\frac{1}{B^{2}}\left(B_{\psi}F_{\alpha}-F_{\psi}B_{\alpha}\right)\right)+\mathcal{J}^{-1}\left(\partial_{\alpha}(\mathcal{J}F_{\psi})-\partial_{\psi}(\mathcal{J}F_{\alpha})\right).$
(C.13)

If the integral constraints (C.9) are satisfied, (C.7) is solvable and we can
rewrite (C.13) as

$\displaystyle\bm{\nabla}\cdot\bm{J}_{\perp}=\bm{B}\cdot\bm{\nabla}\left(\frac{1}{B^{2}}\left(B_{\psi}F_{\alpha}-F_{\psi}B_{\alpha}\right)+\left(\partial_{\alpha}B_{\psi}-\partial_{\psi}B_{\alpha}\right)\right),$
(C.14)

which satisfies (C.11). The MDE (C.10) can then be solved for the Pfirsch-
Schlüter currents.

Therefore, the integral constraints (C.9) are sufficient for the current
$\bm{J}$ to satisfy force-balance and be divergence-free.

##  Appendix D Derivation of the flow-modified GGS

Since the steps in deriving the flow-modified GGS are a bit tricky and one can
easily get lost in the plethora of terms, we give a step-by-step derivation in
this Appendix. Starting with the flow $\bm{u}$ given by (5.10) we obtain the
following expression for the vorticity $\bm{\omega}$ and
$\bm{u}\times\bm{\omega}$

$\displaystyle\bm{\omega}\equiv\bm{\nabla}\times\bm{u}=$
$\displaystyle\frac{\Lambda}{\rho}\bm{J}+\bm{\nabla}\left(\frac{\Lambda}{\rho}\right)\times\bm{B}-\Phi^{\prime\prime}(\psi)\bm{\nabla}\psi\times\bm{Q}-\Phi^{\prime}(\psi)\bm{\nabla}\times\bm{Q}.$
(D.1a) $\displaystyle\bm{u}\times\bm{\omega}=$
$\displaystyle-\frac{\Lambda}{\rho}\bm{J}\times\bm{u}-\bm{u}\cdot\bm{\nabla}\left(\frac{\Lambda}{\rho}\right)\bm{B}+\left(\bm{u}\cdot\bm{B}\right)\bm{\nabla}\left(\frac{\Lambda}{\rho}\right)$
(D.1b)
$\displaystyle-\Phi^{\prime\prime}(\psi)\left(\bm{u}\cdot\bm{Q}\right)\bm{\nabla}\psi-\Phi^{\prime}(\psi)\bm{u}\times(\bm{\nabla}\times\bm{Q}).$

With the help of (D.1) we obtain the right hand side of $\bm{\nabla}\psi$
dotted with the vorticity equation (5.11)

$\displaystyle\bm{\nabla}\psi\cdot\left(\bm{u}\times\bm{\omega}+\frac{1}{\rho}\bm{J}\times\bm{B}\right)=\bm{\nabla}\psi\cdot$
$\displaystyle\left(\frac{1}{\rho}\bm{J}\times\left(\bm{B}-\Lambda\bm{u}\right)+\left(\bm{u}\cdot\bm{B}\right)\bm{\nabla}\left(\frac{\Lambda}{\rho}\right)\right.$
(D.2)
$\displaystyle-\bigg{.}\Phi^{\prime\prime}(\psi)\left(\bm{u}\cdot\bm{Q}\right)\bm{\nabla}\psi\bigg{)}+\Phi^{\prime}(\psi)\left(\bm{u}\times\bm{\nabla}\psi\right)\cdot\left(\bm{\nabla}\times\bm{Q}\right)$

Now, using (5.14a), the left side of $\bm{\nabla}\psi$ dotted with the
vorticity equation yields

$\displaystyle\bm{\nabla}\psi\cdot\left(\frac{1}{2}u^{2}+h\right)=\bm{\nabla}\psi\cdot\left(\bm{\nabla}\psi\left(\frac{dH(\psi)}{d\psi}-\Phi^{\prime\prime}(\psi)\left(\bm{u}\cdot\bm{Q}\right)\right)+\bm{\nabla}\mathcal{H}-\Phi^{\prime}(\psi)\bm{\nabla}\left(\bm{u}\cdot\bm{Q}\right)\right)$
(D.3)

We can see that the $\Phi^{\prime\prime}(\psi)$ term cancels when we equate
the left (D.3) to the right side (D.2).

Let us first evaluate the first term on the right hand side of (D.2) by
replacing $\bm{\nabla}\psi$ by $\bm{B}\times\bm{Q}$ such that

$\displaystyle\bm{\nabla}\psi\cdot\bm{J}\times\left(\bm{B}-\Lambda\bm{u}\right)=(\bm{J}\cdot\bm{Q})(B^{2}-\Lambda\left(\bm{u}\cdot\bm{B}\right))-(\bm{J}\cdot\bm{Q})(\left(\bm{B}\cdot\bm{Q}\right)-\Lambda\left(\bm{u}\cdot\bm{Q}\right))$
(D.4)

The components of $\bm{J}$ can be calculated as follows

$\displaystyle\bm{J}\cdot\bm{Q}$
$\displaystyle=\bm{\nabla}\cdot\left(\bm{B}\times\bm{Q}\right)+\bm{\nabla}\cdot\left(\bm{Q}\times\bm{B}\right)=\nabla^{2}\psi+\frac{1}{Q^{2}}\left(\left(\bm{B}\cdot\bm{Q}\right)\bm{Q}+\bm{Q}\times\bm{\nabla}\psi\right)\cdot(\bm{\nabla}\times\bm{Q})$
$\displaystyle=Q^{2}\bm{\nabla}\cdot\left(\frac{1}{Q^{2}}\bm{\nabla}\psi\right)+\frac{1}{Q^{2}}\left(\bm{w}\cdot\bm{\nabla}\psi+\left(\bm{B}\cdot\bm{Q}\right)\bm{Q}\cdot\bm{\nabla}\times\bm{Q}\right)$
(D.5a) $\displaystyle\bm{J}\cdot\bm{B}$
$\displaystyle=\bm{J}\cdot\frac{1}{Q^{2}}\left(\left(\bm{B}\cdot\bm{Q}\right)\bm{Q}+\bm{Q}\times\bm{\nabla}\psi\right)=\frac{\left(\bm{B}\cdot\bm{Q}\right)}{Q^{2}}\bm{J}\cdot\bm{Q}+\left(\bm{F_{Q}}-\bm{\nabla}\left(\bm{B}\cdot\bm{Q}\right)\right)\cdot\frac{\bm{\nabla}\psi}{Q^{2}}$
(D.5b)

Substituting (D.5) in (D.4) and using

$\displaystyle
B^{2}=\frac{1}{Q^{2}}\left(\left(\bm{B}\cdot\bm{Q}\right)^{2}+|\bm{\nabla}\psi|^{2}\right),\quad\left(\left(\bm{u}\cdot\bm{B}\right)-\frac{\left(\bm{B}\cdot\bm{Q}\right)\left(\bm{u}\cdot\bm{Q}\right)}{Q^{2}}\right)=\frac{\Lambda}{\rho}|\bm{\nabla}\psi|^{2},$
(D.6)

we get

$\displaystyle\frac{\bm{\nabla}\psi}{\rho}\cdot\bm{J}\times\left(\bm{B}-\Lambda\bm{u}\right)=$
$\displaystyle(\bm{J}\cdot\bm{Q})\frac{|\bm{\nabla}\psi|^{2}}{\rho
Q^{2}}\left(1-\frac{\Lambda^{2}}{\rho}\right)+\frac{(\bm{B}-\Lambda\bm{u})\cdot\bm{Q}}{\rho
Q^{2}}\bm{\nabla}\psi\cdot\left(\bm{\nabla}\left(\bm{B}\cdot\bm{Q}\right)+\bm{F_{Q}}\right)$
$\displaystyle=$
$\displaystyle\frac{|\bm{\nabla}\psi|^{2}}{\rho}\bm{\nabla}\cdot\left(\frac{1}{Q^{2}}\left(1-\frac{\Lambda^{2}}{\rho}\right)\bm{\nabla}\psi\right)+Y_{1}$
(D.7) $\displaystyle+\frac{\bm{\nabla}\psi}{\rho
Q^{2}}\cdot\left(|\bm{\nabla}\psi|^{2}\bm{\nabla}\left(\frac{\Lambda^{2}}{\rho}\right)+\left(\left(\bm{B}\cdot\bm{Q}\right)-\Lambda\left(\bm{u}\cdot\bm{Q}\right)\right)\bm{\nabla}\left(\bm{B}\cdot\bm{Q}\right)\right)$

where,

$\displaystyle Y_{1}=\frac{1}{\rho
Q^{2}}\left(\frac{|\bm{\nabla}\psi|^{2}}{Q^{2}}\left(1-\frac{\Lambda^{2}}{\rho}\right)\left(\bm{w}\cdot\bm{\nabla}\psi+\left(\bm{B}\cdot\bm{Q}\right)\bm{Q}\cdot\bm{\nabla}\times\bm{Q}\right)+\left(\left(\bm{B}\cdot\bm{Q}\right)-\Lambda\left(\bm{u}\cdot\bm{Q}\right)\right)\bm{\nabla}\psi\cdot\bm{F_{Q}}\right)$
(D.8)

In simplifying the last term term of (D.2) we use instead of
$\bm{u}\times\bm{\nabla}\psi$

$\displaystyle\bm{u}\times(\bm{B}\times\bm{Q})=\left(\bm{u}\cdot\bm{Q}\right)\bm{B}-\left(\bm{u}\cdot\bm{B}\right)\bm{Q}=$
$\displaystyle\frac{1}{Q^{2}}\left(-|\bm{\nabla}\psi|^{2}\frac{\Lambda}{Q}\bm{Q}+\left(\bm{u}\cdot\bm{Q}\right)\bm{Q}\times\bm{\nabla}\psi\right)$
$\displaystyle\therefore\quad\Phi^{\prime}(\psi)\left(\bm{u}\times\bm{\nabla}\psi\right)\cdot\left(\bm{\nabla}\times\bm{Q}\right)=$
$\displaystyle-\Phi^{\prime}(\psi)\frac{\bm{\nabla}\psi\cdot\bm{\nabla}Q^{2}}{Q^{2}}\left(\bm{u}\cdot\bm{Q}\right)-Y_{2}$
(D.9) $\displaystyle Y_{2}=$
$\displaystyle-\Phi^{\prime}(\psi)\left(\frac{\left(\bm{u}\cdot\bm{Q}\right)}{Q^{2}}\bm{w}\cdot\bm{\nabla}\psi-\frac{\Lambda}{\rho}\frac{\bm{Q}\cdot\bm{\nabla}\times\bm{Q}}{Q^{2}}|\bm{\nabla}\psi|^{2}\right)$

We now rearrange all the terms in the $\bm{\nabla}\psi\cdot$ component of the
vorticity equation in the form

$\displaystyle|\bm{\nabla}\psi|^{2}\left(\rho\frac{dH(\psi)}{d\psi}+\bm{\nabla}\cdot\left(\frac{1}{Q^{2}}\left(1-\frac{\Lambda^{2}}{\rho}\right)\bm{\nabla}\psi\right)\right)+\rho\bm{\nabla}\psi\cdot\bm{\nabla}\mathcal{H}+\rho(Y_{1}+Y_{2}+Y_{3})=0,$
(D.10)

where $Y_{3}$ is given by

$\displaystyle Y_{3}=-\bm{\nabla}\psi\cdot\left[-\frac{1}{\rho
Q^{2}}\left(|\bm{\nabla}\psi|^{2}\bm{\nabla}\left(\frac{\Lambda^{2}}{\rho}\right)+\left(C-\Lambda\left(\bm{u}\cdot\bm{Q}\right)\right)\bm{\nabla}C\right)+\left(\bm{u}\cdot\bm{B}\right)\bm{\nabla}\left(\frac{\Lambda}{\rho}\right)\right.$
$\displaystyle\left.+\Phi^{\prime}(\psi)\left(\bm{\nabla}\left(\bm{u}\cdot\bm{Q}\right)-\frac{\bm{\nabla}Q^{2}}{Q^{2}}\left(\bm{u}\cdot\bm{Q}\right)\right)\right],$
(D.11)

with $C=\left(\bm{B}\cdot\bm{Q}\right)$. We simplify $Y_{3}$ by using
$\bm{\nabla}(\Lambda^{2}/\rho)=\Lambda\bm{\nabla}(\Lambda/\rho)+(\Lambda/\rho)\bm{\nabla}\Lambda$
and (D.6)

$\displaystyle\frac{1}{\rho}|\bm{\nabla}\psi|^{2}\bm{\nabla}\left(\frac{\Lambda^{2}}{\rho}\right)=\left(\left(\bm{u}\cdot\bm{B}\right)-\frac{\left(\bm{u}\cdot\bm{Q}\right)C}{Q^{2}}\right)\left(\bm{\nabla}\left(\frac{\Lambda}{\rho}\right)+\frac{\bm{\nabla}\Lambda}{\rho}\right).$
(D.12)

Substituting (D.12) into $Y_{3}$ and simplifying we get

$\displaystyle
Y_{3}=\frac{|\bm{\nabla}\psi|^{2}}{\rho}\left(\bm{u}\cdot\bm{B}\right)\frac{d\Lambda}{d\psi}+$
$\displaystyle\frac{1}{\rho}\bm{\nabla}\psi\cdot\left[-\frac{C}{\rho
Q^{2}}\bm{\nabla}C+\frac{\left(\bm{u}\cdot\bm{Q}\right)}{Q^{2}}\left(C\bm{\nabla}\left(\frac{\Lambda}{\rho}\right)+\frac{C}{\rho}\bm{\nabla}\Lambda+\frac{\Lambda}{\rho}\bm{\nabla}C\right)\right.$
.
$\displaystyle+\Phi^{\prime}(\psi)\left(\bm{\nabla}\left(\bm{u}\cdot\bm{Q}\right)-\frac{\bm{\nabla}Q^{2}}{Q^{2}}\left(\bm{u}\cdot\bm{Q}\right)\right)\bigg{]}$
(D.13)

The final step is to use the generalized momentum equation (5.14b) to replace
the term, $C\bm{\nabla}C$ in (D.13) by
$C\bm{\nabla}\left(I(\psi)+\mathcal{I}+\Lambda\left(\bm{u}\cdot\bm{Q}\right)\right)$.
After some straightforward algebra we get

$\displaystyle\rho
Y_{3}=|\bm{\nabla}\psi|^{2}\left(\left(\bm{u}\cdot\bm{B}\right)\frac{d\Lambda}{d\psi}+\frac{\left(\bm{B}\cdot\bm{Q}\right)}{Q^{2}}\frac{dI}{d\psi}-\rho\frac{d\Phi^{\prime}(\psi)}{d\psi}\left(\bm{u}\cdot\bm{Q}\right)\right)+\left(\bm{B}\cdot\bm{Q}\right)\bm{\nabla}\psi\cdot\bm{\nabla}\mathcal{I}$
(D.14)

Substituting (D.14) into (D.10) we get (5.20).

##  Appendix E Derivation of the final Hamada condition

To derive the final Hamada condition (5.25), it is convenient to set up flux-
coordinates $(\psi,\alpha,\beta)$ (discussed in (Hameiri, 1998)) such that

$\displaystyle\bm{B}=\bm{\nabla}\psi\times\bm{\nabla}\alpha,\quad\bm{Q}=\bm{\nabla}\psi\times\bm{\nabla}\beta.$
(E.1)

For closed field line systems with zero magnetic shear, $\bm{\nabla}\alpha$ is
a single-valued quantity. Since, $\bm{B}$ and $\bm{Q}$ must be single-valued,
$\bm{\nabla}\beta$ in general is either single-valued or at least the multi-
valued part of $\bm{\nabla}\beta$ is in the $\bm{\nabla}\psi$ direction
(Hameiri, 1998). For closed $\bm{Q}$ lines, $\beta$ is single-valued.

Although we use the $(\psi,\alpha,\beta)$ coordinates to derive (5.25), the
result is independent of the choice and details of the flux coordinates. In
particular, any flux-coordinates that allow radial currents, e.g., generalized
Boozer coordinates (Rodriguez et al., 2020) can be used.

The condition $\bm{B}\times\bm{Q}=\bm{\nabla}\psi$ implies that the Jacobian
$\mathcal{J}^{-1}=\bm{\nabla}\psi\times\bm{\nabla}\alpha\cdot\bm{\nabla}\beta$
is unity. Therefore,

$\displaystyle\bm{B}\cdot\bm{\nabla}=\partial_{\beta},\quad\bm{Q}\cdot\bm{\nabla}=-\partial_{\alpha}.$
(E.2)

In these coordinates, the divergence of a vector of the form
$f\bm{\nabla}\psi$ is given by

$\displaystyle\bm{\nabla}\cdot\left(f\bm{\nabla}\psi\right)=\partial_{\psi}\left(f|\bm{\nabla}\psi|^{2}\right)-\bm{Q}\cdot\bm{\nabla}\left(f\bm{\nabla}\psi\cdot\bm{\nabla}\alpha\right)+\bm{B}\cdot\bm{\nabla}\left(f\bm{\nabla}\psi\cdot\bm{\nabla}\beta\right).$
(E.3)

Averaging (E.3) along the closed magnetic field line we get

$\displaystyle\langle\bm{\nabla}\cdot\left(f\bm{\nabla}\psi\right)\rangle=\partial_{\psi}\langle
f|\bm{\nabla}\psi|^{2}\rangle-\bm{Q}\cdot\bm{\nabla}\langle
f\bm{\nabla}\psi\cdot\bm{\nabla}\alpha\rangle+f|\bm{\nabla}\psi|^{2}[\beta],$
(E.4)

where, $[\beta]$ denotes the jump in $\beta$ along the closed field line. If
$\bm{Q}$ lines are closed $[\beta]=0$.

We obtain (5.25) for $f=(1-\Lambda^{2}/\rho)/Q^{2}$ upon using (5.14),(5.24)
and the identity

$\displaystyle\left(1-\frac{\Lambda^{2}}{\rho}\right)\frac{1}{Q^{2}}|\bm{\nabla}\psi|^{2}=B^{2}-\Lambda\left(\bm{u}\cdot\bm{B}\right)-\frac{\left(\bm{B}\cdot\bm{Q}\right)}{Q^{2}}(I+\mathcal{I}).$
(E.5)

## References

  * Abel et al. (2013) Abel, IG, Plunk, GG, Wang, E, Barnes, M, Cowley, SC, Dorland, W & Schekochihin, AA 2013 Multiscale gyrokinetics for rotating tokamak plasmas: fluctuations, transport and energy flows. Reports on Progress in Physics 76 (11), 116201.
  * Andreussi et al. (2012) Andreussi, Tommaso, Morrison, Philip J & Pegoraro, Francesco 2012 Hamiltonian magnetohydrodynamics: Helically symmetric formulation, casimir invariants, and equilibrium variational principles. Physics of Plasmas 19 (5), 052102.
  * Beskin (1997) Beskin, Vasilii S 1997 Axisymmetric stationary flows in compact astrophysical objects. Physics-Uspekhi 40 (7), 659\.
  * Beskin (2009) Beskin, Vasily S 2009 MHD Flows in Compact Astrophysical Objects: Accretion, Winds and Jets. Springer Science & Business Media.
  * Bogoyavlenskij (2000a) Bogoyavlenskij, Oleg I 2000a Astrophysical jets as exact plasma equilibria. Physical review letters 84 (9), 1914.
  * Bogoyavlenskij (2000b) Bogoyavlenskij, Oleg I 2000b Helically symmetric astrophysical jets. Physical Review E 62 (6), 8616.
  * Boozer (1981) Boozer, Allen H 1981 Plasma equilibrium with rational magnetic surfaces. The Physics of Fluids 24 (11), 1999–2003.
  * Burby et al. (2020a) Burby, JW, Kallinikos, N & MacKay, RS 2020a Some mathematics for quasi-symmetry. Journal of Mathematical Physics 61 (9), 093503.
  * Burby et al. (2020b) Burby, Joshua William, Kallinikos, Nikos & MacKay, RS 2020b Generalized grad–shafranov equation for non-axisymmetric mhd equilibria. Physics of Plasmas 27 (10), 102504\.
  * Burrell (1997) Burrell, KH 1997 Effects of E$\times$B velocity shear and magnetic shear on turbulence and transport in magnetic confinement devices. Physics of Plasmas 4 (5), 1499–1518.
  * Burrell (1999) Burrell, Keith H 1999 Tests of causality: Experimental evidence that sheared e$\times$ b flow alters turbulence and transport in tokamaks. Physics of Plasmas 6 (12), 4418–4435.
  * Canik et al. (2007) Canik, JM, Anderson, DT, Anderson, FSB, Likin, KM, Talmadge, JN & Zhai, K 2007 Experimental demonstration of improved neoclassical transport with quasihelical symmetry. Physical review letters 98 (8), 085002.
  * Chu et al. (1995) Chu, MS, Greene, JM, Jensen, TH, Miller, RL, Bondeson, A, Johnson, RW & Mauel, ME 1995 Effect of toroidal plasma flow and flow shear on global magnetohydrodynamic mhd modes. Physics of Plasmas 2 (6), 2236–2241.
  * Constantin et al. (2020) Constantin, Peter, Drivas, Theodore D & Ginsberg, Daniel 2020 On quasisymmetric plasma equilibria sustained by small force. arXiv preprint arXiv:2009.08860 .
  * Courant & Friedrichs (1999) Courant, Richard & Friedrichs, Kurt Otto 1999 Supersonic flow and shock waves, , vol. 21. Springer Science & Business Media.
  * Davis & Tchekhovskoy (2020) Davis, Shane W & Tchekhovskoy, Alexander 2020 Magnetohydrodynamics simulations of active galactic nucleus disks and jets. Annual Review of Astronomy and Astrophysics 58, 407–439.
  * D’haeseleer et al. (1991) D’haeseleer, William D, Hitchon, William NG, Callen, James D & Shohet, J Leon 1991 Flux coordinates and magnetic field structure: a guide to a fundamental tool of plasma theory. Springer-Verlag Berlin Heidelberg.
  * Freidberg (2014) Freidberg, Jeffrey P 2014 ideal MHD. Cambridge University Press.
  * Fujisawa (2008) Fujisawa, Akihide 2008 A review of zonal flow experiments. Nuclear Fusion 49 (1), 013001.
  * Garren & Boozer (1991) Garren, DA & Boozer, Allen H 1991 Existence of quasihelically symmetric stellarators. Physics of Fluids B: Plasma Physics 3 (10), 2822–2834.
  * Goedbloed & Lifschitz (1997) Goedbloed, JiP & Lifschitz, A 1997 Stationary symmetric magnetohydrodynamic flows. Physics of Plasmas 4 (10), 3544–3564.
  * Goedbloed et al. (2004) Goedbloed, JP Hans, Goedbloed, JP & Poedts, Stefaan 2004 Principles of magnetohydrodynamics: with applications to laboratory and astrophysical plasmas. Cambridge university press.
  * Grad (1967) Grad, Harold 1967 Toroidal containment of a plasma. The Physics of Fluids 10 (1), 137–154.
  * Grad (1971) Grad, Harold 1971 Plasma containment in closed line systems. In Plasma Physics and Controlled Nuclear Fusion Research 1971. Vol. III. Proceedings of the Fourth International Conference on Plasma Physics and Controlled Nuclear Fusion Research.
  * Greene & Karlson (1969) Greene, John M & Karlson, Erik T 1969 Variational principle for stationary magnetohydrodynamic equilibria. The Physics of Fluids 12 (3), 561–567.
  * Groebner et al. (1990) Groebner, RJ, Burrell, KH & Seraydarian, RP 1990 Role of edge electric field and poloidal rotation in the L-H transition. Physical review letters 64 (25), 3015.
  * Guazzotto & Betti (2005) Guazzotto, L & Betti, R 2005 Magnetohydrodynamics equilibria with toroidal and poloidal flow. Physics of plasmas 12 (5), 056107.
  * Hamada (1962) Hamada, Shigeo 1962 Hydromagnetic equilibria and their proper coordinates. Nuclear Fusion 2 (1-2), 23.
  * Hameiri (1981) Hameiri, Eliezer 1981 Spectral estimates, stability conditions, and the rotating screw-pinch. Journal of Mathematical Physics 22 (9), 2080–2088.
  * Hameiri (1983) Hameiri, Eliezer 1983 The equilibrium and stability of rotating plasmas. The Physics of Fluids 26 (1), 230–237.
  * Hameiri (1998) Hameiri, Eliezer 1998 Variational principles for equilibrium states with plasma flow. Physics of Plasmas 5 (9), 3270–3281.
  * Hassam (1996) Hassam, AB 1996 Poloidal rotation of tokamak plasmas at super poloidal sonic speeds. Nuclear fusion 36 (6), 707.
  * Hegna (2012) Hegna, CC 2012 Plasma flow healing of magnetic islands in stellarators. Physics of Plasmas 19 (5), 056101.
  * Helander (1998) Helander, Per 1998 Neoclassical transport in a rotating impure plasma. Physics of Plasmas 5 (4), 1209–1211.
  * Helander (2007) Helander, P 2007 On rapid plasma rotation. Physics of Plasmas 14 (10), 104501.
  * Helander (2014) Helander, Per 2014 Theory of plasma confinement in non-axisymmetric magnetic fields. Reports on Progress in Physics 77 (8), 087001.
  * Helander et al. (2012) Helander, Per, Beidler, CD, Bird, TM, Drevlak, M, Feng, Y, Hatzky, R, Jenko, F, Kleiber, R, Proll, JHE, Turkin, Yu & others 2012 Stellarator and tokamak plasmas: a comparison. Plasma Physics and Controlled Fusion 54 (12), 124009.
  * Helander & Simakov (2008) Helander, P & Simakov, AN 2008 Intrinsic ambipolarity and rotation in stellarators. Physical review letters 101 (14), 145003.
  * Hinton & Wong (1985) Hinton, FL & Wong, SK 1985 Neoclassical ion transport in rotating axisymmetric plasmas. The Physics of fluids 28 (10), 3082–3098.
  * Igumenshchev (2008) Igumenshchev, Igor V 2008 Magnetically arrested disks and the origin of poynting jets: a numerical study. The Astrophysical Journal 677 (1), 317.
  * Igumenshchev et al. (2003) Igumenshchev, Igor V, Narayan, Ramesh & Abramowicz, Marek A 2003 Three-dimensional magnetohydrodynamic simulations of radiatively inefficient accretion flows. The Astrophysical Journal 592 (2), 1042.
  * Ilgisonis & Pastukhov (2000) Ilgisonis, VI & Pastukhov, VP 2000 Variational approaches to the problems of plasma stability and of nonlinear plasma dynamics. Journal of Experimental and Theoretical Physics Letters 72 (10), 530–540.
  * Istomin & Pariev (1996) Istomin, Ya N & Pariev, VI 1996 Stability of a relativistic rotating electron-positron jet: non-axisymmetric perturbations. Monthly Notices of the Royal Astronomical Society 281 (1), 1–26.
  * Kamchatnov (1982) Kamchatnov, AM 1982 Topological solitons in magnetohydrodynamics. Zh. Eksp. Teor. Fiz 82, 117–124.
  * Kovrizhnykh & Shchepetov (1980) Kovrizhnykh, LM & Shchepetov, SV 1980 Description of a stellarator plasma with averaged mhd equations. Sov. J. Plasma Phys.(Engl. Transl.);(United States) 6 (5).
  * Kovrizhnykh & Shchepetov (1989) Kovrizhnykh, LM & Shchepetov, SV 1989 Stellarator equilibria with steady flow. Nuclear fusion 29 (4), 667.
  * Landreman & Sengupta (2019) Landreman, Matt & Sengupta, Wrick 2019 Constructing stellarators with quasisymmetry to high order. Journal of Plasma Physics 85 (6).
  * Lifschitz & Goedbloed (1997) Lifschitz, A & Goedbloed, JP 1997 Transonic magnetohydrodynamic flows. Journal of plasma physics 58 (1), 61–99.
  * Lin et al. (1998) Lin, Zhihong, Hahm, Taik Soo, Lee, WW, Tang, William M & White, Roscoe B 1998 Turbulent transport reduction by zonal flows: Massively parallel simulations. Science 281 (5384), 1835–1837.
  * Loizu et al. (2015a) Loizu, J, Hudson, S, Bhattacharjee, A & Helander, P 2015a Magnetic islands and singular currents at rational surfaces in three-dimensional magnetohydrodynamic equilibria. Physics of Plasmas 22 (2), 022501.
  * Loizu et al. (2015b) Loizu, J, Hudson, SR, Bhattacharjee, A, Lazerson, Sam & Helander, P 2015b Existence of three-dimensional ideal-magnetohydrodynamic equilibria with current sheets. Physics of Plasmas 22 (9), 090704.
  * McClements & Hole (2010) McClements, KG & Hole, MJ 2010 On steady poloidal and toroidal flows in tokamak plasmas. Physics of Plasmas 17 (8), 082509.
  * McKinney & Blandford (2009) McKinney, Jonathan C & Blandford, Roger D 2009 Stability of relativistic jets from rotating, accreting black holes via fully three-dimensional magnetohydrodynamic simulations. Monthly Notices of the Royal Astronomical Society: Letters 394 (1), L126–L130.
  * Morawetz (1985) Morawetz, Cathleen S 1985 On a weak solution for a transonic flow problem. Communications on pure and applied mathematics 38 (6), 797–817.
  * Morozov & Solov’ev (1980) Morozov, AI & Solov’ev, LS 1980 Steady-state plasma flow in a magnetic field. Reviews of Plasma Physics/Voprosy Teorii Plazmy pp. 1–103.
  * Newcomb (1959) Newcomb, William A 1959 Magnetic differential equations. The Physics of Fluids 2 (4), 362–365.
  * Rees (1922) Rees, EL 1922 Graphical discussion of the roots of a quartic equation. The American Mathematical Monthly 29 (2), 51–55.
  * Rodriguez et al. (2020) Rodriguez, Eduardo, Helander, Per & Bhattacharjee, Amitava 2020 Necessary and sufficient conditions for quasisymmetry. Physics of Plasmas 27 (6), 062501.
  * Rodriguez et al. (2021) Rodriguez, Eduardo, Sengupta, Wrick & Bhattacharjee, Amitava 2021 Generalized boozer coordinates (in preparation). Journal of Plasma Physics .
  * Simakov & Helander (2009) Simakov, Andrei N & Helander, Per 2009 Neoclassical momentum transport in a collisional stellarator and a rippled tokamak. Physics of Plasmas 16 (4), 042503.
  * Simakov & Helander (2011) Simakov, Andrei N & Helander, Per 2011 Plasma rotation in a quasi-symmetric stellarator. Plasma Physics and Controlled Fusion 53 (2), 024005.
  * Solov’ev & Shafranov (1970) Solov’ev, LS & Shafranov, VD 1970 Plasma confinement in closed magnetic systems. Reviews of Plasma Physics pp. 1–247.
  * Spong (2005) Spong, DA 2005 Generation and damping of neoclassical plasma flows in stellarators. Physics of plasmas 12 (5), 056114\.
  * Stroth et al. (2011) Stroth, Ulrich, Manz, P & Ramisch, M 2011 On the interaction of turbulence and flows in toroidal plasmas. Plasma Physics and Controlled Fusion 53 (2), 024006.
  * Sugama et al. (2011) Sugama, H, Watanabe, T-H, Nunami, M & Nishimura, S 2011 Quasisymmetric toroidal plasmas with large mean flows. Physics of Plasmas 18 (8), 082505.
  * Tasso & Throumoulopoulos (1998) Tasso, Henri & Throumoulopoulos, George N 1998 Axisymmetric ideal magnetohydrodynamic equilibria with incompressible flows. Physics of plasmas 5 (6), 2378–2383.
  * Terry (2000) Terry, PW 2000 Suppression of turbulence and transport by sheared flow. Reviews of Modern Physics 72 (1), 109.
  * Tessarotto et al. (1996) Tessarotto, M., Johnson, J. L., White, R. B. & Zheng, L.‐J. 1996 Quasi‐helical magnetohydrodynamic equilibria in the presence of flow. Physics of Plasmas 3 (7), 2653–2663.
  * Throumoulopoulos et al. (2006) Throumoulopoulos, GN, Weitzner, Harold & Tasso, Henri 2006 On nonexistence of tokamak equilibria with purely poloidal flow. Physics of plasmas 13 (12), 122501.
  * Villata & Ferrari (1994) Villata, M & Ferrari, A 1994 Exact solutions for helical mhd equilibria of astrophysical jets. Astronomy and Astrophysics 284, 663–678.
  * Villata & Tsinganos (1993) Villata, Massimo & Tsinganos, Kanaris 1993 Exact solutions for helical magnetohydrodynamic equilibria. Physics of Fluids B: Plasma Physics 5 (7), 2153–2164.
  * Vladimirov & Moffatt (1995) Vladimirov, Vladimir A & Moffatt, HK 1995 On general transformations and variational principles for the magnetohydrodynamics of ideal fluids. part 1. fundamental principles. Journal of fluid mechanics 283, 125–139.
  * Weitzner (2014) Weitzner, Harold 2014 Ideal magnetohydrodynamic equilibrium in a non-symmetric topological torus. Physics of Plasmas 21 (2), 022515.
  * Weitzner (2016) Weitzner, Harold 2016 Expansions of non-symmetric toroidal magnetohydrodynamic equilibria. Physics of Plasmas 23 (6), 062512\.
  * White et al. (2019) White, Christopher J, Stone, James M & Quataert, Eliot 2019 A resolution study of magnetically arrested disks. The Astrophysical Journal 874 (2), 168.