AppendixAppendix References

# Taking control of compressible modes: bulk viscosity and the turbulent
dynamo

James R. Beattie${}^{\href https://orcid.org/0000-0001-9199-7771\,1,2,3,4}$,
Christoph Federrath${}^{\href https://orcid.org/0000-0002-0706-2306\,1,5}$,
Neco Kriel${}^{\href https://orcid.org/0000-0002-3558-3926\,1}$, Justin Kin
Jun Hew${}^{\href https://orcid.org/0000-0002-5238-6115\,1,6,7}$ & Amitava
Bhattacharjee${}^{\href https://orcid.org/0000-0001-6411-0178\,2,8}$
1Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia
2Department of Astrophysical Sciences, Peyton Hall, Princeton University,
Princeton, NJ 08544, USA
3Canadian Institute for Theoretical Astrophysics, University of Toronto, 60
St. George Street, Toronto, ON M5S 3H8, Canada
4Department of Astronomy and Astrophysics, University of California, Santa
Cruz, 1156 High Street, Santa Cruz, CA 96054
5Australian Research Council Centre of Excellence in All Sky Astrophysics
(ASTRO3D), Canberra, ACT 2611, Australia
6Space Plasma Power and Propulsion Laboratory, Department of Nuclear Physics
and Accelerator Applications,
Research School of Physics, Australian National University, ACT 2601,
Canberra, Australia
7Mathematical Sciences Institute, Australian National University, Canberra,
ACT, 2601, Australia
8Princeton Plasma Physics Laboratory, Princeton University, Princeton, NJ
08544, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Many polyatomic astrophysical plasmas are compressible and out of chemical and
thermal equilibrium, and yet, due to Stokes’ hypothesis, a means to carefully
control the decay of compressible modes in these systems has largely been
neglected. This is especially important for small-scale, turbulent dynamo
processes, which are known to be sensitive to the effects of compression. To
control the viscous properties of the compressible modes, we perform
supersonic, visco-resistive dynamo simulations with additional bulk viscosity
$\nu_{\text{bulk}}$, deriving a new $\nu_{\text{bulk}}$ Reynolds number
$\text{Re}_{\text{bulk}}$, and viscous Prandtl number
$\text{P}\nu\equiv\text{Re}_{\text{bulk}}/\text{Re}_{\text{shear}}$, where
$\text{Re}_{\text{shear}}$ is the shear viscosity Reynolds number. For
$10^{-3}\leq\text{P}\nu\leq\infty$, we explore a broad range of statistics
critical to the dynamo problem, including the integral and spectral energy
ratios, growth rates, and the magnetic $E_{\rm mag}(k)$ and kinetic
$E_{\text{kin}}(k)$ energy spectrum. We derive a general framework for
decomposing $E_{\rm mag}$ growth rates into incompressible and compressible
terms via orthogonal tensor decompositions of
$\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$, where
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ is the fluid velocity. We find that
compressible modes play a dual role, growing and decaying $E_{\rm mag}$, and
that field-line stretching is the main driver of growth, even in supersonic
dynamos. In the absence of $\nu_{\text{bulk}}$ ($\text{P}\nu\to\infty$),
compressible modes pile up on small-scales, creating an apparent spectral
bottleneck, which disappears for $\text{P}\nu\approx 1$. As $\text{P}\nu$
decreases, compressible modes are dissipated at increasingly larger scales, in
turn suppressing incompressible modes through a coupling between viscosity
operators. We emphasise the importance of further understanding the role of
$\nu_{\text{bulk}}$ in compressible astrophysical plasmas and direct numerical
simulations that include compressibility.

###### keywords:

MHD – turbulence – dynamo – magnetic fields – plasmas

††pubyear: 2023††pagerange: Taking control of compressible modes: bulk
viscosity and the turbulent dynamo–D

##  1 Introduction

Compressibility changes the nature of turbulence, from the energy cascade
(e.g, Kritsuk et al., 2007; Lithwick & Goldreich, 2001; Federrath, 2013; Wang
et al., 2018a; Ferrand et al., 2020; Federrath et al., 2021) to the magnetic
field fluctuations embedded within the turbulent medium (e.g., Mee &
Brandenburg, 2006; Federrath, 2016; White et al., 2019; Mandal et al., 2020;
Beattie & Federrath, 2020; Skalidis et al., 2021; Seta & Federrath, 2021; Bott
et al., 2021; Beattie et al., 2022b; Beattie et al., 2022a). Naturally, one
way of parameterising the compressibility of a medium is through the root-
mean-squared (rms) sonic Mach number
$\mathcal{M}=\left\langle{(v/c_{s})^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}$
where $v$ is the fluid velocity, $\operatorname{\mathcal{V}}$ is the system
volume and $c_{s}$ is the sound speed within the volume. Directly through the
Rankine-Hugoniot shock jump relations, the rms $\mathcal{M}$ tells us both the
amplitude of a characteristic shock in a supersonic $\mathcal{M}>1$ medium
(Lehmann et al., 2016; Park & Ryu, 2019), and by how much the shocks will
compress the gas structure (Landau & Lifshitz, 1959; Padoan & Nordlund, 2011;
Molina et al., 2012; Nolan et al., 2015; Federrath & Banerjee, 2015; Beattie
et al., 2020; Hew & Federrath, 2023; Davidovits et al., 2022), which directly
changes the statistics of all the fluid variables through flux freezing.
Regardless of whether a plasma is compressible or not, in the presence of
turbulent motions, magnetic fields will be grown and nourished through a
small-scale, turbulent dynamo (also called fluctuation dynamo, or just small-
scale dynamo (SSD); for reviews, Brandenburg & Subramanian, 2005; Rincon,
2019; Tobias, 2021).

Broadly speaking111Note we do not discuss the diffusion-free stage, which
happens before the kinematic dynamo, where the magnetic field is successively
folded to smaller-scales, until it reaches the resistive dissipation scale
(Schekochihin et al., 2002b, 2004c; St-Onge et al., 2020)., the turbulent
dynamo converts turbulent kinetic energy in the plasma into magnetic energy
through a three-stage process. First, while
$|\mathrm{{\mn@boldsymbol{\mathit{b}}}}|\ll|\mathrm{{\mn@boldsymbol{\mathit{v}}}}|$
– where $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$ is the magnetic field
permeating through the fluid – the induction equation becomes (approximately)
linear, leading to exponential growth in the volume-averaged magnetic energy,
$\left\langle{E_{\rm mag}}\right\rangle_{\operatorname{\mathcal{V}}}$,

$\displaystyle\frac{\operatorname{d}\\!{\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}}{\operatorname{d}\\!{t}}$
$\displaystyle=2\gamma_{1}\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}.$ (kinematic) (1) Considering
the standard incompressible phenomenology, once the Reynolds stress becomes
comparable to the magnetic tension
$\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})\sim\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})$,
the magnetic field exerts a back-reaction on the momentum equation (through
the Lorentz force) in turn making the system of equations nonlinear, and
suppressing the exponential growth until it becomes algebraic (specifically,
linear) in time $\displaystyle\frac{\operatorname{d}\\!{\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}}{\operatorname{d}\\!{t}}$
$\displaystyle=2\gamma_{2}.$ (nonlinear) (2) Finally, with the field growth
suppressed, the plasma resistivity balances with the growth terms and the
velocity field maintains $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ in a statistically stationary
state where $\displaystyle\frac{\operatorname{d}\\!{\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}}{\operatorname{d}\\!{t}}$
$\displaystyle=0,$ (saturated) (3)

where $\gamma_{1}$ and $\gamma_{2}$ are the growth rates of the respective
stages (e.g., Schober et al., 2015; Xu & Lazarian, 2016; Seta & Federrath,
2020; Galishnikova et al., 2022; Kriel et al., 2023). Compared to
incompressible, solenoidal turbulence, both the kinematic growth rate
$\gamma_{1}$ and $t\rightarrow\infty$ ratio of the volume-averaged energy
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
are reduced if the medium is highly compressible (Federrath et al., 2011;
Schober et al., 2015; Seta & Federrath, 2021; Sur & Subramanian, 2023; Kriel
et al., 2023) – that is, it appears more difficult to grow and maintain a
strong magnetic field in the presence of continuously driven shocks than in
their absence. In Section 3 we will discuss in more detail why this may be the
case.

Considering the effects of compression on the turbulent dynamo is important
for many astrophysical phenomena. For example, the interstellar medium in our
Galaxy is turbulent and trans-sonic in volume-average (Gaensler et al., 2011),
and turbulent and supersonic by mass-average (e.g., Klessen, 2011; Federrath
et al., 2016; Orkisz et al., 2017; Beattie et al., 2019). Hence, the coupling
and flow of energy between the turbulence and magnetic fields that are
embedded within the compressible interstellar medium are maintained with a
super-to-trans-sonic dynamo (e.g., Beck & Wielebinski, 2013; Xu & Lazarian,
2016; McKee et al., 2020; Gent et al., 2023; Beattie et al., 2023). This means
that all phenomena that either depend upon the magnetic field structure or
amplitude in the interstellar medium, e.g., cosmic rays trapped on magnetic
field lines, transported through the medium (Krumholz et al., 2020; Beattie et
al., 2022a; Kempski et al., 2023; Lemoine, 2023; Sampson et al., 2023), or
stars forming in magnetised H2 clouds (Federrath & Klessen, 2012; Krumholz &
Federrath, 2019; McKee et al., 2020; Mathew & Federrath, 2021), must have
initial magnetic fields that are seeded and maintained by a supersonic
turbulent dynamo (supersonic dynamo hereafter). Hence, all of these processes
require an understanding of not just a dynamo, but a supersonic dynamo.

The supersonic dynamo is ubiquitous across many scales in an active galaxy. It
has been shown that galaxy mergers excite these dynamos on kpc scales in the
interstellar medium, (Pakmor et al., 2014, 2017; Brzycki & ZuHone, 2019;
Whittingham et al., 2021, 2023; Pfrommer et al., 2022), but they can also be
found down at km scales in the plasma environment between merging neutrons
stars (Aguilera-Miret et al., 2023a; Chabanov et al., 2023) and ubiquitously
in quark gluon plasmas (e.g., Czajka et al., 2019). Of course, many of these
phenomena contain more than a single dynamo operating in tandem. For example
the magnetorotational instabilitiy dynamo, the $\alpha$ dynamo, shear current
effect and $\omega$ dynamo (see Brandenburg & Ntormousi, 2023, for recent
review) have been found to all operate at the same time at least for some
parameterisations of black hole accretion disks (Jacquemin-Ide et al., 2023).
The turbulent dynamo fits into these dynamo ecosystems by efficiently building
and maintaining the energy reservoir through the electromotive force,
$\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}}}\right\rangle_{\operatorname{\mathcal{V}}}$,
that in turn fuels the operation of the other dynamos (Rincon, 2019;
Brandenburg & Ntormousi, 2023).

Even though compressibility has been shown to have a strong effect on the
statistics of the turbulent dynamo, previous supersonic dynamo experiments
have solely focused on the shear viscosity, $\nu_{\text{shear}}$, and the
corresponding shear Reynolds number $\text{Re}_{\text{shear}}$ to control the
viscous dissipative properties of the plasma. Changing
$\text{Re}_{\text{shear}}$ directly controls the decay timescale for the
solenoidal modes e.g.,
$t_{\nu_{\text{shear}}}\sim(k_{\nu_{\text{shear}}}^{2}\nu_{\text{shear}})^{-1}$
(where $k_{\nu_{\text{shear}}}$ is the dissipation wavemode, with
corresponding length scale $\ell_{\nu_{\text{shear}}}$), and the growth rate
in the kinematic stage, $\gamma_{1}\sim
v_{\nu_{\text{shear}}}/\ell_{\nu_{\text{shear}}}\sim(\left\langle{v^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2})\text{Re}_{\text{shear}}^{1/2}$
where the growth of the magnetic field is set by the dynamical timescale on
$\ell_{\nu_{\text{shear}}}$ (Batchelor, 1950; Yousef et al., 2007;
Schekochihin et al., 2004c; Xu & Lazarian, 2016; Galishnikova et al., 2022).
However, the decay timescale (and corresponding length scale) for compressible
modes,
$t_{\nu_{\text{bulk}}}\sim(k_{\nu_{\text{bulk}}}^{2}\nu_{\text{bulk}})^{-1}$
(and $\ell_{\nu_{\text{bulk}}}$), have been left to decay on the same
timescale as the incompressible modes, where $\nu_{\text{bulk}}$ is the bulk
(volume or longitudinal) viscosity (e.g., Pan & Johnsen, 2017; Chen et al.,
2019; Sharma & Kumar, 2023), which we discuss in much more detail in Section
2. Because the dynamical timescale of compressible turbulent modes is known to
be fast, $t\sim\ell/v_{\ell}\sim\ell/\ell^{1/2}\sim\ell^{1/2}$ (assuming
Burgers 1948-type turbulence), compared to incompressible modes
$t\sim\ell/\ell^{1/3}\sim\ell^{2/3}$ (assuming Kolmogorov 1941-type
turbulence), this means that the diffusion timescale $t_{\nu_{\text{shear}}}$
may be insufficient to decay the compressible modes within a simulation
domain. Therefore, even in direct numerical simulations that explicitly set
$\text{Re}_{\text{shear}}$, the decay of the compressible modes may be
governed completely by numerical dissipation, which is controlled by the
numerical discretisation scheme for the fluxes in the simulation code, as well
as the number of resolution elements in the simulation domain (Grete et al.,
2023; Malvadi Shivakumar & Federrath, 2023).

Theoretically, under the Stokes (1850) hypothesis, $\nu_{\text{bulk}}=0$ is
justified, and exact for mono-atomic gases (Tisza, 1942), however, for di-
atomic gases, such as H2, or tri-atomic gases such as CO2, which are not in
thermal or chemical equilibrium, the magnitude of the bulk viscosity can vary
between tens and thousands of times larger than the shear viscosity (Pan &
Johnsen, 2017; Sharma & Kumar, 2019, 2023). This is certainly the case for the
plasma entrapped by merging neutron stars (Camelio et al., 2023), which
support dynamos excited by Kelvin Helmholtz instabilities (Chabanov et al.,
2023), where strong bulk viscosities
($\nu_{\text{bulk}}/\nu_{\text{shear}}\gtrsim 100$; Most et al. 2022a; Most et
al. 2022b) may play a role in the gravitational wave signatures of the merger
events(Chabanov & Rezzolla, 2023). Even beyond the theory of gases, the exact
role of bulk viscosity on the numerics of supersonic turbulent dynamo
simulations have remained largely unexplored.

In this study, we take the first steps toward describing the appropriate
framework for bringing bulk viscosity into the context of turbulent dynamos.
This includes defining new dimensionless plasma parameters, the bulk viscous
Reynolds number, and the viscous Prandtl number, and making a number of
important measurements that probe the role of this type of viscosity for the
dynamo. Our results are not only broadly applicable to general poly-atomic
plasmas, but also highly important for doing controlled studies of supersonic,
magnetised gases, where spatial and temporal properties of compressible modes
are controlled through an explicit viscosity operator. Our study is organised
in the following manner. In Section 2 we introduce the bulk viscosity through
orthogonal decompositions of the velocity gradient tensor, following a
somewhat pedagogical approach. In Section 3 we go on to provide the theory for
understanding compressibility in the context of key dynamo statistics – the
growth rate and the saturated energy ratio. In Section 4 we describe the
visco(both shear and bulk)-resistive magnetohydrodynamic (MHD) simulation
setup that we use to explore the role of bulk viscosity in turbulent dynamos.
In Section 5 we discuss the results of the saturation, dynamo growth rates,
and in general, the volume integral quantities of the plasma, including
decompositions into compressible and incompressible velocity modes. In Section
6 we explore the total, compressible and incompressible kinetic energy
spectra, highlighting the effects of bulk viscosity, not just on the
compressible spectrum, but also on the coupling between the compressible and
incompressible modes through the mixing of the modes and viscosity operators.
In Section 7 we focus on the magnetic spectra, the evolution of the integral
scales. In Section 8 we combine the decomposed kinetic and magnetic energy
spectra to explore how the dynamo saturates, on a mode-by-mode basis. Finally,
in Section 9 we summarise and conclude.

##  2 Fluid Viscosity & the velocity gradient tensor

###  2.1 Bulk and shear viscosity

The viscous momentum equation for a fluid is,

$\displaystyle\frac{\partial\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}}{\partial
t}-\nabla\cdot\mathbb{F}=\nabla\cdot\bm{\sigma},$ (4)

where $\mathbb{F}$ is a tensor of conserved momentum fluxes, and $\bm{\sigma}$
is viscous stress tensor. Before kinetic theory showed us the way to describe
$\bm{\sigma}$ through perturbations away from thermal equilibrium in the form
of the Maxwell-Boltzmann distribution function it was an empirical fact that
for Newtonian fluids $\bm{\sigma}$ was proportional to the velocity gradient
tensor $\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$, as,

$\displaystyle\|\bm{\sigma}\|\sim\|\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}\|,$
(5)

where $\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}=\partial_{i}v_{j}$
is the tensor product. By decomposing
$\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ into the rate of
expansion $\mathbb{B}$, the traceless, symmetric, rate of strain $\mathbb{S}$,
and the traceless, antisymmetric, rate of rotation $\mathbb{A}$, tensor
components, one can gain some insight into the structure of this type of
viscosity,

$\displaystyle\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
$\displaystyle=\mathbb{B}+\mathbb{S}+\mathbb{A},$ (6)
$\displaystyle\mathbb{B}$
$\displaystyle=\frac{1}{3}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I},$
(7) $\displaystyle\mathbb{S}$
$\displaystyle=\frac{1}{2}\left(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}+[\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}]^{T}\right)-\frac{1}{3}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I},$
(8) $\displaystyle\mathbb{A}$
$\displaystyle=\frac{1}{2}\left(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}-[\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}]^{T}\right).$
(9)

Because $\mathbb{S}$ is traceless,
$\text{tr}(\mathbb{S})=\mathbb{S}_{xx}+\mathbb{S}_{yy}+\mathbb{S}_{zz}=0$, it
describes volume-preserving (incompressible and isochoric) transformations of
fluid elements, and likewise for $\mathbb{A}$; where
$\text{tr}\left\\{\ldots\right\\}$ is the trace operator. The difference is
that all diagonal elements of $\mathbb{A}$ are identically zero, hence there
are no components of $\mathbb{A}$ orthogonal to the fluid element, whereas
$\text{tr}(\mathbb{S})$ sums to zero, hence the orthogonal forces exist, but
cancel out over the whole element. Furthermore, $\mathbb{S}$ is symmetric so
$\mathbb{S}=\mathbb{S}^{T}$, and $\mathbb{A}$ is anti-symmetric,
$\mathbb{A}=-\mathbb{A}^{T}$, which constrains the types of transformations
that are allowable for each tensor.

From this decomposition, it is clear that only $\mathbb{B}$ and $\mathbb{S}$,
contribute to the fluid viscosity, since

$\displaystyle\mathbb{A}=\frac{1}{2}(\partial_{j}v_{i}-\partial_{i}v_{j})=-\frac{1}{2}\epsilon_{ijk}\omega_{k}=\frac{1}{2}\begin{pmatrix}0&-\omega_{z}&\omega_{y}\\\
\omega_{z}&0&-\omega_{x}\\\ -\omega_{y}&\omega_{x}&0\end{pmatrix},$ (10)

where $\omega=\nabla\times\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ is the fluid
vorticity, and $\epsilon_{ijk}$ is the Levi-Civita tensor. This is exactly the
rate of rotation tensor, hence $\mathbb{A}$ cannot be a source of fluid
viscosity (see our derivation in Appendix A; note that for
$\sqrt{\mathbb{S}:\mathbb{S}}=0$, $\mathbb{A}$ describes the rigid body
rotation of a fluid element along its path – formerly, antisymmetric matrices
form a representation of the $\text{SO}(3)$ group, describing infinitesimal
rotations of the fluid elements; Choudhuri 1998). We may therefore write the
viscous stress tensor as,

$\displaystyle\bm{\sigma}=$ $\displaystyle
2\nu_{\text{shear}}\left[\frac{1}{2}\left(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}+[\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}]^{T}\right)-\frac{1}{3}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}\right]$
$\displaystyle\quad+\nu_{\text{bulk}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I},$
(11)

where $\nu_{\text{shear}}$ and $\nu_{\text{bulk}}$ are the coefficients of
bulk viscosity and the shear viscosity, respectively222Note that based on this
definition $\nu_{\text{bulk}}=\lambda+2\nu_{\text{shear}}/3$, where $\lambda$
is the second viscosity coefficient (e.g., Emanuel, 1990; Sharma & Kumar,
2019), but we prefer to parameterise our simulations in terms of
$\nu_{\text{bulk}}$, which is equivalent to changing $\lambda$ for fixed
$\nu_{\text{shear}}$.. Macroscopically, $\nu_{\text{bulk}}\mathbb{B}$ is the
normal viscous stress associated with non-thermal-equilibrium changes to the
volume of the fluid. To get insight into the microscopic description of
$\nu_{\text{bulk}}$, we turn to the two standard procedures for calculating it
(which is more challenging than the shear viscosity; Sharma & Kumar 2023):
Boltzmann equation-based relaxation time approximation and linear response
theory (Sharma & Kumar, 2019). Using the latter, for a dilute gas the bulk
viscosity can be calculated,

$\displaystyle\nu_{\text{bulk}}=\rho_{\rm
eq}c_{s}^{2}\frac{\gamma-1}{\gamma}\sum_{i}^{N}\frac{c_{v,i}}{c_{v}}\tau_{i},$
(12)

where $\rho_{\rm eq}$ is the equilibrium gas density, $\gamma$ is the ratio of
specific heats at equilibrium, $c_{v,i}$ is the heat capacity of the $i^{\rm
th}$ internal mode, $c_{v}$ is the total heat capacity and $\tau_{i}$ is the
relaxation time for each mode (where it assumed $\omega\ll\tau_{i}^{-1}$ and
$\omega$ is the sound wave frequency), summed over $N$ rotational and
vibrational degrees of freedom (Tisza, 1942). It is then easy to see that
$\nu_{\text{bulk}}$ is associated with the summed time lags $\tau_{i}$ that
each internal mode takes to relax to thermal equilibrium (indeed, each mode
may be decomposed into its own associated $\nu_{\text{bulk}}$; Shang et al.
2020). For mono-atomic gases, there are no such modes, so
$\nu_{\text{bulk}}=0$ (Stokes 1850 hypothesis), but for any other gas, this is
no longer the case.

As Sharma & Kumar (2023) suggest in their review of bulk viscosity in fluid
dynamics, the theory for bulk viscosity has largely been neglected in the
literature. Here we summarise the few studies that are relevant to bulk
viscosity in the context of fluid plasma turbulence. Chen et al. (2019)
studied forced homogeneous isotropic (and hydrodynamical) turbulence (HIT) for
moderate sonic Mach numbers $\mathcal{M}\leq 0.6$ and Taylor miscroscale
Reynolds numbers $\text{Re}_{\text{shear}}\approx 100$. Through the kinetic
energy spectrum, they showed that for large
$\nu_{\text{bulk}}/\nu_{\text{shear}}$ ratios the turbulence almost returned
to an incompressible state, i.e., the bulk viscosity acted primarily to decay
the compressible modes. Pan & Johnsen (2017) studied the effect of bulk
viscosity in decaying hydrodynamical turbulence, over a the same range of
$\nu_{\text{bulk}}/\nu_{\text{shear}}$ that we do in our study, finding that
high bulk viscosities enhanced the total decay rate of the kinetic energy by
both suppressing the compressible modes in magnitude, and through the
compressible mode-pressure and the compressible-incompressible mode
interactions. However, they found that the incompressible modes were
reasonably insensitive to bulk viscosity. Touber (2019) studied the eigen-
modes of the Navier-Stokes equation with the addition of bulk viscosity,
showing that in the Stokes (viscous) regime, significant bulk viscosity may be
able to disrupt the small-scale, incompressible structures (the building
blocks for enstrophy, $\omega^{2}$) in the turbulence, enhancing the
dissipation of the incompressible modes on small-scales. However, no studies
exist in the MHD turbulence regime, nor in the presence of the turbulent
dynamo, but if there is a strong coupling of the bulk viscosity at
$\ell_{\nu_{\text{shear}}}$ (undoubtedly transitioning into the Stokes
regime), then we should expect there to be an effect on the turbulent dynamo,
since $\ell_{\nu_{\text{shear}}}$ is responsible for the magnetic field
growth. To situate the bulk viscosity into a regular SSD formalism, we now
define the dimensionless numbers that we can use to parameterise the
additional bulk viscous fluxes in the momentum equation, and thus, our
simulations.

###  2.2 Dimensionless plasma numbers

####  2.2.1 Viscous plasma numbers

The bulk viscosity, as in Section 2.1, can be added simply as an additional
viscous flux in the momentum equation of the MHD fluid model, which we write
for an isothermal fluid,

$\displaystyle\frac{\partial\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}}{\partial
t}$
$\displaystyle-\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}}}}-\left(c_{s}^{2}\rho+\frac{b^{2}}{8\pi}\right)\mathbb{I}$
$\displaystyle+2\nu_{\text{shear}}\rho\Big{\\{}\frac{1}{2}\left(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}+[\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}]^{T}\right)-\frac{1}{3}\left(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}\right)\mathbb{I}\Big{\\}}$
$\displaystyle+\nu_{\text{bulk}}\rho\left(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}\right)\mathbb{I}\Bigg{]}=\rho\mathrm{{\mn@boldsymbol{\mathit{S}}}},$
(13)

where $\rho$ is the gas density, and
$\rho\mathrm{{\mn@boldsymbol{\mathit{S}}}}$ is an arbitrary momentum source.
One way of understanding the magnitude of each term in the momentum equation
is to make it dimensionless to see how each term scales with dimensionless
parameters. In light of this, we scale all length scales by $L$, all
velocities by $V$, all times by $T=L/V$, all magnetic fields by $B$, and all
densities by $\rho_{0}$, such that, e.g.,
$\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}=\mathrm{{\mn@boldsymbol{\mathit{v}}}}/V$
(dimensionless quantities are denoted with hat-notation). As a simple example,
the time-derivative term of the momentum becomes,

$\displaystyle\frac{\partial\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}}{\partial
t}\rightarrow\frac{\rho_{0}V}{T}\frac{\partial\hat{\rho}\mathrm{{\mn<EMAIL_ADDRESS>(14)

Let us non-dimensionalise the momentum equation, rearranging the magnetic
tension and pressure,

$\displaystyle\frac{\rho_{0}V^{2}}{L}\frac{\partial\hat{\rho}\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}}{\partial\hat{t}}$
$\displaystyle-\frac{1}{L}\widehat{\nabla}\cdot\Bigg{[}\frac{B^{2}}{4\pi}\left(\mathrm{{\mn@boldsymbol{\mathit{\hat{b}}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{b}}}}}-\frac{\hat{b}^{2}}{2}\mathbb{I}\right)-V^{2}\rho_{0}\hat{\rho}\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}$
$\displaystyle-
c_{s}^{2}\rho_{0}\hat{\rho}\mathbb{I}+\frac{2\nu_{\text{shear}}\rho_{0}\hat{\rho}V}{L}\left[\frac{1}{2}\left(\widehat{\nabla}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}+[\widehat{\nabla}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}]^{T}\right)\right.$
$\displaystyle\left.-\frac{1}{3}\left(\widehat{\nabla}\cdot\mathrm{{\mn<EMAIL_ADDRESS>(15)

We can now divide the whole equation by $\rho_{0}V^{2}/L$, appearing as the
coefficient for the time derivative of the momentum,

$\displaystyle\frac{\partial\hat{\rho}\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}}{\partial\hat{t}}$
$\displaystyle-\widehat{\nabla}\cdot\Bigg{[}\frac{B^{2}}{4\pi\rho_{0}V^{2}}\left(\mathrm{{\mn@boldsymbol{\mathit{\hat{b}}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{b}}}}}-\frac{\hat{b}^{2}}{2}\mathbb{I}\right)-\hat{\rho}\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}-\frac{c_{s}^{2}}{V^{2}}\hat{\rho}\mathbb{I}$
$\displaystyle+2\frac{\nu_{\text{shear}}\hat{\rho}}{VL}\left[\frac{1}{2}\left(\widehat{\nabla}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}+[\widehat{\nabla}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}]^{T}\right)-\frac{1}{3}\left(\widehat{\nabla}\cdot\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}\right)\mathbb{I}\right]$
<EMAIL_ADDRESS>(16)

Now we can make a number of observations. First, for the Lorentz force,
$V_{A}^{2}=B^{2}/(4\pi\rho_{0})$, and
$V_{A}^{2}/V^{2}=1/\mathcal{M}^{2}_{\text{A}}$ is the dimensionless Alfvén
Mach number. Likewise, for the thermal pressure,
$c_{s}^{2}/V^{2}=1/\mathcal{M}^{2}$ is the sonic Mach number. For the strain
rate dissipation model $\text{Re}_{\text{shear}}\sim VL/\nu_{\text{shear}}$
and for the bulk dissipation model $\text{Re}_{\text{bulk}}\sim
VL/\nu_{\text{bulk}}$. Hence, we can write the dimensionless momentum equation
as,

$\displaystyle\frac{\partial\hat{\rho}\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}}{\partial\hat{t}}$
$\displaystyle-\widehat{\nabla}\cdot\Bigg{[}\frac{1}{\mathcal{M}^{2}_{\text{A}}}\left(\mathrm{{\mn@boldsymbol{\mathit{\hat{b}}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{b}}}}}-\frac{\hat{b}^{2}}{2}\mathbb{I}\right)-\hat{\rho}\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}-\frac{\hat{\rho}}{\mathcal{M}^{2}}\mathbb{I}$
$\displaystyle+\frac{2\hat{\rho}}{\text{Re}_{\text{shear}}}\Bigg{(}\frac{1}{2}\Big{[}\widehat{\nabla}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}+\left(\widehat{\nabla}\otimes\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}\right)^{T}\Big{]}-\frac{1}{3}\left(\widehat{\nabla}\cdot\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}\right)\mathbb{I}\Bigg{)}$
$\displaystyle+\frac{\hat{\rho}}{\text{Re}_{\text{bulk}}}\left(\widehat{\nabla}\cdot\hat{\mathit{v}}\right)\mathbb{I}\Bigg{]}=\hat{\rho}\mathrm{{\mn@boldsymbol{\mathit{\hat{S}}}}},$
(17)

which demonstrates how each of the terms in the momentum equation are related
to each of the dimensionless constants333Previous studies have adopted
$\text{Re}_{\text{shear}}=2\pi\left\langle{v^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}/(k_{0}\nu_{\text{shear}})$,
where $k_{0}=\ell_{0}/(2\pi)$ is the turbulent driving scale in wavenumber
space. This leads to a boost in Reynolds number by a factor of $2\pi$, which
seems confounding, since all we have done is substituted $\ell_{0}$ for
$k_{0}$ in the definition. However, note that real-space and wavenumber space
statistics are linked by integrals, e.g.,
$\left\langle{v^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}=\int_{\forall
k}\operatorname{d}\\!{\mathrm{{\mn@boldsymbol{\mathit{k}}}}}\,|\mathrm{{\mn@boldsymbol{\mathit{v}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})\mathrm{{\mn@boldsymbol{\mathit{v}}}}^{*}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|,$
(18) Parsevel’s theorem, where
$|\mathrm{{\mn@boldsymbol{\mathit{v}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})\mathrm{{\mn@boldsymbol{\mathit{v}}}}^{*}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|$
is the power spectrum of $\mathrm{{\mn@boldsymbol{\mathit{v}}}}$, and hence
the velocity dispersion on $k_{0}$ is actually only a fraction of
$\left\langle{v^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}$,
$\left\langle{v^{2}}\right\rangle_{k_{0}+\operatorname{d}\\!{k}}^{1/2}=\int_{k_{0}}^{k_{0}+\operatorname{d}\\!{k}}\operatorname{d}\\!{\mathrm{{\mn@boldsymbol{\mathit{k}}}}}\,|\mathrm{{\mn@boldsymbol{\mathit{v}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})\mathrm{{\mn@boldsymbol{\mathit{v}}}}^{*}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|$.
This means, to use a wavenumber definition of $\text{Re}_{\text{shear}}$, one
has to also reduce
$\left\langle{v^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}$ such
that one is associating the correct velocities with the correct wavemodes.
This is unnecessarily complicated, so we always chose to use real-space
Reynolds number definitions, which come directly from the momentum equation,
as we have shown. we use to parameterise our flash simulations, including the
new definition of $\text{Re}_{\text{bulk}}$.

More precisely than simple scaling arguments, the shear viscosity is
controlled by changing the shear Reynolds number,

$\displaystyle\text{Re}_{\text{shear}}=\frac{|\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})|}{|2\nu_{\text{shear}}\nabla\cdot(\rho\mathbb{S})|}\sim\frac{\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}\ell_{0}}{\nu_{\text{shear}}},$
(19)

and the bulk viscosity via the bulk Reynolds number,

$\displaystyle\text{Re}_{\text{bulk}}=\frac{|\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})|}{|\nu_{\text{bulk}}\nabla\cdot(\rho\mathbb{B})|}\sim\frac{\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}\ell_{0}}{\nu_{\text{bulk}}},$
(20)

which in ratio (keeping $\nu_{\text{shear}}$ in the numerator, consistent with
the definition of the magnetic Prandtl number) defines our viscous Prandtl
number,

$\displaystyle\text{P}\nu=\frac{\nu_{\text{shear}}}{\nu_{\text{bulk}}}=\frac{\text{Re}_{\text{bulk}}}{\text{Re}_{\text{shear}}}\sim\frac{t_{\nu_{\text{bulk}}}}{t_{\nu_{\text{shear}}}}.$
(21)

With this definition, in the $\text{P}\nu\leq 1$ regime, bulk viscosity
dominates over the shear viscosity ($t_{\nu_{\text{bulk}}}\leq
t_{\nu_{\text{shear}}}$), and we therefore might expect that the compressible
modes decay on larger scales than the incompressible modes, and vice versa for
the $\text{P}\nu>1$ regime ($t_{\nu_{\text{bulk}}}>t_{\nu_{\text{shear}}}$).
We go through the natural implications of these two regimes on the energy
spectra in the following Section 3. But before that, we turn our attention to
the dimensionless parameters for the magnetic field.

Figure 1: Schematics for the energy spectrum in the kinematic stage of the
supersonic ($\mathcal{M}>1$), turbulent dynamo with
$\text{Pm}_{\text{shear}}=1$ (Equation 25), varying $\text{P}\nu$ (Equation
21), including the incompressible ($E_{\text{kin},\perp}$; solid, red),
compressible ($E_{\text{kin},\parallel}$; dotted, red) kinetic and magnetic
spectra ($E_{\rm mag}$; blue). We annotate the energy injection scale $k_{0}$,
peak magnetic energy scale $k_{\rm p}$, dissipation scale for the compressible
modes $k_{\nu_{\text{bulk}}}$, dissipation scale for the incompressible modes
$k_{\nu_{\text{shear}}}$ and dissipation scale for the magnetic field
$k_{\eta}$ on the $\log k$ axis. Note that for both the left and right panel,
we have adopted a key result from Kriel et al. (2023), who showed that $k_{\rm
p}<k_{\nu_{\text{shear}}}$ for the compressible turbulent dynamo, as opposed
to the incompressible case where $k_{\rm p}\sim k_{\eta}$ (Schekochihin et
al., 2002c, 2004c). Left: The $\text{P}\nu<1$ case, where
$k_{\nu_{\text{shear}}}>k_{\nu_{\text{bulk}}}$, decaying the compressible
modes faster than the incompressible modes and creating a sub-bulk-viscous
range of scales $k_{\nu_{\text{bulk}}}<k<k_{\nu_{\text{shear}}}$, where
$E_{\text{kin},\perp}(k)\gg E_{\text{kin},\parallel}(k)$. Right: The
$\text{P}\nu>1$ case, where $k_{\nu_{\text{bulk}}}>k_{\nu_{\text{shear}}}$,
creating a sub-shear-viscous range in kinetic energies
$k_{\nu_{\text{shear}}}<k<k_{\nu_{\text{bulk}}}$, where
$E_{\text{kin},\parallel}(k)\gg E_{\text{kin},\perp}(k)$, and for
$\text{Pm}_{\text{shear}}=1$, $E_{\text{kin},\parallel}(k)\gg E_{\rm mag}(k)$.

####  2.2.2 Resistive plasma numbers

Performing the same analysis as above for fluid the induction equation with
Ohmic resistivity gives,

$\displaystyle\frac{\partial\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}}{\partial\hat{t}}-\widehat{\nabla}\times(\hat{\mathrm{{\mn@boldsymbol{\mathit{v}}}}}\times\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}})=\frac{1}{\text{Rm}}\widehat{\nabla}^{2}\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}},$
(22)

where

$\displaystyle\text{Rm}=\frac{|\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})|}{|\eta\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{b}}}}|}\sim\frac{\left\langle{v^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}\ell_{0}}{\eta},$
(23)

is the magnetic Reynolds number for the Ohmic resistivity and $\eta$ is the
Ohmic resistivity coefficient. Unlike $\text{Re}_{\text{shear}}$ and
$\ell_{\nu_{\text{shear}}}$, the dimensionless number Rm alone is not
sufficient to define the resistive scale in the turbulence (Schekochihin et
al., 2004c; Galishnikova et al., 2022; Kriel et al., 2022; Brandenburg et al.,
2023; Kriel et al., 2023). Instead, by balancing the stretching rate at the
shear viscous scale, $v_{\ell_{\nu_{\text{shear}}}}/\ell_{\nu_{\text{shear}}}$
with the resistive dissipation rate, $\eta/\ell_{\eta}^{2}$,
$v_{\ell_{\nu_{\text{shear}}}}/\ell_{\nu_{\text{shear}}}\sim\eta/\ell_{\eta}^{2}$,
we can rearrange to find,

$\displaystyle\frac{\ell_{\nu_{\text{shear}}}}{\ell_{\eta}}\sim\text{Pm}_{\text{shear}}^{1/2},$
(24)

where

$\displaystyle\text{Pm}_{\text{shear}}=\frac{\nu_{\text{shear}}}{\eta}=\frac{\text{Rm}}{\text{Re}_{\text{shear}}},$
(25)

is the shear magnetic Prandtl number. Equation 24 is independent from the
effects of compressibility (Kriel et al., 2023). $\text{Pm}_{\text{shear}}$
parameterises the range of scales,
$\ell_{\nu_{\text{shear}}}>\ell>\ell_{\eta}$ (the classical sub-viscous range;
Schekochihin et al. 2002a; Schekochihin et al. 2004c). We can furthermore
define a bulk magnetic Prandtl number,

$\displaystyle\text{Pm}_{\text{bulk}}=\frac{\nu_{\text{bulk}}}{\eta}=\frac{\text{Rm}}{\text{Re}_{\text{bulk}}},$
(26)

and it follows that,

$\displaystyle\text{P}\nu=\frac{\text{Pm}_{\text{shear}}}{\text{Pm}_{\text{bulk}}},$
(27)

which relates all three of the Prandtl numbers that we define in this study.
$\text{Pm}_{\text{bulk}}$, like $\text{Pm}_{\text{shear}}$ may be related to a
scale separation
$\ell_{\nu_{\text{bulk}}}/\ell_{\eta}\sim\text{Pm}_{\text{bulk}}^{\alpha}$, if
the
$v_{\ell_{\nu_{\text{bulk}}}}\sim\ell_{\nu_{\text{bulk}}}/t_{\ell_{\nu_{\text{bulk}}}}$
modes are responsible changing the resistive properties of the plasma, but
because Kriel et al. (2023) found that Equation 24 is universal for high and
low $\mathcal{M}$, this is unlikely. We will constrain ourselves to a subset
of the total parameter range in this study, keeping $\text{Re}_{\text{shear}}$
and $\text{Pm}_{\text{shear}}$ fixed, with $\text{Pm}_{\text{shear}}=1$, but
changing $\text{P}\nu$, which means $\text{Pm}_{\text{bulk}}$ is solely
determined by $\text{P}\nu$, $\text{Pm}_{\text{bulk}}=\text{P}\nu^{-1}$. This
is because we are solely focusing on the role of bulk viscosity in the dynamo,
and keeping $\text{Pm}_{\text{shear}}$ and $\text{Re}_{\text{shear}}$ fixed
makes for a simple and informative experiment.

##  3 Compressibility in the turbulent dynamo

In this section we: (1) construct a general picture for what a compressible
small-scale dynamo with bulk viscosity looks like in $k$ space, using our
parameterisations, $\text{P}\nu$, $\text{Pm}_{\text{shear}}$ and
$\mathcal{M}$; and (2), highlight the theoretical aspects of the effects from
compressible modes on the growth rate and saturation of the turbulent dynamo.
We adopt an approach that focuses on the magnetic energy equation, as done
before in Schekochihin et al. (2004c), Seta & Federrath (2021) and Sur &
Subramanian (2023), with some important differences444Schekochihin et al.
(2002b) considered only the incompressible magnetic energy equation, and both
Seta & Federrath (2021) and Sur & Subramanian (2023) did not decompose the
velocity gradient tensor to completely separate the incompressible and
compressible tensors, $\mathbb{A}$, $\mathbb{S}$ and $\mathbb{B}$,
respectively, whilst trying to directly measure the effect of incompressible
(stretching) and compressible effects in the dynamo. They only consider the
effect of projecting the symmetric component of the velocity gradient tensor –
the rate of strain tensor – the first term in Equation 35, onto the magnetic
field, critically neglecting $\mathbb{B}$. Furthermore, Sur & Subramanian
(2023) do not make their stretching tensor traceless, just using
$(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}+[\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}]^{T})/2$,
which no longer makes $\mathbb{S}$ volume-preserving, and hence is
contaminated by isotropic compression. Neglecting the additional compression
within $\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ gives rise to
overestimating the impact of compressions by a factor of $3$ in Seta &
Federrath (2021) and $6$ in Sur & Subramanian (2023), which may significantly
modify the $1^{\rm st}$ moment interpretations and conclusions presented in
Sur & Subramanian (2023). Steinwandel et al. (2023) also neglect the
$\mathbb{B}$ contribution from within their
$\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}\otimes\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}:\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
term in the magnetic field energy equation for their magnetised intracluster
medium simulations, which impacts their decomposition in a similar, negligible
way as Seta & Federrath (2021).. Let us begin with schematics of the dynamo
energy spectrum.

###  3.1 Setting the scene for the compressible dynamo

We schematically represent the $\text{P}\nu<1$ and $\text{P}\nu>1$ regimes in
Figure 1, with $\text{P}\nu<1$ in the left panel, and $\text{P}\nu>1$ in the
right panel, respectively. Both panels correspond to
$\text{Pm}_{\text{shear}}=1\implies\ell_{\nu_{\text{shear}}}\sim\ell_{\eta}$,
and $\mathcal{M}>1$ in both panels. We further annotate the fundamental scales
for the dynamo problem: the integral or energy injection scale $k_{0}$, the
peak magnetic energy scale $k_{\rm p}$, following the results from supersonic
dynamo experiments in Kriel et al. (2023), the viscous dissipation scales,
$k_{\nu_{\text{shear}}}$ and $k_{\nu_{\text{bulk}}}$, and the magnetic energy
dissipation scale $k_{\eta}$. We also annotate a Kolmogorov (1941)-type energy
spectrum $\sim k^{-5/3}$ for the incompressible mode spectrum
$E_{\text{kin},\perp}(k)$, giving rise to
$\ell_{0}\sim\text{Re}^{3/4}\ell_{\nu_{\text{shear}}}$, and a Burgers
(1948)-type spectrum $\sim k^{-2}$ for the compressible mode spectrum
$E_{\text{kin},\parallel}(k)$, which gives
$\ell_{0}\sim\text{Re}^{2/3}\ell_{\nu_{\text{bulk}}}$. This is also the
acoustic wave spectrum for $\mathcal{M}<1$. For the magnetic field, we
annotate a Kazantsev (1968)-type spectrum $\sim k^{3/2}$. For $\text{P}\nu<1$
(left panel), where the compressible modes decay on shorter timescales than
the incompressible modes, a new sub-bulk-viscous range of scales

$\displaystyle\ell_{0}>\ell_{\nu_{\text{bulk}}}>\ell>\ell_{\nu_{\text{shear}}}$
(sub-bulk-viscous range) (28)

emerges, where $E_{\text{kin},\perp}(k)\gg E_{\text{kin},\parallel}(k)$. This
regime is only achievable with explicit control over $k_{\nu_{\text{bulk}}}$,
which we expect may lead to the small-scales becoming fundamentally
incompressible. Conversely, for $\text{P}\nu>1$ (right panel), a sub-shear-
viscous range of scales emerges in the kinetic energy,

$\displaystyle\ell_{0}>\ell_{\nu_{\text{shear}}}>\ell>\ell_{\nu_{\text{bulk}}}$
(sub-shear-viscous range) (29)

where compressible modes are dominant, $E_{\text{kin},\parallel}(k)\gg
E_{\text{kin},\perp}(k)$. For $\text{Pm}=1$, the compressible kinetic energy
on these scales will also be beyond the magnetic energy,
$E_{\text{kin},\parallel}(k)\gg E_{\rm mag}(k)$, making for an interesting
range of scales that is naively fundamentally hydrodynamical, and dominated by
compressions. If Figure 1 is representative of how bulk viscosity works, then
the $\text{P}\nu>1$ regime ought to be the regime that all $\mathcal{M}>1$
simulations are in without any explicit bulk viscosity. It is not obvious if
this is desirable, but what is desirable will clearly depend strongly on the
$\text{P}\nu$ regime that exists in the turbulent plasma of interest. We aim
to test the predictions outlined in this schematic in the remainder of the
study.

Figure 2: The three growth terms in Equation 35 visualised as orthogonal
operations on the three 2D fluid elements (grey) that are evolving from the
elements with solid lines to the elements with dashed lines. Key of the
properties of the tensor operations are listed below the fluid elements. Each
panel represents a geometrical projection between the blue (magnetic) vectors
and the red vectors, which corresponds to the traceless rate of strain tensor
$\mathbb{S}$ (Equation 8; leftmost panel), the traceless rate of rotation
tensor $\mathbb{A}$ (Equation 9; central panel) and the rate of expansion
tensor $\mathbb{B}$ (Equation 7; rightmost panel). The $\mathbb{S}$ and
$\mathbb{A}$ tensors enforce volume-preserving, incompressible transformations
of the fluid element, that is $\text{tr}(\mathbb{S})=\text{tr}(\mathbb{A})=0$
(orthogonal forces sum to zero). Note that we have represented only one
possible stretching mode for $\mathbb{S}$, and $\mathbb{S}$ also contains all
shearing modes. $\mathbb{A}$ only has components in the plane of the fluid
element, with all orthogonal components being exactly zero, and all off-
diagonal components being equal to the components of the fluid vorticity
(Equation 10). Because $\mathbb{A}=-\mathbb{A}^{T}$, $\mathbb{A}$ corresponds
to shear-free rotations of the fluid element. The final
$(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}$ tensor
operation only has components in the orthogonal directions to the fluid
element face. Each force must always point either away or into the fluid
faces, because all components share the sign of
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}$. Hence
$(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}$ is either a
isotropic contraction, $\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}<0$,
as shown in the schematic, or dilation
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}>0$, changing the volume of
the fluid element.

Moreover, we further note that in the literature, previous studies have used
the ratio $\nu_{\text{bulk}}/\nu_{\text{shear}}$ to parameterise the different
plasma regimes, (e.g., Chen et al., 2019), but we prefer our definition555Note
that because $\nu_{\text{bulk}}=(2/3)\nu_{\text{shear}}+\lambda$,
$\text{P}\nu=3\nu_{\text{shear}}/(2\nu_{\text{shear}}+3\lambda)$, hence when
$\text{P}\nu=1$, $\nu_{\text{shear}}/\lambda=3$, as expected from Equation 7.
What this means is that at $\text{P}\nu=1/3$ the dissipation timescale of the
bulk and shear viscosity is the same, which we will later find is consistent
with our simulation results. , which (1) explicitly defines the extra
viscosity through a bulk viscosity Reynolds number, and (2) puts the bulk
viscosity in a framework similar to the resistivity and magnetic Prandtl
number.

###  3.2 Turbulent dynamo growth

The $1^{\rm st}$ moment of the $E_{\rm mag}$ equation can be found by taking
the scalar product of the magnetic field and induction equation and then
averaging over space (similiar to Schekochihin et al., 2002b; Seta &
Federrath, 2021, but retaining the full velocity gradient tensor and the
effects of compressibility; see Appendix B for our derivation). Considering
only Ohmic dissipation as a means of magnetic diffusion, it is,

$\displaystyle\frac{\operatorname{d}\\!{\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}}{\operatorname{d}\\!{t}}$
$\displaystyle=2\gamma\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}-2\eta k_{\rm
u}^{2}\left\langle{E_{\rm mag}}\right\rangle_{\operatorname{\mathcal{V}}},$
(30) $\displaystyle\gamma$
$\displaystyle=\frac{\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):[(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})-\frac{1}{2}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}]}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}},$
(31) $\displaystyle k_{\rm u}^{2}$
$\displaystyle=\frac{\left\langle{(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}},$
(32)

noting that for the incompressible, kinematic case,

$\displaystyle\gamma
t_{0}\sim\left\langle{\hat{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}}\otimes\hat{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})}\right\rangle_{\operatorname{\mathcal{V}}}t_{0}\sim
t_{0}/t_{\nu_{\text{shear}}}\sim\text{Re}_{\text{shear}}^{1/2},$ (33)

since $\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}\sim\mathbb{S}\sim
t_{\nu_{\text{shear}}}^{-1}$, as expected and discussed in Section 1 for
Kolmogorov (1941)-type turbulence. We know from Section 2.1 that we can
decompose $\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ into tensor
components for the different modes, $\mathbb{S}$ and $\mathbb{A}$ for the
incompressible modes, and $\mathbb{B}$ for the compressible modes. Hence we
may write the $\gamma$ equation as

$\displaystyle\gamma$
$\displaystyle=\overbrace{\frac{\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\mathbb{S}+\mathbb{A})}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}}^{\text{incompressible}}+\underbrace{\frac{\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):[\mathbb{B}-\frac{1}{2}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}]}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}}_{\text{compressible}}.$
(34)

Splitting the incompressible term and using $\mathbb{B}$ we find that the
growth rate can be understood as the combination of fundamental fluid element
deformations,

$\displaystyle\gamma$
$\displaystyle=\overbrace{\frac{\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}}^{\text{stretching}}+\underbrace{\frac{\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{A}}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}}_{\text{winding}}-\overbrace{\frac{\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}}{6\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}}^{\text{compression}}.$
(35)

The first term in Equation 35 is the projection of the magnetic field onto the
velocity gradient field that stretches it. One may expect that this term must
dominate the growth of the magnetic field, since field-line-stretching is
well-known to be the key growth mechanism in the turbulent dynamo (Moffatt &
Saffman, 1964; Zel’dovich et al., 1984; Schekochihin et al., 2004c; Seta &
Federrath, 2021; Kriel et al., 2022; St-Onge et al., 2020). Next is the
winding term, which is the projection of the magnetic field onto the rate of
rotation tensor, Equation 10. Consider
$\mathbb{M}=\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}$,
and note that $\mathbb{M}=\mathbb{M}^{T}$ is symmetric
$(b_{i}b_{j}=b_{j}b_{i})$. Then it follows from the definition of
$\mathbb{A}$,

$\displaystyle\begin{split}\mathbb{M}:\mathbb{A}&=\text{tr}\left\\{\mathbb{M}\mathbb{A}^{T}\right\\}=-\text{tr}\left\\{\mathbb{M}\mathbb{A}\right\\},\\\
&=-\frac{1}{2}\text{tr}\left\\{\mathbb{M}(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})\right\\}+\frac{1}{2}\text{tr}\left\\{\mathbb{M}(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})\right\\}\\\
&=0,\end{split}$ (36)

which means we expect
$\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{A}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}\equiv
0$. More broadly, this means turbulent magnetic fields are prevented from
undergoing any kind of local winding to enhance or decay the field, quite
unlike large-scale toroidal fields, which can undergo linear growth via
winding (Aguilera-Miret et al., 2023b). The final term is the compression,
which is the projection of the magnetic field onto the dilatational components
of the fluid element. This term has a negative sign, with the total sign
contribution of the term solely determined by the divergence, because

$\displaystyle\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}=\text{tr}\left\\{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\right\\}\sim
b^{2}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}).$ (37)

Hence, for $\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}<0$, e.g., shocked
fluid, the compressions enhance the growth rate by a local factor
$(1/6)|\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}|$. However, for
rarefactions this term reduces the growth rate by the same factor. Therefore,
one can ask, on average, in a volume, is the fluid dilated
$\left\langle{\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}}\right\rangle_{\operatorname{\mathcal{V}}}>0$
or compressed
$\left\langle{\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}}\right\rangle_{\operatorname{\mathcal{V}}}<0$?
It is well known that for supersonic compressible turbulence, most of the
fluid in the volume is
$\left\langle{\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}}\right\rangle_{\operatorname{\mathcal{V}}}>0$,
with only a very small amount of the fluid contained with the
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}<0$ structures (Passot &
Vázquez-Semadeni 1998; Schmidt et al. 2008; Hopkins 2013; Squire & Hopkins
2017; Mocz & Burkhart 2019; Beattie et al. 2022c and also see the
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ probability density
function in Chen et al. 2019 and in Appendix A1 and A2 in Sur & Subramanian
2023). Hence as the compressions become stronger (for e.g., larger
$\mathcal{M}$), this term, on average, ought to act to reduce the growth rate
of the volume-averaged magnetic energy, which has been shown in numerical
simulations (Federrath et al., 2011), however, not from evolving the bulk
viscosity, but rather changing $\mathcal{M}$.

In the nonlinear induction stage Equation 2, $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}\sim E_{\text{kin}}(\ell_{\rm
eq})$, on energy equipartition scale $\ell_{\rm eq}$ (St-Onge et al., 2020;
Galishnikova et al., 2022; Beattie et al., 2023). Hence the growth term in
Equation 30, $2\gamma\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$, becomes

$\displaystyle 2\gamma t_{0}\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}\sim\frac{\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}E_{\text{kin}}^{1/2}(\ell_{\rm
eq})}{\ell_{\rm eq}}t_{0}\sim\frac{E_{\text{kin}}^{3/2}(\ell_{\rm
eq})}{\ell_{\rm eq}}t_{0},$ (38)

which is approximately constant, giving rise to linear growth as described in
Section 1.

###  3.3 Dynamo saturation

Let us explore Equation 30 in the saturated stage. Setting
$\partial_{t}\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}=0$ for statistically
stationary saturation,

$\displaystyle\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}-$
$\displaystyle\frac{1}{6}\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}=$
$\displaystyle\quad\quad\quad\quad\quad\quad\eta\left\langle{(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})}\right\rangle_{\operatorname{\mathcal{V}}},$
(39)

previously derived in Yousef et al. (2007) for incompressible turbulence. Very
simply put, when the growth terms balance with the dissipation, the dynamo
saturates. However, these may balance in at least three different ways: (1)
the growth terms may be suppressed such that they eventually find their way to
$\eta\left\langle{(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})}\right\rangle_{\operatorname{\mathcal{V}}}$,
or (2)
$\eta\left\langle{(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})}\right\rangle_{\operatorname{\mathcal{V}}}$
may grow, eventually balancing with the growth terms, or (3) there is a
mixture of both. For $\mathcal{M}\approx 0.2$ turbulence, Sur & Subramanian
(2023) showed that
$\eta\left\langle{(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})}\right\rangle_{\operatorname{\mathcal{V}}}$
reduced by a small amount as the dynamo saturates666Note that for
$\mathcal{M}\approx 10$ turbulence Seta & Federrath (2021) showed that the
Ohmic dissipation probability density function (the current; not simply the
$1^{\rm st}$ moment, that we are commenting upon here) reduced in the
saturation. This is also shown for the supersonic runs in Sur & Subramanian
(2023)., hinting that it is the $1^{\rm st}$ moment growth terms that need to
shrink to eventually reach a saturation. Hence, the incompressible growth
term,
$\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}$,
is suppressed (presumably by the magnetic tension
$\sim\nabla\cdot[\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}]$).
Schekochihin et al. (2004b) showed that this could be achieved through
anisotropising the velocity field along folded field structures in
$\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}$.
Furthermore, Seta & Federrath (2021) showed that in the saturation,
$|\mathrm{{\mn@boldsymbol{\mathit{u}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}}|\approx
0$ for both $\mathcal{M}=0.1$ and $\mathcal{M}=10$ dynamos, which suggests
that all tensor projection growth terms weaken in the saturation.

Moreover, as shown in Sur & Subramanian (2023) for the $\mathcal{M}>1$ dynamo,
the compressible growth term,
$(1/6)\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}$,
is suppressed by the magnetic field pressure gradients $\nabla b^{2}/(8\pi)$
that are preferentially parallel to $\nabla\rho$, opposing compression, since
shocked gas can generate large $\nabla b^{2}$, organising the magnetic field
into filamentary structures that are supported by magnetic pressure (Molina et
al., 2012; Beattie et al., 2022c; Kriel et al., 2023). This is relevant in
regions of high $\nabla b^{2}$ (corresponding to high mass), but in the
magnetic voids,
$(1/6)\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}$,
which are the more volume-filling state of the gas, these terms act to reduce
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$,
which ought to decrease final saturation of the dynamo. We will revisit these
volume-averaged terms, computing them through the temporal evolution of the
simulations in both the kinematic and saturated stage in Section 5.

Table 1: Main simulation parameters and derived quantities.

Sim. ID | $\text{Re}_{\text{shear}}$ | $\text{Re}_{\text{bulk}}$ | $\text{P}\nu$ | $\nu_{\text{shear}}t_{0}/\ell_{0}^{2}$ | $\nu_{\text{bulk}}t_{0}/\ell_{0}^{2}$ | $\mathcal{M}$ | $2\gamma_{1}t_{0}$ | $\left(\frac{E_{\rm mag}}{E_{\text{kin}}}\right)_{\text{sat}}$ | $(\operatorname{\mathcal{M}_{\text{A}}})_{\text{sat}}$ | $N_{\text{grid}}^{3}$
---|---|---|---|---|---|---|---|---|---|---
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11)
P$\nu\infty$ | 1000 | $\infty$ | $\infty$ | $1\times 10^{-3}$ | - | 0.96 $\pm$ 0.03 | 1.48 $\pm$ 0.03 | 0.17 $\pm$ 0.03 | 2.49 $\pm$ 0.29 | $512^{3}$
P$\nu$1_1024 | 1000 | 1000 | 1 | $1\times 10^{-3}$ | $1\times 10^{-3}$ | 0.95 $\pm$ 0.03 | 1.42 $\pm$ 0.01 | 0.15 $\pm$ 0.04 | 2.62 $\pm$ 0.40 | $1024^{3}$
P$\nu$1 | 1000 | 1000 | 1 | $1\times 10^{-3}$ | $1\times 10^{-3}$ | 0.95 $\pm$ 0.03 | 1.42 $\pm$ 0.02 | 0.14 $\pm$ 0.03 | 2.71 $\pm$ 0.30 | $512^{3}$
P$\nu$1_256 | 1000 | 1000 | 1 | $1\times 10^{-3}$ | $1\times 10^{-3}$ | 0.95 $\pm$ 0.03 | 1.36 $\pm$ 0.04 | 0.15 $\pm$ 0.02 | 2.61 $\pm$ 0.22 | $256^{3}$
P$\nu$1_128 | 1000 | 1000 | 1 | $1\times 10^{-3}$ | $1\times 10^{-3}$ | 0.94 $\pm$ 0.04 | 1.36 $\pm$ 0.03 | 0.14 $\pm$ 0.04 | 2.76 $\pm$ 0.38 | $128^{3}$
P$\nu$1e-1 | 1000 | 100 | 0.1 | $1\times 10^{-3}$ | $1\times 10^{-2}$ | 0.94 $\pm$ 0.04 | 1.47 $\pm$ 0.05 | 0.20 $\pm$ 0.03 | 2.26 $\pm$ 0.30 | $512^{3}$
P$\nu$1e-1_256 | 1000 | 100 | 0.1 | $1\times 10^{-3}$ | $1\times 10^{-2}$ | 0.92 $\pm$ 0.03 | 1.44 $\pm$ 0.04 | 0.18 $\pm$ 0.05 | 2.44 $\pm$ 0.33 | $256^{3}$
P$\nu$1e-1_128 | 1000 | 100 | 0.1 | $1\times 10^{-3}$ | $1\times 10^{-2}$ | 0.93 $\pm$ 0.03 | 1.35 $\pm$ 0.05 | 0.12 $\pm$ 0.04 | 2.88 $\pm$ 0.36 | $128^{3}$
P$\nu$1e-2 | 1000 | 10 | 0.01 | $1\times 10^{-3}$ | $1\times 10^{-1}$ | 0.90 $\pm$ 0.02 | 1.46 $\pm$ 0.03 | 0.19 $\pm$ 0.07 | 2.38 $\pm$ 0.39 | $256^{3}$
P$\nu$1e-3 | 1000 | 1 | 0.001 | $1\times 10^{-3}$ | $1\times 10^{0\,\,\,\,}$ | 0.93 $\pm$ 0.04 | 1.28 $\pm$ 0.02 | - | - | $256^{3}$

  * Notes. All simulations have magnetic Reynolds number $\text{Rm}=1000$, corresponding to Ohmic resistivity coefficient $\eta=1\times 10^{-3}\ell_{0}^{2}/t_{0}$, a $\text{Pm}_{\text{shear}}=1$ (Equation 25) and a $\text{Pm}_{\text{bulk}}=\text{P}\nu^{-1}$ (Equation 26), hence we list only the viscous material properties in the table. P$\nu$1e-3 was not run to saturation due to the computational cost (extremely viscous, very small diffusion integration time step over long integration times), so we do not report saturated statistics for that simulation. Column (1): the unique simulation ID. Column (2): the shear Reynolds number of the plasma. Column (3): the same as column (2), but for the bulk viscosity. Column (4): the viscous Prandtl number. Column (5): the coefficient of shear viscosity in units of turbulent correlation times on the driving scale of the turbulence. Column (6): the same as column (5) but for the coefficient of bulk viscosity. Column (7): the turbulent Mach number in the kinematic stage of the dynamo. Column (8): the growth rate of the kinematic stage in units of turbulent correlation times. Column (9): the ratio of the volume-averaged magnetic energy and volume-averaged kinetic energy in the saturated stage of the dynamo. Column (10): the same as column (9) but for the Alfvénic Mach number. Column (11): the number of grid points used in the numerical discretisation. We perform a convergence study in Appendix C for both the P$\nu$1e-1 and P$\nu$1 simulations.

##  4 Numerical simulations & Methods

###  4.1 Visco-resistive magnetohydrodynamic fluid model

We use a modified version of the finite volume MHD code flash (Fryxell et al.,
2000; Dubey et al., 2008), utilising a second-order conservative MUSCL-Hancock
5-wave approximate Riemann scheme (Bouchut et al., 2010; Waagan et al., 2011),
and $\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{b}}}}$ parabolic diffusion
flux cleaning (Marder, 1987) to solve the 3D, visco-resistive, isothermal,
compressible MHD equations with a stochastic acceleration field acting to
drive the turbulence,

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}})$ $\displaystyle=0,$
(40)
$\displaystyle\frac{\partial\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}}{\partial
t}-\nabla\cdot\left[\frac{1}{4\pi}\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}-\rho\mathrm{{\mn<EMAIL_ADDRESS>50.00008pt\ignorespaces$
$\displaystyle\left.-\left(c_{s}^{2}\rho+\frac{b^{2}}{8\pi}\right)\mathbb{I}+\rho\bm{\sigma}\right]$
$\displaystyle=\rho\mathrm{{\mn@boldsymbol{\mathit{f}}}},$ (41)
$\displaystyle\frac{\partial\mathrm{{\mn@boldsymbol{\mathit{b}}}}}{\partial
t}-\nabla\times\left[\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}}-\eta\mathrm{{\mn@boldsymbol{\mathit{j}}}}\right]$
$\displaystyle=0,$ (42)
$\displaystyle\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{b}}}}$
$\displaystyle=0,$ (43)

where the viscous stress tensor is

$\displaystyle\bm{\sigma}$ $\displaystyle=2\nu_{\rm shear}\mathbb{S}+\nu_{\rm
bulk}\mathbb{B},$ (44)

using the same definitions as in Section 2.1. We solve the equations on a
triply periodic domain of dimension $\operatorname{\mathcal{V}}\equiv L^{3}$,
with uniform grids of $128^{3}-1024^{3}$ cells (we show a detailed convergence
study in Appendix C), where
$\left\langle{\mathrm{{\mn@boldsymbol{\mathit{b}}}}(t)}\right\rangle_{\operatorname{\mathcal{V}}}=\left\langle{\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}(t)}\right\rangle_{\operatorname{\mathcal{V}}}=0$
(i.e., no net magnetic flux over $\operatorname{\mathcal{V}}$, in the zero-
momentum frame of $\operatorname{\mathcal{V}}$), $c_{s}$ is the sound speed,
$\mathrm{{\mn@boldsymbol{\mathit{j}}}}=\nabla\times\mathrm{{\mn@boldsymbol{\mathit{b}}}}/(4\pi)$
is the current density, and $\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ is the
stochastic turbulent acceleration source term that drives the turbulence.

###  4.2 Turbulent driving

The turbulent source term $\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ follows an
Uhlenbeck & Ornstein (1930) process with finite $e$-fold correlation time,
$t_{0}=\ell_{0}/\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}$
(Schmidt et al., 2009; Federrath et al., 2010). The field
$\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ is constructed in Fourier space with
energy isotropically injected on the peak scale
$|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|=2$ (equivalently,
$\ell_{0}=L/2$) and falls off to zero with a parabolic spectrum within
$1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$. On $\ell_{0}$ we
use $t_{0}$ and the energy injection rate to control the rms velocity, which
we set to
$\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}/c_{s}=\mathcal{M}=1.0$
in the kinematic stage of the dynamo, ensuring that we have compressible modes
for the bulk viscosity to control. Furthermore, we inject energy such that
$\mathrm{{\mn@boldsymbol{\mathit{f}}}}_{\perp}^{2}/\mathrm{{\mn@boldsymbol{\mathit{f}}}}_{\parallel}^{2}=2$,
corresponding to 2:1 energy in incompressible
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{f}}}}_{\perp}=0$ and compressible
$\nabla\times\mathrm{{\mn@boldsymbol{\mathit{f}}}}_{\parallel}=0$ mode
components of $\mathrm{{\mn@boldsymbol{\mathit{f}}}}$, ensuring that the
compressible modes are being replenished on large-scales. See Federrath et al.
(2010) and Federrath et al. (2022) for more details about the TurbGen driving
source. We integrate each of the experiments from $t/t_{0}=0$ to
$t/t_{0}=100$, writing the 3D field variables to disk every $t/t_{0}=0.1$ for
post-processing.

###  4.3 Setting Reynolds and Prandtl numbers

Following Section 2, we parameterise the visco-resistive properties of our
simulations with five numbers. Throughout the study we fix
$\text{Re}_{\text{shear}}=1000$, corresponding to setting
$\nu_{\text{shear}}=1\times 10^{-3}\ell_{0}^{2}/t_{0}$ in Equation 20. This
ensures that the incompressible modes are in a state of turbulence in our
simulations (Frisch, 1995). We then vary $\text{P}\nu$ between
$\text{P}\nu=\text{Re}_{\text{bulk}}=\infty$, where $\nu_{\text{bulk}}=0$ (the
dissipation is controlled purely by numerical means) to $\text{P}\nu=0.001$,
where the $\nu_{\text{bulk}}=10^{3}\nu_{\text{shear}}$ and
$\text{Re}_{\text{bulk}}=1$, which corresponds to $t_{\nu_{\text{bulk}}}\sim
t_{0}$ (i.e., the diffusion is as fast as the dynamical time on $\ell_{0}$;
see column 6 in Table 1 where we write the viscous coefficients in units of
$t_{0}$). For $\text{P}\nu>1$, the simulations develop a sub-shear-viscous
range, Equation 29, where the compressible modes ought to dominate at
high-$k$; and for $\text{P}\nu<1$, they develop a sub-bulk-viscous range,
Equation 28, where compressible modes are suppressed whilst incompressible
modes are left to dynamically evolve with the turbulence. The two ranges of
scales are shown schematically in Figure 1.

We control the resistivity by setting the magnetic Reynolds number, Rm
Equation 23. We choose $\eta=1\times 10^{-3}\ell_{0}^{2}/t_{0}$ corresponding
to $\text{Rm}=1000$. The extent of the (classical) sub-viscous range is set
through the shear magnetic Prandtl number, $\text{Pm}_{\text{shear}}$
(Equation 23; the regular Pm that has been studied in the small-scale dynamo;
Schekochihin et al. 2004c; Federrath et al. 2014; Galishnikova et al. 2022;
Kriel et al. 2022; Brandenburg et al. 2023; Kriel et al. 2023). Since
$\text{Re}_{\text{shear}}=\text{Rm}=1000$, $\text{Pm}_{\text{shear}}=1$ for
all of the simulations in our study, hence we decay the incompressible and
magnetic modes similarly, without the development of a sub-viscous range,
$\ell_{\eta}\leq\ell\leq\ell_{\nu_{\text{shear}}}$. We leave the study of the
bulk viscosity in the context of evolving $\text{Pm}_{\text{shear}}$ for
future studies. Based on Equation 26, the corresponding bulk magnetic Prandtl
number, $\text{Pm}_{\text{bulk}}$, is completely controlled by changing
$\text{P}\nu$, because $\text{Pm}_{\text{bulk}}=\text{P}\nu^{-1}$ for
$\text{Pm}_{\text{shear}}=1$. Hence, for the remainder of the study, we focus
only on the principle parameters that are evolved across the present
simulations, $\text{Re}_{\text{bulk}}$ and $\text{P}\nu$. All of our diffusion
coefficients, $\nu_{\text{shear}}$, $\nu_{\text{bulk}}$, $\eta$ are constant
in time and space.

Figure 3: Two-dimensional slices of the total ($v$; left), incompressible
($v_{\perp}$; middle) and compressible ($v_{\parallel}$; right) velocity modes
in the kinematic stage ($t/t_{0}=30$) of the turbulent dynamo. All velocities
are scaled by the rms, which determines the colour of the plots, with the
corresponding vector field slice of the modes overlaid. The top row
corresponds to the $\text{P}\nu=\infty$ simulation, where there is no explicit
bulk viscosity, and the bottom row is for $\text{P}\nu=0.001$, where the the
decay timescale of the bulk viscosity $t_{\nu_{\text{bulk}}}$ is comparable to
the dynamical timescale of the turbulence on the outer scale $t_{0}$. Both the
$v/\left\langle{v^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}$ and
$v_{\perp}/\left\langle{v_{\perp}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}$
panels are aesthetically similar between the simulations, but the
$v_{\parallel}/\langle
v_{\parallel}^{2}\rangle_{\operatorname{\mathcal{V}}}^{1/2}$ panels are vastly
different, with only low-$k$ $v_{\parallel}$ modes left in the plasma when the
bulk viscosity is strong. Figure 4: Two-dimensional slices in the kinematic
stage of the dynamo ($t/t_{0}=30$) for $E_{\rm mag}/E_{\text{kin}}$ (left
column), $\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ (middle column)
and $|\nabla\times\mathrm{{\mn@boldsymbol{\mathit{v}}}}|$. Each of the
quantities are normalised by the root-mean-squared value of the field to focus
on the structure. The least viscous simulations are in the top row, where
$\text{P}\nu=\infty$ and $\text{Re}_{\text{bulk}}=\infty$, and the dissipation
of the compressible modes is controlled by the numerical viscosity (Malvadi
Shivakumar & Federrath, 2023). The most bulk viscous simulations are shown in
the bottom row where $\text{P}\nu=0.001$ (Equation 21) and
$\text{Re}_{\text{bulk}}=1$ (Equation 20). Magnetic field lines are shown in
blue in the $E_{\rm mag}/E_{\text{kin}}$ panels. The
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ panel show the strongest
effect as the bulk viscous increases down the page. Thin shocks (yellow
structures) in the $\text{P}\nu=\infty$ simulation (top row, middle column),
become wider, diffusing away until they are completely absent in the
$\text{P}\nu=0.001$ simulation (bottom row, middle column). There is a lesser
effect on the vorticity, which looses some of the small-scale structure as
$\text{P}\nu$ gets smaller. A similar effect is shown for $E_{\rm
mag}/E_{\text{kin}}$ and the magnetic field structure that can be seen in the
streamlines.

###  4.4 Initial field conditions

Following, e.g., Kriel et al. (2022), we initialise the velocity field to
$|\mathrm{{\mn@boldsymbol{\mathit{v}}}}(x,y,z,t=0)|/c_{s}=0$, with units
$c_{s}=1$ and the density field $\rho(x,y,z,t=0)=\rho_{0}$, where the density
has units $\rho_{0}=1$. Likewise, $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$ has
units of $c_{s}\rho_{0}^{1/2}=1$, and is initialised as a Gaussian field in
$k$-space with a parabolic spectral power distribution matching the turbulent
forcing function. The total energy in the turbulent fluctuations of the field
is set such that
$\left\langle{E_{\text{mag},0}}\right\rangle_{\operatorname{\mathcal{V}}}/\rho_{0}c_{s}^{2}=10^{-10}=\beta^{-1}$,
where
$\left\langle{E_{\text{mag},0}}\right\rangle_{\operatorname{\mathcal{V}}}$ is
the initial volume integral energy and $\beta$ is the plasma beta. This
corresponds to an incredibly weak seed magnetic field compared to the constant
thermal energy, allowing for the simulations to spend a significant time in
the kinematic stage of the dynamo.

###  4.5 Helmholtz decomposition of velocity field

Any vector field can be written as a vector sum of curl-free and divergence-
free vectors – the Helmholtz decomposition of a vector field. Because we are
seeking to understand the role of bulk viscosity on the compressible modes in
the supersonic dynamo, we specifically decompose our
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ into orthogonal compressible
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$
($\nabla\times\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}=0$) and
incompressible $\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}$
($\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}=0$) modes, such
that,

$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{v}}}}=\mathrm{{\mn<EMAIL_ADDRESS>(45)

This is easily done in Fourier space, where

$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}(\bm{\ell})=\mathfrak{Re}\left\\{\int\operatorname{d}\\!{\mathrm{{\mn@boldsymbol{\mathit{k}}}}}\,[\mathrm{{\mn@boldsymbol{\mathit{k}}}}\cdot\tilde{\mathrm{{\mn@boldsymbol{\mathit{v}}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})]\mathrm{{\mn@boldsymbol{\mathit{k}}}}/k^{2}e^{2\pi
i\mathrm{{\mn@boldsymbol{\mathit{k}}}}\cdot\bm{\ell}}\right\\},$ (46)

is the real $\mathfrak{Re}$ component of the inverse Fourier transform of the
Fourier transformed velocity
$\tilde{\mathrm{{\mn@boldsymbol{\mathit{v}}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})$
projected along $\mathrm{{\mn@boldsymbol{\mathit{k}}}}$. Then
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}=\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$,
ensuring that Equation 45 is maintained. Throughout this study we utilise the
1D shell-integrated (isotropic) power spectra for the magnetic and kinetic
energy,

$\displaystyle E_{\rm mag}(k)$
$\displaystyle=\frac{1}{8\pi}\int_{\Omega_{k}}\operatorname{d}\\!{\Omega_{k}}\,4\pi
k^{2}|\tilde{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}(\mathrm{{\mn@boldsymbol{\mathit{\mathrm{{\mn@boldsymbol{\mathit{k}}}}}}}})|^{2},$
(47) $\displaystyle E_{\text{kin}}(k)$
$\displaystyle=\frac{\rho_{0}}{2}\int_{\Omega_{k}}\operatorname{d}\\!{\Omega_{k}}\,4\pi
k^{2}|\tilde{\mathrm{{\mn@boldsymbol{\mathit{v}}}}}(\mathrm{{\mn@boldsymbol{\mathit{\mathrm{{\mn@boldsymbol{\mathit{k}}}}}}}})|^{2},$
(48)

where $\operatorname{d}\\!{\Omega_{k}}$ is the shell angle at radius $k$ and
$\tilde{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}(\mathrm{{\mn@boldsymbol{\mathit{\mathrm{{\mn@boldsymbol{\mathit{k}}}}}}}})$
and
$\tilde{\mathrm{{\mn@boldsymbol{\mathit{v}}}}}(\mathrm{{\mn@boldsymbol{\mathit{\mathrm{{\mn@boldsymbol{\mathit{k}}}}}}}})$
are the Fourier transforms of the magnetic and velocity field, respectively.
We likewise create the $E_{\text{kin},\parallel}(k)$ and
$E_{\text{kin},\perp}(k)$ spectrum, based on our definitions in Equation 45,
for exploring the Helmholtz decomposed kinetic energies.

We show examples of two-dimensional slices of the Helmholtz decomposed
velocity fields in Figure 3 in our $\text{P}\nu=\infty$ (top row) and
$\text{P}\nu=0.001$ simulations (bottom row). The first column corresponds to
the total velocity field, the second column the incompressible
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}$ field, and the last column the
compressible $\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$ fields, all
which are coloured by the magnitude normalised by the rms. Both the
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ and
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$ fields qualitatively
remained unchanged by the strong bulk viscosity, however the
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$ field loses all of its
high-$k$ structure, and becomes dominated by low-$k$ modes that are being
replenished by the forcing.

We show how this translates into changing
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ (middle panel) and
$\nabla\times\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ (right panel) in Figure 4,
with $\text{P}\nu=\infty$ (top row), $\text{P}\nu=0.1$ (middle row), and
$\text{P}\nu=0.001$ (bottom row) simulations, along with the structure of
$E_{\rm mag}/E_{\text{kin}}$ (left column). Similar to Figure 3, the largest
effect that changing $\text{P}\nu$ has is shown in the
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ field, where strong, thin
compressions (yellow) that form in the $\text{P}\nu=\infty$ simulations are
diffusive and thick in the $\text{P}\nu=0.1$ simulations, until they are
completely dissipated away in the $\text{P}\nu=0.001$ simulation. From this
plot, it is also clear both
$\nabla\times\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ and $E_{\rm
mag}/E_{\text{kin}}$ are also losing high-$k$ mode structure as $\text{P}\nu$
decreases. To understand this process quantitatively, we turn now to the
integral energy statistics.

Figure 5: Top: the integral magnetic to kinetic energy ratio
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
as a function of time in units of turbulent correlation times $t/t_{0}$ on the
driving scale $\ell_{0}$, coloured by $\text{P}\nu$. The purple region
indicates where we define the kinematic stage of the dynamo, $10\leq
t/t_{0}\leq 30$, which we use for the entirety of our study. This is where
$\left\langle{E_{\rm mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ exhibits
exponential growth (Equation 1), which we fit an exponential model to shown
with the black-dashed lines. We tabulate growth rates from these fits
$2\gamma_{1}$ for each curve in Table 1. The grey region indicates the
saturated stage of the dynamo (Equation 3). Similarly, we tabulate the
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
ratios averaged over this region as well as the
$\operatorname{\mathcal{M}_{\text{A}}}$ in Table 1. We find that for fixed
$\text{Re}_{\text{shear}}$ and Rm, moderate $\text{P}\nu$
($\nu_{\text{bulk}}\sim\nu_{\text{shear}}$) does not have a strong influence
on the growth rates or final $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
saturation, however for low $\text{P}\nu$
($\nu_{\text{bulk}}>\nu_{\text{shear}}$), the viscosity begins to have a
strong, negative effect on both statistics. Bottom: The same as the top plot
but for the $\mathcal{M}$ evolution of the simulations, which moves between
stationary states for the kinematic and saturated dynamo stages as the
magnetic flux modifies the kinetic energy spectrum.

##  5 Integral energies and growth rates

We start the exploration for the role of the bulk viscosity on the turbulent
dynamo via the integral quantities of energy, growth rates and growth rate
terms
($\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}$
tensor projections) that we derived in Section 3.2.

###  5.1 Volume integrated energies

We define the volume integral energies as,

$\displaystyle\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}=\frac{\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{8\pi},\quad\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}=\frac{\left\langle{\rho_{0}v^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{2},$
(49)

which we will use throughout this entire study. We plot the integral energy
ratio $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
(top) and $\mathcal{M}$ (bottom) in Figure 5, colouring by $\text{P}\nu$.
These colours will be maintained throughout the entire study. In log-linear
space, we immediately identify the kinematic (exponential growth; highlighted
in dark violet) and saturated (stationarity; highlighted with grey) stages of
the dynamo in the $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
plot, which all experiments undertake (i.e., all simulations are above the
critical Rm for small-scale dynamo growth; Ruzmaikin & Sokolov 1981; Haugen et
al. 2004b; Schekochihin et al. 2004a; Federrath et al. 2014; Seta et al.
2015). In log-linear scale, the nonlinear stage (Equation 2) between the two
regimes is quite obfuscated. However, it is clearly observed in the
$\mathcal{M}(t/t_{0})$ plot, shown in the bottom panel of Figure 5. The
transition happens between $30\leq t/t_{0}\leq 70$ between kinematic and
saturated stages. This corresponds to the backreacting magnetic field
steepening the kinetic energy cascade as
$\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})\sim\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})+\nabla
b^{2}$ throughout the $k$ modes in the cascade (see Fig. B1 in Beattie et al.
2023 and theory in Schekochihin et al. 2002a).

We tabulate all of the $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
saturation statistics in Table 1, averaged over the saturated
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
and $\operatorname{\mathcal{M}_{\text{A}}}$ in the grey region, $t/t_{0}\geq
75$, all of which indicate a quite small effect throughout
$0.01\leq\text{P}\nu\leq\infty$ range, always maintaining $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}\approx
0.15$, corresponding to $\operatorname{\mathcal{M}_{\text{A}}}\approx 2.5$,
consistent within $1\sigma$ with $\text{P}\nu=\infty$ calculations in
Federrath et al. (2011). This indicates that the bulk viscosity does not have
a large effect on the low-$k$ modes, which dominate the integral energies
(Beattie et al., 2022a). This is suggestive that previous $\mathcal{M}$
dependencies upon small-scale dynamo saturation, e.g., Federrath et al.
(2011); Seta & Federrath (2021); Sur & Subramanian (2023), are being
controlled (more specifically the magnetic field being suppressed) by the
lowest $k$ modes. Performing the same calculation with instead the Helmholtz
decomposed integral kinetic energies,

$\displaystyle\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}=\frac{\left\langle{\rho_{0}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{2},\quad\left\langle{E_{\text{kin},\perp}}\right\rangle_{\operatorname{\mathcal{V}}}=\frac{\left\langle{\rho_{0}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{2},$
(50)

provides a clearer picture for how the bulk viscosity effects the integral
plasma quantities.

In Figure 6 we show the $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}$
(solid lines) and $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin},\perp}}\right\rangle_{\operatorname{\mathcal{V}}}$
(dashed lines) in the top panel; and
$\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
(solid lines) and
$\left\langle{E_{\text{kin},\perp}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
(dashed lines) ratios in the bottom panel, coloured in the same fashion as
Figure 5. We turn our attention to the bottom panel. For all $\text{P}\nu$, we
find the turbulence is dominated by the incompressible motions, even though we
are replenishing the source on compressible modes with the same amount of
energy as the incompressible modes in $\mathrm{{\mn@boldsymbol{\mathit{f}}}}$.
Federrath et al. (2010) showed this is the case previously for hydrodynamical
turbulence – i.e., the forcing modes and the momentum modes cannot be treated
analogously. For $\text{P}\nu=\infty$, we find that
$\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}\approx
10^{-1}$ (consistent with what has been measured previously in Federrath et
al. 2011), down to $\text{P}\nu=0.001$, where
$\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}\approx
10^{-3}$, decreasing somewhat linearly with $\text{P}\nu$,

$\displaystyle\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}\sim\text{P}\nu,$
(51)

for $\text{P}\nu<1$. This illustrates that as $\text{P}\nu$ shrinks there is a
strong suppression of $\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$
modes, even for moderate $\text{P}\nu$, that we were not able to see before in
Figure 5.

In the top panel of Figure 6 we show how the kinetic energy in these modes
balance with the magnetic energy. As we expect from the top panel of Figure 5,
the final saturation of the
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$ modes (dashed lines) is
insensitive to the bulk viscosity, with $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin},\perp}}\right\rangle_{\operatorname{\mathcal{V}}}\approx
0.1$ for all $\text{P}\nu$. However, the same is not true for the
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$ modes, which saturate to
$1\leq\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}\leq
10$ (equipartition and beyond). Hence, the turbulence is sub-to-trans-Alfvénic
with respect to $\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$ modes,
which means that, for example, shocked gas will not be able to significantly
change the morphology of the magnetic field due to the fact that
$|\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})+\nabla
b^{2}|\gg|\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel})|$.
This has been recently shown to be the case even in the kinematic stage of the
supersonic dynamo (Kriel et al., 2023), but we find it happens here in a more
extreme manner with decreasing $\text{P}\nu$, highlighting the importance of
treating the $\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$ and
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}$ modes separately in the
supersonic dynamo problem.

###  5.2 Volume integrated growth terms

We compute the growth rate $2\gamma_{1}$ in units of $t_{0}^{-1}$ by fitting a
log-linear model in the dark violet region of the $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
data, which we show in Figure 5. We show the fits with a grey dashed line.
This corresponds to the kinematic stage of the dynamo, which we conservatively
define with the range $10\leq t/t_{0}\leq 30$. We tabulate all of these values
in Table 1, along with the $\mathcal{M}$ averaged over the kinematic regime.

Similar to the $t\rightarrow\infty$ saturation, for moderate $\text{P}\nu$,
$\gamma_{1}$ is quite insensitive to the bulk viscosity. The values range
between $2\gamma_{1}t_{0}=1.40\pm 0.02-1.48\pm 0.03$, i.e., every $1.4t_{0}$
($e$-fold time), the field grows by $\approx
e^{n}\left\langle{E_{\text{mag},0}}\right\rangle_{\operatorname{\mathcal{V}}}$,
where $n$ is the number of $e$-fold times. This suggests that the high-$k$
compressible modes that are being decayed do not have a strong influence on
small-scale shear viscous velocities, which are responsible for growing the
magnetic field (Schekochihin et al., 2004c; Kriel et al., 2022; Galishnikova
et al., 2022).

Figure 6: Top: the same as the top panel in Figure 5, but separated into
incompressible (solid) and compressible (dashed) components of
$E_{\text{kin}}$. The final saturated state for the incompressible modes is
below equipartition
$\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}<1$, whereas it is
significantly above equipartition for the compressible modes
$\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}>1$. Bottom: the integral
compressible kinetic energy
$\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}$
ratios with respect to the total kinetic energy
$\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$. We
plot only
$\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
because
$\left\langle{E_{\text{kin},\perp}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}\approx
1$ for all simulations. As $\text{P}\nu$ decreases, the total energy contained
within the compressible component of the kinetic energy significantly
decreases from
$\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}\approx
10^{-1}$ to
$\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}\approx
10^{-3}$, roughly following
$\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}\sim\text{P}\nu$
for $\text{P}\nu<1$.

However, once the simulations become significantly bulk viscous, such that the
decay time from our input $\nu_{\text{bulk}}$ becomes faster than the
dynamical timescale on the turbulent outer-scale
$t_{\nu_{\text{bulk}}}\sim(\ell_{0}^{2}/\nu_{\text{bulk}})^{-1}\sim t_{0}$,
which happens at $\text{Re}_{\text{bulk}}\leq 10,\;\text{P}\nu\leq 0.01$ at
our $\mathcal{M}$, then $\gamma_{1}$ becomes significantly affected. It
directly follows that it is the low-$k$ compressible modes that affect the
dynamo process, since the decay timescale has to become comparable to
$\ell_{0}$ to have an impact. The fact that there is an impact at all suggests
that either $\nu_{\text{bulk}}$ plays a direct role in the exponential stage
of the dynamo processes, directly by suppressing the growth processes, or that
it plays an indirect role by coupling to viscous incompressible modes (as
hypothesised in Touber, 2019), which are the modes responsible for growing the
field in $\text{Pm}_{\text{shear}}\geq 1$ turbulent plasmas (e.g., Brandenburg
et al., 1995; Schekochihin et al., 2004c; Federrath et al., 2011; Kriel et
al., 2022, 2023). Curiously, this is the opposite of the effect that one
observes when increasing $\mathcal{M}$, or enhancing the injection of
compressibility in the forcing function (Federrath et al., 2011; Achikanath
Chirakkara et al., 2021; Seta & Federrath, 2021; Sur & Subramanian, 2023;
Kriel et al., 2023). We will show in Section 6.3 that we can explain this
phenomena through $k_{\nu_{\text{bulk}}}$-$k_{\nu_{\text{shear}}}$ coupling,
and that changing $\nu_{\text{bulk}}$ (modifying the high-$k$ modes) and
changing $\mathcal{M}$ (modifying the low-$k$ modes) is not equivalent, hence
there is no tension with these studies. Now we turn to the growth terms
derived in Section 3.

In Figure 7 we plot the temporal probability distribution functions of each of
the three volume-averaged growth terms in Equation 35, indicated by changing
the line style of the distribution function. In the top plot, we show the
distributions in the kinematic stage, and the bottom plot, the saturated
stage, maintaining the same time window as we did in the previous section and
shown in Figure 5. Furthermore, we show the complementary cumulative density
function (CDFs),

$\displaystyle\mathbb{P}(x>X)=1-\int_{-\infty}^{X}\operatorname{d}\\!{x}\;p(x),$
(52)

for the
$x=\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$
(top rows) and
$x=\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}/(6\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}})$
(bottom rows) for each simulation in the kinematic stage (left columns) and
saturated stage (right columns) in Figure 8.

Let us begin in the kinematic stage, focusing on the top panel of Figure 7.
Firstly, both simulations show that the
$\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{A}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$
distribution function (dotted line) is identically zero, as we predicted from
Equation 36. This means that the tensors are on average orthogonal to one
another. Because the contribution from the winding term is therefore zero, we
ignore it for the remainder of the study and focus on the other two tensor
projections.

Figure 7: The temporal probability distribution functions of the stretching,
winding and compression terms in the magnetic growth, Equation 35 (shown with
different line styles), in the kinematic stage of the dynamo (top), and the
saturated stage (bottom), coloured by $\text{P}\nu$, showing only the two
$\text{P}\nu$ extremes ($\text{P}\nu=0.001$ and $\text{P}\nu=\infty$). Each
distribution function is constructed through kernel density estimation using a
Gaussian function basis. For both stages, the winding term is zero (as
predicted by Equation 36, i.e., the rate of rotation Equation 10, an
antisymmetric tensor, is everywhere orthogonal to the magnetic field, a
symmetric tensor). The growth rate is therefore an interplay between the
compression term
($\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}/(6\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}})$;
dashed line style; which is suppressed for large bulk viscosities) and the
stretching term
($\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$;
solid line style). In both stages (kinematic and saturated) the stretching
term is on average positive (net growth) and dominant over the compression.
The compression approximately opposes and enhances growth in equal amounts,
with a small tendency towards opposing growth, which is shown in more detail
for all simulations in Figure 8. Figure 8: The complementary cumulative
density functions Equation 52 for the stretching growth term
($\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$;
top panels) and the compression growth term
($\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}/(6\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}})$;
bottom panels), respectively, throughout the kinematic stage (left column) and
saturated state (right column) of the dynamo, binned by the $t/t_{0}$ ranges
in Figure 5. We use vertical, dashed, black lines to separate the values into
growth and decay states, as we did in Figure 7, as well as adding horizontal
lines to highlight the $16^{\rm th}$, $50^{\rm th}$ and $84^{\rm th}$
percentiles of a standard Gaussian distribution function. We find that the
stretching terms dominate the growth states, growing the field for all
quantiles above the $16^{\rm th}$ percentile, while the compression terms are
slightly more predisposed to decay states. Both growth terms contract at the
onset of saturation (become closer to Heaviside functions), and the stretching
term shifts towards $X=0$ by a factor of $\approx 2$, recovering the result
that growth terms are suppressed during the onset of saturation (e.g.,
Schekochihin et al., 2004b; Seta & Federrath, 2021; Sur & Subramanian, 2023).

The next largest contribution to the growth rate is
$\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}/(6\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}})$,
which is almost symmetric about zero (but very slightly negatively skewed; all
CDFs cross zero at $\mathbb{P}<0.5$ in the bottom panel of Figure 8). This
means, on average, in the kinematic stage, compressions make small, negative
contributions to the growth rate at our $\mathcal{M}$777A natural hypothesis
to explain the growth rates and $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
saturation in Federrath et al. (2011) and more recently Sur & Subramanian
(2023) and Kriel et al. (2023), at least from a mathematical perspective, is
that changing $\mathcal{M}$ shifts
$\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}/(6\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}})$
towards more negative values due to the increasing skewness,
$\left\langle{(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})^{3}}\right\rangle_{\operatorname{\mathcal{V}}}$,
as the strongest, positive divergence structures becomes less and less space-
filling, further inhibiting the growth and saturation, on average. In our
study, we focus on a single $\mathcal{M}$ and the effects of bulk viscosity,
so we leave testing this hypothesis for future studies that utilise our
decomposition.. However, temporal fluctuations, which are approximately
Gaussian, ensure that the compressions make contributions to both growing the
$\left\langle{E_{\rm mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ and
decaying $\left\langle{E_{\rm mag}}\right\rangle_{\operatorname{\mathcal{V}}}$
in approximately equal amounts. This highlights the dual role that
compressible modes have in the dynamo problem – almost equally decaying and
growing the field at $\mathcal{M}\approx 1$. As $\text{P}\nu$ gets smaller
(and $\text{Re}_{\text{bulk}}$ gets smaller), the fluctuations of the
compressions reduce (the CDFs become closer to Heaviside functions in Figure
8), but the time-average is maintained close to zero. This describes a process
similar to what Chen et al. (2019) found for the effects of bulk viscosity on
the spatial $\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ distribution
function – a somewhat symmetric contraction of the
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ distribution function that
maintains the $1^{\rm st}$ moment. Hence, as the bulk viscosity increases, any
fluctuations in the growth rate caused by either strong compression or
rarefaction events diminish, but the average contribution of the compressions
is maintained, hence the relatively weak dependence upon the growth rate for
moderate $\text{Re}_{\text{bulk}}$. This is a key result from this analysis.

To understand why at very low $\text{Re}_{\text{bulk}}$ we find that
$\gamma_{1}$ changes quite dramatically, we turn our attention to the growth
rate associated with stretching
$\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$.
Because the distribution function (solid line in Figure 7) is on average
positive, the interaction between the magnetic field and the $\mathbb{S}$
structures on average grow the field exponentially from Equation 30. As is
evident in Figure 8, its not only the average, but as low as the $16^{\rm th}$
percentile of
$\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$
that is able to grow the magnetic field. This is completely consistent with
the key idea behind kinematic turbulent dynamo theories – stretching plays the
dominant growth role in the kinematic stage of the dynamo (e.g., Moffatt &
Saffman, 1964; Haugen et al., 2004a; Zel’dovich et al., 1984; Brandenburg et
al., 1995), which we show is even the case for the supersonic dynamo. This is
why Kriel et al. (2023) was able to show numerically that the
$k_{\eta}/k_{\nu_{\text{shear}}}\sim\text{Pm}_{\text{shear}}^{1/2}$ relation
holds even in $\mathcal{M}>1$ fluids.

Finally, we turn to the bottom panel of Figure 7 and the right column in
Figure 8 to explore the differences between the kinematic and saturated stages
of the dynamo under these tensorial decompositions. Compared to the kinematic
stage, the
$\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$
terms both shift on average, with the $50^{\rm th}$ percentiles moving towards
zero, reducing by $\approx 2$ across all $\text{P}\nu$, and contract, best
shown in the bottom panel of Figure 7. The
$\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}}\right\rangle_{\operatorname{\mathcal{V}}}/(6\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}})$
term similarly contracts in the saturated stage. Seta & Federrath (2021) found
a similar trend, but their $\mathcal{M}$ was significantly higher than ours
and the difference between the saturated and kinematic stage statistics were
not as significant (we find a reduction of roughly a factor of two in the
$1^{\rm st}$ moments for
$\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):\mathbb{S}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{b^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$).
In their decomposition, at $\mathcal{M}=3$, Sur & Subramanian (2023) also find
a factor of $\approx 2$ reduction in the $1^{\rm st}$ moments between the
kinematic and the saturated stages. This highlights that our decomposition
retains the well-known result that the growth terms are suppressed in the
dynamo saturation.

Figure 9: The shell-integrated kinetic energy spectra averaged over the
kinematic stage of the dynamo. Top: The total kinetic energy spectrum,
coloured by $\text{P}\nu$, averaged over the kinematic regime (see Figure 5),
with inset showing flattening of the spectra under the Burgers (1948) scaling,
$E_{\text{kin}}(k)\sim k^{-2}$. The bulk viscosity acts to suppress the
high-$k$ mode structure in the total energy, suggestive that compressible
modes occupy significant energies on the smallest scales, contributing to a
bottleneck effect. Middle: The compressible (longitudinal) kinetic energy
spectrum, demonstrating that for $\text{P}\nu>1$, the dissipation in the
compressible modes is not resolved or controlled. Bottom: The incompressible
(transverse) kinetic energy spectrum with the ratio of the
$E_{\text{kin},\parallel}(k)$ and $E_{\text{kin},\perp}(k)$ in the inset. For
fixed $\text{Re}_{\text{shear}}$, reducing $\text{P}\nu$ suppresses the
high-$k$ compressible modes, acting more strongly on the compressible spectrum
for lower $\text{P}\nu$. The incompressible spectrum is also influenced, with
much less of an effect, with coupling starting deep in the viscous scales, but
slowly moving to lower $k$ modes as the decay timescale,
$t_{\nu_{\text{bulk}}}$, becomes shorter.

##  6 Kinetic Energy spectra

###  6.1 Total spectra

We plot the total kinetic energy spectra in the top panel of Figure 9,
coloured in the same fashion as the previous plots. At low-$k$, all
simulations have similar spectral energy contained on those scales, which
flattens under a $E_{\text{kin}}(k)\sim k^{-2}$ Burgers (1948) scaling as
shown in the inset panel, as expected for supersonic turbulence dominated by
singularities (Federrath, 2013; Verma, 2017; Federrath et al., 2021; Kriel et
al., 2023). However, as expected for a viscosity-like process (a process that
is scale-separated from the large-scales and starts by influencing the small-
scale dynamics), there is a significant difference between simulations in the
high-$k$ modes. Because $2^{\rm nd}$ moments have integrals dominated by
low-$k$ modes (Beattie et al., 2022a), this highlights why the
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\text{kin}}}\right\rangle_{\operatorname{\mathcal{V}}}$
did not have a strong impact from the bulk viscosity in Section 5.

For both the $\text{P}\nu=1.0$ and $\text{P}\nu=\infty$ simulation, we observe
an excess amount of power piling up at high-$k$, which could be qualitatively
described as a bottleneck effect. This is what we would call the sub-shear-
viscous range Equation 29, which ought to appear when $\text{P}\nu>1$, as in
the right panel of Figure 1. However, as $\text{P}\nu$ decreases to $\leq
0.1$, the sub-shear-viscous range and bottleneck completely disappear and we
are left with a smooth, rather featureless decay into dissipatively-dominated,
low-power Fourier energy modes, as expected from the left panel of Figure 1.
That is, once the bulk viscosity is large enough ($\text{Re}_{\text{bulk}}$
and $\text{P}\nu$ are low enough, which ought to be $\text{P}\nu\approx 1/3$
based on our definition; see footnote 5), the physical dissipation of them
allows them to be removed faster than they are able to move down the cascade
(which is fast $E_{\text{kin}}(k)\sim k^{-2}$; Seta & Federrath 2021;
Federrath 2013; Federrath et al. 2021; Kriel et al. 2023), and hence the
bottleneck and sub-shear-viscous range is destroyed.

###  6.2 Helmholtz decomposed spectra

As we did with the integral statistics, to truly observe the $\text{P}\nu$
effect on the plasma we plot the Helmholtz decomposed spectra in the middle
(compressible) and bottom (incompressible) panels in Figure 9, also showing
the ratio between the two spectra in the inset in the bottom panel. Confirming
what we discussed in the previous section, for $\text{P}\nu>1$,
$E_{\text{kin},\parallel}(k)/E_{\text{kin},\perp}(k)$ varies from
$E_{\text{kin},\parallel}(k)/E_{\text{kin},\perp}(k)\approx 10^{-1}$ on large-
scales, to $E_{\text{kin},\parallel}(k)/E_{\text{kin},\perp}(k)\approx 10^{2}$
on small-scales, highlighting how energised the compressible modes become in
the sub-shear-viscous range when $\text{P}\nu$ is large. Once $\text{P}\nu<1$,
the apparent bottleneck and sub-shear-viscous range disappears, and
$E_{\text{kin},\parallel}(k)/E_{\text{kin},\perp}(k)$ becomes relatively scale
invariant, with some systematic shallowing from large to small-scales,
indicating a faster cascade in the compressible modes. Now we turn to the
individual $E_{\text{kin},\parallel}(k)$ and $E_{\text{kin},\perp}(k)$
spectra.

As expected, the $E_{\text{kin},\parallel}(k)$ spectrum (middle panel of
Figure 9) is strongly influenced by $\text{P}\nu$, with the dissipative
effects beginning at high-$k$ and moving towards low-$k$ as the $\text{P}\nu$
decreases. This is expected since
$\nu_{\text{bulk}}\nabla\cdot(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}=\nu_{\text{bulk}}\nabla\cdot(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel})\mathbb{I}=\nu_{\text{bulk}}\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$
ends up having the local properties of a Laplacian operator on the
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$ modes in the momentum
flux. As discussed in Section 3.2, in the lowest $\text{P}\nu$, where
$t_{\nu_{\text{bulk}}}\sim t_{0}$, the low-$k$ modes are significantly
suppressed by roughly $10^{2}$ in total power compared to the
$\text{P}\nu=\infty$ simulation. Because these are the simulations that
actually modify the growth and saturation (see Figure 5 and Figure 6), it is
clear that as we hypothesised earlier, it is predominately the low-$k$
compressible modes that change the characteristics of the dynamo, in stark
contrast with the high-$k$ incompressible modes that are responsible for all
of the dynamo action via stretching.

In Figure 7 and Figure 8, we saw that the $2^{\rm nd}$ moments of both the
shear and compression growth terms, Equation 35, contracted and the CDFs
steepened under the impact of the bulk viscosity. We hypothesised that this
may be due to the correlation between compressible and incompressible modes.
We now see another signature of mode correlation via $E_{\text{kin},\perp}(k)$
(bottom panel of Figure 9), where a Laplacian-like viscosity (starting at
high-$k$) encroaches upon the spectra, and more strongly as $\text{P}\nu$
decreases. This corresponds to a growing sub-bulk-viscous range of $k$, as we
illustrated in Figure 1. It is suggestive that in this $\text{P}\nu<1$ regime
$k_{\nu_{\text{shear}}}\sim k_{\nu_{\text{bulk}}}$, and888In that
$k_{\nu_{\text{shear}}}$ depends upon $k_{\nu_{\text{bulk}}}$, and the
dependence does not need to be linear. However, we leave the determination of
such relations for future works. that the most viscous
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}(k)$ and
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}(k)$ modes are becoming
increasingly coupled as the sub-bulk-viscous range increases in size. This is
similar to the results found in two dimensions in Touber (2019), where they
found that in the viscous limit bulk viscosity would suppress vorticity
production. We explore this correlation in more detail in the following
section.

###  6.3 The coupling of incompressible and compressible modes

####  6.3.1 Measuring mode coupling

There are a number of different probes for measuring mode coupling. For
example, the bispectrum measures the coupling of the mode phases (Burkhart et
al., 2009; Burkhart, 2021) and the turbulence transfer functions measure the
energy flux (or any kind of flux) moving between different modes (e.g.,
Mininni et al., 2005; Grete et al., 2017; Schmidt & Grete, 2019; Grete et al.,
2020; Squire et al., 2023). We choose a simple measure to probe the linear
coupling between the Fourier amplitudes, namely the cross-correlation
$\mathcal{C}$ between each wave mode pair, $(k,k^{\prime})$, across the time-
evolving $E_{\text{kin},\parallel}(k,t)$ and $E_{\text{kin},\perp}(k,t)$
spectra. We define it,

$\displaystyle\mathcal{C}(k,k^{\prime},t^{\prime})=\int^{\infty}_{-\infty}\operatorname{d}\\!{t}\,\overline{E_{\text{kin},\perp}}(k,t)\overline{E_{\text{kin},\parallel}}(k^{\prime},t^{\prime}),$
(53)

where $t^{\prime}=t+\Delta t$ such that $\Delta t$ is an arbitrary time lag,
and the overline indicates a transform of

$\displaystyle\overline{X}=\frac{X-\left\langle
X\right\rangle_{t/t_{0}}}{\left\langle
X\right\rangle_{t/t_{0}}^{2}-\left\langle X^{2}\right\rangle_{t/t_{0}}},$ (54)

to normalise each $k$ mode by both the mean energy (a function of where the
$k$ modes exists in the cascade) and the size-scale of the temporal
fluctuations, which we may change by orders of magnitude between each of the
spectra ($E_{\text{kin},\parallel}(k)$ temporal intermittency need not be the
same as $E_{\text{kin},\perp}(k)$ temporal intermittency), as demonstrated by
comparing the $1\sigma$ bands for the spectra in Figure 9.

We plot the $\mathcal{C}(k,k^{\prime},\underline{t}^{\prime})$ surfaces for
$\text{P}\nu=\infty$ (top panel), $\text{P}\nu=0.1$ (middle panel) and
$\text{P}\nu=0.001$ (bottom panel) in Figure 10, showing only $t^{\prime}$
that maximises the $|\mathcal{C}(k,k^{\prime},t^{\prime})|$ function,
$\underline{t}^{\prime}$, i.e., maximises the strength of the mode coupling.
In general, across the viscous modes we find that $\underline{t}^{\prime}=0$
maximises the correlation, which is a repercussion of the parabolic, non-
causal diffusion operators. Here $k$ on the $y$ axis corresponds to
$E_{\text{kin},\perp}(k)$ modes, and $k^{\prime}$ on the $x$ axis corresponds
to $E_{\text{kin},\parallel}(k)$ modes, as indicated by Equation 53.

Figure 10: The cross correlation function
$\mathcal{C}(k,k^{\prime},t^{\prime})$ given by Equation 53, between the
incompressible $E_{\text{kin},\perp}(k)$ modes and the compressible
$E_{\text{kin},\parallel}(k^{\prime})$ modes, shown at
$t^{\prime}=\underline{t}^{\prime}$, where $\underline{t}^{\prime}$ is the lag
that maximises the magnitude of the correlation function, i.e., the maximum
linear mode coupling. Each panel shows a different simulation, as annotated in
the panel, where $\text{P}\nu=\infty$ (top panel) has no bulk viscosity,
$\text{P}\nu=0.1$ (middle panel) has enough bulk viscosity to $\sim$ balance
the shear viscosity and $\text{P}\nu=0.001$ is dominated by bulk viscosity
($\text{Re}_{\text{bulk}}=1$). We show the regions of strong correlation with
the red, dashed rectangles, highlighting that as the bulk viscosity increases,
the incompressible modes become increasingly coupled within the expanding sub-
bulk-viscous range (see Figure 1).

Firstly, a generic feature for all simulations is that the lowest $k$ modes
have the weakest coupling with all other modes in the simulations. These are,
of course, the driving modes, where we inject Gaussian increments into
turbulence that are correlated in time. These modes are the most orthogonal,
in the sense that the fluctuations there are not strongly coupled on large-
scales. However, as we move to higher $k$ modes, we begin to see stronger
correlations develop, which we highlight with the red quadrilaterals over the
most strongly correlated modes. Let us now explore why this is the case.

####  6.3.2 Viscously coupled modes

The strongest correlations are at the highest $k$ modes, where the plasma is
dominated by both types of viscosities (yellow regions in Figure 10,
highlighted with red rectangles), i.e., in the expanding sub-bulk-viscous
range, as illustrated in Figure 1. This region starts small, but as
$\text{P}\nu$ decreases, more and more incompressible modes become coupled
with the compressible models, until at $\text{P}\nu=0.001$, the correlation is
strong and persistent up to $k\lesssim 10$, consistent with the shallowing
incompressible spectrum in Figure 9. Let us consider what form this coupling
may take in the kinetic energy equations.

In the Stokes regime, where
$\text{Re}_{\text{shear}}(k)\sim\text{Re}_{\text{bulk}}(k)\rightarrow 0$,
relevant to the most dissipative $k$ modes in Figure 9, we can construct the
dimensionless (following Section 2.2) $1^{\rm st}$ moment kinetic energy
evolution equation for each of the modes, noting that
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}=\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}\cdot(\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp})=\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}=v_{\parallel}^{2}$,
by construction, and similarly
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}=v_{\perp}^{2}$.
They are,

$\displaystyle\frac{\partial\left\langle{E_{\text{kin},\perp}}\right\rangle_{\operatorname{\mathcal{V}}}}{\partial
t}=$
$\displaystyle\overbrace{\text{Re}_{\text{shear}}^{-1}\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}}\right\rangle_{\operatorname{\mathcal{V}}}}^{\begin{subarray}{c}\text{incomp.
viscosity}\end{subarray}}$
$\displaystyle\quad+\underbrace{\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}}}}_{\begin{subarray}{c}\text{comp.
mode \&}\\\ \text{incomp. viscosity coupling}\end{subarray}},$ (55)
$\displaystyle\frac{\partial\left\langle{E_{\text{kin},\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}}{\partial
t}=$
$\displaystyle\overbrace{\text{Re}_{\text{shear}}^{-1}\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}}\right\rangle_{\operatorname{\mathcal{V}}}}^{\begin{subarray}{c}\text{comp.
mode \&}\\\ \text{incomp. viscosity coupling}\end{subarray}}$
$\displaystyle\quad+\underbrace{\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}}}}_{\text{shear
enhanced bulk viscosity}},$ (56)

similar to Equations 4.15-4.16 in Wang et al. (2018b) but we have expanded the
viscous term $\left\langle{D}\right\rangle_{\operatorname{\mathcal{V}}}$ and
included bulk viscosity. Hence, not only do both energy equations support
Laplacian-like dissipation operators for each of the modes projected onto
either the compressible or incompressible velocity, but through the definition
of the bulk viscosity (see footnote 5), the
$\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}$ terms are also
scaled by $\text{Re}_{\text{shear}}^{-1}$. Through both the projections, and
the mixing of Reynolds numbers, clearly within the Stokes regime
(corresponding to the sub-bulk-viscous range of scales) the modes become
coupled.

Because we are keeping $\text{Re}_{\text{shear}}$ constant, and we are
interested in the correlations from the compressible modes in the Equation 55,
we are able to determine that the correlation of interest is the second term
of Equation 55. The strength of this term relative to the regular
incompressible diffusion, which is captured by the first term, is

$\displaystyle\frac{\left(\frac{4}{3}\text{Re}_{\text{shear}}^{-1}+\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}}}}{\text{Re}_{\text{shear}}^{-1}\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}}\right\rangle_{\operatorname{\mathcal{V}}}}$
$\displaystyle=$ $\displaystyle\left(\frac{4}{3}+\frac{1}{\text{P}\nu}\right)$
$\displaystyle\frac{\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\cdot\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}}\right\rangle_{\operatorname{\mathcal{V}}}}\sim\frac{1}{\text{P}\nu},$
(57)

which shows that the strength of the coupling relative to the incompressible
diffusion scales with $\text{P}\nu^{-1}$. hence, as we show in Figure 10, the
viscous mode coupling increases with decreasing $\text{P}\nu$, which is the
key result from this section.

Figure 11: The same as Figure 9, but for the magnetic energy spectrum
normalised by time-dependent integral magnetic energy. Top: The normalised
magnetic energy spectra shown in the kinematic stage of the dynamo. We show
the Kazantsev (1968) $E_{\rm mag}(k)\sim k^{3/2}$ power law in the top panel,
which all simulations, regardless of $\text{P}\nu$, exhibit up until the peak
of the energy spectra. Bottom: The same as the top panel but for the saturated
stage of the dynamo. In the inset panel we show the ratio between the
normalised saturated and kinematic magnetic field spectra, highlighting the
different spectral distributions of energy in each of the regimes (more large-
scale in the saturated stage, and vice-versa in the kinematic stage). With a
red dashed, vertical line, we show the extent of the apparent compressible
mode bottleneck (the sub-shear-viscous range) from the inset in Figure 9 (only
present in simulations with $\text{P}\nu>0.1$). This shows that the
suppression of the bottleneck and the sub-shear-viscous range through the bulk
viscosity results in a suppression of high-$k$ modes in $E_{\rm mag}(k)$.
Figure 12: The integral scale of the magnetic field $\ell_{\text{cor},b}$ in
units of the driving scale $\ell_{0}$, coloured by $\text{P}\nu$. As has been
shown before, the kinematic stage of the dynamo quickly builds a small-scale
field ($\ell_{\text{cor},b}>\ell_{0}$), likely though magnetic field folding,
which is present in both $\mathcal{M}<1$ and $\mathcal{M}>1$ dynamos
(Schekochihin et al., 2004c; Kriel et al., 2023). However, as the nonlinear
stage approaches $\sim 20t_{0}$ the integral scale shifts to larger scales
before becoming stationary again at these larger scales in the saturated
stage. The effect of the bulk viscosity suppresses the small-scale field in
the kinematic stage of the dynamo, possibly through disrupting the viscous
solenoidal modes $k\sim k_{\nu_{\text{shear}}}$ responsible for growing the
magnetic field.

####  6.3.3 Coupling in the cascade

In principle,
$\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})$
will provide further correlations in the cascade, since
$\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})=\nabla\cdot(\rho[\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}]\otimes[\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}])$
has interaction terms in the cascade
$\sim\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}+\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel})$
which facilitate the conversion of compressible mode momentum flux into
incompressible modes
$\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp})$,
and vice versa
$\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\parallel})$.
However, our simulations do not have a particularly extended cascade, with
$\text{Re}_{\text{shear}}=1000$, and
$\ell_{0}\sim\text{Re}_{\text{shear}}^{2/3}\ell_{\nu_{\text{shear}}}\sim
0.1(1000)^{2/3}\ell_{\nu_{\text{shear}}}$, only supporting a factor $\sim 10$
separation999Proportionality constant from Figure 6 in Kriel et al. (2023) and
the $\text{Re}_{\text{shear}}^{2/3}$ relation is for Burgers (1948) turbulence
with $E_{\text{kin}}(k)\sim k^{-2}$ (Schober et al., 2012), which surprisingly
works for the total spectrum and not just the decomposed spectrum, as shown in
Kriel et al. (2023). between $\ell_{0}$ and $\ell_{\nu_{\text{shear}}}$. This
means $2\lesssim k\lesssim 20$ is approximately where the cascade should be,
noting that there can be contamination from the driving modes over a much
larger bandwidth than just exactly where the energy is being injected
(Federrath et al., 2021; Grete et al., 2023). We therefore are unable to probe
the coupling within the cascade, in detail, in our simulations. However, we
can see faint stripes through $2\lesssim k\lesssim 20$ in both
$\text{P}\nu=\infty$ and $\text{P}\nu=0.1$ simulations, where weakly-
correlated, non-local fluctuations move across both the incompressible $k$ and
compressible $k^{\prime}$ modes.

Figure 13: The same as the left panel of the spectrum schematic in Figure 1,
but showing that in the $\text{P}\nu<1$ regime, as $\text{P}\nu$ decreases,
$k_{\nu_{\text{bulk}}}$ shifts to lower $k$ modes, and through the viscosity
correlations in the sub-bulk-viscous range, $k_{\nu_{\text{bulk}}}$ drags
$k_{\nu_{\text{shear}}}$ with it, which then shifts $k_{\eta}$, since
$k_{\eta}$ is depends upon the exact placement of $k_{\nu_{\text{shear}}}$,
$k_{\eta}\sim\text{Pm}_{\text{shear}}^{1/2}k_{\nu_{\text{shear}}}$ (for latest