# Thermal damping of Weak Magnetosonic Turbulence in the Interstellar Medium

Kedron Silsbee1, Alexei V. Ivlev1, Munan Gong1 [ [ 1Max-Planck-Institut für
Extraterrestrische Physik, 85748 Garching, Germany

###### Abstract

We present a generic mechanism for the thermal damping of compressive waves in
the interstellar medium (ISM), occurring due to radiative cooling. We solve
for the dispersion relation of magnetosonic waves in a two-fluid (ion-neutral)
system in which density- and temperature-dependent heating and cooling
mechanisms are present. We use this dispersion relation, in addition to an
analytic approximation for the nonlinear turbulent cascade, to model
dissipation of weak magnetosonic turbulence. We show that in some ISM
conditions, the cutoff wavelength for magnetosonic turbulence becomes tens to
hundreds of times larger when the thermal damping is added to the regular ion-
neutral damping. We also run numerical simulations which confirm that this
effect has a dramatic impact on cascade of compressive wave modes.

e-mail<EMAIL_ADDRESS>e-mail<EMAIL_ADDRESS>
##  I. Introduction

Turbulence is ubiquitous in the interstellar medium (ISM), and is extended
over an extremely wide range of scales (Chepurnov & Lazarian 2010; Armstrong
et al. 1995; Minter & Spangler 1996; Marchal et al. 2021). Turbulence on large
scales plays an essential role in the theory of star formation (McKee &
Ostriker 2007). The origin of this turbulence is debated (Mac Low & Klessen
2004; Glazebrook 2013; Krumholz & Burkhart 2016). Two main classes of driving
mechanisms are gravitational instabilities and supernova blast waves. Both
input energy at scales of parsecs or larger. The turbulent cascade then brings
this energy to small scales, generating density and velocity structure which
profoundly impacts the star formation process.

The ISM is composed of gas with a wide range of ionization degrees, ranging
from almost completely ionized near massive stars (Strömgren 1939) to lower
than $10^{-8}$ in dense molecular cloud cores (Caselli et al. 2002). It is
currently thought that MHD turbulence in partially ionized gas in the ISM is
damped primarily by ion-neutral friction (Kulsrud & Pearce 1969; Balsara
1996). Xu et al. (2016) provide an analytic calculation of the scales at which
different wave modes in the gas are damped. Li & Houde (2008) study the line
widths of both charged and neutral species in a molecular cloud, showing that
the ion widths are systematically narrower. They suggest that this occurs
because the ion motions are damped more strongly than the neutral motions at
the scale of ion-neutral decoupling. Interestingly though, Pineda et al.
(2021) find in a different cloud that the ions have a systematically higher
velocity dispersion than the neutrals, indicating the need for a more detailed
analysis of the turbulent cascade at small scales.

Understanding the turbulent cascade at small scales is important for a number
of astrophysical applications, particularly in molecular clouds. Turbulence at
scales of an AU or less affects the transport of the sub-GeV cosmic rays (CRs)
which determine the heating and ionization of molecular gas (Yan & Lazarian
2004; Xu & Yan 2013; Xu et al. 2016; Silsbee & Ivlev 2019). Turbulent eddies
on small scales (Yan et al. 2004; Ormel & Cuzzi 2007) are also likely the
dominant source of relative motion between dust grains in many environments,
thus determining the grain coagulation (Ossenkopf 1993; Hirashita & Yan 2009;
Gong et al. 2020; Silsbee et al. 2020).

In this paper we show that compressive (magnetosonic) modes in weak MHD
turbulence are strongly damped due to radiative energy losses. This thermal
damping can occur at much larger scales than those on which ion-neutral
friction becomes important. The damping mechanism is as follows. The steady-
state gas temperature is determined by a balance of global heating and
cooling. As a compressive disturbance propagates, the gas undergoes adiabatic
heating and cooling in response to the pressure perturbation, while the global
processes bring it back to an equilibrium state on some characteristic
timescale, which we call the cooling time $\tau_{\rm c}$. If $\tau_{\rm c}$ is
large compared to $\omega^{-1}$, then an adiabatic equation of state is
appropriate and thermal damping is negligible. If $\tau_{\rm c}$ is small,
then the wave speed may be modified (by the different equilibrium temperatures
in the compressed and rarefied parts of the wave), but the wave amplitude also
remains constant. However, if the cooling time is comparable to $\omega^{-1}$,
then the restoring pressure force (which governs the wave propagation) is
reduced, thus damping the wave. This damping effect has been studied
previously in the context of waves in the solar atmosphere (e.g. Souffrin
1972; Mihalas & Mihalas 1983; Bunte & Bogdan 1994). Tilley & Balsara (2011)
derive a very general dispersion relation for two-fluid MHD, considering both
cooling and self-gravity. They point out that the cooling can increase the
damping of certain wave modes in molecular clouds, but do not explore under
which circumstances this will substantially affect the turbulent cascade. In
the present paper, we focus on the effect of such damping on the turbulent
cascade in the ISM and explore the regimes under which its effect is
significant.

In Section II, we derive the dispersion relation for magnetosonic waves in a
two-fluid ion-neutral plasma in the presence of thermal damping. We plot the
dispersion relation for realistic astrophysical environments in Section III,
and in Section IV discuss in detail the conditions required in order for the
effect to be large. In Section V, we present a simplified model for the
turbulent cascade, and show that under reasonable assumptions, thermal damping
may play a dominant role in damping turbulence. In Section VI, we describe the
results of numerical simulations illustrating this effect. Several potential
implications of thermal damping are outlined in Section VII, and our
conclusions are presented in Section VIII.

##  II. Dispersion relation

In this section we derive the dispersion relation for compressive MHD waves
propagating in a medium in which there are gas heating and cooling mechanisms.
Assuming small perturbations of density, $\rho=\rho_{0}+\delta\rho$, and
temperature, $T=T_{0}+\delta T$, we can write a linearized energy conservation
equation (see, e.g., Landau & Lifshitz 1960), where the thermal and magnetic
parts of the energy remain decoupled. The resulting equation for the thermal
energy balance reads:

$C_{v}\left(\rho_{0}\frac{\partial\delta T}{\partial
t}+T_{0}\frac{\partial\delta\rho}{\partial
t}\right)+C_{p}\rho_{0}T_{0}\nabla\cdot{\bf v}=\delta q.$ (1)

Here, $\bf v$ is the fluid velocity, $t$ is time,
$\rho\approx\rho_{i}+\rho_{n}$ is the total gas density, $C_{v}$ and $C_{p}$
are the gas heat capacities (per unit mass) at constant volume and pressure,
and $q(\rho,T)\equiv\Lambda_{\rm h}(\rho,T)-\Lambda_{\rm c}(\rho,T)$ is the
net energy deposition (per unit volume) from all heating and cooling
processes. In addition, we have the linearized continuity equation:

$\frac{\partial\delta\rho}{\partial t}+\rho_{0}\nabla\cdot{\bf v}=0.$ (2)

Equations (1) and (2) yield the relation

$\frac{\partial\delta T}{\partial
t}-\left(\gamma-1\right)\frac{T_{0}}{\rho_{0}}\frac{\partial\delta\rho}{\partial
t}=\frac{q^{\prime}_{\rho}\delta\rho+q^{\prime}_{T}\delta
T}{\rho_{0}C_{v}}\,,$ (3)

where $\gamma\equiv C_{p}/C_{v}$, $q^{\prime}_{\rho}\equiv\partial
q/\partial\rho$ and $q^{\prime}_{T}\equiv\partial q/\partial T$. Now, let us
assume harmonic perturbations proportional to $\exp{(-i\omega t)}$. Then,
defining

$\tau_{\rm c}=\frac{C_{v}\rho_{0}}{-q^{\prime}_{T}}\quad{\rm
and}\quad\Phi=\frac{\rho_{0}q^{\prime}_{\rho}}{T_{0}q^{\prime}_{T}}\,,$ (4)

we find from Equation (3),

$\left(1-\frac{1}{i\omega\tau_{\rm c}}\right)\frac{\delta
T}{T_{0}}=\left(\gamma-1+\frac{\Phi}{i\omega\tau_{\rm
c}}\right)\frac{\delta\rho}{\rho_{0}}\,.$ (5)

We note that $\tau_{\rm c}$ is always positive because $q^{\prime}_{T}$ must
be negative to ensure thermal stability. We find in Section III.2 that the
thermal damping effect is only relevant for $k$ such that ions and neutrals
experience many collisions per $\omega^{-1}$. This allows us to consider the
plasma as one fluid, assuming the ions and neutrals to be perfectly coupled,
both thermally and mechanically.

In what follows, we introduce the effect of thermal damping into the standard
framework of the two-fluid approach. For this, we consider the separate motion
of ions and neutrals, still assuming that ions and neutrals are perfectly
thermally coupled. Following Equations 1–3 in Soler et al. (2013), we then
write a linearized equation of motion for ion-electron fluid:

$\rho_{i}\frac{\partial{\bf v}_{i}}{\partial t}=-\nabla\delta
p_{i}+\frac{1}{4\pi}\left(\nabla\times\delta{\bf B}\right)\times{\bf
B}_{0}-\xi\left({\bf v}_{i}-{\bf v}_{n}\right).$ (6)

Here ${\bf B}_{0}$ and $\delta{\bf B}$ are the mean magnetic field and its
perturbation, ${\bf v}_{i}$ and ${\bf v}_{n}$ are the ion-electron and neutral
velocities, respectively, and $\delta p_{i}$ is the pressure perturbation of
the ion-electron fluid. $\delta{\bf B}$ obeys the linearized equation

$\frac{\partial\delta{\bf B}}{\partial t}=\nabla\times({\bf v}_{i}\times{\bf
B}_{0}).$ (7)

$\xi$ is a constant describing the magnitude of the ion-neutral friction,
given by

$\xi=\frac{\rho_{i}\rho_{n}\langle\sigma v\rangle_{in}}{m_{i}+m_{n}}\,,$ (8)

where $\langle\sigma v\rangle_{in}=4\pi
a_{0}^{2}\sqrt{(\alpha/a_{0}^{3})I_{\rm H}/\mu}$ is expressed via the Bohr
radius $a_{0}$, the Rydberg energy $I_{\rm H}$, and
$\mu=m_{i}m_{n}/(m_{i}+m_{n})$. The dimensionless polarizability of the
neutral species, $\alpha/a_{0}^{3}$, is 4.5 for atomic hydrogen and 5.52 for
molecular hydrogen (Raizer et al. 2011).

Similarly, the equation of motion for the neutrals reads

$\rho_{n}\frac{\partial{\bf v}_{n}}{\partial t}=-\nabla\delta
p_{n}-\xi\left({\bf v}_{n}-{\bf v}_{i}\right).$ (9)

By writing $\delta p_{i,n}=n_{i,n}k_{\rm B}\delta T+k_{\rm B}T_{0}\delta
n_{i,n}$, where $k_{\rm B}$ is the Boltzmann constant, $n_{i}\approx
2\rho_{i}/m_{i}$ and $n_{n}=\rho_{n}/m_{n}$, we obtain from Equation (5) for
harmonic perturbations

$\delta p_{i,n}=c_{i,n}^{2}f(\omega)\delta\rho_{i,n}\,,$ (10)

where

$c_{i}=\sqrt{\frac{2\gamma k_{\rm B}T_{0}}{m_{i}}}\,,$ (11)
$c_{n}=\sqrt{\frac{\gamma k_{\rm B}T_{0}}{m_{n}}}\,,$ (12)

are the sound speeds of ions and neutrals, respectively,

$f(\omega)=\frac{i\omega\tau_{\rm c}-\chi}{i\omega\tau_{\rm c}-1}\,,$ (13)

and

$\chi=\frac{1-\Phi}{\gamma}\,.$ (14)

We note that $\chi>0$ is the requirement for a stable equilibrium to exist. As
shown in Appendix A, the imaginary part of the frequency in the long-
wavelength regime becomes positive when $\chi>1$. This would imply a thermal
instability. In principle, this could be realized if the concentration of a
dominant coolant were strongly suppressed at higher density. We do not
consider such a hypothetical situation in this paper, and thus assume
$0<\chi<1$.

Combining Equation (10) with Equation (2), we obtain

$i\omega\delta p_{i,n}=c_{i,n}^{2}\rho_{i,n}f(\omega)\nabla\cdot{\bf
v}_{i,n}\,.$ (15)

Note that we are are interested in magnetosonic perturbations, as non-
compressive Alfvénic perturbations are not radiatively damped. Hence,
following Soler et al. (2013), we take the time derivative of the divergence
of Equation (6). Using Equations (7) and (15) and writing
$\Delta_{i,n}=\nabla\cdot{\bf v}_{i,n}$ for brevity, we find

$\displaystyle-\omega^{2}\rho_{i}\Delta_{i}=c_{i}^{2}\rho_{i}f(\omega)\nabla^{2}\Delta_{i}$
$\displaystyle+\frac{1}{4\pi}\nabla\cdot\left\\{\left[\nabla\times\nabla\times({\bf
v}_{i}\times{\bf B}_{0})\right]\times{\bf
B}_{0}\right\\}+i\omega\xi\left(\Delta_{i}-\Delta_{n}\right).$ (16)

For ${\bf B}_{0}=B_{0}\hat{\bf z}$, some manipulation of the middle term in
Equation (16) allows us to write

$\displaystyle-\omega^{2}\rho_{i}\Delta_{i}=c_{i}^{2}\rho_{i}f(\omega)\nabla^{2}\Delta_{i}$
$\displaystyle+\frac{B_{0}^{2}}{4\pi}\nabla^{2}\left(\Delta_{i}-\Delta_{i}^{\parallel}\right)+i\omega\xi\left(\Delta_{i}-\Delta_{n}\right),$
(17)

where $\Delta_{i}^{\parallel}\equiv dv_{i}^{z}/dz$ is the parallel component
of the divergence. Taking the time derivative of the parallel divergence of
Equation (6) yields

$-\omega^{2}\rho_{i}\Delta_{i}^{\parallel}=c_{i}^{2}\rho_{i}f(\omega)\frac{\partial^{2}\Delta_{i}}{\partial
z^{2}}+i\omega\xi\left(\Delta_{i}^{\parallel}-\Delta_{n}^{\parallel}\right).$
(18)

A complementary pair of equations for $\Delta_{n}$ and
$\Delta_{n}^{\parallel}$ are derived from Equations (9) and (15):

$-\omega^{2}\rho_{n}\Delta_{n}=\rho_{n}c_{n}^{2}f(\omega)\nabla^{2}\Delta_{n}+i\omega\xi(\Delta_{n}-\Delta_{i}),$
(19)
$-\omega^{2}\rho_{n}\Delta_{n}^{\parallel}=c_{n}^{2}\rho_{n}f(\omega)\frac{\partial^{2}\Delta_{n}}{\partial
z^{2}}+i\omega\xi(\Delta_{n}^{\parallel}-\Delta_{i}^{\parallel}).$ (20)

Defining the Alfvén speed,

$c_{\rm A}=\frac{B_{0}}{\sqrt{4\pi\rho_{i}}}\,,$ (21)

the ion-neutral and neutral-ion momentum transfer frequencies,

$\nu_{in}=\frac{\xi}{\rho_{i}}\quad{\rm
and}\quad\nu_{ni}=\frac{\xi}{\rho_{n}}\,,$ (22)

and assuming all perturbed quantities to be proportional to $\exp(-i\omega
t+i{\bf k}\cdot{\bf r})$, we can write the following 4 coupled equations

$\left[\omega^{2}-k^{2}c_{i}^{2}f(\omega)\right]\Delta_{i}+i\omega\nu_{in}(\Delta_{i}-\Delta_{n})=k^{2}c_{\rm
A}^{2}\left(\Delta_{i}-\Delta_{i}^{\parallel}\right),$ (23)
$\left[\omega^{2}-k^{2}c_{n}^{2}f(\omega)\right]\Delta_{n}+i\omega\nu_{ni}(\Delta_{n}-\Delta_{i})=0,$
(24)
$\omega^{2}\Delta_{i,n}^{\parallel}-k_{z}^{2}c_{i,n}^{2}f(\omega)\Delta_{i,n}\pm
i\omega\nu_{in,ni}\left(\Delta_{i}^{\parallel}-\Delta_{n}^{\parallel}\right)=0.$
(25)

Equations (23) - (25) yield a 9th order polynomial dispersion relation
describing the modes present in a partially ionized medium with thermal
damping.

##  III. Magnetosonic modes in astrophysical environments

Figure 1.— Illustration of the dispersion relation of magnetosonic modes
derived from Equations (23)–(25) for two phases of ISM. The solid curves show
the real and imaginary parts of the dispersion relation plotted as a function
of the wavelength $\lambda$. The dashed curves show the result in the limit
$\tau_{\rm c}\rightarrow\infty$ studied in Soler et al. (2013). Here, we have
assumed an angle $\theta$ of 45∘ between ${\bf k}$ and ${\bf B}$, and that the
loaded Alfvén speed $c_{\rm LA}=1.4c_{n}$ (corresponding to plasma
$\beta=0.6$). The shaded gray region indicates the decoupling gap – the
wavelength range over which the decoupling of ions and neutrals occurs. The
blue and green vertical lines show the wavelengths where Re $\omega=\tau_{\rm
c}^{-1}$, near where the thermal damping peaks for the corresponding mode. The
decoupled, loaded (L), and modified loaded (M) regimes are indicated above the
top left panel (see Section III.2 for a detailed discussion).

In this section, we numerically calculate the dispersion relation on the basis
of Equations (23)–(25), with and without thermal damping, and analyze
important features introduced by the damping. The results are illustrated for
characteristic astrophysical environments. In this paper, we have introduced a
number of parameters, including in particular speeds of wave modes for three
different propagation regimes discussed in this section (see Figure 1). The
most important of these parameters are listed in Table 1.

Table 1Notation for the key parameters Parameter | Meaning | Equation
---|---|---
$\tau_{\rm c}$ | gas cooling time | (4)
$\nu_{in}$ | ion-neutral collision frequency | (22)
$\nu_{ni}$ | neutral-ion collision frequency | (22)
$c_{\rm A}$ | Alfvén speed | (21)
$c_{\rm LA}$ | loaded Alfvén speed | (26)
$c_{i}$ | ion sound speed | (11)
$c_{n}$ | neutral sound speed | (12)
$c_{{\rm M}n}$ | modified sound speed | (29)
$c_{\rm f,s}$ | fast/slow speed | (28)
$c_{\rm Lf,s}$ | loaded fast/slow speed | (28) with
|  |  $c_{\rm A}\rightarrow c_{\rm LA}$,
|  | $c_{i}\rightarrow c_{n}$
$c_{\rm Mf,s}$ | modified loaded fast/slow speed | (28) with
|  |  $c_{\rm A}\rightarrow c_{\rm LA}$,
|  | $c_{i}\rightarrow c_{{\rm M}n}$
$\chi$ | ratio $(C_{{\rm M}n}/c_{n})^{2}$ | (14) & (29)
$\gamma$ | adiabatic index | $C_{p}/C_{v}$
$\beta$ | thermal to magnetic pressure ratio | (27)
$\beta_{\rm M}$ |  $\beta$ in the modified loaded regime | (31)
$\mathcal{M}_{n}$ | sonic Mach number | $v_{0}/c_{n}$
$\mathcal{M}_{\rm LA}$ | (loaded) Alfvénic Mach number | $v_{0}/c_{\rm LA}$

###  III.1. Astrophysical Environments

In our analysis, we study the effect in detail in two phases of the ISM: the
cold neutral medium (CNM), and the envelope of a molecular cloud. We discuss
the expected damping in other environments in Section IV.

####  III.1.1 CNM

The CNM is made up of atomic hydrogen, for which we assume $n_{n}=20$ cm-3.
The ion is C+, and the ionization fraction is $1.6\times 10^{-4}$. The
presence of helium or other elements is ignored. Heating is dominated by the
photoelectric effect on dust, and given by $\Lambda_{\rm h}=1.4\times
10^{-26}(n_{n}/{\rm cm}^{-3})$ erg cm-3 s-1 (Draine 2011). The cooling is from
fine-structure lines of C+, and we use the optically thin cooling prescription
considering collisions with H atoms and electrons from Gong et al. (2017).
This gives $T_{0}=60$ K, $\tau_{\rm c}=6\times 10^{11}$ s, and $\chi=0.22$.

####  III.1.2 Cloud Envelope

The neutral constituent is assumed to be molecular hydrogen, the ion is C+,
and the ionization fraction is 0.032% (implying, as in the CNM, a C/H ratio of
0.016%). Goldsmith (2001) suggests that $\Lambda_{\rm c}=4.4\times
10^{-27}(T/{\rm K})^{2.4}$ erg cm-3 s-1 for $n_{n}=10^{3}$ cm-3. We assume
heating to be from a combination of CR and photoelectric heating, with a rate
given by $\Lambda_{\rm h}=1.0\times 10^{-26}(n_{n}/{\rm cm}^{-3})$ erg cm-3
s-1. Setting $n_{n}=10^{3}$ cm-3, we obtain $T_{0}=25$ K. At this temperature,
the logarithmic slope of the cooling function with density, estimated from the
data in Table 2 from Goldsmith (2001), is approximately 1.1. Applying
Equations (4) and (14), we find $\tau_{\rm c}=2.2\times 10^{11}$ s, and
$\chi=0.58$.

###  III.2. Effect of Thermal Damping

Figure 1 illustrates the dispersion relation for magnetosonic waves, showing a
comparison of the cases with and without the thermal damping effect. For all
panels, we set an angle $\theta$ between ${\bf k}$ and ${\bf B}_{0}$ of 45∘
and assume that what we call the loaded Alfvén speed, given by

$c_{\rm LA}=\frac{B_{0}}{\sqrt{4\pi\rho}}\equiv c_{\rm
A}\sqrt{\frac{\rho_{i}}{\rho}}\,.$ (26)

is equal to 1.4 times the neutral sound speed. This is the speed of Alfvén
waves in the regime where collisions are so frequent that the neutrals are
dragged along by the ions, i.e., everything to the right of the gray region in
Figure 1. The coefficient 1.4 was derived from the respective magnetic and
thermal energies indicated in Table 1.5 of Draine (2011). This corresponds to
a plasma $\beta$ of 0.6, where $\beta$ is the ratio of thermal to magnetic
pressure:

$\beta=\frac{8\pi nk_{\rm B}T_{0}}{B_{0}^{2}}\equiv\frac{2c_{n}^{2}}{\gamma
c_{\rm LA}^{2}}\,.$ (27)

The dispersion relation is plotted as a function of wavelength
$\lambda=2\pi/k$. We have only plotted propagating modes (those with a non-
zero real part). In each case, there is a corresponding mode (not shown) with
a negative real part of equal magnitude.

The solid lines in Figure 1 show the dispersion relation of magnetosonic modes
derived from Equations (23)–(25). The dashed lines depict the results in the
adiabatic limit $\tau_{\rm c}\rightarrow\infty$ [where $f(\omega)\rightarrow
1$ in Equation (10)]. In all illustrated cases, there is one mode, shown in
green, which is continuously present at all $\lambda$. For the chosen
parameters, this is the loaded fast mode in the long-wavelength regime, which
becomes the neutral acoustic mode in the short-wavelength (decoupled) regime
(see Section III.4 and Appendix B for general analysis). Also plotted are the
loaded slow mode (blue), the decoupled fast and slow modes (red and orange),
as well as a low-frequency mode (violet), present at intermediate wavelengths.

The speeds of the decoupled fast and slow modes, $c_{\rm f,s}$, can directly
obtained from Equations (23) and (25) by neglecting the frictional coupling
terms. This yields the classical expression for magnetosonic modes (Landau &
Lifshitz 1960)

$c_{\rm f,s}^{2}=\frac{1}{2}\left[c_{\rm A}^{2}+c_{i}^{2}\pm\sqrt{(c_{\rm
A}^{2}+c_{i}^{2})^{2}-4c_{\rm A}^{2}c_{i}^{2}\cos^{2}{\theta}}\>\right].$ (28)

In the loaded regime, where ions and neutrals are moving together, the speeds
$c_{\rm Lf,s}$ of loaded fast and slow modes are given by Equation (28) where
$c_{\rm A}$ is replaced by $c_{\rm LA}$, and $c_{i}$ is replaced by $c_{n}$.
The range of wavelengths over which ion-neutral decoupling occurs – referred
below as the decoupling gap – is indicated in each panel of Figure 1 by the
shaded gray region. As shown in Appendix B, the right edge of the decoupling
gap is approximately described by the condition $kc_{\rm LA}\approx
2\nu_{ni}\cos\theta$ for $c_{\rm LA}\ll c_{n}$ (where $c_{\rm Lf}\approx
c_{n}$ and $c_{\rm Ls}\approx c_{\rm LA}\cos\theta$), and by $kc_{\rm
LA}\approx 2\nu_{ni}$ in the opposite limit (where $c_{\rm Lf}\approx c_{\rm
LA}$ and $c_{\rm Ls}\approx c_{n}\cos{\theta}$). In Figure 1, the edges of the
decoupling gap were determined numerically. In Appendix B we present a
detailed analysis of wave modes in the loaded regime and their connection to
the decoupling gap.

The blue and green vertical lines on the right side from the decoupling gap
correspond to the condition Re $\omega=\tau_{\rm c}^{-1}$, where $\omega$ is
the frequency represented by the curve of the same color. These lines
approximately match the peak of the damping associated with radiative cooling,
where the real parts of the loaded fast and slow modes change to their
modified loaded values. These are given by Equation (28), with $c_{\rm A}$
replaced by $c_{\rm LA}$ (as in the loaded regime), while $c_{i}$ is now
replaced by

$c_{{\rm M}n}=c_{n}\sqrt{\chi}\,,$ (29)

as follows from Equations (10) and (13) in the low-frequency limit.

We see that thermal damping is only important in the modified loaded regime,
at $\lambda$ much larger than the scale of ion-neutral decoupling. This
justifies our initial assumption that ions and neutrals are thermally coupled.
For the chosen parameters, thermal damping has a rather moderate effect on the
real part of both modes, but a dramatic effect – ranging from a factor of
several in the CNM to several hundred in the case of the cloud envelope – on
the imaginary part.

###  III.3. Thermal Damping versus Propagation Angle

Figure 1 suggests that thermal damping can increase the imaginary part of wave
modes in the modified loaded regime by over two orders of magnitude. For that
reason, in this section we focus on the damping rate and explore how that
depends on propagation angle. In Appendix A we show that the imaginary part of
the modified loaded fast and slow modes is approximately given by

$-{\rm Im}~{}\omega=\frac{\chi^{-1}-1}{4}\,\tau_{\rm c}c_{{\rm
M}n}^{2}\mathscr{F}_{\rm f,s}(\beta_{\rm M},\theta)k^{2},$ (30)

where $\mathscr{F}_{\rm f,s}$ is given by Equation (A7) and

$\beta_{\rm M}=\frac{c_{{\rm M}n}^{2}}{c_{\rm
LA}^{2}}\equiv\frac{\chi\gamma}{2}\beta,$ (31)

is the plasma $\beta$ in the modified regime.

The strength of thermal damping is related to the degree to which the mode
results in compression of the gas. Therefore, we expect the imaginary part of
the dispersion relation in the modified loaded regime to be a function of
$({\bf k}\cdot{\bf v})^{2}$, where ${\bf v}$ is the velocity eigenvector. The
components $v_{\perp}$ and $v_{\parallel}$ of the eigenvector, perpendicular
and parallel to ${\bf k}$, are related via (e.g., Landau & Lifshitz 1960)

$\frac{v_{\perp}}{v_{\parallel}}=\frac{\cos{2\theta}-\beta_{\rm
M}\pm\sqrt{1+\beta_{\rm M}^{2}-2\beta_{\rm
M}\cos{2\theta}}}{\sin{2\theta}}\,,$ (32)

for the fast ($+$) and slow ($-$) mode. Using Equations (32) and (A7), we find
that

$\mathscr{F}=2(\hat{\bf k}\cdot\hat{\bf v})^{2}.$ (33)

This dependence allows us to easily understand the behavior of
$\mathscr{F}_{\rm f,s}$, plotted in Figure 2 as a function of $\theta$ for
different values of $\beta_{\rm M}$. For large $\beta_{\rm M}$ (i.e., weak
magnetization) the fast mode is approximately the (modified) neutral sound
mode, with ${\bf v}$ almost parallel to ${\bf k}$, while ${\bf v}$ of the slow
mode (with $c_{\rm Ms}\approx c_{\rm LA}\cos\theta$) is almost perpendicular
to ${\bf k}$, as follows from Equation (32). Therefore, the damping of the
fast mode is strong and almost independent of $\theta$ ($\mathscr{F}_{\rm
f}\approx 2$), while the damping of the slow mode is weak. On the contrary,
for small $\beta_{\rm M}$ (strong magnetization) we have $c_{\rm Mf}\approx
c_{\rm LA}$ with ${\bf v}$ almost perpendicular to ${\bf B}_{0}$, and $c_{\rm
Ms}\approx c_{{\rm M}n}\cos\theta$ with ${\bf v}$ parallel to ${\bf B}_{0}$.
In this case, we have $\mathscr{F}_{\rm f}\approx 2\sin^{2}\theta$ and
$\mathscr{F}_{\rm s}\approx 2\cos^{2}\theta$. We note that the damping of the
fast and slow modes is symmetrically opposite, generally following
$\mathscr{F}_{\rm f}+\mathscr{F}_{\rm s}=2$.

Figure 2.— Prefactor $\mathscr{F}_{\rm f,s}$ of the thermal damping rate,
given by Equation (30). Curves are potted as a function of $\theta$ for
different values of $\beta_{\rm M}=\chi\gamma\beta/2$. The solid lines
correspond to the modified fast loaded mode $\mathscr{F}_{\rm f}$, and the
dashed lines to the modified slow loaded mode, $\mathscr{F}_{\rm
s}=2-\mathscr{F}_{\rm f}$.

###  III.4. Thermal Damping of the Continuous Mode

Let us determine a condition for the existence of the continuous mode. In all
environments considered in this paper, $\rho_{i}\ll\rho_{n}$. In this limit,
if there is a continuous mode, it is the neutral sound mode in the decoupled
regime and throughout the decoupling gap. In Appendix B we derive that,
depending on parameters, it can switch either to the loaded fast mode (as
illustrated in Figure 1) or to the loaded slow mode at the right edge of the
decoupling gap. Figure 3 shows that the continuity is determined by two
parameters, $c_{n}/c_{\rm LA}$ and $\theta$. For $\theta\lesssim 64^{\circ}$,
the neutral sound branch switches between the loaded fast and slow branches at
the solid line, described by Equation (B8). However, the neutral sound mode
disappears completely within the decoupling gap for sufficiently small
$c_{n}/c_{\rm LA}$ and large $\theta$, i.e., for the parameters bound in the
right bottom corner of Figure 3 by the dashed lines, no mode is continuous for
all $k$.

Figure 3.— The parameter space of $c_{n}/c_{\rm LA}\equiv\sqrt{\beta_{\rm
M}/\chi}$ and $\theta$, showing where either the fast mode, the slow mode, or
neither mode are continuous with the neutral sound mode in the decoupled
regime (see the left panel of Figure 1). The solid line is the analytical
condition $c_{n}/c_{\rm LA}=\cos\theta$, Equation (B8), the dotted lines are
computed numerically.

Now we can evaluate the damping rate of the continuous mode in the modified
loaded regime. From Figures 2 and 3 we conclude that for large $c_{n}/c_{\rm
LA}\equiv\sqrt{\beta_{\rm M}/\chi}$, the continuous mode in the modified
loaded regime is represented by the fast mode, whose damping is strong and
almost independent of $\theta$. For small $c_{n}/c_{\rm LA}$, the continuous
mode (for $\theta$ where it exists) is the slow mode, which is strongly damped
too (although the damping decreases with $\theta$). In the range of
$0.45\lesssim c_{n}/c_{\rm LA}<1$, the continuous mode switches from the slow
to the fast mode as $\theta$ increases. The damping is then described by the
dashed green or blue lines in Figure 2 at smaller $\theta$, switching to the
respective solid lines at larger $\theta$.

Something qualitatively different occurs when $1<c_{n}/c_{\rm
LA}<1/\sqrt{\chi}$. In this range, the continuous mode is the fast mode, but
we see in Figure 2 (where the corresponding range is $\chi<\beta_{\rm M}<1$)
that the fast mode is then weakly damped at small $\theta$. This implies that
the resulting turbulent spectrum is expected to be strongly anisotropic for
this range of $c_{n}/c_{\rm LA}$.

We conclude that thermal damping of the continuous mode – as long as it exists
– is generally strong. The damping rate can vary with $\theta$, depending on
the value of $c_{n}/c_{\rm LA}$, but this variation is rather insignificant
outside the relatively narrow range of $1<c_{n}/c_{\rm LA}<1/\sqrt{\chi}$. In
the next section, where we analyze the impact of thermal damping on the
turbulent cascade, the damping rate is therefore assumed to be independent of
$\theta$.

##  IV. Characteristics of Thermal Damping in Different Environments

In assessing the importance of thermal damping, we consider two criteria.
First is the peak value of $\tilde{\omega}_{I}\equiv{\rm|Im\,\omega|}/{\rm
Re\,\omega}$ due to thermal damping, providing a measure of the overall
strength of the effect. Second, we consider the length scale at which the
effect is strongest, and compare this with the scale corresponding to the
right edge of the decoupling gap. Magnetosonic modes are generally affected by
both thermal damping and ion neutral friction, so it is of interest to ask
which mechanism will damp the turbulent cascade at larger scales. To address
this question, in this section we also derive the ratio of thermal damping
rate to damping rate due to ion-neutral friction in the long-wavelength limit.

Equation (A1) describes the dispersion relation in the absence of ion-neutral
friction. We numerically solved Equation (A1) to determine the peak value of
$\tilde{\omega}_{I}$ as a function of $\beta$ and $\chi$. This is shown for
both fast and slow modes in Figure 4. We note that the peak value of
$\tilde{\omega}_{I}$ does not depend on $\tau_{\rm c}$. To see this, define
$\omega(k,\tau_{\rm c})$ as a solution of Equation (A1) for given parameters
$\beta$ and $\chi$. We note that $\omega(k/x,x\tau_{\rm
c})=x^{-1}\omega(k,\tau_{\rm c}$) for an arbitrary constant $x$. This means
that the ratio $\tilde{\omega}_{I}$ satisfies
$\tilde{\omega}_{I}(k/x,x\tau_{\rm c})=\tilde{\omega}_{I}(k,\tau_{\rm c})$.
Therefore, by changing $\tau_{\rm c}$, the location of the peak value of
$\tilde{\omega}_{I}$ is shifted to different $k$ such that $k\tau_{\rm c}$
remains constant, but the peak value of $\tilde{\omega}_{I}$ is unchanged.

As discussed in Section III.3, $\tilde{\omega}_{I}$ depends on the mode
compressibility, and this behavior is reflected in Figure 4. Modes which are
mostly compressive have higher peak values of $\tilde{\omega}_{I}$. These are
the fast modes in the limit of high $\beta$ and the slow modes in the limit of
low $\beta$ (see Figure 2). Figure 4 also shows that the effect is maximized
for small values of $\chi$. This is consistent with Equation (30), which shows
that the damping rate of fast and slow modes is proportional to the deviation
of $\chi^{-1}$ from unity.

The question of the wavelength $\lambda_{\rm peak}$ at which damping peaks
relative to the right edge of the decoupling gap, $\lambda_{\rm dec}$, is more
difficult to address in general. The physics which determines the cooling is
complicated and varies throughout the different ISM phases. In addition to the
two environments analyzed above in detail (the CNM and cloud envelope), in
this section we also briefly discuss the degree of thermal damping in the warm
neutral medium (WNM), and the inner and outer parts of a prestellar core.

Figure 4.— Peak value of $\tilde{\omega}_{I}={\rm|Im\,\omega|}/{\rm
Re\,\omega}$ due to thermal damping, plotted in the plane of $\beta$ and
$\chi$. As expected from the considerations discussed in Sections III.3 and
III.4, the peak of $\tilde{\omega}_{I}$ is higher for the continuous mode,
which is represented by the fast mode if $\beta\gg 1$, and by the slow mode if
$\beta\ll 1$. In both cases, the effect is maximized for small $\chi$.

###  IV.1. Different additional environments

In the WNM, the ion is H+ and the neutral component is atomic hydrogen. The
heating in the WNM is dominated by the photoelectric effect on dust, and we
adopt a constant heating rate of $\Lambda_{\rm h}=1.4\times
10^{-26}(n_{n}/{\rm cm}^{-3})$ erg cm-3 s-1 from Draine (2011). The cooling
function $\Lambda_{\rm c}$ is estimated from Figure 30.1 of Draine (2011),
where an ionization fraction of 1.7% is assumed, and we adopt this value, too.
We set the gas number density $n_{n}=0.4$ cm-3, and calculate an equilibrium
gas temperature $T_{0}=8.7\times 10^{3}$ K. This yields a gas cooling time
$\tau_{\rm c}=5\times 10^{13}$ s, and $\chi=0.37$.

The right edge of the decoupling gap is approximated by $kc_{\rm LA}\approx
2\nu_{ni}$ (see the end of Appendix B). This gives the decoupling scale
$\lambda_{\rm dec}=\pi c_{\rm LA}/\nu_{ni}$, which corresponds to $2\times
10^{4}$ AU in the WNM. The peak of the thermal damping is estimated from the
criterion $\omega\approx\tau_{\rm c}^{-1}$. Since the fast mode is the one
which is continuous over all $\lambda$ for our fiducial values of $\beta=0.6$
and $\theta=45^{\circ}$, we obtain $\lambda_{\rm peak}\approx 2\pi\tau_{\rm
c}c_{\rm Lf}$, which corresponds to a scale of $4\times 10^{7}$ AU in the WNM.

The continuous mode in the WNM is also affected by collisionless damping,
because the thermal speed of ions in this environment is about $c_{n}$.
Therefore, were the medium fully ionized, the rate of collisionless damping
due to ions would be of order $\omega$ for most values of $\theta$ (Ginzburg
1970). However, since the wave energy is proportional to the mass of neutrals,
but only ions contribute to the damping, the resulting damping rate (in units
of $\omega$) is bounded by the ionization fraction, which we assume to be
1.7%.

Let us also consider different regions of a prestellar core. Adjacent to the
cloud envelope is what we call an outer core – a region where the interstellar
UV radiation is already attenuated, so that the ionization fraction is
controlled by CRs, but the gas density is still sufficiently low to ensure the
cooling dominated by gas species. A denser region, where the cooling is due to
dust, is referred to as an inner core.

For the outer core, we set $n_{n}=10^{4}$ cm-3. The heating is assumed to be
dominated by CRs, with $\Lambda_{\rm h}=1.0\times 10^{-27}\,(n_{n}/{\rm
cm}^{-3})$ erg cm-3 s-1 (i.e. an order of magnitude smaller than in Section
III.1.2 for a given value of $n_{n}$). Again, we use the cooling function
taken from Table 2 of Goldsmith (2001), appropriate for a density of $10^{4}$
cm-3, which yields $T_{0}=12.4$ K, $\tau_{c}=1.0\times 10^{12}$ s, and
$\chi=0.67$. HCO+ is assumed to be the principal ion with fractional abundance
of $10^{-7}$ relative to H2 (Williams et al. 1998). We obtain $\lambda_{\rm
peak}\approx 2\times 10^{4}$ AU and $\lambda_{\rm dec}\approx 4\times 10^{3}$
AU. As before, $\lambda_{\rm peak}$ and $\lambda_{\rm dec}$ are evaluated for
the fast mode, assuming $\beta=0.6$ and $\theta=\pi/4$.

The inner core corresponds to gas densities above $10^{5}$ cm-3, where the
equilibrium thermodynamics is determined by CR heating and cooling due to
collisions with dust grains. As discussed in Ivlev et al. (2019), the rate of
cooling generally depends on the gas density and dust size distribution. In
particular, there is a critical grain radius $A$, given by

$A=\frac{1}{\sqrt{18\pi}}\frac{n_{n}v_{n}k_{\rm B}}{q_{\rm abs}\sigma
T_{d0}^{5}}\,,$ (34)

where $v_{n}=\sqrt{k_{\rm B}T/m_{n}}$ is the scale of gas thermal velocity for
the mass $m_{n}=2m_{\rm H}$, $q_{\rm abs}$ is a material-dependent factor
taken to be 0.13 K-2 cm-1, $\sigma$ is the Stefan-Boltzmann constant, and
$T_{\rm d0}$ is the equilibrium dust temperature, set by the local radiation
field in the absence of coupling to the gas.

Table 2Characteristics of thermal dampinga
Environment $c_{\rm LA}/c_{n}=1.4~{}(\beta=0.6)$ $c_{\rm
LA}/c_{n}=0.5~{}(\beta=4.8)$

|  $T$, K |  $n_{n}$, cm-3 | $n_{i}/n_{n}$ |  $\lambda_{\rm dec}$, AU |  $\lambda_{\rm peak}$, AU | $\mathscr{R}$ |  $\lambda_{\rm dec}$, AU |  $\lambda_{\rm peak}$, AU | $\mathscr{R}$
---|---|---|---|---|---|---|---|---|---
WNM | 8700 | 0.4 | 0.017 | $2\times 10^{4}$ | $4\times 10^{7}$ | 140 | 9000 | $2\times 10^{7}$ | 2400
CNM | 60 | 20 | $1.6\times 10^{-4}$ | 3000 | $4\times 10^{4}$ | 1.2 | 1200 | $2\times 10^{4}$ | 28
Cloud envelope | 25 | 1000 | $3.2\times 10^{-4}$ | 20 | $6\times 10^{3}$ | 23 | 7 | 4000 | 300
Outer core | 12.4 | $10^{4}$ | $10^{-7}$ | $4\times 10^{3}$ | $2\times 10^{4}$ | 0.3 | 1500 | $1.3\times 10^{4}$ | 3
Inner core | 6 | $\gtrsim 10^{5}$ | Eq. (36) | $1300/n_{6}^{0.44}$ | $160\,a_{0.1}/n_{6}$ | $0.01\,a_{0.1}/n_{6}^{0.56}$ | $450/n_{6}^{0.44}$ | $100\,a_{0.1}/n_{6}$ | $0.1\,a_{0.1}/n_{6}^{0.56}$
($n_{6}/a_{0.1}\ll 2$) |  |  |  |  |  |  |  |  |
Inner core | 6 | $\gtrsim 10^{5}$ | Eq. (36) | $1300/n_{6}^{0.44}$ | 80 | $0.004\,n_{6}^{0.44}$ | $450/n_{6}^{0.44}$ | 50 | $0.06\,n_{6}^{0.44}$
($n_{6}/a_{0.1}\gg 2$) |  |  |  |  |  |  |  |  |

aWavelengths of the decoupling and damping peaks are calculated as
$\lambda_{\rm dec}=\pi c_{\rm LA}/\nu_{ni}$ and $\lambda_{\rm
peak}=2\pi\tau_{\rm c}c_{\rm Lf}$, the damping rate ratio $\mathscr{R}$ (for
the continuous mode) is given by Equation (37); $n_{6}$ is the gas density
$n_{n}$ in units of $10^{6}$ cm-3, $a_{0.1}$ is the effective grain radius
$a_{\rm eff}$ in units of 0.1 microns.

If the effective grain radius $a_{\rm eff}$ — see Equation 16 of Ivlev et al.
(2019) — is much larger than $A$, then the grain temperature is independent of
the gas temperature. This regime occurs when $n_{6}/a_{0.1}\ll 2$, where
$n_{6}$ is $n_{n}$ in units of $10^{6}$ cm-3 and $a_{0.1}$ is $a_{\rm eff}$ in
units of 0.1 microns. In this regime, we estimate the cooling from Equation 18
of Ivlev et al. (2019). Using Equations (4) and (35), assuming $T-T_{d0}\ll
T_{d0}$, and keeping in mind that $C_{v}=3k_{\rm B}/(4m_{\rm H})$ for
molecular hydrogen, this leads to a cooling time of $\tau_{\rm
c}\approx\sqrt{\pi/8}\>a_{\rm eff}\rho_{d}/(f_{d}v_{n}\rho_{n})=1.2\times
10^{10}n_{6}^{-1}a_{0.1}$ s. In this expression, $f_{d}$ is the dust to gas
mass ratio, $\rho_{n}$ the gas mass density, and $\rho_{d}$ the material
density of the dust grains. In the numerical estimate, we have taken
$f_{d}=0.02$ and $\rho_{d}=2$ g cm-3. We have chosen such values of $f_{d}$
and $\rho_{d}$ to account for the fact that the grains are likely covered with
icy mantles which increase their mass fraction and lower their mean density.
Molecular cloud cores show a wide range of $\beta$ (Crutcher et al. 2010). For
$\beta=0.6$, as in Figure 1, the corresponding peak damping scale is
$\lambda_{\rm peak}=2\pi\tau_{\rm c}c_{\rm LF}\approx 160n_{6}^{-1}a_{0.1}$
AU. Assuming $T-T_{d0}\ll T_{d0}$, Equation (14) yields
$\chi\approx\gamma^{-1}=3/5$.

In the opposite limit, where $n_{6}/a_{0.1}\gg 2$, the cooling can be obtained
from Equation (19) of Ivlev et al. (2019):

$\Lambda_{\rm c}(\rho_{n},T)=\frac{3f_{d}\rho_{n}q_{\rm
abs}\sigma}{\rho_{d}}\left(T^{6}-T_{\rm d0}^{6}\right),$ (35)

while $\Lambda_{\rm h}(\rho_{n},T)\propto\rho_{n}$. From this, we calculate
$\tau_{\rm c}=k_{\rm B}\rho_{d}/(24f_{d}m_{\rm H}q_{\rm abs}\sigma
T^{5})=6.0\times 10^{9}$ s, assuming $T_{\rm d0}=6$ K. Setting, as before,
$\beta=0.6$, this results in $\lambda_{\rm peak}=2\pi\tau_{\rm c}c_{\rm
Lf}\approx 80$ AU (independent of $n_{n}$). Equation (14) yields
$\chi\approx\gamma^{-1}=3/5$ in this limit, too.

To estimate the right edge of the decoupling gap in the inner core, we must
have an estimate for the density and ionization fraction. We use the relation
given for the ionization fraction as a function of density in Caselli et al.
(2002) (their model 3),

$n_{i}/n_{n}=2.3\times 10^{-9}n_{6}^{-0.56},$ (36)

and assume the dominant ion to be HCO+. The criterion $kc_{\rm LA}=2\nu_{ni}$
yields the scale of the right edge of the decoupling gap in terms of density:
$\lambda_{\rm dec}=1.3\times 10^{3}$ $n_{6}^{-0.44}$ AU.

We conclude that, although there are significant sources of uncertainty, the
peak of the thermal damping in the inner core occurs at scales substantially
smaller than the right edge of the decoupling gap.

###  IV.2. Thermal Damping versus Ion-neutral Friction

One quantity of particular interest in studying the turbulent cascade is the
relative magnitude of the thermal damping, as compared to the damping due to
ion-neutral friction. The thermal damping rate in the modified loaded regime
is given by Equation (30). The damping rate due to ion-neutral friction
follows from the expression derived for the loaded regime, Equation (B4): in
the modified regime, we have to replace $c_{n}$ with $c_{{\rm M}n}$, which is
equivalent to setting $\chi=1$ in the expression for $\mathscr{G}_{\rm
f,s}(\beta_{\rm M}/\chi,\theta)$, Equation (B5).

Thus, the ratio of thermal damping to damping from ion-neutral friction in the
modified regime is

$\mathscr{R}_{\rm f,s}=(\chi^{-1}-1)\tau_{\rm c}\nu_{ni}\beta_{\rm
M}\frac{\mathscr{F}_{\rm f,s}(\beta_{\rm M},\theta)}{\mathscr{G}_{\rm
f,s}(\beta_{\rm M},\theta)}\,,$ (37)

where $\mathscr{F}_{\rm f,s}$ and $\mathscr{G}_{\rm f,s}$ are given by
Equations (A7) and (B5) respectively. For $\beta_{\rm M}\gg 1$, where the fast
mode is continuous, we have $\beta_{\rm M}\mathscr{F}_{\rm f}/\mathscr{G}_{\rm
f}\rightarrow\beta_{\rm M}^{2}/\sin^{2}{\theta}$; for $\beta_{\rm M}\ll 1$,
the slow mode is continuous, and $\beta_{\rm M}\mathscr{F}_{\rm
s}/\mathscr{G}_{\rm s}\rightarrow\cot^{2}{\theta}$. We conclude that the
magnitude of $\mathscr{R}_{\rm f,s}$ for the continuous mode (where it exists,
see Figure 3) is, generally, well described by the product $\tau_{\rm
c}\nu_{ni}$; for large $\beta_{\rm M}$ and/or small $\theta$, it is further
enhanced.

The ratio $\mathscr{R}$, evaluated at $\theta=\pi/4$ for the continuous mode,
as well as the wavelengths of the decoupling and damping peak, $\lambda_{\rm
dec}$ and $\lambda_{\rm peak}$, are given in Table 2 for different
environments considered above. The left group of columns corresponds to our
standard value of $c_{\rm LA}/c_{n}=1.4$. The right group shows that the
thermal damping becomes much more important in a higher-$\beta$ regime. We
conclude that the effect of thermal damping on the cascade at large scales is
very important for the cloud envelope, somewhat relevant for the CNM and outer
core environments, and probably unimportant for the inner core. In cases where
the thermal damping is not important in the loaded regimes, it may still play
a significant role in the independent cascade of neutrals, which occurs at
smaller scales. The large difference in $\mathscr{R}$ between the cloud
envelope and the outer core phases is largely attributable to the sharp drop-
off in ionization fraction, assumed to occur between those densities.
$\mathscr{R}$ is large in the WNM, although we note that collisionless damping
may play a significant role in the WNM as well, depending on the ionization
fraction.

##  V. Turbulent cascade

In order to assess the importance of the thermal damping mechanism, we include
it in a simple model of the turbulent cascade. We consider both a Kolmogorov
and a Kraichnan cascade. A model equation for the steady-state spectral energy
density of waves with the amplitude damping rate $\Gamma$ is given in Ptuskin
et al. (2006):

$\frac{d}{dk}\left(\frac{kW}{\tau_{\rm nl}}\right)=-2\Gamma W.$ (38)

Here, $\Gamma(k)=-{\rm Im}~{}\omega(k)$ and $W(k)$ is the turbulent energy
spectrum of magnetic field fluctuations. We assume equipartition between the
energy spectrum of magnetic fluctuations and the kinetic energy density of the
turbulence, such that

$kW(k)=\frac{1}{2}\rho v^{2}(k).$ (39)

$\tau_{\rm nl}(W,k)$ is the characteristic timescale of the nonlinear cascade.
This is given in Ptuskin et al. (2006) as

$\tau_{\rm nl}^{\rm Kol}=C_{\rm Kol}\sqrt{\frac{\rho}{k^{3}W}}\,,$ (40)

and

$\tau_{\rm nl}^{\rm Kr}=C_{\rm Kr}\frac{\rho c_{\rm LA}}{k^{2}W}\,,$ (41)

where $C_{\rm Kol}\approx 12$ and $C_{\rm Kr}\approx 1$.

Figure 5 shows the steady-state turbulent spectrum provided by Equation (38)
with the Kolmogorov cascade time given by Equation (40). We took $\Gamma(k)$
appropriate for the continuous mode (depicted by the green line in Figure 1).
The left panels are for the same parameters as considered in Figure 1. In the
right panels, the magnetic field strength has been reduced by a factor of
$\approx 2.8$, so that $c_{\rm LA}/c_{n}=0.5$ (corresponding to $\beta=4.8$).
Different color curves correspond to different turbulent velocities at the
injection scale (10 pc for the CNM, and 1 pc for the cloud envelope). In all
cases, we have used the wave modes for a propagation angle of $45^{\circ}$ as
representative. As shown in Section III.4, it is only for a narrow range of
$\beta$ that there is a strong variation of the damping of the continuous mode
with propagation direction.111Of the four cases shown in Figures 5 and 6, the
damping strongly depends on $\theta$ only for high $\beta$ and the CNM
environment.

It is evident that the thermal damping has a profound effect on the turbulent
spectrum. Comparing the right and left panels in Figure 5, it is apparent that
even a relatively modest decrease in the magnetic field (below the value which
we consider in our “standard” case depicted in Figure 1) allows the cascade to
proceed to much smaller $\lambda$ if the thermal damping is not considered. On
the contrary, when thermal damping is included, the cascade is cut off at
slightly larger $\lambda$ as the magnetic field is reduced. This is because
the restoring force at lower magnetization is more due to gas pressure and is
thus more heavily influenced by the cooling.

Figure 5.— The turbulent energy spectrum described by Equation (38), with
$\tau_{\rm nl}$ given by Equation (40) appropriate for a Kolmogorov cascade.
Different colors correspond to different injection velocities $v_{0}$, as
shown in the legend in terms of the sonic Mach number
$\mathcal{M}_{n}=v_{0}/c_{n}$. Thermal damping is included for the solid
lines, but omitted for the dotted lines. Turbulence is injected at a scale
$\lambda_{0}$ of 1 pc for the cloud envelope, and 10 pc for the CNM. In the
left panels, the magnetic field is chosen such that $c_{\rm LA}=1.4c_{n}$, as
in Figure 1; in the right panels, $c_{\rm LA}=0.5c_{n}$. The colored vertical
lines correspond to our analytic estimate for the damping scale $\lambda_{\rm
damp}$ provided in Equation (43). The black vertical lines indicate the
condition Re $\omega=\tau_{\rm c}^{-1}$. Figure 6.— Analogous to Figure 5, but
$\tau_{\rm nl}$ was taken from Equation (41) appropriate for a Kraichnan
cascade.

The black vertical lines in Figure 5 show the locations of the peak of the
damping, $\lambda_{\rm peak}$, approximated by the condition ${\rm
Re}~{}\omega=\tau_{\rm c}^{-1}$ (and indicated in Figure 1 by the green
vertical line for the continuous mode). The turbulence is cut off near or
slightly above these values of $\lambda$. For stronger turbulence at the
injection scale, the cascade is damped at shorter $\lambda$, since $\tau_{\rm
nl}$ is smaller for larger $W$.

In the modified loaded regime $\omega\ll\tau_{\rm c}^{-1}$, the thermal
damping is much stronger than the damping due to ion-neutral friction. We
substitute the damping rate $\Gamma(k)=Fk^{2}$, where $F$ is evaluated from
Equation (A6) at $\theta=45^{\circ}$. Then, Equation (38) with $\tau_{\rm nl}$
from Equation (40) has the solution

$\frac{k^{\frac{5}{3}}W(k)}{k_{0}^{\frac{5}{3}}W_{0}}=1-\frac{\sqrt{2}\pi
C_{\rm
Kol}F}{\mathcal{M}_{n}c_{n}\lambda_{0}}\left[\left(\frac{k}{k_{0}}\right)^{\frac{4}{3}}-1\right],$
(42)

where $\mathcal{M}_{n}=v_{0}/c_{n}$ is the sonic Mach number, introduced via
the relation $k_{0}W_{0}=\frac{1}{2}\rho v_{0}^{2}$, with the subscript 0
referring to the value of the quantity at the injection scale. We can solve
Equation (42) for $k_{\rm damp}$ such that $W(k_{\rm damp})=0$, which gives a
damping wavelength $\lambda_{\rm damp}\equiv 2\pi/k_{\rm damp}$ of

$\lambda_{\rm damp}^{\rm
Kol}=\lambda_{0}\left(1+\frac{\mathcal{M}_{n}c_{n}\lambda_{0}}{\sqrt{2}\pi
C_{\rm Kol}F}\right)^{-\frac{3}{4}}.$ (43)

These values for different $\mathcal{M}_{n}$ are shown by the vertical lines
in Figure 5.

Figure 6 is the same as Figure 5, but $\tau_{\rm nl}$ is given by Equation
(41), appropriate for a Kraichnan cascade. The differences are very minor.
Using Equation (39), one can write the ratio of the cascade times given by
Equations (40) and (41) as $\tau_{\rm nl}^{\rm Kol}/\tau_{\rm nl}^{\rm
Kr}\approx 8.4v(k)/c_{\rm LA}$. This ratio is of order unity for most of our
examples, so the cutoff occurs at similar $k$.

Under the same assumptions as for the Kolmogorov case, by using $\tau_{\rm
nl}^{\rm Kr}$ in Equation (38) we find an analytic solution:

$\frac{k^{\frac{3}{2}}W(k)}{k_{0}^{\frac{3}{2}}W_{0}}=1-\frac{8\pi C_{\rm
Kr}F}{3\mathcal{M}_{\rm LA}^{2}c_{\rm
LA}\lambda_{0}}\left[\left(\frac{k}{k_{0}}\right)^{\frac{3}{2}}-1\right],$
(44)

where $\mathcal{M}_{\rm LA}=v_{0}/c_{\rm LA}$ is the Alfvénic Mach number.
Setting $W(k)$ to zero, we solve for the damping wavelength

$\lambda_{\rm damp}^{\rm Kr}=\lambda_{0}\left(1+\frac{3\mathcal{M}_{\rm
LA}^{2}c_{\rm LA}\lambda_{0}}{8\pi C_{\rm Kr}F}\right)^{-\frac{2}{3}},$ (45)

shown by the vertical lines in Figure 6. We point out that if $\lambda_{\rm
damp}$ were significantly smaller than the wavelength of peak damping
$\lambda_{\rm peak}$, then the above analysis (assuming $\Gamma\propto k^{2}$)
would break down. In this case the turbulence would continue to cascade to
smaller wavelength (where $\Gamma$ is constant) until it is eventually damped
by viscosity. For the environments we consider, however, this would only
happen if $\mathcal{M}_{n}>1$, at which point our analysis is not valid
anyway.

##  VI. MHD Simulations with cooling

Figure 7.— Kinetic energy spectra $W_{\rm kin}(k)$ from MHD simulations with
gas heating and cooling (black), as well as with isothermal (green) and
adiabatic (blue) equations of state. The vertical axis is in arbitrary units
and the curves in each simulation are shifted vertically for clarity. In each
case, the solid line shows the time-averaged spectrum at steady state for the
total kinetic energy and the dotted line represents the compressive modes. The
shaded area indicates the $1\sigma$ dispersion among the different time
snapshots. The gray dash-dotted line indicates the slope for the Kolmogorov
spectrum of $k^{-5/3}$. The vertical gray dashed line shows the location for
the fastest damped mode at $k_{\rm peak}/2\pi=1/\lambda_{\rm peak}$ from the
linear theory. In the isothermal or adiabatic cases, the energy spectrum of
the compressive modes has a similar shape to the total energy spectrum and
only slightly less power. In the case with gas heating and cooling, the
compressive mode is damped strongly, with a much steeper spectral slope and
much less power.

In addition to the analytic approach followed in the previous sections, we
also ran MHD simulations with continuously driven turbulence. We solve the
ideal MHD equations (see, e.g., Stone et al. 2008),222In Equations (47)–(50),
the $1/\sqrt{4\pi}$ pre-factor is absorbed in the unit of $\mathbf{B}$.

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0,$ (46)
$\frac{\partial\rho\mathbf{v}}{\partial
t}+\nabla\cdot(\rho\mathbf{vv}-\mathbf{BB})+\nabla\left(p+\frac{1}{2}B^{2}\right)=0,$
(47) $\frac{\partial E}{\partial
t}+\nabla\cdot\left[(p+E)\mathbf{v}-(\mathbf{B\cdot
v})\mathbf{B}\right]=\Lambda_{\rm h}-\Lambda_{\rm c}\,,$ (48)
$E=\frac{p}{\gamma-1}+\frac{1}{2}\rho v^{2}+\frac{1}{2}B^{2}\,,$ (49)
$\frac{\partial\mathbf{B}}{\partial
t}-\nabla\times(\mathbf{v}\times\mathbf{B})=0,$ (50)

where $\rho$ is the mass density of gas, $\mathbf{v}$ is the gas velocity, $p$
is the gas pressure, $\mathbf{B}$ is the magnetic field, $E$ is the total
energy density, and $\Lambda_{\rm h}$ and $\Lambda_{\rm c}$ are the gas
heating and cooling rates. The adiabatic index $\gamma$ is set to $5/3$.

The values of $\Lambda_{\rm h}$ and $\Lambda_{\rm c}$ in our simulations
correspond to the CNM conditions, as described in Section III.1.1. Here we
also include helium with abundance $x_{\rm He}=0.1$. The mean molecular weight
of the gas is then 1.4 accounting for helium. The specific heat capacity per H
atom at constant volume is given by $C_{v}^{\rm H}=1.5k_{\rm
B}(1+x_{\mathrm{He}})$, ignoring the minor contribution from metals and
electrons. In addition, we run two simulations with an isothermal and an
adiabatic equation of state. For the adiabatic simulation, $\Lambda_{\rm
h}=\Lambda_{\rm c}=0$; for the isothermal simulation, Equations (48) and (49)
are replaced by an isothermal equation of state, $p/\rho=c_{n}^{2}/\gamma$,
and the sound speed $c_{n}$ is calculated at the equilibrium CNM temperature
of 60 K (see Section III.1.1).

We perform the simulations using the MHD code Athena++ (Stone et al. 2020).
The simulation set-up is very similar to the driving turbulence simulations in
Gong et al. (2020), and we refer the readers to that paper for more details.
We adopt a Cartesian coordinate system with periodic boundary conditions on
all sides. We use the Harten-Lax-van Leer-Discontinuities (HLLD) Riemann
solver (Miyoshi & Kusano 2005), with the third-order Runge–Kutta integrator
and third-order spatial reconstruction. For each simulation, turbulence was
driven with an energy spectrum of $E(k)\propto k^{-3}$ at integer values of
$[\tilde{k}_{x},\tilde{k}_{y},\tilde{k}_{z}]$ such that
$0<\tilde{k}=\sqrt{\tilde{k}_{x}^{2}+\tilde{k}_{y}^{2}+\tilde{k}_{z}^{2}}<2$.
$\tilde{k}\equiv k/k_{L}$, where $k_{L}$ is the value corresponding to the box
size. The turbulence is driven with a fully compressive velocity field, but
there is a small solenoidal component in the driving due to the finite forcing
correlation time (Grete et al. 2018). The simulations are run with a box-size
of $L_{x}=L_{y}=L_{z}=1~{}\mathrm{pc}$, and a resolution of
$N_{x}=N_{y}=N_{z}=180$. For the fast mode, with $\theta=45^{\circ}$, thermal
damping is peaked at $\lambda_{\rm peak}=0.13$ pc, and we assume this to be
the scale at which the damping effect is strongest. With this resolution, the
scale $\lambda_{\rm peak}$ is well resolved by 23 cells. In the isothermal
simulation, we adopt the sound speed at the equilibrium temperature, 0.62
km/s. In the adiabatic simulation, we use the adiabatic index of 5/3 and start
with the gas at the equilibrium temperature. The gas heats up gradually (the
temperature increases by about 10% at the end of the simulation) as kinetic
energy is injected by turbulence driving. In the simulation, we include the
heating and cooling as explicit source terms in the energy equation. We limit
the simulation time step so that the gas temperature in each cell changes less
than 10% during each time step.

We run the simulations for about 1000 Myr, more than 10 turbulence crossing
times, and analyse the results at steady state between 600 – 1000 Myr. The
steady state sonic turbulent Mach number $\mathcal{M}_{n}$ at a scale of 1 pc
is 0.02–0.08. We pick this low value so as to avoid the introduction of non-
linear effects which might complicate the analysis. We start with a plasma
$\beta=2$ in all simulations, and its value does not vary significantly in
these cases due to the weak turbulence. We experimented with a lower
resolution of $N_{x}=N_{y}=N_{z}=100$, a higher steady-state Mach number
around 0.2, and a higher $\beta$ around 10, and found very similar results to
the simulations shown here.

The energy spectra of the turbulence from the simulations are shown in Figure
7. The lines are vertically offset for clarity. In the isothermal and
adiabatic simulations, the energy spectrum in the compressive modes are
similar to the total energy spectrum, with only slightly less power and
slightly steeper slopes. With realistic cooling, the compressive modes decay
much more rapidly with $k$ than the solenoidal ones, resulting in a much
steeper spectrum and much less power in the compressive modes.

The agreement with the analytic model is at this stage only qualitative. Due
to the limited numerical resolution, we cannot precisely describe the
turbulent spectrum. At high $k$, the spectrum is affected by numerical
damping, and at low $k$ the cascade is not fully developed. In the analytic
case, the damping cuts off the spectrum sharply at a particular scale. In the
simulation, there is not a sharp cutoff, but rather a substantially steeper
slope over a wide range of $k$. Despite of these limitations, however, the
simulations clearly show that the compressive modes of the turbulence are
damped efficiently by gas cooling.

##  VII. Implications

The detailed implications of thermal damping of compressive modes require more
realistic simulations which will better determine the magnitude of the effect
on the turbulent cascade in realistic conditions. However, in this section we
outline a few different possible area in which thermal damping may be
important.

###  VII.1. CR propagation

The CR propagation in the ISM is often modelled as a diffusion process (Strong
& Moskalenko 1998; Strong et al. 2007). The scattering which determines the
diffusion coefficient primarily arises from perturbations of the magnetic
field which are resonant with the gyroradius of CRs (Berezinskii et al. 1984;
Strong et al. 2007). It was shown in Yan & Lazarian (2004) that in realistic
ISM conditions, fast modes dominate the scattering of CRs. Here we demonstrate
that, depending on the plasma $\beta$, such modes are strongly damped in the
WNM and molecular cloud environments, which thus can significantly affect CR
propagation.

###  VII.2. Dust growth

Dust coagulation in molecular clouds is sensitive to the degree of turbulence
on very small scales (Voelk et al. 1980; Yan et al. 2004; Ormel & Cuzzi 2007),
and can begin at densities as low as $10^{4}$ cm-3 (Hirashita & Yan 2009;
Silsbee et al. 2020). If the turbulence is damped at larger scales, the
evolution of the dust population may proceed more slowly. This can have
implications for the gas freeze-out on the surface of dust grains in dark
clouds (which also occurs around densities of $10^{4}$ cm-3). Additionally,
although the volume filling factor of such gas is small, coagulation at
somewhat lower densities ($10^{3}$ cm-3) can affect the grain size
distribution in the ISM as a whole, if such material is mixed back into the
more diffuse phases (see, e.g., Hirashita & Aoyama 2019).

###  VII.3. Observations of small-scale motions

Although small-scale motions cannot typically be directly resolved in
observations, there are strong indications that turbulence on small scales
does exist. This follows both from scintillation observations of electron
density fluctuations in the diffuse ISM (Armstrong et al. 1995) and
observations showing different velocity dispersions of ionized and neutral
tracers in molecular clouds (Li & Houde 2008; Pineda et al. 2021). A
quantitative understanding of the relevant damping processes is therefore
essential for the development of a theoretical framework which could explain
these observations.

##  VIII. Conclusion

In this work we introduced and analyzed the effect of thermal damping,
associated with radiative cooling of gas, on compressive MHD turbulence in the
ISM. For a particular wave mode, the scale at which the damping is maximized
depends on the cooling timescale and the plasma $\beta$. In some environments
the damping peaks at scales hundreds of times the ambipolar damping scale –
the characteristic wavelength at which the ions and neutrals decouple. We
confirm, using a simple analytic model of the turbulent cascade as well as MHD
simulations including cooling, that thermal damping has a substantial effect
on the turbulent spectrum of compressive waves for parameters appropriate for
several characteristic phases of the ISM. We find, for typical turbulent
velocities, that the turbulent cascade is cut off at a scale comparable to the
scale at which the thermal damping peaks. The spatial scale of this peak
varies from environment to environment. In typical WNM conditions, the peak
occurs at scales of around 100 pc, whereas in the CNM or in low-density
molecular gas, the effect peaks around 0.1 pc. At densities greater than
$10^{5}$ cm-3, typical for the inner parts of molecular cloud cores, the peak
becomes on the order of a milliparsec.

We derived Equation (37) for the strength $\mathscr{R}$ of the thermal damping
relative to damping from ion-neutral friction. Its magnitude is determined by
the product $\tau_{\rm c}\nu_{ni}$ of the cooling timescale and the neutral-
ion collision frequency. For the conditions we considered, thermal damping is
most efficient in low-density molecular gas, but can also play a significant
role in the CNM and outer parts of a prestellar core, (as well as in the WNM,
although collisionless damping may be important there as well).

This effect is important for a number of reasons. An open question in the
study of ISM turbulence is how the turbulent energy is transferred from large
scales, where ions and neutrals move together, to small scales, where they
move separately (see Tilley & Balsara 2011; Burkhart et al. 2015). This
process may involve compressive modes, since Alfvén waves in low-ionization
media have a relatively wide range of $k$ in which they cannot propagate. For
this reason, it is important to have a complete understanding of the
mechanisms which impact the survival of compressive waves. Additionally, the
damping described in this paper has potential implications for dust evolution
in the ISM, propagation of CRs, and the interpretation of observed line widths
of ionized versus neutral gas tracers in molecular clouds.

##  Appendix A Appendix A: Rate of Thermal Damping

Thermal damping generally dominates in the modified loaded regime
$\omega\tau_{\rm c}\ll 1$. This allows us to ignore ion-neutral friction and
derive an analytic formula for the damping rate, considering the limit of
small ionization fraction for simplicity. The dispersion relation of fast and
slow modes in the loaded regimes (on the right from the decoupling gap in
Figure 1) is generally given by the standard expression (Landau & Lifshitz
1960) where $c_{\rm A}$ is replaced, according to Equation (15), with $c_{\rm
LA}$, while $c_{i}$ is replaced with $c_{n}f(\omega)$.

$\omega^{4}-\omega^{2}k^{2}c_{n}^{2}f(\omega)-\omega^{2}k^{2}c_{\rm
LA}^{2}+k^{4}\cos^{2}{\theta}c_{n}^{2}c_{\rm LA}^{2}f(\omega)=0.$ (A1)

In the modified loaded regime, where $\omega\tau_{\rm c}\ll 1$, we can
approximate Equation (13) as $f(\omega)=\chi(1-\eta\omega\tau_{\rm c})$, where

$\eta=\chi^{-1}-1$ (A2)

is a positive number. In this case, Equation (A1) becomes

$\omega^{4}-\omega^{2}k^{2}\left[c_{\rm LA}^{2}+c_{{\rm
M}n}^{2}(1-i\eta\omega\tau_{\rm c})\right]+k^{4}c_{\rm LA}^{2}c_{{\rm
M}n}^{2}(1-i\eta\omega\tau_{\rm c})\cos^{2}{\theta}=0,$ (A3)

where $c_{{\rm M}n}$ is given by Equation (29).

We can present the solution for the modified loaded (fast or slow) mode as

$\omega=kc_{\rm Mf,s}+i\>{\rm Im}~{}\omega,$ (A4)

where $c_{\rm Mf,s}$ is the fast or slow mode speed, given by Equation (28)
with $c_{\rm A}$ replaced with $c_{\rm LA}$, and $c_{i}$ replaced with
$c_{{\rm M}n}$. Substituting Equation (A4) into Equation (A3), and keeping
only first-order terms in ${\rm Im}~{}\omega$ gives

$-{\rm Im}~{}\omega=\frac{\eta\tau_{\rm c}c_{{\rm
M}n}^{2}}{2}\left(\frac{c_{\rm Mf,s}^{2}-c_{\rm
LA}^{2}\cos^{2}{\theta}}{2c_{\rm Mf,s}^{2}-c_{\rm LA}^{2}-c_{{\rm
M}n}^{2}}\right)k^{2}\equiv F_{\rm f,s}k^{2}.$ (A5)

After some manipulation, we obtain

$F_{\rm f,s}=\frac{\eta\tau_{\rm c}c_{{\rm M}n}^{2}}{4}\,\mathscr{F_{\rm
f,s}}(\beta_{\rm M},\theta),$ (A6)

where $\beta_{\rm M}=c_{{\rm M}n}^{2}/c_{\rm LA}^{2}\equiv\gamma\chi\beta/2$,
and

$\mathscr{F}_{\rm f,s}(\beta_{\rm
M},\theta)=1\mp\frac{\cos{2\theta}-\beta_{\rm
M}}{\sqrt{(\cos{2\theta}-\beta_{\rm M})^{2}+\sin^{2}{2\theta}}}\,.$ (A7)

We note that $\mathscr{F}_{\rm s}+\mathscr{F}_{\rm f}=2$.

##  Appendix B Appendix B: Loaded regime and decoupling gap

A general dispersion relation given by Equations (23)–(25) can be
substantially simplified for wave modes sustained in the loaded regime and in
the decoupling gap (see Section III.2). For this, we require the ratio of the
ion-to-neutral densities, $\rho_{i}/\rho_{n}$, to be a sufficiently small
number, so that the relative width of the decoupling gap is large,
$\sim\sqrt{\rho_{n}/\rho_{i}}$ (e.g., Kulsrud & Pearce 1969). As damping near
the decoupling gap is dominated by ion-neutral friction, here we can neglect
thermal damping and set $f(\omega)=1$.

Let us identically rewrite Equation (24), by replacing in the first term
$\Delta_{n}$ with $\Delta_{i}$ and subtracting from the resulting equation the
corresponding difference $\propto(\Delta_{i}-\Delta_{n})$. Equation (25) for
neutrals is rewritten in a similar way, by replacing $\Delta_{n}^{\|}$ with
$\Delta_{i}^{\|}$ in the first term, and $\Delta_{n}$ with $\Delta_{i}$ in the
second term. Then, excluding $\Delta_{i}-\Delta_{n}$ and
$\Delta_{i}^{\|}-\Delta_{n}^{\|}$ from the resulting four equations, we
obtain:

$\left[\omega^{2}-k^{2}(c_{\rm
A}^{2}+c_{i}^{2})+\frac{i\omega\nu_{in}}{\omega^{2}+i\omega\nu_{ni}-k^{2}c_{n}^{2}}\,(\omega^{2}-k^{2}c_{n}^{2})\right]\Delta_{i}+k^{2}c_{\rm
A}^{2}\,\Delta_{i}^{\|}=0,$ (B1)
$\omega^{2}\,\Delta_{i}^{\|}-\frac{i\omega\nu_{ni}}{\omega^{2}+i\omega\nu_{ni}-k^{2}c_{n}^{2}}\,k_{z}^{2}c_{n}^{2}\,\Delta_{i}=0.$
(B2)

We are interested in wave modes with $|\omega|\sim\nu_{ni}$, corresponding to
the right edge of the decoupling gap (see below), and therefore in deriving
Equation (B2), we neglected terms
$\propto\nu_{ni}/\nu_{in}=\rho_{i}/\rho_{n}$, and terms
$\propto\omega/\nu_{in}\equiv(\rho_{i}/\rho_{n})(\omega/\nu_{ni})$. For the
same reason, the third term in the brackets in Equation (B1) is estimated as
$\sim(\nu_{in}/\nu_{ni})\omega^{2}$, and therefore the first term $\omega^{2}$
can be omitted. We can also neglect $c_{i}^{2}$ next to $c_{\rm A}^{2}$,
because $c_{\rm LA}\equiv c_{\rm A}\sqrt{\rho_{i}/\rho_{n}}$ in our analysis
is assumed to be comparable to $c_{n}\sim c_{i}$. Substituting
$\nu_{in}=(\rho_{n}/\rho_{i})\nu_{ni}$, after some manipulation we obtain the
following dispersion relation:

$k^{2}c_{\rm
LA}^{2}\omega(\omega^{2}-k^{2}c_{n}^{2})=i\nu_{ni}(\omega^{2}-k^{2}c_{\rm
Lf}^{2})(\omega^{2}-k^{2}c_{\rm Ls}^{2}),$ (B3)

where $c_{\rm Lf,s}^{2}$ are given by Equation (28) with $c_{\rm A}$ replaced
by $c_{\rm LA}$, and $c_{i}$ replaced by $c_{n}$.

Equation (B3) has an intuitive structure. For sufficiently small $k$,
representing wavelengths $\lambda$ to the right of the decoupling gap (see
Section III.2 and Figure 1), the LHS is small (since $\omega\ll\nu_{ni}$) and
thus we recover the loaded fast and slow modes, Re $\omega\approx kc_{\rm
Lf,s}$ and

$-{\rm Im}~{}\omega=\frac{c_{\rm LA}^{2}}{2\nu_{ni}}\frac{|c_{\rm
Lf,s}^{2}-c_{n}^{2}|}{c_{\rm Lf}^{2}-c_{\rm Ls}^{2}}k^{2}\equiv\frac{c_{\rm
LA}^{2}}{4\nu_{\rm ni}}\mathscr{G}_{\rm f,s}(\beta_{\rm M}/\chi,\theta)k^{2},$
(B4)

where

$\mathscr{G}_{\rm f,s}(\beta_{\rm M}/\chi,\theta)=1\pm\frac{1-\beta_{\rm
M}/\chi}{\sqrt{\left(1-\beta_{\rm M}/\chi\right)^{2}+4(\beta_{\rm
M}/\chi)\sin^{2}{\theta}}}\,.$ (B5)

Note that the effect of ion-neutral friction in the modified loaded regime is
obtained by following the same derivation for $f(\omega)=\chi$: this yields Re
$\omega\approx kc_{\rm Mf,s}$ and ${\rm Im}~{}\omega$ given by Equation (B4)
with $c_{\rm Lf,s}\to c_{\rm Mf,s}$ and $c_{n}\to c_{{\rm M}n}$, which is
equivalent to $\chi\to 1$ in Equation (B5). For large $k$, corresponding to
the decoupling gap, we recover the neutral sound mode with Re $\omega\approx
kc_{n}$ and $-{\rm Im}~{}\omega\propto\nu_{ni}$. Equation (B3) is not valid on
the left side of the decoupling gap, because the assumption we made that
$\omega\ll\nu_{in}$ breaks down in this regime.

In order to obtain Figure 3, we introduce dimensionless wave speed
$\tilde{c}$, fast and slow speeds $\tilde{c}_{\rm Lf,s}$, and wavelength
$\tilde{\lambda}$, defined as

$\tilde{c}=\frac{\omega}{kc_{n}},\quad\tilde{c}_{\rm Lf,s}=\frac{c_{\rm
Lf,s}}{c_{n}}\quad{\rm and}\quad\tilde{\lambda}=\frac{c_{n}}{c_{\rm
LA}}\frac{\nu_{ni}}{kc_{\rm LA}}\,.$ (B6)

Then Equation (B3) becomes

$\tilde{c}(\tilde{c}^{2}-1)=i\tilde{\lambda}\left(\tilde{c}^{2}-\tilde{c}_{\rm
Lf}^{2}\right)\left(\tilde{c}^{2}-\tilde{c}_{\rm Ls}^{2}\right).$ (B7)

A condition for the neutral sound mode to switch between the loaded fast and
slow modes (depicted by the solid line in Figure 3) is easily derived by
substituting $\tilde{c}=-1/\tilde{c}_{*}$ in Equation (B7). This transforms
the fast mode into the conjugate slow mode (with the real part of the opposite
sign) and vice versa. The resulting equation for $\tilde{c}_{*}$ is then
reduced to the form of Equation (B7) by setting $c_{\rm Lf}^{2}c_{\rm
Ls}^{2}=c_{n}^{4}$, which yields the sought condition. Using Equation (28) for
$c_{\rm Lf,s}^{2}$, we obtain that the switch between the fast and slow modes
occurs at

$\frac{c_{n}}{c_{\rm LA}}=\cos\theta.$ (B8)

We now can derive an approximate condition for the right edge of the
decoupling gap. Substituting $x=i\tilde{c}$ transforms Equation (B7) into an
equation with all real coefficients. In the case that $c_{\rm Lf}\gg c_{\rm
Ls}$, we can find approximate solutions for the loaded fast mode by discarding
the linear and constant terms in the transformed equation, and for the loaded
slow mode – by discarding the 3rd- and 4th-order terms. This results in the
following equations:

$\displaystyle{\rm Loaded~{}fast:}$ $\displaystyle\quad
x^{2}-\frac{x}{\tilde{\lambda}}+\tilde{c}_{\rm Lf}^{2}+\tilde{c}_{\rm
Ls}^{2}=0,$ (B9) $\displaystyle{\rm Loaded~{}slow:}$
$\displaystyle\quad\left(\tilde{c}_{\rm Lf}^{2}+\tilde{c}_{\rm
Ls}^{2}\right)x^{2}-\frac{x}{\tilde{\lambda}}+\tilde{c}_{\rm
Lf}^{2}\tilde{c}_{\rm Ls}^{2}=0.$ (B10)

We keep in mind that if Equation (B9) or (B10) has a pair of real solutions,
then the corresponding mode has no real solutions (i.e. the mode is
discontinuous). For $c_{n}\gg c_{\rm LA}$, we expect the loaded slow mode to
be cut off at the right edge of the decoupling gap and the loaded fast mode to
be continuous with the neutral sound mode (see Section III.4). Hence, by
setting the discriminant of Equation (B10) to zero, we can solve for the right
edge in this limit: $k=2(\nu_{ni}/c_{\rm LA})\cos{\theta}$. Conversely, for
$c_{n}\ll c_{\rm LA}$ we expect the loaded fast mode to be cut off; setting
the discriminant of Equation (B9) to zero yields $k=2\nu_{ni}/c_{\rm LA}$.

A condition that the continuous mode ceases to exist (depicted by the dashed
lines in Figure 3) is equivalent to that Equation (B7) written for
$x=i\tilde{c}$ has four real roots at $0<x<\infty$. In principle, this
condition could be derived analytically, using Sturm’s theorem for the number
of real roots of a polynomial, but we found it more convenient to compute
these boundaries numerically.

## References

  * Armstrong et al. (1995) Armstrong, J. W., Rickett, B. J., & Spangler, S. R. 1995, ApJ, 443, 209
  * Balsara (1996) Balsara, D. S. 1996, ApJ, 465, 775
  * Berezinskii et al. (1984) Berezinskii, V. S., Bulanov, S. V., Ginzburg, V. L., Dogel, V. A., & Ptuskin, V. S. 1984, Astrophysics of cosmic rays.
  * Bunte & Bogdan (1994) Bunte, M., & Bogdan, T. J. 1994, A&A, 283, 642
  * Burkhart et al. (2015) Burkhart, B., Lazarian, A., Balsara, D., Meyer, C., & Cho, J. 2015, ApJ, 805, 118
  * Caselli et al. (2002) Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002, ApJ, 565, 344
  * Chepurnov & Lazarian (2010) Chepurnov, A., & Lazarian, A. 2010, ApJ, 710, 853
  * Crutcher et al. (2010) Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., & Troland, T. H. 2010, ApJ, 725, 466
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium
  * Ginzburg (1970) Ginzburg, V. L. 1970, The propagation of electromagnetic waves in plasmas
  * Glazebrook (2013) Glazebrook, K. 2013, PASA, 30, e056
  * Goldsmith (2001) Goldsmith, P. F. 2001, ApJ, 557, 736
  * Gong et al. (2020) Gong, M., Ivlev, A. V., Zhao, B., & Caselli, P. 2020, ApJ, 891, 172
  * Gong et al. (2017) Gong, M., Ostriker, E. C., & Wolfire, M. G. 2017, ApJ, 843, 38
  * Grete et al. (2018) Grete, P., O’Shea, B. W., & Beckwith, K. 2018, ApJ, 858, L19
  * Hirashita & Aoyama (2019) Hirashita, H., & Aoyama, S. 2019, MNRAS, 482, 2555
  * Hirashita & Yan (2009) Hirashita, H., & Yan, H. 2009, MNRAS, 394, 1061
  * Ivlev et al. (2019) Ivlev, A. V., Silsbee, K., Sipilä, O., & Caselli, P. 2019, ApJ, 884, 176
  * Krumholz & Burkhart (2016) Krumholz, M. R., & Burkhart, B. 2016, MNRAS, 458, 1671
  * Kulsrud & Pearce (1969) Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445
  * Landau & Lifshitz (1960) Landau, L. D., & Lifshitz, E. M. 1960, Electrodynamics of continuous media
  * Li & Houde (2008) Li, H.-b., & Houde, M. 2008, ApJ, 677, 1151
  * Mac Low & Klessen (2004) Mac Low, M.-M., & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125
  * Marchal et al. (2021) Marchal, A., Martin, P. G., & Gong, M. 2021, arXiv e-prints, arXiv:2106.12683
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565
  * Mihalas & Mihalas (1983) Mihalas, D., & Mihalas, B. W. 1983, ApJ, 273, 355
  * Minter & Spangler (1996) Minter, A. H., & Spangler, S. R. 1996, ApJ, 458, 194
  * Miyoshi & Kusano (2005) Miyoshi, T., & Kusano, K. 2005, Journal of Computational Physics, 208, 315
  * Ormel & Cuzzi (2007) Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413
  * Ossenkopf (1993) Ossenkopf, V. 1993, A&A, 280, 617
  * Pineda et al. (2021) Pineda, J. E., Schmiedeke, A., Caselli, P., et al. 2021, ApJ, 912, 7
  * Ptuskin et al. (2006) Ptuskin, V. S., Moskalenko, I. V., Jones, F. C., Strong, A. W., & Zirakashvili, V. N. 2006, ApJ, 642, 902
  * Raizer et al. (2011) Raizer, Y., Kisin, V., & Allen, J. 2011, Gas Discharge Physics (Springer Berlin Heidelberg)
  * Silsbee & Ivlev (2019) Silsbee, K., & Ivlev, A. V. 2019, ApJ, 879, 14
  * Silsbee et al. (2020) Silsbee, K., Ivlev, A. V., Sipilä, O., Caselli, P., & Zhao, B. 2020, A&A, 641, A39
  * Soler et al. (2013) Soler, R., Carbonell, M., & Ballester, J. L. 2013, ApJS, 209, 16
  * Souffrin (1972) Souffrin, P. 1972, A&A, 17, 458
  * Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
  * Stone et al. (2020) Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4
  * Strömgren (1939) Strömgren, B. 1939, ApJ, 89, 526
  * Strong & Moskalenko (1998) Strong, A. W., & Moskalenko, I. V. 1998, ApJ, 509, 212
  * Strong et al. (2007) Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Annual Review of Nuclear and Particle Science, 57, 285
  * Tilley & Balsara (2011) Tilley, D. A., & Balsara, D. S. 2011, MNRAS, 415, 3681
  * Voelk et al. (1980) Voelk, H. J., Jones, F. C., Morfill, G. E., & Roeser, S. 1980, A&A, 85, 316
  * Williams et al. (1998) Williams, J. P., Bergin, E. A., Caselli, P., Myers, P. C., & Plume, R. 1998, ApJ, 503, 689
  * Xu & Yan (2013) Xu, S., & Yan, H. 2013, ApJ, 779, 140
  * Xu et al. (2016) Xu, S., Yan, H., & Lazarian, A. 2016, ApJ, 826, 166
  * Yan & Lazarian (2004) Yan, H., & Lazarian, A. 2004, ApJ, 614, 757
  * Yan et al. (2004) Yan, H., Lazarian, A., & Draine, B. T. 2004, ApJ, 616, 895


t-r\rangle^{-\frac{1}{2}}\int_{|\omega|=1}|\int_{t-1}^{\infty}\frac{d}{d\rho}\left(\langle
x+\rho\omega\rangle^{\frac{1}{2}}(\omega\cdot\nabla)v_{0}(x+\rho\omega)\sqrt{2\rho+1}\right)d\rho|d\omega$
(B.25) $\displaystyle\lesssim$ $\displaystyle\langle
t-r\rangle^{-\frac{1}{2}}\int_{|\omega|=1}|\int_{t-1}^{\infty}\langle
x+\rho\omega\rangle^{-\frac{1}{2}}|(\omega\cdot\nabla)v_{0}(x+\rho\omega)|\frac{\sqrt{2\rho+1}}{\rho}\rho
d\rho|d\omega$ (B.26) $\displaystyle+\langle
t-r\rangle^{-\frac{1}{2}}\int_{|\omega|=1}|\int_{t-1}^{\infty}\langle
x+\rho\omega\rangle^{\frac{1}{2}}\left((\omega\cdot\nabla)^{2}v_{0}(x+\rho\omega)\sqrt{2\rho+1}+(\omega\cdot\nabla)v_{0}(x+\rho\omega)\frac{1}{\sqrt{2\rho+1}}\right)d\rho|d\omega$
(B.27) $\displaystyle\lesssim$ $\displaystyle\ t^{-\frac{1}{2}}\langle
t-r\rangle^{-\frac{1}{2}}\|\nabla v_{0}\|_{L^{1}}+t^{-\frac{1}{2}}\langle
t-r\rangle^{-\frac{1}{2}}\|\langle
y\rangle^{\frac{1}{2}}\nabla^{2}v_{0}\|_{L^{1}}+t^{-\frac{3}{2}}\langle
t-r\rangle^{-\frac{1}{2}}\|\langle y\rangle^{\frac{1}{2}}\nabla
v_{0}\|_{L^{1}}.$ (B.28)

$J_{4}$ can be handled similarly:

$\displaystyle J_{4}\lesssim t^{-\frac{1}{2}}\langle
t-r\rangle^{-\frac{1}{2}}\|\langle
y\rangle^{\frac{1}{2}}\langle\nabla\rangle^{2}v_{0}\|_{L^{1}}.$ (B.29)

Therefore we conclude (B.1) ∎

##  Appendix C Verification of the smallness conditions of the Initial data

In this section we verify the smallness condition (1.4) is consistent with the
vector fields:

######  Proposition C.1.

Assume $(\Phi_{0},\Phi_{1},N_{0},N_{1})$ satisfy the condition (1.4), then for
any $J\leq K$ we have

$\displaystyle\|\Gamma^{\leq J}\phi|_{t=0}\|_{2}\lesssim\big{(}\|\langle
x\rangle^{J}\langle\nabla\rangle^{J}\Phi_{0}\|_{2}+\|\langle
x\rangle^{J}\langle\nabla\rangle^{J-1}\Phi_{1}\|_{2}\big{)}\big{(}1+\sum_{j=0}^{J-2}\|\langle
x\rangle^{j+1}\nabla^{j}(\nabla N_{0},N_{1})\|_{2}\big{)},$ (C.1)
$\displaystyle\|\Gamma^{\leq J}n|_{t=0}\|_{2}\lesssim\big{(}\|\langle
x\rangle^{J-1}\langle\nabla\rangle^{J-1}\Phi_{0}\|_{2}+\|\langle
x\rangle^{J-1}\langle\nabla\rangle^{J-2}\Phi_{1}\|_{2}\big{)}\big{(}1+\sum_{j=0}^{J-1}\|\langle
x\rangle^{j+1}\nabla^{j}(\nabla N_{0},N_{1})\|_{2}\big{)}.$ (C.2)

In particular, the following also holds:

$\displaystyle\|\partial\Gamma^{\leq
J}\phi|_{t=0}\|_{2}\lesssim\big{(}\|\langle
x\rangle^{J}\langle\nabla\rangle^{J+1}\Phi_{0}\|_{2}+\|\langle
x\rangle^{J}\langle\nabla\rangle^{J}\Phi_{1}\|_{2}\big{)}\big{(}1+\sum_{j=0}^{J-1}\|\langle
x\rangle^{j+1}\nabla^{j}(\nabla N_{0},N_{1})\|_{2}\big{)},$ (C.3)
$\displaystyle\|\widetilde{\Gamma}\Gamma^{\leq
J-1}n|_{t=0}\|_{2}\lesssim\big{(}\|\langle
x\rangle^{J-1}\langle\nabla\rangle^{J-1}\Phi_{0}\|_{2}+\|\langle
x\rangle^{J-1}\langle\nabla\rangle^{J-2}\Phi_{1}\|_{2}\big{)}\big{(}1+\sum_{j=0}^{J-1}\|\langle
x\rangle^{j+1}\nabla^{j}(\nabla N_{0},N_{1})\|_{2}\big{)}.$ (C.4)

Similarly, we also have

$\displaystyle\|\langle x\rangle\log(\langle x\rangle)\Gamma^{\leq
K}\phi|_{t=0}\|_{2}$ (C.5) $\displaystyle\lesssim$
$\displaystyle\big{(}\|\langle x\rangle^{K+1}\log(\langle
x\rangle)\langle\nabla\rangle^{K}\Phi_{0}\|_{2}+\|\langle
x\rangle^{K+1}\log(\langle
x\rangle)\langle\nabla\rangle^{K-1}\Phi_{1}\|_{2}\big{)}\big{(}1+\sum_{j=0}^{K-2}\|\langle
x\rangle^{j+1}\nabla^{j}(\nabla N_{0},N_{1})\|_{2}\big{)},$

and

$\displaystyle\|(\langle
x\rangle^{\frac{3}{2}}\langle\nabla\rangle^{3}\widetilde{\Gamma}^{\leq
K-4}n_{0})|_{t=0}\|_{1}+\|(\langle
x\rangle^{\frac{3}{2}}\langle\nabla\rangle^{2}\partial_{t}\widetilde{\Gamma}^{\leq
K-4}n_{0})|_{t=0}\|_{1}$ (C.6) $\displaystyle\lesssim$
$\displaystyle\sum_{j=0}^{K-5}(\|\langle
x\rangle^{j+\frac{5}{2}}\nabla^{j+3}(\nabla N_{0},N_{1})\|_{1}+\|\langle
x\rangle^{j+\frac{5}{2}}\nabla^{j}(\nabla N_{0},N_{1})\|_{1}).$

All the constants in the inequalities are exact and computable; for simplicity
we identify all quantities related to the initial conditions by
$C_{i}\varepsilon$, where $C_{i}$ is a fixed exact positive constant that can
be computed directly.

###### Proof.

The proof relies on an induction that is very similar to [6], therefore we
only sketch the idea. In fact the wave filed $n$ is simpler to deal with
because of the decomposition $n=n_{0}+A^{\alpha}\partial_{\alpha}n$ so that
the initial conditions $(N_{0},N_{1})$ can be taken care of essentially
through the homogeneous equation up to some higher order perturbation. For the
Klein-Gordon field $\phi$, with no loss, we only consider the “worst” scenario
as follows:

$\displaystyle\|\langle
x\rangle^{J}\partial_{t}^{J}\phi|_{t=0}\|_{2}\lesssim\|\langle
x\rangle^{J}\partial_{t}^{J-2}\Delta\phi|_{t=0}\|_{2}+\|\langle
x\rangle^{J}\partial_{t}^{J-2}\phi|_{t=0}\|_{2}+\|\langle
x\rangle^{J}\partial_{t}^{J-2}(\phi\partial n)|_{t=0}\|_{2}.$ (C.7)

Note that by induction assumption we shall have

$\displaystyle\|\langle
x\rangle^{J}\partial_{t}^{J-2}\Delta\phi|_{t=0}\|_{2}+\|\langle
x\rangle^{J}\partial_{t}^{J-2}\phi|_{t=0}\|_{2}\lesssim\|\langle
x\rangle^{J}\langle\nabla\rangle^{J}\Phi_{0}\|_{2}+\|\langle
x\rangle^{J}\langle\nabla\rangle^{J-1}\Phi_{1}\|_{2}.$ (C.8)

For the nonlinear part, we also only present the worst case:

$\displaystyle\|\langle x\rangle^{J}\partial_{t}^{J-2}(\phi\partial
n)|_{t=0}\|_{2}\lesssim$ $\displaystyle\|\langle
x\rangle^{J}(\partial_{t}^{\leq
J-2}\phi)(\partial_{t}^{\leq\frac{J}{2}}n)|_{t=0}\|_{2}+\|\langle
x\rangle^{J}(\partial_{t}^{\leq
J-1}n)(\partial_{t}^{\leq\frac{J-2}{2}}\phi)|_{t=0}\|_{2}$ (C.9)
$\displaystyle\lesssim$ $\displaystyle\|\langle x\rangle^{J}\partial_{t}^{\leq
J-2}\phi|_{t=0}\|_{2}\|\partial_{t}^{\leq\frac{J}{2}}n|_{t=0}\|_{\infty}+\|\langle
x\rangle^{J}\partial_{t}^{\leq\frac{J-2}{2}}\phi|_{t=0}\|_{\infty}\|\partial_{t}^{\leq{J-1}}n|_{t=0}\|_{2}.$
(C.10)

Then by induction and the standard $H^{2}(\mathbb{R}^{2})\hookrightarrow
L^{\infty}(\mathbb{R}^{2})$ embedding we have

$\displaystyle\|\langle x\rangle^{J}\partial_{t}^{\leq
J-2}\phi|_{t=0}\|_{2}\|\partial_{t}^{\leq\frac{J}{2}}n|_{t=0}\|_{\infty}$
(C.11) $\displaystyle\lesssim$ $\displaystyle(\|\langle
x\rangle^{J}\langle\nabla\rangle^{J-2}\Phi_{0}\|_{2}+\|\langle
x\rangle^{J}\langle\nabla\rangle^{J-3}\Phi_{1}\|_{2})\sum_{j=0}^{\frac{J}{2}+1}\|\langle
x\rangle^{j+1}\nabla^{j}(\nabla N_{0},N_{1})\|_{2}$ (C.12)

and

$\displaystyle\|\langle
x\rangle^{J}\partial_{t}^{\leq\frac{J-2}{2}}\phi|_{t=0}\|_{\infty}\|\partial_{t}^{\leq{J-1}}n|_{t=0}\|_{2}$
(C.13) $\displaystyle\lesssim$ $\displaystyle(\|\langle
x\rangle^{J}\langle\nabla\rangle^{\frac{J+2}{2}}\Phi_{0}\|_{2}+\|\langle
x\rangle^{J}\langle\nabla\rangle^{\frac{J}{2}}\Phi_{1}\|_{2})\sum_{j=0}^{J-2}\|\langle
x\rangle^{j+1}\nabla^{j}(\nabla N_{0},N_{1})\|_{2}.$ (C.14)

The other cases follow from similar arguments and we omit the details.

∎

## Acknowledgement

X. Cheng was partially supported by the Shanghai “Super Postdoc” Incentive
Plan (No. 2021014), the International Postdoctoral Exchange Fellowship Program
(No. YJ20220071) and China Postdoctoral Science Foundation (Grant No.
2022M710796, 2022T150139).

## References

  * [1] L. E. Abbrescia, Y. Chen. Global stability of some totally geodesic wave maps. J. Differential Equations, 284: 219-252, 2021.
  * [2] S. Alinhac. The null condition for quasilinear wave equations in two space dimensions I. Invent. Math., 145: 597-618, 2001.
  * [3] S. Alinhac. The null condition for quasilinear wave equations in two space dimensions II. Amer. J. Math., 123: 1071-1101, 2001.
  * [4] A. Bachelot. Problème de Cauchy global pour des systèmes de Dirac-Klein-Gordon. Ann. Inst. Henri Poincaré, 48: 387–422, 1988.
  * [5] D. Christodoulou. Global solutions of nonlinear hyperbolic equations for small initial data. Comm. Pure Appl. Math., 39: 267-282, 1986.
  * [6] X. Cheng, J. Xu. On the global well-posedness and scattering of the 3D Klein-Gordon-Zakharov system. Preprint, arXiv: 2210.13786.
  * [7] R. O. Dendy. Plasma Dynamics., Oxford University Press, Oxford, UK, 1990.
  * [8] S. Dong. Global solution to the Klein-Gordon-Zakharov equations with uniform energy bounds. SIAM J. Math. Anal., 54(1): 595–615, 2022.
  * [9] S. Dong. Asymptotic behavior of the solution to the Klein-Gordon-Zakharov model in dimension two, Comm. Math. Phys., 384: 587-607, 2021.
  * [10] S. Dong, P. LeFloch, Z. Lei. The top-order energy of quasilinear wave equations in two space dimensions is uniformly bounded., Fundamental Research, 2022.
  * [11] S. Dong, Y. Ma. Global Existence and Scattering of the Klein-Gordon-Zakharov System in Two Space Dimensions, to appear in Peking. Math. J., 2023.
  * [12] S. Dong, Z. Wyatt. Stability of a coupled wave-Klein-Gordon system with quadratic nonlinearities, J. Differential Equations, 269: 7470-7497, 2020.
  * [13] S. Duan, Y. Ma. Global Solutions of Wave-Klein–Gordon Systems in 2+1 Dimensional Space-Time with Strong Couplings in Divergence Form. SIAM J. Math. Anal., 54(3): 2691-2726, 2022.
  * [14] V. Georgiev. Global solution of the system of wave and Klein-Gordon equations. Math. Z., 203: 683-698, 1990.
  * [15] V. Georgiev. Decay estimates for the Klein-Gordon equation. Comm. Partial Differential Equations, 17(7-8): 1111–1139, 1992.
  * [16] Z. Guo, K. Nakanishi, S. Wang. Small energy scattering for the Klein-Gordon-Zakharov system with radial symmetry. Math. Res. Lett., 21: 733-755, 2014.
  * [17] Z. Guo, K. Nakanishi, S. Wang. Global dynamics below the ground state energy for the Klein-Gordon-Zakharov system in the 3D radial case. Comm. Partial Differential Equations, 39: 1158-1184, 2014.
  * [18] B. Guo, G. Yuan. Global smooth solution for the Klein-Gordon-Zakharov equations. J. Math. Phys., 36, 1995.
  * [19] Z. Hani, F. Pusateri, J. Shatah. Scattering for the Zakharov system in 3 dimensions. Comm. Math. Phys., 322: 731-753, 2013.
  * [20] L. Hörmander. Lectures on Nonlinear Hyperbolic Differential Equations Mathématiques & Applications (Berlin) [Mathematics & Applications], vol. 26, Springer-Verlag, Berlin, 1997.
  * [21] F. John. Blow-up for quasi-linear wave equations in three space dimensions. Comm. Pure Appl. Math., 34(1): 29-51, 1981.
  * [22] S. Katayama. Global existence for coupled systems of nonlinear wave and Klein-Gordon equations in three space dimensions. Math. Z., 270: 487-513, 2012.
  * [23] J. Krieger, Global regularity and singularity development for wave maps, Surveys in differential geometry 12 (1) (2007) 167–202.
  * [24] S. Klainerman. Global existence of small amplitude solutions to nonlinear Klein-Gordon equations in four spacetime dimensions. Comm. Pure Appl. Math., 38: 631-641, 1985.
  * [25] S. Klainerman, Q. Wang, S. Yang. Global solution for massive Maxwell -Klein-Gordon equations. Commun. Pure Appl. Math., 73: 63-109, 2020.
  * [26] S. Klainerman. The null condition and global existence to nonlinear wave equations. Nonlinear systems of partial differential equations in applied mathematics, Part 1 (Santa Fe, N.M., 1984), 293–326, Lectures in Appl. Math., 23, Amer. Math. Soc., Providence, RI, 1986.
  * [27] P. G. LeFloch, Y. Ma. The global nonlinear stability of Minkowski space for self-gravitating massive fields. The wave-Klein-Gordon model. Comm. Math. Phys., 346: 603-665, 2016.
  * [28] D. Li. Uniform estimates for 2D quasilinear wave. Advances in Mathematics, 428: 109157, 2023.
  * [29] D. Li. Global well-posedness of hedgehog solutions for the (3+1) Skyrme model. Duke Mathematical Journal, 170(7): 1377-1418, 2021.
  * [30] D. Li, Y. Wu. The Cauchy problem for the two dimensional Euler-Poisson system. J. Eur. Math. Soc. 16(10): 2211–2266, 2014.
  * [31] N. Masmoudi, K. Nakanishi. Energy convergence for singular limits of Zakharov type systems. Invent. Math., 172: 535–583, 2008.
  * [32] J. Metcalfe, C.D. Sogge. Global existence of null-form wave equations in exterior domains. Math. Z., 256: 521-549, 2007.
  * [33] J. Metcalfe, T. Rhoads. Long-time existence for systems of quasilinear wave equations. La Matematica, 1-48, 2023.
  * [34] T. Ozawa, K. Tsutaya, Y. Tsutsumi. Normal form and global solutions for the Klein-Gordon-Zakharov equations, Ann. Inst. Henri Poincaré, 12: 459-503, 1995.
  * [35] T. Ozawa, K. Tsutaya, Y. Tsutsumi. Well-posedness in energy space for the Cauchy problem of the Klein-Gordon-Zakharov equations with different propagation speeds in three space dimensions. Math. Ann., 313: 127-140, 1999.
  * [36] M. Psarelli. Asymptotic behavior of the solutions of Maxwell-Klein-Gordon field equations in 4-dimensional Minkowski space. Comm. Partial Differential Equations, 24: 223–272, 1999.
  * [37] M. Psarelli. Time decay of Maxwell-Klein-Gordon equations in 4-dimensional Minkowski space. Comm. Partial Differential Equations, 24: 273-282, 1999.
  * [38] J. Shatah. Normal forms and quadratic nonlinear Klein-Gordon equations. Comm. Pure Appl. Math., 38: 685-696, 1985.
  * [39] J. Shatah, M. Struwe, Geometric wave equations, Vol. Courant Lecture Notes in Mathematics, vol. 2, American Mathematical Society, Providence, RI, 1998.
  * [40] Q. Shi, S. Wang. Klein-Gordon-Zakharov system in energy space: Blow-up profile and subsonic limit. Math. Methods Appl. Sci., 42: 3211-3221, 2019.
  * [41] T.C. Sideris, S.-Y. Tu. Global existence for systems of nonlinear wave equations in 3D with multiple speeds. SIAM J. Math. Anal., 33: 477-488, 2001.
  * [42] T. H. R. Skyrme. A non-linear field theory, Proc. Roy. Soc. London Ser. A 260: 121-138, 1961.
  * [43] A. Stingo. Global existence of small amplitude solutions for a model quadratic quasi-linear coupled wave-Klein-Gordon system in two space dimension, with mildly decaying Cauchy data, Preprint, arXiv:1810.10235, to appear in Memoir AMS.
  * [44] B. Texier. Derivation of the Zakharov equations. Arch. Ration. Mech. Anal., 184: 121-183, 2007.
  * [45] J. Vilms. Totally geodesic maps, J. Differential Geom., 4: 73–79, 1970.
  * [46] C. Wang, X. Yu. Global existence of null-form wave equations on small asymptotically Euclidean manifolds. J. Funct. Anal., 266: 5676-5708, 2014.
  * [47] Q. Wang. An intrinsic hyperboloid approach for Einstein Klein-Gordon equations. J. Diff. Geom., 115: 27-109, 2020.
  * [48] S. Yang. Pointwise decay for semilinear wave equations in $\mathbb{R}^{1+3}$. J. Funct. Anal., 283(2), Paper No. 109486, 2022.
  * [49] V. E. Zakharov. Collapse of Langmuir waves. Sov. Phys. JETP, 35: 908-914, 1972.
  * [50] D. Zha. Global stability of solutions to two-dimension and one-dimension systems of semilinear wave equations. J. Funct. Anal., 282(1), Paper No. 1092219, 2022.