Universal Properties of Dense Clumps in Magnetized Molecular Clouds Formed through Shock Compression of Two-phase Atomic Gases

The basic equations are the same as those used in Paper I except that self-gravity is taken into account in this work. The basic equations are given by

and

where $\rho$ is the gas density, $v_{i}$ is the gas velocity, $B_{i}$ is the magnetic field, $E=\rho v^{2}/2+P/(\gamma-1)+B^{2}/8\pi$ is the total energy density, $\phi$ is the gravitational potential, $T$ is the gas temperature, $\kappa$ is the thermal conductivity, and ${\cal L}$ is a net cooling rate per unit volume. The thermal conductivity for neutral hydrogen is given by $\kappa(T)=2.5\times 10^{3}T^{1/2}$ cm^-1 K^-1 s^-1 (Parker, 1953).

The Possion equation (Equation (5)) is solved by the multigrid method that has been implemented in Athena++ by Tomida (2022, in preparation). We have implemented the chemical reactions, radiative cooling/heating, ray-tracing of FUV photons into a public MHD code of Athena++ (Stone et al., 2020) in Paper I.

We do not use the adaptive mesh refinement technique but use a uniformly spaced grid. The simulations are conducted with resolution of $2048\times 1024^{2}$, which is two times higher than that in Paper I. The corresponding mesh size is $\Delta x\sim 0.02~{}$pc. Although dense cores with $<0.1$ pc are not resolved sufficiently, the surrounding low-density regions, which is important for the MC formation and driving turbulence, are well resolved (Kobayashi et al., 2020).

To avoid gravitational collapse and numerical fragmentation, the cooling and heating processes are switched off if the local thermal Jeans length $\lambda_{\mathrm{J}}\equiv c_{\mathrm{s}}\sqrt{\pi/G\rho}$ is shorter than $4\Delta x$ (Truelove et al., 1997). Since the temperature of the dense gas is about 20 K in our results, the Truelove condition sets the maximum density $n_{\mathrm{max}}\equiv 1.4\times 10^{5}~{}\mathrm{cm}^{-3}(N_{x}/2048)^{2}(c_{\mathrm{s}}/0.3~{}\mathrm{km~{}s^{-1}})^{-2}$ above which the gas is artificially treated as the adiabatic fluid with $\gamma=5/3$.

The maximum density $n_{\mathrm{max}}$ can be reached only when self-gravity becomes important. As will be explained in Section 2.3, we consider colliding flows of the two-phase atomic gas with collision speeds of $V_{0}=5~{}\mathrm{km~{}s^{-1}}$ and $V_{0}=10~{}\mathrm{km~{}s^{-1}}$. The cold phase of the atomic gas is as dense as $n_{\mathrm{cold}}\sim 50~{}\mathrm{cm}^{-3}$. The maximum densities obtained by the shock compression are estimated to be $n_{\mathrm{cold}}(V_{0}/0.3~{}\mathrm{km~{}s^{-1}})^{2}\sim 10^{4}~{}\mathrm{cm}^{-3}$ for $V_{0}=5~{}\mathrm{km~{}s^{-1}}$ and $\sim 4\times 10^{4}~{}\mathrm{cm}^{-3}$ for $V_{0}=10~{}\mathrm{km~{}s^{-1}}$, both of which are smaller than $n_{\mathrm{max}}$.

Paper I took into account 9 chemical species (H, H₂, H⁺, He, He⁺, C, C⁺, CO, $e$) following Inoue & Inutsuka (2012). We improved the chemical networks and thermal processes following Gong et al. (2017). We consider 15 chemical species (H, H₂, H⁺, H${}_{2}^{+}$, H${}_{3}^{+}$, He, He⁺, C, C⁺, CH_x, CO, O, OH_x, HCO⁺, $e$) which are similar to those in Nelson & Langer (1999). We calculate the far-ultraviolet (FUV) flux, which causes photo-chemistry and related thermal processes, is required to be calculated at any cells. In the same way as in Paper I and Inoue & Inutsuka (2012), FUV radiation is divided into two rays which propagate along the $x$-direction from the left and right $x$-boundaries because the compression region has a sheet-like configuration, and the column densities along the $x$-direction are smaller than those in the $y$- and $z$-directions. The two-ray approximation is also adapted for calculating the escape probabilities of the cooling photons. The detailed description is found in Paper I.

The cosmic-ray ionization rate per hydrogen nucleus $\xi_{\mathrm{H}}$ varies significanly depending on regions. In this paper, we adopt the recently updated one in diffuse molecular clouds $\xi_{\mathrm{H}}=2\times 10^{-16}~{}$s^-1 (Indriolo et al., 2007, 2015; Neufeld & Wolfire, 2017), which is about ten times larger than that adopted in Paper I. We should note that the cosmic ray flux is attenuated toward the deep interior of molecular clouds, and Neufeld & Wolfire (2017) showed marginal evidence that $\xi_{\mathrm{H}}$ decreases with $A_{\mathrm{V}}$ in proportional to $A_{\mathrm{V}}^{-1}$ for $A_{\mathrm{V}}>0.5$. Our adopted value of $\xi_{\mathrm{H}}=2\times 10^{-16}~{}$s^-1 corresponds to an upper limit of the cosmic-ray ionization rate, and the gas temperature at dense regions may be slightly overestimated.

As a solver of the stiff chemical equations, we adopt the Rosenbrock method, which is a semi-implicit 4th-order Runge-Kutta scheme (Press et al., 1986). The time width $\Delta t_{\mathrm{chem}}$ required to solve the chemical reactions with sufficiently small errors can be much smaller than that of the MHD part. We integrate the chemical reactions with adaptive substepping, and the number of the substeps is determined so that a error estimated from 3rd-order and 4th-order solutions becomes sufficiently small.

Before presenting the models in this study, we here review the results of Paper I that investigated the early stage of molecular cloud formation from atomic gas by ignoring self-gravity. They performed a survey in a parameter space of ($\langle n_{0}\rangle$, $V_{0}$, $B_{0}$, $\theta$), where $\langle n_{0}\rangle$ is the mean density of the pre-shock atomic gas, $V_{0}$ is the collision velocity, and $B_{0}$ is the field strength, and $\theta$ is the angle between the collision velocity and magnetic field of the atomic gas. At a given parameter set of $(\langle n_{0}\rangle,B_{0},V_{0})$, they found that there is a critical angle $\theta_{\mathrm{cr}}$ which characterizes the physical properties of the post-shock layers. They obtained an analytic formula of the critical angle as follows:

• •

  For an angle smaller than $\theta_{\mathrm{cr}}$, anisotropic super-Alfvénic turbulence is maintained by the accretion of the two-phase atomic gas. The ram pressure dominates over the magnetic stress in the post-shock layers. As a result, a greatly-extended turbulence-dominated layer is generated.

• •

  For an angle comparable to $\theta_{\mathrm{cr}}$, the shock-amplified magnetic field weakens the post-shock turbulence, making the post-shock layer denser. The mean post-shock density has a peak around $\theta\sim\theta_{\mathrm{cr}}$ as a function of $\theta$.

• •

  For an angle larger than $\theta_{\mathrm{cr}}$, the magnetic pressure plays an important role in the post-shock layer. Their mean densities decrease with $\theta$ and are determined by the balance between the post-shock magnetic pressure and pre-shock ram pressure (see also Inoue & Inutsuka, 2012).

We consider colliding flows of the two-phase atomic gas whose mean density is $\langle n_{0}\rangle$ in the same manner as in Paper I. The collision speed is $V_{0}$, magnetic field strength is $B_{0}$ and the angle between the colliding flow and magnetic field is $\theta$. The colliding flows move in the $\pm x$ directions, and the magnetic fields are tilted toward the $y$-axis.

In order to generate the initial conditions for performing the colliding-flow simulations, the upstream two-phase atomic gas is generated using the same manner as in Paper I. We prepare a thermally unstable gas in thermal equilibrium with a mean density of $\langle n_{0}\rangle$ and a uniform magnetic field of $\mathbf{B}_{0}$ in a cubic simulation box of $-10~{}\mathrm{pc}\leq x,y,z\leq 10~{}\mathrm{pc}$ with the periodic boundary conditions in all the directions. As the initial condition, we add a velocity dispersion having a Kolmogrov spectrum of $0.64~{}\mathrm{km~{}s^{-1}}~{}(L/~{}\mathrm{pc})^{0.37}$, where $L$ is the size of regions. The initial unstable gas evolves into the two thermally stable states, or the CNM and WNM, in about one cooling timescale. We terminate the simulations at $t=8~{}$Myr or $\sim 4$ cooling timescale of the initial unstable gas. The density distribution of our upstream gas is highly inhomogeneous, where the dense CNM clumps are embedded in the diffuse WNM (see Figure 1 in Paper I).

In the colliding flow simulations, the simulation box is doubled in the $x$-direction, or $-20~{}\mathrm{pc}\leq x\leq 20\mathrm{pc}$, $-10~{}\mathrm{pc}\leq y,z\leq 10~{}\mathrm{pc}$. In the $y$- and $z$-directions, periodic boundary conditions are imposed. The $x$-boundary conditions at $x=\pm 20~{}$pc are imposed so that the initial distribution is periodically injected into the computation region from both the $x$-boundaries with a constant speed of $V_{0}$.

The model parameters are listed in Table 1. In all the models, $\langle n_{0}\rangle=10~{}\mathrm{cm}^{-3}$ is adopted. As in Paper I, we here consider the formation of MCs through accretion of HI clouds with a mean density of $\sim 10~{}\mathrm{cm}^{-3}$ (Blitz et al., 2007; Fukui et al., 2009; Inoue & Inutsuka, 2012; Fukui et al., 2017) rather than through accretion of the warm neutral medium. For all the models, the field strength $B_{0}$ is set to $3~{}\mu$G, which is roughly consistent with a typical field strength in the atomic gas (Heiles & Troland, 2005; Crutcher et al., 2010).

As a fiducial model in this paper, we consider $V_{0}=5~{}\mathrm{km~{}s^{-1}}$. In order to investigate the dependence of the results on the collision velocity, we perform additional simulations with $V_{0}=10~{}\mathrm{km~{}s^{-1}}$. Our simulations model accumulation of the atomic gas by expansion of HI shells or spiral shock waves (McCray & Kafatos, 1987; Wada et al., 2011). The critical angle for the fiducial model is $\sin\theta_{\mathrm{cr}}=0.1$ for $V_{0}=5~{}\mathrm{km~{}s^{-1}}$ and $0.2$ for $V_{0}=10~{}\mathrm{km~{}s^{-1}}$ (Equation (7)).

For each $V_{0}$, we consider cases with $\theta=0.25\theta_{\mathrm{cr}}$ and $\theta=2\theta_{\mathrm{cr}}$ as representative models of the turbulence-dominated and magnetic-pressure-dominated layers, respectively (Section 2.2). We name a model by attaching “$\Theta_{\mathrm{cr}}$” in front of a value of $\theta/\theta_{\mathrm{cr}}$ followed by “V” in front of a value of $V_{0}/1~{}\mathrm{km~{}s^{-1}}$ (Table 1).

It is useful to estimate global mass-to-flux ratios $\mu_{\mathrm{layer}}$ of the post-shock layers (Strittmatter, 1966; Mouschovias & Spitzer, 1976; Nakano & Nakamura, 1978; Tomisaka et al., 1988; Tomisaka, 2014). The total mass of the post-shock layer increases with time, following $M_{\mathrm{layer}}=2\langle\rho_{0}\rangle V_{0}L^{2}t$, where $L=20~{}\mathrm{pc}$ is the box size along the $y$- and $z-$axes. The magnetic flux held in the post-shock layer depends on direction. In the plane perpendicular to the $x$-axis ($y$-axis), the magnetic fluxes $\Phi_{x}$ ($\Phi_{y}$) is expressed as $B_{0}L^{2}\cos\theta$ ($2B_{0}V_{0}Lt\sin\theta$). The global mass-to-flux ratio is estimated as $\mu_{\mathrm{layer}}=\min(M_{\mathrm{layer}}/\Phi_{x},M_{\mathrm{layer}}/\Phi_{y})$. It increases with time as $M_{\mathrm{layer}}/\Phi_{x}\propto t$, but does not exceed $M_{\mathrm{layer}}/\Phi_{y}=\langle\rho_{0}\rangle L/(B_{0}\sin\theta)$. The global gravitational instability of the post-shock layers occurs if $\mu_{\mathrm{layer}}$ exceeds a critical value of $\mu_{\mathrm{cr}}=(2\pi\sqrt{G})^{-1}$. The ratio is given by

where $\langle n_{0}\rangle_{\mathrm{1}}=\langle n_{0}\rangle/10~{}\mathrm{cm}^{-3}$, $B_{3}=B_{0}/3~{}\mu\mathrm{G}$, $V_{5}=V_{0}/5~{}$km s^-1, $L_{20}=L/20~{}$pc, and $t_{10}=t/10~{}\mathrm{Myr}$. The mass-to-flux ratios of the post-shock layers increase in proportion to time and exceed unity at $t=2.5V_{5}^{-1}~{}\mathrm{Myr}$. They continue to increase at least until $t=10~{}$Myr for the $\theta=0.25\theta_{\mathrm{cr}}$ models while they reach the maximum values of $2$ at $t=5$ Myr for model $\Theta_{\mathrm{cr}}2$V5 and $4$ at $t=2.5~{}$Myr for model $\Theta_{\mathrm{cr}}2$V10. The star formation will be suppressed for largely-inclined compressions with $\sin\theta>0.8\langle n_{0}\rangle_{10}B_{3}L_{20}$ because $M_{\mathrm{layer}}/\Phi_{y}<\mu_{\mathrm{cr}}$ unless magnetic dissipation processes, such as ambipolar diffusion and/or turbulent reconnection work.

The early time evolution of the post-shock layers exhibits similar behaviors found in Paper I. Figure 1 shows the column densities integrated along the $z$-axis at $t=5$ Myr. The accretion of the CNM clumps through the shock fronts disturbs the post-shock layers significantly. For $\theta=0.25\theta_{\mathrm{cr}}$, the CNM clumps are not decelerated by the Lorentz force significantly after passing through the shock fronts because the compression is almost parallel to the magnetic field (Figures 1a and 1c). Figure 2a shows the time evolution of the mean (volume-averaged) kinetic, magnetic, and gravitational energy densities, ${E}_{\mathrm{kin}}$, ${E}_{\mathrm{mag}}$, ${E}_{\mathrm{grv}}$ of the post-shock layer¹¹1 In a similar way as Paper I, the post-shock layer is defined as the region confined by the two shock fronts whose positions correspond to the minimum and maximum of the $x$ coordinates where $P(x,y,z)\geq P_{\mathrm{th}}$ is satisfied, where $P_{\mathrm{th}}$ is a threshold pressure. In this paper, we adopt a value of $P_{\mathrm{th}}=5.8\times 10^{3}~{}$K cm^-3. for $\theta=0.25\theta_{\mathrm{cr}}$. We first focus on the results with $V_{0}=5~{}\mathrm{km~{}s^{-1}}$. The kinetic energies are always larger than the magnetic energies, indicating that the magnetic field lines are easily bent by the gas motion. This can be seen in the magnetic field directions that are random rather than organized as shown in Figures 1a and 1c.

An increase in $V_{0}$ from $5~{}\mathrm{km~{}s^{-1}}$ to $10~{}\mathrm{km~{}s^{-1}}$ does not change the physical properties of the post-shock layers for $\theta=0.25\theta_{\mathrm{cr}}$. Figure 2a clearly shows that all the energy densities are roughly proportional to $V_{0}^{2}$, keeping the relative ratios constant. This can be understood by the fact that the dynamics in the post-shock layers are driven by the ram pressure that is proportional to $V_{0}^{2}$.

Figure 2c shows that $\delta v_{x,y,z}/V_{0}$ does not depend on $V_{0}$ significantly, indicating that the velocity dispersions in the post-shock layers are proportional to the collision speed. For both the models, the degree of anisotropy of the velocity dispersion, which is defined as $\delta v_{x}/\sqrt{(\delta v_{y}^{2}+\delta v_{z}^{2})/2}$, is as large as $\sim 3$ for $V_{0}=5~{}\mathrm{km~{}s^{-1}}$ and $\sim 4$ for $V_{0}=10~{}\mathrm{km~{}s^{-1}}$ at $t=5$ Myr. In addition, $\delta v_{z}$ is almost identical to $\delta v_{y}$ for each model. This suggests that there is no preferential directions in the $(y,z)$-plane.

When $\theta$ is increased from $0.25\theta_{\mathrm{cr}}$ to $2\theta_{\mathrm{cr}}$ for $V_{0}=5~{}\mathrm{km~{}s^{-1}}$ (model $\Theta_{\mathrm{cr}}2\mathrm{V5}$), the physical properties of the post-shock layer drastically changes. Comparison between Figures 1a and 1b shows that the post-shock layer is thinner and denser for $\theta=2\theta_{\mathrm{cr}}$ because the Lorentz force owing to the shock-amplified magnetic field decelerates the CNM clumps efficiently for $\theta=2\theta_{\mathrm{cr}}$. The shock compression amplifies the magnetic energy that dominates over the kinetic energy at $t>2~{}$Myr (Figure 2b). Figure 1b shows that the magnetic field directions in the post-shock layer are parallel to the $y$-axis. At the early phase ($t<0.5~{}$Myr), Figures 2c and 2d show that $\delta v_{x}$ for model $\Theta_{\mathrm{cr}}2\mathrm{V5}$ is as large as that for model $\Theta_{\mathrm{cr}}0.25\mathrm{V5}$ because the magnetic field has not been amplified yet. In the subsequent evolution, $\delta v_{x}$ decreases in proportion to $t^{-1/2}$ (Paper I), and all the three components $\delta v_{x,y,z}$ become of similar magnitudes around $t\sim 5$ Myr.

Unlike model $\Theta_{\mathrm{cr}}0.25\mathrm{V5}$, there is anisotropy in the velocity dispersion in the $(y,z)$ plane for model $\Theta_{\mathrm{cr}}2\mathrm{V5}$. The velocity dispersion along the $y$-axis increases with time while $\delta v_{z}$ is constant with time. The continuous increase in $\delta v_{y}$ corresponds to the formation of filamentary structure by self-gravity, which will be shown in Figure 4.

Unlike the models with $\theta=0.25\theta_{\mathrm{cr}}$, the time evolution of the physical quantities depends on the collision speed for the $\theta=2\theta_{\mathrm{cr}}$ models. The relative importance of $E_{\mathrm{mag}}$ to the total energy is more significant for $V_{0}=10~{}\mathrm{km~{}s^{-1}}$ than for $V_{0}=5~{}\mathrm{km~{}s^{-1}}$. After the early evolution ($t<1~{}\mathrm{Myr}$) in which the velocity dispersions $\delta v_{x,y,z}$ increase in proportion to $V_{0}$, $\delta v_{x}$ begins to decrease earlier for $V_{0}=10~{}\mathrm{km~{}s^{-1}}$ than for $V_{0}=5~{}\mathrm{km~{}s^{-1}}$. This behavior is consistent with the results of Paper I where they found that $\delta v_{x}/V_{0}=0.3\tau^{-1/2}$, where $\tau=t\left(\langle n_{0}\rangle/5~{}\mathrm{cm}^{-3}\right)\left(V_{0}/20~{}\mathrm{km~{}s^{-1}}\right)/(1~{}\mathrm{Myr})$ ²²2 The results are better reproduced by multiplying $\delta v_{x}/V_{0}=0.3\tau^{-1/2}$ by 1.5. This difference comes from the fact that it originally has a relatively large dispersion as shown in the bottom panel of Figure 17 in Paper I. The difference of cooling/heating processes between this paper and Paper I can be another reason. The formula suggests that $\delta v_{x}/V_{0}$ is independent of $\langle\rho_{0}\rangle$ and $V_{0}$ if $\delta v_{x}$ is measured at the time of reaching the same column density. The effect of the rapid decrease in $\delta v_{x}$ can be seen in the fact that the shock fronts are flatter for $V_{0}=10~{}\mathrm{km~{}s^{-1}}$ (Figure 1b and 1d). The anisotropy of the velocity dispersions in the ($y,z$) plane exists also for the $V_{0}=10~{}\mathrm{km~{}s^{-1}}$. Interestingly, the epochs ($t\sim 5~{}\mathrm{Myr}$) when $\delta v_{x}$ becomes equal to $\delta v_{y}$ appear to be independent of $V_{0}$.

There are two main mechanisms to drive the post-shock turbulence. One is the accretion of the upstream CNN clumps as mentioned above. Since the momentum flux of the CNM clumps is much larger than the averaged post-shock momentum flux $\langle n_{0}\rangle V_{0}^{2}$, the accretion of the CNM clumps significantly disturbs the post-shock layers (Paper I). The other arises from deformation of the shock fronts that are generated by the interaction between the shock fronts and inhomogeneous upstream gas (Kobayashi et al., 2020). The post-shock velocity can be super-sonic if the shock normal is tilted more than $\sim 30$ degrees (Landau & Lifshitz, 1959).

The effect of the thermal instability on the post-shock turbulence is minor. Kobayashi et al. (2020) investigated how the post-shock turbulence depends on the amplitude of initial density perturbations. They found that the thermal instability contributes to drive turbulence only when $\delta n_{0}/\langle n_{0}\rangle$ is less than 10%, where $\delta n_{0}$ is the upstream density perturbation. In such situations, the energy conversion efficiency $\epsilon$ from the upstream kinetic energy $\langle\rho_{0}\rangle V_{0}^{2}/2$ to the post-shock turbulence energy is as low as a few percents. When $\delta n_{0}/n_{0}>10~{}\%$, the interaction between shock fronts and density inhomogeneity drives strong turbulence where $\epsilon$ is as high as 10%, which is roughly consistent with the values obtained in this paper, which is $\epsilon\sim(\delta v_{x}/V_{0})^{2}\sim 16\%$ (Paper I).

Figures 2a and 2b show the time evolution of the gravitational energy densities of the post-shock layers that are estimated by

where the integration is done over the post-shock layer³³3 We here estimate the gravitational energy due to the $x$-component of the self-gravitational force only because the periodic boundary conditions are imposed in $y$ and $z$. . The continuous accumulation of the atomic gas deepens the gravitational potential well, and $-E_{\mathrm{grv}}$ increases with time. The total energies of the post-shock layers however are positive, indicating that the post-shock layers remain unbound until at least the end of the simulations.

Assuming that the gas density is uniform inside the post-shock layer, one can derive the gravitational energy density analytically as follows:

where $\Sigma_{0}(t)=2\langle\rho_{0}\rangle V_{0}t$ is the mean column density of the post-shock layer. Figures 2a and 2b show that the time evolution of ${E}_{\mathrm{grv}}$ is consistent with that of ${E}_{\mathrm{grv,ana}}$ for all the models, indicating that the gravitational energy densities increase in proportion to $V_{0}^{2}$. The deviations of $E_{\mathrm{grv}}$ from $E_{\mathrm{grv,ana}}$ come from the non-uniform density distributions.

As an indicator of star formation, we measure the mass fraction of the dense gas whose density is larger than $n_{\mathrm{max}}$ above which the cooling/heating processes are switched off to prevent gravitational collapse (Section 2.3).

The time evolution of the mass fraction of the dense gas with $n>n_{\mathrm{max}}$ does not depend on $V_{0}$ significantly for each $\theta$ (Figure 3). This behavior can be understood by the fact that the ratios of the thermal, kinetic, and magnetic energies to the gravitational energy are independent of $V_{0}$ (Figure 2). We expect that an increase of $\langle n_{0}\rangle$ accelerates the formation of the dense gas with $n>n_{\mathrm{max}}$ because the gravitational energy density is proportional to $\langle n_{0}\rangle^{2}$ while the ram pressure of the pre-shock gas is proportional to $\langle n_{0}\rangle$.

As is expected in Paper I, the formation of the dense gas with $n>n_{\mathrm{max}}$ is delayed by compressions with smaller $\theta$ because anisotropic super-Alfénic turbulence is driven for $\theta$ less than $\theta_{\mathrm{cr}}$.

Another interesting feature is that the difference in $V_{0}$ affects the formation of the dense gas with $n>n_{\mathrm{max}}$ in a different way at different $\theta$. As $V_{0}$ increases, the time evolution of the dense gas fraction does not change significantly for $\theta=2\theta_{\mathrm{cr}}$ while the formation of the dense gas is suppressed slightly for $\theta=0.25\theta_{\mathrm{cr}}$. This is due to the fact that the $\theta$-dependence of the velocity field differs depending on the collision velocity. For $\theta=0.25\theta_{\mathrm{cr}}$, the anisotropy of the velocity dispersion becomes more significant with increasing $V_{0}$ (Figure 2c).

We define $t_{\mathrm{first}}$, the time when the first bound clump is formed, as the earliest time at which the total energy of a clump with $n>n_{\mathrm{max}}$ becomes negative, where the clumps are identified by a friends-of-friends method. Table 1 lists $t_{\mathrm{first}}$, which are almost identical to the epochs when the mass with $n>n_{\mathrm{max}}$ increases rapidly as shown in Figure 3.

The column density distributions along the collision direction strongly depend on $\theta$ as shown in Figure 4. For the models with $\theta=2\theta_{\mathrm{cr}}$, prominent filamentary structures develop, and their elongations are roughly perpendicular to the mean magnetic field (Figures 4b and 4d). This is because the gas is accumulated by the self-gravity preferentially along with the magnetic fields. By contrast, the column density maps for the models with $\theta=0.25\theta_{\mathrm{cr}}$ do not show filamentary structures clearly. The projected magnetic fields are randomized. The detailed structure of filamentary clouds is investigated in forthcoming papers.

The turbulence structures in the post-shock layers are reflected in the field strength-density ($B$-$n$) relation. Figures 5a and 5b show the mean field strengths as a function of density at $t_{\mathrm{first}}+0.5~{}\mathrm{Myr}$ for the models with $V_{0}=5~{}\mathrm{km~{}s^{-1}}$ and $V_{0}=10~{}\mathrm{km~{}s^{-1}}$, respectively.

We first consider the models with $V_{0}=5~{}\mathrm{km~{}s^{-1}}$ (Figure 5a). In low densities less than $10^{3}~{}\mathrm{cm}^{-3}$, model $\Theta_{\mathrm{cr}}2$V5 shows that the mean field strength is determined by the shock compression that amplifies the field strength up to

which is derived by balance between the pre-shock ram pressure and post-shock magnetic pressure. The independence of $B\sim B_{\mathrm{sh}}$ on $n$ for model $\Theta_{\mathrm{cr}}2$V5 reflects gas condensation along the magnetic field, corresponding to the formation of filaments. By contrast, for model $\Theta_{\mathrm{cr}}0.25$V5, the mean field strength ($\sim 6~{}\mu$G) is amplified by shear motion owing to the anisotropic super-Alfvénic turbulence. Assuming that the energy transfer efficiency from the kinetic energy to the magnetic energy $\epsilon$ is 20%, the following estimated field strength $B_{\mathrm{tur}}$ is consistent with the simulation results,

Once the gas density becomes larger than $\sim 10^{3}~{}\mathrm{cm}^{-3}$, the $\theta$-dependence of the mean field strengths disappears, and the mean field strengths begin to increase with the gas density, following a similar $B$-$n$ relation. Interestingly, the mean field strengths in the high density range are consistent with the field strength $B_{\beta=1}$,

which is given by $\beta=1$ at $T=20~{}$K that is a typical temperature of the dense regions, where $\beta$ is the plasma beta. This is where the self-gravity in the clumps becomes significant and the field amplification mode changes.

The $B$-$n$ relations of the models with $V_{0}=10~{}\mathrm{km~{}s^{-1}}$ behave similarly to those of the models with $V_{0}=5~{}\mathrm{km~{}s^{-1}}$. At each $\theta/\theta_{\mathrm{cr}}$, the field strength increases by a factor of two owing to the increase in $V_{0}$. This is consistent with the predictions from Equations (11) and (12). Similar to the models with $V_{0}=5~{}\mathrm{km~{}s^{-1}}$, in the high density regions, the field strengths follow $B_{\beta=1}$ although the density range is limited. The critical density where the transition between the low and high density regimes occurs is shifted toward high densities when $V_{0}$ increases from $5~{}\mathrm{km~{}s^{-1}}$ to $10~{}\mathrm{km~{}s^{-1}}$.

For the models with $\theta=2\theta_{\mathrm{cr}}$, the mean field strengths take a constant value of $B_{\mathrm{sh}}$ for lower densities and it increases following $B_{\beta=1}\propto n^{1/2}$ for higher densities (Figure 5b). From these facts, a critical density $n_{\mathrm{grv}}$ can be derived by equating $B_{\mathrm{sh}}$ and $B_{\beta=1}$ as follows:

When the density exceeds $n_{\mathrm{grv}}$, the self-gravitational force amplifies the magnetic fields.

Figure 5c is the same as Figures 5a and 5b but for the gas density and field strength are normalized by $n_{\mathrm{grv}}$ and $B_{\mathrm{sh}}$, respectively. The $B$-$n$ relation is characterized by $n_{\mathrm{grv}}$ and $B_{\mathrm{sh}}$ reasonably well for the models with $\theta=2\theta_{\mathrm{cr}}$. The mean field strength for the models is well approximated by

where $\beta_{\mathrm{d}}$ is a typical plasma $\beta$ for higher densities, and set to $1$ in Figure 5. Equation (15) is similar to that found by Tomisaka et al. (1988) who investigate the magnetohydrostatic structure of axisymmetric clouds.

From the virial theorem without the surface terms (Equation (19)), a stability criterion can be derived by $\ddot{I}<0$, which is rewritten as

The virial parameter $\alpha_{\mathrm{tot}}$ is related to that often used in observations as a diagnose of the stability of dense clumps and cores (Bertoldi & McKee, 1992), but it contains the effect of support due to magnetic fields.

Before showing the dependence of $\alpha_{\mathrm{tot}}$ on $n_{\mathrm{th}}$, $M_{\mathrm{cl}}$, $V_{0}$, and $\theta$, we investigate the velocity structure of the dense clumps that are related to ${\cal E}_{\mathrm{kin}}$ in Sections 4.1.1 and 4.1.2.