3.2.1. Continuity of electrostatic potential

One might think that an analogous ‘stepped potential equilibrium’ could provide a solution to the problem of finding a non-trivial but tractable solution of $\boldsymbol {\nabla }_{\parallel }\varPhi \stackrel {\textrm {def}}{=}(\boldsymbol {B}/B)\boldsymbol {\cdot }\boldsymbol {\nabla }\varPhi =0$ in a chaotic magnetic-field-line flow. Unfortunately, however, we must restrict $\boldsymbol {\nabla }_{\perp }\varPhi \stackrel {\textrm {def}}{=}({\boldsymbol{\mathsf{I}}}-\boldsymbol {B}\boldsymbol {B}/B^{2})\boldsymbol {\cdot }\boldsymbol {\nabla }\varPhi$ to square-integrable functions in order to keep the $\boldsymbol {E}{\times }\boldsymbol {B}$ drift (3.3) from acquiring a $\delta$-function component.

This rules out having steps in $\varPhi$ because $\delta$-functions are not square integrable, so stepped potentials would make the kinetic energy integral infinite. However, this does not necessarily imply $\varPhi$ is constant in weakly chaotic regions with a finite measure of Kolmogorov–Arnold–Moser (KAM) surfaces – perhaps weak KAM theory (see e.g. Fathi Reference Fathi2009) would allow fractal potential profiles having finite kinetic energy associated with them.

As a way to handle non-constant $\varPhi$ computationally we propose using a penalty or augmented Lagrangian method (see e.g. Nocedal & Wright Reference Nocedal and Wright2006). That is, we treat Hamilton's principle as a constrained saddle-point optimization problem and add a penalty functional to the Hamiltonian, which regularizes the variational problem by approaching the (perhaps fractal) IMHD ‘feasible region’ of configuration space from outside, in the less-constrained space on which RxMHD is defined (which is smoother, see figure 1 of Dewar et al. (Reference Dewar, Burby, Qu, Sato and Hole2020)).

Another approach might be a time-evolution code with added dissipation such that the long-time solution is attracted to one having chaotic regions of constant pressure interspersed with integrable regions with changing pressure. This can be viewed as a steepest-descents solution of the same optimization problem.

4.1.1. The IOL as a softly constrained optimization problem

In implementing the parallel IOL constraint we propose to adapt the augmented Lagrangian method from finite-dimensional optimization theory, as described by Nocedal & Wright (Reference Nocedal and Wright2006, § 17.3), or for Banach spaces (see e.g. Kanzow et al. (Reference Kanzow, Steck and Wachsmuth2018) and references therein). This is a hybrid numerical method that combines two constraint approaches: the Lagrange multiplier method and the penalty function method sketched above.

As well as adapting notation and methods from optimization theory we have borrowed the term ‘feasible region’, meaning the range in which the vector $X$ of variables to be solved for is such that a set of equality constraints $c_{i}[X]=0$ are satisfied (also inequality constraints $c_{i}[X]>0$, but we do not consider this case). The infeasible region, is its complement, where one or more constraints are violated. By hard constraint we mean one where $X$ must be in the feasible region, and by soft constraint we mean one where $X$ need only be in some neighbourhood of the feasible region, which is useful both practically and for regularizing when, as in MHD, defining the boundary between feasible and infeasible is complicated by the possibility of singular behaviour like current sheets and reconnection points.

We now formulate two related physical tasks, the simpler one being (task 1) the equilibrium problem. In toroidal plasma confinement theory the most physically desirable states are stable, time-independent equilibria, i.e. minima of a Hamiltonian functional $\mathscr {H}[X]$, kinetic plus potential energy within a static boundary $\partial \varOmega$. We seek a numerical algorithm that starts from an initial guess for the physical fields $X$ and iterates to extremize (minimize if seeking a stable equilibrium) a Hamiltonian, under the IOL equality constraint, (1.2). Finding a stable equilibrium can be stated concisely as the optimization problem.

(4.2)\begin{align} Equilibrium:\quad& \min_{X}\mathscr{H}[X]\quad\text{subject to } \boldsymbol{C}[X](\boldsymbol{x})]=0,\notag\\ & \forall\,\boldsymbol{x}\in\varOmega \text{ and boundary conditions}\quad \forall\,\boldsymbol{x}\in\partial\varOmega,\end{align}

where

(4.3)\begin{equation} \boldsymbol{C}\stackrel{\text{def}}{=}\boldsymbol{E}+\boldsymbol{v}{\times}\boldsymbol{B},\end{equation}

is the generalization of the finite set of equality constraint functions $\{ c_{i}\}$ of Nocedal & Wright (Reference Nocedal and Wright2006) (as a 3-vector it is finite-dimensional but as a function of $\boldsymbol {x}$ it is infinite-dimensional).

For the purposes of the present paper $\mathscr {H}$ is the non-canonical version, $H_{\textrm {nc}}^{\textrm {MHD}}$, of the Hamiltonian, $H^{\textrm {MHD}}$ defined in (2.1) plus the global constraint terms described in the next subsection, § 4.2. The ideal boundary conditions are $\boldsymbol {E}+\boldsymbol {v}{\times }\boldsymbol {B}=0$, $\boldsymbol {{n}}\boldsymbol {\cdot }\boldsymbol {v}=0$, $\boldsymbol {{n}}\boldsymbol {\cdot }\boldsymbol {B}=0$ on $\partial \varOmega$ and $\varPhi ={\mathrm {const}.}$ on each disjoint component of $\partial \varOmega$ (think plates of a capacitor or electrodes of a vacuum tube).

To treat the implementation of the IOL in Hamilton's principle, a constrained saddle point optimization problem, we shall use the set of values of the components of $\boldsymbol {E}'(\boldsymbol {{x}}\in \varOmega )$; note the identities

(4.4a–c)\begin{equation} \frac{\partial\boldsymbol{C}}{\partial\boldsymbol{E}}={\boldsymbol{\mathsf{I}}}, \quad\frac{\partial\boldsymbol{C}}{\partial\boldsymbol{v}}={\boldsymbol{\mathsf{I}}}{\times}\boldsymbol{B} \quad\text{and}\quad\frac{\partial\boldsymbol{C}}{\partial\boldsymbol{B}} ={-}{\boldsymbol{\mathsf{I}}}{\times}\boldsymbol{v}.\end{equation}

We seek a soft form of the equilibrium constraint, i.e. a formulation such that $\boldsymbol {C}\to 0, \forall \,\boldsymbol {x}\in \varOmega$, where $\to$ denotes a limiting process whereby $X$ moves from the infeasible class of states where $\boldsymbol {C}\neq 0$ towards the feasible class defined pointwise as $\boldsymbol {C}=0,\,\forall \,\boldsymbol {x}\in \varOmega (t)$, or, in a weak form, as $||\boldsymbol {C}||=0$. Such a soft constraint procedure is provided by the augmented Lagrangian (or, rather, Hamiltonian in the equilibrium problem) as prescribed by Nocedal & Wright (Reference Nocedal and Wright2006, § 17.3),

(4.5)\begin{equation} \mathscr{H}_{{A}}\stackrel{\text{def}}{=}\int_{\varOmega}\left[\mathcal{H}- \boldsymbol{\lambda}\boldsymbol{\cdot}\boldsymbol{C}+\tfrac{1}{2}\mu^{P}_{\varOmega}\boldsymbol{C}^{2}\right]\,\mathrm{d} V,\end{equation}

where $\boldsymbol {\lambda }$ is a Lagrange multiplier and the spatial constant $\mu ^{P}_{\varOmega }\geq 0$ is a penalty multiplier of the non-negative quadratic penalty ${\frac {1}{2}}\int _{\varOmega }\boldsymbol {C}^{2}\,\mathrm {d} V$.

Nocedal & Wright (Reference Nocedal and Wright2006, § 17.3) give an iterative algorithmic framework that combines the advantages of both the Lagrange multiplier and penalty function methods. In their algorithm the user provides an increasing sequence $\{ \mu ^{P}_{\varOmega }|^{n},:n=0,1,2,\ldots \}$ of penalty multipliers and adjusts $\boldsymbol {\lambda }$ to solve for the minima of Hamiltonians with the Lagrange multiplier and penalty terms. The iteration update rule for the sequence of Lagrange multipliers and corresponding constraint residuals $\{ \boldsymbol {\lambda }^{n},\boldsymbol {C}^{n}:n=0,1,2,\ldots \}$ is

(4.6)\begin{equation} {{\boldsymbol{\lambda}^{n+1}}} =\boldsymbol{\lambda}^{n}-\mu^{P}_{\varOmega}|^{n}\boldsymbol{C}^{n}.\end{equation}

In the following sections we shall take the iteration index $n$ as implicit unless needed for clarity, with the updated $\boldsymbol {\lambda }$ as given by the right-hand side of (4.6) denoted by

(4.7)\begin{equation} \boldsymbol{\lambda}_{*}\stackrel{\text{def}}{=}\frac{\partial}{\partial\boldsymbol{C}}\left(\boldsymbol{\lambda}\boldsymbol{\cdot}\boldsymbol{C}- \frac{\mu_{\varOmega}^{{P}}}{2}\boldsymbol{C}^{2}\right)=\boldsymbol{\lambda}-\mu_{\varOmega}^{{P}} \boldsymbol{C},\end{equation}

which is a ‘best estimate’ of the optimum Lagrange multiplier given the current estimate and penalty multiplier.

The more difficult second physical task (task 2) is the time evolution problem. This is similar to task 1 except we seek an evolution, a dynamical path in space-time $\varOmega _{t}\times [t_{1},t_{2}]$ given a time-dependent boundary $\partial \varOmega _{t}$, the objective function for extremization now being an action integral. The stable minima of the Hamiltonian now becoming saddle points of the corresponding action functional $\mathscr {S}[X]$, kinetic minus potential energy. This task can be summarized as the pseudo optimization problem.

(4.8)\begin{equation} Dynamics:\quad\underset{X}{\text{extr}}\,\mathscr{S}[X] \quad\text{subject to }\boldsymbol{C}[X](\boldsymbol{x},t)=0, \forall\, \boldsymbol{x}\in\varOmega_{t},\,t\in[t_{1},t_{2}],\end{equation}

under the same boundary conditions as for equilibrium at each time $t$.

Although sometimes called the ‘principle of least action’, Hamilton's principle is often not an optimization problem but rather a saddle point problem, where the stationary point of $\mathscr {S}[X]$ cannot be found by a descent algorithm. (This is well known in nonlinear Hamiltonian dynamics (Meiss Reference Meiss1992) where periodic orbits are classified as (action) minimizing orbits, which are hyperbolic (unstable), or as minimax orbits, which are elliptic (stable).) Although ‘extremum’ or ‘extremization’ is not quite correct either, as extremum strictly means ‘maximum or minimum’, it is convenient to use the abbreviation ‘extr’ as an abbreviation for these words and add the rider ‘depending on direction of traversal’ (implying also the existence of neutral directions between max and min), so as to include saddle points.

To find saddle points requires some form of Newton method, needing at least estimates of the second variation (Hessian matrix) rather than a descent method. The augmented Lagrangian method still works if we solve (4.8) at each iteration, so here again we adopt it to solve for a stationary point of the augmented phase-space action functional (2.9),

(4.9)\begin{equation} \mathscr{S}_{{\rm ph}}^{{\rm A}}\stackrel{\text{def}}{=}\iint_{\varOmega} (\rho\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{v}-\mathcal{H}+\mathcal{L}_\varOmega^{C})\mathrm{d} V\mathrm{d} t,\end{equation}

where the augmented penalty constraint density $\mathcal {L}_\varOmega ^{C}$ is defined by

(4.10)\begin{equation} \mathcal{L}_\varOmega^{C}\stackrel{\text{def}}{=}\boldsymbol{\lambda}\boldsymbol{\cdot}\boldsymbol{C}-\frac{\mu^{P}_{\varOmega}\,\boldsymbol{C}^{2}}{2},\end{equation}

where $\boldsymbol {\lambda }$ and $\mu ^{P}_{\varOmega }$ are taken as external parameters in the application of Hamilton's principle at each iteration, giving a sequence of regularized magnetofluid models. When $\mu ^{P}_{\varOmega }=0$, the pure Lagrange multiplier method, feasible critical points of $\mathscr {H}_{{A}}$ might be saddle points with descending directions in the infeasible sector even if they are physically stable ideal equilibria where the IMHD Hamiltonian is minimized. When $\boldsymbol {\lambda }=0$, the pure penalty function method, feasible stable equilibria could be approximated arbitrarily well in the limit as $\mu ^{P}_{\varOmega }$ tends to infinity, but this becomes an increasingly ill-posed optimization problem. (It does, however, have the attractive feature of providing a continuous family of RxMHD models running from the RxMHD of Dewar et al. (Reference Dewar, Burby, Qu, Sato and Hole2020) when $\mu ^{P}_{\varOmega }=0$ to a subset of weak IMHD when $\mu ^{P}_{\varOmega }\to +\infty$.)

Remarks:

(i)  task 1 can be treated as a subclass of task 2 in which time derivatives are set to zero and $t$ is taken to be an irrelevant constant, but the Hamiltonian is more appropriate than the Lagrangian for treating it as an optimization problem;

(ii) the iteration method for implementing constraints is implicit, meaning that the state variables in the $n$th iteration need to be found by solving EL equations, taking it for granted the EL equations can be solved and any subiterations required have converged.

We shall not discuss detailed implementation issues here, except to remark that time evolution over a large time interval can be implemented numerically in an outer time stepping loop in which a large time interval is split into multiple short time intervals (time steps) $[t_{i},t_{i+1}]$, within each of which constraint iterations are repeated until converged to the required accuracy. Thus the evolutions required in implementing the constraint iterations are over short time intervals, with each initial guess being the converged evolution from the previous time step and the evolution representable to sufficient accuracy on a low-dimensional interpolation basis (e.g. dimension two for piecewise-linear representation of the full evolution) – the increase in difficulty in going from task 1 to task 2 may not be as great as at first it appears to be.

5.6.1. Electric current

We can also identify the electric current, $\boldsymbol {j}\stackrel {\textrm {def}}{=}\boldsymbol {\nabla }{\times }\boldsymbol {B}/\mu _0$, so (5.31) can be written as

(5.32)\begin{equation} \boldsymbol{j} =\frac{\mu_\varOmega }{\mu_0}\boldsymbol{B}+\frac{\nu_\varOmega}{\mu_0}\boldsymbol{\omega}- \boldsymbol{\nabla}{\times}(\boldsymbol{v}{\times}\boldsymbol{\lambda}_{*})+\frac{\partial\boldsymbol{\lambda}_{*}}{\partial t},\end{equation}

where $\boldsymbol {\omega }\stackrel {\textrm {def}}{=}\boldsymbol {\nabla }{\times }\boldsymbol {u}$ is the fluid vorticity. The first term on the right-hand side of (5.31) is the usual parallel electric current term of the linear-force-free (Beltrami) magnetic field model; the second term is a vorticity-driven current Yokoi (Reference Yokoi2013) term; while the last two terms form a new IOL constraint current which (taking into account the EL equation $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\lambda }_{*}=0$) maintains the divergence-free nature of $\boldsymbol {j}$ as required to maintain quasineutrality.

5.6.2. Physical interpretation of estimated Lagrange multiplier

In the special case $\mu _\varOmega =\nu _\varOmega =0$, $\boldsymbol {v}=\boldsymbol {u}$, if we make the identification $\boldsymbol {\lambda }_{*} = \boldsymbol {P}$ (5.32) becomes identical with the representation of $\boldsymbol {j}$ in terms of the electrostatic dipole moment per unit volume or polarization vector $\boldsymbol {P}$ (see e.g. §§ 1–10 of Panofsky & Phillips Reference Panofsky and Phillips1962). This representation is as given in (12) of Calkin (Reference Calkin1963) and (1.2) of Webb & Anco (Reference Webb and Anco2017), specialized to the MHD case of a quasineutral moving medium, where $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {P}=0$ (consistently with (5.20)).

Calkin (Reference Calkin1963) goes on to derive an IMHD action principle in terms of Clebsch potentials, but these are not globally defined in a 3-D plasma with non-integrable magnetic fields. Our derivation shows the Clebsch representation is not needed to apply this polarization representation for $\boldsymbol {j}$ in an action principle if we apply the Lagrangian variational approach of Newcomb (Reference Newcomb1962). (See also Webb & Anco (Reference Webb and Anco2019) for a discussion of the equivalence of Lagrangian and Eulerian variational approaches.)

5.7.1. Equations of motion

From (4.4a–c) $\partial _{\boldsymbol {v}}\boldsymbol {C}= {\boldsymbol{\mathsf{I}}} {\times }\boldsymbol {B}$, thus

(5.33)\begin{equation} \boldsymbol{\varPi} =\rho\boldsymbol{u}+\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*}.\end{equation}

The EL equation obtained from setting $\delta \mathscr {S}_\varOmega ^{A}/\delta \boldsymbol {x}=0$ in (5.8) thus becomes

(5.34)\begin{align} & \partial_{t}(\rho\boldsymbol{u}+\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*}) +\boldsymbol{\nabla}\boldsymbol{\cdot}[\boldsymbol{v}(\rho\boldsymbol{u}+\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*})] +(\boldsymbol{\nabla}\boldsymbol{v})\boldsymbol{\cdot}(\rho\boldsymbol{u}+\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*}) \nonumber\\ & \quad =\rho\boldsymbol{\nabla}\left(\frac{\partial\mathcal{L}_\varOmega^{A}}{\partial\rho}\right)= \rho\boldsymbol{\nabla}(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{v}-h_{\varOmega}), \end{align}

by (5.18). Cancelling the $\rho (\boldsymbol {\nabla }\boldsymbol {v})\boldsymbol {\cdot }\boldsymbol {u}$ occurring on both sides and rearranging, we have

(5.35)\begin{align} & \rho\partial_{t}\boldsymbol{u}+\rho\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}-\rho (\boldsymbol{\nabla}\boldsymbol{u})\boldsymbol{\cdot}\boldsymbol{v} \nonumber\\ & \quad ={-}\rho\boldsymbol{\nabla} h_{\varOmega}-\partial_{t} (\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*})-\boldsymbol{\nabla}\boldsymbol{\cdot}[\boldsymbol{v}(\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*})] -(\boldsymbol{\nabla}\boldsymbol{v})\boldsymbol{\cdot}(\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*}), \end{align}

where we used (2.4), $\partial _{t}\rho +\boldsymbol {\nabla }\boldsymbol {\cdot }(\rho \boldsymbol {v})=0$, to cancel all derivatives of $\rho$. Thus, dividing both sides by $\rho$ we have the compact Bernoulli-like form

(5.36)\begin{equation} \partial_{t}\boldsymbol{u}+\boldsymbol{\omega}{\times}\boldsymbol{v} ={-}\boldsymbol{\nabla} h_{\varOmega}-\boldsymbol{a}_{\boldsymbol{\lambda}},\end{equation}

the residual acceleration term containing $\boldsymbol {\lambda }_{*}$ being

(5.37)\begin{align} \boldsymbol{a}_{\boldsymbol{\lambda}}& \stackrel{\text{def}}{=} \rho^{{-}1}[\partial_{t} (\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*})+\boldsymbol{\nabla}\boldsymbol{\cdot}[\boldsymbol{v}(\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*})] +(\boldsymbol{\nabla}\boldsymbol{v})\boldsymbol{\cdot}(\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*})] \nonumber\\ & \quad = \partial_{t}\boldsymbol{w}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{w} +(\boldsymbol{\nabla}\boldsymbol{v})\boldsymbol{\cdot}\boldsymbol{w}, \end{align}

where

(5.38)\begin{equation} \boldsymbol{w}\stackrel{\text{def}}{=} \frac{\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*}}{\rho},\end{equation}

again using $\partial _{t}\rho +\boldsymbol {\nabla }\boldsymbol {\cdot }(\rho \boldsymbol {v})=0$.

In the special case $\mu _\varOmega =\nu _\varOmega =0$, $\boldsymbol {v}=\boldsymbol {u}$ we can use the identification in § 5.6.2 of $\boldsymbol {\lambda }_{*}$ as the polarization field $\boldsymbol {P}$ to write $\boldsymbol {w} = \boldsymbol {B}{\times }\boldsymbol {P}/{\rho }$. We can then recognize (5.36) as the Eulerian equation of motion, (23) of Calkin (Reference Calkin1963), thus providing a physical interpretation of our equations of motion in terms of a Lagrange multiplier field.

5.7.2. Conservation form

Now consider the conservation form (5.12) where, from (4.4a–c), (5.9) and (4.10)

(5.39)\begin{equation} \boldsymbol{\varPi}+\frac{\partial\mathcal{L}_\varOmega^{{A}_{-}}}{\partial\boldsymbol{E}}{\times}\boldsymbol{B} =\rho\boldsymbol{u}.\end{equation}

Starting with the coefficient of ${\boldsymbol{\mathsf{I}}}$ in the tensor $ {\boldsymbol{\mathsf{T}}} $, (5.13), and referring to (4.19), (4.20) and (4.21) we find

(5.40)\begin{equation} \mathcal{L}_\varOmega^{{A}_{-}}-\rho\frac{\partial\mathcal{L}_\varOmega^{{A}_{-}}}{\partial\rho} =p-\frac{B^{2}}{2\mu_0}+\nu_\varOmega\frac{\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{B}}{\mu_0}+ \boldsymbol{\lambda}\boldsymbol{\cdot}\boldsymbol{C}-\frac{\mu^{P}_{\varOmega} \boldsymbol{C}^{2}}{2}. \end{equation}

The penultimate term in $ {\boldsymbol{\mathsf{T}}} $ is

(5.41)\begin{equation} -\frac{\partial\mathcal{L}_\varOmega^{C}}{\partial\boldsymbol{E}}\boldsymbol{E} ={-}\boldsymbol{\lambda}_{*}\boldsymbol{E}, \end{equation}

which consists of a symmetric $\boldsymbol {E}'\boldsymbol {E}'$ term and a non-symmetric $\boldsymbol {E}'\boldsymbol {u}{\times }\boldsymbol {B}$ term.

The preceding term of ${\boldsymbol{\mathsf{T}}}$ is, using (5.28) and (B8),

(5.42)\begin{equation} \frac{\partial\mathcal{L}_\varOmega^{{A}_{-}}}{\partial\boldsymbol{B}}{\times}{\boldsymbol{\mathsf{I}}}{\times}\boldsymbol{B} -\frac{{\boldsymbol{\mathsf{I}}}}{\mu_0}({-}B^{2}+\nu_\varOmega\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{B}- \mu_0\boldsymbol{v}{\times}\boldsymbol{\lambda}_{*}\boldsymbol{\cdot}\boldsymbol{B}), \end{equation}

and the remaining first term is

(5.43)\begin{equation} \boldsymbol{v}\boldsymbol{\varPi}=\left(\boldsymbol{u}-\frac{\nu_\varOmega\boldsymbol{B}}{\mu_0\rho}\right) (\rho\boldsymbol{u}+\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*}). \end{equation}

Thus, combining all terms, (5.12) becomes

(5.44)\begin{equation} \partial_{t}(\rho\boldsymbol{u})+\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{\mathsf{T}}}_{{\rm MHD}}+{\boldsymbol{\mathsf{T}}}_{{\rm Res}}) =(\boldsymbol{\nabla}\boldsymbol{\lambda})\boldsymbol{\cdot}\boldsymbol{C},\end{equation}

where ${\boldsymbol{\mathsf{T}}}_{\textrm {MHD}}$ is the momentum transport plus stress tensor for both IMHD and RxMHD (Dewar et al. Reference Dewar, Burby, Qu, Sato and Hole2020),

(5.45)\begin{equation} {\boldsymbol{\mathsf{T}}}_{{\rm MHD}}=\rho\boldsymbol{u}\boldsymbol{u}+ \left(p+\frac{B^{2}}{2\mu_0}\right){\boldsymbol{\mathsf{I}}}-\frac{\boldsymbol{B}\boldsymbol{B}}{\mu_0}, \end{equation}

the terms in $\nu _\varOmega$ that might have contributed to ${\boldsymbol{\mathsf{T}}}_{\textrm {MHD}}$ in the RxMHD case having cancelled.

The new term ${\boldsymbol{\mathsf{T}}}_{\textrm {Res}}$ is the ‘internal’ residual stress contribution arising when action-extremizing solutions are infeasible, i.e. when the IOL constraint is not satisfied exactly,

(5.46)\begin{align} {\boldsymbol{\mathsf{T}}}_{{\rm Res}} & \stackrel{\text{def}}{=}\left(\boldsymbol{\lambda}\boldsymbol{\cdot}\boldsymbol{C}-\frac{\mu^{P}_{\varOmega} \boldsymbol{C}^{2}}{2}-\boldsymbol{\lambda}_{*}\boldsymbol{\cdot}\boldsymbol{v}{\times}\boldsymbol{B}\right){\boldsymbol{\mathsf{I}}} \nonumber\\ & \quad-\boldsymbol{B}\boldsymbol{v}{\times}\boldsymbol{\lambda}_{*}+\boldsymbol{v}\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*}-\boldsymbol{\lambda}_{*}\boldsymbol{E} \nonumber\\ & =\left(\boldsymbol{\lambda}_{*}\boldsymbol{\cdot}\boldsymbol{C}-\boldsymbol{\lambda}_{*}\boldsymbol{\cdot}\boldsymbol{u}{\times}\boldsymbol{B}+ \frac{\mu^{P}_{\varOmega}\,\boldsymbol{C}^{2}}{2}\right){\boldsymbol{\mathsf{I}}} \nonumber\\ & \quad+\boldsymbol{B}\boldsymbol{\lambda}_{*}{\times}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*}+ \boldsymbol{\lambda}_{*}\boldsymbol{u}{\times}\boldsymbol{B}-\boldsymbol{\lambda}_{*}\boldsymbol{C}. \end{align}

(Interestingly, the $\nu _\varOmega$ cancellation also occurred in deriving ${\boldsymbol{\mathsf{T}}}_{\textrm {Res}}$ when $\boldsymbol {v}$ was replaced by $\boldsymbol {u}-\nu _\varOmega \boldsymbol {B}/\mu _0\rho$.)

The ‘external’ residual force on the right-hand side of (5.44) obviously vanishes for feasible solutions. It is not obvious that ${\boldsymbol{\mathsf{T}}}_{\textrm {Res}}$ vanishes when $\boldsymbol {C}=0$ as it involves the unknown converged Lagrange multiplier $\boldsymbol {\lambda }|^{\infty }$ ($=\boldsymbol {\lambda }_{*}|^{\infty }$ as $\boldsymbol {C}|^{\infty }=0$). However, it is easy to verify that both the diagonal and off-diagonal terms of ${\boldsymbol{\mathsf{T}}}_{\textrm {Res}}$ not involving $\boldsymbol {C}$ explicitly do vanish if $\boldsymbol {\lambda }_{*}$ is proportional to $\boldsymbol {B}$ pointwise, implying at least in this case ${\boldsymbol{\mathsf{T}}}_{\textrm {Res}}=0$ if and only if $\boldsymbol {C}=0$ (assuming $\mu ^{P}_{\varOmega }\neq 0$).

5.7.3. Momentum and angular momentum conservation

When a trial solution is IOL-infeasible, i.e. $\boldsymbol {C}\neq 0$, ${\boldsymbol{\mathsf{T}}}_{\textrm {Res}}$ is not a symmetric tensor, indicating it imparts both an isotropic pressure force and a torque on the plasma, presumably tending to change $\boldsymbol {u}$ in such a way as to ‘bend’ the flow towards conformity with the IOL. There is a cyclic symmetry in ${\boldsymbol{\mathsf{T}}}_{\textrm {Res}}$  in the sum of the three terms in $\boldsymbol {\lambda }_{*},\boldsymbol {u},\boldsymbol {B}$ that indicates that the magnetic field is coupled to $\boldsymbol {E}'$ in a similar fashion to $\boldsymbol {u}$, and indeed we see from (5.32) that that there is a ‘dynamo’ term depending on $\boldsymbol {E}'$ in $\boldsymbol {j}$ that modifies $\boldsymbol {B}$, by Ampére's law.

6.1.1. Linearization of Lagrange multiplier determination

Focusing first on the novel part of the calculation we list the linearizations of immediate relevance to the augmented Lagrangian determination of the updated Lagrange multiplier field $\boldsymbol {\lambda }_{*}$.

From (1.5) and (1.6), $\boldsymbol {B}^{(1)}=\boldsymbol {\nabla }{\times }\boldsymbol {A}^{(1)}$ and $\boldsymbol {E}^{(1)}=-\boldsymbol {\nabla }\varPhi ^{(1)}-\partial _{t}\boldsymbol {A}^{(1)}$, so (4.3) becomes

(6.1)\begin{equation} \boldsymbol{C}^{(1)}={-}\boldsymbol{\nabla}\varPhi^{(1)}-\partial_{t}\boldsymbol{A}^{(1)}+\boldsymbol{u}{\times} (\boldsymbol{\nabla}{\times}\boldsymbol{A}^{(1)})+\boldsymbol{u}^{(1)}{\times}\boldsymbol{B},\end{equation}

with $\varPhi ^{(1)}$ to be determined from (5.27), which used the EL equation from the $\varPhi$ variation in its derivation. This becomes

(6.2)\begin{equation} \nabla^{2}\varPhi^{(1)} =\boldsymbol{\nabla}\boldsymbol{\cdot}[\boldsymbol{u}^{(1)}{\times}\boldsymbol{B}+ \boldsymbol{u}{\times}(\boldsymbol{\nabla}{\times}\boldsymbol{A}^{(1)})].\end{equation}

When $\boldsymbol {C}^{(1)}$ is found, the updated Lagrange multiplier is determined from the linearization of (4.7), $\boldsymbol {\lambda }_{*}^{(1)}=\boldsymbol {\lambda }^{(1)}-\mu ^{P}_{\varOmega }\boldsymbol {C}^{(1)}$.

While one occurrence of $\boldsymbol {A}^{(1)}$ in (6.2) has been eliminated by assuming Coulomb gauge, $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {A}^{(1)}=0$, it still arises in the $\boldsymbol {\nabla }{\times }\boldsymbol {A}^{(1)}$ term arising from $\boldsymbol {B}^{(1)}$. Thus we also need the linearization of the $\delta \boldsymbol {A}$ EL equation to give us $\boldsymbol {A}^{(1)}$. We use the modified Beltrami form (5.31)

(6.3)\begin{align} \boldsymbol{\nabla}{\times}(\boldsymbol{\nabla}{\times}\boldsymbol{A}^{(1)}) & =\mu_{\varOmega}\boldsymbol{\nabla}{\times}\boldsymbol{A}^{(1)}+ \nu_\varOmega\boldsymbol{\nabla}{\times}\boldsymbol{u}^{(1)}-\mu_0\boldsymbol{\nabla}{\times}(\boldsymbol{v}^{(1)}{\times}\boldsymbol{\lambda}) \nonumber\\ & \quad+\mu_0[\partial_{t}\boldsymbol{\lambda}_{*}^{(1)}-\boldsymbol{\nabla}{\times} (\boldsymbol{v}{\times}\boldsymbol{\lambda}_{*}^{(1)})], \end{align}

where

(6.4)\begin{equation} \boldsymbol{v}^{(1)}=\boldsymbol{u}^{(1)}-\frac{\nu_\varOmega}{\mu_0} \left(\frac{\boldsymbol{\nabla}{\times}\boldsymbol{A}^{(1)}}{\rho}-\frac{\rho^{(1)}}{\rho} \frac{\boldsymbol{B}}{\rho}\right),\end{equation}

with $\rho ^{(1)}$ to be determined from

(6.5)\begin{equation} \partial_{t}\rho^{(1)}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{u}^{(1)}+\rho^{(1)}\boldsymbol{u})=0,\end{equation}

and $\boldsymbol {u}^{(1)}$ to be treated as the one unknown in terms of which all other physical perturbations are to be expressed.

6.2.1. Eikonal ansatz and natural basis vectors

For short-wavelength high-frequency velocity perturbations we use the eikonal ansatz

(6.6)\begin{equation} \boldsymbol{u}^{(1)}=\tilde{\boldsymbol{u}}(\boldsymbol{x},t)\exp \left(\frac{{\rm i}\varphi(\boldsymbol{x},t)}{\varepsilon}\right),\end{equation}

with similar notations for linear perturbations of other quantities, $\varepsilon$ being the WKB (local plane wave) expansion parameter. The instantaneous local values of wavevector and frequency as seen in the laboratory frame are then defined as $\boldsymbol {k}\stackrel {\textrm {def}}{=}\boldsymbol {\nabla }\varphi$ and $\omega (\boldsymbol {x},t)\stackrel {\textrm {def}}{=}-\partial _{t}\varphi$.

In the following development we shall also encounter two ‘shifted’ frequencies:

(6.7)\begin{equation} \text{(1)}\quad \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} \stackrel{\text{def}}{=}\omega-\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{u},\end{equation}

the Doppler-shifted frequency of the wave as seen in the local rest frame of a fluid element, velocity $\boldsymbol {u}(\boldsymbol {x},t)$, and

(6.8)\begin{equation} \text{(2)}\quad \omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}\stackrel{\text{def}}{=} \omega-\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v},\end{equation}

the same as frequency (1) except with $\boldsymbol {u}$ replaced by the relative velocity $\boldsymbol {v}=\boldsymbol {u}-\boldsymbol {u}^\textrm {Rx}\equiv \boldsymbol {u}-\nu _\varOmega \boldsymbol {B/}\mu _0\rho$.

Taking $\varphi$ and equilibrium quantities to vary $O(1)$ spatial and temporal scales, $\omega$, $\boldsymbol {k}$, $\partial _{t}\boldsymbol {u}$, $\boldsymbol {\nabla }\boldsymbol {u}$, $\mu _{\varOmega }$, $\nu _\varOmega$, etc. are $O(1)$, but $\partial _{t}\boldsymbol {u}^{(1)}$, $\boldsymbol {\nabla }\boldsymbol {u}^{(1)}$, etc. are large, $O(\alpha \varepsilon ^{-1})$, relative to $\boldsymbol {u}^{(1)}=O(\alpha )\varepsilon ^{0}$, and similarly for spatio-temporal derivatives of $\rho ^{(1)}$ and $\boldsymbol {B}^{(1)}$.

In order for $\boldsymbol {B}^{(1)}$ to be the same order as $\boldsymbol {u}^{(1)}$, the potentials $\varPhi ^{(1)}$ and$\boldsymbol {A}^{(1)}$ must be $O(\alpha \varepsilon )$, so we write

(6.9a,b)\begin{equation} \varPhi^{(1)}=\frac{\varepsilon}{{\rm i}}\tilde{\varPhi}(\boldsymbol{x},t) \exp\left(\frac{{\rm i}\varphi(\boldsymbol{x},t)}{\varepsilon}\right) \quad \text{and} \quad \boldsymbol{A}^{(1)}=\frac{\varepsilon}{{\rm i}} \tilde{\boldsymbol{A}}(\boldsymbol{x},t)\exp\left(\frac{{\rm i} \varphi(\boldsymbol{x},t)}{\varepsilon}\right);\end{equation}

(1.5) and (1.6) giving then $\tilde {\boldsymbol {B}}=\boldsymbol {k}{\times }\tilde {\boldsymbol {A}}$ and $\tilde {\boldsymbol {E}}=-\boldsymbol {k} \tilde {\varPhi }+\omega \tilde {\boldsymbol {A}}$.

As in Dewar et al. (Reference Dewar, Burby, Qu, Sato and Hole2020) our strategy is to express all perturbations in terms of $\tilde {\boldsymbol {u}}$, in order to find a $3\times 3$ matrix eigenvalue equation whose roots give the dispersion relations for the three propagating wave branches. We shall also include a forcing term of the form similar to (6.6) in the equation of motion for $\boldsymbol {u}^{(1)}$ so that this $3\times 3$ matrix appears also as a response function, with the dispersion relations giving the location of its poles.

We shall find it useful to expand vectors and dyadics in the orthonormal MHD-wave basis,

(6.10a–c)\begin{equation} \boldsymbol{e}_{1}\stackrel{\text{def}}{=}\frac{\boldsymbol{k}_{{\perp}}}{k_{{\perp}}},\quad \boldsymbol{e}_{2}\stackrel{\text{def}}{=}\frac{\boldsymbol{B}}{B}\equiv\boldsymbol{e}_{\boldsymbol{B}}\quad \text{and}\quad \boldsymbol{e}_{3}\stackrel{\text{def}}{=}\equiv\boldsymbol{e}_{1}{\times}\boldsymbol{e}_{2}\equiv\frac{\boldsymbol{k} {\times}\boldsymbol{B}}{|\boldsymbol{k}{\times}\boldsymbol{B}|},\end{equation}

where $\boldsymbol {k}_{\perp }\stackrel {\textrm {def}}{=} {\boldsymbol{\mathsf{P}}} _{\perp }\boldsymbol {\cdot }\boldsymbol {k}$, so $\boldsymbol {{k}}=k_{\perp }\boldsymbol {e}_{1}+k_{\parallel }\boldsymbol {e}_{2}$ (where $k_{\perp }\stackrel {\textrm {def}}{=}|\boldsymbol {k}_{\perp }|$, $k_{\parallel }\stackrel {\textrm {def}}{=}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {e}_{\boldsymbol {B}}$), $\boldsymbol {u}^\textrm {Rx}=u^{\textrm {Rx}}\boldsymbol {e}_{2}$ (where $u^{\textrm {Rx}}\stackrel {\textrm {def}}{=}|\boldsymbol {u}^\textrm {Rx}|$), $\boldsymbol {k}{\times }\boldsymbol {B}_{{A}}=k_{\perp }c_{{A}}\boldsymbol {e}_{3}$, $ {\boldsymbol{\mathsf{P}}} _{\perp }=\boldsymbol {e}_{1}\boldsymbol {e}_{1}+\boldsymbol {e}_{3}\boldsymbol {e}_{3}$ and $ {\boldsymbol{\mathsf{I}}} =\boldsymbol {e}_{1}\boldsymbol {e}_{1}+\boldsymbol {e}_{2}\boldsymbol {e}_{2}+\boldsymbol {e}_{3}\boldsymbol {e}_{3}$. (There is of course a problem if $\boldsymbol {k}_{\perp }=0$, but we are interested in low-frequency MHD waves around $k_{\parallel }=0$ where $|k_{\perp }|$is maximal.)

6.2.2. IOL constraint in WKB approximation

We now use (6.6) and (6.9a,b) in the linearizations in § 6.1.1, working to leading order in $\varepsilon$ (for instance $\mu _{\varOmega }\boldsymbol {B}^{(1)}$ will be dropped as higher order in $\varepsilon$ than other terms in (6.3)). Then (6.1) becomes

(6.11)\begin{equation} \tilde{\boldsymbol{C}} ={-}\boldsymbol{k}(\tilde{\varPhi}-\boldsymbol{u} \boldsymbol{\cdot}\tilde{\boldsymbol{A}})+\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} \tilde{\boldsymbol{A}}+\tilde{\boldsymbol{u}}{\times}\boldsymbol{B}.\end{equation}

Also (6.2) becomes

(6.12)\begin{equation} \boldsymbol{k}^{2}\tilde{\varPhi} =\boldsymbol{k}^{2}\boldsymbol{u}\boldsymbol{\cdot}\tilde{\boldsymbol{A}}+\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}{\times}\boldsymbol{B},\end{equation}

assuming Coulomb gauge, $\boldsymbol {k}\boldsymbol {\cdot }\tilde {\boldsymbol {A}}=0$. Inserting (6.12) in (6.11) gives

(6.13)\begin{equation} \tilde{\boldsymbol{C}}=\left({\boldsymbol{\mathsf{I}}}-\frac{\boldsymbol{k}\boldsymbol{k}}{\boldsymbol{k}^{2}}\right) \boldsymbol{\cdot}\tilde{\boldsymbol{u}}{\times}\boldsymbol{B}+\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}\tilde{\boldsymbol{A}}.\end{equation}

Next, (6.3) becomes

(6.14)\begin{equation} \boldsymbol{k}^{2}\tilde{\boldsymbol{A}} ={-}\nu_\varOmega\boldsymbol{k}{\times}\tilde{\boldsymbol{u}}-\mu_0\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{v}} \boldsymbol{\lambda}+\mu_0\omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}} \tilde{\boldsymbol{\lambda}_{*}},\end{equation}

where we used $\boldsymbol {v}=\boldsymbol {u}-\boldsymbol {u}^\textrm {Rx}\equiv \boldsymbol {u}-\nu _\varOmega \boldsymbol {B/}\mu _0\rho$.

Finally, (6.5) and (6.4) become

(6.15)\begin{gather} \frac{\tilde{\rho}}{\rho} =\frac{\boldsymbol{k}\boldsymbol{\cdot} \tilde{\boldsymbol{u}}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}}, \end{gather}

(6.16)\begin{gather}\tilde{\boldsymbol{v}} =\tilde{\boldsymbol{u}}+\frac{\nu_\varOmega}{\mu_0\rho} \left(\frac{\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}} \boldsymbol{B}-\boldsymbol{k}{\times}\tilde{\boldsymbol{A}}\right), \end{gather}

hence

(6.17)\begin{equation} \boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{v}} =\left(\frac{\omega_{\boldsymbol{k}}^{\prime \boldsymbol{v}}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}}\right) \boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}.\end{equation}

Substituting (6.17) in (6.14) gives

(6.18)\begin{equation} \boldsymbol{k}^{2}\tilde{\boldsymbol{A}} ={-}\nu_\varOmega\boldsymbol{k}{\times}\tilde{\boldsymbol{u}}-\mu_0 \left(\frac{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}}{\omega_{\boldsymbol{k}}^{\prime \boldsymbol{u}}}\right)\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}\,\boldsymbol{\lambda}+\mu_0\omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}} \tilde{\boldsymbol{\lambda}_{*}},\end{equation}

which in (6.13) then gives

(6.19)\begin{equation} \tilde{\boldsymbol{C}} =\left({\boldsymbol{\mathsf{I}}}-\frac{\boldsymbol{k}\boldsymbol{k}}{\boldsymbol{k}^{2}} \right)\boldsymbol{\cdot}\tilde{\boldsymbol{u}}{\times}\boldsymbol{B}-\nu_\varOmega\omega_{\boldsymbol{k}}^{\prime \boldsymbol{u}}\frac{\boldsymbol{k}{\times}\tilde{\boldsymbol{u}}}{\boldsymbol{k}^{2}}+\mu_0 \frac{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}}{\boldsymbol{k}^{2}} (\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}\tilde{\boldsymbol{\lambda}_{*}}- \boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}\boldsymbol{\lambda}).\end{equation}

Treating $\tilde {\boldsymbol {C}}$ for the moment as a known and solving for $\tilde {\boldsymbol {\lambda }_{*}}$ we have

(6.20)\begin{equation} \tilde{\boldsymbol{\lambda}_{*}} =\frac{\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}}{\omega_{\boldsymbol{k}}^{\prime \boldsymbol{u}}}\boldsymbol{\lambda}+\frac{\nu_\varOmega}{\mu_0} \frac{\boldsymbol{k}{\times} \tilde{\boldsymbol{u}}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}}+ \frac{(\boldsymbol{k}^{2}{\boldsymbol{\mathsf{I}}}-\boldsymbol{k}\boldsymbol{k}) {\times}\boldsymbol{B}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}+\boldsymbol{k}^{2}\tilde{\boldsymbol{C}}}{\mu_0 \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} \omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}}.\end{equation}

Normally we do not need to know the residual IOL error term $\tilde {\boldsymbol {C}}$ exactly, but to convince ourselves it can be made arbitrarily small by iteration, replace $\tilde {\boldsymbol {\lambda }_{*}}$ with its explicit form from the linearization of (4.7), $\tilde {\boldsymbol {\lambda }}-\mu ^{P}_{\varOmega }\tilde {\boldsymbol {C}}$, and collect both $\tilde {\boldsymbol {C}}$ terms on the left,

(6.21)\begin{align} \left(1+\mu_0\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}\omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}} \frac{\mu^{P}_{\varOmega}}{\boldsymbol{k}^{2}}\right)\tilde{\boldsymbol{C}} & ={-}\left({\boldsymbol{\mathsf{I}}}- \frac{\boldsymbol{k}\boldsymbol{k}}{\boldsymbol{k}^{2}}\right){\times}\boldsymbol{B}\boldsymbol{\cdot} \tilde{\boldsymbol{u}}-\nu_\varOmega\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} \frac{\boldsymbol{k}{\times}\tilde{\boldsymbol{u}}}{\boldsymbol{k}^{2}} \nonumber\\ & \quad+\mu_0\frac{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}}{\boldsymbol{k}^{2}} (\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}\tilde{\boldsymbol{\lambda}}-\boldsymbol{k} \boldsymbol{\cdot}\tilde{\boldsymbol{u}}\,\boldsymbol{\lambda}), \end{align}

which confirms the implication in subsection (3.3) that we have enough equations to determine $\tilde {\boldsymbol {C}}$, and hence $\tilde {\varPhi }$, $\tilde {\boldsymbol {A}}$ and $\tilde {\boldsymbol {\lambda }_{*}}$, in terms of $\tilde {\boldsymbol {u}}$.

Dividing both sides of (6.21) by the large penalty multiplier $\mu ^{P}_{\varOmega }$ we see that $\tilde {\boldsymbol {C}}$ is smaller than the other terms and the previous iterate of $\tilde {\boldsymbol {C}}$ by an $O(1/\mu ^{P}_{\varOmega })$ factor. Thus the sequence $\{ \ldots,\tilde {\boldsymbol {C}}^{n},\tilde {\boldsymbol {C}}^{n+1},\tilde {\boldsymbol {C}}^{n+2},\ldots \}$ will converge exponentially towards $0$, or superexponentially if $\mu ^{P}_{\varOmega }$ is increased appropriately at each step. Also $\tilde {\boldsymbol {\lambda }_{*}}$ will converge to $\tilde {\boldsymbol {\lambda }}|^{\infty }$.

However, this linearized calculation is sufficiently simple that we do not actually need to carry out the iteration as we can find $\tilde {\boldsymbol {\lambda }}|^{\infty }$analytically from 6.21 by setting its left-hand side to zero and solving for $\tilde {\boldsymbol {\lambda }}=\tilde {\boldsymbol {\lambda }}|^{\infty }$. Or, if we want to investigate the hypothesis that terminating the iteration at finite $n$, so that $\tilde {\boldsymbol {C}}\neq 0,$ will regularize MHD, we can prescribe $\tilde {\boldsymbol {C}}$ and use (6.20) to give $\tilde {\boldsymbol {\lambda }_{*}}$. To provide a continuous sequence of dynamical fluid models running between the unconstrained RxMHD perturbation dynamics of Dewar et al. (Reference Dewar, Burby, Qu, Sato and Hole2020) to the converged, $\tilde {\boldsymbol {C}}=0$, present model we choose

(6.22)\begin{equation} \tilde{\boldsymbol{C}}={\varepsilon}_{\mathrm{Rx}}\tilde{\boldsymbol{C}}_{0},\end{equation}

with the ‘relaxedness’ parameter ${\varepsilon }_{\mathrm {Rx}}$ running from 0 (IMHD, IOL-compliant) to 1 (RxMHD, may be IOL-infeasible). We shall later also have use of the complementary ‘ideality’ parameter ${\varepsilon }_{{I}}\stackrel {\textrm {def}}{=}1-{\varepsilon }_{\mathrm {Rx}}$.

Here $\tilde {\boldsymbol {C}}_{0}$ is the IOL error for unconstrained, ${\varepsilon }_{{I}}=0, {\varepsilon }_{\mathrm {Rx}}=1$, RxMHD perturbations, which we can find by setting $\tilde {\boldsymbol {\lambda }}=0$, $\mu ^{P}_{\varOmega }=0$, and thus $\tilde {\boldsymbol {\lambda }_{*}}=0$, in (6.19) to give

(6.23)\begin{equation} \widetilde{\boldsymbol{C}_{0}} =\left({\boldsymbol{\mathsf{I}}}-\frac{\boldsymbol{k}\boldsymbol{k}}{\boldsymbol{k}^{2}}\right) \boldsymbol{\cdot}\tilde{\boldsymbol{u}}{\times}\boldsymbol{B}-\nu_\varOmega\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} \frac{\boldsymbol{k}{\times}\tilde{\boldsymbol{u}}}{\boldsymbol{k}^{2}}-\mu_0\boldsymbol{\lambda} \omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}\frac{\boldsymbol{k} \boldsymbol{\cdot}\tilde{\boldsymbol{u}}}{\boldsymbol{k}^{2}},\end{equation}

which is a linear tensor function of the form $\widetilde {\boldsymbol {C}_{0}}=B {\boldsymbol{\mathsf{C}}} _{0}(\boldsymbol {k},\boldsymbol {B})\boldsymbol {\cdot }\tilde {\boldsymbol {u}}$, where the factor $B$ is taken out to make $ {\boldsymbol{\mathsf{C}}} _{0}$ dimensionless. By inspection of (6.23),

(6.24)\begin{equation} {\boldsymbol{\mathsf{C}}}_{0}\stackrel{\text{def}}{=}-\frac{1}{\boldsymbol{k}^{2}B}[\mu_0 \omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}\boldsymbol{\lambda}\boldsymbol{k}+\nu_\varOmega \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}\boldsymbol{k}{\times}{\boldsymbol{\mathsf{I}}}+ (\boldsymbol{k}^{2}{\boldsymbol{\mathsf{I}}}-\boldsymbol{k}\boldsymbol{k}){\times}\boldsymbol{B}].\end{equation}

The linear tensor form of $\widetilde {\boldsymbol {C}_{0}}$ implies $\tilde {\boldsymbol {C}}$, $\tilde {\boldsymbol {\lambda }}$ and $\tilde {\boldsymbol {\lambda }_{*}}$ are of similar form,

(6.25a–c)\begin{equation} \tilde{\boldsymbol{C}}=B {\boldsymbol{\mathsf{C}}}(\boldsymbol{k},\boldsymbol{B})\boldsymbol{\cdot}\tilde{\boldsymbol{u}}, \quad\tilde{\boldsymbol{\lambda}}=\boldsymbol{\varLambda}(\boldsymbol{k},\boldsymbol{B})\boldsymbol{\cdot} \tilde{\boldsymbol{u}},\quad \tilde{\boldsymbol{\lambda}_{*}}=\boldsymbol{\varLambda}_{*}(\boldsymbol{k}, \boldsymbol{B})\boldsymbol{\cdot}\tilde{\boldsymbol{u}},\end{equation}

where from $\boldsymbol {\varLambda }=\boldsymbol {\varLambda }_{*}+\mu ^{P}_{\varOmega } B {\boldsymbol{\mathsf{C}}} $, and from (6.20) and (6.25a–c),

(6.26)\begin{equation} \boldsymbol{\varLambda}_{*} ={\varepsilon}_{{I}}\left[\frac{\boldsymbol{\lambda} \boldsymbol{k}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}} +\frac{\nu_\varOmega}{\mu_0} \frac{\boldsymbol{k}{\times}{\boldsymbol{\mathsf{I}}}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}} +\frac{(\boldsymbol{k}^{2}{\boldsymbol{\mathsf{I}}}-\boldsymbol{k}\boldsymbol{k}){\times}\boldsymbol{B}}{\mu_0 \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} \omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}}\right] ,\end{equation}

using $ {\boldsymbol{\mathsf{C}}} ={\varepsilon }_{\mathrm {Rx}} {\boldsymbol{\mathsf{C}}} _{0}$, 6.22.

6.2.3. Short-wavelength dynamical RxMHD equations

We now consider the linearized equation of motion with the forcing term mentioned in § 6.2.1, a specific force (i.e. force/mass density) we denote as $\boldsymbol {f}^{(1)}$. Thus the linearized (5.36) with forcing term becomes

(6.27)\begin{equation} \partial_{t}\boldsymbol{u}^{(1)}+\boldsymbol{\omega}^{(1)}{\times}\boldsymbol{v} +\boldsymbol{\omega}{\times}\boldsymbol{v}^{(1)}={-}\boldsymbol{a}_{\boldsymbol{\lambda}}^{(1)}- \boldsymbol{\nabla} h_{\varOmega}^{(1)}+\boldsymbol{f}^{(1)},\end{equation}

where, from (5.19),

(6.28)\begin{equation} h_{\varOmega}^{(1)} =\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{u}^{(1)}+{{\tau_{\varOmega}}} \frac{\rho^{(1)}}{\rho}, \end{equation}

from (5.37),

(6.29)\begin{align} \boldsymbol{a}_{\boldsymbol{\lambda}}^{(1)}& = \partial_{t}\boldsymbol{w}^{(1)}+\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{w}^{(1)}+(\boldsymbol{\nabla}\boldsymbol{v})\boldsymbol{\cdot}\boldsymbol{w}^{(1)} \nonumber\\ & \quad +\boldsymbol{v}^{(1)}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{w}+(\boldsymbol{\nabla}\boldsymbol{v}{}^{(1)}) \boldsymbol{\cdot}\boldsymbol{w}, \end{align}

and, from (5.38),

(6.30)\begin{equation} \boldsymbol{w}^{(1)}= \frac{\boldsymbol{B}^{(1)}{\times}\boldsymbol{\lambda}}{\rho}+ \frac{\boldsymbol{B}{\times}\boldsymbol{\lambda}_{*}^{(1)}}{\rho}-\frac{\rho^{(1)} \boldsymbol{B}{\times}\boldsymbol{\lambda}}{\rho^{2}}.\end{equation}

Using the WKB representations in (6.9a,b) and in (6.6), and analogous representations for $\rho ^{(1)}$, $\boldsymbol {B}^{(1)}$, $\boldsymbol {w}^{(1)}$, $\boldsymbol {a}_{\boldsymbol {\lambda }}^{(1)}$ and $\boldsymbol {f}^{(1)}$ in the linearizations above, and working to leading order in $\varepsilon$ as before we have

(6.31)\begin{equation} -\!\omega\tilde{\boldsymbol{u}}+(\boldsymbol{k}{\times}\tilde{\boldsymbol{u}}) {\times}\boldsymbol{v} ={-}\boldsymbol{k}\left(\boldsymbol{u}+\tau_{\varOmega} \frac{\boldsymbol{k}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}}\right)\boldsymbol{\cdot} \tilde{\boldsymbol{u}}-\tilde{\boldsymbol{a}}_{\boldsymbol{\lambda}} +\tilde{\boldsymbol{f}},\end{equation}

which we shall show can be written as

(6.32)\begin{equation} {\boldsymbol{\mathsf{D}}}(\omega,\boldsymbol{k})\boldsymbol{\cdot}\tilde{\boldsymbol{u}}={-}\tilde{\boldsymbol{f}},\end{equation}

where, noting cancelling of $\boldsymbol {k}\boldsymbol {u}\boldsymbol {\cdot }\tilde {\boldsymbol {u}}$ terms,

(6.33)\begin{equation} {\boldsymbol{\mathsf{D}}} \boldsymbol{\cdot}\tilde{\boldsymbol{u}} =\omega\tilde{\boldsymbol{u}}- \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{u}\tilde{\boldsymbol{u}}-\boldsymbol{u}^{\rm Rx}{\times} (\boldsymbol{k}{\times}\tilde{\boldsymbol{u}})-\tilde{\boldsymbol{a}}_{\boldsymbol{\lambda}} -{{\tau_{\varOmega}}}\frac{\boldsymbol{k}\boldsymbol{k}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}} \boldsymbol{\cdot}\tilde{\boldsymbol{u}}.\end{equation}

From (6.29),

(6.34)\begin{equation} \tilde{\boldsymbol{a}}_{\boldsymbol{\lambda}}={-}\omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}} \tilde{\boldsymbol{w}}+\boldsymbol{k}\tilde{\boldsymbol{v}}\boldsymbol{\cdot}\boldsymbol{w},\end{equation}

and, from (6.30), (6.15) and (6.18),

(6.35)\begin{align} \rho\tilde{\boldsymbol{w}}& = \left[\frac{\boldsymbol{k}}{\boldsymbol{k}^{2}}{\times} \left(-\nu_\varOmega\boldsymbol{k}{\times}\tilde{\boldsymbol{u}}-\mu_0\left(\frac{\omega_{\boldsymbol{k}}^{\prime \boldsymbol{v}}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}}\right)\boldsymbol{k}\boldsymbol{\cdot} \tilde{\boldsymbol{u}}\boldsymbol{\lambda}\right)-\frac{\boldsymbol{k}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}}{\omega_{\boldsymbol{k}}^{\prime \boldsymbol{u}}}\boldsymbol{B}\right]{\times}\boldsymbol{\lambda} \nonumber\\ & \quad + \left[\frac{\mu_0\omega_{\boldsymbol{k}}^{\prime\boldsymbol{v}}}{\boldsymbol{k}^{2}} (\boldsymbol{\lambda}\boldsymbol{\cdot}\boldsymbol{k}\,{\boldsymbol{\mathsf{I}}}-\boldsymbol{k}\boldsymbol{\lambda})+\boldsymbol{B}{\times}{\boldsymbol{\mathsf{I}}} \right]\boldsymbol{\cdot}\boldsymbol{\varLambda}_{*}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}, \end{align}

where $\boldsymbol {\varLambda }_{*}$ was given in (6.26).

6.2.4. WKB RxMHD response matrix

To simplify the calculation of the response we expand around a relaxed equilibrium with $\boldsymbol {\lambda }=0$, such as the axisymmetric tokamak equilibrium in Dewar et al. (Reference Dewar, Burby, Qu, Sato and Hole2020), which has a steady flow field $\boldsymbol {u}$ that is the vector sum of an arbitrary rigid toroidal rotation carried by $\boldsymbol {v}$ and an axisymmetric magnetic-field-aligned flow $\boldsymbol {u}^\textrm {Rx}$ proportional to $\nu _\varOmega$, which equilibrium was shown to satisfy the IOL without needing a Lagrange multiplier.

In such a case $\boldsymbol {w}=0$, and (6.35) becomes $\tilde {\boldsymbol {w}}=\boldsymbol {B}{\times }\boldsymbol {\varLambda }_{*}\boldsymbol {\cdot }\tilde {\boldsymbol {u}}/\rho$. Then (6.34) becomes

(6.36)\begin{equation} \tilde{\boldsymbol{a}}_{\boldsymbol{\lambda}}={-}{\varepsilon}_{{I}}\boldsymbol{u}^{\rm Rx}{\times}(\boldsymbol{k} {\times}\tilde{\boldsymbol{u}})-\frac{{\varepsilon}_{{I}}\boldsymbol{B}{\times}(k^{2}{\boldsymbol{\mathsf{I}}}- \boldsymbol{k}\boldsymbol{k}){\times}\boldsymbol{B}}{\mu_0\rho \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}.\end{equation}

In (6.32) the $\boldsymbol {u}^\textrm {Rx}$ terms not involving ${\varepsilon }_{\mathrm {Rx}}$ cancel, giving

(6.37)\begin{align} {\boldsymbol{\mathsf{D}}} & =(\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}+{\varepsilon}_{\mathrm{Rx}}\boldsymbol{k} \boldsymbol{\cdot}\boldsymbol{u}^{\rm Rx}){\boldsymbol{\mathsf{I}}}-{\varepsilon}_{\mathrm{Rx}}\boldsymbol{k}\boldsymbol{u}^{\rm Rx}+{\varepsilon}_{{I}}\frac{\boldsymbol{k}_{{\perp}} {\times}\boldsymbol{c}_{{A}}\boldsymbol{k}_{{\perp}}{\times}\boldsymbol{c}_{{A}} -k^{2}c_{{A}}^{2}{\boldsymbol{\mathsf{P}}}_{{\perp}}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}} \nonumber\\ & \quad-c_{{s}}^{2}\frac{\boldsymbol{k}\boldsymbol{k}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}}, \end{align}

where $\boldsymbol {u}^\textrm {Rx}=\nu _\varOmega \boldsymbol {B}/\mu _0\rho$ is defined in (1.7), $\boldsymbol {c}_{{A}}\stackrel {\textrm {def}}{=}\boldsymbol {B}/(\mu _0\rho )^{1/2}$ is the Alfvén velocity, $c_{{s}}={{\tau _{\varOmega }^{1/2}}}$ is the isothermal sound speed, and we have used (B8) to write $\boldsymbol {c}_{{A}}{\times } {\boldsymbol{\mathsf{I}}} {\times }\boldsymbol {c}_{{A}}=\boldsymbol {c}_{{A}}\boldsymbol {c}_{{A}}-c_{{A}}^{2} {\boldsymbol{\mathsf{I}}} \equiv -c_{{A}}^{2} {\boldsymbol{\mathsf{P}}} _{\perp }$.

To represent $ {\boldsymbol{\mathsf{D}}} $ as a matrix we project onto the orthonormal basis (6.10a–c), which can be written as $\boldsymbol {e}_{1}=\boldsymbol {k}_{\perp }/k_{\perp },\boldsymbol {e}_{2}=\boldsymbol {c}_{{A}}/c_{{A}}\textrm { and }\boldsymbol {e}_{3}=\boldsymbol {k}_{\perp }{\times }\boldsymbol {c}_{{A}}/(k_{\perp }c_{{A}}).$ We thus have

(6.38)\begin{align} {\boldsymbol{\mathsf{D}}} & =(\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}+ {\varepsilon}_{\mathrm{Rx}} k_{{\parallel}}u^{{\rm Rx}}){\boldsymbol{\mathsf{I}}}-{\varepsilon}_{\mathrm{Rx}} u^{{\rm Rx}} (k_{{\perp}}\boldsymbol{e}_{1}+k_{{\parallel}}\boldsymbol{e}_{2})\boldsymbol{e}_{2}- {\varepsilon}_{{I}}\frac{c_{{A}}^{2}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}} (k^{2}\boldsymbol{e}_{1}\boldsymbol{e}_{1}+k_{{\parallel}}^{2}\boldsymbol{e}_{3}\boldsymbol{e}_{3}) \nonumber\\ & \quad-\frac{c_{{s}}^{2}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}} (k_{{\perp}}\boldsymbol{e}_{1}+k_{{\parallel}}\boldsymbol{e}_{2}) (k_{{\perp}}\boldsymbol{e}_{1}+k_{{\parallel}}\boldsymbol{e}_{2}), \end{align}

which can be represented as the block-diagonal matrix

(6.39)\begin{equation} {\boldsymbol{\mathsf{D}}}=\begin{bmatrix} {\boldsymbol{\mathsf{D}}}_{{\rm MS}} & 0\\ 0 & \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}+{\varepsilon}_{\mathrm{Rx}} k_{{\parallel}}u^{{\rm Rx}}- {\varepsilon}_{{I}} k_{{\parallel}}^{2}c_{{A}}^{2}/\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} \end{bmatrix},\end{equation}

with the $1\times 1$ Alfvén block on the lower right and the $2\times 2$ magnetosonic (MS) block,

(6.40)\begin{equation} {\boldsymbol{\mathsf{D}}}_{{\rm MS}}=\begin{bmatrix} \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} +{\varepsilon}_{\mathrm{Rx}} k_{{\parallel}}u^{{\rm Rx}}-({\varepsilon}_{{I}} k^{2}c_{{A}}^{2} +k_{{\perp}}^{2}c_{{s}}^{2})/\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} & -k_{{\perp}}({\varepsilon}_{\mathrm{Rx}} u^{{\rm Rx}}+k_{{\parallel}}c_{{s}}^{2}/ \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}})\\ -k_{{\parallel}}k_{{\perp}}c_{{s}}^{2}/\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} & \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}-k_{{\parallel}}^{2}c_{{s}}^{2}/ \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}} \end{bmatrix},\end{equation}

upper left.

6.2.5. Limiting cases

Consider first the ideal, fully converged case $\boldsymbol {C}=0$ (${\varepsilon }_{{I}}=1, {\varepsilon }_{\mathrm {Rx}}=0$) and use (B8) to write $k^{2} {\boldsymbol{\mathsf{I}}} -\boldsymbol {k}\boldsymbol {k}=-\boldsymbol {k}{\times } {\boldsymbol{\mathsf{I}}} {\times }\boldsymbol {k}$ so

(6.41)\begin{align} {\boldsymbol{\mathsf{D}}} & =\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}{\boldsymbol{\mathsf{I}}}- \frac{\boldsymbol{B}{\times}(\boldsymbol{k}{\times}{\boldsymbol{\mathsf{I}}}{\times}\boldsymbol{k}) {\times}\boldsymbol{B}}{\mu_0\rho \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}} -{{\tau_{\varOmega}}}\frac{\boldsymbol{k}\boldsymbol{k}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}}\nonumber\\ & =\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}{\boldsymbol{\mathsf{I}}}- \frac{(\boldsymbol{k}\boldsymbol{B}-\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{B}{\boldsymbol{\mathsf{I}}})\boldsymbol{\cdot} (\boldsymbol{B}\boldsymbol{k}-\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{B}{\boldsymbol{\mathsf{I}}})}{\mu_0\rho \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}}-{{\tau_{\varOmega}}} \frac{\boldsymbol{k}\boldsymbol{k}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}}, \end{align}

which, apart from the definitions of $ {\boldsymbol{\mathsf{D}}} $ differing by a factor of $\rho \omega _{\boldsymbol {k}}^{\prime \boldsymbol {u}}$, agrees with the IMHD form, (75), of Dewar et al. (Reference Dewar, Burby, Qu, Sato and Hole2020).

In the pure RxMHD case ${\varepsilon }_{{I}}=0, {\varepsilon }_{\mathrm {Rx}}=1$,

(6.42)\begin{equation} {\boldsymbol{\mathsf{D}}}=\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}{\boldsymbol{\mathsf{I}}}-(\boldsymbol{k}\boldsymbol{u}^{\rm Rx}- \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{u}^{\rm Rx}{\boldsymbol{\mathsf{I}}})-{{\tau_{\varOmega}}} \frac{\boldsymbol{k}\boldsymbol{k}}{\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}}. \end{equation}

Apart from the definitions of $ {\boldsymbol{\mathsf{D}}} $ again differing by a factor of $\rho \omega _{\boldsymbol {k}}^{\prime \boldsymbol {u}}$, this agrees with the RxMHD form, (88), of Dewar et al. (Reference Dewar, Burby, Qu, Sato and Hole2020).

Thus ${\varepsilon }_{\mathrm {Rx}}$ parameterizes a continuous interpolation between RxMHD and IMHD.

6.2.6. Dispersion relations – Alfvén branches

Multiplying the first factor of the determinant

(6.43)\begin{equation} \det {\boldsymbol{\mathsf{D}}}=(\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}+ {\varepsilon}_{\mathrm{Rx}} k_{{\parallel}}u^{{\rm Rx}}-{\varepsilon}_{{I}} k_{{\parallel}}^{2}c_{{A}}^{2}/ \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}})\det {\boldsymbol{\mathsf{D}}}_{{\rm MS}},\end{equation}

by $\omega _{\boldsymbol {k}}^{\prime \boldsymbol {u}}$ gives the dispersion relation for the Alfvén-wave branch(es) as the quadratic equation

(6.44)\begin{equation} (\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}})^{2}+{\varepsilon}_{\mathrm{Rx}} k_{{\parallel}} u^{{\rm Rx}}\omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}-(1-{\varepsilon}_{\mathrm{Rx}}) k_{{\parallel}}^{2}c_{{A}}^{2}=0.\end{equation}

The general solution of the quadratic equation is

(6.45)\begin{equation} \omega_{\boldsymbol{k}}^{\prime\boldsymbol{u}}=\frac{k_{{\parallel}}}{2} [-{\varepsilon}_{\mathrm{Rx}} u^{{\rm Rx}}\pm(4{\varepsilon}_{{I}} c_{{A}}^{2}+ ({\varepsilon}_{\mathrm{Rx}} u^{{\rm Rx}})^{2})^{1/2}].\end{equation}

Qualitative analysis is more informative: as ${\varepsilon }_{\mathrm {Rx}}\to 0$ the dispersion relations for the two branches approach the Doppler-shifted Alfvén-wave dispersion relations $\omega -\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {u}=\pm k_{\parallel }c_{{A}}$. Also, inspection shows that $\omega -\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {u}\to 0$ as $k_{\parallel }\to 0$ quite generally, and when $|{\varepsilon }_{\mathrm {Rx}}|\ll 1$ the modification of the dispersion departure from the standard Alfvén-wave dispersion relation is essentially determined by the product ${\varepsilon }_{\mathrm {Rx}} u^{\textrm {Rx}}$. Thus, when $|{\varepsilon }_{\mathrm {Rx}}|\ll 1$ and the parallel flow parameter is at most Alfvénic, $u^{\textrm {Rx}}/c_{{A}}\leq O(1)$, ${\varepsilon }_{\mathrm {Rx}}$ and $u^{\textrm {Rx}}$ will have little effect on the Alfvén-wave branches.

The plots in figure 2 give a visualization of the dependence of the dispersion relation on ${\varepsilon }_{\mathrm {Rx}}$. The figure is for a case where $\boldsymbol {u}^\textrm {Rx}=0$, when (6.45) simplifies to $\omega _{\boldsymbol {k}}^{\prime \boldsymbol {u}}=\pm \sqrt {{\varepsilon }_{{I}}}\,k_{\parallel }c_{{A}}$. (We call the $+$ solution the principal branch.) The square root term $\sqrt {{\varepsilon }_{{I}}}=\sqrt {1-{\varepsilon }_{\mathrm {Rx}}}$ gives rise to a singular dependence on ${\varepsilon }_{\mathrm {Rx}}$ at ${\varepsilon }_{\mathrm {Rx}}=1$ but the vicinity of IMHD is regular.

Figure 2.

This figure shows the transition of the Alfvén-branch dispersion relation: $\omega$ versus $k_{\parallel }$; IMHD (${\varepsilon }_{\mathrm {Rx}}=0$) to RxMHD (${\varepsilon }_{\mathrm {Rx}}=1$; we have used ${\varepsilon }_{\mathrm {Rx}}=0.99$ for clarity as the ${\varepsilon }_{\mathrm {Rx}}=1$ line coincides with the $k_{\parallel }$ axis). Fixed parameters (in arbitrary units) are $k=1$, $\boldsymbol {u}=0$, $\boldsymbol {u}^\textrm {Rx}=0$, $c_{{s}}=1$ and $c_{{A}}=5$. (Colour online – the vertical ordering of the lines in the $k_{\parallel }>0$ half-plane coincides with that of ${\varepsilon }_{\mathrm {Rx}}$ in the legend.)

6.2.7. Dispersion relations – magnetosonic branches

The magnetosonic dispersion relations are obtained by setting $\det {\boldsymbol{\mathsf{D}}} _{\textrm {MS}}=0$, where

(6.46)\begin{equation} \det {\boldsymbol{\mathsf{D}}}_{{\rm MS}}=\frac{c_{{A}}^{2}{\varepsilon}_{{I}} k^{2} (k_{{\parallel}}^{2}c_{{s}}^{2}-\omega^{2})-\omega (c_{{s}}^{2}k^{2}-\omega^{2})({\varepsilon}_{\mathrm{Rx}} k_{{\parallel}}u^{{\rm Rx}}+\omega)}{\omega^{2}}. \end{equation}

The solution of the quartic equation $\omega ^{2}\det {\boldsymbol{\mathsf{D}}} _{\textrm {MS}}=0$ is extremely complicated but figures 3 and 4 give an overview of the ${\varepsilon }_{\mathrm {Rx}}$ dependence. Again, the limit ${\varepsilon }_{\mathrm {Rx}}\to 1$ is clearly singular in the slow magnetosonic case but not ${\varepsilon }_{\mathrm {Rx}}\to 0$. This regularity around IMHD means our dispersion relation analysis is too crude to reveal the potential regularizing effect of softening the IOL constraint.

Figure 3.

This figure shows the transition of the slow-magnetosonic-branch dispersion relation: $\omega$ versus $k_{\parallel }$; from IMHD (${\varepsilon }_{\mathrm {Rx}}=0$) to RxMHD (${\varepsilon }_{\mathrm {Rx}}=1$). Fixed parameters (in arbitrary units) are $k=1$, $\boldsymbol {u}=0$, $\boldsymbol {u}^\textrm {Rx}=0$, $c_{{s}}=1$ and $c_{{A}}=5$. (Colour online – the vertical ordering of the lines in the $k_{\parallel }>0$ half-plane coincides with that of ${\varepsilon }_{\mathrm {Rx}}$ in the legend.)

Figure 4.

This figure shows the transition of the fast-magnetosonic-branch dispersion relation: $\omega$ versus $k_{\parallel }$; from IMHD (${\varepsilon }_{\mathrm {Rx}}=0$) to RxMHD (${\varepsilon }_{\mathrm {Rx}}=1$). Fixed parameters (in arbitrary units) are $k=1$, $\boldsymbol {u}=0$, $\boldsymbol {u}^\textrm {Rx}=0$, $c_{{s}}=1$ and $c_{{A}}=5$. (Colour online – the vertical ordering of the lines in the $k_{\parallel }>0$ half-plane coincides with that of ${\varepsilon }_{\mathrm {Rx}}$ in the legend.)