Issue 
A&A
Volume 623, March 2019



Article Number  A32  
Number of page(s)  10  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201834936  
Published online  28 February 2019 
Transverse waves in coronal flux tubes with thick boundaries: The effect of longitudinal flows
^{1}
Departament de Física, Universitat de les Illes Balears, 07122 Palma de Mallorca, Spain
email: roberto.soler@uib.es
^{2}
Institut d’Aplicacions Computacionals de Codi Comunitari (IAC 3), Universitat de les Illes Balears, 07122 Palma de Mallorca, Spain
Received:
20
December
2018
Accepted:
15
January
2019
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én 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 semianalytic 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én resonance, while the numerical part of the method consists of solving iteratively the transcendental dispersion relation together with the equation for the Alfvén 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. With the present method, 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.
Key words: magnetohydrodynamics (MHD) / Sun: atmosphere / Sun: corona / Sun: oscillations / waves
© ESO 2019
1. Introduction
Highresolution 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, longwavelength 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^{−1} (e.g., Del Zanna 2008; Winebarger et al. 2013), although velocities of ∼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 ∼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 ∼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 (≳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; TerraHomem 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 socalled 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 semianalytic 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 quasimode displacement field in a tube with a wide boundary; Arregui et al. (2015) perform Bayesian inference, model comparison, and modelaveraging techniques to infer the crossfield density structuring in coronal waveguides; Soler & Terradas (2015) compare the solution provided by the Frobeniusbased 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,
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 fieldaligned 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
where U_{i} and U_{e} are internal and external flow velocities and U_{tr}(r) denotes the spatiallydependent 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
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
From here on we adopt the socalled quasimode 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 Fourieranalyzed the perturbations along the uniform φ and zdirections, so that the perturbations are put proportional to exp(imφ + 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 Eq. (5) that the relation between the components of the velocity perturbation and the Lagrangian displacement are
where Ω(r)=ω − k_{z}U(r) is the spatiallydependent Dopplershifted frequency. We note that because of the presence of the spatiallydependent 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
where is the square of the spatiallydependent 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
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
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
where A_{i} is a constant, J_{m} is the Bessel function of the first kind of order m, and
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
where again A_{e} is a constant, K_{m} is the modified Bessel function of the first kind of order m, and
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 cutoff 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 superAlfvénic. Since the vast majority of observed flows in coronal flux tubes are subAlfvénic (see Reale 2014), for simplicity we restricted our analysis to subAlfvé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 outgoing 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 subAlfvé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 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
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 (quasimodes) 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 nonideal process such as, resistivity or viscosity. The present approach, which is based on ideal quasimodes, 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
where the functions h(r), p(r), and q(r) are defined as
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
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
with f_{0} = 0 and
where ρ_{0} = ρ(r_{A}), Ω_{0} = ω − k_{z}U(r_{A}), and
when k ≥ 1.
The general solution to Eq. (18) is
where the subscript “tr” denotes the transitional layer, A_{0} and S_{0} are constants, and and 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
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
For sausage waves, meaning when m = 0, 𝒞 = 0, the singular logarithmic term is dropped from . 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 socalled thin boundary (TB) approximation. In that simple scenario, we find
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_{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
The quantities 𝒢_{±}, ℱ_{±}, Ξ_{±}, and Γ_{±} are defined as
with . 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 zdirection (forward propagation), while propagation is in the opposite direction when ω_{R} < 0 (backward propagation). The phase velocity is computed as
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)
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 Eq. (34) is to recover the dispersion relation previously obtained in the literature in the thin tube (TT, k_{z}R ≪ 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_{1}, and
In order to make further analytic progress, we assumed that the magnetic tube is thin and use the firstorder 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
After long but straightforward algebraic manipulations, the TT and TB approximation of Eq. (34) can be cast as
with
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
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
Although Eqs. (47) and (48) are only strictly valid in the limits k_{z}R ≪ 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 socalled 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_{z}R. Typical values of v_{A,i} and R in coronal loops are v_{A,i} ∼ 1000 km s^{−1} 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_{z}R, 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_{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 ≠ 0 the phase velocity tends to a value somewhat larger than v_{A,i} when k_{z}R increases. This effect is due to the transverse nonuniformity.
Fig. 1. Results in the absence of flow. Contour plots of v_{ph}/v_{A,i} (panel a) and L_{D}/R (panel b) as functions of l/R and k_{z}R. We note that L_{D}/R is given in logarithmic scale. The dashed purple line in panel a denotes v_{ph} = v_{k}. 

Open with DEXTER 
On the other hand, l/R and k_{z}R 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_{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 Eqs. (47) and (48), we define the normalized errors of v_{ph} and L_{D} associated with the use of the approximations as
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 Frobeniusbased method. Figure 2 shows the normalized errors as functions of l/R and k_{z}R. Very small errors are obtained when l/R ≪ 1 and k_{z}R ≪ 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_{z}R ≪ 1 or when l/R ≪ 1 and k_{z}R ≫ 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 err(L_{D}) are ∼70% and are obtained when k_{z}R ≫ 1, while the error when k_{z}R ≪ 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 av_{ph} and panel bL_{D} as functions of l/R and k_{z}R. The dashed purple line in both panels denotes the contour of zero error. 

Open with DEXTER 
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_{z}R = 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 v_{ph}/v_{A,i} (panels a and b) and L_{D}/R (panels c and d) as functions of l/R and U_{i}/v_{A,i}. Panels a and c: forward propagating wave, while panels b and d are for the backward propagating wave. We note that L_{D}/R is given in logarithmic scale. We have used k_{z}R = 0.1. 

Open with DEXTER 
Regarding the damping length, we find that for fixed k_{z}R, 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 v_{ph} (panels a and b) and L_{D} (panels c and d) as functions of l/R and U_{i}/v_{A,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 k_{z}R = 0.1. 

Open with DEXTER 
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 semianalytic 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 semianalytic 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 Frobeniusbased 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 semianalytic Frobeniusbased 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 AYA201785465P (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 20132017) for the Ramón y Cajal grant RYC201414970.
References
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (New York: Dover) [Google Scholar]
 Andries, J., Goossens, M., Hollweg, J. V., Arregui, I., & Van Doorsselaere, T. 2005, A&A, 430, 1109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arregui, I., Van Doorsselaere, T., Andries, J., Goossens, M., & Kimpe, D. 2005, A&A, 441, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arregui, I., Soler, R., & Asensio Ramos, A. 2015, ApJ, 811, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Aschwanden, M. J., Nightingale, R. W., Andries, J., Goossens, M., & Van Doorsselaere, T. 2003, ApJ, 598, 1375 [NASA ADS] [CrossRef] [Google Scholar]
 Brekke, P., KjeldsethMoe, O., & Harrison, R. A. 1997, Sol. Phys., 175, 511 [NASA ADS] [CrossRef] [Google Scholar]
 Cally, P. S. 1986, Sol. Phys., 103, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Cally, P. S. 2003, Sol. Phys., 217, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Del Zanna, G. 2008, A&A, 481, L49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Doschek, G. A., Warren, H. P., Mariska, J. T., et al. 2008, ApJ, 686, 1362 [NASA ADS] [CrossRef] [Google Scholar]
 Edwin, P. M., & Roberts, B. 1983, Sol. Phys., 88, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Erdélyi, R., & Fedun, V. 2007, Science, 318, 1572 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Goddard, C. R., Pascoe, D. J., Anfinogentov, S., & Nakariakov, V. M. 2017, A&A, 605, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goddard, C. R., Antolin, P., & Pascoe, D. J. 2018, ApJ, 863, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Goossens, M., Hollweg, J. V., & Sakurai, T. 1992, Sol. Phys., 138, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Goossens, M., Andries, J., Soler, R., et al. 2012, ApJ, 753, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Goossens, M., Soler, R., Terradas, J., Van Doorsselaere, T., & Verth, G. 2014, ApJ, 788, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Guo, M.Z., Chen, S.X., Li, B., Xia, L.D., & Yu, H. 2016, Sol. Phys., 291, 877 [NASA ADS] [CrossRef] [Google Scholar]
 Holzwarth, V., Schmitt, D., & Schüssler, M. 2007, A&A, 469, 11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Innes, D. E., McKenzie, D. E., & Wang, T. 2003, Sol. Phys., 217, 267 [NASA ADS] [CrossRef] [Google Scholar]
 Ionson, J. A. 1978, ApJ, 226, 650 [NASA ADS] [CrossRef] [Google Scholar]
 Jess, D. B., Morton, R. J., Verth, G., et al. 2015, Space Sci. Rev., 190, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, M. A., & Roberts, B. 1986, ApJ, 301, 430 [NASA ADS] [CrossRef] [Google Scholar]
 Mathioudakis, M., Jess, D. B., & Erdélyi, R. 2013, Space Sci. Rev., 175, 1 [NASA ADS] [CrossRef] [Google Scholar]
 McIntosh, S. W., de Pontieu, B., Carlsson, M., et al. 2011, Nature, 475, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Morton, R. J., Tomczyk, S., & Pinto, R. F. 2016, ApJ, 828, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Nakariakov, V. M., & Roberts, B. 1995, Sol. Phys., 159, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Nakariakov, V. M., Hornsey, C., & Melnikov, V. F. 2012, ApJ, 761, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Nitta, S., Imada, S., & Yamamoto, T. T. 2012, Sol. Phys., 276, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Ofman, L., & Wang, T. J. 2008, A&A, 482, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pascoe, D. J., Wright, A. N., & De Moortel, I. 2010, ApJ, 711, 990 [NASA ADS] [CrossRef] [Google Scholar]
 Pascoe, D. J., Hood, A. W., de Moortel, I., & Wright, A. N. 2012, A&A, 539, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poedts, S., Goossens, M., & Kerner, W. 1989, Sol. Phys., 123, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Reale, F. 2014, Sol. Phys., 11, 4 [Google Scholar]
 Ruderman, M. S. 2010, Sol. Phys., 267, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. S., & Terradas, J. 2015, A&A, 580, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ryutova, M., Berger, T., Frank, Z., Tarbell, T., & Title, A. 2010, Sol. Phys., 267, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, R. 2017, ApJ, 850, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, R., & Terradas, J. 2015, ApJ, 803, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, R., Terradas, J., & Goossens, M. 2011, ApJ, 734, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, R., Goossens, M., Terradas, J., & Oliver, R. 2013, ApJ, 777, 158 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, R., Goossens, M., Terradas, J., & Oliver, R. 2014, ApJ, 781, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Stenuit, H., Tirry, W. J., Keppens, R., & Goossens, M. 1999, A&A, 342, 863 [NASA ADS] [Google Scholar]
 Tagger, M., Falgarone, E., & Shukurov, A. 1995, A&A, 299, 940 [NASA ADS] [Google Scholar]
 TerraHomem, M., Erdélyi, R., & Ballai, I. 2003, Sol. Phys., 217, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Terradas, J., & Goossens, M. 2012, A&A, 548, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Terradas, J., Oliver, R., & Ballester, J. L. 2006, ApJ, 642, 533 [NASA ADS] [CrossRef] [Google Scholar]
 Terradas, J., Goossens, M., & Ballai, I. 2010, A&A, 515, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tomczyk, S., & McIntosh, S. W. 2009, ApJ, 697, 1384 [NASA ADS] [CrossRef] [Google Scholar]
 Tomczyk, S., McIntosh, S. W., Keil, S. L., et al. 2007, Science, 317, 1192 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Van Doorsselaere, T., Andries, J., Poedts, S., & Goossens, M. 2004, ApJ, 606, 1223 [NASA ADS] [CrossRef] [Google Scholar]
 Van Doorsselaere, T., Nakariakov, V. M., & Verwichte, E. 2008, ApJ, 676, L73 [NASA ADS] [CrossRef] [Google Scholar]
 Winebarger, A., Tripathi, D., Mason, H. E., & Del Zanna, G. 2013, ApJ, 767, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Zhelyazkov, I. 2015, JApA, 36, 233 [NASA ADS] [Google Scholar]
Appendix A: Expressions of the Frobenius series coefficients
The expressions of the Frobenius series coefficients a_{k} and s_{k} are as follows:
In the absence of flow, f_{k} = ω^{2}ρ_{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 v_{ph}/v_{A,i} (panel a) and L_{D}/R (panel b) as functions of l/R and k_{z}R. We note that L_{D}/R is given in logarithmic scale. The dashed purple line in panel a denotes v_{ph} = v_{k}. 

Open with DEXTER  
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 av_{ph} and panel bL_{D} as functions of l/R and k_{z}R. The dashed purple line in both panels denotes the contour of zero error. 

Open with DEXTER  
In the text 
Fig. 3. Results in the presence of flow. Contour plots of v_{ph}/v_{A,i} (panels a and b) and L_{D}/R (panels c and d) as functions of l/R and U_{i}/v_{A,i}. Panels a and c: forward propagating wave, while panels b and d are for the backward propagating wave. We note that L_{D}/R is given in logarithmic scale. We have used k_{z}R = 0.1. 

Open with DEXTER  
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 v_{ph} (panels a and b) and L_{D} (panels c and d) as functions of l/R and U_{i}/v_{A,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 k_{z}R = 0.1. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.