Generalized phase mixing: Turbulence-like behaviour from unidirectionally propagating MHD waves

We present the results of three-dimensional (3D) ideal magnetohydrodynamics (MHD) simulations on the dynamics of a perpendicularly inhomogeneous plasma disturbed by propagating Alfvénic waves. Simpler versions of this scenario have been extensively studied as the phenomenon of phase mixing. We show that, by generalizing the textbook version of phase mixing, interesting phenomena are obtained, such as turbulence-like behavior and complex current-sheet structure, a novelty in longitudinally homogeneous plasma excited by unidirectionally propagating waves. This study is in the setting of a coronal hole. However, it constitutes an important finding for turbulence-related phenomena in astrophysics in general, relaxing the conditions that have to be fulfilled in order to generate turbulent behavior.

phenomenology for isothermal compressible MHD turbulence 20 , but needs generalization. This implies that the generation of turbulence is no longer restricted to counterpropagating Alfvén waves in compressible MHD. In the compressible case, by using compressible Elsässer variables (defined as above but with varying density) we can still keep the Elsässer formalism, leading to different equations 16  is the speed of sound. For the sake of simplicity in this introductory section, the equations above were derived using a polytropic equation of state  ρ ρ = γ p p ( ) 0 0 , where γ is the polytropic index. Note that by letting density variations vanish in Eqs. 2, we readily recover the incompressible case of Eqs. 1. Also note that, by considering a weakly compressible case, nonlinearity arises essentially from the same nonlinear advective terms as in the incompressible case, involving both ± z . The crucial difference is, however, that generally we can no longer interpret the Elsässer variables as representing strictly parallel and anti-parallel propagating Alfvén waves, due to contributions from the compressible modes, or traduced in the fully coupled (inhomogeneous) case, due to the mixed compressible/Alfvén properties of the waves. In the incompressible case, starting from a spectrum of purely outgoing Alfvén waves (e.g. as in a coronal hole), a source of incoming waves is necessary for nonlinear interactions to sustain MHD turbulence. This source can be the inhomogeneities, so-called large-scale gradients, parallel to the mean magnetic field in the plasma (e.g. gravitational stratification in a coronal hole), causing reflections of the waves 2,21-23 . The parametric decay instability, a nonlinear instability of Alfvén waves, is also capable of generating backward propagating Alfvén waves, being more effective in the high-amplitude and/or high-frequency regime 24,25 . However, if we consider compressible modes, as stated previously, nonlinear interactions ≠ ± z ( 0) are no longer restricted to counterpropagating Alfvén waves: in an inhomogeneous medium, propagating MHD waves have mixed properties, and appear differently (predominantly Alfvén or fast characteristics) in different regions of the plasma, dictated by the local inhomogeneity 18,26 . This can be viewed, in the Elsässer formalism, as a unidirectionally propagating wave presenting both + z and − z fields to a varying degree, depending on the local inhomogeneity. That is, a wave with mixed properties is necessarily described using both ± z . The + z and − z fields of waves with mixed properties no longer represent Alfvén waves propagating in opposite directions: they propagate in the same direction and with the same phase speed: they describe a unidirectionally propagating wave. In this sense, nonlinear terms represent self-deformation of waves. In this paper we aim to verify these claims and study a perpendicularly inhomogeneous medium perturbed by unidirectionally propagating Alfvénic waves, by employing numerical simulations. Here we would like to point out the meaning of ' Alfvénic' , as formulated by: 26,27 it describes waves which have largely Alfvén characteristics, however, due to plasma inhomogeneity they are not pure Alfvén waves, as compression is also present. Alfvénic waves are an example of MHD waves with mixed properties.
In the following, we describe the numerical setup used in our simulations under the Methods section, followed by the Results. The main conclusions are drawn in the Conclusions section.

Methods
We employ 3D ideal MHD numerical simulations to test the previously described scenario of unidirectionally propagating Alfvénic waves in a perpendicularly inhomogeneous medium, using MPI-AMRVAC 28,29 , which solves the fully nonlinear, ideal MHD equations in 3D Cartesian geometry: is the total energy density. We set the adiabatic index, γ to 5

3
. We supplement the MHD equations with the ideal gas law, where the average mass per particle (in units of hydrogen atom mass m H ) is µ = . 0 6 for coronal abundances. Here, the magnetic field is measured in units for which the magnetic permeability is 1. Note that, unlike in the introductory section, a full energy equation is used, and the equations are solved for the usual variables ρ p v B ( , , , ), from which the Elsässer variables are calculated for analysis. The finite volume method uses the TVD second-order accurate solver and Woodward slope limiter. The solenoidal constraint on the magnetic field ∇ ⋅ = B ( 0 ) is enforced using Powell's scheme. The numerical domain (see Fig. 1), aimed to represent a thin elongated section of a coronal hole is × × 1 1 20 Mm in size, discretized uniformly with × × 512 512 128 numerical cells. This translates in cell sizes in the x-y plane of ≈2 km, putting heavy constraints on the CFL timestep, and limiting the achievable resolution to its present value. We consider a straight, homogeneous magnetic field of = B 5G 0 directed along the z-axis. Gravity is neglected, thus there is no stratification along the magnetic field. In this way, we eliminate possible reflections of unidirectionally propagating Alfvénic waves on equilibrium gradients along the magnetic field. Besides, having in mind the typical density scale height in a 1 MK corona (50 Mm), it is a good approximation for the present study. The plasma β = ≈ . random (in amplitude, position, and width) Gaussian density enhancements in the x-y plane, uniform along the z-axis: are amplitude, position, and Gaussian width, respectively, randomly chosen (uniform distribution) within their respective limits. These density enhancements aim to represent random perpendicular density variations in a coronal hole, in equilibrium. These kind of inhomogeneities, in reality, might arise for example from small localized heating events, leading to chromospheric evaporation. The structure is periodic in the x-y plane, as we set periodic boundaries laterally. Periodicity is ensured by letting the contribution of each Gaussian density enhancement 'pass' through the closest x-axis and y-axis periodic boundaries. The resulting mean density is ≈ . ⋅ − − 1 2 10 kgm 12 3 , while the peak value is ≈ ⋅ − − 3 10 kgm 12 3 . At the top boundary, we use a Neumann-type zero-gradient 'open' boundary condition for all variables. This is essential, as we want to exclude or minimize as much as possible the generation of reflected, counterpropagating waves. Tests with homogeneous density runs show maximum . 0 5% reflection of the incident wave energy on the top boundary. At the bottom, we employ a wave driver with properties intending to mimic observed Alfvénic waves in coronal holes 23,30 . In this sense, we use a superposition of 10 sinusoidal waves: each with a definite angular frequency ω i , velocity amplitudes U i and V i , obtained randomly from the observed log-normal distributions 23 , and random direction in the − x y plane, resulting from the different random velocity amplitudes for x and y-axis components, leading to transverse, Alfvénic waves propagating upwards on the z -axis. Note that the driver is time but not space-dependent, i.e. the whole boundary is driven equally. The resulting rms (root mean-square) velocity amplitude is ≈

Results
We run the simulation for 1000 s. During this time, perturbations originating from the boundary driver propagate upwards (on the z-axis) and interact with the inhomogeneities. The initial (equilibrium) transverse structure is quickly destroyed by the propagating Alfvénic waves, and transformed into a cross-section presenting structures on a large range of scales, reminiscent of turbulence (see Fig. 2). This observation was at the basis of a previous study which doubts the existence of packed thin individual and dynamically independent magnetic elements in the solar corona, so-called 'coronal strands' or filaments, when perturbed by propagating Alfvénic waves 31 . The current density shows the same turbulent structure in the cross-section with numerous current sheets forming (see Fig. 3) and getting dissipated by numerical resistivity. Note how the current in individual strong current sheets is more than an order of magnitude stronger than the average current density. The peak values in the current density depend on the resolution and numerical resistivity (intrinsic to the numerical scheme, not set), but should ultimately saturate at sufficiently high resolutions. We estimate the Lundquist number in the simulation to be of the order − 10 10 5 6 . In test runs with half resolution, the peak values in the current density drop approximately to half as compared to the present setup, indicating that we did not reach the saturated state. The dissipation (heating) leads to an increase in the internal energy of the plasma. However, this represents only a ≈ .
2 2% rise with respect to the equilibrium value of the internal energy at the end of the simulation time, being unable to change the plasma temperature significantly. The available energy flux from the driver is, on average, 1 that would be required to balance radiative losses 32 . As we do not focus on the energetics (processes mainly parallel to the main magnetic field) but rather on the turbulent behavior (perpendicular to the main magnetic field) in this study, heat sources or sinks, thermal conduction or thin radiative losses are not included in the energy equation. The time and space-averaged transverse magnetic field fluctuation magnitude is ≈ .
3 0% of B 0 . We calculated the 1D power spectra of magnetic and kinetic energy, density, pressure (shown in Fig. 4), and ± z (not shown, for clarity) in a plane perpendicular to the mean magnetic field B 0 . Note that the small scale generation or cascade is purely perpendicular to the mean magnetic field B 0 . We use the fft routine of IDL to obtain , , . The magnetic and kinetic energy spectrum is then x y x y x y while the calculation is similar for the other variables. Note, however, that E KM is not the total energy: in compressible MHD the total energy contains also the internal energy 16 . The time-averaged E KM spectrum ( ⊥ − . k 2 54 ) is steeper than all well-known MHD inertial range spectra (e.g. ⊥ − k 5 3  or ⊥ − k 2 in strong or weak incompressible MHD turbulence, respectively). Lower resolution test simulation runs show similar spectral slopes, however, the inertial  range is shorter, due to enhanced numerical dissipation. However, the present inertial range length (less than a decade in wavenumber) and the high variability of its slope in time suggests that the slope values should be treated carefully. Still, the obtained inertial range slopes suggest new dynamics. Indeed, our dynamics resembles an unsteady, imbalanced weak MHD turbulence, with tempo-spatial varying energy ratio of the ± z fields. We use the term 'weak' based on the observation that the cascade is essentially perpendicular to the background B 0 field. Currently, there are no analytical or even phenomenological models describing this scenario in compressible MHD. No attempt is made here to develop such a model. The closest incompressible MHD turbulence description to the present scenario is that of imbalanced weak MHD turbulence 33,34 . Note, however, that there are still fundamental differences with respect to the latter: The nonlinear cascade in our system is not driven by counterpropagating Alfvén waves, but by the mixed properties of unidirectionally propagating MHD waves in an inhomogeneous plasma, presenting self-deformation. The inertial range spectra changes in time due to the varying boundary driver, and the varying presence of + z (see Fig. 5), thus the aforementioned varying  ± ⊥ E k z ( ) 2 ) Elsässer energy ratio. In Fig. 4 we also show the energy spectrum at = t 500s, which presents an inertial range slope close to that in weak MHD turbulence 9,35 , − k 2 . In imbalanced weak MHD turbulence 36 where α ≠ 0 and corresponds to different inertial range slopes for Elsässer energies. We get ≈ − . 1 85 and − .
2 01 for the + E and − E inertial range slopes, respectively, at = t s 500 . The time-average slopes are ≈ − .
2 8, respectively. The pinning effect (converging ± E spectra at the dissipation range, see 33 ) of the Elsässer energies is observed, around = .
⊥ − k 0 02 km 1 . There is a striking difference in the inertial range slopes of pressure and density: for a polytropic law, dimensional analysis leads us to similar  − k 7 3 slopes for both variables, which should be valid even in weakly compressible cases, as the present one 9 . The deviation from this law comes from the use of a full energy equation instead of a polytopic relation, which allows us to set up a moderately inhomogeneous density profile while setting the pressure constant initially. In Fig. 5

Conclusion
We performed a first, full 3D MHD simulation of a generalized phase mixing scenario, with the varying Alfvén speed in the domain given by a random density profile perpendicular to the uniform and straight equilibrium magnetic field. The equilibrium is perturbed by upward-travelling Alfvénic waves (along the magnetic field), which originate from a driver at the lower z-axis boundary, modeled after the observed properties of waves in coronal holes. The waves leave the domain at the top z-axis boundary, and as we consider an initially uniform plasma along the magnetic field, no significant wave reflections occur. Hence we study unidirectonally propagating Alfvénic waves in a perpendicularly inhomogeneous plasma. The resulting dynamics are complex and turbulent in the cross-section, however the solution remains smooth parallel to the z-axis. Current sheets form in the process, displaying ribbon-like appearance along the main magnetic field, when viewed in 3D. These current sheets are dissipated due to numerical resistivity, and this leads to a small rise in the internal energy. This heating (at least in this setup) does not cause significant changes to the average temperature in the domain. The driven, propagating Alfvénic wave has mixed properties due to the inhomogeneities, which manifests as simultaneously nonzero ± z fields, presenting nonlinear self-deformation, in this sense. This nature of the MHD waves with mixed properties is at the origin of the complex, turbulent behaviour observed in the cross section, allowing for nonlinear cascade of wave energy to smaller scales of unidirectionally propagating waves. Let us emphasize the importance of the last sentence, as it is the main conclusion: nonlinear interactions, at the core of turbulent behaviour, are no longer restricted to counterpropagating waves in the scenario described in this paper. This study constitutes a first numerical demonstration of the recently realized 20 fact that the Alfvén effect is no longer valid in its current form in compressible MHD, but needs generalization. Thus, we simulated a potentially new type of MHD turbulence. The perpendicular spectrum of energy displays a power-law like appearance. This spectrum is highly variable in time, and its average slope for the inertial range is not in accordance with any presently available theory, being steeper than the well-known power laws of MHD turbulence. However, there is a clear separation of inertial and dissipation ranges, and the pinning effect, observed in numerical simulations of imbalanced weak MHD turbulence, is present. These findings might relax the criteria for the existence of turbulence-like cascade in plasma with a mean magnetic field and large-scale inhomogeneities perpendicularly to it. As this scenario is likely frequently present in astrophysical context, it might make turbulent behavior nearly ubiquitous.