Phase Mixing of Kink MHD Waves in the Solar Corona: Viscous Dissipation and Heating

Aschwanden et al. (1999) and Nakariakov et al. (1999) were the first to identify the spatially resolved kink oscillations of coronal loops using the Transitional Region And Coronal Explorer (TRACE) observations in the 171 ${\AA}$ Fe IX emission lines. Because of the large number of observations of kink oscillations in coronal loops these oscillations are a great seismological tool to estimate the parameters of the solar coronal plasma such as the magnetic field, the plasma density, and the transport coefficients.

An interesting characteristic of kink oscillations in coronal loops is that they damp fast usually within 3-5 periods. Nakariakov et al. (1999) speculated that an anomalously high viscosity or resistivity in the coronal plasma could be responsible for this rapid damping. However, among the proposed mechanisms to justify the rapid damping of the kink waves (see e.g. Goossens et al. 2002; Ofman & Aschwanden 2002; Ruderman & Roberts 2002; Ofman 2005, 2009; Morton & Erdélyi 2009, 2010) resonant absorption is the strongest candidate, since, unlike other theories, resonant absorption is the only proposed mechanism that is able to explain short damping times (of the order of a few periods of the kink oscillation) without invoking anomalous processes or values of the dissipative coefficients many orders of magnitude larger than those expected in the corona. Resonant absorption is an ideal process that does not need strong diffusion to work. However, there is still no observational evidence for resonant absorption. This mechanism was first proposed by Ionson (1978) as a heating mechanism in coronal loops. Since then, many studies have developed the theory of resonant absorption (see e.g. Davila 1987; Sakurai et al. 1991a, 1991b; Goossens et al. 1995; Goossens & Ruderman 1995; Erdélyi 1997; Cally & Andries 2010 among many others). The necessary condition in this mechanism is that the wave frequency lies in the local Alfvén and/or slow frequency continuum. In this situation the energy of the global mode oscillation transfers to the local perturbations in the inhomogeneous regions of the magnetic flux tube. As a result, the amplitude of the perturbations grows at the resonance point and the dissipation mechanisms become important in the resonance layer, where the oscillations make large gradients. Ruderman & Roberts (2002) studied damping of kink oscillations in coronal loops. Considering the effect of viscosity in their analysis, they confirmed the previously obtained numerical result (see e.g. Poedts & Kerner 1991; Tirry & Goossens 1996) that the decay rate of the transverse oscillation is independent of the Reynolds number $R_{v}$ when $R_{v}\gg 1$. They concluded that Reynolds number affects only the perturbations in the resonance layer so that it is not possible to obtain the value of viscosity from the observations of decaying kink oscillations in coronal loops. Resonant absorption has been studied for various complex configurations of the magnetic flux tubes including curvature of the flux tube (Terradas et al. 2006), longitudinal density stratification (Andries et al. 2005; Karami & Asvar 2007; Soler et al. 2011), twisted magnetic field (Karami & Bahari 2010; Ebrahimi & Karami 2016; Ebrahimi & Bahari 2019) and magnetic field expansion (Shukhobodskiy et al. 2018; Howson et al. 2019). For a review on the theory of resonant absorption, see Goossens et al. (2011).

Another consequence of existing a continuum of Alfvén frequencies in the flux tubes may be the phase mixing of the Alfvén waves (Heyvaerts & Priest 1983; Ireland & Priest 1997; Karami & Ebrahimi 2009; Prokopyszyn & Hood 2019). In this mechanism due to inhomogeneity of the local Alfvén phase speed across the background magnetic field the perturbations on different magnetic surfaces become out of phase as travel in the case of propagating wave or oscillate in the case of standing wave. In the developed stage of phase mixing even with a small amount of viscosity or resistivity the dissipation mechanisms become important and could transform the wave energy to heat. Ofman & Aschwanden (2002) suggested that the loop oscillations are dissipated by phase mixing with viscosity of the order $\nu=10^{5.3\pm 3.5}~{}m^{2}s^{-1}$ that is anomalously many order of magnitudes higher than the classical coronal value of the shear viscosity, $\nu=1~{}m^{2}s^{-1}$ (Ofman et al. 1994). It is believed that some small-scale turbulence and structure enhance the viscosity in the coronal loop plasma. Ofman et al. (1994) investigated the effect of viscous stress tensor on the heating of the corona by the resonant absorption and showed that the shear viscosity has the dominant role in the heating process but the compressive viscosity does not have a significant contribution (see also Erdélyi & Goossens 1995, 1996).

Goossens et al. (2014) investigated the nature of kink waves and stated that kink waves do not only involve purely transverse motions of solar magnetic flux tubes, but the velocity field is a spatially and temporally varying sum of both transverse and rotational motion. In an axisymmetric cylindrical flux tube, wave modes can be classified according to the value of the azimuthal wavenumber, $m$. In this paper we study modes with $m=1$, i.e., kink modes. The global kink mode is the only mode that is able to displace transversely (i.e., laterally) the axis of the cylinder. The global mode with $m=1$ is resonantly coupled to Alfvén modes with $m=1$ in the nonuniform layer of the tube. These modes have both radial and azimuthal components of the displacement. Goossens et al. (2014) called the Alfvén modes with $m\neq 0$ as ”rotational modes”. The reason for calling these modes ”rotational” is that their streamlines follow a closed curve, so that a fluid element flowing along one of those streamlines would describe a ”rotation” around a certain point. In the case of $m=1$, the center of the rotation is not located on the axis of the cylinder (as happens for torsional motions with ”m=0”), but at some place in between the axis and the boundary of the tube. Following Goossens et al. (2014), we use the term ”transverse motion” to refer the lateral displacement of the flux tube caused by the global kink mode and the term ”rotational motion” to refer the local Alfvén perturbations inside the tube. Soler & Terradas (2015, hereafter ST2015) investigated the resonant absorption of the kink MHD wave and phase mixing of its perturbations in coronal flux tubes. Using a modal expansion method they showed that the energy of the global kink oscillation of the flux tube is transformed into small-scale rotational motions in the nonuniform boundary of the tube which are eventually subject to the simultaneously occurring phase mixing process. However, ST2015 used the ideal MHD equations and did not consider the effect of dissipation terms in their analysis. At the developed stage of the phase mixing process where the perturbations are highly phase mixed, the dissipation mechanisms could suppress the rate of generating small scales in the system by coupling the perturbations on the neighboring magnetic surface and finally transform the kink wave energy to heat.

In this paper, our aim is to investigate the effect of viscosity on the kink MHD waves and show how viscosity modifies the previous results obtained by ST2015. Section 2 presents the model and the governing MHD equations of motion. In section 3 we apply and extend the mathematical method used by ST2015 and give a solution to the equation of motion. The results are discussed in section 4. Finally we conclude the paper in section 5.

We model a typical coronal loop by a straight cylinder that has a circular cross section of radius $R$. The background plasma density in cylindrical coordinates ($r$, $\varphi$, $z$) is assumed to be as follows

$$\rho(r)=\left\{\begin{array}[]{lll}\rho_{{\rm i}},&r\leqslant r_{1},\\ \frac{\rho_{{\rm i}}}{2}\left[\left(1+\frac{\rho_{{\rm e}}}{\rho_{{\rm i}}}\right)-\left(1-\frac{\rho_{{\rm e}}}{\rho_{{\rm i}}}\right)\sin\left(\frac{\pi}{l}(r-R)\right)\right],&r_{1}<r<r_{2},\\ \rho_{{\rm e}},&r\geqslant r_{2},\end{array}\right.$$

(1)

where $r_{1}=R-l/2$ and $r_{2}=R+l/2$. Here, $l=r_{2}-r_{1}$ is the width of the inhomogeneous region. The background magnetic field is assumed to be constant and aligned with the flux tube axis everywhere, i.e. $\mathbf{B}=B_{{\rm 0}}\hat{z}$ where $B_{0}$ is constant.

The linearized MHD equations for an incompressible plasma with viscosity are as follows

$$\rho(r)\frac{\partial^{2}\boldsymbol{\xi}}{\partial t^{2}}=-\nabla p^{\prime}+\frac{1}{\mu_{0}}(\nabla\times\mathbf{B^{\prime}})\times{\mathbf{B}}+\rho(r)\nu\nabla^{2}\frac{\partial}{\partial t}\boldsymbol{\xi},$$

(2)

$$\mathbf{B^{\prime}}=\nabla\times(\boldsymbol{\xi}\times\mathbf{B}),$$

(3)

$$\nabla\cdot\boldsymbol{\xi}=0,$$

(4)

where $\boldsymbol{\xi}$ is the lagrangian displacement of the plasma, $\mathbf{B^{\prime}}$ and $p^{\prime}$ are the Eulerian perturbations of the magnetic field and plasma pressure, respectively. Here $\mu_{0}$ is the magnetic permeability of the free space and $\nu$ is coefficient of the kinematic shear viscosity which is assumed to be uniform. Using Eq. (4), we can rewrite Eq. (3) as

$$\mathbf{B^{\prime}}=(\mathbf{B}\cdot\nabla)\boldsymbol{\xi}.$$

(5)

Substituting $\mathbf{B^{\prime}}$ from Eq. (5) into (2) yields

$$\rho(r)\frac{\partial^{2}\boldsymbol{\xi}}{\partial t^{2}}=-\nabla P^{\prime}+\frac{1}{\mu_{0}}B_{0}^{2}\frac{\partial^{2}}{\partial z^{2}}\boldsymbol{\xi}+\rho(r)\nu\nabla^{2}\frac{\partial}{\partial t}\boldsymbol{\xi},$$

(6)

or in a more compact form

$$\mathcal{L}_{A}\boldsymbol{\xi}-\rho(r)\nu\nabla^{2}\frac{\partial}{\partial t}\boldsymbol{\xi}=-\nabla P^{\prime},$$

(7)

where $P^{\prime}=p^{\prime}+(\mathbf{B^{\prime}}\cdot\mathbf{B})/\mu_{0}$ is the Eulerian perturbation of the total (gas plus magnetic) pressure and

$$\mathcal{L}_{A}\equiv\rho(r)\frac{\partial^{2}}{\partial t^{2}}-\frac{1}{\mu_{0}}B_{0}^{2}\frac{\partial^{2}}{\partial z^{2}},$$

(8)

is the Alfvén operator. In the absence of viscosity, Eq. (7) reduces to equation (8) of ST2015. Since the equilibrium quantities are only functions of $r$, the perturbations can be Fourier-analyzed with respect to the $\varphi$ and $z$ coordinates. Hence,

$$\begin{split}P^{\prime}=&P^{\prime}(r,t)~{}e^{i(m\varphi+k_{z}z)},\\ \boldsymbol{\xi}=&\boldsymbol{\xi}(r,t)~{}e^{i(m\varphi+k_{z}z)},\end{split}$$

(9)

where $m$ is the azimuthal mode number and $k_{z}$ is the axial wave number. The three components of Eq. (7) and the incompressibility condition (Eq. 4) form a system of 4 independent equations for $\xi_{r}$, $\xi_{\varphi}$, $\xi_{z}$ and $P^{\prime}$ as follows

	$\displaystyle\mathcal{L}_{A}\xi_{r}-\rho(r)\nu\frac{\partial}{\partial t}\left(\nabla^{2}\xi_{r}-\frac{\xi_{r}}{r^{2}}-\frac{2}{r^{2}}\frac{\partial}{\partial\varphi}\xi_{\varphi}\right)=-\frac{\partial}{\partial r}P^{\prime},$		(10)
	$\displaystyle\mathcal{L}_{A}\xi_{\varphi}-\rho(r)\nu\frac{\partial}{\partial t}\left(\nabla^{2}\xi_{\varphi}-\frac{\xi_{\varphi}}{r^{2}}+\frac{2}{r^{2}}\frac{\partial}{\partial\varphi}\xi_{r}\right)=-\frac{1}{r}\frac{\partial}{\partial\varphi}P^{\prime},$		(11)
	$\displaystyle\mathcal{L}_{A}\xi_{z}-\rho(r)\nu\frac{\partial}{\partial t}\nabla^{2}\xi_{z}=-\frac{\partial}{\partial z}P^{\prime},$		(12)
	$\displaystyle\frac{1}{r}\frac{\partial}{\partial r}\left(r\xi_{r}\right)+\frac{1}{r}\frac{\partial}{\partial\varphi}\xi_{\varphi}+\frac{\partial}{\partial z}\xi_{z}=0.$		(13)

Here, we apply $\frac{1}{r}\frac{\partial}{\partial\varphi}$ and $\frac{\partial}{\partial z}$ from left on Eqs. (11) and (12), respectively and add the resulting equations together. After that with the help of Eqs. (9) and (13) we obtain $P^{\prime}$ in terms of $\xi_{r}$ and $\xi_{\varphi}$ as

$$\begin{split}P^{\prime}&=\frac{-1}{m^{2}/r^{2}+k_{z}^{2}}\left[\mathcal{L}_{A}\left(\frac{1}{r}\frac{\partial(r\xi_{r})}{\partial r}\right)-\right.\\ &\left.\rho(r)\nu\frac{\partial}{\partial t}\left(\nabla^{2}\left(\frac{1}{r}\frac{\partial}{\partial r}(r\xi_{r})\right)+\left(\frac{2}{r^{3}}-\frac{2}{r^{2}}\frac{\partial}{\partial r}\right)\frac{\partial}{\partial\varphi}\xi_{\varphi}-\frac{2}{r^{3}}\frac{\partial^{2}}{\partial\varphi^{2}}\xi_{r}\right)\right].\end{split}$$

(14)

Now we use the thin tube (TT) approximation in which the wavelength of the kink waves, $\lambda$, is much larger than the radius of the cross section of the flux tube, $R=(r_{1}+r_{2})/2$, i.e. $Rk_{z}\ll 1$ or $R\ll\lambda$. The first two terms of Eq. (13) have magnitudes of the order $\xi_{0}/R$ where $\xi_{0}$ is a typical value for the Lagrangian displacement. Since the characteristic length scale in the $z$ direction is $\lambda$, the order of the magnitude of the third term of Eq. (13) is equal to $\xi_{0}/\lambda$. Hence, in TT approximation neglecting the third term of Eq. (13) with respect to the first two terms, yields

$$\frac{\partial}{\partial\varphi}\xi_{\varphi}\simeq-\frac{\partial}{\partial r}(r\xi_{r}).$$

(15)

We use Eq. (15) to obtain a single equation for $\xi_{r}$. Since we use this approximation, it is not possible to obtain $\xi_{z}$ with respect to $\xi_{r}$ from Eqs. (10)-(13). However, in the TT approximation, Goossens et al. (2009) showed that the longitudinal component of the displacement, $\xi_{z}$, is always much smaller than the other components and the dominant motion is in the horizontal plane normal to the background magnetic field. Substituting Eq. (15) in Eqs. (10) and (14) and eliminating $P^{\prime}$ from the resulting equations, gives the equation for $\xi_{r}$ as follows

$$\mathcal{L}_{A}\mathcal{L}_{s}\xi_{r}+\left(\frac{m^{2}}{r^{2}}+k_{z}^{2}\right)\frac{d\rho(r)}{dr}\frac{\partial^{2}}{\partial t^{2}}\frac{1}{r}\frac{\partial}{\partial r}(r\xi_{r})=\nu\frac{\partial}{\partial t}\mathcal{L}_{d}\xi_{r},$$

(16)

where $\mathcal{L}_{s}$ is the surface wave operator (see ST2015 for more details) and $\mathcal{L}_{d}$ is the viscous damping operator which are defined as follows

$$\mathcal{L}_{s}\equiv\left(k_{z}^{2}+\frac{m^{2}}{r^{2}}\right)\frac{\partial^{2}}{\partial r^{2}}+\frac{1}{r}\left(k_{z}^{2}+\frac{3m^{2}}{r^{2}}\right)\frac{\partial}{\partial r}-\frac{1}{r^{2}}\left(k_{z}^{2}-\frac{m^{2}}{r^{2}}\right)-\left(k_{z}^{2}+\frac{m^{2}}{r^{2}}\right)^{2},$$

(17)

$$\begin{split}\mathcal{L}_{d}\equiv&\left(\frac{m^{2}}{r^{2}}+k_{z}^{2}\right)^{2}\left[\rho(r)\left(\nabla^{2}-\frac{1}{r^{2}}+\frac{2}{r^{2}}\left(1+r\frac{\partial}{\partial r}\right)\right)\right.\\ &\left.-\frac{\partial}{\partial r}\left(\frac{\rho(r)}{m^{2}/r^{2}+k_{z}^{2}}\left[\frac{\partial^{3}}{\partial r^{3}}+\frac{4}{r}\frac{\partial^{2}}{\partial r^{2}}-\left(\frac{m^{2}}{r^{2}}+k_{z}^{2}\right)\frac{\partial}{\partial r}-\frac{1}{r}\left(\frac{m^{2}+1}{r^{2}}+k_{z}^{2}\right)\right]\right)\right].\end{split}$$

(18)

In the absence of viscosity, Eq. (16) consistently reverts to Eq. (16) of ST2015. From Eq. (15) we find that $\xi_{\varphi}$ is related to $\xi_{r}$ as

$$i\xi_{\varphi}\simeq-\frac{1}{m}\frac{\partial}{\partial r}(r\xi_{r}).$$

(19)

In this relation, the factor $i$ accounts for a phase difference of $\pi/2$ between $\xi_{\varphi}$ and $\xi_{r}$. So, for convenience, in order to avoid imaginary terms in the calculations, in the rest of the paper we redefine $i\xi_{\varphi}$ as $\xi_{\varphi}$.

In order to solve Eq. (41), the density ratio of the flux tube is considered to be $\rho_{\rm i}/\rho_{\rm e}=5$. The observational values of this parameter has been reported to be in the range [2, 10] (Aschwanden et al. 2003). The azimuthal mode number representing the kink waves is $m=\pm 1$. For the model considered in this paper, the results for both $m=+1$ and $m=-1$ are the same. Here, we take $m=+1$. Since we use the TT approximation, the longitudinal wavenumber is assumed to be $k_{z}=\frac{\pi}{100R}$. For the thickness of the inhomogeneous layer we consider two cases $l/R=0.2$ and $l/R=1$ that represent a thin and a thick layer, respectively. To consider the effect of viscosity, it is appropriate to use a dimensionless quantity, the viscous Reynolds number which is defined as

$$R_{v}\equiv\frac{lv_{\rm Ai}}{\nu}.$$

(54)

Here, the Alfvén speed $v_{\rm Ai}=B_{0}/\sqrt{\mu_{0}\rho_{\rm i}}$ is the characteristic speed for the propagation of kink waves in the flux tube. Traditional value of the coefficient of shear viscosity in the solar corona is of the order $\nu=1~{}m^{2}/s$ (see Ofman et al. 1994 and references therein). With $l=10^{6}m$ and $v_{\rm Ai}=10^{6}m/s$, the corresponding Reynolds number is $R_{v}=10^{12}$. Due to computational limitations, we are forced to consider smaller $R_{v}$ in our results. However, the conclusions we obtain can be easily generalized to the case of larger $R_{v}$. When appropriate, we shall stress what differences would appear if more realistic values of $R_{v}$ were considered.

Solving the generalized eigenvalue problem (41) results to 2N eigenvalues, $\sigma\equiv-i\tilde{\omega}$ where $\tilde{\omega}\equiv\omega+i\gamma$ is the complex eigenfrequency. These eigenvalues are real or come in pairs $(\sigma,\sigma^{*})$ where $\sigma^{*}$ is the complex conjugate of $\sigma$. Figure 1 shows the $\omega_{n}\geqslant 0$ part of the spectrum of the complex eigenvalues $\sigma_{n}\equiv-i\omega_{n}+\gamma_{n}$ with $N=101$ for $l/R=0.2$ and $l/R=1$ and Reynolds numbers $R_{v}=10^{6},10^{7}$. Here, $\omega_{n}$ and $\gamma_{n}$ are the frequency and the damping rate of the $n$’th eigenmode, respectively. In the figure, $\omega_{\rm Ai}=k_{z}B_{0}/\sqrt{\mu_{0}\rho_{\rm i}}$ and $\omega_{\rm Ae}=k_{z}B_{0}/\sqrt{\mu_{0}\rho_{\rm e}}$ are the Alfvén frequencies inside and outside of the tube, respectively. Note that the frequencies and damping rates are in units of the internal Alfvén frequency, $\omega_{\rm Ai}$. The complex spectrum has the typical three-branch structure found in previous papers that computed the resistive Alfvén spectrum in similar configurations (see Poedts & Kerner 1991). It is clear from Figure 1 that the smaller the value of the Reynolds number (larger coefficient of viscosity) the more number of modes with $\omega=0$ are present in the complex spectrum. This result is in well agreement with previous results obtained in resistive MHD (see e.g. Van Doorsselaere and Poedts 2007). An interesting result is that similar to the resistive MHD analysis (see e.g. Poedts and Kerner 1991) in viscous MHD, one of the eigenvalues of the complex spectrum could be identified as the damped quasi-mode solution (global mode) of ideal MHD. The real part of the frequency of this solution is the kink mode frequency, and the imaginary part is its corresponding damping rate due to the resonant absorption mechanism. Following Ruderman and Roberts (2002) and Goossens et al. (2002), the quasi-mode frequency $\omega_{\rm qm}$ and damping rate $\gamma_{\rm qm}$ for MHD kink waves in the TT and thin boundary (TB) ($l/R\ll 1$) approximations are

$$\omega_{\rm qm}\simeq\omega_{\rm k}=\sqrt{\frac{2\rho_{\rm i}}{\rho_{\rm i}+\rho_{\rm e}}}\omega_{\rm Ai},$$

(55)

$$\gamma_{\rm qm}=-\omega_{\rm k}\frac{l}{4R}\frac{\rho_{\rm i}-\rho_{\rm e}}{\rho_{\rm i}+\rho_{\rm e}}.$$

(56)

From these equations for $\rho_{\rm i}/\rho_{\rm e}=5$ the quasi-mode frequency is obtained as $\omega_{qm}=1.29\omega_{\rm Ai}$. Also the quasi-mode damping rates for $l/R=0.2$ and $l/R=1$ are $\gamma_{\rm qm}=-0.043\omega_{\rm Ai}$ and $\gamma_{\rm qm}=-0.215\omega_{\rm Ai}$, respectively. As illustrated in Figure 1, the singled out mode matches the quasi-mode solution especially for the thin transitional layer case, $l/R=0.2$, and larger values of the Reynolds number. For instance, for $l/R=0.2$ and $R_{v}=10^{7}$ the singled out eigenfrequency is $\tilde{\omega}=1.29-0.044i$ that is in well agreement with the result obtained from the above analytic approximations. For $l/R=1$ and $R_{v}=10^{7}$ the corresponding eigenfrequency is $\tilde{\omega}=1.33-0.269i$. Hence, in the case of thick transitional layer, the numerical results and the analytic approximations deviate more from each other. This is consistent with the fact that the analytic approximations are only strictly valid when $l/R\ll 1$. For a study of the validity of the approximations beyond its theoretical range of applicability see Soler et al. (2014). These results show that for the kink waves, the quasi-mode solution of the ideal MHD can be identified as an eigenmode of the viscous MHD spectrum. This correspondence is very clear in the case of thin nonuniform layers and not very large Reynolds numbers. However, as the thickness of the layer or the Reynolds numbers increase, the identification of the quasi-mode in the dissipative spectrum is more confusing since the quasi-mode gets embedded in one of the branches of the spectrum and becomes indistinguishable from an ordinary dissipative Alfvén mode (see discussions on this issue in Van Doorsselaere & Poedts 2007 and Soler et al. 2013). A detailed analysis of the peculiar behaviour of the quasi-mode in the complex spectrum is beyond the purpose of the present paper.

Once the dissipative spectrum is computed, the time-dependent behavior of the perturbations is obtained by the superposition of all the modes in the spectrum according to the prescribed initial condition, which represents a transverse, i.e., kink displacement of the whole tube (see Section 3.3). In short, the evolution of the subsequent global kink oscillation is determined by two simultaneously working mechanisms: resonant absorption and phase mixing. On the one hand, resonant absorption is responsible for a radial flux of energy towards the nonuniform layer of the flux tube, and its net effect is producing the damping of the global kink oscillation. As a result, the amplitude of the displacement at the tube axis decreases in time. On the other hand, the energy accumulated at the nonuniform layer because of resonant absorption drives local Alfvén waves with the same azimuthal symmetry as the original kink oscillation, i.e., $m=1$. Although these Alfvén waves have both radial and azimuthal components of the displacement, they are largely polarized in the azimuthal direction. The Alfvén waves undergo the process of phase mixing because of the spatially-dependent Alfvén velocity. This causes the building up of small scales in the nonuniform layer and the subsequent energy cascade to these small scales. Then, viscous dissipation becomes important. In the following paragraphs, we analyze the dynamics we have just summarized.

Figures 2 and 3 show the evolution of $\xi_{r}$ and $\xi_{\varphi}$, respectively. Time is in units of the period of the kink oscillation in TTTB approximations, $P_{\rm k}=2\pi/\omega_{\rm k}$. In Figures 2 and 3 the solid black curves represent the results previously obtained by ST2015 in the absence of viscosity, i.e. $R_{v}=\infty$. The blue dashed and red dashed-dotted curves are for $R_{v}=10^{7}$ and $R_{v}=10^{6}$, respectively. Figures reveal that in the presence of viscosity, unlike the results obtained in ideal MHD (see Soler and Terradas 2015; Ebrahimi et al. 2017) perturbations are not allowed to be phase mixed indefinitely. The smaller the Reynolds number, the quicker the dissipative solution departs from the ideal solution. The existence of viscosity cause to coupling of the perturbations on the neighboring magnetic surfaces. This effect suppresses the phase mixing when a certain spatial scale is reached and transforms the total (kinetic plus magnetic) energy of the perturbations to heat via dissipation.

To illustrate the flux of the total energy of the kink wave from the internal and external regions to the inhomogeneous region where it is finally dissipated, we calculate the total energy density of the perturbations as

$$E(r,t)=\frac{1}{2}\left(\rho\left|\frac{\partial\boldsymbol{\xi}}{\partial t}\right|^{2}+\frac{1}{\mu}\left|\mathbf{B^{\prime}}\right|^{2}\right),$$

(57)

where $\mathbf{B^{\prime}}$ is obtained from Eq. (5). Figure 4 illustrates the evolution of the total energy density as a function of $r$. In the absence of viscosity after a few time periods the whole energy of the kink wave concentrates in a narrow layer in the inhomogeneous region of the flux tube. The role of viscosity is to dissipate this concentrated energy. Although, these two processes i.e. flow of energy to the inhomogeneous layer and the dissipation occur simultaneously but they are caused by independent and physically different mechanisms. In order to have a better illustration of the flow of energy of the kink wave from internal and external regions to the inhomogeneous layer and its dissipation, it is appropriate to calculate the integrated total energy of the perturbations in the interior, transitional layer and exterior of the flux tube as a function of time, respectively, as follows

	$\displaystyle E_{\rm{in}}=\int_{0}^{r_{1}}\frac{1}{2}\left(\rho\left|\frac{\partial\boldsymbol{\xi}}{\partial t}\right|^{2}+\frac{1}{\mu}\left|\mathbf{B^{\prime}}\right|^{2}\right)r~{}{\rm d}r,$
	$\displaystyle E_{\rm{tr}}=\int_{r_{1}}^{r_{2}}\frac{1}{2}\left(\rho\left|\frac{\partial\boldsymbol{\xi}}{\partial t}\right|^{2}+\frac{1}{\mu}\left|\mathbf{B^{\prime}}\right|^{2}\right)r~{}{\rm d}r,$		(58)
	$\displaystyle E_{\rm{ex}}=\int_{r_{2}}^{\infty}\frac{1}{2}\left(\rho\left|\frac{\partial\boldsymbol{\xi}}{\partial t}\right|^{2}+\frac{1}{\mu}\left|\mathbf{B^{\prime}}\right|^{2}\right)r~{}{\rm d}r.$

Note that in order to compute the third integral of Eq. (4) we must replace the upper limit of the integral, $\infty$, with a sufficiently large radius where the amplitudes of the perturbations are already negligible. Hence, we check the convergence of the total energy by taking a series of large radii as the upper limit of the third integral in Eq. (4). Results show that at $r=20R$ the total energy converges to the desired level of accuracy. Figure 5 shows the integrated energy in these three regions as a function of time for $R_{v}=\infty,~{}10^{7}~{}\&~{}10^{6}$. The solid, dashed and dotted lines denote the integrated energies of the internal, inhomogeneous and external regions, respectively. The black, blue and red colors represent $R_{v}=\infty$ (ideal MHD), $R_{v}=10^{7}$ and $R_{v}=10^{6}$, respectively. Note that in the figure, the plots of the internal and external energies (i.e. the solid and dotted lines) for $R_{v}=\infty$ and $R_{v}=10^{7}\&10^{6}$ coincide with each other which reveals that for large Reynolds numbers the existence of viscosity does not affect the energy flow from internal and external regions to the inhomogeneous layer. In other words, the resonant absorption process is not affected by viscosity. Figure 5 also shows that in the presence of viscosity, the energy in the inhomogeneous layer increases in the initial stage of the evolution and reaches a maximum after a few time periods. After that, the energy decreases, since phase mixing in the inhomogeneous layer enhances the viscous dissipation mechanism. The temporal behaviour of the total energy of the kink wave has been illustrated in Figure 6 for $l/R=0.2~{}\&~{}1$ and $R_{v}=\infty,~{}10^{7}~{}\&~{}10^{6}$. As figure shows, the energy of the kink wave in the absence of viscosity (black line in the figure) is conserved. The results show that, considering a flux tube with thin transitional layer ($l/R=0.2$), for $R_{v}=10^{7}~{}\&~{}10^{6}$ the energy of the kink wave decreases to $1/e$ of its initial energy after $t\simeq 5.8P_{k}$ and $t\simeq 4.4P_{k}$, respectively. For thick transitional layer ($l/R=1$) the corresponding values are $t\simeq 5.2P_{k}$ and $t\simeq 2.7P_{k}$.

Figure 7 shows the radial component of the displacement at the axis of the flux tube. Interestingly the results for $R_{v}=\infty,10^{7}\&10^{6}$ are the same. Note that in the observations of kink waves in coronal flux tubes, we measure the displacement of the axis but detecting the rotational motions related to the kink waves in the inhomogeneous boundary of the flux tubes is not an easy task. Comparing the results illustrated in Figures 6 and 7, it is clear that although the temporal behaviour of the displacement of the flux tube axis is the same in both cases of ideal and viscous MHD (with large Reynolds number), the total energy of the kink waves which is conserved for $R_{v}=\infty$, decays in a few periods for $R_{v}=10^{7}$ and $10^{6}$. Hence, one can conclude that assuming a large Reynolds number for the coronal plasma the global damping of kink waves reported in the observations is due to converting the transverse motion to rotational perturbations in the inhomogeneous layer of the flux tube, i.e. the resonant absorption process, and is not related to the existence of viscosity.

In order to illustrate the rate of plasma heating in the inhomogeneous region due to viscosity we use the viscous dissipation function that for an incompressible plasma with constant viscosity is as follows (White 1991)

$$\Phi=\rho\nu\left[2\left(\varepsilon_{rr}^{2}+\varepsilon_{\varphi\varphi}^{2}+\varepsilon_{zz}^{2}\right)+\varepsilon_{r\varphi}^{2}+\varepsilon_{rz}^{2}+\varepsilon_{\varphi z}^{2}\right],$$

(59)

where $\varepsilon_{ij}$ with $i,j=r,\varphi,z$ are components of the strain rate tensor defined as follows

	$\displaystyle\varepsilon_{rr}$	$\displaystyle=$	$\displaystyle\frac{\partial^{2}\xi_{r}}{\partial r\partial t},$
	$\displaystyle\varepsilon_{\varphi\varphi}$	$\displaystyle=$	$\displaystyle\frac{1}{r}\frac{\partial^{2}\xi_{\varphi}}{\partial\varphi\partial t}+\frac{1}{r}\frac{\partial\xi_{r}}{\partial t},$
	$\displaystyle\varepsilon_{zz}$	$\displaystyle=$	$\displaystyle\frac{\partial^{2}\xi_{z}}{\partial z\partial t},$
	$\displaystyle\varepsilon_{r\varphi}$	$\displaystyle=$	$\displaystyle\frac{1}{r}\frac{\partial^{2}\xi_{r}}{\partial\varphi\partial t}+\frac{\partial^{2}\xi_{\varphi}}{\partial r\partial t}-\frac{1}{r}\frac{\partial\xi_{\varphi}}{\partial t},$
	$\displaystyle\varepsilon_{rz}$	$\displaystyle=$	$\displaystyle\frac{\partial^{2}\xi_{r}}{\partial z\partial t}+\frac{\partial^{2}\xi_{z}}{\partial r\partial t},$
	$\displaystyle\varepsilon_{\varphi z}$	$\displaystyle=$	$\displaystyle\frac{1}{r}\frac{\partial^{2}\xi_{z}}{\partial\varphi\partial t}+\frac{\partial^{2}\xi_{\varphi}}{\partial z\partial t}.$

The dissipation function, $\Phi$, is the work done by the viscous stresses on an element of plasma per unit volume per unit time. Here, we remove the $\varphi$ and $z$ dependency of the dissipation function by integrating $\Phi$ in these directions, i.e.

$$\Phi^{\prime}=\int_{0}^{\lambda}\int_{0}^{2\pi}\Phi rd\varphi dz=\pi\lambda r\Phi,$$

(60)

where, $\lambda=2\pi/k_{z}$ is the wavelength. Figure 8 shows the contour plot of $\Phi^{\prime}$ in $r-t$ plane for $l/R=0.2$ and $l/R=1$ with $R_{v}=10^{6}$ and $R_{v}=10^{7}$. As figure shows, the heating rate has an oblique oscillatory pattern. The obliquity of the contours of $\Phi^{\prime}$ is due to the phase mixing of the perturbations. In order to obtain the time that the dissipation in the inhomogeneous region reaches its maximum we have calculated the integral of $\Phi^{\prime}$ over the range $[r_{1},r_{2}]$ that is a function of time. For $R_{v}=10^{6}$, considering $l/R=0.2$ and $l/R=1$ we obtain the peak time of the integrated dissipation as $t_{peak}=2.75P_{k}$ and $t_{peak}=2.1P_{k}$, respectively. The corresponding values for $R_{v}=10^{7}$ are $t_{peak}=4.5P_{k}$ and $t_{peak}=4.1P_{k}$. For $l/R=0.2$ and $R_{v}=10^{8}$ we get $t_{peak}=8P_{k}$. So, for the larger value of the Reynolds number the system needs to be more phase mixed for viscosity to reach its maximum efficiency as a heating mechanism for the plasma. Obviously, if the Reynolds number is further increased to the expected value in the solar corona, the peak efficiency of dissipation happens in a much later time corresponding to many periods of the original kink oscillation. Since the observationally reported damping time of kink oscillations corresponds to a few periods (as the theory of resonant absorption correctly predicts), we conclude that viscous heating becomes of significance only after the global kink oscillation is damped, i.e., no heating is expected during the damping of the global kink oscillation in realistic coronal conditions.

Here, we investigated the effect of viscosity on the evolution of MHD kink waves in coronal flux tubes. We modelled a magnetic flux tube by a straight magnetic cylinder. Plasma density inside and outside the tube is constant with different values. The interior and exterior of the tube are connected by an inhomogeneous transitional layer in which the plasma density varies smoothly with a sinusoidal profile from the internal value to the external one. The background magnetic field is aligned with tube axis and has constant magnitude everywhere. We neglected the role of viscosity in the constant density regions since in the limit of small viscosities the dissipation is only important in the inhomogeneous region where the phase mixing process is at work. Using the modal expansion technique (Cally 1991; Soler & Terradas 2015) we solved the viscous MHD equations of motion in thin tube approximation and obtained the spatio-temporal behaviour of the perturbations in the flux tube. We considered both the cases of thin and thick inhomogeneous layers in our analysis.

We obtained the spectrum of the complex eigenfrequencies of the Alfvén discrete modes in the inhomogeneous layer. In the spectrum, one of the eigenfrequencies could be identified as the quasi-mode solution of the kink waves by the resonant absorption mechanism.

To investigate the effect of viscosity on the kink waves, we obtained the temporal and spatial behaviour of the perturbations in the interior, inhomogeneous region and exterior of the flux tube. Our results showed that for both cases of thin and thick inhomogeneous layers, considering the viscosity in the system does not affect the transverse motion of the flux tube axis and its decay rate for large Reynolds numbers ($R_{v}\geq 10^{6}$) confirming the previously results obtained by e.g. Ruderman & Roberts (2002) and Goossens et al. (2002). This result confirms that the fast damping of the kink oscillations in coronal loops could be a consequence of resonant absorption mechanism which despite the existence of viscosity naturally results to changing the behaviour of kink mode from mainly transverse motion to rotational motion (Goossens et al. 2014). However, even for small viscosities (relevant for the coronal plasma) the viscous dissipation is important in the developed stage of phase mixing of the perturbations in the inhomogeneous region. Our results show that viscosity eventually suppresses the rate of phase mixing of the perturbations by coupling the neighboring magnetic surfaces and transforming their energy to heat.

In order to investigate the effect of viscosity on cascading the total (kinetic plus magnetic) energy of the kink wave to the inhomogeneous layer of the flux tube, we obtained the total energy density as well as the integrated total energy as a function of $r$ and $t$. As in the ideal MHD case, in presence of viscosity the energy of the kink wave tends to concentrate in a narrow region in the inhomogeneous layer of the flux tube but it is not allowed to achieve its maximum value obtained in the ideal case since the dissipation mechanism is at work and decreases the energy in time. However, the small amount of viscosity considered in our analysis does not affect the energy flow from interior and exterior of the tube to the inhomogeneous region. Temporal behavior of the total energy showed that for both the thin and thick inhomogeneous region with Reynolds numbers of the order $10^{6}$, $10^{7}$ and $10^{8}$ the energy of the kink wave decays to heat the plasma within a few periods. However, if larger and more realistic values of the Reynolds numbers were used, heating would happen much after the observable kink oscillation is completely damped. These conclusions agree with the simple estimations provided in the research note by Terradas & Arregui (2018), although these authors considered resistive heating instead of viscous heating.

We also studied the efficiency of heating duo to viscosity by calculating the dissipation function in the inhomogeneous region. The obtained results show that at initial stage of the evolution the dissipation increases with time and reaches a maximum level after a few periods (2-8 periods for $R_{v}=10^{6}-10^{8}$ but much later for realistic $R_{v}$) in a narrow layer near the boundary of the flux tube. After that the dissipation decreases since the energy budget provided by the initial value problem considered in this paper is finite.

In summary, we showed that viscosity, even in a small amount, can have a significant impact on the later evolution of phase mixing by the suppressing the generation of small scales and transforming the energy of the wave to heat. Reynolds number larger than the values considered in this paper needs more massive numerical computation with the mathematical approach presented here to obtain the correct set of complex eigenvalues of the damped Alfvén discrete modes which could be a subject of future work.

We finally note that this work is based on linear theory. Nonlinear effects may modify somehow these results, since presumably important ingredients as, e.g., the triggering of Kelvin-Helmholtz instabilities are absent from our study (see e.g. Terradas et al. 2008). The nonlinear evolution should be necessarily investigated with high-resolution dissipative MHD simulations.

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