RESONANT ABSORPTION OF AXISYMMETRIC MODES IN TWISTED MAGNETIC FLUX TUBES

Inhomogeneities, such as a density variation across a magnetic flux tube, produce a continuous spectrum of eigenfrequencies. For instance, consider a straight magnetic flux tube of radius r_e and constant temperature, where the density varies smoothly from its center to its boundary, such that cylindrical surfaces have constant density. This means that the sound and Alfvén speeds within every cylindrical surface are also constant. These concentric cylindrical sheaths comprise the flux tube. Because of the difference in characteristic speeds, every surface will have its own eigenfrequency. This results in an infinite set of eigenfrequencies, a continuum. One of the consequences of this continuum in driven systems is resonant absorption, assuming the driving frequency is within the continuum.

Given that inhomogeneities are the rule rather than the exception in the solar atmosphere, resonant absorption is bound to occur there. This has long been recognized, from the first suggestion by Ionson (1978) to subsequent studies motivated by advances in solar observations; see, for example, the following works: Poedts et al. (1989, 1990), Goossens et al. (2002, 2009), Ruderman & Roberts (2002), Andries et al. (2005), Van Doorsselaere et al. (2009), Terradas et al. (2010), Verth et al. (2010), Antolin et al. (2015), and Okamoto et al. (2015), to name but a few. In general, resonant absorption in magnetohydrodynamic (MHD) modes is important for the solar atmosphere. Some of the many reasons for this are the following. Resonant damping of Alfvén waves is a natural and efficient mechanism for energy dissipation of MHD waves in inhomogeneous plasmas (Ionson 1978, 1985; Hollweg & Yang 1988). It can also provide an explanation for the observed loss of power of acoustic modes in sunspots (Hollweg 1988; Sakurai et al. 1991a, 1991b; Goossens & Poedts 1992; Keppens et al. 1994), and it has been shown that it is of importance in transverse oscillations (kink mode); see, for example, Aschwanden et al. (1999), Nakariakov et al. (1999), Goossens et al. (2002), and Ruderman & Roberts (2002). Resonant Alfvén waves can be an energy conduit between photospheric motions at the footpoints of coronal loops (see, for example, De Groof & Goossens 2000, 2002; De Groof et al. 2002), and resonant dissipation plays an important role in the observed damped oscillations in prominences (see Terradas et al. 2008; Arregui et al. 2012). For an in-depth review of resonant absorption in the solar atmosphere, see Goossens et al. (2011).

Since 1999, when the first postflare, standing mode transverse oscillations were detected using the Transition Region and Coronal Explorer (TRACE) (Aschwanden et al. 1999; Nakariakov et al. 1999), there has been a growth in studies of resonant absorption for the kink mode. Ruderman & Roberts (2002) produced relations describing the expected damping for coronal loops using the long wavelength and pressureless plasma⁴ approximations, a result that was previously obtained by Goossens et al. (1992) using the connection formulae derived by Sakurai et al. (1991a, 1991b) for the driven problem and by Tirry & Goossens (1996) for the eigenvalue problem. Later, Goossens et al. (2002) and Aschwanden et al. (2003) used these results and calculated the expected damping times for a sequence of observed parameters for coronal flux tubes. Goossens et al. (2002) concluded that, for the parameter sample used, resonant absorption can explain the observed damping times well, provided that the density contrast is allowed to vary from loop to loop. Another important result in this work is that the observed damping does not require modification of the order-of-magnitude estimates of the Reynolds number (10¹⁴), as suggested by Nakariakov et al. (1999). Aschwanden et al. (2003) also arrived at the conclusion that, on average, the theoretical predictions of the damping rate derived by Goossens et al. (1992) and Ruderman & Roberts (2002) are consistent with observations and suggested that damping times of coronal loops can be used to infer their density contrast with the surrounding plasma. Coronal flux tubes tend to deform in their middle section because of buoyancy, effectively resulting in cross sections that are approximately elliptical. Ruderman (2003) studied the damping of the kink mode in flux tubes with an elliptical cross section and found that, for moderate ratios of the minor to major semiaxis, the difference of the damping rate for resonant absorption compared with flux tubes with circular cross section is not very large. Another deviation from the ideal straight magnetic flux tube is axial curvature. Van Doorsselaere et al. (2004) studied the effect of this curvature and also found that the longitudinal curvature of flux tubes does not significantly alter the damping time of kink modes. The theoretical models for kink oscillations have become progressively more elaborate; for example, Andries et al. (2005) considered longitudinal density stratification. Also, methods for kink wave excitation have been studied; see, for example, Terradas (2009). The larger body of observations of kink waves allowed Verwichte et al. (2013) to perform a statistical study to constrain the free parameters present in theoretical models of resonant absorption in kink modes.

In contrast to this avalanche of theoretical and observational advances related to the kink mode, the resonant absorption for axisymmetric modes has not received much attention. One reason for this is that it was believed that the sausage mode had a long wavelength cutoff (e.g., Edwin & Roberts 1983), which suggested that observation of the sausage mode would be quite challenging. Furthermore, it was correctly believed that, for a straight magnetic field, axisymmetric modes could not be resonantly damped. However, it is apparent, even in early works in resonant absorption (see, for example, Sakurai et al. 1991a, 1991b; Goossens et al. 1992), that, for a weakly twisted magnetic field, axisymmetric modes can and are resonantly damped. What was not known until recently, however, was that the long wavelength cutoff for these modes is also removed in the presence of a weak magnetic twist (Giagkiozis et al. 2015). Therefore, these modes can freely propagate for all wavelengths. And so, at least in principle, these modes should be observable. Additionally, recent works suggest that magnetic twist and axisymmetric modes are ubiquitous throughout the solar atmosphere. Therefore, the study of these modes has become quite relevant and important. Some examples of magnetic twist in the solar atmosphere are flux tubes emerging from the convection zone (see, for example, Hood et al. 2009; Luoni et al. 2011) and sunspot rotation resulting in twisted magnetic fields (Brown et al. 2003; Yan & Qu 2007; Kazachenko et al. 2009). Spicules are observed to have twist (De Pontieu et al. 2012; Sekse et al. 2013), as well as solar tornadoes (Wedemeyer-Böhm et al. 2012). Lastly, observations of axisymmetric modes have been recently reported in Morton et al. (2012) and Grant et al. (2015).

In this work, we focus on the resonant absorption of axisymmetric MHD modes in weakly twisted magnetic flux tubes. Axisymmetric modes correspond to modes with azimuthal wavenumber $m=0$. We accomplish this using the following sequence. First, we recall recent results for axisymmetric modes in magnetic flux tubes with weak twist (Giagkiozis et al. 2015). In that work, the longitudinal component of the magnetic field and the density were discontinuous across the flux tube boundary. This choice was intentional as it avoids the MHD continua and simplifies the analysis. However, this also left out relevant physics. Then having as a starting point the setup in Giagkiozis et al. (2015), we introduce an intermediate layer about the flux tube boundary. Within this layer, we allow the magnetic field and density to vary smoothly, resulting in an overall continuous profile for the longitudinal magnetic field and density. This in turn allows for the existence of the two MHD continua, the slow and the Alfvén continuum. Next, we assume that the layer that connects the internal and external quantities is thin; specifically we assume that ${\ell }\ll {r}_{e}$ where ℓ is the width of the layer and r_e is the flux tube radius. Then we use the conservation laws, and the resulting jump conditions, for the Alfvén continuum by Sakurai et al. (1991a), and we derive the resulting complex dispersion relation. We then solve this dispersion relation numerically. Lastly, to better understand the predicted damping times, we apply the long wavelength limit approximation to the resulting complex dispersion relation. These simpler relations allow us to compare our results with the expected damping for the kink mode predicted using the results by Goossens et al. (1992) and Ruderman & Roberts (2002). We conclude this investigation with a statistical analysis of the resulting approximations to further understand the necessary conditions for the observation of resonantly damped axisymmetric modes. The main contributions of this work can be summarized as follows:

• 1.  
  For the first time, we uncover a dispersion relation for axisymmetric modes in magnetic flux tubes with internal and external twist, including the resonance with the Alfvén continuum. We produce simplified expressions for the frequency and damping time in the long wavelength limit, for which the axisymmetric modes are no longer leaky.
• 2.  
  Given that there are four parameters required for the evaluation of the aforementioned relation, namely density contrast, magnetic field contrast, thickness of the inhomogeneous layer, and magnetic twist, we present a statistical framework to infer what can be drawn from observations.
• 3.  
  We use this statistical framework and show that the predictions of our theoretical model are in agreement with observed damping times that are in agreement with the observed damping times of quasiperiodic pulsations (QPPs). QPPs are interpreted as axisymmetric modes (sausage modes) (Kolotkov et al. 2015).

The plan of this paper is as follows. In Section 2 we present the model and include prior theoretical results required for the derivation of the dispersion relation leading to resonant absorption. In Section 3, using the jump relations in Sakurai et al. (1991a), we derive a dispersion equation. In Section 4 we use the dispersion relation derived in Section 3 to obtain an expression for the damping rate in the long wavelength limit, and then in Section 5 we elaborate on the significance of the results in this work for the observation of axisymmetric modes in the solar atmosphere. Lastly, in Section 6 we summarize and conclude this work.

In this work we assume an idealized, cylindrically symmetric magnetic flux tube in static equilibrium. We employ cylindrical coordinates $r,\varphi$, and z, with the z coordinate along the axis of symmetry of the flux tube. The linearized ideal MHD equations are

$$\begin{eqnarray}\rho \displaystyle \frac{{\partial }^{2}{\boldsymbol{\xi }}}{\partial {t}^{2}} & + & {\rm{\nabla }}{p}^{\prime }\\ & + & \displaystyle \frac{1}{{\mu }_{0}}({\boldsymbol{B}}^{\prime} \times ({\rm{\nabla }}\times {\boldsymbol{B}})+{\boldsymbol{B}}\times ({\rm{\nabla }}\times {\boldsymbol{B}}^{\prime} ))=0,\qquad \end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&{p}^{\prime }+{\boldsymbol{\xi }}\cdot {\rm{\nabla }}p+\gamma p{\rm{\nabla }}\cdot {\boldsymbol{\xi }}=0,\qquad \end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray}&&{\boldsymbol{B}}^{\prime} +{\rm{\nabla }}\times ({\boldsymbol{B}}\times {\boldsymbol{\xi }})=0,\qquad \end{eqnarray} \tag{ 1c }$$

where $\rho ,p$, and ${\boldsymbol{B}}$ are the density, plasma kinetic pressure, and magnetic field, respectively, at equilibrium, ${\boldsymbol{\xi }}$ is the Lagrangian displacement, ${p}^{\prime }$ and ${\boldsymbol{B}}^{\prime}$ are the Eulerian variations of the pressure and magnetic field, γ is the ratio of specific heats (taken to be 5/3 in this work), and ${\mu }_{0}$ is the permeability of free space. In what follows, an index, i, indicates quantities inside the flux tube ($r\lt {r}_{i}$), while variables indexed by e refer to the environment outside the flux tube ($r\gt {r}_{e}$). The inhomogeneous layer has a width equal to ${\ell }={r}_{e}-{r}_{i}$, and it is assumed that ${\ell }\ll {r}_{e}$. Note that in Giagkiozis et al. (2015) r_a was used to denote the tube radius; this is equivalent to r_e in this work. The model configuration is illustrated in Figure 1 when ${B}_{\varphi e}\propto 1/r$. The quantities $\rho ,p$, and ${\boldsymbol{B}}$ are assumed to have only an r dependence; therefore, the following balance equation must be satisfied when ${\ell }=0$:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dr}}\left(p+\displaystyle \frac{{B}_{\varphi }^{2}+{B}_{z}^{2}}{2{\mu }_{0}}\right)=-\displaystyle \frac{{B}_{\varphi }^{2}}{{\mu }_{0}r}.\end{eqnarray} \tag{ 2 }$$

The equilibrium magnetic field is taken to be ${\boldsymbol{B}}=(0,{B}_{\varphi },{B}_{z})$, with ${B}_{\varphi i}={Sr}$, ${B}_{\varphi e}={r}_{e}^{1+\kappa }S/{r}^{\kappa }$, and ${B}_{{zi}},{B}_{{ze}}$ constant. By substituting ${B}_{\varphi i}$ and ${B}_{\varphi e}$ into Equation (2) and defining ${B}_{\varphi A}={B}_{\varphi }({r}_{e})={{Sr}}_{e}$, we obtain

$$\begin{eqnarray}p(r)=\left\{\begin{array}{ll}\displaystyle \frac{{B}_{\varphi A}^{2}}{{\mu }_{0}}\left(1-\displaystyle \frac{{r}^{2}}{{r}_{e}^{2}}\right)+{p}_{e} & {\rm{for}}\ r\leqslant {r}_{e},\\ \displaystyle \frac{{r}_{e}^{2\kappa }{B}_{\varphi A}^{2}(1-\kappa )}{2{\mu }_{0}\kappa }\left(\displaystyle \frac{1}{{r}^{2\kappa }}-\displaystyle \frac{1}{{r}_{e}^{2\kappa }}\right)+{p}_{e} & {\rm{for}}\ r\gt {r}_{e},\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where p_e is the pressure at the boundary of the magnetic flux tube, and the parameter $\kappa \to 1$ corresponds to an external twist proportional to $1/r$ while $\kappa \to 0$ to a constant external twist. Note that although p(r) is continuous, for solar atmospheric conditions and for weak magnetic twist ($\mathrm{sup}({B}_{\varphi }^{2}/{B}_{z}^{2})\ll 1$), its variation is much smaller than p_e and therefore can be assumed to be constant (Giagkiozis et al. 2015). However, in the model used by Giagkiozis et al. (2015), the equilibrium density and the z component of the magnetic field are discontinuous, so the Alfvén continuum was avoided. Note that in Giagkiozis et al. (2015) the equivalent to Equation (3) had a typographical error: $(1-2\kappa )$ should read $(1-\kappa )$.

Zoom In Zoom Out Reset image size

In the present investigation, both the density and the magnetic field are continuous; see Figure 2, which introduces the slow and fast continua into our model. Specifically, the density is assumed to be a piecewise linear function of the following form:

$$\begin{eqnarray}\rho (r)=\left\{\begin{array}{ll}{\rho }_{i} & {\rm{for}}\ r\lt {r}_{i},\\ {\rho }_{i}+\displaystyle \frac{r-{r}_{i}}{{\ell }}({\rho }_{e}-{\rho }_{i}) & {\rm{for}}\ {r}_{i}\leqslant r\leqslant {r}_{e},\\ {\rho }_{e} & {\rm{for}}\ r\gt {r}_{e},\end{array}\right.\end{eqnarray} \tag{ 4 }$$

and a similar form for the variation in the longitudinal component of the magnetic field is assumed:

$$\begin{eqnarray}{B}_{z}(r)=\left\{\begin{array}{ll}{B}_{{zi}} & {\rm{for}}\ r\lt {r}_{i},\\ {B}_{{zi}}+\displaystyle \frac{r-{r}_{i}}{{\ell }}({B}_{{ze}}-{B}_{{zi}}) & {\rm{for}}\ {r}_{i}\leqslant r\leqslant {r}_{e},\\ {B}_{{ze}} & {\rm{for}}\ r\gt {r}_{e}.\end{array}\right.\end{eqnarray} \tag{ 5 }$$

Note that the assumption here is that ${\ell }\ll {r}_{e}$, so pressure balance is maintained (see Equation (2)). Also note that allowing both the density and the magnetic field to vary results in a nonmonotonic variation in the Alfvén frequency across the inhomogeneous layer, as seen in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The equilibrium quantities depend only on r, and therefore the perturbed quantities can be Fourier-analyzed with respect to the φ and z coordinates:

$$\begin{eqnarray}&&{\boldsymbol{\xi }},{p}_{T}^{\prime }\propto {e}^{i(m\varphi +{k}_{z}z-\omega t)}.\end{eqnarray} \tag{ 6 }$$

Here, ω is the angular frequency, $m$ is the azimuthal wavenumber, ${k}_{z}$ is the longitudinal wavenumber, and ${p}_{T}^{\prime }$ is the Eulerian total pressure perturbation, defined as ${p}^{\prime }+{\boldsymbol{B}}{\boldsymbol{B}}^{\prime} /{\mu }_{0}$. Our focus is on axisymmetric modes (sausage waves), so the azimuthal wavenumber is taken to be $m=0$. The Lagrangian displacement vector in flux coordinates is ${\boldsymbol{\xi }}=({\xi }_{r},{\xi }_{\perp },{\xi }_{\parallel })$ where

$$\begin{eqnarray}&&{\xi }_{\perp }=\displaystyle \frac{{B}_{z}{\xi }_{\varphi }-{B}_{\varphi }{\xi }_{z}}{| {\boldsymbol{B}}| },{\xi }_{\parallel }=\displaystyle \frac{{B}_{\varphi }{\xi }_{\varphi }+{B}_{z}{\xi }_{z}}{| {\boldsymbol{B}}| },\end{eqnarray} \tag{ 7 }$$

assuming ${B}_{r}=0$. Using Equation (6), Equation (1) can be transformed into the following two coupled first-order differential equations:

$$\begin{eqnarray}&&D\displaystyle \frac{d(r{\xi }_{r})}{{dr}}={C}_{1}(r{\xi }_{r})-{{rC}}_{2}{p}_{T}^{\prime },\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&D\displaystyle \frac{{{dp}}_{T}^{\prime }}{{dr}}=\displaystyle \frac{1}{r}{C}_{3}(r{\xi }_{r})-{C}_{1}{p}_{T}^{\prime }.\end{eqnarray} \tag{ 8b }$$

$$\begin{eqnarray}&&\rho ({\omega }^{2}-{\omega }_{A}^{2}){\xi }_{\perp }=\displaystyle \frac{{\imath }}{| {\boldsymbol{B}}| }{C}_{A},\end{eqnarray} \tag{ 8c }$$

$$\begin{eqnarray}&&\rho ({\omega }^{2}-{\omega }_{c}^{2}){\xi }_{\parallel }=\displaystyle \frac{{\imath }{f}_{B}}{| {\boldsymbol{B}}| }\displaystyle \frac{{v}_{s}^{2}}{{v}_{s}^{2}+{v}_{A}^{2}}{C}_{S},\end{eqnarray} \tag{ 8d }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{\xi }}=-\displaystyle \frac{{\omega }^{2}{C}_{S}}{\rho ({v}_{s}^{2}+{v}_{A}^{2})({\omega }^{2}-{\omega }_{c}^{2})}\end{eqnarray} \tag{ 8e }$$

and

$$\begin{eqnarray}&&D=\rho ({\omega }^{2}-{\omega }_{A}^{2}){C}_{4},\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&{C}_{1}=\displaystyle \frac{2{B}_{\varphi }}{{\mu }_{0}r}\left({\omega }^{4}{B}_{\varphi }-\displaystyle \frac{m}{r}{f}_{B}{C}_{4}\right),\end{eqnarray} \tag{ 9b }$$

$$\begin{eqnarray}&&{C}_{2}={\omega }^{4}-\left({k}_{z}^{2}+\displaystyle \frac{{m}^{2}}{{r}^{2}}\right){C}_{4},\end{eqnarray} \tag{ 9c }$$

$$\begin{eqnarray}{C}_{3} & = & \rho D\left[{\omega }^{2}-{\omega }_{A}^{2}+\displaystyle \frac{2{B}_{\varphi }}{{\mu }_{0}\rho }\displaystyle \frac{d}{{dr}}\left(\displaystyle \frac{{B}_{\varphi }}{r}\right)\right]\\ & & +4{\omega }^{4}\displaystyle \frac{{B}_{\varphi }^{4}}{{\mu }_{0}^{2}{r}^{2}}-\rho {C}_{4}\displaystyle \frac{4{B}_{\varphi }^{2}{\omega }_{A}^{2}}{{\mu }_{0}{r}^{2}},\end{eqnarray} \tag{ 9d }$$

$$\begin{eqnarray}&&{C}_{4}=({v}_{s}^{2}+{v}_{A}^{2})({\omega }^{2}-{\omega }_{c}^{2}),\end{eqnarray} \tag{ 9e }$$

$$\begin{eqnarray}&&{C}_{A}={g}_{B}{p}_{T}^{\prime }-2\displaystyle \frac{{f}_{B}{B}_{\varphi }{B}_{z}{\xi }_{r}}{{\mu }_{0}r},\ {C}_{S}={p}_{T}^{\prime }-2\displaystyle \frac{{B}_{\varphi }^{2}{\xi }_{r}}{{\mu }_{0}r}\end{eqnarray} \tag{ 9f }$$

where

$$\begin{eqnarray*}{v}_{s}^{2} & = & \gamma \displaystyle \frac{p}{\rho },\ \ {v}_{A}^{2}=\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{{\mu }_{0}\rho },\\ {\omega }_{c}^{2} & = & \displaystyle \frac{{v}_{s}^{2}}{{v}_{A}^{2}+{v}_{s}^{2}}{\omega }_{A}^{2},\ \ {\omega }_{A}^{2}=\displaystyle \frac{{f}_{B}^{2}}{{\mu }_{0}\rho },\\ {f}_{B} & = & {\boldsymbol{k}}\cdot {\boldsymbol{B}}=\displaystyle \frac{m}{r}{B}_{\varphi }+{k}_{z}{B}_{z},\\ {g}_{B} & = & {({\boldsymbol{k}}\times {\boldsymbol{B}})}_{r}=\displaystyle \frac{m}{r}{B}_{z}-{k}_{z}{B}_{\varphi }.\end{eqnarray*}$$

Here, ${\boldsymbol{k}}=(0,m/r,{k}_{z})$ is the wavevector, C_A and C_S are the coupling functions, ${v}_{s}$ is the sound speed, ${v}_{A}$ is the Alfvén speed, ${\omega }_{c}$ is the cusp angular frequency, and ${\omega }_{A}$ is the Alfvén angular frequency. Equation (8) was initially derived by Hain & Lust (1958) and later by Goedbloed (1971) and Sakurai et al. (1991a). The first-order coupled ordinary differential equation (ODEs) in Equation (8) can be reduced to a single second-order ODE for ${\xi }_{r}$:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dr}}\left[\displaystyle \frac{D}{{{rC}}_{2}}\displaystyle \frac{d}{{dr}}(r{\xi }_{r})\right]\\ &&+\left[\displaystyle \frac{1}{D}\left({C}_{3}-\displaystyle \frac{{C}_{1}^{2}}{{C}_{2}}\right)-r\displaystyle \frac{d}{{dr}}\left(\displaystyle \frac{{C}_{1}}{{{rC}}_{2}}\right)\right]{\xi }_{r}=0.\end{eqnarray} \tag{ 10 }$$

The assumption of axisymmetry ($m=0$) leads to

$$\begin{eqnarray}&&{f}_{B}={k}_{z}{B}_{z},{g}_{B}=-{k}_{z}{B}_{\varphi },{C}_{A}=-{k}_{z}{B}_{\varphi }\left({p}_{T}^{\prime }+2\displaystyle \frac{{B}_{z}^{2}}{{\mu }_{0}r}{\xi }_{r}\right).\end{eqnarray} \tag{ 11 }$$

Therefore,

$$\begin{eqnarray}&&\rho ({\omega }^{2}-{\omega }_{A}^{2}){\xi }_{\perp }=-{\imath }\displaystyle \frac{{k}_{z}{B}_{\varphi }}{| {\boldsymbol{B}}| }\left({p}_{T}^{\prime }+2\displaystyle \frac{{B}_{z}^{2}}{{\mu }_{0}r}{\xi }_{r}\right),\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&\rho ({\omega }^{2}-{\omega }_{c}^{2}){\xi }_{\parallel }={\imath }\displaystyle \frac{{k}_{z}{B}_{z}}{| {\boldsymbol{B}}| }\displaystyle \frac{{v}_{s}^{2}}{{v}_{s}^{2}+{v}_{A}^{2}}\left({p}_{T}^{\prime }-2\displaystyle \frac{{B}_{\varphi }^{2}}{{\mu }_{0}r}{\xi }_{r}\right).\end{eqnarray} \tag{ 12b }$$

Note that Equation (12) suggests that the solutions for the components of the Lagrangian displacement vector are coupled. Coupled is in the sense that elimination of one component, such as by setting it to be identical to zero, has direct implications for the remaining components. To see this, consider a solution for which ${\xi }_{r}=0$; then by Equation (10), ${p}_{T}^{\prime }$ must also be equal to zero. As a consequence of Equations (12a) and (12b), it follows immediately that ${\xi }_{\perp }$ and ${\xi }_{\parallel }$ must also be identically equal to zero: setting ${\xi }_{r}=0$ leads to the trivial solution. Alternatively, let us assume that ${\xi }_{\perp }=0$. In this case, by Equation (12a) the following relation must hold:

$$\begin{eqnarray}&&{p}_{T}^{\prime }=-2\displaystyle \frac{{B}_{z}^{2}}{{\mu }_{0}r}{\xi }_{r}.\end{eqnarray} \tag{ 13 }$$

This in turn implies

$$\begin{eqnarray}&&\rho ({\omega }^{2}-{\omega }_{c}^{2}){\xi }_{\parallel }=-2{\imath }\displaystyle \frac{{k}_{z}{B}_{z}| {\boldsymbol{B}}| }{{\mu }_{0}r}\displaystyle \frac{{v}_{s}^{2}}{{v}_{s}^{2}+{v}_{A}^{2}}{\xi }_{r},\end{eqnarray} \tag{ 14 }$$

which in general is nonzero. Now, if we assume ${\xi }_{\parallel }=0$, then

$$\begin{eqnarray}&&{p}_{T}^{\prime }=2\displaystyle \frac{{B}_{\varphi }^{2}}{{\mu }_{0}r}{\xi }_{r},\end{eqnarray} \tag{ 15 }$$

which leads to

$$\begin{eqnarray}&&\rho ({\omega }^{2}-{\omega }_{A}^{2}){\xi }_{\perp }=-2{\imath }\displaystyle \frac{{k}_{z}{B}_{\varphi }| {\boldsymbol{B}}| }{{\mu }_{0}r}{\xi }_{r}.\end{eqnarray} \tag{ 16 }$$

In the case where ${B}_{\varphi }=0$, ${\xi }_{\perp }$ decouples from ${\xi }_{r}$ and ${\xi }_{\parallel }$. At this point it is instructive to mention the interpretation of the three components of ${\boldsymbol{\xi }}$ in flux coordinates by Goossens et al. (2011). Goossens et al. (2011) suggest that ${\xi }_{\perp }$ is the dominant component for Alfvén waves, and for low plasma-β the slow and fast magnetoacoustic waves ${\xi }_{\parallel }$ and ${\xi }_{r}$ are the dominant components, respectively. A quick check, by setting ${B}_{\varphi }=0$ in Equation (7), renders ${\xi }_{\perp }$ equivalent to ${\xi }_{\varphi }$. This illuminates the connection of ${\xi }_{\perp }$ with torsional Alfvén waves.

Giagkiozis et al. (2015) solved Equation (10) for weak internal and external magnetic twists, albeit with the density profile assumed to be piecewise constant. With the help of the conservation relations for the Alfvén continuum derived by Sakurai et al. (1991a), these solutions, which are for ideal MHD, can be used to produce a dispersion relation for MHD waves that undergo damping in the continuum. The solutions by Giagkiozis et al. (2015) are as follows:

$$\begin{eqnarray}&&{\xi }_{{ri}}(s)={A}_{i}\displaystyle \frac{{s}^{1/2}}{{E}^{1/4}}{e}^{-s/2}M(a,b;s),\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}{p}_{{Ti}}^{\prime }(s) & = & {A}_{i}\displaystyle \frac{{k}_{a}{D}_{i}}{{n}_{i}^{2}-{k}_{z}^{2}}{e}^{-s/2}\left[\displaystyle \frac{{n}_{i}+{k}_{z}}{{k}_{z}}{sM}(a,b;s)\right.\\ & & -2M(a,b-1;s)],\end{eqnarray} \tag{ 17b }$$

and

$$\begin{eqnarray}&&{\xi }_{{re}}(r)={A}_{e}{K}_{\nu }({k}_{{re}}r),\end{eqnarray} \tag{ 18a }$$

$$\begin{eqnarray}{p}_{{Te}}^{\prime } & = & {A}_{e}\left(\displaystyle \frac{{\mu }_{0}(1-\nu ){D}_{e}-2{B}_{\varphi A}^{2}{n}_{e}^{2}}{{\mu }_{0}r({k}_{z}^{2}-{n}_{e}^{2})}{K}_{\nu }({k}_{{re}}r)\right.\\ & & \left.-\displaystyle \frac{{D}_{e}}{{k}_{{re}}}{K}_{\nu -1}({k}_{{re}}r)\right).\end{eqnarray} \tag{ 18b }$$

where $M(\cdot )$ is the Kummer function, and $K(\cdot )$ is the modified Bessel function of the second kind (Abramowitz & Stegun 2012). The solutions in Equations (17a) and (17b) were initially derived by Erdélyi & Fedun (2007). The parameters in Equations (17) and (18) are

$$\begin{eqnarray}&&a=1+\displaystyle \frac{{k}_{{ri}}^{2}}{4{k}_{z}^{2}{E}^{1/2}},\ b=2,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{k}_{a}={k}_{z}{(1-{\alpha }^{2})}^{1/2},\ {\alpha }^{2}=\displaystyle \frac{4{B}_{\varphi A}^{2}{\omega }_{{Ai}}^{2}}{{\mu }_{0}{r}_{e}^{2}{\rho }_{i}{({\omega }^{2}-{\omega }_{{Ai}}^{2})}^{2}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&s={k}_{a}^{2}{E}^{1/2}{r}^{2},\ E=\displaystyle \frac{4{B}_{\varphi A}^{4}{n}_{i}^{2}}{{\mu }_{0}^{2}{r}_{e}^{4}{D}_{i}^{2}{k}_{z}^{2}{(1-{\alpha }^{2})}^{2}},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{k}_{r}^{2}={k}_{z}^{2}\left(1-\displaystyle \frac{{n}^{2}}{{k}_{z}^{2}}\right),\ {k}_{r}^{2}=\displaystyle \frac{({k}_{z}^{2}{v}_{s}^{2}-{\omega }^{2})({k}_{z}^{2}{v}_{A}^{2}-{\omega }^{2})}{({v}_{A}^{2}+{v}_{s}^{2})({k}_{z}^{2}{v}_{T}^{2}-{\omega }^{2})},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{n}^{2}={k}_{z}^{2}\displaystyle \frac{{\omega }^{4}}{({\omega }_{s}^{2}+{\omega }_{A}^{2})({\omega }^{2}-{\omega }_{c}^{2})},\ {v}_{T}^{2}=\displaystyle \frac{{v}_{A}^{2}{v}_{s}^{2}}{{v}_{A}^{2}+{v}_{s}^{2}},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{D}_{i}={\rho }_{i}({\omega }^{2}-{\omega }_{{Ai}}^{2}),\ {D}_{e}={\rho }_{e}({\omega }^{2}-{\omega }_{{Ae}}^{2}).\end{eqnarray} \tag{ 24 }$$

and ν is

$$\begin{eqnarray}&&{\nu }^{2}(\kappa ;r)\\ &&\;=1+2\displaystyle \frac{{r}_{e}^{2\kappa }{B}_{\varphi A}^{2}}{{\mu }_{0}^{2}{D}_{e}^{2}{r}^{2\kappa }}\left\{2\displaystyle \frac{{r}_{e}^{2\kappa }{B}_{\varphi A}^{2}{n}_{e}^{2}{k}_{z}^{2}}{{r}^{2\kappa }}+{\mu }_{0}{\rho }_{e}[{\omega }_{{Ae}}^{2}({n}_{e}^{2}(3+\kappa )\right.\\ &&\quad -{k}_{z}^{2}(1-\kappa ))-({n}_{e}^{2}+{k}_{z}^{2})(1+\kappa ){\omega }^{2}]\}.\end{eqnarray} \tag{ 25 }$$

This function in Giagkiozis et al. (2015) is evaluated for $\kappa =0$, resulting in an exact solution for constant twist outside the flux tube that is also a zero-order approximation for the external solution when the magnetic twist is proportional to $1/r$:

$$\begin{eqnarray}{\nu }^{2}(0;r) & = & 1+2\displaystyle \frac{{B}_{\varphi A}^{2}}{{\mu }_{0}^{2}{D}_{e}^{2}}\{2{B}_{\varphi A}^{2}{n}_{e}^{2}{k}_{z}^{2}\\ & & +{\mu }_{0}{\rho }_{e}[{\omega }_{{Ae}}^{2}(3{n}_{e}^{2}-{k}_{z}^{2})-{\omega }^{2}({n}_{e}^{2}+{k}_{z}^{2})]\}.\end{eqnarray} \tag{ 26 }$$

Using $\nu =\nu (0;r)$, i.e., a constant external magnetic twist, results in solutions, namely Equations (18a) and (18b), that have approximately 5% root-mean-square error when compared with the exact solution corresponding to $\nu =\nu (1;r)$, which corresponds to external magnetic twist $\sim 1/r$. For more details, see Giagkiozis et al. (2015).

Imposing continuity for the Lagrangian displacement in the radial direction and total pressure continuity across the flux tube,

$$\begin{eqnarray}&&{{\xi }_{{ri}}| }_{r={r}_{e}}={{\xi }_{{re}}| }_{r={r}_{e}},\end{eqnarray} \tag{ 27a }$$

$$\begin{eqnarray}&&{p}_{{Ti}}^{\prime }-{\left.\displaystyle \frac{{B}_{\varphi i}^{2}}{{\mu }_{0}r}{\xi }_{{ri}}\right|}_{r={r}_{e}}={p}_{{Te}}^{\prime }-{\left.\displaystyle \frac{{B}_{\varphi e}^{2}}{{\mu }_{0}r}{\xi }_{{re}}\right|}_{r={r}_{e}},\end{eqnarray} \tag{ 27b }$$

the following dispersion relation was derived:

$$\begin{eqnarray}\displaystyle \frac{{r}_{e}{D}_{e}}{{k}_{{re}}}\displaystyle \frac{{K}_{\nu -1}({k}_{{re}}{r}_{e})}{{K}_{\nu }({k}_{{re}}{r}_{e})} & = & {\rho }_{i}{v}_{A\varphi i}^{2}\left[\displaystyle \frac{1}{{k}_{{ri}}^{2}}{({n}_{i}+{k}_{z})}^{2}-\displaystyle \frac{1}{{k}_{{re}}^{2}}({n}_{e}^{2}+{k}_{z}^{2})\right]\\ & & +\displaystyle \frac{(1-\nu ){D}_{e}}{{k}_{{re}}^{2}}-2\displaystyle \frac{{D}_{i}}{{k}_{{ri}}^{2}}\displaystyle \frac{M(a,b-1;s)}{M(a,b;s)},\end{eqnarray} \tag{ 28 }$$

where ${v}_{A\varphi i}^{2}={B}_{\varphi A}^{2}/{\mu }_{0}{\rho }_{i}$ and ${\omega }_{A\varphi i}^{2}={k}_{z}^{2}{B}_{\varphi A}^{2}/{\mu }_{0}{\rho }_{i}$.

2.1. Long Wavelength Limit

The long wavelength limit of Equation (28) is needed for the approximation of the location of the resonant point used in subsequent sections and is obtained as follows. From Equation (13.5.5) in Abramowitz & Stegun (2012) we have

$$\begin{eqnarray}&&\underset{\epsilon \to 0}{\mathrm{lim}}\displaystyle \frac{M(a,b-1;s)}{M(a,b;s)}=1,\end{eqnarray} \tag{ 29 }$$

and furthermore, rewriting ${\nu }^{2}(0;r)$ as

$$\begin{eqnarray}{\nu }^{2}(0;r) & = & 1+2\displaystyle \frac{{\rho }_{e}{B}_{\varphi A}^{2}}{{\mu }_{0}{D}_{e}^{2}}\left({\omega }_{{Ae}}^{2}\left(3\displaystyle \frac{{n}_{e}^{2}}{{k}_{z}^{2}}-1\right)-{\omega }^{2}\left(\displaystyle \frac{{n}_{e}^{2}}{{k}_{z}^{2}}+1\right)\right){k}_{z}^{2}\\ & & +4\displaystyle \frac{{B}_{\varphi A}^{4}}{{\mu }_{0}^{2}{D}_{e}^{2}}\displaystyle \frac{{n}_{e}^{2}}{{k}_{z}^{2}}{k}_{z}^{4},\end{eqnarray} \tag{ 30 }$$

it becomes apparent that $\nu =1+{ \mathcal O }({\epsilon }^{2})$, where $\epsilon ={r}_{e}{k}_{z}$. Therefore, using Equations (9.6.8) and (9.6.9) in Abramowitz & Stegun (2012), we obtain that

$$\begin{eqnarray}&&\underset{\epsilon \to 0}{\mathrm{lim}}\displaystyle \frac{{K}_{0}({k}_{{re}}{r}_{e})}{{K}_{1}({k}_{{re}}{r}_{e})}=0.\end{eqnarray} \tag{ 31 }$$

Using Equations (29) and (31) in Equation (28) we have

$$\begin{eqnarray}&&2({\omega }^{2}-{\omega }_{{Ai}}^{2})={v}_{A\varphi i}^{2}\left[{({n}_{i}+{k}_{z})}^{2}-\displaystyle \frac{{k}_{{ri}}^{2}}{{k}_{{re}}^{2}}({n}_{e}^{2}-{k}_{z}^{2})\right].\end{eqnarray} \tag{ 32 }$$

Expanding the part in square brackets on the right-hand side of this equation about $\epsilon =0$ leads to

$$\begin{eqnarray}&&{({n}_{i}+{k}_{z})}^{2}-\displaystyle \frac{{k}_{{ri}}^{2}}{{k}_{{re}}^{2}}({n}_{e}^{2}-{k}_{z}^{2})=2\displaystyle \frac{\omega }{{({v}_{{Ai}}^{2}+{v}_{{si}}^{2})}^{1/2}}{k}_{z}+{ \mathcal O }({\epsilon }^{2}).\end{eqnarray} \tag{ 33 }$$

Using this approximation in Equation (32), the positive solution of the dispersion relation Equation (28) in the long wavelength limit to first order is

$$\begin{eqnarray}&&\omega =\displaystyle \frac{1}{2}\left[\displaystyle \frac{{\omega }_{A\varphi i}^{2}}{{({\omega }_{{Ai}}^{2}+{\omega }_{{si}}^{2})}^{1/2}}+{\left(\displaystyle \frac{{\omega }_{A\varphi i}^{4}}{{\omega }_{{Ai}}^{2}+{\omega }_{{si}}^{2}}+4{\omega }_{{Ai}}^{2}\right)}^{1/2}\right].\end{eqnarray} \tag{ 34 }$$

For notational convenience, Equation (34) is rewritten as follows:

$$\begin{eqnarray}&&\omega ={\omega }_{{Ai}}h,\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&h=\displaystyle \frac{1}{2}\left[\displaystyle \frac{{q}_{i}^{2}}{{(1+{d}^{2})}^{1/2}}+{\left(4+\displaystyle \frac{{q}_{i}^{4}}{1+{d}^{2}}\right)}^{1/2}\right]\end{eqnarray} \tag{ 36 }$$

where ${q}_{i}={B}_{\varphi A}/{B}_{{zi}}$ and $d={v}_{{si}}/{v}_{{Ai}}$. This ω is used as an approximation to the resonance frequency, ${\omega }_{0}$, in Section 4. Lastly, we should note that, given this value for ${\omega }_{0}$, although the variation of the Alfvén speed across the inhomogeneity in the flux tube is quadratic (see Figure 3), since ${\omega }_{{Ai}}\lt {\omega }_{0}\lt {\omega }_{{Ae}}$, there will only be a single resonance point.

For an equilibrium with magnetic twist, such as the model used in this work, the total pressure perturbation is no longer a conserved quantity, so Equations (27a) and (27b) require modification. Sakurai et al. (1991a) derived new conserved quantities for the Alfvén and slow continua. Specifically, for the Alfvén continuum the conserved quantity is

$$\begin{eqnarray}&&{C}_{A}={g}_{B}{p}_{T}^{\prime }-2{f}_{B}{B}_{\varphi }{B}_{z}\displaystyle \frac{{\xi }_{r}}{{\mu }_{0}r}.\end{eqnarray} \tag{ 37 }$$

Using this conserved quantity, they derived jump conditions for ${\xi }_{r}$ and ${p}_{T}^{\prime }$, namely a prescription on how the radial component of the Lagrangian displacement and the total pressure perturbation can vary across the inhomogeneous layer connecting the internal with the external solutions. This prescription then implies that the following conditions must be satisfied:

$$\begin{eqnarray}&&{{\xi }_{{ri}}(r)| }_{r={r}_{i}}+[\kern-2pt[ {\xi }_{r}(r)]\kern-2pt] ={{\xi }_{{re}}(r)| }_{r={r}_{e}}\end{eqnarray} \tag{ 38 }$$

and

$$\begin{eqnarray}&&{{p}_{{Ti}}^{\prime }(r)| }_{r={r}_{i}}+[\kern-2pt[ {p}_{T}^{\prime }(r)]\kern-2pt] ={{p}_{{Te}}^{\prime }(r)| }_{r={r}_{e}},\end{eqnarray} \tag{ 39 }$$

where $[\kern-2pt[ {\xi }_{r}]\kern-2pt]$ and $[\kern-2pt[ {p}_{T}^{\prime }]\kern-2pt]$ are the jump conditions across the resonant layer in the inhomogeneous section of the flux tube, in the radial displacement and total pressure perturbation (Sakurai et al. 1991a). They are given by

$$\begin{eqnarray}&&[\kern-2pt[ {\xi }_{r}]\kern-2pt] =-\displaystyle \frac{{\imath }\pi }{| {{\rm{\Delta }}}_{A}| }\displaystyle \frac{{g}_{B}}{{\mu }_{0}{\rho }^{2}{v}_{A}^{2}}{C}_{A},\end{eqnarray} \tag{ 40 }$$

and

$$\begin{eqnarray}&&[\kern-2pt[ {p}_{T}^{\prime }]\kern-2pt] =-\displaystyle \frac{{\imath }\pi }{| {{\rm{\Delta }}}_{A}| }\displaystyle \frac{2{{TB}}_{z}}{{\mu }_{0}{\rho }^{2}{v}_{A}^{2}r}{C}_{A},\end{eqnarray} \tag{ 41 }$$

where $T={f}_{B}{B}_{\varphi }/{\mu }_{0}={k}_{z}{B}_{z}{B}_{\varphi }/{\mu }_{0}$ and

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{A}=\displaystyle \frac{d}{{dr}}({\omega }^{2}-{\omega }_{A}^{2}(r)).\end{eqnarray} \tag{ 42 }$$

Taking into account that $m=0$ and ${B}_{\varphi }\ne 0$ and Equation (11), the jump conditions, Equations (40) and (41), can be written as

$$\begin{eqnarray}[\kern-2pt[ {\xi }_{r}]\kern-2pt] & = & \displaystyle \frac{{\imath }\pi }{| {{\rm{\Delta }}}_{A}| }\displaystyle \frac{{k}_{z}{B}_{\varphi }}{{\mu }_{0}{\rho }^{2}{v}_{A}^{2}}{C}_{A}\\ & = & {\left.-\displaystyle \frac{{\imath }\pi }{| {{\rm{\Delta }}}_{A}| }\displaystyle \frac{{k}_{z}^{2}{B}_{\varphi }^{2}}{{\mu }_{0}{\rho }^{2}{v}_{A}^{2}}\right|}_{r={r}_{A}}\cdot {\left.\left\{{p}_{T}^{\prime }+\displaystyle \frac{2{B}_{z}^{2}}{{\mu }_{0}}\displaystyle \frac{{\xi }_{r}}{r}\right\}\right|}_{r={r}_{i}},\end{eqnarray} \tag{ 43 }$$

and

$$\begin{eqnarray}[\kern-2pt[ {p}_{T}^{\prime }]\kern-2pt] & = & -\displaystyle \frac{{\imath }\pi }{| {{\rm{\Delta }}}_{A}| }\displaystyle \frac{2{k}_{z}{B}_{\varphi }{B}_{z}^{2}}{{\mu }_{0}^{2}{\rho }^{2}{v}_{A}^{2}r}{C}_{A}\\ & = & {\left.\displaystyle \frac{2{\imath }\pi }{r| {{\rm{\Delta }}}_{A}| }{\left(\displaystyle \frac{{k}_{z}{B}_{\varphi }{B}_{z}}{{\mu }_{0}\rho {v}_{A}}\right)}^{2}\right|}_{r={r}_{A}}\cdot {\left.\left\{{p}_{T}^{\prime }+2\displaystyle \frac{{B}_{z}^{2}}{{\mu }_{0}}\displaystyle \frac{{\xi }_{r}}{r}\right\}\right|}_{r={r}_{i}}.\end{eqnarray} \tag{ 44 }$$

Given that ${B}_{z}={B}_{z}(r)$ and $\rho =\rho (r)$ in the inhomogeneous layer (see Equations (4), (5), and Figure 2), we have

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{A}={\omega }_{A}{(r)}^{2}{{\rm{\Delta }}}_{{AF}}={\omega }_{A}{(r)}^{2}\left[\displaystyle \frac{1}{\rho }\displaystyle \frac{d\rho }{{dr}}-2\displaystyle \frac{1}{{B}_{z}}\displaystyle \frac{{{dB}}_{z}}{{dr}}\right].\end{eqnarray} \tag{ 45 }$$

Obviously, when ${B}_{z}$ is constant across the inhomogeneous layer,

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{A}={\omega }_{A}^{2}(r)\displaystyle \frac{1}{\rho }\displaystyle \frac{d\rho }{{dr}}.\end{eqnarray} \tag{ 46 }$$

Substituting Equations (43)–(45) into Equations (38) and (39), we obtain the dispersion relation for axisymmetric MHD waves that undergo resonant absorption in the Alfvén continuum of frequencies due to the twist in the magnetic field:

$$\begin{eqnarray}&&{{ \mathcal D }}_{{AR}}(\omega ,{k}_{z})+{\imath }{{ \mathcal D }}_{{AI}}(\omega ,{k}_{z})=0,\end{eqnarray} \tag{ 47 }$$

where

$$\begin{eqnarray}{{ \mathcal D }}_{{AR}} & = & 2\displaystyle \frac{{D}_{i}}{{k}_{{ri}}^{2}}\displaystyle \frac{M(a,b-1;{s}_{i})}{M(a,b;{s}_{i})}-2{\rho }_{i}{n}_{i}({n}_{i}+{k}_{z})\displaystyle \frac{{v}_{A\varphi i}^{2}}{{k}_{{ri}}^{2}}\\ & & +\displaystyle \frac{{r}_{i}{D}_{e}}{{k}_{{re}}}\displaystyle \frac{{K}_{\nu -1}({k}_{{re}}{r}_{e})}{{K}_{\nu }({k}_{{re}}{r}_{e})}-\displaystyle \frac{{r}_{i}}{{r}_{e}{k}_{{re}}^{2}}\{(1-\nu ){D}_{e}\\ & & -2{\rho }_{e}{n}_{e}^{2}{v}_{A\varphi e}^{2}\},\end{eqnarray} \tag{ 48 }$$

and

$$\begin{eqnarray}&&{{ \mathcal D }}_{{AI}}={\left.\displaystyle \frac{\pi }{\rho | {{\rm{\Delta }}}_{{AF}}| }\displaystyle \frac{{v}_{A\varphi }^{2}}{{v}_{A}^{4}}\right|}_{r={r}_{A}}\\ &&\times \left[\displaystyle \frac{2}{{k}_{{ri}}^{2}}\left({D}_{i}\displaystyle \frac{M(a,b-1;{s}_{i})}{M(a,b;{s}_{i})}-{n}_{i}({n}_{i}+{k}_{z}){\rho }_{i}{v}_{A\varphi i}^{2}\right)+2{\rho }_{i}{v}_{{Ai}}^{2}\right]\\ &&\times \left[2{\left.\displaystyle \frac{{B}_{z}^{2}}{{\mu }_{0}r}\right|}_{{r}_{A}}+\displaystyle \frac{1}{{r}_{e}{k}_{{re}}^{2}}\{(1-\nu ){D}_{e}-2{\rho }_{e}{n}_{e}^{2}{v}_{A\varphi e}^{2}\}\right.\\ &&-\left.\displaystyle \frac{{D}_{e}}{{k}_{{re}}}\displaystyle \frac{{K}_{\nu -1}({k}_{{re}}{r}_{e})}{{K}_{\nu }({k}_{{re}}{r}_{e})}\right].\end{eqnarray} \tag{ 49 }$$

In these equations the following definitions were used:

$$\begin{eqnarray}&&{v}_{A\varphi A}^{2}=\displaystyle \frac{{B}_{\varphi A}^{2}}{{\mu }_{0}{\rho }_{A}},{v}_{{AA}}^{2}({r}_{A})=\displaystyle \frac{{B}_{{zA}}^{2}}{{\mu }_{0}{\rho }_{A}},\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&| {{\rm{\Delta }}}_{{AF}}({r}_{A})| =\displaystyle \frac{1}{{\ell }}\left|\displaystyle \frac{{\rho }_{e}-{\rho }_{i}}{{\rho }_{A}}-2\displaystyle \frac{{B}_{{ze}}-{B}_{{zi}}}{{B}_{z}({r}_{A})}\right|,\end{eqnarray} \tag{ 51 }$$

where ${\rho }_{A}=\rho ({r}_{A})$, ${v}_{{AA}}={v}_{A}({r}_{A})$, and ${B}_{{zA}}={B}_{z}({r}_{A})$. To find the radius at the resonance point, namely the radius where ${v}_{A}(r)={v}_{0}(={\omega }_{0}/{k}_{z})$,⁵ we can express r_A as a convex combination of the radius r_i and the width of the inhomogeneous layer ℓ since r_A must be within the interval $({r}_{i},{r}_{e})$. Therefore we can write ${r}_{A}={r}_{i}+w{\ell }$, where $w\in (0,1)$. Now we have transformed the problem of solving for r_A into a problem where we have to solve for w, the convex combination parameter. Given this formulation for r_A and Equations (4) and (5), we can write ${B}_{{zA}}={B}_{z}({r}_{A})={B}_{{zi}}+w({B}_{{ze}}-{B}_{{zi}})$ and ${\rho }_{A}=\rho ({r}_{A})={\rho }_{i}+w({\rho }_{e}-{\rho }_{i})$. Equipped with these definitions, the equation that we need to solve to find w becomes

$$\begin{eqnarray}&&{v}_{{AA}}^{2}=\displaystyle \frac{{({B}_{{zi}}+w({B}_{{ze}}-{B}_{{zi}}))}^{2}}{{\mu }_{0}({\rho }_{e}+w({\rho }_{e}-{\rho }_{i}))}={v}_{0}^{2}\end{eqnarray} \tag{ 52 }$$

Using the definitions $\chi ={\rho }_{e}/{\rho }_{i}$ and $\zeta ={B}_{{ze}}/{B}_{{zi}}$, Equation (52) simplifies to

$$\begin{eqnarray}&&{v}_{{Ai}}^{2}\displaystyle \frac{{(1+w(\zeta -1))}^{2}}{1+w(\chi -1)}={v}_{0}^{2}.\end{eqnarray} \tag{ 53 }$$

This equation is solved for w in the next section.

Taking the long wavelength limit, $\epsilon \ll 1$, of Equations (48) and (49) and using Equation (29) and Equation (31), then Equations (48) and (49) reduce to

$$\begin{eqnarray}&&{{ \mathcal D }}_{{AR}}=\displaystyle \frac{{\omega }^{2}-{\omega }_{{Ai}}^{2}}{{k}_{{ri}}^{2}}-\displaystyle \frac{{n}_{i}^{2}}{{k}_{{ri}}^{2}}{v}_{A\varphi i}^{2}+\chi \displaystyle \frac{{n}_{e}^{2}}{{k}_{{re}}^{2}}{v}_{A\varphi e}^{2},\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&{{ \mathcal D }}_{{AI}}=\displaystyle \frac{\pi }{{r}_{e}}\left[{\left.\displaystyle \frac{\pi }{{{\rm{\Delta }}}_{{AF}}}\displaystyle \frac{{v}_{A\varphi }^{2}}{{v}_{A}^{2}}\right|}_{r={r}_{A}}\right]\left[1-\displaystyle \frac{{\rho }_{e}}{{\rho }_{A}}\displaystyle \frac{{n}_{e}^{2}}{{k}_{{re}}^{2}}\displaystyle \frac{{v}_{A\varphi e}^{2}}{{v}_{{AA}}^{2}}\right]\\ &&\quad \times \left[\displaystyle \frac{{\omega }^{2}-{\omega }_{{Ai}}^{2}}{{k}_{{ri}}^{2}}-\displaystyle \frac{{n}_{i}^{2}}{{k}_{{ri}}^{2}}{v}_{A\varphi i}^{2}+{v}_{{Ai}}^{2}\right].\end{eqnarray} \tag{ 55 }$$

These equations can be solved if we allow a complex frequency $\omega ={\omega }_{r}+{\imath }{\gamma }_{A}$, and when ${\gamma }_{A}\ll {\omega }_{r}$ we can obtain the damping rate; then ${\gamma }_{A}$ in the Alfvén continuum frequencies (Goossens et al. 1992) to second order is given by

$$\begin{eqnarray}&&{\gamma }_{A}=-{{ \mathcal D }}_{{AI}}({\omega }_{0}){\left({\left.\displaystyle \frac{\partial {{ \mathcal D }}_{{AR}}}{\partial \omega }\right|}_{\omega ={\omega }_{0}}\right)}^{-1}.\end{eqnarray} \tag{ 56 }$$

This equation results in an expression that is difficult to interpret. For this reason, given that we seek an expression for the damping rate in the long wavelength limit, we expand it in a series about $\varepsilon =0$ where $\varepsilon ={r}_{e}{k}_{z}$. This expansion results in

$$\begin{eqnarray}&&{\gamma }_{A}={\omega }_{0}\displaystyle \frac{\pi }{Z}\displaystyle \frac{{\ell }}{{r}_{e}}\displaystyle \frac{{\rho }_{A}}{{\rho }_{i}}\displaystyle \frac{{B}_{\varphi A}^{2}}{{B}_{{zA}}^{2}}\left(1+\displaystyle \frac{{B}_{\varphi A}^{2}}{{B}_{{zA}}^{2}}\right)\left(1+\displaystyle \frac{{B}_{\varphi A}^{2}}{{B}_{{zi}}^{2}}\right)+{ \mathcal O }({\varepsilon }^{2}),\end{eqnarray} \tag{ 57 }$$

where

$$\begin{eqnarray}&&Z=\left|(\chi -1)-2\displaystyle \frac{{\rho }_{A}}{{\rho }_{i}}\displaystyle \frac{{B}_{{zi}}}{{B}_{{zA}}}(\zeta -1)\right|.\end{eqnarray} \tag{ 58 }$$

Now, in this investigation we assume a weak magnetic twist ($q={B}_{\varphi A}/{B}_{{zA}}\ll 1$), and therefore

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{\varphi A}^{2}}{{B}_{{zA}}^{2}}\left(1+\displaystyle \frac{{B}_{\varphi A}^{2}}{{B}_{{zA}}^{2}}\right)\left(1+\displaystyle \frac{{B}_{\varphi A}^{2}}{{B}_{{zi}}^{2}}\right)=\displaystyle \frac{{B}_{\varphi A}^{2}}{{B}_{{zA}}^{2}}+{ \mathcal O }({q}^{4}),\end{eqnarray} \tag{ 59 }$$

and Equation (57) can be simplified to

$$\begin{eqnarray}&&{\gamma }_{A}={\omega }_{0}\displaystyle \frac{\pi }{Z}\displaystyle \frac{{\ell }}{{r}_{e}}\displaystyle \frac{{\rho }_{A}}{{\rho }_{i}}\displaystyle \frac{{B}_{\varphi A}^{2}}{{B}_{{zA}}^{2}}.\end{eqnarray} \tag{ 60 }$$

Here ${\omega }_{0}$ is approximated by Equation (34), i.e., ${\omega }_{0}\approx {\omega }_{{Ai}}h$, and the radius at the resonance point, r_A, is obtained analogously to Equation (53) by solving

$$\begin{eqnarray}&&\displaystyle \frac{{(1+w(\zeta -1))}^{2}}{1+w(\chi -1)}={h}^{2}.\end{eqnarray} \tag{ 61 }$$

There are two cases to be considered. First, when $\zeta =1$, that is ${B}_{{zi}}={B}_{{ze}}$, and assuming $\chi \in (0,1/{h}^{2})$, the solution for w is

$$\begin{eqnarray}&&w=\displaystyle \frac{{h}^{2}-1}{{h}^{2}(1-\chi )}.\end{eqnarray} \tag{ 62 }$$

When $1/{h}^{2}\lt \chi \lt 1$ in this case, the resonant point is outside of the continuum, and there is no resonant absorption. In the limit $\chi \to 1$, the external and internal Alfvén speeds become equal, and there are no propagating waves either. The second case is for values of $\zeta \in (0,1)$ and $\chi \in (0,{\zeta }^{2}/{h}^{2})$, for which the admissible solution is

$$\begin{eqnarray}w\;= & \displaystyle \frac{2(1-\zeta )+h{(4(\zeta -1)(\zeta -\chi )+{h}^{2}{(\chi -1)}^{2})}^{1/2}+{h}^{2}(\chi -1)}{2{(1-\zeta )}^{2}}.\end{eqnarray} \tag{ 63 }$$

When $\chi \gt {\zeta }^{2}/{h}^{2}$, similarly to the first case there is no resonant absorption since the resonance frequency (${\omega }_{0}$) is outside the continuum. For ${\zeta }^{2}/{h}^{2}\lt \chi \lt {\zeta }^{2}$ there exist undamped propagating waves, but when $\chi \gt {\zeta }^{2}$ the external Alfvén speed is smaller than the internal, and no waves propagate. Lastly note that in this investigation we assume that ${B}_{{zi}}\geqslant {B}_{{ze}}$, so ${\rho }_{i}={\rho }_{e}$ has no admissible solution for w when ${B}_{{zi}}={B}_{{ze}}$.

Now, when ${B}_{z}$ is assumed to be constant, i.e., ${B}_{{zi}}={B}_{{ze}}={B}_{z}$, using Equation (62), ${\rho }_{A}/{\rho }_{i}$, ${B}_{{zA}}/{B}_{{zi}}$, and Z simplify to

$$\begin{eqnarray}&&\displaystyle \frac{{\rho }_{A}}{{\rho }_{i}}=\displaystyle \frac{1}{{h}^{2}},\displaystyle \frac{{B}_{{zA}}}{{B}_{{zi}}}=1,Z=| 1-\chi | ,\end{eqnarray} \tag{ 64 }$$

resulting in ${\gamma }_{A}$ (Equation (60)):

$$\begin{eqnarray}&&{\gamma }_{A}=\displaystyle \frac{{\omega }_{0}}{{h}^{2}}\displaystyle \frac{\pi }{| 1-\chi | }\displaystyle \frac{{\ell }}{{r}_{e}}\displaystyle \frac{{B}_{\varphi A}^{2}}{{B}_{{zA}}^{2}}=\displaystyle \frac{{\omega }_{{Ai}}}{h}\displaystyle \frac{\pi }{| 1-\chi | }\displaystyle \frac{{\ell }}{{r}_{e}}\displaystyle \frac{{B}_{\varphi A}^{2}}{{B}_{{zA}}^{2}}.\end{eqnarray} \tag{ 65 }$$

To obtain the damping time normalized by the period of the wave, we use a typical wavelength ${k}_{z}=\pi /L$, where L is the characteristic length of the tube, and the associated period is $\tau =2\;L/{{hv}}_{{Ai}}$ (see Equation (35)). The damping time ($1/{\gamma }_{A}$) for modes in the continuum as a multiple of the wave period is

$$\begin{eqnarray}&&{\tau }_{d}=\displaystyle \frac{Z}{2{\pi }^{2}}\displaystyle \frac{{r}_{e}}{{\ell }}\displaystyle \frac{{\rho }_{i}}{{\rho }_{A}}\displaystyle \frac{{B}_{{zA}}^{2}}{{B}_{\varphi A}^{2}}\tau .\end{eqnarray} \tag{ 66 }$$

A contour map of this equation for $\zeta \in (0,1)$ and $\chi \in (0,{\zeta }^{2}/{h}^{2})$ can be seen in Figure 4. When ${B}_{{zi}}={B}_{{ze}}$, the damping time becomes

$$\begin{eqnarray}&&{\tau }_{d}={h}^{2}\displaystyle \frac{| 1-\chi | }{2{\pi }^{2}}\displaystyle \frac{{r}_{e}}{{\ell }}\displaystyle \frac{{B}_{{zA}}^{2}}{{B}_{\varphi A}^{2}}\tau .\end{eqnarray} \tag{ 67 }$$

The long wavelength limit approximation of the damping rate ${\gamma }_{A}$ in Equation (57) is accurate to $\approx {10}^{-6}$ at ${k}_{z}{r}_{e}=1$ when compared with the numerical solution of the dispersion relation in Equation (47). This accuracy is better than 10⁻⁶ for ${k}_{z}{r}_{e}\lt 0.1$ and is calculated using the maximum of the root-mean-square error (RMSE):

$$\begin{eqnarray}&&{\rm{RMS}}\;{\rm{Error}}={\left(\displaystyle \frac{1}{N-1}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{\gamma }_{A}-{\hat{\gamma }}_{A}}{{\gamma }_{A}}\right)}^{2}\right)}^{1/2}.\end{eqnarray} \tag{ 68 }$$

In this equation, ${\gamma }_{A}$ is the numerically calculated damping rate, ${\hat{\gamma }}_{A}$ is the theoretical approximation in Equation (57), and N is the number of samples. For this error estimate we used 10⁴ samples in the parameter space ($\chi ,\zeta ,{\ell }/{r}_{e},{B}_{\varphi }/{B}_{z}$), uniformly distributed.⁶

Zoom In Zoom Out Reset image size

Works investigating resonant absorption in the context of solar atmospheric conditions tend to consider solely a radial nonuniformity in either the magnetic field or density. However, accounting for radial variation in both the magnetic field and density can lead to significant variation in the estimated damping times. The ratio of Equation (66) over Equation (67) is

$$\begin{eqnarray}&&\displaystyle \frac{{\tau }_{d}(\chi ,\zeta )}{{\tau }_{d}(\chi ,1)}=\displaystyle \frac{Z}{| 1-\chi | }\displaystyle \frac{{\rho }_{i}}{{\rho }_{A}}\displaystyle \frac{1}{{h}^{2}},\end{eqnarray} \tag{ 69 }$$

and in Figure 5 a contour map is shown for $\zeta \in (0,1)$ and $\chi \in (0,{\zeta }^{2}/{h}^{2})$.

Zoom In Zoom Out Reset image size

It can be seen from Figure 4 that the behavior of the damping rate with respect to changes in the density contrast is in some regions exactly the opposite of that for the kink mode (Goossens & Poedts 1992, for example). Namely, in a roughly triangular region in Figure 4, the damping rate is proportional to $\sim 1/\chi$, in contrast to the kink mode, where the damping rate is proportional to $\sim \chi$. Similar behavior has been been shown to exist in the leaky regime for sausage modes (Vasheghani Farahani et al. 2014). The factor in Equation (66) that determines this behavior is $Z{\rho }_{i}/{\rho }_{A}$. We approximate the local minimum in the χ direction by evaluating the partial derivative of Z with respect to χ:

$$\begin{eqnarray}&&\displaystyle \frac{\partial Z}{\partial \chi }=\displaystyle \frac{-2\zeta +{h}^{2}(\chi -1)+2}{h\sqrt{4(\zeta -1)(\zeta -\chi )+{h}^{2}{(\chi -1)}^{2}}},\end{eqnarray} \tag{ 70 }$$

which is subsequently equated to zero. From this we obtain a relation $\chi =(2/{h}^{2})\zeta +b$, and b is identified by noting that at $\zeta =1$ the maximum value for χ is $1/{h}^{2}$, so the approximation is

$$\begin{eqnarray}&&\chi =\displaystyle \frac{2}{{h}^{2}}\zeta -\displaystyle \frac{1}{{h}^{2}}.\end{eqnarray} \tag{ 71 }$$

Because the remaining terms in Equation (66) do not vary with χ and ζ (note the ratio ${B}_{\varphi }$/${B}_{z}$ is held fixed), this approximation holds for all valid parameters. This approximation allows us to estimate in which regime a specific parameter combination exists. Namely, for parameter combinations that are below the line described by Equation (71), for increasing density contrast ($\chi \downarrow$), damping will be slower (${\tau }_{d}\uparrow$). For parameter combinations that result in points above this line, increasing the density contrast ($\chi \downarrow$) results in decreasing damping time (${\tau }_{d}\downarrow$); thus waves will decay faster. This is illustrated in Figure 4 as a yellow line (Equation (71)), and the exact inflection points are marked with a green line.

Given the form of Equation (66), and especially that of Equation (67), a comparison with previous results for the kink mode is in order, particularly the expression for the damping rate obtained by Goossens et al. (1992) and later by Ruderman & Roberts (2002). In Ruderman & Roberts (2002), and Equation (73) in that work, using the notation in this work, the expression reads as follows:

$$\begin{eqnarray}&&{\tau }_{d}=\displaystyle \frac{2}{\pi }\displaystyle \frac{{r}_{e}}{{\ell }}\displaystyle \frac{1-\chi }{1+\chi }\tau .\end{eqnarray} \tag{ 72 }$$

The relative magnitude of the damping time shown in Equations (72) and (67) is

$$\begin{eqnarray}&&\displaystyle \frac{{\tau }_{d,\mathrm{Axisymmetric}}}{{\tau }_{d,\mathrm{Kink}}}=\displaystyle \frac{{h}^{2}}{4\pi }\displaystyle \frac{{(1-\chi )}^{2}}{1+\chi }\displaystyle \frac{{B}_{{zA}}^{2}}{{B}_{\varphi A}^{2}}.\end{eqnarray} \tag{ 73 }$$

It is evident that there exists a region in the parameter space of $(\chi ,{B}_{\varphi }/{B}_{z})$ for which ${\tau }_{d,\mathrm{Axisymmetric}}$ is smaller than ${\tau }_{d,\mathrm{Kink}}$, but this comparison is given here just as a reference. Caution should be exercised in its interpretation since the damping, ${\tau }_{d,\mathrm{Kink}}$, in Ruderman & Roberts (2002) was calculated for the kink mode without a magnetic twist. It is possible that a magnetic twist amplifies dissipation in the kink mode, and therefore dissipation for the kink mode may be larger than that of axisymmetric modes.

4.1. Numerical Solution of Dispersion Relation in the Alfvén Continuum

We have solved Equation (47) numerically using ${\omega }_{0}$ obtained in Equation (34) as an initial point in the solver. Additionally, by means of investigating whether another solution exists, we solved the dispersion relation again with a random ${\omega }_{0}$ in the range $({v}_{{Ai}},{v}_{{Ae}})$. The solutions and their associated damping rates can be seen in Figure 6. It is interesting that there exists another solution in the long wavelength limit that we could not obtain from our analysis in Section 4. However, given that for this solution ${\tau }_{d}\ll \tau$, it is unlikely that this mode will be observed.

Zoom In Zoom Out Reset image size

Now it has been shown that the singularity about the resonance point at r_A is logarithmic for ${\xi }_{r}$ ($\mathrm{ln}(| r-{r}_{A}| )$) and $1/(r-{r}_{A})$ for ${\xi }_{\perp }$, so the dynamics will be governed by ${\xi }_{\perp }$ since ${\xi }_{r}/{\xi }_{\perp }\to 0$ as $r\to {r}_{A}$, and therefore ${\xi }_{\perp }\gg {\xi }_{{re}}$ in the neighborhood of the resonant point (Poedts et al. 1989; Sakurai et al. 1991a). Also, the ${\xi }_{r}$ component provides its energy to the resonant layer (Goossens et al. 2011), so the characteristic expansion and contraction of axisymmetric modes will be reduced. These facts, along with the proximity of the solution corresponding to the long wavelength limit approximation in Section 4 to the internal Afvén speed, suggest that these waves would appear in observations to have properties similar to Alfvén waves. Given that pure Alfvén waves require ${\rm{\nabla }}\cdot {\boldsymbol{\xi }}$ to be identically zero and a driving mechanism that is solely torsional, we argue that observed Alfvén waves are much more likely to be axisymmetric waves because they do not have these strict requirements. In Figure 6, panels (a) through (d) show solutions for different values of χ, while in panels (e) through (h) solutions are shown when q is allowed to vary. The damping time for the solution for which we have an analytical approximation (see panels (c)) and (d)) increases (${\tau }_{d}\uparrow$) for increasing density contrast ($\chi \downarrow$), while the other solution exhibits the opposite behavior (see panels (a) and (b)), namely ${\tau }_{d}\downarrow$ for $\chi \downarrow$. However, the damping time for both solutions decreases (${\tau }_{d}\downarrow$) for increasing magnetic twist ($q\uparrow$). The bottom panel of Figure 6 shows a different view of the damping times as a function of the magnetic twist (q), shown in panels (f) and (h) at ${k}_{z}{r}_{e}=\{0.1,0.15,0.2,0.25,0.3\}$. From this view it can be seen that the solutions in (e) are much more sensitive to variations in the magnetic twist when compared with the solutions in panel (g). This sensitivity, in combination with the fact that for an extremely small twist the sausage cutoff is reintroduced (Giagkiozis et al. 2015), means that this mode will be observable for a very small interval of magnetic twist. The mode shown in panels (c) and (d) does not present this difficulty, and therefore we expect that observation of this mode is more likely. In both cases, the solution corresponding to the analytic approximation remains very close to the internal Alfvén speed, which is equal to 1 in Figure 6. Since ${\omega }_{0}$ from Equation (35) depends on q, B_zA, and the internal sound speed, these modes will appear to have a strong Alfvén character for virtually all valid parameter combinations. Lastly, ${k}_{r}$ can be likened to the wavenumber in the radial direction, and, since in the long wavelength limit ${k}_{r}$ is proportional to ${k}_{z}$, as ${k}_{z}$ increases, the wavelength in the radial direction decreases and couples with the thin inhomogeneous layer more closely. Therefore more energy per wavelength is absorbed, and thus the damping time is reduced (see Figure 6).

Reports of observations of axisymmetric modes (sausage modes) are increasing in frequency in the recent literature. For example, QPPs in solar flares are believed to be associated with the kink and sausage mode (see, for example, Nakariakov & Zimovets 2011; Van Doorsselaere et al. 2011; De Moortel & Nakariakov 2012; Nakariakov 2012; Kolotkov et al. 2015). Even more interestingly, some of these pulsations appear to have periods in the interval $(15,100)$ s, which could be consistent with the results in the present investigation if the length-scale of these pulsations is on the same order as the length-scale of coronal flux tubes, $\approx 100\;\mathrm{Mm}$. Furthermore, the results by Morton et al. (2012) suggest that axisymmetric modes are ubiquitous and that they appear to coexist with kink modes. This coexistence further supports the argument by Arregui et al. (2015), Arregui & Soler (2015), and Arregui (2015) that Bayesian analysis is an essential tool for the identification of the likely wave modes present in observations, as well as a more systematic method of appropriate model selection. The uncertainty in determining the parameters for the kink mode led Verwichte et al. (2013) to perform a statistical analysis as a way to narrow the range of their values. This departure from certainty and convergence toward probabilistic inference models for solar observations is, in our view, long overdue.

However, despite this increase in interest in axisymmetric modes, the relation that approximates their expected damping rate (see Equation (66)) requires knowledge of four parameters: the density, the magnetic field contrast, the relative magnetic twist, and the ratio of the thickness of the inhomogeneous layer versus the tube radius, or $(\chi ,\zeta ,q={B}_{\varphi A}/{B}_{{zA}},{\ell }/{r}_{e})$. In contrast to the large body of observational evidence for the kink mode, observations of sausage waves are relatively scarce. This makes impossible an analysis similar to that of Verwichte et al. (2013) for these modes. Therefore, we adopt a different approach, a probabilistic approach that is related to the use of Bayesian inference suggested by Arregui et al. (2015).

As a first step toward improving this situation, we provide a way to estimate the probability that an observed sausage wave has a damping rate within a specified range, given that one or more of the four parameters in Equation (57) are known. The assumptions required for the validity of this estimate are the following:

• 1.  
  The four parameters in Equation (57) are independent; that is, no parameter is a function of the others.
• 2.  
  The likelihood of any combination in the parameter space is the same. That is to say that there exists no preferred combination of parameters.

These assumptions are difficult to prove, especially given that there exist no statistical analyses of the properties of sausage waves or reliable estimates of all four parameters. Since we do not know if there is, in fact, a set of preferred parameters, these assumptions are required for an unbiased estimate. Acknowledging these uncertainties, we make a first attempt in identifying the probability predicted by our model that a wave with the characteristics described in this investigation is resonantly damped in the long wavelength limit with a damping rate given by Equation (66) for a given parameter combination.

The aforementioned probability can be estimated as follows. First, we identify the parameters for which reasonably good estimates are available. These parameters we refer to as free parameters, denoted by f. The remaining parameters we refer to as integration parameters and are denoted by i. Subsequently, a domain is defined for the integration parameters. Then the probability of the damping rate being within the open interval (a, b) is given by

$$\begin{eqnarray}&&P(a,b;{f}_{1},...,{f}_{n})\\ &&\quad =\displaystyle \frac{{\displaystyle \int }_{C}{{di}}_{1}...{{di}}_{4-n}w({i}_{1},...,{i}_{4-n}){I}_{{\bar{\tau }}_{d}\gt a,{\bar{\tau }}_{d}\lt b}[\cdot ]}{{\displaystyle \int }_{C}{{di}}_{1}...{{di}}_{4-n}w({i}_{1},...,{i}_{4-n}){I}_{1}[\cdot ]},\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}&&{I}_{{\bar{\tau }}_{d}\gt a,{\bar{\tau }}_{d}\lt b}[\cdot ]={I}_{{\bar{\tau }}_{d}\gt a,{\bar{\tau }}_{d}\lt b}[{\bar{\tau }}_{d}({i}_{1},...,{i}_{4-n};{f}_{1},...,{f}_{n})],\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}&&{I}_{1}[\cdot ]={I}_{1}[{\bar{\tau }}_{d}({i}_{1},...,{i}_{4-n};{f}_{1},...,{f}_{n})],\end{eqnarray} \tag{ 76 }$$

where C is the domain of integration, defined as the set of all elements in the integration parameter space that are valid according to the analysis in this work, and ${\bar{\tau }}_{d}={\tau }_{d}/\tau$. The function $I(\cdot )$ is an indicator function and $n=\{1,2,3\}$; that is, an estimate for at least one parameter is necessary. When the indicator function is subscripted with 1, it simply returns 1 when the parameter combination is valid. That is, the integral in the denominator of Equation (74) simply returns the area where ${\bar{\tau }}_{d}$ is defined. The indicator function in the numerator is defined as follows:

$$\begin{eqnarray}&&{I}_{{\bar{\tau }}_{d}\gt a,{\bar{\tau }}_{d}\lt b}[{\bar{\tau }}_{d}({i}_{1},...,{i}_{4-n};{f}_{1},...,{f}_{n})]\\ &&\quad =\;\left\{\begin{array}{l}1\;a\lt {\bar{\tau }}_{d}\lt b,\\ 0\ {\rm{otherwise}},\end{array}\right.\end{eqnarray} \tag{ 77 }$$

so this function returns 1 when the normalized damping rate (${\bar{\tau }}_{d}$) is within the open interval (a, b). Therefore the numerator of Equation (74) returns the area in C for which the normalized damping rate is within the interval (a, b). The function $w({i}_{1},{i}_{2})$ is a weighting function that is nonnegative, and its integral over C is equal to 1. Because we have assumed that every combination in the integration parameter space $({i}_{1},{i}_{2})$ is equally likely, this function is simply a constant and simplifies out from the integrals. The effect of this function is similar to the prior information in Bayesian inference. Therefore if relevant information of a specific preference in parameter space is present in the solar atmosphere, this can be taken into account by appropriately modifying $w(\cdot )$.

For the contour maps in Figures 7 and 8, it is assumed that the free parameters in Equation (74) are $(\zeta ,q)$ and $(\zeta ,\chi )$, respectively. For this case, Equation (74) becomes

$$\begin{eqnarray}&&P(a,b;{f}_{1},{f}_{2})=\displaystyle \frac{{\int }_{C}{{di}}_{1}{{di}}_{2}w({i}_{1},{i}_{2}){I}_{{\bar{\tau }}_{d}\gt a,{\bar{\tau }}_{d}\lt b}[{\bar{\tau }}_{d}({i}_{1},{i}_{2};{f}_{1},{f}_{2})]}{{\int }_{C}{{di}}_{1}{{di}}_{2}w({i}_{1},{i}_{2}){I}_{1}[{\bar{\tau }}_{d}({i}_{1},{i}_{2};{f}_{1},{f}_{2})]}.\end{eqnarray} \tag{ 78 }$$

Our rationale for selecting the limits for the integration parameters in Figure 7 is based on the values for the parameters $(\chi ,{\ell }/{r}_{e})$ and the normalized damping rate reported in Aschwanden et al. (2003) for the kink mode. These authors used 11 cases of observed damping kink oscillations, and their estimates for these parameters are as follows: $\chi \approx 0.1$ and ${\ell }/{r}_{e}\in (0.1,0.5)$, and the observed normalized damping rates were in the interval $(1,3)$. The q parameter has been selected in an interval that ensures that the magnetic twist is small. As can be seen in both Figures 7 and 8, the probability for the resonantly absorbed axisymmetric modes for a wide range in parameters is significantly high. Also, as seen in Figure 7, normalized damping times in the interval $(1,3)$ are possible even for extremely small magnetic twist ($\approx 0.02$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the case where more information is available, Monte Carlo simulation can be used to estimate the probability density function (PDF) of the normalized damping time. We illustrate this with two examples. First we use the estimates from Morton et al. (2012). In that work the magnetic field is assumed to be constant inside and outside the flux tube, but this assumption is unlikely to be identically satisfied, so we allow a small variation in ζ in the interval $\zeta \in (0.95,1)$. The density contrast is taken to be in the interval $({10}^{-2},{10}^{-1})$. Since Morton et al. (2012) do not provide an estimate for the width of the inhomogenous layer, we allow it to vary uniformly in $(0.1,0.5)$, an interval that is in line with estimates in Goossens et al. (2002) and Aschwanden et al. (2003). In both examples, we assume the magnetic twist is within the interval $(0.1,0.2)$. Using these intervals and assuming a uniform distribution, we sample Equation (66) 10⁶ times. The estimated PDF for this set of parameters is the blue curve in Figure 9. The blue vertical line is the expectation value, which is equal to $E[{\tau }_{d}/\tau ]=7.49$.

Zoom In Zoom Out Reset image size

In the second example, we use parameter estimates from Van Doorsselaere et al. (2011). Assuming a H plasma, $\rho ={{Nm}}_{p}$, $p={{Nk}}_{B}T$ where N is the number density, m_p is the proton mass, and k_B is the Boltzmann constant. With the plasma-β equal to $\beta =2{\mu }_{0}p/{B}^{2}$ and the assumption that ${\beta }_{e}\ll {\beta }_{i}$ (Van Doorsselaere et al. 2011), we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{\beta }_{i}}{{\beta }_{e}}=\displaystyle \frac{{\zeta }^{2}}{\chi }\displaystyle \frac{{T}_{i}}{{T}_{e}}\gg 1.\end{eqnarray} \tag{ 79 }$$

Assuming a lower limit for ${\beta }_{i}/{\beta }_{e}\geqslant 100$ and a hot flux tube, ${T}_{i}/{T}_{e}=10$, we can restrict ζ and χ to

$$\begin{eqnarray}&&\displaystyle \frac{1}{200}\leqslant \chi \leqslant \displaystyle \frac{{\beta }_{e}}{{\beta }_{i}}\displaystyle \frac{{T}_{i}}{{T}_{e}},\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&{\left(\chi \displaystyle \frac{{T}_{e}}{{T}_{i}}\displaystyle \frac{{\beta }_{i}}{{\beta }_{e}}\right)}^{1/2}\leqslant \zeta \leqslant 1.\end{eqnarray} \tag{ 81 }$$

The lower limit for χ is considered as a minimum contrast in Van Doorsselaere et al. (2011) to avoid the sausage cutoff. However, in the presence of a very weak magnetic twist, this cutoff is removed, so we do not need to assume extreme values for the density contrast. The upper limit for χ and lower limit for ζ are taken so that Equation (79) is satisfied. The resulting PDF can be seen in Figure 9. It is interesting that the expected value for the damping time in this case is 5.58, which is very close to the observed damping (${\tau }_{d}/\tau =6$) of a mode that is believed to be a fast sausage mode (Kolotkov et al. 2015). It is apparent from Figure 9 that the PDFs cannot be approximated well using a normal distribution, so their use for obtaining estimates of the damping time from results like Equation (66) in this work and similar equations (e.g., Goossens et al. 1992; Ruderman & Roberts 2002) can be misleading. In contrast, Monte Carlo simulation and nonparametric density estimation can be quite useful tools for exploring this type of problem.

Theoretically, it has been known for some time that, in the presence of a weak magnetic twist, axisymmetric modes will be resonantly damped (see, for example, Goossens et al. 1992). In this work, we have calculated for the first time a dispersion relation for resonantly damped axisymmetric modes in the spectrum of the Alfvén continuum and derived an approximation of the damping time in the long wavelength limit. We have shown that the damping time can be comparable to that observed for the kink mode in the case where there is no magnetic twist. Furthermore, we solved the resulting equation (see Equations (47) and (57)) analytically, and (1) we confirmed the validity of our approximation, and (2) we found an additional solution that decays much faster in comparison. The resulting approximation in the long wavelength limit shows that the damping time is proportional to the magnetic twist and inversely proportional to the density contrast. It is interesting to note that Vasheghani Farahani et al. (2014), who investigated the damping of fast sausage modes in the leaky regime, found a similar relation between the damping time and the density contrast. However, in that work a very large density contrast is required to allow observation of the sausage mode. This is not the case for one of the results of this work, which for even modest density contrasts (see Figures 7 and 8) the damping time is within one to three periods of the wave.

Of the two solutions that we have uncovered, only the one with a phase velocity close to the internal Alfvén speed has, for some parameter combinations, damping times that would allow observation. The other solution is found to be damped on timescales $\approx {10}^{-2}\mbox{--}{10}^{-1}$ of the wave period, as seen in Figure 6, which means that its observation would be extremely challenging. However, the predicted damping times for the solution whose phase speed is close to the Alfvén speed are large enough to allow observation. Also, the fact that its phase speed is so close to the internal Alfvén speed, along with the dominance of the ${\xi }_{\perp }$ component in the wave dynamics, means that the character of this wave will be predominantly Alfvénic (Goossens et al. 2011). Because of this, we argue that it is possible that resonantly damped sausage waves have already been observed, albeit in the guise of Alfvén waves; see, for example, Jess et al. (2009).

Lastly, we estimated the damping time for the parameters presented by Morton et al. (2012) and Van Doorsselaere et al. (2011). Interestingly, the expected damping time is very close to the observed damping in QPPs by Kolotkov et al. (2015) that are believed to be fast sausage waves. We find, subject to certain assumptions, that axisymmetric modes appear to be quite important conduits for energy transfer in the solar atmosphere, perhaps even more important than pure Alfvén waves, given that the excitation mechanism for sausage modes in weakly twisted magnetic flux tubes appears to be more readily available than the purely torsional drivers required for Alfvén waves (Giagkiozis et al. 2015).

I.G. would like to acknowledge the Faculty of Science of the University of Sheffield for the SHINE studentship. I.G. also thanks T.V.D. for financial support during his visit to KU Leuven. M.G. is grateful to Belspo's IAP P7/08 CHARM and KU Leuven GOA-2015-014. G.V. and V.F. would like to acknowledge the STFC for funding received (Grant number ST/M000826/1). T.V.D. has received funding from the Odysseus program of the FWO-Vlaanderen and would also like to acknowledge the framework of Belspo's IAP P7/08 CHARM and the GOA-2015-014 of the Research Council of the KU Leuven.