The effect of a twisted magnetic field on the phase mixing of the kink magnetohydrodynamic waves in coronal loops
Zanyar Ebrahimi Kayoomars Karami Roberto Soler




Zanyar Ebrahimi
z.ebrahimi@uok.ac.ir


Department of Physics, University of Kurdistan, Pasdaran Street, P.O. Box 66177-15175, Sanandaj, Iran


Department of Physics, University of Kurdistan, Pasdaran Street, P.O. Box 66177-15175, Sanandaj, Iran


Departament de Física, Universitat de les Illes Balears, E-07122, Palma de Mallorca, Spain
Institut d'Aplicacions Computacionals de Codi Comunitari (IAC^3), Universitat de les Illes Balears, E-07122, Palma de Mallorca, Spain




There are observational evidences for the existence of twisted magnetic field in the solar corona. This inspires us to investigate the effect of a twisted magnetic field on the evolution of magnetohydrodynamic (MHD) kink waves in coronal loops. To this aim, we solve the incompressible linearized MHD equations in a magnetically twisted nonuniform coronal flux tube in the limit of long wavelengths. Our results show that a twisted magnetic field can enhance or diminish the rate of phase-mixing of the Alfvén continuum modes and the decay rate of the global kink oscillation depending on the twist model and the sign of the longitudinal (k_z) and azimuthal (m) wavenumbers. Also our results confirm that in the presence of a twisted magnetic field, when the sign of one of the two wavenumbers m and k_z is changed, the symmetry with respect to the propagation direction is broken. Even a small amount of twist can have an important impact on the process of energy cascade to small scales.







§ INTRODUCTION
 
Transverse oscillations of the solar coronal loops are one of the greatest seismological tools to extract or approximate the unknown parameters of the solar corona such as the magnetic field, the plasma density and the transport coefficients. Aschwanden et al. (1999) and Nakariakov et al. (1999) were first to report the observation of the transverse oscillations in the coronal loops using the Transition Region and Coronal Explorer (TRACE) telescope on 1998 July 14 in the 171-Å Fe IX emission lines. Nakariakov et al. (1999) indicated that the oscillations are strongly damped and the ratio of the damping time to the period of the oscillation is around 3-5. This observation was identified as a standing kink MHD wave in a magnetic flux tube. However, the mechanism proposed in Nakariakov et al. (1999) to explain the damping involved the assumption of unrealistically large diffusion coefficients. A more satisfactory physical interpretation of the damping was given by Ruderman & Roberts (2002).

 Among the suggested mechanisms responsible for the strong damping of the coronal loop oscillations (e.g. Ruderman and Roberts 2002; Ofman 2005, 2009; Morton and Erdélyi 2009), resonant absorption of the MHD waves, that was established first by Ionson (1978), is a strong candidate. Several works developed this theory (e.g. Davila 1987; Sakurai, Goossens & Hollweg 1991a,b; Goossens et al. 1995; Goossens & Ruderman 1995; Erdélyi 1997; Cally & Andries 2010). The necessary condition for the resonant absorption is a continuum of Alfvén or slow frequency across the loop (Ionson 1978; Hollweg 1984, 1987; Davila 1987; Sakurai, Goossens & Hollweg 1991a). Resonant absorption occurs when the frequency of the global MHD mode matches at least with one of the frequencies of the background Alfvén or slow continuum at a location called resonance point. As a result, the energy of the global MHD mode transfers to the local Alfvén modes in a layer around the resonance point, named resonance layer (Lee & Roberts 1986; see also Goossens et al. 2013; Soler & Terradas 2015). In the absence of dissipation mechanisms, the amplitude of the oscillations diverges at the resonance point. Dissipation is important in the resonance layer where the oscillations make large gradients. The background Alfvén or slow continuum can be due to the variation of the plasma density (e.g. Davila 1987; Ofman, Davila & Steinolfson 1994; Ruderman & Roberts 2002; Terradas, Oliver & Ballester 2006; Soler & Terradas 2015), twisted magnetic field (Ebrahimi & Karami 2016) or both of them together (Karami & Bahari 2010; Giagkiozis et al. 2016). There are a variety of theoretical works related to the damping of the coronal loop oscillations based on the theory of resonant absorption of MHD waves (e.g. Ruderman & Roberts 2002; Goossens, Andries & Aschwanden 2002; Van Doorsselaere et al. 2004; Andries et al. 2005; Terradas, Oliver & Ballester 2006; Goossens et al. 2009; Karami, Nasiri & Amiri 2009; Karami & Bahari 2010; Soler et al. 2013; Soler & Terradas 2015; Ebrahimi & Karami 2016; Jung Yu & Van Doorsselaere 2016; Giagkiozis et al. 2016). For a good review about the theory of resonant absorption, see also Goossens et al. (2011).

Ruderman & Roberts (2002) studied the resonant absorption of kink waves in coronal loops. They suggested that only the loops with transverse density inhomogeneities on a small scale compared to the loop thickness are able to support coherent oscillations and consequently become observable. Safari et al. (2006) investigated the resonant absorption of MHD waves in coronal loops and found that as the longitudinal mode number increases, the maximum value of the wave amplitude moves away from the inhomogeneous region towards the loop axis and as a result the efficiency of the process of resonant absorption decreases.

 Goossens et al. (2014) investigated the transverse and torsional motions of MHD kink waves in coronal loops. They showed that the kink waves are not just transverse motions of coronal loops, but the velocity field of the kink waves involves both the transverse and torsional motions. Soler & Terradas (2015) (hereafter ST2015) investigated the evolution of the MHD kink wave in a coronal loop by solving an initial value problem. Inspired by Cally (1991), they showed that the MHD kink wave can be expressed as a superposition of Alfvén continuum modes. They showed that in the presence of an Alfvén frequency continuum made by the variation of the plasma density across the loop, the energy of the global kink wave transfers to the phase mixed azimuthal perturbations of the local Alfvén waves in the inhomogeneous layer.

An interesting property of the coronal structures is that they can have a twisted magnetic field. Chae et al. (2000) stated that in order to have torsional motions in coronal loops, the magnetic field of the loop should be twisted around the loop axis. Chae & Moon (2005) assumed that the
constriction of plasma (i.e. ∂ p/∂ r≠ 0, where p is the plasma pressure) is due to the magnetic tension of the
azimuthal component of the magnetic field. Using this, they found that for a specific observed coronal loop the magnetic twist on the loop axis is about 1.5 π. The existence of magnetic field twist in coronal structures has been reported in several observations (e.g. Kwon & Chae 2008; Aschwanden et al. 2012; Thalmann et al. 2014; Wang et al. 2015). For instance, Kwon & Chae (2008) using the TRACE 171 Å observations in several coronal loops,  reported that the number of twist turns, N_twist, have values in the range [0.11, 0.87]. These values, are small enough to let a typical coronal loop to be kink stable (for more details, see section <ref>).

There are ample theoretical works on the role of the magnetic twist in the MHD oscillations of the coronal loops (e.g. Bennett, Roberts & Narain 1999; Erdélyi & Carter 2006; Erdélyi & Fedun 2006, 2007, 2010; Carter & Erdélyi 2008; Ruderman 2007, 2015; Karami & Barin 2009; Karami & Bahari 2010, 2012; Terradas & Goossens 2012; Ruderman & Terradas 2015; Ebrahimi & Karami 2016).

Sakurai, Goossens & Hollweg (1991a) investigated resonant absorption in twisted flux tubes and obtained jump conditions of the perturbations across the resonance layer. Using the jump conditions, there is no need to solve dissipative MHD equations in the resonance layer and one can connect ideal MHD solutions of the left and the right sides of the resonance layer (see also Sakurai, Goossens & Hollweg 1991b; Goossens, Hollweg & Sakurai 1992). Karami & Bahari (2010) studied the effect of a twisted magnetic field on the resonant absorption of MHD waves in a coronal flux tube. They showed that when the amount of the magnetic twist is increased, the frequency, the damping rate and the ratio of the frequency to the damping rate increase and the period ratio of the fundamental mode to the first overtone mode decreases from its canonical value. Terradas & Goossens (2012) investigated the MHD kink oscillations of the coronal loops in the presence of magnetic field twist. Solving the MHD equations numerically, they showed that in the presence of magnetic twist for a given value of longitudinal wavenumber, k_z, the quasi-mode frequency of MHD kink waves has different values for different signs of the azimuthal mode number (m=±1). They found that when k_z>0, the frequency for m=1 should be larger than that of the solution without twist. Conversely, the frequency for m=-1 should be smaller. They also showed that for a given value of the azimuthal mode number (m=-1 or m=+1), the frequency for different signs of k_z has different values. Hence, the magnetic twist breaks the symmetry of the phase speed of the MHD kink wave with respect to the propagation direction. As a result, in the presence of magnetic twist, the standing MHD kink oscillation with line-tying boundary conditions at the footpoints of the loop cannot be Fourier-analyzed in azimuthal and longitudinal directions. Ruderman (2015) called the modes corresponding with m=+1, -1 in a twisted flux tube, accelerated and decelerated kink wave, respectively. Ruderman & Terradas (2015) investigated the standing MHD kink oscillations of thin twisted magnetic tubes. They found that depending on the value of the plasma density ratio of the interior and exterior of the loop, the period ratio of the first overtone to the fundamental kink mode can be increased or decreased by increasing the magnetic twist in the loop. They also showed that in the presence of magnetic twist, in general, the eigenmodes of the MHD kink oscillations have elliptical polarization. Recently, in the thin tube thin boundary (TTTB) approximation, Ebrahimi & Karami (2016) analytically showed that the resonant absorption of kink MHD wave in a coronal flux tube with constant densities inside and outside the loop, can occur owing to the existence of a twisted magnetic field around the loop axis. They showed that when the ratio of the azimuthal to axial component of the background magnetic field increases, the frequency and the damping rate of the kink waves increase and the ratio of the frequency to the damping rate decreases. They found that with magnetic twist values in the range of observational values, the ratio of the damping time to the period of the oscillation is in good agreement with the observations.

Another consequence of existence of Alfvén continuum across the loop is phase-mixing of Alfvén waves, in which the oscillations of neighboring field lines become rapidly out of phase. This phenomenon leads to enhanced viscous and ohmic dissipations (Heyvaerts & Priest 1983).
Phase-mixing may occur either spatially in a propagating wave or in time in a standing wave. In both cases, an Alfvén wave is excited on each field line, which has an independent oscillation from its neighbors with a frequency in the Alfvén continuum.
Phase-mixing is an essential ingredient of resonant absorption (Poedts 2002) that causes a cascade of energy to small length scales, where the dissipation mechanisms become more efficient. By studying an initial-value problem we can gain insight into some of the interesting features that a nonuniform medium brings. However, we do not consider dissipation in our work and only are interested in the initial stage of the phase-mixing before the dissipation becomes important.

The main goal of the current paper is to explore the temporal evolution of kink waves in twisted flux tubes. To do so, we add an azimuthal component to the background magnetic field of the model of ST2015 and compare the results obtained in the presence of the magnetic twist with the results of ST2015 (no magnetic twist). To achieve this aim, in section <ref> we introduce the equations of motion and the flux tube model. In section <ref>, we solve an initial-value problem for the obtained equation of motion using the technique developed by ST2015. In section <ref>, we present numerical results. Finally, section <ref> is devoted to our conclusions.






§ EQUATIONS OF MOTION AND MODEL


The linearized ideal MHD equations for an incompressible plasma read

    ρ(r)∂^2 ξ/∂ t^2=-∇δ p+1/μ_0{(∇×δ𝐁)×𝐁 +(∇×𝐁)×δ𝐁},


    δ 𝐁=∇×(ξ×𝐁),


    ∇·ξ=0,

where ξ is the Lagrangian displacement of the
plasma, δ 𝐁 and δ p are the Eulerian perturbations of the
magnetic field and plasma pressure, respectively. Here μ_0 is the magnetic permeability
of the free space. Note that Eq. (<ref>) shows the incompressibility condition, which we adopt for the simplicity of calculations. Although the solar corona in general is a compressible medium, Goossens et al. (2009) elaborated that in the thin tube approximation that is applicable to the problem of long-wavelength transverse oscillations of coronal loops, kink waves are almost incompressible to a high degree of accuracy. They showed that the compressibility of the kink mode is proportional to (k_z R)^2. Hence, in the long-wavelength limit the frequency and damping rate of the kink mode are the same in both compressible and incompressible cases. The same result of Goossens et al. (2009) was previously obtained by Edwin & Roberts (1983), who explained that the kink mode behaves as an incompressible wave in the slender tube limit. Therefore, in order to apply the results of incompressible kink waves to the corona in what follows, we restrict our calculations to the limit of long-wavelength kink modes.

Using Eq. (<ref>), we can rewrite Eq. (<ref>) as

    δ 𝐁=(𝐁·∇)ξ-(ξ·∇)𝐁.

Putting Eq. (<ref>) into (<ref>) and doing some algebra yields

    ρ(r)∂^2 ξ/∂ t^2=-∇δ P+1/μ_0[(𝐁·∇)(𝐁·∇)ξ-(ξ·∇)(𝐁·∇)𝐁],

where δ P=δ p+(δ𝐁·𝐁)/μ_0 is the Eulerian perturbation of the total (gas plus magnetic) pressure.

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, φ, z) is
assumed to be as follows

    ρ(r)={[                                             ρ_ i,                                           r⩽ r_1,; ρ_ i/2[(1+ρ_ e/ρ_ i)-(1-ρ_ e/ρ_ i)sin(π/l(r-R))],                                      r_1< r< r_2,;                                             ρ_ e,                                           r⩾ r_2, ].

where r_1=R-l/2 and r_2=R+l/2. Here, l=r_2-r_1 is the characteristic length of
the radial variation of the background plasma density. The density ratio ρ_ i/ρ_ e is very difficult to estimate from observations. Typical values of this parameter are believed to be in the range ρ_i/ρ_e = 2-10 (Aschwanden et al. 2003). The background magnetic field is assumed to be constant and aligned with the flux tube axis in the regions r⩽ r_1 and r⩾ r_2, but it is twisted around the z-axis in the region r_1<r<r_2,


    𝐁(r)={[           B_0zẑ,          r⩽ r_1,; B_0φ(r)ϕ̂+B_0zẑ,     r_1< r< r_2,;           B_0zẑ,          r⩾ r_2,;                  ].

where B_0z is constant. We should note that in order to satisfy the magnetohydrostatic equation of motion,

    1/μ_0(∇×𝐁)×𝐁-∇ p=0,

the background magnetic field must be non-force-free, i.e. 
    (∇×𝐁)×𝐁=r̂B_0φ/r∂/∂ r(r B_0φ)≠ 0.
 This yields

    B_0φ≠C/r,

where C is a constant of integration. Therefore, in the model presented here, we can consider any profile for the azimuthal component of the background magnetic field other than B_0φ∝ r^-1.

Since the equilibrium quantities are only functions of r, the perturbations can be Fourier-analyzed with respect to the φ and z coordinates. Hence,

    δ P=δ P(r,t) e^i(mφ+k_z z),
    ξ=ξ(r,t) e^i(mφ+k_z
    z),

where m and k_z are the azimuthal and axial wavenumbers, respectively. Therefore, we are not considering the case of standing oscillations line-tied at the ends of the tube (see Terradas & Goossens 2012). Instead, we are implicitly considering propagating waves with fixed values of k_z and m. The case of line-lied oscillations is more difficult to tackle analytically, and a fully numerical approach is generally required. That is beyond the aim of the present work. Inserting perturbations (<ref>) into Eqs. (<ref>) and (<ref>) and eliminating ξ_φ and ξ_z, gives δ P in terms of ξ_r as

    δ P=1/k_z^2+m^2/r^2ℒ(-1/r∂(rξ_r)/∂ r)+m/r/k_z^2+m^2/r^2f(r)ξ_r,

where

    ℒ   ≡ρ(r)∂^2/∂ t^2+1/μ_0(m/rB_0φ+k_z B_0z)^2
       ≡ρ(r)(∂^2/∂ t^2+ω_A^2(r)),


    f(r)≡2/μ_0(mB_0φ^2/r^2+k_zB_0φ B_0z/r).

Here, the operator ℒ is the generalization of the Alfvén operator
ℒ_A=ρ(r)∂^2/∂ t^2+k_z^2
B_0z^2/μ_0, Eq. (9) in ST2015, in the presence of magnetic
twist and

    ω_A(r)≡1/√(μ_0 ρ(r))(m/rB_0φ(r)+k_zB_0z),

is the background Alfvén frequency.
Substituting Eq. (<ref>) in the φ component of Eq. (<ref>) and using Eqs. (<ref>) and (<ref>) one can relate ξ_φ to ξ_r as follows

    ℒξ_φ=i m/r/k_z^2+m^2/r^2ℒ(1/r∂(rξ_r)/∂ r)
    
        +2i/μ_0k_z^2/k_z^2+m^2/r^2(mB_0φ^2/r^2+k_z B_0φ B_0z/r)ξ_r.


Defining α≡ B_0φ(R)/B_0z as the twist parameter and ϵ≡ k_z R, one can obtain the orders of magnitude of the first and second terms on the right-hand side of Eq. (<ref>) denoted by T_0 and T_1, respectively, as

    T_0   ≡m/r/k^2+m^2/r^2ℒ(1/r∂(r ξ_r)/∂ r)
        ≃(B_0z^2/μ_0 R^2)(α+ϵ)^2/1+ϵ^2ξ_r,
       T_1   ≡k_z^2/μ_0/k^2+m^2/r^2(2mB_0φ^2/r^2+2k_zB_0φ B_0z/r)ξ_r
        ≃(B_0z^2/μ_0 R^2)ϵ^2α^2+αϵ/1+ϵ^2ξ_r.

Therefore,

    T_1/T_0∼ϵ^2α/α+ϵ∼{[ O(ϵ^2),    α≳ϵ,;  O(αϵ),    α≪ϵ.         ].

It is clear that in the limit of long-wavelength (ϵ≪1), we can ignore T_1 against T_0 in Eq. (<ref>). For instance, for α=0.01 and ϵ=0.03 we have T_1/T_0∼ 10^-4. Therefore, in the long wavelength limit (ϵ≪1), Eq. (<ref>) takes the form

    ξ_φ=i m/r/k_z^2+m^2/r^21/r∂ (rξ_r)/∂ r.

Putting Eq. (<ref>) into (<ref>) gives ξ_z in
terms of ξ_r as

    ξ_z=i k_z/k_z^2+m^2/r^21/r∂ (rξ_r)/∂ r.

Equations (<ref>) and (<ref>) show that for ϵ≪1, ξ_φ and ξ_z are not explicit functions of
the magnetic twist. However, the magnetic twist indirectly affects
ξ_φ and ξ_z by modifying the equation for ξ_r.

Eliminating δ P, ξ_φ and ξ_z from Eqs.
(<ref>) and (<ref>), we obtain the following differential
equation for ξ_r in the long wavelength limit (ϵ≪ 1)

    ℒℒ_sξ_r+[(k_z^2+m^2/r^2)∂ℒ/∂ r(1/r+∂/∂ r)+Φ(r)]ξ_r=0,

which is the generalized form of Eq. (16) in ST2015, in the presence
of a twisted magnetic field in a thin flux tube (ϵ≪ 1). Here, ℒ_s is the surface wave operator defined as

    ℒ_s≡(k_z^2+m^2/r^2)∂^2/∂ r^2+1/r(k_z^2+3m^2/r^2)∂/∂ r
        -1/r^2(k_z^2-m^2/r^2)-(k_z^2+m^2/r^2)^2,

and

    Φ(r)≡ g(r)(k_z^2+m^2/r^2)^2+2mk_z^2/r^2f(r)
        -m/r(k_z^2+m^2/r^2) df(r)/ d r,

where

    g(r)≡2/μ_0(-B_0φ^2/r^2+B_0φ/r d B_0φ/ d r).

Note that in the case of untwisted magnetic field (i.e. B_0φ=0), we have
f(r)=g(r)=Φ(r)=0 and ℒ=ℒ_A. In this case,
Eq. (<ref>) reduces to Eq. (16) in ST2015.

Here, we should note that in the case of twisted magnetic tubes, to avoid the kink instability  the twist value defined as ϕ_ twist=(L/R)(B_φ/B_z)=2π N_ twist must not exceed a critical value ϕ_ c (see e.g. Shafranov 1957; Kruskal et al. 1958; Hood & Priest 1979; Furno et al. 2006; Lapenta et al. 2006). Here, N_ twist is the number of twist turns in the tube and ϕ_ twsit is the angle of rotation (in radians) of twisted magnetic field per length L along the tube axis. According to the Kruskal-Shafranov analysis, for the particular case of a laboratory torus of major radius R_0, two points located an axial distance L = 2π R_0 apart refer to the same location on the torus, and k equals 2π/L. The kink instability is present in such a torus when ϕ_ twist>ϕ_c=2π. Hood & Priest (1979) considered the effect of line-tying at the ends of a flux tube and showed that for force-free magnetic fields of uniform twist, a magnetic twist larger than ϕ_c= 3.3π leads to the kink instability. Furno et al. (2006) investigated the kink instability in a flux tube that is tied at one end and free at the other end and showed that for ϕ_ twist>ϕ_c=π the flux tube is kink unstable (see also Lapenta 2006). Since we are interested in investigating the propagating waves, to avoid the kink instability, following Furno et al. (2006) we consider ϕ_c= π and take the length scale L=2π/k_z in the longitudinal direction of the tube. Therefore, the constraint

    ϕ_ twist=2π/k_z RB_φ/B_z=2πα/ϵ<ϕ_ c=π,

yields an upper limit for the twist parameter α_max=ϵ/2 in our model. However, to derive the real threshold for instability we need to investigate a stability analysis of the partial differential equations of motion which is beyond the scope of the present work.




§ SOLUTION


Solutions of ξ_r representing the surface kink waves in the
constant density and untwisted regions r<r_1 and r>r_2 in the
TT approximation (ϵ≪ 1) have been obtained by ST2015
as follows

    ξ_r  i(r,t)≈ A_ i(t),       r⩽ r_1,
          ξ_r  e(r,t)≈ A_ e(t)r^-2,  r⩾ r_2,

where A_ i(t) and A_ e(t) are the time-dependent amplitudes. In the nonuniform region r_1<r<r_2, following ST2015, we perform a modal expansion of the radial component of the Lagrangian displacement ξ_r(r,t) as

    ξ_r(r,t)=∑_n=1^∞a_n(t)ψ_n(r),

where the eigenfunctions ψ_n(r) satisfy the regular Sturm-Liouville system defined by the Bessel differential equation

    d^2ψ/ dr^2+1/r dψ/ dr+(λ^2-1/r^2)ψ=0.

Also, the functions ψ_n(r) have the following orthonormality relation

    1/l∫_r_1^r_2ψ_n(r)ψ_n'(r)r  dr=δ_nn'.

Following ST2015, using Eqs. (<ref>)-(<ref>) and applying the continuity of ξ_r and its derivative with respect to r at r=r_1 and r=r_2 one can obtain the boundary conditions governing ψ_n(r) as follows

    .d ψ/d r|_r=r_1=0,
    .(2/rψ+d ψ/d r)|_r=r_2=0.

The coefficient a_n(t) is computed by solving the following generalized eigenvalue problem which is obtained by inserting Eq. (<ref>) into (<ref>) (see Cally 1991 and ST2015)

    ℍ a = ω^2 𝕄 a.

Here ω^2 and a are the eigenvalue and the eigenvector, respectively, and the square matrices ℍ and 𝕄 are as follows

    H_nn'=1/l∫_r_1^r_2[1/μ_0(m/rB_0φ+k_z B_0z)^2ℒ_sψ_n'(r)+Φ(r)ψ_n'(r).
       .+2/μ_0(k_z^2+m^2/r^2)(m/rB_0φ+k_zB_0z)(-m/r^2B_0φ+m/r d B_0φ/ d r)(ψ_n'(r)/r+ dψ_n'(r)/ dr)]ψ_n(r) r  dr,
        
    
          M_nn'=1/l∫_r_1^r_2[ρ(r)ℒ_sψ_n'(r)+ dρ/ dr(k_z^2+m^2/r^2)
        (ψ_n'(r)/r+ dψ_n'(r)/ dr)]ψ_n(r) r  dr.

Following ST2015, the coefficients a_n(t) are obtained as

    a_n(t)=∑_n'=1^∞β_nn'[c_n'cos(ω_n't)+d_n'sin(ω_n't)],

where β_nn' is the nth component of the n'th eigenvector and ω_n' is the n'th eigenvalue. Also, the coefficients c_n and d_n are obtained with the help of suitable initial conditions.
From Eqs. (<ref>) and (<ref>), the expression for ξ_r(r,t) in the region r_1<r<r_2 takes the form

    ξ_r(r,t)=∑_n=1^∞∑_n'=1^∞β_nn'[c_n'cos(ω_n't)+d_n'sin(ω_n't)]ψ_n(r).

Equation (<ref>) can be recast in the following form

    ξ_r(r,t)=∑_n=1^∞[c_ncos(ω_nt)+d_nsin(ω_nt)]ϕ_n(r),

where

    ϕ_n(r)≡∑_n'=1^∞β_n'nψ_n'(r),

is the n's eigenfunction of Alfvén discrete modes (see Cally 1991 and ST2015). These modes are a discretized version of the Alfvén continuum. Hence, in the formalism of ST2015 the kink wave is not a global mode, but instead, it is built up as a superposition of Alfvén continuum modes. Like ST2015, we take the following initial conditions

    ξ_r(r,t=0)={[                           ξ_0,                        r⩽ r_1,;            ξ_0ψ_1(r)/ψ_1(r_1),                    r_1<r< r_2,; ξ_0ψ_1(r_2)/ψ_1(r_1)(r_2/r)^2,                        r⩾ r_2, ].
          ∂ξ_r/∂ t|_(r,t=0)=0,

where ξ_0 is a constant. Using Eqs. (<ref>), (<ref>),
(<ref>), (<ref>) and (<ref>) one can get
 
    c_n=β^-1_n,1ξ_0/ψ_1(r_1),
          d_n=0.

To solve Eq. (<ref>) numerically, we must truncate the infinite series of
Eq. (<ref>) to a finite number N of terms. This means that the Alfvén continuum is discretized in N different discrete modes. The bigger
the value of N the larger the evolution time that we are allowed to
proceed before the energy in the Nth Fourier mode becomes
significant and to the modal expansion starts to become inaccurate (for more details see Cally 1991). We shall consider a sufficiently large N to make sure that the number of terms in the modal expansion is enough for the considered duration of the temporal evolution.




§ NUMERICAL RESULTS 


In subsections <ref> and <ref>, we consider two types of twisted magnetic field to see how different twist profiles affect the phase-mixing of kink MHD waves in comparison with the results of ST2015. To do so, following ST2015, we set ρ_ i/ρ_ e=5, l/R=0.2,1 and k_zR=π/100. Time is in units of the period of the kink oscillation in a thin and untwisted loop, P_ kink=2π/ω_ kink, where

    ω_ kink=k_z√(ρ_ iv_A_i^2+ρ_ ev_A_e^2/ρ_ i+ρ_ e),

is the so-called kink frequency. Here v_A_i=B_0z/√(μ_0ρ_ i) and v_A_e=B_0z/√(μ_0ρ_ e) are the interior and exterior Alfvén speeds, respectively.




 §.§ Model I: discontinuous magnetic field


Following Ebrahimi & Karami (2016), we consider the azimuthal component of the magnetic field in the annulus region (i.e. r_1<r<r_2) as

    B_0φ(r)=Ar(r-r_1).

Note that from Eq. (<ref>), due to having a rotational discontinuity of the background
magnetic field at the location r=r_2, we have a delta-function
current sheet at r=r_2 in the axial direction. Note that in the presence of resistivity (which is absent in our model), the tearing mode instability can occur in this current sheet when
the driving force of the inflow exceeds the opposing Lorentz force. However, Ebrahimi & Karami (2016) showed that even in the presence of resistivity by choosing an appropriate thickness for the current sheet, tearing
mode instability can be avoided in the model (<ref>) during the kink oscillations.

The magnetohydrostatic equilibrium equation
takes the form

    d/ dr(p+B_0φ^2+B_0z^2/2μ_0)=-B_0φ^2/μ_0 r,

where p is the gas pressure. Using Eqs. (<ref>), (<ref>) and
continuity of the total (magnetic plus gas) pressure across r=r_1 and r=r_2, we obtain the gas pressure as

    p(r)={[                                                        p_0,                                                     r⩽ r_1,; ;           p_0-A^2/μ_0(3/4r^4-5/3r^3r_1+r^2r_1^2-1/12r_1^4),                                                 r_1<r< r_2,; ; p_0+A^2/μ_0(-1/4r_2^4-1/2r_1^2r_2^2+2/3r_1r_2^3+1/12r_1^4),                                                     r⩾ r_2, ].

where p_0 is a constant.

Here, we solve Eq. (<ref>) for α=0, 10^-4, 10^-2. The twist parameters considered here are not large enough to allow the flux tube to be kink unstable (see Eq. <ref>). The results for k_z R=π/100, m=±1, l/R=0.2 (thin layer) and l/R=1 (thick layer) are plotted in Figs. <ref>-<ref>. In the case of m=-1 and α=10^-2 we set N=300, but in other cases we set N=101. Figure <ref> shows the background Alfvén frequency, ω_A(r), and the corresponding discrete eigenfrequencies, ω_n, in the region r_1<r<r_2 for three values of α=0, 10^-4 and 10^-2. Note that the results for α=0 and 10^-4 overlap with each other. Here, we consider the small value of the twist parameter, α=10^-4, to show that how the results for the twisted magnetic field converge to the results of the untwisted magnetic field. The case with α = 0 corresponds to the model used by ST2015.




















Figure <ref> illustrates the normalized values of |c_n| versus their corresponding eigenfrequencies, ω_n. The figure reveals that for m=+1 and l/R=0.2, 1, by increasing the twist parameter, α, the peaks of the diagrams become wider and shift toward the larger frequencies. The situation for m=-1 is different. As illustrated in the bottom panels of Fig. <ref>, for m=-1, the peaks become narrower and shift to the smaller frequencies. Since |c_n| is the amplitude of the nth Alfvén discrete mode, the increase/decrease in the width of the frequency distribution implies that with increasing the twist parameter, the number of Alfvén discrete modes that have main contribution to the total displacements increases/decreases. As shown later, the wider frequency distribution results in an enhanced efficiency and a narrower frequency distribution results in a reduced efficiency of the phase-mixing process compared to the case with no twist.

Figures <ref> to <ref> display the temporal evolution of different components of the Lagrangian displacement for (i) m=+1 & l/R=0.2, (ii) m=+1 & l/R=1, (iii) m=-1 & l/R=0.2 and (iv) m=-1 & l/R=1, respectively. As illustrated in Figs. <ref> and <ref>, corresponding to positive m, in the presence of magnetic twist, at a given time, the perturbations are more phase-mixed than those in the case of untwisted magnetic field (α=0). Figures <ref> and <ref>, corresponding to negative m, show that at a given time, the perturbations are less phase-mixed than in the case of untwisted magnetic field. Thus, in model I for k_zR=π/100 and m=+1/-1, the small spatial scales due to phase-mixing in the nonuniform layer develop faster/slower than in the case of a straight field.

Figure <ref> shows the temporal evolution of ξ_r/ξ_0 at r=0 for m=±1 and l/R=0.2, 1. As shown in this figure, for m=+1/-1, in the presence of twisted magnetic field, ξ_r decays faster/slower than that of in the case of no twist for both l/R=0.2 and 1. Note that in Fig. <ref> for a given l/R, the frequency distribution of the Alfvén continuum modes for m=+1/-1 is wider/narrower for larger values of the twist parameter. This is consistent with the results of ST2015 who showed that when the distribution of the frequencies of the Alfvén continuum modes has a wider peak, there should be a larger damping rate.

Figure <ref> shows the power spectrum of ξ_r(r=0) in the time interval t∈[0,10 P_k] for α=0, 10^-4 and 10^-2. As illustrated in the figure, for l/R=0.2 and 1, by increasing the twist parameter, the peak frequency of the power spectrum increases for m=+1 and decreases for m=-1. This behavior is consistent with the results of Terradas & Goossens (2012) and Ruderman (2015) who showed that for k_z>0 and m=+1/-1, by increasing the magnetic twist in a coronal flux tube, the MHD kink frequency increases/decreases.

To illustrate the flux of the total (kinetic plus magnetic) energy from the internal and external regions to the inhomogeneous region, we calculate the integrated energy in each region as
 
    E_in=∫_0^r_11/2(ρ|∂ξ/∂ t|^2+1/μ|δ𝐁|^2)r  dr,
    
        E_nh=∫_r_1^r_21/2(ρ|∂ξ/∂ t|^2+1/μ|δ𝐁|^2)r  dr,
     E_ex=∫_r_2^∞1/2(ρ|∂ξ/∂ t|^2+1/μ|δ𝐁|^2)r  dr,

where δ𝐁 is calculated from Eq. (<ref>). Here, E_in, E_nh and E_ex are the integrated energies in the internal, nonhomogeneous and external regions, respectively. Note that in order to compute the third integral of Eq. (<ref>) we must replace the upper limit of the integral, ∞, with a sufficiently large radius (here r/R=20) where the amplitudes of the perturbations are approximately zero. Figure <ref> shows the integrated energy in these three regions as a function of time for α=0, 10^-4, 10^-2. It is clear from this figure that for both l/R=0.2, 1, when the twist parameter increases, for m=+1/-1 the rate of energy transfer from the internal and external regions to the inhomogeneous region increases/decreases. Hence, for m=+1/-1 the efficiency of the resonant absorption increases/decreases as the twist parameter becomes larger.






 §.§ Model II: continuous magnetic field


Following Terradas & Goossens (2012), in order to have a continuous magnetic field we consider a parabolic profile for the azimuthal component of the magnetic field (in the region r_1<r<r_2) as

    B_0φ(r)=A(r-r_1)(r_2-r).

The corresponding gas pressure can be obtained in a similar way to the one followed for the model (<ref>) as

    p(r)={[                                                                            p_0,                                                                         r⩽ r_1,;                                            p_0+A^2/2μ_0(3/2r^4-10/3(r_1+r_2)r^3;                                                    .+2(r_1^2+r_2^2+4r_1r_2)r^2.; .-6(r_1^2r_2+r_1r_2^2)r-1/6r_1^4+4/3r_2r_1^3+4r_1^2r_2^2+2r_1^2r_2^2ln(r/r_1)),                                                                      r_1<r<r_2,; ;                          p_0-A^2/2μ_0(1/6(r_2^4-r_1^4)+4/3(r_2r_1^3-r_1r_2^3)),                                                                         r⩾ r_2. ].





















Here, we solve Eq. (<ref>) for α=0, 10^-4, 10^-2, and N=101. Figures <ref>-<ref> show the results for m=±1 and l/R=0.2, 1. As illustrated in Fig. <ref>, the background Alfvén frequency for model II is not a monotonic function of r/R for both m=±1. As a result, for the ω_ns values that match the background Alfvén frequency at two positions, there should be two singularities in the corresponding ϕ_n(r) function. For instance, the dotted line in the bottom right panel of Fig. <ref> displays the value of ω_50 for m=-1, l/R=1 and α=10^-2 that matches the corresponding background Alfvén frequency at two points (big crosses). Figure <ref> shows that the ϕ_50(r) eigenfunction of the Alfvén discrete  mode becomes singular at two locations where the ω_50 matches the background Alfvén frequency.

Figure <ref> shows the normalized values of |c_n| versus their corresponding eigenfrequencies, ω_n. The figure presents that for the both l/R=0.2, 1, when m=+1, with increasing the twist parameter, the peaks of diagrams shift slightly toward larger frequencies and the frequency distribution of the Alfvén continuum modes becomes narrower. For m=-1 and α=10^-2 there are two peaks in the diagrams: a wide peak that is shifted to the larger frequencies and a narrow peak located at the beginning of the plots that is shifted to the smaller frequencies. Note that the frequency of the narrow peak corresponds to the minimum value of the Alfvén frequency in the nonuniform region (see the bottom panel of Fig. <ref>). Although the narrow peaks are not well resolved in the plots and are illustrated by only one value of c_n's (look at the single discontinuous red asterisk on the left side of the plots), the effect of this peak can be seen in the power spectrum of ξ_r(r=0) (see the bottom panel of Fig. <ref>). Note that the height of the narrow peak is smaller and larger than that of the wide peak for l/R=0.2, 1, respectively.

Figures <ref>-<ref> show the temporal evolution of different components of the Lagrangian displacement for (i) m=+1 & l/R=0.2, (ii) m=+1 & l/R=1, (iii) m=-1 & l/R=0.2 and (iv) m=-1 & l/R=1, respectively. As illustrated in these figures, in the presence of magnetic twist, at a given t/P_k, for m=+1/-1, the perturbations are less/more phase-mixed than those in the case of untwisted magnetic field (α=0).

Figure <ref> exhibits the temporal evolution of the radial component of the displacement on the loop axis for m=±1 and l/R=0.2, 1. The figure shows that for m=+1/-1, when the twist parameter increases, the decay rate of the perturbations decreases/increases. As illustrated in the bottom panel of Fig. <ref>, for m=-1, when the twist parameter increases there are two phases of oscillations for both l/R=0.2, 1. In the first phase the perturbations have slightly larger frequencies and decay faster than those in the case of untwisted magnetic field. In the second phase, the perturbations are almost decayless and have larger frequencies than those of in the case of untwisted magnetic field. These phases correspond to the wide and narrow peaks in the c_n's distribution (see the bottom panel of Fig. <ref>). Figure <ref> shows the power spectrum of ξ_r(r=0) for m=±1 and l/R=0.2, 1. It is clear in the figure that for m=-1 and l/R=0.2, 1, when the twist parameter increases, there are two frequencies corresponding to two oscillating phases in Fig. <ref> and two peaks of the c_n's distribution in Fig. <ref>. Figures <ref> and <ref> show that for m=+1 and l/R=0.2, 1, when the twist parameter increases, the loop axis oscillates with a frequency higher than the case with no twist.

Figure <ref> illustrates the integrated total (kinetic plus magnetic) energy calculated by Eq. (<ref>) in the internal, nonhomogeneous and external regions as a function of time for m=±1 and l/R=0.2, 1. It is clear from the figure that, for m=-1 and l/R=0.2 (bottom left panel), the magnetic twist enhances the rate of the energy flux toward the nonhomogeneous region but for m=-1 and l/R=1 (bottom right panel), there is not a big difference between the twisted and untwisted case except for the internal energy. As illustrated in the figure, for m=+1, the rate of energy transfer to the nonhomogeneous region decreases for both l/R=0.2 and 1.




 §.§ Comparison with the results of Terradas & Goossens (2012)


Terradas & Goossens (2012) investigated the effect of twisted magnetic field, Eq. (<ref>), on the MHD kink waves in a coronal flux tube. They solved the linearized MHD eigenvalue problem numerically and showed that in a magnetically twisted flux tube, the frequency of kink waves depends on the propagation direction. For instance, for a given value of the longitudinal wavenumber, k_z, the quasi-mode frequency of the kink waves has different values depending on the sign of the azimuthal mode number m=±1. Terradas & Goossens (2012) found that for m=1 and m=-1, the corresponding frequencies are larger and smaller than that in the case without twist, respectively. Now, we are interested in recovering this result in our work. To this aim, following Terradas & Goossens (2012), we take k_zR=π/50 and ρ_ i/ρ_ e=3. We also set r_1=0.5 and r_2=1.5 in Eq. (<ref>). In the model of Terradas & Goossens (2012), because of the existence of the discontinuous piecewise plasma density, the background Alfvén frequency is discontinuous at r=R. Hence, to avoid of this discontinuity, we replace and approximate the piecewise density profile, Eq. (1) in Terradas & Goossens (2012), with a sinusoidal density profile that has a thin transitional layer l/R=0.1.









In the following, we solve Eq. (<ref>) for m=±1, α=0, 0.01, and N=101 and compare our results with those obtained by Terradas & Goossens (2012). Figure <ref> shows the background Alfvén frequency and the corresponding eigenfrequencies of the Alfvén modes. Figure <ref> exhibits the temporal evolution of the radial component of the displacement on the loop axis. As illustrated in Fig. <ref>, the decay rate of the perturbations in the presence of magnetic twist for m=-1 and m=+1 is higher and lower than that of in the case without twist, respectively. Figure <ref> represents the power spectrum of ξ_r(r=0) in the time interval t∈[0,10P_k] for α=0, 0.01 and m=±1. As shown in this figure, in the presence of magnetic twist, the peak frequency of the power spectrum for m=1 and m=-1 is larger and smaller than that in the case of untwisted magnetic field, respectively. Note that the middle vertical dotted line represents the so-called kink frequency, Eq. (<ref>), obtained for an untwisted thin magnetic flux tube with a piecewise plasma density profile. The left and right vertical dotted lines represent the frequencies of the MHD kink waves obtained by Terradas & Goossens (2012) for m=-1 and m=+1, respectively, with the twist parameter α=0.01. Therefore, Fig. <ref> clearly shows that the result of Terradas & Goossens (2012) is completely recovered in our work.




§ CONCLUSIONS


Here, we investigated the effect of twisted magnetic field on the phase-mixing and resonant absorption of the propagating MHD kink waves in coronal flux tubes. The mathematical approach used in this paper is based on the work of Cally (1991) (in Cartesian coordinates) and ST2015 (in cylindrical coordinates). We solved an initial-value problem using the linear ideal MHD equations. Hence, our results cannot be extended to the large times when the strong phase-mixing develops in the system. The reason is that in this limit, due to the strong phase-mixing of the Alfvén waves, the viscous and resistive dissipation mechanisms become significant even though the dissipation coefficients are small in the corona (see, e.g. Heyvaerts & Priest 1983; Karami & Ebrahimi 2009).

Following ST2015, in order to find the temporal and spatial behaviour of the kink perturbations, using a modal expansion technique, we solved the linear incompressible MHD equations in the nonuniform region of a coronal flux tube that has both radial density variation and magnetic field twist. In order to simplify the MHD equations, we used the thin tube approximation (i.e. k_z R≪ 1) in our analysis. Thus, our results are only applicable to the propagating kink MHD waves in the limit of long wavelengths, i.e., when the wavelengths of the waves are much larger than the thickness of the loop.

We considered two types of twisted magnetic fields containing the discontinuous and continuous ones to investigate how different magnetic field profiles affect the resonant absorption and phase-mixing of the MHD kink waves in coronal loops. Also, we examined the effect of magnetic twist on the kink waves in the cases of both thin and thick nonuniform layers. In order to prevent the kink instability, the amount of the magnetic twist in the flux tube must be restricted with a maximum value that is obtained from a stability analysis. Investigating a stability analysis for the models considered here is beyond the scope of the present work. Instead, we considered the twist parameters small enough in order to be in the range of stability obtained in previous works.

One of the interesting effects of the twisted magnetic fields in coronal flux tubes is the asymmetry of the phase speed of the MHD kink waves with respect to the propagation direction (see, e.g., Terradas & Goossens 2012; Ruderman 2015). Hence, in order to investigate this effect, we considered two propagation directions (i) k_z R>0 & m=+1 and (ii) k_z R>0 & m=-1 for the kink MHD wave in both continuous and discontinuous magnetic field models.


For the model I (discontinuous magnetic field), we found the following:




  * By increasing the magnetic twist parameter in the loop, the decay rate of the radial component of the Lagrangian displacement on the axis of the flux tube increases/decreases for m=+1/-1.


  * The power spectrum of ξ_r(r=0) shows that by increasing the twist parameter, the effective frequency of the kink wave increases/decreases when m=+1/-1. Hence, in the presence of a twisted magnetic field, the frequency of the kink wave is asymmetric with respect to the propagation direction. This is in agreement with the result obtained by Terradas & Goossens (2012) and Ruderman (2015).


  * When the twist parameter increases, for both l/R=0.2 and 1, a wider/narrower range of the Alfvén continuum modes contributes to the total displacement of the kink waves for m=+1/-1.


  * The rate of phase-mixing of the perturbations increases/decreases for m=+1/-1 as the twist parameter increases in the loop. The reason is that when the twist parameter increases, for m=+1/-1 the slope of the profile of the Alfvén frequency increases/decreases.


  * The rate of energy flux from the interior and exterior regions of the loop toward the nonuniform region increases/decreases for m=+1/-1 when the twist parameter increases. As the energy of the kink wave transfers to the nonuniform region, the amplitude of perturbations inside and outside the loop decreases. The energy mostly transfers to the azimuthal component of the perturbations in the nonuniform region which is subjected to phase-mixing owing to the existence of an inhomogeneous background Alfvén frequency across the loop.


For the model II (continuous magnetic field), the results show the following:




  * With increasing the twist parameter, for m=+1 and l/R=0.2,1, a narrower range of the Alfvén continuum modes contributes to the total displacement of the kink waves. But for m=-1 and l/R=0.2,1 there are two peaks in the distribution profile of the Alfvén continuum modes. In comparison with the case of untwisted magnetic field, one of these peaks is wider and shifts to higher frequencies and another one is narrower and shifts to lower frequencies.


  * When the twist parameter increases, the decay rate of the radial component of the Lagrangian displacement on the axis of the flux tube decreases/increases for m=+1/-1.


  * For m=+1 and l/R=0.2,1, when the magnetic twist increases, the power spectrum of ξ_r(r=0) becomes narrower and shifts toward higher frequencies. When m=-1, for both l/R=0.2, 1, the power spectrum splits into two peaks: a wider peak in higher frequencies and a narrower peak in lower frequencies with respect to the case of untwisted magnetic field.


  * The rate of phase-mixing of the perturbations decreases/increases for m=+1/-1 as the twist parameter increases in the loop.


  * The rate of energy flux from the interior and exterior regions of the loop toward the nonuniform region decreases/increases for m=+1/-1 when the twist parameter increases.


In the case of the discontinuous magnetic field model, by increasing the twist parameter for the both cases of thin and thick nonuniform layers, the peak frequency of the power spectrum increases/decreases for m=+1/-1. This is in agreement with that obtained using the quasi-mode approach (see, e.g., Karami & Bahari 2010; Terradas & Goossens 2012; Ebrahimi & Karami 2016). However, for the continuous magnetic field model, the situation is different. In this case, when the twist parameter increases, for m=-1 the single peak of the power spectrum for l/R=0.2 splits into two peaks located at higher and lower frequencies than that of in the case of no twist. The reason is that the oscillation of the radial component of the Lagrangian displacement has two phases for m=-1. The first and second phases, respectively, have smaller and larger frequencies than those in the case of untwisted magnetic field. This result also holds in the case of thick nonuniform layer l/R=1. In this case, when the twist is absent, the power spectrum has two peaks. When the twist parameter increases, these peaks get away from each other and move toward higher and lower frequencies. These results point out the important effect that the particular twist model has on the behavior of kink waves. The ignorance of the actual twist profile in coronal loops turns out to be very important in this regard.

Applying the modal expansion approach to the model of Terradas & Goossens (2012), we found that when the twist parameter increases, the peak frequency of the power spectrum for m=1 and m=-1 shifts toward the higher and lower frequencies, respectively. This is in well agreement with the result obtained by Terradas & Goossens (2012).

As illustrated by Terradas & Goossens (2012), in the presence of magnetic field twist, quasi-mode frequencies of MHD kink waves obtained for k_z>0 & m=±1 are the same as for k_z<0 & m=∓1. To investigate this symmetry in our work, we obtained the results for k_z<0 & m=±1 in both models. We have not included these results here for the sake of simplicity. We found that, besides the effective frequency of the kink waves, the whole properties of the evolution of the propagating kink waves are symmetric if we change the signs of k_z and m, simultaneously. It is not straightforward to translate the present results for propagating waves to the case of standing oscillations line-tied at the ends of the tube. In the case of standing oscillations, the perturbations necessarily contain the two possible signs of k_z and m. Hence, in the case of standing waves it is not simple to deduce the net effect that the effect of twist would have on the process of phase-mixing. However, the present results suggest that the effect of twist can be relevant for standing waves as well. Further investigation in this direction is needed.

It is worth to mentioning that during kink oscillations of coronal loops the Kelvin-Helmholtz instability (KHI) can occur around the boundary of the flux tube (Heyvaerts & Priest 1983). The torsional motions, which are amplified in the inhomogeneous region of the flux tube, introduce velocity shears that are liable to be unstable to KHI. Since the observation of KHI has not been reported to date in coronal flux tubes, it is believed that some mechanism is able to suppress it. It is known that the existence of a component of a magnetic field aligned with the direction of the velocity shears has a stabilizing effect and can restrain the KHI (e.g. Chandrasekhar 1961). Soler et al. (2010) showed that a very small amount of magnetic twist, which is very likely and realistic in coronal flux tubes, can suppress the KHI in a cylindrical flux tube. Therefore, the lack of KHI can be one of the possible indirect confirmations of the existence of magnetic twist in coronal loops.

Twist of the magnetic field, even in a small amount, can have a significant impact on the generation of small scales and the energy cascade from the global kink motion to the small scales. This has implications concerning the efficiency of the phase-mixing process and the ability of the process to feed energy to the dissipative scales, where plasma heating takes place. In this paper we have considered a simple scenario to investigate the effect of twist, in order to pave the way for future works that should tackle the full nonlinear 3D problem.





§ ACKNOWLEDGEMENTS


The authors thank the anonymous referee for very valuable comments. R.S. acknowledges the support from grant AYA2014-54485-P (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.




























[Andries et al. (2005)]andries05 Andries, J., Goossens, M., Hollweg, J. V., Arregui, I., & Van Doorsselaere, T. 2005, , 430, 1109
[Aschwanden et al. (1999)]aschwanden99 Aschwanden, M. J., Fletcher, L., Schrijver, C. J., & Alexander D. 1999, , 520, 880
[Aschwanden et al. (2003)]aschwanden03 Aschwanden, M. J., Nightingale, R. W., Andries, J., et al. 2003, , 598, 1375
[Aschwanden et al. (2012)]aschwanden12 Aschwanden, M. J., Wuelser, J. P., Nariaki, V. N., Lemen, J. R., DeRosa, M. L., & Malanushenko, A. 2012, , 756, 124
[Bennett et al. (1999)]bennett99 Bennett, K., Roberts, B., & Narain, U. 1999, , 185, 41
[Cally (1991)]cally91 Cally, P. S. 1991, JPlPh, 45, 453
[Cally & Andries (2010)]cally10 Cally, P. S., & Andries, J. 2010, , 266, 17
[Carter & Erdélyi (2008)]carter08 Carter, B. K., & Erdélyi R. 2008, , 481, 239
[Chae et al. (2000)]chae00 Chae, J., Wang, H., Qiu, J., & Goode, P. R. 2000, , 533, 535
[Chae & Moon (2005)]chae05 Chae, J., & Moon, Y. J. 2005, , 629, 1110
[Chandrasekhar (1961)]chandrasekhar61 Chandrasekhar, S. 1961,Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon)
[Davila (1987)]davila87 Davila, J. M. 1987, , 317, 514
[Ebrahimi & Karami (2016)]ebrahimi16 Ebrahimi, Z., & Karami, K. 2016, , 462, 1002
[Edwin & Roberts (1983)]edwin83 Edwin, P. M., & Roberts, B. 1983, , 88, 179
[Erdélyi (1997)]erdelyi97 Erdélyi, R. 1997, , 171, 49
[Erdélyi & Carter (2006)]erdelyi06a Erdélyi, R., & Carter, B. K. 2006, , 455, 361
[Erdélyi & Fedun (2006)]erdelyi06b Erdélyi, R., & Fedun, V. 2006, , 238, 41
[Erdélyi & Fedun (2007)]erdelyi07 Erdélyi, R., & Fedun, V. 2007, ., 246, 101
[Erdélyi & Fedun (2010)]erdelyi10 Erdélyi, R., & Fedun, V. 2010, , 263, 63
[Furno et al. (2006)]Furno06 Furno, I., Intrator, T. P., Ryutov, D. D., Abbate, S., Madziwa-Nussinov, T., Light, A., Dorf, L., & Lapenta, G. 2006, PRL, 97, 015002
[Giagkiozis et al. (2016)]Giagkiozis16 Giagkiozis, I., Goossens, M., Verth, G., Fedun, V., & Van Doorsselaere, T. 2016, , 823, 71
[Goossens et al. (2002)]goossens02 Goossens, M., Andries, J., & Aschwanden, M. J. 2002, , 394, L39
[Goossens et al. (1992)]goossens92 Goossens, M., Hollweg, J. V., & Sakurai, T. 1992, , 138, 233
[Goossens & Ruderman (1995)]goossens95a Goossens, M., & Ruderman, M. S. 1995, PhSc, T60, 171
[Goossens et al. (1995)]goossens95b Goossens, M., Ruderman, M. S., & Hollweg, J. V. 1995, , 157, 75
[Goossens et al. (2009)]goossens09 Goossens, M., Terradas, J., Andries, J., Arregui, I., & Ballester, J. L. 2009, , 503, 213
[Goossens et al. (2011)]goossens11 Goossens, M., Erdélyi, R., & Ruderman, M. S. 2011, , 158, 289
[Goossens et al. (2013)]goossens13 Goossens, M., Van Doorsselaere, T., Soler, R., & Verth, G. 2013, , 768, 191
[Goossens et al. (2014)]goossens14 Goossens, M., Soler, R., Terradas, J., Van Doorsselaere, T., & Verth, G. 2014, , 788, 9
[Heyvaerts & Priest (1983)]heyvaerts83 Heyvaerts, J., & Priest, E. R. 1983, , 117, 220
[Hollweg (1984)]hollweg84 Hollweg, J. V. 1984, , 277, 392
[Hollweg (1987)]hollweg87 Hollweg, J. V. 1987, , 320, 875
[Hood & Pries (1979)]hood79 Hood A. W., & Priest E. R. 1979, , 64, 303
[Ionson (1978)]ionson78 Ionson, J. A. 1978, , 226, 650
[Jung & Van Doorsselaere (2016)]jung16 Jung Yu, D., & Van Doorsselaere, T. 2016, , 831, 30
[Karami & Bahari (2010)]Karami10 Karami, K., & Bahari, K. 2010, , 263, 87
[Karami & Bahari (2012)]Karami12 Karami, K., & Bahari, K. 2012, , 757, 186
[Karami & Barin (2009)]Karami09a Karami, K., & Barin, M. 2009, , 394, 521
[Karami & Ebrahimi (2009)]karami09b Karami, K., & Ebrahimi, Z. 2009, , 26, 448
[Karami et al. (2009)]karami09c Karami, K., Nasiri, S., & Amiri, S. 2009, , 394, 1973
[Kruskal et al. (1958)]kriskal58 Kruskal M. D., Johnson J. L., Gottlieb M. B., & Goldman L. M. 1958, Phys. Fluids, 1, 421
[Kwon & Chae (2008)]kwon08 Kwon, R. Y., & Chae, J. 2008, , 677, L141
[Lapenta et al. (2006)]lapenta06 Lapenta, G., Furno, I., Intrator, T., & Delzanno, G. L. 2006, JGR, 111, A12S06
[Lee & Roberts (1986)]lee86 Lee, M. A., & Roberts, B. 1986, , 301, 430
[Morton & Erdélyi (2009)]morton09 Morton, R. J., & Erdélyi, R. 2009, , 707, 750
[Nakariakov et al. (1999)]nakariakov99 Nakariakov, V. M., Ofman, L., DeLuca, E. E., Roberts, B., & Davila, J. M. 1999, Sci., 285, 862
[Ofman et al. (1994)]ofman94 Ofman, L., Davila, J. M., & Steinolfson, R. S. 1994, , 421, 360
[Ofman (2005)]ofman05 Ofman, L. 2005, Adv. Space Res., 36, 1572
[Ofman (2009)]ofman09 Ofman, L. 2009, , 694, 502
[Poedts (2002)]poedts02 Poedts, S. 2002, in IAU Coll. 188, SOLMAG 2002, Proc. Magnetic Coupling of the Solar Atmosphere Euroconference, ed. H. Sawaya-Lacoste (Noordwijk: ESA), 273
[Ruderman & Roberts (2002)]ruderman02 Ruderman, M. S., & Roberts, B. 2002, , 577, 475
[Ruderman (2007)]ruderman07 Ruderman, M. S. 2007, , 246, 119
[Ruderman (2015)]ruderman15a Ruderman, M. S. 2015, , 575, A130
[Ruderman & Terradas (2015)]ruderman15b Ruderman, M. S., & Terradas, J. 2015, , 580, A57
[Safari et al. (2006)]safari06 Safari, H., Nasiri, S., Karami, K., & Sobouti, Y. 2006, , 448, 375
[Sakurai et al. (1991a)]sakurai91a Sakurai, T., Goossens, M., & Hollweg J. V. 1991a, , 133, 227
[Sakurai et al. (1991b) ()]sakurai91b Sakurai, T., Goossens, M., & Hollweg J. V. 1991b, , 133, 247
[Soler et al. (2010)]soler10 Soler, R., Terradas, J., Oliver, R., Ballester, J. L., & Goossens, M. 2010, , 712, 875
[Soler et al. (2013)]Soler13 Soler, R., Goossens, M., Terradas, J., & Oliver, R. 2013, , 777, 158
[Soler & Terradas (2015)]soler15 Soler, R., & Terradas, J. 2015, , 803, 43
[Shafranov (1957)]shafranov57 Shafranov V. D. 1957, J. Nucl. Energy II, 5, 86
[Thalmann et al. (2014)]thalmann14 Thalmann, J. K., Tiwari, S. K., & Wiegelmann, T. 2014, , 780, 102
[Terradas & Goossens (2012)]terradas12 Terradas, J., & Goossens, M. 2012, , 548, A112
[Terradas et al. (2006)]terradas16 Terradas, J., Oliver, R., & Ballester, J. L. 2006, , 642, 533
[Van Doorsselaere et al. (2004)]VanDoorsselaere04 Van Doorsselaere, T., Andries, J., Poedts, S., & Goossens, M. 2004, , 606, 1223
[Wang et al. (2015)]Wang2015 Wang, H., Cao, W., Liu, C., Xu, Y., Liu, R., Zeng, Z., Chae, J., Ji, H. 2015, Nature Communications, 6, 70086












