Issue 
A&A
Volume 665, September 2022



Article Number  A113  
Number of page(s)  12  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202244175  
Published online  15 September 2022 
Transition to turbulence in nonuniform coronal loops driven by torsional Alfvén waves
II. Extended analysis and effect of magnetic twist^{⋆}
^{1}
Departament de Física, Universitat de les Illes Balears, 07122 Palma de Mallorca, Spain
email: s.diaz@uib.es
^{2}
Institute of Applied Computing and Community Code (IAC3), Universitat de les Illes Balears, 07122 Palma de Mallorca, Spain
Received:
2
June
2022
Accepted:
8
July
2022
It has been shown in a previous work that torsional Alfvén waves can drive turbulence in nonuniform coronal loops with a purely axial magnetic field. Here we explore the role of the magnetic twist. We modeled a coronal loop as a transversely nonuniform straight flux tube, anchored in the photosphere, and embedded in a uniform coronal environment. We considered that the magnetic field is twisted and control the strength of magnetic twist by a free parameter of the model. We excited the longitudinally fundamental mode of standing torsional Alfvén waves, whose temporal evolution was obtained by means of highresolution threedimensional ideal magnetohydrodynamic numerical simulations. We find that phase mixing of torsional Alfvén waves creates velocity shear in the direction perpendicular to the magnetic field lines. The velocity shear eventually triggers the Kelvin–Helmholtz instability (KHi). In weakly twisted magnetic tubes, the KHi is able to grow nonlinearly, and subsequently, turbulence is driven in the coronal loop in a similar manner as in the untwisted case. When the magnetic twist remains weak, it delays the onset of the KHi and slows the development of turbulence down. In contrast, magnetic tension can suppress the nonlinear growth of the KHi when the magnetic twist is strong enough, even when the KHi has locally been excited by the phasemixing shear. Thus, turbulence is not generated in strongly twisted loops.
Key words: magnetohydrodynamics (MHD) / methods: numerical / Sun: atmosphere / Sun: oscillations / waves
Movie associated to Fig. 5 is available at https://www.aanda.org
© S. DíazSuárez and R. Soler 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
Highresolution and highcadence observations have shown the ubiquity of magnetohydrodynamic (MHD) waves throughout the solar atmosphere (see, e.g., Aschwanden et al. 1999; Nakariakov et al. 1999; Tomczyk et al. 2007; De Pontieu et al. 2007; McIntosh et al. 2011; Morton et al. 2015; Jafarzadeh et al. 2017; Srivastava & Dwivedi 2017). The dissipation of MHD waves could play an important role in the heating of the solar corona (see, e.g., Hollweg 1978; Cranmer & van Ballegooijen 2005; Cargill & de Moortel 2011; Mathioudakis et al. 2013; Soler et al. 2019; Van Doorsselaere et al. 2020; Nakariakov & Kolotkov 2020), and the acceleration of the solar wind (see, e.g., Charbonneau & MacGregor 1995; Cranmer 2009; Matsumoto & Suzuki 2012; Shoda et al. 2018).
Torsional Alfvén waves are a subtype of axisymmetric Alfvén waves in cylindrical flux tubes. They are nearly incompressible, the restoring force is the magnetic tension, and their direction of polarization is perpendicular to the magnetic field lines. This type of waves produces axisymmetric perturbations in the perpendicular components of the velocity and magnetic field. Furthermore, unlike Alfvén waves with other azimuthal symmetries, torsional Alfvén waves do not couple with magnetoacoustic waves when the magnetic field is straight (Goossens et al. 2011). Torsional Alfvén waves have been reported in bright points (Jess et al. 2009) and in a coronal active region structure (Kohutova et al. 2020). The torsional motions found in the chromosphere and transition region by De Pontieu et al. (2012, 2014) can also be interpreted as torsional Alfvén waves. Aschwanden & Wang (2020) interpreted oscillations in the magnetic energy during solar flares as torsional Alfvén waves. This type of waves has also been reported in the solar wind during the interaction between coronal mass ejections (Raghav & Kule 2018).
In coronal loops, which are closed magnetic structures with two footpoints anchored in the solar photosphere, there is no report of such waves. Nonetheless, torsional Alfvén waves can be excited at their footpoints through vortex or twisting plasma motions in the solar photosphere (see, e.g., Fedun et al. 2011; Shelyag et al. 2011, 2012; WedemeyerBöhm et al. 2012; Mumford et al. 2015; Srivastava et al. 2017). Recently, Soler et al. (2021) have investigated the propagation of torsional Alfvén waves from the photosphere to a coronal loop. As the waves reach the coronal loop, they can resonate with standing modes of the loop and drive global torsional oscillations (see also Hollweg 1984b,a). Soler et al. (2021) found that a large amount of energy can be transmitted at coronal heights despite the chromospheric filtering. Particularly, they reported that the energy flux is channeled mostly through the fundamental standing torsional mode of the loop.
Inspired by these results, in DíazSuárez & Soler (2021b; hereafter Paper I), we investigated the nonlinear evolution of torsional Alfvén waves in a straight coronal flux tube with a constant axial magnetic field. In Paper I (see also the torsional model in Guo et al. 2019), we showed that the phase mixing of the torsional Alfvén waves generates azimuthal shear flows. These flows eventually excite the Kelvin–Helmholtz instability (KHi), as Heyvaerts & Priest (1983) and Browning & Priest (1984) predicted analytically. Phase mixing is a linear phenomenon that occurs when waves in adjacent radial positions oscillate with different frequencies. The cause of phase mixing is a frequency continuum, whose origin lies in the inhomogeneities in density and/or magnetic field (see, e.g., Heyvaerts & Priest 1983; Nocera et al. 1984; De Moortel et al. 2000; Smith et al. 2007; Prokopyszyn et al. 2019). Phase mixing increases the values of vorticity and current density and generates small scales, although at a relatively slow pace.
The KHi triggered by phase mixing can evolve nonlinearly by forming large eddies that break into smaller eddies. This process initiates and drives turbulence. Turbulence accelerates the energy cascade to small scales, which might enhance the efficiency of wave energy dissipation (Hillier et al. 2020). There is evidence that coronal loops may be in a turbulent state (De Moortel et al. 2014; Hahn & Savin 2014; Liu et al. 2014; Xie et al. 2017). The generation of turbulence is also a result obtained by numerous studies of numerical simulations of kink oscillations of coronal loops (see, e.g., Terradas et al. 2008, 2018; Antolin et al. 2015; Magyar & Van Doorsselaere 2016; Howson et al. 2017a; Karampelas et al. 2019; Antolin & Van Doorsselaere 2019).
We here extend the investigation of Paper I by replacing the constant axial magnetic field by a twisted magnetic field. The existence of twisted coronal loops has been confirmed by observations (see, e.g., Kwon & Chae 2008; Aschwanden et al. 2012; Thalmann et al. 2014; Aschwanden 2019). The behavior of linear MHD waves in twisted flux tubes has extensively been investigated analytically or semianalytically (see, e.g., Bennett et al. 1999; Erdélyi & Fedun 2006, 2007, 2010; Ruderman 2007; Karami & Barin 2009; Karami & Bahari 2010; Terradas & Goossens 2012; Ebrahimi et al. 2017; Ebrahimi & Soler 2022). The nonlinear evolution of kink MHD waves in twisted flux tubes has also been investigated numerically (Howson et al. 2017a; Terradas et al. 2018). Here, we explore how magnetic twist affects the nonlinear evolution of torsional Alfvén waves.
2. Numerical setup
2.1. Initial configuration
As in Paper I, we use the standard coronal loop model (see, e.g., Edwin & Roberts 1983). The model is made of an overdense cylindrical flux tube of radius, R, and length, L, embedded in a lowβ uniform coronal environment. We set L/R = 10. We used a shorter loop length than what is reported from observations, typically L/R ∼ 100. The main reason for this is that we wish to speed up the simulation times. The periods of the standing torsional oscillations are proportional to the loop length. Thus, if a longer loop were considered, the periods and simulation times would be longer, but the dynamics would be essentially the same. The effect of considering larger loop lengths was explored in Paper I.
The loop footpoints were fixed at two rigid walls representing the solar photosphere. We neglected the curvature of the coronal loop and the thin chromospheric layer at the loop feet. In Paper I we assumed a uniform magnetic field along the cylinder axis. Here, we improved the model by considering a twisted magnetic field, which is a more realistic representation of the magnetic field in coronal loops, namely
where B_{φ} and B_{z} are the azimuthal and longitudinal components of the background magnetic field in a cylindrical coordinate system denoted by r, φ, and z. In our model, the zaxis coincides with the flux tube axis. We considered the same forcefree twist model as was used in Terradas et al. (2018), namely
where B_{0} is a constant that corresponds to the magnetic field strength at the tube axis, and c is a dimensionless parameter that controls the amount of magnetic twist.
The highest value of the magnetic twist parameter used here is c = 0.4, which was chosen according to observational constraints. Aschwanden (2013) computed the twist angle, θ, in coronal loops as , where q_{free} is the free energy ratio. The calculation of q_{free} can be achieved from observations through lineofsight magnetograms, the threedimensional (3D) position of the coronal loop, and a nonlinear forcefree field code. Table 3 in Aschwanden et al. (2014, 2016) show that the free energy ratio ranges from 0 to ∼0.25. Therefore, twist angles smaller than ∼25° are expected. The twist angle of our model can be calculated as
where is the modulus of the magnetic field. For c = 0.4, the twist angle is θ = 21.8° at r = R, which is slightly smaller than the largest angle inferred by Aschwanden et al. (2014, 2016). On the other hand, the number of turns of the magnetic field lines over the cylinder length can be computed as
where the factor LB_{φ}(r)/rB_{z}(r) ≡ Φ(r) is the absolute twist. Kwon & Chae (2008) considered a sample of 14 loops from Transition Region And Coronal Explorer (TRACE; Handy et al. 1999) observations in the extreme ultraviolet and determined that their absolute twist values ranged from 0.22π to 1.73π. This result implies that 0.11 < N < 0.865. At r = R in our model, where B_{φ} is maximum, we obtain N = 0.64 for c = 0.4 and L/R = 10, so that the maximum twist considered here is below the upper limit of the interval inferred by Kwon & Chae (2008).
A twisted magnetic field is prone to develop a kink instability. In the case of uniform twist, Hood & Priest (1979) found that the kink instability can develop when the absolute twist is larger than 3.3π. With an improved analysis, Hood & Priest (1981) reduced the threshold value to 2.49π. Using the definition of N in Eq. (5), we find that the kink instability may appear in our model if N ≥ 1.65 according to the critical twist of Hood & Priest (1979) and if N ≥ 1.245 according to the critical twist of Hood & Priest (1981). Clearly, we are sufficiently far away from the critical twist for the kink instability because our maximum twist (N = 0.64) is still relatively weak.
We note that we used a preexisting twisted magnetic field in our model. A different way to proceed would be to create twist in the model by dynamically rotating the magnetic field lines anchored in the photosphere, as done in Ofman (2009). This alternative approach imposes a more complex scenario and is not considered here for simplicity.
Since the magnetic field is force free, we considered in the model a background uniform gas pressure, p_{0}. In turn, the equilibrium density, ρ_{0}, was the same as in Paper I. Thus, we used
In Eq. (6), ρ_{i} is the density in the uniform inner core, ρ_{e} is the uniform external density, and ρ_{tr}(r) is the density in a nonuniform transitional layer of thickness l that continuously connects the two uniform regions as
where l may range from 0 (abrupt transition) to 2R (fully nonuniform loop). We set ρ_{i}/ρ_{e} = 2 and l = 0.4R. A scheme of our model is depicted in Fig. 1.
Fig. 1. Sketch of the coronal fluxtube model. The isosurface of density with value equal to (ρ_{i} + ρ_{e})/2 is shown in orange. Magnetic field lines located at r = 2R, r = R and r = 0 are drawn in green, blue, and red, respectively, for the model with c = 0.4. The bottom and top gray planes represent the solar photosphere. 
The equilibrium Alfvén, v_{A}(r), and sound, v_{s}(r), speeds are
where γ is the adiabatic constant and μ_{0} is the vacuum magnetic permeability. Figure 2 shows the radial profiles of both speeds for a particular set of parameters. The Alfvén speed and sound speed both vary inside the transition region because the density is nonuniform. However, when c ≠ 0, the Alfvén speed is also position dependent in the inner core and in the external medium, where the density is uniform. This happens because the magnetic field is nonuniform everywhere when twist is present.
Fig. 2. Equilibrium radial profiles of the Alfvén speed, v_{A}(r), and sound speed, v_{s}(r), for the twist parameter ranging from c = 0 to c = 0.4. The sound speed is the same in all cases. The region between the vertical dashed lines denotes the inhomogeneous region in density. The values are normalized with respect to the Alfvén speed at r → ∞. We used ρ_{i}/ρ_{e} = 2 and l/R = 0.4. 
According to Halberstadt & Goedbloed (1993), the linetied Alfvén frequency continuum is
where n is the harmonic number. A continuous spectrum is present because of the radial variation of the zcomponent of the magnetic field, B_{z}(r), and the density, ρ_{0}(r). The density is only nonuniform in the transition region. However, B_{z}(r) is also nonuniform in the inner core of the loop when twist is present. Because of twist, the Alfvén continuum also extends into the uniform core of the loop, as we show in Fig. 3.
Fig. 3. Radial profile of the linetied Alfvén frequency continuum for the twist parameter ranging from c = 0 to c = 0.4. The region between the vertical dashed green lines corresponds to the nonuniform region in density. The values are normalized to the external Alfvén frequency. We used ρ_{i}/ρ_{e} = 2 and l/R = 0.4. 
Throughout the paper, length and density are normalized with respect to the radius of the flux tube, R, and the external density, ρ_{e}, respectively. Velocities are normalized with respect to the external Alfvén speed at r → ∞, namely v_{A, ∞} = v_{A}(r→∞). The normalized time is , with τ_{A} = R/v_{A, ∞}.
2.2. Numerical code
As in Paper I, the numerical simulations were performed with the PLUTO code (Mignone et al. 2007). We solved the 3D ideal MHD equations with a finitevolume shockcapturing spatial discretization on a structured mesh. The equations are as follows:
In Eqs. (11)–(14), denotes the total derivative, ρ is the mass density, v is the velocity, B is the magnetic field, and p is the gas pressure. Gravity and nonideal terms are neglected.
The PLUTO code solves Eqs. (11)–(14) in Cartesian coordinates. As in Paper I, the code uses a RoeRiemann solver (Roe 1981) to compute the numerical fluxes, a piecewise parabolic scheme for spatial reconstruction, and a secondorder RungeKutta algorithm for temporal evolution. In order to maintain the solenoidal constraint, the generalized Lagrange multiplier (GLM) method (Dedner et al. 2002) is used. The background magnetic field is not currentfree because of magnetic twist, and therefore the backgroundfieldsplitting technique cannot be used in the present simulations, in contrast to Paper I. For comparison purposes, the backgroundfieldsplitting technique is not used for the straight magnetic field case either. The code therefore advances the full magnetic field vector here.
We used the same base computational box as in Paper I, namely a numerical resolution of 100 × 100 × 100 cells distributed uniformly from −3R to 3R in the x and ydirections, and from −L/2 to L/2 in the zdirection. The code implements the adaptive mesh refinement (AMR) strategy (Mignone et al. 2012) by which the cells split into two, and so the spatial resolution doubles if a certain criterion is satisfied. The refinement criterion is based on the second derivative error norm (Lohner 1987) of a parameter closely related with the kinetic energy density because the goal of the AMR strategy in the present simulations is to conserve wave energy as much as possible. The use of AMR decreases the computational cost by keeping a low resolution where the dynamics is not relevant and generating a high resolution in the regions of interest. We included four levels of refinement, so that the maximum effective resolution was 1600 × 1600 × 1600. For typical values of a coronal loop radius, the effective transverse resolution was ∼10 km.
Although the AMR scheme allowed us to deal with extremely fine scales, the dynamics developed in the flux tube eventually cascaded the energy to scales below the maximum effective resolution. As in Paper I, to determine the simulation time at which this unphysical energy loss starts to happen, we monitored the kinetic, magnetic, and internal energies. In this way, we were able to determine the time at which numerical dissipation is beginning to play an important role. This time also coincides with the time when the total integrated vorticity saturates and decreases (see later).
We set the same set of boundary conditions as in Paper I. We considered outflow conditions, that is, zero gradient for pressure, density, and the x and ycomponents of the magnetic field at all boundaries. We also set the outflow condition for the zcomponent of the magnetic field, B_{z}, at the lateral boundaries. To mimic the linetying of the magnetic field lines at the solar photosphere, B_{z} was fixed to the equilibrium value through Eq. (3) at the top and bottom boundaries. The three velocity components vanish at all boundaries to avoid energy injection from outside the domain.
2.3. Initial perturbation
In the solar atmosphere, torsional Alfvén waves can be driven by photospheric plasma motions (see, e.g., Fedun et al. 2011; Shelyag et al. 2011, 2012; WedemeyerBöhm et al. 2012; Mumford et al. 2015; Srivastava et al. 2017). The waves propagate to the corona through the chromosphere, where reflection and dissipation heavily reduces the upward wave transmission (see Soler et al. 2019). When a broadband spectrum of waves driven at the photosphere is transmitted into a closed coronal structure such as a coronal loop, Soler et al. (2021) showed that most of the energy is deposited in the longitudinally fundamental mode. In view of this result and as in Paper I, we excited the longitudinally fundamental mode of standing torsional Alfvén waves at t = 0. We did this by perturbing the components of velocity. To prevent the perturbation from also initially exciting slow MHD waves, we imposed an initial velocity field in the same manner as in DíazSuárez & Soler (2021a), so that the velocity was purely perpendicular to the magnetic field lines, namely
where v_{⊥} is the perpendicular velocity to the magnetic field lines. We considered the same form as that used in Paper I, namely
where v_{0} is the maximum velocity amplitude and A(r) contains the radial profile (see Paper I). We set v_{0} = 0.1v_{A, ∞}. The radial profile of v_{φ} and v_{z} at z = 0 is shown in Fig. 4 for different values of the twist parameter, c.
Fig. 4. Radial profiles of the azimuthal (top) and longitudinal (bottom) velocity components at t = 0 for the twist parameter ranging from c = 0 to c = 0.4. The inset in the top panel shows a closeup view of the azimuthal velocity around its maximum. Both velocity components are normalized to the maximum amplitude of the perpendicular velocity perturbation, v_{0}. 
We recover the same initial condition as in Paper I for the untwisted case, c = 0. When we increased the magnetic twist, the maximum of v_{φ} slightly decreased, but its profile reamined mostly unaltered. The longitudinal component of velocity, v_{z}, increased significantly when the twist parameter grew. This is a direct effect of the relative variation of the azimuthal and longitudinal components of the magnetic field.
Although the initial perturbation will mostly excite torsional Alfvén waves, other types of MHD can also appear in the flux tube. The inclusion of a twisted magnetic field introduces linear coupling between torsional Alfvén waves and fast magnetoacoustic sausage waves (Sakurai et al. 1991; Goossens et al. 1992, 2011; Giagkiozis et al. 2015). The waves excited by the initial perturbation when c ≠ 0 will develop mixed properties of both torsional Alfvén waves and fast magnetoacoustic sausage waves. Nevertheless, because the twist is weak, the torsional Alfvén wave character will be dominant. Later in the evolution, the ponderomotive force will additionally nonlinearly drive slow MHD waves (Hollweg 1971; Rankin et al. 1994; Tikhonchuk et al. 1995; Terradas & Ofman 2004).
3. Revisiting the straight magnetic field case
The main goal of this paper is to illustrate the effects of a twisted magnetic field in the nonlinear evolution of the standing torsional Alfvén waves. In particular, we investigate its influence on the triggering of the KHi and the associated turbulence generation. However, before we investigate the twist, we revisit the straight magnetic field case of Paper I and perform an extended study.
We reran the simulation corresponding to the thinlayer case of Paper I but for longer times and deeper into the turbulent phase. The simulation was continued beyond the capacity of the AMR scheme in this way to correctly describe the smallest spatial scales that develop in the computational box. Although the dynamics of the smallest scales will then be affected by numerical dissipation, we can assume that the evolution at spatial scales larger than the maximum effective resolution will remain valid. Having this in mind, we can study how turbulence evolves globally at longer times despite the unreliable information at the smallest scales.
We refer to Paper I for detailed explanations of the simulation evolution. We give a brief summary. Standing torsional Alfvén waves excited in the flux tube develop the process of phase mixing (see, e.g., Heyvaerts & Priest 1983; Nocera et al. 1984; De Moortel et al. 2000; Smith et al. 2007; Prokopyszyn et al. 2019; Ebrahimi et al. 2020). In the transition region, where density is nonuniform, waves in adjacent magnetic surfaces oscillate with a slightly different frequency (see Eq. (10)) and become out of phase as time increases. The effect of phase mixing is to generate an azimuthal velocity shear, which eventually triggers the KHi (see, e.g., Heyvaerts & Priest 1983; Soler et al. 2010; Zaqarashvili et al. 2015; Barbulescu et al. 2019). Initially, the KHi develops locally in the transition region, which means that it does not globally perturb the flux tube. Nonetheless, the KHi evolves nonlinearly as time increases. The nonlinear evolution of the KHi breaks large eddies into smaller eddies. The instability naturally drives turbulence in the flux tube, which develops perpendicularly to the background magnetic field. Turbulence extends beyond the boundaries of the initial nonuniform transition, and further accelerates the energy cascade toward small scales. Because of turbulence, the internal and external plasmas mix. Thus, the dynamics of torsional Alfvén oscillations is similar to that of linearly polarized kink (see, e.g., Terradas et al. 2008; Antolin et al. 2015; Magyar & Van Doorsselaere 2016; Karampelas et al. 2019; Antolin & Van Doorsselaere 2019; Pascoe et al. 2020; MartínezGómez et al. 2022) or circularly polarized kink oscillations (Magyar et al. 2022). Nonetheless, the flux tube is not displaced laterally in the case of torsional oscillations, as occurs with the excitation of kink oscillations, and the resonant absorption process is absent for torsional Alfvén waves (see, e.g., Goossens et al. 2011).
As shown in Paper I, a parameter that helps us understand the different phases of the evolution is the vorticity. We investigated the vorticity squared, ω^{2}(r, t) = ∇×v(r,t)^{2}, and its value integrated over the whole computational domain, Ω^{2}, namely
Figure 5 shows a crosssectional cut at the tube center, z = 0, of ω^{2}(r,t) in logarithmic scale. Only a subdomain of the complete numerical domain is shown, where the relevant dynamics occurs. The inset to the right shows the integrated vorticity squared, Ω^{2}(t), as a function of the computational time. The blue dot tracks the position on the curve of Ω^{2} as time varies. The complete temporal evolution is available as a movie.
Fig. 5. Crosssectional cut at z = 0 of the logarithm with base 10 of the vorticity squared, ω^{2}(r,t), in arbitrary units when . The small panel to the right shows the temporal evolution of ω^{2}(r,t) integrated over the whole computational domain, Ω^{2}(t), normalized with respect to the value at . The blue dot marks the position in the curve for the considered simulation time. The complete temporal evolution is available online as a movie. 
In the small panel of Fig. 5, the temporal evolution of Ω^{2} undergoes three phases (see Paper I). In the first phase for , Ω^{2} oscillates and slightly increases owing to the linear development of phase mixing. In this phase, the spatial distribution of ω^{2}(r,t) is mostly in the form of concentric rings. The concentric rings deform after the KHi is locally excited. In the second phase for , the KHi has grown enough to behave nonlinearly, and it induces turbulence. The KHi vortices break into increasingly smaller vortices. As a result, there is a dramatic increase in vorticity. The spatial distribution of ω^{2}(r,t) clearly reveals the formation of very fine scales owing to turbulence. Finally, in the third phase for , the AMR scheme can no longer fully describe the smallest scales that appear in the evolution. Thus, vorticity decreases unphysically because of numerical dissipation. As a reference, the periods of the internal and external Alfvén oscillations are and 20 normalized time units, respectively.
The particular snapshot displayed in the large panel of Fig. 5 shows the spatial distribution of ω^{2}(r,t) at the time for which Ω^{2} is maximum. The loop boundary is clearly fully turbulent, and extremely fine structures in vorticity have already been generated. Several authors have studied the evolution of vorticity in numerical simulations where the kink mode is excited (see, e.g., Terradas et al. 2008; Antolin et al. 2015; Howson et al. 2017b, 2020, 2021; Karampelas et al. 2017; Karampelas & Van Doorsselaere 2018; Guo et al. 2019), and some of their results are comparable with ours. Similar vorticity structures as those displayed in Fig. 5 can also be seen in previouss works of nonlinear kink oscillations (see, e.g., Antolin et al. 2014; Howson et al. 2017a; Antolin & Van Doorsselaere 2019). Nonetheless, the spatial distribution is different due to the different azimuthal symmetry of the torsional and kink oscillations.
The turbulent evolution of the flux tube is clearly highlighted by the formation of extremely fine vorticity structures captured by the high resolution of the AMR scheme. A consequence of turbulence is the mixing of the internal and external plasmas, which heavily modifies the transverse density profile adopted initially. To illustrate this effect, Fig. 6 displays an azimuthal average of the transverse density profile at z = 0 for various times. Figure 6 is similar to the density panels from Fig. 5 of Karampelas & Van Doorsselaere (2018) for nonlinear kink waves. On average, we find that the width of the nonuniform region increases with time. For long times in our simulation, part of the initially uniform core becomes turbulent, and thus nonuniform. Consequently, at sufficiently long times beyond those simulated here, the whole flux tube should become turbulent. It is likely that for kink waves, which are global oscillations of the whole tube, the fully turbulent regime is reached faster than for torsional Alfvén waves, which are local waves in nature. This would require further analysis. We note that the density variations are smoothed because of the azimuthal average. Along specific directions, we can find density variations that are far stronger than those displayed in Fig. 6. In addition, the full temporal evolution of the results displayed in Fig. 6 also shows periodic density increases in the loop core caused by the ponderomotive force (see, e.g., Hollweg 1971).
Fig. 6. Azimuthally averaged density at the tube center, z = 0, for different simulation times. The vertical dashed purple lines show the limits of the nonuniform region at . Density is normalized to the external density. 
4. Exploring the effect of magnetic twist
Here we investigate the influence of magnetic twist. In addition to the simulation without magnetic twist (c = 0) discussed before, we performed four additional simulations including magnetic twist with c = 0.1, 0.2, 0.3, and 0.4.
Before the onset of the KHi, when the dynamics is governed by the development of phase mixing, the results of all simulations are quite similar, regardless of the strength of magnetic twist. The reason is that the shape of the Alfvén frequency continuum in the nonuniform layer is very similar in all cases (see Fig. 3). Therefore, we focus on later times. In Fig. 7 we show crosssectional cuts of the density at the tube center, z = 0, at three different simulation times corresponding to already advanced stages of the simulations. The results for c = 0 are shown in the first row. The twist parameter increases from c = 0.1 to c = 0.4 from top to bottom in the remaining rows. The results for c = 0.3 are not shown as they are visually indistinguishable from those for c = 0.4. Figure 7 suggests that the inclusion of the magnetic twist in our model delays the onset of KHi. For the same time, the KHi vortices appear to be progressively less developed when the twist increases.
Fig. 7. Crosssectional cuts of density at the tube center, z = 0, for different simulation times: (left column), (mid column), and (right column). The first row corresponds to the simulation without magnetic twist, c = 0. In descending order, the remaining rows correspond to the simulations with magnetic twist parameters of c = 0.1, 0.2, and 0.4. Density is normalized to the external density. 
On the other hand, when magnetic twist is included, the KHi vortices are not strictly perpendicular to the flux tube axis and can also develop in the zdirection (Terradas et al. 2018). To illustrate this fact, we show in Fig. 8 a longitudinal density cut of the simulation with c = 0.2. This cut was made for the same time as in the third column of Fig. 7.
Fig. 8. Longitudinal density cut at y = 0 for the simulation with magnetic twist parameter, c = 0.2 when , i.e., the same time as in the right column of Fig. 7. Density is normalized to the external density. 
4.1. Delay or suppression of the KHi
We aim to quantify how magnetic twist postpone the KHi onset time. Theoretically, a component of the magnetic field along the direction of the shear flow has a stabilizing effect on the KHi (see, e.g., Chandrasekhar 1961). Significant effort has been made to understand this problem for the KHi driven by transverse MHD waves. However, obtaining the onset time of the KHi in a cylindrical tube with timevarying flows under the presence of a twisted magnetic field remains a challenge analytically. Browning & Priest (1984) analytically obtained an expression for the onset time of the KHi for timevarying flows in Cartesian geometry and with a straight field. Soler et al. (2010) considered the KHi at the boundary of a cylindrical tube, but with constant azimuthal flows. Soler et al. (2010) studied the effect of magnetic twist by considering a Cartesian analog with an inclined magnetic field and found that twist can decrease the growth rates and even suppress the instability. Zaqarashvili et al. (2015) studied the instability of a cylindrical, twisted, and rotating jet. They found that jets are unstable to the KHi only when the kinetic rotation energy is higher than the magnetic energy of the twist. Hillier et al. (2019) investigated the linear stability of a discontinuous oscillatory shear flow in Cartesian geometry in the presence of a uniform magnetic field. They showed that parametric instabilities can also grow in addition to the KHi. Barbulescu et al. (2019) also considered a Cartesian model of timevarying flows, but including an inclined magnetic field on one side of the interface to mimic the effect of twist, as was previously done by Soler et al. (2010) for constant flows. Barbulescu et al. (2019) concluded that the magnetic shear may reduce the instability growth rate, but it cannot completely stabilize the interface, in contrast to the constant flow case of Soler et al. (2010).
To understand the effect of twist on the triggering of the KHi, we followed the method introduced by Terradas et al. (2018). Their method is based on studying the excitation of different azimuthal wave numbers. This approach was also used in Antolin & Van Doorsselaere (2019) and in Paper I. We proceeded as follows. In a crosssectional cut at the tube center, z = 0, we calculated the velocity component perpendicular to the magnetic field lines, v_{⊥}, in the middle of the transition region, r = R, as a function of the azimuthal angle, φ, from φ = 0 to φ = 2π. After this, we computed the discrete Fourier transform to the data using the fast Fourier transform (FFT) algorithm with the Scipy module (Virtanen et al. 2020). Following the notation of Terradas et al. (2018), the discrete Fourier transform is
where N is the number of sample points, g(k) is the angular sampling of v_{⊥}, and p is an integer that plays the role of the azimuthal wavenumber. Unlike in Paper I, we considered positive and negative values of p, namely p = 0, ±1, ±2, …, ±(N − 1), since twist breaks the degeneracy in the sign of the azimuthal wavenumber. p = 0 is the torsional plus the sausage mode. p = ±1 are the kink modes, and p≥2 are the fluting modes (see, e.g., Roberts 2019). In the presence of magnetic twist, the torsional mode and the fast sausage mode are coupled linearly (see, e.g., Goossens et al. 2011). However, the contribution of the fast sausage mode is weak because the magnetic twist considered here is relatively weak, even for c = 0.4. Thus, for p = 0, the torsional mode dominates the fast sausage mode.
Figure 9 shows the results of the azimuthal Fourier analysis. To better visualize the results, we plot the temporal evolution of the p = 0 mode alone in the top panel. In the bottom panel, we plot the temporal evolution of the sum of the first 20 positive modes and the first 20 negative modes, excluding the p = 0 mode. We checked that adding more modes to the sum does not affect the results.
Fig. 9. Temporal evolution of the Fourier coefficients at the tube center, z = 0, and, in the middle of the inhomogeneous radial layer, r = R. Top: the Fourier coefficient associated with the torsional mode, p = 0. Bottom: same as the top panel, but for the sum of the first 20 positive and the first 20 negative Fourier coefficients excluding p = 0. In both panels, the twist parameter ranges from c = 0 to c = 0.4. Arbitrary units are used. 
As expected, we obtain that the p = 0 mode is dominant during the linear evolution of phase mixing regardless of the magnetic twist. The dominance of the p = 0 mode lasts until the onset of KHi, which excites higher azimuthal wave numbers. In this moment, the sum of the p ≠ 0 Fourier coefficients begins to display significant variations. These results are consistent with those of Terradas et al. (2018) and Antolin & Van Doorsselaere (2019) for the standing kink mode and also agree with those of Paper I. During the development of KHi and the subsequent turbulence, the initially low amplitudes of the sum of the p ≠ 0 Fourier coefficients become comparable with or even larger than the amplitude of the p = 0 mode. Moreover, we find that the p = 0 mode loses its periodicity.
Similarly to the results of Paper I, regardless the presence or absence of magnetic twist, we verified that the excitation of highorder azimuthal wave numbers initially occurs in the nonuniform transition region. This analysis is not shown here for simplicity. We also verified the persistence of the p = 0 mode in the core where the KHi cannot develop. However, when the turbulence expands into the core, highorder azimuthal wave numbers are also excited in the core. This is consistent with the results from Fig. 6 for the untwisted case.
Therefore, the azimuthal Fourier analysis confirms what was already visualized in Fig. 7: The inclusion of the magnetic twist delays the development of the KHi. If the magnetic twist is sufficiently weak, the KHi is delayed, but the dynamics is similar to that in the absence of twist, although it evolves at a slower pace. In a sense, this conclusion is not completely new, although this is the first time that the role of twist in the nonlinear evolution of standing torsional Alfvén waves has been investigated. Similar results can be found in simulations of nonlinear kink waves in twisted flux tubes (see, e.g., Howson et al. 2017a; Terradas et al. 2018). For the different twist profiles considered by Howson et al. (2017a), the authors reported that the KHi vortices grow less as the magnetic twist is increased. Howson et al. (2017a) also reported that the location of the maximum twist may affect at the development of KHi vortices. Terradas et al. (2018) also studied nonlinear kink waves considering the same twist profile as we used here. Similarly to Howson et al. (2017a), Terradas et al. (2018) reported that the KHi vortices grow less when the magnetic twist increases.
When the magnetic twist is strong enough, however, the physical scenario can be different. The simulations with strong magnetic twist suggest that the KHi can be suppressed altogether even if the phasemixinggenerated shear flows should remain unstable according to linear theory (see Barbulescu et al. 2019). Although the KHi itself can linearly be triggered in a local fashion, its perturbations are not allowed to fully grow into the nonlinear regime when the tension force provided by the twist is strong enough (DíazSuárez & Soler 2021a). The important role of magnetic tension in nonlinearly preventing the growth of the KHi vortices has been discussed by Galinsky & Sonnerup (1994), Ryu et al. (2000), and Hillier (2019), for instance. In our model, the nonlinear suppression of the KHi seems to occur for c ≳ 0.3.
Because the KHi suppression for strong twist is a nonlinear phenomenon, we cannot exclude that the initial amplitude of the excited torsional wave may affect the subsequent stability. For sufficiently large amplitudes, the KHi is expected to break the stabilizing effect of the magnetic tension even if the twist is strong. We therefore speculate that the ability of the KHi to grow might depend on the relation between the kinetic energy of the shear flow that drives the KHi and the magnetic energy of the twisted magnetic field, in a similar fashion as in the model by Zaqarashvili et al. (2015). The effect of the initial amplitude is not investigated here and is left for forthcoming works.
4.2. Turbulence and vorticity generation
The delay or suppression of the KHi has important consequences in the dynamics of the flux tube. In our model, the generation of turbulence, which results from the KHi evolution, can just be delayed or be absent, depending on the strength of magnetic twist. In addition, magnetic twist may affect the subsequent evolution of turbulence. To explore this, we resort again to the vorticity to study how magnetic twist affects the evolution of turbulence.
Figure 10 shows crosssectional cuts of ω^{2}(r,t) at the tube center, z = 0, for the same three times as in Fig. 7. In each row, the strength of the magnetic twist is different and sorted in increasing order from top to bottom. For comparison purposes, we also include the results from the simulation without magnetic twist in the first row. Unlike in Fig. 7, we find slightly different vorticity patterns for the cases with a magnetic twist of c = 0.3 and c = 0.4. This is highlighted below. Figure 10 shows that magnetic twist not only plays a role in delaying or suppressing the KHi, but also heavily affects the value and spatial distribution of vorticity when turbulence is present. Another plot that helps us in this discussion is Fig. 11. There, we show the temporal evolution of the vorticity squared integrated in the whole computational domain, Ω^{2}, as a function of the computational time for all simulations. The curve for c = 0 plotted in Fig. 11 is the same as was shown in the small panel of Fig. 5.
Fig. 10. Crosssectional cuts at z = 0 of the logarithm with base 10 of the vorticity squared, ω^{2}(r,t), for different simulation times: (left column), (middle column), and (right column). The first row corresponds to the simulation without magnetic twist, c = 0. In descending order, the remaining rows correspond to the simulations with a magnetic twist parameter of c = 0.1, 0.2, 0,3, and 0.4. Arbitrary units are used. 
We start by focusing on comparing the simulations with c = 0 (no twist) and c = 0.1 (weakest magnetic twist considered). Figure 11 indicates that the onset of the KHi occurs practically at the same time in both simulations because the sharp increase in vorticity in the two cases is almost simultaneous. Thus, the delaying effect of twist when c = 0.1 is still not very relevant. However, slightly lower values of integrated vorticity are found when c = 0.1 compared with those of c = 0. The first two rows in Fig. 10 show that the spatial distribution of vorticity when turbulence sets in is very similar in the two cases. Nevertheless, smaller vorticity structures appear when magnetic twist is absent. This effect is clearly visible in the two simulations when , a time for which turbulence is already well established.
Fig. 11. Temporal evolution of the vorticity squared integrated in the whole computational domain, Ω^{2}, for the twist parameter ranging from c = 0 to c = 0.4. The curve for c = 0 (dotteddashed red line) is the same as in the small panel of Fig. 5. All curves are normalized to their values at t = 0. 
The effect of magnetic twist becomes more pronounced in the simulation with c = 0.2. Figure 11 shows that the KHi onset is significantly delayed with respect to the cases with c = 0 and c = 0.1 discussed before. Figure 10 plainly shows that the dynamics advances at a slower pace when c = 0.2. For instance, the spatial distribution of ω^{2}(r,t) at is similar to that found at in the simulation with c = 0.1. The case with c = 0.2 corresponds to an intermediate situation in the sense that twist already has an important impact in delaying the KHi onset and in the appearance and evolution of turbulence. However, magnetic twist is still not strong enough to avoid the KHi growth and so prevent the generation of turbulence.
The simulations with stronger magnetic twist, c = 0.3 and c = 0.4, represent a completely different scenario in which turbulence is not driven because the KHi is not allowed to grow. Even for the longest considered simulation times, the shape of ω^{2}(r,t) in the two simulations is essentially in the form of concentric rings, which is the typical structure caused by phase mixing. However, when the shape of ω^{2}(r,t) in the case with c = 0.3 reveals that some of the phasemixingdriven concentric rings are slightly deformed in a wavy fashion. No such deformation is seen in the simulation with c = 0.4 at that time. This deformation in the vorticity rings when c = 0.3 is caused by the very local triggering of the KHi within the nonuniform transition layer. This local inception of the KHi is not appreciable as discernible deformations in density. Nonetheless, vorticity, which is computed from the velocity field, is much more sensitive to the very initial appearance of the KHi.
We ran the simulation with c = 0.3 up to longer times than those displayed in Fig. 10 and confirmed that the KHi does not grow and no vortices form. As explained before, the reason is the nonlinear suppresion of the KHi by magnetic tension (see, e.g., Galinsky & Sonnerup 1994; Ryu et al. 2000; Hillier 2019; DíazSuárez & Soler 2021a). Thus, turbulence is not driven. In the simulation with c = 0.4, the local inception of the KHi is eventually observed at later times, but the perturbations are not allowed to grow and turbulence is equally absent. However, Fig. 11 shows that vorticity continues to increase for both c = 0.3 and c = 0.4 at long times. The reason is that the linear phase mixing remains at work for long times and causes the increase in vorticity. Nevertheless, the increase in vorticity owing to phase mixing alone occurs at a much slower pace than when the KHi is present.
Again, our results about the effect of magnetic twist on the vorticity are comparable with those obtained from nonlinear kink wave simulations (Howson et al. 2017a; Terradas et al. 2018). Howson et al. (2017a) showed crosssectional cuts of the magnitude of vorticity at the loop apex for a twisted and an untwisted case at two different simulation times. In both cases, they reported large and small vortices due to KHi in their untwisted case, but their twisted case did not show small vortices. Thus, the magnetic twist suppresses the small vortices associated with KHi. Howson et al. (2017a) also found that the increase in vorticity due to KHi occurs earlier when twist is absent.
4.3. Kinetic energy cascade to small scales
In Paper I we showed that the onset of turbulence in the flux tube accelerates the cascade of energy from large to small spatial scales. This process is slower before the appearance of turbulence, when phase mixing is the only working mechanism. The delay or even suppression of the turbulent regime when magnetic twist is included should affect the rate at which small scales are generated.
We investigated how kinetic energy is transported to small scales depending on the strength of twist. For this purpose, we only considered the kinetic energy associated with motions perpendicular to the magnetic field lines, namely . We excluded the contribution from the other components of velocity because we are mainly interested in the Alfvénic part of the energy. The presence of twist linearly couples torsional Alfvén waves and fast magnetoacoustic sausage waves (see, e.g., Sakurai et al. 1991). Fast magnetoacoustic sausage waves would predominantly perturb the radial velocity component. In turn, Alfvén waves and slow magnetoacoustic waves are nonlinearly coupled through the ponderomotive force (see, e.g., Hollweg 1971). Slow waves would mainly affect the longitudinal velocity component.
We considered a crosssectional cut at the tube center, z = 0, and averaged the Alfvénic kinetic energy in the azimuthal direction using 512 radial cuts. Thus, the angular resolution is approximately ∼0.7°. We denote the averaged energy by , which is a function of r and t. Then, we applied the continuous Fourier transform in the radial direction, discretized due to the limited numerical resolution. Following Paper I, it can be defined as
where N is the number of samples, k is the radial wave number, Δr is the spatial resolution, and r_{0} is the upper limit of the radial domain. The summatory in Eq. (21) is exactly the 1D discrete Fourier transform, which we computed using the FFT module of Scipy (Virtanen et al. 2020) based on the FFT algorithm (Cooley & Tukey 1965). Finally, we calculated the modulus of and normalized it with respect to the maximum value in the spectrum at each time.
Figure 12 displays the results of the Fourier transform for all simulations at the last common simulation time for which the vorticity still increases in all simulations, . The reason for choosing this time is that for longer times, the simulations with c = 0 and c = 0.1 are already significantly affected by numerical dissipation (see the decrease in integrated vorticity in Fig. 11), which may also affect the results of the Fourier transform at small scales.
Fig. 12. Azimuthally averaged amplitude spectrum of the Alfvénic (perpendicular to the magnetic field) part of the kinetic energy for the twist parameter ranging from c = 0 to c = 0.4. We considered the last common time for which vorticity increases in all simulations, . Arbitrary units are used. 
At large scales, kR < 10^{0}, the amplitude spectrum is independent on the wave number. For intermediate scales, 10^{0} < kR < 10^{2}, we do not find that the values of the amplitude spectrum are ordered in any particular way regarding the twist parameter, c. The temporal evolution of the amplitude spectrum (not shown here) shows oscillations by which the various curves with different c exchange their relative positions depending on the considered time. The reason for this behavior is unclear to us. Understanding the behavior of the amplitude spectrum in these intermediate scales requires additional research.
We are more interested in the spectrum at small scales, 10^{2} < kR < 10^{3}, where the amplitude spectrum decreases as the wavenumber increases. In the absence of twist, the values of the amplitude spectrum are higher than those obtained when twist is considered. In addition, we find a decreasing trend with respect to the twist parameter, c, so that the stronger the magnetic twist, the lower the values of the amplitude spectrum. These results again confirm the relevant role of magnetic twist in suppressing the generation of small scales and in the rate at which Alfvénic energy is deposited into those small scales. We intentionally avoided kR > 10^{3} in our analysis owing to the possible role of numerical dissipation at the smallest scales.
5. Concluding remarks
We concluded the investigation started in Paper I into the ability of torsional Alfvén waves to drive turbulence in coronal loops. Because of plasma and/or magnetic nonuniformity, standing torsional Alfvén waves undergo phase mixing, generating shear flows perpendicular to the magnetic field direction in adjacent radial positions as time increases. Eventually, the shear flows trigger the KHi, as Heyvaerts & Priest (1983) and Browning & Priest (1984) first predicted.
In the absence of magnetic twist (the case studied in Paper I), the KHi can grow nonlinearly without opposition. The nonlinear evolution of the KHi naturally induces turbulence. First, large KHi eddies are formed, which later break into smaller eddies in a cascadelike dynamics. Mixing of plasma occurs. An important increase in vorticity is found, and the generation of small scales speeds up compared with the initial phase, which is dominated by phase mixing alone. Turbulence evolves perpendicularly to the magnetic field, so that for a straight magnetic field, turbulence is pseudo2D.
We have explored the role of magnetic twist. While the linear evolution of phase mixing is similar in the presence or absence of twist, the dynamics of the subsequent KHi growth and turbulence generation strongly depends on the strength of the twist.
If magnetic twist is sufficiently weak, the onset and growth of the KHi is delayed, but the dynamics is similar to that in the case of a straight magnetic field. In this regime, the stronger the magnetic twist, the longer the delay. Although the vorticity still increases dramatically during the development of turbulence, the increase is smaller than in absence of magnetic twist. Small scales are still quickly generated by turbulence, although at a slower pace than when twist is absent. Turbulence still evolves perpendicularly to the magnetic field lines, but because the field is twisted, turbulence is no longer confined to perpendicular planes to the tube axis.
Conversely, when the magnetic twist is strong enough, the scenario is completely different. If the strength of magnetic twist surpasses a critical value or a critical twist angle, no KHi vortices can grow. The KHi itself is still locally excited by the phasemixinggenerated shear flows, but the tension of the twisted magnetic field now prevents the nonlinear development of the instability (see, e.g., Galinsky & Sonnerup 1994; Ryu et al. 2000; Hillier 2019; DíazSuárez & Soler 2021a). As a consequence, there is no enhancement of vorticity, nor is turbulence driven. Thus, the generation of small scales is not accelerated and continues to evolve slowly at the rate dictated by phase mixing.
Although the effect of magnetic twist may oppose the process, the nonlinear evolution of torsional Alfvén waves remains a viable mechanism to induce turbulence as long as coronal loops are weakly twisted. Other effects should be explored in the future. For instance, the tension of the background field in curved coronal loops may also affect the triggering and evolution of the KHi and associated turbulence. It might be also interesting to investigate turbulence generation in multithreaded loops (see, e.g., Ofman 2009).
Movie
Movie 1 associated with Fig. 5 (fig_5_vortyc0) Access here
Acknowledgments
This publication is part of the R+D+i project PID2020112791GBI00, financed by MCIN/AEI/10.13039/501100011033. S.D. acknowledges the financial support from MCIN/AEI/10.13039/501100011033 and European Social Funds for the predoctoral FPI fellowship PRE2018084223. We also thank the UIB for the use of the Foner cluster. For the simulation data analysis we have used VisIT (Childs et al. 2012) and Python 3.6. In particular, we have used Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020), and Numpy (Harris et al. 2020). We are thankful to B. Vaidya and his contributors for the tool pyPLUTO.
References
 Antolin, P., & Van Doorsselaere, T. 2019, Front. Phys., 7, 85 [Google Scholar]
 Antolin, P., Yokoyama, T., & Van Doorsselaere, T. 2014, ApJ, 787, L22 [Google Scholar]
 Antolin, P., Okamoto, T. J., De Pontieu, B., et al. 2015, ApJ, 809, 72 [Google Scholar]
 Aschwanden, M. J. 2013, Sol. Phys., 287, 369 [NASA ADS] [CrossRef] [Google Scholar]
 Aschwanden, M. J. 2019, ApJ, 874, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Aschwanden, M. J., & Wang, T. 2020, ApJ, 891, 99 [Google Scholar]
 Aschwanden, M. J., Fletcher, L., Schrijver, C. J., & Alexander, D. 1999, ApJ, 520, 880 [Google Scholar]
 Aschwanden, M. J., Wuelser, J.P., Nitta, N. V., et al. 2012, ApJ, 756, 124 [NASA ADS] [CrossRef] [Google Scholar]
 Aschwanden, M. J., Xu, Y., & Jing, J. 2014, ApJ, 797, 50 [Google Scholar]
 Aschwanden, M. J., Reardon, K., & Jess, D. B. 2016, ApJ, 826, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Barbulescu, M., Ruderman, M. S., Van Doorsselaere, T., & Erdélyi, R. 2019, ApJ, 870, 108 [Google Scholar]
 Bennett, K., Roberts, B., & Narain, U. 1999, Sol. Phys., 185, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Browning, P. K., & Priest, E. R. 1984, A&A, 131, 283 [NASA ADS] [Google Scholar]
 Cargill, P., & de Moortel, I. 2011, Nature, 475, 463 [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon) [Google Scholar]
 Charbonneau, P., & MacGregor, K. B. 1995, ApJ, 454, 901 [Google Scholar]
 Childs, H., Brugger, E., Whitlock, B., et al. 2012, in High Performance VisualizationEnabling ExtremeScale Scientific Insight, 357 [Google Scholar]
 Cooley, J. W., & Tukey, J. W. 1965, Math. Comput., 19, 297 [Google Scholar]
 Cranmer, S. R. 2009, Liv. Rev. Sol. Phys., 6, 3 [Google Scholar]
 Cranmer, S. R., & van Ballegooijen, A. A. 2005, ApJS, 156, 265 [Google Scholar]
 De Moortel, I., Hood, A. W., & Arber, T. D. 2000, A&A, 354, 334 [NASA ADS] [Google Scholar]
 De Moortel, I., McIntosh, S. W., Threlfall, J., Bethge, C., & Liu, J. 2014, ApJ, 782, L34 [Google Scholar]
 De Pontieu, B., McIntosh, S. W., Carlsson, M., et al. 2007, Science, 318, 1574 [Google Scholar]
 De Pontieu, B., Carlsson, M., Rouppe van der Voort, L. H. M., et al. 2012, ApJ, 752, L12 [Google Scholar]
 De Pontieu, B., Rouppe van der Voort, L., McIntosh, S. W., et al. 2014, Science, 346, 1255732 [Google Scholar]
 Dedner, A., Kemm, F., Kröner, D., et al. 2002, J. Comput. Phys., 175, 645 [Google Scholar]
 DíazSuárez, S., & Soler, R. 2021a, ApJ, 922, L26 [CrossRef] [Google Scholar]
 DíazSuárez, S., & Soler, R. 2021b, A&A, 648, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ebrahimi, Z., & Soler, R. 2022, MNRAS, 511, 3477 [NASA ADS] [CrossRef] [Google Scholar]
 Ebrahimi, Z., Karami, K., & Soler, R. 2017, ApJ, 845, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Ebrahimi, Z., Soler, R., & Karami, K. 2020, ApJ, 893, 157 [Google Scholar]
 Edwin, P. M., & Roberts, B. 1983, Sol. Phys., 88, 179 [Google Scholar]
 Erdélyi, R., & Fedun, V. 2006, Sol. Phys., 238, 41 [CrossRef] [Google Scholar]
 Erdélyi, R., & Fedun, V. 2007, Sol. Phys., 246, 101 [CrossRef] [Google Scholar]
 Erdélyi, R., & Fedun, V. 2010, Sol. Phys., 263, 63 [CrossRef] [Google Scholar]
 Fedun, V., Shelyag, S., Verth, G., Mathioudakis, M., & Erdélyi, R. 2011, Ann. Geophys., 29, 1029 [Google Scholar]
 Galinsky, V. L., & Sonnerup, B. U. Ö. 1994, Geophys. Res. Lett., 21, 2247 [NASA ADS] [CrossRef] [Google Scholar]
 Giagkiozis, I., Fedun, V., Erdélyi, R., & Verth, G. 2015, ApJ, 810, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Goossens, M., Hollweg, J. V., & Sakurai, T. 1992, Sol. Phys., 138, 233 [Google Scholar]
 Goossens, M., Erdélyi, R., & Ruderman, M. S. 2011, Space Sci. Rev., 158, 289 [Google Scholar]
 Guo, M., Van Doorsselaere, T., Karampelas, K., et al. 2019, ApJ, 870, 55 [Google Scholar]
 Hahn, M., & Savin, D. W. 2014, ApJ, 795, 111 [Google Scholar]
 Halberstadt, G., & Goedbloed, J. P. 1993, A&A, 280, 647 [NASA ADS] [Google Scholar]
 Handy, B. N., Acton, L. W., Kankelborg, C. C., et al. 1999, Sol. Phys., 187, 229 [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Heyvaerts, J., & Priest, E. R. 1983, A&A, 117, 220 [NASA ADS] [Google Scholar]
 Hillier, A. 2019, Phys. Plasmas, 26, 16 [Google Scholar]
 Hillier, A., Barker, A., Arregui, I., & Latter, H. 2019, MNRAS, 482, 1143 [Google Scholar]
 Hillier, A., Van Doorsselaere, T., & Karampelas, K. 2020, ApJ, 897, L13 [Google Scholar]
 Hollweg, J. V. 1971, J. Geophys. Res., 76, 5155 [Google Scholar]
 Hollweg, J. V. 1978, Sol. Phys., 56, 305 [Google Scholar]
 Hollweg, J. V. 1984a, Sol. Phys., 91, 269 [CrossRef] [Google Scholar]
 Hollweg, J. V. 1984b, ApJ, 277, 392 [Google Scholar]
 Hood, A. W., & Priest, E. R. 1979, Sol. Phys., 64, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Hood, A. W., & Priest, E. R. 1981, Geophys. Astrophys. Fluid Dyn., 17, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Howson, T. A., De Moortel, I., & Antolin, P. 2017a, A&A, 607, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Howson, T. A., De Moortel, I., & Antolin, P. 2017b, A&A, 602, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Howson, T. A., De Moortel, I., & Reid, J. 2020, A&A, 636, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Howson, T. A., De Moortel, I., & Pontin, D. I. 2021, A&A, 656, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Jafarzadeh, S., Solanki, S. K., Gafeira, R., et al. 2017, ApJS, 229, 9 [Google Scholar]
 Jess, D. B., Mathioudakis, M., Erdélyi, R., et al. 2009, Science, 323, 1582 [Google Scholar]
 Karami, K., & Bahari, K. 2010, Sol. Phys., 263, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Karami, K., & Barin, M. 2009, MNRAS, 394, 521 [NASA ADS] [CrossRef] [Google Scholar]
 Karampelas, K., & Van Doorsselaere, T. 2018, A&A, 610, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Karampelas, K., Van Doorsselaere, T., & Antolin, P. 2017, A&A, 604, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Karampelas, K., Van Doorsselaere, T., & Guo, M. 2019, A&A, 623, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kohutova, P., Verwichte, E., & Froment, C. 2020, A&A, 633, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kwon, R. Y., & Chae, J. 2008, ApJ, 677, L141 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, J., McIntosh, S. W., De Moortel, I., Threlfall, J., & Bethge, C. 2014, ApJ, 797, 7 [Google Scholar]
 Lohner, R. 1987, Comput. Methods Appl. Mech. Eng., 61, 323 [Google Scholar]
 Magyar, N., & Van Doorsselaere, T. 2016, A&A, 595, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Magyar, N., Duckenfield, T., Van Doorsselaere, T., & Nakariakov, V. M. 2022, A&A, 659, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MartínezGómez, D., Soler, R., Terradas, J., & Khomenko, E. 2022, A&A, 658, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathioudakis, M., Jess, D. B., & Erdélyi, R. 2013, Space Sci. Rev., 175, 1 [Google Scholar]
 Matsumoto, T., & Suzuki, T. K. 2012, ApJ, 749, 8 [Google Scholar]
 McIntosh, S. W., de Pontieu, B., Carlsson, M., et al. 2011, Nature, 475, 477 [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7 [Google Scholar]
 Morton, R. J., Tomczyk, S., & Pinto, R. 2015, Nat. Commun., 6, 7813 [Google Scholar]
 Mumford, S. J., Fedun, V., & Erdélyi, R. 2015, ApJ, 799, 6 [Google Scholar]
 Nakariakov, V. M., & Kolotkov, D. Y. 2020, ARA&A, 58, 441 [Google Scholar]
 Nakariakov, V. M., Ofman, L., Deluca, E. E., Roberts, B., & Davila, J. M. 1999, Science, 285, 862 [Google Scholar]
 Nocera, L., Priest, E. R., & Leroy, B. 1984, A&A, 133, 387 [Google Scholar]
 Ofman, L. 2009, ApJ, 694, 502 [NASA ADS] [CrossRef] [Google Scholar]
 Pascoe, D. J., Goddard, C. R., & Van Doorsselaere, T. 2020, Front. Astron. Space Sci., 7, 61 [Google Scholar]
 Prokopyszyn, A. P. K., Hood, A. W., & De Moortel, I. 2019, A&A, 624, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Raghav, A. N., & Kule, A. 2018, MNRAS, 476, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Rankin, R., Frycz, P., Tikhonchuk, V. T., & Samson, J. C. 1994, J. Geophys. Res., 99, 21291 [Google Scholar]
 Roberts, B. 2019, MHD Waves in the Solar Atmosphere (Cambridge: Cambridge University Press) [Google Scholar]
 Roe, P. L. 1981, J. Comput. Phys., 43, 357 [Google Scholar]
 Ruderman, M. S. 2007, Sol. Phys., 246, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Ryu, D., Jones, T. W., & Frank, A. 2000, ApJ, 545, 475 [Google Scholar]
 Sakurai, T., Goossens, M., & Hollweg, J. V. 1991, Sol. Phys., 133, 227 [Google Scholar]
 Shelyag, S., Fedun, V., Keenan, F. P., Erdélyi, R., & Mathioudakis, M. 2011, Ann. Geophys., 29, 883 [Google Scholar]
 Shelyag, S., Mathioudakis, M., & Keenan, F. P. 2012, ApJ, 753, L22 [Google Scholar]
 Shoda, M., Yokoyama, T., & Suzuki, T. K. 2018, ApJ, 853, 190 [Google Scholar]
 Smith, P. D., Tsiklauri, D., & Ruderman, M. S. 2007, A&A, 475, 1111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Soler, R., Terradas, J., Oliver, R., Ballester, J. L., & Goossens, M. 2010, ApJ, 712, 875 [Google Scholar]
 Soler, R., Terradas, J., Oliver, R., & Ballester, J. L. 2019, ApJ, 871, 3 [Google Scholar]
 Soler, R., Terradas, J., Oliver, R., & Ballester, J. L. 2021, ApJ, 909, 190 [Google Scholar]
 Srivastava, A. K., & Dwivedi, B. N. 2017, A&A, 38, 61 [Google Scholar]
 Srivastava, A. K., Shetye, J., Murawski, K., et al. 2017, Sci. Rep., 7, 43147 [Google Scholar]
 Terradas, J., & Goossens, M. 2012, A&A, 548, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Terradas, J., & Ofman, L. 2004, ApJ, 610, 523 [Google Scholar]
 Terradas, J., Andries, J., Goossens, M., et al. 2008, ApJ, 687, L115 [Google Scholar]
 Terradas, J., Magyar, N., & Van Doorsselaere, T. 2018, ApJ, 853, 35 [Google Scholar]
 Thalmann, J. K., Tiwari, S. K., & Wiegelmann, T. 2014, ApJ, 780, 102 [Google Scholar]
 Tikhonchuk, V. T., Rankin, R., Frycz, P., & Samson, J. C. 1995, Phys. Plasmas, 2, 501 [Google Scholar]
 Tomczyk, S., McIntosh, S. W., Keil, S. L., et al. 2007, Science, 317, 1192 [Google Scholar]
 Van Doorsselaere, T., Srivastava, A. K., Antolin, P., et al. 2020, Space Sci. Rev., 216, 140 [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 WedemeyerBöhm, S., Scullion, E., Steiner, O., et al. 2012, Nature, 486, 505 [Google Scholar]
 Xie, H., Madjarska, M. S., Li, B., et al. 2017, ApJ, 842, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Zaqarashvili, T. V., Zhelyazkov, I., & Ofman, L. 2015, ApJ, 813, 123 [Google Scholar]
All Figures
Fig. 1. Sketch of the coronal fluxtube model. The isosurface of density with value equal to (ρ_{i} + ρ_{e})/2 is shown in orange. Magnetic field lines located at r = 2R, r = R and r = 0 are drawn in green, blue, and red, respectively, for the model with c = 0.4. The bottom and top gray planes represent the solar photosphere. 

In the text 
Fig. 2. Equilibrium radial profiles of the Alfvén speed, v_{A}(r), and sound speed, v_{s}(r), for the twist parameter ranging from c = 0 to c = 0.4. The sound speed is the same in all cases. The region between the vertical dashed lines denotes the inhomogeneous region in density. The values are normalized with respect to the Alfvén speed at r → ∞. We used ρ_{i}/ρ_{e} = 2 and l/R = 0.4. 

In the text 
Fig. 3. Radial profile of the linetied Alfvén frequency continuum for the twist parameter ranging from c = 0 to c = 0.4. The region between the vertical dashed green lines corresponds to the nonuniform region in density. The values are normalized to the external Alfvén frequency. We used ρ_{i}/ρ_{e} = 2 and l/R = 0.4. 

In the text 
Fig. 4. Radial profiles of the azimuthal (top) and longitudinal (bottom) velocity components at t = 0 for the twist parameter ranging from c = 0 to c = 0.4. The inset in the top panel shows a closeup view of the azimuthal velocity around its maximum. Both velocity components are normalized to the maximum amplitude of the perpendicular velocity perturbation, v_{0}. 

In the text 
Fig. 5. Crosssectional cut at z = 0 of the logarithm with base 10 of the vorticity squared, ω^{2}(r,t), in arbitrary units when . The small panel to the right shows the temporal evolution of ω^{2}(r,t) integrated over the whole computational domain, Ω^{2}(t), normalized with respect to the value at . The blue dot marks the position in the curve for the considered simulation time. The complete temporal evolution is available online as a movie. 

In the text 
Fig. 6. Azimuthally averaged density at the tube center, z = 0, for different simulation times. The vertical dashed purple lines show the limits of the nonuniform region at . Density is normalized to the external density. 

In the text 
Fig. 7. Crosssectional cuts of density at the tube center, z = 0, for different simulation times: (left column), (mid column), and (right column). The first row corresponds to the simulation without magnetic twist, c = 0. In descending order, the remaining rows correspond to the simulations with magnetic twist parameters of c = 0.1, 0.2, and 0.4. Density is normalized to the external density. 

In the text 
Fig. 8. Longitudinal density cut at y = 0 for the simulation with magnetic twist parameter, c = 0.2 when , i.e., the same time as in the right column of Fig. 7. Density is normalized to the external density. 

In the text 
Fig. 9. Temporal evolution of the Fourier coefficients at the tube center, z = 0, and, in the middle of the inhomogeneous radial layer, r = R. Top: the Fourier coefficient associated with the torsional mode, p = 0. Bottom: same as the top panel, but for the sum of the first 20 positive and the first 20 negative Fourier coefficients excluding p = 0. In both panels, the twist parameter ranges from c = 0 to c = 0.4. Arbitrary units are used. 

In the text 
Fig. 10. Crosssectional cuts at z = 0 of the logarithm with base 10 of the vorticity squared, ω^{2}(r,t), for different simulation times: (left column), (middle column), and (right column). The first row corresponds to the simulation without magnetic twist, c = 0. In descending order, the remaining rows correspond to the simulations with a magnetic twist parameter of c = 0.1, 0.2, 0,3, and 0.4. Arbitrary units are used. 

In the text 
Fig. 11. Temporal evolution of the vorticity squared integrated in the whole computational domain, Ω^{2}, for the twist parameter ranging from c = 0 to c = 0.4. The curve for c = 0 (dotteddashed red line) is the same as in the small panel of Fig. 5. All curves are normalized to their values at t = 0. 

In the text 
Fig. 12. Azimuthally averaged amplitude spectrum of the Alfvénic (perpendicular to the magnetic field) part of the kinetic energy for the twist parameter ranging from c = 0 to c = 0.4. We considered the last common time for which vorticity increases in all simulations, . Arbitrary units are used. 

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