Transverse waves in coronal flux tubes with thick boundaries: The effect of longitudinal flows

TL;DR

This paper develops a semi-analytic method to study how longitudinal flows influence resonant damping of transverse MHD waves in coronal flux tubes with thick nonuniform boundaries, extending previous thin boundary models.

Contribution

It introduces a semi-analytic technique that incorporates flow effects into the analysis of resonant MHD waves in thick-boundary coronal flux tubes, expanding beyond thin boundary assumptions.

Findings

Flow affects phase velocity and damping length of kink waves.

Results recover thin boundary approximation in the appropriate limit.

Extended analysis to thick boundary regimes, assessing approximation errors.

Abstract

Observations show that transverse magnetohydrodynamic (MHD) waves and flows are often simultaneously present in magnetic loops of the solar corona. The waves are resonantly damped in the Alfv\'en continuum because of plasma and/or magnetic field nonuniformity across the loop. The resonant damping is relevant in the context of coronal heating, since it provides a mechanism to cascade energy down to the dissipative scales. It has been theoretically shown that the presence of flow affects the waves propagation and damping, but most of the studies rely on the unjustified assumption that the transverse nonuniformity is confined to a boundary layer much thinner than the radius of the loop. Here we present a semi-analytic technique to explore the effect of flow on resonant MHD waves in coronal flux tubes with thick nonuniform boundaries. We extend a published method, which was originally developed for a static plasma, in order to incorporate the effect of flow. We allowed the flow velocity to continuously vary within the nonuniform boundary from the internal velocity to the external velocity. The analytic part of the method is based on expressing the wave perturbations in the thick nonuniform boundary of the loop as a Frobenius series that contains a singular term accounting for the Alfv\'en resonance, while the numerical part of the method consists of solving iteratively the transcendental dispersion relation together with the equation for the Alfv\'en resonance position. As an application of this method, we investigated the impact of flow on the phase velocity and resonant damping length of MHD kink waves. We consistently recover results in the thin boundary approximation obtained in previous studies. We have extended those results to the case of thick boundaries. We also explored the error associated with the use of the thin boundary approximation beyond its regime of applicability.

Equations

$$\rho(r)=\left\{\begin{array}[]{lll}\rho_{\rm i},&\textrm{if}&r\leq R-l/2,\\ \rho_{\rm tr}(r),&\textrm{if}&R-l/2<r<R+l/2,\\ \rho_{\rm e},&\textrm{if}&r\geq R+l/2,\end{array}\right.$$

$$U(r)=\left\{\begin{array}[]{lll}U_{\rm i},&\textrm{if}&r\leq R-l/2,\\ U_{\rm tr}(r),&\textrm{if}&R-l/2<r<R+l/2,\\ U_{\rm e},&\textrm{if}&r\geq R+l/2,\end{array}\right.$$

$$\rho \left(\frac{\partial \mathbf{v}}{\partial t}+\mathbf{U}\cdot \nabla \mathbf{v}+\mathbf{v}\cdot \nabla \mathbf{U}\right)$$

$$\frac{\partial \mathbf{b}}{\partial t}-\nabla \times \left(\mathbf{U}\times \mathbf{b}\right)$$

$${\bf v}=\frac{\partial\mbox{\boldmath{$\xi$}}}{\partial t}+\left({\bf U}\cdot\nabla\right)\mbox{\boldmath{$\xi$}}-\left(\mbox{\boldmath{$\xi$}}\cdot\nabla\right){\bf U}.$$

$$v_r$$

$$v_{\phi }$$

$$v_z$$

$$\frac{{\partial }^2{P}^{\prime }}{\partial r^2}$$

$${\xi }_r$$

$${\xi }_{\phi }$$

$$\frac{{\mathrm{d}}^2{P}^{\prime }}{\mathrm{d}{r}^2}+\frac{1}{r}\frac{\mathrm{d}{P}^{\prime }}{\mathrm{d}r}+\left(\frac{{\Omega }^2-k_z^2{v}_{\mathrm{A}}^2}{v_{\mathrm{A}}^2}-\frac{m^2}{r^2}\right)P^{\prime }=0,$$

$$P_{\mathrm{i}}^{\prime }=A_{\mathrm{i}}{J}_m\left(k_{\perp ,\mathrm{i}}r\right),$$

$$k_{\perp ,\mathrm{i}}^2=\frac{{\Omega }_{\mathrm{i}}^2-k_z^2{v}_{\mathrm{A},\mathrm{i}}^2}{v_{\mathrm{A},\mathrm{i}}^2}.$$

$$P_{\mathrm{e}}^{\prime }=A_{\mathrm{e}}{K}_m\left(k_{\perp ,\mathrm{e}}r\right),$$

$$k_{\perp ,\mathrm{e}}^2=-\frac{{\Omega }_{\mathrm{e}}^2-k_z^2{v}_{\mathrm{A},\mathrm{e}}^2}{v_{\mathrm{A},\mathrm{e}}^2}.$$

$${\omega }^2-2k_{z}U(r_{\mathrm{A}})\omega +k_z^2\left(U^2(r_{\mathrm{A}})-v_{\mathrm{A}}^2(r_{\mathrm{A}})\right)=0.$$

$${\left(r-r_{\mathrm{A}}\right)}^{2}h(r)\frac{{\partial }^2{P}^{\prime }}{\partial r^2}+\left(r-r_{\mathrm{A}}\right)p(r)\frac{\partial P^{\prime }}{\partial r}+q(r)P^{\prime }=0,$$

$$h(r)$$

$$p(r)$$

$$q(r)$$

$$f(r)=\rho (r)\left({\Omega }^2(r)-k_z^2{v}_{\mathrm{A}}^2(r)\right)={\Omega }^2(r)\rho (r)-k_z^2\frac{B^2}{\mu }.$$

$$f(r)=\sum _{k=0}^{\infty }{f}_k{\left(r-r_{\mathrm{A}}\right)}^k,$$

$$f_k=\sum _{n=0}^k{\rho }_n\sum _{l=0}^{k-n}{\Omega }_l{\Omega }_{k-n-l}\text{for}k\ge 1,$$

$${\rho }_k$$

$${\Omega }_k$$

$$P_{\mathrm{tr}}^{\prime }(r)=A_0{P}_1^{\prime }(r)+S_0{P}_2^{\prime }(r),$$

$$P_1^{\prime }(r)$$

$$P_2^{\prime }(r)$$

$$\mathcal{C}=\frac{m^2}{2r_{\mathrm{A}}^2}.$$

$$P_{\mathrm{tr}}^{\prime }(r)$$

$${\xi }_{r,\mathrm{tr}}(r)$$

$${\xi }_{\phi ,\mathrm{tr}}(r)$$

$$\frac{\frac{k_{\perp ,\mathrm{e}}}{{\rho }_{\mathrm{e}}\left({\Omega }_{\mathrm{e}}^2-k_z^2{v}_{\mathrm{A},\mathrm{e}}^2\right)}\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}}_{+}-{\Xi }_{+}}{\frac{k_{\perp ,\mathrm{e}}}{{\rho }_{\mathrm{e}}\left({\Omega }_{\mathrm{e}}^2-k_z^2{v}_{\mathrm{A},\mathrm{e}}^2\right)}\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}}_{+}-{\Gamma }_{+}}$$

$${\mathcal{G}}_{\pm }$$

Full text

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 $\sim 5-30$ km s*-1* (e.g., Del Zanna, 2008; Winebarger et al., 2013), although velocities of $\sim 50$ km s*-1* have also been reported frequently (e.g., Brekke et al., 1997; Doschek et al., 2008), and even larger velocities up to $\sim 72-123$ km s*-1* have been observed (Ofman & Wang, 2008). Considering that the expected value of the Alfvén velocity in the solar corona is $\sim 1000$ km s*-1*, these observed flow speeds correspond to small fractions of the coronal Alfvén velocity. Observations of flows with Alfvénic speeds ($\gtrsim 500-1000$ km s*-1*) 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 Section 3. Finally, some concluding remarks and prospects for future studies are given in Section 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$, $\varphi$, 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 ${\bf B}=B{\bf 1}_{z}$. The mass density, $\rho$, is uniform in the azimuthal and longitudinal directions and nonuniform in the radial direction, so that $\rho=\rho(r)$. We considered the following radial dependence for the density,

[TABLE]

where $\rho_{\rm i}$ and $\rho_{\rm e}$ are internal and external constant densities and $\rho_{\rm tr}(r)$ is the continuous density profile that connects the internal plasma to the external plasma. We considered $\rho_{\rm i}>\rho_{\rm 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 ${\bf U}=U{\bf 1}_{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

[TABLE]

where $U_{\rm i}$ and $U_{\rm e}$ are internal and external flow velocities and $U_{\rm 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 $\beta=0$ approximation, where $\beta$ refers to the ratio of the thermal pressure to the magnetic pressure. Therefore, the basic equations used in the present work are

[TABLE]

where ${\bf v}=(v_{r},v_{\varphi},v_{z})$ is the velocity perturbation, ${\bf b}=(b_{r},b_{\varphi},b_{z})$ is the magnetic field perturbation, and $\mu$ is the magnetic permeability. In addition, the plasma Lagrangian displacement, \mbox{\boldmath{\xi}}=(\xi_{r},\xi_{\varphi},\xi_{z}), is related to the velocity perturbation and background flow by

[TABLE]

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\omega t)$, where $\omega$ the global mode frequency. In addition, we Fourier-analyzed the perturbations along the uniform $\varphi$- and $z$-directions, so that the perturbations are put proportional to $\exp(im\varphi+ik_{z}z)$, 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 Equation (5) that the relation between the components of the velocity perturbation and the Lagrangian displacement are

[TABLE]

where $\Omega(r)=\omega-k_{z}U(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 $\beta=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^{\prime}={\bf B}\cdot{\bf b}/\mu$, as our main variable. We combined Equations (3) and (4) and, after some algebraic manipulations, we obtain a differential equation involving $P^{\prime}$ alone, namely

[TABLE]

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

The components of the Lagrangian displacement are related to $P^{\prime}$ as

[TABLE]

along with $\xi_{z}=0$ because of the $\beta=0$ approximation. In the absence of flow, i.e., for $U(r)=0$ so that $\Omega(r)=\omega$, Equations (9)–(11) consistently revert to Equations (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, Equation (9) simplifies to

[TABLE]

where now $v_{\mathrm{A}}$ and $\Omega$ are constant. Equation (12) is the Bessel Equation and applies both in the internal ($r\leq R-l/2$) and external ($r\geq 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^{\prime}$ must be regular at $r=0$. Thus, the physical solution of Equation (12) is

[TABLE]

where $A_{\rm i}$ is a constant, $J_{m}$ is the Bessel function of the first kind of order $m$, and

[TABLE]

In the external plasma, we required that $P^{\prime}$ vanishes when $r\to\infty$. This is the condition for the wave to be trapped. The physical solution to Equation (12) is then

[TABLE]

where again $A_{\rm e}$ is a constant, $K_{m}$ is the modified Bessel function of the first kind of order $m$, and

[TABLE]

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=\pm 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 (Equation (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\neq 0$, Equation (9) is singular at the specific position, $r=r_{\mathrm{A}}$, where the resonant condition $\Omega^{2}(r_{\mathrm{A}})=k_{z}^{2}v_{\mathrm{A}}^{2}(r_{\mathrm{A}})$ is satisfied. Here, $r_{\mathrm{A}}$ denotes the Alfvén resonance position, which is a regular singular point of Equation (9). We can expand the resonant condition to find an implicit equation whose solution is the resonant position, namely

[TABLE]

For arbitrary density and flow profiles, Equation (17) has to be solved numerically to find $r_{\mathrm{A}}$. An analytic expression of $r_{\mathrm{A}}$ is only possible for very specific and simple profiles. We note that Equation (17) depends upon $\omega$, which shall be obtained from the dispersion relation that, in turn, requires $r_{\mathrm{A}}$ to be known. So, Equation (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 Equation (9) as an infinite power series expansion around the resonance position $r=r_{\mathrm{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 Equation (9) as

[TABLE]

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

[TABLE]

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

[TABLE]

We assumed that both the density and the flow velocity are analytic functions at $r=r_{\mathrm{A}}$. We performed a Taylor series of $f(r)$ around the location of the resonance as

[TABLE]

with $f_{0}=0$ and

[TABLE]

where $\rho_{0}=\rho(r_{\mathrm{A}})$, $\Omega_{0}=\omega-k_{z}U(r_{\mathrm{A}})$, and

[TABLE]

when $k\geq 1$.

The general solution to Equation (18) is

[TABLE]

where the subscript ‘tr’ denotes the transitional layer, $A_{0}$ and $S_{0}$ are constants, and $P^{\prime}_{1}(r)$ and $P^{\prime}_{2}(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_{\mathrm{A}}$. The indicial equation is obtained from the coefficient of the lowest power of $\left(r-r_{\mathrm{A}}\right)$, and it shows that zero and two are the two possible indices of the expansion. Then, the two linearly independent solutions are

[TABLE]

where $\mathcal{C}$ 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\left(r-r_{\mathrm{A}}\right)=\ln\left(r_{\mathrm{A}}-r\right)\pm i\pi$ if $r<r_{\mathrm{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, $\mathcal{C}$, is independent of the density and flow profiles and is given by

[TABLE]

For sausage waves, meaning when $m=0$, $\mathcal{C}=0,$ the singular logarithmic term is dropped from $P^{\prime}_{2}(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 $\mathcal{O}(l/R)$ in the Frobenius series. This is the so-called thin boundary (TB) approximation. In that simple scenario, we find

[TABLE]

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_{\mathrm{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^{\prime}$, and the radial component of the Lagrangian displacement, $\xi_{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_{\rm i}$, $A_{\rm e}$, $A_{0}$, and $S_{0}$. 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

[TABLE]

The quantities $\mathcal{G}_{\pm}$, $\mathcal{F}_{\pm}$, $\Xi_{\pm}$, and $\Gamma_{\pm}$ are defined as

[TABLE]

with $\zeta_{\pm}=R\pm\frac{l}{2}-r_{\mathrm{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_{\mathrm{A}}$, also needs to be found numerically from Equation (17). The transverse density profiles used by Soler et al. (2013) allow analytic expressions for $r_{\mathrm{A}}$, which is not possible here because of the presence of flow. Hence, the present method iteratively solves Equations (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 Equation (34) in the case of fixed, real, and positive $k_{z}$. So, we impose a particular wavelength of the perturbations, namely $\lambda=2\pi/k_{z}$. Then, the solution is a complex frequency, namely $\omega=\omega_{\rm R}+i\omega_{\rm I}$, where the subscripts R and I denote the real and imaginary parts, respectively. The sign of $\omega_{\rm R}$ indicates the direction of wave propagation. If $\omega_{\rm R}>0$ the wave propagates toward the positive $z$-direction (forward propagation), while propagation is in the opposite direction when $\omega_{\rm R}<0$ (backward propagation). The phase velocity is computed as

[TABLE]

where the sign of $v_{ph}$ follows the same rules as the sign of $\omega_{\rm R}$. On the other hand, $\omega_{\rm I}$ is related to the damping rate of the waves, so that $\omega_{\rm I}<0$ for both directions of propagation. If the waves are not overdamped, i.e., if $\left|\omega_{\rm I}\right|<\left|\omega_{\rm R}\right|$, we can define the exponential damping length of the waves, $L_{\rm 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)

[TABLE]

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_{z}R$. An appropriate way to check the validity of Equation (34) is to recover the dispersion relation previously obtained in the literature in the thin tube (TT, $k_{z}R\ll 1$) and thin boundary (TB, $l/R\ll 1$) approximations (e.g., Goossens et al., 1992; Terradas et al., 2010; Soler et al., 2011). To do so, we kept terms up to $\mathcal{O}(l/R)$ in the expressions of $\mathcal{G}_{\pm}$, $\mathcal{F}_{\pm}$, $\Xi_{\pm}$, and $\Gamma_{\pm}$, and approximate $r_{\mathrm{A}}\approx 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 $\mathcal{G}_{\pm}\approx 0$, $\mathcal{F}_{\pm}\approx 1$, $\Xi_{\pm}\approx 2/f_{1}$, and

[TABLE]

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\neq 0$ of the Bessel functions in Equation (34) (see, e.g., Abramowitz & Stegun, 1972) and also take $R-l/2\approx R+l/2\approx R$, namely

[TABLE]

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

[TABLE]

with

[TABLE]

Equation (45) agrees with Equation (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 Equation (34). We refer readers to those works to study the solutions of the TT and TB dispersion relation.

Soler et al. (2011) considered Equation (45) and obtained approximate expressions of $v_{ph}$ and $L_{\rm D}$ with $U_{\rm e}=0$ and in the case of slow flow (SF, $U_{\rm i}/v_{\mathrm{A,i}}\ll 1$), namely

[TABLE]

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

[TABLE]

Although Equations (47) and (48) are only strictly valid in the limits $k_{z}R\ll 1$, $l/R\ll 1$, and $U_{\rm i}/v_{\mathrm{A,i}}\ll 1$ they provide useful information regarding the effects of the various parameters. We note that Equations (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_{\rm D}$. We compare these approximations with the solutions of the general Equation (34), which remains valid beyond the range of applicability of Equations (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_{\rm e}=0$. In addition, we set the ratio of the internal density to the external density to $\rho_{\rm i}/\rho_{\rm e}=3$, which is representative of a coronal loop. For this density contrast, $v_{k}\approx 1.22v_{\mathrm{A,i}}$.

We investigated the behavior of $v_{ph}$ and $L_{\rm D}$ as functions of three dimensionless quantities: the ratio of the internal flow velocity to the internal Alfvén velocity, $U_{\rm i}/v_{\mathrm{A,i}}$, the relative thickness of the boundary layer, $l/R$, and dimensionless longitudinal wavenumber, $k_{z}R$. Typical values of $v_{\mathrm{A,i}}$ and $R$ in coronal loops are $v_{\mathrm{A,i}}\sim$ 1,000 km s*-1* and $R\sim$ 3,000 km.

3.1 Case without flow

To start with, we considered the case with $U_{\rm 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_{z}R$, which was not done in Soler et al. (2013). Here we have studied the impact of this parameter. Figure 1 shows surface plots of $v_{ph}/v_{\mathrm{A,i}}$ and $L_{\rm D}/R$ as functions of $l/R$ and $k_{z}R$. Concerning the phase velocity, we find that regardless the value of $k_{z}R$, 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_{z}R$. Conversely, the effect of increasing $k_{z}R$ is the opposite, in other words, the phase velocity decreases. Therefore, the parameters $l/R$ and $k_{z}R$ have competing effects on the phase velocity. When $l=0$, the phase velocity approaches the internal Alfvén velocity when $k_{z}R$ increases (see Edwin & Roberts, 1983), but when $l\neq 0$ the phase velocity tends to a value somewhat larger than $v_{\mathrm{A,i}}$ when $k_{z}R$ increases. This effect is due to the transverse nonuniformity.

On the other hand, $l/R$ and $k_{z}R$ have the same effect on the damping length, since $L_{\rm D}/R$ decreases when any of the two parameters increases. This result agrees qualitatively with the approximate analytic dependence of Equation (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_{z}R$ is to produce shorter wavelengths (i.e., higher frequencies) that are more efficiently damped.

In order to test the accuracy of the approximate Equations (47) and (48), we define the normalized errors of $v_{ph}$ and $L_{\rm D}$ associated with the use of the approximations as

[TABLE]

where $v_{ph,\rm app}$ and $L_{\rm D,app}$ denote the analytic approximations given by Equations (47) and (48), respectively, and $v_{ph,\rm Frob}$ and $L_{\rm 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_{z}R$. Very small errors are obtained when $l/R\ll 1$ and $k_{z}R\ll 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 $\sim 15\%$ and are obtained either when $l/R\to 2$ and $k_{z}R\ll 1$ or when $l/R\ll 1$ and $k_{z}R\gg 1$. On the contrary, the error in the case of the damping length seems to be more dominated by $k_{z}R$. The largest values of ${\rm err}\left(L_{\rm D}\right)$ are $\sim 70\%$ and are obtained when $k_{z}R\gg 1$, while the error when $k_{z}R\ll 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_{\rm D}$ for some particular combination of parameters far beyond the regime of applicability of the approximations.

3.2 Effect of flow

We took flow into account and exploited the new extension to the method. Figure 3 shows surface plots of $v_{ph}/v_{\mathrm{A,i}}$ and $L_{\rm D}/R$ as functions of $l/R$ and $U_{\rm i}/v_{\mathrm{A,i}}$ for $k_{z}R=0.1$. Because of the flow, different results are obtained for forward and backward waves, so Figure 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_{\rm i}/v_{\mathrm{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.

Regarding the damping length, we find that for fixed $k_{z}R$, the behavior of $L_{\rm 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_{\rm i}=0$. Hence, the backward wave is more efficiently damped than the forward wave, that is, the backward wave has a shorter $L_{\rm D}$. Qualitatively, this behavior agrees with the approximate analytic dependence of Equation (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 Figure 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 Equation (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 $\sim 13\%$ for the forward wave and $\sim 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.

Next we considered to the error of the damping length. As in the case without flow, ${\rm err}(L_{\rm 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 Equation (48) when $l/R$ is large and some specific values of the flow velocity are considered. On the other hand, the dependence of ${\rm err}(L_{\rm D})$ with $U_{\rm i}/v_{\mathrm{A,i}}$ shows a much simpler relation, namely ${\rm err}(L_{\rm D})$ typically grows when $U_{\rm i}/v_{\mathrm{A,i}}$ increases. The largest values of ${\rm err}(L_{\rm D})$ plotted in Figure 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_{\rm i}/v_{\mathrm{A,i}}$ predicted by the analytic Equations (47) and (48) as a consequence of adopting SF approximation is less accurate for $L_{\rm 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\ll 1$. We note that we computed the results for a fixed value of the density ratio $\rho_{\rm i}/\rho_{\rm 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.

Acknowledgements.

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.

Appendix A Expressions of the Frobenius series coefficients

The expressions of the Frobenius series coefficients $a_{k}$ and $s_{k}$ are as follows:

[TABLE]

In the absence of flow, $f_{k}=\omega^{2}\rho_{k}$ and the coefficients consistently revert to those given in Soler et al. (2013).