© The Authors 2024

1. Introduction

The dissipation of waves in the solar corona is considered one of the viable mechanisms that can explain the heating of the coronal plasma, a long-standing, not fully solved problem (Nakariakov & Kolotkov 2020). Photospheric motions would generate such waves, which propagate upward to the corona along the magnetic field connecting different layers of the solar atmosphere, eventually dissipating within the corona.

There is much observational evidence of the presence of different wave modes in the corona, starting from early observations of nonthermal broadening of coronal lines (Feldman et al. 1988; Dere & Mason 1993), and the detection of slow magnetohydrodynamic (MHD) waves (Chae et al. 1998; Ofman et al. 1999) and transverse oscillations in post-flare loops (Nakariakov et al. 1999; Schrijver et al. 1999; Aschwanden et al. 1999). More recently, the excitation of standing kink oscillations in coronal loops has been studied (e.g., Zimovets & Nakariakov 2015; Goddard et al. 2016) that could be excited by coronal impulsive events through a variety of mechanisms (see also Nakariakov et al. 2021 for a recent review). While the above oscillations are subject to damping within a few wave periods, other oscillations have been detected mainly in quiet loops, with no apparent damping for many periods or even with temporary growing amplitude (Wang et al. 2012; Tian et al. 2012; Nisticò et al. 2013). Wave period scaling linearly with the loop length indicates that the decayless oscillations are standing waves (Anfinogentov et al. 2015). Decayless oscillations could be driven by the interaction between loops and quasi-steady flows (Nakariakov et al. 2016), or by continuous footpoint driving (Karampelas et al. 2017, 2019; Afanasyev et al. 2020). Moreover, it has been proposed that they could be related to Kelvin-Helmholtz (KH) rolls induced by transverse loop oscillations (Antolin et al. 2016). The interplay between torsional waves and turbulent dissipation was studied in detail by Guo et al. (2019). These authors adopt various drivers located at the foot-points of the loop, and find that the heating is more pronounced for the simulation containing the mixed kink and Alfvén driver, compared to the two single drivers. Moreover, KH vortices can be induced by phase-mixed standing Alfvén modes, where the KH instability acts as an agent to dissipate wave energy. In addition, the ubiquitous presence of propagating, mainly transverse waves has been revealed by both ground-based (CoMP, Tomczyk et al. 2007; Tomczyk & McIntosh 2009) and space-based (AIA, on board SDO McIntosh et al. 2011; Lemen et al. 2012) observations, with a dominance of outward to inward wave power (Tomczyk & McIntosh 2009). These oscillations have been interpreted as Alfvén waves, though it has also been proposed that they could be fast-mode kink waves propagating in cylindrical structures (Van Doorsselaere et al. 2008). Their ubiquity seems to indicate a possible role in coronal heating.

From a theoretical point of view, wave dissipation in the coronal plasma represents a nontrivial problem. The Reynolds and Lundquist numbers are estimated to be extremely large in the solar corona. This implies that a wave can be efficiently dissipated only if its energy is transferred to fluctuations at very small spatial scales. Such a process could require a very long dissipative time (Terradas & Arregui 2018), which could be much longer than the coronal cooling timescale. The presence of inhomogeneities in the background structures can speed up the formation of small scales. In particular, when a transverse gradient of the Alfvén velocity is present, Alfvén waves are subject to phase mixing (Heyvaerts & Priest 1983; McLaughlin et al. 2011). During phase mixing, wavefronts are bent due to differences in their propagation velocity on nearby field lines. This results in a progressive generation of increasingly small scales in the direction transverse to the background magnetic field. Phase mixing has been extensively studied both by following a normal-mode approach (Steinolfson 1985; Califano et al. 1990, 1992) and by considering the evolution of an initial disturbance (Lee & Roberts 1986; Malara et al. 1992, 1996). Effects of density stratification and magnetic line divergence have also been considered (Ruderman et al. 1998), as well as nonlinear coupling with compressive modes (Nakariakov et al. 1997, 1998). Phase mixing in 3D configurations in the small wavelength limit has been studied (Petkaki et al. 1998; Malara et al. 2000), also in simplified models of quiet-Sun regions (Malara et al. 2003, 2005, 2007) or open fieldline (Malara 2013; Pucci et al. 2014) regions. It has also been shown that phase mixing is active in loops with a radial density inhomogeneity when azimuthally polarized Alfvén waves are excited by coupling with kink modes (Pagano & De Moortel 2017). Also, phase mixing in complex coronal plasma has been studied in Howson et al. (2020), showing that the wave energy dissipation is enhanced in more complex background magnetic fields.

Torsional Alfvén waves represent one possible oscillating mode in a cylindrically symmetric structure with an axial magnetic field; they represent particular azimuthally polarized waves where fluctuations do not depend on the azimuthal angle θ. They are non compressive and are analogous to shear Alfvén waves in slab geometry. Their evolution is only due to phase mixing, except for large amplitudes when nonlinear effects come into play. In this latter case, nonlinear evolution of torsional Alfvén waves can trigger the KH instability and higher azimuthal wave numbers with respect to the torsional mode are excited in the system, thus driving the flux tube to a turbulent state (Díaz-Suárez & Soler 2021, 2022). Therefore, in studying phase mixing, torsional Alfvén waves represent the most suitable case to be considered. Such waves could be excited by torsional motions at the loop bases. This kind of motion has been revealed in the lower layers of the solar atmosphere. In the photosphere, vortical motions seem to be related to convection (Brandt et al. 1988), mainly located at the downdrafts where the plasma returns to the solar interior after cooling down (Bonet et al. 2008, 2010). Small-scale swirl events have been revealed in the quiet-Sun chromosphere (Wedemeyer-Böhm 2009; Tziotziou et al. 2018), as well as rapidly rotating magnetic structures, ubiquitously distributed in the transition region (Wedemeyer-Böhm et al. 2012). Automated detection of chromospheric swirls has recently been performed (Dakanalis et al. 2022); a mean lifetime 3.4 of min with typical diameter of 1–1.5 Mm was found, with about 40% of the analyzed surface covered by swirls. Torsional waves have been detected in magnetic funnels, with an oscillation period of 126–700 s (Jess et al. 2009). An investigation of a small-scale chromospheric tornado (Tziotziou et al. 2020) has suggested the existence of waves propagating upward with phase speeds of ∼20 − 30 km s⁻¹, in the form of both kink modes and localized Alfvénic torsional waves. Numerical simulations of torsional waves in the magnetic field of a chromospheric funnel have been performed to study frequency filtering (Fedun et al. 2011). Moreover, the presence of Alfvén waves has been detected in simulations of photospheric vortexes (Shelyag et al. 2013).

In addition to phase mixing, turbulence represents an alternative way of generating small scales in a magnetofluid. Turbulence in the solar wind has been characterized using a huge number of in situ measurements for many decades (see, e.g., Bruno & Carbone 2013 for a review). Recent direct measurements of turbulence were performed by the Parker Solar Probe spacecraft in the most external solar corona where the wind is still subsonic or sub-Alfvénic (Kasper et al. 2021; Bandyopadhyay et al. 2022; Zhao et al. 2022). Moreover, indirect indications exist of turbulent fluctuations in the corona. In particular, the nonthermal broadening of coronal spectral lines could be indicative of turbulent fluctuations (Banerjee et al. 1998; Singh et al. 2006; Hahn & Savin 2013, 2014), as well as a f⁻¹ frequency spectrum found in CoMP observations of loops (Morton et al. 2016, 2019; Magyar & Van Doorsselaere 2022). Both phase mixing and turbulence can generate small scales through different mechanisms. In the first case, mode couples with the inhomogeneity of the background structure, while in the second one, the generation of small scales is due to nonlinear coupling between fluctuations. This typically generates power-law spectra in the wave number space. Similar to phase mixing, the magnetohydrodynamic (MHD) turbulent cascade preferentially generates small scales in the direction perpendicular to the background magnetic field. This especially holds in the case of a strong magnetic field (low plasma β), as in the case of the coronal plasma (Shebalin et al. 1983). In Alfvénic turbulence, nonlinear couplings occur between fluctuations propagating in opposite senses. This situation can be easily envisioned in closed magnetic structures like loops. However, even in open-fieldline regions, partial wave reflection due to longitudinal inhomogeneities can give rise to oppositely propagating fluctuations and a consequent activation of a nonlinear cascade. Models of coronal heating based on turbulence have been formulated for coronal loops (Van Ballegooijen et al. 2011, 2017; Downs et al. 2016; Van Ballegooijen & Asgari-Targhi 2018; Rappazzo et al. 2017), and for open structures (Moriyasu et al. 2004; Verdini et al. 2009; Antolin & Shibata 2010; Perez & Chandran 2013; Woolsey & Cranmer 2015; Chandran & Perez 2019). Hybrid shell models are simplified turbulence models based on reduced MHD and hold for a low-β plasma dominated by nearly transverse noncompressive fluctuations. They have been applied to the heating of coronal loops (Nigro et al. 2004, 2005; Reale et al. 2005; Buchlin & Velli 2007), and to characterize the nonlinear energy spectral flux (Nigro et al. 2004; Malara et al. 2010). Such models reproduce power-law frequency spectra (Nigro et al. 2020) similar to those observed in the corona (Tomczyk & McIntosh 2009; McIntosh et al. 2011; Morton et al. 2016)

A situation where both phase mixing and turbulence are active is when an initial standing kink mode resonantly couples with Alfvénic oscillations located at the interface between the loop interior and exterior (see, e.g. Davila 1987 for the resonant absorption heating in the solar corona.) The spatial variation of the Alfvén velocity at the interface generates phase mixing in the Alfvénic azimuthal oscillation, giving origin to a KH instability. The final result is a turbulent state in the interface region that eventually leads to wave energy dissipation. This scenario has been studied in great detail (Terradas et al. 2008a; Antolin et al. 2014; Magyar et al. 2015; Antolin & Van Doorsselaere 2019), also in comparison with observations (Antolin et al. 2017), pointing out the role played by KH rolls in enhancing dissipation, as well as in the formation of fine strand-like structures (Antolin et al. 2014). Finally, it has been shown that in a plasma with transverse inhomogeneities phase mixing can result in a phenomenology resembling that of turbulence also in the case in which uni-directional waves are present (Magyar et al. 2017, 2019).

The present paper is devoted to studying the interplay between phase mixing and turbulence in the generation of small scales and consequent wave dissipation in a simple model of a coronal loop. Since we are interested in selecting the effects of phase mixing, we consider torsional modes with a possible superposition of turbulent perturbation here, excluding other more complex situations such as kink modes. We show that the simultaneous presence of phase mixing and turbulent cascade can act in a synergistic way, so as to reinforce both effects. This can lead to an efficient dissipation, even in situations where small amplitude waves propagate in weakly dissipative plasma, as for the solar corona.

The plan of the paper is the following. An overview of the model is presented in Sect. 2; we report some relevant scaling laws in Sect. 2.1 and the physical model we describe in Sects. 2.3–2.5. In Sect. 3 we report our numerical results. We first focus on the purely Alfvén wave perturbations in Sect. 3.1 and on the turbulence in a homogeneous background in Sect. 3.2, moving to the combined case in Sect. 3.3, both for propagating and standing waves. Finally, in Sect. 4 we discuss our results and its implications.

2. Model

We aim to model a situation where transverse fluctuations propagate and evolve inside a coronal loop characterized by a transverse inhomogeneity. We are interested in studying how the coupling between background inhomogeneity and fluctuations and nonlinear couplings among fluctuations lead to the formation of small-scale fluctuations and eventually to their dissipation. As described below, the loop is modeled as a straight magnetic flux tube with a cylindrical symmetry (at the initial time) where the density in the inner part of the loop is larger than in the outer part, while the magnetic field is axially directed and nearly uniform, due to a low plasma-β. This results in a transverse modulation of the Alfvén velocity $c_{\mathrm{A}}(r)=B(r)/\sqrt{4\pi \rho (r)}$, where r is the radial coordinate with respect to the loop axis, as it can be appreciated in Fig. 1. Different kinds of fluctuations will be considered, namely, (i) torsional Alfvén waves, (ii) turbulent transverse fluctuations, and (iii) a superposition of (i) and (ii).

Fig. 1. 1D profiles of the equilibrium quantities B0, cA, P0, and ρ0 as functions of y, for x = L/2 and z = ΛL/2. In the inset we report a zoom-in for the pressure profile described by Eq. (13) (dashed line) and for the correction described by Eq. (14) (solid line), near the boundary y = 2π. The correction acts in the same way at y = 0 for symmetry.

2.1. Phenomenology and scaling laws

Prior to describing the model and its results, we discuss some general results concerning wave dynamics in the considered situation and their implication for our model.

When an Alfvénic fluctuation propagates in a plasma with a transverse inhomogeneous Alfvén velocity, it is subject to phase mixing that progressively generates small scales perpendicularly to the background magnetic field B₀ (e.g., Heyvaerts & Priest 1983). We indicate by k_⊥ and k_|| the wavevector components perpendicular and parallel to B₀, respectively. During phase mixing k_|| remains constant: k_|| = k_||0, while k_⊥ linearly increases with time t, according to

$$\begin{array}{c}\hfill k_{\perp }(t)\sim k_{||0}\frac{\mathrm{d}{c}_{\mathrm{A}}}{\mathrm{d}r}t,\end{array}$$(1)

where dc_A/dr is a typical value for the transverse gradient of the Alfvén velocity (Kaneko et al. 2015). For simplicity, in Eq. (1) we have assumed that k_⊥(t = 0) = 0. Following Eq. (1), we define the phase-mixing dynamical time t_PM(k_⊥) as

$$\begin{array}{c}\hfill t_{\mathrm{PM}}(k_{\perp })\sim \frac{k_{\perp }}{k_{||0}}{\left(\frac{\mathrm{d}{c}_{\mathrm{A}}}{\mathrm{d}r}\right)}^{-1},\end{array}$$(2)

which is also the time that is needed for a given perpendicular wavevector k_⊥ to double. The stronger the transverse gradient of the Alfvén speed, the faster phase mixing. Moreover, t_PM increases with increasing k_⊥. Therefore, as a mechanism able to generate small scales, phase mixing becomes progressively less efficient with increasing the fluctuation wavevector.

Indicating by ${\ell }_{\mathrm{d}}\sim k_{\perp d}^{-1}$ the dissipation length and with k_⊥d the corresponding perpendicular wavevector, we define the dissipative time $t_{\mathrm{d}}^{\mathrm{PM}}$ as the time phase mixing takes to increase k_⊥(t) up to k_⊥d. Using Eq. (1) we obtain

$$\begin{array}{c}\hfill t_{\mathrm{d}}^{\mathrm{PM}}\sim \frac{k_{\perp d}}{k_{||0}}{\left(\frac{\mathrm{d}{c}_{\mathrm{A}}}{\mathrm{d}r}\right)}^{-1}\sim \frac{L_{||}}{{\ell }_{\mathrm{d}}}{\left(\frac{\mathrm{d}{c}_{\mathrm{A}}}{\mathrm{d}r}\right)}^{-1},\end{array}$$(3)

where we assume that the initial parallel wavelength is of the order on the loop length L_||. The above expression can be normalized to a typical timescale t_A = L_⊥/c_A0, where L_⊥ is the loop width and c_A0 is the mean Alfvén velocity. Our definition of t_A is based on the fact that we use the same L_⊥ for all the simulations and vary L_∥, but other definitions of t_A have been used in the literature. We obtain

$$\begin{array}{c}\hfill \frac{t_{\mathrm{d}}^{\mathrm{PM}}}{t_{\mathrm{A}}}\sim \frac{c_{\mathrm{A}0}}{{\ell }_{\mathrm{d}}}{\left(\frac{\mathrm{d}{c}_{\mathrm{A}}}{\mathrm{d}r}\right)}^{-1}\frac{L_{||}}{L_{\perp }},\end{array}$$(4)

where Λ = L_||/L_⊥ is the loop aspect ratio. To give an order-of-magnitude estimate for the ratio Eq. (4), we assume that the density inside the loop is a factor 2 larger than outside and the magnetic field is nearly uniform as set in our numerical model (see Fig. 1). Therefore, the Alfvén velocity has a relative variation of about 30% from inside to outside the loop. We also assume that c_A varies over a scale Δr ∼ L_⊥/10, with L_⊥ ∼ 10³ km and a loop aspect ratio L_||/L_⊥ ∼ 30. In the coronal plasma, the Lundquist number is extremely high and dissipation is probably due to kinetic effects; therefore, we consider a dissipation length on the order of the ion skin depth: ℓ_d ∼ 0.1 km. Using those values, we obtain an estimation for the normalized phase-mixing dissipative time: $t_{\mathrm{d}}^{\mathrm{PM}}/t_{\mathrm{A}}\sim {10}^5$. Such a large value indicates that phase mixing alone is inefficient in dissipating Alfvén waves, at least in the simple situation here considered.

Another possible dissipative mechanism is represented by a turbulent cascade, where fluctuating energy is transferred to smaller scales by nonlinear couplings. In MHD, the turbulent cascade mainly takes place perpendicular to the background magnetic field (e.g., Shebalin et al. 1983), and this effect is even stronger in the low-β coronal plasma. Since the turbulent cascade can virtually generate indefinitely small scales in a finite time, the dissipative time is larger but on the order of the nonlinear time t_NL(ℓ_⊥0), where ℓ_⊥0 is the large energy-containing scale of fluctuations (Matthaeus et al. 2014), namely

$$\begin{array}{c}\hfill t_{\mathrm{d}}^{\mathrm{turb}}\sim \mathrm{\Gamma }{t}_{\mathrm{NL}}({\ell }_{\perp 0})=\mathrm{\Gamma }\frac{{\ell }_{\perp 0}}{c_{\mathrm{A}0}}\frac{B_0}{\delta B({\ell }_{\perp 0})},\end{array}$$(5)

where δB(ℓ_⊥0) is the amplitude of magnetic field fluctuation at the scale ℓ_⊥0, c_A0 is the mean Alfvén speed and Γ is a numerical factor of the order of a few units, which we numerically estimate to identify the dissipative time scale of the specific turbulence fluctuations we inject in the plasma, as discussed below. Normalized to t_A, it is

$$\begin{array}{c}\hfill \frac{t_{\mathrm{d}}^{\mathrm{turb}}}{t_{\mathrm{A}}}\sim \mathrm{\Gamma }\frac{{\ell }_{\perp 0}}{L_{\perp }}\frac{B_0}{\delta B({\ell }_{\perp 0})}.\end{array}$$(6)

We notice that the dissipative time increases with decreasing fluctuation amplitudes δB. Assuming that the largest scale in fluctuation is on the order of the loop width: ℓ_⊥0 ∼ L_⊥ and a normalized fluctuation amplitude δB/B₀ ∼ 10⁻², we obtain for the normalized dissipative time: $t_{\mathrm{d}}^{\mathrm{turb}}/t_{\mathrm{A}}\sim 100\mathrm{\Gamma }$. This value is many orders of magnitude smaller than the above estimation made for the phase-mixing case. Therefore, in a high Reynolds–Lundquist plasma, as in the corona, turbulence is more efficient than phase mixing to dissipate fluctuations, even for small-amplitude fluctuations.

However, the above picture could change in a more realistic case where turbulence evolves inside an inhomogeneous background, such as a coronal loop. In such a case, besides nonlinear effects driving the turbulent cascade, the coupling between the background inhomogeneity and fluctuations (that is at the bases of phase mixing) also contributes to generating small scales. Therefore, it can be expected that in such a configuration phase mixing and nonlinear effects work together to promote small-scale generation and speed-up dissipation. In particular, let us imagine a situation where, at the initial time, the fluctuation amplitude at large scale δB(ℓ_⊥0) is small enough to have $t_{\mathrm{PM}}({\ell }_{\perp 0})<t_{\mathrm{d}}^{\mathrm{turb}}$. Using Eqs. (2) and (5) this condition can be expressed by

$$\begin{array}{c}\hfill \frac{\delta B}{B_0}<\frac{\mathrm{d}{c}_{\mathrm{A}}}{\mathrm{d}r}\frac{{\ell }_{\perp 0}^2}{L_{\Vert }}\frac{\mathrm{\Gamma }}{c_{\mathrm{A}0}}.\end{array}$$(7)

where δB ≡ δB(ℓ_⊥0). In such a case, at an early stage of the time evolution phase, mixing would proceed faster than nonlinear effects to generate small scales. Moreover, phase mixing leaves the fluctuation amplitude δB(ℓ_⊥) unchanged as ℓ_⊥ decreases. Therefore, as ℓ_⊥ decreases in time, t_NL(ℓ_⊥) will proportionally decrease, too. Hence, at a time t^* and at a corresponding scale ${\mathcal{l}}_{\perp }^{*}$, the situation will be reversed, namely, $t_{\mathrm{d}}^{\mathrm{turb}}({\ell }_{\perp })<t_{\mathrm{PM}}({\ell }_{\perp })$ for ${\mathcal{l}}_{\perp }<{\mathcal{l}}_{\perp }^{*}$. From that time on, nonlinear effects become the faster mechanism producing small scales and the cascade will bring energy to dissipative scales in a time $\gtrsim t_{\mathrm{d}}^{\mathrm{turb}}({\ell }_{\perp }^{\ast })$, which is smaller than the initial $t_{\mathrm{d}}^{\mathrm{turb}}$. A similar mechanism has been recently discussed by Díaz-Suárez & Soler (2021, 2022) showing that turbulence – driven by the nonlinear evolution of phase-mixed torsional Alfvén waves – is able to transport fluctuations’ energy down to dissipative scales much faster compared to what predicted by the theory of linear phase mixing. In this scenario, the two considered mechanisms work in a synergistic way: phase mixing dominates at early times (t ≲ t^*) and nonlinear cascade at later times (t ≳ t^*). The result is a dissipative time that is shorter than those of phase mixing and of the nonlinear cascade when individually considered. In what follows, we will explore the above ideas from a quantitative point of view by means of numerical simulations.

2.2. MHD equations and the numerical method

To describe the evolution of perturbations in the coronal plasma we use the compressible MHD equations that are written in the following form, using dimensionless quantities

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}=-\mathbf{\nabla }·\left(\rho \mathit{v}\right),\hfill \end{array}$$(8)

$$\begin{array}{cc}& \frac{\partial \mathit{v}}{\partial t}=-(\mathit{v}·\mathbf{\nabla })\mathit{v}+\frac{1}{\rho }[(\mathbf{\nabla }\times \mathit{B})\times \mathit{B}]\hfill \\ \hfill & -\frac{\beta }{2\rho }\mathbf{\nabla }\left(\rho T\right)-{\nu }_4{\mathbf{\nabla }}^4\mathit{v},\hfill \end{array}$$(9)

$$\begin{array}{cc}& \frac{\partial \mathit{A}}{\partial t}=\mathit{v}\times \mathit{B}-{\eta }_4{\mathbf{\nabla }}^4\mathit{A},\hfill \end{array}$$(10)

$$\begin{array}{cc}& \frac{\partial T}{\partial t}=-(\mathit{v}·\mathbf{\nabla })T-(\gamma -1)(\mathbf{\nabla }·\mathit{v})T-{\chi }_4{\mathbf{\nabla }}^{4}T.\hfill \end{array}$$(11)

Here, ρ is the density normalized to a typical value $\stackrel{\sim }{\rho }$; B is the magnetic field normalized to a typical value $\stackrel{\sim }{B}$; A is the vector potential normalized to $\stackrel{\sim }{\ell }\stackrel{\sim }{B}$, with $\stackrel{\sim }{\ell }$ a typical length; A is related to B through B = ∇ × A; v is the magnetofluid velocity normalized to the typical Alfvén velocity ${\stackrel{\sim }{c}}_{\mathrm{A}}=\stackrel{\sim }{B}/\left[{(4\pi \stackrel{\sim }{\rho })}^{1/2}\right]$; T is the temperature normalized to a typical value $\stackrel{\sim }{T}$; γ is the adiabatic index. Spatial coordinates x, y and z are normalized to $\stackrel{\sim }{\ell }$ and time is normalized to the Alfvén time $\stackrel{\sim }{\ell }/{\stackrel{\sim }{c}}_{\mathrm{A}}$. The coefficient β is the plasma beta, defined as the ratio between the typical kinetic and the magnetic pressure, $\beta =8\pi {\kappa }_{\mathrm{B}}\stackrel{\sim }{\rho }\stackrel{\sim }{T}/(\mu m_{\mathrm{p}}{\stackrel{\sim }{B}}^2)$, with κ_B the Boltzmann constant, μ the mean atomic weight and m_p the proton mass. In the coronal plasma, it is typically β ≪ 1; we used the value β = 0.05. In code units, the plasma pressure is given by P = βρT/2.

Compressibility is included in the model. As discussed, for instance in Malara et al. (1996), the production of small scales necessary for the occurrence of dissipation is more effective in a compressible medium rather than in an incompressible one. The energy Eq. (11) expresses the adiabaticity condition (except for the thermal hyper-diffusivity, see below). In fact, since we are interested in describing the route to energy dissipation rather than dissipation itself, we use the simplified energy Eq. (11) where effects such as radiative losses are neglected. Integrating Eq. (10), which describes the time evolution of the vector potential A instead of the magnetic field B, guarantees ∇ ⋅ B = 0. Moreover, a logarithmic regularization to the density is used; namely, we set ρ = e^g, and solve the equivalent equation for g with the purpose of a better description of possible discontinuities and shocks. Hyper-dissipation is implemented through the last terms in Eqs. (9)–(11) representing hyper-viscosity, hyper-resistivity and hyper thermal diffusivity, respectively, with ν₄, η₄ and χ₄ the corresponding coefficients. These terms have been introduced to dissipate the turbulent cascade at small scales and control numerical stability; they do not mimic any physical process. With respect to standard dissipation, hyper-dissipation allows one to obtain more extended fluctuation spectra, keeping numerical stability.

Equations (8)–(11) are solved in a 3D Cartesian domain with periodicity boundary conditions along the three space directions x, y, and z. In particular, z represents the direction of the background magnetic field B₀. The 3D Cartesian domain D in normalized units is defined as D = {x, y, z}=[0 : L]×[0 : L]×[0 : ΛL], with L = 2π. Throughout the paper, we indicate the domain size perpendicular to B₀ by L, and the parallel size by ΛL = L_∥, where Λ = L_∥/L is a free parameter that represents the domain aspect ratio.

To solve Eqs. (8)–(11) in 3D configurations, we used the “COmpressible Hall Magnetohydrodynamics simulator for Plasma Astrophysics” (COHMPA) algorithm (Pezzi et al. 2024) with the Hall term switched off. A pseudo-spectral algorithm is used, which adopts a second-order Runge-Kutta scheme to advance in time the MHD variables. In the physical domain D all quantities are calculated on a regular grid formed by N_⊥ × N_⊥ × N_|| grid points: (x_i,y_j,z_l) = (iL/N_⊥,jL/N_⊥,lΛL/N_||), with i, j, l integers, 0 ≤ i, j ≤ N_⊥ − 1 and 0 ≤ l ≤ N_|| − 1. In all numerical runs, we used N_⊥ ≥ N_||. Correspondingly, a discrete set of wavevectors is defined in the spectral space: k = (2πn_x/L,2πn_y/L,2πn_z/(ΛL)), with n_x, n_y, n_z integers, 0 ≤ n_x, n_y ≤ N_⊥/2 − 1 and 0 ≤ n_z ≤ N_||/2 − 1. A spectral filter is employed to suppress numerical artefacts due to aliasing (Orszag 1971; Meringolo & Servidio 2021), and a 2/3 rule is adopted where the spectral coefficients of all quantities are set to zero at k = |k|> 2/3k_N, with k_N = (2π/L)(N_⊥/2 − 1) corresponding to the Nyquist frequency (Canuto et al. 2007). A 2.5D version of this algorithm has already been adopted in literature (Vásconez et al. 2015; Perri et al. 2017; Pezzi et al. 2017). More recently, a fully 3D version of the algorithm has been exploited by Pezzi et al. (2024).

The initial condition is given by the superposition of an ideal (nondissipative) MHD equilibrium and a perturbation. The equilibrium represents a simplified model for a coronal loop. The perturbation represents a torsional propagating or standing Alfvén wave, or a turbulent perturbation, or a combination of both. The explicit forms of the equilibrium and the different kinds of perturbation are given in the following.

2.3. Equilibrium: a cylindrical magnetic flux tube

The equilibrium is a structure with cylindrical symmetry around an axis parallel to the z-axis. All quantities relative to the equilibrium are indicated by the index 0. The magnetic field ${\mathit{B}}_0=B_0(x,\mathit{y})\widehat{\mathit{z}}$ is directed along z, where $\widehat{\mathit{z}}$ is the unit vector in the z direction. Magnetic lines have no curvature, and the Lorenz force per volume unit reduces to the gradient of magnetic pressure. Therefore, the equilibrium is sustained by the balance between magnetic and plasma pressure, expressed by the following relation

$$\begin{array}{c}\hfill \frac{B_0^2(x,y)}{2}+P_0(x,y)=\mathrm{const}.,\end{array}$$(12)

where P₀(x, y) = βρ₀(x, y)T₀/2 is the equilibrium plasma pressure in dimensionless code units, ρ₀(x, y) is the inhomogeneous density and T₀ is the temperature that is assumed to be uniform in space. The inhomogeneities of equilibrium are transverse to the guide field B₀. We have chosen the following functional form for the plasma pressure

$$\begin{array}{c}\hfill P_0(x,y)=\frac{(P_{\mathrm{int}}-P_{\mathrm{ext}})}{2}[1-tanh\left(\frac{r-r_0}{\mathrm{\Delta }r}\right)]+P_{\mathrm{ext}}.\end{array}$$(13)

Here and in what follows, subscripts int and ext indicate values of quantities in regions internal or external to the flux tube, respectively. Moreover, (x₀, y₀) = (L/2, L/2) = (π, π) is the position of the symmetry axis, $r=\sqrt{{(x-x_0)}^2+{(\mathit{y}-{\mathit{y}}_0)}^2}$ is the distance of a point from the symmetry axis, r₀ = L/4 = π/2 is the radius of the flux tube and Δr = L/16 = π/8 is a parameter controlling the width of the shear region separating the interior from the exterior of the flux tube.

Equation (13) describes a cylindrically symmetric function. On the other hand, the domain D is a parallelepiped with periodic boundary conditions. Even though the expression (13) for P₀(x, y) is periodic on the boundaries of D, its derivatives normal to the domain boundaries are not periodic, which could provoke numerical problems. Therefore, we have modified the expression for P₀(x, y) so that all its first-order partial derivatives are periodic. Indicating the original form of P₀ (Eq. (13)) by P_old and the corrected form by P_corr, the two are related by

$$\begin{array}{c}\hfill P_{\mathrm{corr}}(x,y)=P_{\mathrm{old}}(x,y)+a(y){(x-\frac{L}{2})}^2+b(x){(y-\frac{L}{2})}^2.\end{array}$$(14)

The functions a(y) and b(x) are chosen such as the first-order normal derivatives of P_corr are vanishing all along the domain boundaries. This condition is verified by the following choice:

$$\begin{array}{cc}& a(y)=\frac{1}{L}\frac{\partial P_{\mathrm{old}}}{\partial x}(0,y)+\frac{1}{2L^2}{(y-\frac{L}{2})}^2\frac{{\partial }^2{P}_{\mathrm{old}}}{\partial x\partial y}(0,0),\hfill \end{array}$$(15)

$$\begin{array}{cc}& b(x)=\frac{1}{L}\frac{\partial P_{\mathrm{old}}}{\partial y}(x,0)+\frac{1}{2L^2}{(x-\frac{L}{2})}^2\frac{{\partial }^2{P}_{\mathrm{old}}}{\partial x\partial y}(0,0).\hfill \end{array}$$(16)

Taking into account that

$$\begin{array}{cc}& \frac{\partial P_{\mathrm{old}}}{\partial x}(0,y)=-\frac{\partial P_{\mathrm{old}}}{\partial x}(L,y),\frac{\partial P_{\mathrm{old}}}{\partial y}(x,0)=-\frac{\partial P_{\mathrm{old}}}{\partial y}(x,L),\hfill \end{array}$$(17)

$$\begin{array}{cc}& \frac{{\partial }^2{P}_{\mathrm{old}}}{\partial y\partial x}(0,y)=-\frac{{\partial }^2{P}_{\mathrm{old}}}{\partial y\partial x}(L,y),\frac{{\partial }^2{P}_{\mathrm{old}}}{\partial x\partial y}(x,0)=-\frac{{\partial }^2{P}_{\mathrm{old}}}{\partial x\partial y}(x,L),\hfill \end{array}$$(18)

it can be easily checked that (∂P_corr/∂x)(0, y) = (∂P_corr/∂x)(L, y) and (∂P_corr/∂y)(x, 0) = (∂P_corr/∂y)(x, L). We note that, for the values of parameters used in the model, the corrective terms in Eq. (14) are much smaller than the original expression (13). We verified that this procedure essentially removes numerical problems related to the fulfillment of periodicity. From now on, the expression of the equilibrium pressure will be given by P_corr(x, y), even though we continue to indicate it by P₀(x, y) to simplify the notation. A similar procedure was used in Vásconez et al. (2015).

The magnetic field is obtained from Eq. (12) as a function of P₀:

$$\begin{array}{c}\hfill B_0(x,y)=\sqrt{B_{\mathrm{ext}}^2+2P_{\mathrm{ext}}-2P_0(x,y)}.\end{array}$$(19)

Finally, the equilibrium density is proportional to P₀(x, y) according to

$$\begin{array}{c}\hfill {\rho }_0(x,y)=\frac{2P_0(x,y)}{\beta T_0}.\end{array}$$(20)

In particular, it is P_ext = βρ_extT₀/2 and P_int = βρ_intT₀/2.

To completely define the equilibrium model, we have to specify the values of parameters ρ_ext, ρ_int, B_ext, and T₀. We used the values: ρ_ext = B_ext = T₀ = 1, and ρ_int = 2 so that the density inside the flux tube is assumed to be a factor 2 larger than the external density. For β = 0.05, this gives the values P_ext = 0.025 and P_int = 0.05 for the plasma pressure outside and inside the flux tube, respectively. Another relevant quantity is the Alfvén velocity, which is defined by $c_{\mathrm{A}}(x,\mathit{y})=B_0(x,\mathit{y})/\sqrt{{\rho }_0(x,\mathit{y})}$ in the code units, and varies perpedicularly to B₀.

Figure 1 shows 1D profiles of the equilibrium quantities, B₀(L/2, y), c_A(L/2, y), P₀(L/2, y), and ρ₀(L/2, y) as functions of y. Those profiles are taken along a line that crosses the flux tube through the symmetry axis. The variations of density and pressure across the flux tube are clearly visible. In the inset we report a zoom-in for the pressure profile described by Eq. (13) (dashed line) and for the corrected profile described by Eq. (14) (solid line), near the boundary y = 2π. The two profiles are very similar (the relative difference is of the order of 10⁻³), but the corrected profile has a vanishing normal derivative at the boundaries. Due to the low value of β, magnetic pressure dominates plasma pressure, and the magnetic field profile is nearly constant, though B₀ is slightly less intense inside the flux tube. The Alfvén velocity inside the flux tube is lower than outside by a ∼0.7 factor. This variation is at the base of the phase-mixing phenomenon.

Finally, we notice that the unit time (in code units) corresponds to the Alfvén time d_FT/(πc_A, ext) = 1, where d_FT = π is on the order of the flux tube diameter (Fig. 1).

2.4. Alfvénic torsional wave perturbation

The initial perturbation considered at first is a torsional Alfvén wave. This wave is polarized in the azimuthal direction with respect to the flux tube axis. We define the azimuthal angle θ = tan⁻¹[(y−y₀)/(x−x₀)] and the corresponding unit vector $\widehat{\mathit{\theta }}=-sin\theta \widehat{\mathit{x}}+cos\theta \widehat{\mathit{y}}$. The perturbation involves only the transverse components of the magnetic and velocity fields. The magnetic field perturbation reads

$$\begin{array}{c}\hfill \delta {\mathit{B}}^{\mathrm{Aw}}=A_{1}f(r)cos(k_{z}z)\widehat{\mathit{\theta }},\end{array}$$(21)

where the superscript “Aw” indicates quantities related to the torsional Alfvén wave, A₁ is a parameter used to fix the amplitude of the perturbation, k_z = 2π/(ΛL) represents the wavevector in the parallel z direction, and f(r) is a function defining the radial profile of the perturbation. In particular, we choose f(r) as a Tukey window function, that has the form

$$\begin{array}{c}\hfill \{\begin{array}{cc}{\displaystyle f(r)=\frac{1}{2}[1-cos\left(\frac{2r}{a}\right)]}\hfill & 0<r<a\pi /2,\hfill \\ f(r)=1\hfill & a\pi /2\le r\le \pi /2,\hfill \end{array}\end{array}$$(22)

with a = 0.6 and f(π − r) = f(r) for r < π/2. This function guarantees that the perturbation smoothly goes to zero at the domain’s edge and on the loop axis r = 0. We notice that the considered torsional wave has a parallel wavelength that coincides with the parallel domain size: λ_|| = ΛL and does not depend on the azimuthal angle θ (i.e., it is a m = 0 mode). The magnetic field perturbation can be expressed in terms of Cartesian components in the form

$$\begin{array}{c}\hfill \delta {\mathit{B}}^{\mathrm{Aw}}=\delta B_x^{\mathrm{Aw}}\widehat{\mathit{x}}+\delta B_y^{\mathrm{Aw}}\widehat{\mathit{y}}=A_{1}f(r)cos(k_{z}z)(-sin\theta \widehat{\mathit{x}}+cos\theta \widehat{\mathit{y}}).\end{array}$$(23)

The initial perturbation is an Alfvén wave propagating in the positive z direction. Therefore, the velocity field has the form

$$\begin{array}{c}\hfill \delta {\mathit{v}}_{x,y}^{\mathrm{Aw}}(x,y)=-\frac{c_{\mathrm{A}}(r)}{B_0(r)}\delta {\mathit{B}}_{x,y}^{\mathrm{Aw}}(x,y)=-\frac{\delta {\mathit{B}}_{x,y}^{\mathrm{Aw}}(x,y)}{\sqrt{{\rho }_0(r)}}.\end{array}$$(24)

Figure 2 shows 1D profiles of the $\delta {\mathit{B}}_{\mathit{y}}^{\mathrm{Aw}}(x,L/2)$ and $\delta {\mathit{v}}_{\mathit{y}}^{\mathrm{Aw}}(x,L/2)$ components taken along a line that crosses the flux tube axis parallel to the x axis (in Fig. 2 we set A₁ = 1). Moreover, a 2D plot of $\delta {\mathit{B}}_{\mathit{y}}^{\mathrm{Aw}}$ in a plane perpendicular to B₀ is represented in the left panel of Fig. 3. Both magnetic field and velocity perturbations are mainly localized across the shear region, but they are also present in the remaining part of the domain, compatible with the Cartesian geometry of D. The value of the parameter A₁ is determined by imposing that the RMS value of the initial magnetic field perturbation is unitary: $\delta B_{\perp ,\mathrm{rms}}\equiv \sqrt{⟨{(B_x-⟨B_x⟩)}^2+{(B_{\mathit{y}}-⟨B_{\mathit{y}}⟩)}^2⟩}=1$, where ⟨…⟩ indicates a volume average. After this normalization, we multiply the perturbation by a specific value of amplitude that varies for different runs (see Table 1 for details). Note that the perturbation does not affect the remaining physical quantities: δρ = δv_z = δB_z = δT =  0.

Fig. 2. 1D cut representing the components δBy and δvy as functions of x, for the single Alfvén wave perturbation. The section is computed along a straight line parallel to the x-axis, at y = L/2 and z = ΛL/2.

Fig. 3. 2D perpendicular section of the By component for the evolution of Alfvén waves at times t = 0 (left), t = 150 (middle), and t = 300 (right). Here we used a perturbation amplitude A = 10−2, a hyper-dissipation coefficient η = 10−8, and an aspect ratio Λ = 4 (Run 15).

Equations (23) and (24) correspond to a single Alfvén wave propagating along the background magnetic field. However, we also consider the case of a standing Alfvén wave. This can be obtained starting from an initial condition where the magnetic field perturbation has the same form Eq. (23), but the velocity perturbation is initially set to zero. This situation corresponds to two counter-propagating Alfvén waves initially having opposite velocity fields that cancel each other. While the propagating wave perturbation involves both magnetic and kinetic energy (that are equal for an Alfvén wave), the standing wave perturbation initially involves only magnetic energy.

2.5. Turbulent perturbation

Another form of initial perturbation considered here corresponds to a turbulent perturbation. The corresponding magnetic field δB^turb and velocity δv^turb perturbations are written as a superposition of transverse Fourier modes, each polarized in the xy-plane and characterized by a wavevector k = (k_x, k_y, k_z). Both δB^turb and δv^turb are divergence-free, the latter condition to ensure that the initial perturbation is noncompressive, as for Alfvénic perturbations. We describe a procedure to construct δB^turb and δv^turb with the aforementioned characteristics.

To comply with the divergence-free condition, the turbulent magnetic perturbation is expressed as δB^turb = ∇ × A, with A the corresponding vector potential (see, e.g., Servidio et al. 2009).

Since we are interested in fluctuations polarized in the direction perpendicular to B₀, we chose a vector potential having the form $\mathit{A}=\mathit{A}(x,\mathit{y},z)=A_x(z)\widehat{\mathit{x}}+A_{\mathit{y}}(z)\widehat{\mathit{y}}+A_z(x,\mathit{y},z)\widehat{\mathit{z}}$, in which the A_x and A_y components account for modes with k parallel to B₀ and the A_z component accounts for oblique and perpendicular modes. Periodicity allows us to write the vector potential in terms of a Fourier series as $A_x={\displaystyle \sum _{k_z}}{\widehat{A}}_x(k_z)e^{i{k}_{z}z}$, $A_{\mathit{y}}={\displaystyle \sum _{k_z}}{\widehat{A}}_{\mathit{y}}(k_z)e^{i{k}_{z}z}$, $A_z={\displaystyle \sum _{k_x,k_{\mathit{y}},k_z}}{\widehat{A}}_z(k_x,k_{\mathit{y}})e^{i\mathit{k}·\mathit{x}}$, where the wavevector is k = [(2π/L)n_x, (2π/L)n_y, (2π/ΛL)n_z], n_x, n_y and n_z are integers, and (n_x, n_y, n_z) is the wavenumber. The corresponding magnetic field has the form

$$\begin{array}{cc}& \delta B_x^{\mathrm{turb}}={\displaystyle \sum _{k_x,k_y,k_z}}i{k}_y{\widehat{A}}_z(\mathit{k})e^{i\mathit{k}·\mathit{x}}-{\displaystyle \sum _{k_z}}i{k}_z{\widehat{A}}_y(k_z)e^{i{k}_{z}z},\hfill \\ \hfill & \delta B_y^{\mathrm{turb}}={\displaystyle \sum _{k_z}}i{k}_z{\widehat{A}}_x(k_z)e^{i{k}_{z}z}-{\displaystyle \sum _{k_x,k_y,k_z}}i{k}_x{\widehat{A}}_z(\mathit{k})e^{i\mathit{k}·\mathit{x}},\hfill & \hfill \delta B_z^{\mathrm{turb}}=0.\end{array}$$

The Fourier coefficients ${\widehat{A}}_x,{\widehat{A}}_{\mathit{y}},{\widehat{A}}_z$ are chosen such that the spectrum of the initial fluctuations is isotropic in the wavenumber space (see, e.g., Wan et al. 2011). In particular, for $\sqrt{n_x^2+n_{\mathit{y}}^2+n_z^2}\le n_{\max }$, it is $|{\widehat{A}}_x(k_z)|=|{\widehat{A}}_{\mathit{y}}(k_z)|\propto |k_z{|}^{-a}$ and $|{\widehat{A}}_z(k_x,k_{\mathit{y}})|\propto k_{\perp }^{-a}$, with $k_{\perp }=\sqrt{k_x^2+k_{\mathit{y}}^2}$. Otherwise, for $\sqrt{n_x^2+n_{\mathit{y}}^2+n_z^2}>n_{\max }$ it is ${\widehat{A}}_x(k_z)={\widehat{A}}_{\mathit{y}}(k_z)={\widehat{A}}_z(k_x,k_{\mathit{y}})=0$. We used the values a = 1.5 and n_max = 3, giving an initial perturbation localized at the largest spatial scales. The phases of complex Fourier coefficients ${\widehat{A}}_x$, ${\widehat{A}}_{\mathit{y}}$, ${\widehat{A}}_z$ are randomly chosen.

The same procedure is applied to obtain the velocity fluctuation δv^turb. However, the phases of the velocity Fourier coefficients are chosen differently with respect to the magnetic field. In particular, phases are chosen such as the cross helicity H_c = ∫δv^turb ⋅ δB^turb dV of the initial perturbation is small (∼1%) with respect to the fluctuating magnetic energy $E_{\mathrm{M}}={\displaystyle \int {\left(\delta {\mathit{B}}^{\mathrm{turb}}\right)}^2}/2\mathrm{d}V$. This condition is necessary in order to trigger a nonlinear cascade, at least in the case of a uniform background. Finally, as for the Alfvénic torsional wave, both magnetic and velocity fluctuation are normalized to their respective RMS values.

The pattern of the above-described turbulent perturbation is represented in the upper panel of Fig. 4, where a 3D visualization of the B_y component at the initial time t = 0 is drawn.

Fig. 4. 3D visualization of the By component at times t = 0 (top) and t = td (bottom), for Run 14.

3. Numerical results

In this section, we present results from numerical simulations. In order to quantify the evolution of the dissipation rate in our model, we monitored the evolution of the quadratic quantity W = ⟨J²⟩+⟨ω²⟩ in time, where J = ∇ × B and ω = ∇ × v are the rms values of the current density and vorticity, respectively. W represents a proxy for the dissipated power. By summing J² and ω², we are implicitly assuming that viscosity and resistivity have the same values in code units. For each run, we define the dissipation time t_d as the time that corresponds to the maximum dissipation rate, namely W_max = W(t = t_d). To properly quantify the value of t_d, we carried out all simulations up to a time t_final when W becomes definitely smaller than W_max. Each run has been performed using the same value η for the hyper-dissipative coefficients: η = ν₄ = η₄ = χ₄, even though η can be different in different runs. In the following, we refer to the fields in the xy-plane with the subscript ⊥, while ∥ represents the parallel component to the z-axis. We also recall that the Alfvén time is unitary in our normalisation units. All the simulations performed are summarized in Table 1 where we report the values of various parameters.

3.1. Torsional Alfvén wave

We performed a series of runs to characterize the time evolution of the torsional Alfvén wave propagating in the inhomogeneous equilibrium, described in Sects. 2.3 and 2.4, respectively. The main purpose is to verify the adherence of scaling laws characterizing phase mixing by the numerical model. We initialize the fields as

$$\begin{array}{cc}& B_{\Vert }=B_0(x,y),\hfill \end{array}$$(25)

$$\begin{array}{cc}& {\mathit{B}}_{\perp }=A_B\delta {\mathit{B}}_{\perp }^{\mathrm{Aw}}(x,y,z),\hfill \end{array}$$(26)

$$\begin{array}{cc}& v_{\Vert }=0,\hfill \end{array}$$(27)

$$\begin{array}{cc}& {\mathit{v}}_{\perp }=A_{\mathrm{v}}\delta {\mathit{v}}_{\perp }^{\mathrm{Aw}}(x,y,z),\hfill \end{array}$$(28)

$$\begin{array}{cc}& \rho ={\rho }_0(x,y).\hfill \end{array}$$(29)

In the present case of a propagating wave we use A_B = A_v. Due to the normalization of $\delta {\mathit{B}}_{\perp }^{\text{Aw}}$ described in Sect. 2.4, the parameter A = A_B = A_v represents the rms amplitude of the fluctuating magnetic field.

In Fig. 3 we report a 2D perpendicular section at z = ΛL/2 of the B_y component for the high-resolution Run 15, where only the torsional wave is initially present (Table 1), at time t = 0 (left), t = 150 (middle), and t = 300 (right). The last time approximately corresponds to the dissipation time t_d for this test. As the simulation proceeds, phase mixing occurs, generating increasingly smaller length scales in the radial direction (across the background magnetic field). This gradually increases W, until dissipative effects due to hyper-dissipation terms come into play, eroding the fluctuation and reducing W. This happens all across the boundary of the flux tube, where the inhomogeneity is localized.

In order to quantitatively characterize the phenomenon of phase mixing, we can derive an estimate for the dissipative time t_d in numerical runs. At odds with Eq. (3) where t_d has been expressed in terms of the dissipative scale ℓ_d, in the numerical model, such a scale is not fixed a priori. Instead, dissipation is determined by the hyper-dissipative terms included in the model. In particular, we can define the characteristic time t_η associated with hyper-dissipation by balancing the term ∂A/∂t with η∇⁴A in Eq. (10) (or, equivalently, balancing ∂v/∂t with η∇⁴v in Eq. (9)). This gives

$$\begin{array}{c}\hfill t_{\eta }(k_{\perp })\sim \frac{1}{\eta k_{\perp }^4}.\end{array}$$(30)

The dynamical time associated with phase mixing is given by the expression (2) where, in our case, k_||0 = 2π/(ΛL) = Λ⁻¹. This gives

$$\begin{array}{c}\hfill t_{\mathrm{PM}}(k_{\perp })\sim k_{\perp }\mathrm{\Lambda }{\left(\frac{\mathrm{d}{c}_{\mathrm{A}}}{\mathrm{d}r}\right)}^{-1}.\end{array}$$(31)

Dissipation becomes effective at a wavevector k_⊥d where t_PM(k_⊥) becomes on the order of t_η(k_⊥). Therefore, matching Eqs. (30) and (31) we derive

$$\begin{array}{c}\hfill k_{\perp \mathrm{d}}\sim \frac{1}{{(\eta \mathrm{\Lambda })}^{1/5}}{\left(\frac{\mathrm{d}{c}_{\mathrm{A}}}{\mathrm{d}r}\right)}^{1/5}.\end{array}$$(32)

Inserting this expression into Eq. (31) we derive an estimation for the dissipative time for the torsional Alfvén wave case

$$\begin{array}{c}\hfill t_{\mathrm{d}}\sim {\left(\frac{\mathrm{d}{c}_{\mathrm{A}}}{\mathrm{d}r}\right)}^{-4/5}{\eta }^{-1/5}{\mathrm{\Lambda }}^{4/5}.\end{array}$$(33)

Moreover, estimating the current density as J ∼ k_⊥B_⊥, and assuming that phase mixing does not substantially modify the initial wave amplitude A as long as k_⊥ ≲ k_⊥d, we can derive an estimate for the maximum dissipation rate

$$\begin{array}{c}\hfill J_{\max }^2\sim A^2{\left(\frac{1}{\eta \mathrm{\Lambda }}\frac{\mathrm{d}{c}_{\mathrm{A}}}{\mathrm{d}r}\right)}^{2/5}.\end{array}$$(34)

Equations (33) and (34) can be compared with the results of numerical runs in the case of the torsional Alfvén wave. In particular, the dependence of the dissipation on the wave amplitude has been studied by performing some numerical tests (Run 1 to Run 4) where different values of A have been used, varying A over one order of magnitude: A = 0.02, 0.05, 0.1, 0.2. Since a high resolution is not required for these runs, we used a grid with N_⊥ = N_|| = 64 mesh points for each direction. The domain D is a cube (Λ = 1) of size L = 2π. Due to the low number of gridpoints, the hyper-dissipative coefficients have the relatively high value η = 5 × 10⁻⁶. In Fig. 5, we report the time history of W/A² for those runs (Run 1–4). It appears that the four curves are almost superposed, indicating that the dissipation rate W is proportional to A², in accordance with Eq. (34). Moreover, the dissipative time t_d, corresponding to the time of maximum W, is independent of the wave amplitude A, as predicted by Eq. (33).

Fig. 5. Time history of W/A2 for the propagating Alfvén wave, with different values of amplitude A, for Runs 1–4. For these simulations the perturbation amplitude refers both to the magnetic and the velocity fields, namely A = AB = Av.

In Fig. 5 the normalized profile W(t)/A² corresponding to the lowest amplitude A = 0.02 is slightly higher than the others. This small discrepancy can be attributed to the power developed by the dissipation of the current associated with the equilibrium magnetic field B₀(x, y). Of course, for low values of the wave amplitude A this effect becomes more visible in comparison with wave dissipation rate. However, even for low wave amplitudes the dissipation of the equilibrium current remains sufficiently small to be ignored. We also observe that the dissipation of the equilibrium current is completely negligible in the high-resolution runs (Run 14–26) described in Sects. 3.2–3.4, where a much lower value for the dissipative coefficient η has been used.

Another set of low-resolution runs has been performed, where the aspect ratio Λ of the spatial domain has been increased from Λ = 1 up to Λ = 10. For all those runs, we set an amplitude A = 0.02, the coefficient η = 5 × 10⁻⁶, and we used N = 64 mesh points for each direction. These simulations are reported in Table 1 as Run 1 and Run 5–13. For each run, we calculated the dissipative time t_d as the time of maximum W. In Fig. 6, we report the dissipative time as a function of the aspect ratio Λ (dotted line). It appears that the dissipative time follows the law t_d ∝ Λ^4/5 (dotted-dashed line) reasonably well, in accordance with what is predicted by Eq. (33).

Fig. 6. Dissipation time td for different aspect ratios, namely Λ ∈ [1, 10]. The dash-dotted line represents a power law with slope 4/5, according to our estimation in Eq. (33). These simulations are reported in Table 1 as Run 1 and Runs 5–13.

A closer examination of Fig. 6 reveals that the dissipative time increases with Λ slightly faster than the prediction t_d ∝ Λ^4/5. Such a small discrepancy can be explained by taking into account that Eq. (33) has been derived assuming a radial profile of the Alfvén velocity that remains unchanged during the wave evolution. Instead, dissipation also affects the equilibrium magnetic field slightly reducing the gradient of the Alfvén velocity in time, as we explicitly verified. According to Eq. (33), a decrease of dc_A/dr corresponds to an increase of t_d. This effect is more pronounced when t_d is larger (i.e., for larger Λ) since, in this case, dissipation has more time to smooth the profile of c_A(r). This corresponds to a deviation from the scaling law t_d ∝ Λ^4/5 that is larger for larger Λ.

Summarizing, the above results show that, in the configuration considered in our model, the time evolution of the torsional Alfvén wave is fully dominated by phase mixing, and the corresponding properties are well reproduced by the numerical model. On the contrary, no relevant nonlinear effects have been observed on transverse fluctuations characterizing this wave, even for relatively large values of the wave amplitude (A = 0.2).

3.2. Turbulence in homogeneous background

In this section, we describe the time evolution of the turbulent perturbation, defined in Sect. 2.5, in a homogeneous background (Run 14). Though the main focus of this paper is on inhomogeneous structures, here, the main purpose is to fix a reference case. In Run 14, we used a low value of the hyper-dissipative coefficient (η = 10⁻⁸) and a high spatial resolution in the perpendicular direction, trying to obtain a spectrum with a wide range of scales compatible with numerical limitations. The results of Run 14 will be used in Sect. 3.3 for a comparison with the behavior of waves and turbulence in the inhomogeneous equilibrium (Runs 15–20). Therefore, in Run 14 we use the same values for the parameters A_B = A_v, Λ, and η as in Runs 15–20. In particular, we set a low value for the fluctuation amplitude (δB/B₀ ≃ A = A_B = A_v = 0.01). Such a choice will allow us to satisfy the condition (7) for Runs 15–20 of Sect. 3.3. Of course, lower values for A would better fulfill condition (7), but this would be exceedingly costly due to long simulation times. We also take into account that a situation characterized by small fluctuation amplitude δB/B₀ and low dissipation can be considered to be more representative of the corona.

In Run 14, physical quantities have been initialized in the following way: B_|| = B₀ = 1, ρ = ρ₀ ≃ 1.21, which is the average value of ρ₀ in the inhomogeneous cases, and c_A ≃ 0.91. Moreover, ${\mathit{B}}_{\perp }=A_B\delta {\mathit{B}}_{\perp }^{\text{turb}}$, ${\mathit{v}}_{\perp }=A_{\text{v}}\delta {\mathit{v}}_{\perp }^{\text{turb}}$ and v_|| = 0, where the turbulent perturbations $\delta {\mathit{B}}_{\perp }^{\text{turb}}$ and $\delta {\mathit{v}}_{\perp }^{\text{turb}}$ have been specified in Sect. 2.5. A 3D visualization of B_y corresponding to the above initial condition is plotted in Fig. 4, upper panel. It can be observed that the typical transverse scale associated with fluctuations is ℓ_⊥0 ∼ 1 (in code units).

In Fig. 7, the quantity W is plotted as a function of time for Run 14. The value t_d ≃ 600 can be derived for the dissipative time in this case. Using such a value in Eq. (5), as well as c_A0 = 1, ℓ_⊥0 ∼ 1 and δB(ℓ_⊥0)/B₀ = 0.01, we obtain Γ ≃ 6. This value will be used in Sect. 3.3 to identify the regime of turbulence in an inhomogeneous background with respect to the condition (7).

Fig. 7. Time history of W for the turbulence in a homogeneous background with low amplitude A and low hyper-dissipative coefficient η (Run 14).

In Fig. 4, lower panel, the B_y component is represented at the time t = t_d. Comparing with the image at the initial time (upper panel), the presence of structures at small scales is clearly visible at t = t_d. Moreover, such structures are dominated by perpendicular wavevectors, as is expected for turbulence with a strong background magnetic field.

3.3. Phase mixing and turbulence in an inhomogeneous background

In this section, we aim to investigate the synergy between phase mixing and nonlinear effects in the generation of small scales within an inhomogeneous equilibrium, like the flux tube we considered. As previously seen, the evolution of the torsional Alfvén wave can be completely described in terms of phase mixing. However, a fluctuation exactly corresponding to an Alfvén torsional wave is unlikely to be present in a loop. More realistically, one can think of a torsional wave that is more or less distorted. Such a distortion can be represented as a superposed turbulent component on the wave. In this case, important points are investigating whether phase mixing is still effective for a distorted torsional wave and how dissipation depends on the relative amplitude of the turbulent component with respect to the wave. Moreover, considering the turbulent component itself, the inhomogeneity of the background structure acting on it could produce small scales through a mechanism similar to phase mixing. In this case, both phase mixing and nonlinear cascade should contribute to the generation of small scales. This combined mechanism should work even in a case where the amplitude of the turbulent component is comparable with or larger than the torsional wave.

The above speculations have been investigated by considering a perturbation given by the superposition of a turbulent fluctuation on a torsional Alfvén wave and making the whole perturbation evolve in the inhomogeneous equilibrium. Therefore, the initial condition is given by

$$\begin{array}{cc}\hfill B_{\Vert }& =B_0(x,y),\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill {\mathit{B}}_{\perp }& =A_B\delta {\mathit{B}}_{\perp }^{\mathrm{WT}}(x,y,z)\hfill \\ \hfill & =A_B[(1-\alpha )\delta {\mathit{B}}_{\perp }^{\mathrm{Aw}}(x,y,z)+\alpha \delta {\mathit{B}}_{\perp }^{\mathrm{turb}}(x,y,z)],\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill v_{\Vert }& =0,\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill {\mathit{v}}_{\perp }& =A_{\mathrm{v}}\delta {\mathit{v}}_{\perp }^{\mathrm{WT}}(x,y,z)\hfill \\ \hfill & =A_{\mathrm{v}}[(1-\alpha )\delta {\mathit{v}}_{\perp }^{\mathrm{Aw}}(x,y,z)+\alpha \delta {\mathit{v}}_{\perp }^{\mathrm{turb}}(x,y,z)],\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill \rho & ={\rho }_0(x,y),\hfill \end{array}$$(39)

where α ∈ [0, 1] is a free parameter that allows us to tune the amount of turbulence superposed on the torsional wave. Here, WT stands for wave+turbulence. For each value of α, magnetic field $\delta {\mathit{B}}_{\perp }^{\text{WT}}$ and velocity $\delta {\mathit{v}}_{\perp }^{\text{WT}}$ perturbations are normalized to the rms value of $\delta {\mathit{B}}_{\perp }^{\text{WT}}$. The case α = 0 corresponds to a purely torsional Alfvén wave while increasing α the turbulence level increases. We performed direct numerical tests by varying α in two cases: (i) propagating Alfvén wave plus turbulence and (ii) standing Alfvén wave plus turbulence. The forms used for propagating and standing Alfvén waves have been specified in Sect. 2.4. For all the simulations described in the Sects. 3.3.1 and 3.3.2 we used a low perturbation amplitude. Specifically, in the case (i) (Runs 15–20) we used A_B = A_v = 0.01. In the case (ii) (Runs 21–28), the standing Alfvénic perturbation is obtained by setting a null initial velocity perturbation A_v = 0. On the other hand, for the sake of comparison with case (i), it is appropriate that the two kinds of perturbations (propagating and standing) have the same amount of fluctuating energy. Therefore, in Runs 21–26 the amplitude A_B is multiplied by a factor $\sqrt{2}$ with respect to the value used in Runs 15–20. For all those runs we used an aspect ratio Λ = 4, a low hyper-dissipative coefficient η = 10⁻⁸, and N_⊥ = 512, N_∥ = 16 mesh points. As already remarked, these values are the same as those used in the low dissipative, low-amplitude Run 14, where we considered only a turbulent perturbation in a homogeneous background.

An important point is checking whether the regime of parameters we are considering (for the turbulence-dominated case α ∼ 1) is consistent with the inequality (7). As discussed in Sect. 2.1, if the inequality (7) is satisfied by the initial condition, then phase mixing dominates during the initial stage and nonlinear cascade during a later stage. In the present case we estimate terms in the inequality (7) by using the following values: dc_A/dr ≃ 0.3 (Fig. 1); c_A ≃ 1; L_|| = ΛL = 8π; ℓ_⊥0 ∼ 1 and Γ ≃ 6 (Sect. 3.2). This gives the value ≃0.07 for the R.H.S of Eq. (7), while the L.H.S. is δB/B₀ = A = 0.01. Therefore, the choice of parameters corresponds to a regime identified by the condition (7), when α ∼ 1. The same regime holds even more for lower values of α, since in cases when the torsional wave has a larger weight with respect to the turbulent component, the nonlinear cascade should play a less relevant role with respect to phase mixing.

3.3.1. Propagating Alfvén wave plus turbulence

The case of a single propagating torsional Alfvén wave plus turbulence has been considered in Runs 15–20, where the value of the parameter α has been varied between α = 0 and α = 1, thus changing the relative weight of the two components.

The time history of the dissipation rate W for these cases is reported in the top panel of Fig. 8. Each curve corresponds to a different value of α. It can be observed that the dissipative time t_d decreases with increasing α. In particular, t_d is reduced by a factor ≃2 going from α = 0 (torsional wave only) to α = 1 (turbulence only). As expected, the presence of a turbulent component superposed on the torsional Alfvén wave does not inhibit the formation of small scales. On the contrary, the larger the turbulent component the more efficient the dissipation is. As we shall see, the Alfvén wave is still subject to phase mixing even when it is distorted by a turbulent component.

Fig. 8. Time history of W for the case of a single Alfvén wave (top panel for Runs 15–20) and two Alfvén waves (bottom panel for Runs 21–26), and different values of α. The dissipative time decreases for increasing α.

What is more, the fastest dissipation is found when the turbulent component dominates (α ≃ 1). In that case the perturbation simultaneously undergoes phase mixing and nonlinear cascade. The combined effect of these two mechanisms leads to the most efficient dissipation. The case of Run 20, corresponding to α = 1 (Fig. 8, upper panel), can be compared with that of a turbulent perturbation in a homogeneous background (Run 14), which is illustrated in Fig. 7. Both those cases correspond to a perturbation with the same amplitude and the same hyper-dissipative coefficient, the only difference being in the background (inhomogeneous in Run 20 and homogeneous in Run 14). A comparison shows that the presence of inhomogeneity in the background structure strongly reduces the dissipative time. However, as we will see in Sect. 3.4, when turbulence evolves in the inhomogeneous background, phase-mixing develops, localized in high density gradient regions. This mechanism plays a relevant role in the generation of small scales and dissipation and clearly illustrates that phase mixing and nonlinear effects jointly act to promote dissipation.

A further effect that accelerates dissipation is related to modifications of the background density in the transverse directions induced by the velocity field δv_⊥ associated with fluctuations. In the case of a pure torsional wave it is $\delta {\mathit{v}}_{\perp }=\delta v_{\perp }(r,z)\widehat{\mathit{\theta }}$; moreover, the background density is initially is independent of θ: ρ₀ = ρ₀(r). In this case it is ∇ ⋅ (ρ₀δv_⊥) = 0 and Eq. (8) implies ∂ρ₀/∂t = 0. Therefore, velocity fluctuations do not modify the background density, and the cylindrical symmetry of the equilibrium is preserved. In contrast, when the turbulent component is present in the velocity fluctuations (α ≠ 0) the cylindrical symmetry is violated and the velocity δv_⊥ can modify the background density structure. In spite of the small amplitude of fluctuations (A_B = A_v = 0.01), this effect accumulates, and deformations in the density structure become more and more relevant with increasing time. Notice that such a modification of the background density can happen even in the case of an incompressible velocity field. In fact, in the incompressible limit ∇ ⋅ v = 0 the continuity equation reads dρ/dt = ∂ρ/∂t + (v ⋅ ∇)ρ = 0, which implies that fluid elements with different density are advected by the velocity field.

In Fig. 9 we report a 3D visualization of the density field at times t = 200 (top), t = 450 (center), and t = 700 (bottom), for the Run 20 where only the turbulent component is present in the fluctuations (α = 1). It can be seen how the initial cylindrical symmetry is gradually lost, and increasingly larger deformations are produced in the density until, at later times, the monolithic structure of the loop is completely destroyed. Of course, this effect is stronger when the turbulent component dominates fluctuations (α ∼ 1). Such effect is similar to as that reported in Magyar et al. (2017), where density initial inhomogenities are deformed by wave propagation.

Fig. 9. 3D visualization for the density field at time t = 200 (top), t = 450 (center), and t = 700 (bottom), for a case of a turbulent perturbation propagating in the inhomogeneous equilibrium (α = 1, Run 20).

Another feature visible in Fig. 9 is the gradual formation of small scales in the background density, transverse to B₀, which is driven by small-scale generation in the velocity perturbation δv_⊥. Small scales in the density correspond to small scales in the spatial distribution of the Alfvén velocity c_A(x, y). This leads to the formation of regions where the Alfvén velocity gradient is much larger than in the initial equilibrium. Since the phase-mixing dynamical time is proportional to |∇c_A|⁻¹ (Eq. (2)), in those regions, the generation of small scales proceeds with a rate larger than at the initial time. This effect, which becomes more relevant with increasing α, contributes to speeding up the dissipation.

However, deformations in the density structure are also observed for lower values of α, though smaller than for α = 1. We report a 2D section of the density field for the propagating wave (left panel of Fig. 10) at z = L_∥/2 and time t = 300. This figure corresponds to Run 16, where we used the value α = 0.23. In this case, the percentage of turbulence in the initial perturbation is not enough to completely destroy the loop, but it is able to generate relevant changes in the morphology of the background density. In the right panel of the same figure, we report the B_y field component at the same time and at the same value of z. The simulation illustrated in Fig. 10 corresponds to a case of a torsional Alfvén wave plus a moderate level of turbulence, and can be then compared with Fig. 3 where the case α = 0 (Run 15) is illustrated. In Fig. 3 it is visible how phase mixing produces wavefronts in the form of perfect concentric shells. Differently, for α = 0.23 (Fig. 10) the cylindrical symmetry is lost; nevertheless, phase mixing still takes place, locally generating small scales in the form of closely packed wavefronts. The wave vector associated with those wavefronts is locally directed parallel to the background density gradient. This shows that, even in the presence of deformations in the density and wave patterns, phase mixing proceeds with properties similar to what found in the cylindrically symmetric configuration and that the presence of turbulence makes dissipation more efficient. This is in accordance with other works, in which resonant absorption is observed (Pascoe et al. 2011). In particular, when density structures are present, mode coupling in inhomogeneous regions leads to faster dissipation, and to a transfer of energy from kink to Alfvén modes (Terradas et al. 2008b).

Fig. 10. 2D perpendicular section of the density field (left panel) and the By component (right panel) for Run 16 at t ∼ td = 300 and z = L∥/2. The simulation corresponds to a propagating Alfvén wave with a moderate level of turbulence superposed (α = 0.23).

3.3.2. Standing Alfvén wave plus turbulence

A coronal loop with a finite length L_|| supports standing waves. Torsional Alfvén waves are among the possible standing modes. Therefore, in addition to the previous ones, we also performed simulations of a standing torsional Alfvén wave (formed by two counter-propagating waves) with a variable amount of turbulence superposed on it. In order to make a comparison, we adopted the same physical and numerical parameters used in the runs of the propagating waves, except for the amplitudes A_B and A_v. Different tests have been performed with various values of the parameter α; those runs are indicated as Run 21–26 in Table 1. Nevertheless, we indicate this same run with two different numbers to more easily indicate how results vary when varying α in the propagating and standing wave cases.

The time evolution of the W is illustrated in the bottom panel of Fig. 8 for this set of runs. Similar as in the propagating wave case, in the present situation the dissipative time decreases with increasing the relative turbulence level (increasing α). At the same time, the maximum value of W decreases with increasing α. For α ≳ 0.66 the time evolution of W becomes essentially independent of α. Comparing the two panels of Fig. 8, it appears that the time behavior of the dissipation rate W in the standing wave cases is qualitatively similar as for the propagating wave, though with some small differences.

In Fig. 11, we report a 2D section illustrating the density profile (left panel) and the B_y component of the magnetic field (right panel) at z = L_∥/2 and t = 300, for Run 22 where a moderate amount of turbulence (α = 0.23) is initially present. This can be compared with the corresponding Run 15 (Fig. 10) relative to a propagating wave case. The comparison shows net differences in morphology. In particular, in the standing wave case, the flux tube boundary appears to be fringed with curly filaments. Such small-scale structures are generated by rolls related to the formation of a KH instability that develops at the loop boundary as a consequence of phase mixing of the standing Alfvén wave. A similar phenomenon also takes place when an initial standing kink mode couples with Alfvénic oscillations that undergo phase mixing (see, e.g., Antolin et al. 2016). In our case, where a torsional Alfvénic standing wave is considered, a certain amount of turbulence (α ≠ 0) must be present; otherwise, the KH instability is not observed.

Fig. 11. 2D perpendicular section of the density field (left panel) and the By component (right panel) for Run 22 at t ∼ td = 300 and z = L∥/2. Here we evolve a standing Alfvén wave with a moderate level of turbulence superposed (α = 0.23).

Similar as in the propagating wave case, we observe the formation of packed wavefronts in the B_y magnetic field, located across the loop boundary. Such wavefronts appear to be distorted by the underlying irregular shape of the background density (left panel of the same figure). Comparing with the corresponding pattern obtained in Run 16 at the same time, we see that, in the propagating wave case, B_y contains scales that are on average smaller than in the standing wave case. This is in accordance with the profiles of W shown in Fig. 8 for α = 0.23 in the two cases. In fact, for t ∼ t_d the dissipation rate W for the propagating wave is slightly larger than for the standing wave (orange lines). In both cases, it can be seen that small scales generated by phase mixing are localized across the region of the density gradient, which corresponds to the location of the strongest Alfvén velocity gradient.

To complement the comparison between the cases of propagating and standing waves, we computed the perpendicular power spectra of the magnetic field transverse components, namely $P(k_{\perp })=|\widehat{B_x}(k_{\perp }){|}^2+{|\widehat{B_{\mathit{y}}}(k_{\perp })|}^2$, where the hat refers to the Fourier coefficients of the fields. Such spectra are calculated by integrating the 3D Fourier spectra both on concentric 2D shells located perpendicularly to the k_z direction and along k_z. The result is a spectrum that depends only on k_⊥. In Fig. 12, we report the magnetic field power spectra for Run 16 and for Run 22 as a function of k_⊥ and t ∼ t_d = 300. A spectrum proportional to $k_{\perp }^{-5/3}$ (black dashed line) is plotted for reference. It can be seen that the two spectra at α = 0.23 are similar. These exhibit a shallower spectrum with respect to the Kolmogorov prediction $\propto k_{\perp }^{-5/3}$, and scale with $\sim k_{\perp }^{-1}.$ Indeed, for such a moderate value of α, phase mixing is mainly responsible for transferring energy toward smaller scales with respect to turbulent couplings. We notice that other authors have found different slopes for the standing wave case (see, e.g., Díaz-Suárez & Soler 2021, 2022), probably depending on the KH instability, which can enhance the development of small scales, driving the flux tube to a turbulent state. Also, differences between the numerical setups can play a role, such as the perturbation amplitude that is ten times higher in Díaz-Suárez & Soler (2021, 2022) with respect to our case.

Fig. 12. $B_{\perp }^2(k_{\perp })$ power spectra for different cases: the propagating Alfvén wave (Run 16), the standing Alfvén wave (Run 22), the turbulence in homogeneous background (Run 14), and the turbulence in inhomogeneous background (Run 20). All the spectra are reported at the dissipative time of the corresponding run. We also report a k−5/3 slope for comparison.

While in a turbulence the amplitude of fluctuations decreases with increasing the wavevector, phase mixing increases the perturbation wavevector keeping the amplitude unchanged, at least in the region where the background inhomogeneity is located. As a result, the Fourier spectra calculated over the entire spatial domain for Runs 16 and 22 are shallower than a Kolmogorov spectrum. At the smallest scales (k_⊥ ≳ 30), the spectrum of the standing wave case have more energy than for the propagating wave. However, these differences are at energies too small to be appreciated in the perturbation morphology.

To conclude, from the point of view of wave dissipation, there are no relevant differences between the propagating and standing waves. In particular, both the KH rolls, which we find in the standing wave, and the distorted torsional waves we are considering, lead to similar dissipation levels, even through different mechanisms.

3.4. Turbulence in inhomogeneous background

The case where the initial perturbation is formed only by the turbulent component (α = 1) has been considered in Run 20 (A_B = A_v = 0.01) and Run 26 ($A_B=\sqrt{2}\times 0.01$, A_v = 0). Results relative to Run 20 are illustrated in Fig. 13, where the density ρ and the B_y component of the magnetic field are reported, at t = t_d = 200, similar as in Figs. 10 and 11. Even in this case, the pattern of B_y(x, y) is reminiscent of what is generated by phase mixing, namely, small-scale structures appear, formed by closely packed wavefronts. Those structures are essentially localized in regions of inhomogeneous density. Therefore, the coupling between the turbulent perturbation and the inhomogeneous background leads to a phenomenon that is qualitatively similar to phase mixing.

Fig. 13. xy−plane for the density field (left panel) and the By component (right panel) for Run 20 at t ∼ td = 200 and z = L∥/2. Here we used the maximum level of relative turbulence (α = 1.0).

However, in Fig. 13 some regions are visible (see, e.g., the area 1 ≲ x ≲ 3, 1 ≲ y ≲ 3) where the wavevector distribution appears to be more isotropic, also in comparison with cases characterized by lower values of α (e.g., α = 0.23, Fig. 10). Those are the regions where the smallest scales are localized. This feature can be interpreted as a manifestation of a nonlinear cascade that forms locally following the development of phase mixing. We also observe that in those regions, the density structure is particularly distorted and presents a large gradient (Fig. 13, left panel). This contributes to the formation of small scales in the perturbation.

To better investigate the role played by phase mixing in the dissipation of the turbulent perturbation we compared Run 26 with Runs 27–28 where A_v = 0,  α = 1, while the amplitude of the initial magnetic perturbation is decreased: $A_B=0.01\times \sqrt{2}$ (Run 26), $A_B=0.005\times \sqrt{2}$ (Run 27) and $A_B=0.0023\times \sqrt{2}$ (Run 28). For those runs, the time dependence of the normalized dissipation rate $W/A_B^2$ is plotted in Fig. 14. It appears that the dissipative time t_d ∼ 150 is essentially independent of the perturbation amplitude A_B. This feature is typical of phase mixing and it had been also found in the case of a pure torsional wave (Runs 1–4, Fig. 5). This confirms that, even in the case of turbulent initial perturbation, phase mixing represents the main mechanism responsible for dissipation, at least for the amplitudes considered in Runs 26–28. However, at odds with Runs 1–4 where the curves of W/A² are almost perfectly superposed (Fig. 5), in the cases illustrated in Fig. 14 the three plots progressively depart from one-another as the amplitude is increased. Therefore, in the case of a moderate-amplitude turbulent perturbation ($A_B=0.01\times \sqrt{2}$, Run 26), other effects play a role beside the dominant phase mixing. We can identify such effects as the nonlinear cascade, which works together with phase mixing in promoting the dissipation of perturbation.

Fig. 14. Time history of W for the case of turbulence in a non-homogeneous background. The three curves refer to Run 26 ($A_B=0.01\times \sqrt{2})$, Run 27 $(A_B=0.005\times \sqrt{2})$, and Run 28 ($A_B=0.0023\times \sqrt{2})$.

A further characterization can be done by considering the perpendicular spectra of fluctuations. In Fig. 12 the spectrum P(k_⊥) of transverse magnetic fluctuations B_⊥ is plotted for the case of turbulence in the flux tube (Run 20, green curve), as well as in the case of turbulence in an homogeneous background (Run 14, black curve). All the spectra plotted in Fig. 12 are calculated at the dissipative time of the corresponding run, namely the time at which W is maximum, in order to have a spectrum as developed as possible. In the case of Run 14 we observe the presence of an inertial range that approximately follows the slope ∝k^−5/3, indicating that the spectrum has been formed by the nonlinear cascade. In contrast, in the case of turbulence in the flux tube (Run 20) a shallower spectrum is observed that is close to the spectra of phase-mixing-dominated Runs 16 and 22. This further demonstrates that, when the turbulence evolves in the inhomogeneous flux tube, phase-mixing plays a relevant role in the generation of small scales and dissipation. At small scales k_⊥ > 10 − 20 the spectrum of Run 20 becomes less steep than those of Runs 16 and 22, causing to a more efficient dissipation. This feature could be interpreted as a sign of other phenomena (nonlinear cascade, deformation of the background density) that cooperate with phase mixing in the generation of small scales.

4. Summary and conclusions

In the present paper, we investigated the interplay between phase mixing and turbulence in generating small scales and consequent wave dissipation in a simple model of a coronal loop. Both effects tend to generate small scales transverse to the background magnetic field in Alfvénic fluctuations, and therefore, they can jointly act to promote the dissipation of fluctuations. We considered a configuration where an initial torsional Alfvén wave and a variable amount of turbulent Alfvénic transverse fluctuations are present. We chose a torsional wave because its time evolution is essentially determined by phase mixing, which is one of the two mechanisms we investigated. There is observational evidence of torsional motions, both in the photosphere (Brandt et al. 1988; Bonet et al. 2008, 2010) and in the chromosphere (Wedemeyer-Böhm 2009; Wedemeyer-Böhm et al. 2012; Tziotziou et al. 2018; Dakanalis et al. 2022), which could induce torsional waves in coronal loops (Jess et al. 2009).

In this study, we used the COHMPA algorithm, which solves the 3D compressible MHD equations with periodic boundary conditions. We considered an equilibrium structure with a low β value, representing an overdense cylindrical magnetic flux tube. The Alfvén speed gradient across the flux tube boundary is responsible for the phase mixing of Alfvénic perturbations. We characterized the time evolution of a pure torsional Alfvén wave by varying parameters such as wave amplitude and parallel wavelength. We verified that the properties of the dissipation rate in the numerical runs follow what can be predicted for the phase mixing to a good extent. However, the dynamical time associated with phase mixing increases with the decreasing spatial scale ℓ_⊥ of perturbations, and it becomes exceedingly long for realistic values of the dissipative scale in the coronal plasma. In this respect, the turbulent cascade has the advantage of moving the fluctuation energy down to small dissipative scales in a time that is on the order of the large-scale nonlinear time. For that reason MHD turbulence has been considered in several coronal heating models (Nigro et al. 2004; Malara et al. 2010; Van Ballegooijen et al. 2017; Rappazzo et al. 2017; Van Ballegooijen & Asgari-Targhi 2018). In addition, there are indications of the presence of turbulent oscillations in the corona (Banerjee et al. 1998; Singh et al. 2006; Hahn & Savin 2013, 2014; Morton et al. 2016, 2019)

A torsional wave propagating in a cylindrical flux tube represents a peculiar configuration that is completely independent of the azimuthal coordinate θ. It is natural to ask to what extent the phenomenon of phase mixing is related to such a particular symmetry and whether phase mixing is still present when the cylindrical symmetry is violated. In our simulations, we investigated this question by distorting the pure torsional wave with the addition of a certain amount of turbulent perturbation. Our results have shown that, for moderate amplitudes of the turbulent component (small α), phase mixing persists even when the cylindrical symmetry is not preserved: closely packed wave fronts form, with a wave vector that is locally parallel to the density gradient in planes perpendicular to B₀ (Fig. 10). In this respect, phase mixing appears to be a robust mechanism active in configurations more complex than was originally conceived (Heyvaerts & Priest 1983).

We also considered the case of standing Alfvén waves with a variable amount of turbulence superposed on it. In this case, the evolution of the background structure suggests that rolls form at the flux tube boundary, probably generated by a mechanism similar to a KH instability. This phenomenon was previously observed in another situation, where an initial kink mode of an overdense flux tube couples with azimuthal Alfvénic perturbations which in turn undergo phase mixing (e.g., Antolin et al. 2016). In our case, where the initial perturbation is a torsional Alfvén wave, this phenomenon requires a certain amount of turbulent distortion in the initial perturbation (α ≠ 0), otherwise the cylindrical symmetry prevents the formation of the KH instability. Moreover, KH rolls are not observed in the case of the propagating wave, even in the presence of a turbulent component. The velocity perturbation pattern, with its crests and troughs, must probably remain stationary to allow for the development of the KH instability. However, comparing the propagating and standing waves, we did not find relevant differences in the dissipative time for the two cases (Fig. 8). Therefore, the presence of KH rolls seems not to enhance the efficiency of small-scale formation, at least for the considered configuration.

Increasing the parameter α up to α = 1, the initial perturbation becomes dominated by the turbulent component. The case α = 1 corresponds to the shortest dissipative time. Even in the case where no torsional wave is initially present, the time evolution of the turbulent perturbation reveals features qualitatively similar to what is produced by phase mixing, namely closely packed wavefronts localized at the flux tube boundary, with a wave vector locally parallel to the density gradient (Fig. 13). This indicates that a mechanism similar to phase mixing acts, even with a turbulent Alfvénic perturbation, driven by transverse inhomogeneities in the Alfvén velocity.

The case of a turbulent perturbation evolving in a flux tube was studied in a regime expressed by Eq. (7). In this case, the amplitude of the initial turbulent perturbation was small enough to make the nonlinear time, calculated for the large, energy-containing scales, longer than the corresponding phase-mixing dynamical time. In these conditions, the small-scale formation within the turbulent perturbation is initially dominated by phase mixing rather than by the nonlinear cascade. However, while phase mixing proceeds, the perturbation energy moves toward small scales, thus reducing the effective nonlinear time. Eventually, the nonlinear cascade becomes faster than phase mixing and builds up small scales in a time that is essentially independent of the dissipative coefficients. This aspect is relevant for the coronal plasma, where such coefficients are expected to be very small. In the considered regime, our simulations have shown that the corresponding dissipative time is shorter than is found in a pure phase mixing case (represented by the case α = 0) and for a standard turbulent cascade (Run 14). Therefore, phase mixing and turbulence work in a synergistic way, speeding up the dissipation of the perturbation energy.

We refer to the above situation as a phase mixing + nonlinear cascade regime. The opposite case corresponds to a perturbation amplitude large enough to violate condition (7). In this case, the nonlinear time is shorter at large scales than the phase-mixing dynamical time. Therefore, the nonlinear cascade dominates the generation of small scales, with phase mixing playing a minor role. We refer to this latter situation as a nonlinear cascade regime. To better illustrate these two regimes, in the condition (7), we express the gradient of the Alfvén velocity as dc_A/dr ∼ Δc_A/Δr, where Δc_A = c_A, ext − c_A, int is the variation of the Alfvén velocity across the flux tube boundary, which has a width Δr, and $c_{\mathrm{A},\mathrm{ext}}=B_{\mathrm{ext}}/\sqrt{4\pi {\rho }_{\mathrm{ext}}}$ and $c_{\mathrm{A},\mathrm{int}}=B_{\mathrm{int}}/\sqrt{4\pi {\rho }_{\mathrm{int}}}$. Similarly, we estimate the Alfvén velocity as c_A ∼ (c_A, ext + c_A, int)/2. For small values of the plasma β, we can assume that B_ext ≃ B_int ≃ B₀. Therefore, it is

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{c}_{\mathrm{A}}}{c_{\mathrm{A}}}\sim 2\frac{\sqrt{{\rho }_{\mathrm{int}}/{\rho }_{\mathrm{ext}}}-1}{\sqrt{{\rho }_{\mathrm{int}}/{\rho }_{\mathrm{ext}}}+1.}\end{array}$$(40)

We can re-express the condition (7) in terms of the dimensionless ratios ρ_int/ρ_ext, L_||/ℓ_⊥0, ℓ_⊥0/Δr and δB/B₀ characterizing the problem in the following form:

$$\begin{array}{c}\hfill F(\frac{{\rho }_{\mathrm{int}}}{{\rho }_{\mathrm{ext}}},\frac{L_{||}}{{\ell }_{\perp 0}},\frac{\delta B}{B_0},\frac{l_{\perp 0}}{\mathrm{\Delta }r})\equiv \frac{L_{||}}{{\ell }_{\perp 0}}-2\mathrm{\Gamma }\frac{(\sqrt{{\displaystyle \frac{{\rho }_{\mathrm{int}}}{{\rho }_{\mathrm{ext}}}}}-1)}{(\sqrt{{\displaystyle \frac{{\rho }_{\mathrm{int}}}{{\rho }_{\mathrm{ext}}}}}+1)}\frac{\left({\displaystyle \frac{{\ell }_{\perp 0}}{\mathrm{\Delta }r}}\right)}{\left({\displaystyle \frac{\delta B}{B_0}}\right)}<0.\end{array}$$(41)

Condition (41) is represented in Fig. 15 where lines corresponding to F = 0 are plotted in the (ρ_int/ρ_ext, L_||/ℓ_⊥0)-plane, for three values of the ratio δB/B₀ (δB/B₀ = 0.01, 0.02, 0.03), and choosing ℓ_⊥0/Δr = 1 (as in our simulations). The portion of the plane below each curve corresponds to configurations satisfying condition (41) (or, equivalently, condition (7)). Therefore, it corresponds to the regime denoted as phase mixing + nonlinear cascade. This is verified for high values of the density contrast ρ_int/ρ_ext and for a longitudinal scale L_|| that is not much larger than the perpendicular scale ℓ_⊥0. The configuration considered in Run 20 is indicated by a black dot, which falls into this regime. The opposite regime (nonlinear cascade) corresponds to low-density contrast and/or high values for the ratio L_||/ℓ_⊥0. In this case, phase mixing is slower than nonlinear couplings, and the process of small-scale generation is dominated by the nonlinear cascade. Since the nonlinear timescale decreases with increasing the perturbation amplitude, the portion of the plane corresponding to this latter regime becomes larger when the ratio δB/B₀ is increased. The above considerations indicate that the effects of phase mixing in the dynamics and dissipation of turbulent perturbations are more relevant in coronal loops with a high-density contrast and short length. A detailed study of wave damping has been done in Morton et al. (2021), where the authors find that higher density contrasts lead to higher observed damping rates of waves in corona. Even though their study focuses on kink waves, we found similar results when comparing simulations of turbulence in a homogeneous (Run 14) and in an inhomogeneous background (Runs 20, 26–28).

Fig. 15. Curves F(ρint/ρext, L||/ℓ⊥0, δB/B0, l⊥0/Δr) = 0 are represented in the (ρint/ρext, L||/ℓ⊥0)-plane, for δB/B0 = 0.01 (purple line), δB/B0 = 0.02 (green line), and δB/B0 = 0.03 (orange line); each curve is calculated for ℓ⊥0/Δr = 1. The portion of the plane below a given curve corresponds to the regime defined by condition (7) (phase mixing + nonlinear cascade regime), while the portion above the curve corresponds to the opposite condition (nonlinear cascade regime). The black dot indicates the configuration considered in Run 20.

Despite the small amplitude of the considered perturbations, their effects on the background structures are nonnegligible when considered for a sufficiently long time. In particular, the initial density structure is gradually distorted, losing its cylindrical symmetry. This effect becomes more relevant at large α, and it can completely destroy the initial structure of the flux tube after a certain time. This leads to the formation of regions where the gradient of the Alfvén velocity is particularly large. In those regions, the generation of small scales in the perturbation proceeds with a rate higher than elsewhere (Fig. 13). This effect further contributes to accelerate dissipation. Therefore, it can be expected that such distortion of the background structure in a coronal loop corresponds to locations of greater heating.

In previous work, it was shown that multiple density patches coupled with a large-scale wave eventually result in a turbulent distribution of the density (Magyar et al. 2017, 2019). Here, we investigated the case of a single inhomogeneity perturbed by a spectrum fluctuations at large scale, focusing our interest on the interplay between generalized phase mixing (following the definition by Magyar et al. 2017) and nonlinear effects.

In summary, phase mixing appears to be a robust mechanism that works both in the less idealized case of a distorted torsional wave and on a fully turbulent perturbation. We conclude that phase mixing, turbulent cascade, and perturbation-induced background distortions jointly work to promote wave dissipation. The present results are relevant for the solar corona, where the plasma is characterized by extremely high Reynolds and Lundquist number values. In that ambient, continuous injection of fluctuations from photospheric motions is expected, which demands forced rather than decaying simulations, such as those presented in this paper. Such a study is left for future work.

Acknowledgments

The present work has been supported by the Italian Space Agency (ASI) within the Project “Partecipazione italiana alla missione NASA/MUSE” (CUP: F83C22001920005), related to the NASA mission “Multi Slit Solar Explorer (MUSE)”. The simulations have been performed at the Newton cluster at University of Calabria and the work is supported by “Progetto STAR 2-PIR01 00008” (Italian Ministry of University and Research). SS acknowledges supercomputing resources and support from ICSC–Centro Nazionale di Ricerca in High-Performance Computing, Big Data and Quantum Computing–and hosting entity, funded by European Union–NextGenerationEU. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high-performance computing resources and support. FP acknowledges support from the Research Foundation - Flanders (FWO) Junior research project on fundamental research - G020224N;OP is partially supported by the PRIN 2022 project 2022KL38BK - The ULtimate fate of TuRbulence from space to laboratory plAsmas (ULTRA) “(CUP B53D23004850006) by the Italian Ministry of University and Research, Finanziato dall’Unione europea – Next Generation EU’ PIANO NAZIONALE DI RIPRESA E RESILIENZA (PNRR) Missione 4 Istruzione e Ricerca” - Componente C2 Investimento 1.1, Fondo per il Programma Nazionale di Ricerca e Progetti di Rilevante Interesse Nazionale (PRIN)’ Settore PE9.

References

All Tables

All Figures

	Fig. 1. 1D profiles of the equilibrium quantities B0, cA, P0, and ρ0 as functions of y, for x = L/2 and z = ΛL/2. In the inset we report a zoom-in for the pressure profile described by Eq. (13) (dashed line) and for the correction described by Eq. (14) (solid line), near the boundary y = 2π. The correction acts in the same way at y = 0 for symmetry.
In the text

	Fig. 2. 1D cut representing the components δBy and δvy as functions of x, for the single Alfvén wave perturbation. The section is computed along a straight line parallel to the x-axis, at y = L/2 and z = ΛL/2.
In the text

	Fig. 3. 2D perpendicular section of the By component for the evolution of Alfvén waves at times t = 0 (left), t = 150 (middle), and t = 300 (right). Here we used a perturbation amplitude A = 10−2, a hyper-dissipation coefficient η = 10−8, and an aspect ratio Λ = 4 (Run 15).
In the text

	Fig. 4. 3D visualization of the By component at times t = 0 (top) and t = td (bottom), for Run 14.
In the text

	Fig. 5. Time history of W/A2 for the propagating Alfvén wave, with different values of amplitude A, for Runs 1–4. For these simulations the perturbation amplitude refers both to the magnetic and the velocity fields, namely A = AB = Av.
In the text

	Fig. 6. Dissipation time td for different aspect ratios, namely Λ ∈ [1, 10]. The dash-dotted line represents a power law with slope 4/5, according to our estimation in Eq. (33). These simulations are reported in Table 1 as Run 1 and Runs 5–13.
In the text

	Fig. 7. Time history of W for the turbulence in a homogeneous background with low amplitude A and low hyper-dissipative coefficient η (Run 14).
In the text

	Fig. 8. Time history of W for the case of a single Alfvén wave (top panel for Runs 15–20) and two Alfvén waves (bottom panel for Runs 21–26), and different values of α. The dissipative time decreases for increasing α.
In the text

	Fig. 9. 3D visualization for the density field at time t = 200 (top), t = 450 (center), and t = 700 (bottom), for a case of a turbulent perturbation propagating in the inhomogeneous equilibrium (α = 1, Run 20).
In the text

	Fig. 10. 2D perpendicular section of the density field (left panel) and the By component (right panel) for Run 16 at t ∼ td = 300 and z = L∥/2. The simulation corresponds to a propagating Alfvén wave with a moderate level of turbulence superposed (α = 0.23).
In the text

	Fig. 11. 2D perpendicular section of the density field (left panel) and the By component (right panel) for Run 22 at t ∼ td = 300 and z = L∥/2. Here we evolve a standing Alfvén wave with a moderate level of turbulence superposed (α = 0.23).
In the text

	Fig. 12. $B_{\perp }^2(k_{\perp })$ power spectra for different cases: the propagating Alfvén wave (Run 16), the standing Alfvén wave (Run 22), the turbulence in homogeneous background (Run 14), and the turbulence in inhomogeneous background (Run 20). All the spectra are reported at the dissipative time of the corresponding run. We also report a k−5/3 slope for comparison.
In the text

	Fig. 13. xy−plane for the density field (left panel) and the By component (right panel) for Run 20 at t ∼ td = 200 and z = L∥/2. Here we used the maximum level of relative turbulence (α = 1.0).
In the text

	Fig. 14. Time history of W for the case of turbulence in a non-homogeneous background. The three curves refer to Run 26 ($A_B=0.01\times \sqrt{2})$, Run 27 $(A_B=0.005\times \sqrt{2})$, and Run 28 ($A_B=0.0023\times \sqrt{2})$.
In the text

Fig. 15. Curves F(ρint/ρext, L||/ℓ⊥0, δB/B0, l⊥0/Δr) = 0 are represented in the (ρint/ρext, L||/ℓ⊥0)-plane, for δB/B0 = 0.01 (purple line), δB/B0 = 0.02 (green line), and δB/B0 = 0.03 (orange line); each curve is calculated for ℓ⊥0/Δr = 1. The portion of the plane below a given curve corresponds to the regime defined by condition (7) (phase mixing + nonlinear cascade regime), while the portion above the curve corresponds to the opposite condition (nonlinear cascade regime). The black dot indicates the configuration considered in Run 20.

In the text