Legolas: A Modern Tool for Magnetohydrodynamic Spectroscopy

Niels Claes; Jordi De Jonghe; Rony Keppens

doi:10.3847/1538-4365/abc5c4

1. Introduction

The study of stability for plasmas and fluids alike has been a major topic of research over the last century. Understanding how and why a given medium reacts to a linear perturbation is of central importance to many astrophysical phenomena. In incompressible or compressible fluids, governed by hydrodynamic equations, notable instabilities are the Kelvin–Helmholtz instability (KHI), which arises due to a velocity shear at the interface of two fluids; and the Rayleigh–Taylor instability (RTI), where gravitational stratification can lead to an unstable configuration of layered fluids of different densities (Chandrasekhar 2013; Choudhuri 1998). In plasmas, governed by the magnetohydrodynamic (MHD) equations, the study of waves and instabilities becomes much richer due to the inclusion of magnetic fields, with a modern overview provided in Goedbloed et al. (2019). Magnetic fields modify the two aforementioned instabilities in various ways, and the combination of flow, magnetic fields, and pressure gradients introduces many new modes, e.g., the magnetorotational instability (Balbus & Hawley 1991) relevant for (weakly magnetized) accretion disks or the trans-slow-Alfvén continuum modes in disks of arbitrary magnetization (Goedbloed et al. 2004). In the highly magnetized solar corona, observed coronal loop oscillations (periods and damping times) are routinely used to infer loop parameters like their field strength (Nakariakov & Ofman 2001). Embedded in the hot solar corona, we find stable and long-lived quiescent prominences, with internal dynamics due to KHI (Hillier & Polito 2018) and RTI (Hillier 2018). The formation of prominences is due to the thermal instability (TI), as demonstrated in direct observations by, e.g., Berger et al. (2012) or in simulations by Xia & Keppens (2016) and Claes et al. (2020). Together with categorizing all instabilities, knowing the stable eigenoscillations, such as p modes or g modes in stratified atmospheres or stellar interiors, is of prime importance to link theoretical understanding with observed periodic phenomena. In all of these cases, we need to compute the eigenoscillations and corresponding eigenfunctions from the linearized set of governing equations. Linear MHD spectroscopy, which encompasses the entirety of helio- and asteroseismology but incorporates laboratory fusion plasma MHD spectroscopy (Goedbloed et al. 1993), MHD spectroscopy of accretion disks (Keppens et al. 2002) and jets, as well as solar coronal seismology (Roberts 2019), is thus a powerful tool for studying many astrophysical processes.

Since the advent of more powerful computational resources, the main focus of computational astrophysical research has gradually shifted toward solving the fully nonlinear MHD equations, where many nonadiabatic/nonideal effects are incorporated, depending on the application at hand. While this approach successfully reproduced many physical phenomena, especially for realistic solar setups (e.g., sunspots, Rempel 2012; flares, Ruan et al. 2019; or prominences, Xia & Keppens 2016), it usually fails to answer which specific perturbation produces the complex evolution as witnessed. At the same time, theoretical insight showed that MHD spectral theory actually governs the stability of flowing, (self-)gravitating single-fluid evolutions of nonlinear, time-dependent plasmas, and this at any time during their nonlinear evolution (Demaerel & Keppens 2016). Hence, in order to predict the reaction of a certain physical state to perturbations, we should really quantify all its waves and instabilities using linear theory. This has been recognized fully in laboratory fusion plasmas, where MHD spectroscopy is very successful in identifying waves and stability aspects of a given toroidal Grad–Shafranov equilibrium. That this can meaningfully be done for states that include important nonadiabatic effects, like optically thin radiative losses, is important for investigations into prominences and their intriguing fine structure, as revealed by means of direct observations (Engvold 1998; Ballester 2006; Mackay et al. 2010) or through numerical simulations (Xia & Keppens 2016; Xia et al. 2017; Claes et al. 2020). In that context, early analytical work by Van der Linden & Goossens (1991a) based on linear MHD suggests the hypothesis that finite perpendicular thermal conduction induces fine structure in unstable linear eigenmodes. Since this pioneering work of Van der Linden & Goossens (1991a), not much research has been done regarding the full MHD spectrum when nonadiabatic effects are at play, for the simple reason that to date, there existed no numerical tool to solve the full system of linearized MHD equations with all the physical effects included. This is why we developed the new and open-source Legolas solver.

Legolas builds on the heritage of early numerical codes, most notably LEDA (Kerner et al. 1985), which allowed studies of the ideal or resistive MHD spectrum for laboratory plasmas, approximated by a diffuse cylindrical plasma column (or flux tube) and CASTOR (Kerner et al. 1998), which applied to resistive spectra of general tokamak configurations. The latter has follow-up codes such as FINESSE (Beliën et al. 2002) and PHOENIX (Blokland et al. 2007b), extending it to stationary and axisymmetric truly 2D configurations. LEDA was later extended in Van der Linden et al. (1992), where nonadiabatic effects like anisotropic thermal conduction and optically thin radiative losses were added to the equations using a simple analytic function to treat radiative cooling effects. A different branch of LEDA, called LEDAFLOW (Nijboer et al. 1997), was developed to investigate the resistive MHD spectrum, augmented with gravitational and flow effects, but omitting those nonadiabatic terms. Because these codes were developed decades ago and focus has shifted away from linear MHD, their further development was stalled, although in laboratory fusion context, tools to compute multidimensional equilibria and their linear modes are very important for diagnosing experiments. The original codes, like LEDA(FLOW), were not flexible in the sense that adding different equilibria or accounting for additional terms in the equations would be a major undertaking, as parts were hard-coded to (limited) computational resources of that time. Furthermore, programming languages and numerical tools like LAPACK (Anderson et al. 1999) to solve eigenvalue problems have come a long way. This prompted us to develop a brand new, modern MHD spectral code which we named Legolas, short for "Large Eigensystem Generator for One-dimensional pLASmas." The Legolas code is able to handle both Cartesian and cylindrical geometries, and introduces many new features, e.g., selecting between modern cooling curves that treat optically thin radiative cooling effects. Furthermore, every aspect of the code is modularized, making it ready to be extended with additional physics or modern algorithmic requirements (such as mesh refinement). The main goal of this paper is to present the new code in terms of its implementation details and to validate it against a plethora of test cases that ensure a correct treatment of the governing equations.

These tests include eigenmode quantifications of ideal, static MHD configurations under adiabatic conditions, where the static (that is, no equilibrium flow) and adiabatic linear MHD equations make the problem self-adjoint. When performing a standard Fourier analysis in the ignorable directions, the resulting eigenvalue problem is then Hermitian, meaning that all eigenfrequencies will be either fully real (stable waves) or fully complex (pure damped or unstable modes), hence they are found on the real or imaginary axis of the complex eigenfrequency plane, and the full MHD spectrum will be both left–right and up–down symmetric. However, in nature, physical conditions may be far from ideal. The inclusion of nonideal effects like resistivity or thermal conduction lifts the self-adjointness of the eigenvalue problem, allowing the eigenmodes to move away from the axes into the complex plane, and the up–down symmetry gets broken. As long as the equilibrium configuration is static, all (adiabatic or nonadiabatic) modes will still have a complementary mode that lies mirrored around the imaginary axis, making the entire spectrum left–right symmetric. This is related to the forward and backward propagating mode symmetry, or the equivalent statement on the parity–time (PT) symmetry. However, for typical astrophysical plasmas, the conditions are far from static: tokamak plasmas, astrophysical jets, solar coronal loops, accretion disks, etc., all have equilibrium flows. The inclusion of a background flow breaks the left–right symmetry of the MHD spectrum, resulting in an even more complicated structure. However, the study of the ideal, linear MHD spectrum of flowing plasmas is still governed by a pair of self-adjoint operators (Goedbloed 2011; Goedbloed et al. 2019), and it leaves the up–down symmetry of the spectrum intact (where every overstable mode has an equivalent damped counterpart at the same frequency). The combination of flow and nonadiabatic effects, where both left–right and up–down eigenfrequency symmetries are broken, has never been explored in earnest. All of the above makes it clear that a numerical approach becomes essential, especially when the equilibrium is no longer homogeneous. Because in reality, virtually no astrophysical configuration is spatially homogeneous, we need a flexible numerical tool to explore the spectrum systematically.

Legolas solves the linearized MHD equations including various nonadiabatic effects, resistivity, and gravity, assuming a one-dimensional (1D) equilibrium profile with the possibility of background flow. A standard Fourier analysis for the perturbations is combined with a finite element representation of the eigenfunctions in the important coordinate. We transform the original system into an eigensystem for the complex eigenfrequencies ω. We use a general formalism to include two kinds of geometries, a plane Cartesian stratified slab or a (possibly also stratified) cylinder, through the inclusion of a scale factor originating from the divergence, gradient, and curl operators. Legolas can handle the hydrodynamic limit (where all equilibrium magnetic field components are set to zero), enabling us to investigate the stability of hydrodynamic static and stationary equilibria. The resulting system of equations is solved in weak form, transforming the original system into a non-Hermitian complex eigenvalue problem, which is solved using the QR algorithm. This results in a calculation of all eigenfrequencies and corresponding eigenfunctions of the system, such that a detailed analysis of mode stability, but also the entire overview on all supported linear wave modes, becomes possible.

Section 2 introduces the system of equations, along with the linearization procedure and Fourier mode representation. The treatment of the final system using the finite element method is given in the Appendix where we explain the basic mathematical formalism behind the FEM, with a complete treatment of the finite element matrix assembly process and the boundary conditions. Legolas is tested against a large amount of spectra found in the literature in Section 3, which is subdivided into multiple categories. First, we treat ideal MHD with only a gravitational term included, which has as main advantage that solutions can be obtained analytically. The first case handles gravito-MHD waves in a Cartesian slab, after which we quantify quasi-Parker instabilities in stratified atmospheres. Results obtained with Legolas are compared with results found in various modern textbooks, including results for cylindrical geometries, discussing modes for ideal flux tubes, and tokamak current equilibria where we show liability to interchanges. For the inclusion of flow into the equations, we consider a nontrivial case related to astrophysical jet stability, looking at Kelvin–Helmholtz and current-driven instabilities. We demonstrate that we can compute Suydam cluster modes, originating from a surface where the wavevector is locally perpendicular to the magnetic field. We then transition to the inclusion of nonideal effects such as resistivity, looking first at the resistive spectrum for a homogeneous plasma. We then extend to recover the resistive quasi-mode in a nonhomogeneous case. This is further complemented by adding an inhomogeneous medium and background flow in such a way as to give rise to resistive tearing modes, all of which are in a Cartesian geometry. Additionally, by allowing for resistivity and current variation, we show that we can also compute resistive rippling modes. Lastly, we treat nonadiabatic effects, including optically thin radiative losses and thermal conduction into the equations, revisiting some pioneering results of discrete Alfvén waves and magnetothermal instabilities. Along the way, we find interesting extensions to the original published works, due to the much higher resolutions employed here. Because Legolas is the first modern linear MHD code to investigate realistic astrophysical plasmas, this opens the door to further in-depth studies of nonideal equilibrium configurations at high resolutions, ranging from loops, jets, or accretion disks. The Legolas code is available on GitHub at https://github.com/n-claes/legolas.

2. Problem Description and Model Equations

The MHD equations with nonadiabatic effects, resistivity, and gravity included can be written in a (normalized) Eulerian representation as

$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}),\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&\rho \displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}=-{\rm{\nabla }}p-\rho {\boldsymbol{v}}\cdot {\rm{\nabla }}{\boldsymbol{v}}+({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}+\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}\begin{array}{rcl}\rho \displaystyle \frac{\partial T}{\partial t} & = & -\rho {\boldsymbol{v}}\cdot {\rm{\nabla }}T-(\gamma -1)p{\rm{\nabla }}\cdot {\boldsymbol{v}}-(\gamma -1)\rho {\mathscr{L}}\\ & & +\ (\gamma -1){\rm{\nabla }}\cdot ({\boldsymbol{\kappa }}\cdot {\rm{\nabla }}T)+(\gamma -1)\eta {\left({\rm{\nabla }}\times {\boldsymbol{B}}\right)}^{2},\end{array}\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}})-{\rm{\nabla }}\times (\eta {\rm{\nabla }}\times {\boldsymbol{B}}),\end{eqnarray} \tag{ 4 }$

where ρ is the plasma density, ${\boldsymbol{v}}$ is the velocity, T is the temperature, p is the pressure, ${\boldsymbol{B}}$ is the magnetic field (satisfying ${\rm{\nabla }}\cdot {\boldsymbol{B}}=0$ ), η is the resistivity, and ${\boldsymbol{g}}$ the gravitational acceleration. To close the system, the (normalized) ideal gas law p = ρT is used, while γ denotes the ratio of specific heats, taken to be 5/3. The symbol ${\mathscr{L}}$ in Equation (3) represents the heat-loss function, defined as energy losses minus energy gains due to optically thin radiative cooling effects (Parker 1953) and is given by

$\begin{eqnarray}&&{\mathscr{L}}=\rho {\rm{\Lambda }}(T)-{ \mathcal H },\end{eqnarray} \tag{ 5 }$

where ${ \mathcal H }$ represents the total energy gains. In general, ${ \mathcal H }$ can be anything (for example, heating through dissipative Alfvén waves; van der Holst et al. 2014), but because there is still no well-defined parameterization for coronal heating to date, this term is assumed to be as convenient as possible, that is, constant in time but possibly varying in space to ensure thermodynamic balance. Note that ${ \mathcal H }$ should always be consistent with a given equilibrium profile as to exactly balance out the radiative losses (and possibly the thermal conduction effects and/or ohmic heating effects), reaching a thermal equilibrium state. This indirectly implies that the heating term is not necessarily independent of location, but dependent on the connection between the radiative losses and equilibrium temperature profile, which, in general, are both spatially dependent. The first term in Equation (5) denotes the radiative losses, dependent on the cooling curve ${\rm{\Lambda }}(T)$ . These curves are tabulated sets resulting from detailed calculations, and can hence be interpolated to high temperature resolutions. Legolas has multiple cooling curves implemented, most notably those by Colgan et al. (2008) and Schure et al. (2009), where the latter is extended to the low-temperature limit using Dalgarno & McCray (1972). In addition, we also implemented a piecewise power law as described by Rosner et al. (1978), which is an explicit (piecewise) function over the entire temperature domain. In the solar and astrophysical literature, these cooling curves collect detailed knowledge on radiative processes, which are all assumed to be in the optically thin regime. It is worth noting that the inclusion of time-dependent background heating or flows can influence the spectrum in itself, as shown in, for example, Barbulescu et al. (2019) and Hillier et al. (2019), but due to the time dependence of those effects, the assumption of a stationary background state is no longer applicable. A possibility however may be varying the background heating or flow in subsequent runs, provided the background state is known at every snapshot.

Thermal conduction in magnetized plasmas is highly anisotropic, as the effect is a few orders of magnitude stronger along the field lines than across. We hence use a tensor representation to model this anisotropy, denoting the thermal conductivity tensor ${\boldsymbol{\kappa }}$ by

$\begin{eqnarray}&&{\boldsymbol{\kappa }}={\kappa }_{\parallel }{{\boldsymbol{e}}}_{B}{{\boldsymbol{e}}}_{B}+{\kappa }_{\perp }({\boldsymbol{I}}-{{\boldsymbol{e}}}_{B}{{\boldsymbol{e}}}_{B}),\end{eqnarray} \tag{ 6 }$

where ${\boldsymbol{I}}$ denotes the unit tensor and ${{\boldsymbol{e}}}_{B}={\boldsymbol{B}}/B$ is a unit vector along the magnetic field. The coefficients κ_∥ and κ_⊥ denote the conductivity coefficients parallel and perpendicular to the local direction of the magnetic field. For typical astrophysical applications, the Spitzer conductivity is used, given by

$\begin{eqnarray}\begin{array}{rcl}{\kappa }_{\parallel } & \approx & 8\times {10}^{-7}{T}^{5/2}\,\mathrm{erg}\,{\mathrm{cm}}^{-1}\,{{\rm{s}}}^{-1}\,{{\rm{K}}}^{-1},\\ {\kappa }_{\perp } & \approx & 4\times {10}^{-10}{n}^{2}{B}^{-2}{T}^{-3}{\kappa }_{\parallel },\\ {\kappa }_{\parallel } & \approx & 8\times {10}^{-12}{T}^{5/2}\,{\rm{W}}\,{{\rm{m}}}^{-1}\,{{\rm{K}}}^{-1},\\ {\kappa }_{\perp } & \approx & 4\times {10}^{-30}{n}^{2}{B}^{-2}{T}^{-3}{\kappa }_{\parallel },\end{array}\end{eqnarray} \tag{ 7 }$

where n denotes the number density, given by ρ = nm_p with m_p the proton mass. The first and second rows in Equation (7) give the thermal conduction coefficients in cgs and mks units, respectively. For the solar corona, we typically find that κ_∥ is about 12 orders of magnitude larger than κ_⊥ (Priest 2014), so perpendicular thermal conduction is usually ignored. Nevertheless, in Legolas, both parallel and perpendicular thermal conduction are implemented. For the resistivity η, one can in principle take any profile. We implemented the Spitzer resistivity,

$\begin{eqnarray}&&\eta =\displaystyle \frac{4\sqrt{2\pi }}{3}\displaystyle \frac{{Z}_{\mathrm{ion}}{e}^{2}\sqrt{{m}_{e}}\mathrm{ln}(\lambda )}{{\left(4\pi {\epsilon }_{0}\right)}^{2}{\left({k}_{{\rm{B}}}T\right)}^{3/2}},\end{eqnarray} \tag{ 8 }$

where Z_ion denotes the ionization taken to be unity, e and m_e denote the electron charge and mass, respectively, and ₀ and k_B are the electrical permittivity and Boltzmann constant. The Coulomb logarithm is given by $\mathrm{ln}(\lambda )$ and is approximately equal to 22 for solar coronal conditions (Goedbloed et al. 2019). It is important to emphasize that, because we will further linearize the governing nonlinear equations, we can adopt fully realistic values for all the nonideal coefficients, such as the resistivity or thermal conduction coefficients. This is in contrast to fully nonlinear computations, which are severely restrained in reaching magnetic Reynolds numbers beyond 10⁴–10⁵.

**Figure 1.** Unit vectors and examples of ${{\boldsymbol{B}}}_{0}$ and ${{\boldsymbol{v}}}_{0}$ for the Cartesian (left) and cylindrical case (right).
Download figure:
Standard image High-resolution image

**Figure 1.** Unit vectors and examples of ${{\boldsymbol{B}}}_{0}$ and ${{\boldsymbol{v}}}_{0}$ for the Cartesian (left) and cylindrical case (right).
Download figure:
Standard image High-resolution image

2.1. Equilibrium Conditions

We consider a general coordinate system denoted by (u₁, u₂, u₃), corresponding to three orthogonal basis vectors. The main advantage of this approach is that it allows us to include two different geometries with only one basic formalism (and implementation). First, we consider a standard plane slab geometry in Cartesian coordinates, that is, a plasma which is confined in height, and considered to be bounded by two horizontal, perfectly conducting walls at a fixed distance apart, extending outwards to infinity in the other two ignorable coordinates. This case also approximates the limit of a fully infinite free space when the walls are moved off to infinity. In Cartesian geometry, the coordinate system can be written as (x, y, z) and the vectors { u ₁, u ₂, u ₃} are the standard Cartesian triad $\{{\hat{{\boldsymbol{e}}}}_{x},{\hat{{\boldsymbol{e}}}}_{y},{\hat{{\boldsymbol{e}}}}_{z}\}$ along the axes. This makes it quite convenient to include, for example, gravitational effects that will induce an equilibrium stratification in the u₁ coordinate. The second geometry is that of an infinitely long plasma cylinder encased by a solid wall at a certain distance away from the cylinder axis, for which the coordinate system can be defined as (r, θ, z). At each point, the vectors { u ₁, u ₂, u ₃} are defined as the triad of tangent vectors, $\{{\hat{{\boldsymbol{e}}}}_{r},{\hat{{\boldsymbol{e}}}}_{\theta },{\hat{{\boldsymbol{e}}}}_{z}\}$ , with ${\hat{{\boldsymbol{e}}}}_{r}$ along the radial direction, ${\hat{{\boldsymbol{e}}}}_{z}$ in the direction of the cylinder axis, and ${\hat{{\boldsymbol{e}}}}_{\theta }={\hat{{\boldsymbol{e}}}}_{z}\times {\hat{{\boldsymbol{e}}}}_{r}$ tangent to the cylinder. A detailed view of both geometries and their corresponding coordinate systems is shown in Figure 1. The basic operators present in Equations (1)–(4), that is, the divergence, gradient, and curl, introduce a scale factor $\varepsilon =r$ for cylindrical geometries, which is reduced to $\varepsilon =1$ for a Cartesian coordinate system. Hence, exploiting this scale factor in the mathematical formalism allows for one implementation, where one can conveniently switch between both cases. We note that the cylindrical setup is also applicable to the so-called cylindrical accretion disk limit, as, for example, exploited to study MHD instabilities in disks by Blokland et al. (2007a).

Linearization involves splitting variables into two parts: a time-independent part, usually denoted by subscript 0, and a perturbed part, denoted by a subscript 1. Legolas handles one-dimensional equilibria which depend only on u₁, or, more specifically, time-independent equilibria of the form

$\begin{eqnarray}\begin{array}{rcl}{\rho }_{0} & = & {\rho }_{0}({u}_{1}),\quad {{\boldsymbol{v}}}_{0}={v}_{02}({u}_{1}){{\boldsymbol{u}}}_{2}+{v}_{03}({u}_{1}){{\boldsymbol{u}}}_{3},\\ {p}_{0} & = & {p}_{0}({u}_{1}),\quad {{\boldsymbol{B}}}_{0}={B}_{02}({u}_{1}){{\boldsymbol{u}}}_{2}+{B}_{03}({u}_{1}){{\boldsymbol{u}}}_{3}.\\ {T}_{0} & = & {T}_{0}({u}_{1}),\end{array}\end{eqnarray} \tag{ 9 }$

In general, we have ${\boldsymbol{g}}=-g({u}_{1}){{\boldsymbol{u}}}_{1}$ , where the Cartesian case is for a stratified atmosphere or layer, and the cylindrical case can also allow for gravitational stratification of an accretion disk situated for u₁ = r ∈ [1, R]. In the case of a cylinder where u₁ = r ∈ [0, R], this gravitational term is absent. Using these equations in combination with Equations (1)–(4) yields two conditions for the time-independent parts, given by

$\begin{eqnarray}&&\left({p}_{0}+\displaystyle \frac{1}{2}{{\boldsymbol{B}}}_{0}^{2}\right)^{\prime} +{\rho }_{0}g-\displaystyle \frac{\varepsilon ^{\prime} }{\varepsilon }\left({\rho }_{0}{v}_{02}^{2}-{B}_{02}^{2}\right)=0,\end{eqnarray} \tag{ 10 }$

$\begin{eqnarray}&&\displaystyle \frac{1}{\varepsilon }\left(\varepsilon {\kappa }_{\perp }{T}_{0}^{{\prime} }\right)^{\prime} -{\rho }_{0}{{\mathscr{L}}}_{0}=0,\end{eqnarray} \tag{ 11 }$

where the first condition originates from the momentum Equation (2) and should always be satisfied as it expresses a force-balanced state. The second condition originates from the nonadiabatic terms in the energy equation and should be accounted for if these terms are included; the prime denotes the derivative with respect to u₁. It should be noted that resistive terms are not considered here, which is justified by considering that the timescales on which the magnetic fields decay due to resistivity is much, much larger than the timescales of resistive modes. This is a consequence of large magnetic Reynolds numbers R_m in typical astrophysical cases, yielding magnetic decay timescales of τ ∼ R_m τ_A (with τ_A the typical Alfvén time in ideal MHD) compared to the much faster resistive mode timescales of $\tau \sim {R}_{m}^{\alpha }$ (where typically 0 < α < 1). We can hence consider the equilibrium itself to be independent of resistivity, which removes some stringent extra conditions on the energy and induction equations. Also note that the third term in Equation (10) is only included for a cylinder, as $\varepsilon ^{\prime} =0$ in a Cartesian geometry. This translates to the well-known fact that the centrifugal and tensional parts of the Lorentz force are absent for a Cartesian slab. Furthermore, a cylindrical equilibrium profile should satisfy on-axis regularity conditions, meaning that ${v}_{02},{v}_{03}^{{\prime} },{B}_{02}$ , and ${B}_{03}^{{\prime} }$ all have to be equal to zero at r = 0. When considering an accretion disk in the cylindrical limit, the inner edge of the disk is at r = 1, so lengths are then expressed in this inner disk radius and no regularity conditions apply then.

2.2. Linearized Equations

Now, we linearize Equations (1)–(4) around the equilibrium specified in Equation (9), where the unperturbed time-independent parts are denoted by a subscript 0 and the perturbed time-dependent parts are denoted by a subscript 1. It follows from the adopted equilibrium configuration that ${\rm{\nabla }}\cdot {{\boldsymbol{v}}}_{0}=0$ and ${\rm{\nabla }}\cdot {{\boldsymbol{B}}}_{0}=0$ , such that the divergence-free condition on the magnetic field is fulfilled and that the equilibrium flow field is incompressible. However, the perturbed quantities can represent both incompressible or compressible eigenoscillations. The no-monopole condition should also be taken into account for the perturbed magnetic field. Therefore, we adopt a vector potential to write ${{\boldsymbol{B}}}_{1}={\rm{\nabla }}\times {{\boldsymbol{A}}}_{1}$ such that ${\rm{\nabla }}\cdot {{\boldsymbol{B}}}_{1}=0$ is automatically satisfied. The system of linearized equations is thus given by

$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{1}}{\partial t}=-{\rho }_{0}{\rm{\nabla }}\cdot {{\boldsymbol{v}}}_{1}-{{\boldsymbol{v}}}_{0}\cdot {\rm{\nabla }}{\rho }_{1},\end{eqnarray} \tag{ 12 }$

$\begin{eqnarray}\begin{array}{rcl}{\rho }_{0}\displaystyle \frac{\partial {{\boldsymbol{v}}}_{1}}{\partial t} & = & -{\rho }_{0}{{\boldsymbol{v}}}_{0}\cdot {\rm{\nabla }}{{\boldsymbol{v}}}_{1}-{\rho }_{0}{{\boldsymbol{v}}}_{1}\cdot {\rm{\nabla }}{{\boldsymbol{v}}}_{0}-{\rho }_{1}{{\boldsymbol{v}}}_{0}\cdot {\rm{\nabla }}{{\boldsymbol{v}}}_{0}\\ & & -\ {\rm{\nabla }}\left({\rho }_{1}{T}_{0}+{\rho }_{0}{T}_{1}\right)+({\rm{\nabla }}\times {{\boldsymbol{B}}}_{0})\times ({\rm{\nabla }}\times {{\boldsymbol{A}}}_{1})\\ & & +\ \left[{\rm{\nabla }}\times ({\rm{\nabla }}\times {{\boldsymbol{A}}}_{1})\right]\times {{\boldsymbol{B}}}_{0}+{\rho }_{1}{\boldsymbol{g}},\end{array}\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}&&\begin{array}{l}{\rho }_{0}\displaystyle \frac{\partial {T}_{1}}{\partial t}=-{\rho }_{0}{{\boldsymbol{v}}}_{1}\cdot {\rm{\nabla }}{T}_{0}-{\rho }_{0}{{\boldsymbol{v}}}_{0}\cdot {\rm{\nabla }}{T}_{1}-(\gamma -1){\rho }_{0}{T}_{0}{\rm{\nabla }}\cdot {{\boldsymbol{v}}}_{1}\\ \quad \ -\ (\gamma -1){\rho }_{1}{{\mathscr{L}}}_{0}-(\gamma -1){\rho }_{0}\left({{\mathscr{L}}}_{T}{T}_{1}+{{\mathscr{L}}}_{\rho }{\rho }_{1}\right)\\ \quad \ +\ (\gamma -1){\rm{\nabla }}\cdot ({{\boldsymbol{\kappa }}}_{0}\cdot {\rm{\nabla }}{T}_{1})+(\gamma -1){\rm{\nabla }}\cdot ({{\boldsymbol{\kappa }}}_{1}\cdot {\rm{\nabla }}{T}_{0})\\ \quad \ +\ 2(\gamma -1)\eta ({\rm{\nabla }}\times {{\boldsymbol{B}}}_{0})\cdot \left[{\rm{\nabla }}\times ({\rm{\nabla }}\times {{\boldsymbol{A}}}_{1})\right]\\ \quad \ +\ (\gamma -1)\displaystyle \frac{d\eta }{{dT}}{T}_{1}{\left({\rm{\nabla }}\times {{\boldsymbol{B}}}_{0}\right)}^{2},\end{array}\end{eqnarray} \tag{ 14 }$

$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {{\boldsymbol{A}}}_{1}}{\partial t} & = & {{\boldsymbol{v}}}_{1}\times {{\boldsymbol{B}}}_{0}+{{\boldsymbol{v}}}_{0}\times ({\rm{\nabla }}\times {{\boldsymbol{A}}}_{1})-\eta {\rm{\nabla }}\\ & & \times \ ({\rm{\nabla }}\times {{\boldsymbol{A}}}_{1})-\displaystyle \frac{d\eta }{{dT}}{T}_{1}({\rm{\nabla }}\times {{\boldsymbol{B}}}_{0}),\end{array}\end{eqnarray} \tag{ 15 }$

where p₁ is replaced by ρ₁ T₀ + ρ₀ T₁ resulting from the linearized ideal gas law. The perturbation ${{\boldsymbol{\kappa }}}_{1}$ of the thermal conduction tensor is obtained by linearizing the expressions (6)–(7), while the derivatives of the heat-loss function with respect to density and temperature are given by

$\begin{eqnarray}&&{{\mathscr{L}}}_{\rho }={\left.\displaystyle \frac{\partial {\mathscr{L}}}{\partial \rho }\right|}_{T}={\rm{\Lambda }}({T}_{0}),\quad \ {{\mathscr{L}}}_{T}={\left.\displaystyle \frac{\partial {\mathscr{L}}}{\partial T}\right|}_{\rho }={\rho }_{0}\displaystyle \frac{d{\rm{\Lambda }}({T}_{0})}{{dT}},\end{eqnarray} \tag{ 16 }$

which should be evaluated using the equilibrium quantities. The terms containing ∂η/∂T follow from a linearization of the resistivity parameter, which can be written in terms of the variable T₁ using the temperature dependence originating from the assumed Spitzer resistivity in Equation (8), that is, η₁ = T₁∂η/∂T. In addition, Legolas allows for an anomalous resistivity prescription in which we typically have η( u ₁, j ( u ₁, t)), that is, a resistivity profile that is spatiotemporal in general, but for a fixed time, depends on position and current profile. This in turn implies that a total derivative should be used for the resistivity, given by $\eta ^{\prime} =\partial \eta /\partial {u}_{1}+{T}_{0}^{{\prime} }\partial \eta /\partial T$ . In most use cases a resistivity profile η(T( u ₁)) is sufficient, which is only temperature dependent (and hence indirectly spatially varying as well for inhomogeneous temperature profiles). Also note that these linear equations (as also the nonlinear set above) assumed an external gravitational field, so we did not need to linearize the gravity term ${\boldsymbol{g}}$ (this is the so-called Cowling approximation). In the future, we can extend the set of equations with the Poisson equation and also allow for self-gravity-driven Jeans instabilities.

Next, we perform a Fourier analysis with an exponential time dependence, imposing standard Fourier modes on the u₂ and u₃ coordinates of a perturbed quantity f₁, given by

$\begin{eqnarray}&&{f}_{1}=\widehat{{f}_{1}}({u}_{1})\exp \left[i({k}_{2}{u}_{2}+{k}_{3}{u}_{3}-\omega t)\right].\end{eqnarray} \tag{ 17 }$

In this form, the wavenumbers k₂ and k₃ correspond to k_y and k_z in Cartesian geometry, and to m and k in cylindrical geometry, respectively. Note that m is quantified to integer values, because the θ direction is periodic. Additionally, we apply the following transformation to the perturbed quantities:

$\begin{eqnarray}\begin{array}{rcl}\varepsilon {\widehat{\rho }}_{1} & = & {\widetilde{\rho }}_{1},\quad i\varepsilon {\widehat{v}}_{1}={\widetilde{v}}_{1},\quad {\widehat{v}}_{2}={\widetilde{v}}_{2},\quad \varepsilon {\widehat{v}}_{3}={\widetilde{v}}_{3},\\ \varepsilon {\widehat{T}}_{1} & = & {\widetilde{T}}_{1},\quad i{\widehat{A}}_{1}={\widetilde{a}}_{1},\quad \varepsilon {\widehat{A}}_{2}={\widetilde{a}}_{2},\quad {\widehat{A}}_{3}={\widetilde{a}}_{3},\end{array}\end{eqnarray} \tag{ 18 }$

This particular transformation simplifies the resulting set of equations and has as an additional effect that all terms are real except for the nonadiabatic and resistive contributions, such that we are only dealing with imaginary terms when these physical effects are included. This is in analogy to the fact that the purely adiabatic case is governed by self-adjoint operators: one for the case without flow and two for the case with flow included (Goedbloed 2018a, 2018b). The final set of Fourier-analyzed linearized equations is given below, where the tilde notation in Equations (18) is dropped for the sake of simplicity. From now on, tildes will no longer be written explicitly because there is no confusion possible:

$\begin{eqnarray}\begin{array}{rcl}\omega \displaystyle \frac{1}{\varepsilon }{\rho }_{1} & = & -\displaystyle \frac{1}{\varepsilon }{\rho }_{0}^{{\prime} }{v}_{1}-\displaystyle \frac{1}{\varepsilon }{\rho }_{0}\left({v}_{1}^{{\prime} }-{k}_{2}{v}_{2}-{k}_{3}{v}_{3}\right)\\ & & +\ \displaystyle \frac{1}{\varepsilon }\left(\displaystyle \frac{1}{\varepsilon }{k}_{2}{v}_{02}+{k}_{3}{v}_{03}\right){\rho }_{1},\end{array}\end{eqnarray} \tag{ 19 }$

$\begin{eqnarray}&&\begin{array}{l}\omega {\rho }_{0}\displaystyle \frac{1}{\varepsilon }{v}_{1}=\left(\displaystyle \frac{{\rho }_{1}{T}_{0}+{\rho }_{0}{T}_{1}}{\varepsilon }\right)^{\prime} +\displaystyle \frac{1}{\varepsilon }{B}_{03}^{{\prime} }\left({a}_{2}^{{\prime} }-{k}_{2}{a}_{1}\right)\\ \quad \ -\ \displaystyle \frac{(\varepsilon {B}_{02})^{\prime} }{\varepsilon }\left({a}_{3}^{{\prime} }-{k}_{3}{a}_{1}\right)\\ \quad \ -\ {B}_{03}\left\{\displaystyle \frac{{k}_{3}^{2}}{\varepsilon }{a}_{2}-\displaystyle \frac{{k}_{2}{k}_{3}}{\varepsilon }{a}_{3}-\left[\displaystyle \frac{1}{\varepsilon }\left({a}_{2}^{{\prime} }-{k}_{2}{a}_{1}\right)\right]^{\prime} \right\}\\ \quad \ +\ \displaystyle \frac{1}{\varepsilon }\left(\displaystyle \frac{{k}_{2}{v}_{02}}{\varepsilon }+{k}_{3}{v}_{03}\right){\rho }_{0}{v}_{1}\\ \quad \ +\ {B}_{02}\left\{\displaystyle \frac{{k}_{2}^{2}}{{\varepsilon }^{2}}{a}_{3}-\displaystyle \frac{{k}_{2}{k}_{3}}{{\varepsilon }^{2}}{a}_{2}-\displaystyle \frac{1}{\varepsilon }\left[\varepsilon \left({a}_{3}^{{\prime} }-{k}_{3}{a}_{1}\right)\right]^{\prime} \right\}\\ \quad \ -\ 2\displaystyle \frac{\varepsilon ^{\prime} }{\varepsilon }{\rho }_{0}{v}_{02}{v}_{2}-\displaystyle \frac{\varepsilon ^{\prime} }{{\varepsilon }^{2}}{v}_{02}^{2}{\rho }_{1}+\displaystyle \frac{1}{\varepsilon }g{\rho }_{1},\end{array}\end{eqnarray} \tag{ 20 }$

$\begin{eqnarray}&&\begin{array}{l}\omega {\rho }_{0}{v}_{2}=\displaystyle \frac{1}{{\varepsilon }^{2}}\left({\rho }_{1}{T}_{0}+{\rho }_{0}{T}_{1}\right){k}_{2}+{\rho }_{0}\left(\displaystyle \frac{{k}_{2}{v}_{02}}{\varepsilon }+{k}_{3}{v}_{03}\right){v}_{2}\\ \quad \ +\ \displaystyle \frac{1}{{\varepsilon }^{2}}(\varepsilon {B}_{02})^{\prime} ({k}_{3}{a}_{2}-{k}_{2}{a}_{3})\\ \ +\ {B}_{03}\left[-\left({k}_{3}^{2}+\displaystyle \frac{{k}_{2}^{2}}{{\varepsilon }^{2}}\right){a}_{1}+\displaystyle \frac{{k}_{2}}{{\varepsilon }^{2}}{a}_{2}^{{\prime} }+{k}_{3}{a}_{3}^{{\prime} }\right]-\displaystyle \frac{1}{{\varepsilon }^{2}}{\rho }_{0}(\varepsilon {v}_{02})^{\prime} {v}_{1},\end{array}\end{eqnarray} \tag{ 21 }$

$\begin{eqnarray}&&\begin{array}{l}\omega {\rho }_{0}\displaystyle \frac{1}{\varepsilon }{v}_{3}=\displaystyle \frac{1}{\varepsilon }({\rho }_{1}{T}_{0}+{\rho }_{0}{T}_{1}){k}_{3}+\displaystyle \frac{1}{\varepsilon }{\rho }_{0}\left(\displaystyle \frac{{k}_{2}{v}_{02}}{\varepsilon }+{k}_{3}{v}_{03}\right){v}_{3}\\ \quad \ +\displaystyle \frac{1}{\varepsilon }{B}_{03}^{{\prime} }({k}_{3}{a}_{2}-{k}_{2}{a}_{3})\\ \ -\ {B}_{02}\left[-\left({k}_{3}^{2}+\displaystyle \frac{{k}_{2}^{2}}{{\varepsilon }^{2}}\right){a}_{1}+\displaystyle \frac{{k}_{2}}{{\varepsilon }^{2}}{a}_{2}^{{\prime} }+{k}_{3}{a}_{3}^{{\prime} }\right]-\displaystyle \frac{1}{\varepsilon }{\rho }_{0}{v}_{03}^{{\prime} }{v}_{1},\end{array}\end{eqnarray} \tag{ 22 }$

$\begin{eqnarray}&&\begin{array}{l}\omega {\rho }_{0}\displaystyle \frac{1}{\varepsilon }{T}_{1}=-\displaystyle \frac{1}{\varepsilon }{\rho }_{0}{T}_{0}^{{\prime} }{v}_{1}+\displaystyle \frac{1}{\varepsilon }{\rho }_{0}\left(\displaystyle \frac{{k}_{2}{v}_{02}}{\varepsilon }+{k}_{3}{v}_{03}\right){T}_{1}\\ \quad \ -\ (\gamma -1)\displaystyle \frac{1}{\varepsilon }{\rho }_{0}{T}_{0}\left({v}_{1}^{{\prime} }-{k}_{2}{v}_{2}-{k}_{3}{v}_{3}\right)\\ \ -\ i(\gamma -1)\displaystyle \frac{\left({\kappa }_{\parallel ,0}-{\kappa }_{\perp ,0}\right)}{\varepsilon }\displaystyle \frac{1}{{B}^{2}}{\left(\displaystyle \frac{{k}_{2}{B}_{02}}{\varepsilon }+{k}_{3}{B}_{03}\right)}^{2}{T}_{1}\\ \quad \ +\ i(\gamma -1)\displaystyle \frac{1}{\varepsilon }\left[\varepsilon {\kappa }_{\perp ,0}\left(\displaystyle \frac{{T}_{1}}{\varepsilon }\right)^{\prime} \right]^{\prime} \\ \ -\ i(\gamma -1)\displaystyle \frac{{\kappa }_{\perp ,0}}{\varepsilon }\left(\displaystyle \frac{{k}_{2}^{2}}{{\varepsilon }^{2}}+{k}_{3}^{2}\right){T}_{1}+i(\gamma -1)\displaystyle \frac{1}{\varepsilon }\left(\varepsilon {\kappa }_{\perp ,1}{T}_{0}^{{\prime} }\right)^{\prime} \\ \ +\ i(\gamma -1)\displaystyle \frac{\left({\kappa }_{\parallel ,0}-{\kappa }_{\perp ,0}\right)}{\varepsilon }\displaystyle \frac{1}{{B}^{2}}{T}_{0}^{{\prime} }\left[\left(\displaystyle \frac{{k}_{2}{k}_{3}}{\varepsilon }{B}_{02}+{k}_{3}^{2}{B}_{03}\right){a}_{2}\right.\\ \quad \ -\ \left.\left(\displaystyle \frac{{k}_{2}^{2}}{\varepsilon }{B}_{02}+{k}_{2}{k}_{3}{B}_{03}\right){a}_{3}\right]\\ \ -\ i(\gamma -1){{\mathscr{L}}}_{0}\displaystyle \frac{1}{\varepsilon }{\rho }_{1}-i(\gamma -1)\displaystyle \frac{1}{\varepsilon }{\rho }_{0}\left({{\mathscr{L}}}_{T}{T}_{1}+{{\mathscr{L}}}_{\rho }{\rho }_{1}\right)\\ \ +\ 2i(\gamma -1)\eta \left\{{B}_{03}^{{\prime} }\left[\left(\displaystyle \frac{1}{\varepsilon }({a}_{2}^{{\prime} }-{k}_{2}{a}_{1})\right)^{\prime} +\displaystyle \frac{{k}_{2}{k}_{3}}{\varepsilon }{a}_{3}-\displaystyle \frac{{k}_{3}^{2}}{\varepsilon }{a}_{2}\right]\right.\\ \ -\ \left.\displaystyle \frac{(\varepsilon {B}_{02})^{\prime} }{{\varepsilon }^{2}}\left[\left(\varepsilon ({a}_{3}^{{\prime} }-{k}_{3}{a}_{1})\right)^{\prime} +\displaystyle \frac{{k}_{2}{k}_{3}}{\varepsilon }{a}_{2}-\displaystyle \frac{{k}_{2}^{2}}{\varepsilon }{a}_{3}\right].\right\}\\ \ +\ i(\gamma -1)\displaystyle \frac{d\eta }{{dT}}\displaystyle \frac{1}{\varepsilon }\left[{\left({B}_{02}^{{\prime} }\right)}^{2}+{\left({B}_{03}^{{\prime} }\right)}^{2}\right.\\ \quad \ +\ \left.2\displaystyle \frac{\varepsilon ^{\prime} }{\varepsilon }{B}_{02}{B}_{02}^{{\prime} }+{\left(\displaystyle \frac{\varepsilon ^{\prime} }{\varepsilon }{B}_{02}\right)}^{2}\right]{T}_{1},\end{array}\end{eqnarray} \tag{ 23 }$

$\begin{eqnarray}&&\begin{array}{l}\omega {a}_{1}=-{B}_{03}{v}_{2}+\displaystyle \frac{1}{\varepsilon }{B}_{02}{v}_{3}-\displaystyle \frac{1}{\varepsilon }{v}_{02}{a}_{2}^{{\prime} }-{v}_{03}{a}_{3}^{{\prime} }\\ \ +\ \left(\displaystyle \frac{{k}_{2}{v}_{02}}{\varepsilon }+{k}_{3}{v}_{03}\right){a}_{1}\\ +i\eta \left[-\left(\displaystyle \frac{{k}_{2}^{2}}{{\varepsilon }^{2}}+{k}_{3}^{2}\right){a}_{1}+\displaystyle \frac{{k}_{2}}{{\varepsilon }^{2}}{a}_{2}^{{\prime} }+{k}_{3}{a}_{3}^{{\prime} }\right],\end{array}\end{eqnarray} \tag{ 24 }$

$\begin{eqnarray}&&\begin{array}{l}\omega \displaystyle \frac{1}{\varepsilon }{a}_{2}=-\displaystyle \frac{1}{\varepsilon }{B}_{03}{v}_{1}+\displaystyle \frac{{k}_{3}{v}_{03}}{\varepsilon }{a}_{2}-\displaystyle \frac{{k}_{2}{v}_{03}}{\varepsilon }{a}_{3}\\ \ +\ i\eta \left[\left(\displaystyle \frac{1}{\varepsilon }\left({a}_{2}^{{\prime} }-{k}_{2}{a}_{1}\right)\right)^{\prime} +\displaystyle \frac{{k}_{2}{k}_{3}}{\varepsilon }{a}_{3}-\displaystyle \frac{{k}_{3}^{2}}{\varepsilon }{a}_{2}\right]\\ +i\displaystyle \frac{d\eta }{{dT}}\displaystyle \frac{1}{\varepsilon }{B}_{03}^{{\prime} }{T}_{1},\end{array}\end{eqnarray} \tag{ 25 }$

$\begin{eqnarray}\begin{array}{rcl}\omega {a}_{3} & = & \displaystyle \frac{1}{\varepsilon }{B}_{02}{v}_{1}-\displaystyle \frac{{k}_{3}{v}_{02}}{\varepsilon }{a}_{2}+\displaystyle \frac{{k}_{2}{v}_{02}}{\varepsilon }{a}_{3}+i\displaystyle \frac{1}{\varepsilon }\eta \left[\left(\varepsilon \left({a}_{3}^{{\prime} }-{k}_{3}{a}_{1}\right)\right)^{\prime} \right.\\ & & +\ \left.\displaystyle \frac{{k}_{2}{k}_{3}}{\varepsilon }{a}_{2}-\displaystyle \frac{{k}_{2}^{2}}{\varepsilon }{a}_{3}\right]-i\displaystyle \frac{d\eta }{{dT}}\displaystyle \frac{1}{{\varepsilon }^{2}}\left(\varepsilon {B}_{02}\right)^{\prime} {T}_{1}.\end{array}\end{eqnarray} \tag{ 26 }$

The perturbed thermal conductivity tensor κ_⊥,1 in Equation (23) written in terms of the perturbed variables is given by

$\begin{eqnarray}\begin{array}{rcl}\varepsilon {\kappa }_{\perp ,1} & = & \displaystyle \frac{\partial {\kappa }_{\perp ,0}}{\partial T}{T}_{1}+\displaystyle \frac{\partial {\kappa }_{\perp ,0}}{\partial \rho }{\rho }_{1}-2\varepsilon {B}_{02}\left({a}_{3}^{{\prime} }-{k}_{3}{a}_{1}\right)\displaystyle \frac{\partial {\kappa }_{\perp ,0}}{\partial \left({B}^{2}\right)}\\ & & +\ 2{B}_{03}\left({a}_{2}^{{\prime} }-{k}_{2}{a}_{1}\right)\displaystyle \frac{\partial {\kappa }_{\perp ,0}}{\partial \left({B}^{2}\right)}.\end{array}\end{eqnarray} \tag{ 27 }$

An interesting side note is that κ_∥,1 does not appear in the equations, which is due to the fact that this term is accompanied by a ${{\boldsymbol{B}}}_{0}\cdot {\rm{\nabla }}{T}_{0}$ contribution, and this is zero due to the equilibrium profile in Equation (9) (Van der Linden & Goossens 1991a, 1991b). We now have a system of eight ordinary differential equations in u₁ for the perturbed quantities ρ₁, v₁, v₂, v₃, T₁, a₁, a₂, and a₃.

2.3. Boundary Conditions

The above system of differential Equations (19)–(26) has to be complemented by a set of boundary conditions on both sides of the domain. For a Cartesian geometry, we look at a domain enclosed by two conducting walls. Clearly, the velocity component perpendicular to the walls has to be zero because there cannot be any propagation into a solid boundary. Mathematically, this translates into ${\boldsymbol{v}}\cdot {\hat{{\boldsymbol{e}}}}_{n}=0$ , where ${\hat{{\boldsymbol{e}}}}_{n}$ represents the normal vector to the wall. Following the same reasoning, we also require that ${\boldsymbol{B}}\cdot {\hat{{\boldsymbol{e}}}}_{n}=0$ , so in terms of a vector potential, this implies ${\hat{{\boldsymbol{e}}}}_{n}\times {\boldsymbol{A}}=0$ . Hence, applying this to the set of linearized equations means that for the Cartesian case v₁, a₂, and a₃, all have to be zero on the boundaries. Furthermore, because we are dealing with a perfectly conducting wall, one has to take care when thermal conduction is included. In that case, the rigid wall directly influences the temperature because it acts as an energy reservoir essentially eliminating the temperature perturbation. Hence, if and only if perpendicular thermal conduction is taken into account, we have to supplement the boundary conditions by the additional condition T₁ = 0 at the boundary. In theory there is a second possibility, which is treating the wall as a perfect insulator instead of a perfect conductor. In that case, there is no heat flux, which translates to the boundary condition ${T}_{1}^{{\prime} }$ = 0 instead of T₁ = 0. For now, we only consider the latter condition, that is, the one corresponding to a perfectly conducting wall.

In cylindrical geometry, we have the exact same boundary conditions as for the Cartesian case at the outer wall r = R, or at the outer edge of the accretion disk at R. The same is true at the inner disk edge, but for a flux tube extending to r = 0 we have to take the regularity conditions into account when treating the cylinder axis r = 0, which comes down to the fact that rv_r should go to zero when approaching r = 0. Looking back at the transformations (18) we applied, it follows that this condition is equivalent to v₁ = 0. Analogously, the same holds true for a₂ and a₃ such that these conditions are identical to the ones we applied for the Cartesian case, which is convenient implementation-wise. We again have to consider an additional condition if perpendicular thermal conduction is taken into account, because then $r{\widehat{T}}_{1}=0$ should also hold on the cylinder axis, which, similarly to v₁, translates into T₁ = 0 at r = 0.

In the case of confinement by a perfectly conducting wall, we thus have straightforward boundary conditions, that is, v₁ = 0, a₂ = 0, a₃ = 0, and T₁ = 0 for both the Cartesian and cylindrical geometries on both sides. This latter boundary condition should only be taken into account if and only if perpendicular thermal conduction is included.

2.4. Solving the Equations

The system of Equations (19)–(26) is solved through usage of a finite element discretization. Applying a weak Galerkin formalism turns this system of equations into a generalized matrix eigenvalue problem. A detailed explanation on how this is done can be found in the Appendix, where we describe the structure of the finite element approach and the matrix assembly process, along with a detailed treatment of how the boundary conditions are handled.

3. Results

As is common practice when developing a new numerical code, we tested Legolas against various results previously obtained in the literature. We divided this section into four subsections, each of which handles different physical effects. To begin with, we discuss results for adiabatic equilibria where only gravity is included in Section 3.1. In this case, we can compare numerical spectra obtained through Legolas with analytical solutions acquired by solving dispersion relations; here, we focus on stratified atmospheres containing p and g modes. We then move on to cylindrical geometries in Section 3.2 where we first look at adiabatic flux tubes, followed by the inclusion of flow effects by considering equilibria with KHIs and Suydam cluster modes. Next, the focus shifts to nonadiabatic effects in Section 3.3 by looking at a resistive MHD computation for a case without gravity, where a quasi-mode is known analytically. Resistive tearing modes are also discussed, combining the effects of flow and resistivity. Finally, Section 3.4 treats the inclusion of thermal conduction and optically thin radiative cooling effects, where we look at nonadiabatic discrete Alfvén waves and magnetothermal modes.

3.1. Cartesian Cases: Waves in Stratified Atmospheres

First of all, we discuss multiple theoretical results for adiabatic equilibria in a Cartesian geometry, where only gravity is included. We consider p and g modes in stratified layers and pay special attention to specific unstable branches.

3.1.1. Gravito-MHD Waves

The first test case covers gravito-MHD waves as discussed in Goedbloed et al. (2019, Figure 7.9), which handles an exponentially stratified atmosphere with constant sound and Alfvén speeds. This magnetized atmosphere contains the generalization of the p and g modes of an unmagnetized layer, and the constancy of the sound and Alfvén speed renders it analytically tractable, because the slow and Alfvén continua collapse to points. The geometry is Cartesian, with x ∈ [0, 1] and an equilibrium configuration given by

$\begin{eqnarray}\begin{array}{rcl}{\rho }_{0} & = & {\rho }_{c}{e}^{-\alpha x},\quad {p}_{0}={p}_{c}{e}^{-\alpha x},\quad {{\boldsymbol{B}}}_{0}={B}_{c}{e}^{-\tfrac{1}{2}\alpha x}{\hat{{\boldsymbol{e}}}}_{z},\\ \alpha & = & \displaystyle \frac{{\rho }_{c}g}{{p}_{c}+\tfrac{1}{2}{B}_{c}},\end{array}\end{eqnarray} \tag{ 28 }$

where p_c and B_c are taken to be 0.5 and 1, respectively, as to yield a plasma beta equal to unity. The parameter α is taken to be 20, which, together with g = 20, is used to constrain the value for the constant ρ_c. These four equations completely determine the equilibrium configuration, because the temperature is T₀ = p₀/ρ₀, following the ideal gas law. The spectrum discussed in Goedbloed et al. (2019) is actually the solution to the analytic dispersion relation for gravito-MHD waves, which shows the squared eigenvalue as a function of wavenumber for a fixed angle θ = π/4 between the wavevector k ₀ and the magnetic field ${{\boldsymbol{B}}}_{0}$ . However, the spectrum as calculated by Legolas corresponds to one single equilibrium configuration, meaning one value for k_y and k_z. In order to reproduce Figure 7.9 from Goedbloed et al. (2019) and compare the results, we performed 100 different runs where the equilibrium parameters in Equation (28) remained unchanged, but k_y and k_z took on 100 different values between 0 and $\sqrt{250}$ as to yield a wavenumber range for ${k}_{0}^{2}$ between 0 and 500. Because the magnetic field is purely aligned with the z-axis, we can write ${k}_{\parallel }={k}_{z}={k}_{0}\cos (\theta )$ and ${k}_{\perp }={k}_{y}={k}_{0}\sin (\theta )$ . All runs were performed using 351 grid points, yielding a matrix size of 5616 × 5616.

Our results are shown in Figure 2, where every vertical collection of points at the same ${k}_{0}^{2}$ value represents one single Legolas run. Because we are in an MHD regime with β = 1, the three MHD subspectra can be clearly distinguished, showing the fast p modes (top-left branches), Alfvén g modes (middle branches), and slow g modes (bottom branches). The inset shows a zoom-in near the marginal frequency of the spectrum, showing unstable (ω² < 0) slow MHD modes. These long-wavelength unstable modes are related to the Parker instabilities, due to magnetic buoyancy, as we will show in Section 3.1.2. Note that because this case is adiabatic and fully self-adjoint, every individual MHD spectrum is left–right and up–down symmetric in the complex eigenfrequency plane, but this aspect is hidden from the ${\omega }^{2}-{k}_{0}^{2}$ view shown here.

3.1.2. Quasi-Parker Instabilities

Next we discuss a modified case of the gravito-MHD waves, namely a spectrum showing quasi-Parker instabilities as done in Goedbloed et al. (2019, Figure 12.2). The difference with the previous case is that a fully analytic description is no longer possible, because the introduction of magnetic shear leads to continuous ranges in the MHD spectrum. Instead of showing the spectrum for one single value for θ, we now vary the direction of the wavevector k ₀ between 0 and π. The equilibrium configuration is similar to the one in Section 3.1.1, given in Cartesian geometry by

$\begin{eqnarray}\begin{array}{rcl}{\rho }_{0} & = & {\rho }_{c}{e}^{-\alpha x},\quad {p}_{0}={p}_{c}{e}^{-\alpha x},\quad \alpha =\displaystyle \frac{{\rho }_{c}g}{{p}_{c}+\tfrac{1}{2}{B}_{c}},\\ {B}_{02} & = & {B}_{c}{e}^{-\tfrac{1}{2}\alpha x}\sin (\lambda x),\quad {B}_{03}={B}_{c}{e}^{-\tfrac{1}{2}\alpha x}\cos (\lambda x),\end{array}\end{eqnarray} \tag{ 29 }$

where magnetic shear was introduced through the parameter λ.

The quantities α and B_c are assigned the same values as in Equation (28), except that g = 0.5 and p_c = 0.25, which yield a plasma beta β = 0.5. The wavevectors are given by ${k}_{y}=\pi \sin (\theta )$ and ${k}_{z}=\pi \cos (\theta )$ , such that ${k}_{0}^{2}\approx 10$ . The angle θ was varied between 0 and π for a total of 100 runs at 351 grid points each, shown in Figure 3.

**Figure 3.** Spectrum showing Parker and quasi-Parker modes without (left) and with (right) magnetic shear. The slow and Alfvén continua are shown in red and cyan, respectively, where the insets zoom into the region of quasi-interchange modes. The bottom row of panels show the eigenfrequency view for the single case θ = 0.3π. The continua are again annotated in the figures, visualizing the collapsed single point values (left) as well as the genuine continuum ranges (right). The gray dashed line in the top two panels denotes ω = 0.
Download figure:
Standard image High-resolution image

The left panels handle the case without magnetic shear, that is, λ = 0, which basically reduces to the one from the previous subsection. In this case, the slow and Alfvén continua collapse into single point values, denoted in red and cyan, respectively. The right panels show the same configuration where λ = 0.3 was taken, introducing magnetic shear, which introduces genuine continua seen as bands. These continua affect the overall stability and organize the entire MHD spectrum: all discrete modes are fully aware of the essential spectrum formed by these (slow and Alfvén) continua and the (fast) accumulation points at infinite frequency. All features of the original figure in Goedbloed et al. (2019) are reproduced. The inset zooms into the region where both continua overlap, showing quasi-interchange and interchange instabilities. Once more, each run, shown here collectively in Figure 3, actually has a spectrum that is left–right and up–down symmetric in the eigenfrequency plane. This is depicted on the bottom two panels, which show the eigenfrequency view for one single case (θ = 0.3π). The continuum ranges separate nicely: the collapsed single point values are denoted by cyan (Alfvén) and red (slow) points on the left panel, and the genuine continua are shown with cyan and red bands on the right panel. The instabilities themselves are situated on the (positive) imaginary axis, due to the self-adjointness of the eigenvalue problem mentioned earlier.

As explained in Goedbloed et al. (2019), we see from this eigenmode computation that the Parker instability, which is there for k ₀ parallel to ${{\boldsymbol{B}}}_{0}$ , becomes a quasi-Parker instability away from perfect alignment and connects smoothly to well-known quasi-interchange instabilities that occur here (marginally) away from perpendicular orientation. Quantifying how the equilibrium parameters influence the growth rates of these unstable branches can only be done numerically, e.g., with Legolas.

3.2. Adiabatic, Cylindrical Cases

Next we move on to cylindrical configurations, which provide tests for the scale factor $\varepsilon$ in the equations. Analytical results from the literature are again well reproduced. Furthermore, we look at different spectra previously obtained by the LEDA code, discussed in various papers, and compare those with the new spectra from Legolas.

3.2.1. Magnetic Flux Tubes

The first case that we describe in this subsection is a magnetic flux tube embedded in a uniform magnetic environment, discussed in Roberts (2019). The equilibrium configuration is simple, in the sense that we have a uniform magnetic field aligned with the z-axis both inside and outside of the flux tube, with a similar structure for the other equilibrium parameters:

$\begin{eqnarray}{B}_{0}(r),{\rho }_{0}(r),{p}_{0}(r),{T}_{0}(r)=\left\{\begin{array}{ll}{B}_{0},{\rho }_{0},{p}_{0},{T}_{0}, & r\lt a\\ {B}_{e},{\rho }_{e},{p}_{e},{T}_{e}, & r\gt a\end{array}\right.\end{eqnarray} \tag{ 30 }$

where the subscripts 0 and e refer to values inside the tube and for the environment, respectively. The outer radius of the tube is denoted by a and hence represents a discontinuous interface between the tube itself and the environment. Because total pressure balance should be preserved across the boundary, which is something that follows from Equation (10), this yields a relation between pressures and magnetic field components inside and outside of the tube, which in turn implies a connection between the plasma densities, sound speeds, and Alfvén speeds across the boundary:

$\begin{eqnarray}&&{p}_{0}+\displaystyle \frac{1}{2}{B}_{0}^{2}={p}_{e}+\displaystyle \frac{1}{2}{B}_{e}^{2},\displaystyle \frac{{\rho }_{e}}{{\rho }_{0}}=\displaystyle \frac{{c}_{s}^{2}+\tfrac{1}{2}\gamma {c}_{A}^{2}}{{c}_{{se}}^{2}+\tfrac{1}{2}\gamma {c}_{{Ae}}^{2}},\end{eqnarray} \tag{ 31 }$

where ${c}_{s}^{2}=\gamma {p}_{0}/{\rho }_{0}$ and ${c}_{A}^{2}={B}_{0}^{2}/{\rho }_{0}$ denote the sound speed and Alfvén speed, respectively, in which the values outside of the flux tube are used if there is a subscript e present.

It should be noted that this extremely simple equilibrium configuration is the standard case used in many solar coronal loop seismology efforts. Because it simply has two uniform media (one inside the tube and one in its exterior), it has no continuous spectra (they reduce to point values), but the interface makes it possible to again have surface modes that would be affected by true radial variation. Also note that these flux tubes have only stable waves, but we can distinguish between body and surface waves, depending on the variation of the eigenfunctions within the flux tube. In the exterior of the flux tube all eigenfunctions are exponentially varying.

We should also clarify here that the original dispersion relation as given in Roberts (2019) assumes a flux tube embedded in an environment extending toward infinity, while Legolas on the other hand assumes a fixed wall boundary at the outer edge of the domain. Hence, we assume here that the domain is situated in r ∈ [0, 10] with the inner flux tube wall at r = 1 in order to minimize the outer wall influence. However, this introduces an additional computational challenge, in the sense that we are (mainly) interested in the behavior of the inner modes, because we know that the outer modes all have exponentially varying eigenfunctions which decay to infinity (or toward our far-away outer wall). Hence, in order to resolve those inner waves huge resolutions are needed due to the 1:10 ratio. In order to circumvent this issue we used a simple prescription for mesh refinement, that is, a 60−30−10 division of the initial nodes. This means that 60% of the grid points are used for the inner tube region r ∈ [0, 0.95], 30% of the grid points are located near the transition region r ∈ [0.95, 1.05], and the remaining 10% are used for the environment r ∈ [1.05, 10].

(a) Photospheric flux tube. First, we look at a flux tube under photospheric conditions, that is, an equilibrium for which c_Ae < c_s < c_se < c_A. More specifically, we take c_Ae = c_s/2, c_se = 3c_s/2, and c_A = 2c_s following Roberts (2019, Figure 6.5). The relations between the inner and outer regions of the flux tube follow straightforwardly from Equation (31), and hence we only have two degrees of freedom, namely ρ₀ and p₀, which are both taken to be unity, with γ = 5/3. This results in ρ₀ ≈ 0.567 ρ_e, c_t ≈ 0.89 c_s, c_k ≈ 0.63 c_A ≈ 1.27 c_s. Here we also introduced the tube and kink speeds, given by

$\begin{eqnarray}&&{c}_{t}^{2}=\displaystyle \frac{{c}_{s}^{2}{c}_{A}^{2}}{{c}_{s}^{2}+{c}_{A}^{2}},{c}_{k}^{2}=\displaystyle \frac{{\rho }_{0}{c}_{A}^{2}+{\rho }_{e}{c}_{{Ae}}^{2}}{{\rho }_{0}+{\rho }_{e}},\end{eqnarray} \tag{ 32 }$

using the same notation as in Equation (31). As described before, we take r ∈ [0, 10] and place the flux tube boundary at a = 1. Next we perform 40 runs at 300 grid points each for four azimuthal wavenumbers m = 0 to m = 3. For the wavenumber k₃ = k_z, we take 40 values in such a way that the dimensionless wavenumber k_z a has values in [0, 6.2]. The spectrum showing the dispersion relation, where the dimensionless phase speed ω/k_z c_s is plotted as a function of k_z a, is depicted in Figure 4. All speeds indicated in the figure are normalized to the internal sound speed c_s. Panel (a) clearly shows the fast surface waves, including the sausage (m = 0), kink (m = 1), and first two fluting modes (m = 2 and m = 3). The eigenfunctions in the top-right panel (c) correspond to the modes annotated with a transparent circle in panel (a), with colors indicating the mode number m. The left and right sides of panels (c)–(f) show the eigenfunctions for the inner and outer regions of the flux tube, respectively. All eigenfunctions are normalized to their maximum value, and all eigenvalues having a normalized phase speed larger than 1.5 are not shown. It should be noted that the first three runs, that is, the first three dots according to the k_z a axis, were done using 1001 grid points. The reason for this is that when we divide the eigenvalues ω by k_z, small errors are increased selectively for small k_z, explaining why those modes seem slightly scattered and hence why we have to employ such high resolutions in order to minimize said error.

Panel (b) of Figure 4 zooms in between the tube speed (c_t) and internal sound speed (c_s), showing a clear representation of the various body waves that accumulate to the tube speed at long wavelengths. Panel (e) shows the first three modes of the m = 0 and m = 1 sequences in solid, dashed, and dotted lines, respectively; that is, the first mode in the sequence (solid line) corresponds to the mode annotated in panel (b). The next two modes in that sequence are the next two blue dots moving vertically downwards (same k_z a value), these are not annotated to avoid cluttering the figure. Analogously, panel (f) shows the first three modes of the m = 2 and m = 3 sequences, with everything color-coded according to the legend. As indicated before, all eigenfunctions in the outer region are exponentially varying. Panels (a) and (b) in Figure 4 reproduce the analytical results in Roberts (2019, Figure 6.5).

One mode has not yet been discussed, and that is the horizontal line of modes between the internal sound speed and kink speed. The eigenfunctions are shown in panel (d) and correspond to the annotated squares in the top-left panel (a). These modes are not present in the original work, and it is not a priori clear how to interpret these. However, what we do know is that their phase speed is approximately 0.5(c_k + c_s) and that these modes are degenerate, meaning that their position does not change when the mode number m or wavenumber k_z changes. The position of the outer wall also does not seem to have any influence on their value. This, together with the fact that the eigenfunctions seem to indicate that these are surface waves, strengthens the belief that these solutions could be actual waves and not some numerical remnant.

(b) Coronal flux tube. The second application of the magnetic flux tube is one under coronal conditions, that is, an equilibrium for which c_se < c_s < c_A < c_Ae. More specifically, we take c_Ae = 5c_s, c_se = c_s/2, and c_A = 2c_s. Analogous to the previous case, the relations between the equilibrium values inside and outside of the flux tube follow from Equation (31), where we again take p₀ and ρ₀ to be equal to one. This results in ρ₀ ≈ 4.86 ρ_e, c_t ≈ 0.89 c_s, and c_k ≈ 1.38 c_A ≈ 0.55 c_Ae, and we use the same values as for the photospheric case for k_z and the flux tube and outer wall boundaries. Similar to case (a), we perform 40 runs at 300 grid points each for four values of k₂ = m (where again the first three runs have 1000 grid points) and plot the spectrum showing the dispersion relation in Figure 5. Again, all speeds indicated in the figure are normalized to the sound speed c_s. The top-left panel (a) focuses on the fast body waves, the bottom-left panel (b) on the slow body waves. Panel (c) shows the eigenfunctions corresponding to the four body waves indicated with circles in panel (a). Panel (d) depicts the eigenfunctions of the modes annotated with squares between the internal Alfvén (c_A) and kink (c_k) speeds, representing the kink m = 1 and first two fluting (m = 2 and m = 3) modes.

**Figure 5.** Panel (a): spectrum showing the dispersion relation for a flux tube under coronal conditions. The dimensionless phase speed ω/k_z c_s is displayed as a function of the dimensionless wavenumber k_z a for four values of the azimuthal wavenumber m = 0 (blue dots), m = 1 (red crosses), m = 2 (green squares), and m = 3 (yellow triangles). Both Alfvén speeds, together with the kink, tube, and sound speeds are denoted by gray horizontal lines and are all normalized to c_s. Panel (b) zooms in on the tube-speed-related body mode sequences. Panel (b): zoom-in of panel (a) between the tube speed c_t and internal sound speed c_s. Panels (c) and (d): eigenfunctions of the modes annotated with circles (c) and squares (d) on the top-left panel (a). Panel (e): first three eigenfunctions of the m = 0 and m = 1 body wave sequence in the bottom-left panel (b). Panel (f): first three eigenfunctions of the m = 2 and m = 3 body wave sequence. For both panels (e) and f the n = 1, n = 2 and n = 3 modes are shown in solid, dashed and dotted lines, respectively, only the first mode is annotated on the bottom-left panel (b). Every eigenfunction in panels (c)–(f) is normalized to its maximum absolute value in that particular grid interval.
Download figure:
Standard image High-resolution image

Similar to Figure 4, panels (e) and (f) represent the first three body modes in the m = 0 and m = 1 sequences (panel (e)), of which the first mode is indicated in panel (b). Panel (f) shows the first three modes of the m = 2 and m = 3 sequences, with the first, second, and third modes indicated with a solid, dashed and dotted line, respectively. The colors of panels (c)–(f) are consistent with the legend. Figure 5 reproduces the analytical results in Roberts (2019, Figure 6.7), which are based on the analytic dispersion relation containing Bessel functions.

3.2.2. Tokamak Constant Current

Next we discuss an example initially given in Kerner et al. (1985), which shows the ideal MHD spectrum in the presence of an unstable m = −2 interchange mode in a cylindrical geometry. We start from a so-called tokamak current profile, in which an axial current density of the form ${{\boldsymbol{j}}}_{0}=(0,0,{j}_{0}{\left(1-{r}^{2}\right)}^{\nu })$ is assumed, with j₀ a given constant. This yields a twisted magnetic profile in which the longitudinal component B_0z is uniform and equals one, while the poloidal component B_0θ is given by

$\begin{eqnarray}&&{B}_{0\theta }=\displaystyle \frac{{j}_{0}}{2r(\nu +1)}\left[1-{\left(1-{r}^{2}\right)}^{\nu +1}\right],\end{eqnarray} \tag{ 33 }$

for a given value of ν. For the equilibrium considered here, we take ν = 0, making the current profile constant over the flux tube. This means that B_0θ has a linear profile in r such that the magnetic field lines have a constant pitch and the current is distributed equally in the plasma. An expression for the pressure (and hence temperature) can be found by integrating Equation (10) and assuming that, for example, the pressure vanishes at the outer boundary, resulting in a parabolic pressure profile. Hence, for a cylindrical geometry in which r ∈ [0, 1], this yields the following equilibrium configuration:

$\begin{eqnarray}\begin{array}{rcl}{\rho }_{0} & = & 1,\quad {p}_{0}=\displaystyle \frac{1}{4}{j}_{0}^{2}\left(1-{r}^{2}\right),\\ {B}_{02} & = & \displaystyle \frac{1}{2}{j}_{0}r,\quad {B}_{03}=1,\end{array}\end{eqnarray} \tag{ 34 }$

where we assumed a uniform density. We introduce an additional parameter q, called the safety factor, given by

$\begin{eqnarray}&&q(r)=\displaystyle \frac{{{rk}}_{z}{B}_{0z}}{{B}_{0\theta }}=\displaystyle \frac{2{k}_{z}}{{j}_{0}}.\end{eqnarray} \tag{ 35 }$

We performed 39 runs, varying the q factor between 1.9 and 2.1 in order to probe the regime containing the unstable m = −2 interchange mode, which is associated with the vanishing of the factor m + kq, that is, the ${\boldsymbol{k}}\cdot {\boldsymbol{B}}$ product. This implicitly constrains the value for j₀, and we assigned k₂ = m = −2 and k₃ = k = 0.2 for all runs. It should be noted that this particular equilibrium configuration requires a high resolution near q = 2 to correctly resolve the unstable modes. We thus used 501 grid points for all runs. The complete spectrum is shown in Figure 6, where the squared eigenvalues are plotted as a function of the safety factor. The three main branches, that is, fast, Alfvén, and slow, are denoted on the right side of the figure, as well as the region where the modes become unstable (ω² < 0). We see that in the region near the m = −2 interchange instabilities, the slow and Alfvén modes collapse to zero, which is due to the vanishing of the combination F = mB_θ/r + kB_z in the ideal MHD equations (Goedbloed et al. 2019). The slow and Alfvén continua are annotated in the figure in red and cyan, respectively. The Alfvén continuum is collapsed to a single point for this equilibrium configuration, while the slow continuum covers a range in frequency. Note in particular how the full spectrum ranges over many orders of magnitude in the ω² view shown here: an intrinsic property and challenge posed by MHD spectral theory.

**Figure 6.** Complete MHD spectrum for a tokamak current profile in the presence of the m = −2 interchange modes. The squared eigenvalues are plotted against the safety factor q, instabilities correspond to ω² < 0. The various branches are indicated on the right side of the figure.
Download figure:
Standard image High-resolution image

3.2.3. KH and CD Instabilities

As a first test for the inclusion of flow into the equations, we look at the interaction between the KHIs and current-driven (CD) instabilities in a magnetized astrophysical jet, following Baty & Keppens (2002). This model uses a cylindrical jet with a supersonic background flow aligned with the axis and sheared in the radial direction. The equilibrium configuration is taken such that KH surface modes can develop and is generally given by

$\begin{eqnarray}\begin{array}{rcl}{\rho }_{0} & = & 1,{v}_{02}=0,\quad {v}_{03}=\displaystyle \frac{1}{2}V\tanh \left(\displaystyle \frac{{r}_{j}-r}{a}\right),\\ {B}_{02} & = & {B}_{\theta 0}\displaystyle \frac{{{rr}}_{c}}{{r}_{c}^{2}+{r}^{2}},\quad {B}_{03}={B}_{z0},\\ {T}_{0} & = & \displaystyle \frac{{p}_{0}}{{\rho }_{0}}-\displaystyle \frac{{B}_{\theta 0}^{2}}{2{\rho }_{0}}\left(1-\displaystyle \frac{{r}_{c}^{4}}{{\left({r}_{c}^{2}+{r}^{2}\right)}^{2}}\right).\end{array}\end{eqnarray} \tag{ 36 }$

Here, r_j denotes the jet radius and r_c quantifies the radial variation, and they are taken to be r_j = 1 and r_c = 0.5, with r ∈ [0, 2]. The parameter V represents the amplitude of the velocity, given by V = 1.63, while the chosen velocity profile ensures that the shear layer is situated at the jet radius with a radial width given by a = 0.1r_j. Both the density ρ₀ and the pressure on axis p₀ are chosen to be equal to unity. The parameters B_θ0 and B_z0 control the amplitude and twist of the magnetic field, respectively, as explained in Baty & Keppens (2002). The original work discusses three different magnetic field configurations, we choose the profile with

$\begin{eqnarray}&&{B}_{\theta 0}=0.4\left(\displaystyle \frac{{r}_{c}^{2}+{r}_{j}^{2}}{{r}_{j}{r}_{c}}\right),{B}_{z0}=0.25,\end{eqnarray} \tag{ 37 }$

such that B₀₂ = 0.4 at the jet radius. Furthermore, the azimuthal and longitudinal wavenumbers are taken to be k₂ = m = −1 and k₃ = k = π.

The entire spectrum is calculated at high resolution using 501 grid points and shown in Figure 7. Panel (a) depicts the full slow and Alfvén spectra with the KH and first three CD unstable modes (ω_I > 0) denoted in the figure itself. The real and imaginary parts of the ρ, rv_r, and v_z eigenfunctions for each of these modes are shown on the subsequent panels. Those for the KH mode (panels (c)–(d)) are localized around the jet radius r = 1, which is also the point where the sonic and Alfvénic Mach numbers drop to zero (panel (b)). Panels (e) through (j) depict the eigenfunctions of the first three CD modes, with the first, second, and third shown in the first, second, and third rows of the bottom panels, respectively. The left column shows the real part, the right column the imaginary part of the eigenfunction. The CD modes have an increasing number of nodes on r ∈ [0, 1], which is most clearly visible by looking at the rv_r eigenfunction (green): no nodes for the first CD mode (panels (e)–(f)), one node for the second CD mode (panels (g)–(h)), and two nodes for the third CD mode (panels (i)–(j)).

Note that here, we are still adiabatic such that the up–down symmetry (relating to time reversal) is still present in the eigenfrequency plane, but the introduction of equilibrium flow caused left–right symmetry breaking between the forwards and backwards propagating modes. As pointed out in Goedbloed (2018a), the study of MHD spectra of stationary (with flow) equilibria is still governed by two self-adjoint operators, but as seen in Figure 7, modes can enter the complex plane at various locations (identified by the spectral web; Goedbloed 2018b). The correspondence with the original figure in Baty & Keppens (2002) is one to one for the KH and CD modes, but here we have a lot more detail near the axes due to the higher resolution. We will discuss this resolution aspect in Section 4.

3.2.4. Suydam Cluster Modes

Next we look at Suydam cluster modes in a cylindrical geometry, which arise from the presence of a Suydam surface in the equilibrium configuration, that is, a location where ${\boldsymbol{k}}\cdot {{\boldsymbol{B}}}_{0}=0$ . Shear flow effects are included, and the equilibrium is given by

$\begin{eqnarray}\begin{array}{rcl}{\rho }_{0} & = & 1,\quad {v}_{02}=0,\quad {v}_{03}={v}_{z0}\left(1-{r}^{2}\right),\\ {B}_{02} & = & {J}_{1}(\alpha r),\quad {B}_{03}=\sqrt{1-{P}_{1}}{J}_{0}(\alpha r),\\ {p}_{0} & = & {P}_{0}+\displaystyle \frac{1}{2}{P}_{1}{J}_{0}^{2}(\alpha r),\end{array}\end{eqnarray} \tag{ 38 }$

where α = 2, P₀ = 0.05, P₁ = 0.1, and v_z0 = 0.14, and the functions J₀ and J₁ denote the Bessel functions of the first kind. The wavenumbers were chosen to be k₂ = m = 1 and k₃ = k = −1.2, ensuring a Suydam surface at r ≈ 0.483.

The spectrum is calculated using 501 grid points for r ∈ [0, 1] and is shown in Figure 8. The resulting locations of the various off-axis outer modes are in agreement with the results given in Nijboer et al. (1997). However, because this spectrum is calculated using a five times higher resolution, we have much more intricate detail near the Suydam surface, revealing even more off-axis modes (inset in panel (a)). The top two panels on the right side of Figure 8 show the sonic and Alfvénic Mach numbers (panel (b)), together with the ${\boldsymbol{k}}\cdot {{\boldsymbol{B}}}_{0}={B}_{02}{k}_{2}/r\,+{B}_{03}{k}_{3}$ and ${\boldsymbol{k}}\cdot {{\boldsymbol{v}}}_{0}={B}_{02}{v}_{02}/r+{B}_{03}{v}_{03}$ profiles as a function of radius (panel (c)), respectively. The location of the Suydam surface at r ≈ 0.483 is denoted with a red cross. The bottom row of panels (d)–(g) shows the real part of the ρ and rv_r eigenfunctions, for the four modes annotated in panel (a) with the location of the Suydam surface annotated with a vertical red line. ω₁ and ω₃ correspond to modes on the left side of the Suydam surface, the other two correspond to modes on the right side. All eigenfunctions shown here show the specific variation associated with their location relative to the Suydam surface: ω₁ and ω₃ have their localized behavior on the left side of the red dashed line, while it is vice versa for the other two modes. For visual purposes the horizontal axis in panels (e) and (f) is different from the one in panels (c) and (d), because the former correspond to the next modes in the Suydam sequence, showing stronger radial variation. Note that the original Suydam criterion (Goedbloed et al. 2019) is related to static equilibria; the generalization of these Suydam cluster criteria is presented in Wang et al. (2004).

**Figure 8.** Panel (a): Suydam cluster spectrum for the equilibrium given in (38). Inset: zoom-in near the Suydam surface, revealing more detail. Panel (b) depicts the sonic and Alfvénic Mach numbers as a function of radius, panel (c) shows the ${\boldsymbol{k}}\cdot {{\boldsymbol{B}}}_{0}$ and ${\boldsymbol{k}}\cdot {{\boldsymbol{v}}}_{0}$ curves where the red cross denotes the Suydam surface. The bottom row of panels shows the real part of the ρ and rv_r eigenfunctions, for the four modes in the Suydam sequence annotated in panel (a). The dotted red line denotes the location of the Suydam surface.
Download figure:
Standard image High-resolution image

3.3. Resistive, Cartesian Cases

All cases discussed up to now handled an adiabatic equilibrium configuration with or without the inclusion of flow. The up–down symmetry of all the MHD spectra shown so far is perfectly maintained, related to the fact that these cases are in essence time-reversible. Every instability (or overstability in the case with flow) has a damped counterpart. We now move on to include additional effects. Hence, we now compute spectra for time-irreversible cases, where either resistivity or other nonadiabatic effects enter, which will break the up–down symmetry. Here we first focus on the inclusion of resistivity.

3.3.1. Resistive Homogeneous Plasma

First we look at the simplest configuration, that is, a homogeneous plasma in a Cartesian geometry with resistivity included. The uniform equilibrium is given by

$\begin{eqnarray}&&{\rho }_{0}=1,\quad {B}_{02}=0,\quad {B}_{03}=1,\quad {T}_{0}=\displaystyle \frac{1}{2}\beta {B}_{0}^{2},\end{eqnarray} \tag{ 39 }$

where we take a plasma beta of β = 0.25, k₂ = k_y = 0, k₃ = k_z = 1, and x ∈ [0, 1]. The value for the resistivity is assumed to be constant and given by η = 10⁻³, as described in Goedbloed et al. (2019). This spectrum is calculated using 1001 grid points; the result is shown in Figure 9. In ideal MHD the fast modes form a Sturmian sequence (that is, the oscillation of the eigenfrequencies increases when the number of modes increases, which in turn implies a larger real part of ω) of stable fast magnetoacoustic waves with frequencies accumulating to infinite frequency (related to the p modes in our stratified example). The slow modes have an anti-Sturmian sequence toward their accumulation point (in essence the collapsed slow continuum) and the Alfvén modes are degenerate (Goedbloed et al. 2019). When resistivity is included, the fast modes become damped and the Alfvén and slow modes trace out semicircles on the bottom-half of the complex plane. The semicircles and initial fast-mode sequence shown in Figure 9 are in perfect agreement with the spectrum depicted in Goedbloed et al. (2019, Figure 14.6). These semicircles can be quantified analytically; their radius does not depend on the resistivity. The magnetic Reynolds number, calculated as R_m = (x₁ − x₀)c_A/η, is equal to 1000.

**Figure 9.** Panel (a): spectrum for a homogeneous medium with constant resistivity. Panel (b) zooms in near the origin, showing the start of the fast-mode sequence. Panel (c) zooms in further, revealing the semicircles traced out by the Alfvén and slow modes (outer and inner semicircles, respectively).
Download figure:
Standard image High-resolution image

Due to the rather extreme resolution employed here, we trace much further into the fast-mode sequence, where we see something interesting: the fast modes appear on curves that loop around, breaking the purely Sturmian behavior for a moment, after which they continue again toward infinity. This implies that initially, the fast modes become more damped at higher mode frequencies up to a certain turning point at which they achieve maximal damping (the bottom of the loop). After passing this turning point, the oscillation frequency of the modes increases again and the damping frequency seems to converge toward one single value.

Of course, the strong damping for the fast modes as we go further into the fast-mode sequence must have consequences for the original uniform (ideal) equilibrium state. In what follows, we will adopt the common practice of computing resistive MHD spectra about an ideal MHD state, which itself will evolve when resistivity is acting.

3.3.2. Quasi-modes in Resistive MHD

We now turn to a nonadiabatic case, where resistivity is important. We will compute the resistive MHD spectrum for a case where the equilibrium varies across an interface, which gives rise to so-called quasi-modes. These are essentially surface waves undergoing damping, due to the fact that the global quasi-mode overlaps in frequency with the continuum range, causing resonant absorption. Quasi-modes are quite important in solar physics, as they can be indirectly related to the coronal heating problem as discussed in, for example, Poedts et al. (1989) and Poedts & Kerner (1991). A detailed analytical treatment of quasi-modes including theoretical growth rates is given in Priest (2014), where they start from an inhomogeneous layer of width l connecting two regions of uniform plasma. This can be reproduced by introducing a linear density profile between two homogeneous regions in a Cartesian geometry; however, this would mean that the density derivative shows rather strong discontinuities near the edges of the transition layer. We therefore opt for a smooth profile by introducing a sine dependence such that

$\begin{eqnarray}{\rho }_{0}(x)=\left\{\begin{array}{ll}{\rho }_{1} & \mathrm{if}\,{x}_{0}\leqslant x\lt s-\displaystyle \frac{1}{2}l,\\ \displaystyle \frac{1}{2}{\rho }_{1}\left[1+\displaystyle \frac{{\rho }_{2}}{{\rho }_{1}}-\left(1-\displaystyle \frac{{\rho }_{2}}{{\rho }_{1}}\right)\sin \left(\displaystyle \frac{\pi (x-s)}{l}\right)\right] & \mathrm{if}\,s-\displaystyle \frac{1}{2}l\leqslant x\leqslant s+\displaystyle \frac{1}{2}l,\\ {\rho }_{2} & \mathrm{if}\,s+\displaystyle \frac{1}{2}l\lt x\lt {x}_{1},\end{array}\right.\end{eqnarray} \tag{ 40 }$

where s denotes the midpoint between x₀ and x₁, representing the left and right edges of the Cartesian grid, respectively, and equal to 0 and 1 such that s = 0.5. The width of the transition region l = 0.1 is taken to be small as to reduce the influence of the walls. If $l\to 0$ , the inhomogeneous region disappears, and we simply have a discontinuous jump in the plasma. Furthermore, we take ρ₁ = 0.9, ρ₂ = 0.1, B₀₂ = 0, B₀₃ = 1, and T₀ = 0. The magnetic field is unidirectional along the z-axis, and setting the pressure (and thus temperature) to zero provides an additional test on the handling of zero rows in the matrix. This zero-temperature case is frequently encountered in fully nonlinear MHD simulations which artificially adopt a zero plasma beta. It has as an important consequence the slow continuum collapsing to marginal frequency, eliminating many interesting modes from the spectrum. Additionally, we adopt wavenumbers k₃ = 0.05 and k₂ = 1.0, such that k₃ ≪ k₂. As discussed in Priest (2014), the quasi-modes are damped, such that they move away from the real eigenfrequency axis with complex eigenvalues ω = ω_R + iω_I given by

$\begin{eqnarray}&&{\omega }_{R}=\pm \sqrt{\displaystyle \frac{{\rho }_{1}{\omega }_{A1}^{2}+{\rho }_{2}{\omega }_{A2}^{2}}{{\rho }_{1}+{\rho }_{2}}},{\omega }_{I}=-\displaystyle \frac{\pi {k}_{z}l\left({\omega }_{A2}^{2}-{\omega }_{A1}^{2}\right)}{8{\omega }_{R}},\end{eqnarray} \tag{ 41 }$

with ${\omega }_{A}^{2}={k}_{z}^{2}{v}_{A}^{2}$ . This analytic result originates from a complex analysis of the linearized initial value problem. However, there is one caveat: it is impossible to obtain complex eigenvalues away from the real or imaginary axes in an ideal plasma, due to the matrix operator being Hermitian in ideal MHD. Because the (ideal) quasi-mode is damped, we therefore include a (small) value for the resistivity, η = 10⁻⁴, which is sufficient to make the matrix operator non-Hermitian and allows for complex eigenvalues away from the horizontal axis.

The magnetic Reynolds number in this case varies between R_m ≈ 10⁴ and R_m ≈ 3 × 10⁴. Figure 10 shows the MHD spectrum in panel (a), for 501 grid points, with the theoretical prediction for the global quasi-mode location annotated with a red dot. The actual quasi-mode is also denoted in the figure and is slightly more damped than its theoretical counterpart. This is to be expected, because the inclusion of resistivity imposes additional damping on the eigenvalues, shifting them downwards. It should be noted that increasing the width l of the transition layer moves the quasi-mode farther down and to the right in the figure, such that it eventually merges with the damped modes found on the branch immediately to its right. This phenomenon, where the quasi-modes merge with the resistive branches, is discussed in more detail in Van Doorsselaere & Poedts (2007). The inset in panel (a) depicts the ρ (blue) and v_x (orange) eigenfunctions associated with the quasi-mode, showing localized variation near the transition region as expected. Note that the spectrum as shown in Figure 10, left panel (a), is still left–right symmetric, but that resistivity has broken the up–down symmetry, with all modes here found in the stable (damped) half plane.

3.3.3. Resistive Tearing Modes

Next, we move on to an inhomogeneous medium with the inclusion of resistivity and an optional linear flow profile, discussed in Goedbloed et al. (2019). The geometry is Cartesian, with x ∈ [−0.5, 0.5] and an equilibrium configuration given by

$\begin{eqnarray}\begin{array}{rcl}{\rho }_{0} & = & 1,\quad {B}_{02}=\sin (\alpha x),\quad {B}_{03}=\cos (\alpha x),\\ {T}_{0} & = & \displaystyle \frac{1}{2}\beta {B}_{0}^{2},\quad {v}_{02}={Vx},\quad {v}_{03}=0,\end{array}\end{eqnarray} \tag{ 42 }$

where k₃ = k_z = 0, β = 0.15, and α = 4.73884, with a constant resistivity value of η = 10⁻⁴. The magnetic configuration chosen here is a linear force-free field with shear, with α the proportionality constant between the current and magnetic field, that is, ${\boldsymbol{j}}=\alpha {{\boldsymbol{B}}}_{0}$ . This specific choice of parameters violates the tearing mode stability criterion (Goedbloed et al. 2019), which results in an isolated, unstable tearing mode. Spectra are calculated both for V = 0 (no flow) and V = 0.15, the inclusion of the linear velocity profile in the latter introduces a Doppler shift in the slow and Alfvén continuum, and the results are shown for 501 grid points in Figure 11. Panels (a) and (b) show spectra for V = 0, k_y = 0.49 (no flow) and V = 0.15, k_y = 1.5 (linear flow profile), respectively, revealing the intricate behavior of the damped slow and Alfvén sequences. The purely imaginary unstable tearing mode is annotated with an arrow. The resistivity η is equal to 10⁻⁴ for both cases, yielding a magnetic Reynolds number of ≈10⁴. These results are in perfect agreement with the original spectra depicted in Goedbloed et al. (2019, Figures 14.7–14.9). Note that these cases still have left–right symmetry maintained, despite the presence of equilibrium flow. However, this is purely because the flow profile and domain happen to be chosen in a symmetric way, such that forward and backward modes behave symmetrically.

**Figure 11.** Resistive spectra of the equilibrium given in Equation (42), unstable to the tearing mode (annotated with an arrow). Panels (a) and (b) show the spectrum without and with a linear flow profile, respectively, zoomed in to reveal the complex structure of the slow and Alfvén modes (fast modes are not shown). Panels (c) and (d): growth rate of the tearing mode vs. fixed k_y, varying η and fixed η, varying k_y, respectively. Both bottom panels are cases without flow, with 64 runs of 351 grid points each, per panel. The red line in panels (c) and (d) denotes the theoretical growth rate prediction given in Equation (43).
Download figure:
Standard image High-resolution image

For an excellent discussion on tearing modes, we refer to Goedbloed et al. (2019), where they perform a detailed analysis on the resistive MHD equations to derive an analytical expression for the growth rate of the tearing mode, given by

$\begin{eqnarray}\begin{array}{rcl}\mathrm{Im}(\omega ) & = & {R}_{m}^{-3/5}{\left(\mathrm{KH}\right)}^{2/5}{\left(\displaystyle \frac{{\rm{\Delta }}^{\prime} }{C}\right)}^{4/5}{v}_{A}/a,\mathrm{in}\ \mathrm{which}\\ {\rm{\Delta }}^{\prime} & = & -2\sqrt{{H}^{2}-{K}^{2}}\cot \left(\displaystyle \frac{1}{2}\sqrt{{H}^{2}-{K}^{2}}\right),\end{array}\end{eqnarray} \tag{ 43 }$

with R_m = av_A/η the magnetic Reynolds number, a the width of the slab, and v_A the Alfvén velocity. The other parameters in this equation are variables introduced during the analysis, given by H = αa, K = k₀ a, and C = 2πΓ(3/4)/Γ(1/4) ≈ 2.1236. Panels (c) and (d) of Figure 11 show a comparison of the tearing mode growth rate between the Legolas results (using the no-flow case) and Equation (43). Panel (c) holds k_y = 1.5 constant but varies the resistivity, reaching Reynolds numbers of ≈10⁸, a feat that is nearly impossible to achieve if fully nonlinear codes would be used instead. The correspondence between the theoretical and numerical growth rates is nearly one to one, except for large resistivity values as the theoretical approximation starts to break down in those regimes. In panel (d), we keep η = 10⁻⁶ (magnetic Reynolds number of ≈10⁶) constant and vary k_y between 0.1 and 3.5, which again yields excellent agreement between theory and numerical results. For both panels (c) and (d), we performed 64 runs of 351 grid points each.

3.3.4. Resistive Rippling Modes

As an extension to the tearing modes described above, we take a look at resistive rippling modes, which originate whenever there is a spatially varying resistivity profile as discussed in, for example, Priest (2014). To that extent, we take the same equilibrium and parameters as Equation (42) (with V = 0, so without flow and k_y = 0.49), but impose a hyperbolic tangent profile on the resistivity to let it smoothly drop down to zero near the edges of the domain. This spatial variation in resistivity will excite (unstable) rippling modes in a current-carrying equilibrium that is also unstable to the tearing mode. The resulting spectrum is shown in Figure 12 for 1001 grid points, which is again left–right symmetric. The adopted η(x) profile is given explicitly in Equation (44) and depicted in Figure 12, panel (c). Here, η₀ = 10⁻⁴ denotes the constant resistivity value, and s_L = −0.4 and s_R = 0.4 represent the center of the left and right transition region, respectively, having a width of l = 0.1

$\begin{eqnarray}\eta (x)=\left\{\begin{array}{ll}\displaystyle \frac{{\eta }_{0}}{2}+\displaystyle \frac{{\eta }_{0}}{2\tanh \pi }\tanh \left(\displaystyle \frac{2\pi (x-{s}_{L})}{l}\right) & \mathrm{if}\,{s}_{L}-\displaystyle \frac{l}{2}\leqslant x\leqslant {s}_{L}+\displaystyle \frac{l}{2},\\ {\eta }_{0} & \mathrm{if}\,{s}_{L}+\displaystyle \frac{l}{2}\lt x\lt {s}_{R}-\displaystyle \frac{l}{2},\\ \displaystyle \frac{{\eta }_{0}}{2}+\displaystyle \frac{{\eta }_{0}}{2\tanh \pi }\tanh \left(\displaystyle \frac{2\pi ({s}_{R}-x)}{l}\right) & \mathrm{if}\,{s}_{R}-\displaystyle \frac{l}{2}\leqslant x\leqslant {s}_{R}+\displaystyle \frac{l}{2},\\ 0 & \mathrm{elsewhere}\end{array}\right.\end{eqnarray} \tag{ 44 }$

**Figure 12.** Resistive spectrum using the same parameters as panel (a) in Figure 11, but with a spatially varying η(x) profile (panel (c)). Panel (b): full spectrum with modified fast-mode sequences. Panel (a): zoom in on the slow and Alfvén modes, the inset zooms in on the tearing and unstable rippling modes. Panel (d): v_x (solid) and ρ (dotted) eigenfunctions of the tearing mode (inset, blue). Panels (e) and (g): v_x eigenfunctions of the first four modes associated with the left part of the η(x) profile (inset, red circles). Panels (f) and (h): v_x eigenfunctions of the first four modes associated with the right part of the η(x) profile (inset, green squares).
Download figure:
Standard image High-resolution image

The inclusion of this relatively simple η(x) profile has a major influence on the resulting spectrum: the fast-mode sequences trace out intricate patterns in the complex eigenvalue plane. The semicircles traced out by the slow and Alfvén modes are still vaguely present but are much more scattered than their constant-η counterparts. The inset in panel (a) of Figure 12 zooms in near the origin, revealing that the tearing mode is still present, annotated with an arrow and a large blue dot. Panel (d) shows the v_x and ρ eigenfunctions of the tearing mode, with relatively sharp transitions near the x = 0 point and a (minor) influence from the wings of the η(x) profile. The rippling modes visible in the inset come in two branches, annotated with red circles and green squares, which intersect near the top of the branches. The red circles correspond to the influence of the left wing of the η(x) profile, visible in the eigenfunctions which are shown in panels (e) (first and second dots) and (g) (third and fourth dots), and are very localized near the η transition region as can be seen by the grid indication on the horizontal axis. Analogously, the branch annotated by green squares corresponds to the right wing of the η(x) profile, with v_x eigenfunctions shown in panels (f) and (h).

The fact that the rippling modes appear more unstable than the tearing mode indicates that their importance should not be underestimated and that they will most likely be prominently present when fully realistic η profiles are used. Indeed, in actual plasmas, the temperature variation in the equilibrium can itself already cause a spatially varying resistivity profile, and this is accounted for in Legolas. This allows for future systematic studies of rippling versus tearing mode dominance in slabs or loop-like settings.

3.4. Nonadiabatic, Cylindrical Cases

Because Legolas is the first 1D code to simultaneously include nonadiabatic effects, resistivity, and flow, there are no known previously calculated spectra where all these physical effects are included. Nevertheless, there exist some spectra in the literature where solely nonadiabatic effects are included (and where the equilibrium is static and no resistivity is incorporated), so we will use these as a base comparison for testing the nonadiabatic terms in the implementation.

3.4.1. Nonadiabatic Discrete Alfvén Waves

The first spectrum we will look at is that of nonadiabatic discrete Alfvén waves, described in Keppens et al. (1993). The basic cylindrical equilibrium represents a solar coronal loop, and nonadiabatic effects included were optically thin radiative losses and parallel thermal conduction. It was pointed out how a cluster sequence of discrete Alfvén waves may become unstable and could lead to disruptions or oscillatory behavior in loops of prominences.

This particular equilibrium uses an axial current profile in a cylindrical geometry, yielding the same magnetic configuration as given in Equation (33) although this time with ν = 2 such that a current distribution is present throughout the loop. The pressure profile can be obtained through integration of the equation for magnetostatic equilibrium (10) without flow. As a boundary condition, we impose that the pressure vanish at the plasma boundary r = R, which can be used to constrain the integration constant. Furthermore, a parabolic density profile is used, yielding an equilibrium configuration given by

$\begin{eqnarray}\begin{array}{rcl}{\rho }_{0} & = & 1-\left(1-d\right){\left(\displaystyle \frac{r}{R}\right)}^{2},{B}_{02}=\displaystyle \frac{1}{6}{j}_{0}r\left({r}^{4}-3{r}^{2}+3\right),{B}_{03}=1,\\ {p}_{0} & = & \displaystyle \frac{{j}_{0}^{2}}{720}\left[12\left({R}^{10}-{r}^{10}\right)-75\left({R}^{8}-{r}^{8}\right)+200\left({R}^{6}-{r}^{6}\right)\right.\\ & & -\ \left.270\left({R}^{4}-{r}^{4}\right)+180\left({R}^{2}-{r}^{2}\right)\right],\end{array}\end{eqnarray} \tag{ 45 }$

in which d = 0.2 and j₀ = 0.125. The cylinder wall is taken at R = 1. The equilibrium temperature profile can be derived using T₀ = p₀/ρ₀ and ${T}_{0}^{{\prime} }=\left({p}_{0}^{{\prime} }{\rho }_{0}-{\rho }_{0}^{{\prime} }{p}_{0}\right)/{\rho }_{0}^{2}$ , where the prime denotes the derivative with respect to r. Only conduction parallel to the magnetic field lines is taken into account, because cross-field thermal conduction acts on too long a timescale in this case (Keppens et al. 1993).

The wavenumbers k₂ = m and k₃ = k are taken to be 1 and 0.05, respectively. Optically thin radiative losses are included as described in Equation (5), in which we assume that the heating is such that it exactly balances out the cooling terms in the equilibrium state. We use the cooling curve introduced by Rosner et al. (1978), which represents a piecewise cooling law with predetermined coefficients. Because the inclusion of nonadiabatic effects requires us to specify unit normalizations in order to look up the dimensional values in the cooling tables, we take reference values of 50 G for the magnetic field, 1.5 × 10⁻¹⁵ g cm⁻³ for the density and 10¹⁰ cm as a length scale. This automatically constrains all other normalizations through the ideal gas law, assuming a fully ionized plasma.

The spectrum is calculated using 501 grid points, Figure 13 shows the discrete Alfvén spectrum (panel (a)) along with a plot of the Alfvén continuum (panel (b)). In total, six discrete modes (that is, modes that fall outside of the Alfvén continuum) were found in contrast with the five discrete modes in the original paper (Keppens et al. 1993). The discrete modes are annotated according to their overtones, where ω₁ represents the fundamental mode (FM) and ω₆ the fifth overtone (OT). This last overtone does not show up for resolutions below ∼200 grid points, which explains why it is not present in the original work. The imaginary part of the rv_r eigenfunction corresponding to each of these six modes is shown in the bottom-right panels of Figure 13, where the location of the minimum in the Alfvén continuum is indicated by a red dotted line. Note that the number of eigenfunction nodes increases when considering modes further in the Alfvén sequence: the eigenfunction corresponding to ω₁ has no nodes, and ω₂ through ω₆ have one, two, three, four, and five nodes, respectively. Table 1 shows a detailed comparison between the discrete modes found by Legolas and the ones from Keppens et al. (1993). We multiplied these latter by i, because they use the convention $\exp ({st})$ rather than $\exp (-i\omega t)$ in the Fourier analysis (thus, ω = is).

Table 1. Comparison of Discrete Alfvén Eigenvalues between Legolas and the Original Work

Eigenvalue	`Legolas`	Keppens et al. (1993)
FM ω₁	−0.1044165837 + 1.7260866 × 10⁻⁰⁷ i	−0.10436875 + 1.7 × 10⁻⁰⁷ i
1st OT ω₂	−0.1090793809 + 6.5883040 × 10⁻⁰⁸ i	−0.10907061 + 6.7 × 10⁻⁰⁸ i
2nd OT ω₃	−0.1096033464 + 8.3813211 × 10⁻⁰⁹ i	−0.10960184 + 8.4 × 10⁻⁰⁹ i
3rd OT ω₄	−0.1096779374 + 2.5325093 × 10⁻⁰⁹ i	−0.10967774 + 2.5 × 10⁻⁰⁹ i
4th OT ω₅	−0.1096871498 + 3.7770023 × 10⁻¹⁰ i	−0.10968708 + 4.0 × 10⁻¹⁰ i
5th OT ω₆	−0.1096883410 + 7.4335660 × 10⁻¹¹ i	⋯

Download table as: ASCII Typeset image

In reality, the mode sequence is expected to be an infinite sequence accumulating at the local minimum in the Alfvén continuum, as indicated in panel (b) on Figure 13. Both the real and imaginary parts of the eigenvalues are in excellent agreement. It should be noted that when the nonadiabatic effects are omitted, the imaginary parts of these discrete modes become zero, such that they lie on the real axis, representing stable waves. The inclusion of nonadiabatic effects hence has almost no influence on their oscillation frequency and solely pushes these modes into the unstable part of the imaginary plane. In Keppens et al. (1993), it was pointed out how these discrete Alfvén mode sequences can be studied by means of a WKB analysis. However, this only correctly predicted the damped or overstable nature of the higher-order modes: to determine whether the most global modes of the sequence are overstable or damped requires a full numerical computation, possible with general tools like Legolas.

3.4.2. Magnetothermal Instabilities

As a final test, we look at magnetothermal instabilities, originally depicted in Van der Linden et al. (1992). The geometry is cylindrical with r ∈ [0, 1] and an isothermal equilibrium profile given by

$\begin{eqnarray}\begin{array}{rcl}{T}_{0} & = & 1,\quad {B}_{02}=\displaystyle \frac{r}{1+{r}^{2}},\quad {B}_{03}=0,\\ {p}_{0} & = & \displaystyle \frac{1}{2{\left(1+{r}^{2}\right)}^{2}},\quad {\rho }_{0}=\displaystyle \frac{{p}_{0}}{{T}_{0}},\end{array}\end{eqnarray} \tag{ 46 }$

with k₂ = m = 0 and k₃ = k = 1. There is no B_z, such that this configuration actually represents a z-pinch (which is very unstable), and we are looking at axisymmetric (sausage) m = 0 modes here. Field-aligned thermal conduction is included, cross-field thermal conduction is omitted. Optically thin radiative losses are accounted for, and we use the same cooling curve and heating assumptions as for the nonadiabatic discrete Alfvén waves in Equation (45). Reference values are taken to be 2.6 MK for the temperature, 10 G for the magnetic field, and 10⁸ cm for the length scale. As before, this automatically constrains all other normalizations as well. These parameters are representative for solar coronal loops and arcades.

Both the Alfvén and the slow continuum collapse into marginal (zero) frequency. However, including finite parallel thermal conduction and radiative losses introduces the thermal continuum. Together with the marginal slow-Alfvén frequencies, this organizes the modes such that thermal instabilities merge with magnetic modes, introducing magnetothermal branches. Figure 14 shows the MHD spectrum for 1001 grid points in the region of the magnetothermal branches, where the bottom-right panel (c) zooms in further near the origin. The modes denoted by I₊₁ and T₁ represent fundamental magnetic and thermal modes, respectively, meaning they are not coalesced. The overtones of these modes on the other hand are coalesced; these are denoted by (I₊₂, T₂)⁻ and (I₊₁₅, T₁₅)⁻ for the 1st and 14th overtones, respectively. The minus sign in superscript means that these are located in the negative part of the real axis. There are corresponding modes on the positive part of the real axis, as without flow, we maintain the left–right symmetry of the spectrum. The first 14 overtones are in excellent agreement with those described in Van der Linden et al. (1992) and are encircled in red in the left panel (a). However, the original work only displays solutions up to the 14th overtone. The resolution used here is much higher than the originally published results, allowing us to probe the region near the origin as well. Panel (c) zooms in near the origin, denoted by the dotted rectangle in panel (a), revealing a complex mixture of different branches and scattered modes. It seems that the original magnetothermal branches split, and a dense region covered with many magnetothermal modes appears. The analytically known thermal continuum is shaded in green in the figure, and is represented by a dense, continuous range of discrete eigenvalues on the imaginary axis. Because this is the first time that high resolution is possible for computing the magnetothermal modes, it is left to future work to clarify how the MHD spectrum allows for such complex mode interactions, revealing the possible presence of areas covered by modes in the complex eigenfrequency plane.

**Figure 14.** Magnetothermal instabilities for m = 0 and k = 1. The fundamental modes and 14 overtones discussed in the original work are encircled by red and denoted in the left panel (a). Panel (b) shows the thermal continuum as a function of radius. Panel (c) zooms into the region denoted by dashed lines in panel (a), revealing various other modes forming intricate structures.
Download figure:
Standard image High-resolution image

4. Convergence

As discussed in this paper, increasing the resolution can have a major influence on the spectrum, depending on whether or not all modes are sufficiently resolved. Generally speaking, the further one goes in a specific sequence, that is, looking at larger mode numbers that represent overtones having more and more nodes in their eigenfunctions, the higher the resolution that is required in order to resolve the mode completely. For eigenfunctions, it is usually immediately clear if a mode is resolved or not, because higher mode numbers translate into more oscillations in the eigenfunction. Hence, once the number of oscillations approaches the amount of grid points, the eigenfunction is no longer resolved.

For the eigenfrequencies it is not a priori clear when a specific ω is resolved. One way to constrain an eigenvalue is to do multiple runs, each time increasing the resolution. Once an eigenvalue no longer shifts in the complex plane, it can be considered resolved, and increasing the resolution even further will not have (much) effect on its value. Now, the question of how many grid points are typically needed to be able to speak of a "resolved" spectrum naturally arises. This will strongly depend on the type of equilibrium considered: for smooth equilibria without sharp transitions, one can get away with a few dozen grid points and already reach an acceptable accuracy for most modes. However, in the case of equilibria with large gradients, or even with localized discontinuities (interfaces), one has to make sure that a sufficient number of grid points are taken in order to sufficiently resolve that jump.

Furthermore, the amount of eigenvalues in the spectrum increases with resolution, because there are as many eigenfrequencies as the dimension of the matrices (which is 16 times the number of grid points). It is therefore entirely possible that one starts probing parts of the spectrum that were initially not visible at lower resolutions, simply because there were no eigenvalues in that region for that amount of grid points. An example is given here, where we look back at the Kelvin–Helmholtz equilibrium discussed in Section 3.2.3. Figure 15 shows this equilibrium for the same values as employed earlier, but every panel shows the spectrum for a different resolution, with the amount of grid points given in the top-left corner. The first three panels start out at low resolution, increasing from 26 to 76 grid points. Three modes are annotated in the figure, corresponding to the KH and CD modes discussed earlier. These modes are already decently resolved even at the lowest resolutions and can be considered completely resolved at around 100 grid points. However, the region near the horizontal axis is another matter entirely. At low resolutions (first three panels), we see that most of the modes here are almost randomly scattered in the complex plane. At around 100 grid points, they start lining up and form a more intricate "elliptic" pattern. This is the point where the original paper by Baty & Keppens (2002) stopped, because it was not really feasible to run these codes at higher resolutions, and they were only interested in the most unstable KH and CD modes that were resolved, and that determined the early evolution of a full nonlinear MHD run that they performed.

**Figure 15.** Spectrum of the Kelvin–Helmholtz and current-driven equilibrium as shown earlier in Figure 7, using the same parameters as in Section 3.2.3, at various resolutions (indicated on the top-left of each panel). The inset on the top two rows of panels zooms near the horizontal axis, for vertical axis values between −0.1 and 0.1. The insets on the bottom row of panels have slightly different vertical bounds, from −0.03 to 0.03.
Download figure:
Standard image High-resolution image

When the resolution is increased even further, we see that the near-axis modes shift even closer to the real eigenfrequency axis, and as visible on the insets on each panel, these modes start to form intriguing patterns at high resolutions. The amount of "scattered" modes decreases considerably, and a helix-like structure is formed at very high resolutions (1001 grid points, lower-right panel). These complex and almost geometric patterns may arise in various different equilibria as well and call for extensive research in this topic. A complementary means to identify which modes are actually resolved is to exploit the spectral web; see Goedbloed (2018a, 2018b) and Goedbloed et al. (2019, ch. 12–13), who locate eigenmodes on specific curves and their intersections. These curves relate to the two self-adjoint operators at play in stationary, adiabatic MHD. Only a combined approach using high-resolution Legolas runs, modern linear algebra solvers, and physical insight in MHD spectroscopy will in time reveal the true importance of these modes.

Based on the conclusions drawn here, it is clear that the resolution required depends on the part of the spectrum that is to be investigated. For isolated modes, about 100 grid points is more than sufficient to resolve most of them, although this depends on the specified equilibrium. For large-scale surveys of the spectrum, large resolutions should be employed in order to reveal possible regions of interest.

5. Conclusion and Outlook

In this paper, we introduced the novel 1D finite element code Legolas, meant to tackle the complete MHD spectrum including all kinds of physical effects in order to be able to look at various realistic equilibrium configurations. We performed a Fourier mode analysis on the linearized MHD equations, after which the resulting matrix eigenvalue problem is solved using a finite element approach. A general formalism was introduced to handle both Cartesian and cylindrical geometries. Legolas is the first spectral code to combine the effects of flow, resistivity, gravity, radiative cooling, and anisotropic thermal conduction. This opens the door to novel, in-depth studies of the complete MHD spectrum, which was up to now impossible with existing numerical tools.

We tested Legolas against various previously established results from the literature, looking at the comparison to analytical results as well as to spectra previously obtained by similar numerical codes. Correspondence with existing results is in most cases one to one, greatly increasing the confidence in this new tool. Some cases were run in higher resolutions than their original counterparts, revealing interesting features and additional structure in the spectra. The resistive homogeneous case in Figure 9 is in perfect agreement with results from the literature near the origin, revealing the semicircles traced out by Alfvén and slow modes. However, far along the fast-mode sequence, an intriguing loop appears, locally breaking the Sturmian behavior of the fast modes. Even this case calls for further investigation, as this phenomenon has not been described before as no extremely high-resolution studies of that part of the spectrum have been done to date.

Something similar can be seen when looking at the instabilities of a cylindrical magnetized jet flow in Section 3.2.3, where the outer KH and CD modes are in perfect correspondence with the original results. However, the high resolution revealed complex structures that were not probed before. This becomes even clearer in the convergence study of Section 4, where we used the same setup but increased the resolution even further. It is clear that the spectrum evolves drastically at higher resolutions, when more and more modes are becoming properly resolved. Because a complete knowledge of the MHD spectrum of flowing equilibria is lacking, tools like Legolas and careful converged studies will be essential to unravel the role of as yet unresolved spectral structure. We speculate that these spectral structures have a physical meaning, which can be corroborated in adiabatic flowing cases with a complementary spectral web approach. Our speculation extends to cases that address MHD modes in cylindrical accretion disks (Keppens et al. 2002), where it is becoming clear that many more modes than the celebrated MRI instability enrich the spectral structure.

Also in the nonadiabatic cases, we found something interesting: the outer modes of the magnetothermal sequence are again in perfect correspondence with modes found in the original work. However, near the origin—a region that has not been investigated until now—there are indications of a splitting of the magnetothermal sequence in its thermal and magnetic counterparts with various modes scattered in between. An in-depth study of magnetothermal instabilities is called for, a topic of research that essentially stopped progressing for the last two decades. Now that high-resolution MHD simulations start to reveal the complexity of thermal-instability-driven evolutions (Xia & Keppens 2016; Claes et al. 2020), a revival of this topic is urgently needed.

As a side note, it can be argued that a standard Fourier decomposition as done in all cases discussed in this paper misses out on nonmodal growth, that is, so-called transient modes, which have been shown to be of significant importance in resistive systems (MacTaggart 2018). However, it is the authors' opinion that a calculation of the complete spectrum is all that is needed. Transient growth may well be a consequence of solving an actual initial value problem using a Laplace transformation while taking all discrete and continuous modes into account. A detailed proof is far from evident; however, in adiabatic cases we have self-adjoint operators (one for a static case, two if flow is included), such that it seems plausible that a full decomposition in terms of a complete basis of eigenfunctions is possible. The inclusion of flow makes these eigenfunctions nonorthogonal, while in ideal, static plasmas, the eigenfunctions are orthogonal (Goedbloed et al. 2019, ch. 12). For nonideal, stationary plasmas, nonorthogonality can again be expected. How this decomposition would be affected by the inclusion of nonadiabatic effects such as resistivity or viscosity is not a priori clear. The discussion of transient growth is usually done in the context of nonmodal analysis and uses the concept of pseudo-spectra, that is, also allowing for modes that are nearly eigenvalues.

All of the above is a clear indication that a thorough investigation of the entire spectrum is in order to yield new insights in (M)HD instabilities. For solar applications alone, the effect of thermal conduction on the linear modes as indicated by Van der Linden & Goossens (1991a) has never been investigated in fully realistic setups, possibly allowing for a deeper understanding of fine structure in solar prominences. However, as Legolas is such a versatile code, the possibilities are endless, from various hydrodynamic configurations to more magnetically oriented astronomical cases like accretion disks or astrophysical jets. The discussion of the more "advanced" cases in this paper alone reveals how much of the MHD spectrum is still not thoroughly investigated. We plan to tackle various new cases in future work, combining linear results from this code with state-of-the-art nonlinear numerical simulations. Including additional physical effects is also planned, with an extension to viscous (M)HD (Navier–Stokes), Hall MHD, and ambipolar terms. Additionally, we plan to include "real" atmospheric models, that is, numerically generated hydrostatic equilibria based on a given temperature profile. The possibilities are endless, and modern MHD spectroscopy will no doubt shed new light on various physical phenomena.

The authors would like to thank Hans Goedbloed for useful discussions and suggestions. This work is supported by funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program, grant agreement No. 833251 PROMINENT ERC-ADG 2018; by the VSC (Flemish Supercomputer Center), funded by the Research Foundation—Flanders (FWO) and the Flemish Government—department EWI; and by internal funds KU Leuven, project C14/19/089 TRACESpace.

Appendix: Numerical Approach

A.1. Finite Element Analysis

The system of Equations (19)–(26) actually defines a generalized eigenvalue problem in ω of the form ${ \mathcal A }{\boldsymbol{X}}=\omega { \mathcal B }{\boldsymbol{X}}$ in which the state vector ${\boldsymbol{X}}$ is given by

$\begin{eqnarray}&&{\boldsymbol{X}}\,=\,({\rho }_{1},{v}_{1},{v}_{2},{v}_{3},{T}_{1},{a}_{1},{a}_{2},{a}_{3})\,=\,({x}^{1},{x}^{2},{x}^{3},{x}^{4},{x}^{5},{x}^{6},{x}^{7},{x}^{8}),\end{eqnarray} \tag{ A1 }$

where each superscript on the right-hand side denotes the index of the unknown variable in the state vector. The matrices ${ \mathcal A }$ and ${ \mathcal B }$ contain the various equilibrium quantities and differential operators with respect to u₁. The eigenvalue problem is solved by using the finite element method (FEM). The basic idea is to discretize the domain in (not necessarily equally spaced) subdomains that are separated by nodes, that is, prechosen fixed points x_j with j the node number (0 to N) on the interval x (in the Cartesian case), after which a linear combination of basis functions is used to approximate the unknown variable xⁱ. These basis functions are represented by local piecewise polynomials on every subdomain and vanish outside this subdomain. As shown in, for example, Rappaz (1977), using the same basis functions for all variables leads to spectral pollution, while using a combination of linear and constant finite elements can yield inaccurate eigenvalues (Kerner et al. 1985). Hence, Legolas uses a combination of higher-order finite elements, namely quadratic basis functions for the variables ρ₁, v₂, v₃, T₁, a₁, and cubic Hermite basis functions for the variables v₁, a₂, a₃. This mixture of high-order finite elements ensures that the eigenfunctions can represent certain physical properties of specific eigenfunctions exactly, and can, for example, allow for truly incompressible modes that have ${\rm{\nabla }}\cdot {v}_{1}=0$ everywhere on the domain. The trade-off in using higher-order basis functions is that this not only increases the accuracy of the spectrum, but also increases the size of the final matrices by a factor of 2 compared to linear finite elements. Hence, the approximation of a state variable xⁱ for a domain divided into N subdomains (nodes) is given by

$\begin{eqnarray}&&{x}^{i}({u}_{1})\approx \displaystyle \sum _{j=0}^{N}\left({H}_{j}^{1}({u}_{1}){x}_{j1}^{i}+{H}_{j}^{2}({u}_{1}){x}_{j2}^{i}\right),\end{eqnarray} \tag{ A2 }$

where ${H}_{j}^{1}({u}_{1})$ and ${H}_{j}^{2}({u}_{1})$ denote the quadratic or cubic elements at node j. The Hermite (C_j) and quadratic (Q_j) basis functions themselves are given by

$\begin{eqnarray}\begin{array}{rcl}{Q}_{j}^{1}(x) & = & \left\{\begin{array}{l}\displaystyle \frac{4\left(x-{x}_{j-1}\right)\left({x}_{j}-x\right)}{{\left({x}_{j}-{x}_{j-1}\right)}^{2}},\\ 0,\\ 0,\end{array}\right.\\ {C}_{j}^{1}(x) & = & \left\{\begin{array}{ll}{\left(\displaystyle \frac{x-{x}_{j-1}}{{x}_{j}-{x}_{j-1}}\right)}^{2}\left(3-2\displaystyle \frac{x-{x}_{j-1}}{{x}_{j}-{x}_{j-1}}\right) & \mathrm{for}\,{x}_{j-1}\leqslant x\leqslant {x}_{j},\\ {\left(\displaystyle \frac{{x}_{j+1}-x}{{x}_{j+1}-{x}_{j}}\right)}^{2}\left(3-2\displaystyle \frac{{x}_{j+1}-x}{{x}_{j+1}-{x}_{j}}\right) & \mathrm{for}\,{x}_{j}\leqslant x\leqslant {x}_{j+1},\\ 0 & \mathrm{elsewhere},\end{array}\right.\\ {Q}_{j}^{2}(x) & = & \left\{\begin{array}{l}\displaystyle \frac{\left(2x-{x}_{j}-{x}_{j-1}\right)\left(x-{x}_{j-1}\right)}{{\left({x}_{j}-{x}_{j-1}\right)}^{2}},\\ \displaystyle \frac{\left(2x-{x}_{j+1}-{x}_{j}\right)\left(x-{x}_{j+1}\right)}{{\left({x}_{j+1}-{x}_{j}\right)}^{2}},\\ 0,\end{array}\right.\\ {C}_{j}^{2}(x) & = & \left\{\begin{array}{ll}\left(x-{x}_{j}\right){\left(\displaystyle \frac{x-{x}_{j-1}}{{x}_{j}-{x}_{j-1}}\right)}^{2} & \mathrm{for}\,{x}_{j-1}\leqslant x\leqslant {x}_{j},\\ \left(x-{x}_{j}\right){\left(\displaystyle \frac{{x}_{j+1}-x}{{x}_{j+1}-{x}_{j}}\right)}^{2} & \mathrm{for}\,{x}_{j}\leqslant x\leqslant {x}_{j+1},\\ 0 & \mathrm{elsewhere},\end{array}\right.\end{array}\end{eqnarray} \tag{ A3 }$

as in, for example, Kerner et al. (1998) and Goedbloed et al. (2019), and are shown graphically along with their derivatives in Figure 16. Because we have two basis functions per subdomain, this allows for the approximation of both the original variable and its derivative by differentiating Equation (A2).

To turn the problem algebraic in nature, we use the Galerkin method such that the eigenvalue problem is written as a set of integral equations by multiplying each of the eight equations by an appropriate element of the chosen basis, denoted by h_j, and integrating over the relevant domain. Mathematically, this can be written as

$\begin{eqnarray}&&{\int }_{{\rm{\Omega }}}\left({ \mathcal A }{\boldsymbol{X}}-\omega { \mathcal B }{\boldsymbol{X}}\right){h}_{j}{{du}}_{1}=0.\end{eqnarray} \tag{ A4 }$

However, the matrix ${ \mathcal A }$ contains second-order derivatives with respect to u₁. In order to reduce these to derivatives of first order, hence simplifying the integrals, we make use of the Galerkin weak formulation. This is achieved by performing integration by parts, which introduces additional surface terms that have to be evaluated at the boundaries. These surface terms can in turn be exploited to enforce boundary conditions, which will be discussed in Appendix A.3. Because there are eight unknowns in our eigenvalue system, together with two basis functions and N subdomains, this implies that the final matrix eigenvalue problem will have 16N equations, resulting in a 16N × 16N size matrix. By extension, it becomes clear that using basis functions of even higher order will increase the size of this eigensystem considerably.

A.2. Matrix Assembly

The actual expression for the matrix elements can be obtained by applying Equation (A4) to the system of differential Equations (19)–(26). As an example, we will look at the v₁ component of the linearized continuity Equation (19), corresponding to element (1, 2) in the matrices. This "number" links to the indices of the state vector ${\boldsymbol{X}}$ , because the continuity equation is associated with ρ₁, or index 1 in the state vector (A1), and the v₁ component is associated with index 2. How these indices translate to the actual position in each matrix will be discussed further in this section.

In the finite element representation adopted here, the (ρ₁, v₁) contribution can be expanded as

$\begin{eqnarray}&&\displaystyle \sum _{j=0}^{N}\omega \displaystyle \int \displaystyle \frac{1}{\varepsilon }{h}_{j}^{1}{h}_{k}^{1}{{du}}_{1}=-\displaystyle \sum _{j=0}^{N}\displaystyle \int \displaystyle \frac{1}{\varepsilon }{\rho }_{0}^{{\prime} }{h}_{j}^{1}{h}_{k}^{2}{{du}}_{1}-\ \displaystyle \sum _{j=0}^{N}\displaystyle \int \displaystyle \frac{1}{\varepsilon }{\rho }_{0}{h}_{j}^{1}\displaystyle \frac{{{dh}}_{k}^{2}}{{{du}}_{1}}{{du}}_{1},\end{eqnarray} \tag{ A5 }$

where the u₁ dependence of the equilibrium density ρ₀ and the basis functions h¹ and h² are implied. In this case, h¹ is quadratic, because ρ₁ is associated with a quadratic basis function; analogously h² is cubic. The matrix elements for this particular contribution are hence given by

$\begin{eqnarray}&&{{ \mathcal B }}_{{jk}}(1,1)=\displaystyle \int \displaystyle \frac{1}{\varepsilon }{h}_{j}^{1}{h}_{k}^{1}{{du}}_{1},{{ \mathcal A }}_{{jk}}(1,2)\,=\,-\displaystyle \int \left(\displaystyle \frac{1}{\varepsilon }{\rho }_{0}^{{\prime} }\right){h}_{j}^{1}{h}_{k}^{2}{{du}}_{1}-\displaystyle \int \left(\displaystyle \frac{1}{\varepsilon }{\rho }_{0}\right){h}_{j}^{1}\displaystyle \frac{{{dh}}_{k}^{2}}{{{du}}_{1}}{{du}}_{1}.\end{eqnarray} \tag{ A6 }$

If this reasoning is applied to all equations in the linear system, it follows that ${ \mathcal B }$ will only have equal-number elements, implying that ${ \mathcal B }$ is fully symmetric and real. ${ \mathcal A }$ , on the other hand, will have cross-term elements such that it is, in general, not symmetric. Furthermore, it might be complex, depending on the included physical effects. Terms that contain derivatives of the state vector components, as for example (v₁, ρ₁), which corresponds to element (2, 1), will be integrated by parts. It is exactly this integration that gives rise to the surface terms, that is, terms that do not contain an integral that translates into the natural boundary conditions. These are discussed in the next subsection.

The actual assembly of both matrices ${ \mathcal A }$ and ${ \mathcal B }$ in Legolas is done by sequential iteration over the grid points. From Equation (A3), we see that finite elements have a localized nature around every grid point j, meaning that only the elements associated with a certain region (that is, the grid point j itself and its neighbors j − 1 and j + 1) yield a nonzero contribution. However, it is actually easier implementation-wise to loop over the elements in the interval [x_j−1, x_j] instead of over those in [x_j−1, x_j+1], because then the actual integration of the matrix elements can be done in the same way, independent of whether the basis functions are cubic or quadratic. If only the interval [x_j−1, x_j] is considered, we end up with 16 possible combinations of the basis functions for every grid point. This translates into a 4 × 4 submatrix for every variable, where every one of the 16 elements corresponds to one specific combination of the shape functions. As there are eight variables in the state vector ,this implies a 32 × 32 matrix block for every grid point, hereafter dubbed a "quadblock" as it consists of four 16 × 16 blocks (hereafter called "subblocks"). Every one of those subblocks inside a quadblock corresponds to a quarter section of the aforementioned submatrix, which is a 2 × 2 block in every subblock. Hence, to recap, a 2 × 2 block times eight variables represents a 16 × 16 subblock, of which four combined form a 32 × 32 quadblock for every grid point.

Of course, every matrix element contains one or more integrals that still have to be calculated. Because the coefficient functions are in general complicated expressions depending on u₁, this is done numerically using a four-point Gaussian quadrature for which an integral in the interval [x_j−1, x_j] can be expressed as

$\begin{eqnarray}&&{\displaystyle \int }_{{x}_{j-1}}^{{x}_{j}}f({u}_{1}){{du}}_{1}\approx \displaystyle \frac{1}{2}\left({x}_{j}-{x}_{j-1}\right)\,\displaystyle \sum _{i=1}^{4}{w}_{i}f\left(\displaystyle \frac{1}{2}({x}_{j}-{x}_{j-1}){\xi }_{i}+\displaystyle \frac{1}{2}({x}_{j-1}+{x}_{j})\right),\end{eqnarray} \tag{ A7 }$

where ξ_i and w_i are the evaluation points and weights of the Gaussian integration. These values can be found in various textbooks, as, for example, given in Goedbloed et al. (2019). The function f denotes the integral coefficients, which are essentially the equilibrium quantities and basis functions evaluated in the various evaluation points. This actually implies that every grid interval [x_j−1, x_j] is subdivided into four points, meaning that the equilibrium expressions are probed using 4(N − 1) points rather than N points. The way the matrices are then assembled is thus done on a double-loop basis, where the outer loop iterates over the intervals [x_j−1, x_j] and the inner loop iterates over the four Gaussian points. This inner loop will calculate the basis functions and matrix elements at every point, then multiply the coefficients with the Gaussian weights and finally add them all together in a consistent manner.

Because every quadblock corresponds to the interval [x_j−1, x_j], we still have to account for the contribution of the x_j+1 grid point. This is done by partially overlapping the quadblock in the next grid point with the one from the previous grid point. Figure 17 shows a visual representation of the structure and assembly process for the ${ \mathcal A }$ matrix, using the KH and CD equilibrium discussed in Section 3.2.3 with only six grid points here for the purpose of illustration. The left panel shows the general matrix, highlighting the block-tridiagonal structure. The dashed gray lines denote the 32 × 32 quadblocks of the matrix, and every dot represents a nonzero value. In total, five quadblocks can be distinguished for six grid points, one for every grid interval. The middle panel shows a zoom-in of the quadblock corresponding to the second grid interval as annotated on the left panel, where it can be seen that the top-left corner of this quadblock overlaps with the bottom-right corner of the previous quadblock corresponding to the first grid interval.

**Figure 17.** General assembly and structure of the finite element (complex) ${ \mathcal A }$ matrix. Left: example of a full matrix for six grid points, showing the block-tridiagonal structure where dots represent nonzero values. The middle panel zooms in on one quadblock, showing the dependence of different subblock positions with respect to the variables. The 2 × 2 blocks corresponding to ${ \mathcal A }(2,5)$ (cubic, quadratic) are highlighted. Right panel: 2 × 2 block assembly for a general ${ \mathcal A }$ (cubic, quadratic) matrix element. Cubic elements are shown in blue, quadratic ones in orange. The dotted gray line denotes zero.
Download figure:
Standard image High-resolution image

**Figure 17.** General assembly and structure of the finite element (complex) ${ \mathcal A }$ matrix. Left: example of a full matrix for six grid points, showing the block-tridiagonal structure where dots represent nonzero values. The middle panel zooms in on one quadblock, showing the dependence of different subblock positions with respect to the variables. The 2 × 2 blocks corresponding to ${ \mathcal A }(2,5)$ (cubic, quadratic) are highlighted. Right panel: 2 × 2 block assembly for a general ${ \mathcal A }$ (cubic, quadratic) matrix element. Cubic elements are shown in blue, quadratic ones in orange. The dotted gray line denotes zero.
Download figure:
Standard image High-resolution image

A single quadblock is further divided into four subblocks, with the location of the different state variables annotated on the middle panel of Figure 17. The matrix element ${ \mathcal A }(2,5)$ is highlighted in every subblock, which corresponds to the (v₁, T₁) contribution, which are cubic (v₁) and quadratic (T₁) variables. The 2 × 2 block in the top-left corner of the quadblock corresponds to the top-left 2 × 2 corner of the right panel, as indicated by the background colors. The right panel shows the various combinations of the regular basis functions for the ${ \mathcal A }(2,5)$ contribution, where the blue curve corresponds to the cubic basis functions and the orange curves to the quadratic ones. It should be noted that the right panel shows one specific case, that is, the regular ${h}_{j}^{2}{h}_{k}^{5}$ terms of the matrix element. If, for example, the ${h}_{j}^{2}{h}_{k}^{7}$ element is calculated, corresponding to the a₂ term in Equation (20), the quadratic basis functions have to be replaced by their cubic counterparts, because a₂ is also a cubic variable. A similar reasoning can be made for the other matrix elements.

The term ${{pH}}_{j}^{\alpha \beta }{H}_{j}^{\alpha \beta }$ at the top of every subpanel on the right of Figure 17 denotes which combination of the basis functions should be used, where p stands for the integral coefficients. The exponent αβ refers to the expressions for the basis functions in Equation (A3), where α = 1 stands for ${Q}_{j}^{1}$ or ${C}_{j}^{1}$ and α = 2 stands for ${Q}_{j}^{2}$ or ${C}_{j}^{2}$ , depending on the variable under consideration. β stands for which part of the basis function should be taken, β = 1 means the first equation in cases, while β = 2 means the second one. As an example, we can look at ${H}_{j}^{12}{H}_{j}^{21}$ : because ${ \mathcal A }(2,5)$ represents a cubic and quadratic variable, this translates into ${C}_{j}^{12}{Q}_{j}^{21}$ . For the cubic part, we therefore take the second equation of ${C}_{j}^{1}$ , corresponding to x_j ≤ x ≤ x_j+1. The quadratic part on the other hand is given by ${Q}_{j}^{21}$ , which means that we take the first equation of ${Q}_{j}^{2}$ , corresponding to x_j−1 ≤ x ≤ x_j. Boundary conditions are imposed after matrix assembly is completed.

A.3. Implementation of Boundary Conditions

Integration by parts on Equation (A4) gives rise to additional surface terms, which should be evaluated at the boundaries. These kinds of conditions are called natural boundary conditions, as they emerge in a natural way by rewriting the eigenvalue problem. The regularity and fixed wall conditions considered earlier on the other hand are called essential boundary conditions and have to be handled explicitly. Because the additional surface terms originate from reducing second-order derivatives to first-order derivatives, we only have these terms for the variables v₁, T₁, a₂, and a₃, as these are the only equations that contain derivatives of higher order. Hence, for the momentum equation, the additional surface terms can be written as

$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{S}_{{v}_{1}} & = & \left[\displaystyle \frac{{T}_{0}}{\varepsilon }{h}_{j}^{2}{h}_{k}^{1}+\displaystyle \frac{{\rho }_{0}}{\varepsilon }{h}_{j}^{2}{h}_{k}^{5}+\left({B}_{02}{k}_{3}-\displaystyle \frac{1}{\varepsilon }{B}_{03}{k}_{2}\right){h}_{j}^{2}{h}_{k}^{6}\right.+{\left.\displaystyle \frac{{B}_{03}}{\varepsilon }{h}_{j}^{2}\displaystyle \frac{{dh}_{j}^{7}}{{du}_{1}}-{B}_{02}{h}_{j}^{2}\displaystyle \frac{{dh}_{k}^{8}}{{du}_{1}}\right]}_{{\rm{\partial }}{\rm{\Omega }}}\\ & = & \,{v}_{1}\left[\displaystyle \frac{{T}_{0}}{\varepsilon }{\rho }_{1}+\displaystyle \frac{{\rho }_{0}}{\varepsilon }{T}_{1}+\left({B}_{02}{k}_{3}-\displaystyle \frac{1}{\varepsilon }{B}_{03}{k}_{2}\right){a}_{1}\right.+{\left.\displaystyle \frac{{B}_{03}}{\varepsilon }{a}_{2}^{{\rm{{\prime} }}}-{B}_{02}{a}_{3}^{{\rm{{\prime} }}}\right]}_{{\rm{\partial }}{\rm{\Omega }}},\end{array}\end{array}\end{eqnarray} \tag{ A8 }$

where the number in the superscript on h_j/k denotes the index of the variable in the state vector ${\boldsymbol{X}}$ . The subscript ∂Ω means that these terms should be evaluated at the left or inner edge, as well as at the right or outer edge. In a similar manner, the surface terms for the energy and induction equations are given by

$\begin{eqnarray}&&\begin{array}{l}{S}_{{T}_{1}}=\displaystyle \frac{i(\gamma -1)}{\varepsilon }{T}_{1}\left[{T}_{0}^{{\prime} }\displaystyle \frac{\partial {\kappa }_{\perp }}{\partial \rho }{\rho }_{1}+\left({T}_{0}^{{\prime} }\displaystyle \frac{\partial {\kappa }_{\perp }}{\partial T}-\displaystyle \frac{\varepsilon ^{\prime} }{\varepsilon }{\kappa }_{\perp }\right){T}_{1}+{\kappa }_{\perp }{T}_{1}^{{\prime} }\right.\\\,+\ 2\left({T}_{0}^{{\prime} }\left(\varepsilon {B}_{02}{k}_{3}-{B}_{03}{k}_{2}\right)\displaystyle \frac{\partial {\kappa }_{\perp }}{\partial ({B}^{2})}-\eta {B}_{03}^{{\prime} }{k}_{2}+\eta {k}_{3}(\varepsilon {B}_{02})^{\prime} \right){a}_{1}\\ \,+\ 2\left({T}_{0}^{{\prime} }{B}_{03}\displaystyle \frac{\partial {\kappa }_{\perp }}{\partial ({B}^{2})}+\eta {B}_{03}^{{\prime} }\right){a}_{2}^{{\prime} }-\ \left.2\left(\varepsilon {T}_{0}^{{\prime} }{B}_{02}\displaystyle \frac{\partial {\kappa }_{\perp }}{\partial ({B}^{2})}+\eta (\varepsilon {B}_{02})^{\prime} \right){a}_{3}^{{\prime} }\right],\end{array}\end{eqnarray} \tag{ A9 }$

$\begin{eqnarray}&&{S}_{{a}_{2}}=\displaystyle \frac{i\eta }{\varepsilon }{a}_{2}\left(-{k}_{2}{a}_{1}+{a}_{2}^{{\prime} }\right),\end{eqnarray} \tag{ A10 }$

$\begin{eqnarray}&&{S}_{{a}_{3}}=i\eta {a}_{3}\left(-{k}_{3}{a}_{1}+{a}_{3}^{{\prime} }\right),\end{eqnarray} \tag{ A11 }$

which should all be evaluated at the boundaries. For the case of a solid wall, we see that the natural boundaries simplify considerably, because if v₁, a₂, and a₃ are all zero at the wall, ${S}_{{v}_{1}}$ , ${S}_{{a}_{2}}$ and ${S}_{{a}_{3}}$ are also zero and have to be omitted. The natural boundary condition on T₁ is only relevant if resistivity or perpendicular thermal conduction is included. However, in the case of the latter, the additional essential boundary condition requires that T₁ = 0, in which case ${S}_{{T}_{1}}$ drops out as well. The only combination in which the surface terms for the energy equation are nonzero is when resistivity is included but perpendicular thermal conduction is omitted. In that case, the resistive heating terms should be included in the calculation, which is done by adding the appropriate terms to the matrix elements A(5, 6), A(5, 7) and A(5, 8).

Currently, only fixed wall boundary conditions are implemented in Legolas. However, the surface terms described here can be used to impose other types of boundary conditions as well. In the case of a plasma–vacuum–wall transition, we have, for example, Bessel functions at the outer boundary of a cylindrical geometry (Roberts 2019), which encode the analytic vacuum solution for the electromagnetic field in the outer vacuum region. These expressions can then be used to rewrite the surface terms (A8)–(A11) in an appropriate way such that they can be added to their respective subblock positions in the matrix. This is a planned extension to be included in future versions of Legolas. This functionality was previously available in some LEDA versions (Van der Linden et al. 1992).

The essential boundary conditions as described in Section 2.3 have to be implemented explicitly. This is done by omitting the relevant basis functions that do not satisfy the boundary conditions on the edges. Consider as an example the variable v₁, which is associated with a cubic basis function. From Figure 16, we see that the only cubic element that is nonzero at the left boundary is ${C}_{j}^{12}$ , which implies that the matrix elements where it appears should be zeroed out. Looking back at how the quadblock is composed in Figure 17, the boundary condition v₁ = 0 corresponds to forcing the odd rows and columns of the v₁ contribution to zero, for subblocks 1, 2, and 3 on the left boundary (that is, the first node j) because these correspond to the 2 × 2 blocks in blue, green and red on the right panel. Similarly, on the right side (viz. the last node j), only ${C}_{j}^{11}$ is nonzero, which implies that the odd rows and columns of subblocks 2, 3, and 4 have to be zeroed out (corresponding to green, red, and yellow). Extending this reasoning to the essential boundary condition on T₁, we see that in this case the even rows and columns should be handled because T₁ is associated with a quadratic element. This is done for both matrices.

Of course, "just" zeroing out rows and columns in a matrix has the unpleasant side effect that the matrices become singular. For the ${ \mathcal A }$ matrix, this is not necessarily a problem; however, the ${ \mathcal B }$ matrix can never be singular because it is inverted when solving the general eigenvalue problem using the QR algorithm. Therefore, we introduce a one on ${ \mathcal B }$ 's diagonal at the location that was zeroed out, and an element δ on ${ \mathcal A }$ 's diagonal. The effect of this is that one essentially "forces" the boundary condition, because this implies that $\delta {x}_{j}^{i}=\omega {x}_{j}^{i}$ . By extension, if δ is taken to be a large number (we take δ = 10²⁰), this means that ${x}_{j}^{i}$ = 0, which corresponds to the essential boundary condition we wanted to impose. The only side effect of this approach is that it introduces eigenvalues equal to δ. However, because δ is taken to be large, these will not influence the spectrum in any way, and they can be easily filtered out during postprocessing. This method thus provides a relatively easy and straightforward way to impose Dirichlet boundary conditions at the edges. The imposed boundary conditions are noticeable on the left panel of Figure 17, especially for the first node. The odd rows and columns that were zeroed out can clearly be seen, together with the large numbers introduced on the main diagonal.

Legolas: A Modern Tool for Magnetohydrodynamic Spectroscopy

Article metrics

Permissions

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Problem Description and Model Equations

2.1. Equilibrium Conditions

2.2. Linearized Equations

2.3. Boundary Conditions

2.4. Solving the Equations

3. Results

3.1. Cartesian Cases: Waves in Stratified Atmospheres

3.1.1. Gravito-MHD Waves

3.1.2. Quasi-Parker Instabilities

3.2. Adiabatic, Cylindrical Cases

3.2.1. Magnetic Flux Tubes

3.2.2. Tokamak Constant Current

3.2.3. KH and CD Instabilities

3.2.4. Suydam Cluster Modes

3.3. Resistive, Cartesian Cases

3.3.1. Resistive Homogeneous Plasma

3.3.2. Quasi-modes in Resistive MHD

3.3.3. Resistive Tearing Modes

3.3.4. Resistive Rippling Modes

3.4. Nonadiabatic, Cylindrical Cases

3.4.1. Nonadiabatic Discrete Alfvén Waves

3.4.2. Magnetothermal Instabilities

4. Convergence

5. Conclusion and Outlook

Appendix: Numerical Approach

A.1. Finite Element Analysis

A.2. Matrix Assembly

A.3. Implementation of Boundary Conditions

Legolas: A Modern Tool for Magnetohydrodynamic Spectroscopy

Article metrics

Permissions

Share this article

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Problem Description and Model Equations

2.1. Equilibrium Conditions

2.2. Linearized Equations

2.3. Boundary Conditions

2.4. Solving the Equations

3. Results

3.1. Cartesian Cases: Waves in Stratified Atmospheres

3.1.1. Gravito-MHD Waves

3.1.2. Quasi-Parker Instabilities

3.2. Adiabatic, Cylindrical Cases

3.2.1. Magnetic Flux Tubes

3.2.2. Tokamak Constant Current

3.2.3. KH and CD Instabilities

3.2.4. Suydam Cluster Modes

3.3. Resistive, Cartesian Cases

3.3.1. Resistive Homogeneous Plasma

3.3.2. Quasi-modes in Resistive MHD

3.3.3. Resistive Tearing Modes

3.3.4. Resistive Rippling Modes

3.4. Nonadiabatic, Cylindrical Cases

3.4.1. Nonadiabatic Discrete Alfvén Waves

3.4.2. Magnetothermal Instabilities

4. Convergence

5. Conclusion and Outlook

Appendix: Numerical Approach

A.1. Finite Element Analysis

A.2. Matrix Assembly

A.3. Implementation of Boundary Conditions