# The Drag Instability in a 2D Isothermal C-shock

Pin-Gao Gu  Institute of Astronomy & Astrophysics, Academia Sinica, Taipei
10617, Taiwan

###### Abstract

We extend the linear analysis of the drag instability in a 1D perpendicular
isothermal C-shock by Gu & Chen to 2D perpendicular and oblique C-shocks in
the typical environment of star-forming clouds. Simplified dispersion
relations are derived for the unstable modes. We find that the mode property
of the drag instability generally depends on the ratio of the transverse
(normal to the shock flow) to longitudinal (along the shock flow) wavenumber.
For the transversely large-scale mode, the growth rate and wave frequency of
the drag instability in a 2D shock resemble those in a 1D shock. For the
transversely small-scale mode, the drag instability is characterized by an
unstable mode coupled with an acoustic mode primarily along the transverse
direction. When the shock is perpendicular or less oblique, there exists a
slowly propagating mode, which can potentially grow into a nonlinear regime
and contribute to the maximum growth of the instability. In contrast, when the
shock is more oblique, this slowly propagating unstable mode disappears, and
the maximum growth of the drag instability is likely contributed from the
transversely large-scale mode (i.e., almost 1D mode). In all cases that we
consider, the magnitude of the density perturbations is significantly larger
than that of the velocity and magnetic field perturbations, implying that the
density enhancement governs the dynamics in the linear regime of the
instability. A few issues in the linear analysis, as well as the possible
astrophysical implications, are also briefly discussed.

##  1 Introduction

Stars form in the cold and dense cores of molecular clouds (e.g., Kennicutt &
Evans, 2012; Hennebelle & Inutsuka, 2019; Girichidis et al., 2020). Many
observations suggest that molecular clouds exhibit highly supersonic
turbulence (e.g., Elmegreen & Scalo, 2004; Ballesteros-Paredes et al., 2007;
Hennebelle & Falgarone, 2012). Additionally, cosmic rays can weakly ionize the
neutral gas to produce ions with a typical ionization fraction $\lesssim
10^{-6}$ (e.g., Draine et al, 1983; Tielens, 2005; Dalgarno, 2006; Indriolo &
McCall, 2012). Since the magnetic fields are normally assumed to be tightly
coupled to the ions, which in turn interact with the neutrals through the drag
force, the ions together with the magnetic fields can systematically drift
relative to the neutrals, a phenomenon known as ambipolar diffusion (Spitzer,
1956; Mestel & Spitzer, 1956; Shu, 1992; Zweibel, 2015).

Of particular interest in this study is the interplay between ambipolar
diffusion and shocks. While jump-type (J-type) shocks exhibit sharp supersonic
(or super-Alfvénic for magnetized shocks) discontinuities in physical
properties, in the case of nonideal magnetohydrodynamics (MHD), the different
signal speeds of the ions and neutrals along with the ambipolar diffusion
broaden the discontinuity, leading to continuous-type (C-type) shocks with a
smooth transition in physical quantities between the pre- and post-shock
regions (Draine et al, 1983; Draine & McKee, 1993). Observationally, various
efforts have been made to detect such features in turbulent molecular clouds
(e.g., Li & Houde, 2008; Hezareh et al., 2010, 2014; Xu & Li, 2016; Tang et
al., 2018), though most of them are indirect measurements and highly dependent
on the adopted dynamical and chemical models (e.g., Flower & Pineau Des
Forêts, 1998, 2010; Gusdorf et al., 2008; Lehmann & Wardle, 2016; Valdivia et
al., 2017, see also the Introduction of Gu & Chen 2020, hereafter GC (20)).
Theoretically, Chen & Ostriker (2012, 2014) analyzed the structures of
perpendicular and oblique C-shocks in exquisite detail in the typical
environment of star-forming clouds and extended their work to colliding flows
to explore the formation of cores and filaments in numerical simulations.

However, a C-shock is not necessarily a dynamically stable structure. Based on
the 1D background state of steady, plane-parallel C-shocks derived by Chen &
Ostriker (2012), GC (20) conducted a Wentzel-Kramers-Brillouin-Jeffreys (WKBJ)
analysis and confirmed the postulation of Gu et al. (2004) that the drag
instability can occur in a 1D isothermal perpendicular C-shock. In an
environment where the ambipolar diffusion is efficient and the
ionization–recombination equilibrium is nearly attained, the drag instability,
a plasma effect discovered by Gu et al. (2004), ensues from the ion-neutral
drag and is a local linear overstability phenomenon associated with an
exponentially growing mode of a propagating wave. On the other hand, it is
well known that C-shocks are also susceptible to the Wardle instability
(Wardle, 1990, 1991). However, the Wardle instability is suppressed by the
fast ionization–recombination process expected to pervade star-forming clouds
(e.g., Smith & Mac Low, 1997; Stone, 1997; Falle et al., 2009).

It was found by GC (20) that the growing wave mode of the drag instability can
only propagate downstream within a 1D shock and subsequently decay in the
post-shock region. The authors demonstrated that the growth of the drag
instability in a 1D perpendicular shock is limited by the short time span for
the unstable wave to stay within a C-shock. Consequently, the maximum growth
of an unstable wave is given by the maximum total growth (MTG), which is
defined by GC (20) as the maximum value of the total growth of an unstable
mode traveling over an entire shock width before it is damped in the post-
shock region. Specifically, MTG$=\exp\int_{\rm shock\
width}\Gamma_{grow}dx/v_{ph,x}$, where $\Gamma_{grow}$ is the growth rate and
$v_{ph,x}$ is the longitudinal phase velocity. The authors estimated the MTG
for such an unstable wave in typical environments of star-forming clouds. They
found that a stronger shock with a larger shock width favors a more
appreciable growth of the 1D drag instability. Most importantly, the linear
analysis suggests that the density enhancement induced by the drag instability
dominates over the magnetic field and velocity enhancements in the dynamical
growth of the unstable mode.

One of the most important questions for shocks in this context is whether a
shock instability can lead to fragmentation that subsequently undergoes
gravitational collapse and eventually induces star formation. The predominant
density enhancement of the drag instability provides a possible mechanism to
form supercritical clumps/cores within a C-shock in the typical environments
of star-forming clouds. As a preliminary activity, the analysis by GC (20)
focused on the basic behavior of the drag instability within the steady-state
profiles of C-shocks, which is 1D, linear, and non-self-gravitating. Hence,
how exactly this plasma instability plays a role in the observational evidence
of prestellar cores and filamentary networks in relation to magnetic field
morphology (e.g., André et al., 2014; Li et al., 2014) is not clear and has
yet to be addressed in terms of the current theoretical framework in 1D.
Obviously, the theory is still in its infancy and requires further development
toward more realistic physical applications.

Therefore, as a natural consequence of an ongoing effort, we extend the 1D
analysis of the drag instability by GC (20) and conduct a 2D linear analysis
for both perpendicular and oblique shocks in this study. The aim is to develop
a more comprehensive analytical work to provide useful information for probing
and characterizing the drag instability in future numerical simulations and
proceed to a more realistic case for future astrophysical applications.

The contents of the paper are structured as follows. In Section 2, we begin
with a linear analysis of the drag instability in a 2D perpendicular shock and
identify the mode with the maximal growth in the fiducial model of a C-shock.
Based on that experience, we proceed to the linear analysis of the drag
instability in an oblique shock in Section 3 and study the instability
properties in both the fiducial model and a model referred to as model V06 for
a stronger and wider shock for comparison. Finally, the summary and brief
discussions are presented in Section 4.

##  2 Linear analysis: isothermal perpendicular shocks

In cold molecular clouds and their substructures, the dynamical evolution of
ions and neutrals is governed by their individual continuity and momentum
equations, which include cosmic-ray ionization, ion–electron recombination in
the gas phase, mutual collisional drag force, the Lorentz force on ions, and
the pressure force with the isothermal equation of state. Additionally, the
evolution of magnetic fields is governed by the induction equation for ions.
The entire set of equations then reads as follows (e.g., Draine, 1980; Shu,
1992; Chen & Ostriker, 2012; GC, 20):

$\displaystyle{\partial\rho_{n}\over\partial t}+\nabla\cdot(\rho_{n}{\bf
v_{n}})=0,$ (1) $\displaystyle{\partial\rho_{i}\over\partial
t}+\nabla\cdot(\rho_{i}{\bf
v_{i}})=-\beta\rho_{i}^{2}+\xi_{\mathrm{CR}}\rho_{n},$ (2)
$\displaystyle\rho_{n}\left[{\partial{\bf v_{n}}\over\partial t}+({\bf
v_{n}}\cdot\nabla){\bf v_{n}}\right]+\nabla p_{n}={\bf f_{d}},$ (3)
$\displaystyle\rho_{i}\left[{\partial{\bf v_{i}}\over\partial t}+({\bf
v_{i}}\cdot\nabla){\bf v_{i}}\right]+\nabla p_{i}-{1\over
4\pi}(\nabla\times{\bf B})\times{\bf B}=-{\bf f_{d}},$ (4)
$\displaystyle{\partial{\bf B}\over\partial t}+\nabla\times({\bf B}\times{\bf
v_{i}})=0,$ (5)

where ${\bf v}$ is the velocity, $\rho$ is the density, ${\bf B}$ is the
magnetic field, $p=\rho c_{s}^{2}$ is the gas pressure when the isothermal
sound speed $c_{s}$ is 0.2 km/s at a temperature of $\sim 10$ K (Fukui &
Kawamura, 2010), and the subscripts $i$ and $n$ denote the ion and neutral
species, respectively. Note that the neutrals and ions are coupled by the
collisional drag force ${\bf f_{d}}\equiv\gamma\rho_{i}\rho_{n}{\bf
v_{d}}=\gamma\rho_{i}\rho_{n}({\bf v_{i}}-{\bf v_{n}})$, where $\gamma\approx
3.5\times 10^{13}$ cm3 s-1 g-1 is the drag force coefficient (Draine et al,
1983). The evolution of ion number density is controlled by the cosmic-ray
ionization rate $\xi_{\mathrm{CR}}$ and the ion recombination in the gas phase
$\beta$ (see e.g., Chen & Ostriker, 2012). We define the ionization parameter
$\chi_{i0}\equiv 10^{6}\sqrt{\xi_{\mathrm{CR}}(m_{n}/m_{i})/(\beta m_{i})}$.
as done by Chen & Ostriker (2012). We thus adopt $\beta\approx 10^{-7}$ cm3
s${}^{-1}/m_{i}$ and $\xi_{\mathrm{CR}}\approx 10^{-17}$ s-1 ($m_{i}/m_{n}$)
in this study (see, e.g., Shu, 1992; Tielens, 2005), where $m_{n}=2.3\times$
and $m_{i}=30\times$ the hydrogen mass are considered. Indeed, $\chi_{i0}=10$
falls in the typical range of $\chi_{i0}$ observed in star-forming regions
($\sim 1-20$; see, e.g., McKee et al. 2010).

###  2.1 Background States and Linearized Equations

Since there is no background structure parallel to the shock front in a 2D
plane-parallel, perpendicular shock, no fluid motion is normal to the shock
flow. Resultantly, the background states of a 2D perpendicular shock are the
same as those in a 1D perpendicular shock. Therefore, a detailed description
of the background states can be found in GC (20). For clarity, here we simply
summarize the shock model to smoothly connect to the linear analysis to be
presented later in the paper. We consider the gas flow toward the $+x$
direction across the shock with the magnetic field in the $y$-direction. The
equilibrium equations for a plane-parallel shock (i.e. $\partial/\partial
t=\partial/\partial y=0$) are derived from Equations(1)–(5), which are given
by Equations(6)–(10) in GC (20) to provide the background state of our linear
analysis in the shock frame. In these equilibrium equations, GC (20) also
employed a strong coupling approximation under which the ion-neutral drag is
balanced by the magnetic pressure gradient of the ions, i.e.,
$\gamma\rho_{n}L_{B}V_{d}=V_{A,i}^{2}$ where $L_{B}\equiv(-d\ln B/dx)^{-1}$
and $V_{A,i}$ is the Alfvén speed of the ions. Additionally, the equilibrium
between cosmic ionization and recombination is assumed; namely,
$\beta\rho_{i}^{2}=\xi_{\mathrm{CR}}\rho_{n}$ (see Chen & Ostriker (2012) for
justifications of such a choice).

Applying the zero-gradient boundary conditions ($d/dx=0$) far upstream and
downstream (i.e., no structures in the steady pre- and post-shock regions),
Chen & Ostriker (2012) derived the 1D structure equation of a C-shock.
Together with the pre-shock conditions described by $n_{0}$ (neutral number
density), $v_{0}=V_{i,0}=V_{n,0}$ (shock velocity), $B_{0}$, and $\chi_{i0}$
(here and throughout this paper, we use the subscript 0 to denote a physical
quantity in the pre-shock region), the field compression ratio $r_{B}\equiv
B/B_{0}=V_{i,0}/V_{i}$ as well as the neutral compression ratio
$r_{n}\equiv\rho_{n}/\rho_{n,0}=V_{n,0}/V_{n}$ can be solved, and $\rho_{i}$,
$\rho_{n}$, $B$, $V_{i}$, and $V_{n}$ can be subsequently obtained throughout
the shock width. Following GC (20), we place the shock front at $x=0$ pc for
convenience to plot and refer to shock properties as a function of $x$. In
this setup of the problem, the background drift velocity $V_{d}=V_{i}-V_{n}<0$
inside the C-shock (i.e. within the smooth shock transition). As in GC (20),
we adopt the 1D steady C-shock model shown in Figure 3 of Chen & Ostriker
(2012) as the fiducial model for the background state in the shock frame, with
the pre-shock parameters $n_{0}=500$ cm-3, $v_{0}=5$ km/s, $B_{0}=5\mu$G, and
$\chi_{i0}=10$. The left panel of Figure 1 shows $r_{n}/r_{B}$ in the fiducial
model, which is the same as Figure 1 in GC (20).

We now consider the perturbation eigenvector given by
$U(\omega,k_{x},k_{y})\equiv(\delta\rho_{i},\delta v_{x,i},\delta
v_{y,i},\delta B_{x},\delta B_{y},\delta\rho_{n},\delta v_{x,n},\delta
v_{y,n})^{T}$ multiplied by $\exp[\mathrm{i}(k_{x}x+k_{y}y+\omega t)]$ under
the WKBJ approximation, where $k_{x}$ is the longitudinal wavenumber, $k_{y}$
is the transverse wavenumber, and $\omega$ is the eigenvalue evaluated in the
shock frame. By substituting these perturbations and the background states
into the equations (1)–(5), the following linearized equations are obtained
(see Gu et al. (2004) and GC (20) in the 1D case):

$PU=\mathrm{i}\omega U,$ (6)

where

$\displaystyle
P=\left[\begin{array}[]{cccccccc}-\mathrm{i}k_{x}V_{i}-2\beta\rho_{i}&-\mathrm{i}k_{x}\rho_{i}&-\mathrm{i}k_{y}\rho_{i}&0&0&\xi_{\mathrm{CR}}&0&0\\\
-\mathrm{i}k_{x}\frac{c^{2}_{s}}{\rho_{i}}&-\mathrm{i}k_{x}V_{i}-\gamma\rho_{n}&0&\mathrm{i}k_{y}\frac{V^{2}_{A,i}}{B_{y}}&-\mathrm{i}k_{x}\frac{V^{2}_{A,i}}{B_{y}}&-\gamma
V_{d}&\gamma\rho_{n}&0\\\
-\mathrm{i}k_{y}\frac{c^{2}_{s}}{\rho_{i}}&0&-\mathrm{i}k_{x}V_{i}-\gamma\rho_{n}&0&0&0&0&\gamma\rho_{n}\\\
0&\mathrm{i}k_{y}B_{y}&0&-\mathrm{i}k_{x}V_{i}&0&0&0&0\\\
0&-\mathrm{i}k_{x}B_{y}&0&0&-\mathrm{i}k_{x}V_{i}&0&0&0\\\
0&0&0&0&0&-\mathrm{i}k_{x}V_{n}&-\mathrm{i}k_{x}\rho_{n}&-\mathrm{i}k_{y}\rho_{n}\\\
\gamma
V_{d}&\gamma\rho_{i}&0&0&0&-\mathrm{i}k_{x}\frac{c_{s}^{2}}{\rho_{n}}&-\mathrm{i}k_{x}V_{n}-\gamma\rho_{i}&0\\\
0&0&\gamma\rho_{i}&0&0&-\mathrm{i}k_{y}\frac{c_{s}^{2}}{\rho_{n}}&0&-\mathrm{i}k_{x}V_{n}-\gamma\rho_{i}\end{array}\right].$
(15)

Without the loss of generality, we adopt a constant wavenumber $k_{x}$ of
$1/0.015$ pc${}^{-1}(\equiv k_{fid}$), as used in GC (20). The choice of
$k_{x}$ is made to satisfy the WKBJ approximation, as demonstrated in the
right panel of Figure 1 where $k_{x}L_{p}$ and $k_{x}L_{B}$ are shown to be
all much larger than unity (also refer to Figure 1 in GC, 20).

Figure 1: The $r_{n}/r_{B}$ ratio of the background state throughout the
C-shock in the fiducial model (left panel) and a test for the WKBJ
approximation when $k_{x}=k_{fid}\equiv 1/0.015$ pc-1 (right panel). Here
$L_{p}\equiv|d\ln p/dx|^{-1}$ and $L_{B}$ are the scale heights for the gas
pressure and magnetic field of the background states, respectively.

Figure 2: Growth rate $\Gamma_{grow}$ (left panel) and its corresponding wave
frequency $\omega_{wave}$ (right panel) of the drag instability as a function
of $k_{y}/k_{x}$ in the four cases of $k_{x}$ at $x\approx 0.2$ pc. Note that
the wave frequency undergoes a steep change from a negative to a positive
value at $k_{y}/k_{x}\approx 16.84$ regardless of the value of $k_{x}$. The
wave frequency is not presented when the growth rate is zero, i.e., the green
curve in the range of $k_{y}/k_{x}\sim 1-10$ and the red curve in the range of
$k_{y}/k_{x}\lesssim 10$.

Figure 2 shows the growth rate $\Gamma_{grow}$ ($=$Im[$\omega]<0$, left panel)
and its corresponding wave frequency $\omega_{wave}$ ($=$Re[$\omega]$, right
panel) of the unstable wave mode as a function of $k_{y}/k_{x}$ in the shock
frame, at the location of $x\approx 0.2$ pc in the middle of the C-shock. It
is the only unstable mode among all eigenmodes computed from Equation(6). When
the longitudinal wavelengths are not exceptionally short (i.e.,
$k_{x}=k_{fid}$ and $5k_{fid}$ in the figure), the growth rate decreases
approximately with increasing $k_{y}/k_{x}$. However, when $k_{x}$ is as large
as $20k_{fid}$ (green curve) and $50k_{fid}$ (red curve), the growth rate
decreases and even drops to zero, except when $k_{y}/k_{x}$ is sufficiently
large in this particular example, which corresponds to a jump of the wave
frequency around $k_{y}/k_{x}\gtrsim 12.7$, as shown in the right panel of
Figure 2. We explore the underlying physics for the mode properties with
simplified dispersion relations in the following subsection.

###  2.2 Simplified Dispersion Relations

The eigenvalue problem for Equation(6) amounts to a complicated equation
containing a polynomial of $\omega$ to the eighth power. To make the physical
analysis of the unstable mode extractable, we attempt to numerically solve
Equation(6) by removing as many terms as possible in the matrix $P$ but still
obtain the same growth rate and wave frequency as those shown in Figure 2.
After this exercise, we realize that Equation(6) can be further reduced to

$\displaystyle(\Gamma_{i}+2\beta\rho_{i}){\delta\rho_{i}\over\rho_{i}}+\mathrm{i}k_{y}\delta
v_{i,y}-\xi_{CR}{\rho_{n}\over\rho_{i}}{\delta\rho_{n}\over\rho_{n}}=0,$ (16)
$\displaystyle(\Gamma_{i}+\gamma\rho_{n})\delta v_{i,y}-\gamma\rho_{n}\delta
v_{n,y}=0,$ (17)
$\displaystyle\Gamma_{n}{\delta\rho_{n}\over\rho_{n}}+\mathrm{i}k_{x}\delta
v_{n,x}+\mathrm{i}k_{y}\delta v_{n,y}=0,$ (18)
$\displaystyle-\gamma\rho_{i}V_{d}{\delta\rho_{i}\over\rho_{i}}+\mathrm{i}k_{x}c_{s}^{2}{\delta\rho_{n}\over\rho_{n}}+(\Gamma_{n}+\gamma\rho_{i})\delta
v_{n,x}=0,$ (19) $\displaystyle-\gamma\rho_{i}\delta
v_{i,y}+\mathrm{i}k_{y}c_{s}^{2}{\delta\rho_{n}\over\rho_{n}}+(\Gamma_{n}+\gamma\rho_{i})\delta
v_{n,y}=0,$ (20)

where $\Gamma_{i}=\mathrm{i}(\omega+k_{x}V_{i})$ is the eigenvalue in the
comoving frame of the ions and $\Gamma_{n}=\mathrm{i}(\omega+k_{x}V_{n})$ is
the eigenvalue in the comoving frame of the neutrals. Since the collisional
rate for an ion with the ambient neutrals $\gamma\rho_{n}$ is tremendously
large in a weakly ionized cloud, Equation(17) implies that $\delta
v_{i,y}\sim\delta v_{n,y}$, which allows Equation(20) to be further simplified
to

$\mathrm{i}k_{y}c_{s}^{2}{\delta\rho_{n}\over\rho_{n}}+\Gamma_{n}\delta
v_{n,y}\approx 0,$ (21)

and consequently, Equation (16) can be approximated to

$(\Gamma_{i}+2\beta\rho_{i}){\delta\rho_{i}\over\rho_{i}}+{k_{y}^{2}c_{s}^{2}\over\Gamma_{n}}{\delta\rho_{n}\over\rho_{n}}-\xi_{CR}{\rho_{n}\over\rho_{i}}{\delta\rho_{n}\over\rho_{n}}\approx
0.$ (22)

The absence of the induction equations for $\delta{\bf B}$ and the momentum
equation for $\delta v_{i,x}$ in the above reduced set of linear equations
(16)–(20) indicates that the disturbance of the ion flow along the shock
direction and that of magnetic fields plays a minor role in the dynamical
evolution of the 2D drag instability.

In a 2D perpendicular shock, the linearized equations should approach the 1D
result when $k_{y}\ll k_{x}$ (i.e., a transversely large-scale mode). The drag
instability can occur in a 1D steady C-shock (GC, 20). Specifically, for the
definite occurrence of the instability in 1D, the rate of the mode
$\Gamma_{n}$ observed in the comoving frame of the neutrals is considerably
smaller than both the recombination rate $2\beta\rho_{i}$ and the ion-neutral
drift rate across a distance of one wavelength (i.e., $k|V_{d}|$), whereas it
is considerably larger than the neutral collision rate with the ions (i.e.,
$\gamma\rho_{i}$) and the sound-crossing rate over one wavelength (i.e.,
$k_{x}c_{s}$). When the aforementioned conditions are satisfied, the 2D
linearized equations (i.e., Equations(18), (19), & (22)) are reduced to those
for a 1D perpendicular shock, which yield the simplified dispersion relation
(Gu et al., 2004; GC, 20)

$\Gamma_{n}\approx\pm{(1+\mathrm{i})\over
2}\sqrt{{k_{x}V_{d}}\gamma\rho_{i}}=\pm{(1+\mathrm{i})\over 2}\sqrt{k_{x}\over
L_{B}}V_{A,n},$ (23)

where the positive sign corresponds to a growing wave. As pointed out by GC
(20), the wave frequency in the shock frame $\omega_{wave}$ is dominated by
$k_{x}V_{n}$ instead of the imaginary part of $\Gamma_{n}$ due to the fast
background flow across the shock width downstream.

On the other hand, in the regime where $k_{y}\gg k_{x}$ (i.e., a transversely
small-scale mode) such that $k_{y}c_{s}$ is much larger than $2\beta\rho_{i}$,
$k_{x}V_{i}$, and $k_{x}V_{n}$, Equation(22) suggests that the terms
associated with ionization, recombination, and the background shock flows are
less significant, revealing a mode with the wave frequency
Re[$\omega$]$\approx k_{y}c_{s}$ much larger than Im[$\omega]$. Therefore,
Equation(22) suggests the relation
$\delta\rho_{i}/\rho_{i}\approx\delta\rho_{n}/\rho_{n}$. In this regime, the
in-phase relation between $\delta\rho_{i}$ and $\delta\rho_{n}$ is no longer
maintained by the ionization–recombination equilibrium in the regime of a
small $k_{y}/k_{x}$ but is a consequence of the fast acoustic wave along the
background magnetic fields,111In comparison,
$\delta\rho_{i}/\rho_{i}=(1/2)\delta\rho_{n}/\rho_{n}$ when the
ionization–recombination equilibrium is attained. The factor $1/2$ does not
appear in this regime when the slow mode is faster than the recombination
process. i.e., the so-called slow mode, because $c_{s}<V_{A,n}$ in our
fiducial model (e.g., Shu, 1992). This result, along with Equations(18), (19),
and (21), yield the following simplified dispersion relation:

$\Gamma_{n}\approx\pm\left(-\mathrm{i}k_{x}|V_{d}|\gamma\rho_{i}-k^{2}c_{s}^{2}\right)^{1/2}\approx\pm\left({1\over
2}{k_{x}|V_{d}|\gamma\rho_{i}\over
k_{y}^{2}c_{s}^{2}}+\mathrm{i}\right)k_{y}c_{s},$ (24)

where $k^{2}\equiv k_{x}^{2}+k_{y}^{2}$ and moreover, we have used
$\Gamma_{n}^{2}\gamma\rho_{i}\approx k_{y}^{2}c_{s}^{2}\gamma\rho_{i}$ from
Equation(22) and $k_{y}c_{s}\gg\sqrt{k_{x}|V_{d}|\gamma\rho_{i}}$ in deriving
the above equation. The inequality
$k_{y}c_{s}\gg\sqrt{k_{x}|V_{d}|\gamma\rho_{i}}$ holds in the regime where
$k_{y}c_{s}\gg 2\beta\rho_{i}$ is being considered here, along with
$2\beta\rho_{i}\gg\sqrt{k_{x}|V_{d}|\gamma\rho_{i}}$ in the fiducial model
(see the left panel of Figure 3 in GC, 20).222Recall that in a 1D C-shock,
$\sqrt{k_{x}|V_{d}|\gamma\rho_{i}}$ is about the growth rate of the drag
instability (see Equation 23), which is smaller than $2\beta\rho_{i}$ such
that the ionization–recombination equilibrium can be attained for the
perturbations to allow the instability to occur.

The positive sign of Equation(24) corresponds to a growing wave with the
growth rate related to the ion-neutral drag and the wave frequency
$k_{y}c_{s}$ associated with the slow mode. More specifically, in the comoving
frame of the neutrals, the phase velocity of the unstable wave is about
$-k_{y}c_{s}/k_{x}\hat{x}-c_{s}\hat{y}$ and the group velocity is
approximately $-(k_{x}/k_{y})c_{s}\hat{x}-c_{s}\hat{y}$. Hence, in this
regime, while both the phase and signal of the unstable wave propagate mainly
along the background field lines at the sound speed, the signal also
propagates upstream slowly at a speed much smaller than the sound speed. In
the shock frame, the phase velocity of the unstable wave is about
$[-k_{y}c_{s}/k_{x}+V_{n}]\hat{x}+[-c_{s}+(k_{x}/k_{y})V_{n}]\hat{y}\approx-
k_{y}c_{s}/k_{x}\hat{x}-c_{s}\hat{y}$.

Based on the above simplified dispersion relations for unstable modes on
different transverse scales, we may be able to interpret the results
illustrated in Figure 2. When $k_{y}/k_{x}\rightarrow 0$, both the growth rate
and wave frequency converge to the 1D result for the drag instability (GC,
20), as illustrated by the flat curves for $k_{y}/k_{x}<1$ in Figure 2.
However, as $k_{y}/k_{x}$ increases, a new property of the growing mode
associated with 2D emerges when the acoustic wave along the background field
line, i.e. the slow mode, becomes important. This happens when $c_{s}k_{y}\sim
k_{x}V_{n}$ as suggested by the transition of the dispersion relation from
Equation(23) to Equation(24). Physically, it occurs when the slow-mode rate
$c_{s}k_{y}$ introduced by the additional dimension along the background field
in the $y$-direction becomes comparable to the shock-crossing rate
$k_{x}V_{n}$ in the $x$-axis through the C-shock. Owing to this transition at
the large $k_{y}$, the jump of the wave frequency from a negative to a
positive value occurs at $k_{y}/k_{x}\sim
V_{n}/c_{s}\equiv\hat{k}_{jump}\approx 16.84$ at $x\approx 0.2$ pc, which
agrees with the right panel of Figure 2. While we show the value of
$\hat{k}_{jump}$ at one particular $x$, the range of $\hat{k}_{jump}$ within a
C-shock can be estimated by the pre-shock conditions. The $V_{n}$ decreases
from $v_{0}$ at the beginning of a shock to $v_{0}/r_{f}$ at the end of a
shock, where $r_{f}$ is the final compression ratio of the neutral density in
the post-shock region given approximately by $\sqrt{2}v_{0}/v_{A,n,0}$ (Chen &
Ostriker, 2012). Hence, $\hat{k}_{jump}$ decreases from $v_{0}/c_{s}$ to
$(v_{0}/r_{f})/c_{s}\approx B_{0}/(\sqrt{4\pi\rho_{n,0}}c_{s})$. In our
fiducial model, $\hat{k}_{jump}$ changes from about 25 to 3 across the shock
width. In the typical environments of star-forming clouds, $v_{0}\sim$ 1-6
km/s, $n_{0}\sim$ 100-1000 1/cm3, $B_{0}\sim 5-10$ $\mu$G, and the temperature
is about 10 K (e.g., see Table 1 in Chen & Ostriker, 2012). Therefore,
$\hat{k}_{jump}$ can decrease from $\sim$ 10-30 at the beginning of a C-shock
to $\sim$ 1.5-5 at the end of a C-shock.

Figure 3: Blowup of the region where $k_{y}/k_{x}>20$ in Figure 2 to compare
the exact behavior of the growing mode with the approximate result from
Equation(24). The curves of $k_{y}c_{s}$ corresponding to various $k_{x}$ are
also plotted on the right panel for comparison.

As $k_{y}/k_{x}$ further increases to an even larger value, the behavior of
the growing mode follows the dispersion relation described by Equation(24).
This is demonstrated in Figure 3, which zooms in to the region where
$k_{y}/k_{x}>20$ in Figure 2. The dashed lines present the growth rate
Re[$\Gamma_{n}$] (left panel) and the wave frequency
Im[$\Gamma_{n}-k_{x}V_{n}$] (right panel) calculated from Equation(24). As
$k_{y}/k_{x}$ increases, the exact growth rates of the unstable mode for
various $k_{x}$ computed from Equation(6) converge to Re[$\Gamma_{n}$] while
the exact wave frequencies $\omega_{wave}$ overlap with
Im[$\Gamma_{n}-kV_{n}$] in the entire zoomed-in region. The slow-mode
frequency $k_{y}c_{s}$, as illustrated by the dashed-dotted lines, is also
plotted for each value of $k_{x}$. As $k_{y}/k_{x}$ increases, the wave
frequency $\omega_{wave}$ in all cases converges toward $k_{y}c_{s}$,
confirming that the slow mode emerges and predominates on small transverse
scales for the drag instability. Note that the positive wave frequency for
$k_{y}/k_{x}>\hat{k}_{jump}$ indicates a growing wave propagating upstream
within the 2D perpendicular C-shock, in stark contrast to the growing mode for
$k_{y}/k_{x}<\hat{k}_{jump}$, which is more dynamically dominated by the wave
convected with the downstream flow across the shock width, similar to the
result for a 1D C-shock.

In the case of a moderate value of $k_{x}$ (i.e., $k_{x}=k_{fid}$ and
$5k_{fid}$), the growth rate decreases almost monotonically with increasing
$k_{y}$ when $k_{y}/k_{x}\gtrsim\hat{k}_{jump}$, in an approximate agreement
with the trend described by Equation(24). On the other hand, in the regime
where $k_{y}/k_{x}<1$, the growth rate is small or nearly zero when $k_{x}$ is
large enough (i.e., $k_{x}=50$ and $100k_{fid}$ in Figure 2) for the small-
scale gas pressure perturbation along the shock direction to suppress the
growth (Gu et al., 2004; GC, 20). However, a growing mode exists when
$k_{y}/k_{x}\gtrsim 12.7$. Its growth rate rises significantly during the wave
frequency transition and then declines, according to Equation(24), as
$k_{y}/k_{x}$ increases to a large value. We are unable to derive a simple
dispersion relation to describe the rise of the growth rate during the wave
frequency transition. It is because more terms in Equations(16)–(20) become
comparably important and thus cannot be neglected as a proper approximation
for deriving Equation(24). Nevertheless, this complexity of the distinct
behaviors of the growth rate at the wave frequency transition can be unraveled
by studying the phase differences between perturbed quantities, which will be
presented in the next subsection.

###  2.3 An auxiliary analysis based on phase differences between
perturbations

We illustrate the phase differences in Figure 4 corresponding to the growth
rate and wave frequency shown in Figure 2. Figure 4 includes the phase
difference between the density perturbations, as well as those between the
neutral velocity and density perturbations. The flat part of the curves in the
range of $k_{y}/k_{x}\lesssim 1$ corresponds to the flat part of the curves in
Figure 2, resembling the 1D drag instability for the 2D unstable modes. In
contrast to the cases for the small $k_{x}$ (i.e., blue and orange curves),
the phase difference between $\delta\rho_{i}$ and $\delta\rho_{n}$ in the
cases for the large $k_{x}$ (green and red curves in panel (a) of Figure 4)
exhibit a dramatic change when $k_{y}/k_{x}\approx\hat{k}_{jump}=16.8$,
corresponding to the jump of the wave frequency. Since the wave frequency of
the mode $\omega_{wave}$ is small around this jump transition, the Doppler-
shift frequency $k_{x}V_{i}$ with a large $k_{x}$ becomes nonnegligible
compared to the ionization rate in the continuity equation of the ions.
Consequently, the phase difference between the ion and neutral density
perturbations does not stay small by the ionization–recombination equilibrium
but becomes large. The phase difference decreases steeply with increasing
$k_{y}/k_{x}$ because the rate of the mode is predominated quickly by the slow
mode $\approx k_{y}c_{s}$, resulting in the small phase difference between the
ion and neutral density perturbations for a large $k_{y}/k_{x}$. It is in
accordance with the result that
$\delta\rho_{i}/\rho_{i}\approx\delta\rho_{n}/\rho_{n}$ when deriving
Equation(24).

The drastic change in $\phi_{\delta\rho_{i}}-\phi_{\delta\rho_{n}}$ also
causes the steep change in $\phi_{\delta v_{n,x}}-\phi_{\delta\rho_{n}}$ from
$\pi$ to $\approx 0.8\pi$ around the frequency transition (see the green and
red curves in panel (b) of Figure 4) through the linearized continuity and
momentum equations, leading to a bump in the growth rate of the drag
instability near the frequency transition, as shown in the left panel of
Figure 2.333Recall that the drag instability in a 1D case is driven by the
canonical phase difference $3\pi/4$ between $\phi_{\delta v_{n,x}}$ and
$\phi_{\delta\rho_{n}}$ (Gu et al., 2004; GC, 20). Since an acoustic mode can
be characterized by the phase difference $\pi$ between density and velocity
perturbations, the nearly out-of-phase difference between $\delta v_{n,y}$ and
$\delta\rho_{n}$ shown in panel (c) of Figure 4 for $k_{y}/k_{x}\gtrsim 16.8$
reconfirms the right panel of Figure 3, i.e., the emergence of the slow mode
when $k_{y}/k_{x}\gtrsim\hat{k}_{jump}$. It also agrees with Equation(21) on
the relation $\Gamma_{n}\sim\mathrm{i}k_{y}c_{s}$ for
$k_{y}/k_{x}>\hat{k}_{jump}$, which leads to
$\delta\rho_{n}/\rho_{n}\sim-\delta v_{n,y}/c_{s}$, i.e., out of phase between
$\delta\rho_{n}$ and $\delta v_{n,y}$ for a slow mode.

Figure 4: Phase differences of unstable modes between $\delta\rho_{i}$ and
$\delta\rho_{n}$ (panel (a)), $\delta v_{n,x}$ and $\delta\rho_{n}$ (panel
(b)), and $\delta v_{n,y}$ and $\delta\rho_{n}$ (panel (c)) as a function of
$k_{y}/k_{x}$ at the location of $x\approx 0.2$ pc in the fiducial model of
the perpendicular C-shock. Note that there is a steep change in the green and
red curves during the wavelength transition at
$k_{y}/k_{x}\approx\hat{k}_{jump}=16.84$ in panels (a) and (b). The phase
differences are not presented when the growth rates are zero for some values
of $k_{y}/k_{x}$ (refer to the left panel of Figure 2).

###  2.4 Substantial growth of slowly propagating waves

It was shown by GC (20) that the growth of the drag instability in a 1D
perpendicular shock is limited by the short time span for an unstable wave to
stay within a C-shock as the wave is convected downstream by a fast flow
across the shock; specifically, the local growth rate $\Gamma_{grow}$ is much
smaller than the wave frequency in the shock frame $\omega_{wave}$.
Consequently, a stronger shock with a larger shock width favors a more
appreciable growth of the 1D drag instability in the typical environments of
star-forming clouds investigated by the authors. However, for a 2D
perpendicular shock, as has been studied in the preceding subsection, the
transition of the wave frequency from a negative to a positive value allows
for the presence of a slowly traveling mode with an arbitrarily small wave
frequency in the shock frame when $k_{y}/k_{x}\sim\hat{k}_{jump}$, thereby
allowing a tremendous amount of time for the growth of the drag instability
within a shock width in this transition regime.

We investigate this expectation by simply considering a mode with a small wave
frequency $\omega_{wave}=-10^{-14}$ 1/s $\ll$ the flow crossing time over one
longitudinal wavelength $k_{x}V_{n}$. According to the eigenvalue problem, a
family of unstable modes should exist with the proper combination of the
wavenumbers $k_{x}$ and $k_{y}$ as a function of $x$. For the purpose of
computational convenience without loss of generality, we first consider the
uniform $k_{x}$ given by the fiducial wavenumber $k_{x}=k_{fid}=$ 1/(0.015
pc), while allowing the $k_{y}$ of the mode to vary with $x$. The results are
shown in Figure 5. The resulting growth rate $\Gamma_{grow}$ is indeed larger
than $\omega_{wave}$ in most of the region within the C-shock, as illustrated
in the left panel of Figure 5. The corresponding profile of $k_{y}/k_{x}$
decreases from 25 to 3 with increasing $x$, in agreement with the transition
wavenumber ratio $k_{y}/k_{x}\approx\hat{k}_{jump}=V_{n}/c_{s}$, where $V_{n}$
decreases across the shock width due to density compression, as has been
explained in Section 2.2. Next, we consider the uniform $k_{y}$ given by
$15k_{fid}$ (i.e. a transversely small-scale mode) while allowing the $k_{x}$
of the mode to vary with $x$. The results of the growth rate and wave
frequency are considerably similar to those shown in Figure 5, affirming the
slow propagation of the unstable modes with $k_{y}/k_{x}\approx\hat{k}_{jump}$
in a 2D perpendicular C-shock. We also study the “counterpart” mode with the
small positive wave frequency $\omega_{wave}=10^{-14}$ 1/s for the same case
of either the uniform $k_{x}$ or uniform $k_{y}$. We find that the profiles of
the growth rate and $k_{y}/k_{x}$ are almost identical to those presented in
Figure 5, except that the unstable wave slowly propagates upstream rather than
downstream due to the sign change of $\omega_{wave}$.

As the unstable wave propagates slowly at the phase velocity $v_{ph,x}$, the
unstable mode is nearly comoving with the shock during the cloud lifetime,
which is typically about tens of millions of years (Engargiola et al., 2003;
Blitz et al., 2007; Kawamura et al., 2009; Murray, 2011; Miura et al., 2012;
Meidt et al., 2015; Jeffreson & Kruijssen, 2018). Consequently, the MTG is no
longer a proper measure for the growth of the slowly propagating modes in the
shock frame. For instance, it takes about 82 million years for the unstable
mode that we have considered with $\omega_{wave}=-10^{-14}$ 1/s and uniform
$k_{x}=k_{fid}$ to travel across the entire shock width,444The wave-crossing
time through the shock width is given by $\int_{\rm shock\
width}dx/|v_{ph,x}|=\int dx/(|\omega_{wave}|/k_{x})$. which is $\gtrsim$ the
cloud lifetime. It is evident from Figure 2 that there are modes that are
almost/exactly stationary in the shock frame (i.e., $\omega_{wave}\approx 0$)
and have $\Gamma_{grow}\sim 10^{-13}$ 1/s $\sim 1/0.3$ 1/Myr. Hence, the
growth of these modes is approximately given by $\exp(\Gamma_{grow}t)$, which
can be substantially larger than unity if $t$ is some fraction of the cloud
lifetime. In contrast, the maximum growth of the unstable mode in the 1D case
is limited by the shock width and thus is given by MTG, which is merely about
9.9 in the fiducial model (GC, 20). The 2D mode ought to grow to a nonlinear
phase according to the expectation from the linear theory.

Figure 5: Growth rate $\Gamma_{grow}$ as a function of $x$ for the unstable
mode with the wave frequency $\omega_{wave}=-10^{-14}$ 1/s and $1/k_{x}=0.015$
pc (left panel) and the profile of the corresponding wavenumber ratio
$k_{y}/k_{x}$ (right panel) for the fiducial model of the perpendicular
C-shock.

Figure 6 shows the amplitude of perturbations for the unstable mode presented
in Figure 5. The perturbed quantities are normalized by the $y$-component of
the magnetic field perturbation $|\delta B_{y}|/B$. It is evident from the
figure that within the shock width where the unstable mode exists, the density
perturbations $\delta\rho_{n}/\rho_{n}$ and $\delta\rho_{i}/\rho_{i}$ dominate
over the others, implying that the growth of density perturbations plays a
decisive role in the dynamical evolution of the instability. Although this
outcome of promoting the local density growth of a wave turns out to be the
same as that in the 1D case (GC, 20), the drag instability is coupled with the
slow mode in the regime $k_{y}/k_{x}>\hat{k}_{jump}$ for a 2D perpendicular
shock, as described in the previous section. This can also be realized from
Figure 6, which shows $|\delta v_{n,y}|\gg|\delta v_{n,x}|$ and $|\delta
v_{i,y}|\gg|\delta v_{i,x}|$ as a result of the presence of the slow mode
propagating along the background magnetic fields in the $y$-direction.

It is worth noting that $\delta B_{x}$ and $\delta B_{y}$ are generated from
$B$ due to the gradients of the longitudinal motion $\partial_{y}(\delta
v_{i,x})$ (i.e., magnetic wiggling) and $\partial_{x}(\delta v_{i,x})$ (i.e.,
magnetic compression), respectively, through the induction equation such that
$|\delta B_{x}|/|\delta B_{y}|=k_{y}/k_{x}$, in agreement with the ratios
$|\delta B_{x}|/|\delta B_{y}|$ in Figure 6 and $k_{y}/k_{x}$ in the right
panel of Figure 5. However, given their small amplitudes, magnetic field
perturbations simply grow passively through the induction equation and do not
play any significant role in the dynamics of the drag instability, confirming
that the absence of the induction equation in the reduced set of linearized
equations from Equations(16)–(20) still leads to the same growth rate and wave
frequency of the drag instability. Finally, as a reminder, in the transversely
large-scale regime where $k_{y}/k_{x}\ll V_{n}/c_{s}=\hat{k}_{jump}$ , the
property of the drag instability in 2D resembles that in 1D, according to the
analysis in the previous subsection. Therefore, the density growth also
governs the dynamics of the drag instability for the transversely large-scale
mode (GC, 20).

Figure 6: Amplitude of perturbations normalized by the $y$-component of the
magnetic perturbation $|\delta B_{y}|/B$ for the unstable mode presented in
Figure 5.

###  2.5 Recap

When $k_{y}/k_{x}\lesssim V_{n}/c_{s}\equiv\hat{k}_{jump}$ (i.e., transversely
large-scale modes), the drag instability behaves similarly between 1D and 2D
perpendicular C-shocks. When $k_{y}/k_{x}\gtrsim\hat{k}_{jump}$ (i.e.,
transversely small-scale modes), a new property of the growing mode emerges
for 2D perpendicular C-shocks even for a large $k_{x}$, resulting from the
ion-neutral drag coupled with the slow mode along the background magnetic
field. When $k_{y}/k_{x}\sim\hat{k}_{jump}$, the sign transition of wave
frequency occurs. This leads to an unstable mode of a small wave frequency and
thus enables a slowly propagating mode to grow substantially within a C-shock,
unlimited by the shock width. The linear result suggests that the density
growth dominates the evolution of the perturbation growth driven by the drag
instability for the modes on both transversely large and small scales.

##  3 Linear analysis: isothermal oblique shocks

Figure 7: Oblique angle of $\bf B$ (solid curves), $\bf V_{i}$ (dashed
curves), and $\bf V_{n}$ (dotted curves) relative to the direction of the pre-
shock inflow (i.e., the $+x$ direction) as a function of $x$ for the fiducial
model of oblique C-shocks. The cases for three different initial oblique
angles $\theta_{0}$ (70∘, 45∘, 20∘) are presented. These results provide the
background states of oblique C-shocks for the linear analysis of the drag
instability.

###  3.1 Background States and Linearized Equations

In this section, we study the local stability of a 2D steady oblique shock in
which magnetic fields are not normal to the shock flow. We adopt the
background states constructed by Chen & Ostriker (2012). That is, the pre-
shock flow is still along the $x$-direction, the shock front in the $y$-$z$
plane, and the pre-shock magnetic field in the $x$-$y$ plane (${\bf
B_{0}}=B_{x,0}{\hat{x}}+B_{y,0}{\hat{y}}$), at an initial oblique angle
$\theta_{0}$ to the inflow ($B_{y,0}/B_{x,0}=\tan\theta_{0}$). For a steady,
parallel-plane oblique shock (i.e.,
$\partial_{t}=\partial_{y}=\partial_{z}=0$), Chen & Ostriker (2012) derived
the equilibrium equations for $r_{n}\equiv\rho_{n}/\rho_{n,0}=v_{0}/V_{n,x}$,
$r_{B}\equiv B_{y}/B_{y,0}$ ($B_{x}$ is constant because $\nabla\cdot{\bf
B}=0$), and $r_{ix}\equiv v_{0}/V_{i,x}$ from Equations(1)-(5) based on the
strong coupling approximation and ionization–recombination equilibrium. Note
that $r_{B}$ is not equal to $r_{ix}$ for an oblique shock. By solving the
equilibrium equations with a set of pre-shock conditions, the background
states including $\rho_{i}$, $\rho_{n}$, $V_{i,x}$, $V_{n,x}$, $V_{i,y}$,
$V_{n,y}$, and $B_{y}$ can be subsequently obtained as a function of $x$
within the shock and post-shock regions. Since the field line is not
perpendicular to the incoming shock flow, both the field line and gas inflow
will continuously change their directions relative to their initial direction
as they move across the shock width until they reach the post-shock region
(see Figure 7 in the case of the pre-shock conditions given by the fiducial
model). As expected, $\bf B$ and $\bf V_{i}$ are first tilted toward the shock
front (i.e. a large oblique angle shown in Figure 7) due to the earlier
compression of $\rho_{i}$, followed by the subsequent tilt of $\bf V_{n}$
toward the shock front due to the later compression of $\rho_{n}$ within the
shock. Finally, ${\bf V_{i}}={\bf V_{n}}$ in the post-shock region as in the
pre-shock region, but with a final nonzero oblique angle. The shock width is
given by the region between the pre- and post-shock. It is evident from Figure
7 that as $\theta_{0}$ decreases, the compression ratio increases and shock
width decreases. Besides, $\bf B$ in the three cases of $\theta_{0}$ shown in
Figure 7 all quickly becomes more or less normal to the shock flow by shock
compression. All of these results are anticipated from the analysis in Chen &
Ostriker (2012).

With the background states, the linearized equations are given by

$OU=\mathrm{i}\omega U,$ (25)

where

$\displaystyle
O=\left[\begin{array}[]{cccc}-\mathrm{i}k_{x}V_{i,x}-\mathrm{i}k_{y}V_{i,y}-2\beta\rho_{i}&-\mathrm{i}k_{x}\rho_{i}&-\mathrm{i}k_{y}\rho_{i}&0\\\
-\mathrm{i}k_{x}\frac{c^{2}_{s}}{\rho_{i}}&-\mathrm{i}k_{x}V_{i,x}-\mathrm{i}k_{y}V_{i,y}-\gamma\rho_{n}&0&\mathrm{i}k_{y}\frac{V^{2}_{A,i}}{B_{y}}\\\
-\mathrm{i}k_{y}\frac{c^{2}_{s}}{\rho_{i}}&0&-\mathrm{i}k_{x}V_{i,x}-\mathrm{i}k_{y}V_{i,y}-\gamma\rho_{n}&-\mathrm{i}k_{y}{B_{x}\over
B_{y}}{V_{A,i}^{2}\over B_{y}}\\\
0&\mathrm{i}k_{y}B_{y}&-\mathrm{i}k_{y}B_{x}&-\mathrm{i}k_{x}V_{i,x}-\mathrm{i}k_{y}V_{i,y}\\\
0&-\mathrm{i}k_{x}B_{y}&\mathrm{i}k_{x}B_{x}&0\\\ 0&0&0&0\\\ \gamma
V_{d,x}&\gamma\rho_{i}&0&0\\\ \gamma
V_{d,y}&0&\gamma\rho_{i}&0\end{array}\right.$ (34)

$\displaystyle\left.\begin{array}[]{cccc}0&\xi_{\mathrm{CR}}&0&0\\\
-\mathrm{i}k_{x}\frac{V^{2}_{A,i}}{B_{y}}&-\gamma V_{d,x}&\gamma\rho_{n}&0\\\
\mathrm{i}k_{x}{B_{x}\over B_{y}}{V_{A,i}^{2}\over B_{y}}&-\gamma
V_{d,y}&0&\gamma\rho_{n}\\\ 0&0&0&0\\\
-\mathrm{i}k_{x}V_{i,x}-\mathrm{i}k_{y}V_{i,y}&0&0&0\\\
0&-\mathrm{i}k_{x}V_{n,x}-\mathrm{i}k_{y}V_{n,y}&-\mathrm{i}k_{x}\rho_{n}&-\mathrm{i}k_{y}\rho_{n}\\\
0&-\mathrm{i}k_{x}\frac{c_{s}^{2}}{\rho_{n}}&-\mathrm{i}k_{x}V_{n,x}-\mathrm{i}k_{y}V_{n,y}-\gamma\rho_{i}&0\\\
0&-\mathrm{i}k_{y}\frac{c_{s}^{2}}{\rho_{n}}&0&-\mathrm{i}k_{x}V_{n,x}-\mathrm{i}k_{y}V_{n,y}-\gamma\rho_{i}\end{array}\right],$
(43)

where $V_{A,i}$ is still defined as $B_{y}^{2}/4\pi\rho_{i}$ in terms of the
field component parallel to the shock front, which follows the same definition
in the case of the 2D perpendicular shock. When the $x$-component of the pre-
shock $B$ field is zero, $B_{x}=0$, $V_{i,y}=0$, and $V_{n,y}=0$ (thus
$V_{d,y}=0$) throughout the shock. Thus, the matrix $O$ in the above equation
for an oblique shock is reduced to the matrix $P$ in equation(6) for a
perpendicular shock.

Figure 8: Growth rate $\Gamma_{grow}$ (left panel) and its corresponding wave
frequency in the shock frame $\omega_{wave}$ (right panel) of the drag
instability in the fiducial model for a 2D oblique shock with
$\theta_{0}=45^{\circ}$, plotted as a function of $k_{y}/k_{x}$ in the four
cases of $k_{x}$ at $x\approx 0.2$ pc. The wave frequency is not presented
when the growth rate is zero.

We solve the above eigenvalue problem for the fiducial model of the C-shock
with a pre-shock oblique angle $\theta_{0}<90^{\circ}$. As in the
perpendicular shock, we find that there is one unstable mode among the eight
eigenmodes within a C-shock. Figure 8 shows the growth rate $\Gamma_{grow}$
and wave frequency in the shock frame $\omega_{wave}$ of the unstable mode as
a function of $k_{y}/k_{x}$ in the case of $\theta_{0}=45^{\circ}$ with four
values of $k_{x}$ at $x\approx 0.2$ pc, approximately in the middle of the
oblique C-shock. The figure exhibits two clear general features. First of all,
the growth rates and wave frequencies do not vary significantly but have
values similar to the case of $k_{y}=0$ in the range of $k_{y}/k_{x}<1$, as
shown by the flat part of the curves. The second feature is that regardless of
the difference in $k_{x}$, the growth rates of the unstable mode all drop to
zero when $k_{y}/k_{x}\approx 9.3$, and the corresponding wave frequencies
exhibit a discontinuity at the same wavenumber ratio, which turns out to be
exactly equal to $|V_{d,x}|/V_{d,y}$ at $x\approx 0.2$ pc in the fiducial
model. We investigate the reason in the next subsections.

###  3.2 Simplified Dispersion Relation

Since $V_{n,y}$ and $V_{i,y}$ other than $V_{n,x}$ and $V_{i,x}$ are present
in an oblique shock, the Doppler-shift rates in the comoving frame of the ions
$\Gamma_{i}$ and neutrals $\Gamma_{n}$ are expressed by
$\mathrm{i}(\omega+k_{x}V_{i,x}+k_{y}V_{i,y})$ and
$\mathrm{i}(\omega+k_{x}V_{n,x}+k_{y}V_{n,y})$, respectively. Following the
same approach in Section 2.2 for perpendicular shocks, we focus on the
particular mode with $\Gamma_{n}$ smaller than $2\beta\rho_{i}$ and $kV_{d}$
but larger than $\gamma\rho_{i}$. The linear equations expressed in
Equation(25) may be reduced to the same set of Equations(18), (19), (21), and
(22), except that the additional drag term $\gamma V_{d,y}\delta\rho_{i}$
appears in the $y$-component of the momentum equation for the neutrals, namely
(cf. Equation(21)):

$-\gamma\rho_{i}V_{d,y}{\delta\rho_{i}\over\rho_{i}}+\mathrm{i}k_{y}c_{s}^{2}{\delta\rho_{n}\over\rho_{n}}+\Gamma_{n}\delta
v_{n,y}\approx 0.$ (44)

Moreover, as the recombination rate $2\beta\rho_{i}>\Gamma_{n}$ for this mode,
Equation(22) is reduced to
$\delta\rho_{i}/\rho_{i}\approx(1/2)\delta\rho_{n}/\rho_{n}$ due to the
ionization equilibrium. The resulting dispersion reads

$\Gamma_{n}^{2}=-\mathrm{i}{\gamma\rho_{i}\over
4}(k_{x}V_{d,x}+k_{y}V_{d,y})-k^{2}c_{s}^{2}.$ (45)

In the regime where $k_{y}/k_{x}\ll 1$ and $\Gamma_{n}>k_{x}c_{s}$, the above
equation is reduced to Equation(23), the typical growth rate and wave
frequency for the 1D drag instability (Gu et al., 2004; GC, 20), but with a
background state for an oblique shock. As in the case of a 2D perpendicular
shock discussed in the preceding section, this mode behavior of an oblique
shock is consistent with the flat part of the curves for both growth rate and
wave frequency in Figure 8. As $k_{y}/k_{x}$ increases to the moderate value
$|V_{d,x}|/V_{d,y}\equiv\hat{k}_{discon}=9.3$, Figure 8 shows that the growth
rates decline and become zero around $\hat{k}_{discon}$. Besides, the
corresponding wave frequencies decrease gradually with increasing
$k_{y}/k_{x}$ and then abruptly drop around $\hat{k}_{discon}$. This behavior
of the growth rate and wave frequency can be understood in terms of
Equation(45), which admits a special solution Re$[\Gamma_{n}]=0$ when
$k_{x}V_{d,x}+k_{y}V_{d,y}=0$. The dynamics attributed to the absence of the
growing mode at $k_{y}/k_{x}=\hat{k}_{discon}$ is that the drag force in the
$x$-direction (i.e., $\gamma V_{d,x}\delta\rho_{i}$) acts out of phase with
the drag force in the $y$-direction for the neutrals (i.e., $\gamma
V_{d,y}\delta\rho_{i}$), therefore suppressing the density enhancement in the
neutral continuity equation and thus quenching the drag instability.

To explain the property of the wave frequency shown in the right panel of
Figure 8, we investigate the behavior of Equation(45) near
$k_{y}/k_{x}=\hat{k}_{discon}$. Because $kc_{x}\gg
k_{x}V_{d,x}+k_{y}V_{d,y}\approx 0$ around this wavenumber ratio, the
dispersion relation of the unstable mode described by Equation(45) is
simplified to

$\Gamma_{n}=\mathrm{i}\omega+\mathrm{i}k_{x}V_{n,x}+\mathrm{i}k_{y}V_{n,y}\approx{\gamma\rho_{i}\over
4}{|k_{x}V_{d,x}+k_{y}V_{d,y}|\over kc_{s}}\pm\mathrm{i}kc_{s},$ (46)

where the plus (minus) sign of the imaginary part of the above expression
(i.e., the wave frequency in the comoving frame of the neutrals) is taken when
$k_{x}V_{d,x}+k_{y}V_{d,y}<0$ ($>0$). Due to this sign change, the wave
frequency in the shock frame $\omega_{wave}$ (=Re[$\omega$]) at
$k_{y}/k_{x}=k_{discon}$ exhibits a discontinuity, with the difference given
by $\approx 2kc_{s}$. As $k_{y}/k_{x}$ increases beyond $\hat{k}_{discon}$ for
a given $k_{x}$, Re[$\omega$] ($=-k_{x}V_{n,x}-k_{y}V_{n,y}-kc_{s}$) increases
more negatively; approximately speaking, Re[$\omega$]$\approx-
k_{y}(V_{n,y}+c_{s})\propto-k_{y}$, in close agreement with the right panel of
Figure 8.

###  3.3 An auxiliary analysis based on phase differences between
perturbations

Although we attempt to obtain the simplified dispersion relations for guiding
us in understanding the key features of the growth rate and wave frequency of
an unstable wave, not all of the features can be explained. Figure 8 shows
that the unstable mode is suppressed when $k_{y}/k_{x}$ is sufficiently large,
which is unable to be described by Equation(46), because some neglected terms
in deriving the simplified dispersion relations can become comparably
important for a large $k_{y}/k_{x}$. To gain a more comprehensive insight into
the 2D drag instability, we study the phase differences between perturbation
quantities of the unstable mode as in the preceding section to complement the
limited analysis from the dispersion relations. The results as a function of
$k_{y}/k_{x}$ are illustrated in Figure 9.

As in the case of perpendicular shocks, the flat part of the curves in panels
(a) and (b) of Figure 9 correspond to that in Figure 8 for
$k_{y}/k_{x}\lesssim 1$, arising from the fact that 2D oblique shocks resemble
1D shocks (i.e., $k_{y}=0$) in terms of the dynamics in the $x$-direction. In
the $y$-direction, however, the phase difference between $\delta\rho_{n}$ and
$\delta v_{n,y}$ changes gradually from $\approx 1.75\pi$ (equivalent to
$-\pi/4$) to $\pi$ as $k_{y}/k_{x}$ increases from $\approx 0.01$ to
$\hat{k}_{discon}$ ($\approx 9.3$ at $x=0.2$ pc). It arises because when
$k_{y}/k_{x}\ll 1$, the ionization equilibrium (i.e.
$\delta\rho_{i}/\rho_{i}=(1/2)\delta\rho_{n}/\rho_{n}$) and the ion-neutral
drag dictate the dynamics of Equation(44), yielding the phase difference
$-\pi/4$ in the $y$-direction as well. As $k_{y}/k_{x}$ increases for a given
$k_{x}$, Equation(46) implies that the wave frequency is gradually dominated
by $kc_{s}$ rather than the Doppler-shift frequency
$k_{x}V_{n,x}+k_{y}V_{n,y}$. Additionally, the drag term becomes subdominant
to the pressure term in Equation(44). Thus, Equation(44) leads to the relation
$\delta\rho_{n}/\rho_{n}\approx-\delta v_{n,y}/c_{s}$; i.e., $\delta\rho_{n}$
and $\delta v_{n,y}$ gradually become out of phase as $k_{y}/k_{x}$ increases
to $\hat{k}_{discon}$, as illustrated in panel (c) of Figure 9. Once
$k_{y}/k_{x}$ increases to $\hat{k}_{discon}$, the phase difference between
$\delta\rho_{n}$ and $\delta v_{n,y}$ jumps from $\pi$ to zero because of the
sign change of the wave frequency. Furthermore, panel (c) of Figure 9 depicts
that the phase difference between $\delta\rho_{n}$ and $\delta v_{n,y}$ is
almost zero when $k_{y}/k_{x}\gtrsim\hat{k}_{discon}$. It can be realized from
Equation(44) with $\Gamma_{n}\approx-\mathrm{i}k_{y}c_{s}$ as suggested by
Equation(46) for $k_{y}/k_{x}>\hat{k}_{discon}$, resulting in the relation
that $\delta\rho_{n}/\rho_{n}\approx\delta v_{n,y}/c_{s}$; i.e.,
$\delta\rho_{n}$ and $\delta v_{n,y}$ are almost in phase and are associated
with an acoustic wave propagating in the $y$-direction. The resulting
$v_{ph,y}$ of the acoustic wave points to the positive $y$-direction in the
frame comoving with the neutrals.

Panel (a) of Figure 9 also shows that when $k_{y}/k_{x}>\hat{k}_{discon}$, the
phase difference between $\delta\rho_{i}$ and $\delta\rho_{n}$ increases with
$k_{y}/k_{x}$, which can be realized from Equation(22). Since
$\Gamma_{i}=\Gamma_{n}+\mathrm{i}(k_{x}V_{d,x}+k_{y}V_{d,y})$, the term
$k_{x}V_{d,x}+k_{y}V_{d,y}$ is nearly zero at $k_{y}/k_{x}=k_{discon}$ and
thus is negligible compared to the recombination rate $2\beta\rho_{i}$,
resulting in the ionization–recombination equilibrium. However, as
$k_{y}/k_{x}$ increases from $\hat{k}_{discon}$, the term
$k_{x}V_{d,x}+k_{y}V_{d,y}$, which is neglected for deriving Equation(45),
increases and thus becomes comparable to and even larger than $2\beta\rho_{i}$
at a large value of $k_{y}$ in the term $\Gamma_{i}+2\beta\rho_{i}$ of
Equation(22) for the ion continuity equation. Consequently, the
ionization–recombination equilibrium is poorly attained, and the phase
difference between $\delta\rho_{i}$ and $\delta\rho_{n}$ shifts significantly
away from zero at a large value of $k_{y}/k_{x}$.

As a result of the aforementioned phase shift between $\delta\rho_{i}$ and
$\delta\rho_{n}$, the $x$-component of the linearized momentum equation for
the neutrals, i.e. Equation(19), implies that the phase difference between
$\delta\rho_{n}$ and $\delta v_{n,x}$ shifts accordingly away from
$\sim-\pi/4$ (i.e. a typical phase shift for the 1D drag instability; see GC,
20) to a large negative value as $k_{y}/k_{x}$ increases from
$\hat{k}_{discon}$. Together with the continuity equation for the neutrals
described by Equation(18), the phase shift leads to the change of the growth
rate with $k_{y}/k_{x}$ and yields no growth at a certain large value of
$k_{y}/k_{x}$ ($\sim 10^{2}$) where $\delta\rho_{n}$ and $\delta v_{n,x}$ are
almost out of the phase, as shown in panel (b) of Figure 9. At this point, the
wave is no longer unstable but is transformed into an acoustic wave, with
$v_{ph,x}$ pointing to the negative $x$-direction in the comoving frame of the
neutrals. The $v_{ph,y}$ of the acoustic wave points to the positive
$y$-direction in the comoving frame with the neutrals, as was explained
earlier in this subsection.

Figure 9: Phase differences between perturbations as a function of
$k_{y}/k_{x}$ for the unstable mode in the fiducial model of an oblique shock
with $\theta_{0}=45^{\circ}$. The cases of the four values of $k_{x}$ are
shown. The phase differences are not presented when the growth rate is zero
(refer to the left panel of Figure 8).

###  3.4 Growth of an unstable mode and shock obliquity

The distinct behaviors of wave frequency for a large value of $k_{y}/k_{x}$ –
frequency jump versus discontinuity – between Figure 2 (perpendicular shock
with $\theta_{0}=90^{\circ}$) and Figure 8 (oblique shock with
$\theta_{0}=45^{\circ}$) suggests that a shock with a proper range of the
oblique angle between $45^{\circ}$ and $90^{\circ}$ can allow for the unstable
modes with both behaviors of wave frequency. Figure 10 shows the results for
the same shock model with $\theta_{0}=80^{\circ}$. Indeed, the right panel of
the figure shows that a frequency jump from negative to positive values
happens at $k_{y}/k_{x}\approx
V_{n,x}/(c_{s}-V_{n,y})\equiv\hat{k}_{jump}\approx 27$, and a frequency
discontinuity occurs at an even smaller scale where
$k_{y}/k_{x}=|V_{d,x}|/V_{d,y}\equiv\hat{k}_{discon}\approx 37$. Consequently,
the growth rate of the modes with large $k_{x}$ (i.e., $20k_{fid}$ and
$50k_{fid}$) rises around $\hat{k}_{jump}$ and the growth rates all go to zero
at $\hat{k}_{discon}$, as shown and expected in the left panel of Figure 10.
Note that $\hat{k}_{jump}$ involves $V_{n,y}$ for an oblique shock because the
Doppler-shift frequency is $k_{x}V_{n,x}+k_{y}V_{n,y}$. What happens is that
as the oblique angle $\theta_{0}$ decreases from $90^{\circ}$, $V_{n,y}$ and
thus $V_{d,y}$ start deviating from zero such that
$\hat{k}_{jump}<\hat{k}_{discon}<\infty$. Consequently, the modes with
behaviors of the frequency jump and discontinuity both appear in the case for
$\theta_{0}=80^{\circ}$. As $\theta_{0}$ continues to decrease (i.e. the shock
is more oblique) until $\hat{k}_{discon}$ becomes smaller than
$\hat{k}_{jump}$, the mode with the behavior of a frequency jump disappears at
$k_{y}/k_{x}=\hat{k}_{jump}$, leaving the sole phenomenon of the frequency
discontinuity at $k_{y}/k_{x}=\hat{k}_{discon}$ such as the case for
$\theta_{0}=45^{\circ}$. At the location of $x\approx 0.2$ pc in the fiducial
model, $\hat{k}_{jump}=\hat{k}_{discon}$ occurs between
$\theta_{0}=78^{\circ}$ and 79∘; i.e., it occurs when the shock is mildly
oblique. Analogous to perpendicular shocks, the presence of $\hat{k}_{jump}$
enables substantial growth of the slowly propagating unstable mode (i.e.,
unlimited by the shock width) in a mildly oblique shock. Consequently, the
range of $\hat{k}_{jump}$ is similar to that for the perpendicular shocks in
the typical environments of star-forming clouds.

Figure 10: Same as Figures 2 & 8 but for $\theta=80^{\circ}$. The wave
frequency is not presented when the growth rate is zero.

Because the slowly traveling mode with substantial growth appears only in a
background environment close to a perpendicular shock, the growth of the drag
instability in a moderately or exceedingly oblique shock is still limited by
the short time span for the unstable mode to remain within a shock. We explore
this issue by comparing the C-shock model described by the fiducial model and
by the model V06 from Table 2 in GC (20). The pre-shock conditions of model
V06 are given by $n_{0}=200$ cm-3, $v_{0}=6$ km/s, $B_{0}=10\mu$G, and
$\chi_{i0}=5$. Model V06 is discussed in GC (20) because it exhibits a larger
MTG than that of the fiducial model in a 1D shock due to its broader shock
width. The MTG is easily computed for a 1D shock, whereas the same calculation
involves a more elaborate work for a 2D shock. In general, the growth rate
$\Gamma_{grow}$ of a mode with a particular wave frequency $\omega_{wave}$ is
a function of both $k_{x}$ and $k_{y}$, which vary with $x$ in a 2D shock. In
this work, we do not intend to perform an accurate analysis to identify the
combination of $k_{x}(x)$ and $k_{y}(x)$ that provides the MTG of the unstable
mode. Rather, we keep $k_{y}$ uniform but allow $k_{x}$ to change with $x$
across a shock for a given wave frequency of an unstable mode. We then
estimate the MTG of an oblique shock by varying the wave frequency for a few
constant $k_{y}$, which we refer to as MTG${}_{k_{y}}$ in short. Comparing
MTG${}_{k_{y}}$ for different $k_{y}$ and $\theta_{0}$ would provide the
guidance on how the MTG varies with the initial oblique angle $\theta_{0}$.
The results are enumerated in Table 1 for the fiducial model and Table 2 for
model V06. Given the $k_{y}$ of an unstable mode and $\theta_{0}$ for each
C-shock model, the tables display MTG${}_{k_{y}}$ with the corresponding mode
frequency $\omega_{wave}$ as well as the value of $k_{x}$ (in terms of
$k_{y}/k_{x}$) ranging from the beginning to the end of the C-shock width.

Table 1 shows that in the fiducial shock model, MTG${}_{k_{y}}$ does not vary
significantly with $\theta_{0}$. The smaller the $\theta_{0}$, the larger the
compression of the density and the magnetic field (see Figure 7), leading to a
slightly larger growth rate. However, a smaller $\theta_{0}$ results in a
slightly narrower shock width (see Figure 7 and the parameter $L_{shock}$ in
Table 1). As a result of the competition between these two moderate effects,
the oblique angle does not noticeably affect the overall growth of the drag
instability in the fiducial model. Moreover, when $k_{y}/k_{fid}\lesssim 1$,
MTG${}_{k_{y}}$ is comparable to the MTG for the 1D shock ($=9.9$; refer to
GC, 20) and is smaller by more than a factor of 2 for $k_{y}/k_{kid}=10$. It
can be expected from the range of the resulting $k_{y}/k_{x}$ shown in Table 1
. When $k_{y}/k_{x}\ll 1$, the instability behaves as a 1D mode for which
$k_{y}=0$. On the other hand, when $k_{y}/k_{x}\gtrsim 1$, the overall growth
rate decreases within the shock, and thus MTG${}_{k_{y}}$ is small.

Figure 11: The amplitude of perturbations normalized by the $y$ component of
the magnetic perturbation $|\delta B_{y}|/B$ for the unstable mode with
$k_{y}/k_{fid}=0.1$ and $\omega_{wave}=-2\times 10^{-11}$ s-1 in the fiducial
model (left panel) and model V06 (right panel) of the C-shock with
$\theta_{0}=45^{\circ}$.

Similar to the fiducial model, Table 2 shows that MTG${}_{k_{y}}$ in model V06
does not vary significantly with $\theta_{0}$ as well. Though in the case of
$k_{y}/k_{fid}=0.1$, MTG${}_{k_{y}}$ can increase moderately from 36 for
$\theta_{0}=70^{\circ}$ to 47 for $\theta_{0}=20^{\circ}$. This trend actually
also happens in the fiducial model with $k_{y}/k_{fid}=0.1$ but is less
prominent (see Table 1). For a smaller $\theta_{0}$, the effect of the
slightly larger growth rate predominates more than the effect of the slightly
narrower width, resulting in a larger MTG${}_{k_{y}}$. Evidently, a wider
shock in model V06 (see the parameter $L_{shock}$ in Tables 1 and 2) allows an
unstable mode to grow more and therefore makes this trend clearer. In
contrast, the two competing effects are comparable for $k_{y}/k_{fid}=1$ and
10; hence, the MTG${}_{k_{y}}$ changes less noticeably with $\theta_{0}$.
Similar to the fiducial model, the overall trend of the decrease in
MTG${}_{k_{y}}$ with increasing $k_{y}$ also happens in model V06, for the
same reason. Tables 1 and 2 also show that MTG${}_{k_{y}}$ in model V06 is in
general larger than that in the fiducial model, which is expected because the
shock in model V06 has a broader shock width for the instability to acquire
more time to grow. In summary, for a C-shock with the oblique angle $\lesssim
70^{\circ}$, we expect that the MTG of the drag instability in the fiducial
model is about 10 and that in model V06 is about 30-47. Therefore, we expect
that the unstable mode responsible for the MTG has $k_{y}/k_{x}\lesssim 1$.

Figure 11 shows the magnitude of the perturbations across the shock width for
the unstable mode corresponding to $k_{y}/k_{fid}=0.1$ and
$\theta_{0}=45^{\circ}$ in Tables 1 (left panel) and 2 (right panel). It is
evident from the figure that the magnitude of the density perturbation is
always larger than that of the velocity and magnetic field perturbations of
the unstable mode. In fact, it is true for all cases listed in Tables 1 and 2.
The density enhancement is expected to play a critical role in the dynamics of
the drag instability in an oblique shock as well.

When $\theta_{0}\lesssim 20^{\circ}$ in our fiducial model, $B_{y}$ and
$\rho_{n}$ are compressed extremely quickly near the beginning and the end of
the C-shock width, respectively. Therefore, the WKBJ approximation with
$k_{x}=k_{fid}$ becomes invalid in these shock regions. Except for these
regions, the basic behavior of the drag instability for $\theta_{0}\lesssim
20^{\circ}$ is similar to that for $\theta_{0}=45^{\circ}$, shown in Figure 8
at $x\approx 0.2$ pc. When $\theta_{0}$ is smaller than the critical angle
$\theta_{crit}\approx 6^{\circ}$, the background state admits two additional
solutions with field reversal, referred to as intermediate shocks (e.g.,
Wardle, 1998; Chen & Ostriker, 2014). In this study, we restrict ourselves to
the oblique shocks without field reversal and leave the instability analysis
for intermediate shocks to a future work.

Table 1: MTG${}_{k_{y}}$ for various $k_{y}$ and $\theta_{0}$ in the fiducial
model

$k_{y}/k_{fid}$ | $\theta_{0}$ |  $L_{shock}$ (pc) |  $\omega_{wave}$ (1/s) | $k_{y}/k_{x}$ | MTG${}_{k_{y}}$
---|---|---|---|---|---
0.1 | 70∘ | 0.48 |  $-3$e$-11$ | 0.035–0.005 | 10.0
0.1 | 45∘ | 0.40 |  $-2$e$-11$ | 0.052–0.007 | 10.3
0.1 | 20∘ | 0.24 |  $-2$e$-11$ | 0.052–0.005 | 10.9
1 | 70∘ | 0.48 |  $-3$e$-11$ | 0.35–0.05 | 9.57
1 | 45∘ | 0.40 |  $-3$e$-11$ | 0.35–0.05 | 9.49
1 | 20∘ | 0.24 |  $-3$e$-11$ | 0.35–0.05 | 9.55
10 | 70∘ | 0.48 |  $-5$e$-11$ | 1.97–0.32 | 4.4
10 | 45∘ | 0.40 |  $-5$e$-11$ | 2.4–0.32 | 3.8
10 | 20∘ | 0.24 |  $-5$e$-11$ | 2.4–0.08 | 3.4

Note. — Given a uniform value of $k_{y}$ (in units of $k_{fid}$), the mode
frequency $\omega_{wave}$ and the corresponding range of $k_{x}(x)$ (presented
in terms of $k_{y}/k_{x}$) associated with the MTG at a constant $k_{y}$
(i.e., MTG${}_{k_{y}}$) of the drag instability are listed for the fiducial
model of a steady C-shock with various oblique angles $\theta_{0}$. The shock
width $L_{shock}$ is estimated using Equation(A19) in Chen & Ostriker (2012).

Table 2: Same as Table 2 but for the C-shock model given by model V06.  $k_{y}/k_{fid}$ | $\theta_{0}$ |  $L_{shock}$ (pc) |  $\omega_{wave}$ (1/s) | $k_{y}/k_{x}$ | MTG${}_{k_{y}}$
---|---|---|---|---|---
0.1 | 70∘ | 2.08 |  $-2$e$-11$ | 0.062–0.011 | 36
0.1 | 45∘ | 1.70 |  $-2$e$-11$ | 0.062–0.009 | 40
0.1 | 20∘ | 1.0 |  $-2$e$-11$ | 0.063–0.007 | 47
1 | 70∘ | 2.08 |  $-3$e$-11$ | 0.42–0.07 | 31
1 | 45∘ | 1.70 |  $-3$e$-11$ | 0.42–0.07 | 30
1 | 20∘ | 1.0 |  $-3$e$-11$ | 0.45–0.05 | 31
10 | 70∘ | 2.08 |  $-7$e$-11$ | 1.72–0.33 | 7.3
10 | 45∘ | 1.70 |  $-8$e$-11$ | 1.75–0.28 | 5.9
10 | 20∘ | 1.0 |  $-8$e$-11$ | 1.73–0.24 | 5.7

##  4 Summary and Discussions

In this work, we extend the study of the drag instability in 1D perpendicular
C-shocks by GC (20) to 2D perpendicular and oblique C-shocks. We focus on the
fiducial model for an isothermal steady C-shock with the pre-shock conditions
described by the same fiducial model as GC (20) in the typical environment of
star-forming clouds. The WKBJ linear analyses are subsequently performed based
on the background states in the fiducial model. To understand the underlying
physics for the linear results, we make an attempt to derive simplified
dispersion relations, aided by the auxiliary analysis of phase differences
between perturbation quantities. We observe that the drag instability remains
in a 2D shock, and its behavior in general depends on $k_{y}/k_{x}$. When
$k_{y}/k_{x}\lesssim 1$ (i.e., transversely large-scale modes), the growth
rate $\Gamma_{grow}$ and wave frequency $\omega_{wave}$ of the drag
instability in a 2D shock are similar to those in a 1D shock, which is
insensitive to the initial oblique angle $\theta_{0}$ of the shock. When
$k_{y}/k_{x}\gtrsim 1$ (i.e., transversely small-scale modes), the drag
instability is characterized by an unstable mode coupled with the acoustic
mode primarily along the $y$-direction (note that the acoustic mode is the
slow mode in the case of a perpendicular shock). Additionally, in contrast to
the perpendicular shock, $V_{d,y}$ exists in an oblique shock. Therefore,
there exists a particular mode for an oblique shock with
$k_{y}/k_{x}=\hat{k}_{discon}$ where the growth rate is zero and discontinuity
in the wave frequency appears (see Equation(46)).

When the shock is less oblique (i.e., $\theta_{0}\gtrsim 80^{\circ}$ in the
fiducial model), there exists a jump transition of wave frequency in the shock
frame – from the negative Doppler-shift frequency due to the shock flow to the
positive acoustic wave frequency in the $y$-direction – for a mode with
$k_{y}/k_{x}\sim\hat{k}_{jump}$. Owing to the small wave frequency in the
shock frame, this unstable mode propagates slowly within a shock and thus has
sufficient time to potentially grow to a nonlinear phase, thereby contributing
to the maximum growth. While the density enhancement of the ions approximately
lies in phase with that of the neutrals by means of the ionization equilibrium
for $k_{y}/k_{x}<\hat{k}_{jump}$, the same phase overlap for the density
perturbations is maintained by the fast acoustic wave for
$k_{y}/k_{x}>\hat{k}_{jump}$. For the mode with an exceedingly large $k_{x}$,
there is no growing mode for $k_{y}/k_{x}\lesssim 1$, as the drag instability
is suppressed by the pressure effect (Gu et al., 2004; GC, 20). However,
unstable modes appear for $k_{y}/k_{x}\sim\hat{k}_{jump}$, as $k_{x}$ is large
enough for the Doppler-shift frequency of the slowly traveling wave to
dominate over the ionization rate, which in turn produces a proper phase
difference between $\delta v_{n,x}$ and $\delta\rho_{n}$ for the drag
instability to occur.

On the other hand, when the shock is more oblique (i.e., $\theta_{0}\lesssim
80^{\circ}$ in the fiducial model), this slowly propagating unstable mode
disappears. The maximum growth of the drag instability is limited by the short
time span of an unstable mode to stay within a shock and hence is given by
MTG, as is the case for a 1D perpendicular shock. We compute $MTG_{k_{y}}$,
i.e. the MTG for a constant $k_{y}$, to infer the MTG of the drag instability
for a given $\theta_{0}$. We find that the MTG of the drag instability is
contributed by the mode with $k_{y}/k_{x}\ll 1$ (i.e., almost 1D mode) and is
expected to be about 10 insensitive to the initial oblique angle $\theta_{0}$
in the fiducial model. We also conduct the linear analysis for the C-shock
model V06, which has a larger shock width than that of the fiducial model and
thus exhibits a larger MTG in the 1D shock (GC, 20). We find that the MTG of
the drag instability arises from the mode with $k_{y}/k_{x}\ll 1$ as well and
increases from about 36 to 47 as the initial oblique angle $\theta_{0}$
decreases from 70∘ to 20∘. The overall larger MTG in model V06 is expected for
a shock with a larger width. A larger MTG for a smaller $\theta_{0}$ is
primarily caused by the stronger shock compression. In all the cases that we
consider (see Tables 1 and 2, as well as the case of the perpendicular shock),
the magnitude of the density perturbations is much larger than that of the
velocity and magnetic field perturbations (e.g., Figures 6 & 11), implying
that the density enhancement predominantly governs the dynamics of the
instability in the linear regime.

Self-gravity is not considered in our linear analysis in order to study and
present the basic properties of the drag instability in a clean physical
picture. The background magnetic fields with different initial oblique angles
tend to be approximately parallel to the shock front during shock compression
(see Figure 7). In this work, the minimal value of $k_{y}/k_{fid}$ that we
show for the results is 0.01. This transversely large-scale mode has a
transverse wavelength larger than the Jeans scale ($\sim
c_{s}/\sqrt{G\rho_{n}}$) and therefore may be subject to the gravitational
instability along the compressed field lines primarily in the $y$-direction.
However, it should be kept in mind that the background state of the shock is
not static; thus, applying the Jeans criterion is debatable. Our analysis
based on $MTG_{k_{y}}$ suggests that the transversely large-scale mode would
give the MTG of the drag instability within a shock (see Tables 1 & 2),555The
transversely large-scale mode with $k_{y}/k_{fid}=0.1$ shown in the tables is
marginally gravitationally stable against thermal pressure along the field
direction almost in the $y$-direction due to shock compression. Having said
that, we should bear in mind that the shock is not a static structure for the
Jeans criterion to reasonably apply. except when the shock is mildly oblique.
Considered together, after the most unstable mode of the drag instability sets
in, it could grow to become gravitationally unstable more quickly in the
direction of the field line against thermal pressure as a result of its larger
transverse scale. In this regard for a laminar background, it could be
interesting to investigate the interaction between the drag and gravitational
instabilities by including self-gravity in the linear analysis.

While we assume that C-shocks arise from supersonic, turbulent flows, the
effects of turbulence are not actually modeled in our linear analysis using
any phenomenological approaches, such as turbulent diffusion or turbulent
pressure. The turbulent diffusion, if it exists, may weaken/eliminate
transversely small-scale modes of the drag instability or may promote the drag
instability by enhancing ambipolar diffusion in a shock (e.g., Li & Nakamura,
2004). The supersonic turbulent pressure, if it exists, may support the
shocked gas against gravitational collapse along the compressed field lines on
the scales smaller than the turbulent Jeans scale. On the other hand,
supersonic turbulence in shocks is known to dissipate quickly over one
turbulent crossing time, provided that there is no energy supply from the
turbulent injection scale by feedback processes from protostars (e.g., Mac Low
et al., 1998; Stone et al., 1998; Mac Low, 1999; Nakamura & Li, 2005).
Nevertheless, even large-scale modes within C-shocks would not be affected by
the decay of the large-scale turbulence, as that would only reduce the
strength of C-shocks present in the environment, rather than modifying the
behavior of an individual C-shock during its passage. In a nearly
perpendicular shock, the slowly propagating modes responsible for the maximum
growth are transversely small-scale modes (i.e.,
$k_{y}/k_{x}\sim\hat{k}_{jump}\sim 10$ in our models). Analogous to dense
regions seeded by small-scale turbulence in a shock-compressed layer (Chen &
Ostriker, 2014), the small-scale turbulence may help initiate the perturbation
of these small-scale modes of the drag instability in compressed gas within
shocks. As the modes grow to the nonlinear regime, the nonlinear saturation of
the drag instability by itself could be the prime candidate to drive the
turbulence within a C-shock. Compression of background turbulence is likely to
be negligible in comparison, while the background turbulence driven from
stellar feedback operates on much larger time and length scales from those of
an individual C-shock.

A possible consequence of the turbulence driven by drag-induced instabilities
is clump formation. One of the notable examples is the clump formation of dust
in the turbulence driven by the streaming instability due to the dust–gas drag
in a protoplanetary disk under certain favorable conditions (Johansen &
Youdin, 2007). Moreover, we show that the density perturbation dominates over
the velocity and magnetic field perturbations in the linear regime of the drag
instability, which may hint the dynamical importance of density enhancement.
The observations of nearby molecular clouds with the Herschel Space
Observatory suggested that about 70%-80% of dense cores lie within filaments
(Polychroni et al., 2013; Könyves et al., 2015). The numerical simulation
conducted by Chen & Ostriker (2012) revealed that a pair of steady C-shocks
propagate away from the shocked layer of two colliding flows where prestellar
cores and filaments can form (Chen & Ostriker, 2014). The drag instability in
a steady C-shock might be potentially capable of core formation lying outside
of filaments. Alternatively, the drag instability is likely to nonlinearly
saturate, perhaps at a level as high as the velocity difference across the
shock front. This would substantially limit the ability of the instability to
drive a subsequent gravitational instability beyond what has already been
produced by the jump in density across the shock. In any case, whether the
nonlinear outcome of the density enhancement by the drag instability in a
C-shock can facilitate the clump/core formation is beyond the reach of the
linear analysis. Given the small longitudinal wavelength of the drag
instability within a C-shock, nonideal MHD simulations with high resolutions
are required to explore this possibility (see the discussion section in GC,
20).

Apart from the steady C-shock, a transient C-shock appearing in the shocked
layer of two colliding flows has been modeled to be a promising site for the
major formation of cores and filamentary structures (Nakamura & Li, 2008; Chen
& Ostriker, 2012, 2014). We restrict ourselves to the drag instability in
steady C-shocks because the background state has settled to an equilibrium
state and thus provides an easy test bed for demonstrating the existence of
the drag instability in a linear analysis. Nevertheless, the ion-neutral drift
is expected to be extremely fast in a transient C-shock, perhaps favoring the
occurrence of the drag instability (Gu et al., 2004). Given the time-dependent
nature of a transient C-shock, setting up an appropriate background state for
a perturbation theory is expected to be challenging.

Undoubtedly, the aforementioned issues and possible implications that we have
discussed related to the drag instability are highly speculative and
intriguing, requiring prudent and elaborate studies for further investigation.
After all, the interplay between ambipolar diffusion, turbulence, shocks,
gravity, etc. for the dynamical processes of star formation has been a complex
and broad topic. From a theoretical perspective, the 2D linear analysis
conducted in this work based on the 1D analysis presented in Gu et al. (2004)
and GC (20) would advance our understanding of the instabilities of
astrophysical plasma in general. In practice, our framework provides the basic
properties of the drag instability to be studied in future nonideal MHD
simulations for the confirmation of their existence and for understanding the
nonlinear outcome of the linear instability in C-shocks.

We are grateful to Che-Yu Chen, Min-Kai Lin, Hau-Yu Baobab Liu, and Chien-
Chang Yen for useful discussions. We would also like to thank the referee for
helpful comments that greatly improved the manuscript, especially the contents
of the discussion section. This work has been supported by the Ministry of
Science and Technology in Taiwan through the grant MOST 109-2112-M001-052.

## References

  * André et al. (2014) André, P., Di Francesco, J., Ward-Thompson, D., Inutsuka, S. -I., Pudritz, R. E., Pineda, J. E. 2014, Protostars and Planets VI, Henrik Beuther, Ralf S. Klessen, Cornelis P. Dullemond, and Thomas Henning (eds.), University of Arizona Press, Tucson, 914 pp., p.27-51
  * Ballesteros-Paredes et al. (2007) Ballesteros-Paredes, J., Klessen, R. S., Mac Low, M.-M., et al. 2007, Protostars and Planets V, B. Reipurth, D. Jewitt, and K. Keil (eds.), University of Arizona Press, Tucson, 951 pp., p.63-80
  * Blitz et al. (2007) Blitz, L., Fukui, Y., Kawamura, A., Leroy, A., Mizuno, N., & Rosolowsky, E., 2007, Protostars and Planets V. Univ. Arizona Press, Tucson, AZ, p. 81
  * Chen & Ostriker (2012) Chen, C.-Y, & Ostriker, E. 2012, ApJ, 744, 124
  * Chen & Ostriker (2014) Chen, C.-Y, & Ostriker, E. 2014, ApJ, 785, 69 2, ApJ, 567, 947
  * Dalgarno (2006) Dalgarno, A. 2006, Proceedings of the National Academy of Science, 103, 12269
  * Draine (1980) Draine, B. T. 1980, ApJ, 241, 1021
  * Draine & McKee (1993) Draine, B. T., & McKee, C. F. 1993, ARA&A, 31, 373
  * Draine et al (1983) Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, ApJ, 264, 485
  * Elmegreen & Scalo (2004) Elmegreen, B. G., & Scalo, J. 2004, Annu. Rev. Astron. Astrophys., 42, 211
  * Engargiola et al. (2003) Engargiola, G., Plambeck, R. L., Rosolowsky, E., & Blitz, L., 2003, ApJS, 149, 343
  * Falle et al. (2009) Falle, S.A.E.G., Hartquist, T.W., van Loo, S. In: Pogorlov, N.V., Audit, E., Colella, P., Zank, G.P. 2009, (eds.) Numerical Modeling of Space Plasma Flows, p. 80. Astronomical Society of the Pacific
  * Flower & Pineau Des Forêts (1998) Flower, D. R., & Pineau Des Forêts, G. 1998, MNRAS, 297, 1182
  * Flower & Pineau Des Forêts (2010) Flower, D. R., & Pineau Des Forêts, G. 2010, MNRAS, 406, 1745
  * Fukui & Kawamura (2010) Fukui, Y., & Kawamura, A. 2010, ARA&A, 48, 547
  * Girichidis et al. (2020) Girichidis, P., Offner, S.S.R., Kritsuk, A.G. et al. 2020, Space Sci Rev, 216, 68
  * GC (20) Gu, P.-G., & Chen, C.-Y. 2020, ApJ, 898, 67 (GC20)
  * Gu et al. (2004) Gu, P.-G., Lin, D. N. C., & Vishniac, E. T., Astrophysics & Space Science, 292, 261
  * Gusdorf et al. (2008) Gusdorf, A., Cabrit, S., Flower, D. R., et al. 2008, A&A, 482, 809
  * Hennebelle & Falgarone (2012) Hennebelle, P., & Falgarone, E. 2012, Astron. Astrophys. Rev., 20, 55
  * Hennebelle & Inutsuka (2019) Hennebelle, P., & Inutsuka, S.-I. 2019, Frontiers in Astronomy and Space Sciences, 6, 5
  * Hezareh et al. (2010) Hezareh, T., Houde, M., McCoey, C., et al. 2010, ApJ, 720, 603
  * Hezareh et al. (2014) Hezareh, T., Csengeri, T., Houde, M., et al. 2014, MNRAS, 438, 663
  * Indriolo & McCall (2012) Indriolo, N., & McCall, B. J. 2012, ApJ, 745, 91
  * Jeffreson & Kruijssen (2018) Jeffreson, S. M. R. & Kruijssen, J. M. D. 2018, MNRAS, 476, 3688
  * Johansen & Youdin (2007) Johansen, A., & Youdin, A. 2007 ApJ, 662, 627
  * Kawamura et al. (2009) Kawamura, A. et al., 2009, ApJS, 184, 1
  * Kennicutt & Evans (2012) Kennicutt, R. C., & Evans, N. J. 2012, Annu. Rev. Astron. Astrophys. 50, 531
  * Könyves et al. (2015) Könyves, V., André, P., Menśhchikov, A., Palmeirim, P., Arzoumanian, D., Schneider, N., et al. 2015, A&A, 584, 91
  * Lehmann & Wardle (2016) Lehmann, A., & Wardle, M. 2016, MNRAS, 455, 2066
  * Li & Nakamura (2004) Li, Z.-Y., & Nakamura, F. 2004, ApJ, 609. 83
  * Li et al. (2014) Li, H.-B., Goodman, A., Sridharan, T. K., et al. 2014, Protostars and Planets VI, 101
  * Li & Houde (2008) Li, H.-B., & Houde, M. 2008, ApJ, 677, 1151
  * Mac Low (1999) Mac Low, M.-M. 1999, ApJ, 524, 169
  * Mac Low et al. (1998) Mac Low, M.-M., Klessen, R. S., Burkert A. et al., 1998, Phys. Rev. Lett. 80, 2754
  * McKee et al. (2010) McKee, C. F., Li, P. S., & Klein, R. I. 2010, ApJ, 720, 1612
  * Meidt et al. (2015) Meidt, S. E. et al., 2015, ApJ, 806, 72
  * Mestel & Spitzer (1956) Mestel, L., & Spitzer, L. 1956, MNRAS, 116, 503
  * Miura et al. (2012) Miura, R. E. et al., 2012, ApJ, 761, 37
  * Murray (2011) Murray, N., 2011, ApJ, 729, 133
  * Nakamura & Li (2005) Nakamura, F., & Li, Z.-Y. 2005, ApJ, 631, 411
  * Nakamura & Li (2008) Nakamura, F., & Li, Z.-Y. 2008, ApJ, 687, 354
  * Polychroni et al. (2013) Polychroni, D., Schisano, E., Elia, D., Roy, A., Molinari, S., Martin, P., et al. 2013, ApJ,777, 33
  * Shu (1992) Shu, F. H. 1992, in Physics of Astrophysics, Vol. II, ed. F. H. Shu (Mill Valley, CA: Univ. Science Books)
  * Smith & Mac Low (1997) Smith, M. D., & Mac Low, M.-M. 1997, A&A, 326, 801
  * Spitzer (1956) Spitzer, L. 1956, Physics of Fully Ionized Gases, New York: Interscience Publishers
  * Stone (1997) Stone, J. M. 1997, ApJ, 487, 271
  * Stone et al. (1998) Stone, J. M., Ostriker, E. C., Gammie, C. F. 1998, ApJ508, 99
  * Tang et al. (2018) Tang, K. S., Li, H.-B., & Lee, W.-K. 2018, ApJ, 862, 42
  * Tielens (2005) Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium, Cambridge, UK: Cambridge University Press
  * Valdivia et al. (2017) Valdivia, V., Godard, B., Hennebelle, P., et al. 2017, A&A, 600, A114
  * Wardle (1990) Wardle, M. 1990, ApJ, 246, 98
  * Wardle (1991) Wardle, M. 1991, MNRAS, 251, 119
  * Wardle (1998) Wardle, M. 1998, MNRAS, 298, 507
  * Xu & Li (2016) Xu, D., & Li, D. 2016, ApJ, 833, 90
  * Zweibel (2015) Zweibel, E. 2015, Ambipolar Diffusion, in Magnetic Fields in Diffuse Media, ed. Alexander Lazarian, Elisabete M. de Gouveia Dal PinoClaudio Melioli, 407, 285