© ESO 2019

1. Introduction

Kink oscillations of coronal magnetic loops were among the first magnetohydrodynamic oscillations observed by space missions in the solar atmosphere at the end of the last century. They were first observed by the Transition Region and Coronal Explorer (TRACE) mission in 1998 and reported by Aschwanden et al. (1999) and Nakariakov et al. (1999). Since then these oscillations have been routinely observed by space missions (e.g. Erdélyi & Taroyan 2008; Duckenfield et al. 2018; Su et al. 2018; Abedini 2018, and references therein). Recently, kink oscillations have also been observed in prominence threads (e.g. Arregui et al. 2018).

Even the first observation of a coronal loop kink oscillation showed that it was strongly damped, with a damping time of a few oscillation periods. At present, it is generally accepted that this damping is caused by resonant absorption. Resonant absorption was first suggested for complimentary plasma heating in fusion devises (Tataronis & Grossmann 1973; Grossmann & Tataronis 1973). At the same time resonant absorption was used to explain the mechanism of magnetospheric pulsations (Lanzerotti et al. 1973; Southwood 1974; Chen & Hasegawa 1974a,b).

After Ionson (1978) proposed resonant absorption as a mechanism for heating coronal magnetic loops, it remained very popular in solar physics. At first, the main effort of researchers was to elaborate the idea by Ionson (1978) and study the coronal heating by resonant absorption of magnetohydrodynamic (MHD) waves (e.g. Kuperus et al. 1981; Ionson 1985; Hollweg 1990). Resonant absorption was also used to explain the observed loss of acoustic power in sunspots (e.g. Hollweg 1988; Lou 1990; Stenuit et al. 1993).

The launch of space missions at the end of the last century boosted the interest of theoretical physicists in studying waves in the solar atmosphere in general, and in application of resonant absorption to problems in solar physics in particular. Nakariakov et al. (1999) suggested that the observed damping of coronal loop kink oscillations was caused by resonant absorption. Ruderman & Roberts (2002) elaborated this idea and showed that the observed oscillation period and damping time can be used to obtain information about the transversal density distribution in a coronal magnetic loop. Goossens et al. (2002) then applied the method suggested by Ruderman & Roberts (2002) to 11 events of coronal loop kink oscillations. Now, this method is among the main methods of a new branch of solar physics called coronal seismology (e.g. Goossens et al. 2006, 2011).

Observations by Solar and Heliospheric Observatory SoHO (Brekke et al. 1997; Winebarger et al. 2002; Teriaca et al. 2004), TRACE (Winebarger et al. 2001; Doyle et al. 2006), and more recently Hinode (Teriaca et al. 2004; Tian et al. 2008; Ofman & Wang 2008) and STEREO (Tian et al. 2009) showed that flows are ubiquitous in active region loops. Flows with speeds of between about 5 km s⁻¹ and 30 km s⁻¹ are also often observed in prominence threads (Chae et al. 2008; Terradas et al. 2008). In view of these observations it is important to study the effect of flow on resonant damping of kink waves in magnetic flux tubes. To our knowledge, Goossens et al. (1992) were the first to address this problem. These latter authors considered kink oscillations of a twisted magnetic tube in the presence of helical flow and derived so-called connection formulae. Using these formulae they studied resonant damping of propagating kink waves in the thin tube and thin boundary layer (TTTB) approximation. Their analysis was extended by Terradas et al. (2010) who studied the damping of propagating kink waves in the TTTB approximation analytically, and without using this approximation numerically. Soler et al. (2011) used the TTTB approximation to study the effect of flow on spatial damping of propagating kink waves analytically, and then studied the same problem numerically without using the TTTB approximation. They compared the analytical and numerical solutions to outline the limit of applicability of the TTTB approximation. Bahari (2018) studied the effect of magnetic twist on the spatial damping of kink waves in the presence of flow.

Terradas et al. (2010) also addressed the problem of resonant damping of standing waves in the presence of flow. However, these latter authors only obtained an approximate solution under the assumption that the flow speed is much smaller than the Alfvén speed in a magnetic tube. They even cast doubt on the possibility of obtaining an exact solution describing a standing wave in the presence of flow. Nevertheless, Ruderman (2010) obtained the solution describing a standing kink wave in a magnetic flux tube in the presence of flow in the thin tube approximation. This solution is a superposition of two waves propagating in different directions. These waves have the same frequency, but different wave numbers. A similar approach was then used by Terradas et al. (2011) who considered an application to coronal seismology. Both Ruderman (2010) and Terradas et al. (2011) considered straight magnetic tubes. Later, Bahari (2017) extended the theory of standing waves in the presence of flow to twisted magnetic tubes. However, all these authors considered a tube with the sharp boundary and thus did not study the resonant wave damping. Therefore, up to now, there have been no studies of the flow effect on resonant damping of standing kink waves that is valid for arbitrary flow velocity. This paper aims to fill this gap.

The paper is organised as follows. In the following section we formulate the problem. The damping rate due to resonant absorption is calculated using the regular perturbation method with the ratio of thickness of the transitional layer to the tube radius used as a small parameter. In Sect. 3 we calculate the first-order approximation. In Sect. 4 we calculate the second-order approximation and thus find the expression for the decrement. In Sect. 5 we apply the general theory to a particular equilibrium with the linear density profile in the transitional region. In Sect. 6 we discuss a possible application of the obtained results to the coronal and prominence seismology. Section 7 contains the summary of obtained results and our conclusions.

2. Problem formulation and governing equations

We consider the plasma motion in the cold plasma approximation. The unperturbed magnetic field is straight in cylindrical coordinates r, ϕ, z, and is given by B = Be_z, where e_z is the unit vector in the z-direction and B is a constant. In the unperturbed state there is also the plasma flow with the velocity U = U(r)e_z. The equilibrium density is given by

$$\begin{array}{c}\hfill \rho (r)=\{\begin{array}{cc}{\rho }_i,\hfill & r\le R(1-l/2),\hfill \\ {\rho }_t(r),\hfill & R(1-l/2)\le r\le R(1+l/2),\hfill \\ {\rho }_e,\hfill & r\ge R(1+l/2),\hfill \end{array}\end{array}$$(1)

where R is the tube radius, ρ_i and ρ_e are constants, ρ_e <  ρ_i, ρ_t(r) is a monotonically decreasing function, and ρ(r) is continuous at r = R(1 ± l/2). The domain defined by r ≤ R(1 − l/2) is the core part of the magnetic tube, while R(1 − l/2)≤r ≤ R(1 + l/2) defines the transitional region. The velocity magnitude is given by

$$\begin{array}{c}\hfill U=\{\begin{array}{cc}{U}_i,\hfill & r\le R(1-l/2),\hfill \\ U_t(r),\hfill & R(1-l/2)\le r\le R(1+l/2),\hfill \\ 0,\hfill & r\ge R(1+l/2),\hfill \end{array}\end{array}$$(2)

where U_t(r) is a monotonically decreasing function, and U(r) is continuous at r = R(1 ± l/2). We can always choose the direction of the z-axis in such a way that U_i >  0. Since the magnetohydrodynamic equations are invariant with respect to the substitution B → −B we can assume that B >  0. The equilibrium configuration is shown in Fig. 1.

Fig. 1. Equilibrium configuration. The magnetic field is the same everywhere, while the velocity is zero outside the tube. The tube length is L.

The tube length is L. We assume that the tube is thin, R ≪ L. In addition, we assume that the transitional layer is also thin, l ≪ 1. Hence, we use the TTTB approximation. Ruderman et al. (2017) derived the equation describing kink oscillations of a thin expanding non-stationary magnetic tube in the presence of flow. For a particular equilibrium considered in this paper this equation reduces to

$$\begin{array}{c}\hfill {(\frac{\partial }{\partial t}+U_i\frac{\partial }{\partial z})}^2\eta +\frac{1}{\zeta }\frac{{\partial }^2\eta }{\partial t^2}-2V_{\mathit{Ai}}^2\frac{{\partial }^2\eta }{\partial z^2}=\mathcal{L},\end{array}$$(3)

$$\begin{array}{c}\hfill \mathcal{L}=\frac{\delta P}{{\rho }_{i}R}+V_{\mathit{Ai}}^2\frac{{\partial }^2(l\eta +\delta \eta )}{\partial z^2}-\frac{1}{\zeta }\frac{{\partial }^2(l\eta +\delta \eta )}{\partial t^2}·\end{array}$$(4)

Here ζ = ρ_i/ρ_e, η, which is the radial displacement of the tube core, is independent of r, $V_{Ai}^2=B^2/{\mu }_0{\rho }_i$, and

$$\begin{array}{cc}\hfill \delta P& =P{|}_{r=R(1+l/2)}-P{|}_{r=R(1-l/2)},\hfill \\ \hfill \delta \eta & ={\xi }_r{|}_{r=R(1+l/2)}-{\xi }_r{|}_{r=R(1-l/2)},\hfill \end{array}$$(5)

where ξ_r is the radial component of the plasma displacement and P the magnetic pressure perturbation. The solution to Eq. (3) must satisfy the boundary conditions

$$\begin{array}{c}\hfill \eta =0\text{at}z=0,L.\end{array}$$(6)

Below, we look for the solution in the form of normal modes and take all dependent variables proportional to exp(−iωt). We then use the regular perturbation method and look for the solution in the form

$$\begin{array}{c}\hfill \eta ={\eta }_1+l{\eta }_2,\omega ={\omega }_1+l{\omega }_2.\end{array}$$(7)

The variation of P and η across the transitional layer is of the order of l. In accordance with this we introduce the scaled quantities

$$\begin{array}{c}\hfill \stackrel{\sim }{\delta P}=l^{-1}\delta P,\stackrel{\sim }{\delta \eta }=l^{-1}\delta \eta .\end{array}$$(8)

We then transform Eqs. (3) and (4) to

$$\begin{array}{c}\hfill {(i\omega -U_i\frac{\mathrm{d}}{\mathrm{d}z})}^2\eta -\frac{{\omega }^2\eta }{\zeta }-2V_{\mathit{Ai}}^2\frac{{\mathrm{d}}^2\eta }{\mathrm{d}{z}^2}=l\stackrel{\sim }{\mathcal{L}},\end{array}$$(9)

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathcal{L}}=\frac{\stackrel{\sim }{\delta P}}{{\rho }_{i}R}+V_{\mathit{Ai}}^2\frac{{\mathrm{d}}^2(\eta +\stackrel{\sim }{\delta \eta })}{\mathrm{d}{z}^2}+\frac{{\omega }^2}{\zeta }(\eta +\stackrel{\sim }{\delta \eta }),\end{array}$$(10)

where we use the ordinary derivative instead of the partial because all quantities now only depend on z.

3. First-order approximation

Collecting terms of the order of unity in Eq. (9) we obtain

$$\begin{array}{c}\hfill (C_k^2-\frac{\zeta U_i^2}{\zeta +1})\frac{{\mathrm{d}}^2{\eta }_1}{\mathrm{d}{z}^2}+\frac{2i{\omega }_1\zeta U_i}{\zeta +1}\frac{\mathrm{d}{\eta }_1}{\mathrm{d}z}+{\omega }_1^2{\eta }_1=0,\end{array}$$(11)

where

$$\begin{array}{c}\hfill C_k^2=\frac{2\zeta V_{\mathit{Ai}}^2}{\zeta +1}·\end{array}$$(12)

The solution to Eq. (11) must satisfy the boundary conditions

$$\begin{array}{c}\hfill {\eta }_1=0\text{at}z=0,L.\end{array}$$(13)

The general solution to Eq. (11) is

$$\begin{array}{c}\hfill {\eta }_1=e^{-ihz}({\eta }_{1+}{e}^{\mathit{iqz}}+{\eta }_{1-}{e}^{-iqz}),\end{array}$$(14)

where

$$\begin{array}{c}\hfill h=\frac{{\omega }_1\zeta U_i}{C_k^2(\zeta +1)-\zeta U_i^2},q=\frac{{\omega }_1\sqrt{{(\zeta +1)}^2{C}_k^2-\zeta U_i^2}}{C_k^2(\zeta +1)-\zeta U_i^2}·\end{array}$$(15)

Substituting Eq. (14) in the boundary conditions Eq. (13) yields

$$\begin{array}{c}\hfill {\eta }_{1+}+{\eta }_{1-}=0,{\eta }_{1+}{e}^{\mathit{iqL}}+{\eta }_{1-}{e}^{-iqL}=0.\end{array}$$(16)

It follows from these equations that

$$\begin{array}{c}\hfill {\eta }_{1-}=-{\eta }_{1+},e^{2iqL}=1.\end{array}$$(17)

Since we are solving a linear problem we can choose η₁₊ arbitrarily. To have the correct dimension we take η₁₊ = −iη₀/2, where η₀ is a constant with the dimension of length. The second equation in Eq. (17) implies that qL = πn, where n = 1, 2, …. Subsequently with the aid of Eq. (15) we obtain

$$\begin{array}{c}\hfill {\omega }_1=\frac{\pi n|C_k^2(\zeta +1)-\zeta U_i^2|}{L\sqrt{{(\zeta +1)}^2{C}_k^2-\zeta U_i^2}}·\end{array}$$(18)

Here, n = 1 corresponds to the fundamental mode, and n >  1 to the (n − 1)th overtone. Since we use the thin tube approximation we must restrict our analysis to moderate values of n satisfying the condition n ≤ n₀, where n₀R ≪ L. When $U_i^2>{(\zeta +1)}^2{C}_k^2/\zeta$ the tube is subject to the Kelvin–Helmholtz (KH) instability with respect to kink oscillations. Below we assume that $U_i^2<(\zeta +1)C_k^2/\zeta$ which, in particular, guarantees that kink oscillations are KH-stable. It is then straightforward to show that ω₁ is a monotonically decreasing function of U_i.

Using Eq. (17) we obtain

$$\begin{array}{c}\hfill {\eta }_1={\eta }_0{e}^{-ihz}sin\frac{\pi nz}{L}·\end{array}$$(19)

4. Second-order approximation

Collecting terms of the order of l in Eq. (9) yields

$$\begin{array}{cc}& (C_k^2-\frac{\zeta U_i^2}{\zeta +1})\frac{{\mathrm{d}}^2{\eta }_2}{\mathrm{d}{z}^2}+\frac{2i{\omega }_1\zeta U_i}{\zeta +1}\frac{\mathrm{d}{\eta }_2}{\mathrm{d}z}+{\omega }_1^2{\eta }_2\hfill \\ \hfill & =-\frac{\zeta \stackrel{\sim }{\mathcal{L}}}{\zeta +1}-2{\omega }_2({\omega }_1{\eta }_1+\frac{i\zeta U_i}{\zeta +1}\frac{\mathrm{d}{\eta }_1}{\mathrm{d}z}).\hfill \end{array}$$(20)

We aim to calculate ω₂. For this we need to obtain the expression for $\stackrel{\sim }{\mathcal{L}}$.

4.1. Calculation of $\stackrel{\sim }{\mathcal{L}}$

To obtain the expression for $\stackrel{\sim }{\mathcal{L}}$ we need to calculate $\stackrel{\sim }{\delta P}$ and $\stackrel{\sim }{\delta \eta }$. For this we use the linearised MHD equations. The plasma is assumed to be weakly dissipative. We see below that there is at least one resonant surface in the transitional layer where the frequency of kink oscillation coincides with the Doppler-shifted local Alfvénic frequency. Near this surface there are large gradients of perturbations. This implies that, although we can neglect dissipation almost everywhere, we have to take it into account near the resonant surface. We need dissipation to remove the singularity at the resonant surface. Below we only need the jumps of the radial plasma displacement and the magnetic pressure perturbation across the dissipative layer. The exact form of dissipation and the value of dissipation coefficients only affect the behaviour of solution in the dissipative layer. However, they do not affect the jumps of the radial plasma displacement and the magnetic pressure perturbation (e.g. Ruderman et al. 1995; Goossens et al. 2011). In particular, these jumps would remain the same even if we account for viscosity and neglect resistivity, or account for resistivity and neglect viscosity, or keep both of them. This observation enables us to keep viscosity and neglect resistivity.

Large gradients of perturbations only exist in the radial direction. As a result, we can write the term describing viscosity in the momentum equation in a simplified form as the second derivative with respect to r of the velocity perturbation. Summarising, we write the MHD equations as

$$\begin{array}{c}\hfill \rho (\frac{\partial \mathit{v}}{\partial t}+(\mathit{U}·\mathrm{\nabla })\mathit{v}+(\mathit{v}·\mathrm{\nabla })\mathit{U})=-\mathrm{\nabla }P+\frac{1}{{\mu }_0}(\mathit{B}·\mathrm{\nabla })\mathit{b}+\frac{1}{{\mu }_0}(\mathit{b}·\mathrm{\nabla })\mathit{B}+\rho \nu \frac{{\partial }^2\mathit{v}}{\partial r^2},\end{array}$$(21)

$$\begin{array}{c}\hfill \frac{\partial \mathit{b}}{\partial t}+(\mathit{U}·\mathrm{\nabla })\mathit{b}+(\mathit{v}·\mathrm{\nabla })\mathit{B}=(\mathit{B}·\mathrm{\nabla })\mathit{v}+(\mathit{b}·\mathrm{\nabla })\mathit{U}-\mathit{B}\mathrm{\nabla }·\mathit{v}.\end{array}$$(22)

Here, ρ = ρ_t(r), U = U_t(r)e_z, v = (v_r, v_ϕ, v_z) is the velocity perturbation, b = (b_r, b_ϕ, b_z) the magnetic field perturbation, P = Bb_z/μ₀ the magnetic pressure perturbation, ν the kinematic viscosity, and μ₀ the magnetic permeability of free space. We do not write the mass conservation equation because it is not used below.

We now introduce the plasma displacement ξ = (ξ_r, ξ_ϕ, ξ_z). Its relation with the velocity perturbation can be written in two equivalent forms (e.g. Goossens et al. 1992; Ruderman et al. 2017),

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

where we take into account that ∇ ⋅ U = 0. Ruderman et al. (2017) used this expression for v in terms of ξ to derive the system of equations for P, ξ_r, and ξ_ϕ. In the particular case of a non-expanding tube, this system reduces to

$$\begin{array}{c}\hfill P=-\rho V_A^2(\frac{1}{r}\frac{\partial (r{\xi }_r)}{\partial r}+\frac{1}{r}\frac{\partial {\xi }_{\varphi }}{\partial \varphi }),\end{array}$$(24)

$$\begin{array}{c}\hfill {(\frac{\partial }{\partial t}+U\frac{\partial }{\partial z})}^2{\xi }_r=-\frac{1}{\rho }\frac{\partial P}{\partial r}+V_A^2\frac{{\partial }^2{\xi }_r}{\partial z^2}+\nu (\frac{\partial }{\partial t}+U\frac{\partial }{\partial z})\frac{{\partial }^2{\xi }_r}{\partial r^2},\end{array}$$(25)

$$\begin{array}{c}\hfill {(\frac{\partial }{\partial t}+U\frac{\partial }{\partial z})}^2{\xi }_{\varphi }=-\frac{1}{\rho r}\frac{\partial P}{\partial \varphi }+V_A^2\frac{{\partial }^2{\xi }_{\varphi }}{\partial z^2}+\nu (\frac{\partial }{\partial t}+U\frac{\partial }{\partial z})\frac{{\partial }^2{\xi }_{\varphi }}{\partial r^2}·\end{array}$$(26)

The characteristic scales in the r and z directions are R and L, respectively. The characteristic time is L/V_Ai. Using Eqs. (25) and (26) we then obtain the estimate that $P/\rho V_A^2$ is of the order of ϵ²ξ_r/R, where ϵ = R/L. This estimate implies that the ratio of the left-hand side of Eq. (24) to its right-hand side is of the order of ϵ², meaning that the left-hand side of Eq. (24) can be neglected and it can be reduced to

$$\begin{array}{c}\hfill \frac{\partial (r{\xi }_r)}{\partial r}+\frac{\partial {\xi }_{\varphi }}{\partial \varphi }=0.\end{array}$$(27)

Now, taking the dependent variables proportional to exp(−iωt + iϕ) and eliminating ξ_ϕ we obtain the system of equations for ξ_r and P,

$$\begin{array}{c}\hfill \frac{1}{\rho }\frac{\partial P}{\partial r}=V_A^2\frac{{\partial }^2{\xi }_r}{\partial z^2}-{(U\frac{\partial }{\partial z}-i\omega )}^2{\xi }_r+\nu (U\frac{\partial }{\partial z}-i\omega )\frac{{\partial }^2{\xi }_r}{\partial r^2},\end{array}$$(28)

$$\begin{array}{c}\hfill P=r\rho [V_A^2\frac{{\partial }^2}{\partial z^2}-{(U\frac{\partial }{\partial z}-i\omega )}^2+\nu (U\frac{\partial }{\partial z}-i\omega )\frac{{\partial }^2}{\partial r^2}]\frac{\partial (r{\xi }_r)}{\partial r}·\end{array}$$(29)

The variation of P in the radial direction in the transitional layer is of the order of lP. We need to calculate the jumps of ξ_r and P across the transitional layer in the leading-order approximation with respect to l. In accordance with this we can take P = P_i in Eq. (29), where P_i is the value of P at r = R(1 − l/2). We also take ω = ω₁. We note that Eqs. (28) and (29) are valid everywhere and not only in the transitional layer. Since in the core region ξ_r = η(z) and we can take ν = 0 there, it follows from Eq. (29) that in this region

$$\begin{array}{c}\hfill P=r{\rho }_i[V_{\mathit{Ai}}^2\frac{{\mathrm{d}}^2\eta }{\mathrm{d}{z}^2}-(U_i\frac{\mathrm{d}}{\mathrm{d}z}-i\omega {)}^2\eta ].\end{array}$$(30)

Now, taking into account that R(1 − l/2)≈R in the leading-order approximation with respect to l, we obtain from Eq. (30)

$$\begin{array}{c}\hfill P_i=R{\rho }_i[V_{\mathit{Ai}}^2\frac{{\mathrm{d}}^2{\eta }_1}{\mathrm{d}{z}^2}-{(U_i\frac{\mathrm{d}}{\mathrm{d}z}-i{\omega }_1)}^2{\eta }_1].\end{array}$$(31)

Substituting this result in Eq. (29) and using the approximation r ≈ R in the transitional layer yields

$$\begin{array}{cc}& [V_A^2\frac{{\partial }^2}{\partial z^2}-{(U\frac{\partial }{\partial z}-i\omega )}^2+\nu (U\frac{\partial }{\partial z}-i{\omega }_1)\frac{{\partial }^2}{\partial r^2}]\frac{\partial (r{\xi }_r)}{\partial r}=\hfill \\ \hfill & \frac{{\rho }_i}{{\rho }_t(r)}[V_{\mathit{Ai}}^2\frac{{\mathrm{d}}^2{\eta }_1}{\mathrm{d}{z}^2}-{(U_i\frac{\mathrm{d}}{\mathrm{d}z}-i{\omega }_1)}^2{\eta }_1].\hfill \end{array}$$(32)

We did not substitute ω₁ for ω in the second term in brackets on the left-hand side of this equation because, as we see further below, the small correction to ω₁ can be important in the dissipative layer. Using Eq. (19) and the variable substitution

$$\begin{array}{c}\hfill \frac{\partial (r{\xi }_r)}{\partial r}=Xexp\left(\frac{i\omega Uz}{U^2-V_A^2}\right),\end{array}$$(33)

we transform Eq. (32) to

$$\begin{array}{cc}& {\sigma }^2\frac{{\partial }^{2}X}{\partial z^2}+{\omega }^{2}X+\nu (\frac{U(V_A^2-U^2)}{V_A^2}\frac{\partial }{\partial z}-i{\omega }_1)\frac{{\partial }^{2}X}{\partial r^2}=\hfill \\ \hfill & \frac{{\eta }_0{e}^{\mathit{igz}}(U^2-V_A^2)}{V_{\mathit{Ai}}^2}(i{C}_{1}cos\frac{\pi nz}{L}+C_{2}sin\frac{\pi nz}{L}),\hfill \end{array}$$(34)

where

$$\begin{array}{c}\hfill \sigma =\frac{V_A^2-U^2}{V_A},g=\frac{{\omega }_{1}U}{V_A^2-U^2}-h,\end{array}$$(35)

$$\begin{array}{c}\hfill C_1=\frac{2\pi n}{L}[h(V_{\mathit{Ai}}^2-U_i^2)-U_i{\omega }_1],\end{array}$$(36)

$$\begin{array}{c}\hfill C_2=(h^2+\frac{{\pi }^2{n}^2}{L^2})(V_{\mathit{Ai}}^2-U_i^2)-{\omega }_1({\omega }_1+2h{U}_i).\end{array}$$(37)

When deriving Eq. (34), we neglected terms proportional to the product of ν and derivatives of unperturbed quantities. The solution to Eq. (34) must satisfy the boundary conditions

$$\begin{array}{c}\hfill X=0\text{at}z=0,L.\end{array}$$(38)

4.1.1. Solution far from the resonant surface

Usually the observed flows in coronal loops and prominence threads are sub-Alfvénic. In accordance with this, below we assume that U_i <  V_Ai. Since ρ_t(r) is a monotonically decreasing function, it follows that V_A(r) is monotonically increasing. Since, in addition, U_t(r) is monotonically decreasing it follows from our assumption that U(r) < V_A(r) in the transitional layer.

Now, we consider the solution sufficiently far from the resonant layer. This enables us to neglect the term proportional to ν in Eq. (34) and reduce it to

$$\begin{array}{c}\hfill {\sigma }^2\frac{{\partial }^{2}X}{\partial z^2}+{\omega }^{2}X=\frac{{\eta }_0{e}^{\mathit{igz}}(U^2-V_A^2)}{V_{\mathit{Ai}}^2}(i{C}_{1}cos\frac{\pi nz}{L}+C_{2}sin\frac{\pi nz}{L}).\end{array}$$(39)

Now we consider the differential operator ℳ defined by the differential expression −σ²∂²X/∂z², and the boundary conditions Eq. (38). It is easy to show that the eigenvalues of this operator are equal to ${\mathrm{\Omega }}_n^2$, where

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_n(r)=\frac{\pi n\sigma (r)}{L},n=1,2,\dots \end{array}$$(40)

When r varies from R(1 − l/2) to R(1 + l/2), Ω_n(r) monotonically increases from nΩ_i to nΩ_e, where

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_i=\frac{\pi (V_{\mathit{Ai}}^2-U_i^2)}{L{V}_{\mathit{Ai}}},{\mathrm{\Omega }}_e=\frac{\pi V_{\mathit{Ae}}}{L},\end{array}$$(41)

taking into account that, by assumption, U(R(1 + l/2)) = 0. Hence, the operator defined by the left-hand side of Eq. (39) and the boundary conditions, Eq. (38) is invertible if and only if ω ≠ ±Ω_n(r) for all r ∈ [R(1 − l/2),R(1 + l/2)] and all n = 1, 2, …. This implies that the Alfvén continuum consists of the infinite union of intervals,

$$\begin{array}{c}\hfill \mathfrak{V}={\displaystyle \bigcup _{n=1}^{\infty }([-n{\mathrm{\Omega }}_e,-n{\mathrm{\Omega }}_i]\bigcup [n{\mathrm{\Omega }}_i,n{\mathrm{\Omega }}_e]).}\end{array}$$(42)

The operator defined by the left-hand side of Eq. (39) and the boundary conditions of Eq. (38) can be written as ω² − ℳ. When ω₁ ∉ 𝔙 the solution to the equation (ω² − ℳ)X = f(z) is

$$\begin{array}{cc}\hfill X& ={({\omega }^2-\mathcal{M})}^{-1}f\hfill \\ \hfill & =-\frac{1}{{\sigma }^2}(\frac{sin[\kappa (L-z)]}{\kappa sin(\kappa L)}{\displaystyle {\int }_0^{z}f(u)sin(\kappa u)\mathrm{d}u}\hfill \\ \hfill & -\frac{sin(\kappa z)}{\kappa sin(\kappa L)}{\displaystyle {\int }_z^{L}f(u)sin[(\kappa (L-u)]\mathrm{d}u),}\hfill \end{array}$$(43)

where κ = ω/σ. Substituting f(z) equal to the right-hand side of Eq. (39) we obtain

$$\begin{array}{cc}\hfill X& =\frac{i{\eta }_0{V}_A}{\sigma V_{\mathit{Ai}}^{2}sin(\kappa L)}{\displaystyle \sum _{j=1}^2\frac{A_j}{{\kappa }^2-{\kappa }_j^2}\{sin(\kappa L)[cos(\kappa z)}\hfill \\ \hfill & -cos({\kappa }_{j}z)]-sin(\kappa z)[cos(\kappa L)-cos({\kappa }_{j}L)]\hfill \\ \hfill & +i[sin(\kappa z)sin({\kappa }_{j}L)-sin(\kappa L)sin({\kappa }_{j}z)]\},\hfill \end{array}$$(44)

where

$$\begin{array}{c}\hfill {\kappa }_j=g+\frac{{(-1)}^j\pi n}{L},A_j=\frac{C_1}{2}-\frac{{(-1)}^j{C}_2}{2}·\end{array}$$(45)

This result can be obtained directly solving the ordinary differential equation, Eq. (39), with the boundary conditions of Eq. (38) where r is only present as a parameter.

It is straightforward to see that the expression for X is regular when ${\kappa }^2={\kappa }_j^2$. Hence, the only singularity in this expression appears when sin(κL) = 0. This condition is equivalent to ω = πnσ(r)/L for some positive integer n. We write the second term in the expansion for ω with respect to l as the sum of the real and imaginary part, ω₂ = ω_2r − iγ. We see further below that γ >  0. This implies that the condition ω = πnσ(r)/L is never satisfied. However, sin(kL) = 𝒪(l⁻¹) when ω₁ = πnσ(r)/L, which implies that there are large gradients in the vicinity of point r_A determined by this equation. Hence, in the vicinity of r_A we have to take into account the term proportional to ν on the left-hand side of Eq. (34). The surface defined by the equation r = r_A is called the Alfvén resonant surface.

Since Ω_i = πσ_i/L and Ω_e = πσ_e/L, where σ_i and σ_e are values of function σ(r) at r = R(1 − l/2) and r = R(1 + l/2), respectively, the condition that there is r_A satisfying the equation ω₁ = πnσ(r_A)/L is equivalent to ω₁ ∈ 𝔙. It is easy to show that the condition nΩ_i <  ω₁ <  nΩ_e is always satisfied when U_i <  V_Ai and ζ >  1. Hence, ω₁ ∈ [nΩ_i, nΩ_e] and, consequently, there is at least one resonant surface. Since the intervals constituting the Alfvén continuum can overlap, it is possible that there are a few resonant surfaces.

4.1.2. Connection formulae

We see that there are large gradients in the radial direction at the resonant surface in the solution to ideal MHD equations. Hence, we need to solve dissipative MHD equations in a narrow dissipative layer embracing the resonance surface. To obtain the solution in the transitional layer we match the ideal solution far from the resonant surface and the dissipative solution in two overlap regions to the left and right of the resonant surface. In these overlap regions, both ideal and dissipative solutions are valid. The only quantities related to the dissipative solution that we need to connect to solutions to the ideal MHD equations are the jumps of ξ_r and P across the dissipative layer. The expressions for these jumps are called the connection formulae. The notion of connection formulae was first introduced by Sakurai et al. (1991). Goossens et al. (1992) derived the connection formulae for propagating waves in the presence of flow, and later Goossens et al. (1995) improved the derivation by Sakurai et al. (1991) and expressed the solution in the dissipative layer in terms of F and G functions.

Here it is worth making an important comment. When there is no flow, we can take all dependent variables proportional to sin(πnz/L). As a result, after that we obtain a 1D problem with all variables only depending on r. Probably the most important point is that Eq. (34) becomes an ordinary differential equation. However, when there is flow the situation is more complex. We managed to construct the solution to Eq. (11) in the form of a standing wave as the superposition of two propagating waves with the same frequency and different wave numbers. However this approach does not work for Eq. (34). It is even not clear if it is possible to obtain the solution to Eq. (34) as a sum of a few functions with prescribed dependence on z.

A natural way to obtain the solution to Eq. (34) is to expand X in the Fourier series:

$$\begin{array}{c}\hfill X(r,z)={\displaystyle \sum _{m=1}^{\infty }{X}_m(r)sin\frac{\pi mz}{L}·}\end{array}$$(46)

We also use the expansions

$$\begin{array}{cc}\hfill \frac{\partial X}{\partial z}& =\frac{8}{L}\{{\displaystyle \sum _{m=1}^{\infty }(2m-1)sin\frac{\pi (2m-1)z}{L}}\hfill \\ \hfill & \times {\displaystyle \sum _{k=1}^{\infty }\frac{k{X}_{2k}(r)}{{(2m-1)}^2-4k^2}}\hfill \\ \hfill & +{\displaystyle \sum _{m=1}^{\infty }msin\frac{2\pi mz}{L}\sum _{k=1}^{\infty }\frac{(2k-1)X_{2k-1}(r)}{4m^2-{(2k-1)}^2}\},}\hfill \end{array}$$(47)

$$\begin{array}{c}\hfill e^{\mathit{igz}}(i{C}_{1}cos\frac{\pi nz}{L}+C_{2}sin\frac{\pi nz}{L})=i{\displaystyle \sum _{m=1}^{\infty }{a}_m(r)sin\frac{\pi mz}{L},}\end{array}$$(48)

where

$$\begin{array}{c}\hfill a_m(r)=[{(-1)}^{n+m}{e}^{\mathit{igL}}-1](\frac{\pi (m-n)C_1-gL{C}_2}{g^2{L}^2-{\pi }^2{(n-m)}^2}+\frac{\pi (m+n)C_1-gL{C}_2}{g^2{L}^2-{\pi }^2{(n+m)}^2}).\end{array}$$(49)

Substituting Eqs. (46) and (47) in Eq. (34) we obtain with the aid of Eq. (48)

$$\begin{array}{c}\hfill (\frac{{\pi }^2{m}^2{\sigma }^2}{L^2}-{\omega }^2)X_m-\frac{8\nu U(V_A^2-U^2)}{L{V}_A^2}\frac{{\mathrm{d}}^2{Y}_m}{\mathrm{d}{r}^2}+i{\omega }_1\nu \frac{{\mathrm{d}}^2{X}_m}{\mathrm{d}{r}^2}=\frac{i{\eta }_0{a}_m(V_A^2-U^2)}{V_{\mathit{Ai}}^2},\end{array}$$(50)

where

$$\begin{array}{c}\hfill Y_m(r)=\{\begin{array}{cc}{\displaystyle {\displaystyle {\sum }_{k=1}^{\infty }\frac{k(2m-1)X_{2k}}{{(2m-1)}^2-4k^2},}}\hfill & m\text{odd},\hfill \\ {\displaystyle {\displaystyle {\sum }_{k=1}^{\infty }\frac{m(2k-1)X_{2k-1}}{4m^2-{(2k-1)}^2},}}\hfill & m\text{even}.\hfill \end{array}\end{array}$$(51)

If ω₁ ≠ πmσ(r)/L in the whole transitional layer then the first term on the left-hand side of this equation strongly dominates the terms proportional to ν, which can be neglected. As a result, X_m(r) is given by

$$\begin{array}{c}\hfill X_m(r)=\frac{i{\eta }_0{a}_m{L}^2(V_A^2-U^2)}{V_{\mathit{Ai}}^2({\pi }^2{m}^2{\sigma }^2-L^2{\omega }_1^2)},\end{array}$$(52)

where we substituted ω₁ for ω because we calculate all quantities in the transitional layer in the leading order approximation with respect to l.

Now we solve Eq. (50) for values of m corresponding to the existence of resonant surfaces. As we have already shown, the resonant layer always exists when m = n. We also point out that there is a possibility of existence of more than one resonant surface. Since V_A(r) monotonically increases and U(r) monotonically decreases in the transitional layer, it follows that σ(r) monotonically increases. It then becomes apparent that if there is a resonant surface for m = n₂ >  n, then there is a resonant surface for any m ∈ [n, n₂]. Similarly, if there is a resonant surface for m = n₁ <  n (in the case where n >  1), then there is a resonant surface for any m ∈ [n₁, n]. Where there are resonant surfaces for m ∈ [n₁, n₂], the mth resonant surface is defined by r = r_m, where r_m is determined by the equation

$$\begin{array}{c}\hfill \sigma (r_m)=\frac{L{\omega }_1}{\pi m},m=n_1,n_1+1,\dots n_2.\end{array}$$(53)

Since σ(r) is a monotonically increasing function, it follows that r_n2 <  r_n2 − 1 <  … <  r_n1.

Let m ∈ [n₁, n₂]. First of all, we notice that the expression for Y_m does not contain X_m. In the vicinity of r = r_m the function X_k(r) with k ≠ m is given by Eq. (52) and its gradient is much smaller than the gradient of X_m(r). This implies that the third term on the left-hand side of Eq. (50) strongly dominates the second term, and thus the second term can be neglected. Now, following Goossens et al. (1995) we make the variable substitution s = r − r_m in Eq. (50). Since the dissipative layer is very thin, we can use the approximate relation

$$\begin{array}{c}\hfill \frac{{\pi }^2{m}^2{\sigma }^2}{L^2}-{\omega }^2=\mathrm{\Delta }s+2il{\omega }_1\gamma ,\mathrm{\Delta }=\frac{{\pi }^2{m}^2}{L^2}\frac{\mathrm{d}{\sigma }^2}{\mathrm{d}r}{|}_{s=0},\end{array}$$(54)

where ω₂ = ω_2r − iγ. We neglect the term −2ω₁ω_2r because, in principle, it can be included in the definition of r_m causing a small shift in the position of the resonant surface. It practically does not change the solution in the dissipative layer. In contrast, the account of the term proportional to γ can substantially change this solution as was first shown by Ruderman et al. (1995); see also Tirry & Goossens (1996). We can also substitute r_m for r in the right-hand side of Eq. (50). Again following Goossens et al. (1995) we introduce the dimensionless variable

$$\begin{array}{c}\hfill \tau =\frac{s}{{\delta }_A},{\delta }_A={\left(\frac{{\omega }_1\nu }{\mathrm{\Delta }}\right)}^{1/3}.\end{array}$$(55)

Using this variable and Eq. (54) we transform Eq. (50) to

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{X}_m}{\mathrm{d}{\tau }^2}+(\mathrm{\Lambda }-i\tau )X_m=\frac{{\eta }_0{a}_m}{{\delta }_A\mathrm{\Delta }{V}_{\mathit{Ai}}^2}(V_A^2-U^2),\end{array}$$(56)

where

$$\begin{array}{c}\hfill \mathrm{\Lambda }=\frac{2l{\omega }_1\gamma }{{\delta }_A\mathrm{\Delta }}·\end{array}$$(57)

The solution to this equation was obtained by Tirry & Goossens (1996); see also Goossens et al. (2011); it reads

$$\begin{array}{c}\hfill X_m=-\frac{{\eta }_0{a}_m}{{\delta }_A\mathrm{\Delta }{V}_{\mathit{Ai}}^2}(V_A^2-U^2)F_{\mathrm{\Lambda }}(\tau ),\end{array}$$(58)

where

$$\begin{array}{c}\hfill F_{\mathrm{\Lambda }}(\tau )={\displaystyle {\int }_0^{\infty }exp(-iu\tau -u^3/3+\mathrm{\Lambda }u)\mathrm{d}u.}\end{array}$$(59)

We now use Eq. (33). In the dissipative layer we can neglect the derivatives of unperturbed quantities. In addition, we substitute R for r in the transitional layer. As a result, using Eq. (46) we transform Eq. (33) to

$$\begin{array}{c}\hfill \frac{\partial {\xi }_r}{\partial r}=\frac{1}{R}exp\left(\frac{i\omega Uz}{U^2-V_A^2}\right){\displaystyle \sum _{m=1}^{\infty }{X}_m(r)sin\frac{\pi mz}{L},}\end{array}$$(60)

where U and V_A are calculated at r = r_m. When m ∉ [n₁, n₂] the integral of X_m across the dissipative layer is of the order of the thickness of this layer, which is much smaller than lR. Hence, we can neglect this variation and only leave terms in the sum in Eq. (60) with m ∈ [n₁, n₂]. Moreover, the dissipative layers are separated meaning that we need to leave only one term in this sum when calculating the jump of ξ_r across a particular dissipative layer. Hence, to calculate the jump of ξ_r across the dissipative layer embracing the resonant surface r = r_m, we can use the reduced Eq. (60),

$$\begin{array}{c}\hfill \frac{\partial {\xi }_r}{\partial \tau }=\frac{{\delta }_A}{R}exp\left(\frac{i\omega Uz}{U^2-V_A^2}\right)X_m(\tau )sin\frac{\pi mz}{L}·\end{array}$$(61)

Using Eqs. (58) and (59) we then obtain

$$\begin{array}{cc}\hfill {\xi }_r& =\frac{i{\eta }_0{a}_m(U^2-V_A^2)}{R\mathrm{\Delta }{V}_{\mathit{Ai}}^2}exp\left(\frac{i\omega Uz}{U^2-V_A^2}\right)sin\frac{\pi mz}{L}{G}_{\mathrm{\Lambda }}(\tau )+f_m(z),\hfill \end{array}$$(62)

where f_m(z) is an arbitrary function, and

$$\begin{array}{c}\hfill G_{\mathrm{\Lambda }}(\tau )={\displaystyle {\int }_0^{\infty }\frac{exp(-iu\tau )-1}{u}{e}^{-u^3/3+\mathrm{\Lambda }u}\mathrm{d}u.}\end{array}$$(63)

The functions F_Λ(τ) and G_Λ(τ) were introduced by Goossens et al. (2011). When Λ = 0 they coincide with F(τ) and G(τ), respectively. We note that the expression for G_Λ(τ) is slightly different from that given by Goossens et al. (2011). However, the two expressions are different by a constant. Since ξ_r is defined with the accuracy up to an additive function of z, this difference does not affect the expression for ξ_r.

We define the jump of a function of τ across the dissipative layer as

$$\begin{array}{c}\hfill [f(\tau )]=\underset{\tau \to \infty }{lim}\{f(\tau )-f(-\tau )\}.\end{array}$$(64)

To calculate [G_Λ(τ)] we use the variable substitution u′ = uτ and then drop the prime. As a result we obtain

$$\begin{array}{cc}\hfill [G_{\mathrm{\Lambda }}(\tau )]& =-2i\underset{\tau \to \infty }{lim}{\displaystyle {\int }_0^{\infty }\frac{sin(u\tau )}{u}{e}^{-u^3/3+\mathrm{\Lambda }u}\mathrm{d}u}\hfill \\ \hfill & =-2i\underset{\tau \to \infty }{lim}{\displaystyle {\int }_0^{\infty }\frac{sinu}{u}{e}^{-u^3/3{\tau }^3+\mathrm{\Lambda }u/\tau }\mathrm{d}u=-i\pi .}\hfill \end{array}$$(65)

Using this result we obtain from Eq. (62)

$$\begin{array}{c}\hfill {[{\xi }_r]}_m=\frac{\pi {\eta }_0{a}_m(U^2-V_A^2)}{R\mathrm{\Delta }{V}_{\mathit{Ai}}^2}exp\left(\frac{i{\omega }_{1}Uz}{U^2-V_A^2}\right)sin\frac{\pi mz}{L},\end{array}$$(66)

where all unperturbed quantities are calculated at r = r_m. This expression is called the connection formula, first introduced by Sakurai et al. (1991), as stated above. The subscript m indicates that this expression gives the jump of ξ_r across the dissipative layer embracing the resonant surface r = r_m. We note that, when Λ >  1, which is the case for typical coronal conditions, the thickness of the dissipative layer is of the order of l²R, and inside the dissipative layer the behaviour of ξ_r is strongly oscillatory (Goossens et al. 2011).

We now calculate [P]. It is easy to show that the ratio of the last term to the first term on the right-hand side of in Eq. (28) is of the order of δ_A/lR, but δ_A does not exceed the thickness of the dissipative layer, that is δ_A ≲ l²R. It then follows that the last term is much smaller than the first and can be neglected. Subsequently, we obtain that the jump of P across the dissipative layer is of the order of $\rho V_A^2{\eta }_0/L^2$ times the thickness of the dissipative layer, that is, it is of the order of $(\rho V_A^{2}l{\eta }_0/R){(R/L)}^2$. On the other hand, the variation of P across the transitional layer is of the order of $\rho V_A^{2}l{\eta }_0/R$, which is by far larger than the jump of P across the dissipative layer. Hence, we can neglect the jump of P across the dissipative layer and write

$$\begin{array}{c}\hfill {[P]}_m=0.\end{array}$$(67)

A similar result has been obtained in many studies (e.g. Sakurai et al. 1991; Goossens et al. 1995, 2011).

4.1.3. Calculating $\stackrel{\sim }{\delta P}$ and $\stackrel{\sim }{\delta \eta }$

In the method of matched asymptotic expansions, the solution in the dissipative layer is called the internal solution, and the solution outside of the dissipative layer the external solution. Equation (64) defines the jump of a quantity across the dissipative layer from the point of view of the internal solution. We can also define the jump of a quantity across the dissipation layer from the point of view of the external solution:

$$\begin{array}{c}\hfill [f(r){]}_m=\underset{\epsilon \to +0}{lim}\{f(r_m+\epsilon )-f(r_m-\epsilon )\}.\end{array}$$(68)

The two expressions for the jump must coincide.

The transitional layer with the resonant surfaces removed is the union of n₂ − n₁ + 1 intervals, [R(1 − l/2),r_n1), (r_n1, r_n1 + 1),…(r_n2, R(1 + l/2)]. The solution to the ideal MHD equations in the transitional layer is given by

$$\begin{array}{cc}\hfill {\xi }_r& =\frac{1}{R}{\displaystyle {\int }_{R(1-l/2)}^{r}Xexp\left(\frac{i{\omega }_{1}Uz}{U^2-V_A^2}\right)\mathrm{d}r\prime }\hfill \\ \hfill & +{\xi }_r{|}_{r=R(1-l/2)},r\in [R(1-l/2),r_{n_2}],\hfill \end{array}$$(69)

$$\begin{array}{cc}\hfill {\xi }_r& =\frac{1}{R}{\displaystyle {\int }_{{\overline{r}}_m}^{r}Xexp\left(\frac{i{\omega }_{1}Uz}{U^2-V_A^2}\right)\mathrm{d}r\prime +{\overline{\xi }}_{\mathit{rm}}(z),}\hfill \\ \hfill & r\in (r_m,r_{m-1}),\hfill \end{array}$$(70)

$$\begin{array}{cc}\hfill {\xi }_r& =\frac{1}{R}{\displaystyle {\int }_{R(1+l/2)}^{r}Xexp\left(\frac{i{\omega }_{1}Uz}{U^2-V_A^2}\right)\mathrm{d}r\prime }\hfill \\ \hfill & +{\xi }_r{|}_{r=R(1+l/2)},r\in (r_{n_1},R(1+l/2)],\hfill \end{array}$$(71)

where m = n₁ + 1, n₁ + 2, …, n₂, X(r, z) is given by Eq. (46) with X_m(r) defined by Eq. (52), ${\overline{\xi }}_{\mathit{rm}}(z)$ is an arbitrary function, and ${\overline{r}}_m\in (r_m,r_{m-1})$, but otherwise it is arbitrary. Since we calculate ξ_r in the transitional layer in the leading-order approximation with respect to l, we substituted R for r and ω₁ for ω. The function X(r, z) is given by Eq. (44) with ω₁ substituted for ω. When n₁ = n₂ = n, that is there is only one resonant surface, the ideal solution in the transitional layer is given by Eqs. (69) and (71), while Eq. (70) is redundant.

Using Eqs. (68)–(71) we obtain

$$\begin{array}{cc}\hfill [{\xi }_r{]}_{n_2}& ={\overline{\xi }}_{r{n}_2}(z)-{\xi }_r{|}_{r=R(1-l/2)}\hfill \\ \hfill & -\frac{1}{R}\mathcal{P}{\displaystyle {\int }_{R(1-l/2)}^{{\overline{r}}_{n_2}}Xexp\left(\frac{i{\omega }_{1}Uz}{U^2-V_A^2}\right)\mathrm{d}r\prime ,}\hfill \end{array}$$(72)

$$\begin{array}{cc}\hfill [{\xi }_r{]}_m& ={\overline{\xi }}_{rm-1}(z)-{\overline{\xi }}_{\mathit{rm}}(z)\hfill \\ \hfill & -\frac{1}{R}\mathcal{P}{\displaystyle {\int }_{{\overline{r}}_m}^{{\overline{r}}_{m-1}}Xexp\left(\frac{i{\omega }_{1}Uz}{U^2-V_A^2}\right)\mathrm{d}r\prime ,}\hfill \end{array}$$(73)

$$\begin{array}{cc}\hfill [{\xi }_r{]}_{n_1}& ={\xi }_r{|}_{r=R(1+l/2)}-{\overline{\xi }}_{r{n}_1}(z)\hfill \\ \hfill & -\frac{1}{R}\mathcal{P}{\displaystyle {\int }_{{\overline{r}}_{n_1-1}}^{R(1+l/2)}Xexp\left(\frac{i{\omega }_{1}Uz}{U^2-V_A^2}\right)\mathrm{d}r\prime ,}\hfill \end{array}$$(74)

where again m = n₁ + 1, n₁ + 2, …, n₂, and X(r, z) is given by Eq. (46) with X_m(r) defined by Eq. (52). When n₂ = n₁ + 1, that is there are only two resonant surfaces, the jumps of the ideal solution through these surfaces are given by Eqs. (72) and (74), while Eq. (73) is redundant. Finally, when n₁ = n₂ = n, that is there is only one resonant surface, the jump of the ideal solution through this surface is given by

$$\begin{array}{cc}\hfill [{\xi }_r{]}_n& ={\xi }_r{|}_{r=R(1+l/2)}-{\xi }_r{|}_{r=R(1-l/2)}\hfill \\ \hfill & -\frac{1}{R}\mathcal{P}{\displaystyle {\int }_{R(1-l/2)}^{R(1+l/2)}Xexp\left(\frac{i{\omega }_{1}Uz}{U^2-V_A^2}\right)\mathrm{d}r\prime .}\hfill \end{array}$$(75)

In Eqs. (72)–(75), 𝒫 indicates the Cauchy principal part of an integral. We use this symbol because X has a singularity at r = r_m.

Equations (72)–(74) represent n₂ − n₁ + 1 equations because in Eq. (73), m runs from n₁ + 1 to n₂ − 1. Taking the sum of all these equations and recalling that the internal and external expressions for the jump of ξ_r through the dissipative layer must coincide, we obtain, with the aid of Eq. (5),

$$\begin{array}{c}\hfill l\stackrel{\sim }{\delta \eta }=\frac{1}{R}\mathcal{P}{\displaystyle {\int }_{R(1-l/2)}^{R(1+l/2)}Xexp\left(\frac{i{\omega }_{1}Uz}{U^2-V_A^2}\right)\mathrm{d}r+\sum _{m=n_1}^{n_2}{[{\xi }_r]}_m,}\end{array}$$(76)

where [ξ_r]_m is defined by Eq. (66).

When calculating [P]_m we showed that the last term on the right-hand side of Eq. (28) can be neglected in the dissipative layer. Since it can also be neglected outside of the dissipative layer, we conclude that this term can be neglected everywhere in the transitional layer. Since the variation of ξ_r in the transitional layer is of the order of lη, and ξ_r = η in the core region, we can substitute η₁ for ξ_r in Eq. (28) in the leading-order approximation with respect to l. We then obtain

$$\begin{array}{c}\hfill \frac{1}{\rho }\frac{\partial P}{\partial r}=V_A^2\frac{{\partial }^2{\eta }_1}{\partial z^2}-{(U\frac{\partial }{\partial z}-i{\omega }_1)}^2{\eta }_1.\end{array}$$(77)

Using Eqs. (19) and (67) we obtain from this equation in the leading-order approximation with respect to l:

$$\begin{array}{c}\hfill l\stackrel{\sim }{\delta P}=-{\eta }_0{e}^{-ihz}(i{C}_{1}cos\frac{\pi nz}{L}+C_{2}sin\frac{\pi nz}{L}){\displaystyle {\int }_{R(1-l/2)}^{R(1+l/2)}\rho \mathrm{d}r.}\end{array}$$(78)

Now, using Eqs. (19), (66), (76), and (78) we obtain

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathcal{L}}={\stackrel{\sim }{\mathcal{L}}}_1+{\stackrel{\sim }{\mathcal{L}}}_2+{\stackrel{\sim }{\mathcal{L}}}_3,\end{array}$$(79)

where

$$\begin{array}{cc}\hfill {\stackrel{\sim }{\mathcal{L}}}_1& =(V_{\mathit{Ai}}^2\frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}+\frac{{\omega }_1^2}{\zeta }){\displaystyle \sum _{m=n_1}^{n_2}{[{\xi }_r]}_m}\hfill \\ \hfill & =\frac{\pi {\eta }_0}{\mathit{lR}}{\displaystyle \sum _{m=n_1}^{n_2}\frac{a_m(V_{\mathit{Am}}^2-U_m^2)}{{\mathrm{\Delta }}_m{V}_{\mathit{Ai}}^2}exp\left(\frac{i{\omega }_1{U}_{m}z}{U_m^2-V_{\mathit{Am}}^2}\right)}\hfill \\ \hfill & \times (M_{m}sin\frac{\pi mz}{L}+i{N}_{m}cos\frac{\pi mz}{L}),\hfill \end{array}$$(80)

$$\begin{array}{cc}\hfill {\stackrel{\sim }{\mathcal{L}}}_2& =-{\eta }_0{e}^{-ihz}[(S_{n}sin\frac{\pi nz}{L}+\frac{2\pi inh{V}_{\mathit{Ai}}^2}{L}cos\frac{\pi nz}{L})\hfill \\ \hfill & +\frac{1}{\mathit{lR}}(i{C}_{1}cos\frac{\pi nz}{L}+C_{2}sin\frac{\pi nz}{L}){\displaystyle {\int }_{R(1-l/2)}^{R(1+l/2)}\frac{\rho }{{\rho }_i}\mathrm{d}r],}\hfill \end{array}$$(81)

$$\begin{array}{cc}\hfill {\stackrel{\sim }{\mathcal{L}}}_3& =-\frac{i{\eta }_0{L}^2}{lR{V}_{\mathit{Ai}}^2}\mathcal{P}{\displaystyle {\int }_{R(1-l/2)}^{R(1+l/2)}(V_A^2-U^2)exp\left(\frac{i{\omega }_{1}Uz}{U^2-V_A^2}\right)}\hfill \\ \hfill & \times {\displaystyle \sum _{m=1}^{\infty }\frac{a_m}{{\pi }^2{m}^2{\sigma }^2-L^2{\omega }_1^2}(T_{m}sin\frac{\pi mz}{L}+i{W}_{m}cos\frac{\pi mz}{L})\mathrm{d}r,}\hfill \end{array}$$(82)