Probing the nature of dissipation in compressible MHD turbulence

In the present study, we run a set of simulations of decaying magnetohydrodynamics (MHD) turbulence.

The numerical scheme we use (Godunov) implicitly introduces dissipation to evolve the ideal MHD equations, but as stated above, we incorporate additional explicit physical dissipation terms in our evolution equations. Indeed, it is important to retain some amount of physical viscosity as Godunov schemes do not provide an implicit viscosity in shear layers. Here, we discuss our methods to estimate the fraction of the dissipation which is due to the numerical scheme.

We set values for the viscous and resistive coefficients $\nu$ and $\eta$ identical to those used by Momferratos et al. (2014) in pseudo-spectral simulations with 512³ spectral elements: $\nu=\eta=7\times 10^{-4}$ in the same non-dimensional units. This is motivated by the common belief that spectral codes are approximately twice more efficient as grid based codes. Our study for shocks in appendix B presents a more detailed picture. Figure 21 shows the dissipation bump in a fiducial shock front at various resolutions. For our chosen values for the dissipative coefficients and a resolution of $N=1024$, we see it is effectively spread up by nearly a factor of three, while one would have to increase the resolution by a factor 8 to fully resolve it. A resolution twice less would spread the front by a factor six and thus our current choice is a good compromise between accuracy and CPU efficiency.

In isothermal MHD, the integrated total isothermal generalised mechanical energy $E=<\frac{1}{2}\rho u^{2}+\frac{1}{8\pi}B^{2}+p\log\rho>$ decreases due to all irreversible processes taking place (see equation 51). Because our Godunov time integration scheme is conservative to round-off error, we can use its time derivative to estimate the global budget of dissipated energy:

where $<\varepsilon_{\mathrm{tot}}>$ is the total rate of irreversible heating integrated over the whole computational domain. Appendix B presents and tests a new method to estimate locally the total irreversible heating $\varepsilon_{\mathrm{tot}}$. The chosen method has the additional advantage that it preserves to round-off error the validity of equation (9) when integrated over the whole domain.

We can now decompose the local total heating rate as

where

and

are the local viscous and resistive dissipative heating rates, and $\varepsilon_{\rm num}$ is the dissipation due to the numerical scheme.

We can then estimate the local numerical dissipation rate simply by computing $\varepsilon_{\rm num}=\varepsilon_{\rm tot}-(\varepsilon_{\nu}+\varepsilon_{\eta})$ where we use well centered estimates for equations (11) and (12). If our estimate for the local dissipation were perfect, this quantity would always be positive, because we are performing our simulations with a time step small enough for the scheme to be stable (it is set to 70% of the shortest unstable time step). However, we are subject to truncation errors in both the $\varepsilon_{\rm tot}$ term (see appendix B) and the $\varepsilon_{\nu}+\varepsilon_{\eta}$ terms (where a centered difference is used). The difference between the two terms can hence be negative due to these truncation errors. We thus define $\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}$ as a corrected local total dissipation rate which ensures the resulting estimate for $\varepsilon_{\mathrm{num}}$ is positive. It is equal to the total local dissipation $\varepsilon_{\rm tot}$ where the numerical dissipation is positive (i.e. where $\varepsilon_{\rm tot}>(\varepsilon_{\nu}+\varepsilon_{\eta})$), while it is equal to the total physical dissipation $\varepsilon_{\nu}+\varepsilon_{\eta}$ elsewhere. This ensures the corrected local numerical dissipation rate $\varepsilon_{\mathrm{num}}^{\mathrm{corr}}=\varepsilon_{\rm tot}^{\mathrm{corr}}-(\varepsilon_{\nu}+\varepsilon_{\eta})$ is always positive. In particular, the local corrected total dissipation rate $\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}$ is always greater than $\varepsilon_{\mathrm{tot}}$. It is then shared between resistive and viscous natures in the same proportions as the physical terms we introduced to provide estimates for viscous and resistive dissipations including numerical dissipation:

Figure 1 displays the temporal evolution of various total dissipation rates. Thanks to the equality (9), we can compute the exact total dissipation rate at each time step (blue curves), and we can compare it to the integrated local estimate $<\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}>$ (green curves) which by construction is always greater. The difference between the two gives an estimate of the error we make on the estimation of the dissipation (on the order of a percent at most). It corresponds to the integrated estimated $\varepsilon_{\mathrm{num}}$ in all the pixels where it is negative. The orange curves show the integrated physical dissipative terms $<\varepsilon_{\nu}+\varepsilon_{\eta}>$. They amount to about two-thirds of the total, while the remainder is numerical dissipation by the scheme.

We know local intense dissipation events are caused by strong variations of some of the fluid state variables. Here, we want to identify regions where the fluid state varies strongly and to characterise its variations in each direction.

The fluid state is characterised by the seven (1+3+3) components of $\@vec{W}=\left(\rho,\mathbf{u},\mathbf{B}\right)$, which do not have the same physical dimensions. We want to put the variations of density, velocity and magnetic fields on equal footing. Hence, we need to rescale the gradient of each component of $\@vec{W}$ to make them homogeneous to the same physical dimension. We now choose to define the rescaled gradient of $\@vec{W}$ in a given direction $\@vec{r}$ as

where $\@vec{\hat{r}}=\@vec{r}/r$ is the unit vector in the direction of $\@vec{r}$. This rescaled gradient has the dimension of the inverse of a length scale, which represents the typical length scale over which the state variables vary in the direction $\@vec{\hat{r}}$.

The norm of this gradient will be large whenever there is a rapid change in one or several state variables. Its square can be expressed as

where $\alpha_{ij}=\partial_{i}\@vec{W}\cdot\partial_{j}\@vec{W}$ is a 3$\times$3 matrix (and the dot product applies to the seven components vectors) with coefficients homogeneous to an inverse squared length. It is real symmetric, and therefore diagonal in an orthonormal basis. We can rewrite equation (16) in a more explicit form

where $\ell_{\rm scan}^{2}$, $\ell_{\perp 1}^{2}$, and $\ell_{\perp 2}^{2}$ are the inverse of the eigenvalues associated to the eigenvectors $\@vec{\hat{r}}_{\mathrm{scan}}$, $\@vec{\hat{r}}_{\perp 1}$, and $\@vec{\hat{r}}_{\perp 2}$ of the matrix $\alpha_{ij}$. Equation (17) shows how the gradient of state variables depends on directions. A 3D polar plot of the norm of this gradient takes the form of an ellipsoid whose principal axes are in the three orthogonal eigenvalue directions of the above matrix:

with the three length scales ordered so that $\ell_{\rm scan}\leq\ell_{\perp 1}\leq\ell_{\perp 2}$. These three variation length scales and their associated orthogonal directions characterise the local geometry of the gradients of the fluid state variables.

Figure 2 shows how the aspect ratios between these typical variation length scales are distributed in all cells of a simulation (left panel) and for only highly dissipating ones (four standard deviations over the mean, right panel). It shows that most fluid state variables vary primarily in one direction for extreme dissipation events, whereas aspect ratios span all possibilities if we consider the full simulation domain. We also notice a slight imbalance towards ribbons compared to sheets. When one variation direction is dominant ($\ell_{\rm scan}<<\ell_{\perp 1}\leq\ell_{\perp 2}$), quantities are essentially constant in the direction orthogonal to it and the local situation is hence nearly 1D plane-parallel. We thus define the planarity as the ratio $\ell_{\perp 1}/\ell_{\rm scan}$ which is large whenever $\ell_{\rm scan}<<\ell_{\perp 1}$, i.e. when the local geometry is close to plane-parallel.

This one-dimensional geometry of gradients for intense dissipation regions is consistent with the typical two-dimensional geometry of the structures found in MHD turbulence (Uritsky et al., 2010; Zhdankin et al., 2013; Momferratos et al., 2014). On intense dissipation structures, we should thus be able to capture most fluid variations by browsing those along the maximum gradient direction. In section 3.2 we use $\@vec{\hat{r}_{\mathrm{scan}}}$ as a sampling direction to probe the variation of physical quantities around strong dissipation regions.

In the ideal case where the gradient would be strictly in one direction, the gas dynamics are governed by 1D plane-parallel MHD equations, and we show here how local gradients can be projected onto ideal MHD waves.

We write $x$ the space coordinate along the direction of the gradient and $t$ the time coordinate. The requirement $\@vec{\nabla}.\@vec{B}=0$ implies that $\partial_{x}B_{x}=0$: the corresponding component of $\partial_{x}\@vec{W}$ is thus zero. It turns out that the six non zero components of $\partial_{x}\@vec{W}$ are spanned by the six ideal MHD waves, as we now turn to show.

Wave solutions take the form $\@vec{W}(x,t)=\@vec{F}(x-\varv t)$ where $\varv$ is the traveling speed of the wave. We note that $\partial_{t}\@vec{F}=-\varv\partial_{x}\@vec{F}$ and plug this form into the ideal MHD part of the equations (without the dissipation terms). We arrive at a linear eigenvalue problem for which we can find six eigenvectors $\partial_{x}\@vec{F}$ with eigenvalues $\varv$ corresponding to the six waves of ideal isothermal MHD²²2Note that formally, finding the gradient $\partial_{x}\@vec{F}$ is equivalent to solving the amplitude for the linear wave problem when we identify $\partial_{t}\equiv i\omega$, $\partial_{x}\equiv ik$ and $\varv=\omega/k$ matches the phase velocity.. We label them by their wave type $s,i,f$ for slow, intermediate and fast and their direction of propagation $R,L$ for right (or forward, $\varv>0$) and left (or backward, $\varv<0$). To within a multiplicative constant, the expressions for these eigenvectors are (see section 5.2.3 of Goedbloed et al. (2019) or section 6.5 of Gurnett & Bhattacharjee (2005), for example) :

• •

  for intermediate (Alfvén) waves

  $$\partial_{\@vec{\hat{x}}}\@vec{F}^{R,L}_{i}\propto(0,\epsilon^{R,L}\@vec{a}_{t}^{\perp},-\mathrm{sign}(a_{x})\@vec{a}_{t}^{\perp})$$

  (19)

  where $\epsilon^{R,L}=-1$ for left-going (backward going) waves and $\epsilon^{R,L}=1$ for right-going (forward going) waves, $\@vec{a}=\@vec{B}/\sqrt{4\pi\rho}$ is the Alfvén velocity vector, $\@vec{a}_{t}$ is the transverse component of $\@vec{a}$, $a_{x}$ is its $x$-component and $\@vec{a}_{t}^{\perp}$ is $\@vec{a}_{t}$ rotated by $\pi/2$ in the transverse plane. The first component of this gradient is zero, hence the density is uniform. And the transverse magnetic field has its gradient orthogonal to itself: it rotates along the scanning direction. The corresponding travelling speed is $c_{i}^{R,L}=\epsilon^{R,L}|a_{x}|$.

• •

  for fast and slow magnetosonic waves

  $$\partial_{\@vec{\hat{x}}}\@vec{F}^{R,L}_{s,f}\propto(-\frac{c}{c^{R,L}_{s,f}},\@vec{\hat{x}}-\frac{a_{x}}{d}\@vec{a}_{t},\frac{c^{R,L}_{s,f}}{d}\@vec{a}_{t})$$

  (20)

  where the propagation speed $c^{R,L}_{s,f}$ reads:

  $$c^{R,L}_{s,f}=\epsilon^{R,L}\sqrt{(c^{2}+a^{2})+\epsilon_{f,s}\sqrt{(c^{2}+a^{2})^{2}-4a_{x}^{2}c^{2}}}$$

  (21)

  with $d=(c^{R,L}_{s,f})^{2}-a_{x}^{2}$ and $\epsilon_{f,s}=1$ for fast waves or -1 for slow waves. These waves are compressive (the density gradient is non-zero) and the gradient of transverse magnetic field is aligned with itself. In other words, the transverse magnetic field remains in the same direction, which also happens to be the same direction as the variation of the transverse velocity. Both the velocity and the magnetic field vectors thus remain in the plane defined by the scanning direction $\@vec{\hat{x}}$ and the initial transverse field (a property sometimes referred to as the coplanarity of these waves).

These six gradients form an orthogonal basis which can be easily normalised to make it an orthonormal basis $\@vec{\hat{e}}^{R,L}_{s,i,f}=\partial_{\@vec{\hat{x}}}\@vec{F}^{R,L}_{s,i,f}/||\partial_{\@vec{\hat{x}}}\@vec{F}^{R,L}_{s,i,f}||$ .

Any gradient $\partial_{\@vec{\hat{x}}}\@vec{W}$ can now easily be decomposed into the six waves by computing the scalar product $\alpha^{R,L}_{s,i,f}=\hat{\@vec{e}}^{R,L}_{s,i,f}.\partial_{\@vec{\hat{x}}}\@vec{W}$. Thanks to orthonormality, we have $\sum_{R,L,s,i,f}(\alpha^{R,L}_{s,i,f})^{2}=||\partial_{\@vec{\hat{x}}}\@vec{W}||^{2}$ and each coefficient $(\alpha^{R,L}_{s,i,f})^{2}/||\partial_{\@vec{\hat{x}}}\@vec{W}||^{2}$ can be interpreted as a 0 to 1 coefficient which characterises how similar the gradient $\partial_{\@vec{\hat{x}}}\@vec{W}$ is to the corresponding ideal MHD wave. We call ’most representative wave’ the wave with the largest coefficient in this decomposition. The most representative wave characterizes the local gradient as slow, intermediate or fast, each one in a left (backward) or right (forward) going version depending on the sign of its speed relative to the fluid $c^{R,L}_{s,i,f}$. Note also that this decomposition does not change if we add a constant vector to the velocity: it is independent from the choice of Galilean frame.

Until now, we have only considered wave solutions of the ideal part of the MHD equations (without dissipation), while the gradients in our simulation result from the evolution of fully dissipative MHD. Let us now consider a non-linear wave solution of the 1D fully dissipative MHD $\@vec{F}_{\mathrm{full}}(x-\varv_{\mathrm{full}}t)$, such as the isothermal shocks of appendix A, for example. The profile of this wave continuously joins two uniform states related by the Rankine-Hugoniot relations (see subsection 3.3). These two states are separated by a region where dissipation occurs. Consider the gas state at the local maximum of dissipation: this is where the gradients of the state variables are the largest, and where the gradient of viscous and resistive stresses are likely to be small (because we are close to their maximum). At this position, the 1D dissipative physics behaves like the 1D ideal physics, and we can expect that the measured gradients fall along one of the ideal wave gradients we described above. As a result, the fully dissipative wave speed should be well approximated by its ideal estimate: $\varv_{\mathrm{full}}\simeq u_{x}+c^{R,L}_{s,i,f}$ where $u_{x}$ is the fluid velocity and $c^{R,L}_{s,i,f}$ applies to the most representative wave at the dissipation maximum. We will make use of this fact in the following to estimate the steady-state velocity of the structures we detect (see section 3.3.2). Furthermore, we investigated the gradients of semi-analytic isothermal shock profiles (computed in appendix A) and we noticed that gradients in slow shocks are dominated by slow magnetosonic waves all along their profiles. Similarly fast shocks gradients are dominated by fast magnetosonic waves. This result seems natural but we find it nevertheless surprising that dissipative physics does not affect more the nature of gradients and we did not yet find a satisfactory explanation for this behaviour.

Finally, note that we can always decompose a gradient in a given direction, but it makes less sense if the 3D gradient is not strongly dominated by a single direction. By selecting intense dissipative cells, however, we are more likely to be in a situation where the gradient is well directed (see figure 2 and previous subsection).

In turbulent MHD flows, the bulk dissipation of kinetic and magnetic energy occurs in a small volume compared to the global scale of the flow. Dissipation has been analyzed and observed in several studies (e.g. Uritsky et al., 2010; Zhdankin et al., 2013; Momferratos et al., 2014) to be organized in coherent structures, ribbon-shaped or sheet-like.

Figure 3 shows isocontours of the total dissipation rate $\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}$. The dissipation rate in each cell is computed using the method described in appendices A and B. We follow previous work (Uritsky et al., 2010) and define a connected dissipation structure as a connected set of cells where

with $\lambda$ a parameter we use to tune the detection threshold, $\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}$ the dissipation rate determined by our method and $\sigma_{\mathrm{\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}}}$ the standard deviation of the dissipation rate distribution. We choose $\lambda=4$ because we find that energy transfers are mainly due to events above $4\sigma$: we checked that the bulk of the third order structure function (responsible for energy transfers) is obtained from increments above 3-4 sigma. We also want the structure to be identifiable as clearly as possible and we expect such high dissipation structures to be associated with more intense gradients and a more clear-cut physical nature.

As already hinted at by local gradients (figure 2), we see on figure 3 that extracted dissipation structures are mainly sheets. Another way to see this is to look at a thin slice of the dissipation field in our OT simulation with $\mathcal{P}_{m}=1$ (figure 4) where the trace of the sheets appears as thin ridges. Compared to the same figure for the incompressible runs of Momferratos et al. (2014), the viscous and Ohmic natures of dissipation are now much more entangled and sometimes even overlap. A close eye inspection of this figure (and of similar cuts at other time steps and initial conditions) reveals various sub-layering of Ohmic dissipation sheets (red striations or ohmic dissipation wrapped by shear) or isolated viscous and ohmic heating sheets (purple color, which hints at a mix of compressive viscous heating and ohmic heating). These are not the only situations which occur, but it reveals that the intense dissipation sheets are not always randomly positioned with respect to one another.

A careful inspection of figure 3 allows to witness a few small filament-like structures. Some of them may be traced on figure 2 by the low-probability tail in the bottom right hand corner of the right panel, where the aspect ratios of the gradients are such that $\ell_{\mathrm{scan}}\simeq\ell_{\mathrm{\perp 1}}$ while $\ell_{\mathrm{\perp 1}}>>\ell_{\mathrm{\perp 2}}$. These tube-like structures will unfortunately be missed by our systematic investigation which focuses on locally planar structures, but we checked a posteriori that these structures only account for a very small fraction of the dissipation (less than one percent).

Figure 5 shows how the dissipation is distributed in volume: it gives the volume filling factor of the regions of large dissipation as a function of their dissipation fraction. This figure compiles several time steps up to $t=1.33t_{\mathrm{turnover}}$, where $t_{\mathrm{turnover}}=L/\sqrt{<u^{2}>}=2\pi$ is the initial eddy turnover time. It shows that the intermittency of the dissipation decreases over time, as the r.m.s. sonic Mach number decreases, because we consider decaying turbulence simulations (figure 1 shows the rapid decline of the sonic Mach number). We see here that structures shown on figure 3 (yellow lines for dissipation greater than four standard deviations above the mean) occupy $\simeq 0.8\%$ of the volume while they are at the origin of $\simeq 25\%$ of the total dissipation rate.

As we use decaying simulations, the Mach number decreases rapidly as well as the dissipation. We have considered snapshots at two characteristic times. The first snapshot is at 1/3 turnover time, when the first large dissipative structures form, and shortly after the dissipation peak. It is customary to think of the dissipation peak as a point similar to a steady state because the time derivative of the dissipation is zero. However, as we will show, this epoch still bears a strong imprint from the initial conditions. We therefore also consider a second snapshot at one turnover time. We do not consider much later times, as turbulence quickly decays and r.m.s. Mach numbers become much lower than at the beginning (see figure 1).

We develop here the details of our procedure to identify scans along the connected structures.

Rankine-Hugoniot (RH) relations express jump conditions across discontinuities in their stationary frame (Macquorn Rankine, 1870; Gurnett & Bhattacharjee, 2005). RH relations hold in a very specific situation where the fluid is stationary, with a plane parallel symmetry and homogeneous conditions on either side of a discontinuity. Nothing seems further away than our fully turbulent decaying turbulence simulations. Nevertheless, we want to check if our structure identification allows to recover some of the properties expected from the RH relations. If they hold, it would bring more weight to the selection criteria we have devised, and it would generalise to 3D MHD the results of Lesaffre et al. (2020) who found that in 2D decaying unmagnetised turbulence, 1D steady-state shocks could be used to model the strongest dissipation structures. In 1D steady-state isothermal MHD, conservation of mass, momentum and magnetic field read

where $\left[\right]_{\mathrm{pre}}^{\mathrm{post}}$ denotes the difference between the states at pre- and post- discontinuity. $\@vec{n}$ is the normal to the discontinuity. In our study, we take $\@vec{n}=\hat{\@vec{r}}_{\rm scan}$, which contains most of the gradient for highly dissipative cells (see figure 2). In other words, the planar region hypothesis which subtends RH relations is well verified for the most intense dissipative regions.

Across the discontinuity, the velocity of the fluid transitions from above to under a characteristic speed set by the three MHD linear wave speeds ($c^{R}_{s,i,f}$, see section 2.4). This leads to the traditional MHD velocity regimes classification (Delmont & Keppens, 2011), where numbers designate up- and downstream states, and $u_{n}=\@vec{u.n}$:

• •

  1- superfast $u_{n}\geq c^{R}_{f}$

• •

  2- subfast/super-alfvénic $c^{R}_{i}\leq u_{n}\leq c^{R}_{f}$

• •

  3- sub-Alfvénic/superslow $c^{R}_{s}\leq u_{n}\leq c^{R}_{i}$

• •

  4- subslow $u_{n}\leq c^{R}_{s}$

The discontinuity type is labelled by $i\rightarrow j$, where $i\geq j$. These discontinuity types show different behavior for the transverse magnetic fields :

• •

  $1\rightarrow 2$ are fast shocks. Magnetic field is refracted away from the shock normal, which yields a transverse magnetic field increase. Fast shocks efficiently convert kinetic to transverse magnetic energy.

• •

  $3\rightarrow 4$ are slow shocks. Magnetic field is refracted toward the shock normal, which yields a transverse magnetic field decrease. Slow shocks are efficient at compressing the gas.

• •

  $1\rightarrow 3$, $1\rightarrow 4$, $2\rightarrow 3$, $2\rightarrow 4$ are intermediate shocks. The transverse magnetic field flips across the shock normal.

• •

  $2=3\rightarrow 2=3$ are called Alfvén discontinuities or rotational discontinuities. The norm of the transverse magnetic field is unchanged between pre- and post-discontinuity regions, and only its direction changes in the plane parallel to the discontinuity. Alfvén discontinuities are believed to be efficient at reconnecting the field lines (Zweibel & Brandenburg, 1997; Zhdankin et al., 2013).

Density and total pressure profiles also show different signatures. In the first three cases, these profiles are jumps whose amplitude depends on the parameters of the shock. In the case of Alfvén discontinuities, these quantities must be identical on both sides of the discontinuity.

We consider here the effect of the initial conditions of the simulations on the nature of the structures formed. Figure 12 shows distributions of the identification scans as a function of the mean dissipation rate in the volume probed by each scan, and table 2 summarises these results averaged over all scans. The time step chosen for the two graphs on the left half of this figure (as well as for the left half of table 2) is 1/3 of the initial turnover time, shortly after the dissipation peak, when the first and most intense dissipation structures form. The time step chosen for the right hand side of figure 12 and table 2 is after one turn-over.

As described in section 2.1 we use two types of initial flows, ABC and OT. The main difference between these two flows resides in their magnetic and cross helicities (see subsection 2.1.2). This initially yields very different types of structures. Early time, OT is dominated by shocks (mainly fast shocks) while ABC is dominated by Alfvén discontinuities (rotational discontinuities and Parker sheets). After one turnover time, the impact of the initial conditions on the formation of the dissipation structures seems to be erased. The main dissipation mechanism is then through rotational discontinuities and Parker sheets for both ABC and OT.

Interestingly, for both types of initial conditions, at early and late times, the distribution of physical natures of scans does not seem to depend on their level of dissipation. Intense dissipative scans and weak scans have about the same proportions of each nature.

Table 2 also shows the amount of unidentified and misidentified scans. Between 58% and 78% of intense dissipation is identified by our technique, with a greater success rate for ABC than OT. Early time structures are also better identified than later times.

One critical parameter of dissipation is the magnetic Prandtl number, $\mathcal{P}_{m}=\nu/\eta$, the ratio of kinematic viscosity $\nu$ to magnetic diffusivity $\eta$ (see e.g. Brandenburg & Rempel, 2019). We perform our dissipation structure analysis on simulations with a range of magnetic Prandtl numbers, from $\mathcal{P}_{m}=1$ to $\mathcal{P}_{m}=16$ (see table 1). As dicussed in appendix B, the intrinsic numerical dissipation from the scheme causes the effective Prandtl number to be slightly different from the input one. Thanks to semi-analytical solutions (computed in appendix A) we probe the effective Prandtl number of our scheme in 1D MHD shocks. Figure 24 shows that for moderate velocity shocks (or $u_{0}\leq 1$) our scheme is already converged while the highest shock velocities we find in simulations are at $u_{0}\simeq 5$. We thus remain confident that the effective Prandtl number in our simulations is overall close to the input one, at least as regards shocks.

Figure 13 shows OT and ABC identification distributions for input $\mathcal{P}_{m}=1,4,16$: the magnetic Prandtl does not seem to have any impact on the distribution of structures. The only noticeable difference is a slight increase in the number of rotational discontinuities at the expense of Parker sheets for ABC initial conditions.

It was shown by Brandenburg & Rempel (2019) that an increase in $\mathcal{P}_{\mathrm{m}}$ causes an increase in $<\epsilon_{v}>/<\epsilon_{\eta}>$, a result which we confirm and make more precise here. If there is no statistical difference in high dissipation structures distributions, it must be differences in the internal structure of the dissipation layers that leads to differences in dissipation rates. On figure 14 we show scatter plots of viscous versus Ohmic dissipation rates, where each dot marker expresses the mean of the corresponding dissipation rate within each scan. At $\mathcal{P}_{m}=1$, it is clear that rotational discontinuities and Parker sheets are dominated by magnetic energy dissipation. Fast shocks are an intermediate case, with a more balanced share between Ohmic and viscous dissipation rates. Slow shocks are dominated by viscous dissipation. Each type of structure is more or less characterised by a given slope in these graphs (i.e. the ratio $<\varepsilon_{v}>_{\mathrm{scan}}/<\epsilon_{\eta}>_{\mathrm{scan}}$ is within a more or less well defined sector for each type of structure, regardless of initial conditions OT or ABC). In particular for shocks (both slow and fast), when $\mathcal{P}_{m}$ varies, this ratio simply scales as $\mathcal{P}_{m}$ : the behaviour of dissipation within fast shock scans follows the rule $<\varepsilon_{v}>_{\mathrm{scan}}/<\epsilon_{\eta}>_{\mathrm{scan}}\propto\mathcal{P}_{m}$. This scaling is less clear for the other types of structures: Alfvén discontinuities experience a wider range of ratios at fixed $\mathcal{P}_{m}$ which makes it less easy to assess if such a scaling is present (in fact, the envelope of green and cyan points suggests a different scaling, closer to $\mathcal{P}_{m}^{1/2}$). Provided the global dissipation is reflected by intense events, this could explain why the global average ratio $<\varepsilon_{\nu}>/<\varepsilon_{\eta}>$ is also found to scale approximately as $\mathcal{P}_{m}$ in our simulations. We therefore conclude that in our simulations the increase in the viscous over Ohmic fraction when $\mathcal{P}_{m}$ rises comes not from a difference in the nature of the dissipation structures that form, but from a modification in the way each of these types of structure dissipates internally.

In this section we consider values of the state variables at the pre- and post- positions on either side of the discontinuities, and we look for statistical differences between the various natures of discontinuities. The distributions of entrance sonic Mach numbers (see figure 15) show without surprise that the entrance velocities are very small for Alfvén discontinuities, moderate for slow shocks and on the order of the r.m.s. Mach number for fast shocks, but with a wide spread distribution. Previous work by Lehmann et al. (2016) has shown the distributions of fast and slow shocks display exponential tails at large velocities. We have less statistics but our data is also consistent with this picture. The bottom row of figure 15 displays the distributions for the entrance orthogonal Alfvénic Mach number $\mathcal{M}_{\mathrm{a}}=\tilde{u}_{\mathrm{scan}}/(B_{\mathrm{scan}}/\sqrt{4\pi\rho})=\tilde{u}_{\mathrm{scan}}/c_{i}^{R}$. Naturally, it is above 1 for fast shocks, and below 1 for slow shocks. Its distribution for Alfvén discontinuities is more surprising, though, with wide spread values ranging to values even above 1, while one would expect it to be close to zero. This comes from the fact that the entrance velocities in these discontinuities are on the order or below the sound speed, but the magnetic field happens to be almost transverse, thus yielding very small values for the intermediate speed $c_{i}^{R}$. Finally, we also looked at the statistical distributions of the density on the pre-discontinuity side, and found these distributions were independent of the nature of the discontinuity considered.

The environment of each discontinuity is defined by 7 state variables on either side (pre- and post-) of the discontinuity, a total of 14 independent state variables. We can reduce this number by using the 7 conservation relations (mass, momentum and magnetic field, -7 independent variables), a normalization by the pre-discontinuity density (-1 variable), a rotation of the transverse axes to cancel one component of the pre-discontinuity magnetic field (-1 variable), as well as a transverse boost of the frame to cancel the transverse velocity components pre-discontinuity (-2 variables). The environment can thus be fully characterized by a remainder of 3 independent variables, which can all be considered on the pre-discontinuity side.

For these three independent ’entrance parameters’, we choose the normal $p_{1}$ and transverse $p_{2}$ magnetic pressures as well as the difference $p_{3}$ between the ram pressure and the normal magnetic pressure (all are normalised by the thermal pressure $p=\rho c^{2}$). This choice makes it easy to predict where each discontinuity should lie in the 3D parameter space, according to the sign of $p_{3}$: negative for slow shocks, zero for Alfvén discontinuities and positive for fast shocks, while $p_{1}$ and $p_{2}$ could be arbitrary positive numbers.

Figure 16 displays two projections of this 3D space onto the planes $(p_{3},p_{1})$ and $(p_{3},p_{2})$. In this parameter space the only forbidden region is set by the positivity of the ram pressure $\frac{B_{\mathrm{scan}}^{2}}{4\pi}\geq-(\rho u_{\mathrm{scan}}^{2}-\frac{B_{\mathrm{scan}}^{2}}{4\pi})$, which translates as $p_{1}\geq-p_{3}$. However, we find it is not the only region that is devoid of points: the entrance parameters of our structures do not explore fully the available parameter space. In fact, tight correlations for fast shocks are visible on the right hand panel of figure 16 (where $\rho u_{\mathrm{scan}}^{2}\simeq\frac{B_{\mathrm{scan}}^{2}}{4\pi}+0.5\frac{B_{\perp}^{2}}{4\pi})$ and similarly for the slow shocks on the left hand panel (where $\rho u_{\mathrm{scan}}^{2}<<\frac{B^{2}}{4\pi}$). It also appears that rotational discontinuities and Parker sheets are preferably parallel to the magnetic field ($B_{\mathrm{scan}}<<B_{\perp}$). The latter two discontinuity types are not distinguishable in this parameters space, as we consider here only the pre- and post-discontinuity states without taking into account the internal structure.

These statistical constraints on the three entrance parameters of the dissipation structures reduce to two the number of independent parameters. It will be the subject of future work to understand the origin of these correlations in MHD turbulence.

Molecular line observations probe the radial velocities across the plane of sky, while we have access to the full three dimensional geometry of the intense dissipation regions in our simulations. When projected on the plane of sky, an intense dissipation sheet yields a salient filament-like feature on observable maps where the plane of the sheet is orthogonal to the plane of sky, i.e. where column-density is greatly enhanced through a caustic-like effect. On figure 17 we hence display the velocity difference statistics projected on the two transverse directions $\hat{\@vec{r}}_{\perp 1}$ and $\hat{\@vec{r}}_{\perp 2}$, as a proxy to what an observer would measure for the velocity difference across the projected ridge of an intense dissipation sheet, in the two cases where the line of sight is along $\hat{\@vec{r}}_{\perp 1}$ or $\hat{\@vec{r}}_{\perp 2}$. As expected due to rotation symmetry, rotational discontinuities have no noticeable difference depending on the direction, while fast and slow shocks are co-planar, hence the difference is greater in direction $\hat{\@vec{r}}_{\perp 1}$ than in direction $\hat{\@vec{r}}_{\perp 2}$. What is more surprising though is the fact that Parker sheets follow the same trend as rotational discontinuities: this is an indication that there is more continuity between the rotational discontinuity and Parker sheet classes than what our arbitrary division between the two would suggest. An interesting feature we also observe is the bimodality of the slow shocks compared to the fast shocks, which is linked to the fact that $B_{\mathrm{scan}}/\sqrt{4\pi\rho}$ needs to be greater than the typical velocity to yield a slow shock.

We retained the sign of the velocity differences although they technically cannot be probed by observations, due to the unknown projection angle. In the OT case, the positive and negative velocity differences for slow and fast shocks along $\hat{\@vec{r}}_{\perp 1}$ have markedly different statistics. This is due to the relative orientation between the fluid velocity and the magnetic field direction along the scanning vector. When they have the same orientation (i.e. $B_{\mathrm{scan}}=B_{n}>0$), transverse velocity differences have the same sign as transverse magnetic field differences (see equation 29): in this case fast and slow shocks have respectively positive and negative velocity differences. For the OT initial conditions, the positive cross helicity results from a mean positive alignment between velocity and magnetic fields, hence a more likely positive velocity difference for fast shocks, or negative for slow shocks. For the ABC initial conditions, the mean cross helicity is zero, which yields symmetric statistics.

We performed simulations at half the resolution (512³ pixels) in order to check the stability our results. Note that the dimensions of our scanning cylinders are defined with respect to the pixel size, with 3 pixels lateral radius and a 6 $\ell_{\varepsilon}$ scanning length centered on the detected local maxima of dissipation. The appropriate quantity to consider is hence the average energy dissipation rate per unit surface, or the energy flux through the surface of the discontinuity.

We consider on figure 18 the statistics of these dissipation fluxes nature by nature for two corresponding runs at 512³ and 1024³. They are seen to match perfectly, except for statistical noise which jags the lower resolution results. Indeed, we identify approximately four times less scans at low resolution, which is another indication that our dissipation structures are mainly sheets.

We also find that for both $N=512$ and $N=1024$, $\ell_{\varepsilon}$ is on the order of 1.5 pixels. This is consistent with our findings on 1-D shocks in appendix B that a twice lower resolution would yield an energy deposit twice more wide spread in the scanning direction, while keeping its integrated value constant (thanks to our method to recover numerical dissipation). It is also a hint that the same behaviour holds for Alfvénic discontinuities, which was not obvious.

Finally, this is an other indication that large scale dynamics set up the environmental characteristics of the discontinuities (the values of the state variables of the gas on either side of them), while the microphysics (physical and numerical dissipation) control the internal profiles of these discontinuities. The first evidence for this was uncovered with our Prandtl number study in section 4.2.

Previously (section 2.1.2, for example), we discussed the relative distributions within our scans. However, it is unfortunately difficult to relate it to the global dissipation in the computational domain, because the scans focus on the most intense events.

Nevertheless, we have seen in section 3.1 that strong dissipative pixels considered in this study ($>4\sigma$) represent a large fraction ($\simeq 25\%$ for ABC near the dissipation peak and $\simeq 30\%$ for OT) of the global dissipation rate of the simulation time step. However, the dissipation rate exceeds this threshold only near the peak of each scan, so the dissipation in a given scan also accounts for some dissipation below that threshold. Hence, the dissipation within our scans must amount to more than these global fractions of dissipation.

Furthermore, if we had chosen a lower threshold to identify structures, we would have detected weaker scans closer to the lateral edges of the dissipation structures. Since the proportions in physical natures do not appear to depend much on the strength of the scan, we might expect to find the same proportions in these weaker scans. As a result the fractions we currently measure might apply to a more significant fraction of the global dissipation, although it is difficult to ascertain it (overlap between scans of neighbouring structures and lack of planarity for some of the weaker scans might moderate the above arguments). In particular, we are not able to assess whether there is a diffuse component of dissipation, outside the intense dissipative regions.

Our method identifies a large majority of the scans we probed in each simulation and timesteps studied. The criteria for the identification of the scans are kept in this study voluntarily strict, with the aim to discover the structures basis causing the dissipation in the isothermal compressible MHD regime. Although restricted to the purest structures, we identify a significant fraction of the total dissipation of the simulations. We mentioned in the previous section 5.2 the difficulty to relate total dissipation rates to the scan distribution. Here, we attempt to link the distributions of scans to the distribution of the pixels above a given threshold. To partly remedy to overlapping problems which might occur between scans, we flag cells above a given threshold of the simulation according to the first identification that contain them. The resulting dissipation rate identification is shown on the top plots of figure 19.

The black line shows the fraction of dissipation captured by the threshold four standard deviations over the mean. The colors below show the fraction of each pixels above this threshold identified as each of the four main natures we find in our scans (the white space combines pixels which were never flagged because they always fall in unidentified or unknown scans). This time evolution graph clearly shows how the distributions differ at early time for our two different initial conditions, and how they stabilise after little bit less than one turnover time. This confirms the result that we found for scan distributions in section 4.1.

Furthermore, if we consider only the well identified scans, we observe that about 70% of the related connected structures are identified by a single type of scan and 80% have more than 75% of their scans identified by the same nature (see figure 20). We can therefore consider propagating the dominant type of a connected structure to the remainder of its cells. This allows to avoid the edge effects of structures and increases the identified dissipation fraction. The result is shown on the bottom plots of figure 19. This graph tells us, first of all, that the fully unidentified connected structures, although representing a significant fraction of the studied structures in number, participate very little in the total dissipation of the cube. These are small fragmented events which also comprise the short filament-like structures seen on figure 3. Second, we notice that the unidentified scans often belong to structures dominated by rotational discontinuities, except for the OT simulations at an early time, when they are sometimes part of fast shocks (see figure 20). For ABC runs, until about 0.7 turnover times, the dissipation generated by the Parker sheets decreases slightly to the benefit of that produced by the rotational discontinuities (bottom right panel of figure 19). This implies that despite the uniqueness of the identifications within a related structure, a significant fraction of the Parker sheets are found within structures formed by rotational discontinuities, which relates to our previous remarks on the continuity between Parker-Sheets and rotational discontinuities: our divide between the two is rather arbitrary and these connected structures could probably be gathered into a single Alfénic discontinuity class.

We also tried to decrease the detection threshold of the dissipation structures to $\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}=<\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}>+2\sigma_{\mathrm{\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}}}$ to increase the fraction of the total identified dissipation. In this case, the identification rate decreases only by a few percent compared to a threshold of $\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}=<\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}>+4\sigma_{\mathrm{\varepsilon_{\mathrm{tot}}^{\mathrm{corr}}}}$. The contribution of the different natures to the overall dissipation remains similar and the fraction of dissipation above threshold identified by the scans increases by about 10% overall.

By using two sets of criteria which are rather independent, we have biased our identifications towards more false negative and less false positives. Thus, there remains many scans misidentified because they either do not fit any of our heuristic criteria (unidentified scans) or because the two sets of criteria do not match (misidentified scans). We list here some of the reasons why our identification criteria might miss a significant fraction of the scans.

The main culprits are ”edge” scans. These are scans at the periphery of structures where the main direction of the gradient is less well defined and therefore the scanning direction is less relevant. For instance, the scanning direction is irrelevant in the case of the small filament-like structures observed in figure 3 and 20, where scans probably fall at least in the misidenditified category. Also, when two structures are too close to each other, the heuristic part of the identification is confused, because bumps or jumps are less well defined. Note the wave decomposition suffers less from adjacent structures, because it is sensitive only at the cell scale. Some of the unidentifications could also be due to the presence of intermediate shocks (see section 5.6 below) but we reckon they probably account for a small fraction only.

Given the strong correlation between our two sets of criteria for identification (heuristic and ideal waves), one could suggest to use only ideal wave decomposition to greatly increase the identification rate. We would thus reach 100% of identification, but our results would then be subject to caution, and would be biased toward false positives. Indeed, one should restrict this wave only decomposition to the most planar cells where the gradient approach makes sense. Also, the identification would then rely on the velocity regimes, which we have shown can be subject to caution depending on the method to estimate the travelling speed of the discontinuities.

If intermediate shocks are present in our simulations, our heuristic criteria would voluntarily miss them, and they would fall in the unidentified category. Indeed, these shocks have either a density or a pressure jump, but often display a magnetic field trough. We chose not to add this criterion here because we would not have had an independent criterion on gradients to solidify this heuristic one. This might be a reason why we get less identification in the OT case at early times, which seems more prone to generate shocks. We have in fact attempted to target intermediate shocks, and have found some convincing cases. However, the uncertainty on our estimate for the steady state velocity of these discontinuities makes it difficult to validate the speed regimes of these shocks on a statistically significant population. We hence decided to postpone our investigation on these shocks. In any case, the fraction of unidentification that we publish here puts an upper limit on the fraction of intermediate shocks.

Some astrophysical situations (a solar coronal ejection, a supernova, or a runaway star encountering a cloud, for example) can locally inject mechanical energy in a short amount of time. In these events, well defined initial conditions can suddenly be imposed, which will later develop into turbulence. In these contexts, decaying turbulence is probably a sensible approach. However, one can question the applicability of decaying turbulence in many other realistic astrophysical situations. First, turbulence being ubiquitous, it is very unlikely that an initial set up would be entirely devoid of turbulent perturbations at the onset. Also, the fast decay of turbulence has led Stone et al. (1998) amongst other authors to advocate the necessity of persistent driving at the scales relevant for the interstellar clouds. Our experiments have shown that initial conditions can imprint significant changes in the early phases of the development of turbulence: similarly, we expect that a given forcing may impact at least to some degree the resulting turbulent cascade, and it will do so at all times. One should therefore be careful to properly characterise the properties of the forcing (such as helicity injection) used in the experiments of driven turbulence.

In particular, a lot of attention has been devoted to the importance of solenoidal vs. compressible forcing (Federrath et al., 2010). This has led to techniques to probe observationally the degree of compressibility of the gas: some regions of the ISM near the galactic center were found to be dominated by solenoidal driving (Federrath et al., 2016) while others associated to stellar feedback were found to be more compressively driven (Menon et al., 2021). Intermediate situations were also witnessed in Orion B (see Orkisz et al., 2017, for example). These considerations are important, as the nature of star formation is believed to be strongly influenced by the driving of the turbulence (Federrath & Klessen, 2012). In the present study, we have not attempted to check the influence of compressible vs. solenoidal initial velocity fields (our initial velocity fields are divergence free), but one can imagine that initially more weight would be given to shocks had we started with velocities including a compressive component. We still believe though that the tendency to evolve towards incompressible structures would persist (imagine starting the origin of times after a fraction of the turnover time: this might be an illustration of what could happen if we start the simulation with compressible initial conditions). Finally, an intermediate solution would be to start decaying turbulence with a snapshot of fully established driven turbulence, or to imprint on the initial conditions synthetic models of turbulence such as designed in Durrive et al. (2020).