We use the ideal incompressible MHD equations, ∇ ⋅ v ⃗ = 0 , (1) ρ ∂ v ⃗ ∂ t + ρ v ⃗ ⋅ ∇ v ⃗ = − ∇ p + 1 μ ( ∇ × B ⃗ ) × B ⃗ , (2) ∂ B ⃗ ∂ t = ∇ × ( v ⃗ × B ⃗ ) , (3) where v ⃗ denotes the velocity field of the plasma, B ⃗ is the magnetic field, ρ is the density, μ is the magnetic permeability, and p is the gas pressure.

We consider an equilibrium configuration in a Cartesian coordinate system with a uniform background magnetic field which is directed along with the z axis ( B 0 ⃗ = B 0 e ⃗ z ) and no background flow ( v ⃗ 0 = 0 ). Here and hereafter, we chose a piece-wise constant density (Goossens et al., 2012): ρ ( x ) = ρ l if x ≤ 0 , ρ r if x > 0 , (4) where ρ _l and ρ _r are constant and ρ _l ≠ ρ _r . We use l and r indices to distinguish the density values on the left and right sides of the interface. As a next step, the incompressible MHD equations have been linearised ( B ⃗ = B ⃗ 0 + B ⃗ ′ , v ⃗ = 0 + v ⃗ ′ , ρ = ρ 0 + ρ ′ , p = p 0 + p ′ ): ∇ ⋅ v ⃗ = 0 , ρ 0 ∂ v ⃗ ′ ∂ t = − ∇ p ′ + 1 μ ( B ⃗ 0 ⋅ ∇ ) B ⃗ ′ − 1 μ ∇ B ⃗ 0 ⋅ B ⃗ ′ , ∂ B ⃗ ′ ∂ t − ( B ⃗ 0 ⋅ ∇ ) v ⃗ ′ = 0 , We take all perturbed quantities proportional to f′ = f(x)  exp(i (k _y y + k _z z − ω t)). We Fourier decompose the perturbation with respect to y, z and t, as the equilibrium quantities are only depending on x. Here f represents the physical variables, ω is the wave frequency and k _y , k _z are the wavenumbers. As a next step, we rewrite the equation in terms of displacement ( ξ ⃗ = d v ⃗ ′ / d t ) and total pressure (P′ = p′ + B ₀ B _z /μ). Eliminating the variables (ξ _y , ξ _z , ρ ₀, B _x , B _y ), we can combine the linearized MHD equations into one 2nd order differential equation (ODE). Using the x-component of the equation of motion d P ′ d x = ρ 0 ( ω 2 − ω A 2 ) ξ x , (5) where ω _A = k _z V _A and V A = B 0 / μ ρ ( x ) . We get a 2nd order ODE for total pressure as follows ρ 0 ( ω 2 − ω A 2 ) d d x 1 ρ 0 ( ω 2 − ω A 2 ) d P ′ d x = κ 2 P ′ . (6) We introduce the new variable κ 2 = k y 2 + k z 2 , which is the wavenumber in the perpendicular direction (i.e., perpendicular to the density interface). We consider the unidirectional magnetic field along the z direction, which allows us to take different Alfvén speeds on the left and the right side. In other words, surface waves can exist only when V _Al ≠ V _Ar (indices l and r denotes left and right side of the density medium, see Eq. 4). Now due to constant quantities in each half-space Eq. 6 simplifies into d 2 P ′ d x 2 = κ 2 P ′ , (7) The solutions for Eq. 7 are finite at x = ±∞ if and only if P ′ ( x ) = P l ′ = C 1 e κ x , if x ≤ 0 , P r ′ = C 2 e − κ x , if x > 0 , (8) At this point we implement the continuity conditions for total pressure P ′ x = 0 = 0 and displacement ξ x x = 0 = 0 at the interface x = 0. The braces denote the difference across the interface. This resulted in C ₁ = C ₂ = C and the known dispersion relation for surface waves (Ionson, 1978; Wentzel, 1979b; Roberts, 1981) as follows ω 2 = ρ l ω A l 2 + ρ r ω A r 2 ρ l + ρ r , (9) where ω _Al = k _z V _Al and ω _Ar = k _z V _Ar .

Now we compute the Elsässer variables for surface Alfvén waves in our setup. We show that surface Alfvén waves carry both Elsässer variables even if they propagate only in one direction (inward or outward). This is in clear contrast to Alfvén waves in uniform plasmas. Therefore, we describe the Elsässer variables for perturbed quantities. Z ⃗ ± = v ⃗ ′ ± B ⃗ ′ μ ρ 0 , (10) where the Elsässer variables are the well-known combination of flow and Alfvén speed perturbation. In an attempt to calculate the Elsässer variables, we determine the linear variables as follows v ⃗ ′ = d ξ ⃗ d t = 1 ρ 0 ( ω 2 − ω A 2 ) d d t ( ∇ P ′ ) , B ⃗ ′ = B 0 d ξ ⃗ d z = B 0 ρ 0 ( ω 2 − ω A 2 ) d d z ( ∇ P ′ ) . (11) Here frequency omega is ω = ω _k given in Eq. 9. For a propagating wave, we take total pressure as P′ = P′(x) cos k _y y  cos(k _z z − ω t), where we take P′(x) as in Eq. 8. P′ is harmonic in time, in y and z directions, as the κ contains both k _y and k _z . We substitute Eq. 11 into Eq. 10 in computing the Elsässer variables as follows Z x ± = d P ′ ( x ) d x ω ∓ ω A ρ 0 ( ω 2 − ω A 2 ) cos k y y sin ( k z z − ω t ) , (12) Z y ± = − P ′ ( x ) k y ω ∓ ω A ρ 0 ( ω 2 − ω A 2 ) sin k y y sin ( k z z − ω t ) , (13) Z z ± = P ′ ( x ) k z ω ∓ ω A ρ 0 ( ω 2 − ω A 2 ) cos k y y cos ( k z z − ω t ) , (14)

In this section, we derive the evolution of the wave energy densities. With the use of the definition of Elsässer variables, we obtain an expression for wave energy densities. w ± = ρ 0 ( Z ⃗ ± ) 2 4 = ( ω ∓ ω A ) 2 4 ρ 0 ( ω 2 − ω A 2 ) 2 ∂ P ′ ∂ x 2 ⁡ cos 2 k y y + P ′ 2 k y 2 sin 2 k y y sin 2 ( k z z − ω t ) + + k z 2 P ′ 2 ( ω ∓ ω A ) 2 4 ρ 0 ( ω 2 − ω A 2 ) 2 cos 2 ( k y y ) cos 2 ( k z z − ω t ) . (15) This can be simplified for the solutions (Eq. 8) in the two regions w r ± = C 2 e − 2 κ x 4 ρ r ( ω ± ω A r ) 2 k y 2 sin 2 ( k z z − ω t ) + k z 2 cos 2 ( k y y ) , w l ± = C 2 e 2 κ x 4 ρ l ( ω ± ω A l ) 2 k y 2 sin 2 ( k z z − ω t ) + k z 2 cos 2 ( k y y ) , (16) We calculate the energy density averaged over the cross-section and over the period ⟨ w ⟩ = ∫ x d x k y 2 π ∫ 0 2 π / k y d y ω 2 π ∫ 0 2 π / ω w d t , When we average over the x direction, we use the variable as a predetermined length L in the transverse direction. The integration in the x direction is from − L to + L 〈 w r 〉 = 〈 w r + + w r − 〉 = C 2 κ 8 ( ω 2 + ω A r 2 ) ρ r ( ω 2 − ω A r 2 ) 2 ( 1 − e − 2 κ L ) , 〈 w l 〉 = 〈 w l + + w l − 〉 = C 2 κ 8 ( ω 2 + ω A l 2 ) ρ l ( ω 2 − ω A l 2 ) 2 ( 1 − e − 2 κ L ) , After taking the average wave energy density for two different regions, we take the sum of ⟨w⟩ = ⟨w _l ⟩ + ⟨w _r ⟩. We obtain the energy density averaged over the cross-section as follow ⟨ w ⟩ = C 2 κ 8 ρ r ( ω 2 − ω A r 2 ) 2 ( ω 2 + ω A l 2 ) + ρ l ( ω 2 − ω A l 2 ) 2 ( ω 2 + ω A r 2 ) ρ l ρ r ( ω 2 − ω A l 2 ) 2 ( ω 2 − ω A r 2 ) 2 ( 1 − e − 2 κ L ) , (17) By using the dispersion relation (Eq. 9) it simplifies to ⟨ w ⟩ = C 2 κ 4 ( ρ r + ρ l ) 3 ω 2 ρ r 2 ρ l 2 ( ω A l 2 − ω A r 2 ) 2 ( 1 − e − 2 κ L ) . (18)

The energy cascade rate is computed from the non-linear term in the incompressible MHD equations. Following Van Doorsselaere et al. (2020a), the energy dissipation rate for upward and downward surface Alfvén waves can be expressed by using. ϵ ± = Z ⃗ ∓ ⋅ ∇ w ± (19) We find that ϵ r ± = C 3 e − 3 κ x cos ( k y y ) sin ( k z z − ω t ) 2 ρ r 2 ( ω 2 − ω A r 2 ) ( ω ± ω A r ) k z 4 cos 2 ( k y y ) + k y 4 sin 2 ( k z z − ω t ) + 2 k y 2 k z 2 , x > 0 , (20) ϵ l ± = C 3 e 3 κ x cos ( k y y ) sin ( k z z − ω t ) 2 ρ l 2 ( ω 2 − ω A l 2 ) ( ω ± ω A l ) k z 4 cos 2 ( k y y ) + k y 4 sin 2 ( k z z − ω t ) + 2 k y 2 k z 2 , x ≤ 0 . (21) Since there are the third power of a periodic function in the y direction and time, taking the average of Eqs 20, 21 leads to an energy cascade rate of 0. In solar wind Alfvén turbulence model, this zero average in the time domain is avoided by approximating Z ⃗ ± ⋅ ∇ w ± to ( Z ⃗ ± ) 2 w ± / L , where L is the correlation length that defines the perpendicular scales of the wave energy distributions. Van Doorsselaere et al. (2020a) generalized this method and used the root-mean-square approach while averaging over a period and the φ direction. In order to avoid an energy cascade equal to zero, we use the root-mean-square approximation while averaging over the y direction and period. ⟨ ϵ ⟩ = ∫ x d x k y 2 π ∫ 0 2 π / k y d y ω 2 π ∫ 0 2 π / ω ϵ 2 ⁡ d t 1 / 2 , (22) Finally, 〈 ϵ r 〉 = 〈 ϵ r + + ϵ r − 〉 = 2 24 ω | C 3 | κ 5 k y 8 + 24 k y 6 k z 2 + 41 k y 4 k z 4 + 24 k y 2 k z 6 + 5 k z 8 ( ω 2 − ω A r 2 ) 2 ρ r 2 ( 1 − e − 3 κ L ) , (23) 〈 ϵ l 〉 = 〈 ϵ l + + ϵ l − 〉 = 2 24 ω | C 3 | κ 5 k y 8 + 24 k y 6 k z 2 + 41 k y 4 k z 4 + 24 k y 2 k z 6 + 5 k z 8 ( ω 2 − ω A l 2 ) 2 ρ l 2 ( 1 − e − 3 κ L ) , (24) Once more we use the dispersion relation to reduce the expression as follows ⟨ ϵ ⟩ = 2 12 ω | C 3 | κ 5 k y 8 + 24 k y 6 k z 2 + 41 k y 4 k z 4 + 24 k y 2 k z 6 + 5 k z 8 ( ρ r + ρ l ) 2 ρ r 2 ρ l 2 ( ω A l 2 − ω A r 2 ) 2 ( 1 − e − 3 κ L ) , (25) Taking the ratio τ = ⟨w⟩/⟨ϵ⟩ and considering L → ∞ we get τ = 3 2 C ω ( ρ r + ρ l ) ( k y 2 + k z 2 ) 5 k y 8 + 24 k y 6 k z 2 + 41 k y 4 k z 4 + 24 k y 2 k z 6 + 5 k z 8 . (26)

We know that in solar corona structure, k _z ≫ k _⊥ (length of the structure is much longer than the length scale in the perpendicular direction), which allows us to acknowledge the long-wavelength limit. In order to simplify the last equation, we consider quasi-perpendicular propagation, also known as the long-wavelength limit. We use this assumption to reduce the complexity of the wavenumber expression. We introduce the new parameter δ = k _z /k _y ≪ 1. Taking the perpendicular component v _y from Eq. 11 yields. v y V A ∼ − C V A r k y ω ρ r ( ω 2 − ω A r 2 ) , if x > 0 Considering last expression and rewriting in terms of the density contrast ζ = ρ _r /ρ _l leads Eq. 26 to τ = 3 10 5 V k y ζ + 1 ζ − 1 , (27) We can see that in an inhomogeneous incompressible MHD, the damping time scale is inversely proportional to the oscillation amplitude V, where for the density contrast, we obtained same relation as was discussed by Van Doorsselaere et al. (2020a). They found a similar expression for damping time in the cylindrical configuration: τ = 2 5 π R V ζ + 1 ζ − 1 , (28) Here V is the velocity amplitude, and R is the radius of the coronal loop, which in planar geometry R = 1/k _y can be used. Remarkably, the only difference is the factor, which is approximately 4 times smaller. We show that the fundamental radial mode of kink waves in a cylinder behaves very much like a surface Alfvén waves.

We will be able to vary in the simulations the parameters involved in this formula for the damping by the energy cascade. These surface Alfvén waves can be easily numerically modeled, and we can directly compare the energy dissipation rate to this analytical work we derived.

3D ideal MHD simulations were run for numerical demonstration of the self-cascading model using the code MPI-AMRVAC (Keppens et al., 2012; Xia et al., 2018). In the simulation, flux schemes were adopted with a five-wave variant, HLLD solver, and Woodward slope limiter. The base grid resolution is chosen as 60 × 60 × 40 in x, y, and z directions, respectively. After applying three-level refinement, the grid resolution becomes 240 × 240 × 160. However, the physical system domain is 0.12 × 0.1 × Z ₀, where Z ₀ has been changed with the density contrast as in Table 1, and the lower density contrast requires additional length along the z direction due to weaker damping. All the results are presented in code units. We establish the following normalization values of unit velocity = 1.16 × 10⁴ cm/s, unit length = 1 cm, number density = 1 cm ⁻³, and unit time = 8.6 × 10⁻⁵ s in our code. For a series of simulations the chosen values for amplitudes are V = {0.01, 0.03, 0.06, 0.12, 0.18}, where ζ = 5 considered and similarly we change the density contrast ζ = {1.1, 2, 5, 10, 20} for V = 0.06.

We chose simple initial conditions of our model. The following quantities are in code units. Firstly, the straight, homogeneous magnetic field value is B ₀ = 1.115 along the z-axis, but this value does not play a significant role. Secondly, we employ a step function of density perpendicular to the magnetic field, which is, in our case, at x = 0, where the density changes discontinuously for ρ _l to ρ _r . In the y and z direction, it does not vary. Moreover, in the simulations, we considered a nearly incompressible regime. The minimum Alfvén speed is V _A ∼ 0.7 and sound speed is V _s ∼ 1.78, resulting in a high plasma β ∼ 4.4. These are significantly greater than our velocity perturbation. Furthermore, in the paper by Magyar et al. (2019b), they considered tests with different plasma beta values (from β ≈ 0.02 to 15), indicating that the dynamics perpendicular to the magnetic field is not very sensitive to its value.

Firstly, we check if the numerical dissipation influences the measured damping time driven by cascade. The higher resolution affects the generation of smaller scales at the interface but leads to almost similar results for the damping time. We run the simulation for V = 0.06 and ζ = 5 case with the 400 × 400 × 240 grid resolution and the damping time is τ ≈ 0.85. The calculated damping time is similar with a lower resolution: τ ≈ 0.85.

We implement the following boundary conditions: First, at the bottom of the z direction, we take a velocity along the x direction, which varies sinusoidally in y. v x ( y , t ) = V ⁡ cos ( ω t ) sin ( k y y ) . According to Magyar et al. (2019b), such a driver will result in uniturbulence in our setup. Here V is the initial velocity amplitude, the frequency of driver is ω = 2 π and the wavenumber chosen for one wavelength in the y direction, k _y = 2 π/0.1.

Second, continuous conditions (Neumann style zero gradient conditions) are taken in the top boundary in the z direction and both x directions for all variables. Third, for the y direction, we used periodic boundary conditions.

Figures 1, 2 1constitute a set of three-dimensional fluctuations of density change for various initial velocity amplitudes (V) and density contrasts (ζ), respectively. We take snapshots of the cross-section at half of the z domain, with the relevant Z ₀ values (Table 1). However, all snapshots were taken at four time period (t = 4.0) except V = 0.18 and ζ = 20. These latter snapshots were taken at t = 2.5 due to the high velocity amplitude or high density contrast numerical simulations becoming unstable for reasons we do not fully understand. In any case, we analyze the first five wavefronts in the remaining text, and the t = 2.5 time is admissible to determine the mean damping time of the first five wavefronts avoiding the instability. In Figure 1 we can see that even at low amplitude, the interface changes its sinusoidal shape and shows effects of non-linearity. Increasing density contrast shows the turbulence and energy cascade loss from large-scale motion to small-scale motion leading to numerical energy dissipation. On the other hand, increasing the amplitude shows that the nonlinear generation of smaller scales is increasing drastically on the density surface, showing evidence of uniturbulence. However, in Figure 2 at ζ = 1.1, we can see that the interface keeps its shape for a more extended time period, showing that uniturbulence is delayed. We can also see that the interface moves with much higher amplitude in higher density contrast regimes due to the significant differences of slow and fast Alfvén speed ( B ⃗ / μ ρ l ≫ B ⃗ / μ ρ r ).

In order to see how the energy transfers to smaller scales in the system, we compute the power spectrum of the velocity field by performing FFT of the v _x component along the y direction (Popescu Braileanu et al., 2021). In particular, we consider the Fourier power at k _y = 2π/0.1 to track the wave’s amplitude in the fundamental mode that was injected by the driver. As a next step, we take the average of the Fourier power in the strip between x = X ₀ and x = − X ₀ as in Figure 3. In other words, we integrate over the volume of the strip region in 2D. We calculate the average at each time moment (t _i , i = 0.1‥4.0), and for every cross-section (x − y plane) of z _j , j = 0‥160 (grid resolution in z). Furthermore, for different simulations, the width X ₀ was taken slightly differently (see Table 1). The width increased or decreased for high or small ζ and V values. However, we check the influence of chosen widths for ζ = 5 and V = 0.06 simulation case in the vertical direction. This is illustrated in Figure 4, showing the damping length of different widths for symmetric strip region (red squares) and non-symmetric strip region (blue dots). In a symmetric case X ₀ = | − X ₀|, but for a non-symmetric case, we take X ₀ and − X ₀ strips close to the small structures. We find that the width X ₀ does not have a large impact on the damping length. Thus we consider only the values as in Table 1.

In the following analysis, we take the averaged power amplitudes of the fundamental mode and follow the peak values with a kink speed (V _k ) in the z direction (z _i = V _k t _i ) at each time step t _i (Figure 5). The power amplitudes in the fundamental mode values were calculated with the help of n u m p y p o l y f i t package by fitting a parabola to the power to quantify the power in the fundamental mode. The top of the fitted parabola was used further. We can see that the data points decay exponentially (Figure 6). We then fit these data points with the exponential damping profile (exp( − z/L _d )) with damping length (L _d ). These damping lengths were used to estimate the damping time with the following expression: τ d = L d / V k , where we have assumed a simple propagating wave.

In order to understand the evolution of the wavefront (how it decays and evolves to higher wavenumbers), we need to follow the wavefront as it propagates. If we consider the evolution in terms of time dependence, we can see that there will be a series of waves passing. However, the next wavefront has a similar evolution stage as the first one, except for a slight change in the density structure. This is because it has been smoothed by the mixing of density by the first wavefront. Because of the modification of the wave, we find that turbulence is affected only in a minor way, and there will be a third-order effect on the wave. Therefore, we can conclude that time variation will be third order, and hence the effect of time will be insignificant. Thus, following the wave propagation along the z axis is the crucial point of studying the energy loss.

Figure 7 shows the comparison of our theoretical model with the damping lengths from the simulation for each density contrast and velocity amplitude. The green line represents the theoretical model for the damping time scale Eq. 27. Red squares are the calculated damping time from the simulation for the first wavefront. In comparison, blue dots show the average damping time for the first five wavefronts. The error bars were calculated starting from the Fourier transform by taking the standard deviation σ, where we computed it using n u m p y s t d routine. We apply uncertainties on the polyfit function, and the estimated damping length was taken from the fitted exponential function, including the uncertainties from the covariance matrix. The errors on the blue dots are bigger than the red squares due to the loss of the parabolic shape for subsequent wavefronts. Finally, we use the following equation to find the mean and the error in the damping length when averaging over multiple wavefronts as in the paper by Van Doorsselaere et al. (2007). L ̄ d = ∑ i = 0 N − 1 L d i / σ i 2 ∑ i = 0 N − 1 1 / σ i 2 , (29) where mean error calculated with σ L ̄ d 2 = 1 ∑ i = 0 N − 1 1 / σ i 2 . (30) We run the simulation with a different driver v x ( y , t ) = V ⁡ cos ( ω t ) sin ( k y y ) exp ( κ x ) if x ≤ 0 , exp ( − κ x ) if x > 0 . where we take κ = k _y = 2 π/0.1 and V = 0.12, corresponding to the surface Alfvén wave’s analytical eigenfunction. One of the main reasons was to check if the different driver has an effect on the results. Figure 8 demonstrates the results for the simulation by the various drivers adopting the appropriate theoretical solution of the surface Alfvén wave. The z − = 0.3 and t = 3 are applied to the snapshot. We can see the non-linear deformation of the interface different than the initial simulation, but the calculated damping time matches with the theoretical solution.

We have also simulated with a smaller wavelength k _y , where we take the new wavenumber as k y ∗ = 2 k y and compare with Eq. 27. The predicted damping time for the half-wavelength in Eq. 27 gives τ _d = 0.38. The numerically estimated value for the simulation is τ _d = 0.55 ± 0.16. We believe that the discrepancy is that the nonlinear deformation is much faster and might require higher resolution. Moreover, in k y ∗ case, the first wavefront deform faster, loses its parabolic shape faster, and will have it in a shorter period than the k _y .

Lastly, we compare the wavelength of the Alfvén waves in the non-uniform regions and Surface Alfvén waves at the interface. It is known that the phase speed of surface Alfvén waves lies between V A r 2 and V A l 2 (e.g., Roberts, 1981; Goossens et al., 2012). The purpose here is to highlight the distinction between classic Alfvén waves and surface Alfvén waves. Obviously, V A r 2 v k 2 < 1 < V A l 2 v k 2 , λ A W r < λ S A W < λ A W l . It can easily be seen from the following expression that the Alfvén waves have a shorter wavelength in the upper half than surface Alfvén waves of the non-uniform region. Similarly, they have a longer wavelength in the lower half, which is also confirmed in the numerical simulation (Figure 9). This shows that the uniform driver excites both the surface Alfvén wave and classic Alfvén waves in the top and bottom of the domain. Moreover, the classic Alfvén waves show less damping because they do not have the uniturbulent cascade.