© ESO 2021

1. Introduction

The presence of a complex magnetic field, which makes the plasma highly structured and dynamic, dominates the dynamics in a wide range of regions of the solar atmosphere. Improvements in solar telescope technology over the past two decades have enabled detections of a myriad of periodic perturbations in a wide range of magnetic structures (see, e.g., Banerjee et al. 2007, and Tomczyk et al. 2007), that are often described in terms of magnetohydrodynamic (MHD) waves. MHD wave theory is a powerful tool to be exploited for plasma diagnostics of these structures, known as solar magneto-seismology (SMS; see reviews by Erdélyi 2006a; 2006b; Andries et al. 2009, and Ruderman & Erdélyi 2009). SMS combines the measured characteristics of waves with theoretical MHD modelling, and makes it possible to determine waveguide parameters that are hard to measure directly, like the coronal magnetic field strength (Nakariakov & Ofman 2001; Erdélyi & Taroyan 2008, and Kuridze et al. 2019).

One of the first approaches used for modelling solar atmospheric waveguides was based on representing the magnetic flux as a straight homogeneous magnetic cylinder (see, e.g., Ryutov & Ryutova 1976, and Edwin & Roberts 1983). Since then, more complex and more realistic models have been developed (for recent reviews, see, e.g., Nakariakov & Verwichte 2005, Ruderman & Erdélyi 2009 and Nakariakov et al. 2016).

It is customary to categorise the wave modes in a magnetic flux tube depending on the azimuthal wavenumber m (see, e.g., Edwin & Roberts 1983), where m = 0 corresponds to sausage waves, modes with m = 1 are the kink modes, and modes with m ≥ 2 are called the fluting modes. The properties of sausage and kink modes have been studied in a significant number of articles; on the contrary, fluting modes are not getting sufficient attention which can probably be explained by the absence of observations associated with fluting modes. A reason for the lack of observed detection of these modes may be related to the need of high spatial and temporal resolution, where the former is not that easy to achieve.

Nevertheless, there is a great interest in investigating the fluting modes with a theoretical approach as analysing these modes may provide further insights into the sub-resolution structure of a waveguide, and into their high potentials contributing to heating. An initial insight, popular today, into the properties of the oscillations of the flux tube with fluting modes is given by Edwin & Roberts (1983); a homogeneous flux tube was studied, and it was shown how the sausage, kink, and fluting modes behave under solar atmospheric circumstances. However, Ruderman et al. (2017), among others, has theoretically revealed for the kink modes that the dispersion is rather different for a non-homogeneous flux tube with plasma flows even in the thin tube and thin boundary approximation. Moreover, numerical experiments and contemporary theory suggest that the kink and fluting waves could couple in a non-linear fashion (see, e.g., Ruderman et al. 2010; Ruderman 2017, and Terradas et al. 2018). Another interesting aspect is that high azimuthal wavenumbers could be determinant in the cause of the Kelvin-Helmholtz instability (see, e.g., Terradas et al. 2008, Antolin et al. 2014, and Magyar & Van Doorsselaere 2016). As a result, additional investigation of the behaviour of the fluting oscillations could complement the analysis of the sausage and kink modes. Examples of the latter are Pascoe et al. (2016, 2017), Nelson et al. (2019), and Shukhobodskaia et al. (2021), among others, where significant development has lead to a better understanding of kink oscillations of coronal loops, although contemporary theory is still unable to perfectly match the observed oscillation patterns.

In this work, we further extend the study of oscillations of coronal loops. Coronal loops with elliptical cross-sections and a constant density profile have been studied previously in both cold (Ruderman 2003), and finite-β (Erdélyi & Morton 2009) plasmas. Fluting oscillation, central to the current paper, in the absence of background flows and magnetic flux tube expansion was studied comprehensively by Soler (2017). In our work we now consider fluting modes in expanding cylindrical flux tubes in the presence of background bulk flows and a variable density.

The paper is organised as follows. Section 2 contains the formulation of the problem and recalls the linearised MHD equation for the model that is studied here. In Sect. 3 we derive the governing equation for fluting oscillations in the thin tube approximation in the case of the absence of transitional layer between the tube interior and exterior. In Sect. 4 we conduct a general analysis of the eigenvalue problem describing the fluting perturbations in the presence of a stationary flow. Section 5 consists of the derivation of an adiabatic invariant for the fluting oscillations of an expanding magnetic flux tube with slowly varying density. An analytical and a numerical study of the effect of cooling of the coronal magnetic loops, associated with the invariant, are carried out. Finally, Sect. 5 contains the summary of the obtained results and the main conclusions of the work.

2. Proposed model

In this article, we consider a coronal loop as a thin straight expanding magnetic flux tube with circular cross-section R(z). Figure 1 shows the sketch of the equilibrium configuration of the proposed model. In what follows, the cylindrical coordinates r, ϕ, and z are used. The density inside the tube ρ_i(t, r, z) varies with time along and across the magnetic flux tube; the density in the surrounding plasma ρ_e(t, r, z) and in the transitional layer ρ_t(t, r, z) are considered and are given by

$$\begin{array}{c}\hfill \rho =\{\begin{array}{c}{\rho }_{\mathrm{i}}(t,r,z),0\le r\le R(z)(1-l/2),\\ {\rho }_t(t,r,z),R(z)(1-l/2)\le r\le R(z)(1+l/2),\\ {\rho }_{\mathrm{e}}(t,r,z),r\ge R(z)(1+l/2).\end{array}\end{array}$$(1)

Fig. 1. Equilibrium configuration of the straight expanding two-layered magnetic flux tube in the presence of a background flow.

Here l is a constant determining the thickness of the transitional layer sandwiched between the loops’s interior and exterior. It is assumed that ρ(t, r, z) is a continuous function, and ρ_t(t, r, z) is a monotonically decreasing function of r. The reason why the tube is split into core and transitional layers is because we want to take into account damping due to resonance absorption. The time-independent non-twisted equilibrium magnetic field is B = (B_r(r, z),0, B_z(r, z)). Therefore, the divergence-free condition ∇ ⋅ B = 0 for the magnetic field can be written as

$$\begin{array}{c}\hfill \frac{1}{r}\frac{\partial (r{B}_r)}{\partial r}+\frac{\partial B_z}{\partial z}=0.\end{array}$$(2)

It follows from this equation that B can be expressed in terms of the magnetic flux function ψ as

$$\begin{array}{c}\hfill B_r=-\frac{1}{r}\frac{\partial \psi }{\partial z},B_z=\frac{1}{r}\frac{\partial \psi }{\partial r}.\end{array}$$(3)

In addition, we assume that the boundaries of the transitional layer are magnetic surfaces. Therefore, equations r = R(z)(1 − l/2) and r = R(z)(1 + l/2) can be written as ψ = ψ_i = const and ψ = ψ_e = const, respectively, where indices i and e refer to internal and external values. In what follows we use the cold plasma approximation and the thin tube thin boundary (TTTB) approximation. Therefore, it follows from the equilibrium configuration that the magnetic field must be potential. We have ∇ × B = 0 as a result, which can be rewritten in polar coordinates as

$$\begin{array}{c}\hfill \frac{\partial B_r}{\partial z}=\frac{\partial B_z}{\partial r}.\end{array}$$(4)

Substituting Eqs. (3) into (4), we express the latter equation in terms of magnetic flux function ψ as

$$\begin{array}{c}\hfill r\frac{\partial }{\partial r}(\frac{1}{r}\frac{\partial \psi }{\partial r})+\frac{{\partial }^2\psi }{\partial z^2}=0.\end{array}$$(5)

There is also a time-dependent background flow U = (U_r(t, r, z),0, U_z(t, r, z)) present. It is assumed that the background flow velocity is parallel to the equilibrium magnetic field, U ∥ B. The plasma density and velocity are related by the mass conservation equation

$$\begin{array}{c}\hfill \frac{\partial \rho }{\partial t}+\frac{1}{r}\frac{\partial (r\rho U_r)}{\partial r}+\frac{\partial (\rho U_z)}{\partial z}=0.\end{array}$$(6)

The perturbations of the magnetic field and plasma velocity, b = (b_r, b_ϕ, b_z) and U = (u_r, u_ϕ, u_z), are described by the linearised MHD equations in the cold plasma approximation,

$$\begin{array}{c}\hfill \frac{\partial \mathit{u}}{\partial t}+(\mathit{U}·\mathrm{\nabla })\mathit{u}+(\mathit{u}·\mathrm{\nabla })\mathit{U}=\frac{1}{{\mathrm{\mu }}_0\rho }(\mathrm{\nabla }\times \mathit{b})\times \mathit{B},\end{array}$$(7)

$$\begin{array}{c}\hfill \frac{\partial \mathit{b}}{\partial t}=\mathrm{\nabla }\times (\mathit{u}\times \mathit{B}+\mathit{U}\times \mathit{b}),\end{array}$$(8)

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathit{b}=0,\end{array}$$(9)

where μ₀ is the magnetic permeability of free space.

3. Derivation of governing equation

We first introduce the plasma displacement ξ = x − x₀, where x(t, a) is the trajectory of the plasma element and a is the initial position. The components perpendicular to magnetic field lines of the displacement and velocity are given by

$$\begin{array}{c}\hfill {\xi }_{\perp }={\xi }_r{b}_{0z}-{\xi }_z{b}_{0r}{u}_{\perp }=u_r{b}_{0z}-u_z{b}_{0r}.\end{array}$$(10)

We assume that B_z > 0. Then, since the magnetic flux tube is thin, we can approximate ψ by its first terms of a Taylor’s expansion. In addition, we consider the tube axis as the magnetic field line ψ = const. at r = 0. Therefore, without loss of generality, we assume that ψ = 0 at r = 0, and we have from Eq. (3) that

$$\begin{array}{c}\hfill \psi =\frac{1}{2}{r}^{2}h(z).\end{array}$$(11)

Ruderman et al. (2017) showed that

$$\begin{array}{c}\hfill h(z)R^2(z)=C=\mathrm{const}.\end{array}$$(12)

It is easy to see that the analysis of Ruderman et al. (2017) is still valid up to Eq. (57) therein. Thus, we use Eqs. (53)–(57) of Ruderman et al. (2017) and take perturbations proportional to $e^{i\stackrel{\sim }{m}\varphi }$, where $\stackrel{\sim }{m}$ is the azimuthal wavenumber, to obtain

$$\begin{array}{c}\hfill \mathcal{O}({ϵ}^2)=r{B}^2\frac{\partial w}{\partial \psi }+\frac{i\stackrel{\sim }{m}{B}^2{\xi }_{\varphi }}{r}+B_z\frac{w}{r},\end{array}$$(13)

$$\begin{array}{c}\hfill {\stackrel{\sim }{u}}_{\perp }=\frac{1}{B}\frac{\partial w}{\partial T}+\frac{U_z}{\mathit{rB}}\frac{\partial \mathit{rw}}{\partial Z},\end{array}$$(14)

$$\begin{array}{c}\hfill {\stackrel{\sim }{u}}_{\varphi }=\frac{\partial {\xi }_{\varphi }}{\partial T}+r{U}_z\frac{\partial }{\partial Z}(\frac{{\xi }_{\varphi }}{r}),\end{array}$$(15)

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{u}}_{\perp }}{\partial T}+rU{B}_z\frac{\partial }{\partial Z}(\frac{{\stackrel{\sim }{u}}_{\perp }}{\mathit{rB}})=\frac{rB{B}_z}{{\mathrm{\mu }}_o\rho }\frac{\partial }{\partial Z}(\frac{B_z}{r^2{B}^2}\frac{\partial \mathit{rw}}{\partial Z})-\frac{r{B}^3}{\rho }\frac{\partial Q}{\partial \psi }+\mathcal{O}({ϵ}^2),\end{array}$$(16)

$$\begin{array}{c}\hfill \frac{\partial \stackrel{\sim }{u_{\varphi }}}{\partial T}+\frac{U_z}{r}\frac{\partial r{\stackrel{\sim }{u}}_{\varphi }}{\partial Z}=\frac{B_z}{{\mathrm{\mu }}_0\rho r}\frac{\partial }{\partial Z}[r^2{B}_z\frac{\partial }{\partial Z}(\frac{{\xi }_{\varphi }}{r})]-\frac{i\stackrel{\sim }{m}{B}^{2}Q}{\rho r}.\end{array}$$(17)

Now, eliminating ξ_ϕ, ${\stackrel{\sim }{u}}_{\perp }$, and ${\stackrel{\sim }{u}}_{\varphi }$, and only keeping the leading terms with respect to ϵ, the system of equations with the respect to Q and w transforms to

$$\begin{array}{cc}\hfill Q=& \frac{1}{{\mathrm{\mu }}_0{\stackrel{\sim }{m}}^{2}B}\frac{\partial }{\partial Z}(r^{2}B\frac{{\partial }^{2}W}{\partial \psi \partial Z})\hfill \\ \hfill & -\frac{\rho }{{\stackrel{\sim }{m}}^2{B}^2}(r^2\frac{\partial }{\partial T}+U\frac{\partial }{\partial Z}{r}^2)(\frac{{\partial }^{2}W}{\partial T\partial \psi }+U\frac{{\partial }^{2}W}{\partial \psi \partial Z}),\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill (\frac{\partial }{\partial T}& +r^{2}U{B}^2\frac{\partial }{\partial Z}\frac{1}{r^2{B}^2})(\frac{\partial W}{\partial T}+U\frac{\partial W}{\partial Z})\hfill \\ \hfill & =\frac{r^2{B}^3}{{\mathrm{\mu }}_0\rho }\frac{\partial }{\partial Z}(\frac{1}{r^{2}B}\frac{\partial W}{\partial Z})-\frac{r^2{B}^4}{\rho }\frac{\partial Q}{\partial \psi },\hfill \end{array}$$(19)

where W = rw = Brξ_⊥. Then, since ρ, B, and U are quantities that are independent of r and thus of ψ, employing Eq. (11) we obtain

$$\begin{array}{c}\hfill Q=\frac{2\psi }{{\stackrel{\sim }{m}}^{2}h}\frac{\partial }{\partial \psi }[\frac{1}{{\mathrm{\mu }}_0}\frac{{\partial }^{2}W}{\partial Z^2}-\frac{\rho }{h^2}(\frac{\partial }{\partial T}+hU\frac{\partial }{\partial Z}\frac{1}{h})(\frac{\partial W}{\partial T}+U\frac{\partial W}{\partial Z})],\end{array}$$(20)

$$\begin{array}{c}\hfill (\frac{\partial }{\partial T}+hU\frac{\partial }{\partial Z}\frac{1}{h})(\frac{\partial W}{\partial T}+U\frac{\partial W}{\partial Z})=\frac{h^2}{{\mathrm{\mu }}_0\rho }\frac{{\partial }^{2}W}{\partial Z^2}-\frac{2\psi h^3}{\rho }\frac{\partial Q}{\partial \psi }.\end{array}$$(21)

After differentiating Eq. (21) with respect to ψ and employing Eq. (20), we have

$$\begin{array}{c}\hfill 4\psi \frac{\partial }{\partial \psi }(\psi \frac{\partial Q}{\partial \psi })-{\stackrel{\sim }{m}}^{2}Q=0.\end{array}$$(22)

The solutions to Eq. (22) must decay as ψ → ∞ and be regular at ψ = 0. Therefore, as a result we have

$$\begin{array}{c}\hfill Q=\{\begin{array}{c}{Q}_{\mathrm{i}}(T,Z){\psi }^{\stackrel{\sim }{m}/2},0\le \psi \le {\psi }_{\mathrm{i}},\\ Q_{\mathrm{e}}(T,Z){\psi }^{-\stackrel{\sim }{m}/2},\psi \ge {\psi }_{\mathrm{e}},\end{array},\end{array}$$(23)

where Q_i and Q_e are arbitrary functions. The internal and external boundaries are respectively defined by the equations ψ = ψ_i and ψ = ψ_e, where ψ_i and ψ_e are constants. Next, it follows from Eqs. (11) and (12) that

$$\begin{array}{c}\hfill {\psi }_{\mathrm{i}}=\frac{C}{2}(1-\frac{l}{2}{)}^2,{\psi }_{\mathrm{e}}=\frac{C}{2}(1+\frac{l}{2}{)}^2.\end{array}$$(24)

It follows from Eqs. (21) and (23) that in the core layer of the magnetic flux tube $W={\psi }^{\stackrel{\sim }{m}/2}\stackrel{\sim }{W}(T,Z)$, where $\stackrel{\sim }{W}(T,Z)$ is an arbitrary function. Now, we substitute Eq. (23) in Eq. (21) to obtain

$$\begin{array}{c}\hfill (\frac{\partial }{\partial T}+h{U}_{\mathrm{i}}\frac{\partial }{\partial Z}\frac{1}{h})(\frac{\partial W_{\mathrm{i}}}{\partial T}+U_{\mathrm{i}}\frac{\partial W_{\mathrm{i}}}{\partial Z})-\frac{h^2}{{\mathrm{\mu }}_0{\rho }_{\mathrm{i}}}\frac{{\partial }^2{W}_{\mathrm{i}}}{\partial Z^2}=-\frac{\stackrel{\sim }{m}{\psi }_{\mathrm{i}}^{\stackrel{\sim }{m}/2}{h}^3{Q}_{\mathrm{i}}}{{\rho }_{\mathrm{i}}},\end{array}$$(25)

$$\begin{array}{c}\hfill (\frac{\partial }{\partial T}+h{U}_{\mathrm{e}}\frac{\partial }{\partial Z}\frac{1}{h})(\frac{\partial W_{\mathrm{e}}}{\partial T}+U_{\mathrm{e}}\frac{\partial W_{\mathrm{e}}}{\partial Z})-\frac{h^2}{{\mathrm{\mu }}_0{\rho }_{\mathrm{e}}}\frac{{\partial }^2{W}_{\mathrm{e}}}{\partial Z^2}=\frac{\stackrel{\sim }{m}{h}^3{Q}_{\mathrm{e}}}{{\rho }_{\mathrm{e}}{\psi }_{\mathrm{e}}^{\stackrel{\sim }{m}/2}},\end{array}$$(26)

where W_i and W_e are calculated at ψ = ψ_i and ψ = ψ_e, respectively. We now introduce the new variable

$$\begin{array}{c}\hfill \eta =\frac{1}{R(z)}{\xi }_{\perp }.\end{array}$$(27)

In the case $\stackrel{\sim }{m}=1$, we note that ξ_⊥ is independent of r; therefore, the magnetic flux tube oscillates as a solid. In addition, Eqs. (25) and (26) reduce to the ones obtained by Ruderman et al. (2017). We now introduce jumps to this new variable and the magnetic pressure perturbation across the transitional layer:

$$\begin{array}{c}\hfill \delta \eta =\eta {|}_{\psi ={\psi }_{\mathrm{e}}}-\eta {|}_{\psi ={\psi }_{\mathrm{i}}},\delta P=P{|}_{\psi ={\psi }_{\mathrm{e}}}-P{|}_{\psi ={\psi }_{\mathrm{i}}}.\end{array}$$(28)

We have the estimates $\delta \eta =\mathcal{O}(l)$, $\delta P=\mathcal{O}(l)$. Then, employing Eq. (23) we obtain

$$\begin{array}{cc}& W_{\mathrm{i}}=C(1-l/2)\eta ,W_{\mathrm{e}}=C[(1+l/2)\eta +\delta \eta ],\hfill \\ \hfill & {\displaystyle {\psi }_{\mathrm{i}}^{\stackrel{\sim }{m}/2}{h}^3{Q}_{\mathrm{i}}=\frac{{ϵ}^{-2}C}{R^2}P{|}_{\psi ={\psi }_{\mathrm{i}}},{\psi }_{\mathrm{e}}^{-\stackrel{\sim }{m}/2}{h}^3{Q}_{\mathrm{e}}=\frac{{ϵ}^{-2}C}{R^2}P{|}_{\psi ={\psi }_{\mathrm{e}}}.}\hfill \end{array}$$(29)

Now, let us multiply Eq. (25) by ρ_i, Eq. (26) by ρ_e, add the results, using Eqs. (12), (28), and (29), and return to the original non-scaled independent variables, yielding

$$\begin{array}{cc}& {\rho }_{\mathrm{i}}(\frac{\partial }{\partial t}+\frac{U_{\mathrm{i}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial \eta }{\partial t}+U_{\mathrm{i}}\frac{\partial \eta }{\partial z})\hfill \\ \hfill & x+{\rho }_{\mathrm{e}}(\frac{\partial }{\partial t}+\frac{U_{\mathrm{e}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial \eta }{\partial t}+U_{\mathrm{e}}\frac{\partial \eta }{\partial z})-\frac{2B^2}{{\mathrm{\mu }}_0}\frac{{\partial }^2\eta }{\partial z^2}=\mathcal{L},\hfill \end{array}$$(30)

with

$$\begin{array}{cc}\hfill \mathcal{L}& =\frac{\stackrel{\sim }{m}\delta P}{R^2}+\frac{B^2}{{\mathrm{\mu }}_0}\frac{{\partial }^2\delta \eta }{\partial z^2}\hfill \\ \hfill & -{\rho }_{\mathrm{e}}(\frac{\partial }{\partial t}+\frac{U_{\mathrm{e}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial }{\partial t}+U_{\mathrm{e}}\frac{\partial }{\partial z})(\frac{l\eta }{2}+\delta \eta )\hfill \\ \hfill & +{\rho }_{\mathrm{i}}(\frac{\partial }{\partial t}+\frac{U_{\mathrm{i}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial }{\partial t}+U_{\mathrm{i}}\frac{\partial }{\partial z})\frac{l\eta }{2}.\hfill \end{array}$$(31)

For 0 ≤ ψ ≤ ψ_i we have

$$\begin{array}{cc}& (\frac{\partial }{\partial t}+\frac{U_{\mathrm{i}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial \eta }{\partial t}+U_{\mathrm{i}}\frac{\partial \eta }{\partial z})-\frac{B^2}{{\mathrm{\mu }}_0{\rho }_{\mathrm{i}}}\frac{{\partial }^2\eta }{\partial z^2}=\hfill \\ \hfill & {\left(\frac{\widehat{r}}{1-l/2}\right)}^{\stackrel{\sim }{m}-1}((\frac{\partial }{\partial t}+\frac{U_{\mathrm{i}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial {\eta }_{\mathrm{i}}}{\partial t}+U_{\mathrm{i}}\frac{\partial {\eta }_{\mathrm{i}}}{\partial z})-\frac{B^2}{{\mathrm{\mu }}_0{\rho }_{\mathrm{i}}}\frac{{\partial }^2{\eta }_{\mathrm{i}}}{\partial z^2}),\hfill \end{array}$$(32)

whereas for ψ ≥ ψ_e we have

$$\begin{array}{cc}& (\frac{\partial }{\partial t}+\frac{U_{\mathrm{e}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial \eta }{\partial t}+U_{\mathrm{e}}\frac{\partial \eta }{\partial z})-\frac{B^2}{{\mathrm{\mu }}_0{\rho }_{\mathrm{e}}}\frac{{\partial }^2\eta }{\partial z^2}=\hfill \\ \hfill & {\left(\frac{1+l/2}{\widehat{r}}\right)}^{\stackrel{\sim }{m}-1}((\frac{\partial }{\partial t}+\frac{U_{\mathrm{e}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial {\eta }_{\mathrm{e}}}{\partial t}+U_{\mathrm{e}}\frac{\partial {\eta }_{\mathrm{e}}}{\partial z})-\frac{B^2}{{\mathrm{\mu }}_0{\rho }_{\mathrm{e}}}\frac{{\partial }^2{\eta }_{\mathrm{e}}}{\partial z^2}),\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill & {\xi }_{\varphi }=\frac{\mathit{ir}}{\stackrel{\sim }{m}}\frac{\partial \mathit{rw}}{\partial \psi },\hfill \end{array}$$(34)

where $\widehat{r}=r/R(z)={(\psi /{\psi }_{\mathrm{i}})}^{1/2}$, η_i = η(ψ_i), and η_e = η(ψ_e). Equations (30)–(33) govern oscillations of the magnetic flux tube for $\stackrel{\sim }{m}>0$ in the absence of a transitional layer (i.e. l = 0). This system is valid both for standing waves and propagating waves. In the case l ≠ 0 the system is not closed, and the jumps across the transitional layer δη and δP should be determined to close the governing system. We should also note that, in the case $\stackrel{\sim }{m}=1$, Eq. (32) is fulfilled identically since η is independent of the r-component. In this case the above system of governing equation coincides with that obtained by Ruderman et al. (2017); the detailed investigation carried out is available in Shukhobodskiy & Ruderman (2018), Shukhobodskiy et al. (2018), Ruderman et al. (2019), and Ruderman & Petrukhin (2019), among others.

3.1. Fluting modes in the absence of transitional layer

We assume that the magnetic tube has a sharp boundary, meaning that l = 0. Since the transitional layer is absent in this equilibrium, we use Eq. (30) with the right-hand side equal to 0. Hence, we rewrite Eq. (30) as

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{i}}(\frac{\partial }{\partial t}& +\frac{U_{\mathrm{i}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial \eta }{\partial t}+U_{\mathrm{i}}\frac{\partial \eta }{\partial z})+\hfill \\ \hfill & {\rho }_{\mathrm{e}}(\frac{\partial }{\partial t}+\frac{U_{\mathrm{e}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial \eta }{\partial t}+U_{\mathrm{e}}\frac{\partial \eta }{\partial z})-\frac{2B^2}{{\mathrm{\mu }}_0}\frac{{\partial }^2\eta }{\partial z^2}=0.\hfill \end{array}$$(35)

Using Eq. (32), we reduce it to

$$\begin{array}{cc}& (\frac{\partial }{\partial t}+\frac{U}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial \eta }{\partial t}+U\frac{\partial \eta }{\partial z})-\frac{B^2}{{\mathrm{\mu }}_0\rho }\frac{{\partial }^2\eta }{\partial z^2}\hfill \\ \hfill & ={\widehat{r}}^{\stackrel{\sim }{m}-1}((\frac{\partial }{\partial t}+\frac{U_{\mathrm{i}}}{R^2}\frac{\partial }{\partial z}{R}^2)(\frac{\partial {\eta }_{\mathrm{i}}}{\partial t}+U_{\mathrm{i}}\frac{\partial {\eta }_{\mathrm{i}}}{\partial z})-\frac{B^2}{{\mathrm{\mu }}_0{\rho }_{\mathrm{i}}}\frac{{\partial }^2{\eta }_{\mathrm{i}}}{\partial z^2}).\hfill \end{array}$$(36)

These equations are governing for fluting oscillations and contain general information regarding the system in the absence of transitional layer. Now we employ these equations in order to facilitate the analysis.

3.2. Eigenvalue problem in the presence of stationary flow: General analysis

In this section we assume that the density and flow velocity are both independent of time and that the external plasma is at rest (i.e. U_e = 0). Since we assume that the characteristic scale of variation of the equilibrium quantities in the radial direction is the same as in the axial direction, we can neglect their radial variation inside the tube in the thin tube approximation. Then, it follows from the mass conservation in Eq. (6) that the density and flow speed are related by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{i}}U{R}^2=\mathrm{const},\end{array}$$(37)

where we have dropped the subscript ‘i’ for the internal flow speed U. It follows from the magnetic flux conservation that

$$\begin{array}{c}\hfill B{R}^2=\mathrm{const}.\end{array}$$(38)

Below, we consider standing waves and assume that the magnetic flux tube ends are fixed at the dense photosphere. Thus, we impose the boundary conditions

$$\begin{array}{c}\hfill \eta =0\mathrm{at}z=\pm L/2,\end{array}$$(39)

where L is the tube length. We look for stationary solutions, so we take η ∝ e^−iωt. Then, Eq. (36) reduces to

$$\begin{array}{cc}& (\frac{B^2}{{\mathrm{\mu }}_0}-{\rho }_{\mathrm{i}}{U}_{\mathrm{i}}^2)(\frac{{\partial }^2\eta }{\partial z^2}-{\widehat{r}}^{\stackrel{\sim }{m}-1}\frac{{\partial }^2{\eta }_{\mathrm{i}}}{\partial z^2})+(2i\omega {\rho }_{\mathrm{i}}{U}_{\mathrm{i}}-\frac{{\rho }_{\mathrm{i}}{U}_{\mathrm{i}}}{R^2}\frac{\partial U_{\mathrm{i}}{R}^2}{\partial z})\hfill \\ \hfill & \times (\frac{\partial \eta }{\partial z}-{\widehat{r}}^{\stackrel{\sim }{m}-1}\frac{\partial {\eta }_{\mathrm{i}}}{\partial z})+(\frac{2i\omega {\rho }_{\mathrm{i}}{U}_{\mathrm{i}}}{R}\frac{\partial R}{\partial z}+{\omega }^2{\rho }_{\mathrm{i}})(\eta -{\widehat{r}}^{\stackrel{\sim }{m}-1}{\eta }_{\mathrm{i}})=0.\hfill \end{array}$$(40)

This equation, with the boundary conditions Eq. (39), constitute an eigenvalue problem. Here we now study the general properties of this problem. We are interested in the conditions for which the tube will be stable with respect to long flute perturbations. Thus, we assume that μ₀ρ_iU² < 2B², that is $U^2<2V_{Ai}^2\equiv 2B^2/{\mathrm{\mu }}_0{\rho }_{\text{i}}$, and this condition is satisfied for all z ∈ [ − L/2, L/2]. Then, we make the variable substitution

$$\begin{array}{c}\hfill \eta -{\widehat{r}}^{\stackrel{\sim }{m}-1}{\eta }_{\mathrm{i}}=qexp(-i\omega \sigma (z)),\sigma ={\displaystyle {\int }_{-L/2}^z\frac{{\mathrm{\mu }}_0{\rho }_{\mathrm{i}}U\mathrm{d}z\prime }{B^2-{\mathrm{\mu }}_0{\rho }_{\mathrm{i}}{U}^2}.}\end{array}$$(41)

Now, multiplying the obtained equation by R⁴, we rewrite Eq. (40) as

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}z}\left[R^4(\frac{B^2}{{\mathrm{\mu }}_0}-{\rho }_{\mathrm{i}}{U}^2)\frac{\mathrm{d}q}{\mathrm{d}z}\right]+{\omega }^{2}W(z)q=0,\end{array}$$(42)

where

$$\begin{array}{c}\hfill W(z)=\frac{R^4{\rho }_{\mathrm{i}}{B}^2}{B^2-{\mathrm{\mu }}_0{\rho }_{\mathrm{i}}{U}^2}.\end{array}$$(43)

Equation (42) with the boundary conditions Eq. (39) constitute the classical Sturm-Liouville problem with the eigenvalue ω². It follows from the general theory of the Sturm-Liouville problem (e.g. Coddington & Levinson 1955) that the eigenvalues are real and constitute a monotonically increasing unbounded sequence. The first eigenvalue is the square of the fundamental frequency, and the corresponding eigenfunction has no nodes in ( − L/2, L/2). All other eigenvalues are the squares of frequencies of corresponding overtones. The eigenfunction corresponding to the nth overtone has n − 1 nodes in ( − L/2, L/2).

It is clear that we can always take q to be real. Multiplying Eq. (42) by q, integrating the obtained equation, and using Eq. (43) we obtain

$$\begin{array}{c}\hfill {\omega }^2{\displaystyle {\int }_{-L/2}^{L/2}W(z)q^2\mathrm{d}z={\int }_{-L/2}^{L/2}{R}^4(\frac{B^2}{{\mathrm{\mu }}_0}-{\rho }_{\mathrm{i}}{U}^2){\left(\frac{\mathrm{d}q}{\mathrm{d}z}\right)}^2\mathrm{d}z.}\end{array}$$(44)

Since W(z) > 0, it follows from this equation that ω² > 0, meaning that all the eigenfrequencies are real. This also implies that the inequality $U^2<V_{Ai}^2$ is a sufficient condition for stability inside the tube with respect to long flute perturbations. To obtain an equation describing the boundary, we follow the analysis of Ruderman et al. (2017) to determine whether the stability condition on the boundary is different, with a sufficient condition $U^2<2V_{Ai}^2$. We note that this condition is sufficient but not necessary. On the other hand, this analysis does not prove that the tube is stable if the condition $U^2<V_{Ai}^2$ is satisfied. While it guaranties that the tube is stable with respect to long flute perturbations, it can be unstable with respect to other types of perturbations.

4. Flute oscillations of magnetic flux tube with slowly varying density

Below, we now consider flute oscillations of a magnetic flux tube with slowly varying density. We assume that the magnetic flux tube is fixed at the dense photosphere, which means that the boundary condition in Eq. (40) is satisfied. In addition, we also assume that the characteristic time of density variation, t_ch, is much greater than the characteristic time of oscillations. We note that this is not a very restrictive assumption.

4.1. Derivation of adiabatic invariant

We denote the ratio of t_ch to the characteristic time of oscillations as ν⁻¹, where ν ≪ 1. Hence, the characteristic time of oscillations is νt_ch. On the other hand, it is also of the order of the loop length divided by the characteristic flute speed, which is equal to $L/(B/\sqrt{{\mathrm{\mu }}_0{\rho }_{\mathrm{ch}}})$, where ρ_ch is the characteristic density. We have, as a result,

$$\begin{array}{c}\hfill B\sim {\nu }^{-1}\sqrt{{\mathrm{\mu }}_0{\rho }_{\mathrm{ch}}}\frac{L}{t_{\mathrm{ch}}}.\end{array}$$(45)

Following this estimate, let us now introduce the scaled magnetic field $\stackrel{\sim }{B}=\nu B$. Then, we rewrite Eqs. (33) and (42) as

$$\begin{array}{cc}& (\frac{{\nu }^{-2}{\stackrel{\sim }{B}}^2}{{\mathrm{\mu }}_0}-{\rho }_{\mathrm{i}}{U}_{\mathrm{i}}^2)(\frac{{\partial }^2\eta }{\partial z^2}-{\widehat{r}}^{\stackrel{\sim }{m}-1}\frac{{\partial }^2{\eta }_{\mathrm{i}}}{\partial z^2})+(2i\omega {\rho }_{\mathrm{i}}{U}_{\mathrm{i}}-\frac{{\rho }_{\mathrm{i}}{U}_{\mathrm{i}}}{R^2}\frac{\partial U_{\mathrm{i}}{R}^2}{\partial z})\hfill \\ \hfill & \times (\frac{\partial \eta }{\partial z}-{\widehat{r}}^{\stackrel{\sim }{m}-1}\frac{\partial {\eta }_{\mathrm{i}}}{\partial z})+(\frac{2i\omega {\rho }_{\mathrm{i}}{U}_{\mathrm{i}}}{R}\frac{\partial R}{\partial z}+{\omega }^2{\rho }_{\mathrm{i}})(\eta -{\widehat{r}}^{\stackrel{\sim }{m}-1}{\eta }_{\mathrm{i}})=0,\hfill \end{array}$$(46)

$$\begin{array}{cc}& (\frac{{\nu }^{-2}{\stackrel{\sim }{B}}^2}{{\mathrm{\mu }}_0}-{\rho }_{\mathrm{e}}{U}_{\mathrm{e}}^2)(\frac{{\partial }^2\eta }{\partial z^2}-{\widehat{r}}^{1-\stackrel{\sim }{m}}\frac{{\partial }^2{\eta }_{\mathrm{e}}}{\partial z^2})\hfill \\ \hfill & +(2i\omega {\rho }_{\mathrm{e}}{U}_{\mathrm{e}}-\frac{{\rho }_{\mathrm{e}}{U}_{\mathrm{e}}}{R^2}\frac{\partial U_{\mathrm{e}}{R}^2}{\partial z})(\frac{\partial \eta }{\partial z}-{\widehat{r}}^{1-\stackrel{\sim }{m}}\frac{\partial {\eta }_{\mathrm{i}}}{\partial z})\hfill \\ \hfill & +(\frac{2i\omega {\rho }_{\mathrm{e}}{U}_{\mathrm{e}}}{R}\frac{\partial R}{\partial z}+{\omega }^2{\rho }_{\mathrm{e}})(\eta -{\widehat{r}}^{1-\stackrel{\sim }{m}}{\eta }_{\mathrm{e}})=0.\hfill \end{array}$$(47)

Now, we use the Wentzel-Kramer-Brillouin (WKB) method (e.g. Bender & Orszag 1978) and look for solution to this equation in the form

$$\begin{array}{c}\hfill \eta =S(t,r,z)exp[{\nu }^{-1}\mathrm{\Theta }(t)].\end{array}$$(48)

Then, we expand S in the series

$$\begin{array}{c}\hfill S=S_0+\nu S_1+\cdots \end{array}$$(49)

We substitute Eq. (48) in Eqs. (46) and (47), using Eq. (49), and collecting terms of order ν⁻² we obtain

$$\begin{array}{c}\hfill \frac{{\partial }^2{\stackrel{\sim }{S}}_0}{\partial z^2}+\frac{{\mathrm{\Omega }}^2}{{\stackrel{\sim }{V}}_{\mathit{Ai}}^2}{\stackrel{\sim }{S}}_0=0\end{array}$$(50)

and

$$\begin{array}{c}\hfill \frac{{\partial }^2{\widehat{S}}_0}{\partial z^2}+\frac{{\mathrm{\Omega }}^2}{{\stackrel{\sim }{V}}_{\mathit{Ae}}^2}{\widehat{S}}_0=0,\end{array}$$(51)

where

$$\begin{array}{cc}\hfill {\stackrel{\sim }{S}}_0& =S_0(t,r,z)-{\widehat{r}}^{\stackrel{\sim }{m}-1}{S}_0(t,R(z),z),\hfill \\ \hfill {\widehat{S}}_0& =S_0(t,r,z)-{\widehat{r}}^{1-\stackrel{\sim }{m}}{S}_0(t,R(z),z),\hfill \\ \hfill \mathrm{\Omega }& =\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}t},{\stackrel{\sim }{V}}_{\mathit{Ai}}^2=\frac{{\stackrel{\sim }{B}}^2}{{\mathrm{\mu }}_0{\rho }_{\mathrm{i}}}{\stackrel{\sim }{V}}_{\mathit{Ae}}^2=\frac{{\stackrel{\sim }{B}}^2}{{\mathrm{\mu }}_0{\rho }_{\mathrm{e}}}.\hfill \end{array}$$(52)

Inside the magnetic flux tube ${\widehat{r}}^{\stackrel{\sim }{m}-1}\ll 1$ and outside the flux tube ${\widehat{r}}^{1-\stackrel{\sim }{m}}\ll 1$, hence adding Eqs. (50) and (51), we obtain

$$\begin{array}{c}\hfill \frac{{\partial }^2{S}_0}{\partial z^2}+\frac{{\mathrm{\Omega }}^2}{{\stackrel{\sim }{C}}_k^2}{S}_0=0,\end{array}$$(53)

where

$$\begin{array}{c}\hfill {\stackrel{\sim }{C}}_k^2=\frac{2{\stackrel{\sim }{B}}^2}{{\mathrm{\mu }}_0({\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}})}.\end{array}$$(54)

Since we assume that the magnetic flux tube is fixed in the dense photosphere (i.e. there is line-tying) we now have that

$$\begin{array}{c}\hfill S_0=0\text{at}z=\pm L/2.\end{array}$$(55)

Therefore, Eqs. (53) and (55) constitute the classical Sturm-Liouville problem as a result. An approximation of such order is called the approximation of geometrical optics (e.g. Bender & Orszag 1978). The Sturm-Liouville problem of Eqs. (53) and (55) coincide with a problem obtained by Dymova & Ruderman (2005) to describe oscillations of magnetic flux tubes whose density varies along the tube, and by Ruderman et al. (2008) to describe kink oscillation of magnetic tubes whose density and cross-section radius vary along the tube.

Below, we assume that Ω² is an eigenvalue and S₀ the corresponding eigenfunction. In accordance with the general theory (e.g. Coddington & Levinson 1955), Ω² is real. It is easy to prove that Ω² > 0 by multiplying Eq. (53) by S₀, then integrating by parts with respect to z from −L/2 to L/2 and using the boundary condition (55). We can always assume that S₀ is a real function.

Now, we collect terms of the order ν⁻¹ to obtain

$$\begin{array}{c}\hfill \frac{{\partial }^2{\stackrel{\sim }{S}}_1}{\partial z^2}+\frac{{\mathrm{\Omega }}^2}{{\stackrel{\sim }{V}}_{\mathit{Ai}}^2}{\stackrel{\sim }{S}}_1=\frac{2i\mathrm{\Omega }}{{\stackrel{\sim }{V}}_{\mathit{Ai}}^2}(\frac{\partial {\stackrel{\sim }{S}}_0}{\partial t}+\frac{{\stackrel{\sim }{S}}_0}{2\mathrm{\Omega }}\frac{\partial \mathrm{\Omega }}{\partial t}+\frac{U_{\mathrm{i}}}{R}\frac{\partial R{\stackrel{\sim }{S}}_0}{\partial z})\end{array}$$(56)

and

$$\begin{array}{c}\hfill \frac{{\partial }^2{\widehat{S}}_1}{\partial z^2}+\frac{{\mathrm{\Omega }}^2}{{\stackrel{\sim }{V}}_{\mathit{Ae}}^2}{\stackrel{\sim }{S}}_1=\frac{2i\mathrm{\Omega }}{{\stackrel{\sim }{V}}_{\mathit{Ae}}^2}(\frac{\partial {\widehat{S}}_0}{\partial t}+\frac{{\widehat{S}}_0}{2\mathrm{\Omega }}\frac{\partial \mathrm{\Omega }}{\partial t}+\frac{U_{\mathrm{e}}}{R}\frac{\partial R{\widehat{S}}_0}{\partial z}),\end{array}$$(57)

where

$$\begin{array}{c}\hfill {\stackrel{\sim }{S}}_1=S_1(t,r,z)-{\widehat{r}}^{\stackrel{\sim }{m}-1}{S}_1(t,R(z),z)\end{array}$$(58)

and

$$\begin{array}{c}\hfill {\widehat{S}}_1=S_1(t,r,z)-{\widehat{r}}^{1-\stackrel{\sim }{m}}{S}_1(t,R(z),z).\end{array}$$(59)

Taking into account that inside the magnetic flux tube ${\widehat{r}}^{\stackrel{\sim }{m}-1}\ll 1$ and outside the flux tube ${\widehat{r}}^{1-\stackrel{\sim }{m}}\ll 1$, from Eqs. (56) and (57), we have

$$\begin{array}{c}\hfill \frac{{\partial }^2{S}_1}{\partial z^2}+\frac{{\mathrm{\Omega }}^2}{{\stackrel{\sim }{C}}_k^2}{S}_1=\frac{2i\mathrm{\Omega }}{{\stackrel{\sim }{C}}_k^2}(\frac{\partial S_0}{\partial t}+\frac{S_0}{2\mathrm{\Omega }}\frac{\partial \mathrm{\Omega }}{\partial t}+\frac{U_{\mathrm{e}}}{R}\frac{\partial R{S}_0}{\partial z}).\end{array}$$(60)

This order of approximation is the approximation of the physical optics. S₁ must satisfy the boundary conditions

$$\begin{array}{c}\hfill S_1=0z=\pm L/2.\end{array}$$(61)

The homogeneous counterparts of the boundary value problem constituted by Eqs. (53) and (60), and the boundary condition Eqs. (55) and (61) have the non-trivial solution S₁ = S₀. This implies that the boundary value problem determining S₁ has a solution only if the right-hand side of Eq. (60) satisfies the compatibility condition, which is the condition that it should be orthogonal to S₀. To obtain this condition we multiply Eq. (60) by S₀, then integrate with respect to z from −L/2 to L/2. After some algebra we obtain that

$$\begin{array}{c}\hfill {\displaystyle {\int }_{-L/2}^{L/2}\frac{1}{C_k^2}\frac{\partial \mathrm{\Omega }{S}_0^2}{\partial t}\mathrm{d}z=-\frac{\mathrm{\Omega }}{2}{\int }_{-L/2}^{L/2}\frac{U_{\mathrm{i}}}{R^2{C}_k^2}\frac{\partial (S_0^2{R}^2)}{\partial z}\mathrm{d}z.}\end{array}$$(62)

Using integration by parts, and some algebra, we eventually obtain

$$\begin{array}{c}\hfill \omega {\displaystyle {\int }_{-L/2}^{L/2}\frac{S_0^2}{C_k^2}\mathrm{d}z=\mathrm{const},}\end{array}$$(63)

where

$$\begin{array}{c}\hfill \omega ={\nu }^{-1}\mathrm{\Omega },C_k^2=\frac{2B^2}{{\mathrm{\mu }}_0({\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}})}.\end{array}$$(64)

The left-hand side of Eq. (63) is an adiabatic invariant. Equation (63) states that this invariant is conserved. It is easy to see that the derivation for adiabatic invariant and approximations of the 0th order on the boundary presented by Ruderman et al. (2017) are still valid. It is also worth noting that the adiabatic invariant is the same as that obtained by Ruderman (2011a) for oscillations of magnetic tubes with constant cross-section radius. Hence, the tube expansion only affects the temporal evolution of flute oscillations of magnetic flux tubes with varying density through C_k, V_Ai, V_Ae, and S₀ that all depend on R(z).

4.2. Effect of cooling on flute oscillations of coronal magnetic loops

The main cause of cooling of a coronal loop is radiation. The intensity is proportional to the plasma density squared for optically think plasmas. The energy deposition that can cover the losses of energy in the rarefied external plasma could be too small to cover the energy loss in the internal region of a magnetic flux tube. Therefore, we can assume that the plasma temperature remains constant outside the magnetic flux tube. As was found in Aschwanden & Nightingale (2005), and Morton & Erdélyi (2010), the temperature evolution inside the loop is approximated by exponentially decaying function,

$$\begin{array}{c}\hfill T(t)=T_{0}exp(-t/t_{\mathrm{cool}}),\end{array}$$(65)

where we assume that, at the initial time t = 0, the cooling begins. In addition, we assume that the plasma temperature is the same at this initial time in the external and internal regions, so that the temperature of the external plasma remains equal to T₀. In what follows we consider a magnetic flux tube that has a semi-circular shape. We neglect the effect of flux tube curvature on fluting oscillations. The shape of the loop determines the density variation along the loop only as a result, t. The density of the external plasma is given by the barometric formula

$$\begin{array}{c}\hfill {\rho }_{\mathrm{e}}(z)=\frac{{\rho }_{\mathrm{f}}}{\zeta }exp(\frac{-L}{\pi H_0}cos\frac{\pi z}{L}),\end{array}$$(66)

where ρ_f is the plasma density at the footpoints inside the magnetic flux tube, which is assumed to be constant, and ζ is the ratio of internal to external plasma densities at initial time,

$$\begin{array}{c}\hfill H_0=\frac{k_{\mathrm{B}}{T}_0}{\mathit{mg}},\end{array}$$(67)

k_B is the Boltzmann constant, g is the gravity acceleration, and m is the mean mass per particle equal to one-half of the proton mass for a proton-electron plasma. The flow inside the magnetic flux tube is caused by plasma cooling.

It is easy to show that the analysis of Ruderman (2011b) still holds true even for fluting modes described in the model presented in this manuscript. Thus, the effect of the flow on the plasma density is fairly weak for the typical coronal conditions and the observed cooling times. Therefore, the barometric formula provides a sufficiently good approximation of the plasma density in the internal region of the magnetic flux tube. We have

$$\begin{array}{c}\hfill {\rho }_{\mathrm{i}}(z)={\rho }_{\mathrm{f}}exp(\frac{-L}{\pi H(t)}cos\frac{\pi z}{L}),\end{array}$$(68)

where

$$\begin{array}{c}\hfill H(t)=\frac{k_{\mathrm{B}}T(t)}{\mathit{mg}}.\end{array}$$(69)

Since the background magnetic field is straight and not twisted, we can describe the variation of cross-section radius of the tube with the same expression as in Ruderman et al. (2008),

$$\begin{array}{c}\hfill R(z)=R_f\lambda \sqrt{\frac{cosh(L/2l_c)-1}{cosh(2L/2l_c)-{\lambda }^2+({\lambda }^2-1)cosh(z/l_c)}},\end{array}$$(70)

where R_f is the cross-section radius of the magnetic flux tube at the footpoints, l_c is an arbitrary positive constant with the dimension of length, and λ = R(0)/R_f. Ruderman et al. (2008) showed that for the z-component of the magnetic field to be positive everywhere in the region |z|≤L/2, the expansion factor must satisfy the inequality λ < λ_m, where

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{m}}^2\approx \frac{1.4cosh(L/2l_{\mathrm{c}})}{1+0.4cosh(L/2lc)}.\end{array}$$(71)

Here λ_m is a monotonically increasing function of L/l_c, such that λ_m → 1 as L/l_c → 0, and λ_m → 1.87 as L/l_c → ∞. Since the expansion of coronal loop factor does not exceed 1.5, the careful choice of the ration L/l_c allows us to cover the entire range of the expansion factor variation.

Now, we introduce the following dimensionless variables and parameters

$$\begin{array}{cc}& Z=\frac{2z}{L},\tau =\frac{t}{t_{\mathrm{cool}}},\varpi =\frac{\omega L}{C_{\mathrm{f}}},\kappa =\frac{L}{\pi H_0},\hfill \\ \hfill & \mathrm{\Lambda }=\frac{R}{R_F},S_0\prime =\frac{S_0}{R_{\mathrm{f}}},\hfill \end{array}$$(72)

where the kink speed at the footpoints is

$$\begin{array}{c}\hfill C_f^2=\frac{2\zeta B_{\mathrm{f}}^2}{{\mathrm{\mu }}_0{\rho }_{\mathrm{f}}(\zeta +1)}.\end{array}$$(73)

Therefore, we can reduce Eq. (53) to

$$\begin{array}{cc}\hfill \frac{{\partial }^2{S}_0}{\partial Z^2}& +\frac{{\varpi }^2{\mathrm{\Lambda }}^4(Z)S_0}{4(\zeta +1)}\hfill \\ \hfill & \times [\zeta exp(-\kappa e^{\tau }cos(\pi Z/2)+exp(-\kappa cos(\pi Z/2)]=0.\hfill \end{array}$$(74)

We have dropped the prime in Eq. (74). We should also note that in order to obtain solutions for the plasma displacement in the internal region and in the external region, which both depend on radial component, we first need to find a solution for the plasma displacement at the boundary of Eq. (74). We also know that W = rR(z)Bη, implying that inside the tube $W={\psi }^{\stackrel{\sim }{m}/2}\stackrel{\sim }{W}(t,z)$ and in the external region, in the vicinity of the magnetic flux tube, $W={\psi }^{\stackrel{\sim }{-m}/2}\widehat{W}(t,z)$, where $\widehat{W}(t,z)$ is an arbitrary function. Therefore we obtain

$$\begin{array}{cc}\hfill \eta & =\frac{1}{2}{\left(\frac{{\psi }^{\stackrel{\sim }{m}-1}}{{\psi }^i}\right)}^{\frac{1}{2}}\stackrel{\sim }{W}(t,z)\text{for}0\le r\le R(z)\text{and}\hfill \\ \hfill \eta & =\frac{1}{2}{\left(\frac{1}{{\psi }^{\stackrel{\sim }{m}+1}{\psi }^i}\right)}^{\frac{1}{2}}\stackrel{\sim }{W}(t,z)\text{for}r\ge R(z).\hfill \end{array}$$(75)

The eigenfunction S₀ is accurately determined up to the multiplication of the arbitrary function of τ, and the dependence on ψ is determined by Eq. (75). We can always fix its value to one particular point. We assume X(τ, Z) to be an eigenfunction corresponding to the fundamental mode, which satisfies the condition X(τ, 0) = 1. Then, the general solution to the eigenvalue problem of the fundamental mode is S₀(τ, Z, ψ) = A(τ, ψ)X(τ, Z), where A(τ, ψ) is the oscillation amplitude at Z = 0. With the aid of Eq. (75) we derive

$$\begin{array}{c}\hfill S_0=\frac{1}{2}{\left(\frac{{\psi }^{m-1}}{{\psi }_{\mathrm{i}}}\right)}^{\frac{1}{2}}\widehat{A}(\tau )X(\tau ,Z)\text{for}0\le r\le R(Z)\text{and}\end{array}$$(76)

$$\begin{array}{c}\hfill S_0=\frac{1}{2}{\left(\frac{1}{{\psi }^{m+1}{\psi }_{\mathrm{i}}}\right)}^{\frac{1}{2}}\stackrel{\sim }{A}(\tau )X(\tau ,Z)\text{for}r\ge R(Z),\end{array}$$(77)

where

$$\begin{array}{c}\hfill A(\tau ,\psi )=\frac{1}{2}{\left(\frac{{\psi }^{m-1}}{{\psi }_{\mathrm{i}}}\right)}^{\frac{1}{2}}\widehat{A}(\tau )0\le r\le R(Z)\text{and}\end{array}$$(78)

$$\begin{array}{c}\hfill A(\tau ,\psi )=\frac{1}{2}{\left(\frac{1}{{\psi }^{m+1}{\psi }_{\mathrm{i}}}\right)}^{\frac{1}{2}}\stackrel{\sim }{A}(\tau )r\ge R(Z).\end{array}$$(79)

Therefore, we rewrite Eq. (63) as

$$\begin{array}{cc}\hfill \varpi A^2(\tau ,\psi ){\displaystyle {\int }_{-1}^1}& X^2{\mathrm{\Lambda }}^4[\zeta exp(-\kappa e^{\tau }cos(\pi Z/2))\hfill \\ \hfill & +exp(-\kappa cos(\pi Z/2))]\mathrm{d}Z=\text{const}.\hfill \end{array}$$(80)

We also note that

$$\begin{array}{c}\hfill \frac{A(\tau ,\psi )}{A(0,{\psi }_{\mathrm{i}})}=\frac{\widehat{A}(\tau )}{\widehat{A}(0)}{\left(\frac{r}{R(Z)}\right)}^{m-1}\text{for}0\le r\le R(Z)\text{and}\end{array}$$(81)

$$\begin{array}{c}\hfill \frac{A(\tau ,\psi )}{A(0,{\psi }_{\mathrm{i}})}=\frac{\stackrel{\sim }{A}(\tau )}{\stackrel{\sim }{A}(0)}{\left(\frac{R(Z)}{r}\right)}^{m+1}\text{for}r\ge R(Z).\end{array}$$(82)

Then, we set l_c = 4 so that the entire observed range of possible λ is covered by the expansion model and solve Eqs. (80)–(82) numerically to obtain the dependence of dimensionless amplitude on r for various values of $\stackrel{\sim }{m}$.

First, we note that it follows from Fig. 2 that the higher the value of λ, the lower the amplitude closer to the magnetic flux tube core. The phenomenon that can be spotted is that the amplitude increases between the initial and final time. Furthermore, the higher the expansion factor, the higher the final amplitude, which is consistent with the result obtained by Shukhobodskiy et al. (2018). The other pattern to be noted in Fig. 2 is that when the value of $\stackrel{\sim }{m}$ is higher, more oscillations are localised at the tube boundary. This latter conclusion is in agreement with Soler (2017). Furthermore, we note that the expansion factor λ of the magnetic flux tube is of paramount importance in determining the oscillation intensity. This can be especially observed in the bottom panels of Fig. 2, where oscillations occur close to the tube boundary magnetic field lines to such an extent, that even the region above the boundary in the footpoints would not oscillate at the apex of the tube for even a mildly expanding flux tube.

Fig. 2. Dependence of dimensionless amplitude on r/Rf. Left panels: amplitude variation at the initial time τ = 0, and right panels: amplitude variation at the final time τ = 1. The solid, dash-dotted, and dot-dashed lines correspond to λ = 1, 1.1, and (1.3, 1.5), respectively. From top to bottom, the panels correspond to $\stackrel{\sim }{m}=2$, 3, 5, 10, and 100.

It follows from Fig. 3 that ζ has no effect on the initial amplitude variation in the distance from the core region. On the other hand, at the final time τ = 1 it is possible to observe three separate phenomena. First of all, increasing the wave number $\stackrel{\sim }{m}$ localises the oscillation closer to the boundary. Secondly, when the value of ζ increases, the amplification of oscillations due to cooling become stronger. Finally, the increase in $\stackrel{\sim }{m}$ leads to a reduction in the spread between the oscillation variation across the tube for different values of ζ, preserving the trend that the higher the value of ζ, the higher the amplification due to cooling.

Fig. 3. Dependence of dimensionless amplitude on r/Rf. Left panels: amplitude variation at the initial time τ = 0, and right panels: amplitude variation at final time τ = 1. The solid, dash-dotted, and dot-dashed lines correspond to ζ = 1, 2, and (3, 4), respectively. From top to bottom, the panels correspond to $\stackrel{\sim }{m}=2$, 3, 5, 10, and 100.

Figure 4 represents the evolution of oscillation of magnetic flux tube both through time and distance from the annulus. We find, for both the top and bottom panels, that the amplitude increases with time. Furthermore, higher values of $\stackrel{\sim }{m}$ lead to the localisation of oscillation closer to the boundary. On the other hand, an increasing expansion factor λ not only leads to an increase in the oscillation amplitude, but also shifts the localisation of oscillation of the magnetic flux tube further from the centre. Therefore we conclude that not only does $\stackrel{\sim }{m}$ have an important role in the evolution of the cooling system, but also variations in the expansion λ and density ratio ζ, with the latter having a significant deviation between the initial and final times.

Fig. 4. Dependence of dimensionless amplitude on r/Rf and time τ for ζ = 3. The orange, purple, yellow, and dark purple colour schemes respectively correspond to $\stackrel{\sim }{m}=2$, 3, 4, and 5 in the top panel, and to λ = 1.1, 1.2, 1.3, and 1.5 in bottom panel. The top panel has a fixed value of ζ = 1.5 and the bottom panel has a fixed value of $\stackrel{\sim }{m}=3$.

We also note that the above analysis is valid for any cooling function T(t). For example, if we take

$$\begin{array}{c}\hfill T(t)=T_0(1-\alpha (sin(t/t_{\mathrm{cool}}\pi )),\end{array}$$(83)

where

$$\begin{array}{c}\hfill \alpha <1,\end{array}$$(84)

the loop is first cooled and then heated. The parameter α in Eq. (84) is responsible for the strength of cooling and heating in such temperature distributions. Therefore, as a result, it is possible to rewrite Eqs. (74) and (80) as

$$\begin{array}{cc}& \frac{{\partial }^2{S}_0}{\partial Z^2}+\frac{{\varpi }^2{\mathrm{\Lambda }}^4(Z)S_0}{4(\zeta +1)}[\zeta exp(-\kappa {(1-\alpha sin(\tau \pi ))}^{-1}cos(\pi Z/2))\hfill \\ \hfill & +exp(-\kappa cos(\pi Z/2))]=0\hfill \end{array}$$(85)