Revealing the Drag Instability in One-fluid Nonideal Magnetohydrodynamic Simulations of a 1D Isothermal C-shock

We employ the Athena code to study the presence of the drag instability in a 1D isothermal C-shock. Athena is a grid-based Godunov MHD code designed for astrophysical problems (Stone et al. 2008), which has been widely used for studies of galaxy dynamics, star formation, and accretion disks. In particular, the ambipolar diffusion effect in Athena was developed and tested by Bai & Stone (2011), and has been adopted in star formation simulations studied by Chen & Ostriker (2014). Athena is therefore an ideal numerical tool for evolving Equations (1)–(3) for the environments of star-forming clouds in our study. To study the problem of the drag instability in a C-shock, we make use of the sample test problem provided in Athena, which generates the steady-state C-shock profile following Equations (1)–(3) in Mac Low et al. (1995), and implement the ionization-recombination equilibrium via a power-law dependence of ρ_i on ρ_n , i.e., ${\rho }_{i}\propto {\rho }_{n}^{0.5}$.

Chen & Ostriker (2012) analytically derived the steady-state profile of 1D isothermal C-shocks under the typical environment of star-forming clouds. We use such an equilibrium solution to construct the base-state C-shock profile as the initial condition of our 1D simulations. Since the drag instability is dominated by density perturbation (GC20), we introduce density fluctuations on top of the background state in the incoming inflow of the pre-shock region, and let them propagate downstream at the inflow velocity v₀. We set the density perturbation to be in the form of sinusoidal waves with given amplitude (≪1) and wavenumber k₀. To maintain the amplitude and wavenumber of the perturbation in the pre-shock region, we continuously drive the sinusoidal waves in density with angular frequency v₀ k₀ at the inflow boundary during the simulation. That is, we set the density in the upstream ghost zones to follow the same sinusoidal wave pattern based on the simulation time and the base velocity. The simulation is then evolved with Equations (1)–(3) but without self-gravity, to focus on the instability growth within the C-shock. Note that as the perturbation travels downstream through the C-shock, these sinusoidal waves retain their initial ω_wave but change k due to the background gradient.

It can be difficult to compute the growth rate of the drag instability directly from a nonideal MHD simulation because the unstable mode is expected to travel a few wavelengths over one growth timescale in the shock frame and because its growth rate and the wavenumber vary with x due to the background gradient (GC20). However, we can attempt to extract the local growth rate indirectly from a simulation by using the linear equation. Specifically, the linearized continuity equation for the neutrals reads (i.e., the second row of Equation (4))

$\begin{eqnarray}&&i\omega \delta {\rho }_{n}=-{ik}(\delta {\rho }_{n}{V}_{n}+{\rho }_{n}\delta {v}_{n}).\end{eqnarray} \tag{ 9 }$

We now substitute the identity $\delta {v}_{n}/{V}_{n}=C(x)\delta {\rho }_{n}/{\rho }_{n}\exp (i{\rm{\Delta }}\phi )$ in the above equation, where C(x) =(∣δ v_n ∣/V_n )/(∣δ ρ_n ∣/ρ_n ) and Δϕ is the phase difference between δ v_n and δ ρ_n . Thus, Equation (9) becomes

$\begin{eqnarray}&&\omega =-{{kV}}_{n}(1+C\exp (i{\rm{\Delta }}\phi )),\end{eqnarray} \tag{ 10 }$

which in turn gives the instantaneous growth rate, which is the negative of the imaginary part of ω:

$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{grow}}=-\mathrm{Im}[\omega ]={{CkV}}_{n}\sin ({\rm{\Delta }}\phi ),\end{eqnarray} \tag{ 11 }$

with the wave frequency given by

$\begin{eqnarray}&&{\omega }_{\mathrm{wave}}=\mathrm{Re}[\omega ]=-{{kV}}_{n}[1+C\cos ({\rm{\Delta }}\phi )].\end{eqnarray} \tag{ 12 }$

Similarly, we can also compute the growth rate from the linearized momentum equation for the neutrals (i.e., the third row of Equation (4)):

$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Gamma }}}_{\mathrm{grow}}=-{\omega }_{I}\\ =\displaystyle \frac{{{kV}}_{A,n}^{2}{C}_{B}\sin ({\rm{\Delta }}{\phi }_{B}-{\rm{\Delta }}\phi )-{{kc}}_{s}^{2}\sin ({\rm{\Delta }}\phi )-\gamma {\rho }_{i}{V}_{d}\cos ({\rm{\Delta }}\phi )}{{V}_{n}C},\end{array}\end{eqnarray} \tag{ 13 }$

where we have used the identity $\delta B/B={C}_{B}(\delta {\rho }_{n}/{\rho }_{n})\exp (i{\rm{\Delta }}{\phi }_{B})$ with C_B ≡ (∣δ B∣/B)/(∣δ ρ_n ∣/ρ_n ) and ${\rm{\Delta }}{\phi }_{B}\equiv {\phi }_{\delta B}-{\phi }_{\delta {\rho }_{n}}$.

Due to the different eigenvalues in Equation (4), any initial sinusoidal perturbation is the linear combination of the three linearly independent modes from the eigenvalue problem described by Equation (4). After the damping timescales of the decaying modes inferred from the linear theory, the initial sinusoidal perturbation has propagated into the C-shock and settled to the unstable growing mode in a simulation. At this point, C, C_B , Δϕ, and Δϕ_B of the unstable mode can be directly extracted from the simulation at each x. The background states are also known at each x. Therefore, we can compare the analytical and simulated profiles of the perturbations δ ρ/ρ, δ v_n /V_n , and δ B/B in terms of their relative amplitudes and phase to examine whether they match. If they match reasonably, we can estimate the instantaneous growth rate of the drag instability as a function of x from the simulation at a particular time, using either Equation (11) or (13). The growth rates of density, velocity, and magnetic-field perturbations should be the same during the exponential growth phase, except when the initial perturbations have not yet settled to the pure growing mode driven by the drag instability or when the perturbations have grown significantly to be subject to nonlinear saturation.

In reality, the situation is expected to be more complex because the perturbations are not exactly sinusoidal even over a distance of 2–3 wavelengths (e.g., see Figure 2 later in this paper). A range of wavenumber k, amplitudes (C and C_B ), and phases (Δϕ and Δϕ_B ) could be involved in simulated data, say, due to the truncation of the sinusoidal waves and due to k, amplitudes and phases being a continuous function of x (see Figures 6 and 8 in GC20, and Figures 5 and 6 later in this paper). To compare to sinusoidal modes at a particular x from the linear analysis (see Figure 4 later in this paper), a Fourier transform is performed for the local perturbation from a nonideal MHD simulation in a spatial range covering two wavelengths to obtain the predominant mode in the k space. This k value of the predominant mode is then applied to the linear theory for further investigation of its relevance to the drag instability.

GC20 performed a local linear analysis and found the drag instability in C-shocks. The authors adopted the shock profile described by model Fig3CO12 as a fiducial case, which was presented in Figure 3 in Chen & Ostriker (2012) with the pre-shock conditions n₀ = 500 cm⁻³, v₀ = 5 km s⁻¹, B₀ = 10 μG, and χ_i,0 = 10. The shock width is about 0.4 pc and the total growth is significantly limited by the short time span of the unstable mode within the shock in their fiducial model. Following GC20, we consider this case in this study. GC20 also conducted a parameter study for a range of pre-shock conditions (v₀ ∼ 1–6 km s⁻¹, n₀ ∼ 100–1000 1/cm³, B₀ ∼ 5 − 10 μG, and T = 10 K; see Table 1 in Chen & Ostriker 2012 and Table 2 in GC20), which are typical for star-forming regions in molecular clouds (see e.g., McKee et al. 2010). They found that a stronger shock with a broader shock width promotes the total growth of the drag instability across the shock width. Following the discussion of GC20, we also consider model V06 as a representative case for a wider C-shock, of which the pre-shock conditions are given by n₀ = 200 cm⁻³, v₀ = 6 km s⁻¹, B₀ = 10 μG, and χ_i,0 = 5. In this model, the instability has more total growth in the two-fluid linear analysis because the shock width is about five times larger than that in model Fig3CO12.

We investigate the drag instabilities in C-shocks described by models Fig3CO12 and V06. Figure 1 shows the background states of these two C-shock models, presented in terms of the profile of r_n /r_B , where r_n and r_B are the compression ratio of neutral density and magnetic field, respectively. The simulated C-shock structures generated by Athena following Equations (1)–(3) are highly consistent with those integrated from the ordinary differential equation derived by Chen & Ostriker (2012), thus confirming the valid setup of the background states in these simulations.

Zoom In Zoom Out Reset image size

The simulated C-shocks with density perturbations are presented in Figure 2. The top panel of Figure 2 shows snapshots of the sinusoidal propagation of the perturbations with the gas flow along the +x direction, at simulation time t = 1 Myr for model Fig3CO12 (left) and at t = 1.15 Myr for model V06 (right), respectively. At these moments in both models, the perturbations have propagated across the entire shock width (see Figure 1), where the drag instability is to be revealed. The perturbations (initiated at t = 0 Myr) are maintained by continuously imposing sinusoidal waves to gas density at the inflow boundary (i.e., x = 0) with a relative amplitude ∣δ ρ_n /ρ_n ∣_init = 0.001. Since the drag instability is a local instability, the numerical setup is more consistent with the WKBJ approximation for perturbations with large wavenumbers k₀, especially those that satisfy k₀ ≫ 1/L_shock where L_shock is the C-shock width. We thus choose k₀ = 1/0.003 pc⁻¹ in model Fig3CO12 (L_shock ≈ 0.4 pc) and k₀ = 1/0.015 pc⁻¹ in model V06 (L_shock ≈ 1.5 pc). The grid size of the computation domain is 1/8000 and 4/16,000 pc⁻¹ for models Fig3CO12 and V06, respectively. ⁸

Zoom In Zoom Out Reset image size

As expected from the drag instability in the linear theory, the simulated perturbation grows quickly inside the C-shock (see the top panels of Figure 2) with decreasing wavelength with x due to shock compression (see the middle panels of Figure 2). The amplified perturbation then starts to damp at the beginning of the post-shock region (x ≈ 0.6 pc in model Fig3CO12 and ≈ 2.9 pc in model V06; see Figure 1 and the middle panels of Figure 2). The middle panels of Figure 2 also indicate that the profiles of the velocity perturbation appear to be tilted away from the sinusoidal shape near the end of the shock width in the simulations. It is to be recalled that the phase velocity of the perturbations is close to the background velocity V_n . As the velocity perturbation δ v_n /V_n ≈ δ v_n /v_ph increases with x, the weak nonlinear effect becomes non-negligible. The velocity crest travels even faster than the velocity trough, which slightly steepens the sinusoidal profiles downstream.

Further, the lower panels of Figure 2 show that before the perturbation is amplified within the shock, the wave amplitude goes to zero at particular locations independent of time, i.e., x ≈ 0.125 pc in model Fig3CO12 and x ≈ 0.78 pc in model V06. This simulated result for the nearly zero amplitude can be explained by the destructive interference of two modes, as the initial density perturbation disperses in the linear theory. Solving Equation (4) in the pre-shock region with Re [ω] = ω_wave = −v₀ k₀, we obtain the three decaying modes (i.e., Im [ω] > 0) as illustrated in Figure 3. The first mode, denoted by subscript 1, is predominated by magnetic-field perturbation and is thus subject to magnetic diffusion. This mode is short-lived with the damping timescale given by the ambipolar diffusion timescale 1/(k² D_ambi) ≈ 270 and 850 yr in the linear theory for models Fig3CO12 and V06, respectively. In contrast, the second and the third modes, denoted by the subscripts 2 and 3 in Figure 3, form a pair of modes predominated by density perturbation and are hence subject to the ion-neutral drag, thus being long-lived with the damping timescale $\approx {(\gamma {\rho }_{i}/2)}^{-1}$ (GC20), which are 0.15 and 0.48 Myr for models Fig3CO12 and V06, respectively.

Zoom In Zoom Out Reset image size

The density-only sinusoidal perturbation excited at the left boundary can be expressed by the linear combination of the three aforementioned modes. As explained in detail in Appendix, the pair of the modes, i.e., the second and third modes, are strongly coupled with the initial density perturbation and are therefore predominantly excited with equal amplitude in contrast to the weakly excited first mode. Further, the second and third modes are slowly damped, compared to the short-lived first mode. Consequently, the second and third modes are crucial for the subsequent evolution. This is the reason why a density-only perturbation is excited in the simulations. As has been demonstrated in GC20, one of these two primary modes becomes the growth mode due to the drag instability with the other as its counterpart, yet decaying mode, once the pair of modes travel into the shock, which can be realized in Equation (8) also. ⁹ Owing to the slight difference in the wavelength, the second and third modes produce "wave beats in space" until one of the pair damps and the other grows significantly inside the shock width. In other words, the envelope of the wave beats does not travel with time and its spatial interval λ_beat =2π/(k₂ − k₃) between the locations of complete destructive interference. Using the values of k₂ and k₃ derived from the linear theory as shown in Figure 3, we obtain the beat wavelength λ_beat ≈ 0.24 and 1.43 pc for models Fig3CO12 and V06, respectively. The linear results are consistent with the simulated locations of complete destructive interference in the pre-shock region at x ∼ 0.12 pc in model Fig3CO12 and x ∼ 0.71 pc in model V06 (i.e., ≈ λ_beat/2 from the left boundary), as shown in the bottom panels of Figure 2. In Appendix, we elaborate on the exact locations of zero amplitude of the beat envelope, which match perfectly with the simulation results by also considering the initial phase difference between the mode pair when they are excited at the left boundary. Besides the first zero of the beat envelope, it can be also observed from the bottom panels of Figure 2 that the next minimal amplitude of the wave envelope due to beats occurs inside the shock, where x ≈ 0.32 pc in model Fig3CO2 and 1.75 pc in model V06. The amplitudes are not nil at these locations, as one of the mode pairs has been growing and the other has been decaying within the C-shock.

Figure 4 compares the perturbation profiles of the simulated (solid line) and analytical results (dashed). The comparison is made at the neighborhood with the range of about two wavelengths of three locations within the C-shock in each model, Fig3CO12 (left panels) and V06 (right panels). The three locations are selected in the x range of the simulation, where the perturbations have exhibited clear growth within the C-shock (see Figure 2). In other words, to identify the instability easily, we are looking at the region where the decaying mode has become much weaker than the growing one in the pair of density-dominated modes and thus the total perturbation should be primarily contributed by the growing mode associated with the drag instability. As described earlier, the simulated profile for perturbations in Figure 4 is the predominant Fourier mode extracted from the simulated result. The wavenumber k of the predominant mode obtained from the Fourier transform is applied to the linear theory to find the unstable mode driven by the drag instability. Subsequently, we scale the amplitude of δ ρ_n /ρ_n of the growing mode from the linear theory to overlap that of the predominant mode from the simulation, while keeping the relative amplitude and phase among the perturbations in the linear theory. The procedure yields the comparison plots in Figure 4.

Zoom In Zoom Out Reset image size

Figure 4 illustrates that the simulated curves for perturbations almost match the analytical curves with little perceptible difference, strongly suggesting the presence of the drag instability in the simulation. Hence, we identify C, C_B , Δϕ, and Δϕ_B from the simulated curves in Figure 4 and present the numerical results to compare the analytical solutions in Figures 5 and 6 for models Fig3CO12 and V06, respectively. The amplitudes of δ ρ_n /ρ_n are much larger than those of δ v_n /V_n and δ B/B, i.e., both C and C_B ≪ 1, as has been realized by GC20. Additionally, δ v_n and δ B/B generally lead δ ρ_n by a phase of Δϕ ∼ 0.9π and Δϕ_B ∼ 1.5π, respectively, which yields the instability. Applying Equations (11)–(13), we then obtain ω_wave and Γ_grow at the three different locations in each model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 5 and 6 show that in both models, the analytical and simulated growth rates have a similar trend along the shock flow, though they do not exactly match. The discrepancy is probably caused by the difference between the instantaneous growth rate of the predominant mode and the exact growth rate, i.e., the exact growth rate and wavenumber still slightly vary in the spatial range covering two wavelengths, where the single growth rate and k of the predominant Fourier mode are obtained in our methodology. Nevertheless, the curves for the analytical growth rate are sandwiched by the simulated growth rates evaluated at the three locations using Equations (11) and (13). The wavenumber k, the relative amplitudes (C and C_B ), and phase differences (Δϕ and Δϕ_B ) of perturbations from the simulated results follow the analytical curves reasonably well. The top left panels of Figures 5 and 6 show that ω_wave is about 30–40× larger than Γ_grow for the predominant mode extracted from the simulation, meaning that the instability propagates faster than its growth in accordance with the linear results. Because ω_wave is dominated by the Doppler-shifted frequency kV_n due to the fast downstream flow (GC20), k increases with x due to shock compression as shown in the figures, which has been observed in Figure 2. Further, C slightly increases with x but C_B decreases with x. Overall, the simulated results agree with analytical results in terms of the mode profiles (i.e., the eigenvectors) shown in Figure 4 as well as the mode growth rates and frequencies (i.e., the eigenvalues) plotted in Figures 5 and 6. We confirm the presence of the drag instability in nonideal MHD simulations for 1D isothermal C-shocks.

Besides the confirmation of the linear theory of the drag instability, numerical simulations enable us to directly evaluate the global total growth TG_sim for an initial perturbation across the entire shock width, which cannot be achieved by the local linear theory. Recalling that the perturbations have the form $U\exp ({ikx}+i\omega t)\equiv W$ in the local linear analysis, we may express the amplitude of the growing mode as follows:

$\begin{eqnarray}\begin{array}{rcl}| W(x)| & = & \left(| U(x={x}_{1})| +{\displaystyle \int }_{{x}_{1}}^{x}\displaystyle \frac{d| U| }{{dx}}{dx}\right)\\ & & \times \exp \left({\displaystyle \int }_{{x}_{1}}^{x}{{\rm{\Gamma }}}_{\mathrm{grow}}\displaystyle \frac{{dx}}{{v}_{\mathrm{ph}}}\right),\end{array}\end{eqnarray} \tag{ 14 }$

where x₁ is the location of the beginning of a C-shock and ∣U(x₁)∣ = ∣W(x₁)∣. ∣W(x)∣/∣W(x₁)∣ should equal TG_sim when x is the location of the shock end. However, d∣U∣/dx is unknown because ∣U∣ is arbitrary in a local linear analysis for normal modes. To get a general sense of the total growth from the WKBJ analysis, GC20 set the amplitude of the unstable mode ∣U∣ = 1 everywhere in the shock and thus d∣U∣/dx = 0. The simplification allows the authors to estimate the total growth TG_∣U∣=1 by integrating only the local growth of the mode over the entire shock width, i.e., the exponential term on the right-hand side of Equation (14). The total growth of the unstable mode at x in a shock with ∣U∣ = 1 is then given by

$\begin{eqnarray}&&{\mathrm{TG}}_{| U| =1,x}=\exp \left({\int }_{{x}_{1}}^{x}{{\rm{\Gamma }}}_{\mathrm{grow}}\displaystyle \frac{{dx}}{{v}_{\mathrm{ph}}}\right).\end{eqnarray} \tag{ 15 }$

Therefore, the total growth across the entire shock, denoted by TG_∣U∣=1, is obtained by setting x to be the location of the shock end in the above equation.

GC20 used TG_∣U∣=1 as a proxy of TG_sim to evaluate which shock model favors the growth of the drag instability in their parameter study. As shown above, ∣U∣ should vary with x in the shock and needs to smoothly connect to the initial condition at x = 0. Figure 7 indicates that the profiles of simulated density perturbation δ ρ_n /ρ_n increase faster than TG_∣U∣=1,x , i.e., the total growth from the beginning of the shock based on the linear theory with ∣U∣ = 1. The peak value of TG_∣U∣=1,x in Figure 7 is TG_∣U∣=1. It is to be recalled from Section 4.2 that the amplitude of the unstable mode is made of only one-half of the initial density perturbation ∣δ ρ_n /ρ_n ∣_init = 0.001. The simulated density perturbation increases from 0.0005 to 0.044 in model Fig3CO12 and to 0.057 in model V06, which amount to ${\mathrm{TG}}_{{sim}}\approx 88$ and 114 in model Fig3CO12 and V06, respectively. Comparatively, TG_∣U∣=1 is only about 42.7 in model Fig3CO12 and 58.3 in model V06. Evidently, the simulated perturbations increase faster than those based on the simple estimate by a factor of two in the two models—the local analysis performed at each x cannot account for the growth contributed from d∣U∣/dx in Equation (14). Additionally, the exact growth should be a little larger because the unstable mode has been slowly damped in the pre-shock region before it travels to the C-shock and starts growing.

Zoom In Zoom Out Reset image size

Although the simulated growth rate within the C-shock is generally consistent with the analytical growth rate, they do not exactly match in our study. It is probably because the exact local growth rate and wavenumber still slightly vary in the spatial range covering two wavelengths, where the single growth rate and wavenumber of the predominant Fourier mode are obtained by invoking the instantaneous growth rate. In such a case, the mismatch is expected to only be of order 1/(kL), where L is the gradient length scale. The local analysis yields the dispersion relation described by Equation (8). The analysis can approximately model the effect of a background gradient length scale L by adding a small imaginary part, i/L, to k. The dispersion relation then predicts a correspondingly small shift in ω in the complex plane. However, since ω is mostly real (i.e., ω ≈ kV_n ), a small fractional shift of ω in the complex plane as k is corrected to k + i/L is likely to change the imaginary part of ω by a significant fraction. The dispersion relation is then $\omega \sim {{kV}}_{n}+[(1+i)/2]\sqrt{{{kV}}_{d}\gamma {\rho }_{i}}$, implying a correction to the imaginary part of ω of order iV_n /L. Figure 9 illustrates the error estimate (V_n /L)/Γ for L = the pressure scale length L_p (solid curve) and L = the magnetic-field scale length L_B (dashed curve). They are plotted to compare to the deviation percentage of simulated Γ_grow from that obtained in the complete linear theory as shown in the top left panel of Figures 5 and 6 (i.e., deviation of dots or crosses from the curve for Γ_grow in the top left panel of Figures 5 and 6). We see that the deviation of simulated data from the linear results (∼20%–50%) is in general larger than the error estimate (V_n /L)/Γ ∼ 20% in the second half of the shock regions. The effect of k/L may contribute to part of the deviations. While the error estimate is made based on the dispersion relation Equation (8), it is worth noting that the full linearized equations are more complicated than that. Namely, the result from the complete linear analysis does not exactly follow the dispersion relation due to the presence of other less dominant terms that are dropped to derive the simplified dispersion relation (refer to Section 2).

Zoom In Zoom Out Reset image size

As described in Section 3.2, it can be difficult to directly compare the growth rates based on the local linear analysis to those from numerical simulations of the drag instability, which travels a few wavelengths in one growth timescale. Many of the difficulties in robustly connecting the simulations to the local linear analysis could be avoided by explicitly integrating the linearized equations globally through the steady C-shock background. In contrast to local linear analysis, a global linear analysis preserves information on the causality of propagating waves. Therefore, it would be worthwhile to improve the comparison between the linear analysis and simulated results by performing a global linear analysis of the local, yet fast traveling, overstability.

As has been noted by GC20, the density perturbation dominates over velocity and magnetic-field perturbations in the drag instability, which leads the authors to postulate that clump/core formation could be one of the nonlinear outcomes of the instability in C-shocks. However, our numerical simulations suggest that the growth of the instability in 1D isothermal C-shocks saturates in the weak nonlinear regime where δ ρ_n /ρ_n is only a fraction of unity. Because the corresponding δ v_n /v_ph is even smaller when the nonlinear saturation occurs, we conjecture that the steep pressure gradient of the predominant δ ρ_n /ρ_n would be important in saturating the growth of the drag instability. Analogous to the steepening of 1D ideal MHD waves, the pressure force resulting from the steepened density perturbation is expected to play a role in the nonlinear saturation of the drag instability in 1D C-shocks (Suzuki & Inutsuka 2005; Gu & Suzuki 2009). More numerical simulations for a parameter study with a large initial perturbation should be conducted to fully understand the nonlinear effects.

Larger perturbations can play a more important role in the thermal properties of the modes. The significance of thermodynamics for a perturbation may be characterized by δ v_n /c_s , of which the maximum values are ≈0.2–0.3 for the nonlinear cases presented in Figure 8. Moreover, the radiative cooling rate associated with a mode, Γ_rad, may be approximated as $\delta {\rho }_{n}{c}_{s}^{2}/{\rm{\Lambda }}$, where Λ is the cooling function. Adopting the expressions of Λ from Goldsmith & Langer (1978) due to molecular and atomic line emissions for T ∼ 10 K and n ∼ a few ×10³ cm⁻³, we have Γ_rad ∼ Γ_grow in our nonlinear cases. These all together imply that the adiabatic heating or cooling of the perturbations would be small but non-negligible. In our present study based on isothermal perturbations, heating effects (e.g., due to the ion-neutral drag, nonlinear dissipation) and aforementioned thermal processes are ignored. More self-consistent calculations incorporating thermochemical evolution would be worthwhile for future study of mode saturation.

Although the numerical simulations for the drag instability are conducted for the one-fluid approach in this study, it is worth discussing the difference between the one-fluid approach and the more realistic two-fluid approach from the perspective of local linear theory. As described in Section 4.4, TG_∣U∣=1 = 42.7 and 58.3 in models Fig3CO12 and V06, respectively. However, GC20 identified the maximum-growth mode in the two-fluid approach and showed that the total growth of the maximum-growth mode, i.e., the maximum TG_∣U∣=1 by varying ω_wave, ¹⁰ is only about 10 and 36 for models Fig3CO12 and V06, respectively. Apparently, the growth rate is overestimated in the one-fluid approach. Indeed, we apply the linear analysis in the two-fluid approach and find that TG_∣U∣=1 = 6.4 and 30.3 for model Fig3CO12 and 30.3 for model V06, which are smaller than those in the one-fluid approach by a factor of about 6.7 and 1.9, respectively. We now investigate the underlying physics.

Figure 10 presents the comparison of the growth rate, Δϕ ($={\phi }_{\delta {v}_{n}}-{\phi }_{\delta {\rho }_{n}}$), and ${\phi }_{\delta {\rho }_{i}}-{\phi }_{\delta {\rho }_{n}}$ within the shock width between the one- and two-fluid linear approaches for model Fig3CO12 (left panels) and model V06 (right panels). The curves from the one-fluid theory have been demonstrated to agree with the numerical simulation in Figures 5 and 6. However, the growth rate for the one-fluid approach is larger than that for the two-fluid approach, albeit still in the same order. Given the relatively high ionization and recombination rates in the problem, the primary approximation made for the one-fluid approach is that the ionization-recombination equilibrium is assumed to strictly hold for the density perturbations in the continuity equation for the ions such that $\delta {\rho }_{i}\propto \delta {\rho }_{n}^{1/2}$, which helps simplify the two-fluid model to the one-fluid model. Consequently, the ion and neutral density excesses are always in phase, maximizing the ion-neutral drag, thus promoting the growth due to the drag instability. To further demonstrate the physics more quantitatively, we increase both the ionization and recombination rate by the same factor of 10 and perform the linear analysis. In the process, the background state, which depends only on the ionization fraction, remains unchanged. The dashed lines in Figure 10 show the results for the models with 10 times the ionization and recombination rates. It is evident that in both models, the growth rate in the two-fluid approach increases toward that in the one-fluid approach when both the ionization and recombination rates are enhanced. We can see from Figure 10 that δ ρ_i and δ ρ_n become closer in phase (i.e., ${\phi }_{\delta {\rho }_{i}}-{\phi }_{\delta {\rho }_{n}}\to 0$) due to the even faster processes of the cosmic ionization and ion-neutral recombination in the two-fluid approach. Therefore, the ionization equilibrium is almost attained to increase the phase difference Δϕ and the growth rate closer to those in the one-fluid approach.

Zoom In Zoom Out Reset image size

The importance of the ionization equilibrium in exciting the drag instability can be also realized by ignoring the ionization and recombination effects from the local linearized equations. In the two-fluid approach, the drag instability vanishes in a 1D isothermal C-shock in the absence of ionization and recombination (GC20). In the one-fluid approach, the same physical process can be achieved by using δ ρ_i /ρ_i = δ B/B (frozen-in field condition) instead of δ ρ_i /ρ_i = (1/2)δ ρ_n /ρ_n (ionization equilibrium) in the linearized equations. Consequently, we find that the matrix element (3/2)ikV_d B/ρ_n in C_n becomes ikV_d B/ρ_n in Equation (4). It can then be shown that there is no unstable eigenmode from Equation (4), with the frozen-in background state of 1D isothermal C-shocks where the ion compression ratio equals r_B (Chen & Ostriker 2012).

As k increases in the two-fluid approach, the unstable mode is more akin to an acoustic mode, and the ionization equilibrium becomes less valid. Hence, the drag instability is suppressed by the small-scale gas pressure perturbation along the shock direction (Gu et al. 2004; Gu & Chen 2020; Gu 2021). In contrast, the demise of the instability does not occur for a large k in the one-fluid approach. As the initial k of the perturbation increases in the incoming flow, we find that the total growth TG_∣U=1∣ increases to about 58.15 and 121.06 in models Fig3CO12 and V06, respectively, and then remains almost constant independent of the initial k. It arises because the thermal pressure term ${k}^{2}{c}_{s}^{2}$ becomes important in Equation (6) in the regime of a large k. Retaining this term in Equation (7) yields the simplified dispersion relation for the unstable mode to be Γ ∼ kc_s − i(∣V_d ∣/c_s )γ ρ_i /4. As in the two-fluid approach, the unstable mode of a large k is transformed to an acoustic mode in the one-fluid analysis (i.e., Re [Γ] ∼ kc_s ), but with a growth rate −Im[Γ] independent of k. However, this growth is spurious because the ionization equilibrium breaks down in the limit of a large k, i.e., kc_s > 2β ρ_i , and thus the one-fluid approach fails. Therefore, while the one-fluid analysis is still reasonable to study the general properties of the drag instability for a moderate value of k, it is unable to hint at the maximum-growth modes as has been done by GC20 in their two-fluid analysis.

Adopting the two-fluid approach, Gu (2021) extended the linear analysis of the drag instability in 1D C-shocks to 2D perpendicular and oblique isothermal C-shocks. With one more dimension for fluid motion and magnetic-field component, the author found that the properties of the unstable mode generally depend on the ratio of the transverse (i.e., normal to the shock flow) to longitudinal (i.e., parallel to the shock flow) wavenumber. It follows that the growth rate of the instability is dynamically correlated with ${\phi }_{\delta {\rho }_{i}}-{\phi }_{\delta {\rho }_{n}}$ in different regimes of the wavenumber ratio, as illustrated in panel (a) in Figures 4 and 9 of Gu (2021). Nevertheless, the 2D linear analysis shows that unless the shock is mildly oblique or perpendicular, the maximum growth of the drag instability is probably contributed by the transversely large-scale mode (i.e., almost 1D mode), suggesting that the numerical simulation in the one-fluid approach employed in this study can capture the drag instability in 2D oblique shocks, albeit with overestimated growth rates as obtained in this study. As for the drag instability in mildly oblique and perpendicular shocks, Gu (2021) found a type of unstable modes, which almost comoves with the shock and thus has sufficient time to grow without being limited by the shock width. Since the maximum-growth modes occur when the wavenumber ratio transitions from the transversely large-scale mode (i.e., ionization equilibrium is crucial) to the transversely small-scale mode (i.e., slow modes become essential). The numerical simulation based on the one-fluid approach due to ionization equilibrium may miss this piece of physics and hence fail to capture these modes. More robust studies in the future would shed more light on this issue.

The linear theory, supported by the numerical simulations of this work, suggests that the drag instability imprints the density-banded substructure inside a C-shock in the typical environment of star-forming clouds. The maximum-growth modes have been suggested by GC20 in their two-fluid analysis. That is, for the predominant growing mode, ω_wave ≈ −3 × 10⁻¹¹ s⁻¹ with k ranging from ≈190–1300 pc⁻¹ across the shock for model Fig3CO12, and ω_wave ≈ −2 ×10⁻¹¹ s⁻¹ with k ranging from ≈110–610 pc⁻¹ throughout the shock for model V06. However, the density enhancement in a C-shock is still weaker than the density in the post-shock region, which may imply that the density substructure is probably challenging to probe. For instance, there is no evidence of subcritical clouds in observations (Crutcher 2012), and indeed, Nakano (1998) demonstrated that a subcritical cloud would generally be quite unobservable because it would have too low a column-density contrast with its surroundings. Nevertheless, it is worth mentioning that Hoang & Tram (2019) recently proposed a new probe of the C-shock structure by detecting microwave emissions from fast spinning nanoparticles due to supersonic neutral gas-charged grain drift. In their work, the spin-up timescale of nanoparticles is inversely proportional to the neutral density. While the required resolution and sensitivity would be unrealistically high to resolve the density growth within the C-shock due to the drag instability, the new method could shed new light on observationally probing some plasma properties.