{'fraction_non_alphanumeric': 0.04068803752931978, 'fraction_numerical': 0.02004691164972635, 'mean_word_length': 4.227399051822789, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Piecewise Parabolic Method on a Local Stencil for Magnetized Supersonic Turbulence Simulation

Stable, accurate, divergence-free simulation of magnetized supersonic turbulence is a severe test of numerical MHD schemes and has been surprisingly difficult to achieve due to the range of flow conditions present. Here we present a new, higher order-accurate, low dissipation numerical method which requires no additional dissipation or local"fixes"for stable execution. We describe PPML, a local stencil variant of the popular PPM algorithm for solving the equations of compressible ideal magnetohydrodynamics. The principal difference between PPML and PPM is that cell interface states are evolved rather that reconstructed at every timestep, resulting in a compact stencil. Interface states are evolved using Riemann invariants containing all transverse derivative information. The conservation laws are updated in an unsplit fashion, making the scheme fully multidimensional. Divergence-free evolution of the magnetic field is maintained using the higher order-accurate constrained transport technique of Gardiner and Stone. The accuracy and stability of the scheme is documented against a bank of standard test problems drawn from the literature. The method is applied to numerical simulation of supersonic MHD turbulence, which is important for many problems in astrophysics, including star formation in dark molecular clouds. PPML accurately reproduces in three-dimensions a transition to turbulence in highly compressible isothermal gas in a molecular cloud model. The low dissipation and wide spectral bandwidth of this method make it an ideal candidate for direct turbulence simulations.


Introduction
Piecewise Parabolic Method on a Local Stencil (PPML) [1,2] is a new numerical scheme developed for solving multidimensional compressible Euler (HD) and ideal mag-netohydrodynamic (MHD) equations. The method is based on a piecewise parabolic approximation of variables inside individual grid cells. It is third order-accurate in space and it implies second order temporal accuracy. This method is an improvement over the popular PPM introduced by Colella and Woodward [3,4] for nonmagnetized flows with strong shocks and extended by Dai and Woodward to compressible ideal magnetohydrodynamics [5]. PPM has been widely used in computational practice ever since, and versions of the PPM gas dynamics scheme have been incorporated into a number of codes for astrophysical applications [6].
The main new feature of PPML is the procedure for calculating interface values between adjacent cells. Instead of an interpolation procedure used in the original PPM formulation, which employs a four-point stencil, PPML relies on the information preserved from a previous time step. The required values are obtained via Riemann invariants transferred along the characteristic curves of the equations to cell boundaries using an approximate parabolic solution within a cell [7,8,9,10,11]. To preserve the order of the scheme at local extrema, a monotonicity constraint is applied to these interface values [12,13,14,15]. In a multidimensional case a monotonicity preserving method from [16] is additionally applied. The scheme is multidimensional as it keeps terms containing derivatives with respect to the tangential directions in the equations for wave amplitudes. This approach provides the left and right states for the Riemann problem based on multidimensional reconstruction. For the ideal MHD case in three dimensions, a zero-divergence constraint on the magnetic field is enforced by the use of the Stokes theorem (the so-called constrained transport approach [17]). An updated component of the electric field at a cell boundary is calculated by averaging the quantities obtained from given components of flux-vectors, taking into account a value of the electric field gradient and the information about the sign of the tangential velocity at that boundary [18].
PPML has been tested on a number HD and MHD problems that demonstrated the ability of the method to resolve discontinuous solutions without adding excessive dissipation. PPML provides a very accurate treatment of strong discontinuities, minimizes numerical dissipation of the kinetic and magnetic energy, and substantially improves the spectral bandwidth for compressible turbulence models compared to its predecessors [19].
In this paper we present a comprehensive description of the numerical method for 3D MHD simulations as well as results of numerical tests. The method has been substantially improved and remastered compared to its previous version presented in [1,2]. The new PPML features in this paper include: (i) an improved approach to maintaining zero divergence for the magnetic field following [18], see Section 5; (ii) a monotonicity constraint proposed by Rider et al. [14], as described in Section 6; (iii) an extended suite of test problems that includes a comparison with the FLASH3 MHD solver, see Section 8; Finally, in Section 9 we illustrate the performance of PPML on a threedimensional problem of highly compressible weakly magnetized forced turbulence with a sonic Mach number of 10. The problem of supersonic, super-Alfvénic turbulence with an isothermal equation of state proved to be a challenging regime for numerical modeling due to the presence of strong rarefactions and very high density contrasts in the flow. PPML scheme is perfectly stable numerically on problems of this sort. We also briefly discuss the effects of the weak magnetic field on the spectral properties of MHD turbulence and obtain good correspondence between our numerical results and observed characteristics of supersonic turbulence in star forming molecular clouds.

PPML description
Let us consider a homogeneous one-dimensional grid with the spacing h and a function q(x, t 0 ) ≡ q 0 (x) defined on this grid at t = t 0 . It is assumed that the function q 0 (x) can be approximated by a parabola inside every grid cell ( Fig. 1): where Function q i satisfies a condition Let us consider a linear advection equation where F (x, t) = a q(x, t) and find its solution for the moment t = t 0 + τ . On a discrete grid we have a number of local Riemann problems which lead to some average states q * (x i+1/2 , t) on every interface x i+1/2 between the adjacent cells. The equation (2) has only one characteristic defined by dx/dt = a. To find a value q i+1/2 , for example, on the right boundary of a cell number i at the time step t = t 0 + τ , we suggest to move along the characteristic line from the point x i+1/2 back to the moment t = t 0 and use a value from some point on a previous parabola (Fig. 2). Then for a > 0 we have where All the values in the rhs of (3) must be taken from a previous time step t = t 0 . For a < 0 a value q i+1/2 is defined by a parabola in the cell number i + 1: In non-monotonic regions we must redefine q L i and q R i . If q i is a local maximum or minimum, the interpolation function (1) must be constant, (1) can have an extremum inside the grid cell (it happens when |∆q i | < |q (6) i |). In this case we must move this extremum to the boundary of the cell. These conditions can be written as and If we know the function q(x), we can compute its average for the interval x i+1/2 − y . . . x i+1/2 (for y > 0): For a > 0, the solution of (2) at time t = t 0 + τ can be found by averaging over the interval The flux on the interface can be computed as where a + = max(a, 0) = (a + |a|) /2, a − = min(a, 0) = (a − |a|) /2. We can use an arbitrary value for q L i+1/2 if a < 0, and for q R i+1/2 if a > 0.

The governing equations
Let us consider the ideal MHD equations in 3D in the following form: Here U is a vector of eight conservative variables, F, G and H are the fluxes: where ρ is the density; u, v and w are the velocity components; B x , B y and B z are the magnetic field components; E is the total energy and p is the total pressure: We have included the factor 1/ √ 4π in the definition of B. An equation for the total energy and an ideal gas equation of state are where γ -the adiabatic index, ε -the specific internal energy. If we denote we can write the sound velocity c, Alfvén velocity c a , the fast and the slow magnetoacoustic velocities c f, s as We will also deal with a non-conservative form of the MHD equations: where The matrices A, B and C can be computed using Jacobians of (10) and a transition matrix For example, the Jacoby matrix A is The corresponding eigenvalues are x represent a pair of fast magneto-acoustic waves, λ 2, 7 x -a pair of Alfvén waves, λ 3, 6 xa pair of slow magneto-acoustic waves, λ 4 x -an entropy wave, λ 5 x -a magnetic-flux wave. The eigenvectors of the Jacobians could have singularities at the points of degeneracy of the eigenvalues since the MHD equations are nonstrictly hyperbolic. To avoid those, Brio and Wu [20] suggested a scaled version of the eigenvectors that comes from defining Thus the left and the right eigenvectors are

A numerical scheme
To solve (10) we apply a conservative difference scheme: Half-integer indices such as i + 1/2 denote the boundaries of grid cells, half-integer time index n + 1/2 means that we use the averaged values of the fluxes over a time step τ in order to get a second-order temporal accuracy. The solution inside every grid cell is approximated by a parabola along any Cartesian grid axis. The boundary values for each parabola are determined from a conservation property of Riemann invariants that remain constant along the characteristics of the initial linearized system of equations. Parabolae must be built using the primitive variables (16), so we need to consider a non-conservative form of the MHD equations (15).
We can expand A, B and C in (15) into their eigenvectors. For example, for the x-direction: where R x is a matrix with columns filled by the right eigenvectors r p x (p = 1, . . . , 8), is an inverse matrix, with rows filled by the left eigenvectors l p x . Λ x is a diagonal matrix of the eigenvalues: To construct piecewise parabolae for every time step, one needs to define the states on the cell edges and the states at their centers. For simplicity, let us further consider a 1D case: Inserting (19) into (20) and multiplying by the L matrix from the left, we arrive at Let us expand a vector V(x, t) into the local basis of the right eigenvectors r p , which are fixed in every cell: Inserting (22) into (21), we arrive at The equations (23) mean that the coefficients α p (x, t) in the expansion (22) (the wave amplitudes) must be constant along the characteristics x p (t): i.e. α p (x, t) are Riemann invariants. A value of Riemann invariant on the boundary of a cell (x = x i+1/2 ) at the moment t + τ could be computed using its value at the moment t as Fig. 3 shows two adjacent cells i and i + 1. The characteristics in the cell i have index p 1 , in the cell i + 1 -index p 2 . One of the characteristics x p 1 (t) with the eigenvalue λ p 1 > 0 is shown in the cell i, another one x p 2 (t) with the eigenvalue λ p 2 < 0 is shown in the cell i + 1. According to (24) the amplitude of a wave at point 3, which propagates inside the cell i along the characteristic x p 1 (t) with the eigenvalue λ p 1 , is equal to its value at point 1. In the same way the amplitude of a wave at point 3, which propagates inside the cell i + 1 along the characteristic x p 2 (t), is equal to its value at point 2.
The state at point 3, which is computed according to (22) by summation with respect to all the eigenvectors, fixed in cell i, with λ p 1 > 0 will be on the left side of the interface. Let V L denote this value. Similarly, let V R denote the state at point 3 on the right side of the interface which is computed by summation with respect to all the eigenvectors fixed in cell i + 1, with λ p 2 < 0.
The amplitudes of waves α p (x p , t) at the point (24), which influence the right boundary of cell i (λ p > 0), can be computed by multiplying the expansion (22) by the left eigenvectors, fixed in cell i: where We can use arbitrary values for the wave amplitudes for λ p < 0 because these waves have no effect on the right boundary of the cell and will be omitted in the sum (22). For convenience let them be With the help of Θ-function we can rearrange (25)- (27) as Multiplying (28) by r p and summing over all p such as λ p > 0 according to (22), after some simple manipulations we obtain the boundary value at time t + τ : If we consider cell i + 1 and waves with λ p < 0 we will obtain a similar expression for the value V R on the right side of the interface at time t + τ : 9 The left and the right eigenvectors in (29)- (30) are fixed in every cell. To compute them we can use a state from any point x inside the cell -it has been shown by numerical experiments that this choice has no influence on the solution. We suggest using the values of states in the centres of the cells, i.e. l p = l p (V(x i , t)).
In order to solve the 3D problem we split the initial set (15) by the space variables and solve the 1D equations separately for the x-, y-and z-directions. However, in this case we have an additional change in the quantities because of the fluxes in the orthogonal directions. For example, a flux in the y-direction will affect the quantities at point 1 between cells (i, j) and (i + 1, j) (see Fig. 5). To obtain the correct result, we can solve equation (15) in the x-direction considering the terms B ∂ y V and C ∂ z V as sources.
Then instead of (23) we arrive at where D p are the components of the vector The components of the partial derivatives in (32) can be calculated as The solution of (31) can be obtained from a Taylor series expansion of α p (x, t) near the boundary x = x i+1/2 . It is similar to (24) but has the additional term Then instead of (29)- (30) we arrive at To obtain states V i+1/2 in PPML, we solve the Riemann problem for the V L and V R states on every interface using, e.g., the Roe solver [21] or the HLLD solver [22]: where R is the Riemann solver. Note, that in the original PPM, the states V i+1/2 are obtained through monotonic interpolation [7]. We then apply a monotonicity-and extrema-preserving procedure proposed by Rider et al. [14] to the values of V i+1/2 , as described in Section 6. Finally, we modify the resulting interface states with the standard PPM monotonicity procedure (5)- (6). So far we have obtained the boundary values of piecewise parabolae at time t + τ in each grid cell. Now we need to define the fluxes to compute the new central states. Again we can use the Roe solver [21] or the HLLD solver [22] with V L and V R from (35)- (36) but in this case the numerical scheme will have the first order of temporal accuracy. To design a second order scheme, we must average the amplitudes α p (x, t) over the zones of influence. Fig. 6 shows a set of characteristics corresponding to waves with λ p > 0. The characteristic x 1 (t) has the maximum eigenvalue λ 1 . Point 1 is the point of intersection between this characteristic and the piecewise parabola at time t. Obviously, only the zone between the interface x = x i+1/2 and point 1 affects the left boundary state at point 2. 11 If we consider a wave which propagates along the characteristic with λ p > 0 inside cell i, its averaged amplitude on the interface x = x i+1/2 at time t + τ can be calculated as Multiplying the expansion (22) by the left eigenvectors, fixed in cell i, will yield Inserting (39) into (38) and removing the factor l p from the integral, we arrive at After that we can arrive at which is similar to (35). Here V L, 1 i+1/2 is the averaged by formula (40) solution V(x) at time t over the zone of influence of the wave in cell i, corresponding to the maximum eigenvalue (λ 1 > 0).
For cell i + 1 and the negative eigenvalues we arrive at  (7) and (8), respectively.
The time step τ in (18) is obtained from the Courant condition where σ is the Courant number and c x f , c y f , and c z f are the fast magneto-acoustic velocities along the coordinate directions.

Zero divergence constraint for the magnetic field
Our numerical method must provide numerical solutions that satisfy a condition div B = 0.
There are several approaches to this problem in the literature [18,23,24,25,17,26]. In our numerical scheme we used an unsplit Godunov method for ideal MHD with a constrained transport developed in [18] that is based on the Stokes theorem Calculated from (46) the magnetic field obviously satisfies condition (45). To apply the Stokes theorem in 3D case, we define the magnetic field components on the cell faces and the electric field components on the cell edges. The algorithm exploits the fact that some components of the fluxes F, G and H are actually the components of the electric field.
For simplicity let us consider a 2D case. Thus the sixth component of the flux F is the z-component of the electric field reversed in sign, which correspond to the left and right boundaries of a cell. The fifth component of the flux G is the z-component of electric field, which correspond to the top and bottom boundaries of a cell. Computing Taylor series for these components near the nodes of the mesh and averaging them, we obtain the components of the electric field at the nodes (Fig. 7).
where for example These values of E z i+1/2,j+1/2 are used in a discrete version of (46): Updated components of magnetic field at the center of (i, j)-cell are computed by averaging: The magnetic field computed this way will automatically satisfy (45). To make sure that this is indeed the case, divB should be approximated by the following expression

Monotonicity constraints
A standard PPM monotonicity preserving procedure (5)-(6) is insufficient for the ideal MHD case. We need to use additional procedures to suppress spurious oscillations. Procedure 1. To keep a solution monotonic without reducing the order of the scheme and to preserve all the local extrema in MHD simulations we can follow a number of approaches described in the literature [12,13,14,15]. In our method we rely on a Piecewise Parabolic Accurate Monotonicity-and Extrema-Preserving procedure described in [14].
We compute the wave amplitudes α p i for the central states V i and α p i±1/2 for the interface states V i±1/2 . We then calculate new values for wave amplitudes of the interface states: If α p, * * i±1/2 = α p i±1/2 for all p, the procedure is completed. Otherwise, we compute a set of fifth-order WENO interface values α p i±1/2, * , using Algorithm 2.2.4 from [14], and check for a local extremum. If α p, * * i±1/2 = α p i for all p, we apply Algorithm 2.1.2, 4(b), otherwise the region is monotonic but too steep to be approximated with the values obtained from the Riemann solver (37) and we apply Algorithm 2.1.2, 4(c). Procedure 2. To keep a solution monotonic in a multidimensional case, we employ a method proposed in [16].
Let us consider a 2D case for simplicity. A parabola that approximates a solution along the x-axis for every component V (x) of a state V(x, t) at some point in time can be defined as As a limiting function φ(V ) we can use that described in [16]: where l = i − 2, i − 1, i, i + 1, i + 2, m = j − 2, j − 1, j, j + 1, j + 2 except (l, m) = (i, j). In a 3D case, the limiting function (47) must include all the neighbors of the cell (i, j). 15

A FORTRAN implementation
The algorithm for computing the left and the right boundary values V L and V R (35)-(36) could be implemented in FORTRAN this way: integer n, n2 real dt, dx real VL (8) To obtain a monotonic solution for a 3D MHD problem, we suggest using the following algorithm: 6. Solve Riemann problem between V L and V R to determine the new interface states. 7. Apply procedure 1 in x-direction (section 6). 8. Apply procedure 1 in y-direction (section 6). 9. Apply procedure 1 in z-direction (section 6). 10. Apply procedure 2 (section 6). 11. Apply PPM procedure (5)- (6).
Note that only the interface states are modified by the monotonicity preserving procedures. The solution of the Riemann problem joins the interface values, but the monotonicity procedures split them again into the left and right values.

Riemann problem with multiple weak discontinuities
This is a 1D problem from [27]. An interval x ∈ [0 . . . 1] is divided in two by x = 0.5. The left and the right states at the initial moment are defined as ρ L , u L , v L , w L , B L y , B L z , p L = (1.08, 1.2, 0.01, 0.5, 3.6, 2, 0.95),

Numerical dissipation and decay of Alfvén waves
Numerical calculations on a discrete grid always lead to energy loss due to numerical dissipation. In order to estimate the properties of numerical dissipation of the PPML ideal MHD scheme, we used the test problem from [28] and followed a decay of twodimensional Alfvén wave. We used a standing wave propagating along the grid diagonal with initial conditions δv x = v amp c a sin(k x x + k y y),  Figure 11 shows the envelope for the maxima of z-component of the magnetic field and velocity obtained with PPML and PPM reconstruction procedures as a function of time. While both schemes show very low dissipation, PPML dissipation is even smaller than that of PPM.

Travelling circularly polarized Alfvén wave problem
This problem was suggested in [23] where the solution on the mesh N = 128 regarded as the exact one U E i,j . The rate of convergence was calculated as follows   where δ N is an averaged value: The calculations were carried out up to t = 5 with a Courant number σ = 0.4 and γ = 5/3. We applied periodic boundary conditions. Figure 12 demonstrates the convergence of the numerical solution. It shows the orthogonal component B ⊥ in travelling Alfvén wave for computations on meshes of different size N . Table 1 gives average relative numerical errors and rates of convergence obtained for the PPML scheme.

Rotor problem
The Rotor problem was suggested in [26] and has been widely used to test numerical schemes [e.g., 23,29,30]. It turned out to be a hard nut to crack for many codes due to the appearance of negative pressure values [23].
The computational domain in this case is a square [0, 1]×[0, 1] with a uniform pressure p = 1 and magnetic field components B x = 5/ √ 4π, B y = 0. There is a rotating disk of dense fluid at the center with a radius r 0 = 0.1. For r < r 0 we specify ρ = 10, For r > r 1 = 0.115 the fluid is initially at rest (u = v = 0) with density ρ = 1. In the intermediate zone r 0 < r < r 1 we use a linear interpolation of the variables: ρ = 1 + 9f , r 0 ). In this setup, the initial configuration is imbalanced due to centrifugal forces. The rotating fluid will tend to equilibrate, while the magnetic field holds the oblate shape of the rotor.
The  We compare numerical solutions obtained with PPML with those presented in [23,29,30]. In Table 2 the average relative numerical errors and the rates of convergence are given for PPML and for Flash3 USM-MEC (unsplit staggered mesh algorithm with modified electric field construction introduced in [30]). Both PPML and Flash3 USM-MEC codes use a Roe solver [21] for this test. The relative numerical errors were computed using eq. (48), where for U E i,j we used the highest resolution result (N = 400). The average errorδ N (U ) is defined as the average δ N (U ) for all non-zero variables U . The rate of convergence is estimated as in (49). PPML results are more accurate and have a comparable rate of (self-)convergence with those from the new Flash3 MHD solver.

Orszag-Tang vortex problem
This problem was suggested in [31] and since then has been used in many papers as a standard test problem for numerical codes in 2D MHD. It involves formation and an interaction of multiple shocks and a transition to supersonic turbulence.
In the computational domain [0, 1] × [0, 1], we set a uniform density ρ = 25/(36π) and pressure p = 5/(12π) with γ = 5/3 (in this case the sound velocity c = γp/ρ = 1). The initial velocities and components of the magnetic field are set using harmonic functions: u = − sin 2πy, v = sin 2πx, w = 0, B x = −B 0 sin 2πy, B y = B 0 sin 4πx, B z = 0, where  We carried out computations on N ×N meshes with several values of N using periodic boundary conditions and a Courant number σ = 0.3. Figure 15 demonstrates the pressure distribution at time t = 0.5 for N = 256. In Fig. 16 the pressure distributions along the 23 lines y = 0.3125 (j = 83) and y = 0.4277 (j = 112) are shown to illustrate the accuracy and sharpness of the main flow features. Table 3 contains the average relative numerical errors and the rates of convergence for PPML and Flash3 USM-MEC [30] solvers at t = 0.5. As an "exact" solution U E i,j , we used a solution obtained on the grid with N = 400. The corresponding PPML results are more accurate and demonstrate better convergence.

A compressible turbulence simulation
Magnetized supersonic turbulence plays an important role in statistical star formation theories [35]. This stimulated development of accurate numerical methods suitable for modeling turbulent molecular clouds. One of the motivations behind the design of PPML has been a need for an MHD scheme with low numerical dissipation comparable or better than that of PPM. In this section we illustrate the performance of PPML on a challenging problem of forced super-Alfvénic turbulence. Some numerical methods that successfully pass the tests discussed above turn unstable on this application. Since adding more dissipation where needed -the usual way to cure for "blow ups" caused by numerical instabilities -would ultimately damage the derived statistics of turbulence [33], the issue of inherent stability of numerical methods is crucial for both supersonic turbulence and star formation simulations.
For illustrative purposes, we present here a simulation of weakly magnetized supersonic turbulence. In this experiment, turbulence in a periodic domain of linear dimension L = 1 is driven by a large-scale solenoidal force for 8 flow-crossing times t d ≡ L/2M s . At time t = 0, a uniform gas with density ρ ≡ 1 is permeated by a weak uniform magnetic field B 0 x x x, such that β 0 ≡ 2p/B 2 0 = 20. We apply an initial large-scale velocity field that corresponds to an rms sonic Mach number M s ∼ 10 and assume an isothermal equation of state (c ≡ 1) to mimic the average physical conditions in the dense parts of molecular clouds (n = 10 3 cm −3 , T = 10 K).
The evolution begins with the formation of strong shocks on "caustics" of the initial velocity field. Shock interactions then cascade the initial kinetic energy of large-scale motion of the gas down to smaller and smaller scalesá la Kolmogorov-Richardson, see Fig. 17. The magnetic field gets amplified by a factor of about 50 via the small-scale dynamo action [34]. The large-scale solenoidal force (acceleration) keeps the rms sonic Mach number roughly constant at M s ≈ 10. The evolution of kinetic and magnetic energies is shown in Fig. 18, left panel. Also shown is max(|∇· ∇· ∇·B|) as a function of time during this simulation. The method keeps the absolute value of the divergence 24 of magnetic field below 10 −12 at all times, even after 70,000 integration time steps (if double precision is used). After about 4 crossing times of evolution, the system completes a transition to a fully developed isotropic state with M s ≈ 10 and M A ≈ 10. The right panel of Fig. 18 illustrates spectral characteristics of this saturated state by showing the time-average (over 25 flow snapshots taken between t = 4t d and t = 8t d ) power spectra for the density, velocity, and magnetic field strength.
The velocity spectrum has an extended scaling range with a slope of about −2, as in the Burgers turbulence, similar to the corresponding scaling in non-magnetized flows at high Mach numbers [32]. This is expected, as turbulence here is only weakly magnetized. The density spectrum slope of about −0.7 is again consistent with our previous results for non-magnetized flows obtained with PPM (−1.0 at M s = 6) and with an anticipated trend towards a flat "white noise" spectrum at M s = ∞. The magnetic energy spectrum does not show a clear scaling range, as expected at this modest resolution, assuming the effective magnetic Prandtl number of PPML is of order unity. We also looked at more advanced spectral characteristics for compressible flows, such as the power spectrum of ρ 1/3 v, and found a slope of −1.7. This power spectrum related to the energy transfer rate in wavenumber space is insensitive to the turbulent Mach number and should have a Kolmogorov −5/3 slope [32] in both incompressible and highly compressible regimes, although a steeper scaling does occur due to intermittency [36]. This is true for both non-magnetized and weakly magnetized flows.
We have also carried out two additional simulations of the same kind but with higher degrees of magnetization, β 0 = 2 and 0.2 [37]. While the saturated turbulent state in the β 0 = 2 variant is still super-Alfvénic with M A ≈ 3 and the magnetic energy is about 3 times smaller than the kinetic energy, the trans-Alfvénic case, β 0 = 0.2, reaches an equipartition of kinetic and magnetic energies. In both cases, PPML proved to be perfectly stable at a Courant number σ = 0.2, as in the super-Alfvénic case β 0 = 20 discussed above.
Our approach to handle the stability issues in MHD turbulence simulations with PPML is as follows: (i) we use locally multidimensional reconstruction that improves the quality of Right and Left interface states and helps to avoid numerous well-known pathologies, such as "carbuncles", etc. [38]; (ii) we control the quality of these states before moving forward with the flux calculation; (iii) if we find that the states are not satisfactory, we reduce the order of reconstruction to linear or further skip the whole reconstruction step; (iv) we use only nonlinear Riemann solver in all cases lacking the reconstruction step.
Overall, the derived spectral properties of weakly magnetized highly compressible turbulence demonstrate that low dissipation and wide spectral bandwidth of PPML make it an ideal numerical scheme for large-scale simulations of magnetized supersonic turbulence.

Conclusions
In this paper we presented PPML, a new numerical method for compressible ideal MHD that is based on the piecewise parabolic approximation. Interface values for the interpolation parabolae in every grid cell are defined with the help of the Riemann invariants which remain constant along the characteristics. The monotonicity of the states on the interfaces between adjacent cells is provided by the monotonicity-and extremapreserving procedure from [14]. The scheme is fully multidimensional as it includes the terms corresponding to the tangential directions in the amplitude equations. This helps to avoid numerous well-known pathologies, such as "carbuncles", etc.
The states in the cell centers are defined by the conservative difference scheme (18). To obtain the second-order temporal accuracy we must average the wave amplitudes over the corresponding domains of influence. To define the fluxes we need to solve the Riemann problem between the states at the cell interfaces computed with the averaged amplitudes.
To preserve zero divergence of the magnetic field in three dimensions, we use an unsplit Godunov method based on the constrained transport approach [18]. We use the information about the magnetic field gradients to fulfill the constraint on the magnetic field more accurately.
We tested the performance of PPML on several numerical problems which demonstrated its high accuracy on both smooth and discontinuous solutions. Two-dimensional flow fields generated by PPML are highly resolved without any wiggled contour lines. Our pilot simulations of supersonic magnetized turbulence in three dimensions with PPML show that low dissipation and wide spectral bandwidth of this method make it an ideal candidate for direct turbulence simulations.