An intuitive form of the energy functional for resistive plasma perturbations

Abstract

The energy functional of ideal magnetohydrodynamics (MHD) is extended to the regime of finite resistivity and transformed into an intuitive form, which separates pressure-gradient, current-density drive and stabilising phenomena, analogously to the ideal version of the functional derived by Greene and Johnson (Plasma Phys., vol. 10, issue 8, 1968, pp. 729–745). The resulting resistive energy functional is not in the form of a solvable eigenvalue problem but is intended to analyse the energetic composition of eigenfunctions which were calculated in advance by solving the MHD equations. The new resistive energy functional provides information on the different energetic drives and stabilising phenomena of an eigenfunction. The new energy functional has been implemented in the CASTOR3D code in general 3D curvilinear coordinates and successfully validated across three different coordinate systems for a numerical test case.

1. Introduction

Instabilities of magnetohydrodynamic (MHD) equilibria represent a major limitation of the operational space of magnetically confined fusion plasmas. There are various methods for the analysis of linear MHD stability which grant different levels of insight depending on the complexity of the problem. They range from analytical stability criteria, which are restricted to certain types of instability but reveal the detailed underlying physics, to numerical solutions of eigenvalue problems or variational principles. For a single equilibrium configuration, the eigenvalue solutions grant information on the structure of the instability and its growth rate. Scans over parameters such as resistivity and density or different equilibrium configurations yield additional information. The purpose of this paper is to provide expressions for the different forms of energy which result from a small perturbation of the MHD equilibrium in the framework of resistive MHD. These expressions are meant to be applied to the numerical solutions of the corresponding eigenvalue problem in order to gain additional insight into the mechanisms which drive or stabilise/damp these instabilities. The obtained energy decomposition of the perturbations contributes to their physical interpretation and understanding and allows one to distinguish between, e.g., current-density- and pressure-gradient-driven instabilities.

Resistive MHD is described by the following set of equations (Goedbloed & Poedts 2004):

(1.1)\begin{gather} \partial_t \rho + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{v}) = 0, \end{gather}

(1.2)\begin{gather} \rho (\partial_t + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{v} ={-} \boldsymbol{\nabla} p + \boldsymbol{j} \times \boldsymbol{B}, \end{gather}

(1.3)\begin{gather} \boldsymbol{E} + \boldsymbol{v} \times \boldsymbol{B} = \eta \boldsymbol{j}, \end{gather}

(1.4)\begin{gather} (\partial_t + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}) (p \rho^{-\varGamma}) = (\varGamma-1)\eta \boldsymbol{j}^{2}, \end{gather}

(1.5)\begin{gather} \partial_t \boldsymbol{B} ={-} \boldsymbol{\nabla} \times \boldsymbol{E}, \end{gather}

(1.6a,b)\begin{gather} \boldsymbol{\nabla} \times \boldsymbol{B} = \mu_0 \boldsymbol{j}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{B} = 0 \end{gather}

with plasma density $\rho$, velocity $\boldsymbol {v}$, pressure $p$, magnetic field $\boldsymbol {B}$, current density $\boldsymbol {j}$, electric field $\boldsymbol {E}$, resistivity $\eta$ and adiabatic coefficient $\varGamma$. The Ohmic heating term $(\varGamma -1)\eta \boldsymbol {j}^{2}$ in (1.4) is usually neglected by linear resistive MHD stability codes (Chu et al. 1995; Kerner et al. 1998). In the following, we neglect the Ohmic heating term in order to be consistent with the resistive MHD equations as they are implemented in the CASTOR and CASTOR3D codes (Kerner et al. 1998; Strumberger & Günter 2017, 2019). However, we remark that there are instabilities that are affected by Ohmic heating in the limit of extremely low conductivity or high wave numbers (Furth, Killeen & Rosenbluth 1963; Zhu et al. 2013).

The boundary conditions at the plasma–vacuum interface are given by (Goedbloed & Poedts 2004)

(1.7)\begin{gather} \left[\!\!\left[p + \frac{\boldsymbol{B}^{2}}{2 \mu_0} \right]\!\!\right]= 0, \end{gather}

(1.8)\begin{gather} \boldsymbol{n} \times [\![\boldsymbol{E}]\!] = (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{v}) [\![\boldsymbol{B}]\!], \end{gather}

(1.9a,b)\begin{gather} \boldsymbol{n} \boldsymbol{\cdot} [\![\boldsymbol{B}]\!] = 0, \quad \boldsymbol{n} \times [\![\boldsymbol{B}]\!] = \mu_0 \boldsymbol{K}, \end{gather}

where $[\![\cdots ]\!]$ denotes the jump of a quantity across the plasma–vacuum boundary, $\boldsymbol {n}$ is the surface normal vector of the plasma volume and $\boldsymbol {K}$ is a surface current.

Linearising the set of equations above yields an eigenvalue problem that describes resistive linear MHD stability. For ideal MHD ($\eta \equiv 0$), this eigenvalue problem can be used to derive an energy principle, which is based on minimisation of the energy functional representing the potential energy of a perturbation (Bernstein et al. 1958). Note that, in principle, the definition of a potential energy is not possible for MHD, because the MHD force is not a conservative force and no force potential exists. However, because the restriction to linear perturbations causes the work integral to become path-independent, the MHD energy functional only depends on the displacement. For this reason, the energy functional is named the ‘potential energy’ of the perturbation. As the work integral is also path-independent for resistive MHD and only depends on the spatial and magnetic perturbation, we keep the term potential energy for the energy functional. The energy functional of ideal MHD was transformed to an intuitive form (Greene & Johnson 1968), where the different terms of the functional separate pressure-gradient from current-density drive and stabilising phenomena. There are also energy principles for resistive MHD (Tasso 1977; Tasso & Virtamo 1980) which are meant to determine stability criteria or evaluate stability but do not provide information on the different contributions to the perturbed energy.

While the MHD stability of fusion plasmas is often well-described by ideal MHD, finite resistivity affects the stability of some equilibria significantly. In addition, there are types of instability that only exist in plasmas with finite resistivity. In order to describe the energy contributions to instabilities in a resistive plasma, the energy functional by Greene & Johnson (1968) must be extended to the regime of finite resistivity.

2. Resistive energy functional

Linearising the MHD equations and inserting (1.4) in (1.2) and (1.3) in (1.5) we obtain

(2.1)\begin{gather} \rho_0 \lambda^{2} \boldsymbol{\xi} = \boldsymbol{\nabla} (\boldsymbol{\xi} \boldsymbol{\cdot} \boldsymbol{\nabla} p_0) + \boldsymbol{\nabla} (p_0 \varGamma (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\xi})) + \frac{1}{\mu_0} (\boldsymbol{\nabla} \times \boldsymbol{B}_0) \times \boldsymbol{B}_1 + \frac{1}{\mu_0} (\boldsymbol{\nabla} \times \boldsymbol{B}_1) \times \boldsymbol{B}_0, \end{gather}

(2.2)\begin{gather} \boldsymbol{B}_1 = \boldsymbol{\nabla} \times (\boldsymbol{\xi} \times \boldsymbol{B}_0 ) - \lambda^{{-}1} \boldsymbol{\nabla} \times \left(\frac{\eta}{\mu_0} \boldsymbol{\nabla} \times \boldsymbol{B}_1\right), \end{gather}

where we denote equilibrium quantities by an index 0, perturbed quantities by an index 1 and assumed $f_1(x, t)=f_1(x)\text {e}^{\lambda t}$ for all perturbed quantities. Furthermore, we introduced the eigenvalue $\lambda =\gamma + {\rm i} \omega \in \mathbb {C}$, with growth rate $\gamma$ and oscillation frequency $\omega$, and the plasma displacement $\boldsymbol {\xi }$ defined by $\partial _t \boldsymbol {\xi }=\boldsymbol {v}_1$.

In contrast to ideal MHD, (2.1) can no longer be written purely in terms of $\boldsymbol {\xi }$ by inserting (2.2) into (2.1). The goal of the following derivation is to obtain a set of equations describing the different energetic drives of an eigenfunction similar to the ideal MHD version (Greene & Johnson 1968). The resulting energy functional is not in the form of a solvable eigenvalue problem but is intended to analyse the energetic composition of eigenfunctions that were calculated in advance by solving the MHD equations. In order to get an expression for the plasma energy $\delta W$, we multiply the displacement $\boldsymbol {\xi }^{\ast }/2$ of the perturbation with the restoring force defined by (2.1), which results in an energy density, and integrate the product over the plasma volume $V$. This is the common approach to derive the general complex energy functional $\delta W$. It should be noted that the energy densities related to the complex energy functional are not the physical energy densities of the real perturbation $\boldsymbol {\xi }_r=\mathrm {Re}(\boldsymbol {\xi })$. However, the physical energies and energy densities of $\delta W$ are easily obtained from the complex energy functional by taking the real part of all complex quantities and all complex-conjugate quantities, separately. This is because the physical work done by the real displacement $\boldsymbol {\xi }_r$ is obtained by its multiplication with the physical restoring force $\mathrm {Re}(\lambda ^{2}\boldsymbol {\xi })$. An example for deducing the physical energy density of the complex functional $\delta W$ is given at the end of this section.

To form the resistive energy functional, we apply cyclic rotation and integration by parts to the last term of (2.1), resulting in

(2.3)\begin{equation} F \equiv \int (\boldsymbol{\nabla} \times \boldsymbol{B}_1) \times \boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{\xi}^{{\ast}} \,\text{d}V ={-} \int \boldsymbol{\nabla} \times (\boldsymbol{\xi}^{{\ast}} \times \boldsymbol{B}_0) \boldsymbol{\cdot} \boldsymbol{B}_1 \,\text{d}V - \oint \boldsymbol{B}_1 \times (\boldsymbol{\xi}^{{\ast}} \times \boldsymbol{B}_0) \boldsymbol{\cdot} \boldsymbol{n} \,\text{d}S. \end{equation}

Inserting 2.2, integrating by parts and using some vector algebra results in

(2.4)\begin{align} - \frac{F}{\mu_0} & = \int \frac{1}{\mu_0} \| \boldsymbol{B}_1 \|^{2} + \frac{\eta}{\mu_0^{2}\lambda^{{\ast}}} \| \boldsymbol{\nabla} \times \boldsymbol{B}_1 \|^{2} \,\text{d}V \nonumber\\ & \quad + \oint \frac{\eta}{\mu_0^{2}\lambda^{{\ast}}} (\boldsymbol{\nabla} \times \boldsymbol{B}_1^{{\ast}}) \boldsymbol{\cdot} (\boldsymbol{B}_1 \times \boldsymbol{n}) + \frac{1}{\mu_0} (\boldsymbol{B}_1 \boldsymbol{\cdot} \boldsymbol{B}_0)(\boldsymbol{\xi}^{{\ast}} \boldsymbol{\cdot} \boldsymbol{n}) \,\text{d}S. \end{align}

Finally, integration by parts of the pressure-dependent terms in (2.1) yields an expression for the energy functional similar to Bernstein et al. (1958):

(2.5)\begin{equation} - \int \rho_0 \lambda^{2} \|\boldsymbol{\xi}\|^{2} \,\text{d}V = 2 \delta W = 2 \delta W_V + 2 \delta W_S, \end{equation}

where

(2.6)\begin{gather} 2 \delta W_V = W_\text{MAG} + W_\text{JXB} + W_\text{RCD} + W_\text{SND} + W_{\text{DP},C}, \end{gather}

(2.7)\begin{gather}2 \delta W_S = W_\text{MAG}^{s} + W_\text{RES}^{s} + W_\text{PRE}^{s} + W_\text{DP}^{s} \end{gather}

with the volume energy contributions $\delta W_V$

(2.8a,b)\begin{gather} W_\text{MAG}=\frac{1}{\mu_0} \int \| \boldsymbol{B}_1 \|^{2} \,\text{d}V, \quad W_\text{JXB}={-}\frac{1}{\mu_0}\int (\boldsymbol{\nabla} \times \boldsymbol{B}_0)\times \boldsymbol{B}_1 \boldsymbol{\cdot} \boldsymbol{\xi}^{{\ast}} \,\text{d}V, \end{gather}

(2.9a,b)\begin{gather}W_\text{RCD}=\frac{1}{\mu_0^{2}\lambda^{{\ast}}} \int \eta \| \boldsymbol{\nabla} \times \boldsymbol{B}_1 \|^{2} \,\text{d}V ,\quad W_\text{SND}=\varGamma \int p_0 \vert\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\xi}\vert^{2} \,\text{d}V, \end{gather}

(2.10)\begin{gather}W_{\text{DP},C}=\int (\boldsymbol{\xi} \boldsymbol{\cdot} \boldsymbol{\nabla} p_0)(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\xi}^{{\ast}}) \,\text{d}V \end{gather}

and surface energy contributions $\delta W_S$

(2.11a,b)\begin{gather} W_\text{MAG}^{s}=\frac{1}{\mu_0} \oint (\boldsymbol{B}_1 \boldsymbol{\cdot} \boldsymbol{B}_0)(\boldsymbol{\xi}^{{\ast}} \boldsymbol{\cdot} \boldsymbol{n}) \,\text{d}S, \quad W_\text{RES}^{s}=\frac{1}{\mu_0^{2}\lambda^{{\ast}}}\oint \eta (\boldsymbol{\nabla} \times \boldsymbol{B}_1^{{\ast}}) \boldsymbol{\cdot} (\boldsymbol{B}_1 \times \boldsymbol{n}) \,\text{d}S, \end{gather}

(2.12a,b)\begin{gather}W_\text{PRE}^{s}={-}\varGamma \oint p_0 (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\xi})(\boldsymbol{\xi}^{{\ast}} \boldsymbol{\cdot} \boldsymbol{n}) \,\text{d}S, \quad W_\text{DP}^{s}={-}\oint (\boldsymbol{\nabla} p_0 \boldsymbol{\cdot} \boldsymbol{\xi})(\boldsymbol{\xi}^{{\ast}}\boldsymbol{\cdot} \boldsymbol{n}) \,\text{d}S. \end{gather}

Two resistive energy contributions appear in (2.5), the resistive current diffusion $W_\text {RCD}$, and a surface term $W_\text {RES}^{s}$, which is a correction to $W_\text {MAG}^{s}$. Combining $W_\text {MAG}^{s}$ and $W_\text {RES}^{s}$, we obtain

(2.13)\begin{equation} W_\text{EXB}^{s} = W_\text{MAG}^{s}+W_\text{RES}^{s} = \frac{1}{\mu_0\lambda^{{\ast}}} \oint \boldsymbol{E}_1^{{\ast}} \times \boldsymbol{B}_1 \boldsymbol{\cdot} \boldsymbol{n}\,\text{d}S, \end{equation}

which describes a flow of electromagnetic energy through the plasma surface, the Poynting flux. The ideal volume contributions in (2.5) are related to the magnetic energy $W_\text {MAG}$, the Lorentz force $W_\text {JXB}$, adiabatic compression $W_\text {SND}$ and compression against the pressure-gradient $W_{\text {DP},C}$. The surface contributions always relate to a jump of a quantity at the plasma boundary. In the following, we bring the terms $W_\text {MAG}$, $W_\text {JXB}$ and $W_{\text {DP},C}$ into their intuitive form, which separates current-density from pressure-gradient drive, analogously to Greene & Johnson (1968). Splitting the equilibrium current density into a parallel $\boldsymbol {j}_{0,\|}$ and perpendicular $\boldsymbol {j}_{0,\perp }$ component with respect to the magnetic field $\boldsymbol {B}_0$ results in

(2.14)\begin{gather} - W_\text{JXB} = \int \frac{j_{0,\|}}{\| \boldsymbol{B}_0 \|} (\boldsymbol{\xi}^{{\ast}} \times \boldsymbol{B}_0) \boldsymbol{\cdot} \boldsymbol{B}_1 \,\text{d}V + \int \frac{1}{\boldsymbol{B}_0^{2}}(\boldsymbol{B}_1 \boldsymbol{\cdot} \boldsymbol{B}_0)(\boldsymbol{\nabla} p_0 \boldsymbol{\cdot} \boldsymbol{\xi}^{{\ast}}) \,\text{d}V \nonumber\\ - \int \frac{1}{\boldsymbol{B}_0^{2}}(\boldsymbol{B}_1 \boldsymbol{\cdot} \boldsymbol{\nabla} p_0)(\boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{\xi}^{{\ast}}) \,\text{d}V \end{gather}

with $j_{0,\|}=(\boldsymbol {j}_{0}\boldsymbol {\cdot } \boldsymbol {B}_0)/\|\boldsymbol {B}_0\|$. Analogously, we split the divergence $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\xi }$ in $W_{\text {DP},C}$ into $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\xi }_\|$ and $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\xi }_\perp$. After some vector algebra and integration by parts, we obtain

(2.15)\begin{gather} \boldsymbol{B}_0^{2} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\xi}_\perp{=} \mu_0 \, \boldsymbol{\xi}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} p_0 - 2 \boldsymbol{B}_0^{2} (\boldsymbol{\xi}_\perp \boldsymbol{\cdot} \boldsymbol{\kappa}) - \boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{B}_1 - \boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{B}_R, \end{gather}

(2.16)\begin{gather} - \int (\boldsymbol{\xi} \boldsymbol{\cdot} \boldsymbol{\nabla} p_0)(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\xi}_\|^{{\ast}}) \,\text{d}V = \int \frac{1}{\boldsymbol{B}_0^{2}}(\boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{\xi}^{{\ast}})(\boldsymbol{B}_1\boldsymbol{\cdot}\boldsymbol{\nabla} p_0) + \frac{1}{\boldsymbol{B}_0^{2}}(\boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{\xi}^{{\ast}})(\boldsymbol{B}_R\boldsymbol{\cdot}\boldsymbol{\nabla} p_0) \,\text{d}V, \end{gather}

where we introduce the self-induced ohmic field $\boldsymbol {B}_R=\lambda ^{-1}\boldsymbol {\nabla }\times (\eta \boldsymbol {j}_1)$ and the equilibrium curvature vector $\boldsymbol {\kappa }=(\boldsymbol {b} \boldsymbol {\cdot } \boldsymbol {\nabla }) \boldsymbol {b}$ with $\boldsymbol {b}=\boldsymbol {B}_0/\|\boldsymbol {B}_0\|$. Note that for constant resistivity, we get $\boldsymbol {B}_R=({\eta }/{\mu _0\lambda })\Delta \boldsymbol {B}_1$ which describes resistive diffusion. Equations (2.15) and (2.16) relate the (work done by) parallel and perpendicular divergence of the plasma displacement to the (work done by the) magnetic perturbation using (2.2). In contrast to ideal MHD, the perturbed field generated by the displacement is changed by its self-induced ohmic field $\boldsymbol {B}_R$, which is taken into account by the corrective resistive terms in (2.15) and (2.16). Finally, combining (2.14) to (2.16) yields

(2.17)\begin{align} & W_\text{MAG} + W_\text{JXB} + W_{\text{DP},C} = \int \frac{1}{\mu_0} \left\| \boldsymbol{B}_1 - \boldsymbol{B}_0 \frac{\mu_0 \, \boldsymbol{\xi}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} p_0}{\boldsymbol{B}_0^{2}} \right\|^{2} \nonumber\\ & \qquad - \frac{j_{0,\|}}{\| \boldsymbol{B}_0 \|} (\boldsymbol{\xi}^{{\ast}}_\perp{\times} \boldsymbol{B}_0) \boldsymbol{\cdot} \boldsymbol{B}_1 - 2 (\boldsymbol{\xi}^{{\ast}}_\perp \boldsymbol{\cdot} \boldsymbol{\kappa}) (\boldsymbol{\xi}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} p_0) \nonumber\\ & \qquad - \frac{1}{\boldsymbol{B}_0^{2}} \left[ (\boldsymbol{\xi}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} p_0)(\boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{B}_R^{{\ast}}) + (\boldsymbol{\xi}_\|^{{\ast}} \boldsymbol{\cdot} \boldsymbol{B}_0)(\boldsymbol{\nabla} p_0 \boldsymbol{\cdot} \boldsymbol{B}_R) \right] \,\text{d}V. \end{align}

Then, using (2.15) and (2.17), the intuitive form of the energy functional for resistive MHD reads

(2.18)\begin{align} 2 \delta W_V & =\underbrace{\int \vphantom{\frac{1}{\boldsymbol{B}_0^{2}}} \frac{\| \boldsymbol{B}_{1,\perp} \|^{2}}{\mu_0} \,\text{d}V}_{W_\text{SHA}} + \underbrace{\int \frac{\boldsymbol{B}_0^{2}}{\mu_0} \left\vert \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\xi}_\perp{+} 2 (\boldsymbol{\xi}_\perp \boldsymbol{\cdot} \boldsymbol{\kappa}) + \frac{1}{\boldsymbol{B}_0^{2}}\boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{B}_R \right\vert^{2} \,\text{d}V}_{W_\text{CPA}}\nonumber\\ & \quad \underbrace{+ \varGamma \int p_0 \vert\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\xi}\vert^{2} \,\text{d}V}_{W_\text{SND}} \quad \underbrace{+ \frac{1}{\mu_0^{2} \lambda^{{\ast}}} \int \eta \| \boldsymbol{\nabla} \times \boldsymbol{B}_1 \|^{2} \,\text{d}V}_{W_\text{RCD}}\nonumber\\ & \quad \underbrace{- \int \frac{j_{0,\|}}{\| \boldsymbol{B}_0 \|} (\boldsymbol{\xi}_\perp^{{\ast}} \times \boldsymbol{B}_0) \boldsymbol{\cdot} \boldsymbol{B}_{1,\perp} \,\text{d}V}_{W_\text{CUR}} \quad \underbrace{- \int 2 (\boldsymbol{\xi}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} p_0) (\boldsymbol{\xi}_\perp^{{\ast}} \boldsymbol{\cdot} \boldsymbol{\kappa}) \,\text{d}V}_{W_\text{DP0}}\nonumber\\ & \quad \underbrace{- \int \frac{1}{\boldsymbol{B}_0^{2}} (\boldsymbol{\xi}_\|^{{\ast}} \times \boldsymbol{\nabla} p_0)\boldsymbol{\cdot}(\boldsymbol{B}_0 \times \boldsymbol{B}_R) \,\text{d}V}_{W_{\text{RD,}\|}} \quad \underbrace{- \int \frac{1}{\boldsymbol{B}_0^{2}} (\boldsymbol{\xi}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} p_0)(\boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{B}_R^{{\ast}})\,\text{d}V}_{W_{\text{RD,}\perp}}, \end{align}

where the first four terms, which are generally stabilising, relate to perpendicular magnetic perturbations (shear Alfvén waves) $W_\text {SHA}$, parallel magneto-compressional perturbations (compressional Alfvén waves) $W_\text {CPA}$, adiabatic compression (sound waves) $W_\text {SND}$ and resistive current diffusion $W_\text {RCD}$, respectively. The next two terms describe the parallel current-density $W_\text {CUR}$ and field-line-curvature-dependent pressure-gradient $W_\text {DP0}$ drives, respectively. The last two terms are resistive corrections to the pressure-gradient drive. They account for the effect of resistive diffusion on the work done by the induced magnetic perturbation that is generated by parallel $W_{\text {RD,}\|}$ and perpendicular $W_{\text {RD,}\perp }$ compression of the plasma, respectively. Finally, the resistive correction to the relation between perpendicular compression and its induced magnetic perturbation also appears in the stabilising energy of magneto-compressional perturbations $W_\text {CPA}$. In summary, taking finite resistivity into account, there are three new volume contributions to the energy functional ($W_\text {RCD}$, $W_{\text {RD,}\|}$, $W_{\text {RD,}\perp }$) as well as a correction to the stabilising energy of magneto-compressional perturbations ($W_\text {CPA}$). Note that resistive diffusion can affect the perturbation over a longer time-scale if the mode grows or oscillates slower, which is described by the inverse proportionality of the resistive energy terms with the eigenvalue of the perturbation.

In the following, we extend the energy functional to the vacuum region, analogously to Bernstein et al. (1958). The linearised boundary condition for the pressure balance (1.7) is not dependent on resistivity and remains unchanged:

(2.19)\begin{equation} \boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{B}_1 - \mu_0 \varGamma p_0 (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\xi}) = \frac{1}{2} \boldsymbol{\xi} \boldsymbol{\cdot} (\boldsymbol{\nabla} \hat{\boldsymbol{B}}_0^{2} - \boldsymbol{\nabla} \boldsymbol{B}_0^{2}) + \hat{\boldsymbol{B}}_0 \boldsymbol{\cdot} \hat{\boldsymbol{B}}_1, \end{equation}

where we denote the vacuum quantities by a hat. An additional boundary condition is obtained by combining Ohm's law (1.3) with the boundary condition for the tangential electric field (1.8):

(2.20)\begin{equation} {-}\boldsymbol{n} \times \hat{\boldsymbol{A}}_1 = \frac{\eta}{\lambda} \boldsymbol{n} \times \boldsymbol{j}_1 + \hat{\boldsymbol{B}}_0 (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\xi}) \quad\text{with }\ \hat{\boldsymbol{A}}_1={-}\frac{1}{\lambda}\hat{\boldsymbol{E}}_1 \end{equation}

resulting in a resistive correction on the right-hand side. Applying (2.19) to $W_\text {MAG}^{s}$ and $W_\text {PRE}^{s}$ results in

(2.21)\begin{equation} W_\text{MAG}^{s}+W_\text{PRE}^{s} = \oint \frac{1}{2\mu_0} (\boldsymbol{\xi}^{{\ast}} \boldsymbol{\cdot} \boldsymbol{n}) \boldsymbol{\xi} \boldsymbol{\cdot} [\![\boldsymbol{\nabla} \boldsymbol{B}_0^{2}]\!] \,\text{d}S + \oint \frac{1}{\mu_0} (\boldsymbol{\xi}^{{\ast}} \boldsymbol{\cdot} \boldsymbol{n})(\hat{\boldsymbol{B}}_1 \boldsymbol{\cdot} \hat{\boldsymbol{B}}_0) \,\text{d}S. \end{equation}

Adding the latter surface integral in (2.21) to $W_\text {RES}^{s}$ and using (2.20) yields

(2.22)\begin{align} \oint \frac{1}{\mu_0} (\boldsymbol{\xi}^{{\ast}} \boldsymbol{\cdot} \boldsymbol{n})(\hat{\boldsymbol{B}}_1 \boldsymbol{\cdot} \hat{\boldsymbol{B}}_0) \,\text{d}S + \oint \frac{\eta}{\mu_0 \lambda^{{\ast}}} \boldsymbol{j}_1^{{\ast}} \boldsymbol{\cdot} (\boldsymbol{B}_1 \times \boldsymbol{n}) \,\text{d}S = \int \frac{ \| \hat{\boldsymbol{B}}_1 \|^{2} }{\mu_0}\,\text{d}\hat{{V}} + \oint \frac{\eta}{\lambda^{{\ast}}} \boldsymbol{j}_1^{{\ast}} \boldsymbol{\cdot} \boldsymbol{K}_1 \,\text{d}S, \end{align}

where we introduced the perturbed surface current $\boldsymbol {K}_1=\mu _0^{-1} \boldsymbol {n} \times [\![\boldsymbol {B}_1]\!]$. However, surface currents cannot be maintained in resistive plasmas, i.e. $\eta \neq 0$ implies ${\boldsymbol {K}_1=0}$. This means that all components of the magnetic field must be continuous in resistive plasmas. Finally, combining (2.21) and (2.22) with $W_\text {DP}^{s}$, we obtain

(2.23)\begin{equation} \delta W_S=\delta W_{S'}+\delta W_{\hat{V}} \end{equation}

with

(2.24a,b)\begin{equation} 2\delta W_{S'}=\oint\vert \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\xi}\vert^{2} \boldsymbol{n}\boldsymbol{\cdot} \left[\!\!\left[\boldsymbol{\nabla} \left( p_0 + \frac{\boldsymbol{B}_0^{2}}{2\mu_0} \right)\right]\!\!\right] \,\text{d}S , \quad 2\delta W_{\hat{V}}=\int \frac{\| \hat{\boldsymbol{B}}_1 \|^{2}}{\mu_0} \,\text{d}\hat{{V}}. \end{equation}

As a result, no additional surface or vacuum energy contributions have to be considered in the case of finite resistivity. Note that $\delta W_{S'}=0$ if there are no equilibrium surface currents, which is a common constraint on equilibria and is consistent with the resistive boundary conditions. However, in principle, one could allow for finite surface currents in an ideal equilibrium and test this equilibrium for resistive instabilities. In the case of vanishing surface current, the surface terms in (2.5) are equivalent to the vacuum energy of the perturbation, i.e. $\delta W_S=\delta W_{\hat {V}} > 0$. While the vacuum contribution $\delta W_{\hat {V}}$ provides a simpler, more accessible understanding of the perturbation energy outside the plasma, the surface contributions $\delta W_S$ encode the vacuum energy in a single surface, which can be computationally advantageous. If there is a resistive wall or resistive components in the vacuum region, additional terms in (2.23) are required to take the interaction of the perturbation with the wall into account. These additional terms can be easily obtained from the divergence theorem applied in (2.22).

Finally, we demonstrate how to obtain the physical energies and energy densities from the complex terms of the resistive energy functional derived previously. As discussed in the beginning of this section, one must take the real part of the complex and complex-conjugate quantities in order to obtain the physical energies. For the energy density related to resistive current diffusion $w_\text {RCD}$, we obtain

(2.25)\begin{gather} w_\text{RCD} = \frac{1}{\mu_0^{2}\lambda^{{\ast}}} \eta \| \boldsymbol{\nabla} \times \boldsymbol{B}_1 \|^{2} \ \rightarrow\ \mathcal{P}\left\{w_\text{RCD}\right\} = \frac{\eta}{\mu_0^{2}}\ \mathrm{Re} \left(\lambda^{{-}1} \boldsymbol{\nabla} \times \boldsymbol{B}_1 \right) \boldsymbol{\cdot} \mathrm{Re} \left(\boldsymbol{\nabla} \times \boldsymbol{B}_1\right)\nonumber\\ \stackrel{\text{if }\lambda \in \mathbb{R}}{=} \frac{\eta}{\mu_0^{2}\lambda} \mathrm{Re} \left(\boldsymbol{\nabla} \times \boldsymbol{B}_1 \right)^{2}, \end{gather}

where $\mathcal {P}\{\cdots \}$ denotes the physical energy or energy density and the lowercase $w_{\Box }$ denotes the energy density or integrand of the energy term $W_{\Box }$. Note that one must express the absolute values in terms of complex and complex-conjugate quantities before extracting the physical energy.

The real part of the energy functional and the physical energy functional are related to the kinetic energy by

(2.26a,b)\begin{equation} \delta E_{\text{kin}} + \delta E_{\text{osc}} + \mathrm{Re} \left(\delta W\right) = 0,\quad \mathcal{P}\left\{\delta E_{\text{kin}}\right\} + \mathcal{P}\left\{\delta E_{\text{osc}}\right\} + \mathcal{P}\left\{\delta W\right\} = 0, \end{equation}

where the kinetic energy is $\delta E_{\text {kin}}=\int \tfrac {1}{2} \rho _0 \|\boldsymbol {v}_1\|^{2} \,\text {d}V$ and the energy stored in the oscillation frequency is $\delta E_{\text {osc}}=-\int \omega ^{2} \|\boldsymbol {\xi }\|^{2} \,\text {d}V$; for the physical energies, we have $\mathcal {P}\{\delta E_{\text {kin}}\}=\int \tfrac {1}{2} \rho _0\ \mathrm {Re} (\boldsymbol {v}_1)^{2} \,\text {d}V$ and $\mathcal {P}\{\delta E_{\text {osc}}\}=\tfrac {1}{2} \delta E_{\text {osc}}$. The oscillation energy is negligible if $\gamma \gg \omega$.

3. Numerical validation

The intuitive form of the energy functional for finite resistivity has been implemented as a diagnostic for eigenfunctions in the linear MHD stability code CASTOR3D (Strumberger & Günter 2017, 2019), which is formulated in general 3D curvilinear coordinates. We validate the energy functional by comparison of potential and kinetic energy and of the energy contributions of $\delta W$ over three different coordinate systems: NEMEC coordinates (NEM), 2D straight-field-line coordinates (SFL) and Boozer coordinates (BZR) (Strumberger & Günter 2017). The following validation is based on an $n=8$ edge-localised mode of a simple low-$\beta$ axisymmetric equilibrium surrounded by an infinite vacuum region. The equilibrium was calculated using the NEMEC equilibrium code (Hirshman & Whitson 1983; Hirshman, van RIJ & Merkel 1986).

Figure 1 shows the flux surfaces and lines of constant poloidal coordinate $u$ of the validation case for a toroidal cross-section, i.e. constant toroidal coordinate $v$, in the different coordinate systems. One can see that the lines of constant poloidal coordinates are equal for SFL and BZR coordinates. However, the toroidal cross-section in SFL coordinates is planar in cylindrical coordinates, while in BZR coordinates it is deformed and extends over multiple toroidal angles $\varPhi$. The equilibrium pressure and safety factor profiles (see Puchmayr, Dunne & Strumberger 2022) of the validation case are displayed in figure 2(a). For the resistivity, a polynomial profile that is small in the plasma-core and large at the edge is used:

(3.1)\begin{equation} \eta(\rho_\text{tor})=\eta_0 \times (1+10 \rho_\text{tor}+10 \rho_\text{tor}^{2}+100 \rho_\text{tor}^{8}+100 \rho_\text{tor}^{10}), \end{equation}

where $\rho _\text {tor}$ is the square-root of the normalised toroidal flux. The shape of this resistivity profile is shown in figure 2(b). In the following, all energies and energy densities are normalised with respect to $\int \tau _A^{-2} \rho \| \boldsymbol {\xi } \|^{2} \,\text {d}V$ where $\tau _A={\sqrt {\mu _0\rho _\text {ax}}R_\text {ax}}/{B_\text {ax}}$ is the Alfvén time, $R_\text {ax}$ is the major radius and the subscript ‘ax’ denotes equilibrium quantities at the magnetic axis. For the physical energies and energy densities, the normalisation factor naturally becomes $\int \tau _A^{-2} \rho \,\mathrm {Re}(\boldsymbol {\xi })^{2} \,\text {d}V$. The normalisation of the energy terms is necessary in order to achieve comparability of the energy for different eigenfunctions, since the amplitude of the perturbation is arbitrary in linear MHD.

Figure 1. Toroidal cross-section with flux surfaces (black) and lines of constant poloidal coordinate for NEM (dark blue), SFL (purple) and BZR (green, dashed) coordinates; toroidal angle of the cross-section for BZR coordinates (background colour, blue to red).

Figure 2. (a) Equilibrium pressure profile (solid) and safety factor profile (dashed). (b) Profile shape of the resistivity as defined in (3.1) normalised to $\eta (1)=1$.

Table 1 lists the potential energy contributions of the $n=8$ eigenfunction for a resistivity of $\eta _0=10^{-8}\, \Omega \text {m}$. The Fourier spectra in BZR coordinates of the radial and toroidal velocity perturbation for this eigenfunction are displayed in figure 3. In order to resolve the eigenfunction, different sets of poloidal mode numbers are required for each coordinate system. From the comparison of the last two columns of table 1, one can see that the kinetic energy, i.e. the eigenvalue, eventually converges much faster with respect to the number of poloidal modes than some of the energy terms, especially the resistive current diffusion term $W_\text {RCD}$ which depends on ${\boldsymbol {\nabla } \times \boldsymbol {B}_1}$. This leads to the large difference of kinetic and potential energy for the BZR coordinates with the small set of poloidal mode numbers $m=20\unicode{x2013}70$. While for BZR coordinates, the spectrum $m=20\unicode{x2013}70$ is probably sufficient to solve the numerical eigenvalue problem for the growth rate of the mode, some derivatives of the eigenfunction can still be poorly resolved. Figure 4(a) shows the convergence of the radially resolved energy density $w_\text {RCD}$ for different sets of poloidal mode numbers. The radial resolution of the eigenfunction must also be large enough to resolve the behaviour of the potential energy densities at the resonant flux surfaces. The radially resolved energy densities can vary strongly near these resonant surfaces, as can be seen in figure 4(b).

Figure 3. Fourier spectra ($m=20\unicode{x2013}90$) of the (a) radial and (b) toroidal velocity perturbation $\boldsymbol {v}_1$ in BZR coordinates for the $n=8$ eigenfunction and a resistivity of $\eta _0=10^{-8}\,\Omega \text {m}$.

Figure 4. (a) Normalised energy density $w_\text {RCD}$ for poloidal mode numbers $m=20\unicode{x2013}70$ (red), $m=20\unicode{x2013}75$ (purple), $m=20\unicode{x2013}80$ (blue) and $m=20\unicode{x2013}90$ (black) in BZR coordinates. (b) Normalised kinetic energy density (dark blue) and normalised potential energy densities $w_\text {RCD}$ (red) and $w_\text {SHA}$ (dark red) in BZR coordinates and resonant surfaces (light grey, dashed).

Table 1. Normalised potential energy contributions calculated in the different coordinate systems and numerical deviation of the potential and kinetic energies.

As one can see in table 1, the kinetic and potential energy as well as the different energy contributions match well for all coordinate systems provided that sufficient poloidal resolution has been chosen. The numerical equality of potential and kinetic energy validates that all energy contributions have been considered, while the numerical equality between the coordinate systems validates the implementation in general curvilinear 3D coordinates. Note that in order to get such a high accuracy of the energy functional over different coordinate systems requires not only a sufficient radial and poloidal resolution of the eigenfunction but also a highly resolved force equilibrium.

4. Influence of finite resistivity

In the following, we investigate the influence of increasing resistivity on the potential energy composition of the eigenfunction for the validation case. In order to estimate the effect of the resistive corrections with increasing resistivity, we define the following potential energy proportions:

(4.1a,b)\begin{gather} \chi_\text{RES}^{s}=\frac{\left|W_\text{RES}^{s}\right|}{\left|W_\text{MAG}^{s}\right|+\left|W_\text{RES}^{s}\right|} ,\quad \chi_\text{RCD}=\frac{W_\text{RCD}}{W_\text{RCD}+W_\text{SHA}+W_\text{CPA}+W_\text{SND}}, \end{gather}

(4.2)\begin{gather} \chi_\text{RD}=\frac{W_{\text{RD,}\perp}+W_{\text{RD,}\|}}{W_{\text{RD,}\perp}+W_{\text{RD,}\|}+W_\text{DP0}}, \end{gather}

with the proportion of the resistive correction relative to the surface or vacuum energy $\chi _\text {RES}^{s}$, to the stabilising energy terms $\chi _\text {RCD}$ and to the pressure-gradient drive $\chi _\text {RD}$. We also define the ratio of pressure-gradient drive to the destabilising terms (current-density and pressure-gradient drive):

(4.3)\begin{equation} \chi_\text{DP/CUR}=\frac{W_{\text{RD,}\perp}+W_{\text{RD,}\|}+W_\text{DP0}}{W_{\text{RD,}\perp}+W_{\text{RD,}\|}+W_\text{DP0}+W_\text{CUR}} \end{equation}

The newly defined energy proportions are shown in figure 5(a). One can see that the surface correction $W_\text {RES}^{s}$ (red) and the resistive current diffusion $W_\text {RCD}$ (purple) become important even for small values of resistivity, while the relative correction to the pressure-gradient drive $W_{\text {RD,}\perp }+W_{\text {RD,}\|}$ (blue) only increases slowly with increasing resistivity. Note that the validation case is ideally unstable and the effect of the correction $W_{\text {RD,}\perp }+W_{\text {RD,}\|}$ might be larger for purely resistive modes. Moreover, the relative effect of the corrections on the growth rate is often larger compared with the effect of the corrections on the related energy terms because the growth rate is given by the sum of the energy terms which might be smaller than its summands. Finally, from the ratio of pressure-gradient to current-density drive $\chi _\text {DP/CUR}$ (black curve in figure 5a) it can be seen that the eigenfunction becomes increasingly pressure-gradient driven for increasing resistivity.

Figure 5. (a) Proportions of the resistive corrections with respect to the energy terms $\chi _\text {RES}^{s}$ (red), $\chi _\text {RCD}$ (purple) and $\chi _\text {RD}$ (blue) as well as the ratio of pressure-gradient to current-density drive $\chi _\text {DP/CUR}$ (black); the grey area indicates the range of realistic resistivity values. (b) Magnetic energy density of the perpendicular perturbation $w_\text {SHA}$ for resistivity values of $\eta _0=0\, \Omega \text {m}$ (black), $\eta _0=10^{-8}\, \Omega \text {m}$ (blue), $\eta _0=10^{-7}\, \Omega \text {m}$ (purple) and $\eta _0=10^{-6}\, \Omega \text {m}$ (red).

Figure 5(b) shows the radially resolved perpendicular magnetic energy density $w_\text {SHA}$ for different values of resistivity. The radial structure of the magnetic perturbation becomes less sensitive to resonant surfaces and the maximum of the magnetic perturbation shifts inwards with increasing resistivity. In addition, the poloidal localisation of the potential energy contributions is affected by resistivity as can be seen in figure 6 which shows the energy density corresponding to the pressure-gradient drive. With increasing resistivity, the pressure-gradient drive of the eigenfunction becomes at first slightly less localised at the low-field side compared with the ideal mode, which can be seen from the relative increase of the energy density in the upper, lower and high-field side regions of the plasma (see, e.g., enlarged regions of figure 6a,b). Further increasing the resistivity causes the energy density of the pressure-gradient drive to move to the upper and lower regions of the plasma, away from the midplane. Furthermore, it can be seen that the eigenfunction is weakly stabilised (magenta regions in figure 6) by the pressure-gradient on the high-field side and at some locations on the low-field side. We note here that such trends can also be seen for resistive modes at realistic resistivity values.

Figure 6. Energy density of the pressure-gradient drive including resistive corrections ($w_\text {DP0}+w_{\text {RD,}\perp }+w_{\text {RD,}\|}$) for resistivity values of ${\eta _0=0\, \Omega \text {m}}$ (a), ${\eta _0=10^{-8}\, \Omega \text {m}}$ (b) and ${\eta _0=10^{-6}\, \Omega \text {m}}$ (c). The colour range is normalised with respect to the maximum and minimum energy densities for each plot separately and the colour scale is linear in the intervals $[-1,0]$ and $[0,1]$. The minimum normalised energy densities are ${w_\text {DP}^{\min}=-1.132}$ (a), ${w_\text {DP}^{\min}=-1.029}$ (b) and ${w_\text {DP}^{\min}=-1.427}$ (c); the maximum normalised energy densities are ${w_\text {DP}^{\max}=0.033}$ (a), ${w_\text {DP}^{\max}=0.027}$ (b) and ${w_\text {DP}^{\max}=0.073}$ (c).

5. Conclusion

Expressions for the different stabilising and destabilising energy contributions of plasma perturbations in the framework of resistive MHD have been derived similar to the intuitive form of the ideal energy functional by Greene & Johnson (1968). By including finite resistivity, three additional resistive energy terms ($W_\text {RCD}$, $W_{\text {RD,}\|}$, $W_{\text {RD,}\perp }$) as well as a correction to the stabilising energy of magneto-compressional perturbations ($W_\text {CPA}$) appear in the energy contributions of the plasma volume (2.18) and have been discussed in § 4. We have shown that no corrections have to be taken into account for the energy contributions of the plasma surface and vacuum region (2.23). Furthermore, we have demonstrated how to obtain the physical energy densities from the complex energy functional by inserting real parts for the complex and complex-conjugate quantities, separately.

The terms of the resistive energy functional were implemented in the non-ideal linear stability code CASTOR3D (Strumberger & Günter 2017, 2019) in general curvilinear coordinates. The convergence of the energy terms has been analysed for a simple test equilibrium, showing excellent numerical equality of the potential and kinetic energy and between the different coordinate systems that are implemented in the CASTOR3D code. In order to obtain the numerical equality achieved for the test equilibrium, it is necessary to use a highly accurate equilibrium as well as a sufficiently high number of flux surfaces and poloidal/toroidal mode numbers. It has been shown that the resistive corrections can significantly affect the energy contributions of the ideal energy functional for realistic values of the resistivity. Finally, trends in the localisation of the energy density with increasing resistivity were presented, where, with increasing resistivity, the perturbation moves away from the midplane. In conclusion, comparison of the energy terms and the analysis of the locally resolved physical energy densities using the new energy functional for finite resistivity grants increased insight to the linear stability of resistive modes compared with trends in growth rates and perturbed quantities. The energetic decomposition contributes to the physical understanding of these perturbations and allows to distinguish current-density and pressure-gradient driven instabilities.