© The Authors 2024

1. Introduction

The solar atmosphere is a very dynamic environment. Reconnection is the cause of energetic phenomena such as solar flares and jets. Jets may be represented as fast, rotating outflows (Shen 2021). This transient flow of plasma has been observed all around the solar corona (Chitta et al. 2021; Long et al. 2023). Waves and quasi-periodic flows are also common (De Moortel et al. 2015; Banerjee et al. 2021). Upflows from the transition layer and corona are considered the driving force behind the slow solar wind (Barczynski et al. 2021a, 2023). Furthermore, prominences are involved in violent eruptions, creating shock waves throughout the corona and launching coronal mass ejections (CMEs) into space (Webb 2015). Prominences are host to a variety of flows (Kucera 2015).

A solar tornadoes is a kind of prominence with plasma that appears to be rotating around a vertical axis (Pettit 1932). The phenomenon has recently been the topic of the review paper by Gunár et al. (2023). Many observations were made in the last decade (see e.g. Su et al. 2012, 2014; Li et al. 2012; Wedemeyer et al. 2013; Yang et al. 2018). Line-of-sight Doppler velocities were measured by Orozco Suárez et al. (2012). The red- and blueshifted pattern was interpreted as rotation around a vertical axis. Yang et al. (2018) observed two tornadoes, both with coherent stable red and blue Doppler shifts, in favour of this interpretation. Models of helical magnetic structures with cool plasma flowing helically were constructed by Luna et al. (2015) and Onishchenko et al. (2018). For the tornado studied by Li et al. (2012), the rotational axis is believed to be parallel to the axis of a horizontal flux rope. Such a rotational flow in a prominence cavity has recently been studied in a multidimensional simulation of a flux rope prominence (Liakh & Keppens 2023). Prominences are dense structures that are assumed to be suspended in the corona by the magnetic field. Luna et al. (2015) showed that prominence mass can be supported by a vertical helical field if the magnetic field is highly twisted and/or significant poloidal upflows are present. As the rotating nature of tornado prominences is highly controversial, counterarguments and interpretations are provided (Gunár et al. 2023). Panasenco et al. (2014) argued that counterstreaming flows (Schmieder et al. 1991) and oscillations on horizontal threads cause the illusion of rotation. The projection of the 3D structure of a horizontal prominence onto the 2D plane of view can make the field lines look elliptical (Gunár et al. 2018).

Local thermal instability is the preferred method to form cool condensations in hot media, such as prominences in the solar corona. Parker (1953) found that the increase in radiative losses with decreasing temperature can lead to runaway cooling. The theory of thermal instability was first described by Field (1965). Claes & Keppens (2019), Claes et al. (2020b), and Hermans & Keppens (2021) performed multidimensional simulations of condensation formation by the thermal mode in a coronal volume.

The spectral study of the MHD waves and instabilities is a fundamental methodology to understand plasmas and their evolution. Frequencies of the eigenmodes can be represented as a spectrum on the complex frequency plane. Depending on the chosen wave-form, the modes with a positive imaginary part of the frequency are unstable and lead to instabilities. The general structure of the eigenvalue spectrum of a cylindrical equilibrium consists of two continua (Appert et al. 1974), the Alfvén and the slow continua, due to the inhomogeneity and radial dependence of the background variables. Discrete modes can arise that cluster towards the edges of the continua (Goedbloed & Sakanaka 1974; Goedbloed 1984). Due to the complexity it is challenging to interpret the complete MHD spectrum of waves and instabilities for equilibria with both axial and azimuthal flow. The first systematic studies of the spectrum of cylindrical plasmas with flow were performed by Hameiri (1981); Hameiri (1983). Bondeson et al. (1987) examined the MHD spectrum for cylindrical plasmas with axial flow and obtained local criteria for clustering of discrete modes, whereas Wang et al. (2004) included both azimuthal and axial flow. The effect of rotational flows on magnetoacoustic surface and body modes of an adiabatic photospheric magnetic flux tube was recently studied by Skirvin et al. (2023). The existence of the thermal continuum, analogous to the Alfvén and slow continua, was demonstrated by van der Linden & Goossens (1991). This arises for an equilibrium with radial variation when non-adiabatic terms are included. They also derived approximate expressions to predict the occurance of discrete thermal modes using a Wentzel-Kramers-Brillouin (WKB) analysis, similar to Goedbloed (1984). However, background flow is typically neglected in spectroscopic studies of thermal instability.

In recent years the numerical study of MHD spectroscopy has been revived by the development of the eigenvalue code Legolas (Claes et al. 2020a; De Jonghe et al. 2022; Claes & Keppens 2023). It is an open-source code used to solve the eigenvalue problem of the linearised and Fourier-analysed MHD equations. It provides the eigenfunctions and the corresponding eigenfrequencies, which can be represented as a spectrum on the complex plane. Claes & Keppens (2021) used Legolas to study the eigenmodes of a model solar atmosphere, focusing on the thermal modes. A large-scale systematic investigation of the scaling of the tearing growth rate with resistivity, density, and velocity parameters was recently completed using Legolas (De Jonghe 2023).

In Sect. 2 we extend the work of van der Linden & Goossens (1991) by including a generic background flow in a cylindrical coordinate system. The non-adiabatic MHD equations are linearised, Fourier-analysed, and transformed into two coupled first-order differential equations. These are examined to understand how a background flow changes the continua. The numerical investigation is presented in Sect. 3. We set up a benchmark equilibrium and explore its eigenvalue spectrum using Legolas. The benchmark has a force-free Gold-Hoyle magnetic field and helical flow field. The derived analytic expressions of the previous section are verified. Several equilibrium parameters are varied to study the effect on the eigenmodes. In Sect. 4 discrete thermal modes are the topic of interest. An approximate expression is derived using a WKB analysis. This expression can be used to predict the appearance of discrete thermal modes based on the equilibrium parameters. The modes are then investigated numerically and the previously obtained expression is tested. Lastly, the eigenfunctions of the discrete eigenmodes are visualised. We present a summary of the results, discussion, and conclusions in Sect. 5.

2. Analytic investigation

In this section we derive the expression for the thermal continuum following closely the work by van der Linden & Goossens (1991). The additional physical effect that we take into account is a background flow. The solar tornadoes that we are interested in are vertical column-like structures. A cylindrical coordinate system is thus the most convenient way to describe the plasma. In order to study thermal instability radiative losses need to be included. We also include a background heating to match the losses and obtain a thermal equilibrium to investigate. Since thermal conduction is a relevant and omnipresent effect in the solar corona, it is included. The influence of thermal conduction on thermal instability has been shown in many works (Field 1965; van der Linden & Goossens 1991; Sharma et al. 2010). For simplicity we only include parallel conduction. The inclusion of perpendicular conduction can modify the thermal modes drastically (van der Linden & Goossens 1991). In fact, perpendicular conduction replaces the continuum with a closely packed set of quasi-continuum modes that are in essence discrete. The eigenfunctions of those quasi-continuum modes might have a localised, rapidly spatially varying stucture, providing us with a natural explanation for finestructure in solar prominences. Since the thermal continuum does not strictly exist with perpendicular conduction included, we do not include it here. We also do not include gravity, as it is not in the radial direction when considering a vertical cylinder.

To start, we consider a fully ionised plasma under influence of aforementioned non-adiabatic effects. The following set of dimensionless MHD equations are appropriate to describe this medium:

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}& =-\mathrm{\nabla }·(\mathit{v}\rho ),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \rho \frac{\partial \mathit{v}}{\partial t}& =-\mathrm{\nabla }p-\rho \mathit{v}·\mathrm{\nabla }\mathit{v}+(\mathrm{\nabla }\times \mathit{B})\times \mathit{B},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \rho \frac{\partial T}{\partial t}& =-\rho \mathit{v}·\mathrm{\nabla }T-(\gamma -1)p\mathrm{\nabla }·\mathit{v}-(\gamma -1)\rho \mathcal{L}\hfill \\ \hfill & +(\gamma -1)\mathrm{\nabla }·(\mathit{\kappa }·\mathrm{\nabla }T),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{B}}{\partial t}& =-\mathit{B}(\mathrm{\nabla }·\mathit{v})+(\mathit{B}·\mathrm{\nabla })\mathit{v}-(\mathit{v}·\mathrm{\nabla })\mathit{B},\hfill \end{array}$$(4)

with the magnetic field in units, where μ₀ = 1. The quantities ρ, T, p, v, and B are the density, temperature, pressure, velocity, and magnetic field, respectively. The function ℒ is the net heat-loss function and is defined as

$$\begin{array}{c}\hfill \mathcal{L}=\rho \mathrm{\Lambda }(T)-\mathcal{H}.\end{array}$$(5)

It represents the difference in energy loss between the radiative cooling, for which we assume an optically thin plasma, and an unknown heating function ℋ. The optically thin cooling function or curve, Λ(T), describes radiative energy loss in function of the temperature. A wide variety is used in the literature. The effect of the cooling curve on numerical simulations of thermal instability is discussed in Hermans & Keppens (2021).

The last term of the right-hand side of the energy equation is related to thermal conduction. κ is the thermal conduction tensor. In magnetised plasmas thermal conduction is anisotropic and typically respresented as

$$\begin{array}{c}\hfill \mathit{\kappa }={\kappa }_{\Vert }{\mathit{e}}_B{\mathit{e}}_B+{\kappa }_{\perp }(\mathit{I}-{\mathit{e}}_B{\mathit{e}}_B),\end{array}$$(6)

where I denotes the unit tensor and e_B = B/B is the unit vector along the magnetic field lines. The coefficients κ_∥ and κ_⊥ denote the conductivity coefficients parallel and perpendicular to the local direction of the magnetic field, respectively. Explicit formulae in terms of temperature, density, and magnetic field strength for astrophysical applications can be found in Spitzer (1962). For the solar corona the perpendicular coefficient is typically twelve orders of magnitude smaller than the parallel coefficient (Braginskii 1965). Hence, as in most works we therefore neglect perpendicular thermal conduction (i.e. we take κ_⊥ = 0). This also simplifies the analytic derivations.

We consider an equilibrium around which we linearise the equations. As our main interest is helical flows and magnetic fields, working in a cylindrical coordinate system is most beneficial. We use the convential meaning of the variables, (r, θ, z) being the radial coordinate, azimuthal angle and axial coordinate, respectively. We consider an axisymmetric tube with all background parameters varying in only one dimension, the radial direction. The background quantities are denoted with a subscript ‘0’. The ideal gas law defines the relation between the background pressure, density and, temperature as

$$\begin{array}{c}\hfill p_0(r)={\rho }_0(r)T_0(r)\end{array}$$(7)

in the dimensionless notation adopted in this work.

We do include a background flow. This flow is stationary (i.e. not varying in time). The background profiles of the equilibrium flow are given by

$$\begin{array}{c}\hfill {\mathit{v}}_0=v_{0\theta }(r){\mathit{e}}_{\theta }+v_{0z}(r){\mathit{e}}_z.\end{array}$$(8)

The background magnetic field is given by

$$\begin{array}{c}\hfill {\mathit{B}}_0=B_{0\theta }(r){\mathit{e}}_{\theta }+B_{0z}(r){\mathit{e}}_z.\end{array}$$(9)

The divergence-free condition of the magnetic field, ∇ ⋅ B = 0, is automatically satisfied for our background magnetic field because the non-radial components only depend on the radial coordinate and B_r = 0. Similarly, the equilibrium flow is incompressible, but compressibility effects are fully accounted for in the linearised equations.

Furthermore, the background profiles are at equilibrium, meaning that they have to satisfy the mechanical and thermal equilibrium. In order to obtain mechanical equilibrium the following radial force balance needs to be satisfied for the equilibrium profiles.

$$\begin{array}{c}\hfill D(p_0+\frac{1}{2}{B}_0^2)=-\frac{B_{0\theta }^2}{r}+\frac{{\rho }_0{v}_{0\theta }^2}{r}.\end{array}$$(10)

We denote the derivative in the direction of r (i.e. $\frac{\partial }{\partial r}$) as D. Hence, we follow the notation of van der Linden & Goossens (1991). Compared to static equilibria, such as in the case of van der Linden & Goossens (1991), the background flow introduces an additional term on the right-hand-side. This term corresponds to the centrifugal force and is only apparent in a cylindrical configuration. Equilibria with rotational flow and static equilibria are thus inherently slightly different. The axial flow profile v_0z is of no relevance to the mechanical equilibrium and can be chosen freely.

Since we are interested in thermal instabilities, we require an initial state in thermal equilibrium. Under the assumption of neglected perpendicular conduction, this can be written as

$$\begin{array}{c}\hfill {\rho }_0{\mathcal{L}}_0=0,\end{array}$$(11)

where ℒ₀ is the heat-loss function at the initial state. Intuitively this condition means that the energy-loss by the optically thin radiative cooling needs to be balanced by the background heating at t = 0. As an equation this means

$$\begin{array}{c}\hfill {\mathcal{H}}_0={\rho }_0\mathrm{\Lambda }(T_0).\end{array}$$(12)

We consider the Eulerian perturbations to be small, so we can expand around the equilibrium in the form

$$\begin{array}{c}\hfill f(\mathit{r},t)=f_0(r)+f_1(\mathit{r},t),\end{array}$$(13)

where f represents any of the physical quantities and the subscript ‘1’ is used for the dynamic, perturbed part. The equilibrium is constant in time, such that temporal derivatives vanish. We are only interested in linear perturbations and can ignore all higher-order terms. Substituting these expansions into the MHD equations, Eqs. (1)–(3), yields

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_1}{\partial t}& =-\mathrm{\nabla }·({\rho }_0{\mathit{v}}_1)-\mathrm{\nabla }·({\rho }_1{\mathit{v}}_0),\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill {\rho }_0\frac{\partial {\mathit{v}}_1}{\partial t}& =-\mathrm{\nabla }{p}_1-{\rho }_0{\mathit{v}}_0·\mathrm{\nabla }{\mathit{v}}_1-{\rho }_0{\mathit{v}}_1·\mathrm{\nabla }{\mathit{v}}_0-{\rho }_1{\mathit{v}}_0·\mathrm{\nabla }{\mathit{v}}_0\hfill \\ \hfill & +(\mathrm{\nabla }\times {\mathit{B}}_1)\times {\mathit{B}}_0+(\mathrm{\nabla }\times {\mathit{B}}_0)\times {\mathit{B}}_1,\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill {\rho }_0\frac{\partial T_1}{\partial t}& =-{\rho }_0{\mathit{v}}_0·\mathrm{\nabla }{T}_1-{\rho }_0{\mathit{v}}_1·\mathrm{\nabla }{T}_0-(\gamma -1)p_0\mathrm{\nabla }·{\mathit{v}}_1\hfill \\ \hfill & -(\gamma -1){\rho }_0[{\rho }_1{\mathcal{L}}_{\rho }+T_1{\mathcal{L}}_T]\hfill \\ \hfill & +(\gamma -1){\mathit{B}}_0·\mathrm{\nabla }[\frac{{\kappa }_{\Vert ,0}}{B_0^2}{\mathit{B}}_0·\mathrm{\nabla }{T}_1]\hfill \\ \hfill & +(\gamma -1){\mathit{B}}_0·\mathrm{\nabla }[\frac{{\kappa }_{\Vert ,0}}{B_0^2}{\mathit{B}}_1·\mathrm{\nabla }{T}_0],\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill \frac{\partial {\mathit{B}}_1}{\partial t}& =({\mathit{B}}_0·\mathrm{\nabla }){\mathit{v}}_1+({\mathit{B}}_1·\mathrm{\nabla }){\mathit{v}}_0-{\mathit{B}}_0(\mathrm{\nabla }·{\mathit{v}}_1)-({\mathit{v}}_0·\mathrm{\nabla }){\mathit{B}}_1\hfill \\ \hfill & -({\mathit{v}}_1·\mathrm{\nabla }){\mathit{B}}_0.\hfill \end{array}$$(17)

We used the equilibrium conditions of mechanical and thermal equilibrium, together with the fact that the background is stationary, $\frac{\partial f_0}{\partial t}=0$.

The quantities ℒ_ρ and ℒ_T are the partial derivatives of the heat-loss function with respect to density and temperature, respectively. Hence, they are given by

$$\begin{array}{c}\hfill {\mathcal{L}}_{\rho }={\left(\frac{\partial \mathcal{L}}{\partial \rho }\right)}_T,{\mathcal{L}}_T={\left(\frac{\partial \mathcal{L}}{\partial T}\right)}_{\rho },\end{array}$$(18)

to be evaluated at the equilibrium values. These derivatives are important ingredients in the instability criteria of Field (1965) for the thermal instability. The derivative with respect to the density is usually positive for astrophysical plasmas. This is the case for the optically thin plasma in the solar corona. The derivative with respect to the temperature is determined by the shape of the cooling curve, as the latter is a function of temperature. For most cooling curves discussed in this work the cooling curve decreases with increasing temperature around coronal values (see e.g. Hermans & Keppens 2021). The derivative with respect to temperature is thus negative and this leads to thermal instability. Although, thermal conduction can stabilise the medium.

Neglecting perpendicular thermal conduction simplifies the linearisation of the thermal conduction term. To linearise the thermal conduction term, we made use of the definition Eq. (6) and the vector calculus identity

$$\begin{array}{c}\hfill \mathrm{\nabla }·({\mathit{B}}_{0}f)=f(\mathrm{\nabla }·{\mathit{B}}_0)+{\mathit{B}}_0·(\mathrm{\nabla }f).\end{array}$$(19)

The first term on the right-handside vanishes because of the divergence-free background magnetic field. We also used the fact B_r = 0, which implies that

$$\begin{array}{c}\hfill {\mathit{B}}_0·\mathrm{\nabla }{T}_0=0,\end{array}$$(20)

to have no terms with the perturbed parallel conduction coefficient present. Therefore, we can drop the subscript ‘0’ of the parallel conduction coefficient.

Since the equilibrium quantities only depend on the radial coordinates, we Fourier-analyse the perturbed quantities with respect to the ignorable coordinates, being θ and z. One typically introduces the wavenumber m and k in the azimuthal and axial directions, respectively. The time dependence is assumed exponential with a complex frequency ω. The perturbed quantities can then be written as plane-waves

$$\begin{array}{c}\hfill f_1(\mathit{r},t)=f^{\prime }(r)e^{i(m\theta +kz-\omega t)},\end{array}$$(21)

with the primed f′(r) being the amplitudes of the waves. Note that in the work by van der Linden & Goossens (1991) the notation of the time dependence and frequency are slightly different. They use the exponential as e^st, hence for direct comparison one should use the substitution s = −iω or ω = is. This has as most important effect that the real and complex axis are swapped. Purely imaginary modes in the notation used in this work, such as the thermal modes, are denoted as purely real modes in van der Linden’s work.

We introduce an additional quantity, the total pressure perturbation as

$$\begin{array}{c}\hfill Y=p^{\prime }+B_{0\theta }{B}_{\theta }^{\prime }+B_{0z}{B}_z^{\prime },\end{array}$$(22)

which allows us to reduce the set of eight linear MHD equations to two first order differential equations. Furthermore, we use the linearised ideal gas law to define the pressure perturbation.

$$\begin{array}{c}\hfill \frac{p^{\prime }}{p_0}=\frac{{\rho }^{\prime }}{{\rho }_0}+\frac{T^{\prime }}{T_0}.\end{array}$$(23)

Since equilibrium values are always denoted with a subscript ‘0’, we can drop the prime on the perturbed quantities. The Fourier-analysed MHD equations are given by

$$\begin{array}{cc}\hfill [\omega -{\mathrm{\Omega }}_0]\rho & =-\frac{i}{r}D(r{\rho }_0{v}_r)+{\rho }_0[\frac{m}{r}{v}_{\theta }+k{v}_z],\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\rho }_0{v}_r& =-iD(Y)-F{B}_r-\frac{2i{B}_{0\theta }}{r}{B}_{\theta }+\frac{i{v}_{0\theta }^2}{r}\rho \hfill \\ \hfill & +\frac{2i{\rho }_0{v}_{0\theta }}{r}{v}_{\theta },\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill {\rho }_{0}r{v}_{\theta }& =mY-rF{B}_{\theta }+iD(r{B}_{0\theta })B_r-i{\rho }_{0}D(r{v}_{0\theta })v_r,\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill {\rho }_0{v}_z& =kY-F{B}_z+iD(B_{0z})B_r-i{\rho }_{0}D(v_{0z})v_r,\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill [[\omega -{\mathrm{\Omega }}_0]{\rho }_0& +i(\gamma -1)[\frac{{\kappa }_{\Vert }}{B_0^2}{F}^2+{\rho }_0{\mathcal{L}}_T]]T\hfill \\ \hfill & =-i{\rho }_{0}D(T_0)v_r-i(\gamma -1)\frac{p_0}{r}D(r{v}_r)\hfill \\ \hfill & +(\gamma -1)p_0[\frac{m}{r}{v}_{\theta }+k{v}_z]-i(\gamma -1){\rho }_0{\mathcal{L}}_{\rho }\rho \hfill \\ \hfill & -(\gamma -1)\frac{{\kappa }_{\Vert }}{B_0^2}FD(T_0)B_r,\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill B_r& =-F{v}_r,\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill B_{\theta }& =-iD(B_{0\theta }{v}_r)-k[B_{0z}{v}_{\theta }-B_{0\theta }{v}_z]+irD\left(\frac{v_{0\theta }}{r}\right)B_r,\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill B_z& =-\frac{i}{r}D(r{B}_{0z}{v}_r)+\frac{m}{r}[B_{0z}{v}_{\theta }-B_{0\theta }{v}_z]+iD(v_{0z})B_r,\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill \mathit{rY}& =[\omega -{\mathrm{\Omega }}_0]\frac{r{p}_0}{{\rho }_0}\rho +[\omega -{\mathrm{\Omega }}_0]\frac{r{p}_0}{T_0}T\hfill \\ \hfill & -ir{B}_{0\theta }D(B_{0\theta }{v}_r)-i{B}_{0z}D(r{B}_{0z}{v}_r)\hfill \\ \hfill & +r[\frac{m}{r}{B}_{0z}-k{B}_{0\theta }][B_{0z}{v}_{\theta }-B_{0\theta }{v}_z]\hfill \\ \hfill & +[r{B}_{0\theta }D\left(\frac{v_{0\theta }}{r}\right)+B_{0z}D(v_{0z})]ir{B}_r,\hfill \end{array}$$(32)

where

$$\begin{array}{c}\hfill F=\mathit{k}·{\mathit{B}}_0\text{and}{\mathrm{\Omega }}_0=\mathit{k}·{\mathit{v}}_0.\end{array}$$(33)

The former quantity is the parallel gradient operator and is an important quantity in the study of spectral stability (Goedbloed et al. 2019). At the radius where F vanishes, the Alfvén and slow continua both extend to the marginal frequency ω² = 0. The magnetic tension is minimised at that location, counteracting less the driving forces of instability. Suydam (1958) derived a local criterion for such a point to be a clusterpoint where unstable discrete modes can accumulate at. These Suydam modes are highly localised to the radial surface. The generalisation of Suydam’s criterion for a plasma with background flow was derived by Wang et al. (2004).

The latter quantity in Eq. (33) is the frequency corresponding to the Doppler shift. From these equations it can be seen that every frequency ω becomes Doppler-shifted. This shift is along the real axis in the complex plane because the wavevector and flow are both real quantities. We use the following notation for the shifted frequency

$$\begin{array}{c}\hfill \stackrel{\sim }{\omega }=\omega -{\mathrm{\Omega }}_0.\end{array}$$(34)

Shifting the frequency is not the only alteration to the modes that the background flow might be responsible for. This can be seen from the right-hand-sides of the equations. In most equations, terms which are related to the background flow are present.

We now reduce the set of Fourier-analysed MHD equations, Eqs. (24)–(31), to two coupled first order differential equations; one for the total perturbed pressure and one for the radial velocity perturbation. This is a lengthy derivation, where in the first phase the equations are substituted into each other in favor for the perturbations ρ, T, v_r, and Y. In the second part the equations are reduced to two first-order differential equations, for v_r and Y, by collecting, reshaping and eliminating terms. We finally obtain

$$\begin{array}{cc}\hfill C_{0}D(\stackrel{\sim }{\omega }Y)& =-C_{1+}(\stackrel{\sim }{\omega }Y)+C_3(r{v}_r),\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill C_{0}D(r{v}_r)& =C_2(\stackrel{\sim }{\omega }Y)+C_{1-}(r{v}_r),\hfill \end{array}$$(36)

with the coefficients given by

$$\begin{array}{cc}\hfill C_0& =r[{\rho }_0{\stackrel{\sim }{\omega }}^2-F^2]C_t,\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill C_{1\pm }& =C_1\pm C_0\frac{D({\mathrm{\Omega }}_0)}{\stackrel{\sim }{\omega }}\hfill \\ \hfill & =-\frac{2m}{r}[F{B}_{0\theta }+v_{0\theta }{\rho }_0\stackrel{\sim }{\omega }]C_t\hfill \\ \hfill & +i[\stackrel{\sim }{\omega }{\rho }_0+i(\gamma -1)[\frac{{\kappa }_{\Vert }}{B_0^2}{F}^2+{\rho }_0{\mathcal{L}}_T]]\hfill \\ \hfill & \times [(B_{0\theta }^2-{\rho }_0{v}_{0\theta }^2){\stackrel{\sim }{\omega }}^2+{[B_{0\theta }\stackrel{\sim }{\omega }+v_{0\theta }F]}^2]{\rho }_0^2{\stackrel{\sim }{\omega }}^2\pm C_0\frac{D({\mathrm{\Omega }}_0)}{\stackrel{\sim }{\omega }},\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill C_2& =-i{r}^2[\frac{m^2}{r^2}+k^2]C_t\hfill \\ \hfill & -[\stackrel{\sim }{\omega }{\rho }_0+i(\gamma -1)[\frac{{\kappa }_{\Vert }}{B_0^2}{F}^2+{\rho }_0{\mathcal{L}}_T]]r^2{\rho }_0^2{\stackrel{\sim }{\omega }}^4,\hfill \end{array}$$(39)

$$\begin{array}{cc}\hfill C_3& =i[{[{\rho }_0{\stackrel{\sim }{\omega }}^2-F^2]}^2+[{\rho }_0{\stackrel{\sim }{\omega }}^2-F^2]rD[\frac{B_{0\theta }^2}{r^2}-\frac{{\rho }_0{v}_{0\theta }^2}{r^2}]\hfill \\ \hfill & -\frac{4}{r^2}{[F{B}_{0\theta }+{\rho }_0{v}_{0\theta }\stackrel{\sim }{\omega }]}^2]C_t\hfill \\ \hfill & -\frac{1}{r^2}[\stackrel{\sim }{\omega }{\rho }_0+i(\gamma -1)[\frac{{\kappa }_{\Vert }}{B_0^2}{F}^2+{\rho }_0{\mathcal{L}}_T]]\hfill \\ \hfill & \times {[(B_{0\theta }^2-{\rho }_0{v}_{0\theta }^2){\rho }_0{\stackrel{\sim }{\omega }}^2+{\rho }_0{[B_{0\theta }\stackrel{\sim }{\omega }+v_{0\theta }F]}^2]}^2,\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill C_t& =i{\rho }_0\stackrel{\sim }{\omega }[{\rho }_0{\stackrel{\sim }{\omega }}^2(\gamma p_0+B_0^2)\hfill \\ \hfill & -\gamma p_0{F}^2]+(\gamma -1)[{\rho }_0{\stackrel{\sim }{\omega }}^2-F^2]{\rho }_0^3{\mathcal{L}}_{\rho }\hfill \\ \hfill & -(\gamma -1)[\frac{{\kappa }_{\Vert }{F}^2}{B_0^2}+{\rho }_0{\mathcal{L}}_T][{\rho }_0{\stackrel{\sim }{\omega }}^2(p_0+B_0^2)-p_0{F}^2].\hfill \end{array}$$(41)

In the limiting case of vanishing background flow, the Doppler shift and flow related terms are neglected. Equations (35) and (36) reduces to the set of equations derived by van der Linden & Goossens (1991), which in turn is formally identical to the equations obtained by Appert et al. (1974), without any non-adiabatic terms. Note the difference in notation with van der Linden & Goossens (1991).

2.1. Frieman-Rotenberg formulation

One of the things to note from the previous equations, Eqs. (37)–(40), is the ±-sign in the coefficient C_1±. This coefficient consists out of the expression C₁ and a term due to the background flow. The second term arises only because of the use of a static Eulerian perturbation, v_r.

For spectral studies of equilibria with stationary flow a more appropriate formulation exists, the Frieman-Rotenberg formulation. A full description can be found in Frieman & Rotenberg (1960) and Goedbloed et al. (2019), among others. A plasma displacement vector ξ is defined in a Lagrangian way, to connect the perturbed flow with the unperturbed flow.

In this work we use the substitution from Blokland et al. (2005)

$$\begin{array}{c}\hfill {\mathit{v}}_1=\partial \mathit{\xi }t+{\mathit{v}}_0·\mathrm{\nabla }\mathit{\xi }-\mathit{\xi }·\mathrm{\nabla }{\mathit{v}}_0.\end{array}$$(42)

For the radial component of the velocity perturbation, this yields

$$\begin{array}{c}\hfill v_r=-i\stackrel{\sim }{\omega }{\xi }_r,\end{array}$$(43)

so there is only a contribution due to the Doppler shift. This is because the background flow is only dependent on the radial coordinate.

Using this substitution, the set of two first-order differential equations are given by

$$\begin{array}{cc}\hfill C_{0}D(Y)& =-C_{1}Y-i{C}_3(r{\xi }_r),\hfill \end{array}$$(44)

$$\begin{array}{cc}\hfill C_{0}D(r{\xi }_r)& =i{C}_{2}Y+C_1(r{\xi }_r),\hfill \end{array}$$(45)

The benefit of this notation is that the cumbersome C_1± just becomes C₁, dropping the $\pm C_0\frac{D({\mathrm{\Omega }}_0)}{\stackrel{\sim }{\omega }}$.

These two first-order differential equations can be combined into one second-order differential equation. Substituting Eq. (44) into Eq. (45) yields after some manipulation

$$\begin{array}{c}\hfill D\left[\frac{C_0}{C_2}D(r{\xi }_r)\right]-[\frac{C_1^2+C_2{C}_3}{C_0{C}_2}+D\left[\frac{C_1}{C_2}\right]]r{\xi }_r=0.\end{array}$$(46)

This second-order differential equation is formally the same as obtained by Hain & Lüst (1958). The coefficients are different because they considered a static, adiabatic plasma. This equation is used in Sect. 4. Neglecting the non-adiabatic effects, Eqs. (44) and (45), are formally identical and reduce to the set of equations derived by Bondeson et al. (1987), Goossens et al. (1992) and Wang et al. (2004).

2.2. The continua

The set of first-order differential equations, Eqs. (35) and (36), becomes singular as the coefficient C₀ vanishes, as discussed by Appert et al. (1974). The eigenfunctions of such continuum solutions are non-square integrable. These are of special interest because of their local character, leading in the case of thermal instabilities to rapid, in situ condensation formation. The coefficient C₀ is given by

$$\begin{array}{c}\hfill C_0=r[{\rho }_0{\stackrel{\sim }{\omega }}^2-F^2]C_t=0.\end{array}$$(47)

The first possibility for this coefficient to vanish is when the term in the brackets vanishes. This leads to the forward- and backward-propagating Alfvén continua. The Alfvén modes are not influenced by the non-adiabatic effects. The continuum modes are real and resonate locally, meaning that the Alfvén speed varies with r and at every point you can define a continuum mode as

$$\begin{array}{c}\hfill {\stackrel{\sim }{\omega }}_{A\pm }(r)=\pm \frac{F(r)}{\sqrt{{\rho }_0(r)}}.\end{array}$$(48)

The Alfvén continuum is thus a set of frequencies where each frequency depends on a certain radius. Those modes are altered by the background flow. However, this is just a shift along the real axis and does not alter the stability of the modes.

The second way how the coefficient C₀ can vanish, is when the coefficient C_t vanishes. This is the thermal coefficient and shown here again for clarity:

$$\begin{array}{cc}& C_t(\stackrel{\sim }{\omega },r)=i{\rho }_0\stackrel{\sim }{\omega }[{\rho }_0{\stackrel{\sim }{\omega }}^2(\gamma p_0+B_0^2)-\gamma p_0{F}^2]\hfill \\ \hfill & -(\gamma -1)[\frac{{\kappa }_{\Vert }{F}^2}{B_0^2}+{\rho }_0{\mathcal{L}}_T][{\rho }_0{\stackrel{\sim }{\omega }}^2(p_0+B_0^2)-p_0{F}^2]\hfill \\ \hfill & +(\gamma -1)[{\rho }_0{\stackrel{\sim }{\omega }}^2-F^2]{\rho }_0^3{\mathcal{L}}_{\rho }=0.\hfill \end{array}$$(49)

This equation is a third-order polynomial in frequency and can be solved for every radial coordinate. The three solutions are the thermal continuum and the two, forward- and backward-propagating, slow continua. The slow continua are modified by the non-adiabatic effects and can be unstable or damped based on the physical properties of the background equilibrium. The thermal continuum was first described by van der Linden & Goossens (1991). It is the extension of the thermal mode discussed by Field (1965) for an non-uniform medium. The thermal continuum is quite useful, since it allows one to easily determine the stability of a plasma, disregarding discrete modes for the moment, without solving the set of Fourier-analysed MHD equations numerically.

When the non-adiabatic effects are neglected Eq. (49) reduces to a second-order polynomial, which holds only the slow continua. When also flow is neglected the equation than reduces to the results in for example the review paper by Goedbloed (1984).

Now what is the effect of the background flow on the coefficient C_t and by extension on the slow and thermal continua? The equation has the same functional shape as the equation obtained by van der Linden & Goossens (1991), taking into account the difference in notation. However, it is not to be solved for just the frequency. Here it needs to be solved for the frequency modified by the Doppler shift. Hence the solutions of the third-order polynomial are shifted in the complex plane. Since the Doppler shift is a shift along the real axis, it does not alter the stability. So modes that are unstable in a given equilibrium are also unstable in the analogous, flow-included equilibrium, as long as the flow does not change the equilibrium itself.

3. Numerical investigation of the continua

In the previous section we derived the equations for the continua and investigated the influence of a stationary background flow on the thermal continuum in particular. In this section we confirm our analytic results numerically using the eigenvalue code Legolas (Claes et al. 2020a; De Jonghe et al. 2022; Claes & Keppens 2023.

We start with the description of the numerical set-up which we use to represent a solar tornado in the solar corona. Secondly, we investigate the MHD spectrum of such an equilibrium. Furthermore, we vary several parameters of the equilibrium, which include the azimuthal velocity profile of the background flow and the optically thin radiative cooling curve. Lastly, we study thermal conduction and how different wavenumbers influence the stability of the thermal continuum.

3.1. Numerical set-up

The numerical code used in this work is Legolas (Claes et al. 2020a; De Jonghe et al. 2022; Claes & Keppens 2023). It is an open-source code designed to solve the set of linearised MHD equations to perform spectral analysis of the plasma. While relying on a finite element representation, the Fourier-analysed equations are transformed into a generalised eigenvalue problem of the form

$$\begin{array}{c}\hfill A\mathit{x}=\omega B\mathit{x}\end{array}$$(50)

where A and B are matrices containing equilibrium variations. The vector x denotes the state vector of perturbed quantities. Legolas is written in Fortran and modularised. It therefore has a very large range of physical effects implemented. For the work considered here, only background flow, optically thin radiative cooling, and thermal conduction are needed. However, other effects, such as resistivity and external gravity, are also available. Its post-processing Python package, pylbo, allows for easy access to the eigenvalues and eigenfunctions, and to the details of the equilibrium. Legolas supports 3D cylindrical and cartesian geometries with one dimensional variation.

The equilibrium that we use is a cylindrical plasma with one dimensional variation that represents a solar tornado. Magnetic models of solar tornadoes exist in the literature (Luna et al. 2015, 2018; Onishchenko et al. 2018). Since our main result is independent of the profiles of the equilibrium parameters, we base our model on the simpler and well-known Gold-Hoyle equilibrium, first described by Gold & Hoyle (1960). This model was also used by van der Linden & Goossens (1991), giving us the opportunity for direct comparison.

The magnetic field of a Gold-Hoyle equilibrium is given by

$$\begin{array}{cc}& B_{0\theta }=\frac{B_{0c}\mu r}{1+{\mu }^2{r}^2},\hfill \end{array}$$(51)

$$\begin{array}{cc}& B_{0z}=\frac{B_{0c}}{1+{\mu }^2{r}^2},\hfill \end{array}$$(52)

where B_0c is the total magnetic field strength at the axis. In these equations μ is the inverse pitch and is related to the magnetic field components via

$$\begin{array}{c}\hfill \mu =\frac{B_{0\theta }}{r{B}_{0z}}.\end{array}$$(53)

The magnetic field is helical, just as we want for a magnetic tornado model. The parameter μ determines how twisted the field is. We consider a constant pitch. One of the benefits of using the Gold-Hoyle magnetic field is that it is force-free.

The background flow is also helical with the profiles taken as

$$\begin{array}{cc}& v_{0\theta }=v_{0\theta b}{r}^{\alpha },\hfill \end{array}$$(54)

$$\begin{array}{cc}& v_{0z}=v_{0z},\hfill \end{array}$$(55)

where v_0θb is the azimuthal velocity at the boundary at r = 1 and α is the exponent determining the shape of the azimuthal velocity profile. The axial velocity is taken to be constant throughout the cylinder.

The pressure, temperature, and density are coupled due to the ideal gas law. If we take the density to be constant, the temperature cannot be constant, with a non-constant pressure. We assume a radially dependent density profile in the equilibrium when using Legolas

$$\begin{array}{c}\hfill {\rho }_0={\rho }_{0b}{r}^{\tau },\end{array}$$(56)

with ρ_0b the density at the boundary and τ being the exponent. However, we keep τ fixed to zero to effectively work with a constant density. Unlike the real Gold-Hoyle equilibrium the pressure is not constant. The pressure is fixed by the mechanical equilibrium of our flow-included set-up, Eq. (10). As the magnetic field is force-free, it becomes a relation between the thermal pressure and the azimuthal velocity. The thermal pressure needs to balance the centrifugal force. Using the ideal gas law and the profiles of aforementioned equilibrium quantities the profile of the background temperature is given by

$$\begin{array}{c}\hfill T_0=\frac{v_{0\theta b}^2}{2\alpha +\tau }{r}^{2\alpha }+T_{0c}{r}^{-\tau },\end{array}$$(57)

where T_0c is the background temperature at the axis when the density is constant. In the limiting case of vanishing background flow and constant density, this simplifies to the literature Gold-Hoyle equilibrium of Gold & Hoyle (1960).

We consider a domain for the dimensionless radial coordinate r from 0 to 1. The normalisation used in Legolas is described in Claes et al. (2020a). We chose to set the unit temperature, unit magnetic field, and unit length scale to 10⁶ K, 10 G, and 10⁹ cm, respectively. We consider an electron-proton plasma with a mean molecular weight of 0.5.

The default values for the physical parameters are given in Table 1. We consider a coronal plasma with a background magnetic field strength of 10 G. The plasma beta of the medium is approximately 0.1, which is typical for prominences in the solar corona (Gibson 2018). For the background flow we only consider velocities in the order of a few tens of kilometers, as observed for solar tornados by use of solar spectroscopy (Yang et al. 2018; Barczynski et al. 2021b). We initially use a r-squared profile for the azimuthal velocity. With our choice of values for the parameters, we consider a weak inhomogeneity. This ensures that the continua are well separated in the complex plane. The equilibrium values are chosen such that discrete modes corresponding to well-known instabilities, like Kelvin-Helmholtz ones which are not the topic of this investigation, are not present.

In Fig. 1 several of the background parameters and some profiles of important derived quantities are shown. In panel (A) the constant density is shown. The temperature profile is shown in panel (B). The temperature is not constant, because the pressure balances the centrifugal force. It slightly increases outwards, but remains close to 10⁶ Kelvin. The panels (C) and (D) show the azimuthal and axial magnetic field components, respectively. The magnetic field is mostly vertical, as to be expected for vertical columns. The azimuthal field is an order of magnitude weaker than the axial field. The medium is dominated by the magnetic forces as can be seen from the plasma beta shown in panel (G). The plasma beta is around 0.096 and spatially varies a little bit because of the variations in the pressure and magnetic field. The velocity profiles are shown in panels (E) and (F) of Fig. 1. The azimuthal velocity is quadratically increasing, while the axial velocity is constant. From panels (H) and (I), it can be seen that the azimuthal velocity is always subsonic and sub-Alfvénic, via the corresponding Mach numbers. The sound speed and Alfvén speed are around 165 km s⁻¹ and 587 km s⁻¹ for the given parameters, respectively. Figure 2 depicts the magnetic and velocity field lines of the stationary background in 3D. The magnetic field lines are indeed predominantly vertical, only the outermost have more twist. They are coloured according to the magnitude of the B_0θ component, highlighting the helicity. The velocity field lines are more twisted, whereas the lines at the center are mostly straight. The velocity increases going outwards, due to the increase in azimuthal flow.

Fig. 1. Equilibrium profiles of the benchmark set-up with respect to the radial coordinate. From left to right, top to bottom the quantities shown are density, temperature, azimuthal magnetic field, axial magnetic field, azimuthal velocity, axial velocity, plasma beta, Alfvén mach number, and standard mach number with respect to the sound speed.

Fig. 2. Magnetic and velocity field lines of the background equilibrium are shown left and right, respectively. The magnetic field lines are colour-coded according to the magnitude of the B0θ component. The velocity field lines take the colour of their magnitude.

For all cases, the dimensionless wavenumbers, m and k, are both set to unity, unless stated otherwise. One cannot investigate thermal instability without radiative cooling. Legolas handles this using optically thin radiative cooling tables. Commonly used cooling tables are described by Hermans & Keppens (2021) and in their appendix. In the default case we use the SPEX_DM cooling curve (Schure et al. 2009). The cooling curve is one of the parameters that we vary in the following sections. Thermal equilibrium needs to hold at the initial state. A constant background heating, equal to the radiative heating at t = 0, is set in order to facilitate this. Thermal conduction is not used in the benchmark case. However, we use it later on, but only parallel conduction.

We use Legolas 2.0, the most recent version. This second version was recently discussed in detail by Claes & Keppens (2023). Legolas has several solvers at its disposal. They have become more accurate and faster in the second version. Solvers, such as the default QR-invert, calculate the whole spectrum of eigenmodes. Other solvers, such as the Arnoldi solvers, are tailored to finding a certain amount of specific modes. We typically use the default QR-invert solver, but have used the Arnoldi solvers when in doubt of the convergence of the modes. For details about the different solvers in Legolas we refer to Claes & Keppens (2023). We used 500 gridpoints in the base grid to keep the memory usage and computation time in check. The only available boundary condition is a perfectly conduction wall at r = 1.

3.2. Benchmark

3.2.1. The adiabatic case

We first take a look at the adiabatic spectrum in order to understand some peculiarities that arise due to the background flow, regardless of non-adiabatic effects. The spectrum of the benchmark equilibrium, a Gold-Hoyle magnetic field with a constant axial and r-squared azimuthal flow profile, is given in the top panel of Fig. 3. There are insets showing the two slow continua, the backward-propagating Alfvén continuum, and the Doppler continuum. The continua are spread out along the real axis according to their intrinsic resonances and the Doppler shift. There is a single discrete slow mode next to the upper edge of the forward-propagating slow continuum. No discrete modes appear next to the backward-propagating slow continuum or the backward-propagating Alfvén continuum.

Fig. 3. Top panel: overview of the spectrum of the adiabatic benchmark equilibrium with both azimuthal and axial flow is shown. Insets zooming in on the continua are added. Middle left panel: zoom-in image of the upper edge of the forward-propagating Alfvén continuum. Three discrete Alfvén modes are marked. The corresponding eigenfunctions are shown in the middle right panel. Bottom left panel: forward-propagating Alfvén continuum as a function of radius for the adiabatic benchmark equilibrium, for the static version of the benchmark equilibrium, and for a version of the benchmark equilibrium with only axial flow. The backward-propagating Alfvén continuum is shown in the bottom right panel.

It is interesting to note that the azimuthal background flow facilitates the clustering of a series of discrete Alfvén modes to the upper edge of the forward-propagating Alfvén continuum. The discrete modes do not appear in the static case or when only including axial flow in the background equilibrium. The clustering can be seen in the middle left panel of Fig. 3. The three discrete modes are marked with coloured crosses and their eigenfunctions are shown in the middle right panel in the same colours. It is an anti-Sturmian sequence as the frequency reduces with the number of nodes. Discrete modes are known to cluster towards extrema of continua (Goedbloed & Sakanaka 1974). To explain the appearance of these discrete Alfvén modes the forward-propagating Alfvén continuum is plotted in the bottom left panel together with the Alfvén continuum of a static version of the equilibrium and that of an equilibrium that only includes axial flow. The bottom right panel shows the backward-propagating Alfvén continuum of the general flow case discussed here. From the former panel it can be seen that only in the case of both azimuthal and axial flow there is an internal maximum. In the case with only axial flow the continuum is still monotonously decreasing, but it is increased with a constant Doppler shift since the axial velocity is taken constant. The inclusion of the azimuthal flow, not linearly dependent on the radius, causes the Doppler shift to be non-constant, leading to this maximum. The backward-propagating Alfvén continuum also has no internal extremum because the addition of the Doppler shift to the continuum is not in the correct radial position to create it. Wang et al. (2004) derived a criterion for clustering towards extrema of slow and Alfvén continua using a Frobenius expansion. When a factor dependent on mode numbers and equilibrium and variations, $D_A^{+}$, evaluated at the location of the extremum is larger than 0.25, discrete modes arise. For the Alfvén continuum considered here we have $D_A^{+}(r=0.78)=1.58$. This is larger than 0.25, and hence we would indeed expect discrete Alfvén modes to cluster towards the maximum of the forward-propagating Alfvén continuum.

3.2.2. The non-adiabatic case

We now extend on the literature results by investigating the spectrum of a plasma both influenced by non-adiabatic effects and background flow. We still neglect thermal conduction for the time being.

The spectrum showing the overview of the modes is given in Fig. 4. The black dots are the results of the Legolas run. The blue, green, and red dots are the analytic expressions for the Alfvén, thermal, and slow continuum modes, respectively, and correspond to equations Eqs. (48) and (49). For the used equilibrium the thermal continuum is unstable, while the slow continua are damped. The Alfvén continua remain real and stable, as to be expected. The beginning of the fast branches are visible at the outsides of the Alfvén continua. Two modes are shown on each side. Those branches extend to infinity and are damped.

Fig. 4. Spectrum of the benchmark equilibrium. The black dots are the results obtained using Legolas. The green, red, and blue dots are the thermal, slow and Alfvén continua, respectively.

It is also clear from Fig. 4 that the results obtained with Legolas match perfectly with the analytic expressions for the continua.

In the left panel of Fig. 5 the thermal continuum demonstrates that the numerical result matches perfectly with the analytic green curve. A background flow Doppler shifts the thermal continuum as can be seen from the analytical expression, Eq. (49), derived earlier in this work. This shift is along the real axis. In the right panel of Fig. 5 the density eigenfunctions of two continuum modes, denoted with yellow and blue crosses in the left panel, are depicted. The eigenfunctions are very localised and show a non-square integrable shape, as can be seen most easily in the inset.

Fig. 5. Left panel: zoomed-in image of the thermal continuum of the benchmark equilibrium with the vertical axis as a grey dotted line. Right panel: density eigenfunctions of the orange and blue modes indicated by crosses in the left panel are shown. An inset focused on the localised non-square integrable shape is given in the right panel; vertical grey lines indicate the radial position of the resonance.

The thermal continuum has a curved shape in the complex plane. At first glance, one might think that the shape of the shifted continuum is due to the chosen velocity profiles. However this is not the complete truth. Since we are working in a cylindrical coordinate system, the Doppler shift is not just the squared shape of the azimuthal velocity profile. It is also influenced by the axial velocity profile and the wavenumbers, such that

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_0=\mathit{k}·{\mathit{v}}_0=\frac{m}{r}{v}_{\theta 0}+k{v}_{z0}.\end{array}$$(58)

In the case of a squared profile for the azimuthal velocity and a constant axial velocity, the Doppler shift is linear. This can be seen on Fig. 6, where the real and imaginary parts of the thermal continuum modes are plotted with respect to the radial coordinate. If the shape was dominated by the Doppler shift, one would expect just a straight line in the complex plane. The curved shape of the continuum is composed out of the real and imaginary part. Hence, the distribution of the imaginary part of the continuum in space is of relevance. This also has a curved shape as can be seen by the blue curve in Fig. 6. About half of the continuum modes, the ones corresponding to radii smaller than 5 Mm, have relatively low and similar growth rates. Due to the linear shift of those modes with nearly the same growth rate, an approximately horizontal part of the thermal continuum is obtained. The continuum modes, which are localised farthest outwards, are the most unstable and experience the largest Doppler shift. The fact that the outermost modes are more unstable than the innermost modes, can also be seen in Fig. 5. The equilibrium parameters, set in Sect. 3.1 and shown in Fig. 1, determine how the imaginary part of the thermal continuum varies with radius. This can change drastically with only a minor change to any of the equilibrium profiles, which is discussed further.

Fig. 6. Imaginary and real parts of the thermal continuum of the benchmark equilibrium as a function of radius (blue and orange curves, respectively).

The forward- and backward- propagating slow and Alfvén continua are shown in more detail in the top left panel of Fig. 7, denoted by their typical colours of red and blue, respectively. For both types of continua the numerical results match the analytic expressions. The slow continua are governed by Eq. (49), which is the same equation as the thermal continuum. The Alfvén continuum is spread out according to Eq. (48). Besides the continua there are some more things to note. Just as in the adiabatic case shown in Sect. 3.2.1, the inclusion of a background flow creates a maximum in the forward-propagating Alfvén continuum. Three discrete modes can be seen to cluster towards the right edge in the bottom panels of Fig. 7. The modes are marked with different colours and are very slightly damped, contrary to their completely real nature in the adiabatic case. This is due to the non-adiabatic effects. The damping of discrete Alfvén modes by non-adiabaticity has been shown by Keppens et al. (1993). The eigenfunctions are shown in the bottom right panel. The grey dotted line denotes the location of the maximum in the Alfvén continuum. It can be seen that the number of nodes increases for modes closer to the right edge of the continuum (i.e. with decreasing frequency). This is thus an anti-Sturmian sequence, as to be expected. The backward-propagating Alfvén continuum has no extremum and no discrete modes cluster towards its edges. There also appears a discrete slow mode next to the forward-propagating slow continuum. It is marked by a green cross and its eigenfunction is shown in the top right panel. The fast modes, of which one can be seen on either side in the overview panel in the top left of Fig. 7, are also damped.

Fig. 7. Slow and Alfvén continua of the benchmark equilibrium (red and blue, respectively). The results obtained with Legolas are the black dots. In the top left panel the insets show zoomed-in images of the backward-propagating Alfvén continuum and the two slow continua. The two panels below the main spectrum are zoom-ins onto the forward-progating Alfvén continuum, with another zoom-in towards the right edge. The discrete modes are marked with coloured crosses. The eigenfunction of the discrete slow mode is shown in the top right panel. In the bottom right panel the eigenfunctions of the three discrete Alfvén modes are shown together with a grey dotted line denoting the location where the Alfvén continuum reaches its maximum.

3.3. Azimuthal velocity profile

One of the important and interesting parameters of a background flow is the azimuthal velocity profile. Here we compare two different profiles with the benchmark case discussed previously. We keep the axial velocity profile constant and do not change the magnitudes.

First, we consider the basic case of an azimuthal velocity given by

$$\begin{array}{c}\hfill v_{0\theta }=v_{0\theta b}r,\end{array}$$(59)

which has a simple linear dependency on the radius, with α = 1. This linear dependence on radius means that the Doppler shift, calculated with Eq. (58), is constant. In the left panel of Fig. 8, the constant shift to the right in the complex plane of all the thermal eigenmodes can be seen. Beside the difference in shape compared to Fig. 5, the growth rate is also slightly altered. The most unstable modes have become more unstable. From Eq. (49) we have learned that the inclusion of a background flow does not alter the stability of the modes. However, including an azimuthal velocity or altering it does change other background parameters via the mechanical equilibrium that needs to be maintained. It manifests itself as a change in the pressure gradient or the magnetic field, see Eq. (10). Those quantities do influence the growth rate of the thermal continuum. Modifying the profile in an already flow-included equilibrium does only slighty change the growth rate. There is no change in growth rate when the axial velocity is altered.

Fig. 8. Thermal continuum for the α = 1 and α = 0.5 cases (left and right, respectively). The analytic expressions are shown in green, while the Legolas results are the black dots.

The second profile we consider here is a square root dependence on radius given by

$$\begin{array}{c}\hfill v_{0\theta }=v_{0\theta b}\sqrt{r}.\end{array}$$(60)

We take the parameter α as 0.5. This is a Keplerian rotation profile, in a different physical setting. The shape of the Doppler shifted thermal continuum shown in the right panel of Fig. 8 is again different. The growth rates have increased, becoming even more unstable. The most unstable modes are still the outermost, radially. Due to the different Doppler shift they are now the least shifted along the real axis.

Compared to the spectra for the previous velocity profiles, this spectrum has some kinks in it. There is one at the upper left part of the continuum and one around a growth rate of 0.00093 Hz. To see whether the Doppler shift, the imaginary part of the thermal continuum or the interplay of both is responsible, we look at the real and imaginary components of the thermal continuum. In Fig. 9 the real and imaginary parts of the continuum are shown with respect to the radial coordinate in orange and blue, respectively. We also plot the value of the derivative of the cooling rate for each radius in red. The Doppler shift is as expected from a square root azimuthal velocity profile. It is the imaginary part of the continuum that has the kinks. One can easily see that the kinks correspond to cusps in the derivative of the cooling rate with respect to the temperature. The cooling curve or its derivative do not need to be, and typically are not, smooth functions. This may even lead to discontinuous slow and thermal continua, as for example the slow continua for a solar coronal slab obtained by Claes & Keppens (2021). The different velocity profile in combination with the constraint of mechanical equilibrium translates to a slightly different temperature profile compared to the benchmark case. However, a different temperature profile also alters the cooling rate at every radial position and the derivative of the cooling rate. Hence, the change of such a small parameter may lead to different thermal and slow continua.

Fig. 9. Imaginary and real parts of the thermal continuum of the α = 0.5 case. They are shown as blue and orange curves, respectively. This figure explains how the precise shape of the continuum relates to the temperature derivative of the cooling curve, denoted in red, used in this case.

In all cases it is confirmed that the analytic expression matches the numerical results obtained using Legolas. This is to be expected since the equations are derived from the general MHD equations in a cylindrical coordinate system using a generic background flow profile.

3.4. Optically thin cooling curves

In the benchmark case we used the SPEX_DM cooling curve (Schure et al. 2009), however a wide variety of them is used in the literature. A similar cooling curve suitable for MHD simulations, being a more modern interpolated table, is the Colgan_DM curve (Colgan & Abdallah 2008). The curves are discussed in more detail in Hermans & Keppens (2021). They are plotted in Fig. 10. The corresponding thermal continuum of the equilibrium, but with the Colgan_DM cooling curve used instead, is shown in the left panel of Fig. 11. The only notable difference is the difference of factor two in the growth rate of the modes. The fact that growth and damping rates can differ due to the use of a different cooling curve has also been discussed by Soler et al. (2012) and Hermans & Keppens (2021). The outermost modes are still the most unstable.

Fig. 10. Three optically thin cooling curves used in this work. The cooling rate is plotted as a function of the temperature. The blue, orange, and green curves represent the SPEX_DM, Colgan_DM, and Rosner cooling curves, respectively.

Fig. 11. Thermal continuum for the Colgan_DM and Rosner cases (left and right, respectively). The analytic expressions are shown in green, while the Legolas results are the black dots.

A second interesting cooling curve is the Rosner curve (Rosner et al. 1978; Priest 1982). It is also shown in Fig. 10. It is an older piece-wise power law with a lower radiative cooling rate at 1 million Kelvin, leading to expected longer growth and damping rates. In the right panel of Fig. 11 the thermal continuum using the Rosner curve is shown. The growth rate is indeed significantly lower. However, what is most striking is the shape of the Doppler shifted continuum. It is mirrored with respect to the benchmark and Colgan_DM cases. This is due to the fact that the outermost modes are not the most unstable. The innermost modes are the most unstable when using the Rosner curve because of the physical parameters of the background equilibrium. The innermost modes are shifted the least, creating the flipped shape of the Doppler shifted continuum. It is important to note that the choice of cooling curve can thus alter the location of the most unstable modes. The location of the most unstable modes can be determined from the spectrum, if the velocity profiles are known.

3.5. Thermal conduction and axial wavenumber

A physical effect that we did not yet take into account in the numerical investigation of the previous subsections is thermal conduction. Its influence has been extensively studied in the literature for static equilibria and in multidimensional simulations (Field 1965; Begelman & McKee 1990; Sharma et al. 2010; Claes et al. 2020b; Jennings & Li 2021). Thermal conduction smooths out temperature gradients and small perturbations with wavelengths shorter than a characteristic length scale. The thermal and slow modes might be damped depending on physical parameters of the equilibrium and the wavenumber. The Alfvén modes are not influenced. We restrict ourselves to field-aligned, parallel conduction.

In Fig. 12 the spectrum is shown for the benchmark case, but with parallel conduction. The first thing to note is that all the continua obtained with Legolas still confirm the analytic expressions. Secondly, the thermal continuum is now damped. The modes do not lead to thermal instability in this case. However, the shape in the complex plane is still the same. The slow continuum modes are more damped and the Alfvén modes are unaltered.

Fig. 12. Spectrum of the benchmark equilibrium with parallel thermal conduction included. The black dots are the results obtained with Legolas. The green, red, and blue dots are the thermal, slow and Alfvén continua, respectively.

The previous results show that the plasma is not thermally unstable for a dimensionless axial wavenumber of unity. However, this purely coronal volume of plasma might still be unstable for different values of this wavenumber. The magnitude of the axial wavenumber can be related to the wavelength of the observed tornado structures. In the literature several values are quoted (see e.g. Wedemeyer et al. 2013 and Tziotziou et al. 2023). We take typical heights to be between 5 and 100 Mm. The associated dimensionless axial wavenumbers can be calculated as follows

$$\begin{array}{c}\hfill |k|=\frac{2\pi }{\lambda }{L}_u=\frac{2\pi \times {10}^{\mathbf{9}}\mathrm{cm}}{[5-100]\mathrm{Mm}}=\frac{2\pi \times {10}^9\mathrm{cm}}{[5-100]\times {10}^8\mathrm{cm}}\approx [0.5-12].\end{array}$$(61)

Shorter perturbations with larger axial wavenumbers might represent finestructure. In Fig. 13 the thermal continua for different values of k are shown in the complex plane. We varied k from -4 to 2, where the minus-sign denotes the direction of the wavevector. We keep the azimuthal wavenumber fixed at 1. For every continuum the Legolas results in black are accompanied by the analytic expression in colour. The continuum modes with an axial wavenumber of 0 are the most unstable. They are shifted to the right by a completely positive Doppler shift. Increasing the axial wavenumber dampens the thermal modes. The continua are also more shifted to the right. If the axial wavenumber is decreased, becoming more and more negative, the continua are more damped and shifted to the left. For small, negative axial wavenumbers the shift might not be completely positive, leading to a part of the real thermal continuum that is positive and one that is negative. Some modes of the continuum will be forward-propagating, while others will be backward-propagating.

Fig. 13. Thermal continua for the different wavenumbers given in the legend in the complex plane. The black dots correspond to a run with a given value of k, as denoted by the colour of the overlapping thermal continuum.

In Fig. 14 we show the largest growth rate of the thermal continuum with respect to the axial wavenumber. This mode triggers and determines the thermal instability, if unstable. It can be seen that the growth rate is diminished for large k. Hence small perturbations are smoothed out. Perturbations with wavenumbers less than unity can unstable for this coronal volume of plasma, as can be seen in the inset. This corresponds to perturbations with the size of roughly 60 Mm. Note that the growth rates are different for negative and positive versions with the same k. This is because k also modifies the thermal continuum through F in Eq. (49).

Fig. 14. Largest growth rate of the thermal continuum as a function of the axial wavenumber. The inset showns a to zoom-in on the unstable wavenumber range.

Parallel thermal conduction has the same effect of damping the thermal modes, based on a cut-off wavelength. Small perturbations are damped for the equilibrium parameters considered here. However, a slight change to the temperature or density might be enough to make the medium unstable for certain wavenumbers which were previously damped.

4. Discrete modes

Besides the continua, also discrete solutions to the linearised equations exist. While continuum modes have very localised singular eigenfunctions, the eigenfunctions of discrete modes vary with radius. Discrete modes exist in the form of well-known instabilities, for example in linear resistive MHD one frequently invokes the tearing mode (Furth et al. 1963). However, discrete modes can also appear as a series of modes that cluster towards the edges of continua (Goedbloed & Sakanaka 1974; Goedbloed 1984). The stability of such modes depends on the physical parameters of the equilibrium. This has been studied for Alfvén and slow modes by Goedbloed (1984). van der Linden & Goossens (1991) proved that discrete thermal modes can also exist and argued that they can alter the global, thermal stability of an equilibrium. The influence of a background flow on adiabatic discrete modes was studied by Wang et al. (2004).

In this section, we extend on the literature of discrete modes by including both a background flow and non-adiabatic effects. We look in particular at the discrete thermal and slow modes. In the first subsection we derive an analytic approximation for the frequency of the discrete modes. In the following subsections we investigate the discrete modes numerically using Legolas and visualise them in 3D.

4.1. Analytic investigation

In this section we perform a WKB analysis to derive an approximate dispersion relation near an internal extremum of the thermal or slow continuum. Discrete modes might cluster towards such a point. We follow the methodology of Goedbloed (1984), van der Linden & Goossens (1991), and Keppens et al. (1993).

The starting point is the second-order differential equation Eq. (46), which can be written in the form

$$\begin{array}{c}\hfill D\left[f(r)D(r{\xi }_r)\right]-g(r)r{\xi }_r=0\end{array}$$(62)

with the coefficients f(r) and g(r) dependent on the coefficients C₀ to C_t given by Eqs. (37)–(40). A fundamental asssumption of the WKB approximation is that the coefficients f(r) and g(r) are only weakly varying. We then assume a solution of the form

$$\begin{array}{c}\hfill r{\xi }_r=p(r)e^{i{\displaystyle \int q}(r)dr}\end{array}$$(63)

where p(r) is the amplitude modulation function and q(r) is the local radial wavenumber (Keppens et al. 1993). The exponential part is assumed to rapidly vary compared to the prefactor p(r) and with respect to the characteristic length scale of the equilibrium L: q²L² ≫ 1. The applicability of the WKB results is thus limited to modes of high radial order. The expressions for p(r) and q(r) can be determined by substituting Eq. (63) into Eq. (62) and only considering the leading order terms in the inhomogeneity. We obtain

$$\begin{array}{cc}& p\approx {(-fg)}^{-1/4},\hfill \end{array}$$(64)

$$\begin{array}{cc}& q\approx {(-g/f)}^{1/2}.\hfill \end{array}$$(65)

The latter relation defines a local dispersion relation that connects the local radial wavenumber q to the frequency. Expressing f and g in terms of the coefficients of the first-order differential equations, Eqs. (37)–(40), yields

$$\begin{array}{c}\hfill q^2=-\frac{g}{f}=-\frac{C_2}{C_0}[\frac{C_1^2+C_2{C}_3}{C_0{C}_2}+D\left[\frac{C_1}{C_2}\right]].\end{array}$$(66)

Multiplying both sides with C₀ and substituting its expression, Eq. (37), gives

$$\begin{array}{c}\hfill r[{\rho }_0{\stackrel{\sim }{\omega }}^2-F^2]C_t{q}^2=-C_2[\frac{C_1^2+C_2{C}_3}{C_0{C}_2}+D\left[\frac{C_1}{C_2}\right]].\end{array}$$(67)

We now determine approximate solutions to this local dispersion relation near the extremum of the slow and thermal continua. To that end we suppose that the continua are sufficiently far apart in the complex plane. The slow and thermal continua are strongly connected because of their coupling due to the non-adiabatic, radiative effects. Both kinds of modes are solutions of the same equation Eq. (49). We consider a Doppler-shifted thermal or slow continuum that has an internal extremum at r₀, such that ${\stackrel{\sim }{\omega }}_0=\stackrel{\sim }{\omega }(r_0)$ and $C_t(r_0,{\stackrel{\sim }{\omega }}_0)=0$. Such a continuum can be found or constructed using the equilibrium profiles.

We expand the coefficient C_t in the neighbourhood of the internal extremum

$$\begin{array}{c}\hfill C_t(r_0,{\stackrel{\sim }{\omega }}_0+ϵ)\approx C_t(r_0,{\stackrel{\sim }{\omega }}_0)+\frac{\partial C_t}{\partial \stackrel{\sim }{\omega }}ϵ=ϵ\frac{\partial C_t}{\partial \stackrel{\sim }{\omega }},\end{array}$$(68)

where ϵ is a small, complex perturbation. The left-hand side of the local dispersion relation becomes

$$\begin{array}{c}\hfill r[{\rho }_0{\stackrel{\sim }{\omega }}_0^2-F^2]\frac{\partial C_t}{\partial \stackrel{\sim }{\omega }}ϵq^2,\end{array}$$(69)

where every equilibrium quantity is evaluated at r = r₀. The left-hand side is thus accurate up to first-order of ϵ. The right-hand side of the local dispersion relation is expanded to zeroth-order of ϵ and evaluated at r = r₀. After a lengthy derivation ϵ can be determined. The approximation of the eigenfrequency is finally given by

$$\begin{array}{c}\hfill \stackrel{\sim }{\omega }\approx {\stackrel{\sim }{\omega }}_0+ϵ={\stackrel{\sim }{\omega }}_0+\frac{1}{q^2}\zeta ,\end{array}$$(70)

with ${\stackrel{\sim }{\omega }}_0$ the continuum frequency around which was expanded. The parameter ζ is given by

$$\begin{array}{cc}& \zeta =\frac{i[{\stackrel{\sim }{\omega }}_0{\rho }_0+i(\gamma -1)[\frac{{\kappa }_{\Vert }}{B_0^2}{F}^2+{\rho }_0{\mathcal{L}}_T]]}{[{\rho }_0{\stackrel{\sim }{\omega }}_0^2-F^2]\frac{\partial C_t}{\partial \stackrel{\sim }{\omega }}}\{{\rho }_0^2{\stackrel{\sim }{\omega }}_0^4[{\rho }_0{\stackrel{\sim }{\omega }}_0^2-F^2]\hfill \\ \hfill & -r{\rho }_0^2{\stackrel{\sim }{\omega }}_0^{4}D\left[\frac{B_{0\theta }^2-{\rho }_0{v}_{0\theta }^2}{r^2}\right]-\frac{4{\rho }_0{\stackrel{\sim }{\omega }}_0^2}{r^2}\frac{{[F{B}_{0\theta }+v_{0\theta }{\rho }_0{\stackrel{\sim }{\omega }}_0]}^2}{{\rho }_0{\stackrel{\sim }{\omega }}_0^2-F^2}\hfill \\ \hfill & \times [k^2{B}_0^2+{\rho }_0{\stackrel{\sim }{\omega }}_0^2-F^2]+\frac{4m{\rho }_0{\stackrel{\sim }{\omega }}_0^2}{r^3}[F{B}_{0\theta }+v_{0\theta }{\rho }_0{\stackrel{\sim }{\omega }}_0]\hfill \\ \hfill & \times [2B_{0\theta }^2-{\rho }_0{v}_{0\theta }^2]-\frac{1}{r^2}(\frac{m^2}{r^2}+k^2)[{[2B_{0\theta }^2-{\rho }_0{v}_{0\theta }^2]}^2{\rho }_0{\stackrel{\sim }{\omega }}_0^2\hfill \\ \hfill & -{\rho }_0^2{v}_{0\theta }^2{[2B_{0\theta }{\stackrel{\sim }{\omega }}_0+v_{0\theta }F]}^2]-r{\rho }_0^2{\stackrel{\sim }{\omega }}_0^2[D\left[\frac{v_{0\theta }^2{F}^2}{r^2}\right]\hfill \\ \hfill & +2{\stackrel{\sim }{\omega }}_{0}D\left[\frac{v_{0\theta }{B}_{0\theta }F}{r^2}\right]+{\stackrel{\sim }{\omega }}_0^{2}D\left[\frac{{\rho }_0{v}_{0\theta }^2}{r^2}\right]]\}\hfill \end{array}$$(71)

with

$$\begin{array}{cc}& \frac{\partial C_t}{\partial \stackrel{\sim }{\omega }}=2i{\rho }_0^2{\stackrel{\sim }{\omega }}_0^2[\gamma p_0+B_0^2]-\frac{p_0(\gamma -1)}{{\stackrel{\sim }{\omega }}_0}[\frac{{\kappa }_{\Vert }}{B_0^2}{F}^2+{\rho }_0{\mathcal{L}}_T-\frac{{\rho }_0^2{\mathcal{L}}_{\rho }}{T_0}]\hfill \\ \hfill & \times [{\rho }_0{\stackrel{\sim }{\omega }}_0^2+F^2]-{\rho }_0{\stackrel{\sim }{\omega }}_0(\gamma -1)B_0^2[\frac{{\kappa }_{\Vert }}{B_0^2}{F}^2+{\rho }_0{\mathcal{L}}_T].\hfill \end{array}$$(72)

All parameters are evaluated at the local extremum of the continuum, (r₀, ${\stackrel{\sim }{\omega }}_0$). Furthermore, we used the fact that $D(C_t)(r_0,{\stackrel{\sim }{\omega }}_0)=0$ to simplify the derivation.

When background flow is neglected, our approximated solution reduces to Eqs. (4.10) and (4.11) of van der Linden & Goossens (1991), taking into account the difference in notation. Moreover, in the case of an ideal non-adiabatic plasma, the expressions reduces to those for discrete slow modes found by Goedbloed (1984). A Frobenius expansion around the extremum would give similar results. Wang et al. (2004) derived a local clustering criterium for discrete slow modes using this method. Our expression reduces to their result. Hence, the results obtained in this chapter reduce to all simpler cases.

Equation (70) is not a closed instability criterion. Nonetheless, it can be used to predict the appearance of discrete modes numerically without solving the differential equations, as we show in Sect. 4.2. The parameter ζ, given by Eq. (71), contains many free parameters and profiles. This makes the expression quite unclear and hard to handle. With respect to the background flow, it is important to note that the axial component does not play a role. The azimuthal component is nevertheless omnipresent in the expression.

The sign of ζ is crucial in the determination of the thermal stability of an equilibrium. An equilibrium is susceptible to thermal instability when thermal or slow modes are unstable; in other words, when they have a positive growth rate. The most unstable mode determines the evolution of the system. Discrete modes can alter the stability when they are more unstable. Discrete modes can exist when they lie above a maximum or below a minimum of the continuum. Otherwise, they would be absorbed into it. In the case of a maximum in the continuum ζ needs to be positive, vice versa for a minimum. An interesting case would be an equilibrium where all the continua are damped, but with unstable discrete thermal modes. The condensations of equilibria with the most unstable modes being discrete modes rather than continuum modes also have a different shape, a spread-out profile instead of a very localised perturbation.

Note that a similar expression can be derived for the non-adiabatic Doppler shifted Alfvén continuum. However, we did not pursue it here because it deviates too far from our main topic of thermal modes. Discrete Alfvén modes were studied by Goedbloed (1984) in an ideal plasma. van der Linden & Goossens (1991) and Keppens et al. (1993) included radiative cooling and thermal conduction. The adiabatic case with a background flow was studied by Wang et al. (2004).

4.2. Numerical investigation

The spectra shown in Sect. 3 do not show any discrete thermal modes. This is because of the equilibrium parameters used. The thermal continuum does not have extrema for modes to cluster to. In this section we look at a few equilibria that do have an extremum in their thermal continuum. We verify our expression to predict the existence of discrete thermal modes derived in the previous subsection.

4.2.1. Set-up

We base the three equilibria on the benchmark set-up discussed in Sect. 3.1. We keep most parameters the same. The wavenumber k is for all three cases set to 0.5. The two parameters that are changed are the temperature at the axis and the magnetic inverse pitch, μ. Varying the inverse pitch affects the shape of the magnetic field components. In Fig. 15 the profiles for the background temperature, azimuthal and axial magnetic field components, and the plasma beta are shown. Each row corresponds to a case. The temperature of case 1 and 2 are the same and around 357 000 K, which is about a third of the benchmark case. The optically thin cooling rate is different, as is its derivative. This influences the thermal continuum a lot and can cause extrema. The difference between the first and second case is the inverse pitch, this is set to 5 for the former and to 60 for the latter. The field is predominantly vertical for low inverse pitch. The plasma beta also remains below unity for the whole domain, as can be seen in the right-most panel of the top row of Fig. 15. For case two the magnetic field drops off quickly (i.e. it is more centralised around the axis). In the third case the inverse pitch is set to 50 and the temperature is again lowered. The derivative of the cooling rate with respect to temperature is positive for low temperatures, around 84 000 K for the third case. We therefore expect a damped thermal continuum.

Fig. 15. Background temperature, azimuthal magnetic field, axial magnetic field, and plasma beta profiles of the three equilibria with discrete thermal modes. Each row corresponds to a different case.

4.2.2. Results

The results of the Legolas runs for the three equilibria described previously are shown in Fig. 16. The left panels contain the thermal continua, with the discrete thermal modes marked, in the complex plane for the three cases. For the second case only a part of the continuum is plotted. This is to focus on the discrete modes. The right panels show the density eigenfunctions of the marked discrete modes. A grey dotted line denotes the radial location of the extremum. Each row of Fig. 16 corresponds to one of the three equilibria.

Fig. 16. Thermal continua, discrete thermal modes, and corresponding density eigenfunctions for the three cases. Each row corresponds to a different case. The Legolas results are shown as black dots. The analytic thermal continuum is overlaid in green. The discrete thermal modes are marked by coloured crosses. The colours of the eigenfunctions match the colours of the discrete modes. The grey dotted lines in the right panels denote the radial location of the extremum of the thermal continuum.

First of all, in contrast to the thermal continua discussed in Sect. 3, there is an extremum. In the first and second case it is a maximum. Hence, a different temperature and by extension cooling rate and its derivatives, has a large impact on the continuum. For the case for which the plasma beta is still completely below unity, case 1, there is one discrete mode. The eigenmode for the density is positive as to be expected for the growth of a condensation. For the second equilibrium, shown in the middle row of Fig. 16, there are four discrete modes. The number of nodes increase as the frequency decreases in an anti-Sturmian fashion. The third equilibrium has a much lower temperature. The thermal continuum is completely damped due to the derivative of the cooling rate being positive. The continuum has a minimum and there are four discrete modes that cluster towards that minimum. This is a Sturmian sequence.

For all three cases we can check whether our expression found by the WKB-approximation gives the expected sign. As there is still the unknown parameter q in the expression of ϵ, in Eq. (70), only the sign of ζ can be used to predict the existence of discrete modes. The magnitude of ζ is not important because q is unknown and assumed to be large, ensuring that ϵ is small. When there is a maximum in the thermal continuum for discrete modes to cluster to the imaginary part of ζ has to be positive, and vice versa for a minimum. The discrete modes are on the right side of the extremum for all three cases and hence the real part is expected to be positive. The calculated values are 0.022 + 0.09i, 0.003 + 0.66i, 0.448 − 10.42i for the first, second, and third case, respectively. The signs are all in perfect agreement with expectations. It should be noted that the precision in the numerical values of these numbers is very sensitive to round-off since quantities that are of very similar size are subtracted. Small differences due to normalisation and physical constants used also affects them.

4.3. Visualisation of eigenmodes

By solving the eigenvalue problem for a given equilibrium, you obtain the frequencies and eigenfunctions of the eigenmodes of the system. Claes (2022) describes how Legolas solves the eigenvalue problem and calculates the eigenfunctions. For this to work, we assumed a wave mode representation of the perturbations. In Eq. (21), this was given in cylindrical coordinates by

$$\begin{array}{c}\hfill f_1(\mathit{r},t)=f^{\prime }(r)e^{i(m\theta +kz-\omega t)},\end{array}$$(73)

where the primed f′(r) are the amplitudes of the waves, i.e. the eigenfunctions. A new feature in Legolas 2.0 uses this equation to visualise the modes in multiple dimensions, as presented in Claes & Keppens (2023). Setting values for t, θ, and/or z allows one to calculate the perturbation at that given location using the eigenfunction, wavenumbers, and frequency. Depending on for which parameter a range is supplied, the temporal evolution, spatial slices or full 3D datacubes can be visualised. It is important to note that the views obtained using this method represent the linear evolution of the waves and not the non-linear behaviour. Nevertheless, it is a useful and simple way to gain information without performing full non-linear multidimensional MHD simulation.

We here look at the discrete modes and use the eigenfunctions and frequencies determined in the previous subsection. All the quantities plotted in the figures of the subsection are dimensionless, while the axis are rescaled to represent the actual length scales and timescales. We use a colour scale where darker tones represent higher densities.

Let us first look at the temporal evolution of the least interesting case, the third equilibrium with damped discrete modes. The temporal evolution can be visualised by choosing a fixed θ and z to look at and vary t. We took the default value of 0 for θ and z. We use a superposition of the four discrete modes with an amplitude ratio of unity for the four modes. The evolution is shown in the bottom panel of Fig. 17. The initial perturbation based on the eigenfunctions is excited at t = 0 and also depicted in the top panel. The eigenfunctions are coloured using the same colourscheme as in Fig. 16. Since the growth rates are negative the modes are damped. This is clearly visible. Locations where the amplitude is bigger take more time to diminish.

Fig. 17. Temporal view for discrete modes of the third case, the damped eigenmodes. The eigenfunctions are shown in the top panel and coloured according to Fig. 16.

The discrete modes of the second equilibrium, set up in Sect. 4.2.1, are of more physical relevance. The four discrete modes have positive growth rates and are thus unstable. The temporal evolution is shown in Fig. 18. Indeed, the modes become amplified. At later times it can be seen how the density is spread along the radius. The superposition of the eigenfunctions is negative between 20 Mm and 40 Mm, as can be seen in the top panel. This leads to a decrease in density in this region. In the outer part the density is increased. The discrete modes have a more global behaviour, in contrast to the ultra-localised eigenfunctions of the continuum modes. We did not favour any of the discrete modes over the others; in other words, none of the eigenfunctions was scaled differently.

Fig. 18. Temporal view for discrete modes of the second case, the unstable eigenmodes. The eigenfunctions are shown in the top panel and coloured according to Fig. 16.

In Fig. 19 we present the 2D slice at constant height z = 0 for the four discrete modes of the second equilibrium in the different panels for each mode. The top left panel contains the angular distribution of the most unstable mode, the mode with the largest growth rate and with the least nodes because of the Sturmian character of the sequence. The other modes are visualised in the other three panels, from left to right and top to bottom, the modes are ordered following the Sturmian sequence. We look at t = 0, to negate the effect of the temporal evolution and to focus on the angular distribution.

Fig. 19. Spatial cuts through the cylinder at z = 0 and t = 0, for the four discrete eigenmodes of the second equilibrium. From left to right and top to bottom each panel shows a different mode starting from the mode with the least number of nodes to the one with the most nodes in the eigenfunction.

There is a clear difference between the left and right side of the slices. In the case of the most unstable mode in the top left panel, a condensation forms in the right while the left is vacated. This is due to the azimuthal wavenumber of the mode. The eigenmodes that we considered all have m = 1 and thus are non-axisymmetric, like a ‘kink’ mode. From a mathematical point of view this can be understood by looking at the angular part of Eq. (21), e^imθ.

For m = 0 modes the exponential is constant and independent of the angle. However, for other values of m this is not the case. For m = 1 modes the exponential varies between −1 and 1, taking those values at θ = π and θ = 0 respectively. The magnitude of the eigenfunction is thus opposite on the other side of the axis. When going through the panels, it can also be seen that the modes become increasingly localised as they approach the continuum.

The last 2D slice that we present here is one at constant angle, θ = 0, and at the initial time. Figure 20 shows the variation of the density eigenfunction of the superposition of the four discrete modes of the second equilibrium with height for each radial position. The height is here varied from 0 to 8 Mm in the coordinate z. This corresponds the approximate height for prominence tornado structures. Due to the shape of the eigenfunctions the most dense region is between 40 and 80 Mm in height and 2 and 4 Mm in width. The condensations formed by these eigenmodes are thus more patchy in space.

Fig. 20. Spatial cut at a fixed angle, the z-r plane, for the superposition of the discrete modes of case 2 at time t = 0.

The data cube can also be visualised in 3D. We exported the equilibrium and eigenfunctions for the four unstable discrete modes using the functionality of Legolas. The obtained files were visualised using ParaView¹. In Fig. 21 we show the evolution of the density perturbation in the cylinder of height 80 Mm. Cuts are made for clarity. The left cylinder is at the initial time, the second at t = 3, and the third at t = 5. The density and noted timestamps are dimensionless. The density perturbation grows over time based on the initial shape of the eigenfunctions. The 2D views can be observed very clearly in the 3D images, by considering the cuts that have been made. The half of the bottom circle depicts the variation by the angle, as in Fig. 19. The variation with height of Fig. 20 is clearly seen in the right half of the open cylinder. The large imaginary part of the frequency of the superposition of the discrete modes dominates the temporal behaviour of the modes, whereas there appears little to no influence of the much smaller real part of the frequency.

Fig. 21. 3D visualisations of the density perturbation throughout the cylinder at t = 0, t = 3, and t = 5.

The streamlines of the total velocity field, background field plus perturbation, is shown in 3D in Fig. 22. The field lines are coloured according to the local v_r component of the velocity perturbation. There are three slices with v_r at different heights provided. The top two are made a little bit more transparent in order to see the field lines more clearly. The total velocity field is mostly helical. However, due to the perturbation it is more complex than the equilibrium background flow, as shown in Fig. 2. The added radial component distorts it. The radial component varies with height, as can be seen from the slices. The field lines can be seen moving inwards and outwards depending on their location. It should be noted that we did not rescale the eigenfunctions. Since they are of roughly the same order of magnitude as the background quantities, their influence is expected to be large.

Fig. 22. Field lines of the total flow field. The field lines and slices are coloured according to the vr component of the velocity perturbation.

5. Discussion and conclusions

Flows are of the upmost importance in understanding the evolution of plasma in dynamical environments. Due to the increased complexity that they introduce, they are typically neglected when studying waves via MHD spectroscopy. Solar tornadoes are a class of prominences that appear to be rotating and could be formed by thermal instability. We therefore extended the work of van der Linden & Goossens (1991) by including a helical background flow.

We went through the linearisation and Fourier-analysis of the non-adiabatic MHD equations in a cylindrical coordinate system with a generic helical flow. The set of equations is transformed into two first-order differential equations, one for v_r and one for Y, in the spirit of the work by Appert et al. (1974). The equations differ from those derived by van der Linden & Goossens (1991). Most importantly, the frequency is Doppler shifted. Additional terms appear they are all dependent on the azimuthal background velocity v_0θ. There are no terms that depend on the axial background velocity v_0z. This means that azimuthal flow is of the utmost importance, which was also noticed by Wang et al. (2004) for adiabatic equilibria with rotational flow. Because of this, clustering conditions derived from these equations, such as clustering of slow, Alfvén, and Suydam modes (Wang et al. 2004) and clustering of thermal modes considered in this work, are altered by azimuthal flow. The expressions for the continua were studied in more detail. They are only altered by flow due to the additional Doppler shift to the frequency.

The analytic expressions were verified using the eigenvalue code Legolas. The force-free helical magnetic field based on Gold & Hoyle (1960), helical flow field and coronal parameters mimic the environment to form a solar tornado. The true nature of solar tornadoes is a contested topic (Gunár et al. 2023). It should be noted that we are interested in thermal instability with background flow and only use the tornadoes as examples of condensations formed in the solar corona. The magnetic field of the benchmark case was predominantly vertical, with a smaller azimuthal component. It was set up this way because structures in the solar corona, such as coronal loops and prominences, have been observed or set up in numerical simulations vertically. All analytic expressions match the numerical results. For the benchmark case the thermal continuum is unstable. Slightly damped discrete Alfvén modes accumulate to a maximum in the Doppler shifted Alfvén continuum solely arising due to the background flow. The choice of cooling curve can cause the thermal continuum to take a different shape in the spectral plane. The location of the most unstable modes can vary from the inside to the outside. This has significant implications. The most unstable modes govern the behaviour of the plasma as a whole. Since thermal continuum modes are hyper-localised, the location can thus be altered just by chosing a different cooling curve. Thermal conduction damps out small perturbations. This has been found in the literature (Field 1965) and has been shown in this work.

We derived an approximate expression to check whether discrete thermal and slow modes cluster towards an internal extremum of their respective continuum in the manner of van der Linden & Goossens (1991). Equation (70) is not a closed criterion due to the unknown local wavenumber q that it contains. However, using the sign of ζ, given by Eq. (71) and obtained through the WKB analysis, it can still be used to predict the appearance of the discrete modes. Using Legolas we verified the expression for three equilibria. Two of them have discrete thermal modes above the continuum. For these the formation process of the condensations is governed by the eigenfunctions of the discrete modes. We therefore visualised the eigenfunctions in 2D spatial and temporal slices and in 3D. The distribution of the condensation due to the density perturbation has been shown. It is not a uniform condensation or singular condensation, due to the complicated shape of the eigenfunction and the dependence of the assumed waveform on the spatial coordinates and wavenumbers. Three-dimensional visualisation of the total velocity field shows that the field retains its helical shape but is also heavily influenced by the radial velocity perturbation. The choice of m = 1, is clearly visible in the 2D and 3D plots. It would be interesting to investigate whether there are discrete thermal modes for different wavenumbers and how their eigenfunctions behave in space and time. Skirvin et al. (2023) investigated the behaviour of adiabatic m = 0 and m = ±1 modes in a flux tube with rotational flow and showed that the m = 0 modes are not influenced that much by the flow. The flux tube that they considered is embedded in an external photosphere. Just as in the well-known work by Edwin & Roberts (1983), this kind of set-up leads to surface and body modes, depending on where they propagate. However, we did not consider an external medium as we mainly focused on the thermal modes, which are internal due to their local nature.

Multidimensional MHD simulations are an interesting way to further study the formation of condensations in helical configurations with flow. Due to the 1D radial variation of the equilibria used in MHD spectroscopy, the physical effects depending on other coordinates are neglected. In the context of prominences and solar tornadoes, gravity and the support against it are of the utmost importance. In MHD simulations this can be included, and parametric studies can be performed to see what kinds of flows and magnetic field configurations are needed to support the prominence mass in helical field, futher extending the analytic work by Luna et al. (2015). The formation and appearance of condensations might be more transient. Comparing synthetic images with observations can teach us a lot about the existence of solar tornadoes. A second important physical improvement is the inclusion of a realistic solar atmosphere with density and temperature stratification, as recently used in Jenkins & Keppens (2021) and Jerčić & Keppens (2023), among others. Different prescriptions for background heating can be used. Models varying by height alone or as a function of local density and magnetic field are commonly used. They have been shown to affect the location and morphology of prominence formation (Brughmans et al. 2022). The interplay between fragmentation of the condensations in the far non-linear stage, as shown by Hermans & Keppens (2021), and flow is definitely worth investigating.

It is important to note that the analytic equations derived in this work are valid for more physical environments than solely the solar corona. The solar coronal set-up, which mimicks possible tornado formation, was used as an example. Thermal instability is a general phenomenon based on energy loss by radiation, and flows are present on every scale in the universe. Therefore, there are many more interesting applications, such as the clumpiness of galactic ouflows (Veilleux et al. 2020) and the three-phase nature of the dynamic interstellar medium (Cox 2005).

Acknowledgments

We wish to thank the anonymous referee, for the constructive comments that improved the paper, and the editorial office. J.H. would like to thank N. Claes and J. De Jonghe for the Legolas support and V. Jerčić and N. Brughmans for the useful discussions. The visualisations were achieved using the open source software Python (http://python.org) and ParaView (http://paraview.org). Both authors are supported by the ERC Advanced Grant PROMINENT and a joint FWO-NSFC grant G0E9619N. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 833251 PROMINENT ERC-ADG 2018). This research is further supported by Internal funds KU Leuven, project C14/19/089 TRACESpace.

References

All Tables

All Figures

	Fig. 1. Equilibrium profiles of the benchmark set-up with respect to the radial coordinate. From left to right, top to bottom the quantities shown are density, temperature, azimuthal magnetic field, axial magnetic field, azimuthal velocity, axial velocity, plasma beta, Alfvén mach number, and standard mach number with respect to the sound speed.
In the text

	Fig. 2. Magnetic and velocity field lines of the background equilibrium are shown left and right, respectively. The magnetic field lines are colour-coded according to the magnitude of the B0θ component. The velocity field lines take the colour of their magnitude.
In the text

Fig. 3. Top panel: overview of the spectrum of the adiabatic benchmark equilibrium with both azimuthal and axial flow is shown. Insets zooming in on the continua are added. Middle left panel: zoom-in image of the upper edge of the forward-propagating Alfvén continuum. Three discrete Alfvén modes are marked. The corresponding eigenfunctions are shown in the middle right panel. Bottom left panel: forward-propagating Alfvén continuum as a function of radius for the adiabatic benchmark equilibrium, for the static version of the benchmark equilibrium, and for a version of the benchmark equilibrium with only axial flow. The backward-propagating Alfvén continuum is shown in the bottom right panel.

In the text

	Fig. 4. Spectrum of the benchmark equilibrium. The black dots are the results obtained using Legolas. The green, red, and blue dots are the thermal, slow and Alfvén continua, respectively.
In the text

Fig. 5. Left panel: zoomed-in image of the thermal continuum of the benchmark equilibrium with the vertical axis as a grey dotted line. Right panel: density eigenfunctions of the orange and blue modes indicated by crosses in the left panel are shown. An inset focused on the localised non-square integrable shape is given in the right panel; vertical grey lines indicate the radial position of the resonance.

In the text

	Fig. 6. Imaginary and real parts of the thermal continuum of the benchmark equilibrium as a function of radius (blue and orange curves, respectively).
In the text

Fig. 7. Slow and Alfvén continua of the benchmark equilibrium (red and blue, respectively). The results obtained with Legolas are the black dots. In the top left panel the insets show zoomed-in images of the backward-propagating Alfvén continuum and the two slow continua. The two panels below the main spectrum are zoom-ins onto the forward-progating Alfvén continuum, with another zoom-in towards the right edge. The discrete modes are marked with coloured crosses. The eigenfunction of the discrete slow mode is shown in the top right panel. In the bottom right panel the eigenfunctions of the three discrete Alfvén modes are shown together with a grey dotted line denoting the location where the Alfvén continuum reaches its maximum.

In the text

	Fig. 8. Thermal continuum for the α = 1 and α = 0.5 cases (left and right, respectively). The analytic expressions are shown in green, while the Legolas results are the black dots.
In the text

	Fig. 9. Imaginary and real parts of the thermal continuum of the α = 0.5 case. They are shown as blue and orange curves, respectively. This figure explains how the precise shape of the continuum relates to the temperature derivative of the cooling curve, denoted in red, used in this case.
In the text

	Fig. 10. Three optically thin cooling curves used in this work. The cooling rate is plotted as a function of the temperature. The blue, orange, and green curves represent the SPEX_DM, Colgan_DM, and Rosner cooling curves, respectively.
In the text

	Fig. 11. Thermal continuum for the Colgan_DM and Rosner cases (left and right, respectively). The analytic expressions are shown in green, while the Legolas results are the black dots.
In the text

	Fig. 12. Spectrum of the benchmark equilibrium with parallel thermal conduction included. The black dots are the results obtained with Legolas. The green, red, and blue dots are the thermal, slow and Alfvén continua, respectively.
In the text

	Fig. 13. Thermal continua for the different wavenumbers given in the legend in the complex plane. The black dots correspond to a run with a given value of k, as denoted by the colour of the overlapping thermal continuum.
In the text

	Fig. 14. Largest growth rate of the thermal continuum as a function of the axial wavenumber. The inset showns a to zoom-in on the unstable wavenumber range.
In the text

	Fig. 15. Background temperature, azimuthal magnetic field, axial magnetic field, and plasma beta profiles of the three equilibria with discrete thermal modes. Each row corresponds to a different case.
In the text

Fig. 16. Thermal continua, discrete thermal modes, and corresponding density eigenfunctions for the three cases. Each row corresponds to a different case. The Legolas results are shown as black dots. The analytic thermal continuum is overlaid in green. The discrete thermal modes are marked by coloured crosses. The colours of the eigenfunctions match the colours of the discrete modes. The grey dotted lines in the right panels denote the radial location of the extremum of the thermal continuum.

In the text

	Fig. 17. Temporal view for discrete modes of the third case, the damped eigenmodes. The eigenfunctions are shown in the top panel and coloured according to Fig. 16.
In the text

	Fig. 18. Temporal view for discrete modes of the second case, the unstable eigenmodes. The eigenfunctions are shown in the top panel and coloured according to Fig. 16.
In the text

	Fig. 19. Spatial cuts through the cylinder at z = 0 and t = 0, for the four discrete eigenmodes of the second equilibrium. From left to right and top to bottom each panel shows a different mode starting from the mode with the least number of nodes to the one with the most nodes in the eigenfunction.
In the text

	Fig. 20. Spatial cut at a fixed angle, the z-r plane, for the superposition of the discrete modes of case 2 at time t = 0.
In the text

	Fig. 21. 3D visualisations of the density perturbation throughout the cylinder at t = 0, t = 3, and t = 5.
In the text

	Fig. 22. Field lines of the total flow field. The field lines and slices are coloured according to the vr component of the velocity perturbation.
In the text