# Role of magnetic pressure forces in Fluctuation dynamo saturation

Sharanya Sur,1 Kandaswamy Subramanian2,3
1Indian Institute of Astrophysics, 2nd Block, 100 Feet Road, Koramangala,
Bangalore 560034, INDIA
2IUCAA, Post Bag 4, Ganeshkhind, Pune 411007, INDIA
3Department of Physics, Ashoka University, Rajiv Gandhi Education City, Rai,
Sonipat 131029, Haryana, INDIA E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Using magnetohydrodynamic simulations of Fluctuation dynamos in turbulent
flows with rms Mach numbers $\mathcal{M}\approx 0.2,1.1$ and $3$, we show that
magnetic pressure forces play a crucial role in dynamo saturation in
supersonic flows. Firstly, as expected when pressure forces oppose
compression, an increase in anti-correlation between density and magnetic
field strengths obtains even in subsonic flows with the anti-correlation
arising from the intense but rarer magnetic structures. In supersonic flows,
due to stronger compressive motions density and magnetic field strength
continue to maintain a positive correlation. However, the degree of positive
correlation decreases as the dynamo saturates. Secondly, we find that the unit
vectors of $\nabla\rho$ and $\nabla B^{2}$ are preferentially anti-parallel to
each other in subsonic flows. This is indicative of magnetic pressure opposing
compression. This anti-parallel alignment persists in transonic and supersonic
flows at dynamo saturation. However, compressive motions also lead to the
emergence of a parallel alignment in these flows. Finally, we consider the
work done against the components of the Lorentz force and the different
sources of magnetic energy growth and dissipation. We show that while in
subsonic flows, suppression of field line stretching is dominant in saturating
the dynamo, the picture is different in supersonic flows. Both field line
stretching and compression amplifies the field initially. But the growing
magnetic pressure opposes further compression of magnetic flux which then
dominates the saturation of the dynamo.

###### keywords:

dynamo–(magnetohydrodynamics) MHD – turbulence – methods: numerical

††pubyear: 2024††pagerange: Role of magnetic pressure forces in Fluctuation
dynamo saturation–A

##  1 Introduction

It is now well established that dynamically insignificant seed magnetic fields
embedded in a conducting fluid can be amplified to near equipartition strength
by a three-dimensional random or turbulent flow (Batchelor, 1950; Kazantsev,
1968; Haugen et al., 2004; Schekochihin et al., 2004; Brandenburg &
Subramanian, 2005; Subramanian et al., 2006; Federrath et al., 2011; Bhat &
Subramanian, 2013; Sur et al., 2014a; Porter et al., 2015; Federrath, 2016; Xu
& Lazarian, 2016; Sur et al., 2018; Rincon, 2019; Seta et al., 2020; Xu &
Lazarian, 2020; Seta & Federrath, 2021a; Tobias, 2021). The growth of the
field depends on competition between inductive growth and resistive
dissipation which translates to the requirement that the magnetic Reynolds
number ($\rm Rm$) is above a critical instability threshold ${\rm Rm}_{\rm
cr}$. This process, referred to as Fluctuation dynamo (FD) amplifies seed
fields exponentially fast (on eddy-turnover time-scales) by stretching and
compression of field lines driven by random (in time) motions that is supplied
by a random/turbulent flow, until a saturation process becomes important. FDs
are crucial to explain the possible presence of magnetic fields in elliptical
galaxies (Moss & Shukurov, 1996; Sur, 2019; Seta et al., 2021), in young
galaxies at high redshifts (Sur et al., 2018) and in clusters of galaxies
(Subramanian et al., 2006; Cho & Ryu, 2009; Bhat & Subramanian, 2013; Porter
et al., 2015; Marinacci et al., 2018; Vazza et al., 2018; Donnert et al.,
2018; Sur, 2019; Sur et al., 2021). In addition, they are also likely to
operate in disc galaxies, generated from turbulent fluid motions (Pakmor et
al., 2017; Gent et al., 2023) and provide seed fields for large-scale dynamo
action.

Given the ubiquity of the occurrence of fluctuation dynamos in astrophysical
objects, we revisit the mechanism of fluctuation dynamo saturation in
turbulent flows in this work. The standard picture of the non-linear evolution
of the dynamo is the following. For magnetic Prandtl numbers $\rm Pm>1$ and
fluid Reynolds number (${\rm Re}$) large enough such that the velocity
spectrum consists of multiple scales, the dynamo is excited initially by
random motions on small enough scales $\ell$, where $\rm Rm(\ell)>{\rm
Rm}_{\rm cr}$. The amplification of the field due to such eddies continue till
the local magnetic energy density approaches values comparable to the energy
density of fluid motions on similar scales. Further growth of the field by
such eddies is suppressed due to non-linear back reaction arising from the
Lorentz force. However, turbulent eddies at a neighbouring larger scale can
still continue to grow the field as they have a larger energy. Gradually over
time, scale-by-scale saturation of the field occurs with the peak of the
magnetic energy spectrum approaching a fraction of the outer scale (or the
forcing scale of turbulent motions).

Although the above mechanism has been confirmed in numerical simulations by a
number of authors mentioned earlier, the exact manner in which the approach to
a global non-linear steady state occurs still remains unclear. For this
purpose it is instructive to consider a qualitative picture of FD
amplification. In an ideal plasma, $\bm{{{B}}}/\rho\propto\delta\bm{{{r}}}$,
where $\bm{{{B}}}$ is the magnetic field, $\rho$ the plasma density and
$\delta\bm{{{r}}}$ the infinitesimal separation vector of fluid elements along
the field line (Brandenburg & Subramanian, 2005). If $\delta\bm{{{r}}}$
increases in a mean square sense, due to random stretching of the field, and
if $\rho$ is assumed strictly constant as would obtain in an incompressible
flow, the magnetic energy also increases. Suppression of random stretching of
the magnetic field lines by the growing magnetic tension (of the Lorentz
force) has been recognised as the key mechanism ushering in the saturation of
the dynamo (Schekochihin et al., 2002).

However there is another possibility to limit growth in case we do not impose
strict incompressibility, which obtains due to the effect of magnetic
pressure. When a region containing a flux of field lines (or a flux tube) of
length ’$L$’ is stretched, magnetic pressure opposes the compression of the
flux tube cross-section of area ’$A$’. Consequently, $\rho$ decreases to
conserve mass, and the magnetic field strength $B$ then need not increase.
This mechanism could be important for FD saturation in the interstellar medium
(ISM) of galaxies, where turbulence is expected to be supersonic at the
driving scale. Analysis of the alignment between the magnetic field and the
eigenvectors of the rate of strain tensor ($S_{ij}$) in compressible MHD
simulations by Sur et al. (2014b) further showed that compressive motions are
statistically oriented perpendicular to the direction of stretching. In a
similar vein, magnetic pressure can also oppose any purely compressive motions
and prevent the increase of $\rho$ hence $B$. We examine here how important is
this mode of saturation for the FD in random compressible flows, even if the
compressibility effects are small as in subsonic flows. Note however that the
role of magnetic pressure in saturating the dynamo cannot be captured in
numerical simulations which solves the exactly incompressible MHD equations.

The paper is organised as follows. In Section 2 we present in brief the
initial conditions and the set-up of the simulations. Thereafter, results of
our investigation is presented under various sections. To guide the readers,
we begin by first discussing the evolution of the ratio of magnetic to kinetic
energies and the PDFs of the magnetic energy density in Section 3. With the
overarching goal to elucidate the role of magnetic pressure in FD saturation,
we then proceed to explore the time evolution of the $\rho-B$ correlation and
how regions with varying field strengths (as obtained in the FD) and over
densities contribute to this correlation in Section 4. If magnetic pressure
plays a role as described earlier, we expect this correlation to decrease as
the field becomes dynamically important. Next in Section 5, we study in detail
the nature of the alignment between the magnetic field direction, the gradient
of the density, and the components of the Lorentz force, namely the gradient
of the magnetic pressure and the magnetic tension and explore how these
alignments manifest themselves in the kinematic and non-linear stages of
dynamo evolution. In Section 6, we study the work done against the magnetic
tension and gradient of the magnetic pressure and how it underlines the
importance of magnetic pressure in dynamo saturation. Next, similar to
previous studies (e.g., Seta et al., 2020; Seta & Federrath, 2021a), we also
explore the relative importance of stretching, advection, compression and
dissipation on the growth and saturation of magnetic energy in Section 7.
Finally, in Section 8 we conclude with highlights of the key results and
emphasise the role of magnetic pressure forces in dynamo saturation.

##  2 Numerical method and initial conditions

We focus our analysis on three non-ideal MHD simulations of non-helically
forced fluctuation dynamo simulations in the simplest possible numerical
setup. The simulations were performed with an isothermal equation of state and
periodic boundaries using the FLASH code (Fryxell et al., 2000; Benzi et al.,
2008; Eswaran & Pope, 1988) (version 4.2), which solves the three-dimensional
compressible MHD equations using explicit viscosity and resistivity. In all
the simulations, turbulence is driven solenoidally over a range of wave
numbers between $1\leq|\bm{{{k}}}|L/2\pi\leq 3$ (such that the average forcing
wave number $k_{\rm f}L/2\pi=2$ and ’$L$’ is the length of the box) as a
stochastic Orstein-Ulhenbeck (OU) process with a finite time correlation. To
cover different regimes of compressibility, we analyse the data from a
subsonic, a transonic and one supersonic simulation with steady state rms Mach
numbers $\mathcal{M}_{\rm rms}=u_{\rm rms}/c_{\rm s}\approx 0.2,1.1$ and
$\approx 3.0$, respectively. All simulations have a resolution of $512^{3},\rm
Pm=\rm Rm/{\rm Re}=1$ and are in dimensionless units. The initial conditions
consist of a box of unit length, the density $\rho$ and sound speed $c_{\rm
s}$ set to unity with zero initial velocities. The magnetic field is
initialised as $\bm{{{B}}}_{\rm init}=B_{0}[0,0,\sin(10\pi x)]$, with $B_{0}$
adjusted to a value such that the initial plasma beta $\beta_{0}=p_{\rm
th}/p_{\rm mag}\approx 10^{6}-10^{7}$. Thus, the magnetic field is dynamically
weak to start with in all the simulations. Further, to maintain
$\nabla\cdot\bm{{{B}}}$ to machine precision level, we use the unsplit
staggered mesh algorithm in FLASHv4.2 with a constrained transport scheme (Lee
& Deane, 2009; Lee, 2013) and a Harten-Lax-van Leer discontinuities (HLLD)
Riemann solver (Miyoshi & Kusano, 2005). Table 1 lists the important
parameters of these runs.

Table 1: Key parameters of simulations used in this study. The resolution in each case is $512^{3}$. $k_{\rm f}L/2\pi=2$ is the average forcing wave number and $\mathcal{M}_{\rm rms}$ is the average values of the rms Mach number in the steady state. $\ell_{\rm f}=2\pi/k_{\rm f}$ is the forcing scale and $\rm Pm$ and ${\rm Re}$ are the magnetic Prandtl number and the fluid Reynolds numbers, respectively. Run | $k_{\rm f}L/2\pi$ | $\mathcal{M}_{\rm rms}$ | $\rm Pm$ | ${\rm Re}=u\,\ell_{\rm f}/\nu$
---|---|---|---|---
A | 2.0 | $\approx 0.2$ | 1.0 | 1500
B | 2.0 | $\approx 1.1$ | 1.0 | 1250
C | 2.0 | $\approx 3.0$ | 1.0 | 2250

##  3 Magnetic energy evolution and its probability distribution function

Figure 1: Time evolution of the ratio of the magnetic to kinetic energies for
all the runs.

In Fig. 1 we show the time evolution of the ratio of magnetic to the density-
weighted kinetic energies as a function of the eddy-turnover time. Compared to
$\mathcal{M}_{\rm rms}=1.1$ and $3$, the subsonic run with $\mathcal{M}_{\rm
rms}=0.2$ has a short kinematic phase of growth up to $\approx 9t_{\rm ed}$.
Thereafter, it has a saturated phase extending from $\approx 13-25t_{\rm ed}$.
We refer the phase in between these two phases as an intermediate one. In run
B ($\mathcal{M}_{\rm rms}=1.1$), the kinematic phase extends up to $\approx
12t_{\rm ed}$, followed by a short intermediate phase of slower growth and
then a saturated phase from $\approx 21t_{\rm ed}$. On the other hand, run C
($\mathcal{M}_{\rm rms}=3.0$) shows a kinematic phase up to $\approx 21t_{\rm
ed}$, followed by a period of slower growth up to $\approx 30t_{\rm ed}$.
Beyond this, the dynamo approaches the non-linear phase with $E_{\rm m}/E_{\rm
k}$ still showing a marginal increase till the time the simulation was evolved
keeping computational costs within reasonable limits.

Figure 2: PDFs of $(B/B_{\rm rms})^{3}P(B/B_{\rm rms})$ in the kinematic
(blue, dashed) and saturated phases (black, solid) for runs A (top panel) and
C (bottom panel).

A distinct feature of FD is that it generates both less intense volume-filling
fields together with much rarer but strong fields. To illustrate our point, we
show in Fig. 2, the probability distribution functions (PDFs) of $(B/B_{\rm
rms})^{3}P(B/B_{\rm rms})$ in the kinematic (blue, dotted lines) and saturated
phases (black, solid lines) for run A (subsonic) in the top panel and run C
(supersonic) in the bottom panel. This measures the magnetic energy density
contributed by fields in a logarithmic interval of $B/B_{\rm rms}$ where
$P(B/B_{\rm rms})$ is the PDF of the magnetic field strength. Two features are
evident from these plots. Firstly, in the kinematic phase, magnetic fields of
a wide range of strengths are present and contribute to the total energy.
Especially in the supersonic case, we find fields with strengths $\approx
10B_{\rm rms}$ contribute comparable energy density to the total as the rms
fields. This is not the case with the subsonic run. This feature could be due
to the continuous amplification of fields by density compression in addition
to amplification by stretching. In the nonlinear phase, the peak of the PDFs
in both subsonic and supersonic cases lies at $\approx 2\,B/B_{\rm rms}$.
However, strong field regions with $B/B_{\rm rms}>2$ also contribute to the
magnetic energy with a probability of $10^{-2}$ for $B/B_{\rm rms}=6$ in the
supersonic case. This implies that together with the general sea of volume
filling fields, strong field regions may also play a role in dynamo
saturation, particularly when compressibility effects are significant. We also
note that in the supersonic case, regions with much higher field strengths
($B/B_{\rm rms}>6$) continue to exist albeit with a much lower probability.
This could imply that even though the dynamo suppresses the compressive
motions when the fields become dynamically important, density compression
still manages to amplify fields in the rarer regions.

##  4 Correlation between density and magnetic field strength

In our simulations, density and magnetic field strength are uncorrelated to
start with. To explore the temporal evolution of the (anti)correlation between
them as the FD evolves, we compute the Pearson correlation coefficient
($r_{p}$) defined as,

$\displaystyle r_{p}(\rho,B)=\frac{{\rm
Cov}(\rho,B)}{\sigma_{\rho}\sigma_{B}}$ (1)
$\displaystyle=\frac{\sum_{i,j,k}(\rho_{i,j,k}-\bar{\rho})(B_{i,j,k}-\bar{B})}{\sqrt{\sum_{i,j,k}(\rho_{i,j,k}-\bar{\rho})^{2}}\,\sqrt{\sum_{i,j,k}(B_{i,j,k}-\bar{B})^{2}}},$

where $B=\sqrt{B_{x}^{2}+B_{y}^{2}+B_{z}^{2}}$ is the magnitude of the field,
$\rho_{i,j,k}$ and $B_{i,j,k}$ are the density and the magnetic field strength
at a point $(i,j,k)$ in the simulation volume. $\bar{\rho}$ and $\bar{B}$ are
the mean values of density and $B$, respectively.

Fig. 3 shows that in the subsonic case (run A, blue dotted line with
asterisks) density and magnetic field strength evolves to become anti-
correlated after about one eddy-turnover time ($t_{\rm ed}=l_{\rm f}/u_{\rm
rms}$) and eventually settles at an average value $\langle
r_{p}\rangle\approx-0.3$ in the saturated state of dynamo. Thus, when the
effects of compressibility are weak, high density regions correspond to low
magnetic field strength regions and vice versa. This is a consequence of the
magnetic pressure acting to oppose the compression of the flux tube.

When compressibility effects become important (with the increase in
$\mathcal{M}_{\rm rms}$) the growth of the field is expected to be driven by a
combination of amplification by random stretching and compression in
converging flows which also enhances the density. The former amplifies fields
in regions where the vorticity is strong and the latter in regions where over-
densities are large. However, regions of strong vorticity may not correspond
to regions with large over-densities. In the kinematic phase when fields are
dynamically weak, amplification due to density compressions seem to dominate
the evolution of the correlation co-efficient. This leads to an initial
positive correlation with the degree of correlation strongly dependent on the
rms $\mathcal{M}$ number of the flow as shown also by Yoon et al. (2016) and
Seta & Federrath (2021b).

Indeed, in the transonic case (run B, black line with squares), we observe a
positive correlation with $\langle r_{p}\rangle\approx 0.35$ which then
gradually declines as the dynamo saturates. After about $25t_{\rm ed}$, the
density and the magnetic field strength becomes uncorrelated. Thus, with
exception to the initial positive correlation, the evolution of
$r_{p}(\rho,B)$ in the transonic case bears resemblance to the corresponding
evolution in the subsonic case. The decrease in $r_{p}$ can be understood as a
consequence of magnetic pressure forces gaining in importance, which acts to
further resist the compressive motions resulting in suppression of density
enhancements. This results in regions of strong fields amplified by vortical
motions, that are not necessarily associated with high density regions which
arose from compression.

Figure 3: Evolution of $r_{p}(\rho,B)$ for runs A, B and C. The blue dotted,
red dashed and black solid vertical lines indicate the transition to the non-
linear phase in the respective runs.

On the other hand, in the supersonic case (run C, red dotted lines with
triangles) where $\delta\rho/\bar{\rho}\propto\mathcal{M}^{2}\sim 9$, the
dominance of density compressions amplifying the field pushes the $\langle
r_{p}\rangle\approx 0.6$ up to $\approx 20\,t_{\rm ed}$, which corresponds to
the duration of the kinematic phase. Thereafter, the correlation decreases
only slightly as the field evolves to the non-linear saturated state with
$\langle r_{p}\rangle\approx 0.45$ in the interval $t/t_{\rm ed}=(30-48)$.
This suggests that when compressible effects are strong, field amplification
due to density compression dominates the evolution of $r_{p}$ even when
Lorentz forces are strong enough to be dynamically significant.

###  4.1 Evolution of the correlation coefficient for different ranges in
$B/B_{\rm rms}$ and over densities

A generic feature of FD is that it generates magnetic fields consisting of
rarer, intense field structures embedded in a sea of less intense, volume
filling fields. Therefore, it is worthwhile to distinguish the contribution to
$r_{p}(\rho,B)$ seen in Fig. 3 from regions with differing magnetic field
strengths. Note here that for a given range of $B/B_{\rm rms}$ we focus only
on those regions in the simulation volume where the range condition is
satisfied. Thus, the data in each of the ranges discussed below correspond to
a subset of the full data set and Eqn. 1 is solved for the respective data
sets.

Figure 4: Evolution of $r_{p}(\rho,B)$ for runs A, B and C for different
ranges of $B/B_{\rm rms}$. The dotted vertical lines denote the transition to
the non-linear phase in the respective runs.

In figure 4, we show the evolution of $r_{p}$ for the subsonic (top panel),
transonic (middle) and supersonic (bottom panel) for two different ranges of
$B/B_{\rm rms}$ while the black solid lines in each panel depicts the
variations of $r_{p}$ as in Fig. 3 where no such ranges in $B/B_{\rm rms}$ are
considered. For the subsonic case, it is now clearly evident that the strong
anti-correlation $\langle r_{p}\rangle\approx-0.46$ in the range $t/t_{\rm
ed}=13-25$ arises from the rarer, intense field regions where $B/B_{\rm
rms}>3$ (green dashed line with asterisk). In the absence of strong density
compressions, such intense magnetic field regions solely arise due to random
stretching. Thus, the observed strong anti-correlation implies that magnetic
pressure exerted by these strong fields resist further compression of flux
tubes. At the same time, the sea of less intense volume filling fields
however, remain uncorrelated with the density at all times.

In the transonic case, initial density enhancements due to compression also
enhances the magnetic field. This leads to an initial positive correlation,
although the degree of positive correlation is lesser when regions with fields
$B/B_{\rm rms}>3$ are considered. However, the plots clearly show that once
the dynamo builds up the fields, the degree of positive correlation steadily
decreases even in regions where $B/B_{\rm rms}\leq 1$, while for the strong
field regions, $r_{p}(\rho,B)$ evolves to be anti-correlated. Thus, with
exception to the initial positive correlation, the general trend of
$r_{p}(\rho,B)$ is very similar to the one observed in the subsonic case.
Finally, in the supersonic case $\rho$ and $B$ remain positively correlated
throughout the evolution. However, regions with $B/B_{\rm rms}>3$ only show a
weaker degree of positive correlation compared to regions with $B/B_{\rm
rms}\leq 1$. This again reflects the fact that such strong and intense field
regions arise due to vorticity. However, density enhancements due to strong
compressions continue to dominate over the magnetic field enhancements due to
vortical motions leading to a net albeit weaker positive correlation.

Figure 5: Evolution of $r_{p}(\rho,B)$ for runs A, B and C for over density
ranges segregated by regions : $\delta\rho/\overline{\rho}<\mathcal{M}^{2}$
and $\delta\rho/\overline{\rho}\geq\mathcal{M}^{2}$. As before, the dotted
vertical lines in all the panels indicate the transition to the non-linear
phase in the respective runs.

The fact that we consider three simulations with varying degrees of flow
compressibility also makes it equally important to probe the degree of
(anti)-correlation in different over density ranges. Recall that density
inhomegenities $\delta\rho/\overline{\rho}\propto\mathcal{M}^{2}$, where
$\delta\rho=(\rho-\bar{\rho})$ is the fluctuation in density. While
$\delta\rho/\overline{\rho}$ is negligible in the subsonic case, significant
density fluctuations $\propto\mathcal{M}^{2}$ could influence the correlation
in transonic and supersonic regimes. To explore this, we compute the evolution
of $r(\rho,B)$ in density ranges segregated by regions :
$\delta\rho/\overline{\rho}<\mathcal{M}^{2}$ and
$\delta\rho/\overline{\rho}\geq\mathcal{M}^{2}$. The results are shown in
figure 5 with panel (a) showing the evolution for the subsonic, (b) for
transonic and panel (c) for the supersonic run. We find that in the subsonic
case, the observed anti-correlation arises from regions with
$\delta\rho/\overline{\rho}<0.04$ as the magnetic fields in these regions are
amplified solely due to vortical motions. Regions with
$\delta\rho/\overline{\rho}\geq 0.04$ remain uncorrelated. In the transonic
and supersonic cases, the evolution of $r_{p}$ in regions with
$\delta\rho/\overline{\rho}<\mathcal{M}^{2}$ follow the pattern seen when no
cutoffs are introduced (black solid lines). On the other hand, in both cases,
$r_{p}$ starts out with a much weaker positive correlation when over dense
regions with $\delta\rho/\overline{\rho}\geq\mathcal{M}^{2}$ are considered.
In the transonic case, magnetic field strength eventually becomes uncorrelated
with the density while in the supersonic case, they maintain a weaker positive
correlation. This could be due to the fact that such high density regions are
fewer in number in the simulation domain.

##  5 Alignment angles

The evolution of $r_{p}(\rho,B)$ offers a first glimpse into the possible role
of magnetic pressure in saturating the dynamo. Additionally it is also
worthwhile to investigate the nature of the alignments between the density
gradient, magnetic field and the components of the Lorentz force. Accordingly,
we show in Fig. 6, the PDFs of the cosine of the angles between : the magnetic
field and the gradient of density (left column), the gradient of density and
gradient of the magnetic pressure (middle column) and finally, between the
gradient of the density and the magnetic tension (right column). These PDFs
are computed over a range of eddy-turnover times in the kinematic phase (blue,
dash-dotted) and in the saturated phase (red, dotted).

Figure 6: PDFs of the cosine of the angles for the subsonic
($\mathcal{M}=0.2$) in the top row, transonic run ($\mathcal{M}=1.1$) in the
middle and for the supersonic ($\mathcal{M}=3.0$) run in the bottom row. The
alignments are shown between : $n_{\bm{{{B}}}}\cdot n_{\nabla\rho}$ (left
column), $n_{\nabla\rho}\cdot n_{\nabla B^{2}}$ (middle column) and
$n_{\nabla\rho}\cdot n_{(\bm{{{B}}}\cdot\nabla)\bm{{{B}}}}$ (right column).
The $x$-axis in the middle column range from -1.0 to +1.0 from left to right.

The left column of the above figure shows that the distributions of the unit
vectors of $\nabla\rho$ and $B$ are symmetric about $\cos(\theta)=0$ with the
unit vector of the magnetic field (directed along the field line)
preferentially perpendicular ($90^{\circ}$ or $270^{\circ}$) to the unit
vector of the density gradient. This is consistent with the simple picture of
a flux tube threaded by a field of strength ’$B$’ where the density can be
enhanced or depleted in the flux tube. Moreover, the degree of orthogonal
alignment appears to be stronger in the saturated phase in comparison to the
kinematic phase. On the other hand, the gradient of the magnetic pressure
($\nabla B^{2}$) is directed inwards to the flux tube. The middle column in
the figure shows that in the subsonic case (run A), $\nabla\rho$ and $\nabla
B^{2}$ are anti-parallel (i.e., at $180^{\circ}$) to each other in both
kinematic and saturated phases. This is a manifestation of the fact that as
magnetic pressure gains in importance due to dynamo amplification, the force
due to magnetic pressure opposes the compression of the flux tubes. This
results in density variations that are anti-correlated with variations in
magnetic pressure. Thus, force due to magnetic pressure tends to empty out the
flux tubes rather than allowing the cross-sectional area to decrease. The
kinematic phase in the transonic case (run B) shows an anti-parallel alignment
together with a weak parallel alignment between the two, which evolves to a
stronger anti-parallel alignment in the saturated phase.

However, the picture changes considerably in the supersonic case (run C) where
compressibility effects are strong enough that in the kinematic phase, density
enhancements drive the amplification of the magnetic field. In this case, the
strong positive correlation seen in Fig. 3 results in density variations that
are also correlated with magnetic pressure variations. Consequently,
$\nabla\rho$ and $\nabla B^{2}$ are initially aligned parallel to each other.
However, over time random stretching also amplifies magnetic fields in
addition to field amplification due to compressions. Similar to the subsonic
case, we find that as the dynamo saturates force due to magnetic pressure
opposes further compression of the field lines resulting in the emergence of
an anti-parallel alignment between the two that coexists with the parallel
alignment111The (anti)parallel alignment further implies that the contribution
from the $\nabla\rho\times\nabla B^{2}$ term towards vorticity generation is
statistically insignificant even if magnetic fields are dynamically important
and/or density gradients are large as in compressible turbulence..

The right column shows that magnetic tension does not show any preferred
alignment with the density gradient in the subsonic case. However, with the
increase in flow compressibility, a weak orthogonal alignment between the two
emerges which becomes more prominent at $\mathcal{M}_{\rm rms}=3$. Thus, the
nature of the alignments discussed above provides further appreciation of the
role of magnetic pressure in dynamo saturation.

##  6 Impact of Lorentz forces on magnetic energy evolution

In this section, we discuss how the energy exchange between the velocity and
magnetic fields facilitated by the different components of the Lorentz force
elucidates the role of magnetic pressure and magnetic tension in saturating
the dynamo.

The volume averaged magnetic energy evolution in terms of the work done by the
Lorentz force is,

$\displaystyle\frac{d}{dt}\left\langle{|\bm{{{B}}}|^{2}/2}\right\rangle=$
$\displaystyle-$
$\displaystyle\left\langle\bm{{{U}}}\cdot(\bm{{{J}}}\times\bm{{{B}}})\right\rangle-\left\langle\frac{\bm{{{J}}}^{2}}{\sigma}\right\rangle$
(2) $\displaystyle-$
$\displaystyle\left\langle\nabla\cdot(\bm{{{E}}}\times\bm{{{B}}})\right\rangle$

where $\sigma$ is the conductivity and the angular brackets
$\langle...\rangle$ denotes volume averaging. Substituting
$\bm{{{J}}}\times\bm{{{B}}}=(\bm{{{B}}}\cdot\nabla)\bm{{{B}}}-\nabla(B^{2}/2)$
in the above equation and multiplying both sides by $t_{\rm ed}/B_{\rm
rms}^{2}$, the magnetic energy evolution in dimensionless form is,

$\displaystyle\frac{d}{dt}\left\langle\frac{|\bm{{{B}}}|^{2}}{2}\right\rangle\frac{t_{\rm
ed}}{B_{\rm rms}^{2}}=$ $\displaystyle-$
$\displaystyle\left\langle\bm{{{U}}}\cdot\bm{{{L}}}_{\rm
T}\right\rangle\frac{t_{\rm ed}}{B^{2}_{\rm rms}}$ (3) $\displaystyle+$
$\displaystyle\left\langle\bm{{{U}}}\cdot\bm{{{L}}}_{\rm
p}\right\rangle\frac{t_{\rm ed}}{B^{2}_{\rm
rms}}-\left\langle\frac{\bm{{{J}}}^{2}}{\sigma}\right\rangle\frac{t_{\rm
ed}}{B_{\rm rms}^{2}}.$

The first two terms on the r.h.s correspond to the rate of work done by the
magnetic tension $\bm{{{L}}}_{\rm T}=\tilde{\bm{{{L}}}}_{\rm
T}-(\hat{\bm{{{B}}}}\cdot\tilde{\bm{{{L}}}}_{\rm T})\hat{\bm{{{B}}}}$ and the
gradient of the magnetic pressure $\bm{{{L}}}_{\rm p}=\tilde{\bm{{{L}}}}_{\rm
p}-(\hat{\bm{{{B}}}}\cdot\tilde{\bm{{{L}}}}_{\rm p})\hat{\bm{{{B}}}}$. Here,
$\tilde{\bm{{{L}}}}_{\rm T}=(\bm{{{B}}}\cdot\nabla)\bm{{{B}}}$ and
$\tilde{\bm{{{L}}}}_{\rm p}=\nabla(B^{2}/2)$ and $\hat{\bm{{{B}}}}$ is the
unit vector of the magnetic field. Thus, $\bm{{{L}}}_{\rm T}$ and
$\bm{{{L}}}_{\rm p}$ as defined above only retain components perpendicular to
$\bm{{{B}}}$. This is because the Lorentz force has no component parallel to
$\bm{{{B}}}$. The last term represents the decrease in magnetic energy due to
Joule dissipation.222The term proportional to the Poynting vector
$\bm{{{E}}}\times\bm{{{B}}}$ goes to zero as the surface integral vanishes
over the volume.

Figure 7: Evolution of the first and second term in Equation 3 for runs A (top
panel), B (middle panel) and C (bottom panel).

Using the instantaneous values of $u_{\rm rms}$ and $B_{\rm rms}$ for the
normalisation parameter $t_{\rm ed}/B_{\rm rms}$, we show in Fig. 7 the time
evolution of the average values of the first term (blue dotted line with
diamonds), and the second term (red, dashed line with ’+’ symbol). In all the
three cases corresponding to different degrees of flow compressibility, both
the terms have positive average values which implies that they contribute to
the growth in magnetic energy (see Eqn. 3). In the subsonic and transonic
cases (top and middle panel), the work done against the gradient of the
magnetic pressure is always small compared to the magnetic tension. The plots
also show that both magnetic pressure and magnetic tension are involved in
dynamo saturation. We see that the average values of both the terms decrease
as the dynamo evolved from the kinematic to the saturated phase. In the
subsonic case, both magnetic pressure and tension contributions decrease by
identical factors. However, in the transonic case, the suppression of the
magnetic pressure term is more compared to the decrease due to the magnetic
tension.

Interestingly in the supersonic case (bottom panel), the energy gained from
the velocity due to $\left\langle\bm{{{U}}}\cdot\bm{{{L}}}_{\rm
T}\right\rangle$ and $\left\langle\bm{{{U}}}\cdot\bm{{{L}}}_{\rm
p}\right\rangle$ terms are comparable up to $20t_{\rm ed}$. As can be seen
from Fig. 1 this corresponds to the time by when the intermediate phase of
slower growth of the magnetic energy culminates and the non-linear saturation
phase begins. Similar to the transonic case, once the dynamo evolves to the
non-linear phase, the work done against the gradient of the magnetic pressure
is suppressed more in comparison to the magnetic tension. This again reflects
the fact that the magnetic pressure gains in importance in opposing
compressive motions as the dynamo saturates.

Thus, expressing the equation for magnetic energy evolution in terms of the
rate of work done against the components of the Lorentz force offers a fresh
perspective which highlights the significance of magnetic pressure in the
supersonic case.

##  7 Impact of local stretching, advection, compression and dissipation on
magnetic energy evolution

Another approach to understand how the FD saturates is to analyse the effects
of local stretching, advection, compression and dissipation on the growth or
decay of magnetic energy (Seta et al., 2020; Seta & Federrath, 2021a). In this
context, it is important to clarify if suppression of field line stretching is
the only agent responsible for saturating the dynamo. To this end, we take the
dot product of the magnetic induction equation with $\bm{{{B}}}$ and then
expand the $\bm{{{B}}}\cdot[\nabla\times(\bm{{{U}}}\times\bm{{{B}}})]$ term to
obtain the evolution of the mean magnetic energy. In the following, we present
a detailed derivation of the magnetic energy evolution equation and also
highlight a subtle point.

Assuming $\eta$ to be constant, the $i$-th component of the magnetic induction
equation is,

$\displaystyle\frac{\partial B_{i}}{\partial
t}=[\nabla\times(\bm{{{U}}}\times\bm{{{B}}})]_{i}-\eta(\nabla\times\bm{{{J}}})_{i},$
(4)

where $\bm{{{J}}}=\nabla\times\bm{{{B}}}$ is the current density with
$\mu_{0}=1$.

Expanding the cross-product in the induction term using Levi-Civita symbols
yield,

$\displaystyle\frac{\partial B_{i}}{\partial
t}=\epsilon_{ijk}\frac{\partial}{\partial
x_{j}}(\epsilon_{klm}U_{l}B_{m})-\eta(\nabla\times\bm{{{J}}})_{i}.$ (5)

Now, substituting
$\epsilon_{ijk}\epsilon_{klm}=(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl})$
and utilising the properties of Kronecker Delta $\delta_{ij}$’s gives,

$\displaystyle\frac{\partial B_{i}}{\partial t}=\frac{\partial}{\partial
x_{j}}(U_{i}B_{j})-\frac{\partial}{\partial
x_{j}}(U_{j}B_{i})-\eta(\nabla\times\bm{{{J}}})_{i}$ (6)

Taking a dot product of Eqn. 6 with $B_{i}$, and neglecting the term
$\propto\nabla\cdot\bm{{{B}}}$ we obtain the evolution equation of the
magnetic energy in terms of local stretching, advection, compression and
dissipation terms,

$\displaystyle\frac{\partial}{\partial
t}\left(\frac{B^{2}}{2}\right)=\underbrace{B_{i}B_{j}\frac{\partial
U_{i}}{\partial x_{j}}}_{\text{stretching}}-\underbrace{B_{i}\frac{\partial
U_{j}}{\partial
x_{j}}B_{i}}_{\text{compression}}-\underbrace{B_{i}U_{j}\frac{\partial
B_{i}}{\partial x_{j}}}_{\text{advection}}-\underbrace{\eta
B_{i}(\nabla\times\bm{{{J}}})_{i}}_{\text{dissipation}}.$ (7)

Next, we define the rate of the strain tensor $S_{ij}=(U_{i,j}+U_{j,i})/2$
where $U_{i,j}$ denote partial derivative of $U_{i}$ w.r.t. $j$ and vice
versa. Now, expressing the stretching term in Eqn. 7 in terms of $S_{ij}$ and
recasting the third term $B_{i}(U_{j}\partial_{j})B_{i}$ as
$U_{j}\partial_{j}(B_{i}^{2}/2)$ the volume integrated magnetic energy
evolution is given by,

$\displaystyle\int_{V}\frac{1}{2}\frac{\partial|\bm{{{B}}}|^{2}}{\partial
t}~{}\mathrm{d}V$
$\displaystyle=\int_{V}B_{i}B_{j}S_{ij}~{}\mathrm{d}V-\int_{V}\bm{{{U}}}\cdot\frac{1}{2}\nabla|\bm{{{B}}}|^{2}~{}\mathrm{d}V$
$\displaystyle-\int_{V}|\bm{{{B}}}|^{2}(\nabla\cdot\bm{{{U}}})~{}\mathrm{d}V$
$\displaystyle-\eta\int_{V}\bm{{{B}}}\cdot(\nabla\times\bm{{{J}}})~{}\mathrm{d}V.$
(8)

Note that in deriving Eqn. 8, the term $\propto(U_{i,j}-U_{j,i})$ is anti-
symmetric in $(i,j)$ and hence does not contribute to the magnetic energy.

The first term on the r.h.s of Eqn. 8 represents the stretching and
compression of the magnetic field lines by the turbulent flow which may
increase or decrease the magnetic energy. The second term reflects the effects
of advection of the field lines by the turbulent flow. The third term
representing the effects of compression can also reduce the magnetic energy
depending on the divergence of the velocity field and is generally important
in transonic and supersonic flows. The last term is dissipative in nature
which acts to reduce the magnetic energy. The effects of local stretching,
advection, compression and dissipation on the magnetic energy growth (or
decay) can then be analysed from the PDFs of each of the terms on the r.h.s of
Equation 8.

Before proceeding further, we wish to elucidate a subtle point. We note that
the advection term in Eqn. 8 can be further simplified to
$[\nabla\cdot(\bm{{{U}}}|\bm{{{B}}}^{2}|)-|\bm{{{B}}}^{2}|(\nabla\cdot\bm{{{U}}})]/2$
which together with the compression term leads to a term
$\propto(\nabla\cdot\bm{{{U}}})$ and a surface term
$\propto(\bm{{{U}}}|\bm{{{B}}}^{2}|)$. While for periodic boundary conditions,
the latter term vanishes, it will still contribute when computing the PDFs at
each mesh point. Similar reasoning is applicable to the dissipation term in
Eqn 8 which can be expressed as a combination of a divergence term
$\propto\nabla\cdot(\bm{{{J}}}\times\bm{{{B}}})$ and a term
$\propto|\bm{{{J}}}|^{2}$. Again, the divergence term vanishes for periodic
boundaries but will contribute to the PDFs at each mesh point. Therefore,
while computing the PDFs it is essential to retain the contribution from the
advection and dissipation terms in the form given in Eqn. 8. This is in
contrast to Seta & Federrath (2021a) where the aforesaid surface terms were
neglected.

Similar to in Sec. 6 we express Eqn. 8 in dimensionless form by multiplying
both sides of the equation by $t_{\rm ed}/B_{\rm rms}^{2}$ which results in,

$\displaystyle\int_{V}\frac{\partial}{\partial
t}\left(\frac{|\bm{{{B}}}|^{2}}{2}\right)\,\frac{t_{\rm ed}}{B_{\rm
rms}^{2}}~{}\mathrm{d}V$ $\displaystyle=+\int_{V}S_{ij}B_{i}B_{j}\frac{t_{\rm
ed}}{B_{\rm rms}^{2}}~{}\mathrm{d}V$
$\displaystyle-\int_{V}\bm{{{U}}}\cdot\frac{1}{2}\nabla|\bm{{{B}}}|^{2}\frac{t_{\rm
ed}}{B_{\rm rms}^{2}}~{}\mathrm{d}V$
$\displaystyle-\int_{V}|{\bm{{{B}}}}|^{2}({\nabla\cdot\bm{{{u}}}})\frac{t_{\rm
ed}}{B_{\rm rms}^{2}}~{}\mathrm{d}V$
$\displaystyle-\eta\int_{V}\bm{{{B}}}\cdot(\nabla\times\bm{{{J}}})\frac{t_{\rm
ed}}{B_{\rm rms}^{2}}~{}\mathrm{d}V.$ (9)

While computing the PDFs of the compression, advection and dissipation terms
are fairly straightforward, to compute the PDFs of the local stretching we
first project the local magnetic field along the three eigenvectors of the
rate of strain tensor ($\bm{e}_{1},\bm{e}_{2},\bm{e}_{3}$) (Brandenburg, 1995;
Sur et al., 2014b; Seta et al., 2020; Seta & Federrath, 2021a). Then local
stretching is computed as, $\Sigma_{i=1,2,3}[\lambda_{i}t_{\rm
ed}(\bm{{{B}}}\cdot\bm{e}_{i})^{2}]/B_{\rm rms}^{2}$, where
($\bm{{{B}}}\cdot\bm{e}_{i}$) denotes the component of the local magnetic
field along the direction of the eigenvectors, and $\lambda_{i}$ are the
eigenvalues corresponding to the eigenvectors $\bm{e}_{i}$.

In the kinematic phase, there is a net growth in the magnetic energy, while in
the saturated phase the magnetic energy is expected to remain constant. It is
expected that this is caused by the mutual adjustments between the local
stretching, compression, advection and dissipation terms. In what follows, we
discuss the effects of each of these terms in two regimes of compressibility
focussing on run A where $\mathcal{M}_{\rm rms}=0.2$ and run C where
$\mathcal{M}_{\rm rms}=3.0$. Since we normalised Eqn. 9 using the
instantaneous values of $u_{\rm rms}$ and $B_{\rm rms}$ we refrain from
computing the total PDFs over the range of $t_{\rm ed}$ in the kinematic and
saturated phases. Instead, we show in Fig. 8 the time evolution of the mean
values of the PDFs of local stretching ($\mu_{\rm s}$), advection ($\mu_{\rm
a}$), compression ($\mu_{\rm c}$) and dissipation ($\mu_{\rm d}$) and the sum
total of the mean values ($\mu_{\rm tot}$) for the subsonic run with
$\mathcal{M}_{\rm rms}=0.2$. In the same vein, Fig. 9 shows the evolution for
the supersonic run with $\mathcal{M}_{\rm rms}=3$. For the sake of
completeness, the PDFs of stretching, advection, compression and dissipation
terms at a few representative times in the kinematic and saturated phases in
shown in Figs. 10 and 11 in Appendix A.

Figure 8: Evolution of the mean values ($\mu$) of the PDFs of local stretching
(blue, asterisks), compression (green, filled stars), advection (red, open
diamonds) and dissipation term (purple, ’+’) for the subsonic run with
$\mathcal{M}_{\rm rms}=0.2$. The sum total of the mean values is depicted by
black crosses.

In the subsonic case (run A), Fig. 8 shows that field amplification in the
kinematic phase is predominantly driven by an increase in local stretching and
also somewhat by advection of the field by the turbulent flow. This is
accompanied by an increase in dissipation. However, $\mu_{\rm c}$ remains
negligible in both kinematic and saturated phases. This reconfirms the fact
that compression plays a negligible role in dynamo saturation in subsonic
flows. The mean values of local stretching, advection and dissipation peaks
$\approx 3-4t_{\rm ed}$ beyond which, the growth rate of the magnetic energy
changes (see Fig. 1). In fact, Fig. 10 show that dissipation is relatively
stronger than both field line stretching and advection initially. The change
in growth rate of magnetic energy is manifested by a decrease in the values of
$\mu_{\rm s},\mu_{\rm a}$ and $\mu_{\rm d}$. The net growth of magnetic energy
in the kinematic phase is confirmed from the non-zero value of the total mean.
Of course, as the dynamo approaches the non-linear phase, the net growth
consistently decreases. In the saturated phase, field line stretching and
dissipation compensate each other while advection is totally suppressed. A
quick comparison between the PDFs in Fig. 10 show that both local stretching
and dissipation are reduced in the saturated phase. Note that the rate of
strain tensor in the stretching term in Eqn. 9 is only affected by magnetic
tension. Thus, the reduction in stretching observed in the saturated phase
suggest that magnetic tension is resisting the further stretching of the field
lines.

Figure 9: Same as in Fig. 8 but now for the supersonic run with
$\mathcal{M}_{\rm rms}=3.0$.

In contrast to the subsonic case, Fig. 9 shows that amplification of the
magnetic field in the kinematic phase in supersonic turbulence (run C) is
dominated by compression of magnetic field lines. Growth of magnetic energy
due to stretching of the field lines remains subdominant to the effects of
compression by almost a factor of two in this phase. The compressible nature
of turbulence ($\mathcal{M}_{\rm rms}=3.0$) drives local density enhancements
which further amplifies the magnetic field. On the other hand, the effects of
advection remain marginal compared to compression, stretching and dissipation.
The saturated phase is marked by a decrease in field line stretching,
compression, advection and dissipation. But field amplification by compressive
motions suffer the strongest suppression as $\mu_{\rm c}$ decreases by a
factor $\approx 3.3$ (also see Fig. 11) compared to a factor of $\approx 1.5$
for the reduction in stretching. This implies that in the compressible regime,
in addition to reduction in field line stretching, advection and dissipation,
the suppression of further compression of field lines due to the action of
magnetic pressure plays a major role in saturating the dynamo. Interestingly,
the aforesaid deduction is consistent with Sur et al. (2014b) where it was
shown that compressive motions of the strain (perpendicular to the direction
of stretching) are suppressed in the saturated phase of the fluctuation
dynamo.

##  8 Conclusions

Over the course of the last two decades, a great deal of progress has been
made to understand the manifold features of fluctuation dynamos and it’s
application to magnetic fields in the ISM of galaxies and in the intracluster
medium (ICM) of galaxy clusters. A majority of these studies were based on
direct numerical simulations that dealt with amplification of dynamically
insignificant seed magnetic fields and subsequent saturation in flows
encompassing various degrees of compressibility. However, the precise way in
which the dynamo evolves to the non-linear saturated state and the factors
that play a deciding role in achieving this steady state remained far from
understood, more so in compressible flows. While in a broader sense, the back
reaction due to the Lorentz force is rightly adjudged to be responsible for
saturating the dynamo, it had been argued that the magnetic tension force is
the main ingredient that suppresses Lagrangian chaos via reduction of field
line stretching which in turns results in the saturation of the dynamo.

Although the role of magnetic tension is appreciated for dynamo saturation in
incompressible flows, a better understanding is required in compressible flows
which are likely to occur in the galactic ISM where density fluctuations can
be significant. In particular, two important questions arise - does magnetic
pressure gradient associated with the Lorentz force also play a role in
saturating the dynamo and, is suppression of field line stretching the only
mechanism responsible for dynamo saturation when the condition of strict
incompressibility is relaxed?

With the aforesaid questions in mind, here we have provided a new perspective
on how fluctuation dynamos saturate. To this effect, we performed numerical
simulations of fluctuation dynamos at $\rm Pm=1$ in turbulent flows driven at
half the scale of the box with the amplitude of turbulent driving adjusted
such that the steady state $\mathcal{M}_{\rm rms}=0.2,1.1$ and $3.0$. We then
analysed the data from these simulations using a variety of probes. These
include : (i) exploring the correlation between the density and magnetic field
strength, (ii) the nature of the alignments between the magnetic field, the
gradient of the density and the components of the Lorentz force, (iii) the
work done against different parts of the Lorentz force to grow the magnetic
field and, (iv) the impact of local stretching, advection, compression and
dissipation on the evolution of the magnetic energy. Our main conclusions are
as follows :

  * •

The PDFs of the magnetic energy density shown in Fig. 2, demonstrate in a
novel manner that together with the general sea of less intense, volume
filling fields, strong fields regions are also likely to play a role in dynamo
saturation in supersonic turbulence. In fact, as Fig. 2 shows strong fields of
strength $\approx 10B_{\rm rms}$ in the supersonic case contribute comparable
energy density as the rms strength fields in the kinematic phase. Even after
saturation, fields upto $4B_{\rm rms}$ contribute significantly to the
magnetic energy density.

  * •

Density and magnetic fields strength are anti-correlated in subsonic flows.
Due to the incompressible nature of the flow, densities cannot be changed by a
large factor. Thus, saturation of the dynamo can happen by suppression of
field line stretching. This is also confirmed from Fig. 8.

  * •

As the effects of compressibility become important (e.g., runs B and C)
magnetic fields are amplified by a combination of random stretching and
compression which also enhances the density. The latter ensures an initial
positive value of $r_{p}$ which then decreases as dynamo saturates (see Fig.
3). In transonic flows, $\rho$ and $B$ tend to become uncorrelated in the
saturated state while in the supersonic case, $r_{p}$ shows a lesser degree of
positive correlation in the saturated phase compared to the value in the
kinematic phase.

  * •

A quick glance at Fig. 4 reveals that the anti-correlation between $\rho$ and
$B$ seen in the subsonic case stems from the strong field regions with
$B/B_{\rm rms}>3$. Such rarer and intense fields also lead to negligible
$r_{p}$ in the transonic case and a much weaker positive correlation in the
supersonic case. Interestingly, we find that regions with
$\delta\overline{\rho}/\overline{\rho}<\mathcal{M}^{2}$ (see Fig. 5) lead to
the anti-correlation in the subsonic case. In the transonic case, these over
dense regions tend to be uncorrelated with the magnetic field while in the
supersonic case, they tend to reduce the degree of positive correlation.

  * •

The fact that magnetic pressure forces play a significant role in dynamo
saturation is clearly evident from the PDFs of the alignment presented in Fig.
6. In the subsonic case we find that ${\bf n}_{\nabla\rho}$ is statistically
aligned anti-parallel to ${\bf n}_{\nabla B^{2}}$, with the alignment becoming
stronger in the saturated phase. In the transonic and more prominently in the
supersonic case, apart from the anti-parallel alignment, a strong parallel
alignment between the two emerges in the kinematic phase. Nevertheless, the
anti-parallel alignment due to magnetic pressure forces dominate in the
saturated phase. In contrast, magnetic tension does not show any preferred
alignment with the gradient of the density irrespective of the flow
compressibility.

  * •

An analysis of the impact of Lorentz forces on the growth of magnetic energy
(Fig. 7) reveals that the contribution from the magnetic pressure gradient is
always small compared to the magnetic tension up to $\mathcal{M}_{\rm
rms}\approx 1.1$. However, in the supersonic regime, the two are comparable
till the time the non-linear saturated phase ushers in. Thereafter, the work
done against the gradient of magnetic pressure (aiding field growth) is
suppressed more in comparison to the magnetic tension. This again illustrates
the importance of magnetic pressure gradients in dynamo saturation for
supersonic turbulence.

  * •

The time evolution of the mean values of the PDFs of local stretching,
advection, compression and dissipation terms are shown in Figs. 8 and 9. They
clearly demonstrates firstly that local stretching of the field lines play a
crucial role in amplifying the field in subsonic flows. And the saturation of
the dynamo in such flows occurs via a reduction in field line stretching. In
contrast, in supersonic flows, amplification of magnetic fields is dominated
by compression of the field lines in the kinematic phase with stretching
playing a subdominant role. As the magnetic energy evolves to equipartition
strengths, the dynamo saturates more due to suppression of compressive
amplification than due reduction in field line stretching.

In a nutshell, the results from our work confirm in detail a new route to
dynamo saturation based on the simple schematic picture of magnetic flux tubes
discussed in Section 1. In incompressible flows where magnetic tension plays a
decisive role in dynamo saturation, magnetic pressure forces nevertheless act
simultaneously to oppose the compression of the flux tubes, especially in rare
but intense magnetic structures. The pivotal role of magnetic pressure forces
in saturating the dynamo comes into prominence in transonic and particularly
in supersonic flows.There, growing magnetic pressure forces resist compression
of flux tubes and lead to regions of strong magnetic flux being emptied out
rather than a decrease of their cross-sectional areas. This fact has not been
sufficiently appreciated so far and could be important as many astrophysical
systems host compressible turbulent flows. Future work extending these
considerations also to cases with $\rm Pm\gg 1$, and where turbulence is
inhomogeneous as expected in galaxies, would be of interest.

## Acknowledgements

S. Sur thanks the Science and Engineering Research Board (SERB) of the
Department of Science & Technology (DST), Government of India, for support
through research grant ECR/2017/001535. He also acknowledges computing time
awarded at the CDAC National Param Supercomputing facility (NPSF), India,
under the grant ‘Hydromagnetic-Turbulence-PR’ and the use of the High
Performance Computing (HPC) resources made available by the Computer Centre of
the Indian Institute of Astrophysics (IIA).

## Data Availability

The data supporting this article are available on reasonable request to the
corresponding author.

## References

  * Batchelor (1950) Batchelor G. K., 1950, Proceedings of the Royal Society of London Series A, 201, 405
  * Benzi et al. (2008) Benzi R., Biferale L., Fisher R. T., Kadanoff L. P., Lamb D. Q., Toschi F., 2008, Physical Review Letters, 100, 234503
  * Bhat & Subramanian (2013) Bhat P., Subramanian K., 2013, MNRAS, 429, 2469
  * Brandenburg (1995) Brandenburg A., 1995, Chaos Solitons and Fractals, 5, 2023
  * Brandenburg & Subramanian (2005) Brandenburg A., Subramanian K., 2005, Phys. Rep., 417, 1
  * Cho & Ryu (2009) Cho J., Ryu D., 2009, Astrophys. J. L., 705, L90
  * Donnert et al. (2018) Donnert J., Vazza F., Brüggen M., ZuHone J., 2018, Space Sci. Rev., 214, 122
  * Eswaran & Pope (1988) Eswaran V., Pope S. B., 1988, Physics of Fluids, 31, 506
  * Federrath (2016) Federrath C., 2016, Journal of Plasma Physics, 82, 535820601
  * Federrath et al. (2011) Federrath C., Chabrier G., Schober J., Banerjee R., Klessen R. S., Schleicher D. R. G., 2011, Phys. Rev. Lett., 107, 114504
  * Fryxell et al. (2000) Fryxell B., et al., 2000, Astrophys. J. S., 131, 273
  * Gent et al. (2023) Gent F. A., Mac Low M.-M., Korpi-Lagg M. J., Singh N. K., 2023, Astrophys. J., 943, 176
  * Haugen et al. (2004) Haugen N. E., Brandenburg A., Dobler W., 2004, Phys. Rev. E, 70, 016308
  * Kazantsev (1968) Kazantsev A. P., 1968, Soviet Journal of Experimental and Theoretical Physics, 26, 1031
  * Lee (2013) Lee D., 2013, Journal of Computational Physics, 243, 269
  * Lee & Deane (2009) Lee D., Deane A. E., 2009, Journal of Computational Physics, 228, 952
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * Miyoshi & Kusano (2005) Miyoshi T., Kusano K., 2005, Journal of Computational Physics, 208, 315
  * Moss & Shukurov (1996) Moss D., Shukurov A., 1996, MNRAS, 279, 229
  * Pakmor et al. (2017) Pakmor R., et al., 2017, MNRAS, 469, 3185
  * Porter et al. (2015) Porter D. H., Jones T. W., Ryu D., 2015, Astrophys. J., 810, 93
  * Rincon (2019) Rincon F., 2019, Journal of Plasma Physics, 85, 205850401
  * Schekochihin et al. (2002) Schekochihin A. A., Maron J. L., Cowley S. C., McWilliams J. C., 2002, Astrophys. J., 576, 806
  * Schekochihin et al. (2004) Schekochihin A. A., Cowley S. C., Taylor S. F., Maron J. L., McWilliams J. C., 2004, Astrophys. J., 612, 276
  * Seta & Federrath (2021a) Seta A., Federrath C., 2021a, Physical Review Fluids, 6, 103701
  * Seta & Federrath (2021b) Seta A., Federrath C., 2021b, MNRAS, 502, 2220
  * Seta et al. (2020) Seta A., Bushby P. J., Shukurov A., Wood T. S., 2020, Physical Review Fluids, 5, 043702
  * Seta et al. (2021) Seta A., Rodrigues L. F. S., Federrath C., Hales C. A., 2021, Astrophys. J., 907, 2
  * Subramanian et al. (2006) Subramanian K., Shukurov A., Haugen N. E. L., 2006, MNRAS, 366, 1437
  * Sur (2019) Sur S., 2019, MNRAS, 488, 3439
  * Sur et al. (2014a) Sur S., Pan L., Scannapieco E., 2014a, Astrophys. J., 784, 94
  * Sur et al. (2014b) Sur S., Pan L., Scannapieco E., 2014b, Astrophys. J. L., 790, L9
  * Sur et al. (2018) Sur S., Bhat P., Subramanian K., 2018, MNRAS, 475, L72
  * Sur et al. (2021) Sur S., Basu A., Subramanian K., 2021, MNRAS, 501, 3332
  * Tobias (2021) Tobias S. M., 2021, Journal of Fluid Mechanics, 912, P1
  * Vazza et al. (2018) Vazza F., Brunetti G., Brüggen M., Bonafede A., 2018, MNRAS, 474, 1672
  * Xu & Lazarian (2016) Xu S., Lazarian A., 2016, Astrophys. J., 833, 215
  * Xu & Lazarian (2020) Xu S., Lazarian A., 2020, Astrophys. J., 899, 115
  * Yoon et al. (2016) Yoon H., Cho J., Kim J., 2016, Astrophys. J., 831, 85

##  Appendix A Evolution of the PDFs of stretching, advection, compression and
dissipation in compressible flows

In Figs. 10 and 11 we show the PDFs of local stretching (in panel ’a’),
advection (in panel ’b’), compression (in panel ’c’) and dissipation (in panel
’d’) at five representative times encompassing the kinematic and saturated
phases of the dynamo. It is apparent from these plots that the standard
deviations ($\sigma$) values are much larger compared to the mean ($\mu$)
values. Fig. 10 clearly shows that the compression term plays a negligible
role in turbulent incompressible flows.

Figure 10: Evolution of the PDFs of the stretching term in panel (a), the
advection term in panel (b), the compression term in panel (c) and the
dissipation term in panel (d) for run A with $\mathcal{M}_{\rm rms}=0.2$
(subsonic). The values of $\mu$ and $\sigma$ shown in the legends corresponds
to the mean and standard deviation of the distribution in each case. Figure
11: Same as in Fig. 10 but now for run C with $\mathcal{M}_{\rm rms}=3.0$
(supersonic).