© ESO 2019

1. Introduction

High-resolution observations have shown the ubiquitous presence of nearly incompressible transverse waves propagating in the solar corona (see, e.g., Tomczyk et al. 2007; Tomczyk & McIntosh 2009; McIntosh et al. 2011; Morton et al. 2016). From the theoretical point of view, these waves are interpreted as magnetohydrodynamic (MHD) kink waves in magnetic flux tubes (see, e.g., Erdélyi & Fedun 2007; Van Doorsselaere et al. 2008; Mathioudakis et al. 2013; Jess et al. 2015). In solar coronal conditions, long-wavelength MHD kink waves are almost incompressible and their dominant restoring force is magnetic tension (Goossens et al. 2012). Due to plasma density and/or magnetic field nonuniformity across the waveguide, MHD kink waves undergo the process of resonant absorption in the Alfvén continuum, by which their energy is transferred to localized Alfvén waves that later develop small scales because of phase mixing (see, e.g., Lee & Roberts 1986; Poedts et al. 1989; Pascoe et al. 2012; Goossens et al. 2014; Soler & Terradas 2015). This process naturally brings wave energy down to dissipative scales at which it can be thermalized and heat the coronal plasma (see, e.g., Ionson 1978).

In addition to waves, observations often show the presence of flows along coronal loops (see the review by Reale 2014). Most of the measured Doppler velocities associated to flows are typically in the range ∼5 − 30 km s⁻¹ (e.g., Del Zanna 2008; Winebarger et al. 2013), although velocities of ∼50 km s⁻¹ have also been reported frequently (e.g., Brekke et al. 1997; Doschek et al. 2008), and even larger velocities up to ∼72 − 123 km s⁻¹ have been observed (Ofman & Wang 2008). Considering that the expected value of the Alfvén velocity in the solar corona is ∼1000 km s⁻¹, these observed flow speeds correspond to small fractions of the coronal Alfvén velocity. Observations of flows with Alfvénic speeds (≳500 − 1000 km s⁻¹) are scarce and often related to very energetic or explosive events (e.g., Innes et al. 2003; Nitta et al. 2012).

The presence of flows along the coronal waveguides modifies the behavior and properties of the MHD waves (see, e.g., Nakariakov & Roberts 1995; Terra-Homem et al. 2003; Ofman & Wang 2008; Ruderman 2010, to name a few works). In the case of resonantly damped kink waves, the effect of flow has been studied in a number of papers. Goossens et al. (1992) presented an analytic theory based on the use of jump conditions for the wave perturbations at the Alfvén resonance. Goossens et al. (1992) derived an approximate expression for the damping rate of the waves, which was applicable to the case that the variation of density and flow velocity across the waveguide were strictly confined to a thin layer, in other words, the so-called thin boundary approximation. Terradas et al. (2010) and Soler et al. (2011) use the analytic formalism of Goossens et al. (1992) to study the effect of flow on the resonant damping of standing and propagating kink waves in coronal loops. Terradas et al. (2010) and Soler et al. (2011) also complement their results with the use of fully numerical resistive eigenvalue computations, which show a good agreement with the analytic approximations.

The thin boundary approximation is useful from an analytic point of view, because it allows the derivation of a simple expression for the damping rate. However, there is no observational justification for the use of such an approximation. Indeed, some observations indicate that coronal loops are largely inhomogeneous in the transverse direction (see, e.g., Aschwanden et al. 2003; Goddard et al. 2017), and numerical simulations point out that the waves themselves may contribute to generating wide nonuniform layers owing to nonlinear Kelvin–Helmholtz instabilites (see Goddard et al. 2018). Therefore, it is important to explore the properties of resonant kink waves beyond the thin boundary approximation. Moreover, the role of flows need to be explored beyond the limit of thin transitions.

In the absence of flow the resonant damping of kink waves in loops with thick boundaries is investigated by Van Doorsselaere et al. (2004) and Arregui et al. (2005) using numerical resistive MHD eigenvalue computations. Soler et al. (2013) revisit the same problem in ideal MHD, but with an entirely different semi-analytic approach. The analytic part of the method is based on expressing the wave perturbations in the thick nonuniform boundary of the waveguide as a Frobenius series that contains a singular term accounting for the Alfvén resonance, while the numerical part of the method simply consists of solving a transcendental equation that plays the role of the dispersion relation. Since the technique of Soler et al. (2013) is much faster than the solution of the resistive eigenvalue problem, detailed parameter studies can be tackled. The solutions provided by the method of Soler et al. (2013) have been successfully used in a number of papers: Soler et al. (2014) test the error associated with the use of the thin boundary approximation beyond its theoretical range of applicability; Goossens et al. (2014) discuss the kink quasi-mode displacement field in a tube with a wide boundary; Arregui et al. (2015) perform Bayesian inference, model comparison, and model-averaging techniques to infer the cross-field density structuring in coronal waveguides; Soler & Terradas (2015) compare the solution provided by the Frobenius-based method with the temporal evolution obtained by expressing the kink wave as a superposition of Alfvén continuum modes; and Soler (2017) investigates the behavior of fluting modes in transversely nonuniform tubes. Independently, a similar series expansion approach (but without a resonant term) is used to study sausage waves in nonuniform tubes (see Guo et al. 2016).

The purpose of the present work is to further extend the method of Soler et al. (2013) by incorporating the effect of mass flow along the waveguide. Both the density and the flow velocity are allowed to vary across the flux tube in a nonuniform layer of arbitrary thickness. Section 2 presents the mathematical formalism, which is largely based on that by Soler et al. (2013) but with the appropriate modifications to incorporate the flow. One of the most important differences is that, because of the spatially varying flow velocity, the equation for the radial position of the Alfvén resonance becomes a transcendental equation that has to be solved iteratively along with the dispersion relation. In order to verify the method, approximate results in the thin tube, thin boundary, and slow flow approximations are obtained and compared with previous expressions in the literature. As an application of this method, we explored the effect of flow on the phase velocity and resonant damping length of forward and backward propagating kink waves in coronal tubes with thick boundaries, described in Sect. 3. Finally, some concluding remarks and prospects for future studies are given in Sect. 4.

2. Method

2.1. Background

As the background configuration to represent a coronal waveguide we consider a straight magnetic cylinder of radius R embedded in a uniform and unbounded plasma. Cylindrical coordinates were used, with r, φ, and z denoting the radial, azimuthal, and longitudinal coordinates, respectively. The magnetic field is straight and constant and along the axis of the cylinder, namely B = B1_z. The mass density, ρ, is uniform in the azimuthal and longitudinal directions and nonuniform in the radial direction, so that ρ = ρ(r). We considered the following radial dependence for the density,

$$\begin{array}{c}\hfill \rho (r)=\{\begin{array}{ccc}{\rho }_{\mathrm{i}},\hfill & \text{if}\hfill & r\le R-l/2,\hfill \\ {\rho }_{\mathrm{tr}}(r),\hfill & \text{if}\hfill & R-l/2<r<R+l/2,\hfill \\ {\rho }_{\mathrm{e}},\hfill & \text{if}\hfill & r\ge R+l/2,\hfill \end{array}\end{array}$$(1)

where ρ_i and ρ_e are internal and external constant densities and ρ_tr(r) is the continuous density profile that connects the internal plasma to the external plasma. We considered ρ_i >  ρ_e to represent a tube that is denser than the surrounding plasma. The thickness of the nonuniform boundary layer, l, is arbitrary and can take any value between l = 0 (abrupt jump) and l = 2R (fully inhomogeneous tube).

In addition, we considered the presence of a field-aligned mass flow, namely U = U1_z, where U is the flow velocity. As in the case of the density, we assume the flow velocity to be uniform in the azimuthal and longitudinal directions and nonuniform in the radial direction, so that U = U(r). For simplicity, we assumed that the radial dependence of the flow velocity mimics that of the density, namely

$$\begin{array}{c}\hfill U(r)=\{\begin{array}{ccc}{U}_{\mathrm{i}},\hfill & \text{if}\hfill & r\le R-l/2,\hfill \\ U_{\mathrm{tr}}(r),\hfill & \text{if}\hfill & R-l/2<r<R+l/2,\hfill \\ U_{\mathrm{e}},\hfill & \text{if}\hfill & r\ge R+l/2,\hfill \end{array}\end{array}$$(2)

where U_i and U_e are internal and external flow velocities and U_tr(r) denotes the spatially-dependent flow velocity in the transitional layer.

2.2. Linear perturbations

Linear ideal MHD waves are superimposed on the background state. To study coronal transverse waves, we considered the linearized ideal MHD equations in the β = 0 approximation, where β refers to the ratio of the thermal pressure to the magnetic pressure. Therefore, the basic equations used in the present work are

$$\begin{array}{cc}& \rho {\displaystyle (\frac{\partial \mathbf{v}}{\partial t}+\mathbf{U}·\mathrm{\nabla }\mathbf{v}+\mathbf{v}·\mathrm{\nabla }\mathbf{U})}=\frac{1}{\mu }(\mathrm{\nabla }\times \mathbf{b})\times \mathbf{B},\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathbf{b}}{\partial t}-\mathrm{\nabla }\times (\mathbf{U}\times \mathbf{b})=\mathrm{\nabla }\times (\mathbf{v}\times \mathbf{B}),\hfill \end{array}$$(4)

where v = (v_r, v_φ, v_z) is the velocity perturbation, b = (b_r, b_φ, b_z) is the magnetic field perturbation, and μ is the magnetic permeability. In addition, the plasma Lagrangian displacement, ξ = (ξ_r, ξ_φ, ξ_z), is related to the velocity perturbation and background flow by

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

From here on we adopt the so-called quasi-mode approach, which assumes that the nonuniform waveguide supports global modes (see, e.g., Goossens et al. 2014, for a discussion on the validity of this approach). We expressed the temporal dependence of perturbations as exp(−iωt), where ω the global mode frequency. In addition, we Fourier-analyzed the perturbations along the uniform φ- and z-directions, so that the perturbations are put proportional to exp(imφ + ik_zz), where m and k_z and the azimuthal and longitudinal wavenumbers, respectively. In the linear regime, different azimuthal and longitudinal wavenumbers do not interact with each other. Hence, we only retained the dependence of the perturbations on the radial direction. We find from Eq. (5) that the relation between the components of the velocity perturbation and the Lagrangian displacement are

$$\begin{array}{cc}\hfill v_{\mathrm{r}}& =-i\mathrm{\Omega }(r){\xi }_{\mathrm{r}},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill v_{\phi }& =-i\mathrm{\Omega }(r){\xi }_{\phi },\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill v_z& =-\frac{\partial U(r)}{\partial r}{\xi }_{\mathrm{r}},\hfill \end{array}$$(8)

where Ω(r)=ω − k_zU(r) is the spatially-dependent Doppler-shifted frequency. We note that because of the presence of the spatially-dependent longitudinal flow, v_z is nonzero even in the β = 0 approximation.

The following procedure closely follows that of Soler et al. (2013), with the difference here we have included the effect of flow. We used the total pressure Eulerian perturbation, P′=B ⋅ b/μ, as our main variable. We combined Eqs. (3) and (4) and, after some algebraic manipulations, we obtain a differential equation involving P′ alone, namely

$$\begin{array}{cc}\hfill \frac{{\partial }^2{P}^{\prime }}{\partial r^2}& +[\frac{1}{r}-\frac{\frac{\mathrm{d}}{\mathrm{d}r}\left(\rho (r)({\mathrm{\Omega }}^2(r)-k_z^2{v}_{\mathrm{A}}^2(r))\right)}{\rho (r)({\mathrm{\Omega }}^2(r)-k_z^2{v}_{\mathrm{A}}^2(r))}]\frac{\partial P^{\prime }}{\partial r}\hfill \\ \hfill & +(\frac{\rho (r)({\mathrm{\Omega }}^2(r)-k_z^2{v}_{\mathrm{A}}^2(r))}{B^2/\mu }-\frac{m^2}{r^2})P^{\prime }=0,\hfill \end{array}$$(9)

where $v_A^2(r)=B^2/\mu \rho (r)$ is the square of the spatially-dependent Alfvén velocity. We note that Eq. (9) can also be obtained from Eq. (18) of Goossens et al. (1992) in the β = 0 case and in the absence of magnetic twist.

The components of the Lagrangian displacement are related to P′ as

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{r}}& =\frac{1}{\rho (r)({\mathrm{\Omega }}^2(r)-k_z^2{v}_{\mathrm{A}}^2(r))}\frac{\partial P^{\prime }}{\partial r},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\xi }_{\phi }& =\frac{1}{\rho (r)({\mathrm{\Omega }}^2(r)-k_z^2{v}_{\mathrm{A}}^2(r))}\frac{\mathit{im}}{r}{P}^{\prime }.\hfill \end{array}$$(11)

along with ξ_z = 0 because of the β = 0 approximation. In the absence of flow, i.e., for U(r)=0 so that Ω(r)=ω, Eqs. (9)–(11) consistently revert to Eqs. (4)–(6) of Soler et al. (2013).

2.3. Solution in the internal and external plasmas

In the regions with constant density and flow velocity, Eq. (9) simplifies to

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{P}^{\prime }}{\mathrm{d}{r}^2}+\frac{1}{r}\frac{\mathrm{d}{P}^{\prime }}{\mathrm{d}r}+(\frac{{\mathrm{\Omega }}^2-k_z^2{v}_{\mathrm{A}}^2}{v_{\mathrm{A}}^2}-\frac{m^2}{r^2})P^{\prime }=0,\end{array}$$(12)

where now v_A and Ω are constant. Equation (12) is the Bessel Equation and applies both in the internal (r ≤ R − l/2) and external (r ≥ R + l/2) plasmas. We use the subscripts “i” and “e” to denote quantities related to the internal and external plasmas, respectively.

In the internal plasma, P′ must be regular at r = 0. Thus, the physical solution of Eq. (12) is

$$\begin{array}{c}\hfill P_{\mathrm{i}}^{\prime }=A_{\mathrm{i}}{J}_m\left(k_{\perp ,\mathrm{i}}r\right),\end{array}$$(13)

where A_i is a constant, J_m is the Bessel function of the first kind of order m, and

$$\begin{array}{c}\hfill k_{\perp ,\mathrm{i}}^2=\frac{{\mathrm{\Omega }}_{\mathrm{i}}^2-k_z^2{v}_{\mathrm{A},\mathrm{i}}^2}{v_{\mathrm{A},\mathrm{i}}^2}·\end{array}$$(14)

In the external plasma, we required that P′ vanishes when r → ∞. This is the condition for the wave to be trapped. The physical solution to Eq. (12) is then

$$\begin{array}{c}\hfill P_{\mathrm{e}}^{\prime }=A_{\mathrm{e}}{K}_m\left(k_{\perp ,\mathrm{e}}r\right),\end{array}$$(15)

where again A_e is a constant, K_m is the modified Bessel function of the first kind of order m, and

$$\begin{array}{c}\hfill k_{\perp ,\mathrm{e}}^2=-\frac{{\mathrm{\Omega }}_{\mathrm{e}}^2-k_z^2{v}_{\mathrm{A},\mathrm{e}}^2}{v_{\mathrm{A},\mathrm{e}}^2}·\end{array}$$(16)

As in Soler et al. (2013) we have focussed on trapped waves and discard leaky waves from the present investigation. In the absence of flow, leaky waves have been investigated in detail by, for example, Cally (1986, 2003) in the case of tubes with a piecewise constant density, and by, for example, Stenuit et al. (1999), Nakariakov et al. (2012), Guo et al. (2016) in the case of transversely nonuniform tubes. As shown in Terradas et al. (2010) and Soler et al. (2011), when the flow velocity surpasses a certain threshold the trapped waves are forced to become leaky waves because their frequency is located above the external cut-off frequency. In the case of kink waves, i.e., for m = ±1, the transition to the leaky regime occurs for an internal flow velocity that is super-Alfvénic. Since the vast majority of observed flows in coronal flux tubes are sub-Alfvénic (see Reale 2014), for simplicity we restricted our analysis to sub-Alfvénic flows and so avoid the possibility that the waves become leaky. If leaky waves were to be considered, the solution in the external plasma (Eq. (15)) should be expressed using Hankel functions and the condition of out-going waves should be enforced, that is, that there is no energy input from infinity (see, e.g., Stenuit et al. 1999; Guo et al. 2016, for more details on the treatment of leaky modes). By restricting ourselves to sub-Alfvénic flows, we are also discarding the triggering of the Kelvin–Helmholtz instability due to velocity shear at the boundary of the tube. The Kelvin–Helmholtz instability in flux tubes with longitudinal shear flows has been extensively investigated in the literature (e.g., Holzwarth et al. 2007; Ryutova et al. 2010; Zhelyazkov 2015) and is not the subject of the present study.

2.4. Solution in the nonuniform boundary layer

In the nonuniform transitional layer, and for m ≠ 0, Eq. (9) is singular at the specific position, r = r_A, where the resonant condition ${\mathrm{\Omega }}^2(r_{\text{A}})=k_z^2{v}_A^2(r_{\text{A}})$ is satisfied. Here, r_A denotes the Alfvén resonance position, which is a regular singular point of Eq. (9). We can expand the resonant condition to find an implicit equation whose solution is the resonant position, namely

$$\begin{array}{c}\hfill {\omega }^2-2k_{z}U(r_{\mathrm{A}})\omega +k_z^2(U^2(r_{\mathrm{A}})-v_{\mathrm{A}}^2(r_{\mathrm{A}}))=0.\end{array}$$(17)

For arbitrary density and flow profiles, Eq. (17) has to be solved numerically to find r_A. An analytic expression of r_A is only possible for very specific and simple profiles. We note that Eq. (17) depends upon ω, which shall be obtained from the dispersion relation that, in turn, requires r_A to be known. So, Eq. (17) must be solved iteratively along with the dispersion relation.

Concerning the behavior of the waves, the presence of the resonance causes wave damping due to coupling between the global modes (quasi-modes) and the localized Alfvén waves around the resonance position. The temporal evolution of these coupled modes (see, e.g., Terradas et al. 2006; Pascoe et al. 2010; Soler & Terradas 2015) reveals that the energy of the global mode is ideally transferred to the Alfvén continuum modes in the nonuniform layer. These Alfvén continuum modes develop small length scales due to phase mixing. Finally, the small scales are dissipated by some non-ideal process such as, resistivity or viscosity. The present approach, which is based on ideal quasi-modes, allows us to study the first phase of this process, that is, the damping of the global mode, while the generation of small scales and their dissipation would requite other approaches beyond the aim of the present study.

As in Soler et al. (2013; see also references therein), we used the method of Frobenius to express the solution to Eq. (9) as an infinite power series expansion around the resonance position r = r_A. The method assumes that there is only one resonance. The existence of multiple resonances is only possible when the density and flow profiles have very peculiar shapes that hardly represent the actual conditions in coronal flux tubes. For instance, various resonances may occur if the density within the transitional layer is not monotonic and increases and decreases following an oscillatory pattern, or when the density varies smoothly but the flow velocity varies very abruptly. This last case is analyzed in some detail by Terradas et al. (2010).

We rewrite Eq. (9) as

$$\begin{array}{c}\hfill {(r-r_{\mathrm{A}})}^{2}h(r)\frac{{\partial }^2{P}^{\prime }}{\partial r^2}+(r-r_{\mathrm{A}})p(r)\frac{\partial P^{\prime }}{\partial r}+q(r)P^{\prime }=0,\end{array}$$(18)

where the functions h(r), p(r), and q(r) are defined as

$$\begin{array}{cc}\hfill h(r)& =r^{2}f(r),\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill p(r)& =r(r-r_{\mathrm{A}})(f(r)-r\frac{\partial f(r)}{\partial r}),\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill q(r)& ={(r-r_{\mathrm{A}})}^2(\frac{\mu }{B^2}{r}^{2}f(r)-m^2)f(r).\hfill \end{array}$$(21)

The functions h(r), p(r), and q(r) take the same form as in Soler et al. (2013). The difference resides in the expression of the function f(r), which now contains the effect of flow and is given by

$$\begin{array}{c}\hfill f(r)=\rho (r)({\mathrm{\Omega }}^2(r)-k_z^2{v}_{\mathrm{A}}^2(r))={\mathrm{\Omega }}^2(r)\rho (r)-k_z^2\frac{B^2}{\mu }·\end{array}$$(22)

We assumed that both the density and the flow velocity are analytic functions at r = r_A. We performed a Taylor series of f(r) around the location of the resonance as

$$\begin{array}{c}\hfill f(r)={\displaystyle \sum _{k=0}^{\infty }}{f}_k{(r-r_{\mathrm{A}})}^k,\end{array}$$(23)

with f₀ = 0 and

$$\begin{array}{c}\hfill f_k={\displaystyle \sum _{n=0}^k{\rho }_n\sum _{l=0}^{k-n}}{\mathrm{\Omega }}_l{\mathrm{\Omega }}_{k-n-l}\text{for}k\ge 1,\end{array}$$(24)

where ρ₀ = ρ(r_A), Ω₀ = ω − k_zU(r_A), and

$$\begin{array}{cc}\hfill {\rho }_k& =\frac{1}{k!}{\frac{{\mathrm{d}}^k\rho }{\mathrm{d}{r}^k}|}_{r=r_{\mathrm{A}}},\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill {\mathrm{\Omega }}_k& =\frac{1}{k!}{\frac{{\mathrm{d}}^k\mathrm{\Omega }}{\mathrm{d}{r}^k}|}_{r=r_{\mathrm{A}}}=-\frac{k_z}{k!}{\frac{{\mathrm{d}}^{k}U}{\mathrm{d}{r}^k}|}_{r=r_{\mathrm{A}}},\hfill \end{array}$$(26)

when k ≥ 1.

The general solution to Eq. (18) is

$$\begin{array}{c}\hfill P_{\mathrm{tr}}^{\prime }(r)=A_0{P}_1^{\prime }(r)+S_0{P}_2^{\prime }(r),\end{array}$$(27)

where the subscript “tr” denotes the transitional layer, A₀ and S₀ are constants, and $P_1^{\prime }(r)$ and $P_2^{\prime }(r)$ are two linearly independent solutions. We find the two independent solutions with the help of a Frobenius series expansion around the regular singular point r = r_A. The indicial equation is obtained from the coefficient of the lowest power of (r−r_A), and it shows that zero and two are the two possible indices of the expansion. Then, the two linearly independent solutions are

$$\begin{array}{cc}\hfill P_1^{\prime }(r)& ={(r-r_{\mathrm{A}})}^2{\displaystyle \sum _{k=0}^{\infty }}{a}_k{(r-r_{\mathrm{A}})}^k,\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill P_2^{\prime }(r)& ={\displaystyle \sum _{k=0}^{\infty }}{s}_k{(r-r_{\mathrm{A}})}^k+\mathcal{C}{P}_1^{\prime }(r)ln(r-r_{\mathrm{A}}),\hfill \end{array}$$(29)

where 𝒞 is the coupling constant and a_k and s_k are the series coefficients. We consider that the causal branch of the logarithm is the one where ln(r−r_A) = ln(r_A−r) ± iπ if r <  r_A, where the criterion for choosing either the + sign or the − sign is not arbitrary but based on the physical argument that the effect of the resonance is to produce the damping of the waves. The coupling constant, 𝒞, is independent of the density and flow profiles and is given by

$$\begin{array}{c}\hfill \mathcal{C}=\frac{m^2}{2r_{\mathrm{A}}^2}·\end{array}$$(30)

For sausage waves, meaning when m = 0, 𝒞 = 0, the singular logarithmic term is dropped from $P_2^{\prime }(r)$. As a consequence, no resonant damping occurs if m = 0. Conversely, the coefficients a_k and s_k depend upon the choice of the density and flow profiles. General expressions of the coefficients a_k and s_k are given in the Appendix A.

In the simplified case that the nonuniform layer is thin compared with the tube radius, it suffices to keep terms up to 𝒪(l/R) in the Frobenius series. This is the so-called thin boundary (TB) approximation. In that simple scenario, we find

$$\begin{array}{ccc}\hfill P_{\mathrm{tr}}^{\prime }(r)& \sim \hfill & \hfill S_0,\end{array}$$(31)

$$\begin{array}{ccc}\hfill {\xi }_{\mathrm{r},\mathrm{tr}}(r)& \sim \hfill & \hfill \frac{S_0}{f_1}\frac{m^2}{r_{\mathrm{A}}^2}ln(r-r_{\mathrm{A}}),\end{array}$$(32)

$$\begin{array}{ccc}\hfill {\xi }_{\phi ,\mathrm{t}r}(r)& \sim \hfill & \hfill i\frac{S_0}{f_1}\frac{m}{r_{\mathrm{A}}}\frac{1}{r-r_{\mathrm{A}}}·\end{array}$$(33)

In a thin nonuniform layer the total pressure perturbation is constant, the radial displacement jumps logarithmically, and the azimuthal displacement behaves as 1/(r − r_A). This behavior agrees with that found in previous works (e.g., Goossens et al. 1992). In the present work we have used the full Frobenious solution that contains this fundamental behavior as well as the corrections beyond the TB approximation owing to the larger thickness of the nonuniform layer.

2.5. Dispersion relation

The conditions that the total pressure perturbation, P′, and the radial component of the Lagrangian displacement, ξ_r, are continuous at r = R − l/2 and r = R + l/2 provide us with a system of four algebraic equations for the constants A_i, A_e, A₀, and S₀. The dispersion relation of the waves is then obtained from the requirement that there is a nontrivial solution of the system, in other words, by setting the associated determinant equal to zero. The expression of the general dispersion relation is

$$\begin{array}{cc}& \frac{\frac{k_{\perp ,\mathrm{e}}}{{\rho }_{\mathrm{e}}({\mathrm{\Omega }}_{\mathrm{e}}^2-k_z^2{v}_{\mathrm{A},\mathrm{e}}^2)}\frac{K_m^{\prime }\left[k_{\perp ,\mathrm{e}}(R+l/2)\right]}{K_m\left[k_{\perp ,\mathrm{e}}(R+l/2)\right]}{\mathcal{G}}_{+}-{\mathrm{\Xi }}_{+}}{\frac{k_{\perp ,\mathrm{e}}}{{\rho }_{\mathrm{e}}({\mathrm{\Omega }}_{\mathrm{e}}^2-k_z^2{v}_{\mathrm{A},\mathrm{e}}^2)}\frac{K_m^{\prime }\left[k_{\perp ,\mathrm{e}}(R+l/2)\right]}{K_m\left[k_{\perp ,\mathrm{e}}(R+l/2)\right]}{\mathcal{F}}_{+}-{\mathrm{\Gamma }}_{+}}\hfill \\ \hfill & -\frac{\frac{k_{\perp ,\mathrm{i}}}{{\rho }_{\mathrm{i}}({\mathrm{\Omega }}_{\mathrm{i}}^2-k_z^2{v}_{\mathrm{A},\mathrm{i}}^2)}\frac{J_m^{\prime }\left[k_{\perp ,\mathrm{i}}(R-l/2)\right]}{J_m\left[k_{\perp ,\mathrm{i}}(R-l/2)\right]}{\mathcal{G}}_{-}-{\mathrm{\Xi }}_{-}}{\frac{k_{\perp ,\mathrm{i}}}{{\rho }_{\mathrm{i}}({\mathrm{\Omega }}_{\mathrm{i}}^2-k_z^2{v}_{\mathrm{A},\mathrm{i}}^2)}\frac{J_m^{\prime }\left[k_{\perp ,\mathrm{i}}(R-l/2)\right]}{J_m\left[k_{\perp ,\mathrm{i}}(R-l/2)\right]}{\mathcal{F}}_{-}-{\mathrm{\Gamma }}_{-}}=0.\hfill \end{array}$$(34)

The quantities 𝒢_±, ℱ_±, Ξ_±, and Γ_± are defined as

$$\begin{array}{cc}\hfill {\mathcal{G}}_{\pm }& ={\displaystyle \sum _{k=0}^{\infty }}{a}_k{\zeta }_{\pm }^{k+2},\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill {\mathcal{F}}_{\pm }& ={\displaystyle \sum _{k=0}^{\infty }}(s_k{\zeta }_{\pm }^k+\frac{m^2}{2r_{\mathrm{A}}^2}ln\left({\zeta }_{\pm }\right)a_k{\zeta }_{\pm }^{k+2}),\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill {\mathrm{\Xi }}_{\pm }& =\frac{1}{{\sum }_{k=0}^{\infty }{f}_{k+1}{\zeta }_{\pm }^k}{\displaystyle \sum _{k=0}^{\infty }}(k+2)a_k{\zeta }_{\pm }^k,\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\pm }& =\frac{1}{{\sum }_{k=0}^{\infty }{f}_{k+1}{\zeta }_{\pm }^k}{\displaystyle \sum _{k=0}^{\infty }}(k{s}_k{\zeta }_{\pm }^{k-2}+\frac{m^2}{2r_{\mathrm{A}}^2}{a}_k{\zeta }_{\pm }^k\hfill \\ \hfill & +\frac{m^2}{2r_{\mathrm{A}}^2}ln\left({\zeta }_{\pm }\right)(k+2)a_k{\zeta }_{\pm }^k),\hfill \end{array}$$(38)

with ${\zeta }_{\pm }=R\pm \frac{l}{2}-r_{\text{A}}$. In the absence of flow the dispersion relation reverts to that given in Soler et al. (2013).

Equation (34) is a transcendental equation with roots that are to be found numerically. To do this, we used a numerical routine based on that used by Soler et al. (2013). A difference of this case from the one without flow in Soler et al. (2013) is that the resonance position, r_A, also needs to be found numerically from Eq. (17). The transverse density profiles used by Soler et al. (2013) allow analytic expressions for r_A, which is not possible here because of the presence of flow. Hence, the present method iteratively solves Eqs. (17) and (34) until both solutions converge. We note that the dispersion relation involves series with infinite number of terms. To proceed numerically we must truncate the infinite series so that only the first N terms are accounted for. To make sure that the number of terms considered is large enough for the error to be negligible, we performed convergence tests by increasing N until a good convergence of the solution is obtained. Typically, we consider N = 51.

We shall find the solutions to Eq. (34) in the case of fixed, real, and positive k_z. So, we impose a particular wavelength of the perturbations, namely λ = 2π/k_z. Then, the solution is a complex frequency, namely ω = ω_R + iω_I, where the subscripts R and I denote the real and imaginary parts, respectively. The sign of ω_R indicates the direction of wave propagation. If ω_R >  0 the wave propagates toward the positive z-direction (forward propagation), while propagation is in the opposite direction when ω_R <  0 (backward propagation). The phase velocity is computed as

$$\begin{array}{c}\hfill v_{\mathrm{ph}}=\frac{{\omega }_R}{k_z},\end{array}$$(39)

where the sign of v_ph follows the same rules as the sign of ω_R. On the other hand, ω_I is related to the damping rate of the waves, so that ω_I <  0 for both directions of propagation. If the waves are not overdamped, i.e., if |ω_I| <  |ω_R|, we can define the exponential damping length of the waves, L_D, as the distance the waves need to travel for their amplitude to be reduced by a factor of e, namely (see Tagger et al. 1995)

$$\begin{array}{c}\hfill L_{\mathrm{D}}=\left|\frac{1}{{\omega }_I}\frac{\partial {\omega }_R}{\partial k_z}\right|.\end{array}$$(40)

The damping length so defined is positive for both forward and backward propagating waves.

2.6. Recovering the thin tube, thin boundary, and slow flow approximation

Equation (34) is valid for arbitrary values of l/R and k_zR. An appropriate way to check the validity of Eq. (34) is to recover the dispersion relation previously obtained in the literature in the thin tube (TT, k_zR ≪ 1) and thin boundary (TB, l/R ≪ 1) approximations (e.g., Goossens et al. 1992; Terradas et al. 2010; Soler et al. 2011). To do so, we kept terms up to 𝒪(l/R) in the expressions of 𝒢_±, ℱ_±, Ξ_±, and Γ_±, and approximate r_A ≈ R as consistent with the assumption that the boundary layer is so thin that the resonance position is necessarily close to r = R. Then, we find 𝒢_± ≈ 0, ℱ_± ≈ 1, Ξ_± ≈ 2/f₁, and

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{+}& \approx \frac{m^2/R^2}{f_1}{\displaystyle (ln\frac{l}{2}+\frac{1}{2})},\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{-}& \approx \frac{m^2/R^2}{f_1}{\displaystyle (ln\frac{l}{2}\pm i\pi +\frac{1}{2})}.\hfill \end{array}$$(42)

In order to make further analytic progress, we assumed that the magnetic tube is thin and use the first-order expansion for small arguments and m ≠ 0 of the Bessel functions in Eq. (34) (see, e.g., Abramowitz & Stegun 1972) and also take R − l/2 ≈ R + l/2 ≈ R, namely

$$\begin{array}{cc}\hfill \frac{J_m^{\prime }\left(k_{\perp ,\mathrm{i}}(R-l/2)\right)}{J_m\left(k_{\perp ,\mathrm{i}}(R-l/2)\right)}& \approx \frac{J_m^{\prime }\left(k_{\perp ,\mathrm{i}}R\right)}{J_m\left(k_{\perp ,\mathrm{i}}R\right)}\approx \frac{m}{k_{\perp ,\mathrm{i}}R},\hfill \end{array}$$(43)

$$\begin{array}{cc}\hfill \frac{K_m^{\prime }\left(k_{\perp ,\mathrm{e}}(R+l/2)\right)}{K_m\left(k_{\perp ,\mathrm{e}}(R+l/2)\right)}& \approx \frac{K_m^{\prime }\left(k_{\perp ,\mathrm{e}}R\right)}{K_m\left(k_{\perp ,\mathrm{e}}R\right)}\approx -\frac{m}{k_{\perp ,\mathrm{e}}R}·\hfill \end{array}$$(44)

After long but straightforward algebraic manipulations, the TT and TB approximation of Eq. (34) can be cast as

$$\begin{array}{cc}& {\rho }_{\mathrm{i}}({\mathrm{\Omega }}_{\mathrm{i}}^2-k_z^2{v}_{\mathrm{A},\mathrm{i}}^2)+{\rho }_{\mathrm{e}}({\mathrm{\Omega }}_{\mathrm{e}}^2-k_z^2{v}_{\mathrm{A},\mathrm{e}}^2)\hfill \\ \hfill & =i\pi \frac{m/R}{\left|f_1\right|}{\rho }_{\mathrm{i}}({\mathrm{\Omega }}_{\mathrm{i}}^2-k_z^2{v}_{\mathrm{A},\mathrm{i}}^2){\rho }_{\mathrm{e}}({\mathrm{\Omega }}_{\mathrm{e}}^2-k_z^2{v}_{\mathrm{A},\mathrm{e}}^2),\hfill \end{array}$$(45)

with

$$\begin{array}{c}\hfill \left|f_1\right|={|{\mathrm{\Omega }}^2\frac{\partial \rho }{\partial r}-2k_z\mathrm{\Omega }\rho \frac{\partial U}{\partial r}|}_{r=R}.\end{array}$$(46)

Equation (45) agrees with Eq. (74) of Goossens et al. (1992), which is also used in Terradas et al. (2010) and Soler et al. (2011). Thus, the TT and TB dispersion relation is correctly recovered from the more general Eq. (34). We refer readers to those works to study the solutions of the TT and TB dispersion relation.

Soler et al. (2011) considered Eq. (45) and obtained approximate expressions of v_ph and L_D with U_e = 0 and in the case of slow flow (SF, U_i/v_A,i ≪ 1), namely

$$\begin{array}{cc}\hfill v_{\mathrm{ph}}& \approx \pm v_k+\frac{{\rho }_{\mathrm{i}}}{{\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}}}{U}_{\mathrm{i}},\hfill \end{array}$$(47)

$$\begin{array}{cc}\hfill L_{\mathrm{D}}& \approx \frac{2\pi }{k_z}\frac{F}{m}\frac{R}{l}\frac{{\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}-{\rho }_{\mathrm{e}}}\frac{v_k}{v_{\mathrm{ph}}}(1\pm \frac{2{\rho }_{\mathrm{i}}}{{\rho }_{\mathrm{i}}-{\rho }_{\mathrm{e}}}\frac{U_{\mathrm{i}}}{v_k})\hfill \\ \hfill & \approx \frac{2\pi }{k_z}\frac{F}{m}\frac{R}{l}\frac{{\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}-{\rho }_{\mathrm{e}}}(1\pm \sqrt{\frac{{\rho }_{\mathrm{i}}}{2({\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}})}}\frac{{\rho }_{\mathrm{i}}+3{\rho }_{\mathrm{e}}}{{\rho }_{\mathrm{i}}-{\rho }_{\mathrm{e}}}\frac{U_{\mathrm{i}}}{v_{\mathrm{A},\mathrm{i}}}),\hfill \end{array}$$(48)

where the + and − signs stand for forward and backward propagating waves respectively, F is a numerical factor that depends on the shape of the transitional layer (see Soler et al. 2014), and v_k is the kink velocity given by

$$\begin{array}{c}\hfill v_k=\sqrt{\frac{{\rho }_{\mathrm{i}}{v}_{\mathrm{A},\mathrm{i}}^2+{\rho }_{\mathrm{e}}{v}_{\mathrm{A},\mathrm{e}}^2}{{\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}}}}=\sqrt{\frac{2{\rho }_{\mathrm{i}}}{{\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}}}}{v}_{\mathrm{A},\mathrm{i}}.\end{array}$$(49)

Although Eqs. (47) and (48) are only strictly valid in the limits k_zR ≪ 1, l/R ≪ 1, and U_i/v_A,i ≪ 1 they provide useful information regarding the effects of the various parameters. We note that Eqs (47) and (48) predict that the effect of flow is to produce a linear correction in the flow velocity to both v_ph and L_D. We compare these approximations with the solutions of the general Eq. (34), which remains valid beyond the range of applicability of Eqs. (47) and (48).

3. Application

As a brief application of the method described above, we computed results corresponding to the fundamental radial harmonic of the m = 1 waves, in other words, the so-called kink mode. We note, however, that the dispersion relation is valid for any value of m and for any radial harmonic as long as the modes are not leaky. In all the numerical solutions, we considered a sinusoidal transition for both the density and the flow velocity within the nonuniform boundary layer. For simplicity, we chose the reference frame so that the external plasma is static, meaning that we set U_e = 0. In addition, we set the ratio of the internal density to the external density to ρ_i/ρ_e = 3, which is representative of a coronal loop. For this density contrast, v_k ≈ 1.22v_A,i.

We investigated the behavior of v_ph and L_D as functions of three dimensionless quantities: the ratio of the internal flow velocity to the internal Alfvén velocity, U_i/v_A,i, the relative thickness of the boundary layer, l/R, and dimensionless longitudinal wavenumber, k_zR. Typical values of v_A,i and R in coronal loops are v_A,i ∼ 1000 km s⁻¹ and R ∼ 3000 km.

3.1. Case without flow

To start with, we considered the case with U_i = 0. In this scenario, forward and backward propagating waves are degenerate, so for simplicity, only the forward wave was considered. This corresponds to the situation studied by Soler et al. (2013). First of all, we check that the solutions given in Soler et al. (2013) are correctly recovered with the present numerical routine. We have fully confirmed this for some selected configurations (this analysis is not shown here).

We then focussed on investigating the dependence on k_zR, which was not done in Soler et al. (2013). Here we have studied the impact of this parameter. Figure 1 shows contour plots of v_ph/v_A,i and L_D/R as functions of l/R and k_zR. Concerning the phase velocity, we find that regardless the value of k_zR, there is an increasing trend in the phase velocity as l/R increases, in agreement with the results of Soler et al. (2013) for fixed k_zR. Conversely, the effect of increasing k_zR is the opposite, in other words, the phase velocity decreases. Therefore, the parameters l/R and k_zR have competing effects on the phase velocity. When l = 0, the phase velocity approaches the internal Alfvén velocity when k_zR increases (see Edwin & Roberts 1983), but when l ≠ 0 the phase velocity tends to a value somewhat larger than v_A,i when k_zR increases. This effect is due to the transverse nonuniformity.

Fig. 1. Results in the absence of flow. Contour plots of vph/vA,i (panel a) and LD/R (panel b) as functions of l/R and kzR. We note that LD/R is given in logarithmic scale. The dashed purple line in panel a denotes vph = vk.

On the other hand, l/R and k_zR have the same effect on the damping length, since L_D/R decreases when any of the two parameters increases. This result agrees qualitatively with the approximate analytic dependence of Eq. (48). From the physical point of view, the damping length becomes smaller when l/R increases because the efficiency of resonant damping grows, while the effect of increasing k_zR is to produce shorter wavelengths (i.e., higher frequencies) that are more efficiently damped.

In order to test the accuracy of the approximate Eqs. (47) and (48), we define the normalized errors of v_ph and L_D associated with the use of the approximations as

$$\begin{array}{cc}\hfill \mathrm{err}\left(v_{\mathrm{ph}}\right)& =\frac{v_{\mathrm{ph},\mathrm{F}\mathrm{rob}}-v_{\mathrm{ph},\mathrm{Frob}}}{v_{\mathrm{ph},\mathrm{Frob}}},\hfill \end{array}$$(50)

$$\begin{array}{cc}\hfill \mathrm{err}\left(L_{\mathrm{D}}\right)& =\frac{L_{\mathrm{D},\mathrm{Frob}}-L_{\mathrm{D},\mathrm{app}}}{L_{\mathrm{D},\mathrm{Frob}}},\hfill \end{array}$$(51)

where v_ph, app and L_D, app denote the analytic approximations given by Eqs. (47) and (48), respectively, and v_ph, Frob and L_D, Frob are the actual results obtained with the present Frobenius-based method. Figure 2 shows the normalized errors as functions of l/R and k_zR. Very small errors are obtained when l/R ≪ 1 and k_zR ≪ 1, as consistent with the regime of applicability of the approximations. For larger values of the two parameters, the errors display a nonmonotonic behavior. In the case of the phase velocity, the largest errors are ∼15% and are obtained either when l/R → 2 and k_zR ≪ 1 or when l/R ≪ 1 and k_zR ≫ 1. On the contrary, the error in the case of the damping length seems to be more dominated by k_zR. The largest values of err(L_D) are ∼70% and are obtained when k_zR ≫ 1, while the error when k_zR ≪ 1 remains moderate even for large values of l/R. Remarkably, owing to the peculiar nonmonotonic behavior of the errors, we obtain zero errors of both v_ph and L_D for some particular combination of parameters far beyond the regime of applicability of the approximations.

Fig. 2. Results in the absence of flow. Contour plots of the errors associated with the approximate Eqs. (47) and (48) for panel avph and panel bLD as functions of l/R and kzR. The dashed purple line in both panels denotes the contour of zero error.

3.2. Effect of flow

We took flow into account and exploited the new extension to the method. Figure 3 shows contour plots of v_ph/v_A,i and L_D/R as functions of l/R and U_i/v_A,i for k_zR = 0.1. Because of the flow, different results are obtained for forward and backward waves, so Fig. 3 displays both. Physically, the effect of the flow is to drag the waves toward to direction of the flow, and hence a shift of the phase velocity is produced. This results in larger phase velocities for the forward propagating wave and smaller phase velocities (in absolute value) for the backward propagating wave when compared with the case without flow. As in the absence flow, the effect of increasing l/R is to increase the phase velocity (in absolute value) for both forward and backward waves. Thus, l/R and U_i/v_A,i have the same effect on the phase velocity of forward waves, while these parameters have competing effects on the phase velocity of backward waves.

Fig. 3. Results in the presence of flow. Contour plots of vph/vA,i (panels a and b) and LD/R (panels c and d) as functions of l/R and Ui/vA,i. Panels a and c: forward propagating wave, while panels b and d are for the backward propagating wave. We note that LD/R is given in logarithmic scale. We have used kzR = 0.1.

Regarding the damping length, we find that for fixed k_zR, the behavior of L_D is essentially governed by l/R, while the role of the flow is to produce a small positive (for the forward wave) or negative (for the backward wave) shift with respect to the value for U_i = 0. Hence, the backward wave is more efficiently damped than the forward wave, that is, the backward wave has a shorter L_D. Qualitatively, this behavior agrees with the approximate analytic dependence of Eq. (48) and is also consistent with the resistive MHD results of Soler et al. (2011).

As before, we have computed the errors associated with the use of the approximations. These results are given in Fig. 4. In the case of the phase velocity, we obtain that the error grows with l/R so that the largest errors are found when l/R = 2. Clearly, the departure from the theoretical range of the TB approximation significantly reduces the accuracy of Eq. (47). In this regard, the error of backward wave phase velocity is larger than that of the forward wave. In the considered range of parameters, the maximum error of phase velocity is ∼13% for the forward wave and ∼23% for the backward wave. Conversely, the departure from the theoretical range of the SF approximation produces comparatively small errors. Equation (47) predicts a linear shift of the phase velocity with the flow velocity, and such a dependence remains quite accurate in the whole range of considered flow velocities.

Fig. 4. Results in the presence of flow. Contour plots of the errors associated with the approximate Eqs. (47) and (48) for vph (panels a and b) and LD (panels c and d) as functions of l/R and Ui/vA,i. Panels a and c: forward propagating wave, while panels b and d are for the backward propagating wave. The dashed purple line denotes the contour of zero error. We have used kzR = 0.1.

Next we considered to the error of the damping length. As in the case without flow, err(L_D) displays a nonmonotonic behavior with l/R. This fact makes it is possible to find very small errors (even zero error) associated with the use of Eq. (48) when l/R is large and some specific values of the flow velocity are considered. On the other hand, the dependence of err(L_D) with U_i/v_A,i shows a much simpler relation, namely err(L_D) typically grows when U_i/v_A,i increases. The largest values of err(L_D) plotted in Fig. 4 are obtained when the flow velocity takes the maximum value used in the computations and, as before, the backward wave is the solution for which the approximation performs worst. As a general rule, the linear dependence on U_i/v_A,i predicted by the analytic Eqs. (47) and (48) as a consequence of adopting SF approximation is less accurate for L_D than for v_ph.

4. Concluding remarks

The semi-analytic technique developed by Soler et al. (2013) to compute transverse waves in flux tubes with thick boundary layers has been extended to incorporate the effect of longitudinal mass flows. The flow velocity is allowed to vary within the nonuniform boundary from the internal velocity to the external velocity. In the past, similar configurations have been studied analytically and/or numerically with resistive eigenvalue computations but in the case of tubes with thin transitions (see, e.g., Goossens et al. 1992; Terradas et al. 2010; Soler et al. 2011). Two advantages of the present semi-analytic approach based on the Frobenius method are that nonuniform boundaries of arbitrary width can be considered and that its numerical implementation is much faster than resistive eigenvalue computations. Therefore, detailed studies of the impact of the various model parameters on the wave properties are feasible with the present approach. We foresee future works that could exploit the method.

As an exemplary application, we have performed a parameter study of the effect of flow on the phase velocity and damping length of resonantly damped kink waves. First, we consistently recover the results in the TT, TB, and SF approximations obtained in previous works. Thanks to the present approach, we then extend those results beyond the range of applicability of the approximations. We have focussed on testing the validity of the TT, TB, and SF approximations against the actual results provided by the Frobenius-based method. While the TT and SF approximations perform relatively well for the considered wavenumbers and flow velocities, the use of the TB approximation implies a significant error when used beyond the limit l/R ≪ 1. We note that we computed the results for a fixed value of the density ratio ρ_i/ρ_e and boundary layer of sinusoidal shape, while Soler et al. (2014) showed that both ingredients also affect the accuracy of the approximations. Therefore, it is difficult to deduce a simple universal estimation of the error associated with the approximations, especially the TB approximation, since their accuracy is quite sensitive to the background configuration considered.

We have used the simplification that the density and the flow velocity follow that same radial dependence. Terradas et al. (2010) and Soler et al. (2011) considered the case that the density and flow velocity vary transversely to the tube within transitional layers of different width. They found that this has little impact on the results except in the case that the transitional layer for the flow velocity is much thinner than the corresponding layer for the density. In that case, the forward propagating wave damps faster than the backward propagating wave (see details in Terradas et al. 2010; Soler et al. 2011). This is, however, a very peculiar situation that hardly represents the expected conditions in coronal flux tubes.

Besides the effect of flows, there are still various physical ingredients that could be incorporated to the semi-analytic Frobenius-based method. For instance, two additional improvements of the method would be to consider longitudinal stratification (e.g., Andries et al. 2005; Arregui et al. 2005) and magnetic twist (e.g., Terradas & Goossens 2012; Ruderman & Terradas 2015). These extensions could be tackled in forthcoming works.

Acknowledgments

RS acknowledges the support from grant AYA2017-85465-P (MINECO/AEI/FEDER, UE) and from the Ministerio de Economía, Industria y Competitividad and the Conselleria d’Innovació, Recerca i Turisme del Govern Balear (Pla de ciència, tecnologia, innovació i emprenedoria 2013-2017) for the Ramón y Cajal grant RYC-2014-14970.

References

Appendix A: Expressions of the Frobenius series coefficients

The expressions of the Frobenius series coefficients a_k and s_k are as follows:

$$\begin{array}{cc}\hfill a_0& =1,\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill a_1& =-\frac{2f_1-2f_2{r}_{\mathrm{A}}}{3f_1{r}_{\mathrm{A}}}{a}_0,\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill a_2& =-\frac{9f_1{r}_{\mathrm{A}}{a}_1+(2f_1-2f_2{r}_{\mathrm{A}}-4f_3{r}_{\mathrm{A}}^2-m^2{f}_1)a_0}{8f_1{r}_{\mathrm{A}}^2},\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill a_3& =-\frac{1}{15f_1{r}_{\mathrm{A}}^2}[(4f_2{r}_{\mathrm{A}}^2+20f_1{r}_{\mathrm{A}})a_2+(-3f_3{r}_{\mathrm{A}}^2+3f_2{r}_{\mathrm{A}}+6f_1-m^2{f}_1)a_1\hfill \\ \hfill & +(-6f_4{r}_{\mathrm{A}}^2-6f_3{r}_{\mathrm{A}}+\frac{r_{\mathrm{A}}^2}{B^2/\mu }{f}_1^2-m^2{f}_2)a_0],\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill a_4& =-\frac{1}{24f_1{r}_{\mathrm{A}}^2}{\displaystyle [\sum _{j=0}^3(j+2)(2j-4)r_{\mathrm{A}}^2{f}_{5-j}{a}_j+\sum _{j=0}^3(j+2)(4j-5)r_{\mathrm{A}}{f}_{4-j}{a}_j}\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^2((j+2)(2j-1)-m^2)f_{3-j}{a}_j+\sum _{j=0}^1\sum _{l=0}^{1-j}\frac{r_{\mathrm{A}}^2}{B^2/\mu }{f}_{l+1}{f}_{2-j-l}{a}_j+2\frac{r_{\mathrm{A}}}{B^2/\mu }{f}_1^2{a}_0]},\hfill \\ \hfill a_k& =-\frac{1}{k(k+2)f_1{r}_{\mathrm{A}}^2}{\displaystyle [\sum _{j=0}^{k-1}(j+2)(2j-k)r_{\mathrm{A}}^2{f}_{k-j+1}{a}_j+\sum _{j=0}^{k-1}(j+2)(4j-2k+3)r_{\mathrm{A}}{f}_{k-j}{a}_j}\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{k-2}((j+2)(2j-k+3)-m^2)f_{k-j-1}{a}_j+\sum _{j=0}^{k-3}\sum _{l=0}^{k-j-3}\frac{r_{\mathrm{A}}^2}{B^2/\mu }{f}_{l+1}{f}_{k-j-l-2}{a}_j}\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{k-4}\sum _{l=0}^{k-j-4}\frac{2r_{\mathrm{A}}}{B^2/\mu }{f}_{l+1}{f}_{k-j-l-3}{a}_j+\sum _{j=0}^{k-5}\sum _{l=0}^{k-j-5}\frac{\mu }{B^2}{f}_{l+1}{f}_{k-j-l-4}{a}_j]},\text{for}k\ge 5,\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill s_0& =1,\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill s_1& =0,\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill s_2& =0,\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill s_3& =\frac{1}{3f_1{r}_{\mathrm{A}}^2}[(m^2{f}_2-\frac{r_{\mathrm{A}}^2}{B^2/\mu }{f}_1^2)s_0-\mathcal{C}(4r_{\mathrm{A}}^2{f}_1{a}_1+(r_{\mathrm{A}}^2{f}_2+5r_{\mathrm{A}}{f}_1)a_0)],\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill s_4& =-\frac{1}{8f_1{r}_{\mathrm{A}}^2}[9f_1{r}_{\mathrm{A}}{s}_3+(\frac{2r_{\mathrm{A}}}{B^2/\mu }{f}_1^2+\frac{2r_{\mathrm{A}}^2}{B^2/\mu }{f}_1{f}_2-m^2{f}_3)s_0\hfill \\ \hfill & +\mathcal{C}(6r_{\mathrm{A}}^2{f}_1{a}_2+3r_{\mathrm{A}}^2{f}_2{a}_1+9r_{\mathrm{A}}{f}_1{a}_1+3(r_{\mathrm{A}}{f}_2+f_1)a_0)],\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill s_k& =-\frac{1}{k(k-2)f_1{r}_{\mathrm{A}}^2}{\displaystyle \{\sum _{j=0}^{k-1}j(2j-k-2)r_{\mathrm{A}}^2{f}_{k-j+1}{s}_j+\sum _{j=0}^{k-1}j(4j-2k-1)r_{\mathrm{A}}{f}_{k-j}{s}_j}\hfill \end{array}$$(A.11)

$$\begin{array}{cc}& +{\displaystyle \sum _{j=0}^{k-2}}[(j(2j-k+1)-m^2)f_{k-j-1}{s}_j+\mathcal{C}(3j-k+4)r_{\mathrm{A}}^2{f}_{k-j-1}{a}_j]\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{k-3}[\sum _{l=0}^{k-j-3}\frac{r_{\mathrm{A}}^2}{B^2/\mu }{f}_{l+1}{f}_{k-j-l-2}{s}_j+\mathcal{C}(6j-2k+11)r_{\mathrm{A}}{f}_{k-j-2}{a}_j]}\hfill \\ \hfill & {\displaystyle +\sum _{j=0}^{k-4}[\sum _{l=0}^{k-j-4}\frac{2r_{\mathrm{A}}}{B^2/\mu }{f}_{l+1}{f}_{k-j-l-3}{s}_j+\mathcal{C}(3j-k+7)f_{k-j-3}{a}_j]}\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{k-5}\sum _{l=0}^{k-j-5}\frac{\mu }{B^2}{f}_{l+1}{f}_{k-j-l-4}{s}_j\}},\text{for}k\ge 5.\hfill \end{array}$$(A.12)

In the absence of flow, f_k = ω²ρ_k and the coefficients consistently revert to those given in Soler et al. (2013).

All Figures

	Fig. 1. Results in the absence of flow. Contour plots of vph/vA,i (panel a) and LD/R (panel b) as functions of l/R and kzR. We note that LD/R is given in logarithmic scale. The dashed purple line in panel a denotes vph = vk.
In the text

	Fig. 2. Results in the absence of flow. Contour plots of the errors associated with the approximate Eqs. (47) and (48) for panel avph and panel bLD as functions of l/R and kzR. The dashed purple line in both panels denotes the contour of zero error.
In the text

	Fig. 3. Results in the presence of flow. Contour plots of vph/vA,i (panels a and b) and LD/R (panels c and d) as functions of l/R and Ui/vA,i. Panels a and c: forward propagating wave, while panels b and d are for the backward propagating wave. We note that LD/R is given in logarithmic scale. We have used kzR = 0.1.
In the text

	Fig. 4. Results in the presence of flow. Contour plots of the errors associated with the approximate Eqs. (47) and (48) for vph (panels a and b) and LD (panels c and d) as functions of l/R and Ui/vA,i. Panels a and c: forward propagating wave, while panels b and d are for the backward propagating wave. The dashed purple line denotes the contour of zero error. We have used kzR = 0.1.
In the text