© The Authors 2025

1. Introduction

Magnetohydrodynamic (MHD) waves are fundamental to both experimental and astrophysical plasma physics, and many phenomena observed in the solar atmosphere are attributed to their influence. It is well theorized that these waves play a crucial role in the transfer of energy from the lower to the upper layers of the solar atmosphere, offering a key explanation for coronal heating (see e.g. McIntosh et al. 2011; De Moortel & Browning 2015; Khomenko & Collados 2015; Van Doorsselaere et al. 2020). Moreover, MHD waves serve as natural seismological tools, providing valuable insights into the physical properties of solar plasmas (Liu & Ofman 2014) (for a recent review see e.g. Nakariakov et al. 2024). They are also implicated in the acceleration of solar winds (Cranmer 2012; Vasheghani Farahani et al. 2021), and their non-linear effects may contribute to the geometry of solar phenomena such as collimation of coronal jets (Vasheghani Farahani & Hejazi 2017).

A key manifestation of MHD waves in the solar atmosphere is through magnetoacoustic (MA) modes, classified by their symmetry properties such as sausage and kink waves, which manifest in both fast and slow modes as fundamental examples of MHD wave phenomena. These modes are characterized by distinct oscillation patterns, frequencies, and phase speeds along magnetic structures, and have been widely observed in various solar structures, including coronal loops, jets, and prominences (Nakariakov et al. 1999; Van Doorsselaere et al. 2008; Dorotovič et al. 2014; De Moortel et al. 2015; Moreels et al. 2015a; Nakariakov & Kolotkov 2020). The deviation from thermal equilibrium caused by MHD waves, leading to enhanced plasma heating over cooling in the upper layers of the solar atmosphere through various damping mechanisms, has been extensively studied (for a comprehensive review, see Van Doorsselaere et al. 2020). Conversely, the heating can influence the dynamics of MHD waves, detectable through a broad range of spectroscopic observations, from microwave to X-ray emission lines (Kolotkov et al. 2019). The observation of MA waves in the solar atmosphere serves as an effective seismological tool, where the wave properties “such as propagation speeds, periods, amplitudes, and damping times” enable the determination of essential plasma parameters, including density, temperature, magnetic field strength, and heating function. Several theoretical and observational studies have highlighted the potential of MHD seismology to reveal physical properties of solar atmospheric plasmas (Andries et al. 2009; Verwichte et al. 2009; Kolotkov et al. 2020; Nakariakov et al. 2021).

Slow MA waves, whether propagating or standing, are fundamental components of the MHD wave spectrum observed and theorized in the solar atmosphere. These waves are detected via Doppler shifts or intensity oscillations in microwave, extreme-ultraviolet (EUV), and X-ray emission bands (Nakariakov et al. 2000a, 2019; Wang et al. 2021). Due to their ubiquitous presence and relative ease of detection compared to other MHD waves, they have attracted considerable attention as diagnostic tools for probing solar plasma characteristics (Wang et al. 2003; Nakariakov & Verwichte 2005; de Moortel 2009). One notable example of slow MA waves is the so-called SUMER oscillations, named after observations by the Solar Ultraviolet Measurements of Emitted Radiation (SUMER) instrument on board the Solar and Heliospheric Observatory (SoHO) (Wang et al. 2002; Wang 2011). These oscillations, also often observed with instruments such as the Interface Region Imaging Spectrograph (IRIS) (Li et al. 2017), manifest as rapidly damped, standing, slow-mode oscillations. They are commonly characterized by periodic or quasi-periodic fluctuations in temperature and density within hot coronal loops, providing crucial insights into the physical processes governing solar dynamics.

A particularly interesting aspect of wave behaviour that can significantly influence the propagation characteristics of MA modes, including their damping and amplification within magnetic flux tubes, lies in the role of thermal misbalance, a non-adiabatic process resulting from the disparity between the timescales of heating and radiative cooling perturbations in the plasma. The phenomenon of thermal instability in plasma, first studied by Field (1965), revealed that thermal equilibrium can be inherently unstable in a uniform medium, leading to the formation of condensations of a higher density and lower temperature. This study provided foundation under how thermal instability operates under a wide range of conditions, influencing not only solar phenomena but also broader astrophysical structures such as planetary nebulae and interstellar mediums. Earlier studies on thermal misbalance, such as those done by Nakariakov et al. (2000b), Kumar et al. (2016), Kolotkov et al. (2019), and Zavershinskii et al. (2019), further developed the understanding of wave-induced thermal misbalance under the approximations of infinite magnetic field and weak non-adiabaticity, highlighting its effects on the evolution of slow MA waves.

Recent studies have shown that slow MA waves, particularly in solar coronal loops, are strongly influenced by thermal misbalance. This misbalance causes wave amplification, attenuation, and dispersion, depending on the dominant thermodynamic processes in the plasma. Observations have shown that slow MA waves undergo frequency-dependent damping, which traditional damping theories fail to explain. This opens the possibility for using thermal misbalance to diagnose unknown coronal heating mechanisms through detailed seismological analysis of wave properties (Nakariakov et al. 2017; Kolotkov et al. 2019; Kolotkov & Nakariakov 2022). In both propagating and standing MA waves, heating-cooling misbalance can result in different damping regimes, oscillatory patterns, and even wave amplification, depending on the plasma conditions (Zavershinskii et al. 2019; Belov et al. 2020; Kolotkov et al. 2020; Agapova et al. 2022). This has led to ongoing efforts to better understand how wave properties such as phase shifts, group speeds, and damping lengths are affected by both the plasma’s non-adiabatic processes and its magnetic structuring (Kolotkov et al. 2021; Prasad et al. 2022; Ibañez & Ballester 2022).

However, the solar corona is a complex environment where the infinite magnetic field approximation is not universally applicable. Variations in magnetic field strength, as reported in the works by Afanasyev & Nakariakov (2015) and Nakariakov et al. (2017), suggest that the assumption of an overwhelmingly dominant magnetic pressure may not always hold. Thus, it becomes crucial to explore scenarios where the magnetic field strength is finite, corresponding to non-zero plasma-β conditions, and to assess their impact on wave dynamics. This consideration is particularly relevant since coronal heating, believed to be primarily driven by magnetic fields, can have its efficiency modulated by the field strength itself, consequently affecting the thermal misbalance and the associated behaviour of MHD waves (e.g. Duckenfield et al. 2021). Building on this, Kolotkov et al. (2019, 2020) presented compelling evidence of rapidly decaying slow MA waves due to thermal misbalance observed in hot coronal loops, demonstrating its potential for constraining the coronal heating function through observed damping rates of slow modes. These studies further conclude that coronal plasmas, in all scenarios, should always be regarded as an active medium, powered by thermal misbalance and engaged in a continuous exchange of energy with MHD waves (Kolotkov et al. 2021).

In addition to affecting slow MA waves, thermal misbalance also influences other MHD wave modes, including Alfvén waves. These waves are subject to non-linear effects due to the interaction of torsional and shear Alfvén waves with the surrounding plasma. Studies have shown that the induced longitudinal velocity perturbations in such waves are modulated by thermal misbalance, with the effect being more pronounced in shear waves than torsional ones (Belov et al. 2022). Furthermore, the back-reaction caused by the perturbed thermal equilibrium in the solar corona results in thermal instabilities, which play a key role in the formation of fine thermal structures in coronal loops (Kolotkov et al. 2023). These instabilities are stabilized to some extent by thermal conduction, especially in hot active regions, but the sensitivity of slow wave stability to coronal heating functions presents an opportunity for further seismological diagnostics (Zavershinskii et al. 2023; Ballester et al. 2024). Overall, this demonstrates the effect of thermal misbalance on MHD wave dynamics and highlight its importance in constraining coronal heating mechanisms and explaining solar atmospheric behaviour (Kolotkov & Nakariakov 2022; Riashchikov et al. 2023).

The implications of these findings are significant, potentially resolving long-standing discrepancies in the observed frequency-dependent damping of slow MA waves in various coronal structures, as reported by Mariska (2006), Gupta (2014), Krishna Prasad et al. (2014, 2018), and Kolotkov & Nakariakov (2022). Numerical studies have demonstrated that damping lengths vary with wave periods (Gupta 2014; Mandal et al. 2016). Furthermore, Arregui et al. (2023) highlighted the importance of considering both thermal conduction and thermal misbalance in explaining the damping of slow MA waves in hot coronal loops. The consideration of thermal misbalance as a damping mechanism may reconcile observed phase shifts between density and temperature, as well as the growth of the poly-tropic index, with the theoretical predictions done by Owen et al. (2009), Van Doorsselaere et al. (2011), and Krishna Prasad et al. (2018). Additionally, Kolotkov (2022) explored the use of slow wave seismology to compare the effects of weak and strong thermal conduction, aiming to derive a more accurate value for the effective adiabatic index in coronal plasma.

In this study, we develop an analytical model to investigate the impact of thermal misbalance in none zero plasma-β on a spectrum of MHD modes in coronal magnetic structures, using coronal seismology as a tool to refine the understanding of the heating function. The MHD modes analysed in this research cover a range of phase speeds, from fast to slow MA waves, propagating within solar magnetic cylinders. By expanding the non-adiabatic linear dispersion relation to the long wavelength limit, we aim to provide deeper insight into the damping times of slow MA waves in the solar corona and the dependence of wave damping on the heating and cooling functions. The context of the present study and its associated model are inspired by a series of works by Edwin & Roberts (1983), van der Linden & Goossens (1991), Ireland et al. (1992), and Kolotkov et al. (2019), to advance the field by incorporating the geometrical constraints of a cylindrical configuration with consideration of finite magnetic field, resulting transversal structuring, alongside the influence of thermal misbalance of arbitrary intensity, which are intrinsic to the solar coronal environment.

2. Governing equations and dispersion relation

We consider a magnetic cylinder in a uniform medium characterized by equilibrium finite magnetic field B₀ = (0, 0, B₀), pressure p₀, and density ρ₀. Our analysis employs the ideal MHD equations while also considering the influence of thermal misbalance between heating and optically thin cooling radiation:

$$\begin{array}{cc}& \rho [\frac{\partial \mathbf{v}}{\partial t}+(\mathbf{v}.\mathrm{\nabla })\mathbf{v}]=-\mathrm{\nabla }p+\frac{1}{4\pi }(\mathrm{\nabla }\times \mathit{B})\times \mathit{B},\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }.(\rho \mathbf{v})=0,\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times (\mathbf{v}\times \mathit{B}),\hfill \end{array}$$(3)

$$\begin{array}{cc}& p=\frac{k_{\mathrm{B}}\rho T}{m},\hfill \end{array}$$(4)

$$\begin{array}{cc}& \frac{{\rho }^{\gamma }}{\gamma -1}\frac{\mathrm{d}}{\mathrm{d}t}(\frac{p}{{\rho }^{\gamma }})=-\rho Q(\rho ,T),\hfill \end{array}$$(5)

$$\begin{array}{cc}& \mathrm{\nabla }·\mathit{B}=0,\hfill \end{array}$$(6)

where, Q(ρ, T) represents the discrepancy between the heating and the radiative cooling processes, defined as L(ρ, T)−H(ρ, T) (Parker 1953). The cooling function, L(ρ, T), and the heating function, H(ρ, T), are both expressed as general power-law functions of density and temperature as below, with unspecified exponents, each of which must be independently determined for different regions within the solar atmosphere,

$$\begin{array}{c}\hfill L(\rho ,T)=\chi \rho T^{\beta },H(\rho ,T)=h{\rho }^a{T}^b.\end{array}$$(7)

The optically thin radiative cooling function, L(ρ, T), is readily determinable, thanks to the CHIANTI atomic database (since Dere et al. 1997). In this study, we employed the latest version, 10.1 (Dere et al. 2023), to calculate the temperature exponent, β, and the coefficient χ in the cooling function. Regarding the unspecified heating function, it is necessary to consider the constraints governing the values of the power-law exponents a and b (Kolotkov et al. 2020). Assuming an isothermal initial equilibrium, we set Q(ρ₀, T₀) = 0, which is essential for determining the coefficient of h in the heating function. It is noteworthy that some heating mechanisms may exhibit time-varying components on very short timescales (Reale 2016) or dependence on the magnetic field (Duckenfield et al. 2021). In this study, we adopt a non-specified model for both heating and cooling functions. Our primary objective is to develop a descriptive model to understand how arbitrary forms of thermal misbalance functions, each characterized by its own timescales, influence MHD modes within a magnetic cylindrical waveguide. While this approach may lead to precise forms of the heating function, such details are beyond the scope of this article.

At this point, it is important to note that both the general heating and cooling functions depend on density, and the heating coefficient, h, is directly linked to the background density profile, ρ₀, where under realistic solar atmospheric conditions, should be considered stratified. Consequently, assuming a uniform background plasma may be inaccurate, and the effects of stratification should be incorporated into our mathematical models. However, it is essential to recognize that the MA waves, particularly slow MA waves considered in this study, are characterized by compressive oscillations and represent the MHD modes significantly influenced by thermal non-adiabaticity in the solar corona (Zavershinskii et al. 2019; Kolotkov et al. 2019, 2020; Agapova et al. 2022). These waves exhibit wavelengths ranging from a few thousand to several tens of thousands of kilometres. Notably, this range is much smaller than the length scale of stratification, which reaches a minimum value of approximately 50 Mega-metres (Nakariakov et al. 2000c; Andries et al. 2005). Therefore, the assumption of a uniform medium remains justified.

Upon linearizing the governing equations by introducing small amplitude perturbations to each physical parameter around the initial equilibrium, we focus on the variable that represents compressibility in the system (∇ ⋅ v) (Edwin & Roberts 1983). This procedure results in a linear set of operators that govern the physical characteristics of the system, as outlined below

$$\begin{array}{c}\hfill (\widehat{G}+\widehat{Q})\mathrm{\Delta }=0.\end{array}$$(8)

Let Δ = ∇ ⋅ v, and the subsequent expressions represent each operator:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \widehat{G}& \equiv \frac{{\partial }^2}{\partial t^2}[\frac{{\partial }^2}{\partial t^2}-(V_{\mathrm{A}}^2+C_{\mathrm{s}}^2){\mathrm{\nabla }}^2]+V_{\mathrm{A}}^2{C}_{\mathrm{s}}^2\frac{{\partial }^2}{\partial z^2}{\mathrm{\nabla }}^2,\hfill \\ \hfill \widehat{Q}& \equiv Q_{\rho }\left[{\rho }_0(\gamma -1)\frac{\partial }{\partial t}{\mathrm{\nabla }}^2\right]\hfill \\ \hfill & +Q_{\mathrm{T}}[\frac{1}{C_{\mathrm{V}}}\frac{\partial }{\partial t}(\frac{{\partial }^2}{\partial t^2}-V_{\mathrm{A}}^2{\mathrm{\nabla }}^2+V_{\mathrm{A}}^2\frac{{\partial }^2}{\partial z^2})-(\gamma -1)T_0\frac{\partial }{\partial t}{\mathrm{\nabla }}^2].\hfill \end{array}\end{array}$$(9)

Each differential operator in Eq. (8) corresponds to a specific physical effect within the system. The operator $\widehat{G}$ represents the geometrical constraints accounted for in the dispersion relation, which can incorporate the cylindrical magnetic configuration (see also Eq. (3a) in Edwin & Roberts 1983). The operator $\widehat{Q}$ describes the effect of thermal misbalance between heating and radiative cooling processes as well. For further details on the derivation of the mathematical expressions of the operators in Eq. (9), refer to Appendix A.

In Eq. (9), V_A and C_s represent the Alfvén and acoustic speeds in the system, respectively. They are defined by the relations V_A = B₀(4πρ₀)^−1/2 and $C_{\mathrm{s}}=\sqrt{\gamma k_{\mathrm{B}}{T}_0{m}^{-1}}$, where B₀ is the equilibrium magnetic field strength, ρ₀ is the background density, γ is the adiabatic index, k_B is the Boltzmann constant, T₀ is the equilibrium temperature, and m is the mean particle mass. Additionally, the specific heat capacity at constant volume, denoted as C_V = (γ − 1)⁻¹k_B/m, is determined using the standard adiabatic index (γ = 5/3). It is important to note that the polytropic index of the system deviates from the standard adiabatic index when non-adiabatic terms are included in the energy equation (Eq. (5)). Adjusting this index becomes necessary to account for the effects on the index in subsequent calculations. However, the focus of this analytical study is on the complex influence of thermal misbalance on MHD waves, and therefore, this topic is reserved for future investigation. For more detailed discussions on the effective adiabatic index in coronal plasma, readers are encouraged to refer to Van Doorsselaere et al. (2011), Krishna Prasad et al. (2018), and Kolotkov (2022).

Moreover, in Eq. (9), We express the rate of change of the heating-cooling function with respect to density at constant temperature as Q_ρ, and with respect to temperature at constant density as Q_T, as follows:

$$\begin{array}{c}\hfill Q_{\rho }\equiv \frac{\partial Q}{\partial \rho }{{|}_{\mathrm{T}},Q_{\mathrm{T}}\equiv \frac{\partial Q}{\partial T}|}_{\rho }.\end{array}$$(10)

In the thin flux tube approximation, neglecting parallel thermal conduction reduces Eq. (10) in Kolotkov et al. (2021) to Eq. (8). It is also deducible that neglecting the effect of thermal misbalance in Eq. (8), results in the presence of only the geometrical dispersion operator, reducing it to Eq. (3a) in Edwin & Roberts (1983). At this point, it is noteworthy that if we consider the governing equations in the absence of magnetic effects (i.e. in the hydrodynamic limit) while accounting for thermal conduction, the linearized perturbed equations (Eqs. (A.1, A.2, A.4, A.5) in Appendix A) regenerate the forms of Eqs. (11–14) in Field (1965). This ends up resulting in the dispersion relation given by Eq. (15) in Field (1965), for an exponentially propagating wave solution in an infinite, uniform medium.

By defining the characteristic timescales τ₁ and τ₂ as follows, it becomes possible to quantify the effects of cooling and heating processes on MA waves, for instance, regarding attenuation, amplification, and phase speed changes (Zavershinskii et al. 2019),

$$\begin{array}{c}\hfill {\tau }_1=\gamma C_{\mathrm{V}}/[Q_{\mathrm{T}}-\frac{{\rho }_0}{T_0}{Q}_{\rho }],{\tau }_2=C_{\mathrm{V}}/Q_{\mathrm{T}}.\end{array}$$(11)

Additionally, in the infinite magnetic field approximation–where the background magnetic field is considered infinitely strong, allowing the perturbation of the magnetic fields in the system to be assumed negligible –the geometrical dispersion effect on slow modes can be disregarded. In this scenario, the governing equations are reduced to one-dimensional forms along the $\widehat{z}$ axis. Therefore, incorporating the characteristic timescale for thermal conduction, τ_cond, which is proportional to the square of wavelength (see Zavershinskii et al. 2019), allows for the reduction of Eq. (8) to Eq. (5) in Kolotkov et al. (2019).

In this study, we focus on a segment of the magnetic flux tube located within the coronal plasma, where the τ_cond is significantly longer than the characteristic timescales associated with thermal misbalance (Kolotkov et al. 2019). Moreover, due to the location in the coronal region of the solar atmosphere, we assume a uniform plasma temperature aligned with the background magnetic field, with a negligible background temperature gradient. This assumption, widely adopted in theoretical coronal plasma studies (e.g. Owen et al. 2009; Mandal et al. 2016; Wang et al. 2018; Kolotkov et al. 2020), implies that, compared to thermal misbalance, the effects of thermal conduction can be neglected. This allows us to focus exclusively on the influence of thermal misbalance on wave damping (pure thermal misbalance case in Kolotkov et al. 2019).

In this context, it is where to introduce two dimensionless characteristic timescales, obtained by normalizing τ₁ and τ₂ by the wave period (P = 2π/ω_r, ω_r = ℜ𝔢(ω)), to account for the thermal misbalance. These are defined as $\overline{{\tau }_1}=\frac{{\tau }_1}{P}$ and $\overline{{\tau }_2}=\frac{{\tau }_2}{P}$. The dependence of these characteristic timescales on the wave period is significant, especially in scenarios where thermal misbalance perturbation exhibits a wave-induced nature (Zavershinskii et al. 2019; Kolotkov et al. 2019; Belov et al. 2021). Through introducing these two dimensionless period-dependent timescales, the impact of thermal misbalance on a wave, characterized by physical parameters such as frequency or period, can be directly correlated with the wave’s susceptibility to thermal misbalance.

Figure 1 illustrates the dependence of $\overline{{\tau }_1}$ and $\overline{{\tau }_2}$ on temperature and wave period. The cooling function is derived from the CHIANTI database v. 10.1 for an equilibrium density (ρ₀) of 2.3 × 10⁻¹⁵ g/cm³, while the heating function is parametrized with power-law exponents a = −0.5 and b = −3. It is evident that these characteristic timescales satisfy the conditions $\overline{{\tau }_2}<\overline{{\tau }_1}$ and $\overline{{\tau }_1}-\overline{{\tau }_2}>0$, ensuring that thermal misbalance acts as a damping factor on MHD waves (Kolotkov et al. 2020). Furthermore, the difference between the equilibrium densities of the external and internal media may cause variations in the characteristic timescales of thermal misbalance. However, this study assumes that such discrepancies in equilibrium densities do not significantly affect the characteristic timescales of thermal misbalance within and outside the tube. Moreover, for waves with longer periods and environmental temperatures ranging from 1 to 5 MK, as shown in Fig. 1, the values of $\overline{{\tau }_1}$ and $\overline{{\tau }_2}$ are minimized, and vice versa.

Fig. 1. Dependence of dimensionless characteristic timescales, ${\overline{\tau }}_1$ (a) and ${\overline{\tau }}_2$ (b), on temperature and wave period calculated using Eq. (11), normalized by the wave period. This analysis considers a heating model characterized by power indices a = −0.5 and b = −3 for density and temperature, respectively.

Next, expressing the Laplacian operator in cylindrical coordinates (r, θ, z) within the set of differential operators in Eq. (9) as

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\equiv \frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}+\frac{1}{r^2}\frac{{\partial }^2}{\partial {\theta }^2}+\frac{{\partial }^2}{\partial z^2}·\end{array}$$(12)

We adopt the Fourier exponential form to describe the solution’s dependence on the coordinates and time as:

$$\begin{array}{c}\hfill R(r)exp(ikz+im\theta -i\omega t),\end{array}$$(13)

where k, m, and ω represent the longitudinal wavenumber, azimuthal wavenumber, and wave frequency, respectively. Substituting this form into Eq. (8) leads to a Bessel equation for R(r) as

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^{2}R}{\mathrm{d}{r}^2}+\frac{1}{r}\frac{\mathrm{d}R}{\mathrm{d}r}-({\kappa }^2+\frac{m^2}{r^2})R=0.\end{array}$$(14)

Here, κ represents the interior radial wavenumber and is a complex function of wave frequency (ω) derived from the following expression

$$\begin{array}{c}\hfill {\kappa }^2=\frac{(k^2{C}_{\mathrm{s}}^2-{\omega }^2)(k^2{V}_{\mathrm{A}}^2-{\omega }^2)-\frac{i{\omega }^2}{2\pi \overline{{\tau }_2}}(\frac{\overline{{\tau }_2}}{\overline{{\tau }_1}}{k}^2{C}_{\mathrm{s}}^2-{\omega }^2)}{(k^2{C}_{\mathrm{T}}^2-{\omega }^2)(C_{\mathrm{s}}^2+V_{\mathrm{A}}^2)-\frac{i{\omega }^2}{2\pi \gamma \overline{{\tau }_2}}(C_{\mathrm{s}}^2+\gamma V_{\mathrm{A}}^2)}·\end{array}$$(15)

The tube speed, denoted as C_T and expressed as $\sqrt{\frac{C_{\mathrm{s}}^2{V}_{\mathrm{A}}^2}{C_{\mathrm{s}}^2+V_{\mathrm{A}}^2}}$, is commonly referred to as the cusp speed (V_C) (Roberts & Webb 1978). However, this study demonstrates that the cusp speed deviates from the tube speed due to the influence of non-adiabatic effects within the system, corroborating earlier findings illustrated in Figs. 1 and 3 of Belov et al. (2021), where it is also referred to as C_TQ.

Assuming no energy exchange between the interior and exterior of the magnetic flux tube implies that the oscillations within the tube at the boundary are effectively isolated from external energy sources or sinks. This isolation allows the oscillations inside to be treated as independent of external drivers or damping effects. Consequently, we also arrive at a Bessel equation with distinct Alfvén, sonic, and tube speeds for the exterior region. The solutions to these Bessel equations, while maintaining continuity at r = R (the radius of the tube) and ensuring finite solutions at both r = 0 and r → ∞, are then obtained by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill R(r)=\{\begin{array}{c}{A}_{\mathrm{in}}{J}_m({\kappa }_{\mathrm{in}}r)r<R,\hfill \\ A_{\mathrm{out}}{K}_m({\kappa }_{\mathrm{out}}r)r>R,\hfill \end{array}\end{array}\end{array}$$(16)

where A_in and A_out are complex constants, and under the condition that κ_in² = −κ² > 0 (Edwin & Roberts 1983), J_m and K_m are the Bessel function of the first kind and the modified Bessel function of the second kind, respectively, and κ_out is the external radial wavenumber, given by

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{out}}^2=\frac{(k^2{C}_{\mathrm{se}}^2-{\omega }^2)(k^2{V}_{\mathrm{Ae}}^2-{\omega }^2)-\frac{i{\omega }^2}{2\pi \overline{{\tau }_2}}(\frac{\overline{{\tau }_2}}{\overline{{\tau }_1}}{k}^2{C}_{\mathrm{se}}^2-{\omega }^2)}{(k^2{C}_{\mathrm{Te}}^2-{\omega }^2)(C_{\mathrm{se}}^2+V_{\mathrm{Ae}}^2)-\frac{i{\omega }^2}{2\pi \gamma \overline{{\tau }_2}}(C_{\mathrm{se}}^2+\gamma V_{\mathrm{Ae}}^2)},\end{array}$$(17)

where $C_{\mathrm{Te}}=\sqrt{\frac{C_{\mathrm{se}}^2{V}_{\mathrm{Ae}}^2}{C_{\mathrm{se}}^2+V_{\mathrm{Ae}}^2}}$ represents the exterior tube speed. It is worth noting that the radial wavenumbers in Eqs. (15, 17), as the arguments of the Bessel functions, are complex expressions. therefore, the real part of both Bessel functions Y_m and K_m meet the condition for finite exterior solutions as the distance from the flux tube boundary increases. However, we take the Bessel function K_m as the solution for the external region to remain consistent with the adiabatic limit. This limit corresponds to the characteristic timescales τ₁ and τ₂ approaching infinity, thus recovering the exterior solution in Eq. (7) of Edwin & Roberts (1983). The dispersion relation is then obtained by examining the continuity of the radial velocity component and the total pressure at the tube boundary (r = R), expressed as the following function:

$$\begin{array}{c}\hfill D_m(\omega )=(k^2{V}_{\mathrm{A}}^2-{\omega }^2)-\frac{{\rho }_{\mathrm{e}}}{{\rho }_0}\frac{{\kappa }_{\mathrm{in}}}{{\kappa }_{\mathrm{out}}}(k^2{V}_{\mathrm{Ae}}^2-{\omega }^2)\frac{J_m^{\prime }({\kappa }_{\mathrm{in}}R)K_m({\kappa }_{\mathrm{out}}R)}{J_m({\kappa }_{\mathrm{in}}R)K_m^{\prime }({\kappa }_{\mathrm{out}}R)}=0.\end{array}$$(18)

This equation recovers Eq. (8b) of Edwin & Roberts (1983), in the adiabatic limit. As evident in Eq. (18), the distinction between surface and body modes within the tube is no longer clear, as the radial wavenumbers κ_in and κ_out are complex numbers. This is leading to an intertwined state of surface and body waves, where the solution to the Bessel equation inside the tube (J_m(κ_inr)) becomes a combination of both modes.

Building on the findings of Agapova et al. (2022) concerning the dynamics of fast and slow MA waves in plasma slabs affected by thermal misbalance, and also supported by the results presented in this study, it becomes obvious that the responses to the dispersion relation are primarily dominated by body modes rather than surface modes.

3. Long wavelength limit

For the case $\frac{r}{\lambda }\ll 1$ or, in other words, for small values of the longitudinal wavenumber k, we can implement approximations for Bessel functions with the small arguments as discussed in Section 2. This allows us to derive an analytical expression that describes the effect of thermal misbalance on the MA modes within the flux tube. In this section, mathematical expressions are calculated for the axisymmetric (sausage) and the first azimuthal harmonic (kink) modes.

It is important to note that all the numerical and analytical solutions presented in this section are carried out under typical coronal conditions, where the interior Alfvén speed (V_Ai) is assumed to be 1 Mm/s, the exterior Alfvén speed (V_Ae) is 1.5 Mm/s, the interior sound speed (C_si) is 0.15 Mm/s, and the exterior sound speed (C_se) is 0.13 Mm/s. These parameter values refer to the coronal conditions used in all figures throughout this work.

3.1. Fast sausage modes

In the case of zero azimuthal wavenumber (m = 0), the dispersion relation (Eq. (18)) takes the following form:

$$\begin{array}{c}\hfill D_0(\omega )=(k^2{V}_{\mathrm{A}}^2-{\omega }^2)-\frac{{\rho }_{\mathrm{e}}}{{\rho }_0}\frac{{\kappa }_{\mathrm{in}}}{{\kappa }_{\mathrm{out}}}(k^2{V}_{\mathrm{Ae}}^2-{\omega }^2)\frac{J_1({\kappa }_{\mathrm{in}}R)K_0({\kappa }_{\mathrm{out}}R)}{J_0({\kappa }_{\mathrm{in}}R)K_1({\kappa }_{\mathrm{out}}R)}=0.\end{array}$$(19)

Figure 2 presents a numerical solution computed using MATLAB to find the roots of the dispersion relation given by Eq. (19). The plot depicts the fundamental MHD modes with phase speeds faster than the internal Alfvén speed, which are classified as fast MHD modes. As depicted in Fig. 2a (solid line), the fundamental fast sausage wave exhibits a cut-off in wavenumbers at high frequencies with values exceeding the external Alfvén frequency. This cut-off behaviour is an inevitable characteristic of all fast sausage and fast non-fundamental kink modes in the coronal region (Edwin & Roberts 1983). Consequently, for fast MHD mode responses with frequencies or phase speeds surpassing those of the external Alfvén’s, the solutions to the Bessel equation in the exterior region manifest in cylindrical Hankel function forms rather than Bessel K-function forms. This observation indicates that the waves would exhibit leaky forms of propagation, and the dispersion relation Eq. (18) loses its validity in describing fast MHD modes with phase speeds exceeding the external Alfvén speed (for an extensive discussion, refer to Farahani et al. 2014). Furthermore, in light of a recent study into the cut-off wavenumber observed in the context of sausage oscillations within a slab configuration with varying external density, it is advisable to consult Wang et al. (2023) as a supplementary reference.

Fig. 2. In the coronal condition (VAe = 1.5, VAi = 1, Csi = 0.15, Cse = 0.13 Mm/s) and with $\overline{{\tau }_1}=6.8$, $\overline{{\tau }_2}=2.4$ (corresponding to approximate values of 34 and 12 minutes for τ1 and τ2, respectively, at a temperature of about 5 MK), the study considers a wave with a period of 5 min: (a) numerical solution of the dispersion relation from Eq. (19) for frequencies (phase speed) normalized by the internal acoustic frequency (acoustic speed) (ωr/ωsi ≡ vph/Csi) of the fast sausage fundamental mode without (solid blue line) and with (solid black line) consideration of the thermal misbalance effect, as well as the numerical solution of the dispersion relation from Eq. (27) for the fast kink fundamental mode (m = 1) without (dashed blue line) and with (dashed black line) consideration of the thermal misbalance effect versus the normalized longitudinal wavenumber (kR). (b) Damping rate to frequency ratio (−γ/ωr) of fundamental states of the fast sausage (solid line) and kink (dashed line) modes with consideration of the thermal misbalance effect versus the normalized longitudinal wavenumber (kR).

Hence, in this study, it is crucial to confine our analysis to responses where the phase speed does not surpass the external Alfvén speed to avoid compromising the validity of the dispersion relation, Eq. (18). Subsequent research on the effect of thermal misbalance on leaky wave solutions could be pursued in future research.

In the regime of long-wavelength sausage modes where the argument of the Bessel functions is approaching to small values, the dispersion relation Eq. (19) takes on an approximate form as

$$\begin{array}{c}\hfill \frac{{\rho }_0(k^2{V}_{\mathrm{A}}^2-{\omega }^2)}{{\rho }_{\mathrm{e}}(k^2{V}_{\mathrm{Ae}}^2-{\omega }^2)}\approx \frac{{\kappa }_{\mathrm{in}}^2{R}^2}{2}[-ln\left(\frac{{\kappa }_{\mathrm{out}}R}{2}\right)-{\gamma }_{\mathrm{e}}]\approx -\frac{{\kappa }_{\mathrm{in}}^2{R}^2}{2}ln(kR).\end{array}$$(20)

Here we have taken the long-wavelength limit expansion of the Bessel K-function, where |ln(kR)|≫1. As evident from the simplified dispersion relation (Eq. (20)), as k approaches zero, κ_in² diminishes much faster than ln(kR) approaches infinity. Consequently, the dispersion relation only yields a solution at the frequency equal to the internal Alfvén frequency (ω_Ai), indicating that in the limit of small longitudinal wavenumbers, the dispersion relation has no solutions for frequencies exceeding the internal Alfvén frequency associated with fast sausage modes. Consequently, the dispersion relation, Eq. (19), yields a solution at the frequency corresponding to the internal Alfvén frequency (ω_Ai), indicating that, in the limit of small longitudinal wavenumbers, the dispersion relation lacks solutions for frequencies exceeding the internal Alfvén frequency associated with long wavelength fast sausage modes. Based on analytical results and referring to the numerical solution (Fig. 2, solid lines), the cut-off occurring in the range of small longitudinal wavenumbers and frequencies between internal and external Alfvén frequencies is a phenomenon more pronounced than the ability to find solutions for frequencies even under the influence of thermal misbalance on the fast asymmetric MHD modes, with the cut-off effect remaining dominant. Therefore, in this limit, it is impossible to find an analytical expression describing long wavelength fast sausage modes.

Moreover, as evident from the numerical solution in Fig. 2a (solid line), thermal misbalance exhibits no impact on the position of cut-off wavenumber and minimal impact on the phase speed of fast sausage modes. In Fig. 2b (solid line), the damping rate closely approaches that of the fast kink mode (dashed line) after the cut-off wavenumber, with damping rate values being negligible, thus insignificantly influencing the overall behaviour of fast sausage modes, contrary to the slow sausage modes studied in the next subsection.

3.2. Slow sausage modes

Slow MHD modes arise from the analysis of wave behaviour approaching either the subsonic or sub-Alfvénic frequencies around the internal cusp frequency, denoted as ω_C. When k²C_T² − ω² = ϵ, and ϵ is a very small value, the interior radial wavenumber Eq. (15) displays a singularity implying a resonance condition in the absence of the thermal misbalance effect. At this point, the cusp frequency is equal to the internal tube frequency (ω_Ti = kC_Ti), see Fig. 3a blue solid (slow sausage) and blue dashed (slow kink) lines. However, the presence of thermal misbalance introduces a slight shift in the phase speed values deviated from the values obtained without considering the influence of thermal misbalance. This results in a complex-valued cusp frequency that differs from the tube frequency in its real part, while also introducing an imaginary part directly associated with one of the characteristic timescales ($\overline{{\tau }_2}$) connected to the thermal misbalance effect (e.g. Fig. 3a black solid and dashed lines). Subsequently, the cusp frequency undergoes modification and can be obtained by examining the denominator of the interior radial wavenumber given by Eq. (15) for frequency values that cause it to approach zero. In this study, the modified cusp frequency is denoted by $\omega {\prime }_{\mathrm{C}}$ and expressed as

Fig. 3. In the coronal condition and $\overline{{\tau }_1}=6.8$, $\overline{{\tau }_2}=2.4$: (a) numerical solution of the dispersion relation Eq. (18), computed by MATLAB, for the frequency (phase speed) normalized by the internal acoustic frequency (acoustic speed) (ℜ𝔢(ω)/ωsi ≡ vph/Csi) of the fundamental slow sausage mode without consideration of the thermal misbalance effect (unaffected) (solid blue line) and with consideration of the thermal misbalance effect (affected) (solid black line), as well as the fundamental slow kink mode unaffected (dashed blue line) and affected (dashed black line) by the thermal misbalance versus the normalized longitudinal wavenumber (kR). (b) Damping rate to the frequency ratio (−γ/ωr) of fundamental states of the affected slow sausage (solid line) and kink (dashed line) modes by the thermal misbalance versus the normalized longitudinal wavenumber (kR). (c) Damping time to the wave period ratio (τd/P) of the affected slow sausage (solid line) and kink (dashed line) modes by the thermal misbalance versus the normalized longitudinal wavenumber (kR). For comparison, the analytical results obtained by Eqs. (25, 32) are shown by the black and blue solid curves in panels d, e, f.

$$\begin{array}{c}\hfill {{\omega }_{\mathrm{C}}^{\prime }}^2={\omega }_{\mathrm{T}}^2{[1+\frac{i}{2\pi \gamma \overline{{\tau }_2}}\frac{C_{\mathrm{s}}^2+\gamma V_{\mathrm{A}}^2}{C_{\mathrm{s}}^2+V_{\mathrm{A}}^2}]}^{-1},\end{array}$$(21)

As an illustration, when considering a typical condition associated with coronal active regions and normalizing all system speeds by the internal sound speed, the modified cusp speed ($V{\prime }_{\mathrm{C}}$) takes on a complex form, approximately 0.987339 − 0.0324i for timescales specified by values of $\overline{{\tau }_1}=6.8,\overline{{\tau }_2}=2.4$, whereas the interior tube speed (C_T) under the same conditions is approximately 0.988936. Consequently, the recalibration of the interior radial wavenumber (κ_in) with the utilization of this modified cusp frequency is denoted by $\kappa {\prime }_{\mathrm{in}}$ and rewritten as

$$\begin{array}{c}\hfill {{\kappa }_{\mathrm{in}}^{\prime }}^2=-\frac{({\omega }_{\mathrm{s}}^2-{\omega }^2)({\omega }_{\mathrm{A}}^2-{\omega }^2)-\frac{i{\omega }^2}{2\pi \overline{{\tau }_2}}(\frac{\overline{{\tau }_2}}{\overline{{\tau }_1}}{\omega }_{\mathrm{s}}^2-{\omega }^2)}{({{\omega }_{\mathrm{C}}^{\prime }}^2-{\omega }^2)[C_{\mathrm{s}}^2+V_{\mathrm{A}}^2+\frac{i}{2\pi \gamma \overline{{\tau }_2}}(C_{\mathrm{s}}^2+\gamma V_{\mathrm{A}}^2)]},\end{array}$$(22)

where ω_s² ≡ k²C_s² and ω_A² ≡ k²C_A².

To obtain the solution of the dispersion relation Eq. (19) concerning frequencies in the vicinity of the modified cusp frequency given in Eq. (21) under the long-wavelength limit (kR → 0), we can employ a frequency expansion with the form of ${\omega }^2\approx {\omega {\prime }_{\mathrm{C}}}^2[1+{\delta }_0{k}^2{R}^{2}ln(kR)]$ as derived by Moreels et al. (2015b). This expansion is applied on the case of m = 0, where δ₀ represents a small deviation determined by the solution of the dispersion relation around the modified cusp frequency. Therefore, utilizing the frequency expansion and expressing in term of phase speeds, Eq. (22) transforms to

$$\begin{array}{cc}& {{\kappa }_{0\mathrm{in}}^{\prime }}^2{R}^2=\frac{(C_{\mathrm{s}}^2-{V_{\mathrm{C}}^{\prime }}^2)(V_{\mathrm{A}}^2-{V_{\mathrm{C}}^{\prime }}^2)-\frac{i{V_{\mathrm{C}}^{\prime }}^2}{2\pi \overline{{\tau }_2}}(\frac{\overline{{\tau }_2}}{\overline{{\tau }_1}}{C}_{\mathrm{s}}^2-{V_{\mathrm{C}}^{\prime }}^2)}{{\delta }_0{V_{\mathrm{C}}^{\prime }}^{2}ln(kR)[C_{\mathrm{s}}^2+V_{\mathrm{A}}^2+\frac{i}{2\pi \gamma \overline{{\tau }_2}}(C_{\mathrm{s}}^2+\gamma V_{\mathrm{A}}^2)]},\hfill \\ & {V_{\mathrm{C}}^{\prime }}^2=\frac{{{\omega }_{\mathrm{C}}^{\prime }}^2}{k^2}·\hfill \end{array}$$(23)

where $V{\prime }_{\mathrm{C}}$ is modified cusp speed. By substituting this expression into the dispersion relation Eq. (19) and applying standard approximations for Bessel functions, as outlined in deriving Eq. (20), we obtain

$$\begin{array}{c}\hfill D_0(V_{\mathrm{C}}^{\prime })=(V_{\mathrm{A}}^2-{V_{\mathrm{C}}^{\prime }}^2)+\frac{{\rho }_{\mathrm{e}}}{{\rho }_0}(V_{\mathrm{Ae}}^2-{V_{\mathrm{C}}^{\prime }}^2)\frac{{{\kappa }_{0\mathrm{in}}^{\prime }}^2{R}^{2}ln(kR)}{2}=0,\end{array}$$(24)

By utilizing this expression and employing the frequency expansion used for slow sausage modes, the damping rate is expressed as

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{damp}}=\mathfrak{Im}\left[{\omega }_{\mathrm{C}}^{\prime }(1+\frac{1}{2}{\delta }_0{k}^2{R}^{2}ln(kR))\right],\end{array}$$(25)

Here, δ₀ is obtained by solving Eq. (24) and takes on a mathematical form of

$$\begin{array}{cc}\hfill {\delta }_0& =-\frac{{\rho }_{\mathrm{e}}}{2{\rho }_0{V_{\mathrm{C}}^{\prime }}^2}\left(\frac{V_{\mathrm{Ae}}^2-{V_{\mathrm{C}}^{\prime }}^2}{V_{\mathrm{A}}^2-{V_{\mathrm{C}}^{\prime }}^2}\right)\hfill \\ \hfill & \times \left[\frac{(C_{\mathrm{s}}^2-{V_{\mathrm{C}}^{\prime }}^2)(V_{\mathrm{A}}^2-{V_{\mathrm{C}}^{\prime }}^2)-\frac{i{V_{\mathrm{C}}^{\prime }}^2}{2\pi \overline{{\tau }_2}}(\frac{\overline{{\tau }_2}}{\overline{{\tau }_1}}{C}_{\mathrm{s}}^2-{V_{\mathrm{C}}^{\prime }}^2)}{C_{\mathrm{s}}^2+V_{\mathrm{A}}^2+\frac{i}{2\pi \gamma \overline{{\tau }_2}}(C_{\mathrm{s}}^2+\gamma V_{\mathrm{A}}^2)}\right].\hfill \end{array}$$(26)

In Figs. 3d–f, numerical solutions (dashed black lines) obtained by numerically solving Eq. (18) for the axisymmetric MHD modes with m = 0 are compared with analytical solutions (solid black lines) derived from the frequency expansion ${\omega }^2\approx {\omega {\prime }_{\mathrm{C}}}^2[1+{\delta }_0{k}^2{R}^{2}ln(kR)]$, where the parameters are obtained from Eqs. (21, 26), and Eq. (25) for the fundamental state of slow sausage modes. The comparison reveals a high degree of consistency between the numerical and analytical solutions, particularly evident in the long-wavelength limit (small longitudinal wavenumbers k), affirming the validity of our analytical expressions.

3.3. Fast kink modes

Higher azimuthal harmonics (m > 0) of the dispersion relation in Eq. (18) can take the form of

$$\begin{array}{cc}& D_m(\omega )=(k^2{V}_{\mathrm{A}}^2-{\omega }^2)\hfill \\ \hfill & -\frac{{\rho }_{\mathrm{e}}}{{\rho }_0}\frac{{\kappa }_{\mathrm{in}}}{{\kappa }_{\mathrm{out}}}(k^2{V}_{\mathrm{Ae}}^2-{\omega }^2)\frac{(\frac{m}{{\kappa }_{\mathrm{in}}R}{J}_m({\kappa }_{\mathrm{in}}R)-J_{m+1}({\kappa }_{\mathrm{in}}R))K_m({\kappa }_{\mathrm{out}}R)}{J_m({\kappa }_{\mathrm{in}}R)(\frac{m}{{\kappa }_{\mathrm{out}}R}{K}_m({\kappa }_{\mathrm{out}}R)-K_{m+1}({\kappa }_{\mathrm{out}}R))}\hfill \\ \hfill & =0,m>0.\hfill \end{array}$$(27)

In the limit of long wavelength, the Bessel approximations for small values of the longitudinal wavenumber (k) allow us to rewrite Eq. (27) to a reduced form of

$$\begin{array}{c}\hfill \frac{{\rho }_0(k^2{V}_{\mathrm{A}}^2-{\omega }^2)}{{\rho }_{\mathrm{e}}(k^2{V}_{\mathrm{Ae}}^2-{\omega }^2)}\approx \frac{{\kappa }_{\mathrm{in}}^2{R}^2}{2m(m+1)}-1,m>0.\end{array}$$(28)

Subsequently, in the regime of small longitudinal wavenumbers where the interior radial wavenumber Eq. (15) approaches zero, for non-zero azimuthal harmonic numbers (including m = 1 to describe kink modes), the following expression can be obtained from Eq. (28):

$$\begin{array}{c}\hfill \frac{{\rho }_0(k^2{V}_{\mathrm{A}}^2-{\omega }^2)}{{\rho }_{\mathrm{e}}(k^2{V}_{\mathrm{Ae}}^2-{\omega }^2)}\approx -1,\end{array}$$(29)

which is only solved by the kink frequency of

$$\begin{array}{c}\hfill {\omega }_k^2=\frac{{\rho }_0{V}_{\mathrm{A}}^2+{\rho }_{\mathrm{e}}{V}_{\mathrm{Ae}}^2}{{\rho }_0+{\rho }_{\mathrm{e}}}{k}^2.\end{array}$$

It can be readily inferred from Eq. (28) that thermal misbalance has no significant effect on fast either kink or higher azimuthal harmonic waves (m > 1) in the long wavelength limit. The numerical analysis is also shown in Fig. 2 (dashed lines) for fast kink modes that thermal misbalance has no impact on the fundamental fast kink frequency (phase speed). However, a very small decrement is observable for the damping rate in Fig. 2b (dashed black line) for small values of longitudinal wavenumbers which may be inferred as suppressed damping state of kink oscillations (Agapova et al. 2022), and very interesting for further research. In addition, a slight increase is observed for greater values of longitudinal wavenumbers while the damping rate values are very small and negligible. The conclusion of no significant impact is totally different for slow kink modes.

3.4. Slow kink modes

In the case of slow kink modes, the criterion for finding an analytical expression is similar to what was done in Sect. 3.2 for slow sausage modes. The analytical expression for the damping rate is derived by employing the modified cusp frequency ($\omega {\prime }_{\mathrm{C}}$) defined in Eq. (21) and the rewritten form of the interior radial wavenumber ($\kappa {\prime }_{\mathrm{in}}$) in Eq. (22). By setting the dispersion relation Eq. (27) for the first azimuthal wavenumber (m = 1), and employing the Bessel function approximations under the assumption of kR ≪ 1, an approximate form of the dispersion relation Eq. (27) describing the long wavelength slow kink modes is obtained as follows

$$\begin{array}{c}\hfill D_1(V_{\mathrm{ph}})=(V_{\mathrm{A}}^2-V_{\mathrm{ph}}^2)+\frac{{\rho }_{\mathrm{e}}}{{\rho }_0}(V_{\mathrm{Ae}}^2-V_{\mathrm{ph}}^2)[1-\frac{{{\kappa }_{\mathrm{in}}^{\prime }}^2{R}^2}{4}]=0.\end{array}$$(30)

where $V_{\mathrm{ph}}=\frac{\omega }{k}$. Once again, similar to the slow sausage modes, to analytically study the kink wave around the cusp frequency in the long-wavelength limit, we need to employ a frequency expansion. However, this time, it is specified for the case of m = 1 with the form of ${\omega }^2\approx {\omega {\prime }_{\mathrm{C}}}^2(1+{\delta }_1{k}^2{R}^2)$ (Edwin & Roberts 1983), where the small deviation δ₁ is obtained by solving the dispersion relation Eq. (30) for the modified cusp speed ($V{\prime }_{\mathrm{C}}$).

By utilizing the frequency expansion, the radial wavenumber $\kappa {\prime }_{\mathrm{in}}$ given by Eq. (22) for the modified cusp speed has the following mathematical expression

$$\begin{array}{c}\hfill {{\kappa }_{1\mathrm{in}}^{\prime }}^2{R}^2=\frac{(C_{\mathrm{s}}^2-{V_{\mathrm{C}}^{\prime }}^2)(V_{\mathrm{A}}^2-{V_{\mathrm{C}}^{\prime }}^2)-\frac{i{V_{\mathrm{C}}^{\prime }}^2}{2\pi \overline{{\tau }_2}}(\frac{\overline{{\tau }_2}}{\overline{{\tau }_1}}{C}_{\mathrm{s}}^2-{V_{\mathrm{C}}^{\prime }}^2)}{{\delta }_1{V_{\mathrm{C}}^{\prime }}^2(C_{\mathrm{s}}^2+V_{\mathrm{A}}^2+\frac{i}{2\pi \gamma \overline{{\tau }_2}}(C_{\mathrm{s}}^2+\gamma V_{\mathrm{A}}^2))},\end{array}$$(31)

By replacing this expression into the dispersion relation Eq. (30) instead of ${\kappa {\prime }_{\mathrm{in}}}^2{R}^2$, for the phase speed equal to $V{\prime }_{\mathrm{C}}$, and solving it for the deviation of δ₁, an analytical expression for the damping rate is obtained as follows

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{damp}}=\mathfrak{Im}\left[{\omega }_{\mathrm{C}}^{\prime }(1+\frac{1}{2}{\delta }_1{k}^2{R}^2)\right],\end{array}$$(32)

where includes the deviation of δ₁ defining with the following expression

$$\begin{array}{cc}\hfill {\delta }_1& =\left(\frac{(C_{\mathrm{s}}^2-{V_{\mathrm{C}}^{\prime }}^2)(V_{\mathrm{A}}^2-{V_{\mathrm{C}}^{\prime }}^2)-\frac{i{V_{\mathrm{C}}^{\prime }}^2}{2\pi \overline{{\tau }_2}}(\frac{\overline{{\tau }_2}}{\overline{{\tau }_1}}{C}_{\mathrm{s}}^2-{V_{\mathrm{C}}^{\prime }}^2)}{4{V_{\mathrm{C}}^{\prime }}^2[C_{\mathrm{s}}^2+V_{\mathrm{A}}^2+\frac{i}{2\pi \gamma \overline{{\tau }_2}}(C_{\mathrm{s}}^2+\gamma V_{\mathrm{A}}^2)]}\right)\hfill \\ \hfill & \times {[1+\frac{{\rho }_0}{{\rho }_{\mathrm{e}}}\frac{V_{\mathrm{A}}^2-{V_{\mathrm{C}}^{\prime }}^2}{V_{\mathrm{Ae}}^2-{V_{\mathrm{C}}^{\prime }}^2}]}^{-1}.\hfill \end{array}$$(33)

The comparison in Figs. 3d–f reveals a high degree of consistency between the numerical (dashed blue lines) and analytical (solid blue lines) solutions for the fundamental state of slow kink modes, arisen from Eq. (27) (for transverse harmonic m = 1) and Eq. (33) respectively, particularly evident in the long wavelength limit (small longitudinal wavenumbers k), and slight deviation observed for larger values of longitudinal wavenumbers, affirming the validity of our analytical equation (Eq. (32)).

Finally, in Fig. 4, we numerically analyse the behaviour of slow sausage and kink waves and their associated damping characteristics across varying periods of 3, 5, and 8 minutes. It is observed that longer periods correspond to smaller values of $\overline{{\tau }_1}$ and $\overline{{\tau }_2}$, leading to enhanced damping rates. The analytical expressions in Eqs. (25, 32) also demonstrate this conclusion.

Fig. 4. In the coronal condition, the figure displays: (a), (b), and (c) the normalized numerical solutions of the dispersion relation Eq. (19) by the internal acoustic speed (Cs) for sausage modes (m = 0) with different periods of 8 min ($\overline{{\tau }_1}=4.25$, $\overline{{\tau }_2}=1.5$) (solid line), 5 min ($\overline{{\tau }_1}=6.8$, $\overline{{\tau }_2}=2.4$) (dashed line), and 3 min ($\overline{{\tau }_1}=11.34$, $\overline{{\tau }_2}=4$) (dotted line); and (d), (e), and (f): the normalized numerical solutions of the dispersion relation Eq. (27) by the internal acoustic speed (Cs) for kink modes (m = 1) with different periods of 8 min ($\overline{{\tau }_1}=4.25$, $\overline{{\tau }_2}=1.5$) (solid line), 5 min ($\overline{{\tau }_1}=6.8$, $\overline{{\tau }_2}=2.4$) (dashed line), and 3 min ($\overline{{\tau }_1}=11.34$, $\overline{{\tau }_2}=4$) (dotted line).

4. Conclusions

In this study, we investigated the propagation of MHD waves dispersed in a cylindrical waveguide and affected by non-adiabatic factors, with thermal misbalance at the centre of focus. We characterized the impact of thermal misbalance resulting from the deviation from thermal equilibrium in the coronal region between heating and radiative cooling mechanisms, using two distinct dimensionless characteristic times $\overline{{\tau }_1}$ and $\overline{{\tau }_2}$. These times play a crucial role in either attenuating or amplifying MHD wave amplitudes within a cylindrical magnetic waveguide in a coronal uniform plasma. The characteristic times are directly linked to the heating-cooling function Q(ρ, T), the sole non-adiabatic term considered in this analysis. We deliberately selected positive values for these timescales, ensuring that $\overline{{\tau }_1}-\overline{{\tau }_2}>0$, aligning with constraints consistent with the stronger damping (enhanced damping) regimes identified in Kolotkov et al. (2020). The ensuing results are presented as follows:

1. The dispersion relation governing MHD modes in cylindrical waveguides immersed in a uniform finite magnetic field and subjected to the discrepancy between heating and radiative cooling mechanisms was obtained. The impact of thermal misbalance, as elucidated by two characteristic times linked to the derivative of the heating-cooling function with respect to density perturbation in constant temperature and with respect to temperature perturbation in constant density, initially introduced by Kolotkov et al. (2019), were normalized by the MHD wave period. These normalized characteristic timescales are employed to encapsulate the influence of thermal misbalance properties on the radial wavenumber.

2. Both the analytical and numerical solutions in the limit of fast MHD modes were obtained to show that thermal misbalance exhibits only a marginal impact on these modes which stems from the consideration that, in high-frequency regimes (ω ≫ 1/min{|τ₁|,|τ₂|}) (Nakariakov et al. 2017), the influence of thermal misbalance on the damping rate remains minimal and negligible, in consistency with the assumption established firstly by Zavershinskii et al. (2019). However, the slight decrement observed in the small longitudinal wavenumbers for the fundamental fast kink wave, Fig. 2b (dashed line), may signify a suppressed damping mode, which merits further investigation.

3. In the case of slow MHD modes within a magnetic tube, both analytical and numerical results were obtained. Analytical expressions in Eqs. (25, 32) governing the damping rates in the long wavelength limit were derived for slow sausage and kink modes respectively, demonstrating exemplary agreement with the corresponding numerical results presented in Fig. 3. The equations and figures collectively reveal that thermal misbalance significantly impacts slow modes in comparison to fast modes. Furthermore, our analysis indicates that the influence of thermal misbalance on slow MHD wave propagation in the coronal conditions is more pronounced in asymmetric modes than in axisymmetric ones. This result opposes with the expected behaviour, as MA waves characterized by higher compressibility exhibit greater sensitivity to thermal misbalance. Figures 3b and e illustrate that the damping rates for slow kink modes are higher than those for slow sausage modes with the same period. This finding suggests that thermal misbalance has a more substantial impact on slow asymmetric MA modes than their axisymmetric counterparts. In fast modes, as depicted in Fig. 2, a slight inverse relationship is evident. As a result, we can conclude that the MA waves owning the faster phase speeds are less affected than the slower ones by thermal misbalance.

4. In the investigation of either high-frequency (small period) waves or elevated temperatures, an examination of the surface plot Fig. 1 reveals that both normalized characteristic timescales, denoted as $\overline{{\tau }_1}$ and $\overline{{\tau }_2}$, tend towards their maximum values. This observation suggests that, under the specified conditions which are subject to the heating model and consequently the values considered for the parameters, a and b, in the heating function Eq. (7), the influence of thermal misbalance on the dynamical behaviour of cylindrical MHD modes is minimized. This is because the characteristic timescales, in their respective maximum values, minimize the imaginary part of the radial wavenumbers, as evident from Eqs. (15 and 17). On the contrary, the diminution of the characteristic timescales exerts a pronounced influence, leading to an enhanced damping rate for MHD waves, particularly for slow modes. Moreover, as $\overline{{\tau }_1}$ and $\overline{{\tau }_2}$ exhibit a correlation with the temperature of the active plasma region, as observable in Fig. 1, the damping rates Eqs. (25, 32) consequently demonstrate a correlated behaviour with temperature. In such a way, the damping rate decreases with rising temperature and increases as the temperature decreases to the 1−5 MK range which is evident in Fig. 5 for transverse slow oscillations as well. Importantly, our analysis underscores the significance of $\overline{{\tau }_2}$ over $\overline{{\tau }_1}$. Especially, the heating-cooling function’s rate with respect to temperature, encapsulated in $\overline{{\tau }_2}$, plays a more prominent role in shaping wave dynamics. This is particularly evident as $\overline{{\tau }_2}$ is the exclusive timescale in the denominator of the fourth power of wave frequency in the imaginary parts of the radial wavenumbers (Eqs. (15) and (17)), and it is the only characteristic timescale appearing in the cusp frequency Eq. (21). Consequently, smaller values of $\overline{{\tau }_2}$ signify a more substantial impact of thermal misbalance.

5. As mentioned, the analytical damping rates of slow MHD waves were determined in Eqs. (25) and (32), and the obtained values were found to be comparable to their respective wave periods, aligning well with observed characteristics. To elucidate this relationship, we employed the analytical expression derived from Eq. (32), which describes the damping rate of slow transverse MA modes in a magnetic tube. The results, as illustrated in Fig. 5, demonstrate that the fraction of damping time per the period of the wave for SUMER oscillations falls within the range of approximately 1.55 to 2.85 times the wave period at 6.3 MK and approximately 2.9 to 5.5 times the wave period at 8.9 MK, for a range of normalized longitudinal wavenumbers (kR) between 0 and 2. This finding is consistent with the observed rapid damping of slow modes, especially in the limit of long-wavelength waves, substantiating the validity of the analytical model in capturing the damping behaviour of the studied MHD waves.

6. Adjusting the thermal misbalance characteristic timescales, $\overline{{\tau }_1}$ and $\overline{{\tau }_2}$, to match the observed ratio of rapidly damping or amplifying times per the period of the oscillations in various layers of the solar atmosphere offers seismological insights on plasma characteristics or mechanisms, for example, plasma density or the heating function. It is also examinable alongside other methods, especially in the context of identifying the heating function. As an instance, in this study, we demonstrated that the selected values for these timescales, taken from the work of Kolotkov et al. (2019), show compatibility with the damping times observed for the considered SUMER oscillations, as shown in Fig. 5. This study propose that gathering the values of these timescales consistent with each oscillation provides a set of statistical data to specify the heating function form with the more accurate values of the heating function exponents of a and b in Eq. (7).

Fig. 5. In the coronal condition, for temperatures associated with SUMER oscillations with periods of 15.7 minutes at 6.3 MK ($\overline{{\tau }_1}\approx 2.36$, $\overline{{\tau }_2}\approx 0.76$) (dashed line) and 13.3 minutes at 8.9 MK ($\overline{{\tau }_1}\approx 4.89$, $\overline{{\tau }_2}\approx 1.43$) (solid line), as calculated using the analytical expression in Eq. (32), the figure depicts the ratio of the damping time to the wave period (τd/P) of the affected slow transverse MA modes by the thermal misbalance versus kR.

In conclusion, our current study has provided robust analytical and numerical evaluations of MHD modes in a cylindrical configuration influenced by thermal misbalance. Building upon this foundation, the next crucial steps involve extending our theoretical framework to incorporate other non-adiabatic terms, notably thermal conduction and scalar resistivity. Given that thermal conduction always acts as a damping for slow MA modes, as demonstrated by Duckenfield et al. (2021), incorporating it alongside thermal misbalance promises more precise estimates for damping times, aligning closely with observed values.

Furthermore, the methodology and results presented in this article can be replicated for alternative sets of thermal misbalance characteristic times, shedding light on non-decaying, amplifying, and suppressed damping oscillations, aligning with the states identified by Kolotkov et al. (2020). It is also valuable to consider other physical parameters participating in the heating function such as magnetic field (Duckenfield et al. 2021), and study their effects on MHD modes inside a cylindrical waveguide. These investigations can contribute to refining heating function models governing magnetic structures within the solar atmosphere.

An intriguing avenue for future research involves introducing the concept of a resonant layer at the boundary between the inside and outside of solar tubes. Moreover, a differential operator, similar to those defined for the geometrical constraint and thermal misbalance effect in Section 2, can be introduced to describe the impact of thermal conduction as an additional cooling process alongside radiative cooling in the dispersion relation Eq. (8) (Field 1965), and by incorporating the characteristic timescale (τ_cond) related to thermal conduction (Kolotkov et al. 2019) in the solar atmosphere, we anticipate more accurate predictions of damping times for slow modes observed in this region.

Expanding our focus to include the non-linear evolution of MHD waves in non-uniform plasmas, especially due to gravitational stratification, promises exciting insights. Exploring the impact of thermal misbalance on non-linear phenomena, such as shock waves, becomes particularly relevant in understanding energy release processes in the solar atmosphere, unravelling the heating mechanism. Among these, investigating the compressible aspects of non-linear evolution of Alfvén waves, and their response to thermal misbalance, especially when propagating through inhomogeneous plasma and forming shock waves, or how thermal misbalance influences force equilibrium, specifically its effects on the ponderomotive force, presents opportunities for further investigations.

Finally, it can be the best point to run simulations that play a pivotal role in obtaining a more tangible understanding of how non-adiabatic terms affect MHD waves within magnetic structures in the solar atmosphere. The implications derived from the findings presented in this article and the proposed follow-up investigations significantly enhance our understanding of the intricate dynamics that govern the solar corona, particularly concerning its thermal characteristics.

Acknowledgments

CHIANTI is a collaborative project involving George Mason University, the University of Michigan (USA) and the University of Cambridge (UK). TVD received financial support from the Flemish Government under the long-term structural Methusalem funding program, project SOUL: Stellar evolution in full glory, grant METH/24/012 at KU Leuven, the DynaSun project (number 101131534 of HORIZON-MSCA-2022-SE-01), and also a Senior Research Project (G088021N) of the FWO Vlaanderen. DYK acknowledges funding from the STFC consolidated grant ST/X000915/1 and Latvian Council of Science Project No. lzp2022/1-0017.

References

Appendix A: Derivation of the differential operators in Eq. (9)

Considering the system at an initial magnetohydrodynamical equilibrium and by applying small amplitude perturbations to the governing equations (Eqs. (1-6)), their linearized forms are derived, respectively, as follows

$$\begin{array}{cc}& {\rho }_0\frac{\partial {\mathbf{v}}_1}{\partial t}=-\mathrm{\nabla }{p}_1-\frac{1}{4\pi }[\mathrm{\nabla }({\mathit{B}}_0·{\mathit{B}}_1)-({\mathit{B}}_0·\mathrm{\nabla }){\mathit{B}}_1],\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& \frac{\partial {\rho }_1}{\partial t}+{\rho }_0\mathrm{\nabla }·{\mathbf{v}}_1=0,\hfill \end{array}$$(A.2)

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

$$\begin{array}{cc}& p_1=\frac{k_B}{m}({\rho }_0{T}_1+T_0{\rho }_1),\hfill \end{array}$$(A.4)

$$\begin{array}{cc}& \frac{1}{\gamma -1}(\frac{\partial p_1}{\partial t}-\frac{\gamma p_0}{{\rho }_0}\frac{\partial {\rho }_1}{\partial t})=-{\rho }_0[{\rho }_1{Q}_T+T_1{Q}_{\rho }],\hfill \end{array}$$(A.5)

where, the first-order perturbed parameters are denoted by the subscript 1, while the equilibrium quantities are denoted by the subscript 0. For the definitions of the quantities Q_ρ and Q_T, refer to Eq. (10). By substituting the first order magnetic field perturbation (B₁) in Eq. (A.1), using Eq. (A.3), we obtain

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\mathrm{\nabla }}^2{p}_1={\rho }_0(V_A^2{\mathrm{\nabla }}^2(\mathrm{\nabla }·{\mathbf{v}}_1)-V_A^2\frac{{\partial }^2}{\partial z^2}\mathrm{\nabla }·{\mathbf{v}}_1-\frac{{\partial }^2}{\partial t^2}\mathrm{\nabla }·{\mathbf{v}}_1),\end{array}$$(A.6)

where V_A = B₀(4πρ₀)^−1/2 is Alfvén speed. The first order temperature perturbation (T₁) is also obtained by Eq. (A.4) as

$$\begin{array}{c}\hfill T_1=\frac{1}{{\rho }_0}(\frac{m}{k_B}{p}_1-T_0{\rho }_1),\end{array}$$(A.7)

Replacing the perturbed physical parameters in Eq. (A.5) with their expressions from Eqs. (A.2, A.6, A.7) and rearranging the equation by ∇ ⋅ v₁, we derive the following differential equation (the subscript 1 is omitted):

$$\begin{array}{cc}& \frac{{\partial }^2}{\partial t^2}[\frac{{\partial }^2}{\partial t^2}(\mathrm{\nabla }·\mathbf{v})-(V_A^2+C_s^2){\mathrm{\nabla }}^2(\mathrm{\nabla }·\mathbf{v})]+V_A^2{C}_s^2\frac{{\partial }^2}{\partial z^2}{\mathrm{\nabla }}^2(\mathrm{\nabla }·\mathbf{v})\hfill \\ \hfill & =-{\rho }_0(\gamma -1)Q_{\rho }\frac{\partial }{\partial t}{\mathrm{\nabla }}^2(\mathrm{\nabla }·\mathbf{v})\hfill \\ \hfill & -\frac{Q_T}{C_V}\frac{\partial }{\partial t}(\frac{{\partial }^2}{\partial t^2}(\mathrm{\nabla }·\mathbf{v})-V_A^2{\mathrm{\nabla }}^2(\mathrm{\nabla }·\mathbf{v})+V_A^2\frac{{\partial }^2}{\partial z^2}(\mathrm{\nabla }·\mathbf{v}))\hfill \\ \hfill & +(\gamma -1)T_0{Q}_T\frac{\partial }{\partial t}{\mathrm{\nabla }}^2(\mathrm{\nabla }·\mathbf{v}).\hfill \end{array}$$(A.8)

where $C_s=\sqrt{\gamma k_B{T}_0{m}^{-1}}$ and C_V = (γ − 1)⁻¹k_B/m are acoustic speed and specific heat capacity at constant volume, respectively. It is interesting to note when the context of study is focused only on the plasma with lack of thermal misbalance effect, all the terms on the RHS of Eq. (A.8) disappear, and the outcome aligns with Eq. (3a) of Edwin & Roberts (1983).

All Figures

	Fig. 1. Dependence of dimensionless characteristic timescales, ${\overline{\tau }}_1$ (a) and ${\overline{\tau }}_2$ (b), on temperature and wave period calculated using Eq. (11), normalized by the wave period. This analysis considers a heating model characterized by power indices a = −0.5 and b = −3 for density and temperature, respectively.
In the text

Fig. 2. In the coronal condition (VAe = 1.5, VAi = 1, Csi = 0.15, Cse = 0.13 Mm/s) and with $\overline{{\tau }_1}=6.8$, $\overline{{\tau }_2}=2.4$ (corresponding to approximate values of 34 and 12 minutes for τ1 and τ2, respectively, at a temperature of about 5 MK), the study considers a wave with a period of 5 min: (a) numerical solution of the dispersion relation from Eq. (19) for frequencies (phase speed) normalized by the internal acoustic frequency (acoustic speed) (ωr/ωsi ≡ vph/Csi) of the fast sausage fundamental mode without (solid blue line) and with (solid black line) consideration of the thermal misbalance effect, as well as the numerical solution of the dispersion relation from Eq. (27) for the fast kink fundamental mode (m = 1) without (dashed blue line) and with (dashed black line) consideration of the thermal misbalance effect versus the normalized longitudinal wavenumber (kR). (b) Damping rate to frequency ratio (−γ/ωr) of fundamental states of the fast sausage (solid line) and kink (dashed line) modes with consideration of the thermal misbalance effect versus the normalized longitudinal wavenumber (kR).

In the text

Fig. 3. In the coronal condition and $\overline{{\tau }_1}=6.8$, $\overline{{\tau }_2}=2.4$: (a) numerical solution of the dispersion relation Eq. (18), computed by MATLAB, for the frequency (phase speed) normalized by the internal acoustic frequency (acoustic speed) (ℜ𝔢(ω)/ωsi ≡ vph/Csi) of the fundamental slow sausage mode without consideration of the thermal misbalance effect (unaffected) (solid blue line) and with consideration of the thermal misbalance effect (affected) (solid black line), as well as the fundamental slow kink mode unaffected (dashed blue line) and affected (dashed black line) by the thermal misbalance versus the normalized longitudinal wavenumber (kR). (b) Damping rate to the frequency ratio (−γ/ωr) of fundamental states of the affected slow sausage (solid line) and kink (dashed line) modes by the thermal misbalance versus the normalized longitudinal wavenumber (kR). (c) Damping time to the wave period ratio (τd/P) of the affected slow sausage (solid line) and kink (dashed line) modes by the thermal misbalance versus the normalized longitudinal wavenumber (kR). For comparison, the analytical results obtained by Eqs. (25, 32) are shown by the black and blue solid curves in panels d, e, f.

In the text

Fig. 4. In the coronal condition, the figure displays: (a), (b), and (c) the normalized numerical solutions of the dispersion relation Eq. (19) by the internal acoustic speed (Cs) for sausage modes (m = 0) with different periods of 8 min ($\overline{{\tau }_1}=4.25$, $\overline{{\tau }_2}=1.5$) (solid line), 5 min ($\overline{{\tau }_1}=6.8$, $\overline{{\tau }_2}=2.4$) (dashed line), and 3 min ($\overline{{\tau }_1}=11.34$, $\overline{{\tau }_2}=4$) (dotted line); and (d), (e), and (f): the normalized numerical solutions of the dispersion relation Eq. (27) by the internal acoustic speed (Cs) for kink modes (m = 1) with different periods of 8 min ($\overline{{\tau }_1}=4.25$, $\overline{{\tau }_2}=1.5$) (solid line), 5 min ($\overline{{\tau }_1}=6.8$, $\overline{{\tau }_2}=2.4$) (dashed line), and 3 min ($\overline{{\tau }_1}=11.34$, $\overline{{\tau }_2}=4$) (dotted line).

In the text

Fig. 5. In the coronal condition, for temperatures associated with SUMER oscillations with periods of 15.7 minutes at 6.3 MK ($\overline{{\tau }_1}\approx 2.36$, $\overline{{\tau }_2}\approx 0.76$) (dashed line) and 13.3 minutes at 8.9 MK ($\overline{{\tau }_1}\approx 4.89$, $\overline{{\tau }_2}\approx 1.43$) (solid line), as calculated using the analytical expression in Eq. (32), the figure depicts the ratio of the damping time to the wave period (τd/P) of the affected slow transverse MA modes by the thermal misbalance versus kR.

In the text