Transition to turbulence in nonuniform coronal loops driven by torsional Alfvén waves. II. Extended analysis and effect of magnetic twist

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 ﬁeld. Here we explore the role of the magnetic twist. We modeled a coronal loop as a transversely nonuniform straight ﬂux tube, anchored in the photosphere, and embedded in a uniform coronal environment. We considered that the magnetic ﬁeld 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 high-resolution three-dimensional ideal magnetohydrodynamic numerical simulations. We ﬁnd that phase mixing of torsional Alfvén waves creates velocity shear in the direction perpendicular to the magnetic ﬁeld 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 phase-mixing shear. Thus, turbulence is not generated in strongly twisted loops.

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. (2012Pontieu et al. ( , 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. 2011Shelyag et al. , 2012Wedemeyer-Böhm et al. 2012;Mumford et al. 2015;. 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  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 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. 2008Terradas et al. , 2018Antolin 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 2006Ruderman 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.

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 z-axis coincides with the flux tube axis. We considered the same force-free 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 θ = arctan √ q free , where q free is the free energy ratio. The calculation of q free can be achieved from observations through line-of-sight magnetograms, the three-dimensional (3D) position of the coronal loop, and a nonlinear force-free field code. Tables 3 in Aschwanden et al. (2014) and Aschwanden et al. (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 B(r) = B 2 ϕ (r) + B 2 z (r) 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) and Aschwanden et al. (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. The equilibrium Alfvén, v A (r), and sound, v s (r), speeds are where γ is the adiabatic constant. 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. According to Halberstadt & Goedbloed (1993), the line-tied Alfvén frequency continuum is where n is the harmonic number. A continuous spectrum is present because of the radial variation of the z-component 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. 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 t = t/τ A , with τ A = R/v A,∞ .

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 finite-volume shock-capturing spatial discretization on a structured mesh. The equations are as follows: In Eqs. (11)-(14), D Dt ≡ ∂ ∂t + v · ∇ denotes the total derivative, ρ is the mass density, v is the velocity, B is the magnetic field, p is the gas pressure, and µ 0 is the vacuum magnetic permeability. Gravity and nonideal terms are neglected.
The PLUTO code solves Eqs. (11)- (14) in Cartesian coordinates. As in Paper I, the code uses a Roe-Riemann solver (Roe 1981) to compute the numerical fluxes, a piecewise parabolic scheme for spatial reconstruction, and a second-order Runge-Kutta 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 current-free because of magnetic twist, and therefore the background-field-splitting technique cannot be used in the present simulations, in contrast to Paper I. For comparison purposes, the background-field-splitting 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 100x100x100 cells distributed uniformly from −3R to 3R in the x-and y-directions, and from −L/2 to L/2 in the z-direction. 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 1600x1600x1600. 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 y-components 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 line-tying 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.

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. 2011Shelyag et al. , 2012Wedemeyer-Böhm et al. 2012;Mumford et al. 2015;. The waves propagate to the corona through the chrosmosphere, 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 , so that the velocity was purely perpendicular to the magnetic field lines, namely v r (r, z) = 0, 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. 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. 1992Goossens et al. , 2011Giagkiozis 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).

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 thin-layer 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 os-cillate 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 (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 cross-sectional 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.
In the small panel of Fig. 5, the temporal evolution of Ω 2 undergoes three phases (see Paper I). In the first phase fort 65, Ω 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 65 t 87, 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 fort 87, 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 20 √ 2 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;Karampelas et al. 2017;Karampelas & Van Doorsselaere 2018;Guo et al. 2019;Howson et al. 2020Howson et al. , 2021, and some of their results are comparable with ours. Similar vorticity structures as those displayed in Fig. 5 can also be seen in previous 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).

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 cross-sectional 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.
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 z-direction (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.

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 time-varying 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 time-varying 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 time-varying 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 cross-sectional 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 S. Díaz-Suárez & R. Soler: Torsional Alfvén waves in twisted magnetic fields Fig. 5. Cross-sectional cut at z = 0 of the logarithm with base 10 of the vorticity squared, ω 2 (r, t), in arbitrary units when t = 87. 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 t = 0. The blue dot marks the position in the curve for the considered simulation time. The complete temporal evolution is available as a movie. 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.
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, high-order 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 Fig. 7. Cross-sectional cuts of density at the tube center, z = 0, for different simulation times: t = 76 (left column), t = 82 (mid column), and t = 88 (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.
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 phase-mixing-generated 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íaz-Suárez & Soler 2021). 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. Fig. 8. Longitudinal density cut at y = 0 for the simulation with magnetic twist parameter, c = 0.2 when t = 88, i.e., the same time as in the right column of Fig. 7. Density is normalized to the external density. 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.

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 cross-sectional 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.
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 t = 88, a time for which turbulence is already well established.
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 t = 88 is similar to that found at t = 76 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 t = 88, the shape of ω 2 (r, t) in the case with c = 0.3 reveals that some of the phase-mixing-driven concentric rings are slightly deformed in a wavy fashion. No such defor-mation 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;. 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 cross-sectional 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.

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 E = 1 2 ρv 2 ⊥ . 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 cross-sectional 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 E, which is a function of r and t. Then, we applied the continuous Fourier transform in the radial direction, discretized due to the limited Article number, page 9 of 12 A&A proofs: manuscript no. 44175corr Fig. 10. Cross-sectional cuts at z = 0 of the logarithm with base 10 of the vorticity squared, ω 2 (r, t), for different simulation times: t = 76 (left column), t = 82 (middle column), and t = 88 (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. 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  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 E F (k, t) 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, t = 87. 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.
At large scales, kR < 10 0 , the amplitude spectrum is independent of 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 of 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 inten- tionally avoided kR > 10 3 in our analysis owing to the possible role of numerical dissipation at the smallest scales.

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 cascade-like 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 pseudo-2D.
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 just 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 phasemixing-generated 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;. 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).