Issue 
A&A
Volume 686, June 2024



Article Number  A180  
Number of page(s)  18  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202348337  
Published online  11 June 2024 
A spectroscopic investigation of thermal instability for cylindrical equilibria with background flow
Centre for Mathematical PlasmaAstrophysics, Celestijnenlaan 200B, 3001 Leuven, KU Leuven, Belgium
email: jorishermans96@gmail.com
Received:
20
October
2023
Accepted:
27
February
2024
Context. Flows are omnipresent and govern the dynamics of plasma. Solar tornadoes are a class of apparently rotating prominences that might be formed by thermal instability. In spectroscopic studies on thermal instability, background flow is commonly neglected.
Aims. We here determine the effect of background flow on thermal instability in cylindrical magnetic field configurations. How various parameters affect the distribution of eigenmodes in the magnetohydrodynamic (MHD) spectrum is also explored. We investigate whether discrete thermal modes exist.
Methods. In an analytical study, we extended upon the literature by including a generic background flow in a cylindrical coordinate system. The nonadiabatic MHD equations are linearised, Fourieranalysed, and examined to understand how a background flow changes the continua. An approximate expression for discrete thermal modes is derived using a WentzelKramersBrillouin (WKB) analysis. The analytical results are then verified for a benchmark equilibrium using the eigenvalue code Legolas. The eigenfunctions of discrete thermal modes are visualised in 2D and 3D.
Results. The thermal continuum is Dopplershifted due to the background flow, just like the slow and Alfvén continua. Discrete modes are altered because the governing equations contain flowrelated terms. An approximate expression to predict the appearance of discrete thermal modes based on the equilibrium parameters is derived. All analytical expressions match the numerical results. The distribution of the density perturbations of the discrete thermal modes is not a uniform or singular condensation, due to the shape of the eigenfunctions and the dependence of the assumed waveform on the coordinates and wavenumbers. A 3D visualisation of the total velocity field shows that the helical field is heavily influenced by the radial velocity perturbation.
Conclusions. We derived analytic expressions for nonadiabatic MHD modes of a cylindrical equilibrium with background flow and verified them using a coronal equilibrium. However, the equations are valid for and can be applied in other astrophysical environments.
Key words: instabilities / magnetohydrodynamics (MHD) / methods: analytical / Sun: filaments, prominences
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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 quasiperiodic 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). Lineofsight 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 waveform, 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 nonadiabatic terms are included. They also derived approximate expressions to predict the occurance of discrete thermal modes using a WentzelKramersBrillouin (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 opensource code used to solve the eigenvalue problem of the linearised and Fourieranalysed 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 largescale 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 nonadiabatic MHD equations are linearised, Fourieranalysed, and transformed into two coupled firstorder 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 forcefree GoldHoyle 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 columnlike 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 quasicontinuum modes that are in essence discrete. The eigenfunctions of those quasicontinuum 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 nonadiabatic effects. The following set of dimensionless MHD equations are appropriate to describe this medium:
with the magnetic field in units, where μ_{0} = 1. The quantities ρ, T, p, v, and B are the density, temperature, pressure, velocity, and magnetic field, respectively. The function ℒ is the net heatloss function and is defined as
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 righthand 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
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
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
The background magnetic field is given by
The divergencefree condition of the magnetic field, ∇ ⋅ B = 0, is automatically satisfied for our background magnetic field because the nonradial 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.
We denote the derivative in the direction of r (i.e. ) 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 righthandside. 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
where ℒ_{0} is the heatloss function at the initial state. Intuitively this condition means that the energyloss by the optically thin radiative cooling needs to be balanced by the background heating at t = 0. As an equation this means
We consider the Eulerian perturbations to be small, so we can expand around the equilibrium in the form
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 higherorder terms. Substituting these expansions into the MHD equations, Eqs. (1)–(3), yields
We used the equilibrium conditions of mechanical and thermal equilibrium, together with the fact that the background is stationary, .
The quantities ℒ_{ρ} and ℒ_{T} are the partial derivatives of the heatloss function with respect to density and temperature, respectively. Hence, they are given by
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
The first term on the righthandside vanishes because of the divergencefree background magnetic field. We also used the fact B_{r} = 0, which implies that
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 Fourieranalyse 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 planewaves
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
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.
Since equilibrium values are always denoted with a subscript ‘0’, we can drop the prime on the perturbed quantities. The Fourieranalysed MHD equations are given by
where
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 ω^{2} = 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 Dopplershifted. 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
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 righthandsides of the equations. In most equations, terms which are related to the background flow are present.
We now reduce the set of Fourieranalysed 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 firstorder differential equations, for v_{r} and Y, by collecting, reshaping and eliminating terms. We finally obtain
with the coefficients given by
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 nonadiabatic terms. Note the difference in notation with van der Linden & Goossens (1991).
2.1. FriemanRotenberg 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_{1} 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 FriemanRotenberg 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)
For the radial component of the velocity perturbation, this yields
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 firstorder differential equations are given by
The benefit of this notation is that the cumbersome C_{1±} just becomes C_{1}, dropping the .
These two firstorder differential equations can be combined into one secondorder differential equation. Substituting Eq. (44) into Eq. (45) yields after some manipulation
This secondorder 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 nonadiabatic 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 firstorder differential equations, Eqs. (35) and (36), becomes singular as the coefficient C_{0} vanishes, as discussed by Appert et al. (1974). The eigenfunctions of such continuum solutions are nonsquare 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_{0} is given by
The first possibility for this coefficient to vanish is when the term in the brackets vanishes. This leads to the forward and backwardpropagating Alfvén continua. The Alfvén modes are not influenced by the nonadiabatic 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
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_{0} can vanish, is when the coefficient C_{t} vanishes. This is the thermal coefficient and shown here again for clarity:
This equation is a thirdorder polynomial in frequency and can be solved for every radial coordinate. The three solutions are the thermal continuum and the two, forward and backwardpropagating, slow continua. The slow continua are modified by the nonadiabatic 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 nonuniform 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 Fourieranalysed MHD equations numerically.
When the nonadiabatic effects are neglected Eq. (49) reduces to a secondorder 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 thirdorder 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, flowincluded 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 setup 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 setup
The numerical code used in this work is Legolas (Claes et al. 2020a; De Jonghe et al. 2022; Claes & Keppens 2023). It is an opensource 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 Fourieranalysed equations are transformed into a generalised eigenvalue problem of the form
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 postprocessing 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 wellknown GoldHoyle 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 GoldHoyle equilibrium is given by
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
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 GoldHoyle magnetic field is that it is forcefree.
The background flow is also helical with the profiles taken as
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 nonconstant pressure. We assume a radially dependent density profile in the equilibrium when using Legolas
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 GoldHoyle equilibrium the pressure is not constant. The pressure is fixed by the mechanical equilibrium of our flowincluded setup, Eq. (10). As the magnetic field is forcefree, 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
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 GoldHoyle 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^{6} K, 10 G, and 10^{9} cm, respectively. We consider an electronproton 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 rsquared 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 wellknown instabilities, like KelvinHelmholtz ones which are not the topic of this investigation, are not present.
Values for the physical quantities used in the benchmark case. B_{0c} is the magnetic field strength at the axis; μ is the inverse pitch of the magnetic field; ρ_{0b}, τ, and T_{0c} denote the density at the boundary, the exponent of the density profile, and the temperature at the axis, respectively.
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^{6} 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 subAlfvénic, via the corresponding Mach numbers. The sound speed and Alfvén speed are around 165 km s^{−1} and 587 km s^{−1} 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 setup 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 colourcoded according to the magnitude of the B_{0θ} 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 QRinvert, 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 QRinvert 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 nonadiabatic effects. The spectrum of the benchmark equilibrium, a GoldHoyle magnetic field with a constant axial and rsquared azimuthal flow profile, is given in the top panel of Fig. 3. There are insets showing the two slow continua, the backwardpropagating 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 forwardpropagating slow continuum. No discrete modes appear next to the backwardpropagating slow continuum or the backwardpropagating 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: zoomin image of the upper edge of the forwardpropagating Alfvén continuum. Three discrete Alfvén modes are marked. The corresponding eigenfunctions are shown in the middle right panel. Bottom left panel: forwardpropagating 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 backwardpropagating 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 forwardpropagating 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 antiSturmian 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 forwardpropagating 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 backwardpropagating 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 nonconstant, leading to this maximum. The backwardpropagating 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, , evaluated at the location of the extremum is larger than 0.25, discrete modes arise. For the Alfvén continuum considered here we have . This is larger than 0.25, and hence we would indeed expect discrete Alfvén modes to cluster towards the maximum of the forwardpropagating Alfvén continuum.
3.2.2. The nonadiabatic case
We now extend on the literature results by investigating the spectrum of a plasma both influenced by nonadiabatic 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 nonsquare integrable shape, as can be seen most easily in the inset.
Fig. 5. Left panel: zoomedin 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 nonsquare 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
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 forwardpropagating 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 nonadiabatic effects. The damping of discrete Alfvén modes by nonadiabaticity 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 antiSturmian sequence, as to be expected. The backwardpropagating Alfvén continuum has no extremum and no discrete modes cluster towards its edges. There also appears a discrete slow mode next to the forwardpropagating 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 zoomedin images of the backwardpropagating Alfvén continuum and the two slow continua. The two panels below the main spectrum are zoomins onto the forwardprogating Alfvén continuum, with another zoomin 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
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 flowincluded 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
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 piecewise 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 fieldaligned, 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
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 minussign 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 forwardpropagating, while others will be backwardpropagating.
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 zoomin on the unstable wavenumber range. 
Parallel thermal conduction has the same effect of damping the thermal modes, based on a cutoff 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 wellknown 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 nonadiabatic 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 secondorder differential equation Eq. (46), which can be written in the form
with the coefficients f(r) and g(r) dependent on the coefficients C_{0} 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
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^{2}L^{2} ≫ 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
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 firstorder differential equations, Eqs. (37)–(40), yields
Multiplying both sides with C_{0} and substituting its expression, Eq. (37), gives
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 nonadiabatic, radiative effects. Both kinds of modes are solutions of the same equation Eq. (49). We consider a Dopplershifted thermal or slow continuum that has an internal extremum at r_{0}, such that and . 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
where ϵ is a small, complex perturbation. The lefthand side of the local dispersion relation becomes
where every equilibrium quantity is evaluated at r = r_{0}. The lefthand side is thus accurate up to firstorder of ϵ. The righthand side of the local dispersion relation is expanded to zerothorder of ϵ and evaluated at r = r_{0}. After a lengthy derivation ϵ can be determined. The approximation of the eigenfrequency is finally given by
with the continuum frequency around which was expanded. The parameter ζ is given by
with
All parameters are evaluated at the local extremum of the continuum, (r_{0}, ). Furthermore, we used the fact that 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 nonadiabatic 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 spreadout profile instead of a very localised perturbation.
Note that a similar expression can be derived for the nonadiabatic 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. Setup
We base the three equilibria on the benchmark setup 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 rightmost 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 antiSturmian 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 WKBapproximation 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 roundoff 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
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 nonlinear behaviour. Nevertheless, it is a useful and simple way to gain information without performing full nonlinear 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 ultralocalised 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 nonaxisymmetric, 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 zr 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^{1}. 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 v_{r} 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 Fourieranalysis of the nonadiabatic MHD equations in a cylindrical coordinate system with a generic helical flow. The set of equations is transformed into two firstorder 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 forcefree 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 hyperlocalised, 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. Threedimensional 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 wellknown work by Edwin & Roberts (1983), this kind of setup 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 nonlinear 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 setup, 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 threephase 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 FWONSFC 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 ERCADG 2018). This research is further supported by Internal funds KU Leuven, project C14/19/089 TRACESpace.
References
 Appert, K., Gruber, R., & Vaclavik, J. 1974, Phys. Fluids, 17, 1471 [Google Scholar]
 Banerjee, D., Krishna Prasad, S., Pant, V., et al. 2021, Space Sci. Rev., 217, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Barczynski, K., Harra, L., Kleint, L., Panos, B., & Brooks, D. H. 2021a, A&A, 651, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barczynski, K., Schmieder, B., Peat, A. W., et al. 2021b, A&A, 653, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barczynski, K., Harra, L., Schwanitz, C., et al. 2023, A&A, 673, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Begelman, M. C., & McKee, C. F. 1990, ApJ, 358, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Blokland, J. W. S., van der Swaluw, E., Keppens, R., & Goedbloed, J. P. 2005, A&A, 444, 337 [Google Scholar]
 Bondeson, A., Iacono, R., & Bhattacharjee, A. 1987, Phys. Fluids, 30, 2167 [Google Scholar]
 Braginskii, S. I. 1965, Rev. Plasma Phys., 1, 205 [Google Scholar]
 Brughmans, N., Jenkins, J. M., & Keppens, R. 2022, A&A, 668, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chitta, L. P., Solanki, S. K., Peter, H., et al. 2021, A&A, 656, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claes, N. 2022, Ph.D. thesis, Faculty of Science, KU Leuven, Belgium [Google Scholar]
 Claes, N., & Keppens, R. 2019, A&A, 624, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claes, N., & Keppens, R. 2021, Sol. Phys., 296, 143 [Google Scholar]
 Claes, N., & Keppens, R. 2023, Comput. Phys. Commun., 291, 108856 [Google Scholar]
 Claes, N., De Jonghe, J., & Keppens, R. 2020a, ApJS, 251, 25 [Google Scholar]
 Claes, N., Keppens, R., & Xia, C. 2020b, A&A, 636, A112 [EDP Sciences] [Google Scholar]
 Colgan, J., Abdallah, J. Jr., Sherrill, M. E., et al. 2008, ApJ, 689, 585 [Google Scholar]
 Cox, D. P. 2005, ARA&A, 43, 337 [Google Scholar]
 De Jonghe, J. 2023, Ph.D. thesis, KU Leuven, Belgium [Google Scholar]
 De Jonghe, J., Claes, N., & Keppens, R. 2022, J. Plasma Phys., 88, 905880321 [Google Scholar]
 De Moortel, I., Antolin, P., & Van Doorsselaere, T. 2015, Sol. Phys., 290, 399 [NASA ADS] [CrossRef] [Google Scholar]
 Edwin, P. M., & Roberts, B. 1983, Sol. Phys., 88, 179 [Google Scholar]
 Field, G. B. 1965, ApJ, 142, 531 [Google Scholar]
 Frieman, E., & Rotenberg, M. 1960, Rev. Mod. Phys., 32, 898 [Google Scholar]
 Furth, H. P., Killeen, J., & Rosenbluth, M. N. 1963, Phys. Fluids, 6, 459 [Google Scholar]
 Gibson, S. E. 2018, Liv. Rev. Sol. Phys., 15, 7 [Google Scholar]
 Goedbloed, J. P. 1984, Phys. D Nonlinear Phenomena, 12, 107 [Google Scholar]
 Goedbloed, J. P., & Sakanaka, P. H. 1974, Phys. Fluids, 17, 908 [Google Scholar]
 Goedbloed, J., Keppens, R., & Poedts, S. 2019, Magnetohydrodynamics of Laboratory and Astrophysical Plasmas (Cambridge: Cambridge University Press) [Google Scholar]
 Gold, T., & Hoyle, F. 1960, MNRAS, 120, 89 [NASA ADS] [Google Scholar]
 Goossens, M., Hollweg, J. V., & Sakurai, T. 1992, Sol. Phys., 138, 233 [Google Scholar]
 Gunár, S., Dudík, J., Aulanier, G., Schmieder, B., & Heinzel, P. 2018, ApJ, 867, 115 [Google Scholar]
 Gunár, S., Labrosse, N., Luna, M., et al. 2023, Space Sci. Rev., 219, 33 [Google Scholar]
 Hain, K., & Lüst, R. 1958, Zeitschrift Naturforschung Teil A, 13, 936 [Google Scholar]
 Hameiri, E. 1981, J. Math. Phys., 22, 2080 [Google Scholar]
 Hameiri, E. 1983, Phys. Fluids, 26, 230 [Google Scholar]
 Hermans, J., & Keppens, R. 2021, A&A, 655, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jenkins, J. M., & Keppens, R. 2021, A&A, 646, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jennings, R. M., & Li, Y. 2021, MNRAS, 505, 5238 [NASA ADS] [CrossRef] [Google Scholar]
 Jerčić, V., & Keppens, R. 2023, A&A, 670, A64 [Google Scholar]
 Keppens, R., van der Linden, R. A. M., & Goossens, M. 1993, Sol. Phys., 144, 267 [Google Scholar]
 Kucera, T. A. 2015, in Solar Prominences, eds. J. C. Vial, & O. Engvold, Astrophys. Space Sci. Lib., 415, 79 [Google Scholar]
 Li, X., Morgan, H., Leonard, D., & Jeska, L. 2012, ApJ, 752, L22 [Google Scholar]
 Liakh, V., & Keppens, R. 2023, ApJ, 953, L13 [Google Scholar]
 Long, D. M., Chitta, L. P., Baker, D., et al. 2023, ApJ, 944, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Luna, M., MorenoInsertis, F., & Priest, E. 2015, ApJ, 808, L23 [Google Scholar]
 Luna, M., Priest, E., & MorenoInsertis, F. 2018, ApJ, 863, 147 [Google Scholar]
 Onishchenko, O. G., Fedun, V., Smolyakov, A., et al. 2018, Phys. Plasmas, 25, 054503a [Google Scholar]
 Orozco Suárez, D., Asensio Ramos, A., & Trujillo Bueno, J. 2012, ApJ, 761, L25 [Google Scholar]
 Panasenco, O., Martin, S. F., & Velli, M. 2014, Sol. Phys., 289, 603 [Google Scholar]
 Parker, E. N. 1953, ApJ, 117, 431 [Google Scholar]
 Pettit, E. 1932, ApJ, 76, 9 [Google Scholar]
 Priest, E. R. 1982, Solar Magnetohydrodynamics (D. Reidel Publishing Company) [Google Scholar]
 Rosner, R., Tucker, W. H., & Vaiana, G. S. 1978, ApJ, 220, 643 [Google Scholar]
 Schmieder, B., Raadu, M. A., & Wiik, J. E. 1991, A&A, 252, 353 [NASA ADS] [Google Scholar]
 Schure, K. M., Kosenko, D., Kaastra, J. S., Keppens, R., & Vink, J. 2009, A&A, 508, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sharma, P., Parrish, I. J., & Quataert, E. 2010, ApJ, 720, 652 [NASA ADS] [CrossRef] [Google Scholar]
 Shen, Y. 2021, Proc. R. Soc. London Ser. A, 477, 217 [NASA ADS] [Google Scholar]
 Skirvin, S. J., Fedun, V., Silva, S. S. A., et al. 2023, MNRAS, 518, 6355 [Google Scholar]
 Soler, R., Ballester, J. L., & Parenti, S. 2012, A&A, 540, A7 [Google Scholar]
 Spitzer, L. 1962, Physics of Fully Ionized Gases (Courier Corporation) [Google Scholar]
 Su, Y., Wang, T., Veronig, A., Temmer, M., & Gan, W. 2012, ApJ, 756, L41 [Google Scholar]
 Su, Y., Gömöry, P., Veronig, A., et al. 2014, ApJ, 785, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Suydam, B. R. 1958, Proc. Second U.N. Intern. Conf. Peaceful Uses At. Energy, Geneva, 31, 157 [Google Scholar]
 Tziotziou, K., Scullion, E., Shelyag, S., et al. 2023, Space Sci. Rev., 219, 1 [NASA ADS] [CrossRef] [Google Scholar]
 van der Linden, R. A. M., & Goossens, M. 1991, Sol. Phys., 134, 247 [NASA ADS] [CrossRef] [Google Scholar]
 Veilleux, S., Maiolino, R., Bolatto, A. D., & Aalto, S. 2020, A&A Rev., 28, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, C., Blokland, J. W. S., Keppens, R., & Goedbloed, J. P. 2004, J. Plasma Phys., 70, 651 [Google Scholar]
 Webb, D. F. 2015, in Solar Prominences, eds. J. C. Vial, & O. Engvold, Astrophys. Space Sci. Lib., 415, 411 [Google Scholar]
 Wedemeyer, S., Scullion, E., Rouppe van der Voort, L., Bosnjak, A., & Antolin, P. 2013, ApJ, 774, 123 [Google Scholar]
 Yang, Z., Tian, H., Peter, H., et al. 2018, ApJ, 852, 79 [Google Scholar]
All Tables
Values for the physical quantities used in the benchmark case. B_{0c} is the magnetic field strength at the axis; μ is the inverse pitch of the magnetic field; ρ_{0b}, τ, and T_{0c} denote the density at the boundary, the exponent of the density profile, and the temperature at the axis, respectively.
All Figures
Fig. 1. Equilibrium profiles of the benchmark setup 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 colourcoded according to the magnitude of the B_{0θ} 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: zoomin image of the upper edge of the forwardpropagating Alfvén continuum. Three discrete Alfvén modes are marked. The corresponding eigenfunctions are shown in the middle right panel. Bottom left panel: forwardpropagating 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 backwardpropagating 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: zoomedin 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 nonsquare 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 zoomedin images of the backwardpropagating Alfvén continuum and the two slow continua. The two panels below the main spectrum are zoomins onto the forwardprogating Alfvén continuum, with another zoomin 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 zoomin 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 zr 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 v_{r} component of the velocity perturbation. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.