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 (Reference Gebhardt and Grossmann1993), 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 (Reference Gebhardt and Grossmann1993) are recovered in figure 1(a–d). These spectra feature two different types of modes, as described in Gebhardt & Grossmann (Reference Gebhardt and Grossmann1993), 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 (Reference Gebhardt and Grossmann1993) ($\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. (Reference Dahlburg, Zang, Montgomery and Hussaini1983), 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. (Reference Dahlburg, Zang, Montgomery and Hussaini1983) 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. (Reference Dahlburg, Zang, Montgomery and Hussaini1983), 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. (Reference Dahlburg, Zang, Montgomery and Hussaini1983) ($\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. (Reference Dahlburg, Zang, Montgomery and Hussaini1983), 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. (Reference Dahlburg, Zang, Montgomery and Hussaini1983) 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. (Reference Dahlburg, Zang, Montgomery and Hussaini1983), 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. Reference Goedbloed, Keppens and Poedts2019) 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. (Reference Goedbloed, Keppens and Poedts2019). 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.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 (Reference Hameiri, Ishizawa and Ishida2005), 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. (Reference Hameiri, Ishizawa and Ishida2005) (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 Reference De Jonghe and Keppens2020), 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. Reference Hameiri, Ishizawa and Ishida2005). 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. (Reference Hameiri, Ishizawa and Ishida2005, 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. (Reference Shi, Velli, Pucci, Tenerani and Innocenti2020), 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. (Reference Shi, Velli, Pucci, Tenerani and Innocenti2020) assume incompressibility, so we also use the incompressible approximation (large $\gamma$) in Legolas.

Shi et al. (Reference Shi, Velli, Pucci, Tenerani and Innocenti2020) 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. (Reference Shi, Velli, Pucci, Tenerani and Innocenti2020). 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. (Reference Shi, Velli, Pucci, Tenerani and Innocenti2020) 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. (Reference Shi, Velli, Pucci, Tenerani and Innocenti2020). 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. (Reference Shi, Velli, Pucci, Tenerani and Innocenti2020) 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. (Reference Furth, Killeen and Rosenbluth1963). 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. (Reference Furth, Killeen and Rosenbluth1963) 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 Reference Rheinhardt and Geppert2002). 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.