numerical confirmation, see, e.g., Kriel et al., 2022, 2023). Figure 14: The
magnetic $E_{\rm mag}(k)$, compressible ($E_{\text{kin},\parallel}(k)$; top)
and incompressible ($E_{\text{kin},\perp}(k)$; bottom) kinetic energy ratios
as a function of $k$ for $\text{P}\nu=\infty$ (no bulk viscosity; left column)
and $\text{P}\nu=0.01$ (strong bulk viscosity; right column), coloured by
$t/t_{0}$ (indicated in the top left panel). We start at $t/t_{0}=30$, at the
end of the kinematic stage (black), and evolve to $t/t_{0}=90$, well within
the saturated stage (aqua; see Figure 5 for how $t/t_{0}$ corresponds to the
different stages). Each panel has the equipartition $E_{\rm
mag}(k)=E_{\text{kin}}(k)$ drawn in red. Compressible modes (top row): for
$\text{P}\nu=\infty$, the saturated state for the compressible modes is in
subequipartition ($E_{\rm mag}(k)/E_{\text{kin},\parallel}(k)<1$) at low and
high-$k$. The high-$k$ mode subequipartition energy is associated with the
sub-shear-viscous range, outlined in Figure 1. However, for
$\text{P}\nu=0.01$, the sub-shear-viscous range is destroyed through the bulk
viscosity, and only the low-$k$ modes are left in subequipartition, meaning
that all but the lowest compressible $k$ modes ought to be dominated by the
effects of the magnetic field. Incompressible modes (bottom row): during the
kinematic stage, $E_{\rm mag}(k)\sim k^{2}E_{\text{kin},\perp}(k)$ up to the
peak scale of the spectral energy ratio, previously shown in Beattie et al.
(2023). As the dynamo saturates, the energy ratio flattens to roughly $E_{\rm
mag}(k)\sim kE_{\text{kin},\perp}(k)$ for superequipartition $k$ modes
($E_{\rm mag}(k)/E_{\text{kin},\perp}(k)<1$), potentially due in part to
$\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})$
suppressing the field-line stretching term in Equation 35.

##  7 Magnetic energy spectra

We plot the magnetic spectra normalised by the integral, and averaged over
time, in Figure 11, showing the kinematic stage spectra in the top panel, and
the saturated spectra in the bottom panel, along with the ratio between the
spectra in the inset plot. Firstly, it is apparent that the magnetic fields
for the low $\text{P}\nu$ simulations are becoming more large-scale (high-$k$
modes are being suppressed), during both the kinematic stage, where they
follow the Kazantsev (1968) $k^{3/2}$ power law, and in the saturated stage.
We can quantify this more precisely by computing the magnetic integral (or
correlation) scale. For $E_{\rm mag}(k)$, it is

$\displaystyle\frac{\ell_{\text{cor},b}}{\ell_{0}}=\frac{2}{L}\left(\frac{\displaystyle\int\operatorname{d}\\!{k}\,k^{-1}E_{\rm
mag}(k)}{\displaystyle\int\operatorname{d}\\!{k}\,E_{\rm
mag}(k)}\right)^{-1},$ (58)

which we compute for each time realisation, for each of our simulations, and
plot across time in Figure 12, fitting cubic splines to the data to extract
the general low-frequency trends. We see the well-established trend that the
field is small-scale in the kinematic stage, before becoming correlated on
large-scales, reducing $\ell_{\text{cor},b}/\ell_{0}$ through the non-linear
and finally settling on $\ell_{\text{cor},b}/\ell_{0}\approx 2.5$ in the
saturated stage (Schekochihin et al., 2004c; Rincon, 2019; Seta & Federrath,
2021; Galishnikova et al., 2022; Beattie et al., 2023). The energy
redistribution can be seen on a mode-by-mode basis in the inset plot in the
bottom panel of Figure 11, where there is now a factor of $\approx 2$ more
energy concentrated at low-$k$, and an order of magnitude less at high-$k$.
For $\text{P}\nu<1$, the magnetic fields lose energy on high-$k$, potentially
through the suppression of the sub-shear-viscous modes and compressible mode
bottleneck, which could in principle grow the field through compressions, but
these will be highly-inefficient since flux freezing will be violated strongly
on these scales (resistivity becomes comparable to induction at high-$k$), at
least in our simulations where a sub-shear-viscous range is generated. We
annotate the rough extent of this range with red in Figure 11. Furthermore, as
expected from Figure 11, we find that for $\text{P}\nu<0.01$, where we know
the compressible/incompressible mode coupling is strong, the correlation
scales are pushed closer to $\ell_{0}$. This could be simply because
$\ell_{\nu_{\text{shear}}}$ shifts to lower $k$ modes through the coupling,
which in turn reorganises the magnetic field for a fixed
$\text{Pm}_{\text{shear}}$, since
$k_{\eta}/k_{\nu_{\text{shear}}}\sim\text{Pm}_{\text{shear}}^{1/2}$
(Schekochihin et al., 2002c; Kriel et al., 2022, 2023), and as we have
discussed in Section 6.2, $k_{\nu_{\text{shear}}}\sim k_{\nu_{\text{bulk}}}$
for low $\text{Re}_{\text{bulk}}$. The key takeaway is that the bulk viscosity
suppresses magnetic energy in the high-$k$ modes, either due to the
suppression of the sub-shear-viscous range, where compressible modes dominate
and potentially amplify the field through weak compressions, or through direct
correlations between $k_{\nu_{\text{shear}}}$ and $k_{\nu_{\text{bulk}}}$, in
both the stages of the dynamo. We annotate how this may work in Figure 13,
where we show show $k_{\nu_{\text{bulk}}}$ shifting to lower $k$ modes, which
in turn moves $k_{\nu_{\text{shear}}}$, which in turn shifts $k_{\eta}$,
destroying the high-$k$ magnetic modes. The net effect of this process is that
magnetic fields are significantly more large-scale than the pure shear-viscous
dynamo counterparts.

##  8 Spectral ratio functions

The final statistic we use to probe the role of bulk viscosity in the
compressible turbulent dynamo are the ratios between the magnetic and
compressible and incompressible kinetic energies as a function of time, Figure
14, revealing scale-by-scale information about the transition between the
stages of the dynamo (see Section 1), as previously done in Beattie et al.
(2023). To further build upon Beattie et al. (2023), we split the saturation
into compressible $E_{\text{kin},\parallel}(k)/E_{\rm mag}(k)$ (top panels)
and incompressible $E_{\text{kin},\perp}(k)/E_{\rm mag}(k)$ (bottom panels),
with the columns organised by $\text{P}\nu$, with $\text{P}\nu=\infty$ (left)
and $\text{P}\nu=0.01$ (right). We take the spectral ratios from $t/t_{0}=30$,
which is within the kinematic stage, to $t/t_{0}=90$, which is well and truly
saturated (see Figure 5). The intermediate purple stage is within the
nonlinear, transition stage between the kinematic and saturation. We draw a
dashed, red line at equipartition, $E_{\text{kin}}(k)/E_{\rm mag}(k)=1$, which
defines an energy equipartition scale, $k_{\rm eq}$ such that
$E_{\text{kin}}(k_{\rm eq})/E_{\rm mag}(k_{\rm eq})=1$, which has gone by a
number of names in the literature, e.g., tearing scale (Galishnikova et al.,
2022, and references therein), or stretching scale (Schekochihin et al.,
2004c), to describe the onset of saturation (as shown in Equation 38).

Let us first focus on the left column for the $\text{P}\nu=\infty$
simulations. The bottom panel shows the onset of the $E_{\rm mag}(k)$
saturation w.r.t, the incompressible modes, which we can directly compare with
the incompressible, $\text{Pm}_{\text{shear}}=1$ simulations in Beattie et al.
(2023). As Beattie et al. (2023) showed previously, in the kinematic stage
there is an extended power law $E_{\rm mag}(k)\sim
k^{2}E_{\text{kin},\perp}(k)$ up to the peak of $E_{\rm
mag}(k)/E_{\text{kin},\perp}(k)$. Throughout the nonlinear stage, the modes
$k>k_{\rm eq}$ are suppressed, in that $E_{\rm mag}(k)\sim
k^{2}E_{\text{kin},\perp}(k)$ becomes significantly shallower, close to
$E_{\rm mag}(k)\sim kE_{\text{kin},\perp}(k)$. Based on our analysis in
Section 3.2, this must be when the field growth is suppressed when the
stretching is being suppressed, which is presumably happening on these scales.
As the dynamo saturates, $k_{\rm eq}$ progressively moves from the peak of
$E_{\rm mag}(k)/E_{\text{kin},\perp}(k)$, $k\approx 100$, to $k\approx 10$,
where it settles.
$\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}\otimes\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}:\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
should still be operational on these subequipartition scales, but now the
stretching is slow, since the stretching time scale $t_{k}\sim
1/(kv_{\perp,k})$ becomes larger towards lower $k$ modes (Schekochihin et al.,
2004b). In the top panel, we see the same spectral ratio shape and final
saturated state is not retained in the compressible modes, as we expect basdd
on our integral statistics in Section 5.

Firstly, whilst $E_{\rm mag}(k)/E_{\text{kin},\parallel}(k)$ retains similar
power-law like structure into the peak as $E_{\rm
mag}(k)/E_{\text{kin},\perp}(k)$, the peak of $E_{\rm
mag}(k)/E_{\text{kin},\parallel}(k)$ is shifted to much smaller modes
$k\approx 30$, and now a significant range of scales, the sub-shear-viscous
scales, emerge at high-$k$ with $E_{\rm
mag}(k)/E_{\text{kin},\parallel}(k)\approx 0.1$. Moreover, there is no clear
change in the shape of $E_{\rm mag}(k)/E_{\text{kin},\parallel}(k)$ as the
dynamo saturates, in contrast to the suppression we see in the bottom panel.
Naturally, there are two $k_{\rm eq}$, the low-$k$ $k_{\rm eq}$ corresponding
to the transition between compressible modes that are in subequipartition, and
hence strong enough to influence the field, and then the high-$k$ $k_{\rm
eq}$, corresponding to the onset of the sub-shear-viscous range. Sandwiched in
between are a range of scales that are in superequipartition, making the
scale-by-scale saturation quite complex in the compressible modes of the
$\text{P}\nu\gg 1$ plasmas.

Finally, we turn to the $\text{P}\nu=0.01$ simulations in the right panels of
Figure 14. For $E_{\rm mag}(k)/E_{\text{kin},\perp}(k)$, a similar
$k^{2}\rightarrow k$ transition is found through the evolution from the
kinematic to saturated state. From Figure 9, we showed that
$E_{\text{kin},\perp}(k)$ modes $k\gtrsim 10$ are influenced by viscosity
coupling. Compared to the $\text{P}\nu=\infty$ simulation, these modes are
suppressed, with $E_{\rm mag}(k)/E_{\text{kin},\perp}(k)\approx 2$ at the peak
energy scale in saturation. It is also clear that $k_{\rm eq}$ takes a slower
journey towards the low-$k$ modes, albeit saturating to a similar value, which
is probably the reason why that even at this low $\text{P}\nu$, the integral
statistics shown in Section 5 do not vary significantly from the
$\text{P}\nu=\infty$ case. Because $\text{P}\nu$ is small, the $E_{\rm
mag}(k)/E_{\text{kin},\parallel}(k)$ ratio is notably different from the
incompressible counterpart. No sub-shear-viscous range exists, and hence the
high-$k$ modes are highly magnetised with respect to the compressible modes,
finding values of $E_{\rm mag}(k)/E_{\text{kin},\parallel}(k)\approx 10^{3}$.
Furthermore, $k_{\rm eq}$ saturates to $k_{\rm eq}\approx 1$. Intuitively,
based on our previous findings, this means that all but the lowest
compressible $k$ modes are able to participate in the saturated state. As we
commented upon previously, the key take away is that both the journey towards
and the role in the saturation of the dynamo varies significantly for the
compressible and incompressible modes, with the incompressible modes

##  9 Summary and Conclusions

We perform three-dimensional, visco-resistive direct numerical simulations of
trans-sonic, magnetohydrodynamic turbulence undergoing magnetic field
amplification and saturation through the small-scale dynamo, including viscous
effects from both the shear and bulk viscosity. We list the conclusions below.

###  9.1 Bulk viscosity in the turbulent dynamo

To appropriately place bulk viscosity into a small-scale dynamo framework, we
derive and define two new dimensionless numbers to study the role of bulk
(volume) viscosity, the “bulk viscous Reynolds number"
$\text{Re}_{\text{bulk}}$, Equation 20, which is the ratio between Reynolds
stress and the bulk viscosity, as well as the “viscous Prandtl number"
$\text{P}\nu$, Equation 21, which is the ratio between the bulk and shear
viscosities (or equivalently, Reynolds numbers). For $\text{P}\nu<1$,
$\text{P}\nu$ defines the size of a sub-bulk-viscous range,
$k_{\nu_{\text{shear}}}>k>k_{\nu_{\text{bulk}}}$, and for $\text{P}\nu>1$, the
sub-shear-viscous range, $k_{\nu_{\text{bulk}}}>k>k_{\nu_{\text{shear}}}$,
which we illustrate in Figure 1. Intuitively, the bulk viscosity acts like a
Laplacian-like viscosity operator that primarily effects the compressible
modes, starting at high-$k$ and shifting to low-$k$ as $\text{P}\nu$
decreases.

In the $\text{P}\nu>1$ regime, which simulations without explicit bulk
viscosity will be in ($\text{P}\nu=\infty$), have a strong, high-$k$ build-up
of compressible modes that appear as a bottleneck in the total kinetic energy
spectrum. For $\text{Pm}_{\text{shear}}=1$ plasmas (Equation 25), where
$k_{\nu_{\text{shear}}}=k_{\eta}$, this results in low-levels of magnetisation
and high-levels of compressible kinetic energy building up on these scales.
Whether or not this is desirable behaviour for simulations of astrophysical
plasmas is therefore sensitive to the exact nature of bulk viscosity in the
plasma of interest, but in general, the nature of bulk viscosity is largely
unknown with a very limited amount of exploration in the literature both in
the laboratory and in simulation.

In the $\text{P}\nu<1$ regime, the compressible and incompressible modes
become coupled through the viscosity operators in the sub-bulk-viscous range.
This means that as $k_{\nu_{\text{bulk}}}$ shifts to low-$k$,
$k_{\nu_{\text{shear}}}$ shifts with it, which for fixed
$\text{Pm}_{\text{shear}}$ also brings $k_{\eta}$ to lower $k$ modes. When
$\text{P}\nu$ is very low, such that $t_{\nu_{\text{bulk}}}\sim t_{0}$, there
are significant changes in the shape of all energy spectra, dynamo integral
energy statistics, and of course, an extremely strong suppression of
compressible modes. If $\text{P}\nu$ is closer to one, the impact the bulk
viscosity has on the dynamo is only in the high-$k$ modes, and such, integral
energies, growth rates, and saturations are largely unaffected. Future work is
required to determine the relations for the scaling between
$k_{\nu_{\text{bulk}}}$, $k_{\nu_{\text{shear}}}$ and $k_{\eta}$ in the
$\text{P}\nu<1$ regime, where the coupling is strong.

###  9.2 Compressibility in the turbulent dynamo

To understand the role of compressible modes in the kinematic and saturated
dynamo stages, we derive the $1^{\rm st}$ moment magnetic energy equation that
separates out the effects of compressible and incompressible modes, in the
sense that they are truly orthogonal contributions to the growth rate and
saturation. This includes three main terms: volume-preserving stretching,
volume-preserving winding and volume-changing compressions (Equation 7–9,
respectively). Winding is orthogonal to the magnetic field tensor, hence one
cannot locally wind a turbulent magnetic field to grow it (see Equation 36).
By directly computing the growth terms (Figure 7), we find that stretching is
always the most dominant pathway for growing the magnetic field, even in the
supersonic dynamo, and irrespective of the phase of the dynamo. This explains
why $\ell_{\nu_{\text{shear}}}/\ell_{\eta}\sim\text{Pm}_{\text{shear}}^{1/2}$,
(the stretching rate at $\ell_{\nu_{\text{shear}}}$ balances with resistive
dissipation rate) is true for the compressible regime, which was recently
found in Kriel et al. (2023). This is in tension with some recent results
i.e., Sur & Subramanian (2023), but the tension is removed if they carefully
remove the trace of $\mathbb{S}$ in their decomposition (see footnote 4).
Compression plays a dual, almost symmetric (at least at $\mathcal{M}=1$) role,
in both growing and decaying the magnetic field through compressions and
dilatations. There is however a slight negative skew in the compression growth
term, which means on average the compressions oppose the growth, resulting in
slower growth rates and lower saturations than for the incompressible dynamo
counterparts. By splitting the velocity filed into compressible and
incompressible modes, and analysing the saturation across $k$ space, we come
to the conclusion that at least for our parameter range, for the compressible
modes, it is only the very low-$k$ modes that are strong enough to influence
the field through the late-kinematic stages and saturation. However, it is
also the low-$k$ modes that dominate the integral energies, hence why changing
the compressibility through various, large-scale parameterisations (changing
$\mathcal{M}$, changing the compressible mode energy at the driving scale) is
able to change the integral statistics and saturation energies in the
turbulent dynamo.

###  9.3 Beyond the fluid continuum limit

Our study is embedded within the context of fluid plasma theory, where
$\lambda_{\text{mfp}}/L\ll 1$, and $\lambda_{\text{mfp}}$ is the molecular
mean-free-path101010For the interstellar medium, the cold neutral medium has
an electron mean-free-path $\lambda_{e}\sim 5\,\text{R}_{\oplus}$, the warm
neutral medium $\lambda_{e}\sim 1.3\,\text{AU}$, the warm ionised medium
$\lambda_{e}\sim 0.17\,\text{AU}$ and the hot ionised medium $\lambda_{e}\sim
0.3\,\text{pc}$, with $L$ varying between 10 to 100 pc (Ferrière, 2020), hence
the interstellar medium really is very well-approximated by a fluid.. We have
studied the implications of bulk viscosity
$\sim\nabla\cdot(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}$
in $\mathcal{M}\approx 1$, turbulent dynamos, where shocks are able to
permeate through the fluid medium. However, for collisionless
$\lambda_{\text{mfp}}/L\gg 1$, weakly collisional $\lambda_{\text{mfp}}/L\sim
1$ and Braginskii (1965)-MHD the viscous properties of the plasma are
significantly different from the vanilla MHD fluid form (Section 2.1; St-Onge
et al. 2020; Squire et al. 2023). For example, in Braginskii (1965)-MHD, which
is important for plasmas in the intracluster medium (Berlok et al., 2020),
particle collisions are rare, in turn preventing the pressure tensor from
isotropising. This creates anisotropic viscous stress
$\sim(\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}\otimes\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}):\nabla\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
along $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$, taking the same form as our our
dynamo growth terms Equation 35, making the competition for dynamo growth and
viscous diffusion more fierce in this regime (St-Onge et al., 2020; Squire et
al., 2023). As Squire et al. (2023) points out, whether or not either of these
terms,
$\sim\nabla\cdot\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}\otimes\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}[(\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}\otimes\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}):\nabla\mathrm{{\mn@boldsymbol{\mathit{v}}}}]$
(Braginskii viscosity) or
$\sim\nabla\cdot(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}$
(bulk viscosity), are viscous dissipative or non-dissipative (i.e., akin to
the thermal response of $\nabla p$), depends strongly upon the aforementioned
plasma regime (see footnote on page 10 in Squire et al. 2023 for a longer
discussion). Interestingly, if we use our decomposition in Equation 35 for the
Braginskii-MHD viscosity, we find that there is a bulk-viscous-type term
$\sim(\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}\otimes\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}):\mathbb{B}$,
which will act to preferentially damp fast and slow magnetosonic waves along
$\mathrm{{\mn@boldsymbol{\mathit{b}}}}$. Furthermore, it has been seen in
controlled laboratory astrophysics experiments that kinetic instabilities of
the Weibel (1959)-type can play a critical role in the formation of
collisionless shocks (Fox et al., 2013; Huntington et al., 2015). Whether such
mechanisms can be parameterised at all in terms of viscosity coefficients
remains an open question. We leave the detailed exploration of these bulk
viscous-type operators in the context of these different plasma regimes for
future studies.

## Acknowledgements

We thank Bart Ripperda, Aditya Vijaykumar, Janosz Dewberry, Elias Most, Philip
Kempski, Amit Seta, Mark Krumholz, Chris Thompson and Drummond Fielding for
the helpful discussions. J. R. B. acknowledges financial support from the
Australian National University, via the Deakin PhD and Dean’s Higher Degree
Research (theoretical physics) Scholarships and the Australian Government via
the Australian Government Research Training Program Fee-Offset Scholarship and
the Australian Capital Territory Government funded Fulbright scholarship. C.
F. acknowledges funding provided by the Australian Research Council (Future
Fellowship FT180100495 and Discovery Project DP230102280), and the Australia-
Germany Joint Research Cooperation Scheme (UA-DAAD). A. B. and J. R. B. also
acknowledge partial support from NSF Award 2206756. We further acknowledge
high-performance computing resources provided by the Leibniz Rechenzentrum and
the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale
project 10391), the Australian National Computational Infrastructure (grant
ek9) and the Pawsey Supercomputing Centre (project pawsey0810) in the
framework of the National Computational Merit Allocation Scheme and the ANU
Merit Allocation Scheme.

The simulation software, flash, was in part developed by the flash Centre for
Computational Science at the Department of Physics and Astronomy of the
University of Rochester. We utilise the TurbGen code for the turbulent forcing
function, $\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ Federrath et al. (2022).
Data analysis and visualisation software used in this study: C++ (Stroustrup,
2013), numpy (Oliphant, 2006; Harris et al., 2020), matplotlib (Hunter, 2007),
cython (Behnel et al., 2011), visit (Childs et al., 2012), scipy (Virtanen et
al., 2020), scikit-image (van der Walt et al., 2014), cmasher (van der Velden,
2020).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Achikanath Chirakkara et al. (2021) Achikanath Chirakkara R., Federrath C., Trivedi P., Banerjee R., 2021, Phys. Rev. Lett., 126, 091103
  * Aguilera-Miret et al. (2023a) Aguilera-Miret R., Palenzuela C., Carrasco F., Viganò D., 2023a, arXiv e-prints, p. arXiv:2307.04837
  * Aguilera-Miret et al. (2023b) Aguilera-Miret R., Palenzuela C., Carrasco F., Viganò D., 2023b, Phys. Rev. D, 108, 103001
  * Batchelor (1950) Batchelor G. K., 1950, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 201, 405
  * Beattie & Federrath (2020) Beattie J. R., Federrath C., 2020, MNRAS, 492, 668
  * Beattie et al. (2019) Beattie J. R., Federrath C., Klessen R. S., Schneider N., 2019, MNRAS, 488, 2493
  * Beattie et al. (2020) Beattie J. R., Federrath C., Seta A., 2020, MNRAS, 498, 1593
  * Beattie et al. (2022a) Beattie J. R., Krumholz M. R., Federrath C., Sampson M. L., Crocker R. M., 2022a, Frontiers in Astronomy and Space Sciences, 9, 900900
  * Beattie et al. (2022b) Beattie J. R., Krumholz M. R., Skalidis R., Federrath C., Seta A., Crocker R. M., Mocz P., Kriel N., 2022b, MNRAS, 515, 5267
  * Beattie et al. (2022c) Beattie J. R., Mocz P., Federrath C., Klessen R. S., 2022c, MNRAS, 517, 5003
  * Beattie et al. (2023) Beattie J. R., Federrath C., Kriel N., Mocz P., Seta A., 2023, MNRAS, 524, 3201
  * Beck & Wielebinski (2013) Beck R., Wielebinski R., 2013, Magnetic Fields in Galaxies. p. 641, doi:10.1007/978-94-007-5612-0_13
  * Behnel et al. (2011) Behnel S., Bradshaw R., Citro C., Dalcin L., Seljebotn D. S., Smith K., 2011, Computing in Science & Engineering, 13, 31
  * Berlok et al. (2020) Berlok T., Pakmor R., Pfrommer C., 2020, MNRAS, 491, 2919
  * Bott et al. (2021) Bott A. F. A., et al., 2021, Phys. Rev. Lett., 127, 175002
  * Bouchut et al. (2010) Bouchut F., Klingenberg C., Waagan K., 2010, Numerische Mathematik, 115, 647
  * Braginskii (1965) Braginskii S. I., 1965, Reviews of Plasma Physics, 1, 205
  * Brandenburg & Ntormousi (2023) Brandenburg A., Ntormousi E., 2023, ARA&A, 61, 561
  * Brandenburg & Subramanian (2005) Brandenburg A., Subramanian K., 2005, Phys. Rep., 417, 1
  * Brandenburg et al. (1995) Brandenburg A., Procaccia I., Segel D., 1995, Physics of Plasmas, 2, 1148
  * Brandenburg et al. (2023) Brandenburg A., Rogachevskii I., Schober J., 2023, MNRAS, 518, 6367
  * Brzycki & ZuHone (2019) Brzycki B., ZuHone J., 2019, The Astrophysical Journal, 883, 118
  * Burgers (1948) Burgers J., 1948, Advances in Applied Mechanics, 1, 171
  * Burkhart (2021) Burkhart B., 2021, PASP, 133, 102001
  * Burkhart et al. (2009) Burkhart B., Falceta-Gonçalves D., Kowal G., Lazarian A., 2009, ApJ, 693, 250
  * Camelio et al. (2023) Camelio G., Gavassino L., Antonelli M., Bernuzzi S., Haskell B., 2023, Phys. Rev. D, 107, 103031
  * Chabanov & Rezzolla (2023) Chabanov M., Rezzolla L., 2023, arXiv e-prints, p. arXiv:2307.10464
  * Chabanov et al. (2023) Chabanov M., Tootle S. D., Most E. R., Rezzolla L., 2023, ApJ, 945, L14
  * Chen et al. (2019) Chen S., Wang X., Wang J., Wan M., Li H., Chen S., 2019, Physics of Fluids, 31
  * Childs et al. (2012) Childs H., et al., 2012, in , High Performance Visualization–Enabling Extreme-Scale Scientific Insight. Taylor & Francis, pp 357–372
  * Choudhuri (1998) Choudhuri A. R., 1998, The physics of fluids and plasmas : an introduction for astrophysicists /
  * Czajka et al. (2019) Czajka A., Dasgupta K., Gale C., Jeon S., Misra A., Richard M., Sil K., 2019, Journal of High Energy Physics, 2019, 145
  * Davidovits et al. (2022) Davidovits S., Federrath C., Teyssier R., Raman K. S., Collins D. C., Nagel S. R., 2022, Phys. Rev. E, 105, 065206
  * Dubey et al. (2008) Dubey A., et al., 2008, in Pogorelov N. V., Audit E., Zank G. P., eds, Astronomical Society of the Pacific Conference Series Vol. 385, Numerical Modeling of Space Plasma Flows. p. 145
  * Emanuel (1990) Emanuel G., 1990, Physics of Fluids A: Fluid Dynamics, 2, 2252
  * Federrath (2013) Federrath C., 2013, MNRAS, 436, 1245
  * Federrath (2016) Federrath C., 2016, Journal of Plasma Physics, 82, 535820601
  * Federrath & Banerjee (2015) Federrath C., Banerjee S., 2015, MNRAS, 448, 3297
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761
  * Federrath et al. (2010) Federrath C., Roman-Duval J., Klessen R., Schmidt W., Mac Low M. M., 2010, A&A, 512
  * Federrath et al. (2011) Federrath C., Chabrier G., Schober J., Banerjee R., Klessen R. S., Schleicher D. R. G., 2011, Phys. Rev. Lett., 107, 114504
  * Federrath et al. (2014) Federrath C., Schober J., Bovino S., Schleicher D. R. G., 2014, ApJ, 797, L19
  * Federrath et al. (2016) Federrath C., et al., 2016, ApJ, 832, 143
  * Federrath et al. (2021) Federrath C., Klessen R. S., Iapichino L., Beattie J. R., 2021, Nature Astronomy
  * Federrath et al. (2022) Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2022, TG: Turbulence Generator, Astrophysics Source Code Library, record ascl:2204.001 (ascl:2204.001)
  * Ferrand et al. (2020) Ferrand R., Galtier S., Sahraoui F., Federrath C., 2020, ApJ, 904, 160
  * Ferrière (2020) Ferrière K., 2020, Plasma Physics and Controlled Fusion, 62, 014014
  * Fox et al. (2013) Fox W., Fiksel G., Bhattacharjee A., Chang P.-Y., Germaschewski K., Hu S. X., Nilson P. M., 2013, Phys. Rev. Lett., 111, 225002
  * Frisch (1995) Frisch U., 1995, Turbulence. The legacy of A. N. Kolmogorov.
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Gaensler et al. (2011) Gaensler B. M., et al., 2011, Nature, 478, 214
  * Galishnikova et al. (2022) Galishnikova A. K., Kunz M. W., Schekochihin A. A., 2022, Physical Review X, 12, 041027
  * Gent et al. (2023) Gent F. A., Mac Low M.-M., Korpi-Lagg M. J., Singh N. K., 2023, ApJ, 943, 176
  * Grete et al. (2017) Grete P., O’Shea B. W., Beckwith K., Schmidt W., Christlieb A., 2017, Physics of Plasmas, 24, 092311
  * Grete et al. (2020) Grete P., O’Shea B. W., Beckwith K., 2020, ApJ, 889, 19
  * Grete et al. (2023) Grete P., O’Shea B. W., Beckwith K., 2023, ApJ, 942, L34
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Haugen et al. (2004a) Haugen N. E. L., Brandenburg A., Dobler W., 2004a, Phys. Rev. E, 70, 016308
  * Haugen et al. (2004b) Haugen N. E. L., Brandenburg A., Mee A. J., 2004b, MNRAS, 353, 947
  * Hew & Federrath (2023) Hew J. K. J., Federrath C., 2023, MNRAS, 520, 6268
  * Hopkins (2013) Hopkins P. F., 2013, MNRAS, 430, 1880
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Huntington et al. (2015) Huntington C. M., et al., 2015, Nature Physics, 11, 173
  * Jacquemin-Ide et al. (2023) Jacquemin-Ide J., Rincon F., Tchekhovskoy A., Liska M., 2023, arXiv e-prints, p. arXiv:2311.00034
  * Kazantsev (1968) Kazantsev A. P., 1968, Soviet Journal of Experimental and Theoretical Physics, 26, 1031
  * Kempski et al. (2023) Kempski P., Fielding D. B., Quataert E., Galishnikova A. K., Kunz M. W., Philippov A. A., Ripperda B., 2023, MNRAS, 525, 4985
  * Klessen (2011) Klessen R. S., 2011, in Charbonnel C., Montmerle T., eds, EAS Publications Series Vol. 51, EAS Publications Series. pp 133–167 (arXiv:1109.0467), doi:10.1051/eas/1151009
  * Kolmogorov (1941) Kolmogorov A. N., 1941, Doklady Akademii Nauk Sssr, 30, 301
  * Kriel et al. (2022) Kriel N., Beattie J. R., Seta A., Federrath C., 2022, MNRAS, 513, 2457
  * Kriel et al. (2023) Kriel N., Beattie J. R., Federrath C., Krumholz M. R., Hew J. K. J., 2023, arXiv e-prints, p. arXiv:2310.17036
  * Kritsuk et al. (2007) Kritsuk A. G., Norman M. L., Padoan P., Wagner R., 2007, ApJ, 665, 416
  * Krumholz & Federrath (2019) Krumholz M. R., Federrath C., 2019, Frontiers in Astronomy and Space Sciences, 6, 7
  * Krumholz et al. (2020) Krumholz M. R., Crocker R. M., Xu S., Lazarian A., Rosevear M. T., Bedwell-Wilson J., 2020, MNRAS, 493, 2817
  * Landau & Lifshitz (1959) Landau L., Lifshitz E., 1959, Fluid Mechanics: Landau and Lifshitz: Course of Theoretical Physics. Butterworth-Heinemann
  * Lehmann et al. (2016) Lehmann A., Federrath C., Wardle M., 2016, MNRAS, 463, 1026
  * Lemoine (2023) Lemoine M., 2023, Journal of Plasma Physics, 89, 175890501
  * Lithwick & Goldreich (2001) Lithwick Y., Goldreich P., 2001, ApJ, 562, 279
  * Malvadi Shivakumar & Federrath (2023) Malvadi Shivakumar L., Federrath C., 2023, arXiv e-prints, p. arXiv:2311.10350
  * Mandal et al. (2020) Mandal A., Federrath C., Körtgen B., 2020, Monthly Notices of the Royal Astronomical Society, 493, 3098
  * Marder (1987) Marder B., 1987, Journal of Computational Physics, 68, 48
  * Mathew & Federrath (2021) Mathew S. S., Federrath C., 2021, MNRAS, 507, 2448
  * McKee et al. (2020) McKee C. F., Stacy A., Li P. S., 2020, MNRAS, 496, 5528
  * Mee & Brandenburg (2006) Mee A. J., Brandenburg A., 2006, MNRAS, 370, 415
  * Mininni et al. (2005) Mininni P., Alexakis A., Pouquet A., 2005, Phys. Rev. E, 72, 046302
  * Mocz & Burkhart (2019) Mocz P., Burkhart B., 2019, ApJ, 884, L35
  * Moffatt & Saffman (1964) Moffatt H. K., Saffman P. G., 1964, The Physics of Fluids, 7, 155
  * Molina et al. (2012) Molina F. Z., Glover S. C. O., Federrath C., Klessen R. S., 2012, MNRAS, 423, 2680
  * Most et al. (2022a) Most E. R., Haber A., Harris S. P., Zhang Z., Alford M. G., Noronha J., 2022a, arXiv e-prints, p. arXiv:2207.00442
  * Most et al. (2022b) Most E. R., et al., 2022b, MNRAS, 509, 1096
  * Nolan et al. (2015) Nolan C. A., Federrath C., Sutherland R. S., 2015, MNRAS, 451, 1380
  * Oliphant (2006) Oliphant T., 2006, NumPy: A guide to NumPy, USA: Trelgol Publishing, http://www.numpy.org/
  * Orkisz et al. (2017) Orkisz J. H., et al., 2017, A&A, 599, A99
  * Padoan & Nordlund (2011) Padoan P., Nordlund Å., 2011, ApJ, 730, 40
  * Pakmor et al. (2014) Pakmor R., Marinacci F., Springel V., 2014, The Astrophysical Journal Letters, 783, L20
  * Pakmor et al. (2017) Pakmor R., et al., 2017, Monthly Notices of the Royal Astronomical Society, 469, 3185
  * Pan & Johnsen (2017) Pan S., Johnsen E., 2017, Journal of Fluid Mechanics, 833, 717–744
  * Park & Ryu (2019) Park J., Ryu D., 2019, ApJ, 875, 2
  * Passot & Vázquez-Semadeni (1998) Passot T., Vázquez-Semadeni E., 1998, Phys. Rev. E, 58, 4501
  * Pfrommer et al. (2022) Pfrommer C., Werhahn M., Pakmor R., Girichidis P., Simpson C. M., 2022, Monthly Notices of the Royal Astronomical Society, 515, 4229
  * Rincon (2019) Rincon F., 2019, Journal of Plasma Physics, 85, 205850401
  * Ruzmaikin & Sokolov (1981) Ruzmaikin A. A., Sokolov D. D., 1981, Pisma v Astronomicheskii Zhurnal, 7, 701
  * Sampson et al. (2023) Sampson M. L., Beattie J. R., Krumholz M. R., Crocker R. M., Federrath C., Seta A., 2023, MNRAS, 519, 1503
  * Schekochihin et al. (2002a) Schekochihin A. A., Cowley S. C., Hammett G. W., Maron J. L., McWilliams J. C., 2002a, New Journal of Physics, 4, 84
  * Schekochihin et al. (2002b) Schekochihin A. A., Boldyrev S. A., Kulsrud R. M., 2002b, ApJ, 567, 828
  * Schekochihin et al. (2002c) Schekochihin A. A., Boldyrev S. A., Kulsrud R. M., 2002c, ApJ, 567, 828
  * Schekochihin et al. (2004a) Schekochihin A. A., Cowley S. C., Maron J. L., McWilliams J. C., 2004a, Physical review letters, 92, 054502
  * Schekochihin et al. (2004b) Schekochihin A. A., Cowley S. C., Taylor S. F., Hammett G. W., Maron J. L., McWilliams J. C., 2004b, Phys. Rev. Lett., 92, 084504
  * Schekochihin et al. (2004c) Schekochihin A. A., Cowley S. C., Taylor S. F., Maron J. L., McWilliams J. C., 2004c, The Astrophysical Journal, 612, 276
  * Schmidt & Grete (2019) Schmidt W., Grete P., 2019, Phys. Rev. E, 100, 043116
  * Schmidt et al. (2008) Schmidt W., Federrath C., Klessen R., 2008, Phys. Rev. Lett., 101, 194505
  * Schmidt et al. (2009) Schmidt W., Federrath C., Hupp M., Kern S., Niemeyer J. C., 2009, A&A, 494, 127
  * Schober et al. (2012) Schober J., Schleicher D., Federrath C., Klessen R., Banerjee R., 2012, Phys. Rev. E, 85, 026303
  * Schober et al. (2015) Schober J., Schleicher D. R. G., Federrath C., Bovino S., Klessen R. S., 2015, Phys. Rev. E, 92, 023010
  * Seta & Federrath (2020) Seta A., Federrath C., 2020, MNRAS, 499, 2076
  * Seta & Federrath (2021) Seta A., Federrath C., 2021, Physical Review Fluids, 6, 103701
  * Seta et al. (2015) Seta A., Bhat P., Subramanian K., 2015, Journal of Plasma Physics, 81, 395810503
  * Seta et al. (2020) Seta A., Bushby P. J., Shukurov A., Wood T. S., 2020, Physical Review Fluids, 5, 043702
  * Shang et al. (2020) Shang J., Wu T., Wang H., Xu W., Ye C., Hu R., Tao J., He X., 2020, IEEE Access, 8, 40909
  * Sharma & Kumar (2019) Sharma B., Kumar R., 2019, Phys. Rev. E, 100, 013309
  * Sharma & Kumar (2023) Sharma B., Kumar R., 2023, arXiv e-prints, p. arXiv:2303.08400
  * Skalidis et al. (2021) Skalidis R., Sternberg J., Beattie J. R., Pavlidou V., Tassis K., 2021, A&A, 656, A118
  * Squire & Hopkins (2017) Squire J., Hopkins P. F., 2017, MNRAS, 471, 3753
  * Squire et al. (2023) Squire J., Kunz M. W., Arzamasskiy L., Johnston Z., Quataert E., Schekochihin A. A., 2023, Journal of Plasma Physics, 89, 905890417
  * St-Onge et al. (2020) St-Onge D. A., Kunz M. W., Squire J., Schekochihin A. A., 2020, Journal of Plasma Physics, 86, 905860503
  * Steinwandel et al. (2023) Steinwandel U. P., Dolag K., Böss L., Marin-Gilabert T., 2023, arXiv e-prints, p. arXiv:2306.04692
  * Stokes (1850) Stokes G. G., 1850, On the Effect of the Internal Friction of Fluids on the Motion of Pendulums. Cambridge University Press, p. 1–10, doi:10.1017/CBO9780511702266.002
  * Stroustrup (2013) Stroustrup B., 2013, The C++ Programming Language, 4th edn. Addison-Wesley Professional
  * Sur & Subramanian (2023) Sur S., Subramanian K., 2023, MNRAS,
  * Tisza (1942) Tisza L., 1942, Phys. Rev., 61, 531
  * Tobias (2021) Tobias S., 2021, Journal of Fluid Mechanics, 912, P1
  * Touber (2019) Touber E., 2019, Journal of Fluid Mechanics, 875, 974–1003
  * Uhlenbeck & Ornstein (1930) Uhlenbeck G. E., Ornstein L. S., 1930, Physical Review, 36, 823
  * Verma (2017) Verma M. K., 2017, Reports on Progress in Physics, 80, 087001
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Waagan et al. (2011) Waagan K., Federrath C., Klingenberg C., 2011, Journal of Computational Physics, 230, 3331
  * Wang et al. (2018a) Wang J., Wan M., Chen S., Xie C., Chen S., 2018a, Phys. Rev. E, 97, 043108
  * Wang et al. (2018b) Wang J., Wan M., Chen S., Chen S., 2018b, Journal of Fluid Mechanics, 841, 581–613
  * Weibel (1959) Weibel E. S., 1959, Phys. Rev. Lett., 2, 83
  * White et al. (2019) White T. G., et al., 2019, Nature Communications, 10, 1758
  * Whittingham et al. (2021) Whittingham J., Sparre M., Pfrommer C., Pakmor R., 2021, Monthly Notices of the Royal Astronomical Society, 506, 229
  * Whittingham et al. (2023) Whittingham J., Sparre M., Pfrommer C., Pakmor R., 2023, arXiv preprint arXiv:2301.13208
  * Xu & Lazarian (2016) Xu S., Lazarian A., 2016, ApJ, 833, 215
  * Yousef et al. (2007) Yousef T. A., Rincon F., Schekochihin A. A., 2007, Journal of Fluid Mechanics, 575, 111
  * Zel’dovich et al. (1984) Zel’dovich Y. B., Ruzmaikin A. A., Molchanov S. A., Sokolov D. D., 1984, Journal of Fluid Mechanics, 144, 1
  * van der Velden (2020) van der Velden E., 2020, The Journal of Open Source Software, 5, 2004
  * van der Walt et al. (2014) van der Walt S., et al., 2014, PeerJ, 2, e453

##  Appendix A Non-dissipative property of $\mathbb{A}$, the rate of rotation
tensor

While the rate of strain tensor Equation 8 contributes to shearing deformation
through isochoric transformations of fluid parcels, the rate of rotation
tensor Equation 10, instead represents rigid body rotations about a point. The
latter can be shown to neither grow nor dissipate energy. We show this by
constructing the $E_{\text{kin}}$ equation by left-multiplying
$\mathbf{v}\cdot$ to the strong form of the MHD momentum equation Equation 41
with zero external forcing ($\rho\mathrm{{\mn@boldsymbol{\mathit{f}}}}=0$),

$\displaystyle\frac{\partial E_{\text{kin}}}{\partial t}$
$\displaystyle-\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla\cdot\Bigg{[}\frac{1}{4\pi}\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}-\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}-p\mathbb{I}+\rho\bm{\sigma}\Bigg{]}=0,$
(59)

where $p$ is the total (thermal and magnetic) pressure. We direct our
attention solely to the term involving the fully decomposed viscous stress
tensor with $\mathbb{A}$ (indices explicitly shown for clarity),

$\bm{v}\cdot\nabla\cdot(\rho\mathbb{A})=v_{i}\mathbb{A}_{ij}\partial_{j}\rho+\rho
v_{i}\partial_{j}\mathbb{A}_{ij},$ (60)

where $\rho v_{i}\partial_{j}\mathbb{A}_{ij}=-\rho
v_{i}\partial_{j}\mathbb{A}_{ji}$ and
$v_{i}\mathbb{A}_{ij}\partial_{j}\rho=-v_{i}\mathbb{A}_{ji}\partial_{j}\rho$,
where we have applied the skew-symmetric property of $\mathbb{A}$, i.e.,
$\mathbb{A}_{ij}=-\mathbb{A}_{ji}$. Expanding the sum of either $\rho
v_{i}\partial_{j}\mathbb{A}_{ij}$ or $v_{i}\mathbb{A}_{ij}\partial_{j}\rho$,
one then finds that due to the skew-symmetric property of $\mathbb{A}_{ij}$,
each term has a corresponding term that cancels the other out, e.g.,

$\displaystyle\rho v_{i}\partial_{j}\mathbb{A}_{ij}$
$\displaystyle=\rho\big{(}\underbrace{v_{x}\partial_{y}\mathbb{A}_{xy}+v_{y}\partial_{x}\mathbb{A}_{yx}}_{0}+\ldots\big{)}=0.$
(61)

Therefore, $\mathbb{A}$ has no contribution of the evolution (growth or
dissipation) in the viscous stress term in the kinetic energy evolution, and
hence need not be considered in the isotropic-deviatoric decomposition of the
viscous stress tensor, Section 2.1.

##  Appendix B Expressing the $1^{\rm st}$ moment magnetic energy equation in
terms of tensor operations

In our main study, we present Equation 30, which is a compressible version of
$1^{\rm st}$ moment of the magnetic equation shown in Schekochihin et al.
(2004c). Seta et al. (2020) and Sur & Subramanian (2023) derive a similar
$1^{\rm st}$ moment magnetic energy equation, but we explicitly derive it in
terms of tensor products, unlike either of these studies, so that we can apply
our $\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ decomposition,
Equation 35. The non-relativistic induction equation with Ohmic dissipation
is,

$\displaystyle\frac{\partial\mathrm{{\mn@boldsymbol{\mathit{b}}}}}{\partial
t}=\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn<EMAIL_ADDRESS>(62)

Now taking the scalar product with $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$,

$\displaystyle\frac{\partial b^{2}}{\partial
t}=2\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})+2\eta\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{b}}}},$
(63)

where

$\displaystyle\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})=\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\mathrm{{\mn@boldsymbol{\mathit{b}}}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})-\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}},$
(64)

and,

$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})$
$\displaystyle=\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}:\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
(65)
$\displaystyle\quad-\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}:(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}-\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn<EMAIL_ADDRESS>(66)

The last term can be simplified,

$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}:\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}$
$\displaystyle=\frac{1}{2}\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{v}}}}b^{2})-\frac{1}{2}\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}:(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I},$
(67)

which means

$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})$
$\displaystyle=\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}:\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
(68)
$\displaystyle\quad-\frac{1}{2}\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}:(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}-\frac{1}{2}\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{v}}}}b^{2}).$
(69)

The dissipation term is simply,

$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{b}}}}=\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\cdot(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}).$
(70)

Hence the dynamical equation for $b^{2}$ is

$\displaystyle\frac{\partial b^{2}}{\partial t}$
$\displaystyle=2\left[\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}:\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\frac{1}{2}\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}:(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}-\frac{1}{2}\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{v}}}}b^{2})\right]$
$\displaystyle\quad+2\eta\left[\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\cdot(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})\right],$
(71)

Now we take the $1^{\rm st}$ spatial moments of the equation and divide by
$8\pi$ to put the time-derivative in our energy density unit system,

$\displaystyle\frac{\partial\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}{\partial t}$
$\displaystyle=\frac{2}{8\pi}\left\langle{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}:\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}:\frac{1}{2}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I})}\right\rangle_{\operatorname{\mathcal{V}}}$
$\displaystyle\quad+\frac{2\eta}{8\pi}\left\langle{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\cdot(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})}\right\rangle_{\operatorname{\mathcal{V}}},$
(72)

where
$\left\langle{\frac{1}{2}\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{v}}}}b^{2})}\right\rangle_{\operatorname{\mathcal{V}}}=0$
because there is no net magnetic flux through the triply periodic boundary.
This is in every way equivalent to Seta et al. (2020), but we retain the full
$\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$, instead of just
$\mathbb{S}$ (see Section 3.2 for more discussions on how
$\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ also contains tensorial
compression effects), as we mentioned in the main text. We do not retain the
advection term in the $1^{\rm st}$ moment, as in Sur & Subramanian (2023). Now
we simplify the dissipation term to get to the form expressed in Schekochihin
et al. (2004c). Firstly,

$\displaystyle\left\langle{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\cdot(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})}\right\rangle_{\operatorname{\mathcal{V}}}=\int_{\operatorname{\mathcal{V}}}\operatorname{d}\\!{\operatorname{\mathcal{V}}}\,\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\cdot(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}),$
(73)

which by parts $\int u\operatorname{d}\\!{v}=uv-\int v\operatorname{d}\\!{u}$,
where $u=\mathrm{{\mn@boldsymbol{\mathit{b}}}}$,
$\operatorname{d}\\!{u}=\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}$,
$\operatorname{d}\\!{v}=\nabla\cdot(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})$,
$v=\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}$ is

$\displaystyle\int_{\operatorname{\mathcal{V}}}\operatorname{d}\\!{\operatorname{\mathcal{V}}}\,\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\cdot(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})=$
$\displaystyle\cancelto{0}{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}\bigg{|}^{\infty}_{-\infty}}$
(74)
$\displaystyle\quad-\int_{\operatorname{\mathcal{V}}}\operatorname{d}\\!{\operatorname{\mathcal{V}}}\,(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}),$
(75)

where the first term in the integral goes to zero due to periodicity. Hence,

$\displaystyle\left\langle{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\nabla\cdot(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})}\right\rangle_{\operatorname{\mathcal{V}}}=-\left\langle{(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})}\right\rangle_{\operatorname{\mathcal{V}}}.$
(76)

Finally, we put this expression back into our evolution equation, and multiply
and divide by the integral energy density to get

$\displaystyle\frac{\partial\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}{\partial t}$
$\displaystyle=\frac{2\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}\Bigg{\langle}(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})$
$\displaystyle\quad\quad\quad-(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\left(\frac{1}{2}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}\right)\Bigg{\rangle}_{\operatorname{\mathcal{V}}}$
$\displaystyle\quad-\frac{2\eta\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}\Big{\langle}(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})\Big{\rangle}_{\operatorname{\mathcal{V}}},$
(77)

the $1^{\rm st}$ moment equations shown in Equation 30.

##  Appendix C Resolution study

Figure 15: The same as Figure 5 but instead of colouring by $\text{P}\nu$, we
fix $\text{P}\nu=1$ and change grid resolution of the simulation,
$N_{\text{grid}}$. We list the growth rate $\gamma_{1}$ and the final
saturation of $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
in Table 1, and find good convergence for both $\gamma_{1}$ and
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
saturation at $N_{\text{grid}}\approx 256$. Figure 16: The same as Figure 15
but for $\text{P}\nu=0.1$. At low $\text{P}\nu$, we similar convergence
properties as the $\text{P}\nu=1$ simulation, for both $\gamma_{1}$ and the
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
saturation.

In this section we provide a resolution study for our $\text{P}\nu=1$ and
$\text{P}\nu=0.1$ simulations, showing both the integral energy statistics, as
in Figure 5, in Figure 15 and Figure 16, and the Helmholtz decomposed spectra
in Figure 17 and Figure 18, for $\text{P}\nu=1$ and $\text{P}\nu=0.1$,
respectively. All statistics are coloured by the grid resolution, $N_{\rm
grid}^{3}=128^{3}$, $256^{3}$, $512^{3}$ and (only for $\text{P}\nu=1$)
$1024^{3}$. Beyond $256^{3}$, the integral energy statistics show strong
agreement, both in saturation and in growth rates (see Table 1), with very
little deviation based on the grid resolution. As we have pointed out in the
main text, this is because the energy integral statistics are always dominated
by the low-$k$ modes (Beattie et al., 2022a), so we expect them to converge
quickly, and at low grid resolutions. Similarly for the growth rates, because
we have only moderate $\text{Re}_{\text{shear}}=1000$ the growth rates, which
are sensitive to $\ell_{\nu_{\text{shear}}}$, also converge quickly.

Figure 17: The compressible ($E_{\text{kin},\parallel}(k)$; top) and
incompressible ($E_{\text{kin},\perp}(k)$; bottom) kinetic energy spectra for
the $\text{P}\nu=1.0$ ensemble of simulations (the same as the middle and
bottom panel in Figure 9), tabulated in Table 1. The inset in the top panel
shows the total kinetic energy $E_{\text{kin}}(k)$, and the inset in bottom
panel shows the ratio between compressible and incompressible spectral kinetic
energies, $E_{\text{kin},\parallel}(k)/E_{\text{kin},\perp}(k)$. We find good
convergence in $E_{\text{kin},\perp}(k)$, flattening out well before the
Nyquist frequency, indicating that $k_{\nu_{\text{shear}}}$ is resolved.
However, for $E_{\text{kin},\parallel}(k)$, the dissipation scale appears to
be unresolved, which results in an extended sub-shear-viscous range of scales
that grows with $N_{\text{grid}}$. Figure 18: The same as Figure 17 but for
the $\text{P}\nu=0.1$ ensemble of simulations, which is the $\text{P}\nu$
where the bulk viscosity and shear viscosity are approximately equal. With the
inclusion of bulk viscosity that matches the shear viscosity,
$E_{\text{kin},\parallel}(k)$ now has a resolved dissipation scale that is
(qualitatively) at a similar $k$ mode as the dissipation scale in
$E_{\text{kin},\perp}(k)$.

Now we compare the decomposed kinetic energy spectra in Figure 17 and Figure
18. Firstly, we can confidently say that the $\text{P}\nu=0.1$ simulation is
converged on all scales at $512^{3}$. That is, the physical length scales
associated with the bulk and shear viscosity are properly resolved in the
domain, since both $E_{\text{kin},\parallel}(k)$ and $E_{\text{kin},\perp}(k)$
trace each other almost perfectly through all the $k$ modes, and then truncate
before the Nyquist frequency. Furthermore, as we expect from our definition of
$\mathbb{B}$ (bulk viscosity operator; see Equation 35), the overall structure
of the two spectra are almost identical, i.e., $k_{\nu_{\text{shear}}}\approx
k_{\nu_{\text{bulk}}}$ (which should be exact for $\text{P}\nu=1/3$). However,
for the $\text{P}\nu=1$ simulation, $k_{\nu_{\text{shear}}}$ looks resolved
between $256^{3}$ and $512^{3}$ (as expected from our growth rates), but
$k_{\nu_{\text{bulk}}}$ does not, i.e., at $\text{P}\nu=1$, we do not have
sufficient resolution, even at $1024^{3}$, to resolve the entire sub-shear-
viscous range, which we see growing in size in the
$E_{\text{kin},\parallel}(k)/E_{\text{kin},\perp}(k)$ Figure 17 inset,
reaching values up to
$E_{\text{kin},\parallel}(k)/E_{\text{kin},\perp}(k)=10^{3}$ within this range
of scales. This would be a problem if in our study we focused on measuring
precise viscous dissipation scales, and viscous scale relations. However, we
focus on a mode-by-mode analysis, for the most part, and one can read off from
Figure 17, that even though the physical $k_{\nu_{\text{bulk}}}$ is not
resolved, we get relatively good convergence up to the $k\approx 100$ mode
before there is significant deviation between the spectra due to numerical
(compressible) viscosity. This is completely sufficient for our study. This
also highlights that for no explicit $\mathbb{B}$ viscosity
($\text{P}\nu=\infty$), we should not expect any convergence of the high-$k$
compressible modes and $E_{\text{kin},\parallel}(k)/E_{\text{kin},\perp}(k)$
should grow unbounded as $N_{\rm grid}^{3}$ increases, which exactly what one
expects for an ILES.

##  Appendix D Compressible and incompressible $1^{\rm st}$ moment energy
equations

We stated the Stokes limit for the compressible and incompressible $1^{\rm
st}$ moment kinetic equations in Section 6.3. Here we derive the full
equations and take the Stokes limit to demonstrate how we got to Equation 56
and Equation 55. We start with the dimensionless isothermal momentum equation,

$\displaystyle\frac{\partial\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}}{\partial
t}$
$\displaystyle-\nabla\cdot\left(\frac{1}{\mathcal{M}^{2}_{\text{A}}}\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}-\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\frac{1}{\mathcal{M}^{2}}\rho\mathbb{I}-\frac{1}{\mathcal{M}^{2}_{\text{A}}}b^{2}\mathbb{I}\right)$
$\displaystyle=\nabla\cdot\left(\frac{2\rho}{\text{Re}_{\text{shear}}}\mathbb{S}+\frac{\rho}{\text{Re}_{\text{bulk}}}\mathbb{B}\right)$
(78)

as in Section 2.2, and take the Stokes limit,
$\text{Re}_{\text{bulk}}=\text{Re}_{\text{shear}}\rightarrow 0$, which maps us
to the dynamical equations for the high-$k$, viscous modes in the sub-bulk-
viscous range (see Figure 1),

$\displaystyle\frac{\partial\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}}{\partial
t}=\nabla\cdot\left(\frac{2}{\text{Re}_{\text{shear}}}\mathbb{S}+\frac{1}{\text{Re}_{\text{bulk}}}\mathbb{B}\right)$
(79)

where we have assumed $\nabla\rho$ is small compared to
$\sim\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$. Now taking the dot
product with $\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}$ and taking
ensemble averages,

$\displaystyle\frac{1}{2}\frac{\partial\left\langle{\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{\partial
t}$
$\displaystyle=\frac{\partial\left\langle{E_{\text{kin},\perp}}\right\rangle_{\operatorname{\mathcal{V}}}}{\partial
t}$ (80)
$\displaystyle=\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla\cdot\left(\frac{2}{\text{Re}_{\text{shear}}}\mathbb{S}+\frac{1}{\text{Re}_{\text{bulk}}}\mathbb{B}\right)}\right\rangle_{\operatorname{\mathcal{V}}},$
(81)

utilising $\left\langle{v_{j}\partial_{t}\rho
v_{i}}\right\rangle_{\operatorname{\mathcal{V}}}=\partial_{t}\left\langle{\rho
v_{j}v_{i}}\right\rangle_{\operatorname{\mathcal{V}}}$ (Wang et al., 2018b).
Decomposing the remaining total velocity terms into
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}=\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$,
the first term is simply,

$\displaystyle\frac{2}{\text{Re}_{\text{shear}}}\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla\cdot\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}=$
$\displaystyle\frac{1}{\text{Re}_{\text{shear}}}\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}+\frac{4}{3}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}}\right\rangle_{\operatorname{\mathcal{V}}},$
(82)

and the second term is,

$\displaystyle\frac{1}{\text{Re}_{\text{bulk}}}\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla\cdot\mathbb{B}}\right\rangle_{\operatorname{\mathcal{V}}}=$
$\displaystyle\frac{1}{\text{Re}_{\text{bulk}}}\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}}\right\rangle_{\operatorname{\mathcal{V}}},$
(83)

hence

$\displaystyle\frac{\partial\left\langle{E_{\text{kin},\perp}}\right\rangle_{\operatorname{\mathcal{V}}}}{\partial
t}=$
$\displaystyle\left(\frac{4\text{Re}_{\text{shear}}^{-1}}{3}+\text{Re}_{\text{bulk}}^{-1}\right)\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}$
$\displaystyle+$
$\displaystyle\text{Re}_{\text{shear}}^{-1}\left\langle{\mathrm{{\mn<EMAIL_ADDRESS>(84)

Likewise for the $1^{\rm st}$ moment compressible kinetic energy equation,

$\displaystyle\frac{\partial\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}}{\partial
t}=$
$\displaystyle\left(\frac{4\text{Re}_{\text{shear}}^{-1}}{3}+\text{Re}_{\text{bulk}}^{-1}\right)\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}$
$\displaystyle+$
$\displaystyle\text{Re}_{\text{shear}}^{-1}\left\langle{\mathrm{{\mn<EMAIL_ADDRESS>(85)