Legolas: magnetohydrodynamic spectroscopy with viscosity and Hall current

Abstract

Many linear stability aspects in plasmas are heavily influenced by non-ideal effects beyond the basic ideal magnetohydrodynamics (MHD) description. Here, the extension of the modern open-source MHD spectroscopy code Legolas with viscosity and the Hall current is highlighted and benchmarked on a stringent set of historic and recent findings. The viscosity extension is demonstrated in a cylindrical set-up featuring Taylor–Couette flow and in a viscoresistive plasma slab with a tearing mode. For the Hall extension, we show how the full eigenmode spectrum relates to the analytic dispersion relation in an infinite homogeneous medium. We quantify the Hall term influence on the resistive tearing mode in a Harris current sheet, including the effect of compressibility, which is absent in earlier studies. Furthermore, we illustrate how Legolas mimics the incompressible limit easily to compare with literature results. Going beyond published findings, we emphasise the importance of computing the full eigenmode spectrum, and how elements of the spectrum are modified by compressibility. These extensions allow for future stability studies with Legolas that are relevant to ongoing dynamo experiments, protoplanetary disks or magnetic reconnection.

1. Introduction

The open-source linear magnetohydrodynamics (MHD) spectroscopy code Legolas was first introduced in Claes, De Jonghe & Keppens (2020, see also https://legolas.science), and allows for computation of all linear eigenmodes and their eigenfunctions, for one-dimensionally (1-D) stratified plasmas in a wide range of settings. Such non-homogeneous 3-D plasma states with 1-D variation are very common, e.g. in plane-parallel, gravitationally stratified atmospheres, in cylindrical set-ups for Taylor–Couette experiments or in magnetic flux tubes or loops in the solar atmosphere. MHD spectroscopy is useful to determine the complete stability properties of a given force-balanced equilibrium state, and quantifies how this (in)stability is influenced by specific equilibrium ingredients, such as the magnetic pitch, or the presence of non-trivial background flows. In the original release (Claes et al. 2020), the linearised set of compressible MHD equations included the (possibly combined) effects of flow, external gravity, resistivity, anisotropic thermal conduction and radiative cooling. The code was tested and validated against a wide variety of theoretically known plasma stability results, e.g. those from modern plasma physics textbooks focusing on MHD spectroscopy (Goedbloed, Keppens & Poedts 2019). A typical application can be found in Claes & Keppens (2021), where it was shown that especially the effect of radiative cooling in a realistic, magnetised solar atmosphere model yields MHD eigenspectra dominated by thermal instabilities and magneto-thermal overstabilities.

Ideal MHD stability aspects have been studied extensively in the plasma physics literature, in various applications. For fusion devices such as tokamaks, a solid understanding of instabilities is required to create MHD-stable operation conditions, whilst in solar physics both stable wavemodes and instabilities are of interest to understand the observed periodicities, or the evolution of initially stable coronal loops towards destructive events such as coronal mass ejections. The non-adiabatic effects already included in Claes et al. (2020) allow for investigation of radiatively driven processes such as solar coronal rain or prominence formation. Non-ideal effects like resistivity are known to introduce new paths to instability: the well-known resistive tearing instability (Furth, Killeen & Rosenbluth 1963) has been the subject of many studies and received renewed interest due to the role it plays in triggering magnetic reconnection events. However, physical effects that were previously omitted in Legolas, in particular viscosity and the Hall current, may influence growth rates or modify stability properties in a significant way.

For resistive tearing and the resulting reconnection in particular, nonlinear simulations have shown that the rates at which magnetic field lines reconnect (and convert magnetic energy into kinetic or thermal energy in the process) can become much higher than resistivity can account for on its own (see e.g. the review by Yamada, Kulsrud & Ji 2010). Both viscosity and Hall effects are candidates for modifying the growth rate of the resistive tearing instability. In fact, it has been known for quite some time that viscosity can act as a stabilising mechanism (Coppi, Greene & Johnson 1966; Loureiro, Schekochihin & Uzdensky 2013; Tenerani et al. 2015) whilst the Hall current may introduce destabilising effects, resulting in faster reconnection rates (Terasawa 1983; Pucci, Velli & Tenerani 2017). Similarly, exploration of the relationship between resistivity and viscosity on tearing by Dahlburg et al. (1983) revealed more intricacies of the growth rate as a function of the resistivity and viscosity. More recent work has focused on evaluating the influence of the Hall effect on the tearing instability in current sheets (Shi et al. 2020). We here document the extension of Legolas with viscous and Hall terms, and benchmark our MHD spectroscopy tool on these and some other published results. At the same time, we show how we can easily extend published findings with full spectral knowledge, or with quantifications of how incompressible and compressible regimes differ.

Whilst there are many more applications of Hall-MHD, such as the modification of the magnetorotational instability (MRI) by the Hall current in protoplanetary accretion discs (see e.g. the lecture notes by Lesur 2021) or the influence of the Hall current on different instabilities in dynamo experiments (see e.g. Mishra, Mamatsashvili & Stefani 2022), these subjects are not treated in this current report, but are suitable to be investigated with Legolas in future work. In particular, Legolas's multirun framework is quite apt to perform parametric studies exploring stability in function of the magnetic Prandtl number in viscoresistive set-ups. Such studies could also benefit from the inclusion of ambipolar diffusion, which we intend to implement in a future release of the code.

To demonstrate Legolas's new capabilities, this paper is structured as follows. In § 2 we present a brief overview of the compressible MHD equations implemented in Legolas (contrary to the incompressible equations that are often used for the tearing instability). Section 3 contains the main results of this paper, with the presentation of new diagnostic tools in § 3.1, the viscosity module in § 3.2 and the Hall module in § 3.3. Finally, § 4 concludes the paper with a discussion of the preliminary results using these new Legolas extensions.

2. Problem description and model equations

As presented in Claes et al. (2020), but extended here with previously ignored physical effects, Legolas considers the full set of compressible MHD equations

(2.1)\begin{gather} \frac{\partial \rho}{\partial t} ={-}\boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{v}), \end{gather}

(2.2)\begin{gather}\rho\frac{\partial \boldsymbol{v}}{\partial t} ={-}\boldsymbol{\nabla} p - \rho \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v} + \boldsymbol{J} \times \boldsymbol{B} + \rho\boldsymbol{g} + \boldsymbol{F}_\mathrm{visc}, \end{gather}

(2.3)\begin{gather} \rho\frac{\partial T}{\partial t}={-}\rho \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} T - (\gamma - 1)p\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} - (\gamma - 1)\rho\mathscr{L} + (\gamma - 1)\boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{\kappa} \boldsymbol{\cdot} \boldsymbol{\nabla} T) \nonumber\\ + (\gamma - 1)\eta\boldsymbol{J}^2 + (\gamma-1) H_\mathrm{visc}, \end{gather}

(2.4)\begin{gather}\frac{\partial \boldsymbol{B}}{\partial t} = \boldsymbol{\nabla} \times (\boldsymbol{v} \times \boldsymbol{B}) - \boldsymbol{\nabla} \times (\eta\boldsymbol{J}) - \boldsymbol{\nabla}\times\boldsymbol{J}_\mathrm{Hall}, \end{gather}

with variables $\rho$ density, $\boldsymbol {v}$ velocity, $T$ temperature, $\boldsymbol {B}$ magnetic field, $p = \rho T$ pressure and $\boldsymbol {J} = \boldsymbol {\nabla }\times \boldsymbol {B}$ the current density. We adopt a suitable dimensionalisation, so dimensional factors like the gas constant or the permeability of vacuum no longer appear. The adiabatic index is denoted by $\gamma$, and taken equal to $5/3$ usually. Additionally, $\boldsymbol {g}$ is the (external) gravitational acceleration, $\mathscr {L}$ the heat loss function, defined as energy losses (optically thin radiation) minus energy gains (e.g. heating), $\boldsymbol {\kappa }$ the thermal conduction tensor and $\eta$ the resistivity. In this follow-up work, we focus on the additional viscous and Hall effects. The viscous force term $\boldsymbol {F}_\mathrm {visc}$ in the momentum equation (2.2) and the viscous heating term $H_\mathrm {visc}$ in the energy equation (2.3) are introduced in § 3.2. The Hall term $\boldsymbol {J}_\mathrm {Hall}$ in the induction equation (2.4) is detailed in § 3.3.

These eight nonlinear partial differential equations are linearised around a one-dimensionally varying equilibrium

(2.5)\begin{equation} \left. \begin{aligned} & \rho_0 = \rho_0(u_1), \quad \boldsymbol{v}_0 = v_{02}(u_1)\,\boldsymbol{u}_2 + v_{03}(u_1)\boldsymbol{u}_3 ,\\ & T_0 = T_0(u_1), \quad \boldsymbol{B}_0 = B_{02}(u_1)\,\boldsymbol{u}_2 + B_{03}(u_1)\,\boldsymbol{u}_3, \end{aligned} \right\} \end{equation}

where in Cartesian coordinates $u_1$ is the $x$-coordinate, and $\boldsymbol {u}_2$ and $\boldsymbol {u}_3$ are the unit vectors in the then invariant $y$- and $z$-directions, respectively. In cylindrical coordinates, $u_1$ is the radial coordinate, $\boldsymbol {u}_2$ is then a unit vector in the angular direction and $\boldsymbol {u}_3$ is aligned along the cylinder axis. The dependence of the background equilibrium state on the $u_1$-coordinate is considered on a bounded domain (e.g. a slab of a plane-parallel atmosphere of a given vertical extent, or a flux tube of given radius), whilst the other coordinates are unrestricted (but the $u_2$-coordinate is periodic in the cylindrical case). The linearisation introduces the perturbations $(\rho _1, \boldsymbol {v}_1, T_1, \boldsymbol {B}_1)$, which in principle are fully three-dimensionally structured, time-dependent functions. Note that the indices $0$ and $1$ refer to equilibrium and perturbed quantities, respectively.

Subsequently, after adopting a vector potential $\boldsymbol {A}_1$ to describe the perturbed magnetic field as $\boldsymbol {B}_1 = \boldsymbol {\nabla }\times \boldsymbol {A}_1$, a 3-D Fourier analysis is applied to all perturbed quantities ($\rho _1$, $\boldsymbol {v}_1$, $T_1$, $\boldsymbol {A}_1$) as

(2.6)\begin{equation} f_1(\boldsymbol{r},t) = \hat{f}_1(u_1)\exp\left[ \mathrm{i}\left( k_2 u_2 + k_3 u_3 - \omega t \right) \right], \end{equation}

introducing the wavevector $\boldsymbol {k} = k_2\,\boldsymbol {u}_2 + k_3\,\boldsymbol {u}_3$ and the frequency $\omega$ (note that in the cylindrical case, $k_2$ is an integer usually denoted by $m$, enforcing annular periodicity). In essence, this reduces the problem to a generalised eigenvalue problem

(2.7)\begin{equation} \boldsymbol{\mathsf{A}}\boldsymbol{x} = \omega\boldsymbol{\mathsf{B}}\boldsymbol{x} \end{equation}

for matrices $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$, and the state vector $\boldsymbol {x} = (\rho _1, \boldsymbol {v}_1, T_1, \boldsymbol {A}_1)^\top$. Subsequently, it is transformed using a finite element method, where the domain is discretised using a specified number of grid points $N$ and linear combinations of basis functions are used to approximate all perturbed quantities $\hat {f}_1(u_1)$ in every subinterval. Note that this discretisation and the subsequent construction of the matrices $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$ means that the number of output eigenmodes is directly related to the number of grid points $N$ (for more information, see Claes et al. 2020). The resulting eigenproblem is passed to a QR solver from the LAPACK library (Anderson et al. 1999). The linearised equations of this eigenvalue problem are given in Appendix A with the inclusion of the new viscous and Hall contributions, listed in § 3.2 and § 3.3. The Legolas code then returns couples of eigenvalues $\omega$ and state vectors $\boldsymbol {x} = (\rho _1, \boldsymbol {v}_1, T_1, \boldsymbol {A}_1)^\top$, each of which describes a fundamental linear wave of the system. Any system may have both discrete and continuous solutions, such that all of these eigenmodes in the spectrum either belong to a continuum, or correspond to a discrete solution or overtone thereof (see e.g. Goedbloed et al. 2019).

For the boundaries in the $u_1$-direction we consider perfectly conducting walls, i.e.

(2.8a,b)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{u}_1 = 0, \quad \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{u}_1 = 0 \end{equation}

at the edges ($\boldsymbol {u}_1$ is the normal to the wall). The choice of equilibrium above guarantees that the equilibrium fields automatically satisfy these boundary conditions. For the perturbed quantities, these boundary conditions become

(2.9a,b)\begin{equation} v_1 = 0, \quad k_3 A_2 - k_2 A_3 = 0, \end{equation}

where $v_1$ is the first component of the velocity perturbation $\boldsymbol {v}_1 = (v_1,v_2,v_3)^\top$, and $A_2$ and $A_3$ denote components of the vector potential $\boldsymbol {A}_1 = (A_1, A_2, A_3)^\top$. However, instead of this second equation the more stringent condition $A_2 = A_3 = 0$ is imposed if both $k_2$ and $k_3$ are non-zero. If either wavevector component vanishes, only the corresponding $\boldsymbol {A}_1$-component is set to zero, i.e. if $k_2 = 0$ ($k_3 = 0$), the constraint reduces to $k_3 A_2 = 0$ ($k_2 A_3 = 0$) and only $A_2$ ($A_3$) is set to zero.

3. Code extensions and validation results

In this section, we present a couple of new features of the publicly available, open-source Legolas code. The first extension allows for the computation of physically relevant, derived quantities such as the perturbed magnetic field (as opposed to the auxiliary vector potential $\boldsymbol {A}_1$) or the entropy perturbation (as opposed to the density or temperature eigenfunction). These derived eigenfunctions are useful diagnostics, e.g. to evaluate specific eigenmode changes due to viscosity and Hall extensions. Similar to Claes et al. (2020), these viscosity and Hall modules were tested against a fair selection of recent to historic literature results, as will be demonstrated further. However, whilst the effects of flow, gravity, resistivity, thermal conduction and radiative cooling on full MHD spectra are relatively well documented, the literature contains many fewer fully reproducible tests on the linear effects of viscosity and the Hall current. Also, relevant published results were originally obtained using the incompressible MHD equations, i.e. satisfying $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v} = 0$ whilst typically ignoring the energy equation and/or mass conservation equation. Since the regular Legolas code uses the full set of compressible MHD equations, we want to be able to approximate the incompressible regime. In theory, this corresponds to the $\gamma \rightarrow \infty$ limit. In practice, we discard all terms in the energy equation's right-hand side except for the term $-(\gamma -1)p_0\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}_1$, which relates to the finite pressure perturbation. Then, the perturbed divergence becomes $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}_1 = \omega T_1 / T_0 (\gamma -1)$, whose real and imaginary parts clearly go to zero for $\gamma \rightarrow \infty$. Implementation-wise, $\gamma$ is set to a value of $10^{12}$ such that $|\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}_1|$ remains very small on the whole domain, usually yielding values smaller than $10^{-12}$ in the region of the spectrum near the origin. The value of $|\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}_1|$ is indeed observed to increase for increasingly large eigenvalues.

3.1. Derived eigenfunctions

There are various physical quantities of interest that can be derived from the eight eigenfunctions $(\rho _1, v_1, v_2, v_3, T_1, A_1, A_2, A_3)$ that Legolas computes (every eigenfunction here is now using the changed notation $f_1\equiv \hat {f}_1(u_1)$ and belonging to a specific set $(k_2, k_3, \omega )$). The most evident is the perturbed magnetic field $\boldsymbol {B}_1$, which is a combination of the $A_j$-eigenfunctions ($j=1,2,3$), through

(3.1)\begin{equation} \boldsymbol{B}_1 = \boldsymbol{\nabla}\times\boldsymbol{A}_1 = \mathrm{i}\left( \frac{k_2}{\varepsilon} A_3 - k_3 A_2 \right)\boldsymbol{u}_1 + \left( \mathrm{i} k_3 A_1 - A_3' \right)\boldsymbol{u}_2 + \frac{1}{\varepsilon} [(\varepsilon A_2)' - \mathrm{i} k_2 A_1 ]\boldsymbol{u}_3, \end{equation}

where a prime denotes the derivative with respect to $u_1$ from now on. We can similarly compute its divergence $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B}_1$ to validate that it is numerically zero, and its curl, $\boldsymbol {\nabla }\times \boldsymbol {B}_1$, yielding the perturbed current. Besides these magnetic-field-derived quantities, we can also determine the divergence of the velocity perturbation $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}_1$, which serves as a diagnostic tool when exploiting the incompressible approximation, and the perturbed vorticity $\boldsymbol {\nabla }\times \boldsymbol {v}_1$. Further worth mentioning is the entropy perturbation

(3.2)\begin{equation} S_1 = (p\rho^{-\gamma})_1 = \rho_0^{1-\gamma} T_1 + (1-\gamma) \rho_0^{-\gamma} T_0 \rho_1, \end{equation}

where we used the ideal gas law $p = \rho T$ to write the entropy in terms of density and temperature.

In addition, in the presence of an equilibrium magnetic field $\boldsymbol {B}_0$, all perturbed vector quantities ($\boldsymbol {B}_1$, $\boldsymbol {\nabla }\times \boldsymbol {B}_1$, $\boldsymbol {v}_1$, $\boldsymbol {\nabla }\times \boldsymbol {v}_1$) can be expressed in a reference frame consisting of a component along the equilibrium magnetic field and two perpendicular components. This is of interest to verify or determine the (theoretically expected) polarisations of specific eigenmodes. Note that the unit vector $\boldsymbol {u}_1$ is always perpendicular to the equilibrium magnetic field due to the chosen equilibrium form (2.5).

3.2. Viscosity

In MHD, viscosity appears as a force term $\boldsymbol {F}_\mathrm {visc}$ in the momentum equation (2.2). In its most general form, $\boldsymbol {F}_\mathrm {visc}$ can be written as $\boldsymbol {F}_\mathrm {visc} = -\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\pi }$, where $\boldsymbol {\pi }$ denotes the viscous stress tensor (Braginskii 1965). However, as shown by Erdélyi & Goossens (1995), where the authors used the full viscous stress tensor, only the shear viscosity contributes to resonant absorption, and the compressive and perpendicular components have negligible effects. Hence, for a constant dynamic viscosity $\mu$ the viscous force is in good approximation equal to (see e.g. Goedbloed et al. 2019)

(3.3)\begin{equation} \boldsymbol{F}_\mathrm{visc} = \mu\left[\nabla^2\boldsymbol{v} + \frac{1}{3}\boldsymbol{\nabla}(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v})\right]. \end{equation}

The linearisation of this expression is implemented in Legolas. Sometimes it is assumed that the kinematic viscosity $\nu =\mu /\rho$ is constant rather than the dynamic viscosity. This introduces additional terms in the linearisation, which are not implemented in Legolas. Note, however, that in a Cartesian set-up with constant $\rho _0$ and $\boldsymbol {v}_0$ the additional terms introduced by assuming a constant kinematic viscosity, rather than a constant dynamic viscosity, vanish.

In addition to a contribution in the momentum equation, viscosity also adds a viscous heating term $(\gamma -1) H_\mathrm {visc}$ to the energy equation (2.3). The source term $H_\mathrm {visc}$ is given by $H_\mathrm {visc} = -(\boldsymbol {\pi }\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {\cdot }\boldsymbol {v}$ and is approximately (see e.g. Goedbloed et al. 2019)

(3.4)\begin{equation} H_\mathrm{visc} \approx \mu\left|\boldsymbol{\nabla} \boldsymbol{v} \right|^2. \end{equation}

The linearisation of this approximation in Legolas assumes the Frobenius norm, resulting in the linearised term

(3.5)\begin{equation} H_{\mathrm{visc},1} = 2\mu \sum_{i=1}^3 \sum_{j=1}^3 (\boldsymbol{\nabla}\boldsymbol{v}_0)_{ij} (\boldsymbol{\nabla}\boldsymbol{v}_1)_{ij}. \end{equation}

Note that this contribution vanishes if the equilibrium flow is constant or zero, as is the case for the set-up of § 3.2.2 whilst it introduces two non-zero terms for the set-up in § 3.2.1. However, since both set-ups employ the incompressible approximation, which eliminates the energy equation, this term is not represented in either test case. Note further that the background equilibrium flow $\boldsymbol {v}_0$ is always adopted as a stationary, Eulerian flow, much like one normally computes eigenspectra for an ideal MHD equilibrium with time-independent $\boldsymbol {B}_0$, even in the presence of a finite resistivity. The viscous terms are thus omitted in the equilibrium equations.

The inclusion of viscosity also imposes additional no-slip boundary conditions at a rigid wall. In essence, this implies that the total plasma velocity at the boundary equals the wall's velocity. Implementation-wise, we impose that the velocity perturbation $\boldsymbol {v}_1$ at the boundary is exactly zero,

(3.6)\begin{equation} v_1 = v_2 = v_3 = 0. \end{equation}

As a consequence of the no-slip boundary condition, a non-zero equilibrium velocity at a boundary then simulates a boundary moving at that constant speed. We use this to study a viscous hydrodynamic Taylor–Couette flow below, which serves as a test case for cylindrical geometry by comparing with results of Gebhardt & Grossmann (1993). Additionally, the new viscosity module is tested in an MHD set-up using the results from Dahlburg et al. (1983), where the authors study the tearing mode in a viscoresistive plasma slab (i.e. in Cartesian geometry).

3.2.1. Viscous eigenspectra for Taylor–Couette flow

When considering a viscous fluid confined between two concentric cylinders that both rotate with constant angular velocity, the flow established under no-slip boundary conditions is called Taylor–Couette flow. A hydrodynamic equilibrium of this form, studied spectroscopically under incompressible conditions in Gebhardt & Grossmann (1993), is given by a uniform density $\rho _0$, and temperature and velocity profiles

(3.7a,b)\begin{equation} T_0(r) = \frac{1}{2} \left( A^2 r^2 + 4AB \log(r) - \frac{B^2}{r^2} \right) + C,\quad \boldsymbol{v}_0(r) = \left( Ar + \frac{B}{r} \right)\boldsymbol{u}_2, \end{equation}

where $C$ is an arbitrary constant to guarantee that $T_0$ is positive everywhere and

(3.8a,b)\begin{equation} A = \dfrac{\beta - \left( \dfrac{R_1}{R_2} \right)^2 \alpha}{1 - \left( \dfrac{R_1}{R_2} \right)^2},\quad B ={-}\dfrac{R_1^2 (\beta - \alpha)}{1 - \left( \dfrac{R_1}{R_2} \right)^2}, \end{equation}

with $\alpha$ ($\beta$) the angular speed of the inner (outer) cylinder at radius $R_1$ ($R_2$). Since this is a hydrodynamic test case, there is no equilibrium magnetic field $\boldsymbol {B}_0$.

Using the incompressible approximation in Legolas, four representative eigenspectra from figure 3 in Gebhardt & Grossmann (1993) are recovered in figure 1(a–d). These spectra feature two different types of modes, as described in Gebhardt & Grossmann (1993), namely translational modes with $v_1 = 0 = v_2$ (with only a non-trivial $v_3(u_1)=v_z(r)$ variation) and ‘azimuthal’ modes with non-zero $v_1$ and $v_2$ eigenfunctions and $v_3 = 0$. Admittedly, we recover the azimuthal modes but they do have a non-vanishing $v_3$ eigenfunction in the test cases with Legolas's incompressible approximation. However, we find that the spectra in figure 1(a,d) are almost identical in the compressible case (the spectrum is more heavily influenced by compressibility for smaller radius-to-thickness ratios), and the $v_3$ eigenfunction indeed vanishes in the compressible set-up. Hence, the incompressible approximation results in slightly spurious $v_3$ eigenfunctions for this case, presumably due to the vanishing of the $k_3v_3$ term in $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}$ because $k_3$ is zero. The $v_3$ eigenfunction is shown for a translational mode in figure 1(e), and the $v_2$ and $v_3$ eigenfunctions of an azimuthal mode are shown in figures 1( f) and 1(g), respectively. The corresponding modes are marked in figure 1(d). These eigenmodes in Legolas ($\omega _{L}$) are consistent with those reported by Gebhardt & Grossmann (1993) ($\omega _{G}$), namely $\omega _{L} = 1583.50-1053.70\,\mathrm {i}$ and $\omega _{L} = 2406.82-579.55\,\mathrm {i}$ compared with $\omega _{G} = 1583.56-1053.83\,\mathrm {i}$ and $\omega _{G} = 2406.81-579.54\,\mathrm {i}$. The Legolas eigenfunctions match their eigenfunctions up to a complex factor (this represents the freedom to choose a reference amplitude and phase in a linear eigenvalue problem) when comparing with their figures 6(b) and 8(b).

Figure 1. Parts of the incompressible ($\gamma \rightarrow \infty$) Taylor–Couette spectrum with an inner cylinder at rest ($\alpha = 0$) for $k_3 = 0$, $\rho _0 = 1$, $\mu = 1$ and different parameter choices: (a) $k_2 = 1$, $R_1 = 7/3$, $R_2 = 10/3$, $\beta = 3\times 10^3$; (b) $k_2 = 2$, $R_1 = 1$, $R_2 = 2$, $\beta = 2.5\times 10^3$; (c) $k_2 = 2$, $R_1 = 0.25$, $R_2 = 1.25$, $\beta = 2\times 10^3$; and (d) $k_2 = 3$, $R_1 = 9$, $R_2 = 10$, $\beta = 10^3$. The modes represented by a dot (or $\blacklozenge$ in (d)) are translational modes with $v_1$ and $v_2$ numerically zero whilst the crosses (and $\blacksquare$ in (d)) represent azimuthal modes with non-zero $v_1$ or $v_2$ components. (e) The $v_3$ eigenfunction of the translational (d)-eigenvalue $\omega = 1583.50-1053.70\,\mathrm {i}$ ($\blacklozenge$). ( f,g) The $v_1$ and $v_2$ eigenfunction, respectively, of the azimuthal (d)-eigenvalue $\omega = 2406.82\textrm {-}579.55\,\mathrm {i}$ ($\blacksquare$). Solid lines represent real parts, dotted lines imaginary parts. All runs were performed at $251$ grid points.

Taking the analysis of Taylor–Couette flow one step further using the fully compressible functionality of the code, we take a look at the entropy perturbation in the compressible spectrum, motivated by the observation that the azimuthal modes have a non-zero entropy perturbation whilst the translational modes have no entropy variation. In particular, we compare the entropy perturbation with and without the inclusion of viscous heating in the energy equation. Here, we use the set-up of figure 1(a) again, with parameters $k_2 = 1$, $k_3 = 0$, $R_1 = 7/3$, $R_2 = 10/3$, $\beta = 3\times 10^3$, $\rho _0 = 1$ and $\mu = 1$, but without the incompressible approximation, i.e. $\gamma = 5/3$. The compressible spectrum (without viscous heating) is shown in figure 2(a). It is extremely similar to the corresponding incompressible spectrum in figure 1(a) and is hardly influenced by viscous heating. The entropy perturbations of an azimuthal mode ($\omega = 2529.12 - 485.63\,\mathrm {i}$ without viscous heating; $\omega = 2529.66 - 485.80\,\mathrm {i}$ with viscous heating; marked by $\blacksquare$ in figure 2a) are shown in figure 2(b,c) for the compressible case without and with viscous heating, respectively. The viscous heating introduces a limited but noticeable change in the entropy perturbation $S$.

Figure 2. (a) Part of the compressible Taylor–Couette spectrum (3.7a,b), with parameters $k_2 = 1$, $k_3 = 0$, $R_1 = 7/3$, $R_2 = 10/3$, $\beta = 3\times 10^3$, $\rho _0 = 1$ and $\mu = 1$. Dots represent translational modes, crosses (and $\blacksquare$) are azimuthal modes. (b) Entropy perturbation $S$ of the azimuthal mode $\omega = 2529.12 - 485.63\,\mathrm {i}$ ($\blacksquare$) without the inclusion of viscous heating. (c) Entropy perturbation $S$ of the azimuthal mode $\omega = 2529.66 - 485.80\,\mathrm {i}$ ($\blacksquare$) with the influence of viscous heating. Solid lines represent real parts, dotted lines imaginary parts. Both runs were performed at $251$ grid points.

3.2.2. Viscoresistive plasma slab

As an MHD test case, consider the incompressible, viscoresistive stability analysis of a plane-parallel plasma slab from Dahlburg et al. (1983), with equilibrium magnetic field profile

(3.9)\begin{equation} \boldsymbol{B}_0 = \left( \arctan \alpha x - \frac{\alpha x}{1+\alpha^2} \right)\,\boldsymbol{u}_2, \end{equation}

with parameter $\alpha$, uniform density $\rho _0$ and $T_0$ positive and satisfying the constant total pressure condition $\partial (\rho _0 T_0 + \frac {1}{2}\boldsymbol {B}_0^2)/\partial x = 0$. Note that the field is not force-free, and induces a current

(3.10)\begin{equation} \boldsymbol{J}_0 = \alpha \left( \frac{1}{1+\alpha^2x^2} - \frac{1}{1+\alpha^2} \right)\,\boldsymbol{u}_3, \end{equation}

which vanishes at $x=\pm 1$, where we introduce perfectly conducting walls with a no-slip boundary condition. The simultaneous inclusion of resistivity and viscosity in the linear stability analysis leads to different tearing mode regimes, based on the resistivity $\eta$ and dynamic viscosity $\mu$ coefficients. The formulation in Dahlburg et al. (1983) actually uses the kinematic viscosity $\nu = \mu /\rho$, but since they assume a uniform density and no equilibrium flow, our constant $\mu$ formulation is equivalent. The relation between the resistivity $\eta$ and the kinematic viscosity $\nu$ is often expressed in terms of the magnetic Prandtl number $Pm = \nu /\eta = \mu /\rho _0\eta$. In the remainder of this section, $Pm$ will vary between $10^{-4}$ and $10^4$. Note that $\rho _0 = 1$ in all examples in this section, such that the Prandtl number reduces to $Pm = \mu /\eta$.

In Dahlburg et al. (1983), the authors give numerical values for the purely unstable tearing eigenmode and show the $v_{1x}$ and $B_{1x}$ eigenfunctions of the tearing mode for a few different values of $\eta$ and $\mu$ as well as the evolution of the tearing mode growth rate as a function of $\eta$, $\mu$ and the parallel wavenumber $k_2$. Here, we reproduce these results using Legolas.

First, we recover the eigenvalues and eigenfunctions for three cases: (a) $\eta = \mu = 10^{-3}$, (b) $\eta = 0.1$, $\mu = 10^{-5}$ and (c) $\eta = 10^{-5}$, $\mu = 0.1$, all with $\rho _0 = 1$, $\alpha = 10$ and $\boldsymbol {k} = \boldsymbol {u}_2$. Due to Legolas's incompressible approximation, the tearing modes from Legolas ($\omega _{L}$) deviate slightly from those reported in Dahlburg et al. (1983) ($\omega _{D}$), namely (a) $\omega _{L} = 0.1965\,\mathrm {i}$ compared with $\omega _{D} = 0.19687\,\mathrm {i}$, (b) $\omega _{L} = 0.4393\,\mathrm {i}$ compared with $\omega _{D} = 0.4397\,\mathrm {i}$ and (c) $\omega _{L} = 0.002531\,\mathrm {i}$ compared with $\omega _{D} = 0.002537\,\mathrm {i}$. The $v_{1x}$ and $B_{1x}$ eigenfunctions, defined up to a complex factor and rescaled here for comparison with figure 2(c–e) in Dahlburg et al. (1983), are shown in figure 3(a–c). The arbitrary complex factor is chosen for each case such that all shown eigenfunctions are real.

Figure 3. The $v_{1x}$ and $B_{1x}$ eigenfunctions are shown for the equilibrium (3.9) with $\rho _0 = 1$, $\alpha = 10$, $\boldsymbol {k} = \boldsymbol {u}_2$ and (a) $\eta = \mu = 10^{-3}$, (b) $\eta = 0.1$, $\mu = 10^{-5}$ and (c) $\eta = 10^{-5}$, $\mu = 0.1$. (d) Growth rate as a function of $\eta ^{-1}$ for given values of $\mu$. (e) Growth rate as a function of $\mu ^{-1}$ for given values of $\eta$. ( f) Growth rate as a function of $\boldsymbol {k} = k_2\,\boldsymbol {u}_2$ for (A) $\eta = 10^{-2}$ and $\mu = 10^{-2}$; (B) $\eta = 10^{-3}$ and $\mu = 10^{-2}$; (C) $\eta = 10^{-2}$ and $\mu = 10^{-3}$; (D) $\eta = 2\times 10^{-3}$ and $\mu = 2\times 10^{-3}$. All runs were performed at $201$ grid points.

Next, we reproduce the evolution of the tearing mode growth rate (d) as a function of $\eta ^{-1}$ for fixed values of $\mu$, (e) as a function of $\mu ^{-1}$ for fixed values of $\eta$ and ( f) as a function of $k_2$ for fixed values of $\eta$ and $\mu$. In figure 3(d–f) the tearing growth rate is shown for each configuration of parameters, obtaining the positive growth rates shown in Dahlburg et al. (1983) in their figures 6(a), 7 and 9, respectively. Once again, they agree very well. Note that each marker represents a single Legolas run at $201$ grid points.

Unlike Dahlburg et al. (1983), Legolas not only computes the tearing mode, but the entire spectrum. Hence, we can also compare the purely resistive, purely viscous and truly viscoresistive spectra for the same equilibrium profile (3.9). This is shown in figure 4 for runs at $251$ grid points. In this figure, the left column displays the incompressible limit and the right column shows the fully compressible spectra. All spectra are supplemented with the analytical, ideal MHD slow and Alfvén continua, which correspond to singular solutions of the ordinary differential equation obtained through a reformulation of the ideal MHD equations in terms of the $x$-component of the Lagrangian displacement field. It can be shown for homogeneous backgrounds (see e.g. Goedbloed et al. 2019) that in the presence of resistivity the Alfvén and slow modes trace out semi-circles in the stable part of the spectrum with infinitely degenerate (collapsed) continua. For inhomogeneous resistive spectra the ideal continuum ranges will relocate to collections of discrete modes in the stable half-plane, still resembling the semi-circular curves, as seen in figure 4(a,d), which will have links to extremal or edge values of the ideal continua. Figure 4(b,e) now shows that viscosity exerts a similar influence as resistivity. Finally, figure 4(c, f) represents modified variants of the semi-circle-like curves in the other panels, due to the combined effects of viscosity and resistivity. Since the ideal slow and Alfvén continua partially overlap and are symmetric with respect to the imaginary axis, the slow continuum is only drawn in the left half-plane (red dashed line) and the Alfvén continuum in the right half-plane (cyan solid line). In the left column, the slow continua are eliminated by the incompressible assumption.

Figure 4. Comparison of spectra for equilibrium profile (3.9) for resistive (a,d), viscous (b,e) and viscoresistive (c, f) cases. The left column (a–c) represents the incompressible approximation, the right column (d–f) the compressible equations. All runs use parameters $\boldsymbol {k} = \boldsymbol {u}_2$, $\rho _0 = 1$ and $\alpha = 10$. In case of resistivity (viscosity), the parameter is $\eta = 10^{-3}$ ($\mu = 10^{-3}$). The cyan solid and red dashed lines represent the ideal MHD Alfvén and slow continua, respectively. Both continua are symmetric with respect to the imaginary axis, but only shown in one half-plane to avoid overlap. All runs were performed at $251$ grid points.

The first row of figure 4(a,d) shows the resistive slab with $\eta = 10^{-3}$. This case is well known and discussed in e.g. Goedbloed et al. (2019). In panel (a), the Alfvén modes form a semicircle and the slow (magnetoacoustic) modes are eliminated by the incompressible approximation. In the compressible case of panel (d), the slow modes reappear as the inner semicircle. Finally, both the compressible and incompressible spectra feature a resistive tearing mode, as the only purely unstable eigenmode of this system for the chosen parameters.

The second row of figure 4(b,e) on the other hand shows the viscous case with $\mu = 10^{-3}$. The result looks surprisingly similar to the resistive case, with both the slow and Alfvén modes taking on the same semicircular shape of similar magnitude (note that the axes are scaled identically in panels a,b,d,e). Whilst there are many minute differences with the resistive case in the first row, the key difference is the absence of a tearing mode in the viscous spectra.

Ultimately, the third row (panels c, f) shows the viscoresistive spectrum with $\eta = \mu = 10^{-3}$. Although resistivity and viscosity exert a similar influence on the slow and Alfvén modes when they are the only physical effect in consideration, the combination of both effects reveals new behaviour in both the incompressible (panel c) and compressible case (panel f). Whilst the slow and Alfvén branches still originate in the same point on the real axis, the semicircular structures are replaced by stretched-out curves along the imaginary axis. The resistive tearing mode is still present, but damped by the viscosity. Therefore, the changes are most pronounced on the stable and damped parts of the spectrum, whose physical relevance must also consider the fact that the ideal MHD equilibrium itself will evolve on a specific diffusive time scale when viscoresistive effects are active.

3.3. Hall-MHD

The induction equation, given by the Maxwell–Faraday equation

(3.11)\begin{equation} \frac{\partial\boldsymbol{B}}{\partial t} ={-}\boldsymbol{\nabla}\times\boldsymbol{E}, \end{equation}

can be extended to include the effects of the Hall current, electron pressure and electron inertia by expressing the electric field $\boldsymbol {E}$ using the (dimensionless) generalised Ohm's law

(3.12)\begin{equation} \boldsymbol{E} ={-}\boldsymbol{v}\times\boldsymbol{B} + \eta\boldsymbol{J} + \frac{\eta_{H}}{\rho} \left( \boldsymbol{J}\times\boldsymbol{B} - \boldsymbol{\nabla} p_{e} \right) + \frac{\eta_{e}}{\rho} \frac{\partial\boldsymbol{J}}{\partial t}. \end{equation}

Here, $p_{e}$ denotes the electron pressure and is related to $p$ through the electron fraction $f_{e}$ as $p_{e} = f_{e}p$, with $f_{e} = n_{e} / (n_{e}+n_{p}) = 1/2$ for a charge-neutral electron–proton plasma. Furthermore, $\eta _{H}$ and $\eta _{e}$ are the normalised Hall and electron inertia coefficients,

(3.13a,b)\begin{equation} \eta_{H} = \frac{m_{i}}{e} \frac{V_{R}}{L_{R} B_{R}},\quad \eta_{e} = \frac{m_{e} m_{i}}{e^2} \left( \frac{V_{R}}{L_{R} B_{R}} \right)^2, \end{equation}

respectively. Here, $e$ denotes the fundamental charge, $m_{i}$ and $m_{e}$ are the ion and electron mass, respectively, and $V_{R}$, $L_{R}$ and $B_{R}$ are the reference velocity, length and magnetic field strength. Consequently, the electron inertia coefficient $\eta _{e}$ is several orders of magnitude smaller than the Hall coefficient $\eta _{H}$. Therefore, the effect of electron inertia is usually negligible. Hence, most results in the literature do not include it. Whilst this effect is implemented in Legolas, the reference tests that follow all set $\eta _e = 0$, so its effect is only quantified for one limit case here. It should also be pointed out that any equilibrium of the form (2.5) that satisfies the ideal MHD equilibrium conditions, also satisfies the Hall-MHD equilibrium conditions (neglecting electron inertia) because the Hall term in the induction equation (2.4), given by

(3.14)\begin{equation} \boldsymbol{\nabla}\times\boldsymbol{J}_{\mathrm{Hall},0} = \boldsymbol{\nabla}\times \left[ \frac{\eta_{H}}{\rho_0} \left( \boldsymbol{J}_0\times\boldsymbol{B}_0 - \boldsymbol{\nabla} p_{{e}0} \right) \right], \end{equation}

reduces to zero. For the first term this follows because $\boldsymbol {B}_0$ and $\boldsymbol {J}_0$ both lie in the $\boldsymbol {u}_2\boldsymbol {u}_3$-plane, such that their vector product is proportional to $\boldsymbol {u}_1$, and since they only depend on $u_1$, this implies that $\boldsymbol {\nabla }\times f(u_1)\,\boldsymbol {u}_1 = 0$. The second term is the curl of a gradient, which is always zero. Note that whilst there are many similarities between the resistive and Hall terms, here they differ since the resistive term does not disappear in the induction equation (2.4) for an equilibrium of the form (2.5). As pointed out in Claes et al. (2020), the resistive term is neglected in the equilibrium equations by assuming that the time scales on which magnetic fields decay is much larger than the time scales of resistive modes.

The Hall and electron pressure terms in the induction equation are not implemented directly in Legolas as written in equation (3.12). Instead, $\boldsymbol {J}\times \boldsymbol {B}$ is substituted into this expression using the momentum equation (2.2) as done in e.g. Ahedo & Ramos (2009) because it is observed to be more stable numerically. The result is

(3.15)\begin{align} \boldsymbol{E} & ={-}\boldsymbol{v}\times\boldsymbol{B} + \eta\boldsymbol{J} + \eta_{H} \left\{ \frac{\partial\boldsymbol{v}}{\partial t} + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{v} - \boldsymbol{g} - \frac{\mu}{\rho} \left[ \nabla^2\boldsymbol{v} + \frac{1}{3}\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}) \right] + \frac{\boldsymbol{\nabla} p_{i}}{\rho} \right\}\nonumber\\ & \quad + \frac{\eta_{e}}{\rho} \frac{\partial\boldsymbol{J}}{\partial t}. \end{align}

This equation now features the ion pressure $p_i$ instead, which is related to the total pressure as $p_{i} = (1-f_{e})p$. Note that this expression for the electric field now has two time derivatives on the right-hand side, which will then enter the induction equation. Exploiting $\boldsymbol {A}_1$ instead of $\boldsymbol {B}_1$, the linearised induction equation becomes $\partial \boldsymbol {A}_1/\partial t = -\boldsymbol {E}_1$. Hence, we linearise equation (3.15), allowing for a temperature-dependent Spitzer resistivity $\eta (T)$, which gives

(3.16)\begin{align} \boldsymbol{E}_1 & ={-}\boldsymbol{v}_1\times\boldsymbol{B}_0 - \boldsymbol{v}_0\times\left(\boldsymbol{\nabla}\times\boldsymbol{A}_1\right) + \eta_0 \boldsymbol{\nabla}\times\left(\boldsymbol{\nabla}\times\boldsymbol{A}_1\right) + \frac{\mathrm{d}\eta}{\mathrm{d}T} T_1 \boldsymbol{\nabla}\times\boldsymbol{B}_0\nonumber\\ & \quad +\eta_{H} \left\{ \frac{\partial\boldsymbol{v}_1}{\partial t} + \boldsymbol{v}_1\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{v}_0 + \boldsymbol{v}_0\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{v}_1 - \frac{\mu}{\rho_0} \left[ \nabla^2\boldsymbol{v}_1 + \frac{1}{3}\boldsymbol{\nabla}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}_1\right) \right]\right.\nonumber\\ & \left.\quad+ \mu\frac{\rho_1}{\rho_0^2} \left[ \nabla^2\boldsymbol{v}_0 + \frac{1}{3}\boldsymbol{\nabla}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}_0\right) \right] + \frac{\boldsymbol{\nabla} p_{{i}1}}{\rho_0} - \frac{\rho_1\boldsymbol{\nabla} p_{{i}0}}{\rho_0^2} \right\}\nonumber\\ & \quad + \frac{\eta_{e}}{\rho_0}\frac{\partial}{\partial t} \left[ \boldsymbol{\nabla}\times(\boldsymbol{\nabla}\times\boldsymbol{A}_1) \right] - \eta_{e}\frac{\rho_1}{\rho_0^2} \frac{\partial}{\partial t} \left( \boldsymbol{\nabla}\times\boldsymbol{B}_0 \right). \end{align}

This expression can be simplified by observing that the term $\boldsymbol {\nabla } p_{{i}1}/\rho _0$ can be written as

(3.17)\begin{equation} \frac{\boldsymbol{\nabla} p_{{i}1}}{\rho_0} = \boldsymbol{\nabla}\left( \frac{p_{{i}1}}{\rho_0} \right) + \frac{\boldsymbol{\nabla}\rho_0}{\rho_0^2} p_{{i}1}. \end{equation}

Hence, this term is a pure gradient if $\rho _0$ is uniform. Since the electric field is only defined up to a gradient, we can redefine it as $\widetilde {\boldsymbol {E}}_1 = \boldsymbol {E}_1 - \eta _{H} \boldsymbol {\nabla } (p_{{i}1}/\rho _0)$ (with $\boldsymbol {E}_1$ expression (3.16)) such that after substituting (3.17) $p_{{i}1}$ only appears in $\widetilde {\boldsymbol {E}}_1$ in the term $\eta _{H} p_{{i}1}\boldsymbol {\nabla }\rho _0/\rho _0^2$, and thus only if $\rho _0$ is not uniform. The resulting induction equation $\partial \boldsymbol {A}_1/\partial t = -\widetilde {\boldsymbol {E}}_1$ is implemented in the code. Note that in the generalised eigenvalue problem (2.7) resulting from the Fourier analysis (2.6) the time derivatives in the Hall and electron inertia terms enter in the $\boldsymbol{\mathsf{B}}$-matrix and break its former symmetry.

In what follows, we present a series of stringent test cases to validate our Hall-MHD linear solver.

3.3.1. Hall-MHD waves in a homogeneous plasma slab

For the first test case, consider a homogeneous Cartesian plasma slab confined between two perfectly conducting walls (perpendicular to the $x$-axis). The equilibrium is given by

(3.18a–c)\begin{equation} \rho_0 = 1, \quad T_0 = 1, \quad \boldsymbol{B}_0 = \boldsymbol{u}_3, \end{equation}

and our normalisation is chosen such that $\eta _{H} = 1$. We want to quantify Hall-MHD eigenmodes of this slab, where we can compare with the analytical result of waves for an infinite homogeneous plasma in Hameiri, Ishizawa & Ishida (2005), where they give the dispersion relation (temporarily reintroducing dimensions)

(3.19)\begin{gather} \left( \frac{\omega}{kv_{{A}0}} \right)^6 - \left( \frac{\omega}{kv_{{A}0}} \right)^4 \left[ 1 + \frac{\gamma T_0}{v_{{A}0}^2} + \cos^2\theta \left( 1 + \frac{(k\eta_{H})^2}{\rho_0} \right) \right]\nonumber\\ + \left( \frac{\omega}{kv_{{A}0}} \right)^2 \cos^2\theta \left[ 1 + \frac{\gamma T_0}{v_{{A}0}^2} \left( 2 + \frac{(k\eta_{H})^2}{\rho_0} \right) \right] - \frac{\gamma T_0}{v_{{A}0}^2} \cos^4\theta = 0. \end{gather}

Here, $k$ is the magnitude of the wavevector $\boldsymbol {k}$, $\theta$ is the angle between $\boldsymbol {k}$ and $\boldsymbol {B}_0$ and $v_{{A}0} = |\boldsymbol {B}_0|/\sqrt {\mu _0\rho _0}$ is the equilibrium Alfvén speed, where $\mu _0$ is the vacuum permeability. The inclusion of the Hall term introduces a length scale into the previously scale-independent MHD equations through the ion skin depth $d_{i} = \eta _{H}/\sqrt {\rho _0}$. This makes the Hall-MHD waves dispersive as seen from the above dispersion relation. To simulate an infinite medium, we need to ensure that the ratio of the equilibrium ion skin depth to the system size is sufficiently small. Hence, for the choice of $d_{{i}0} = 1$, we solve in the interval $x\in [0,10^3]$. The exact choice of interval size is largely arbitrary, but it should be kept in mind that when we increase the interval size, we may also be forced to increase the resolution to ensure that Legolas picks up the Hall modes. This is due to the fact that the grid resolution can be directly linked to the dimensions of the $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$ matrices in the eigenvalue problem (2.7), and thus also to the number of eigenvalues returned.

Since the medium in Legolas is bounded in the $x$-direction, each solution of the dispersion relation (3.19) should approximate the first mode in a sequence in the spectrum, which can be verified by the number of nodes in the mode's corresponding eigenfunctions. For given angles $\theta$, the first mode of each sequence is shown in figure 5(a) alongside the theoretical curves and a comparison with the ideal MHD dispersion relation. A full spectrum version is also shown in figure 5(b). As can be seen in (3.19), in the case of perpendicular propagation the dispersion relation is not influenced by the Hall parameter ($\cos \theta = 0$). There, the highest mode reduces to the regular fast MHD mode and the lower two modes (slow and Alfvén) vanish, also visible in figure 5(a,b). The spectrum itself is shown for an angle $\theta \approx 0.564$ in figure 5(c) alongside the analytical infinite-medium solutions, each one corresponding to the start of a sequence, indicated by vertical lines. The sequences themselves, whose modes are much more tightly packed than in the ideal MHD sequences, are shown in the insets of figure 5(c). The smallest sequence displays anti-Sturmian behaviour, similar to the ideal MHD slow modes, whilst the two larger sequences behave in a Sturmian way, like the ideal MHD Alfvén and fast modes.

Figure 5. (a) Comparison of the first mode of each sequence (dots) with the theoretical Hall prediction by Hameiri et al. (2005) (solid lines) and ideal MHD (dashed lines) as a function of the angle $\theta$ between $\boldsymbol {k} = {\rm \pi}(\sin \theta \,\boldsymbol {u}_2 + \cos \theta \,\boldsymbol {u}_3)$ and $\boldsymbol {B}_0 = \boldsymbol {u}_3$ with $\rho _0 = 1$, $T_0 = 1$, $\eta _{H} = 1$ and $x\in [0,10^3]$. (b) Comparison of the full spectrum with MHD and Hall-MHD solutions for the set-up from (a). (The isolated branches are unresolved modes.) (c) Spectrum for an angle $\theta \approx 0.564$ with the three (positive) solutions of the dispersion relation (3.19) as vertical (dotted) lines. (d) The $\rho _1$ eigenfunctions of the first three modes of the smallest solution sequence for $\theta \approx 1.007$. (e) Comparison of the full spectrum with ideal and Hall-MHD predictions for varying wavenumber for $\boldsymbol {k} = k\,(\boldsymbol {u}_2/2 + \sqrt {3}\,\boldsymbol {u}_3/2)$, $\boldsymbol {B}_0 = \boldsymbol {u}_3$, $\rho _0 = 1$, $T_0 = 1$, $\eta _{H} = 1$ and $x\in [0,10^3]$. All runs were performed at $501$ grid points.

Furthermore, the (real) $\rho _1$ eigenfunctions of the first three modes in the smallest sequence of the $\theta \approx 1.007$ spectrum are given in figure 5(d). Contrary to the slow, Alfvén, and fast modes in ideal MHD, the density perturbation vanishes at the edges here. This behaviour is easily derived from the equations in Appendix A for the adiabatic homogeneous set-up considered here. Neglecting equilibrium flow, resistivity and viscosity, and applying the perfectly conducting wall boundary conditions $\tilde {v}_1 = \tilde {a}_2 = \tilde {a}_3 = 0$ (with the tildes indicating the transformed variables A1) reduces the second and third components of the induction equation (A7) and (A8) to $\tilde {v}_2 = 0$ and $\tilde {v}_3 = 0$, respectively, for non-zero frequency. Since these equations vanish altogether for an adiabatic homogeneous set-up in ideal MHD, these emerging no-slip boundary conditions are naturally imposed by the Hall current. Using these newfound conditions alongside the others in the continuity equation (A2), the third component of the momentum equation (A5), and the energy equation (A9) reduces these equations to $\omega \tilde {\rho }_1 = -\rho _0 \tilde {v}_1'$, $k_3 (\tilde {\rho }_1 T_0 + \rho _0 \tilde {T}_1) = 0$ and $\omega \tilde {T}_1 = -(\gamma -1)T_0 \tilde {v}_1'$, respectively, where we also used that $B_{02} = 0$ in our reference frame (for a different reference frame a linear combination of (A4) and (A5) gives similar results for any constant $\boldsymbol {B}_0$). Combining all three implies that $\tilde {\rho }_1 = \tilde {T}_1 = \tilde {v}_1' = 0$ at the wall boundaries. Since this derivation made no assumptions about the waves, this behaviour is present for all modes in the adiabatic homogeneous Hall spectrum. Note, however, that it only holds for oblique angles between $\boldsymbol {k}$ and $\boldsymbol {B}_0$. If $\boldsymbol {k}$ is parallel to $\boldsymbol {B}_0$ and along the $y$- or $z$-axis in a reference frame of choice, either $\tilde {a}_2 = 0$ or $\tilde {a}_3 = 0$ does not hold and therefore the other equations do not reduce in the way described above. Therefore, for parallel propagation the density perturbation is non-zero at the boundaries, just like in the ideal MHD case.

Finally, figure 5(e) shows the dispersion of all three sequences by comparing the full spectrum for different wavenumbers with the theoretical ideal MHD and Hall-MHD predictions (the dispersive nature of the middle sequence is more subtle). This dispersive behaviour identifies the largest sequence as the so-called whistler wave, which derives its name from the property that higher frequencies travel faster such that higher frequencies reach observers at earlier times than lower frequencies. The smallest sequence is sometimes called the ion cyclotron wave because its frequency approaches $\varOmega _{i}\cos \theta$ asymptotically for increasing wavenumber, where $\varOmega _{i} = ZeB_0/m_{i}$ is the ion cyclotron frequency with $Z$ the charge number, $e$ the fundamental charge and $m_{i}$ the ion mass (in Legolas, the ions are protons). Note that the final (intermediate) mode in this panel, which is related to the MHD Alfvén wave and two-fluid A mode (De Jonghe & Keppens 2020), fails to capture the electron cyclotron resonance $\omega \rightarrow \varOmega _{e}\cos \theta$ ($\varOmega _{e} = eB_0/m_{e}$, with $m_{e}$ the electron mass) in the short wavelength (large wavenumber) limit, which is present in the two-fluid description, because $\eta _{e}$ ($\propto m_{e}$) was set to zero (Hameiri et al. 2005). It has been verified that the sequences indeed start at the theoretical Hall-MHD results (up to an error of $10^{-5}$ at $501$ grid points and a ratio of $10^{-3}$ of $\eta _{H}$ to slab thickness), even though it is somewhat unclear in this image due to the large frequency range and the proximity of various sequences.

As a quick test of the electron inertia term, it is shown in figure 6 that the inclusion of this term ($\eta _{e} \neq 0$) captures the behaviour near the electron cyclotron resonance for large wavenumbers.

Figure 6. Comparison with the theoretical Hall-MHD curves from Hameiri et al. (2005, where $\eta _{e} = 0$) and two-fluid resonances ($\varOmega _{e}\cos \theta$, $\varOmega _{i}\cos \theta$) of the frequency $\omega$ as a function of wavenumber $k$ in Hall-MHD (a) without electron inertia ($\eta _{e} = 0$) and (b) with electron inertia ($\eta _{e} \neq 0$). The set-up is identical to the one used in figure 5(e). All runs were performed at $501$ grid points.

3.3.2. Resistive Harris sheet: tearing in Hall-MHD

In Shi et al. (2020), the authors investigate the influence of the Hall current on the resistive tearing mode of a Harris current sheet. The equilibrium profile takes the form

(3.20a–c)\begin{equation} \rho_0 = \tilde{\rho}_0,\quad T_0 = \frac{B_0^2}{2\tilde{\rho}_0}\ \mathrm{sech}^2 \left( \frac{x}{a} \right),\quad \boldsymbol{B}_0 = B_0 \tanh\left( \frac{x}{a} \right)\,\boldsymbol{u}_2 + B_{g}\,\boldsymbol{u}_3, \end{equation}

with $\tilde {\rho }_0 = B_0 = a =1$ and $B_{g}$ a variable guide field parameter. The included physical effects are a constant resistivity $\eta = 10^{-4}$ and a Hall current with coefficient $\eta _{H} = 1$. As explained earlier, in the vector potential formulation in Legolas we include the Hall term and the electron pressure term, but in this test we ignore the electron inertia effect (i.e. $\eta _e=0$). Furthermore, Shi et al. (2020) assume incompressibility, so we also use the incompressible approximation (large $\gamma$) in Legolas.

Shi et al. (2020) solve for the tearing mode on the interval $x\in [-15,15]$ and assume exponential decay of the perturbed quantities outside of that interval since the profile (3.20a–c) is approximately constant there for the chosen parameters. In Legolas, the default boundary settings are conducting wall boundary conditions at a finite distance, which may modify the linear MHD spectrum due to e.g. wall stabilisation effects. However, for $a = 1$ the interval $[-15,15]$ seems large enough such that the stabilising influence of the conducting walls is expected to be negligible. Our results are shown for two different values of $k_2$ in figure 7, recovering figure 4 in Shi et al. (2020). In this figure, we show the growth rate $\mathrm {Im}(\omega )$ and frequency $\mathrm {Re}(\omega )$ in top and bottom panels, respectively. Each marker represents the tearing mode in a Legolas run at $501$ grid points and a Laplace distributed grid was used to focus the grid points around the region of strongest change in equilibrium magnetic field at $x = 0$. Note that the non-zero $\mathrm {Re}(\omega )$ values are due to the inclusion of the Hall terms, which results in spectrum asymmetry with respect to the imaginary axis here, similar to equilibrium flow. For any guide field value $B_{g}$, the growth rate is influenced by the wavevector, with the maximum growth rate depending on both wavevector components, $k_2$ and $k_3$. Whilst the real part of the frequency $\mathrm {Re}(\omega )$ has an extremum as a function of $k_3$ in the presence of a guide field ($B_{g} \neq 0$), $|\mathrm {Re}(\omega )|$ increases linearly with increasing $k_3$ in the absence of a guide field ($B_{g} = 0$), until the tearing mode is fully damped.

Figure 7. The real (b,d) and imaginary (a,c) parts of the tearing mode in a Harris current sheet (3.20a–c) as a function of $k_3$, with $\tilde {\rho }=1$, $a=1$, $B_0=1$, $\eta = 10^{-4}$ and $\eta _{H} = 1$ for different guide field strengths $B_{g}$. (a,b) $k_2 = 0.155$; (c,d) $k_2 = 0.5$. All runs were performed at $501$ grid points.

Besides quantifying the tearing mode complex eigenfrequencies, Shi et al. (2020) also reported on the tearing mode eigenfunctions. Up to a complex factor, the eigenfunctions obtained by Legolas, shown in figure 8, match the results in the first two rows of figure 7 in Shi et al. (2020). For the cases in figures 8(a,d) and 8(b,e), where $B_{g} = 0$, the $B_1$ and $v_1$ eigenfunctions are symmetric and antisymmetric, respectively, with respect to the location of the Harris sheet ($x=0$) whereas the (anti)symmetry is broken with the introduction of a non-zero guide field $B_{g}$ (panel c, f).

Figure 8. The incompressible tearing mode's $v_1$ and $B_1$ eigenfunctions for $\eta = 10^{-4}$, $\eta _{H} = 1$ and $k_2 = 0.5$, with different values of $k_3$ and $B_{g}$: (a,d) $k_3 = 0$, $B_{g} = 0$; (b,e) $k_3 = 0.06$, $B_{g} = 0$; and (c, f) $k_3 = 0.06$, $B_{g} = 5$. Real parts are shown as solid (blue) lines, imaginary parts as dotted (orange) lines. All runs were performed at 501 grid points.

Shi et al. (2020) only quantify incompressible linear eigenmodes, which they justify by stating that the resistive tearing mode has a negligible contribution due to compressibility, based on the reasoning followed by Furth et al. (1963). We can here easily verify that assumption, using the full compressible functionality of Legolas. It appears that the inclusion of Hall terms in the treatment of the tearing mode causes differences in incompressible vs compressible plasma settings. The influence of compressibility is shown in figure 9(a,b), where the compressible growth rate and frequency, respectively, are shown for $k_2 = 0.155$, to be compared with the incompressible case in figure 7(a,b). Although Furth et al. (1963) showed that compressibility has a negligible effect on the resistive tearing mode, which our tests with Legolas also confirm, their treatment did not take the Hall current into account. When the Hall terms are taken into account, the effect of compressibility on the resistive tearing mode growth rate is no longer negligible. In particular, stronger guide fields result in stronger damping of the growth rate, as evidenced by figure 9. Additionally, new unstable modes appear in the spectrum and become more unstable than the tearing mode for sufficiently large $k_3$. These are Hall instabilities, occurring in a Cartesian slab when the magnetic field is sufficiently curved, i.e. if $\partial ^2\boldsymbol {B}_0/\partial x^2$ is non-zero (Rheinhardt & Geppert 2002). The ranges where the largest Hall instability overtakes the tearing instability as the most unstable mode are indicated in figure 9(a) with lines on the horizontal axis. The part of a spectrum containing the tearing mode and the other unstable modes is shown in figure 9(c).

Figure 9. The real (a) and imaginary (b) parts of the compressible tearing mode in a Harris current sheet (3.20a–c) as a function of $k_3$, with $k_2 = 0.155$, $\tilde {\rho }=1$, $a=1$, $B_0=1$, $\eta = 10^{-4}$ and $\eta _{H} = 1$ for different guide field strengths $B_{g}$, to be compared with the incompressible case in figure 7(a,b). The horizontal lines in (a) indicate the ranges where the tearing mode is not the most unstable mode in the spectrum. (c) Spectrum from the $B_{g} = 5$ series with $k_3 \approx 0.015346$. The tearing mode is circled in red. All runs were performed at $501$ grid points.

4. Conclusion and outlook

In this paper, the extension of the MHD spectroscopy code Legolas (Claes et al. 2020) with viscosity and the Hall current was presented and verified using test cases taken from the literature. To validate the implementation of the viscosity module, we first accurately reproduced the spectrum and eigenfunctions of an incompressible, hydrodynamic Taylor–Couette flow in a cylindrical set-up, taken from Gebhardt & Grossmann (1993), with the newly implemented incompressible approximation. As a second test case, we considered the Cartesian MHD equilibrium with finite resistivity from Dahlburg et al. (1983), where we reproduced their results concerning the interplay of viscous and resistive effects on the growth rate of the resistive tearing instability. As an extension of their results, we showed that the full resistive and viscous spectra are extremely similar, with the prominent distinction that the viscous spectrum does not have an unstable tearing mode. The combination of viscous and resistive effects was mostly seen in the stable part of the spectrum, and its role in further nonlinear evolutions warrants further exploration. However, since both viscosity test cases used the incompressible approximation, which eliminates the energy equation, the viscous heating term did not play a role. For one selected Taylor–Couette case it was shown that the viscous heating did not significantly alter the spectrum, but had a limited influence on the entropy perturbation. The flexible Legolas implementation allows for future linear stability studies with or without viscous heating. Such viscoresistive stability studies of magnetised Taylor–Couette set-ups can be very important for aiding the interpretation of dynamo experiments (Willis & Barenghi 2002; Rüdiger et al. 2007), and especially to determine when the MRI is sufficiently suppressed by viscoresistive effects (Eckhardt & Herron 2018) to create stable configurations.

To test the implementation of the Hall module, another two cases were considered. The simplest case considered an ideal, homogeneous, Cartesian plasma slab with a Hall current. For a small ratio of ion inertial length to plate separation this case is comparable to the infinite, homogeneous medium, described by the dispersion relation of Hameiri et al. (2005). The solutions of the infinite medium corresponded to the first modes in several sequences of modes, as evidenced by the eigenfunctions, which is expected when going from an infinite medium to a semi-infinite medium that is bounded in one direction. Whilst the smallest sequence behaves anti-Sturmian, the larger two display Sturmian behaviour. All three wave sequences become dispersive in Hall-MHD, and the smallest and largest sequences are known as ion cyclotron and whistler modes respectively. The middle sequence fails to capture the electron cyclotron resonance because the electron mass was set to zero. If the electron inertia term is included as well, the electron cyclotron resonance is recovered.

The more advanced test case also included resistivity and evaluated Hall effects on the resistive tearing mode growth rate in a Harris sheet set-up described by Shi et al. (2020). The reproduction of these results required an incompressible approximation, but a good match between both the tearing mode and the eigenfunctions was achieved. However, contrary to the assumption of Shi et al. (2020) that compressibility has a negligible effect on the resistive tearing mode, which was shown for the purely resistive MHD case by Furth et al. (1963), a guide field introduces a non-negligible damping effect when both compressibility and the Hall current are considered.

Our Legolas tool can now be used to combine and explore linear eigenmodes and full eigenspectra for cases where we have multiple effects at play, such as the influence of equilibrium flow or non-uniform equilibrium density, the combination of Hall and viscosity and the electron inertia term. However, the effect of the latter is likely negligible in many cases because the electron inertia coefficient $\eta _{e}$ is several orders of magnitude smaller than the Hall coefficient $\eta _{H}$.

The inclusion of viscosity and the Hall current opens up various research avenues, such as the investigation of the influence of viscosity in resistive set-ups, and in particular how the introduction of viscosity affects resistive instabilities like the resistive tearing mode discussed in § 3.2.2. For the Hall current it is now possible to examine its effect on the previously mentioned MRI in accretion discs (Lesur 2021) or to explore instabilities requiring a Hall current, such as the Hall-shear instability (Kunz 2008). In the context of the MRI, a similar future extension of the Legolas code can implement ambipolar diffusion as a proxy for charge-neutral decoupling effects. This would also introduce the ambipolar-diffusion-shear instability (Kunz 2008).

Looking ahead, this extension brings Legolas one step closer to describing realistic 1-D set-ups. We foresee that the tool can even be used to identify the modes responsible for specific evolutions seen in 2-D and 3-D nonlinear simulations, when we can describe the instantaneous multidimensional MHD fields with a (e.g. height and azimuthally averaged) force-balanced 1-D background state, and which physical effects play a relevant role in observed growth rates. Furthermore, since Legolas captures both fundamental modes and overtones, it can act as a tool in the analysis of observed waves and overtones in coronal loop seismology (see e.g. Andries et al. 2009), albeit in the infinite cylinder approximation, and tokamaks (see e.g. Ochoukov et al. 2018; Spong et al. 2018), although toroidal effects are lost in the cylinder set-up. In the future, vacuum or vacuum-wall boundary conditions could be implemented to better model these physical systems.