The Drag Instability in a 1D Isothermal C-shock

Following the simplified scenario discussed in Chen & Ostriker (2012) that the magnetic field is perpendicular to the gas flow toward the +x direction through the shock, the equilibrium equations in this 1D C-shock system are given by

$$\begin{eqnarray}&&{V}_{n}\displaystyle \frac{d{\rho }_{n}}{{dx}}=-{\rho }_{n}\displaystyle \frac{{{dV}}_{n}}{{dx}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{V}_{n}\displaystyle \frac{{{dV}}_{n}}{{dx}}=\gamma {\rho }_{i}{V}_{d}-\displaystyle \frac{{c}_{s}^{2}}{{\rho }_{n}}\displaystyle \frac{d{\rho }_{n}}{{dx}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\beta {\rho }_{i}^{2}={\xi }_{\mathrm{CR}}{\rho }_{n},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\gamma {\rho }_{n}{V}_{d}=-{V}_{A,i}^{2}\displaystyle \frac{d\mathrm{ln}B}{{dx}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{V}_{i}\displaystyle \frac{{dB}}{{dx}}=-B\displaystyle \frac{{{dV}}_{i}}{{dx}}.\end{eqnarray} \tag{ 10 }$$

In the above equilibrium equations, we consider the strong-coupling approximation under which the ion–neutral drag is balanced by the magnetic pressure gradient for the ions; i.e., $\gamma {\rho }_{n}{L}_{B}{V}_{d}={V}_{A,i}^{2}$ where ${L}_{B}\equiv {(-d\mathrm{ln}B/{dx})}^{-1}$. In addition, the equilibrium between cosmic ionization and recombination is assumed; namely, $\beta {\rho }_{i}^{2}={\xi }_{\mathrm{CR}}{\rho }_{n}$ (see CO12 for justifications of this choice). We note that these equilibrium states are consistent with the background states considered by GLV.

By subjecting the equilibrium equations to the zero-gradient boundary conditions (d/dx = 0) far upstream and downstream (i.e., no structures in the steady pre- and postshock regions), CO12 derived the 1D structure equation of a C-shock as follows (here and throughout this paper, we use the subscript 0 to denote a physical quantity in the preshock region):

$$\begin{eqnarray}&&\displaystyle \frac{{{dr}}_{B}}{{dx}}=-\displaystyle \frac{\gamma {\rho }_{i,0}}{{v}_{0}}{{ \mathcal M }}_{A}^{2}\displaystyle \frac{{r}_{n}^{3/2}}{{r}_{B}}\left(\displaystyle \frac{1}{{r}_{B}}-\displaystyle \frac{1}{{r}_{n}}\right),\end{eqnarray} \tag{ 11 }$$

where the field compression ratio is ${r}_{B}\equiv B/{B}_{0}={V}_{i,0}/{V}_{i}$, the neutral compression ratio is ${r}_{n}\equiv {\rho }_{n}/{\rho }_{n,0}={V}_{n,0}/{V}_{n}$, and the Alfvén Mach number ${{ \mathcal M }}_{A}$ for the shock velocity v₀ is defined as v₀/V_A,n,0. Note that r_n can be written as a function of r_B:

$$\begin{eqnarray}&&\begin{array}{l}{r}_{n}=\left[1+2{{ \mathcal M }}_{A}^{2}+{\beta }_{\mathrm{plasma}}-{r}_{B}^{2}\right.\\ -\left.\sqrt{{\left[1+2{{ \mathcal M }}_{A}^{2}+{\beta }_{\mathrm{plasma}}-{r}_{B}^{2}\right]}^{2}-8{\beta }_{\mathrm{plasma}}{{ \mathcal M }}_{A}^{2}}\right]/2{\beta }_{\mathrm{plasma}},\end{array}\end{eqnarray} \tag{ 12 }$$

with the plasma beta value ${\beta }_{\mathrm{plasma}}\equiv 8\pi {\rho }_{0}{c}_{s}^{2}/{B}_{0}^{2}$. By placing the shock front at x = 0, the physical quantities along the C-shock can be obtained by integrating the ordinary different equation (Equation (11)) backward from far downstream, where

$$\begin{eqnarray}&&\begin{array}{l}{r}_{B}={r}_{n}=4{{ \mathcal M }}_{A}^{2}/\left[1+{\beta }_{\mathrm{plasma}}\right.\\ \left.+{\left[{\left(1+{\beta }_{\mathrm{plamsa}}\right)}^{2}+8{{ \mathcal M }}_{A}^{2}\right]}^{1/2}\right]\end{array}\end{eqnarray} \tag{ 13 }$$

to far upstream. 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).

We now consider the perturbations $U(\omega ,k)\equiv (\delta {\rho }_{i},\delta {v}_{i},{\delta B,\delta {\rho }_{n},\delta {v}_{n})}^{T}$ multiplied by $\exp [{\rm{i}}({kx}+\omega t)]$ under the Wentzel–Kramers–Brillouin–Jeffreys (WKBJ) approximation. By substituting these perturbations and the background states into the Equations (1)–(5), the following linearized equations are obtained (GLV):

$$\begin{eqnarray}&&{CU}={\rm{i}}\omega U,\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}\begin{array}{l}C=\left[\begin{array}{ccccc}-{\rm{i}}{{kV}}_{i}-2\beta {\rho }_{i} & -{\rm{i}}k{\rho }_{i} & 0 & {\xi }_{\mathrm{CR}} & 0\\ -{\rm{i}}k\tfrac{{c}_{s}^{2}}{{\rho }_{i}} & -{\rm{i}}{{kV}}_{i}-\gamma {\rho }_{n} & -{\rm{i}}k\tfrac{{V}_{A,i}^{2}}{B} & -\gamma {V}_{d} & \gamma {\rho }_{n}\\ 0 & -{\rm{i}}{kB} & -{\rm{i}}{{kV}}_{i} & 0 & 0\\ 0 & 0 & 0 & -{\rm{i}}{{kV}}_{n} & -{\rm{i}}k{\rho }_{n}\\ \gamma {V}_{d} & \gamma {\rho }_{i} & 0 & -{\rm{i}}k\displaystyle \frac{{c}_{s}^{2}}{{\rho }_{n}} & -{\rm{i}}{{kV}}_{n}-\gamma {\rho }_{i}\end{array}\right].\end{array}\end{eqnarray} \tag{ 15 }$$

Hence, we can solve the above equations as an eigenvalue problem with ω being the eigenvalue and U being the eigenfunctions. The goal is to identify a maximum-growth mode associated with the drag instability within the C-shock.

We adopt the 1D steady-state C-shock profile shown in Figure 3 of CO12 as the background state (in the shock frame) of our fiducial model. The preshock parameters are n₀ = 500 cm⁻³ (neutral number density), V_i,0 = V_n,0 = v₀ = 5 km s⁻¹ (shock velocity), B₀ = 10 μG, and ionization fraction coefficient χ_i0 = 10. Here, the parameter χ_i0 is defined in the expression ${n}_{i}={10}^{-6}{\chi }_{i0}{n}_{n}^{1/2}$ of CO12 assuming ionization-recombination equilibrium, and is therefore given by ${\chi }_{i0}\equiv {10}^{6}\sqrt{{\xi }_{\mathrm{CR}}({m}_{n}/{m}_{i})/(\beta {m}_{i})}$. We thus adopted β ≈ 10⁻⁷ cm³ s⁻¹/m_i and ξ_CR ≈ 10⁻¹⁷ s⁻¹ (m_i/m_n) in this study (see, e.g., Shu 1992; Tielens 2005), where m_n = 2.3× and m_i = 30× the hydrogen mass are considered. Indeed, χ_i0 = 10 falls in the typical range of χ_i0 observed in star-forming regions (∼1−20; see, e.g., McKee et al. 2010). The fiducial model is selected such that it can be reproduced easily by comparison with the results of CO12.

Without loss of generality, we adopt a constant wavenumber (k) of 1/0.015 pc⁻¹ to keep the fiducial model as simple as possible. Figure 1 displays the r_n/r_B ratio of the background state (left panel) and the validity of the WKBJ approximation (right panel). The C-shock transition begins from x = 0 pc and ends at approximately x = 0.4 pc. We obtained the same U-shape for the r_n/r_B ratio as that obtained by CO12 in their Figure 3. The U-shaped profile is a notable feature of the C-shock model, as demonstrated by CO12. The ratio r_n/r_B = 1 in the pre- (x < 0 pc) and postshock (x > 0.4 pc) regions where no background gradients are present. Throughout the C-shock, 1/(kL_B) and 1/(kL_p) are considerably smaller than 1; thus, the WKBJ approximation is justified. The marginally high value of 1/(kL_B) and 1/(kL_p) around x = 0 pc and x = 0.4 pc, respectively, is caused by the initial compression of the ions and thus the magnetic fields, followed by a delayed compression of the neutrals by means of the ion–neutral drag as the gas flows downstream across the steady C-shock.

Zoom In Zoom Out Reset image size

Before solving the eigenvalue problem in Equation (14), we analyze Equation (14) in terms of a couple of simplified dispersion relations, which will provide the basic and clear physics to better understand the exact solutions obtained using the eigenvalue approach. Table 1 summarizes the meanings of the main symbols that we use in the linear analysis.

Table 1.  Summary of the Main Symbols Adopted in the Linear Analysis

Symbol and definition	Meaning
${\rm{\Gamma }}\equiv {\rm{i}}\omega +{\rm{i}}{{kV}}_{n}$	rate of a mode in the comoving frame of the neutrals
ΓGLV	rate of the unstable/decaying mode derived by GLV
${{\rm{\Gamma }}}_{\mathrm{re}}\equiv 2\beta {\rho }_{i}$	recombination rate
${{\rm{\Gamma }}}_{\mathrm{grav}}\equiv \sqrt{G{\rho }_{n}}$	rate of the gravitational instability
${{\rm{\Gamma }}}_{\mathrm{th}}\equiv {{kc}}_{s}$	sound-crossing rate over one wavelength
${{\rm{\Gamma }}}_{{{kV}}_{d}}\equiv k| {V}_{d}|$	ion–neutral drift rate across a distance of wavelength
${{\rm{\Gamma }}}_{\mathrm{alf},i}\equiv {{kV}}_{A,i}$	speed of the Alfvén wave in the ions crossing one wavelength
${{\rm{\Gamma }}}_{\mathrm{alf},n}\equiv {{kV}}_{A,n}$	speed of the Alfvén wave in the neutrals crossing one wavelength
${{\rm{\Gamma }}}_{\mathrm{ambi}}\equiv {k}^{2}{D}_{\mathrm{ambi}}$	ambipolar diffusion rate
${\gamma }_{n}\equiv \gamma {\rho }_{n}$	ion collision rate with the neutrals
${\gamma }_{i}\equiv \gamma {\rho }_{i}$	neutral collision rate with the ions
${{\rm{\Gamma }}}_{\mathrm{grow}}$	growth rate of an unstable mode
${\omega }_{\mathrm{wave},n}$	wave frequency of an unstable mode in the comoving frame of the neutrals
${\omega }_{\mathrm{wave}}\equiv \mathrm{Re}[\omega ]$	wave frequency of a mode
${D}_{\mathrm{ambi}}\equiv {V}_{A,n}^{2}/\gamma {\rho }_{i}$	ambipolar diffusion coefficient
$\gamma \approx 3.5\times {10}^{13}$ cm3 s−1 g−1	drag force coefficient (Draine et al. 1983)

Download table as:  ASCIITypeset image

The analysis begins with a brief review of the drag instability. GLV indicated that when the drift velocity V_d is sufficiently high, a type of overstability, called the drag instability, can occur. Specifically, for the definite occurrence of the instability, the rate of the mode ${\rm{\Gamma }}\equiv {\rm{i}}\omega +{\rm{i}}{{kV}}_{n}$ observed in the comoving frame of the neutrals is considerably lower than both the recombination rate Γ_re ≡ 2βρ_i and the ion–neutral drift rate across a distance of wavelength (i.e., ${{\rm{\Gamma }}}_{{{kV}}_{d}}\equiv k| {V}_{d}|$), whereas it is considerably higher than the neutral collision rate with the ions (i.e., γ_i ≡ γρ_i), the sound-crossing rate over one wavelength (i.e., Γ_th ≡ kc_s), and the rate of the gravitational instability (i.e., ${{\rm{\Gamma }}}_{\mathrm{grav}}\equiv \sqrt{G{\rho }_{n}}$). Although self-gravity is not included in our equations, Γ_grav is still estimated to evaluate its significance. When the aforementioned conditions are satisfied, the linearized equation (Equation (14)) is substantially reduced to (see GLV)

$$\begin{eqnarray}&&{\rm{\Gamma }}\displaystyle \frac{\delta {\rho }_{n}}{{\rho }_{n}}={\rm{i}}k\delta {v}_{n},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\rm{\Gamma }}\delta {v}_{n}\approx \gamma {\rho }_{i}{V}_{d}\displaystyle \frac{\delta {\rho }_{i}}{{\rho }_{i}},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&2\beta {\rho }_{i}\displaystyle \frac{\delta {\rho }_{i}}{{\rho }_{i}}\approx {\xi }_{\mathrm{CR}}\displaystyle \frac{{\rho }_{n}}{{\rho }_{i}}\displaystyle \frac{\delta {\rho }_{n}}{{\rho }_{n}}.\end{eqnarray} \tag{ 18 }$$

Along with the background states, the above equations leads to the dispersion relation,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{GLV}}\approx \pm \displaystyle \frac{(1+{\rm{i}})}{2}\sqrt{{{\rm{\Gamma }}}_{{{kV}}_{d}}{\gamma }_{i}}=\pm \displaystyle \frac{(1+{\rm{i}})}{2}\sqrt{\displaystyle \frac{k}{{L}_{B}}}{V}_{A,n}.\end{eqnarray} \tag{ 19 }$$

In the above dispersion relation, the positive and negative parts correspond to the growing and decaying waves, respectively. The subscript GLV is added to Γ to indicate the growth/damping rate (Re[Γ]) and wave frequency (Im[Γ]) of the unstable/decaying mode derived by GLV, which is compared to the eigenvalue of the growing mode in Section 2.4. The growth/damping rate is lower than the speed of the Alfvén wave in the neutrals crossing one wavelength ${{\rm{\Gamma }}}_{\mathrm{alf},n}\equiv {{kV}}_{A,n}$, which indicates that the unstable/decaying mode is slower than the magneto-acoustic mode of the bulk fluid. It is also evident that the above dispersion relation appears very different from that for the magneto-acoustic mode, and therefore the density perturbation is not caused by magneto-acoustic oscillations.³ The reason why the 1D overdensity/underdensity occurs in the bulk of the fluid (i.e., the neutrals) is that the neutrals experience high drag due to the density clump of the ions (Equation (17)), which is also the density clump of the neutrals due to the rapid ionization equilibrium (Equation (18)).

By substituting the positive part of Equation (19) (for an unstable wave) into Equation (16), we have the relation $\delta {v}_{n}\propto \delta {\rho }_{n}\exp ({\rm{i}}3\pi /4)$ and obtain the phase velocity in the comoving frame of the neutrals given by ${v}_{{ph},n}=-\mathrm{Im}[| {{\rm{\Gamma }}}_{{GLV}}| ]/k\lt 0$. These imply that v_n leads δρ_n by a phase of 3π/4, and the unstable wave travels upstream in the rest frame of the neutrals. As explained in GLV, this specific phase difference between δv_n and δρ_n in the rest frame of neutrals is the physical origin of the drag instability. Figure 2 illustrates how the instability occurs. The ion and neutral density perturbations δρ_i and δρ_n are in phase due to the ionization equilibrium. In the rest frame of the neutrals, the wave and V_i propagate upstream (to the left), and δv_n leads δρ_n by a phase difference 3π/4. Consequently, the peak of the density perturbation δρ_n at x = π/2 continues to increase due to the converging velocity field (i.e., dδv_n/dx < 0), while the trough of the density perturbation at x = −π/2 continues to decrease due to the diverging velocity field (i.e., dδv_n/dx > 0), thereby leading to further growth of the perturbations. We refer to the study of GLV for more detailed descriptions.⁴ Note that the unstable wave in fact propagates downstream in the shock frame at the phase velocity given by v_ph,n + V_n: the fast streaming motion of the background flow brings the growing wave downstream through the shock.

Zoom In Zoom Out Reset image size

We now examine the dispersion relations outside a C-shock where V_i = V_n. Thus, V_d = 0 and ${{\rm{\Gamma }}}_{i}\equiv {\rm{i}}\omega +{\rm{i}}{{kV}}_{i}={\rm{\Gamma }}$ in both the pre- and postshock regions. Hence, the dispersion relation derived from Equation (14) reads (e.g., GLV)

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Gamma }}}_{\mathrm{th}}^{4}{\rm{\Gamma }}+{\rm{\Gamma }}{{\rm{\Gamma }}}_{\mathrm{re}}[{{\rm{\Gamma }}}_{\mathrm{alf},i}^{2}({\rm{\Gamma }}+{\gamma }_{i})+{{\rm{\Gamma }}}^{2}{\gamma }_{n}]\\ +\ {{\rm{\Gamma }}}_{\mathrm{th}}^{2}\left[{{\rm{\Gamma }}}_{\mathrm{alf},i}^{2}{{\rm{\Gamma }}}_{\mathrm{re}}+{\rm{\Gamma }}(2{\rm{\Gamma }}+{\gamma }_{i})\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{re}}}{2}+{\rm{\Gamma }}{{\rm{\Gamma }}}_{\mathrm{re}}{\gamma }_{n}\right]=0,\end{array}\end{eqnarray} \tag{ 20 }$$

where ${{\rm{\Gamma }}}_{\mathrm{alf},i}\equiv {{kV}}_{A,i}$ and ${\gamma }_{n}\equiv \gamma {\rho }_{n}$. The rate Γ_th usually has a low value. It may be even lower than γ_i in the above dispersion relation, which sometimes occurs in the postshock region where the neutral density is compressed and the ion density is subsequently enhanced by ionization. Ignoring the thermal terms associated with Γ_th at this moment, Equation (20) can be reduced as

$$\begin{eqnarray}&&({\rm{\Gamma }}+{{\rm{\Gamma }}}_{\mathrm{re}})[({{\rm{\Gamma }}}_{\mathrm{alf},i}^{2}+{{\rm{\Gamma }}}^{2})({\rm{\Gamma }}+{\gamma }_{i})+{\gamma }_{n}{{\rm{\Gamma }}}^{2}]=0.\end{eqnarray} \tag{ 21 }$$

The first term of the above equation represents a decaying mode with a damping rate of −Γ_re. In the second term of the above equation, we have recovered the well-known dispersion relation for 1D linear Alfvén waves in a weakly ionized plasma (Kulsrud & Pearce 1969).⁵ Two branches of this mode can exist depending on the strength of the coupling between ions and neutrals (i.e., the strong and weak-coupling branches). In the strong-coupling branch, Γ_alf,i ≪ γ_n; thus, the dispersion relation can be further reduced as follows (e.g., Kulsrud & Pearce 1969; McKee et al. 2010):

$$\begin{eqnarray}\begin{array}{rcl}\omega +{{kV}}_{n} & = & \pm {\left({{\rm{\Gamma }}}_{\mathrm{alf},i}^{2}\displaystyle \frac{{\rho }_{i}}{{\rho }_{n}}-\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{alf},i}^{4}}{4{\gamma }_{n}^{2}}\right)}^{1/2}+{\rm{i}}\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{alf},i}^{2}}{2{\gamma }_{n}}\\ & \approx & {\rm{i}}\left[\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{ambi}}}{2}\pm \left(\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{ambi}}}{2}-{\gamma }_{i}\right)\right],\end{array}\end{eqnarray} \tag{ 22 }$$

where Γ_ambi is the ambipolar diffusion rate equal to ${k}^{2}{D}_{\mathrm{ambi}}$, with the ambipolar diffusion coefficient ${D}_{\mathrm{ambi}}\equiv {V}_{A,n}^{2}/\gamma {\rho }_{i}$. We have assumed that Γ_ambi ≫ Γ_alf,n to expand the expression of the square root to further simplify the final result on the right-hand side of Equation (22). On the other hand, the dispersion relation is obtained for the weak-coupling regime (Kulsrud & Pearce 1969; McKee et al. 2010),

$$\begin{eqnarray}\begin{array}{rcl}\omega +{{kV}}_{n} & = & \pm {\left({{\rm{\Gamma }}}_{\mathrm{alf},i}^{2}-\displaystyle \frac{{\gamma }_{n}^{2}}{4}\right)}^{1/2}+{\rm{i}}\displaystyle \frac{{\gamma }_{n}}{2}\\ & \approx & {\rm{i}}\left[\displaystyle \frac{{\gamma }_{n}}{2}\pm \left(\displaystyle \frac{{\gamma }_{n}}{2}-{{\rm{\Gamma }}}_{\mathrm{ambi}}\right)\right],\end{array}\end{eqnarray} \tag{ 23 }$$

where we have assumed that ${\gamma }_{n}^{2}\gg {{\rm{\Gamma }}}_{\mathrm{ambi}}^{2}$, which is equivalent to ${\gamma }_{n}\gg {{\rm{\Gamma }}}_{\mathrm{alf},i}$, to expand the expression of the square root to derive the right-hand side of Equation (23). Thus, Equations(22) and (23) indicate that no waves but decaying modes exist in the frame comoving with the flow. The damping rates in the strong-coupling branch are ∼Γ_ambi and γ_i, and the damping rates in the weak-coupling branch are ∼γ_n and Γ_ambi.

Finally, there exist slow decaying modes associated with the weak thermal effect that we have ignored so far. When we consider that Γ ≪ Γ_alf,i, γ_i and because Γ_th ≪ Γ_ambi, Equation (20) can be reduced to ${{\rm{\Gamma }}}_{\mathrm{th}}^{2}+{{\rm{\Gamma }}}^{2}+{\rm{\Gamma }}{\gamma }_{i}\approx 0$, which has two solutions: ${\rm{\Gamma }}\approx (-{\gamma }_{i}\pm \sqrt{{\gamma }_{i}-4{{\rm{\Gamma }}}_{\mathrm{th}}^{2}})/2$. If ${\gamma }_{i}^{2}\gg 4{{\rm{\Gamma }}}_{\mathrm{th}}^{2}$, the reduced dispersion relation yields the following two decaying modes with no waves in the comoving frame of the flow:

$$\begin{eqnarray}&&\omega +{{kV}}_{n}\approx {\rm{i}}\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{th}}^{2}}{{\gamma }_{i}},{\rm{i}}{\gamma }_{i}.\end{eqnarray} \tag{ 24 }$$

However, if ${\gamma }_{i}^{2}\ll 4{{\rm{\Gamma }}}_{\mathrm{th}}^{2}$, the following two decaying wave modes exist in the comoving frame of the flow:

$$\begin{eqnarray}&&\omega +{{kV}}_{n}\approx {\rm{i}}\displaystyle \frac{{\gamma }_{i}}{2}\pm {{\rm{\Gamma }}}_{\mathrm{th}}.\end{eqnarray} \tag{ 25 }$$

In summary, the dispersion relation in the postshock region (Equation (20)) indicates the presence of five decaying modes with damping rates of γ_n, 2βρ_i, Γ_ambi, γ_i, and ${{\rm{\Gamma }}}_{\mathrm{th}}^{2}/{\gamma }_{i}$ when Γ_th < γ_i or γ_n, 2βρ_i, Γ_ambi, γ_i/2, and γ_i/2 (the same as for the two slowest modes) when Γ_th > γ_i.

Figure 3 illustrates the suitability of the assumptions made in our fiducial C-shock model, where k is set to 1/0.015 pc⁻¹; the assumptions were made to derive the dispersion relations inside and outside the C-shock. The left panel of Figure 3 indicates that the ion–neutral drift rate across one wavelength, ${{\rm{\Gamma }}}_{{{kV}}_{d}}$, becomes high within the C-shock (i.e., between x = 0 and x ≈ 0.4 pc) because of the high drift velocity V_d that is caused by the shock compression. Consequently, ${{\rm{\Gamma }}}_{{{kV}}_{d}}$ is higher than γ_i but still lower than Γ_re and γ_n inside the C-shock. The drag instability thus occurs inside the shock according to the dispersion relation presented in Equation (19), with the growth rate Re[Γ_GLV] being higher than γ_i and Γ_th. Although we plot Γ_GLV beyond x ≈ 0.4 pc in the postshock region, the instability is expected to disappear because ${{\rm{\Gamma }}}_{{{kV}}_{d}}$ decreases quickly outside the C-shock. We also plot the rate of gravitational instability ${{\rm{\Gamma }}}_{\mathrm{grav}}$, which is considerably lower than the other rates and can be reasonably neglected in the linear analysis.

Zoom In Zoom Out Reset image size

The right panel of Figure 3 indicates that ${\gamma }_{n}\gg {{\rm{\Gamma }}}_{\mathrm{alf},i}\,\gg {{\rm{\Gamma }}}_{\mathrm{ambi}}\gg {{\rm{\Gamma }}}_{\mathrm{alf},n}$, which ensures the presence of the decaying modes described by the dispersion relations with V_d = 0, i.e., Equations (22) and (23), in the postshock region. We also see from the figure that Γ_th is the lowest rate of the rates of interest. It can be ignored except when Γ ≲ Γ_th, which results in the presence of two modes of the lowest decaying rate associated with γ_i and Γ_th in the postshock region, as described by Equations (24) and (25). In the fiduical model, the left panel of Figure 3 shows that γ_i > Γ_th in the postshock region. Hence, the decaying mode is expected to follow the dispersion relation described by Equation (24) more closely (see the next subsection).

After identifying the modes and determining their underlying physics through the simple dispersion relations, we study the exact solutions of the entire set of linearized equations in Equation (14). Figure 4 shows the properties of the eigenvalues, which describe the behaviors of the growth/damping rates (left panel) and the wave frequencies for the modes of interest (right panel). The left panel of Figure 4 overplots the five eigenvalues Im[$| \omega | ]$ at four different locations (colored dots), both inside (x = 0.1, 0.2, 0.3 pc) and outside (x = 0.5 pc) the C-shock, with relevant rates in the system. Of the five eigenmodes, only one unstable wave mode exists inside the C-shock (red dots). In contrast, there are only decaying wave modes (blue dots) in the postshock region (the x = 0.5 pc location).

Zoom In Zoom Out Reset image size

To understand the properties of these eigenvalues, the left panel of Figure 4 also shows the rates relevant to the different decaying modes according to the dispersion relations in the previous subsection. We see that the first three largest Im[$| \omega |$] are almost the same as γ_n, Γ_ambi, and Γ_re and therefore correspond to the the damping processes due to ion–neutral collisions, ambipolar diffusion, and recombination, respectively. Although the result is expected from the dispersion relations for the postshock region, the left panel of Figure 4 suggests that the aforementioned three damping modes also exist inside the shock. In addition, Figure 4 also illustrates that the two smallest Im[$| \omega | ]$ at x = 0.1, 0.2, and 0.3 pc are consistent with Re[${{\rm{\Gamma }}}_{{GLV}}]$ for the pair of growing and decaying wave modes.

Furthermore, the left panel of Figure 4 shows that the last two smallest Im[$| \omega |$] are close to γ_i and ${{\rm{\Gamma }}}_{\mathrm{th}}^{2}/{\gamma }_{i}$ at x = 0.5 pc; thus these eigenvalues are consistent with the dispersion relations for γ_i > Γ_th in the postshock region, as expected from the previous subsection. When k is increased such that γ_i ≈ Γ_th or even γ_i ≪ Γ_th, the decaying rate of the modes with the last two smallest Im[ω] becomes close to γ_i/2 in the postshock region (not shown), in accordance with the dispersion relation described by Equation (25) or the more general form ${\rm{\Gamma }}=-{\gamma }_{i}\pm \sqrt{{\gamma }_{i}^{2}-4{{\rm{\Gamma }}}_{\mathrm{th}}^{2}}/2$ shown in Section 2.3. Combined with the fact that these two modes represent the pair of the growing and decaying wave modes within the C-shock (see Equation (19)), this suggests that these two decaying modes in the postshock region replace the overstable mode and its counterpart of decaying mode inside the shock. We trace the evolution of the eigenvalue and eigenmode from the in-shock to the postshock region around x ≈ 0.4 pc. We find that when the overstable mode propagates to the postshock region, it gradually transforms into a decaying mode with the damping rate of $\approx {{\rm{\Gamma }}}_{\mathrm{th}}^{2}/{\gamma }_{i}$ in our fiducial case.

Figure 5 illustrates the physical properties of the unstable mode as a function of x in the comoving frame of the neutrals. The left panel of Figure 5 displays the growth rate Γ_grow ($=-\mathrm{Im}[\omega ]\gt 0$) and wave frequency ${\omega }_{\mathrm{wave},n}$ ($=\mathrm{Re}[\omega ]+{{kV}}_{n}$) of the unstable mode in the comoving frame of the neutrals. The growth rate Re[$| {{\rm{\Gamma }}}_{{GLV}}| ]$ is also plotted in this panel for comparison. According to the simplified dispersion relation in Equation (19), Re$[| {{\rm{\Gamma }}}_{{GLV}}| ]={\omega }_{\mathrm{wave},n}={{\rm{\Gamma }}}_{\mathrm{grow}}$. These parameters are not completely identical for the exact solutions shown in left panel of Figure 5, however; ${\omega }_{\mathrm{wave},n}$ and Re$[| {{\rm{\Gamma }}}_{{GLV}}| ]$ are marginally larger than Γ_grow due to the presence of less dominant terms that are not significantly smaller than Re[$| {{\rm{\Gamma }}}_{{GLV}}|$], such as γ_i and Γ_th (left panel of Figure 3). Owing to the same reason, the phase difference between δv_n and δρ_n of the unstable model is not exactly (3/4)π, as expected from the dispersion relation, but approaches this toward ≈0.85π from the shock boundaries to the middle of the shock width, as depicted in the right panel of Figure 5. The same panel also shows that δρ_n and δρ_i are almost in phase due to the ionization-recombination equilibrium. When we remove the ionization and recombination terms in the linearized equations, the unstable mode almost disappears in the eigenvalue problem. Consequently, the overall results are consistent with the results expected from the simple dispersion relation for the drag instability. The physical picture is that in the comoving frame of the neutrals, the ions drift toward the shock front at x = 0 pc (V_d < 0), and the wave travels toward the shock front as well (ω_wave,n > 0) with a phase of δv_n, which leads δρ_n by approximately (3/4)π.

Zoom In Zoom Out Reset image size

Because the overstable mode propagates inside the shock and subsequently decays in the postshock region, a problem arises regarding whether sufficient time is available for the unstable wave to grow. The right panel of Figure 4 shows the smallest Im[$| \omega | ]$ (dots) and its corresponding wave frequency Re[−ω] (crosses) inside the C-shock at x = 0.1, 0.2, and 0.3 pc and in the postshock region at x = 0.5 pc for our fiducial model with k = 1/(0.015 pc). The ratio Re[−ω]/Im[$| \omega |$] is approximately 10–20, which indicates that unstable waves travel downstream approximately 1020 times faster than their growth rate. This high rate of wave propagation in the shock frame is caused by the advection of waves by the fast downstream motion with a speed of V_n, which dominates v_ph,n of the unstable mode. The aforementioned statement is verified by the result Re[−ω] ≈ kV_n, as presented in the right panel of Figure 4. In particular, Re[−ω] is exactly equal to kV_n in the postshock region, which agrees with the dispersion relations. In the following section, we investigate whether the shock width is sufficient or if any favorable preshock conditions exist for the fast-traveling wave to grow significantly.

In the preceding WKBJ analysis, we simply kept k constant in the fiducial case to study the basic properties of a local unstable/decaying mode. To examine whether an unstable wave mode can grow appreciably over the shock width, we consider a mode of a given wave frequency of ω_wave = Re[ω] in the shock frame (e.g., ∼−1 × 10⁻¹¹ s⁻¹ according to the right panel of Figure 4). Equation (14) is solved for both the growth rate Γ_grow (≡Im[ω] when Im[ω] < 0 or zero otherwise) and the wavenumber k, corresponding to a given ω_wave everywhere in the shock. The perturbation amplitude $| U|$ is arbitrary in a linear analysis for normal modes. For our purpose of evaluating the global growth, we can gain a general sense of the total growth of the unstable model U by setting its norm equal to 1 everywhere in the shock. The local growth of the unstable wave mode is $\exp [{{\rm{\Gamma }}}_{\mathrm{grow}}({dx}/{v}_{\mathrm{ph}})]$, where v_ph is the phase velocity of the wave in the shock frame and is equal to −ω_wave/k. Consequently, the total growth of the mode can be computed by integrating the local growth over the entire shock width (i.e., $\exp ({\int }_{\mathrm{shock}\mathrm{width}}{{\rm{\Gamma }}}_{\mathrm{grow}}{dx}/{v}_{\mathrm{ph}})$). We vary the wave frequency ω_wave and repeat the above procedure to seek the particular mode with the maximum total growth (MTG). MTG=1 when no growth occurs (i.e., Γ_grow = 0 everywhere).

In our fiducial model for the steady C-shock, the mode ω_wave ≈ −3 × 10⁻¹¹ s⁻¹ is responsible for the MTG. The mode properties as a function of x are shown in Figure 6. The figure indicates that the growth rate Γ_grow increases from the shock front around x ≳ 0.027 pc, declines after the midpoint of the shock width, and then drops to zero at the end of the shock width at approximately x = 0.37 pc. The parameter Γ_grow is more than 10 times smaller than ω_wave. Nevertheless, the wavenumber k increases downstream across the shock width due to the gradient of the background states. The result can be explained by the relation Re[−ω] ≈ kV_n, as discussed in the preceding section. Because the wave frequency ω_wave (=Re[ω] ≈ − kV_n) is approximately constant for a given wave mode, V_n decreases and hence k increases with x as the neutrals are compressed and thus decelerated downstream inside the shock.

Zoom In Zoom Out Reset image size

The left panel of Figure 7 indicates that 1/(kL_B) and 1/(kL_p) of the unstable mode are considerably smaller than 1, which justifies the WKBJ approximation for the calculations. The right panel of Figure 7 displays the profile of the exponential exponent of the local growth (${{\rm{\Gamma }}}_{\mathrm{grow}}{dx}/{v}_{\mathrm{ph}}$) across the shock. The profile indicates that the unstable mode gains more growth in the rear of the shock transition because the mode has a larger k farther downstream within the shock and hence propagates more slowly to allow for more growth. This maximum-growth mode caused by the drag instability results in an MTG value of approximately 9.9 within the steady C-shock, which implies that an initial perturbation of finite magnitude (i.e., δρ_n/ρ_n ∼ 1/10) is required for a substantial growth of the mode to the nonlinear regime.

Zoom In Zoom Out Reset image size

In addition to the fiducial case, we also calculate the MTG for the C-shock models listed in Table 1 of CO12, which represent various conditions for C-shocks to form in star-forming clouds. The results of the parameter study are presented in Table 2. In accordance with CO12, the letters N, V, B, and X of the model names denote the variations in the n₀, v₀, B₀, and χ_i0 values of the models, respectively. In this study, the variation in χ_i0 is simply calculated by changing the recombination rate β while maintaining the ionization rate ξ_CR constant. Table 2 is almost identical to Table 1 of CO12, except for the last two columns, which present the MTG and the corresponding ω_wave of the unstable mode in each model.

Table 2.  MTG and Corresponding Mode Frequency ${\omega }_{\mathrm{wave}}$ of the Drag Instability in the Steady C-shock Models of CO12

Model	n0 (cm−3)	v0 (km s−1)	B0 (μG)	${\chi }_{i0}$	Lshock (pc)	ωwave (s−1)	MTG
N01	100	5	10	5	3.15	−2e−11	10.9
N03	300	5	10	5	1.38	−3e−11	19.1
N05	500	5	10	5	0.94	−3e−11	22.0
N08	800	5	10	5	0.66	−4e−11	23.8
N10	1000	5	10	5	0.56	−4e−11	24.7
V04	200	4	10	5	1.68	−3e−11	6.4
V06	200	6	10	5	2.05	−2e−11	36.0
V08	200	8	10	5	2.37	−2e−11	121.2
V10	200	10	10	5	2.65	−4e−11	327.3
V12	200	12	10	5	2.90	−3e−11	816.2
B02	200	5	2	5	0.84	−2e−11	27.7
B04	200	5	4	5	1.18	−2e−11	25.3
B06	200	5	6	5	1.45	−2e−11	22.5
B08	200	5	8	5	1.68	−2e−11	19.5
B10	200	5	10	5	1.87	−2e−11	16.6
B12	200	5	12	5	2.05	−2e−11	13.7
B14	200	5	14	5	2.22	−2e−11	11.2
X01	200	5	10	1	9.37	−2e−11	23.5
X06	200	5	10	6	1.56	−2e−11	14.7
X10	200	5	10	10	0.94	−2e−11	8.9
X15	200	5	10	15	0.62	−2e−11	5.1
X20	200	5	10	20	0.47	−1e−11	3.6
Fig4CO12	200	1	2	10	0.18	−8e−12	1.05
Fig5CO12	500	1	4	10	0.13	N/A	1

Note. The final two rows correspond to the two additional models used for Figures 4 and 5 in CO12, respectively. The parameter L_shock is the shock width estimated using Equation (42) in CO12.

Download table as:  ASCIITypeset image

Table 2 indicates that for the unstable modes with the MTG, ω_wave has a low model dependence and has the value of approximately 1–4 × 10⁻¹¹ s⁻¹. In Table 2, the MTG exhibits an increasing trend with the increasing n₀ and v₀ but decreasing B₀ and χ_i0. A wider shock width L_shock does not necessarily result in a larger MTG. The MTG changes by approximately 10-20 times due to the variations in n₀, B₀, and χ_i0 in the parameter study, resulting in a modest mode growth as in the fiducial model. By comparison, the MTG is more sensitive to the change in v₀. In the parameter study, the stronger shocks characterized by higher shock speeds (v₀ = 8–12 km s⁻¹) with wider shock widths (L_shock ≳ 2 pc), as presented in Models V08, V10, and V12, can boost the MTG to a value of several hundred.

In addition to considering the original models in Table 1 of CO12, we ran two additional models, which are denoted as "Fig4CO12" and "Fig5CO12" in the last two rows of Table 2. These two models correspond to the scenarios shown in Figures 4 and 5 of CO12 for their 1D MHD simulations, which represent the interesting cases of weak C-shocks with v₀ = 1 km s⁻¹. The MTG values of these two cases are almost 1, meaning little to no growth of the initial perturbation for weak C-shocks.

The overall trend of the change in the MTG with the preshock parameters n₀, v₀, B₀, and χ_i0 may be qualitatively but not quantitatively understood as follows: The exponent of the MTG (i.e., ${\int }_{\mathrm{shock}\mathrm{width}}{{\rm{\Gamma }}}_{\mathrm{grow}}{dx}/{v}_{\mathrm{ph}})$) can be approximated to t_ad Γ_grow, where t_ad is the ambipolar drift timescale across the shock width, which is equal to L_shock/V_d. Inside the shock, the ions are compressed prior to the neutrals. Thus, we consider r_n ∼ 1 and r_B ∼ r_f, where r_f is the final compression ratio, which is proportional to v₀/V_A,n,0 (CO12). Using the dispersion relation ${{\rm{\Gamma }}}_{\mathrm{grow}}\sim \sqrt{k/{L}_{B}}{V}_{A,n}\approx \sqrt{{\omega }_{\mathrm{wave}}/({L}_{B}{V}_{n})}{V}_{A,n}$ as well as L_B ∼ L_shock, V_d ∼ V_n ∼ v₀/r_f, and ${V}_{A,n}\propto {{Bn}}_{n}^{-1/2}\propto ({r}_{B}{B}_{0}){({r}_{n}{n}_{0})}^{-1/2}$, we obtain the following relations: V_A,n ∼ v₀ and V_d ∼ V_A,n,0. Because ${L}_{{shock}}\propto {n}_{0}^{-3/8}{v}_{0}^{1/4}{B}_{0}^{1/4}{\chi }_{i0}^{-1/2}$ (CO12), it follows that ${t}_{{ad}}{{\rm{\Gamma }}}_{\mathrm{grow}}\propto {L}_{{shock}}^{1/2}{V}_{d}^{-3/2}{v}_{0}\propto {n}_{0}^{3/8}{v}_{0}^{5/4}{B}_{0}^{-5/4}{\chi }_{i0}^{-1/2}$. The scaling result for the exponent qualitatively gives the trend that MTG increases with n₀ and v₀, but decreases with B₀ and χ_i0. In terms of a rough physical picture, the ion–neutral drift velocity, which is comparable to the neutral Alfvén velocity in the preshock region, is enhanced by strong magnetic fields or a low neutral density, leading to a short ambipolar drift time t_ad for the unstable wave to grow. Hence, a positive correlation exists between the MTG and the neutral density n₀, whereas a negative correlation exists between the MTG and the magnetic fields B₀. In addition, the neutral Alfvén speed in the shock is enhanced in a strong shock to result in a high growth rate; that is, a positive correlation between the MTG and the shock speed v₀. Finally, χ_i0 changes the MTG by varying the shock width. A large χ_i0 leads to a wide shock width, which allows more time for the unstable mode to grow. We discuss the astronomical implication of these results in Section 4.2.

While the model of the drag instability was developed based on the steady-state profile of C-shocks, it is possible that the drag instability could occur in time-dependent simulations of C-shocks. Conceptually, this could occur when the evolving C-shock system is very close but not exactly equal to the steady-state C-shock structure. If the deviation from the steady-state profile happens to satisfy the unstable mode favored by the drag instability, this local perturbation could evolve and grow with time.

In addition to deriving the structure of a steady-state C-shock, CO12 also numerically obtained the C-shock structure by simulating two colliding flows using the Athena code (Stone et al. 2008). Instead of computing two fluids comprising the ions and neutrals, as shown in Equations (1)–(5), CO12 simulated the equations for the neutrals alone under the strong-coupling approximation, that is, ${{\boldsymbol{f}}}_{{\boldsymbol{d}}}={{\boldsymbol{f}}}_{{\boldsymbol{L}}}=(1/4\pi )({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}$. Therefore the two-fluid equations considered in this study can be reduced to the following one-fluid equations for the neutrals (see also, e.g., Mac Low et al. 1995):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{n}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{n}{{\boldsymbol{v}}}_{{\boldsymbol{n}}})=0\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\rho }_{n}\left[\displaystyle \frac{\partial {{\boldsymbol{v}}}_{{\boldsymbol{n}}}}{\partial t}+({{\boldsymbol{v}}}_{{\boldsymbol{n}}}\cdot {\rm{\nabla }}){{\boldsymbol{v}}}_{{\boldsymbol{n}}}\right]+{\rm{\nabla }}{p}_{n}=\displaystyle \frac{1}{4\pi }({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+{\rm{\nabla }}\times ({\boldsymbol{B}}\times {{\boldsymbol{v}}}_{{\boldsymbol{n}}})\\ \quad =\ {\rm{\nabla }}\times \left\{{\boldsymbol{B}}\times \left[\displaystyle \frac{1}{4\pi \gamma {\rho }_{i}{\rho }_{n}}({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}\right]\right\}.\end{array}\end{eqnarray} \tag{ 28 }$$

Note that the momentum equation (Equations (26)) and the mass conservation equation (Equations (27)) for neutrals are identical to that in the ideal MHD limit (see Equations (1) and (3)), but the induction equation (Equations (5)) now has a correction term from the ion–neutral drift (see also Equation (46) in CO12). Although the numerical results of CO12 were consistent with the analytical expectation for the C-shock structure, no instabilities were observed in their simulations. This result contrasts with the result of our linear analysis. We discuss possible explanations below.

Because the drag instability is derived from the two-fluid model in this study, we examine whether the strong-coupling approximation can dismiss the drag instability in the one-fluid model adopted in Mac Low et al. (1995) and CO12, for example. In the two-fluid model, the drag instability arises from the perturbed drag term $\gamma {\rho }_{i}{V}_{d}\delta {\rho }_{i}/{\rho }_{i}$ in the momentum equation for the neutrals (see Equation (17)). In the one-fluid model, this term is replaced by the perturbation of magnetic pressure $-({V}_{A,n}^{2}/B)d\delta B/{dx}=-{ik}({V}_{A,n}^{2}/B)\delta B$, which is in turn linked to the perturbation of the diffusion-corrected induction equation (Equation (28)):

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Gamma }}\delta B+{\rm{i}}k\delta {v}_{n}B & \approx & \displaystyle \frac{{V}_{A,n}^{2}}{\gamma {\rho }_{i}}\left(\displaystyle \frac{{d}^{2}\delta B}{{{dx}}^{2}}-\displaystyle \frac{1}{{\rho }_{n}}\displaystyle \frac{d\delta {\rho }_{n}}{{dx}}\displaystyle \frac{{dB}}{{dx}}\right)\\ & = & -{k}^{2}\displaystyle \frac{{V}_{A,n}^{2}}{\gamma {\rho }_{i}}\delta B+{\rm{i}}{{kV}}_{d}\displaystyle \frac{\delta {\rho }_{n}}{{\rho }_{n}}B.\end{array}\end{eqnarray} \tag{ 29 }$$

Note that we have adopted the WKBJ approximation but still retained the term with dB/dx due to a large V_d (∝dB/dx) for our interest. The two terms from the right-hand side of Equation (29) arise from the perturbation of the ambipolar diffusion. The first term is the typical diffusion term ${k}^{2}{D}_{\mathrm{ambi}}$ with the ambipolar diffusivity ${D}_{\mathrm{ambi}}\equiv {V}_{A,n}^{2}/(\gamma {\rho }_{i})$. Rearranging the above equation, we find

$$\begin{eqnarray}&&-{\rm{i}}{{kV}}_{A,n}^{2}\displaystyle \frac{\delta B}{B}\approx \gamma {\rho }_{i}{V}_{d}\displaystyle \frac{\delta {\rho }_{n}}{{\rho }_{n}}-\gamma {\rho }_{i}\delta {v}_{n}+{\rm{i}}\displaystyle \frac{1}{k}{\rm{\Gamma }}\gamma {\rho }_{i}\displaystyle \frac{\delta B}{B},\end{eqnarray} \tag{ 30 }$$

where the first term on the right-hand side (i.e., the last term on the right-hand side of Equation (29)) is the term required for the drag instability in the two-fluid model. This suggests that the drag instability can occur in the strong-coupling limit.

Indeed, the ion–neutral drift in some of the simulated C-shocks may not be sufficiently strong to either initiate the drag instability (models Fig4CO12 and Fig5CO12 in Table 2) or generate appreciable growth without an initial perturbation from the background structure (i.e., the steady-state solution). Another possible factor of the missing instability is numerical resolution. The spacial resolution adopted in CO12's 1D simulations is 0.01 pc, which could be too coarse to resolve the drag instability. In fact, as shown in Figure 6, the wavenumber k corresponds to the unstable mode in our fiducial model is ≳500 pc⁻¹ within the C-shock, which suggests that the physical scale of the growing perturbation could be smaller than 0.002 pc.

We further note that the drag instability has not been reported in most of the previous studies investigating the 1D C-shock structure (e.g., Smith & Mac Low 1997; Chieze et al. 1998; Ciolek & Roberge 2002; van Loo et al. 2009). Smith & Mac Low (1997) adopted the so-called frozen-in condition (e.g., Wardle 1990) assuming ion conservation. This means that ionizations and recombinations are neglected, and the ion number density through the C-shock is purely determined by the compression of the magnetic field via the ion conservation equation and the induction equation. Because the dependence of the ion density on the neutral density is essential for the drag instability to occur, it is not surprising that the drag instability was suppressed in their simulations. For works that included microphysics and/or the chemistry of the C-shock system (e.g., Chieze et al. 1998; Ciolek & Roberge 2002; van Loo et al. 2009), the ionization and recombination processes became more complicated and could affect the timescale on which the drag instability grew. We also note that the increased gas temperature from shock compression (up to ∼10²–10³ K) may completely prevent the drag instability (which is derived using isothermal equation of state) in these simulations.

In this work, we study the drag instability in C-shocks with conditions that can arise from clump–clump collisions or cloud-scale supersonic turbulent flows in typical star-forming regions. One of the most tantalizing questions for shocks in this context is whether a shock instability can lead to fragmentation that is subject to gravitational collapse and eventually induce star formation. However, our analysis here focuses on the behavior of the drag instability within the steady-state profiles of C-shocks, which is linear and non-self-gravitating. These linear analyses are therefore not applicable to directly address this issue and predict any nonlinear outcomes.

Nonetheless, the linear theory could still provide indications on the possible consequence of the instability. The values of the MTG listed in Table 2 indicate that under preferred circumstances, a small perturbation from the steady-state solution could lead to large (≳100×) growth. Based on the parameter study, the most favorable condition for drag instability to grow significantly is in strong shocks (high inflow velocities; v₀ ≳ 5 km s⁻¹). This condition is consistent with the typical environment in giant molecular clouds or molecular cloud complexes (velocity dispersion σ_v ∼ 1–10 km s⁻¹; see recent observations in, e.g., Heyer et al. 2009; Miura et al. 2012; Evans et al. 2014; García et al. 2014; Nguyen-Luong et al. 2016). The total growth induced by the drag instability can also be further enhanced by efficient ambipolar diffusion, i.e., weak-coupling between neutrals and ions, and/or relatively low ionization fractions (n_i/n_n ≲ 10⁻⁷). Because the efficiency of nonideal MHD diffusivity is highly dependent on chemical composition and microscopic physical processes, this condition could be typical in some molecular clouds permitted by dust grain properties, e.g., regions with larger grains (Nishi et al. 1991; Nakano et al. 2002). Still, we note that these conditions potentially favored by the drag instability to develop gravitationally unstable structures were derived following the guidance from the 1D linear analysis, and thus may not be applicable to more complex systems.

In addition, we investigate the perturbation amplitude induced by the drag instability obtained from our eigenvalue problem (see also, e.g., Equation (16)), which is illustrated in Figure 8. Because the amplitude ratio between any two of the perturbations (density, velocity, or magnetic field) remains the same in the linear regime as the unstable mode grows, we plot the relative amplitude of the perturbations in ρ_n, v_n, ρ_i, and v_i normalized by the perturbation in B. Figure 8 shows these perturbations for the growing mode within the C-shock in model V06, which has a moderate MTG in our parameter study (see Table 2). It is evident from the figure that the density perturbations (δρ_n/ρ_n, δρ_i/ρ_i) are much larger than both the velocity and magnetic field perturbations (δv_n/V_n, δv_i/V_i, δB/B) everywhere in the shock. Similar results are found in other models, which implies that as the perturbations grow due to the drag instability, the density perturbation would reach the nonlinear phase faster than other perturbations. As a result, the density perturbation driven by the drag instability is dynamically significant, and these unstable density enhancements induced by the growing wave mode within C-shocks could become gravitationally important in turbulent molecular clouds. Further examinations in numerical simulations will help clarify in which scenario the drag instability would become effective during the star-forming process.

Zoom In Zoom Out Reset image size