11institutetext: Owens Valley Radio Observatory, California Institute of
Technology, MC 249-17, Pasadena, CA 91125, USA 22institutetext: Department of
Physics & ITCP, University of Crete, GR-70013, Heraklion, Greece
33institutetext: Institute of Astrophysics, Foundation for Research and
Technology-Hellas, Vasilika Vouton, GR-70013 Heraklion, Greece

# Analytic characterization of sub-Alfvénic turbulence energetics

R. Skalidis<EMAIL_ADDRESS>K. Tassis , 2233 V. Pavlidou , 2233

Magnetohydrodynamic (MHD) turbulence is a cross-field process relevant to many
systems. A prerequisite for understanding these systems is to constrain the
role of MHD turbulence, and in particular, the energy exchange between kinetic
and magnetic forms. The energetics of strongly magnetized and compressible
turbulence has so far resisted attempts to understand them. Numerical
simulations reveal that kinetic energy can be orders of magnitude higher than
fluctuating magnetic energy. We solved this lack-of-balance puzzle by
calculating the energetics of compressible and sub-Alfvénic turbulence based
on the dynamics of coherent cylindrical fluid parcels. Using the MHD
Lagrangian, we proved analytically that the bulk of the magnetic energy
transferred to kinetic energy is the energy that is stored in the coupling
between the ordered and fluctuating magnetic field. The analytical relations
are in strikingly good agreement with numerical data, up to second-order
terms.

##  1 Introduction

Magnetohydrodynamic (MHD) turbulence is involved in a plethora of physical
phenomena (Biskamp, 2003; Beresnyak, 2019; Matthaeus & Velli, 2011; Matthaeus,
2021; Schekochihin, 2020). The interplay between kinetic and magnetic energy
is important for understanding these processes (Goldstein et al., 1995; Ciolek
& Basu, 2006; Kirk et al., 2009; Oughton et al., 2013; Matthaeus et al., 1983;
Zweibel & McKee, 1995; Schekochihin et al., 2007; Cho & Lazarian, 2002;
Federrath et al., 2011). It is challenging to understand the energy exchange
between kinetic and magnetic forms because the MHD equations are nonlinear.
For this reason, several assumptions and approximations are usually employed.

A widely employed approximation is the incompressibility of the gas (Sridhar &
Goldreich, 1994; Goldreich & Sridhar, 1995), although this is only applicable
to a limited number of systems. Compressible MHD turbulence is more complex,
and additional energy terms contribute to the energy cascade. One main
difference in the energy cascade rate of incompressible and compressible
turbulence is that in the latter, the background magnetic field
($\@vec{B}_{0}$) appears with leading-order terms (Banerjee & Galtier, 2013;
Andrés & Sahraoui, 2017). In contrast, the incompressible turbulence energy
cascade is dominated by the increments of the magnetic and velocity
fluctuations, and $\@vec{B}_{0}$ only appears in higher-order statistics (Wan
et al., 2012). This result motivated the hypothesis that $\@vec{B}_{0}$ might
also appear in the total (kinetic and magnetic) fluctuating energy of
compressible MHD turbulence (Andrés & Sahraoui, 2017), whereas in
incompressible turbulence, the total fluctuating energy is dominated by the
fluctuating (second-order) kinetic and magnetic energy.

In incompressible and sub-Alfvénic turbulence, the fluctuating magnetic energy
is completely transferred to kinetic energy, and the volume-averaged
quantities are in equilibrium, $\rho\langle u^{2}\rangle/2\sim\langle\delta
B^{2}\rangle/8\pi$, when turbulence is maintained in a steady state. In
contrast, direct numerical simulations of sub-Alfvénic and compressible
turbulence show that the volume-averaged kinetic energy is much higher than
the second-order fluctuating magnetic energy, $\rho\langle
u^{2}\rangle/2\gg\langle\delta B^{2}\rangle/8\pi$ (Heitsch et al., 2001; Li et
al., 2012a, b), and their relative ratio depends on the amplitude of
$\@vec{B}_{0}$ (Andrés et al., 2018; Lim et al., 2020; Skalidis et al., 2021;
Beattie et al., 2022b). The excess of the kinetic energy suggests that
$\@vec{B}_{0}$ might provide additional energy to the fluid.

The role of $\@vec{B}_{0}$ in the energetics can be intuitively understood
when we decompose the total magnetic field into a background and a fluctuating
component. In incompressible turbulence, the fluctuating magnetic energy comes
only from the perturbations of the magnetic field, which are of second order.
However, in compressible turbulence, the background field appears in the the
total fluctuating magnetic energy due to the coupling between the background
field and magnetic perturbations ($\delta\@vec{B}$). The magnetic coupling,
expressed as $\@vec{B}_{0}\cdot\delta\@vec{B}$, can only be realized in
compressible turbulence (Montgomery et al., 1987; Bhattacharjee & Hameiri,
1988; Bhattacharjee et al., 1998; Fujimura & Tsuneta, 2009) and is the
dominant (first-order) term of the fluctuating magnetic energy.

In sub-Alfvénic and compressible turbulence, numerical data show that
$\@vec{B}_{0}\cdot\delta\@vec{B}$ stores most of the magnetic energy, and that
the kinetic energy approximately reaches equipartition with the fluctuations
of the coupling term (Skalidis & Tassis, 2021; Skalidis et al., 2021; Beattie
et al., 2022b, a). Thus, the magnetic coupling holds the key for understanding
the energetics of strongly magnetized and compressible turbulence. However,
there is still a lack of first-principle understanding of the role of
$\@vec{B}_{0}\cdot\delta\@vec{B}$ in MHD turbulence dynamics and how it
contributes to the averaged energetics.

We present an analytical theory of the role of the coupling potential in the
energy exchange of sub-Alfvénic and compressible turbulence, which is
encountered in systems such as tokamaks (Strauss, 1976, 1977; Zocco &
Schekochihin, 2011), in the interstellar medium (Mouschovias et al., 2006;
Panopoulou et al., 2015, 2016; Planck Collaboration et al., 2016; Skalidis et
al., 2022), and the Sun (Verdini & Velli, 2007; Tenerani & Velli, 2017; Kasper
et al., 2021; Zank et al., 2022). We write the Lagrangian of coherent flux
structures (Crowley et al., 2022), which allows us to approximate turbulence
properties in a deterministic manner, and calculate analytically the energy
exchange between kinetic and magnetic forms as a function of the Alfvénic Mach
number ($\mathcal{M_{A}}$). We find remarkable agreement between the
analytically calculated energetics and numerical data. We conclude that the
majority of the fluctuating magnetic energy transferred to kinetic energy is
provided by the coupling between the background and the fluctuating magnetic
field.

Figure 1: Magnetized fluid consisting of multiple coherent cylindrical fluid
parcels. Red arrows show the initial magnetic field morphology. Untwisted
fluid parcels are elongated, $\ell_{\parallel}$ $\gg$ $\ell_{\perp}$, and
their motion is longitudinal along or perpendicular to $\@vec{B}_{0}$. In sub-
Alfvénic turbulence, the motion of these fluid parcels can be decomposed into
two independent velocity components, parallel (black arrows) and perpendicular
(orange arrows) to $\@vec{B}_{0}$.

##  2 Model

We considered a turbulent fluid characterized by the commonly employed
properties: 1) spatial homogeneity, 2) infinite magnetic and kinetic Reynolds
number, and 3) time stationarity. We considered that the fluid consists of
coherent flux tubes (e.g., Fig. 1 in Banerjee & Galtier, 2013) (or fluid
parcels) with coordinates $\left(r(t),\phi(t),z(t)\right)$, as shown in Fig.
1. Cylindrical coordinates are motivated by studies showing that the
properties of strongly magnetized turbulence are axially symmetric, with
$\@vec{B}_{0}$ being the axis of symmetry (Goldreich & Sridhar, 1995; Maron &
Goldreich, 2001). We assumed the following initial conditions: 1) uniform
temperature, 2) uniform density, 3) no bulk velocity, 4) uniform static
magnetic field ($\@vec{B}_{0}$ = $B_{0}\hat{z}$), and 5) no self-gravity. We
henceforth adopt the following notation: $z=\ell_{\parallel}$ and
$r=\ell_{\perp}$, where $\ell_{\parallel}$ and $\ell_{\perp}$ denote parcel
sizes parallel and perpendicular to $\@vec{B}_{0}$ , respectively.

We perturbed the magnetic field of a coherent fluid structure with a length
scale $\ell=\sqrt{\ell^{2}_{\parallel}+\ell_{\perp}^{2}}$ by
$\delta\@vec{B}_{\ell}$ such that $|\@vec{B}_{0}|\gg|\delta\@vec{B}_{\ell}|$,
which applies to sub-Alfvénic turbulence. Magnetic perturbations tend to
redistribute the magnetic flux within a fluid. For ideal-MHD (flux-freezing)
conditions, the magnetic flux is preserved. Thus, the surface of the perturbed
fluid parcel ($\@vec{S_{\ell}}$) follows the magnetic field lines. The motion
of the field lines, and hence of $\@vec{S_{\ell}}$, can be either parallel or
perpendicular to $\@vec{B}_{0}$ (Fig. 1): 1) Squeezing and stretching of
$\@vec{S_{\ell}}$ along $\@vec{B}_{0}$ leads to parallel motions,
$\dot{\ell}_{\parallel}\neq 0$. 2) Fluctuations of $\ell_{\perp}$ lead to
perpendicular motions, $\dot{\ell}_{\perp}\neq 0$. Finally, 3) twisting leads
to rotational motions, $\dot{\phi}\neq 0$. This naturally defines
$\ell_{\parallel}$ and $\ell_{\perp}$ as the coherence lengths of the
perturbed volume parallel and perpendicular to $\@vec{B}_{0}$ , respectively.
We focused on large scales since coherent structures are prominent there (De
Giorgio et al., 2017). We invoke as a boundary condition a local environment
beyond $\ell$ (pressure wall).

The flux freezing theorem is

$\frac{d\@vec{B_{\ell}}}{dt}\cdot\@vec{S_{\ell}}=-\<EMAIL_ADDRESS>(1)

The cross section of the coherent volume perpendicular and parallel to
$\@vec{B}_{0}$ is
$\@vec{S}_{\perp,\ell}=2\pi\ell_{\perp}\ell_{\parallel}\hat{r}$, and
$\@vec{S}_{\parallel,\ell}=\pi\ell^{2}_{\perp}\hat{z}$ , respectively. The
cross section related to the rotational motion is
$\@vec{S}_{\phi,\ell}=\ell_{\parallel}\ell_{\perp}\hat{\phi}$. The total
magnetic field in cylindrical coordinates can be expressed as
$\@vec{B_{\ell}}=\delta B_{\perp r,\ell}\hat{r}+\delta
B_{\perp\phi,\ell}\hat{\phi}+\left(B_{0}+\delta
B_{\parallel,\ell}\right)\hat{z}$. From Eq. 1, we obtain that when
$|\@vec{B}_{0}|\gg|\delta\@vec{B}|$, magnetic perturbations along
$\@vec{S}_{\parallel,\ell}$ are associated with a longitudinal motion such
that

$u_{\perp
r,\ell}\equiv\dot{\ell}_{\perp}(t)=-\frac{\delta\dot{B}_{\parallel,\ell}(t)}{2B_{0}+\delta
B_{\parallel,\ell}}\ell_{\perp}(t)\approx-\frac{\delta\dot{B}_{\parallel,\ell}(t)}{2B_{0}}\ell_{\perp,0},$
(2)

where we have considered that the initial dimension of the perturbed volume
$\ell_{\perp,0}$ is much larger than its perturbations. Along
$\@vec{S}_{\perp,\ell}$ , we find that

$u_{\parallel,\ell}\equiv\dot{\ell}_{\parallel}(t)=-\left(\frac{\delta\dot{B}_{\perp
r,\ell}(t)}{\delta B_{\perp
r,\ell}(t)}-\frac{\delta\dot{B}_{\parallel,\ell}(t)}{2B_{0}}\right)\ell_{\parallel}(t)\approx-\frac{\delta\dot{B}_{\perp
r,\ell}(t)}{\delta B_{\perp r,\ell}(t)}\ell_{\parallel}(t),$ (3)

while the azimuthal velocity along $\@vec{S}_{\phi,\ell}$ is

$u_{\perp\phi,\ell}\equiv\dot{\phi}(t)\ell_{\perp}(t)\approx-\left(\frac{\delta\dot{B}_{\perp
r,\ell}(t)}{\delta B_{\perp
r,\ell}(t)}-\frac{\delta\dot{B}_{\perp\phi,\ell}(t)}{\delta
B_{\perp\phi,\ell}(t)}\right)\ell_{\perp}(t).$ (4)

As a result of assuming $|\@vec{B}_{0}|\gg|\delta\@vec{B}|$, we have obtained
that parallel and perpendicular motions are decoupled. The coupling of
parallel and perpendicular motions becomes inevitable when
$|\@vec{B}_{0}|\sim|\delta\@vec{B}|$ (Eq. 3).

In sub-Alfvénic turbulence, magnetic tension dominates magnetic pressure
(Passot & Vázquez-Semadeni, 2003). The high tension suppresses transverse
oscillations due to the strong restoring torques. Thus, twisting would have
minimum contribution to the dynamics (e.g., Longcope & Klapper, 1997) and
motions would be mostly longitudinal ($\dot{\phi}\,,\delta
B_{\perp\phi,\ell}\approx 0$). Since $u_{\perp\phi,\ell}\rightarrow 0$, then
due to Eqs. 3, and 4, $\ell_{\parallel}\gg\ell_{\perp}$, which implies that
untwisted coherent structures are stretched toward the $\@vec{B}_{0}$ axis,
which is consistent with the anisotropic properties of sub-Alfvénic turbulence
(Shebalin et al., 1983; Higdon, 1984; Oughton et al., 1994; Sridhar &
Goldreich, 1994; Goldreich & Sridhar, 1995; Oughton et al., 2013; Oughton &
Matthaeus, 2020; Cho & Lazarian, 2003; Yang et al., 2018; Makwana & Yan, 2020;
Gan et al., 2022).

For untwisted fluid parcels, the perpendicular component of the magnetic
fluctuations has a dominant radial component such that $\delta
B_{\perp,\ell}\approx\delta B_{\perp r,\ell}$. From Eqs. 2 and 3, we derive

$\displaystyle\delta B_{\parallel,\ell}(t)$ $\displaystyle\propto-
B_{0}\log\ell_{\perp}(t),$ (5) $\displaystyle\delta B_{\perp,\ell}(t)$
$\displaystyle\propto\ell_{\parallel}^{-1}(t).$ (6)

The difference in the scaling is due to the Lorenz force by $\@vec{B}_{0}$,
which affects perpendicular motions, while it has no effect on parallel
motions.

##  3 MHD Lagrangian of coherent structures

We write the Lagrangian for the perturbed volume. We place the reference frame
at the center of mass of the target volume, hence there is no bulk velocity
term in the Lagrangian. Therefore, all the velocity components are due to
internal motions induced by magnetic perturbations. We focus on low plasma-
beta fluids111The relative ratio of the thermal and magnetic pressure is the
plasma beta, which is defined as
$\beta=2\mathcal{M}^{2}_{\mathcal{A}}/\mathcal{M}^{2}_{\mathcal{S}}$., which
for sub-Alfvénic turbulence corresponds to high sonic Mach numbers
($\mathcal{M}_{s}$). The perturbed Lagrangian (Newcomb, 1962; Andreussi et
al., 2016; Kulsrud, 2005) of the coherent cylindrical fluid parcel, with
surface $\@vec{S_{\ell}}$, can be split into a parallel and a perpendicular
term (Appendix A),

$\delta\mathcal{L}=\overbrace{\left(\frac{1}{2}\rho\dot{\ell}_{\parallel}^{2}-\frac{\delta
B_{\perp,\ell}^{2}}{8\pi}\right)}^{\delta\mathcal{L}_{\perp}}+\overbrace{\left(\frac{1}{2}\rho\dot{\ell}_{\perp}^{2}-\frac{B_{0}\delta
B_{\parallel,\ell}}{4\pi}-\frac{\delta
B_{\parallel,\ell}^{2}}{8\pi}\right)}^{\delta\mathcal{L}_{\parallel}}.$ (7)

Due to Eqs. 2 and 3, $\delta B_{\parallel,\ell}$, and $\delta B_{\perp,\ell}$
are generalized coordinates of $\delta\mathcal{L}$ and
$\ell_{\parallel}(t)=C/\delta B_{\perp,\ell}(t)$ (Eq. 6), where $C$ is a
constant determined from the initial conditions. With this expression, we
eliminate $\ell_{\parallel}$ from the Lagrangian, which up to second-order
terms is separable into a parallel and a perpendicular part, and is
analytically solvable,

$\displaystyle\delta\mathcal{L}_{\perp}\left(\delta
B_{\perp,\ell},\delta\dot{B}_{\perp,\ell}\right)\approx\frac{1}{2}\rho
C^{2}\frac{\delta\dot{B}_{\perp,\ell}^{2}}{\delta
B_{\perp,\ell}^{4}}-\frac{\delta B_{\perp,\ell}^{2}}{8\pi},$ (8)
$\displaystyle\delta\mathcal{L}_{\parallel}\left(\delta
B_{\parallel,\ell},\delta\dot{B}_{\parallel,\ell}\right)\approx\frac{1}{8}\rho\frac{\delta\dot{B}_{\parallel,\ell}^{2}}{B_{0}^{2}}\ell_{\perp,0}^{2}-\frac{B_{0}\delta
B_{\parallel,\ell}}{4\pi}-\frac{\delta B_{\parallel,\ell}^{2}}{8\pi}.$ (9)

We solve the Euler-Lagrange equations for $\delta\mathcal{L_{\parallel}}$
(Appendix B) and $\delta\mathcal{L_{\perp}}$ (Appendix C) and derive the
analytical solutions of the velocity ($u_{\parallel,\ell}$(t),
$u_{\perp,\ell}$(t) ) and magnetic fluctuations ($\delta
B_{\parallel}$${}_{,\ell}(t)$, $\delta B_{\perp}$${}_{,\ell}(t)$) of
$\@vec{S_{\ell}}$. We find that $\delta B_{\parallel,\ell}$ $\sim t^{2}$,
$u_{\perp,\ell}$ $\sim t$, and $\delta B_{\perp,\ell}$ $\sim t^{-1}$, while
$u_{\parallel,\ell}$ is set by the initial conditions (free streaming of the
gas). We used these analytical solutions in order to calculate the averaged
energetics of a strongly magnetized and compressible fluid.

##  4 Energetics

The total energy of fully developed turbulence is stationary because energy
diffusion is balanced by injection. Time stationarity enables us to
approximate turbulence energetics with the leading-order solutions that we
obtained because our approximations preserve time symmetry, and thus energy is
conserved. The statistical properties of large-scale coherent structures
accurately approximate the volume-averaged turbulent statistical properties.
Thus, for an ergodic fluid (Monin & I’Aglom, 1971; Galanti & Tsinober, 2004),
averaging the turbulent statistical properties over the volume of the fluid at
a given time step $\left(\langle
f\rangle_{\mathcal{V}}=\int_{\mathcal{V}}f\right)$ is approximately equivalent
to averaging over multiple realizations of a typical large-scale coherent
structure $\left(\langle
f_{\ell}\rangle_{\mathcal{T}}=\int_{\mathcal{T}}f_{\ell}\right)$, hence
$\langle f\rangle_{\mathcal{V}}\sim\langle f_{\ell}\rangle_{\mathcal{T}}$,
where $f$ denotes an energy term, and $\mathcal{T}$ corresponds to the
coherent structure crossing time. We next analytically compute the $\langle
f_{\ell}\rangle_{\mathcal{T}}$ energy contribution of each Lagrangian term
(Eq. 7) and their relative ratios. Since coherent cylindrical parcels are
characterized by two different coherence lengths $\ell_{\parallel}$and
$\ell_{\perp}$, they also have two different crossing times:
$\mathcal{T}_{\parallel}$ and $\mathcal{T}_{\perp}$ , respectively. We compare
the $\langle f_{\ell}\rangle_{\mathcal{T}}$ analytical energy ratios with the
$\langle f\rangle_{\mathcal{V}}$ numerical values. The numerical results
correspond to simulations of ideal isothermal MHD turbulence without self-
gravity, and turbulence is maintained in a quasi-static state by injecting
energy with an external forcing mechanism (Beattie et al., 2022b). These
simulation are forced with a mixture of compressible and incompressible modes,
but the driving modes do not affect the energetics of sub-Alfvénic and
compressible turbulence (Skalidis et al., 2021).

###  4.1 Kinetic energy

The total averaged kinetic energy ($\rm E_{kinetic}$) of the coherent fluid
parcel with scale $\ell$ is

$\rm{E_{kinetic}}\equiv\frac{1}{2}\rho\langle
u_{\ell}^{2}\rangle_{\mathcal{T}}=\frac{1}{2}\rho\left(\langle
u_{\perp,\ell}^{2}\rangle_{\mathcal{T}_{\perp}}+\langle
u_{\parallel,\ell}^{2}\rangle_{\mathcal{T}_{\parallel}}\right)\approx\frac{B_{0}\delta
B_{\parallel,\max}}{6\pi}+\frac{\delta B_{\perp,\max}^{2}}{8\pi}.$ (10)

The kinetic energy is dominated to first order by $u_{\perp,\ell}$. Thus, the
average Alfvénic Mach number to first order is

$\mathcal{M_{A}}\equiv\frac{\sqrt{\langle
u_{\ell}^{2}\rangle_{\mathcal{T}}}}{V_{A}}\approx\sqrt{\frac{4\delta
B_{\parallel,\max}}{3B_{0}}}.$ (11)

###  4.2 Harmonic potential

From Eqs. 22 and 26, we find that $\langle\delta
B_{\parallel,\ell}^{2}\rangle_{\mathcal{T_{\perp}}}=7\delta
B_{\parallel,\max}^{2}/15$, and $\langle\delta
B_{\perp,\ell}^{2}\rangle_{\mathcal{T_{\parallel}}}=\delta
B_{\perp,\max}^{2}/2$. The total time-averaged harmonic potential energy ($\rm
E_{harmonic}$) density is equal to

$\rm{E_{harmonic}}\equiv\frac{\langle\delta
B_{\ell}^{2}\rangle_{\mathcal{T}}}{8\pi}\approx\frac{\delta B_{\parallel,\rm
max}^{2}}{8\pi}\left(\frac{7}{15}+\frac{\zeta^{2}(\mathcal{M}_{A})}{2}\right),$
(12)

where $\zeta=\delta B_{\perp,\max}/\delta B_{\parallel,\max}$. Sub-Alfvénic
turbulence is anisotropic (Shebalin et al., 1983; Higdon, 1984; Oughton et
al., 1994; Goldreich & Sridhar, 1995), with the anisotropy between $\delta
B_{\perp}$ and $\delta B_{\parallel}$ depending on $\mathcal{M_{A}}$ (Beattie
et al., 2020). To account for this property, we assumed that $\zeta$ is a
function of $\mathcal{M_{A}}$. When $\mathcal{M_{A}}\rightarrow 0$,
$\@vec{B}_{0}$ suppresses any bending of the magnetic field lines with the
amplitude of $\delta B_{\parallel}$ being larger than that of $\delta
B_{\perp}$ (Beattie et al., 2020), hence $\zeta\rightarrow 0$; this is also a
consequence of $\@vec{\nabla}\cdot\@vec{B}=0$ for anisotropic fluid parcels
with $\ell_{\parallel}$ $\gg$ $\ell_{\perp}$. For $\mathcal{M_{A}}\rightarrow
1$, fluctuations tend to become more isotropic, and hence
$\zeta\rightarrow\sqrt{2}$. These limiting behaviors are consistent with
numerical simulations (Beattie et al., 2020, 2022b).

###  4.3 Coupling potential

According to Eq. 10, $\@vec{B}_{0}\cdot\delta\@vec{B}$ contributes to $\rm
E_{kinetic}$ since

$\rm{E_{coupling}}\equiv\frac{B_{0}\sqrt{\langle\delta
B_{\parallel,\ell}^{2}\rangle_{\mathcal{T_{\perp}}}}}{4\pi}\approx\frac{B_{0}\delta
B_{\parallel,\max}}{6\pi}\approx\rm{E_{kinetic}},$ (13)

to first order. This equation demonstrates that the energy stored in the
coupling potential is in equipartition with the averaged kinetic energy when
turbulence is sub-Alfvénic.

###  4.4 Energetics ratios

The $\rm E_{kinetic}$$/$$\rm E_{coupling}$ ratio is

$\frac{\rm E_{kinetic}}{\rm E_{coupling}}\approx
1+\frac{9}{16}\mathcal{M}^{2}_{A}\zeta^{2}(\mathcal{M_{A}}).$ (14)

For $\mathcal{M_{A}}\rightarrow 0$, $\rm E_{coupling}$ $\approx$ $\rm
E_{kinetic}$, while for $\mathcal{M_{A}}\rightarrow 1$, $\rm E_{kinetic}$
$\gtrsim$ $\rm E_{coupling}$. $\rm E_{kinetic}$ becomes higher than $\rm
E_{coupling}$ because $u_{\parallel}$,ℓ contributes more to $\rm E_{kinetic}$
as $\mathcal{M_{A}}$ increases. When $\mathcal{M_{A}}\rightarrow 1$,
$\zeta\approx\sqrt{2}$ , so that the $\rm E_{kinetic}$/$\rm E_{coupling}$
ratio in trans-Alfvénic turbulence scales as

$\frac{\rm E_{kinetic}}{\rm E_{coupling}}\approx
1+\frac{9}{8}\mathcal{M_{A}}^{2}.$ (15)

Regarding the $\rm E_{harmonic}$/$\rm E_{coupling}$ ratio, we find that

$\frac{\rm E_{harmonic}}{\rm
E_{coupling}}\approx\frac{3}{8}\sqrt{\frac{15}{7}}\mathcal{M}^{2}_{\mathcal{A}}\left(\frac{7}{15}+\frac{\zeta^{2}(\mathcal{M_{A}})}{2}\right),$
(16)

which for the two limiting cases of $\zeta$ ($\mathcal{M_{A}}$) becomes

$\frac{\rm E_{harmonic}}{\rm
E_{coupling}}\approx\left\\{\begin{array}[]{l}0.25\leavevmode\nobreak\
\mathcal{M}^{2}_{\mathcal{A}},\leavevmode\nobreak\ \mathcal{M_{A}}\rightarrow
0\\\ \\\ 0.80\leavevmode\nobreak\
\mathcal{M}^{2}_{\mathcal{A}},\leavevmode\nobreak\ \mathcal{M_{A}}\rightarrow
1\end{array}\right..$ (17)

###  4.5 Comparison between analytical and numerical results

In Fig. 2 we compare the analytically calculated energy ratios with numerical
results from the literature (Beattie et al., 2022b). The lines correspond to
the analytical relations for $\rm E_{harmonic}$/$\rm E_{coupling}$ (Eq. 17)
and $\rm E_{kinetic}$$/$$\rm E_{harmonic}$ (Eq. 15), while the colored points
correspond to the numerical values. The numerical data behave as predicted by
the analytical relations. The scatter of triangles increases at higher
$\mathcal{M_{A}}$ because thermal pressure starts becoming important there,
and hence the contribution of thermal motions to the kinetic energy increases.
In the limit of $\mathcal{M_{A}}$ $\ll 1$, thermal pressure is subdominant and
$\beta\rightarrow 0$. For $\mathcal{M_{A}}$ $=1$, we obtain that
$\beta\rightarrow 0$ when $\mathcal{M}_{s}$ $\gg 1$, while when
$\mathcal{M}_{s}$ $\lesssim 1$, $\beta\rightarrow 1$. Thus, for trans-Alvfénic
turbulence, thermal pressure becomes important only for low $\mathcal{M}_{s}$,
while at high $\mathcal{M}_{s}$ , it has a minor contribution to the
energetics. In our calculations, we neglected thermal pressure, and for this
reason, at $\mathcal{M_{A}}$ = 1, triangles are consistent with the analytical
ratio of $\rm E_{kinetic}$/$\rm E_{coupling}$ (blue line) when
$\mathcal{M}_{s}$ $\geq 2$, while at lower $\mathcal{M}_{s}$, the deviation
between numerical and analytical results increases because $\beta$, hence the
relative contribution of thermal pressure, increases. For $\mathcal{M_{A}}$ ¡
1, $\beta\ll 1,$ and for this reason, the numerical data agree perfectly with
the analytical ratio (blue line). Finally, when we account for the
contribution from both $\@vec{B}_{0}\cdot\delta\@vec{B}$ and $\delta B^{2}$,
the total energy stored in magnetic fluctuations $\left(\rm
E_{m,total}=E_{coupling}+E_{harmonic}\right)$ is very close to equipartition
with kinetic energy, as shown by the red boxes.

Figure 2: Comparison between analytical and numerical results. The solid and
dashed thick black lines correspond to the $\rm E_{harmonic}$/$\rm
E_{coupling}$ratio obtained analytically for $\mathcal{M_{A}}\rightarrow 0$
($\zeta=0$) and $\mathcal{M_{A}}\rightarrow 1$ ($\zeta=\sqrt{2}$),
respectively. Numerical data are shown with colored dots. The blue line
corresponds to the analytically obtained $\rm E_{kinetic}$/$\rm E_{coupling}$
ratio, while colored triangles show the same quantities calculated from
numerical data. The red boxes correspond to $\rm E_{kinetic}$$/\rm E_{\rm
m,total}$. The thin green line shows the energy terms in equipartition. The
color bar shows the sonic Mach number ($\mathcal{M}_{s}$) of the simulations.

##  5 Discussion and conclusions

Analytical calculations of strongly magnetized and compressible (isothermal)
turbulence show that $\@vec{B}_{0}$ appears in the energy cascade with
leading-order terms (Banerjee & Galtier, 2013; Andrés et al., 2018). This is
in striking contrast to incompressible turbulence, where $\@vec{B}_{0}$
appears in higher-order terms (Wan et al., 2012). In the formalism presented
here, the incompressible limit is approximated when
$\@vec{B}_{0}\cdot\delta\@vec{B}$ = 0. In this case, $\@vec{B}_{0}$ does not
appear in the dominant Lagrangian terms, hence the averaged kinetic energy
would scale linearly with the fluctuating magnetic energy, $\delta
u_{\ell}\sim\delta B_{\ell}$ (or equivalently, $\mathcal{M_{A}}\sim\delta
B_{\ell}$). However, in agreement with previous works (Wan et al., 2012), our
formalism shows that $\@vec{B}_{0}$ appears in the energetics of
incompressible turbulence because of the coupling between $\delta
B_{\parallel,\ell}$and $\delta B_{\perp,\ell}$ (Eqs. 3) when higher-order
terms are considered.

For sub-Alfvénic and compressible turbulence, we find that
$\@vec{B}_{0}\cdot\delta\@vec{B}$ is the leading term in the dynamics, and as
a result, the scaling between velocity and magnetic fluctuations becomes
$\delta u_{\ell}\sim\sqrt{B_{0}\delta B_{\ell}}$, or equivalently,
$\mathcal{M}_{\mathcal{A}}\sim\delta u_{\ell}/V_{A}\sim\sqrt{\delta
B_{\ell}/B_{0}}$, which is supported by numerical data (Beattie et al., 2020).
In compressible and strongly magnetized turbulence, compression and dilatation
of the gas locally changes the energy cascade rate (Banerjee & Galtier, 2013).
These local energy fluctuations can only be realized in compressible
turbulence and might be related to the fluctuations of the
$\@vec{B}_{0}\cdot\delta\@vec{B}$ potential. Our analytical results prove that
the total averaged magnetic energy transferred to kinetic is equal to
$\left(2B_{0}\sqrt{\langle\delta B_{\parallel}^{2}\rangle}+\langle\delta
B^{2}\rangle\right)/8\pi$.

The consistency between our analytical relations and numerical data is
remarkable. It is not the first time that simple analytical arguments agree
quantitatively with numerical simulations of nonlinear problems (e.g.,
Mouschovias et al., 2011). However, an analytical theory is always
advantageous because it allows us to achieve a deeper understanding of complex
problems. For this reason, the formalism we presented might offer new insights
into the energetics of strongly magnetized and compressible turbulence. We
hope that it motivates future works about the role of magnetic couplings in
the energy cascade.

## References

  * Andrés & Sahraoui (2017) Andrés, N. & Sahraoui, F. 2017, Phys. Rev. E, 96, 053205
  * Andrés et al. (2018) Andrés, N., Sahraoui, F., Galtier, S., et al. 2018, Journal of Plasma Physics, 84, 905840404
  * Andreussi et al. (2016) Andreussi, T., Morrison, P. J., & Pegoraro, F. 2016, Physics of Plasmas, 23, 102112
  * Banerjee & Galtier (2013) Banerjee, S. & Galtier, S. 2013, Phys. Rev. E, 87, 013019
  * Basu et al. (2009) Basu, S., Ciolek, G. E., Dapp, W. B., & Wurster, J. 2009, New A, 14, 483
  * Beattie et al. (2020) Beattie, J. R., Federrath, C., & Seta, A. 2020, MNRAS, 498, 1593
  * Beattie et al. (2022a) Beattie, J. R., Krumholz, M. R., Federrath, C., Sampson, M., & Crocker, R. M. 2022a, arXiv e-prints, arXiv:2203.13952
  * Beattie et al. (2022b) Beattie, J. R., Krumholz, M. R., Skalidis, R., et al. 2022b, arXiv e-prints, arXiv:2202.13020
  * Beresnyak (2019) Beresnyak, A. 2019, Living Reviews in Computational Astrophysics, 5, 2
  * Bhattacharjee & Hameiri (1988) Bhattacharjee, A. & Hameiri, E. 1988, Physics of Fluids, 31, 1153
  * Bhattacharjee et al. (1998) Bhattacharjee, A., Ng, C. S., & Spangler, S. R. 1998, ApJ, 494, 409
  * Biskamp (2003) Biskamp, D. 2003, Plasma Physics and Controlled Fusion, 45, 1827
  * Cho et al. (2022) Cho, H., Ryu, D., & Kang, H. 2022, ApJ, 926, 183
  * Cho & Lazarian (2002) Cho, J. & Lazarian, A. 2002, Phys. Rev. Lett., 88, 245001
  * Cho & Lazarian (2003) Cho, J. & Lazarian, A. 2003, MNRAS, 345, 325
  * Ciolek & Basu (2006) Ciolek, G. E. & Basu, S. 2006, ApJ, 652, 442
  * Colman et al. (2022) Colman, T., Robitaille, J.-F., Hennebelle, P., et al. 2022, MNRAS, 514, 3670
  * Crowley et al. (2022) Crowley, C. J., Pughe-Sanford, J. L., Toler, W., et al. 2022, Proceedings of the National Academy of Science, 119, e2120665119
  * De Giorgio et al. (2017) De Giorgio, E., Servidio, S., & Veltri, P. 2017, Scientific Reports, 7, 13849
  * Elmegreen (2009) Elmegreen, B. G. 2009, in The Galaxy Disk in Cosmological Context, ed. J. Andersen, Nordströara, B. m, & J. Bland-Hawthorn, Vol. 254, 289–300
  * Elmegreen & Scalo (2004) Elmegreen, B. G. & Scalo, J. 2004, ARA&A, 42, 211
  * Eswaran & Pope (1988) Eswaran, V. & Pope, S. B. 1988, Computers and Fluids, 16, 257
  * Federrath et al. (2011) Federrath, C., Chabrier, G., Schober, J., et al. 2011, Phys. Rev. Lett., 107, 114504
  * Fujimura & Tsuneta (2009) Fujimura, D. & Tsuneta, S. 2009, ApJ, 702, 1443
  * Galanti & Tsinober (2004) Galanti, B. & Tsinober, A. 2004, Physics Letters A, 330, 173
  * Gan et al. (2022) Gan, Z., Li, H., Fu, X., & Du, S. 2022, ApJ, 926, 222
  * Girichidis et al. (2016) Girichidis, P., Walch, S., Naab, T., et al. 2016, MNRAS, 456, 3432
  * Goldreich & Sridhar (1995) Goldreich, P. & Sridhar, S. 1995, ApJ, 438, 763
  * Goldstein et al. (1995) Goldstein, M. L., Roberts, D. A., & Matthaeus, W. H. 1995, Annual Review of Astronomy and Astrophysics, 33, 283
  * Hanasoge et al. (2020) Hanasoge, S. M., Hotta, H., & Sreenivasan, K. R. 2020, Science Advances, 6, eaba9639
  * Heitsch et al. (2001) Heitsch, F., Zweibel, E. G., Mac Low, M.-M., Li, P., & Norman, M. L. 2001, ApJ, 561, 800
  * Higdon (1984) Higdon, J. C. 1984, ApJ, 285, 109
  * Iffrig, Olivier & Hennebelle, Patrick (2017) Iffrig, Olivier & Hennebelle, Patrick. 2017, A&A, 604, A70
  * Kasper et al. (2021) Kasper, J. C., Klein, K. G., Lichko, E., et al. 2021, Phys. Rev. Lett., 127, 255101
  * Kirk et al. (2009) Kirk, H., Johnstone, D., & Basu, S. 2009, ApJ, 699, 1433
  * Klessen & Hennebelle (2010) Klessen, R. S. & Hennebelle, P. 2010, A&A, 520, A17
  * Kritsuk et al. (2017) Kritsuk, A. G., Ustyugov, S. D., & Norman, M. L. 2017, New Journal of Physics, 19, 065003
  * Krumholz & Burkert (2010) Krumholz, M. & Burkert, A. 2010, ApJ, 724, 895
  * Krumholz & Burkhart (2016) Krumholz, M. R. & Burkhart, B. 2016, MNRAS, 458, 1671
  * Kulsrud (2005) Kulsrud, R. M. 2005, Plasma physics for astrophysics
  * Li et al. (2012a) Li, P. S., McKee, C. F., & Klein, R. I. 2012a, ApJ, 744, 73
  * Li et al. (2012b) Li, P. S., Myers, A., & McKee, C. F. 2012b, ApJ, 760, 33
  * Lim et al. (2020) Lim, J., Cho, J., & Yoon, H. 2020, ApJ, 893, 75
  * Longcope & Klapper (1997) Longcope, D. W. & Klapper, I. 1997, ApJ, 488, 443
  * Mac Low et al. (1998) Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Phys. Rev. Lett., 80, 2754
  * Makwana & Yan (2020) Makwana, K. D. & Yan, H. 2020, Physical Review X, 10, 031021
  * Maron & Goldreich (2001) Maron, J. & Goldreich, P. 2001, ApJ, 554, 1175
  * Matthaeus (2021) Matthaeus, W. H. 2021, Physics of Plasmas, 28, 032306
  * Matthaeus et al. (1983) Matthaeus, W. H., Goldstein, M. L., & Montgomery, D. C. 1983, Phys. Rev. Lett., 51, 1484
  * Matthaeus & Velli (2011) Matthaeus, W. H. & Velli, M. 2011, Space Sci Rev, 160, 145
  * McKee (1989) McKee, C. F. 1989, ApJ, 345, 782
  * Monin & I’Aglom (1971) Monin, A. S. & I’Aglom, A. M. 1971, Statistical fluid mechanics; mechanics of turbulence
  * Montgomery et al. (1987) Montgomery, D., Brown, M. R., & Matthaeus, W. H. 1987, J. Geophys. Res., 92, 282
  * Mouschovias et al. (2011) Mouschovias, T. C., Ciolek, G. E., & Morton, S. A. 2011, MNRAS, 415, 1751
  * Mouschovias et al. (2006) Mouschovias, T. C., Tassis, K., & Kunz, M. W. 2006, ApJ, 646, 1043
  * Newcomb (1962) Newcomb, W. A. 1962, Nuclear Fusion Supplement Part, 2
  * Oughton & Matthaeus (2020) Oughton, S. & Matthaeus, W. H. 2020, ApJ, 897, 37
  * Oughton et al. (1994) Oughton, S., Priest, E. R., & Matthaeus, W. H. 1994, Journal of Fluid Mechanics, 280, 95–117
  * Oughton et al. (2013) Oughton, S., Wan, M., Servidio, S., & Matthaeus, W. H. 2013, ApJ, 768, 10
  * Panopoulou et al. (2015) Panopoulou, G., Tassis, K., Blinov, D., et al. 2015, MNRAS, 452, 715
  * Panopoulou et al. (2016) Panopoulou, G. V., Psaradaki, I., & Tassis, K. 2016, MNRAS, 462, 1517
  * Park & Ryu (2019) Park, J. & Ryu, D. 2019, ApJ, 875, 2
  * Passot & Vázquez-Semadeni (2003) Passot, T. & Vázquez-Semadeni, E. 2003, A&A, 398, 845
  * Piontek & Ostriker (2007) Piontek, R. A. & Ostriker, E. C. 2007, ApJ, 663, 183
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 586, A138
  * Schekochihin (2020) Schekochihin, A. A. 2020, arXiv e-prints, arXiv:2010.00699
  * Schekochihin et al. (2007) Schekochihin, A. A., Iskakov, A. B., Cowley, S. C., et al. 2007, New Journal of Physics, 9, 300
  * Shebalin (2013) Shebalin, J. V. 2013, Geophysical and Astrophysical Fluid Dynamics, 107, 411
  * Shebalin et al. (1983) Shebalin, J. V., Matthaeus, W. H., & Montgomery, D. 1983, Journal of Plasma Physics, 29, 525
  * Skalidis et al. (2021) Skalidis, R., Sternberg, J., Beattie, J. R., Pavlidou, V., & Tassis, K. 2021, A&A, 656, A118
  * Skalidis & Tassis (2021) Skalidis, R. & Tassis, K. 2021, A&A, 647, A186
  * Skalidis et al. (2022) Skalidis, R., Tassis, K., Panopoulou, G. V., et al. 2022, A&A, 665, A77
  * Sridhar & Goldreich (1994) Sridhar, S. & Goldreich, P. 1994, ApJ, 432, 612
  * Strauss (1976) Strauss, H. R. 1976, Physics of Fluids, 19, 134
  * Strauss (1977) Strauss, H. R. 1977, Physics of Fluids, 20, 1354
  * Tenerani & Velli (2017) Tenerani, A. & Velli, M. 2017, ApJ, 843, 26
  * Verdini & Velli (2007) Verdini, A. & Velli, M. 2007, ApJ, 662, 669
  * Wan et al. (2012) Wan, M., Oughton, S., Servidio, S., & Matthaeus, W. H. 2012, Journal of Fluid Mechanics, 697, 296–315
  * Yang et al. (2018) Yang, L., Zhang, L., He, J., et al. 2018, The Astrophysical Journal, 866, 41
  * Yang et al. (2021) Yang, Y., Wan, M., Matthaeus, W. H., & Chen, S. 2021, Journal of Fluid Mechanics, 916, A4
  * Zank et al. (2022) Zank, G. P., Zhao, L. L., Adhikari, L., et al. 2022, ApJ, 926, L16
  * Zocco & Schekochihin (2011) Zocco, A. & Schekochihin, A. A. 2011, Physics of Plasmas, 18, 102309
  * Zweibel & McKee (1995) Zweibel, E. G. & McKee, C. F. 1995, ApJ, 439, 779

###### Acknowledgements.

We would like to thank the anonymous referee for reviewing our manuscript. We
are grateful to T. Ch. Mouschovias, and P. F. Hopkins for stimulating
discussions. We also thank E. N. Economou, V. Pelgrims, E. Ntormousi, A.
Tsouros, and I. Komis for useful suggestions on the manuscript. This work was
supported by NSF grant AST-2109127. We acknowledge support by the European
Research Council under the European Union’s Horizon 2020 research and
innovation programme, grant agreement No. 771282 (RS and KT); by the Hellenic
Foundation for Research and Innovation under the “First Call for H.F.R.I.
Research Projects to support Faculty members and Researchers and the
procurement of high-cost research equipment grant”, Project 1552 CIRCE (VP);
and from the Foundation of Research and Technology - Hellas Synergy Grants
Program (project MagMASim, VP, and project POLAR, KT).

##  Appendix A Lagrangian of coherent cylindrical parcels

We used the MHD Lagrangian as derived by Newcomb (Newcomb 1962) for isothermal
magnetized fluids. The total MHD Lagrangian is the sum of the kinetic and the
total potential energy of all the fluid elements within a volume
$\mathcal{V}$,

$\mathcal{L}=\int_{\mathcal{V}}\left(\frac{1}{2}\rho
u^{2}-P_{s}-\frac{B^{2}}{8\pi}-\frac{1}{2}\rho\Phi\right),$ (18)

where $P_{s}$ is the thermal pressure, and $\Phi$ is the gravitational
potential. The equation of motion for magnetized turbulent fluids is obtained
from the stationary-action principle
$\left(\delta\int_{t_{1}}^{t_{2}}dt\mathcal{L}=0\right)$. We focused on fluids
in which magnetic pressure dominates thermal pressure $\left(B^{2}/8\pi\gg
P_{s}\right)$, and we ignored self-gravity, $\Phi=0$. Therefore, the dominant
potential term in the Lagrangian is magnetic pressure. When we consider that
the ensemble of the fluid elements moves as a coherent cylinder (Fig. 1), then
the integration of the Lagrangian takes place within the volume of the
cylinder. In this case, the integrated Lagrangian terms correspond to the
kinetic and the magnetic energy of the cylinder. The perturbed Lagrangian of a
cylinder is (using Eqs. 2, 3, and 4)

$\begin{split}\delta\mathcal{L}=\frac{1}{2}\rho\ell^{2}\left(\frac{\delta\dot{B}_{\parallel}^{2}}{4B_{0}^{2}}+\frac{\delta\dot{B}_{\perp
r}^{2}}{\delta B_{\perp
r}^{2}}\right)+\frac{1}{2}\rho\ell_{\perp}^{2}\left(\frac{\delta\dot{B}_{\perp\phi}^{2}}{\delta
B_{\perp\phi}^{2}}-2\frac{\delta\dot{B}_{\perp r}}{\delta B_{\perp
r}}\frac{\delta\dot{B}_{\perp\phi}}{\delta B_{\perp\phi}}\right)\\\
-\frac{1}{2}\rho\ell_{\parallel}^{2}\frac{\delta\dot{B}_{\perp r}}{\delta
B_{\perp r}}\frac{\delta\dot{B}_{\parallel}}{B_{0}}-\left(\frac{\delta
B_{\parallel}^{2}}{8\pi}+\frac{\delta B_{\perp r}^{2}}{8\pi}+\frac{\delta
B_{\perp\phi}^{2}}{8\pi}\right)-\frac{B_{0}\delta
B_{\parallel}}{4\pi}.\end{split}$ (19)

For untwisted cylinders, all terms containing a $\phi$ component are zero.
Then, the Lagrangian contains only the parallel and the perpendicular
(longitudinal) components, which are generally coupled due to the
$\delta\dot{B}_{\perp r}\delta\dot{B}_{\parallel}/(\delta B_{\perp
r}B_{0})\ell_{\parallel}^{2}$ term. For sub-Alfvénic turbulence, this is a
higher-order term because —$\@vec{B}_{0}$— $\gg$ —$\delta\@vec{B}$—. By
keeping the dominant (second-order) terms, we derived the total perturbed
Lagrangian of a coherent cylindrical structure, which to leading order, can be
expressed as the sum of two independent parts (parallel and perpendicular to
the background magnetic field, Eqs. 8 and 9).

##  Appendix B Solutions of $\delta\mathcal{L_{\parallel}}$

From the Euler-Lagrange equation of $\delta\mathcal{L_{\parallel}}$ , we
obtain

$\delta\ddot{B}_{\parallel,\ell}(t)=-\left(\delta
B_{\parallel,\ell}(t)+B_{0}\right)\frac{4V_{A}^{2}}{\ell_{\perp,0}^{2}},$ (20)

where $V_{A}=B_{0}/\sqrt{4\pi\rho}$ is the Alfvénic speed.

Initially, we compressed the perturbed volume perpendicular to $\@vec{B}_{0}$,
then released it and allowed the compression to propagate. For the initial
conditions, we considered that $u_{\perp,\ell}(t=0)=0$ and $\delta
B_{\parallel,\ell}(t=0)=\delta B_{\parallel,\rm max}$. We might have initiated
the fluid parcel at $\delta B_{\parallel,\ell}(t=0)=-\delta B_{\parallel,\rm
max}$, but in that case, $u_{\perp,\ell}(t=0)\neq 0$. Solutions of Eq. 20 are
harmonic, but are valid only for early times because at later times, nonlinear
interactions become important and energy is diffused. Below, we consider the
scenario of energy diffusing due to the shock formation because we considered
highly compressible fluids. Without loss of generality, we can consider any
diffusive process.

From the jump conditions, we analytically obtained that when $\mathcal{M}_{s}$
$\gg 1$, an isothermal shock perpendicular to $\@vec{B}_{0}$ forms when

$\delta
B_{\parallel}\lesssim\frac{B_{0}}{2}\left(\mathcal{M}^{2}_{A}-1\right).$ (21)

Thus, in sub-Alfvénic turbulence, $\mathcal{M_{A}}$ $<1$, magnetized shocks
form when $\delta B_{\parallel}$ $<0$, which means that $\delta B_{\parallel}$
will never perform a full harmonic cycle. Keeping the dominant term of the
expansion of the harmonic solutions (Eq. 20), we derive that

$\delta B_{\parallel,\ell}(t)\approx\delta
B_{\parallel,\max}-2B_{0}\frac{V_{A}^{2}}{\ell_{\perp,0}^{2}}t^{2}.$ (22)

The above solution through Eq. 2 yields

$u_{\perp,\ell}(t)\approx\frac{2V^{2}_{A}}{\ell_{\perp,0}}t.$ (23)

From Eqs. 22 and 23, we obtain that as the magnetic field of the perturbed
volume decreases, $u_{\perp,\ell}$ increases. When the shock is formed, the
perturbed volume instantaneously bounces off its environment, which acts as a
pressure wall (Basu et al. 2009). At the post-shock phase, the motion is
reversed and the coherent volume will start contracting until until $\delta
B_{\parallel,\ell}$ reaches a value of $+\delta B_{\parallel,\max,p}$. The
post-shock solutions are obtained from Eq. 20 with initial conditions
$u_{p}(t=0)>0$ and $\delta B_{\parallel,p}(t=0)<0$, where the subscript $p$
denotes post-shock quantities. At the post-shock phase, the solution of
$\delta B_{\parallel}$ is

$\delta B_{\parallel,\ell}(t)\approx-\delta
B_{\parallel,p}+u_{p}(t=0)t-2B_{0}\frac{V_{A}^{2}}{\ell_{\perp,0}^{2}}t^{2}.$
(24)

At the post-shock phase, the magnetic field increases until $+\delta
B_{\parallel,\max,p}$, which is smaller than the initial magnetic field
increase ($+\delta B_{\parallel,\max}$) of the pre-shock phase because energy
has been dissipated by the shock (Park & Ryu 2019; Cho et al. 2022). When the
perturbed volume reaches $+\delta B_{\parallel,\max,p}$, the velocity is zero,
and the motion is reversed. Then, the volume starts expanding until it again
forms a shock. Overall, the perturbed volume would perform damped oscillations
until all the energy is dissipated (Basu et al. 2009; Yang et al. 2021).

Fluids in nature are commonly assumed to be constantly perturbed until
turbulence reaches a steady state (Krumholz & Burkert 2010; Kritsuk et al.
2017; Colman et al. 2022). Various driving mechanisms could maintain turbulent
energy in nature (Eswaran & Pope 1988; McKee 1989; Mac Low et al. 1998;
Piontek & Ostriker 2007; Elmegreen 2009; Krumholz & Burkhart 2016; Girichidis
et al. 2016; Hanasoge et al. 2020; Iffrig, Olivier & Hennebelle, Patrick 2017;
Klessen & Hennebelle 2010; Elmegreen & Scalo 2004). In our model, turbulent
driving is equivalent to adding externally kinetic energy to the perturbed
volume, such that the initial velocity at the post-shock phase, $u_{p}(t=0)$,
is sufficient to compress the perturbed volume until it reaches the maximum
compression it had in the pre-shock phase, $\delta
B_{\parallel,\max,p}\approx+\delta B_{\max,\parallel}$.

We considered an external driver, which ensured that $\delta B_{\parallel}$
fluctuations, and hence energy, were maintained in a quasi-static state. In
addition, we considered that the fluid is ergodic (Monin & I’Aglom 1971;
Galanti & Tsinober 2004). For ergodic fluids, $\delta B_{\parallel,\ell}$ are
characterized by ballistic profiles, $\delta B_{\parallel,\ell}\propto t^{2}$,
and as we argue below, they bounce between $+\delta B_{\parallel,\max}$ and
$-\delta B_{\parallel,\max}$ within a characteristic timescale
$\mathcal{T}_{\perp}\approx 4\ell_{\perp,0}V_{A}^{-1}\sqrt{\delta
B_{\parallel,\max}/2B_{0}}$.

When we initially compressed the magnetic field of the perturbed volume along
$\@vec{B}_{0}$, then $\ell_{\perp}$ decreased due to Eq. 5. This forced the
surface of the environment of the perturbed volume to increase by equal
amounts. Thus, the $+\delta B_{\parallel,\max}$ initial increase of the
magnetic field of the perturbed volume causes the magnetic field of the
environment to decrease by $-\delta B_{\parallel,\max}$ due to flux freezing.
If the fluid is ergodic, then different fluid parcels correspond to different
oscillation phases of the target fluid parcel (Monin & I’Aglom 1971; Galanti &
Tsinober 2004). Therefore, the $-\delta B_{\parallel,\max}$ of the environment
corresponds to the maximum decrease in magnetic field strength of the target
volume. Nonlinear effects can break the symmetry between $+\delta
B_{\parallel,\max}$ and $-\delta B_{\parallel,\max}$, but ergodicity is only
weakly broken when $\@vec{B}_{0}$ $\neq 0$ (Shebalin 2013).

The perturbed volume would spend most of its time in the compressed state
because the velocity is minimum there. On the other hand, the velocity of the
fluid parcel is maximum when $\delta B_{\parallel,\ell}<0$, and hence the
fluid parcel would spend minimum time there. As a result, the majority of
fluid parcels at a given time would be compressed ($\delta
B_{\parallel,\ell}>0$) due to ergodicity, which is verified by numerical
simulations (Beattie et al. 2022b).

##  Appendix C Solutions of $\delta\mathcal{L_{\perp}}$

From the Euler-Lagrange equation of $\delta\mathcal{L_{\perp}}$ , we obtain

$\delta\ddot{B}_{\perp,\ell}(t)\delta
B_{\perp,\ell}(t)-2\delta\dot{B}_{\perp,\ell}^{2}(t)+\frac{\delta
B_{\perp,\ell}^{6}(t)}{4\pi\rho C^{2}}=0.$ (25)

For $|\@vec{B}_{0}|\gg|\delta\@vec{B}|$, the sixth-order term above can be
neglected, and then the solutions are straightforward. The total pressure of
the fluid exerted by $\delta B_{\perp}$ is transferred to parallel motions
(Eq. 3), hence $\rho u_{\parallel,\max}^{2}/2=\delta
B_{\perp,\max}^{2}/(8\pi)$. We derive the following solutions:

$\delta B_{\perp,\ell}(t)\approx\frac{fB_{0}}{1\pm
t/\mathcal{T}_{\parallel}},\leavevmode\nobreak\
u_{\parallel,\ell}(t)\approx\pm fV_{A},$ (26)

where $\ell_{\parallel}(t=0)=\ell_{\parallel,0}$, $f=\delta
B_{\perp,\max}/B_{0}\ll 1$, and
$\mathcal{T}_{\parallel}=fV_{A}^{-1}\ell_{\parallel,0}$. In the above
equations, the signs depend on the initial conditions. Initially, we
considered that $\delta B_{\perp,\ell}(t=0)=\delta B_{\perp,\max}$, and
$u_{\parallel,\ell}(t=0)=u_{\perp,\max}$, which leads to positive signs.

If the initial velocity along $\@vec{B}_{0}$ were zero, then both
$u_{\parallel,\ell}$ and $\delta B_{\perp,\ell}$ would remain static. The
coupling of parallel and perpendicular motions (Eq. 3) would induce parallel
motions when $\delta\dot{B}_{\parallel,\ell}\neq 0$, even if
$u_{\parallel,\ell}(0)=0$. However, because we have neglected the coupling of
motions, we initiated $u_{\parallel,\ell}$ from the initial conditions.

From Eq. 6, we obtain that the free streaming of the perturbed volume causes
$\ell_{\parallel}$ to expand or contract as

$\ell_{\parallel}(t)\approx\ell_{\parallel,0}\left(1\pm\frac{t}{\mathcal{T}_{\parallel}}\right).$
(27)

As the target fluid parcel expands, its environment along the $\@vec{B}_{0}$
axis contracts, provided that the fluid has fixed boundaries. Due to the
expansion of the target volume, the initial velocity of the environment would
be $-u_{\parallel,\max}$, which causes a negative sign in the denominator of
Eq. 26, and hence $\delta B_{\perp,\ell}$ increases in the environment. On the
other hand, $\delta B_{\perp,\ell}$ in the target volume stops increasing when
$t=\mathcal{T}_{\parallel}$ because $\delta B_{\perp,\ell}$ in the environment
becomes infinite. In sub-Alfvénic flows,
$|\@vec{B}_{0}|\gg|\delta\@vec{B}_{\perp}|$, so that this infinity should be
treated as an asymptotic behaviour of $\delta B_{\perp,\ell}$: there is a
physical limit above which $\delta B_{\perp,\ell}$ cannot grow. After
$\mathcal{T}_{\parallel}$, the motion is reversed and the environment starts
expanding along $\@vec{B}_{0}$, causing the target volume to contract with
$\delta B_{\perp,\ell}$ growing as 222This solution is obtained by considering
that the initial conditions in the reversed motion of the fluid parcel are
$\delta B_{\perp,\ell}(0)=\delta B_{\perp,\max}/2$,
$u_{\parallel,\ell}(0)=-u_{\parallel,\max}$, and
$\ell_{\parallel}(0)=2\ell_{\parallel,0}$. These values correspond to the
solutions of Eqs. 26 and 27 for $t=T_{\parallel}$. $\delta
B_{\perp,\ell}(t)\approx fB_{0}/(2-fV_{A}\ell_{\parallel,0}^{-1}t)$ until it
reaches $\delta B_{\perp,\max}$.